87 {
88
92
96 switch (theParticle->
getType()) {
98 if (drnd < 0.3972) {
99
100 nbpart=2;
103 }
104 else if (drnd < 0.7265) {
105
109 }
110 else if (drnd < 0.9575) {
111
115 }
116 else {
117
121 }
122 break;
124 if (drnd < 0.9009) {
125
129 }
130 else if (drnd < 0.9845) {
131
132 nbpart=2;
135 }
136 else {
137
138 nbpart=2;
141 }
142 break;
143 default:
144 INCL_FATAL(
"Unrecognized pion resonance type; type=" << theParticle->
getType() <<
'\n');
145 break;
146 }
147
148 if (nbpart == 2) {
150 sampleAngles(&ctet, &stet, &fi);
151
155
158 if (beta >= 1.0e-10)
159 sal = incidentDirection.
perp()/beta;
160 if (sal >= 1.0e-6) {
165 G4double t1 = ctet+cal*stet*sfi/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
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);
188
189 fs->addModifiedParticle(theParticle);
190 fs->addCreatedParticle(createdParticle);
191
192 }
193 else if (nbpart == 3) {
194
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
210 }
211
212 }
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 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.