Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCL::PionResonanceDecayChannel Class Reference

#include <G4INCLPionResonanceDecayChannel.hh>

+ Inheritance diagram for G4INCL::PionResonanceDecayChannel:

Public Member Functions

 PionResonanceDecayChannel (Particle *, ThreeVector const &)
 
virtual ~PionResonanceDecayChannel ()
 
void fillFinalState (FinalState *fs)
 
- Public Member Functions inherited from G4INCL::IChannel
 IChannel ()
 
virtual ~IChannel ()
 
FinalStategetFinalState ()
 
virtual void fillFinalState (FinalState *fs)=0
 

Static Public Member Functions

static G4double computeDecayTime (Particle *p)
 

Detailed Description

Definition at line 49 of file G4INCLPionResonanceDecayChannel.hh.

Constructor & Destructor Documentation

◆ PionResonanceDecayChannel()

G4INCL::PionResonanceDecayChannel::PionResonanceDecayChannel ( Particle p,
ThreeVector const &  dir 
)

Definition at line 51 of file G4INCLPionResonanceDecayChannel.cc.

52 :theParticle(p), incidentDirection(dir)
53 { }

◆ ~PionResonanceDecayChannel()

G4INCL::PionResonanceDecayChannel::~PionResonanceDecayChannel ( )
virtual

Definition at line 55 of file G4INCLPionResonanceDecayChannel.cc.

55{}

Member Function Documentation

◆ computeDecayTime()

G4double G4INCL::PionResonanceDecayChannel::computeDecayTime ( Particle p)
static

Definition at line 58 of file G4INCLPionResonanceDecayChannel.cc.

58 {
59 const G4double m = p->getMass();
60 const G4double geff = p->getEnergy()/m;
61// const G4double geta = 1.31e-3;
62 const G4double gomega = 8.49;
63 G4double gg=0.;
64 switch (p->getType()) {
65/* case Eta:
66 gg=geta;
67 break;*/
68 case Omega:
69 gg=gomega;
70 break;
71 default:
72 INCL_FATAL("Unrecognized pion resonance type; type=" << p->getType() << '\n');
73 break;
74 }
75 const G4double tpires = -G4INCL::PhysicalConstants::hc/gg*std::log(Random::shoot())*geff;
76 return tpires;
77 }
#define INCL_FATAL(x)
double G4double
Definition: G4Types.hh:83
const G4double hc
[MeV*fm]
G4double shoot()
Definition: G4INCLRandom.cc:93

Referenced by G4INCL::StandardPropagationModel::generateDecays().

◆ fillFinalState()

void G4INCL::PionResonanceDecayChannel::fillFinalState ( FinalState fs)
virtual

Implements G4INCL::IChannel.

Definition at line 87 of file G4INCLPionResonanceDecayChannel.cc.

87 {
88
89 ParticleType createdType;
90 ParticleType pionType1=Neutron; // to avoid forgetting pionType definition when 3 particles are emitted
91 ParticleType pionType2=Neutron;
92
93 const G4double sqrtS = theParticle->getMass();
94 G4int nbpart = 3; // number of emitted particles
96 switch (theParticle->getType()) {
97 case Eta:
98 if (drnd < 0.3972) { // renormalized to the only four decays taken into account here
99// 2 photons
100 nbpart=2;
101 theParticle->setType(Photon);
102 createdType = Photon;
103 }
104 else if (drnd < 0.7265) {
105// 3 pi0
106 theParticle->setType(PiZero);
107 pionType1 = PiZero;
108 pionType2 = PiZero;
109 }
110 else if (drnd < 0.9575) {
111// pi+ pi- pi0
112 theParticle->setType(PiZero);
113 pionType1 = PiPlus;
114 pionType2 = PiMinus;
115 }
116 else {
117// pi+ pi- photon
118 theParticle->setType(Photon);
119 pionType1 = PiPlus;
120 pionType2 = PiMinus;
121 }
122 break;
123 case Omega:
124 if (drnd < 0.9009) { // renormalized to the only three decays taken into account here
125// pi+ pi- pi0
126 theParticle->setType(PiZero);
127 pionType1 = PiPlus;
128 pionType2 = PiMinus;
129 }
130 else if (drnd < 0.9845) {
131// pi0 photon
132 nbpart=2;
133 theParticle->setType(PiZero);
134 createdType = Photon;
135 }
136 else {
137// pi+ pi-
138 nbpart=2;
139 theParticle->setType(PiPlus);
140 createdType = PiMinus;
141 }
142 break;
143 default:
144 INCL_FATAL("Unrecognized pion resonance type; type=" << theParticle->getType() << '\n');
145 break;
146 }
147
148 if (nbpart == 2) {
149 G4double fi, ctet, stet;
150 sampleAngles(&ctet, &stet, &fi);
151
152 G4double cfi = std::cos(fi);
153 G4double sfi = std::sin(fi);
154 G4double beta = incidentDirection.mag();
155
156 G4double q1, q2, q3;
157 G4double sal=0.0;
158 if (beta >= 1.0e-10)
159 sal = incidentDirection.perp()/beta;
160 if (sal >= 1.0e-6) {
161 G4double b1 = incidentDirection.getX();
162 G4double b2 = incidentDirection.getY();
163 G4double b3 = incidentDirection.getZ();
164 G4double cal = b3/beta;
165 G4double t1 = ctet+cal*stet*sfi/sal;
166 G4double t2 = stet/sal;
167 q1=(b1*t1+b2*t2*cfi)/beta;
168 q2=(b2*t1-b1*t2*cfi)/beta;
169 q3=(b3*t1/beta-t2*sfi);
170 } else {
171 q1 = stet*cfi;
172 q2 = stet*sfi;
173 q3 = ctet;
174 }
175
177 theParticle->getMass(),
178 ParticleTable::getINCLMass(createdType));
179 q1 *= xq;
180 q2 *= xq;
181 q3 *= xq;
182
183 ThreeVector createdMomentum(q1, q2, q3);
184 ThreeVector createdPosition(theParticle->getPosition());
185 Particle *createdParticle = new Particle(createdType, createdMomentum, createdPosition);
186 theParticle->setMomentum(-createdMomentum);
187 theParticle->adjustEnergyFromMomentum();
188
189 fs->addModifiedParticle(theParticle);
190 fs->addCreatedParticle(createdParticle);
191
192 }
193 else if (nbpart == 3) {
194// assert(pionType1!=Neutron && pionType2!=Neutron);
195 ParticleList list;
196 list.push_back(theParticle);
197 const ThreeVector &rposdecay = theParticle->getPosition();
198 const ThreeVector zero;
199 Particle *Pion1 = new Particle(pionType1,zero,rposdecay);
200 Particle *Pion2 = new Particle(pionType2,zero,rposdecay);
201 list.push_back(Pion1);
202 list.push_back(Pion2);
203
204 fs->addModifiedParticle(theParticle);
205 fs->addCreatedParticle(Pion1);
206 fs->addCreatedParticle(Pion2);
207
208 // PhaseSpaceGenerator::generateBiased(sqrtS, list, 0, angularSlope); Biasing?
210 }
211
212 }
int G4int
Definition: G4Types.hh:85
const G4INCL::ThreeVector & getPosition() const
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
G4INCL::ParticleType getType() const
void setType(ParticleType t)
G4double getMass() const
Get the cached particle mass.
G4double getY() const
G4double getZ() const
G4double mag() const
G4double perp() const
G4double getX() const
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
G4double getINCLMass(const G4int A, const G4int Z, const G4int S)
Get INCL nuclear mass (in MeV/c^2)
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.

The documentation for this class was generated from the following files: