Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4ReactionDynamics Class Reference

#include <G4ReactionDynamics.hh>

Public Member Functions

 G4ReactionDynamics ()
 
virtual ~G4ReactionDynamics ()
 
virtual G4double FindInelasticity ()
 
virtual G4double FindTimeDelay ()
 
G4bool GenerateXandPt (G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4DynamicParticle *originalTarget, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool leadFlag, G4ReactionProduct &leadingStrangeParticle)
 
void SuppressChargedPions (G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged)
 
G4bool TwoCluster (G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4DynamicParticle *originalTarget, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool leadFlag, G4ReactionProduct &leadingStrangeParticle)
 
void TwoBody (G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4DynamicParticle *originalTarget, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &targetHasChanged)
 
G4int Factorial (G4int n)
 
G4double GenerateNBodyEvent (const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen)
 
void ProduceStrangeParticlePairs (G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, const G4ReactionProduct &modifiedOriginal, const G4DynamicParticle *originalTarget, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged)
 
void NuclearReaction (G4FastVector< G4ReactionProduct, 4 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4Nucleus &aNucleus, const G4double theAtomicMass, const G4double *massVec)
 

Detailed Description

Definition at line 48 of file G4ReactionDynamics.hh.

Constructor & Destructor Documentation

◆ G4ReactionDynamics()

G4ReactionDynamics::G4ReactionDynamics ( )
inline

Definition at line 52 of file G4ReactionDynamics.hh.

52{}

◆ ~G4ReactionDynamics()

virtual G4ReactionDynamics::~G4ReactionDynamics ( )
inlinevirtual

Definition at line 54 of file G4ReactionDynamics.hh.

54{}

Member Function Documentation

◆ Factorial()

G4int G4ReactionDynamics::Factorial ( G4int  n)

Definition at line 2743 of file G4ReactionDynamics.cc.

2744 { // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
2745 G4int mn = std::min(n,10);
2746 G4int result = 1;
2747 if( mn <= 1 )return result;
2748 for( G4int i=2; i<=mn; ++i )result *= i;
2749 return result;
2750 }
G4int
int G4int
Definition: G4Types.hh:66

◆ FindInelasticity()

virtual G4double G4ReactionDynamics::FindInelasticity ( )
inlinevirtual

Definition at line 56 of file G4ReactionDynamics.hh.

57 { return 0.0; }

◆ FindTimeDelay()

virtual G4double G4ReactionDynamics::FindTimeDelay ( )
inlinevirtual

Definition at line 59 of file G4ReactionDynamics.hh.

60 { return 0.0; }

◆ GenerateNBodyEvent()

G4double G4ReactionDynamics::GenerateNBodyEvent ( const G4double  totalEnergy,
const G4bool  constantCrossSection,
G4FastVector< G4ReactionProduct, GHADLISTSIZE > &  vec,
G4int &  vecLen 
)

Definition at line 2491 of file G4ReactionDynamics.cc.

2496 {
2497 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2498 // derived from original FORTRAN code PHASP by H. Fesefeldt (02-Dec-1986)
2499 // Returns the weight of the event
2500 //
2501 G4int i;
2502 const G4double expxu = 82.; // upper bound for arg. of exp
2503 const G4double expxl = -expxu; // lower bound for arg. of exp
2504 if( vecLen < 2 )
2505 {
2506 G4cerr << "*** Error in G4ReactionDynamics::GenerateNBodyEvent" << G4endl;
2507 G4cerr << " number of particles < 2" << G4endl;
2508 G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
2509 return -1.0;
2510 }
2511 G4double mass[18]; // mass of each particle
2512 G4double energy[18]; // total energy of each particle
2513 G4double pcm[3][18]; // pcm is an array with 3 rows and vecLen columns
2514
2515 G4double totalMass = 0.0;
2516 G4double extraMass = 0;
2517 G4double sm[18];
2518
2519 for( i=0; i<vecLen; ++i )
2520 {
2521 mass[i] = vec[i]->GetMass()/GeV;
2522 if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
2523 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
2524 pcm[0][i] = 0.0; // x-momentum of i-th particle
2525 pcm[1][i] = 0.0; // y-momentum of i-th particle
2526 pcm[2][i] = 0.0; // z-momentum of i-th particle
2527 energy[i] = mass[i]; // total energy of i-th particle
2528 totalMass += mass[i];
2529 sm[i] = totalMass;
2530 }
2531 G4double totalE = totalEnergy/GeV;
2532 if( totalMass > totalE )
2533 {
2534 //G4cerr << "*** Error in G4ReactionDynamics::GenerateNBodyEvent" << G4endl;
2535 //G4cerr << " total mass (" << totalMass*GeV << "MeV) > total energy ("
2536 // << totalEnergy << "MeV)" << G4endl;
2537 totalE = totalMass;
2538 return -1.0;
2539 }
2540 G4double kineticEnergy = totalE - totalMass;
2541 G4double emm[18];
2542 //G4double *emm = new G4double [vecLen];
2543 emm[0] = mass[0];
2544 emm[vecLen-1] = totalE;
2545 if( vecLen > 2 ) // the random numbers are sorted
2546 {
2547 G4double ran[18];
2548 for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
2549 for( i=0; i<vecLen-2; ++i )
2550 {
2551 for( G4int j=vecLen-2; j>i; --j )
2552 {
2553 if( ran[i] > ran[j] )
2554 {
2555 G4double temp = ran[i];
2556 ran[i] = ran[j];
2557 ran[j] = temp;
2558 }
2559 }
2560 }
2561 for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
2562 }
2563 // Weight is the sum of logarithms of terms instead of the product of terms
2564 G4bool lzero = true;
2565 G4double wtmax = 0.0;
2566 if( constantCrossSection ) // this is KGENEV=1 in PHASP
2567 {
2568 G4double emmax = kineticEnergy + mass[0];
2569 G4double emmin = 0.0;
2570 for( i=1; i<vecLen; ++i )
2571 {
2572 emmin += mass[i-1];
2573 emmax += mass[i];
2574 G4double wtfc = 0.0;
2575 if( emmax*emmax > 0.0 )
2576 {
2577 G4double arg = emmax*emmax
2578 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
2579 - 2.0*(emmin*emmin+mass[i]*mass[i]);
2580 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
2581 }
2582 if( wtfc == 0.0 )
2583 {
2584 lzero = false;
2585 break;
2586 }
2587 wtmax += std::log( wtfc );
2588 }
2589 if( lzero )
2590 wtmax = -wtmax;
2591 else
2592 wtmax = expxu;
2593 }
2594 else
2595 {
2596 // ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
2597 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
2598 256.3704, 268.4705, 240.9780, 189.2637,
2599 132.1308, 83.0202, 47.4210, 24.8295,
2600 12.0006, 5.3858, 2.2560, 0.8859 };
2601 wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
2602 }
2603 lzero = true;
2604 G4double pd[50];
2605 //G4double *pd = new G4double [vecLen-1];
2606 for( i=0; i<vecLen-1; ++i )
2607 {
2608 pd[i] = 0.0;
2609 if( emm[i+1]*emm[i+1] > 0.0 )
2610 {
2611 G4double arg = emm[i+1]*emm[i+1]
2612 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
2613 /(emm[i+1]*emm[i+1])
2614 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
2615 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
2616 }
2617 if( pd[i] <= 0.0 ) // changed from == on 02 April 98
2618 lzero = false;
2619 else
2620 wtmax += std::log( pd[i] );
2621 }
2622 G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
2623 if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
2624
2625 G4double bang, cb, sb, s0, s1, s2, c, ss, esys, a, b, gama, beta;
2626 pcm[0][0] = 0.0;
2627 pcm[1][0] = pd[0];
2628 pcm[2][0] = 0.0;
2629 for( i=1; i<vecLen; ++i )
2630 {
2631 pcm[0][i] = 0.0;
2632 pcm[1][i] = -pd[i-1];
2633 pcm[2][i] = 0.0;
2634 bang = twopi*G4UniformRand();
2635 cb = std::cos(bang);
2636 sb = std::sin(bang);
2637 c = 2.0*G4UniformRand() - 1.0;
2638 ss = std::sqrt( std::fabs( 1.0-c*c ) );
2639 if( i < vecLen-1 )
2640 {
2641 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
2642 beta = pd[i]/esys;
2643 gama = esys/emm[i];
2644 for( G4int j=0; j<=i; ++j )
2645 {
2646 s0 = pcm[0][j];
2647 s1 = pcm[1][j];
2648 s2 = pcm[2][j];
2649 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
2650 a = s0*c - s1*ss; // rotation
2651 pcm[1][j] = s0*ss + s1*c;
2652 b = pcm[2][j];
2653 pcm[0][j] = a*cb - b*sb;
2654 pcm[2][j] = a*sb + b*cb;
2655 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
2656 }
2657 }
2658 else
2659 {
2660 for( G4int j=0; j<=i; ++j )
2661 {
2662 s0 = pcm[0][j];
2663 s1 = pcm[1][j];
2664 s2 = pcm[2][j];
2665 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
2666 a = s0*c - s1*ss; // rotation
2667 pcm[1][j] = s0*ss + s1*c;
2668 b = pcm[2][j];
2669 pcm[0][j] = a*cb - b*sb;
2670 pcm[2][j] = a*sb + b*cb;
2671 }
2672 }
2673 }
2674 for( i=0; i<vecLen; ++i )
2675 {
2676 vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
2677 vec[i]->SetTotalEnergy( energy[i]*GeV );
2678 }
2679
2680 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2681 return weight;
2682 }
G4double
double G4double
Definition: G4Types.hh:64
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53

Referenced by GenerateXandPt(), NuclearReaction(), and TwoCluster().

◆ GenerateXandPt()

G4bool G4ReactionDynamics::GenerateXandPt ( G4FastVector< G4ReactionProduct, GHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4HadProjectile *  originalIncident,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4DynamicParticle *  originalTarget,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged,
G4bool  leadFlag,
G4ReactionProduct &  leadingStrangeParticle 
)

Definition at line 105 of file G4ReactionDynamics.cc.

118 {
119 //
120 // derived from original FORTRAN code GENXPT by H. Fesefeldt (11-Oct-1987)
121 //
122 // Generation of X- and PT- values for incident, target, and all secondary particles
123 // A simple single variable description E D3S/DP3= F(Q) with
124 // Q^2 = (M*X)^2 + PT^2 is used. Final state kinematic is produced
125 // by an FF-type iterative cascade method
126 //
127 // internal units are GeV
128 //
129 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
130
131 // Protection in case no secondary has been created; cascades down to two-body.
132 if(vecLen == 0) return false;
133
134 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
135 G4ParticleDefinition *aProton = G4Proton::Proton();
136 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
137 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
138 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
139
140 G4int i, l;
141 G4bool veryForward = false;
142
143 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
144 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
145 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
146 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
147 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
148 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
149 targetMass*targetMass +
150 2.0*targetMass*etOriginal ); // GeV
151 G4double currentMass = currentParticle.GetMass()/GeV;
152 targetMass = targetParticle.GetMass()/GeV;
153 //
154 // randomize the order of the secondary particles
155 // note that the current and target particles are not affected
156 //
157 for( i=0; i<vecLen; ++i )
158 {
159 G4int itemp = G4int( G4UniformRand()*vecLen );
160 G4ReactionProduct pTemp = *vec[itemp];
161 *vec[itemp] = *vec[i];
162 *vec[i] = pTemp;
163 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
164 }
165
166 if( currentMass == 0.0 && targetMass == 0.0 )
167 {
168 // Target and projectile have annihilated. Replace them with the first
169 // two secondaries in the list. Current particle KE is maintained.
170
171 G4double ek = currentParticle.GetKineticEnergy();
172 G4ThreeVector mom = currentParticle.GetMomentum();
173 currentParticle = *vec[0];
174 targetParticle = *vec[1];
175 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
176 G4ReactionProduct *temp = vec[vecLen-1];
177 delete temp;
178 temp = vec[vecLen-2];
179 delete temp;
180 vecLen -= 2;
181 currentMass = currentParticle.GetMass()/GeV;
182 targetMass = targetParticle.GetMass()/GeV;
183 incidentHasChanged = true;
184 targetHasChanged = true;
185 currentParticle.SetKineticEnergy( ek );
186 currentParticle.SetMomentum( mom );
187 veryForward = true;
188 }
189 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
190 const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
191 const G4double protonMass = aProton->GetPDGMass()/MeV;
192
193 if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
194 && G4UniformRand() >= 0.7) {
195 G4ReactionProduct temp = currentParticle;
196 currentParticle = targetParticle;
197 targetParticle = temp;
198 incidentHasChanged = true;
199 targetHasChanged = true;
200 currentMass = currentParticle.GetMass()/GeV;
201 targetMass = targetParticle.GetMass()/GeV;
202 }
203 const G4double afc = std::min( 0.75,
204 0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
205 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
206
207 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
208 G4double forwardEnergy = freeEnergy/2.;
209 G4int forwardCount = 1; // number of particles in forward hemisphere
210
211 G4double backwardEnergy = freeEnergy/2.;
212 G4int backwardCount = 1; // number of particles in backward hemisphere
213
214 if(veryForward)
215 {
216 if(currentParticle.GetSide()==-1)
217 {
218 forwardEnergy += currentMass;
219 forwardCount --;
220 backwardEnergy -= currentMass;
221 backwardCount ++;
222 }
223 if(targetParticle.GetSide()!=-1)
224 {
225 backwardEnergy += targetMass;
226 backwardCount --;
227 forwardEnergy -= targetMass;
228 forwardCount ++;
229 }
230 }
231
232 for( i=0; i<vecLen; ++i )
233 {
234 if( vec[i]->GetSide() == -1 )
235 {
236 ++backwardCount;
237 backwardEnergy -= vec[i]->GetMass()/GeV;
238 } else {
239 ++forwardCount;
240 forwardEnergy -= vec[i]->GetMass()/GeV;
241 }
242 }
243 //
244 // Add particles from intranuclear cascade.
245 // nuclearExcitationCount = number of new secondaries produced by nuclear excitation
246 // extraCount = number of nucleons within these new secondaries
247 //
248 G4double xtarg;
249 if( centerofmassEnergy < (2.0+G4UniformRand()) )
250 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
251 else
252 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
253 if( xtarg <= 0.0 )xtarg = 0.01;
254 G4int nuclearExcitationCount = Poisson( xtarg );
255 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
256 G4int extraNucleonCount = 0;
257 G4double extraNucleonMass = 0.0;
258 if( nuclearExcitationCount > 0 )
259 {
260 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
261 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
262 G4int momentumBin = 0;
263 while( (momentumBin < 6) &&
264 (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
265 ++momentumBin;
266 momentumBin = std::min( 5, momentumBin );
267 //
268 // NOTE: in GENXPT, these new particles were given negative codes
269 // here I use NewlyAdded = true instead
270 //
271
272 for( i=0; i<nuclearExcitationCount; ++i )
273 {
274 G4ReactionProduct * pVec = new G4ReactionProduct();
275 if( G4UniformRand() < nucsup[momentumBin] )
276 {
277 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
278 pVec->SetDefinition( aProton );
279 else
280 pVec->SetDefinition( aNeutron );
281 pVec->SetSide( -2 ); // -2 means backside nucleon
282 ++extraNucleonCount;
283 backwardEnergy += pVec->GetMass()/GeV;
284 extraNucleonMass += pVec->GetMass()/GeV;
285 }
286 else
287 {
288 G4double ran = G4UniformRand();
289 if( ran < 0.3181 )
290 pVec->SetDefinition( aPiPlus );
291 else if( ran < 0.6819 )
292 pVec->SetDefinition( aPiZero );
293 else
294 pVec->SetDefinition( aPiMinus );
295 pVec->SetSide( -1 ); // backside particle, but not a nucleon
296 }
297 pVec->SetNewlyAdded( true ); // true is the same as IPA(i)<0
298 vec.SetElement( vecLen++, pVec );
299 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
300 backwardEnergy -= pVec->GetMass()/GeV;
301 ++backwardCount;
302 }
303 }
304 //
305 // assume conservation of kinetic energy in forward & backward hemispheres
306 //
307 G4int is, iskip;
308 while( forwardEnergy <= 0.0 ) // must eliminate a particle from the forward side
309 {
310 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
311 iskip = G4int(G4UniformRand()*forwardCount) + 1; // 1 <= iskip <= forwardCount
312 is = 0;
313 G4int forwardParticlesLeft = 0;
314 for( i=(vecLen-1); i>=0; --i )
315 {
316 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
317 {
318 forwardParticlesLeft = 1;
319 if( ++is == iskip )
320 {
321 forwardEnergy += vec[i]->GetMass()/GeV;
322 for( G4int j=i; j<(vecLen-1); j++ )*vec[j] = *vec[j+1]; // shift up
323 --forwardCount;
324 delete vec[vecLen-1];
325 if( --vecLen == 0 )return false; // all the secondaries have been eliminated
326 break; // --+
327 } // |
328 } // |
329 } // break goes down to here
330 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
331 if( forwardParticlesLeft == 0 )
332 {
333 G4int iremove = -1;
334 for (i = 0; i < vecLen; i++) {
335 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
336 iremove = i;
337 break;
338 }
339 }
340 if (iremove == -1) {
341 for (i = 0; i < vecLen; i++) {
342 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
343 iremove = i;
344 break;
345 }
346 }
347 }
348 if (iremove == -1) iremove = 0;
349
350 forwardEnergy += vec[iremove]->GetMass()/GeV;
351 if (vec[iremove]->GetSide() > 0) --forwardCount;
352
353 for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
354 delete vec[vecLen-1];
355 vecLen--;
356 if (vecLen == 0) return false; // all secondaries have been eliminated
357 break;
358 }
359 } // while
360
361 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
362 while( backwardEnergy <= 0.0 ) // must eliminate a particle from the backward side
363 {
364 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
365 iskip = G4int(G4UniformRand()*backwardCount) + 1; // 1 <= iskip <= backwardCount
366 is = 0;
367 G4int backwardParticlesLeft = 0;
368 for( i=(vecLen-1); i>=0; --i )
369 {
370 if( vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled())
371 {
372 backwardParticlesLeft = 1;
373 if( ++is == iskip ) // eliminate the i'th particle
374 {
375 if( vec[i]->GetSide() == -2 )
376 {
377 --extraNucleonCount;
378 extraNucleonMass -= vec[i]->GetMass()/GeV;
379 backwardEnergy -= vec[i]->GetTotalEnergy()/GeV;
380 }
381 backwardEnergy += vec[i]->GetTotalEnergy()/GeV;
382 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
383 --backwardCount;
384 delete vec[vecLen-1];
385 if( --vecLen == 0 )return false; // all the secondaries have been eliminated
386 break;
387 }
388 }
389 }
390 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
391 if( backwardParticlesLeft == 0 )
392 {
393 G4int iremove = -1;
394 for (i = 0; i < vecLen; i++) {
395 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
396 iremove = i;
397 break;
398 }
399 }
400 if (iremove == -1) {
401 for (i = 0; i < vecLen; i++) {
402 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
403 iremove = i;
404 break;
405 }
406 }
407 }
408 if (iremove == -1) iremove = 0;
409
410 backwardEnergy += vec[iremove]->GetMass()/GeV;
411 if (vec[iremove]->GetSide() > 0) --backwardCount;
412
413 for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
414 delete vec[vecLen-1];
415 vecLen--;
416 if (vecLen == 0) return false; // all secondaries have been eliminated
417 break;
418 }
419 } // while
420
421 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
422 //
423 // define initial state vectors for Lorentz transformations
424 // the pseudoParticles have non-standard masses, hence the "pseudo"
425 //
426 G4ReactionProduct pseudoParticle[10];
427 for( i=0; i<10; ++i )pseudoParticle[i].SetZero();
428
429 pseudoParticle[0].SetMass( mOriginal*GeV );
430 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
431 pseudoParticle[0].SetTotalEnergy(
432 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
433
434 pseudoParticle[1].SetMass( protonMass*MeV ); // this could be targetMass
435 pseudoParticle[1].SetTotalEnergy( protonMass*MeV );
436
437 pseudoParticle[3].SetMass( protonMass*(1+extraNucleonCount)*MeV );
438 pseudoParticle[3].SetTotalEnergy( protonMass*(1+extraNucleonCount)*MeV );
439
440 pseudoParticle[8].SetMomentum( 1.0*GeV, 0.0, 0.0 );
441
442 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
443 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
444
445 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
446 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
447
448 G4double dndl[20];
449 //
450 // main loop for 4-momentum generation
451 // see Pitha-report (Aachen) for a detailed description of the method
452 //
453 G4double aspar, pt, et, x, pp, pp1, rmb, wgt;
454 G4int innerCounter, outerCounter;
455 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
456
457 G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
458 //
459 // process the secondary particles in reverse order
460 // the incident particle is Done after the secondaries
461 // nucleons, including the target, in the backward hemisphere are also Done later
462 //
463 G4double binl[20] = {0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.11,1.25,
464 1.43,1.67,2.0,2.5,3.33,5.00,10.00};
465 G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere
466 G4double totalEnergy, kineticEnergy, vecMass;
467
468 for( i=(vecLen-1); i>=0; --i )
469 {
470 G4double phi = G4UniformRand()*twopi;
471 if( vec[i]->GetNewlyAdded() ) // added from intranuclear cascade
472 {
473 if( vec[i]->GetSide() == -2 ) // is a nucleon
474 {
475 if( backwardNucleonCount < 18 )
476 {
477 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
478 for(G4int j=0; j<vecLen; j++) delete vec[j];
479 vecLen = 0;
480 throw G4HadReentrentException(__FILE__, __LINE__,
481 "G4ReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon");
482 }
483 vec[i]->SetSide( -3 );
484 ++backwardNucleonCount;
485 continue;
486 }
487 }
488 }
489 //
490 // set pt and phi values, they are changed somewhat in the iteration loop
491 // set mass parameter for lambda fragmentation model
492 //
493 vecMass = vec[i]->GetMass()/GeV;
494 G4double ran = -std::log(1.0-G4UniformRand())/3.5;
495 if( vec[i]->GetSide() == -2 ) // backward nucleon
496 {
497 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon" ||
498 vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
499 aspar = 0.75;
500 pt = std::sqrt( std::pow( ran, 1.7 ) );
501 } else { // vec[i] must be a proton, neutron,
502 aspar = 0.20; // lambda, sigma, xsi, or ion
503 pt = std::sqrt( std::pow( ran, 1.2 ) );
504 }
505
506 } else { // not a backward nucleon
507
508 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
509 aspar = 0.75;
510 pt = std::sqrt( std::pow( ran, 1.7 ) );
511 } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
512 aspar = 0.70;
513 pt = std::sqrt( std::pow( ran, 1.7 ) );
514 } else { // vec[i] must be a proton, neutron,
515 aspar = 0.65; // lambda, sigma, xsi, or ion
516 pt = std::sqrt( std::pow( ran, 1.5 ) );
517 }
518 }
519 pt = std::max( 0.001, pt );
520 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
521 for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
522 if( vec[i]->GetSide() > 0 )
523 et = pseudoParticle[0].GetTotalEnergy()/GeV;
524 else
525 et = pseudoParticle[1].GetTotalEnergy()/GeV;
526 dndl[0] = 0.0;
527 //
528 // start of outer iteration loop
529 //
530 outerCounter = 0;
531 eliminateThisParticle = true;
532 resetEnergies = true;
533 while( ++outerCounter < 3 )
534 {
535 for( l=1; l<20; ++l )
536 {
537 x = (binl[l]+binl[l-1])/2.;
538 pt = std::max( 0.001, pt );
539 if( x > 1.0/pt )
540 dndl[l] += dndl[l-1]; // changed from just = on 02 April 98
541 else
542 dndl[l] = et * aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) )
543 * (binl[l]-binl[l-1]) / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass )
544 + dndl[l-1];
545 }
546 innerCounter = 0;
547 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
548 //
549 // start of inner iteration loop
550 //
551 while( ++innerCounter < 7 )
552 {
553 ran = G4UniformRand()*dndl[19];
554 l = 1;
555 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
556 l = std::min( 19, l );
557 x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
558 if( vec[i]->GetSide() < 0 )x *= -1.;
559 vec[i]->SetMomentum( x*et*GeV ); // set the z-momentum
560 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
561 vec[i]->SetTotalEnergy( totalEnergy*GeV );
562 kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
563 if( vec[i]->GetSide() > 0 ) // forward side
564 {
565 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy )
566 {
567 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
568 forwardKinetic += kineticEnergy;
569 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
570 pseudoParticle[6].SetMomentum( 0.0 ); // set the z-momentum
571 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
572 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
573 phi += pi + normal()*pi/12.0;
574 if( phi > twopi )phi -= twopi;
575 if( phi < 0.0 )phi = twopi - phi;
576 outerCounter = 2; // leave outer loop
577 eliminateThisParticle = false; // don't eliminate this particle
578 resetEnergies = false;
579 break; // leave inner loop
580 }
581 if( innerCounter > 5 )break; // leave inner loop
582 if( backwardEnergy >= vecMass ) // switch sides
583 {
584 vec[i]->SetSide( -1 );
585 forwardEnergy += vecMass;
586 backwardEnergy -= vecMass;
587 ++backwardCount;
588 }
589 } else { // backward side
590 if( extraNucleonCount > 19 ) // commented out to duplicate ?bug? in GENXPT
591 x = 0.999;
592 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
593 if( (backwardKinetic+kineticEnergy) < xxx*backwardEnergy )
594 {
595 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
596 backwardKinetic += kineticEnergy;
597 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
598 pseudoParticle[6].SetMomentum( 0.0 ); // set the z-momentum
599 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
600 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
601 phi += pi + normal() * pi / 12.0;
602 if( phi > twopi )phi -= twopi;
603 if( phi < 0.0 )phi = twopi - phi;
604 outerCounter = 2; // leave outer loop
605 eliminateThisParticle = false; // don't eliminate this particle
606 resetEnergies = false;
607 break; // leave inner loop
608 }
609 if( innerCounter > 5 )break; // leave inner loop
610 if( forwardEnergy >= vecMass ) // switch sides
611 {
612 vec[i]->SetSide( 1 );
613 forwardEnergy -= vecMass;
614 backwardEnergy += vecMass;
615 backwardCount--;
616 }
617 }
618 G4ThreeVector momentum = vec[i]->GetMomentum();
619 vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
620 pt *= 0.9;
621 dndl[19] *= 0.9;
622 } // closes inner loop
623 if( resetEnergies )
624 {
625 // if we get to here, the inner loop has been Done 6 Times
626 // reset the kinetic energies of previously Done particles, if they are lighter
627 // than protons and in the forward hemisphere
628 // then continue with outer loop
629 //
630 forwardKinetic = 0.0;
631 backwardKinetic = 0.0;
632 pseudoParticle[4].SetZero();
633 pseudoParticle[5].SetZero();
634 for( l=i+1; l<vecLen; ++l )
635 {
636 if (vec[l]->GetSide() > 0 ||
637 vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
638 vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
639
640 G4double tempMass = vec[l]->GetMass()/MeV;
641 totalEnergy = 0.95*vec[l]->GetTotalEnergy()/MeV + 0.05*tempMass;
642 totalEnergy = std::max( tempMass, totalEnergy );
643 vec[l]->SetTotalEnergy( totalEnergy*MeV );
644 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) );
645 pp1 = vec[l]->GetMomentum().mag()/MeV;
646 if( pp1 < 1.0e-6*GeV )
647 {
648 G4double costheta = 2.*G4UniformRand() - 1.;
649 G4double sintheta = std::sqrt(1. - costheta*costheta);
650 G4double phi2 = twopi*G4UniformRand();
651 vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
652 pp*sintheta*std::sin(phi2)*MeV,
653 pp*costheta*MeV ) ;
654 } else {
655 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
656 }
657 G4double px = vec[l]->GetMomentum().x()/MeV;
658 G4double py = vec[l]->GetMomentum().y()/MeV;
659 pt = std::max( 1.0, std::sqrt( px*px + py*py ) )/GeV;
660 if( vec[l]->GetSide() > 0 )
661 {
662 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
663 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
664 } else {
665 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
666 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
667 }
668 } // if pi, K or forward
669 } // for l
670 } // if resetEnergies
671 } // closes outer loop
672
673 if( eliminateThisParticle && vec[i]->GetMayBeKilled()) // not enough energy, eliminate this particle
674 {
675 if( vec[i]->GetSide() > 0 )
676 {
677 --forwardCount;
678 forwardEnergy += vecMass;
679 } else {
680 if( vec[i]->GetSide() == -2 )
681 {
682 --extraNucleonCount;
683 extraNucleonMass -= vecMass;
684 backwardEnergy -= vecMass;
685 }
686 --backwardCount;
687 backwardEnergy += vecMass;
688 }
689 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
690 G4ReactionProduct *temp = vec[vecLen-1];
691 delete temp;
692 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
693 if( --vecLen == 0 )return false; // all the secondaries have been eliminated
694 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
695 pseudoParticle[6].SetMomentum( 0.0 ); // set z-momentum
696 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
697 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
698 phi += pi + normal() * pi / 12.0;
699 if( phi > twopi )phi -= twopi;
700 if( phi < 0.0 )phi = twopi - phi;
701 }
702 } // closes main for loop
703
704 //
705 // for the incident particle: it was placed in the forward hemisphere
706 // set pt and phi values, they are changed somewhat in the iteration loop
707 // set mass parameter for lambda fragmentation model
708 //
709 G4double phi = G4UniformRand()*twopi;
710 G4double ran = -std::log(1.0-G4UniformRand());
711 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
712 aspar = 0.60;
713 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
714 } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
715 aspar = 0.50;
716 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
717 } else {
718 aspar = 0.40;
719 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
720 }
721
722 for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
723 currentParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
724 et = pseudoParticle[0].GetTotalEnergy()/GeV;
725 dndl[0] = 0.0;
726 vecMass = currentParticle.GetMass()/GeV;
727 for( l=1; l<20; ++l )
728 {
729 x = (binl[l]+binl[l-1])/2.;
730 if( x > 1.0/pt )
731 dndl[l] += dndl[l-1]; // changed from just = on 02 April 98
732 else
733 dndl[l] = aspar/std::sqrt( std::pow((1.+sqr(aspar*x)),3) ) *
734 (binl[l]-binl[l-1]) * et / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass ) +
735 dndl[l-1];
736 }
737 ran = G4UniformRand()*dndl[19];
738 l = 1;
739 while( (ran>dndl[l]) && (l<20) )l++;
740 l = std::min( 19, l );
741 x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
742 currentParticle.SetMomentum( x*et*GeV ); // set the z-momentum
743 if( forwardEnergy < forwardKinetic )
744 totalEnergy = vecMass + 0.04*std::fabs(normal());
745 else
746 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
747 currentParticle.SetTotalEnergy( totalEnergy*GeV );
748 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
749 pp1 = currentParticle.GetMomentum().mag()/MeV;
750 if( pp1 < 1.0e-6*GeV )
751 {
752 G4double costheta = 2.*G4UniformRand() - 1.;
753 G4double sintheta = std::sqrt(1. - costheta*costheta);
754 G4double phi2 = twopi*G4UniformRand();
755 currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
756 pp*sintheta*std::sin(phi2)*MeV,
757 pp*costheta*MeV ) ;
758 } else {
759 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
760 }
761 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
762
763 //
764 // Current particle now finished
765 //
766 // Begin target particle
767 //
768
769 if( backwardNucleonCount < 18 )
770 {
771 targetParticle.SetSide( -3 );
772 ++backwardNucleonCount;
773 }
774 else
775 {
776 // set pt and phi values, they are changed somewhat in the iteration loop
777 // set mass parameter for lambda fragmentation model
778 //
779 vecMass = targetParticle.GetMass()/GeV;
780 ran = -std::log(1.0-G4UniformRand());
781 aspar = 0.40;
782 pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
783 targetParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
784 for( G4int j=0; j<20; ++j )binl[j] = (j-1.)/(19.*pt);
785 et = pseudoParticle[1].GetTotalEnergy()/GeV;
786 dndl[0] = 0.0;
787 outerCounter = 0;
788 eliminateThisParticle = true; // should never eliminate the target particle
789 resetEnergies = true;
790 while( ++outerCounter < 3 ) // start of outer iteration loop
791 {
792 for( l=1; l<20; ++l )
793 {
794 x = (binl[l]+binl[l-1])/2.;
795 if( x > 1.0/pt )
796 dndl[l] += dndl[l-1]; // changed from just = on 02 April 98
797 else
798 dndl[l] = aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) ) *
799 (binl[l]-binl[l-1])*et / std::sqrt( pt*x*et*pt*x*et+pt*pt+vecMass*vecMass ) +
800 dndl[l-1];
801 }
802 innerCounter = 0;
803 while( ++innerCounter < 7 ) // start of inner iteration loop
804 {
805 l = 1;
806 ran = G4UniformRand()*dndl[19];
807 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
808 l = std::min( 19, l );
809 x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
810 if( targetParticle.GetSide() < 0 )x *= -1.;
811 targetParticle.SetMomentum( x*et*GeV ); // set the z-momentum
812 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
813 targetParticle.SetTotalEnergy( totalEnergy*GeV );
814 if( targetParticle.GetSide() < 0 )
815 {
816 if( extraNucleonCount > 19 )x=0.999;
817 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
818 if( (backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy )
819 {
820 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
821 backwardKinetic += totalEnergy - vecMass;
822 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
823 pseudoParticle[6].SetMomentum( 0.0 ); // set z-momentum
824 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
825 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
826 phi += pi + normal() * pi / 12.0;
827 if( phi > twopi )phi -= twopi;
828 if( phi < 0.0 )phi = twopi - phi;
829 outerCounter = 2; // leave outer loop
830 eliminateThisParticle = false; // don't eliminate this particle
831 resetEnergies = false;
832 break; // leave inner loop
833 }
834 if( innerCounter > 5 )break; // leave inner loop
835 if( forwardEnergy >= vecMass ) // switch sides
836 {
837 targetParticle.SetSide( 1 );
838 forwardEnergy -= vecMass;
839 backwardEnergy += vecMass;
840 --backwardCount;
841 }
842 G4ThreeVector momentum = targetParticle.GetMomentum();
843 targetParticle.SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
844 pt *= 0.9;
845 dndl[19] *= 0.9;
846 }
847 else // target has gone to forward side
848 {
849 if( forwardEnergy < forwardKinetic )
850 totalEnergy = vecMass + 0.04*std::fabs(normal());
851 else
852 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
853 targetParticle.SetTotalEnergy( totalEnergy*GeV );
854 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
855 pp1 = targetParticle.GetMomentum().mag()/MeV;
856 if( pp1 < 1.0e-6*GeV )
857 {
858 G4double costheta = 2.*G4UniformRand() - 1.;
859 G4double sintheta = std::sqrt(1. - costheta*costheta);
860 G4double phi2 = twopi*G4UniformRand();
861 targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
862 pp*sintheta*std::sin(phi2)*MeV,
863 pp*costheta*MeV ) ;
864 }
865 else
866 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
867
868 pseudoParticle[4] = pseudoParticle[4] + targetParticle;
869 outerCounter = 2; // leave outer loop
870 eliminateThisParticle = false; // don't eliminate this particle
871 resetEnergies = false;
872 break; // leave inner loop
873 }
874 } // closes inner loop
875 if( resetEnergies )
876 {
877 // if we get to here, the inner loop has been Done 6 Times
878 // reset the kinetic energies of previously Done particles, if they are lighter
879 // than protons and in the forward hemisphere
880 // then continue with outer loop
881
882 forwardKinetic = backwardKinetic = 0.0;
883 pseudoParticle[4].SetZero();
884 pseudoParticle[5].SetZero();
885 for( l=0; l<vecLen; ++l ) // changed from l=1 on 02 April 98
886 {
887 if (vec[l]->GetSide() > 0 ||
888 vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
889 vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
890 G4double tempMass = vec[l]->GetMass()/GeV;
891 totalEnergy =
892 std::max( tempMass, 0.95*vec[l]->GetTotalEnergy()/GeV + 0.05*tempMass );
893 vec[l]->SetTotalEnergy( totalEnergy*GeV );
894 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) )*GeV;
895 pp1 = vec[l]->GetMomentum().mag()/MeV;
896 if( pp1 < 1.0e-6*GeV )
897 {
898 G4double costheta = 2.*G4UniformRand() - 1.;
899 G4double sintheta = std::sqrt(1. - costheta*costheta);
900 G4double phi2 = twopi*G4UniformRand();
901 vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
902 pp*sintheta*std::sin(phi2)*MeV,
903 pp*costheta*MeV ) ;
904 }
905 else
906 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
907
908 pt = std::max( 0.001*GeV, std::sqrt( sqr(vec[l]->GetMomentum().x()/MeV) +
909 sqr(vec[l]->GetMomentum().y()/MeV) ) )/GeV;
910 if( vec[l]->GetSide() > 0)
911 {
912 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
913 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
914 } else {
915 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
916 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
917 }
918 } // if pi, K or forward
919 } // for l
920 } // if (resetEnergies)
921 } // closes outer loop
922 }
923
924 //
925 // Target particle finished.
926 //
927 // Now produce backward nucleons with a cluster model
928 //
929 pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
930 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
931 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
932 if( backwardNucleonCount == 1 ) // target particle is the only backward nucleon
933 {
934 G4double ekin =
935 std::min( backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV );
936
937 if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
938 vecMass = targetParticle.GetMass()/GeV;
939 totalEnergy = ekin+vecMass;
940 targetParticle.SetTotalEnergy( totalEnergy*GeV );
941 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
942 pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
943 if( pp1 < 1.0e-6*GeV )
944 {
945 G4double costheta = 2.*G4UniformRand() - 1.;
946 G4double sintheta = std::sqrt(1. - costheta*costheta);
947 G4double phi2 = twopi*G4UniformRand();
948 targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
949 pp*sintheta*std::sin(phi2)*MeV,
950 pp*costheta*MeV ) ;
951 } else {
952 targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1) );
953 }
954 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
955 }
956 else // more than one backward nucleon
957 {
958 const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
959 const G4double gpar[] = { 2.6, 2.6, 1.80, 1.30, 1.20 };
960 // Replaced the following min function to get correct behaviour on DEC.
961 // G4int tempCount = std::min( 5, backwardNucleonCount ) - 1;
962 G4int tempCount;
963 if (backwardNucleonCount < 5)
964 {
965 tempCount = backwardNucleonCount;
966 }
967 else
968 {
969 tempCount = 5;
970 }
971 tempCount--;
972 //cout << "backwardNucleonCount " << backwardNucleonCount << G4endl;
973 //cout << "tempCount " << tempCount << G4endl;
974 G4double rmb0 = 0.0;
975 if( targetParticle.GetSide() == -3 )
976 rmb0 += targetParticle.GetMass()/GeV;
977 for( i=0; i<vecLen; ++i )
978 {
979 if( vec[i]->GetSide() == -3 )rmb0 += vec[i]->GetMass()/GeV;
980 }
981 rmb = rmb0 + std::pow(-std::log(1.0-G4UniformRand()),cpar[tempCount]) / gpar[tempCount];
982 totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
983 vecMass = std::min( rmb, totalEnergy );
984 pseudoParticle[6].SetMass( vecMass*GeV );
985 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
986 pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
987 if( pp1 < 1.0e-6*GeV )
988 {
989 G4double costheta = 2.*G4UniformRand() - 1.;
990 G4double sintheta = std::sqrt(1. - costheta*costheta);
991 G4double phi2 = twopi*G4UniformRand();
992 pseudoParticle[6].SetMomentum( -pp*sintheta*std::cos(phi2)*MeV,
993 -pp*sintheta*std::sin(phi2)*MeV,
994 -pp*costheta*MeV ) ;
995 }
996 else
997 pseudoParticle[6].SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
998
999 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV; // tempV contains the backward nucleons
1000 tempV.Initialize( backwardNucleonCount );
1001 G4int tempLen = 0;
1002 if( targetParticle.GetSide() == -3 )tempV.SetElement( tempLen++, &targetParticle );
1003 for( i=0; i<vecLen; ++i )
1004 {
1005 if( vec[i]->GetSide() == -3 )tempV.SetElement( tempLen++, vec[i] );
1006 }
1007 if( tempLen != backwardNucleonCount )
1008 {
1009 G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
1010 G4cerr << "tempLen = " << tempLen;
1011 G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
1012 G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
1013 G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
1014 for( i=0; i<vecLen; ++i )
1015 G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
1016 G4Exception("G4ReactionDynamics::GenerateXandPt", "601",
1017 FatalException, "Mismatch in nucleon count");
1018 }
1019 constantCrossSection = true;
1020 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1021 if( tempLen >= 2 )
1022 {
1023 wgt = GenerateNBodyEvent(
1024 pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen );
1025 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1026 if( targetParticle.GetSide() == -3 )
1027 {
1028 targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
1029 // tempV contains the real stuff
1030 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
1031 }
1032 for( i=0; i<vecLen; ++i )
1033 {
1034 if( vec[i]->GetSide() == -3 )
1035 {
1036 vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
1037 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
1038 }
1039 }
1040 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1041 }
1042 }
1043 //
1044 // Lorentz transformation in lab system
1045 //
1046 if( vecLen == 0 )return false; // all the secondaries have been eliminated
1047 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1048
1049 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
1050 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
1051 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[1] );
1052
1053 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1054 //
1055 // leadFlag will be true
1056 // iff original particle is at least as heavy as K+ and not a proton or
1057 // neutron AND if incident particle is at least as heavy as K+ and it is
1058 // not a proton or neutron leadFlag is set to the incident particle
1059 // or
1060 // target particle is at least as heavy as K+ and it is not a proton or
1061 // neutron leadFlag is set to the target particle
1062 //
1063 G4bool leadingStrangeParticleHasChanged = true;
1064 if( leadFlag )
1065 {
1066 if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
1067 leadingStrangeParticleHasChanged = false;
1068 if( leadingStrangeParticleHasChanged &&
1069 ( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() ) )
1070 leadingStrangeParticleHasChanged = false;
1071 if( leadingStrangeParticleHasChanged )
1072 {
1073 for( i=0; i<vecLen; i++ )
1074 {
1075 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
1076 {
1077 leadingStrangeParticleHasChanged = false;
1078 break;
1079 }
1080 }
1081 }
1082 if( leadingStrangeParticleHasChanged )
1083 {
1084 G4bool leadTest =
1085 (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
1086 leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
1087 G4bool targetTest =
1088 (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
1089 targetParticle.GetDefinition()->GetParticleSubType() == "pi");
1090
1091 // following modified by JLC 22-Oct-97
1092
1093 if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
1094 {
1095 targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
1096 targetHasChanged = true;
1097 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1098 }
1099 else
1100 {
1101 currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
1102 incidentHasChanged = false;
1103 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1104 }
1105 }
1106 } // end of if( leadFlag )
1107
1108 // Get number of final state nucleons and nucleons remaining in
1109 // target nucleus
1110
1111 std::pair<G4int, G4int> finalStateNucleons =
1112 GetFinalStateNucleons(originalTarget, vec, vecLen);
1113
1114 G4int protonsInFinalState = finalStateNucleons.first;
1115 G4int neutronsInFinalState = finalStateNucleons.second;
1116
1117 G4int numberofFinalStateNucleons =
1118 protonsInFinalState + neutronsInFinalState;
1119
1120 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
1121 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
1122 originalIncident->GetDefinition()->GetPDGMass() <
1123 G4Lambda::Lambda()->GetPDGMass())
1124 numberofFinalStateNucleons++;
1125
1126 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
1127
1128 G4int PinNucleus = std::max(0,
1129 targetNucleus.GetZ_asInt() - protonsInFinalState);
1130 G4int NinNucleus = std::max(0,
1131 targetNucleus.GetN_asInt() - neutronsInFinalState);
1132
1133 pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
1134 pseudoParticle[3].SetMass( mOriginal*GeV );
1135 pseudoParticle[3].SetTotalEnergy(
1136 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
1137
1138 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
1139 G4int diff = 0;
1140 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
1141 if(numberofFinalStateNucleons == 1) diff = 0;
1142 pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
1143 pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV );
1144 pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV );
1145
1146 G4double theoreticalKinetic =
1147 pseudoParticle[3].GetTotalEnergy()/MeV +
1148 pseudoParticle[4].GetTotalEnergy()/MeV -
1149 currentParticle.GetMass()/MeV -
1150 targetParticle.GetMass()/MeV;
1151
1152 G4double simulatedKinetic =
1153 currentParticle.GetKineticEnergy()/MeV + targetParticle.GetKineticEnergy()/MeV;
1154
1155 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
1156 pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
1157 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
1158
1159 pseudoParticle[7].SetZero();
1160 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
1161 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
1162
1163 for( i=0; i<vecLen; ++i )
1164 {
1165 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
1166 simulatedKinetic += vec[i]->GetKineticEnergy()/MeV;
1167 theoreticalKinetic -= vec[i]->GetMass()/MeV;
1168 }
1169
1170 if( vecLen <= 16 && vecLen > 0 )
1171 {
1172 // must create a new set of ReactionProducts here because GenerateNBody will
1173 // modify the momenta for the particles, and we don't want to do this
1174 //
1175 G4ReactionProduct tempR[130];
1176 tempR[0] = currentParticle;
1177 tempR[1] = targetParticle;
1178 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
1179 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
1180 tempV.Initialize( vecLen+2 );
1181 G4int tempLen = 0;
1182 for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
1183 constantCrossSection = true;
1184
1185 wgt = GenerateNBodyEvent( pseudoParticle[3].GetTotalEnergy()/MeV+
1186 pseudoParticle[4].GetTotalEnergy()/MeV,
1187 constantCrossSection, tempV, tempLen );
1188 if (wgt == -1) {
1189 G4double Qvalue = 0;
1190 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
1191 wgt = GenerateNBodyEvent( Qvalue/MeV,
1192 constantCrossSection, tempV, tempLen );
1193 }
1194 if(wgt>-.5)
1195 {
1196 theoreticalKinetic = 0.0;
1197 for( i=0; i<tempLen; ++i )
1198 {
1199 pseudoParticle[6].Lorentz( *tempV[i], pseudoParticle[4] );
1200 theoreticalKinetic += pseudoParticle[6].GetKineticEnergy()/MeV;
1201 }
1202 }
1203 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1204 }
1205 //
1206 // Make sure, that the kinetic energies are correct
1207 //
1208 if( simulatedKinetic != 0.0 )
1209 {
1210 wgt = (theoreticalKinetic)/simulatedKinetic;
1211 theoreticalKinetic = currentParticle.GetKineticEnergy()/MeV * wgt;
1212 simulatedKinetic = theoreticalKinetic;
1213 currentParticle.SetKineticEnergy( theoreticalKinetic*MeV );
1214 pp = currentParticle.GetTotalMomentum()/MeV;
1215 pp1 = currentParticle.GetMomentum().mag()/MeV;
1216 if( pp1 < 1.0e-6*GeV )
1217 {
1218 G4double costheta = 2.*G4UniformRand() - 1.;
1219 G4double sintheta = std::sqrt(1. - costheta*costheta);
1220 G4double phi2 = twopi*G4UniformRand();
1221 currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
1222 pp*sintheta*std::sin(phi2)*MeV,
1223 pp*costheta*MeV ) ;
1224 }
1225 else
1226 {
1227 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
1228 }
1229 theoreticalKinetic = targetParticle.GetKineticEnergy()/MeV * wgt;
1230 targetParticle.SetKineticEnergy( theoreticalKinetic*MeV );
1231 simulatedKinetic += theoreticalKinetic;
1232 pp = targetParticle.GetTotalMomentum()/MeV;
1233 pp1 = targetParticle.GetMomentum().mag()/MeV;
1234 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1235 if( pp1 < 1.0e-6*GeV )
1236 {
1237 G4double costheta = 2.*G4UniformRand() - 1.;
1238 G4double sintheta = std::sqrt(1. - costheta*costheta);
1239 G4double phi2 = twopi*G4UniformRand();
1240 targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
1241 pp*sintheta*std::sin(phi2)*MeV,
1242 pp*costheta*MeV ) ;
1243 } else {
1244 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
1245 }
1246 for( i=0; i<vecLen; ++i )
1247 {
1248 theoreticalKinetic = vec[i]->GetKineticEnergy()/MeV * wgt;
1249 simulatedKinetic += theoreticalKinetic;
1250 vec[i]->SetKineticEnergy( theoreticalKinetic*MeV );
1251 pp = vec[i]->GetTotalMomentum()/MeV;
1252 pp1 = vec[i]->GetMomentum().mag()/MeV;
1253 if( pp1 < 1.0e-6*GeV )
1254 {
1255 G4double costheta = 2.*G4UniformRand() - 1.;
1256 G4double sintheta = std::sqrt(1. - costheta*costheta);
1257 G4double phi2 = twopi*G4UniformRand();
1258 vec[i]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
1259 pp*sintheta*std::sin(phi2)*MeV,
1260 pp*costheta*MeV ) ;
1261 }
1262 else
1263 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
1264 }
1265 }
1266 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1267
1268 Rotate( numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
1269 modifiedOriginal, originalIncident, targetNucleus,
1270 currentParticle, targetParticle, vec, vecLen );
1271 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1272 //
1273 // add black track particles
1274 // the total number of particles produced is restricted to 198
1275 // this may have influence on very high energies
1276 //
1277 if( atomicWeight >= 1.5 )
1278 {
1279 // npnb is number of proton/neutron black track particles
1280 // ndta is the number of deuterons, tritons, and alphas produced
1281 // epnb is the kinetic energy available for proton/neutron black track particles
1282 // edta is the kinetic energy available for deuteron/triton/alpha particles
1283 //
1284 G4int npnb = 0;
1285 G4int ndta = 0;
1286
1287 G4double epnb, edta;
1288 if (veryForward) {
1289 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
1290 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
1291 } else {
1292 epnb = targetNucleus.GetPNBlackTrackEnergy();
1293 edta = targetNucleus.GetDTABlackTrackEnergy();
1294 }
1295
1296 const G4double pnCutOff = 0.001;
1297 const G4double dtaCutOff = 0.001;
1298 const G4double kineticMinimum = 1.e-6;
1299 const G4double kineticFactor = -0.010;
1300 G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
1301 const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
1302 if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
1303 if( epnb >= pnCutOff )
1304 {
1305 npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
1306 if( numberofFinalStateNucleons + npnb > atomicWeight )
1307 npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
1308 npnb = std::min( npnb, 127-vecLen );
1309 }
1310 if( edta >= dtaCutOff )
1311 {
1312 ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
1313 ndta = std::min( ndta, 127-vecLen );
1314 }
1315 if (npnb == 0 && ndta == 0) npnb = 1;
1316
1317 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1318
1319 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
1320 kineticFactor, modifiedOriginal,
1321 PinNucleus, NinNucleus, targetNucleus,
1322 vec, vecLen);
1323
1324 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1325 }
1326 //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
1327 // MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
1328 //
1329 // calculate time delay for nuclear reactions
1330 //
1331 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
1332 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
1333 else
1334 currentParticle.SetTOF( 1.0 );
1335 return true;
1336 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1337 }
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::mag
double mag() const
G4FastVector::SetElement
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
G4FastVector::Initialize
void Initialize(G4int items)
Definition: G4FastVector.hh:63
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4Nucleus::GetAnnihilationPNBlackTrackEnergy
G4double GetAnnihilationPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:156
G4Nucleus::GetPNBlackTrackEnergy
G4double GetPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:150
G4Nucleus::GetAnnihilationDTABlackTrackEnergy
G4double GetAnnihilationDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:159
G4Nucleus::GetN_asInt
G4int GetN_asInt() const
Definition: G4Nucleus.hh:112
G4Nucleus::GetDTABlackTrackEnergy
G4double GetDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:153
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:137
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4ReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen)
Definition: G4ReactionDynamics.cc:2491
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:166
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:123
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:135
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::SetNewlyAdded
void SetNewlyAdded(const G4bool f)
Definition: G4ReactionProduct.hh:162
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:120
G4ReactionProduct::SetSide
void SetSide(const G4int sid)
Definition: G4ReactionProduct.hh:156
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:200
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:129
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:104
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:155
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:147
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:144
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70
sqr
T sqr(const T &x)
Definition: templates.hh:145

Referenced by G4InelasticInteraction::CalculateMomenta().

◆ NuclearReaction()

void G4ReactionDynamics::NuclearReaction ( G4FastVector< G4ReactionProduct, 4 > &  vec,
G4int &  vecLen,
const G4HadProjectile *  originalIncident,
const G4Nucleus &  aNucleus,
const G4double  theAtomicMass,
const G4double *  massVec 
)

Definition at line 3777 of file G4ReactionDynamics.cc.

3784 {
3785 // derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
3786 //
3787 // Nuclear reaction kinematics at low energies
3788 //
3789 G4ParticleDefinition *aGamma = G4Gamma::Gamma();
3790 G4ParticleDefinition *aProton = G4Proton::Proton();
3791 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3792 G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
3793 G4ParticleDefinition *aTriton = G4Triton::Triton();
3794 G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
3795
3796 const G4double aProtonMass = aProton->GetPDGMass()/MeV;
3797 const G4double aNeutronMass = aNeutron->GetPDGMass()/MeV;
3798 const G4double aDeuteronMass = aDeuteron->GetPDGMass()/MeV;
3799 const G4double aTritonMass = aTriton->GetPDGMass()/MeV;
3800 const G4double anAlphaMass = anAlpha->GetPDGMass()/MeV;
3801
3802 G4ReactionProduct currentParticle;
3803 currentParticle = *originalIncident;
3804 //
3805 // Set beam particle, take kinetic energy of current particle as the
3806 // fundamental quantity. Due to the difficult kinematic, all masses have to
3807 // be assigned the best measured values
3808 //
3809 G4double p = currentParticle.GetTotalMomentum();
3810 G4double pp = currentParticle.GetMomentum().mag();
3811 if( pp <= 0.001*MeV )
3812 {
3813 G4double phinve = twopi*G4UniformRand();
3814 G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
3815 currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3816 p*std::sin(rthnve)*std::sin(phinve),
3817 p*std::cos(rthnve) );
3818 }
3819 else
3820 currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
3821 //
3822 // calculate Q-value of reactions
3823 //
3824 G4double currentKinetic = currentParticle.GetKineticEnergy()/MeV;
3825 G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/MeV;
3826 G4double qv = currentKinetic + theAtomicMass + currentMass;
3827
3828 G4double qval[9];
3829 qval[0] = qv - mass[0];
3830 qval[1] = qv - mass[1] - aNeutronMass;
3831 qval[2] = qv - mass[2] - aProtonMass;
3832 qval[3] = qv - mass[3] - aDeuteronMass;
3833 qval[4] = qv - mass[4] - aTritonMass;
3834 qval[5] = qv - mass[5] - anAlphaMass;
3835 qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
3836 qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
3837 qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
3838
3839 if( currentParticle.GetDefinition() == aNeutron )
3840 {
3841 const G4double A = G4double(targetNucleus.GetA_asInt()); // atomic weight
3842 if( G4UniformRand() > ((A-1.0)/230.0)*((A-1.0)/230.0) )
3843 qval[0] = 0.0;
3844 if( G4UniformRand() >= currentKinetic/7.9254*A )
3845 qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
3846 }
3847 else
3848 qval[0] = 0.0;
3849
3850 G4int i;
3851 qv = 0.0;
3852 for( i=0; i<9; ++i )
3853 {
3854 if( mass[i] < 500.0*MeV )qval[i] = 0.0;
3855 if( qval[i] < 0.0 )qval[i] = 0.0;
3856 qv += qval[i];
3857 }
3858 G4double qv1 = 0.0;
3859 G4double ran = G4UniformRand();
3860 G4int index;
3861 for( index=0; index<9; ++index )
3862 {
3863 if( qval[index] > 0.0 )
3864 {
3865 qv1 += qval[index]/qv;
3866 if( ran <= qv1 )break;
3867 }
3868 }
3869 if( index == 9 ) // loop continued to the end
3870 {
3871 throw G4HadronicException(__FILE__, __LINE__,
3872 "G4ReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
3873 }
3874 G4double ke = currentParticle.GetKineticEnergy()/GeV;
3875 G4int nt = 2;
3876 if( (index>=6) || (G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
3877
3878 G4ReactionProduct **v = new G4ReactionProduct * [3];
3879 v[0] = new G4ReactionProduct;
3880 v[1] = new G4ReactionProduct;
3881 v[2] = new G4ReactionProduct;
3882
3883 v[0]->SetMass( mass[index]*MeV );
3884 switch( index )
3885 {
3886 case 0:
3887 v[1]->SetDefinition( aGamma );
3888 v[2]->SetDefinition( aGamma );
3889 break;
3890 case 1:
3891 v[1]->SetDefinition( aNeutron );
3892 v[2]->SetDefinition( aGamma );
3893 break;
3894 case 2:
3895 v[1]->SetDefinition( aProton );
3896 v[2]->SetDefinition( aGamma );
3897 break;
3898 case 3:
3899 v[1]->SetDefinition( aDeuteron );
3900 v[2]->SetDefinition( aGamma );
3901 break;
3902 case 4:
3903 v[1]->SetDefinition( aTriton );
3904 v[2]->SetDefinition( aGamma );
3905 break;
3906 case 5:
3907 v[1]->SetDefinition( anAlpha );
3908 v[2]->SetDefinition( aGamma );
3909 break;
3910 case 6:
3911 v[1]->SetDefinition( aNeutron );
3912 v[2]->SetDefinition( aNeutron );
3913 break;
3914 case 7:
3915 v[1]->SetDefinition( aNeutron );
3916 v[2]->SetDefinition( aProton );
3917 break;
3918 case 8:
3919 v[1]->SetDefinition( aProton );
3920 v[2]->SetDefinition( aProton );
3921 break;
3922 }
3923 //
3924 // calculate centre of mass energy
3925 //
3926 G4ReactionProduct pseudo1;
3927 pseudo1.SetMass( theAtomicMass*MeV );
3928 pseudo1.SetTotalEnergy( theAtomicMass*MeV );
3929 G4ReactionProduct pseudo2 = currentParticle + pseudo1;
3930 pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
3931 //
3932 // use phase space routine in centre of mass system
3933 //
3934 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
3935 tempV.Initialize( nt );
3936 G4int tempLen = 0;
3937 tempV.SetElement( tempLen++, v[0] );
3938 tempV.SetElement( tempLen++, v[1] );
3939 if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
3940 G4bool constantCrossSection = true;
3941 GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
3942 v[0]->Lorentz( *v[0], pseudo2 );
3943 v[1]->Lorentz( *v[1], pseudo2 );
3944 if( nt == 3 )v[2]->Lorentz( *v[2], pseudo2 );
3945
3946 G4bool particleIsDefined = false;
3947 if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
3948 {
3949 v[0]->SetDefinition( aProton );
3950 particleIsDefined = true;
3951 }
3952 else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
3953 {
3954 v[0]->SetDefinition( aNeutron );
3955 particleIsDefined = true;
3956 }
3957 else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
3958 {
3959 v[0]->SetDefinition( aDeuteron );
3960 particleIsDefined = true;
3961 }
3962 else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
3963 {
3964 v[0]->SetDefinition( aTriton );
3965 particleIsDefined = true;
3966 }
3967 else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
3968 {
3969 v[0]->SetDefinition( anAlpha );
3970 particleIsDefined = true;
3971 }
3972 currentParticle.SetKineticEnergy(
3973 std::max( 0.001, currentParticle.GetKineticEnergy()/MeV ) );
3974 p = currentParticle.GetTotalMomentum();
3975 pp = currentParticle.GetMomentum().mag();
3976 if( pp <= 0.001*MeV )
3977 {
3978 G4double phinve = twopi*G4UniformRand();
3979 G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
3980 currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3981 p*std::sin(rthnve)*std::sin(phinve),
3982 p*std::cos(rthnve) );
3983 }
3984 else
3985 currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
3986
3987 if( particleIsDefined )
3988 {
3989 v[0]->SetKineticEnergy(
3990 std::max( 0.001, 0.5*G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
3991 p = v[0]->GetTotalMomentum();
3992 pp = v[0]->GetMomentum().mag();
3993 if( pp <= 0.001*MeV )
3994 {
3995 G4double phinve = twopi*G4UniformRand();
3996 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
3997 v[0]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3998 p*std::sin(rthnve)*std::sin(phinve),
3999 p*std::cos(rthnve) );
4000 }
4001 else
4002 v[0]->SetMomentum( v[0]->GetMomentum() * (p/pp) );
4003 }
4004 if( (v[1]->GetDefinition() == aDeuteron) ||
4005 (v[1]->GetDefinition() == aTriton) ||
4006 (v[1]->GetDefinition() == anAlpha) )
4007 v[1]->SetKineticEnergy(
4008 std::max( 0.001, 0.5*G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
4009 else
4010 v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/MeV ) );
4011
4012 p = v[1]->GetTotalMomentum();
4013 pp = v[1]->GetMomentum().mag();
4014 if( pp <= 0.001*MeV )
4015 {
4016 G4double phinve = twopi*G4UniformRand();
4017 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
4018 v[1]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
4019 p*std::sin(rthnve)*std::sin(phinve),
4020 p*std::cos(rthnve) );
4021 }
4022 else
4023 v[1]->SetMomentum( v[1]->GetMomentum() * (p/pp) );
4024
4025 if( nt == 3 )
4026 {
4027 if( (v[2]->GetDefinition() == aDeuteron) ||
4028 (v[2]->GetDefinition() == aTriton) ||
4029 (v[2]->GetDefinition() == anAlpha) )
4030 v[2]->SetKineticEnergy(
4031 std::max( 0.001, 0.5*G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
4032 else
4033 v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/MeV ) );
4034
4035 p = v[2]->GetTotalMomentum();
4036 pp = v[2]->GetMomentum().mag();
4037 if( pp <= 0.001*MeV )
4038 {
4039 G4double phinve = twopi*G4UniformRand();
4040 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
4041 v[2]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
4042 p*std::sin(rthnve)*std::sin(phinve),
4043 p*std::cos(rthnve) );
4044 }
4045 else
4046 v[2]->SetMomentum( v[2]->GetMomentum() * (p/pp) );
4047 }
4048 G4int del;
4049 for(del=0; del<vecLen; del++) delete vec[del];
4050 vecLen = 0;
4051 if( particleIsDefined )
4052 {
4053 vec.SetElement( vecLen++, v[0] );
4054 }
4055 else
4056 {
4057 delete v[0];
4058 }
4059 vec.SetElement( vecLen++, v[1] );
4060 if( nt == 3 )
4061 {
4062 vec.SetElement( vecLen++, v[2] );
4063 }
4064 else
4065 {
4066 delete v[2];
4067 }
4068 delete [] v;
4069 return;
4070 }
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95

Referenced by G4LEAlphaInelastic::ApplyYourself(), G4LEDeuteronInelastic::ApplyYourself(), and G4LETritonInelastic::ApplyYourself().

◆ ProduceStrangeParticlePairs()

void G4ReactionDynamics::ProduceStrangeParticlePairs ( G4FastVector< G4ReactionProduct, GHADLISTSIZE > &  vec,
G4int &  vecLen,
const G4ReactionProduct &  modifiedOriginal,
const G4DynamicParticle *  originalTarget,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged 
)

Definition at line 3386 of file G4ReactionDynamics.cc.

3395 {
3396 // derived from original FORTRAN code STPAIR by H. Fesefeldt (16-Dec-1987)
3397 //
3398 // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
3399 // K+ Y0, K0 Y+, K0 Y-
3400 // For antibaryon induced reactions half of the cross sections KB YB
3401 // pairs are produced. Charge is not conserved, no experimental data available
3402 // for exclusive reactions, therefore some average behaviour assumed.
3403 // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
3404 //
3405 if( vecLen == 0 )return;
3406 //
3407 // the following protects against annihilation processes
3408 //
3409 if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return;
3410
3411 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
3412 const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
3413 G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
3414 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
3415 targetMass*targetMass +
3416 2.0*targetMass*etOriginal ); // GeV
3417 G4double currentMass = currentParticle.GetMass()/GeV;
3418 G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
3419 if( availableEnergy <= 1.0 )return;
3420
3421 G4ParticleDefinition *aProton = G4Proton::Proton();
3422 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
3423 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3424 G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
3425 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
3426 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
3427 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
3428 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
3429 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
3430 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
3431 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
3432 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
3433 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
3434 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
3435 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
3436 G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
3437
3438 const G4double protonMass = aProton->GetPDGMass()/GeV;
3439 const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
3440 //
3441 // determine the center of mass energy bin
3442 //
3443 const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
3444
3445 G4int ibin, i3, i4;
3446 G4double avk, avy, avn, ran;
3447 G4int i = 1;
3448 while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
3449 if( i == 12 )
3450 ibin = 11;
3451 else
3452 ibin = i;
3453 //
3454 // the fortran code chooses a random replacement of produced kaons
3455 // but does not take into account charge conservation
3456 //
3457 if( vecLen == 1 ) // we know that vecLen > 0
3458 {
3459 i3 = 0;
3460 i4 = 1; // note that we will be adding a new secondary particle in this case only
3461 }
3462 else // otherwise 0 <= i3,i4 < vecLen
3463 {
3464 ran = G4UniformRand();
3465 while( ran == 1.0 )ran = G4UniformRand();
3466 i4 = i3 = G4int( vecLen*ran );
3467 while( i3 == i4 )
3468 {
3469 ran = G4UniformRand();
3470 while( ran == 1.0 )ran = G4UniformRand();
3471 i4 = G4int( vecLen*ran );
3472 }
3473 }
3474 //
3475 // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
3476 //
3477 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
3478 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
3479 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
3480 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
3481 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
3482 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
3483
3484 avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3485 /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
3486 avk = std::exp(avk);
3487
3488 avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3489 /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
3490 avy = std::exp(avy);
3491
3492 avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3493 /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
3494 avn = std::exp(avn);
3495
3496 if( avk+avy+avn <= 0.0 )return;
3497
3498 if( currentMass < protonMass )avy /= 2.0;
3499 if( targetMass < protonMass )avy = 0.0;
3500 avy += avk+avn;
3501 avk += avn;
3502 ran = G4UniformRand();
3503 if( ran < avn )
3504 {
3505 if( availableEnergy < 2.0 )return;
3506 if( vecLen == 1 ) // add a new secondary
3507 {
3508 G4ReactionProduct *p1 = new G4ReactionProduct;
3509 if( G4UniformRand() < 0.5 )
3510 {
3511 vec[0]->SetDefinition( aNeutron );
3512 p1->SetDefinition( anAntiNeutron );
3513 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
3514 vec[0]->SetMayBeKilled(false);
3515 p1->SetMayBeKilled(false);
3516 }
3517 else
3518 {
3519 vec[0]->SetDefinition( aProton );
3520 p1->SetDefinition( anAntiProton );
3521 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
3522 vec[0]->SetMayBeKilled(false);
3523 p1->SetMayBeKilled(false);
3524 }
3525 vec.SetElement( vecLen++, p1 );
3526 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3527 }
3528 else
3529 { // replace two secondaries
3530 if( G4UniformRand() < 0.5 )
3531 {
3532 vec[i3]->SetDefinition( aNeutron );
3533 vec[i4]->SetDefinition( anAntiNeutron );
3534 vec[i3]->SetMayBeKilled(false);
3535 vec[i4]->SetMayBeKilled(false);
3536 }
3537 else
3538 {
3539 vec[i3]->SetDefinition( aProton );
3540 vec[i4]->SetDefinition( anAntiProton );
3541 vec[i3]->SetMayBeKilled(false);
3542 vec[i4]->SetMayBeKilled(false);
3543 }
3544 }
3545 }
3546 else if( ran < avk )
3547 {
3548 if( availableEnergy < 1.0 )return;
3549
3550 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
3551 0.6875, 0.7500, 0.8750, 1.000 };
3552 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
3553 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
3554 ran = G4UniformRand();
3555 i = 0;
3556 while( (i<9) && (ran>=kkb[i]) )++i;
3557 if( i == 9 )return;
3558 //
3559 // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
3560 // charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 -
3561 //
3562 switch( ipakkb1[i] )
3563 {
3564 case 10:
3565 vec[i3]->SetDefinition( aKaonPlus );
3566 vec[i3]->SetMayBeKilled(false);
3567 break;
3568 case 11:
3569 vec[i3]->SetDefinition( aKaonZS );
3570 vec[i3]->SetMayBeKilled(false);
3571 break;
3572 case 12:
3573 vec[i3]->SetDefinition( aKaonZL );
3574 vec[i3]->SetMayBeKilled(false);
3575 break;
3576 }
3577 if( vecLen == 1 ) // add a secondary
3578 {
3579 G4ReactionProduct *p1 = new G4ReactionProduct;
3580 switch( ipakkb2[i] )
3581 {
3582 case 11:
3583 p1->SetDefinition( aKaonZS );
3584 p1->SetMayBeKilled(false);
3585 break;
3586 case 12:
3587 p1->SetDefinition( aKaonZL );
3588 p1->SetMayBeKilled(false);
3589 break;
3590 case 13:
3591 p1->SetDefinition( aKaonMinus );
3592 p1->SetMayBeKilled(false);
3593 break;
3594 }
3595 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
3596 vec.SetElement( vecLen++, p1 );
3597 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3598 }
3599 else // replace
3600 {
3601 switch( ipakkb2[i] )
3602 {
3603 case 11:
3604 vec[i4]->SetDefinition( aKaonZS );
3605 vec[i4]->SetMayBeKilled(false);
3606 break;
3607 case 12:
3608 vec[i4]->SetDefinition( aKaonZL );
3609 vec[i4]->SetMayBeKilled(false);
3610 break;
3611 case 13:
3612 vec[i4]->SetDefinition( aKaonMinus );
3613 vec[i4]->SetMayBeKilled(false);
3614 break;
3615 }
3616 }
3617 }
3618 else if( ran < avy )
3619 {
3620 if( availableEnergy < 1.6 )return;
3621
3622 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
3623 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
3624 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
3625 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
3626 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
3627 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
3628 ran = G4UniformRand();
3629 i = 0;
3630 while( (i<12) && (ran>ky[i]) )++i;
3631 if( i == 12 )return;
3632 if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
3633 {
3634 // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
3635 // 0 + 0 0 0 0 + + + 0 + 0
3636 //
3637 // 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
3638 // 0 + 0 0 0 0 - + - 0 - 0
3639 switch( ipaky1[i] )
3640 {
3641 case 18:
3642 targetParticle.SetDefinition( aLambda );
3643 break;
3644 case 20:
3645 targetParticle.SetDefinition( aSigmaPlus );
3646 break;
3647 case 21:
3648 targetParticle.SetDefinition( aSigmaZero );
3649 break;
3650 case 22:
3651 targetParticle.SetDefinition( aSigmaMinus );
3652 break;
3653 }
3654 targetHasChanged = true;
3655 switch( ipaky2[i] )
3656 {
3657 case 10:
3658 vec[i3]->SetDefinition( aKaonPlus );
3659 vec[i3]->SetMayBeKilled(false);
3660 break;
3661 case 11:
3662 vec[i3]->SetDefinition( aKaonZS );
3663 vec[i3]->SetMayBeKilled(false);
3664 break;
3665 case 12:
3666 vec[i3]->SetDefinition( aKaonZL );
3667 vec[i3]->SetMayBeKilled(false);
3668 break;
3669 }
3670 }
3671 else // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
3672 {
3673 // ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
3674 // 24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
3675 if( (currentParticle.GetDefinition() == anAntiProton) ||
3676 (currentParticle.GetDefinition() == anAntiNeutron) ||
3677 (currentParticle.GetDefinition() == anAntiLambda) ||
3678 (currentMass > sigmaMinusMass) )
3679 {
3680 switch( ipakyb1[i] )
3681 {
3682 case 19:
3683 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
3684 break;
3685 case 23:
3686 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
3687 break;
3688 case 24:
3689 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
3690 break;
3691 case 25:
3692 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
3693 break;
3694 }
3695 incidentHasChanged = true;
3696 switch( ipakyb2[i] )
3697 {
3698 case 11:
3699 vec[i3]->SetDefinition( aKaonZS );
3700 vec[i3]->SetMayBeKilled(false);
3701 break;
3702 case 12:
3703 vec[i3]->SetDefinition( aKaonZL );
3704 vec[i3]->SetMayBeKilled(false);
3705 break;
3706 case 13:
3707 vec[i3]->SetDefinition( aKaonMinus );
3708 vec[i3]->SetMayBeKilled(false);
3709 break;
3710 }
3711 }
3712 else
3713 {
3714 switch( ipaky1[i] )
3715 {
3716 case 18:
3717 currentParticle.SetDefinitionAndUpdateE( aLambda );
3718 break;
3719 case 20:
3720 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
3721 break;
3722 case 21:
3723 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
3724 break;
3725 case 22:
3726 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
3727 break;
3728 }
3729 incidentHasChanged = true;
3730 switch( ipaky2[i] )
3731 {
3732 case 10:
3733 vec[i3]->SetDefinition( aKaonPlus );
3734 vec[i3]->SetMayBeKilled(false);
3735 break;
3736 case 11:
3737 vec[i3]->SetDefinition( aKaonZS );
3738 vec[i3]->SetMayBeKilled(false);
3739 break;
3740 case 12:
3741 vec[i3]->SetDefinition( aKaonZL );
3742 vec[i3]->SetMayBeKilled(false);
3743 break;
3744 }
3745 }
3746 }
3747 }
3748 else return;
3749 //
3750 // check the available energy
3751 // if there is not enough energy for kkb/ky pair production
3752 // then reduce the number of secondary particles
3753 // NOTE:
3754 // the number of secondaries may have been changed
3755 // the incident and/or target particles may have changed
3756 // charge conservation is ignored (as well as strangness conservation)
3757 //
3758 currentMass = currentParticle.GetMass()/GeV;
3759 targetMass = targetParticle.GetMass()/GeV;
3760
3761 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
3762 for( i=0; i<vecLen; ++i )
3763 {
3764 energyCheck -= vec[i]->GetMass()/GeV;
3765 if( energyCheck < 0.0 ) // chop off the secondary List
3766 {
3767 vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
3768 G4int j;
3769 for(j=i; j<vecLen; j++) delete vec[j];
3770 break;
3771 }
3772 }
3773 return;
3774 }
G4AntiLambda::AntiLambda
static G4AntiLambda * AntiLambda()
Definition: G4AntiLambda.cc:108
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4AntiSigmaMinus::AntiSigmaMinus
static G4AntiSigmaMinus * AntiSigmaMinus()
Definition: G4AntiSigmaMinus.cc:106
G4AntiSigmaPlus::AntiSigmaPlus
static G4AntiSigmaPlus * AntiSigmaPlus()
Definition: G4AntiSigmaPlus.cc:108
G4AntiSigmaZero::AntiSigmaZero
static G4AntiSigmaZero * AntiSigmaZero()
Definition: G4AntiSigmaZero.cc:99
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4KaonMinus::KaonMinus
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4KaonZeroLong::KaonZeroLong
static G4KaonZeroLong * KaonZeroLong()
Definition: G4KaonZeroLong.cc:115
G4KaonZeroShort::KaonZeroShort
static G4KaonZeroShort * KaonZeroShort()
Definition: G4KaonZeroShort.cc:104
G4ReactionProduct::SetMayBeKilled
void SetMayBeKilled(const G4bool f)
Definition: G4ReactionProduct.hh:168
G4ReactionProduct::SetDefinitionAndUpdateE
void SetDefinitionAndUpdateE(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:142
G4SigmaMinus::SigmaMinus
static G4SigmaMinus * SigmaMinus()
Definition: G4SigmaMinus.cc:106
G4SigmaPlus::SigmaPlus
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
G4SigmaZero::SigmaZero
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99

Referenced by G4InelasticInteraction::CalculateMomenta().

◆ SuppressChargedPions()

void G4ReactionDynamics::SuppressChargedPions ( G4FastVector< G4ReactionProduct, GHADLISTSIZE > &  vec,
G4int &  vecLen,
const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged 
)

Definition at line 1339 of file G4ReactionDynamics.cc.

1348 {
1349 // this code was originally in the fortran code TWOCLU
1350 //
1351 // suppress charged pions, for various reasons
1352 //
1353 G4double mOriginal = modifiedOriginal.GetMass()/GeV;
1354 G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
1355 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
1356 G4double cmEnergy = std::sqrt( mOriginal*mOriginal + targetMass*targetMass +
1357 2.0*targetMass*etOriginal );
1358 G4double eAvailable = cmEnergy - mOriginal - targetMass;
1359 for (G4int i = 0; i < vecLen; i++) eAvailable -= vec[i]->GetMass()/GeV;
1360
1361 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
1362 const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
1363 const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/GeV;
1364
1365 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1366 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1367 G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
1368 G4ParticleDefinition *aProton = G4Proton::Proton();
1369 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1370 G4double piMass = aPiPlus->GetPDGMass()/GeV;
1371 G4double nucleonMass = aNeutron->GetPDGMass()/GeV;
1372
1373 const G4bool antiTest =
1374 modifiedOriginal.GetDefinition() != G4AntiProton::AntiProton() &&
1375 modifiedOriginal.GetDefinition() != G4AntiNeutron::AntiNeutron() &&
1376 modifiedOriginal.GetDefinition() != G4AntiLambda::AntiLambda() &&
1377 modifiedOriginal.GetDefinition() != G4AntiSigmaPlus::AntiSigmaPlus() &&
1378 modifiedOriginal.GetDefinition() != G4AntiSigmaMinus::AntiSigmaMinus() &&
1379 modifiedOriginal.GetDefinition() != G4AntiXiZero::AntiXiZero() &&
1380 modifiedOriginal.GetDefinition() != G4AntiXiMinus::AntiXiMinus() &&
1381 modifiedOriginal.GetDefinition() != G4AntiOmegaMinus::AntiOmegaMinus();
1382
1383 if( antiTest && (
1384 currentParticle.GetDefinition() == aPiPlus ||
1385 currentParticle.GetDefinition() == aPiZero ||
1386 currentParticle.GetDefinition() == aPiMinus ) &&
1387 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
1388 ( G4UniformRand() <= atomicWeight/300.0 ) )
1389 {
1390 if (eAvailable > nucleonMass - piMass) {
1391 if( G4UniformRand() > atomicNumber/atomicWeight )
1392 currentParticle.SetDefinitionAndUpdateE( aNeutron );
1393 else
1394 currentParticle.SetDefinitionAndUpdateE( aProton );
1395 incidentHasChanged = true;
1396 }
1397 }
1398 if( antiTest && (
1399 targetParticle.GetDefinition() == aPiPlus ||
1400 targetParticle.GetDefinition() == aPiZero ||
1401 targetParticle.GetDefinition() == aPiMinus ) &&
1402 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
1403 ( G4UniformRand() <= atomicWeight/300.0 ) )
1404 {
1405 if (eAvailable > nucleonMass - piMass) {
1406 if( G4UniformRand() > atomicNumber/atomicWeight )
1407 targetParticle.SetDefinitionAndUpdateE( aNeutron );
1408 else
1409 targetParticle.SetDefinitionAndUpdateE( aProton );
1410 targetHasChanged = true;
1411 }
1412 }
1413 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1414 for( G4int i=0; i<vecLen; ++i )
1415 {
1416 if( antiTest && (
1417 vec[i]->GetDefinition() == aPiPlus ||
1418 vec[i]->GetDefinition() == aPiZero ||
1419 vec[i]->GetDefinition() == aPiMinus ) &&
1420 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
1421 ( G4UniformRand() <= atomicWeight/300.0 ) )
1422 {
1423 if (eAvailable > nucleonMass - piMass) {
1424 if( G4UniformRand() > atomicNumber/atomicWeight )
1425 vec[i]->SetDefinitionAndUpdateE( aNeutron );
1426 else
1427 vec[i]->SetDefinitionAndUpdateE( aProton );
1428 }
1429 }
1430 }
1431 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1432 }
G4AntiOmegaMinus::AntiOmegaMinus
static G4AntiOmegaMinus * AntiOmegaMinus()
Definition: G4AntiOmegaMinus.cc:110
G4AntiXiMinus::AntiXiMinus
static G4AntiXiMinus * AntiXiMinus()
Definition: G4AntiXiMinus.cc:105
G4AntiXiZero::AntiXiZero
static G4AntiXiZero * AntiXiZero()
Definition: G4AntiXiZero.cc:106

Referenced by G4InelasticInteraction::CalculateMomenta().

◆ TwoBody()

void G4ReactionDynamics::TwoBody ( G4FastVector< G4ReactionProduct, GHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4DynamicParticle *  originalTarget,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  targetHasChanged 
)

Definition at line 2243 of file G4ReactionDynamics.cc.

2252 {
2253 //
2254 // derived from original FORTRAN code TWOB by H. Fesefeldt (15-Sep-1987)
2255 //
2256 // Generation of momenta for elastic and quasi-elastic 2 body reactions
2257 //
2258 // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
2259 // The b values are parametrizations from experimental data.
2260 // Not available values are taken from those of similar reactions.
2261 //
2262
2263 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2264 static const G4double expxu = 82.; // upper bound for arg. of exp
2265 static const G4double expxl = -expxu; // lower bound for arg. of exp
2266
2267 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
2268 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
2269 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
2270 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
2271 G4double currentMass = currentParticle.GetMass()/GeV;
2272 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
2273
2274 targetMass = targetParticle.GetMass()/GeV;
2275 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
2276
2277 G4double etCurrent = currentParticle.GetTotalEnergy()/GeV;
2278 G4double pCurrent = currentParticle.GetTotalMomentum()/GeV;
2279
2280 G4double cmEnergy = std::sqrt( currentMass*currentMass +
2281 targetMass*targetMass +
2282 2.0*targetMass*etCurrent ); // in GeV
2283
2284 //if( (pOriginal < 0.1) ||
2285 // (centerofmassEnergy < 0.01) ) // 2-body scattering not possible
2286 // Continue with original particle, but spend the nuclear evaporation energy
2287 // targetParticle.SetMass( 0.0 ); // flag that the target doesn't exist
2288 //else // Two-body scattering is possible
2289
2290 if( (pCurrent < 0.1) || (cmEnergy < 0.01) ) // 2-body scattering not possible
2291 {
2292 targetParticle.SetMass( 0.0 ); // flag that the target particle doesn't exist
2293 }
2294 else
2295 {
2296// moved this if-block to a later stage, i.e. to the assignment of the scattering angle
2297// @@@@@ double-check.
2298// if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
2299// targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
2300// if( G4UniformRand() < 0.5 )
2301// targetParticle.SetDefinitionAndUpdateE( aNeutron );
2302// else
2303// targetParticle.SetDefinitionAndUpdateE( aProton );
2304// targetHasChanged = true;
2305// targetMass = targetParticle.GetMass()/GeV;
2306// }
2307 //
2308 // Set masses and momenta for final state particles
2309 //
2310 G4double pf = cmEnergy*cmEnergy + targetMass*targetMass - currentMass*currentMass;
2311 pf = pf*pf - 4*cmEnergy*cmEnergy*targetMass*targetMass;
2312
2313 if( pf < 0.001 )
2314 {
2315 for(G4int i=0; i<vecLen; i++) delete vec[i];
2316 vecLen = 0;
2317 throw G4HadronicException(__FILE__, __LINE__, "G4ReactionDynamics::TwoBody: pf is too small ");
2318 }
2319
2320 pf = std::sqrt( pf ) / ( 2.0*cmEnergy );
2321 //
2322 // Set beam and target in centre of mass system
2323 //
2324 G4ReactionProduct pseudoParticle[3];
2325
2326 if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
2327 targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
2328 pseudoParticle[0].SetMass( targetMass*GeV );
2329 pseudoParticle[0].SetTotalEnergy( etOriginal*GeV );
2330 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
2331
2332 pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
2333 pseudoParticle[1].SetMass( mOriginal*GeV );
2334 pseudoParticle[1].SetKineticEnergy( 0.0 );
2335
2336 } else {
2337 pseudoParticle[0].SetMass( currentMass*GeV );
2338 pseudoParticle[0].SetTotalEnergy( etCurrent*GeV );
2339 pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV );
2340
2341 pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
2342 pseudoParticle[1].SetMass( targetMass*GeV );
2343 pseudoParticle[1].SetKineticEnergy( 0.0 );
2344 }
2345 //
2346 // Transform into centre of mass system
2347 //
2348 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
2349 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
2350 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
2351 //
2352 // Set final state masses and energies in centre of mass system
2353 //
2354 currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
2355 targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
2356 //
2357 // Set |t| and |tmin|
2358 //
2359 const G4double cb = 0.01;
2360 const G4double b1 = 4.225;
2361 const G4double b2 = 1.795;
2362 //
2363 // Calculate slope b for elastic scattering on proton/neutron
2364 //
2365 G4double b = std::max( cb, b1+b2*std::log(pOriginal) );
2366 G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV;
2367
2368 G4double exindt = -1.0;
2369 exindt += std::exp(std::max(-btrang,expxl));
2370 //
2371 // Calculate sqr(std::sin(teta/2.) and std::cos(teta), set azimuth angle phi
2372 //
2373 G4double ctet = 1.0 + 2*std::log( 1.0+G4UniformRand()*exindt ) / btrang;
2374 if( std::fabs(ctet) > 1.0 )ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
2375 G4double stet = std::sqrt( (1.0-ctet)*(1.0+ctet) );
2376 G4double phi = twopi * G4UniformRand();
2377 //
2378 // Calculate final state momenta in centre of mass system
2379 //
2380 if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
2381 targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
2382
2383 currentParticle.SetMomentum( -pf*stet*std::sin(phi)*GeV,
2384 -pf*stet*std::cos(phi)*GeV,
2385 -pf*ctet*GeV );
2386 } else {
2387
2388 currentParticle.SetMomentum( pf*stet*std::sin(phi)*GeV,
2389 pf*stet*std::cos(phi)*GeV,
2390 pf*ctet*GeV );
2391 }
2392 targetParticle.SetMomentum( currentParticle.GetMomentum() * (-1.0) );
2393 //
2394 // Transform into lab system
2395 //
2396 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
2397 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
2398
2399 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2400
2401 G4double pp, pp1, ekin;
2402 if( atomicWeight >= 1.5 )
2403 {
2404 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
2405 pp1 = currentParticle.GetMomentum().mag()/MeV;
2406 if( pp1 >= 1.0 )
2407 {
2408 ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
2409 ekin = std::max( 0.0001*GeV, ekin );
2410 currentParticle.SetKineticEnergy( ekin*MeV );
2411 pp = currentParticle.GetTotalMomentum()/MeV;
2412 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
2413 }
2414 pp1 = targetParticle.GetMomentum().mag()/MeV;
2415 if( pp1 >= 1.0 )
2416 {
2417 ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV;
2418 ekin = std::max( 0.0001*GeV, ekin );
2419 targetParticle.SetKineticEnergy( ekin*MeV );
2420 pp = targetParticle.GetTotalMomentum()/MeV;
2421 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
2422 }
2423 }
2424 }
2425
2426 // Get number of final state nucleons and nucleons remaining in
2427 // target nucleus
2428
2429 std::pair<G4int, G4int> finalStateNucleons =
2430 GetFinalStateNucleons(originalTarget, vec, vecLen);
2431 G4int protonsInFinalState = finalStateNucleons.first;
2432 G4int neutronsInFinalState = finalStateNucleons.second;
2433
2434 G4int PinNucleus = std::max(0,
2435 targetNucleus.GetZ_asInt() - protonsInFinalState);
2436 G4int NinNucleus = std::max(0,
2437 targetNucleus.GetN_asInt() - neutronsInFinalState);
2438
2439 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2440 if( atomicWeight >= 1.5 )
2441 {
2442 // Add black track particles
2443 // npnb is number of proton/neutron black track particles
2444 // ndta is the number of deuterons, tritons, and alphas produced
2445 // epnb is the kinetic energy available for proton/neutron black track particles
2446 // edta is the kinetic energy available for deuteron/triton/alpha particles
2447 //
2448 G4double epnb, edta;
2449 G4int npnb=0, ndta=0;
2450
2451 epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code
2452 edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code
2453 const G4double pnCutOff = 0.0001; // GeV
2454 const G4double dtaCutOff = 0.0001; // GeV
2455 const G4double kineticMinimum = 0.0001;
2456 const G4double kineticFactor = -0.010;
2457 G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
2458 if( epnb >= pnCutOff )
2459 {
2460 npnb = Poisson( epnb/0.02 );
2461 if( npnb > atomicWeight )npnb = G4int(atomicWeight);
2462 if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
2463 npnb = std::min( npnb, 127-vecLen );
2464 }
2465 if( edta >= dtaCutOff )
2466 {
2467 ndta = G4int(2.0 * std::log(atomicWeight));
2468 ndta = std::min( ndta, 127-vecLen );
2469 }
2470
2471 if (npnb == 0 && ndta == 0) npnb = 1;
2472
2473 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2474
2475 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
2476 kineticFactor, modifiedOriginal,
2477 PinNucleus, NinNucleus, targetNucleus,
2478 vec, vecLen);
2479 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2480 }
2481 //
2482 // calculate time delay for nuclear reactions
2483 //
2484 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
2485 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
2486 else
2487 currentParticle.SetTOF( 1.0 );
2488 return;
2489 }
G4ReactionProduct::SetTOF
void SetTOF(const G4double t)
Definition: G4ReactionProduct.hh:150

Referenced by G4InelasticInteraction::CalculateMomenta().

◆ TwoCluster()

G4bool G4ReactionDynamics::TwoCluster ( G4FastVector< G4ReactionProduct, GHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4HadProjectile *  originalIncident,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4DynamicParticle *  originalTarget,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged,
G4bool  leadFlag,
G4ReactionProduct &  leadingStrangeParticle 
)

Definition at line 1434 of file G4ReactionDynamics.cc.

1447 {
1448 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1449 // derived from original FORTRAN code TWOCLU by H. Fesefeldt (11-Oct-1987)
1450 //
1451 // Generation of X- and PT- values for incident, target, and all secondary particles
1452 //
1453 // A simple two cluster model is used.
1454 // This should be sufficient for low energy interactions.
1455 //
1456
1457 // + debugging
1458 // raise(SIGSEGV);
1459 // - debugging
1460
1461 G4int i;
1462 G4ParticleDefinition *aProton = G4Proton::Proton();
1463 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1464 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1465 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1466 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
1467 G4bool veryForward = false;
1468
1469 const G4double protonMass = aProton->GetPDGMass()/MeV;
1470 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
1471 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
1472 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
1473 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
1474 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
1475 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
1476 targetMass*targetMass +
1477 2.0*targetMass*etOriginal ); // GeV
1478 G4double currentMass = currentParticle.GetMass()/GeV;
1479 targetMass = targetParticle.GetMass()/GeV;
1480
1481 if( currentMass == 0.0 && targetMass == 0.0 )
1482 {
1483 G4double ek = currentParticle.GetKineticEnergy();
1484 G4ThreeVector mom = currentParticle.GetMomentum();
1485 currentParticle = *vec[0];
1486 targetParticle = *vec[1];
1487 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
1488 if(vecLen<2)
1489 {
1490 for(G4int j=0; j<vecLen; j++) delete vec[j];
1491 vecLen = 0;
1492 throw G4HadReentrentException(__FILE__, __LINE__,
1493 "G4ReactionDynamics::TwoCluster: Negative number of particles");
1494 }
1495 delete vec[vecLen-1];
1496 delete vec[vecLen-2];
1497 vecLen -= 2;
1498 currentMass = currentParticle.GetMass()/GeV;
1499 targetMass = targetParticle.GetMass()/GeV;
1500 incidentHasChanged = true;
1501 targetHasChanged = true;
1502 currentParticle.SetKineticEnergy( ek );
1503 currentParticle.SetMomentum( mom );
1504 veryForward = true;
1505 }
1506
1507 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
1508 const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
1509 //
1510 // particles have been distributed in forward and backward hemispheres
1511 // in center of mass system of the hadron nucleon interaction
1512 //
1513 // incident is always in forward hemisphere
1514 G4int forwardCount = 1; // number of particles in forward hemisphere
1515 currentParticle.SetSide( 1 );
1516 G4double forwardMass = currentParticle.GetMass()/GeV;
1517 G4double cMass = forwardMass;
1518
1519 // target is always in backward hemisphere
1520 G4int backwardCount = 1; // number of particles in backward hemisphere
1521 G4int backwardNucleonCount = 1; // number of nucleons in backward hemisphere
1522 targetParticle.SetSide( -1 );
1523 G4double backwardMass = targetParticle.GetMass()/GeV;
1524 G4double bMass = backwardMass;
1525
1526 for( i=0; i<vecLen; ++i )
1527 {
1528 if( vec[i]->GetSide() < 0 )vec[i]->SetSide( -1 ); // added by JLC, 2Jul97
1529 // to take care of the case where vec has been preprocessed by GenerateXandPt
1530 // and some of them have been set to -2 or -3
1531 if( vec[i]->GetSide() == -1 )
1532 {
1533 ++backwardCount;
1534 backwardMass += vec[i]->GetMass()/GeV;
1535 }
1536 else
1537 {
1538 ++forwardCount;
1539 forwardMass += vec[i]->GetMass()/GeV;
1540 }
1541 }
1542 //
1543 // nucleons and some pions from intranuclear cascade
1544 //
1545 G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
1546 if(term1 < 0) term1 = 0.0001; // making sure xtarg<0;
1547 const G4double afc = 0.312 + 0.2 * std::log(term1);
1548 G4double xtarg;
1549 if( centerofmassEnergy < 2.0+G4UniformRand() ) // added +2 below, JLC 4Jul97
1550 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
1551 else
1552 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
1553 if( xtarg <= 0.0 )xtarg = 0.01;
1554 G4int nuclearExcitationCount = Poisson( xtarg );
1555 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
1556 G4int extraNucleonCount = 0;
1557 G4double extraMass = 0.0;
1558 G4double extraNucleonMass = 0.0;
1559 if( nuclearExcitationCount > 0 )
1560 {
1561 G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );
1562 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
1563 //
1564 // NOTE: in TWOCLU, these new particles were given negative codes
1565 // here we use NewlyAdded = true instead
1566 //
1567// G4ReactionProduct *pVec = new G4ReactionProduct [nuclearExcitationCount];
1568 for( i=0; i<nuclearExcitationCount; ++i )
1569 {
1570 G4ReactionProduct* pVec = new G4ReactionProduct();
1571 if( G4UniformRand() < nucsup[momentumBin] ) // add proton or neutron
1572 {
1573 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
1574 pVec->SetDefinition( aProton );
1575 else
1576 pVec->SetDefinition( aNeutron );
1577 ++backwardNucleonCount;
1578 ++extraNucleonCount;
1579 extraNucleonMass += pVec->GetMass()/GeV;
1580 }
1581 else
1582 { // add a pion
1583 G4double ran = G4UniformRand();
1584 if( ran < 0.3181 )
1585 pVec->SetDefinition( aPiPlus );
1586 else if( ran < 0.6819 )
1587 pVec->SetDefinition( aPiZero );
1588 else
1589 pVec->SetDefinition( aPiMinus );
1590 }
1591 pVec->SetSide( -2 ); // backside particle
1592 extraMass += pVec->GetMass()/GeV;
1593 pVec->SetNewlyAdded( true );
1594 vec.SetElement( vecLen++, pVec );
1595 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1596 }
1597 }
1598 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1599 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
1600 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
1601 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
1602 G4bool secondaryDeleted;
1603 G4double pMass;
1604
1605 while( eAvailable <= 0.0 ) // must eliminate a particle
1606 {
1607 secondaryDeleted = false;
1608 for( i=(vecLen-1); i>=0; --i )
1609 {
1610 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
1611 {
1612 pMass = vec[i]->GetMass()/GeV;
1613 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
1614 --forwardCount;
1615 forwardEnergy += pMass;
1616 forwardMass -= pMass;
1617 secondaryDeleted = true;
1618 break;
1619 }
1620 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
1621 {
1622 pMass = vec[i]->GetMass()/GeV;
1623 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
1624 --backwardCount;
1625 backwardEnergy += pMass;
1626 backwardMass -= pMass;
1627 secondaryDeleted = true;
1628 break;
1629 }
1630 } // breaks go down to here
1631 if( secondaryDeleted )
1632 {
1633 delete vec[vecLen-1];
1634 --vecLen;
1635 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1636 }
1637 else
1638 {
1639 if( vecLen == 0 )
1640 {
1641 return false; // all secondaries have been eliminated
1642 }
1643 if( targetParticle.GetSide() == -1 )
1644 {
1645 pMass = targetParticle.GetMass()/GeV;
1646 targetParticle = *vec[0];
1647 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
1648 --backwardCount;
1649 backwardEnergy += pMass;
1650 backwardMass -= pMass;
1651 secondaryDeleted = true;
1652 }
1653 else if( targetParticle.GetSide() == 1 )
1654 {
1655 pMass = targetParticle.GetMass()/GeV;
1656 targetParticle = *vec[0];
1657 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
1658 --forwardCount;
1659 forwardEnergy += pMass;
1660 forwardMass -= pMass;
1661 secondaryDeleted = true;
1662 }
1663 if( secondaryDeleted )
1664 {
1665 delete vec[vecLen-1];
1666 --vecLen;
1667 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1668 }
1669 else
1670 {
1671 if( currentParticle.GetSide() == -1 )
1672 {
1673 pMass = currentParticle.GetMass()/GeV;
1674 currentParticle = *vec[0];
1675 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
1676 --backwardCount;
1677 backwardEnergy += pMass;
1678 backwardMass -= pMass;
1679 secondaryDeleted = true;
1680 }
1681 else if( currentParticle.GetSide() == 1 )
1682 {
1683 pMass = currentParticle.GetMass()/GeV;
1684 currentParticle = *vec[0];
1685 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
1686 --forwardCount;
1687 forwardEnergy += pMass;
1688 forwardMass -= pMass;
1689 secondaryDeleted = true;
1690 }
1691 if( secondaryDeleted )
1692 {
1693 delete vec[vecLen-1];
1694 --vecLen;
1695 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1696 }
1697 else break;
1698 }
1699 }
1700 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
1701 }
1702 //
1703 // This is the start of the TwoCluster function
1704 // Choose masses for the 3 clusters:
1705 // forward cluster
1706 // backward meson cluster
1707 // backward nucleon cluster
1708 //
1709 G4double rmc = 0.0, rmd = 0.0;
1710 const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
1711 const G4double gpar[] = { 2.6, 2.6, 1.8, 1.30, 1.20 };
1712
1713 if (forwardCount <= 0 || backwardCount <= 0) return false; // array bounds protection
1714
1715 if (forwardCount == 1) rmc = forwardMass;
1716 else
1717 {
1718 G4int ntc = std::max(1, std::min(5,forwardCount))-1; // check if offset by 1 @@
1719 rmc = forwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
1720 }
1721
1722 if (backwardCount == 1) rmd = backwardMass;
1723 else
1724 {
1725 G4int ntc = std::max(1, std::min(5,backwardCount)); // check, if offfset by 1 @@
1726 rmd = backwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
1727 }
1728
1729 while( rmc+rmd > centerofmassEnergy )
1730 {
1731 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
1732 {
1733 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
1734 rmc *= temp;
1735 rmd *= temp;
1736 }
1737 else
1738 {
1739 rmc = 0.1*forwardMass + 0.9*rmc;
1740 rmd = 0.1*backwardMass + 0.9*rmd;
1741 }
1742 }
1743
1744 G4ReactionProduct pseudoParticle[8];
1745 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
1746
1747 pseudoParticle[1].SetMass( mOriginal*GeV );
1748 pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
1749 pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
1750
1751 pseudoParticle[2].SetMass( protonMass*MeV );
1752 pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
1753 pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
1754 //
1755 // transform into centre of mass system
1756 //
1757 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
1758 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
1759 pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
1760
1761 const G4double pfMin = 0.0001;
1762 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
1763 pf *= pf;
1764 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
1765 pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
1766 //
1767 // set final state masses and energies in centre of mass system
1768 //
1769 pseudoParticle[3].SetMass( rmc*GeV );
1770 pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
1771
1772 pseudoParticle[4].SetMass( rmd*GeV );
1773 pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
1774 //
1775 // set |T| and |TMIN|
1776 //
1777 const G4double bMin = 0.01;
1778 const G4double b1 = 4.0;
1779 const G4double b2 = 1.6;
1780
1781 G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
1782 G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
1783 G4double factor = 1.0 - std::exp(-dtb);
1784 G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
1785 costheta = std::max(-1.0, std::min(1.0, costheta));
1786 G4double sintheta = std::sqrt(1.0-costheta*costheta);
1787 G4double phi = G4UniformRand() * twopi;
1788
1789 //
1790 // calculate final state momenta in centre of mass system
1791 //
1792 pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
1793 pf*sintheta*std::sin(phi)*GeV,
1794 pf*costheta*GeV );
1795
1796 pseudoParticle[4].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
1797 //
1798 // simulate backward nucleon cluster in lab. system and transform in cms
1799 //
1800 G4double pp, pp1;
1801 if( nuclearExcitationCount > 0 )
1802 {
1803 const G4double ga = 1.2;
1804 G4double ekit1 = 0.04;
1805 G4double ekit2 = 0.6;
1806 if( ekOriginal <= 5.0 )
1807 {
1808 ekit1 *= ekOriginal*ekOriginal/25.0;
1809 ekit2 *= ekOriginal*ekOriginal/25.0;
1810 }
1811 const G4double a = (1.0-ga)/(std::pow(ekit2,(1.0-ga)) - std::pow(ekit1,(1.0-ga)));
1812 for( i=0; i<vecLen; ++i )
1813 {
1814 if( vec[i]->GetSide() == -2 )
1815 {
1816 G4double kineticE =
1817 std::pow( (G4UniformRand()*(1.0-ga)/a+std::pow(ekit1,(1.0-ga))), (1.0/(1.0-ga)) );
1818 vec[i]->SetKineticEnergy( kineticE*GeV );
1819 G4double vMass = vec[i]->GetMass()/MeV;
1820 G4double totalE = kineticE*GeV + vMass;
1821 pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
1822 G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
1823 G4double sint = std::sqrt( std::max( 0.0, (1.0-cost*cost) ) );
1824 phi = twopi*G4UniformRand();
1825 vec[i]->SetMomentum( pp*sint*std::sin(phi)*MeV,
1826 pp*sint*std::cos(phi)*MeV,
1827 pp*cost*MeV );
1828 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
1829 }
1830 }
1831 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1832 }
1833 //
1834 // fragmentation of forward cluster and backward meson cluster
1835 //
1836 currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
1837 currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
1838
1839 targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
1840 targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
1841
1842 pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
1843 pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
1844 pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
1845
1846 pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
1847 pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
1848 pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
1849
1850 G4double wgt;
1851 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1852 if( forwardCount > 1 ) // tempV will contain the forward particles
1853 {
1854 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
1855 tempV.Initialize( forwardCount );
1856 G4bool constantCrossSection = true;
1857 G4int tempLen = 0;
1858 if( currentParticle.GetSide() == 1 )
1859 tempV.SetElement( tempLen++, &currentParticle );
1860 if( targetParticle.GetSide() == 1 )
1861 tempV.SetElement( tempLen++, &targetParticle );
1862 for( i=0; i<vecLen; ++i )
1863 {
1864 if( vec[i]->GetSide() == 1 )
1865 {
1866 if( tempLen < 18 )
1867 tempV.SetElement( tempLen++, vec[i] );
1868 else
1869 {
1870 vec[i]->SetSide( -1 );
1871 continue;
1872 }
1873 }
1874 }
1875 if( tempLen >= 2 )
1876 {
1877 wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
1878 constantCrossSection, tempV, tempLen );
1879 if( currentParticle.GetSide() == 1 )
1880 currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
1881 if( targetParticle.GetSide() == 1 )
1882 targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
1883 for( i=0; i<vecLen; ++i )
1884 {
1885 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
1886 }
1887 }
1888 }
1889 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1890 if( backwardCount > 1 ) // tempV will contain the backward particles,
1891 { // but not those created from the intranuclear cascade
1892 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
1893 tempV.Initialize( backwardCount );
1894 G4bool constantCrossSection = true;
1895 G4int tempLen = 0;
1896 if( currentParticle.GetSide() == -1 )
1897 tempV.SetElement( tempLen++, &currentParticle );
1898 if( targetParticle.GetSide() == -1 )
1899 tempV.SetElement( tempLen++, &targetParticle );
1900 for( i=0; i<vecLen; ++i )
1901 {
1902 if( vec[i]->GetSide() == -1 )
1903 {
1904 if( tempLen < 18 )
1905 tempV.SetElement( tempLen++, vec[i] );
1906 else
1907 {
1908 vec[i]->SetSide( -2 );
1909 vec[i]->SetKineticEnergy( 0.0 );
1910 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
1911 continue;
1912 }
1913 }
1914 }
1915 if( tempLen >= 2 )
1916 {
1917 wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
1918 constantCrossSection, tempV, tempLen );
1919 if( currentParticle.GetSide() == -1 )
1920 currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
1921 if( targetParticle.GetSide() == -1 )
1922 targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
1923 for( i=0; i<vecLen; ++i )
1924 {
1925 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
1926 }
1927 }
1928 }
1929 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1930 //
1931 // Lorentz transformation in lab system
1932 //
1933 currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
1934 targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
1935 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
1936
1937 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1938 //
1939 // sometimes the leading strange particle is lost, set it back
1940 //
1941 G4bool dum = true;
1942 if( leadFlag )
1943 {
1944 // leadFlag will be true
1945 // iff original particle is at least as heavy as K+ and not a proton or
1946 // neutron AND if incident particle is at least as heavy as K+ and it is
1947 // not a proton or neutron leadFlag is set to the incident particle
1948 // or
1949 // target particle is at least as heavy as K+ and it is not a proton or
1950 // neutron leadFlag is set to the target particle
1951 //
1952 if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
1953 dum = false;
1954 else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
1955 dum = false;
1956 else
1957 {
1958 for( i=0; i<vecLen; ++i )
1959 {
1960 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
1961 {
1962 dum = false;
1963 break;
1964 }
1965 }
1966 }
1967 if( dum )
1968 {
1969 G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
1970 G4double ekin;
1971 if( ((leadMass < protonMass) && (targetParticle.GetMass()/MeV < protonMass)) ||
1972 ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
1973 {
1974 ekin = targetParticle.GetKineticEnergy()/GeV;
1975 pp1 = targetParticle.GetMomentum().mag()/MeV; // old momentum
1976 targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
1977 targetParticle.SetKineticEnergy( ekin*GeV );
1978 pp = targetParticle.GetTotalMomentum()/MeV; // new momentum
1979 if( pp1 < 1.0e-3 )
1980 {
1981 G4double costheta2 = 2.*G4UniformRand() - 1.;
1982 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
1983 G4double phi2 = twopi*G4UniformRand();
1984 targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
1985 pp*sintheta2*std::sin(phi2)*MeV,
1986 pp*costheta2*MeV ) ;
1987 }
1988 else
1989 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
1990
1991 targetHasChanged = true;
1992 }
1993 else
1994 {
1995 ekin = currentParticle.GetKineticEnergy()/GeV;
1996 pp1 = currentParticle.GetMomentum().mag()/MeV;
1997 currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
1998 currentParticle.SetKineticEnergy( ekin*GeV );
1999 pp = currentParticle.GetTotalMomentum()/MeV;
2000 if( pp1 < 1.0e-3 )
2001 {
2002 G4double costheta2 = 2.*G4UniformRand() - 1.;
2003 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2004 G4double phi2 = twopi*G4UniformRand();
2005 currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2006 pp*sintheta2*std::sin(phi2)*MeV,
2007 pp*costheta2*MeV ) ;
2008 }
2009 else
2010 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
2011
2012 incidentHasChanged = true;
2013 }
2014 }
2015 } // end of if( leadFlag )
2016
2017 // Get number of final state nucleons and nucleons remaining in
2018 // target nucleus
2019
2020 std::pair<G4int, G4int> finalStateNucleons =
2021 GetFinalStateNucleons(originalTarget, vec, vecLen);
2022
2023 G4int protonsInFinalState = finalStateNucleons.first;
2024 G4int neutronsInFinalState = finalStateNucleons.second;
2025
2026 G4int numberofFinalStateNucleons =
2027 protonsInFinalState + neutronsInFinalState;
2028
2029 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
2030 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
2031 originalIncident->GetDefinition()->GetPDGMass() <
2032 G4Lambda::Lambda()->GetPDGMass())
2033 numberofFinalStateNucleons++;
2034
2035 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
2036
2037 G4int PinNucleus = std::max(0,
2038 targetNucleus.GetZ_asInt() - protonsInFinalState);
2039 G4int NinNucleus = std::max(0,
2040 targetNucleus.GetN_asInt() - neutronsInFinalState);
2041 //
2042 // for various reasons, the energy balance is not sufficient,
2043 // check that, energy balance, angle of final system, etc.
2044 //
2045 pseudoParticle[4].SetMass( mOriginal*GeV );
2046 pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
2047 pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
2048
2049 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
2050 G4int diff = 0;
2051 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
2052 if(numberofFinalStateNucleons == 1) diff = 0;
2053 pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
2054 pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
2055 pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
2056
2057 G4double theoreticalKinetic =
2058 pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
2059
2060 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
2061 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
2062 pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
2063
2064 if( vecLen < 16 )
2065 {
2066 G4ReactionProduct tempR[130];
2067 tempR[0] = currentParticle;
2068 tempR[1] = targetParticle;
2069 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
2070
2071 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
2072 tempV.Initialize( vecLen+2 );
2073 G4bool constantCrossSection = true;
2074 G4int tempLen = 0;
2075 for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
2076
2077 if( tempLen >= 2 )
2078 {
2079 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2080 wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
2081 pseudoParticle[5].GetTotalEnergy()/MeV,
2082 constantCrossSection, tempV, tempLen );
2083 if (wgt == -1) {
2084 G4double Qvalue = 0;
2085 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
2086 wgt = GenerateNBodyEvent( Qvalue/MeV,
2087 constantCrossSection, tempV, tempLen );
2088 }
2089 theoreticalKinetic = 0.0;
2090 for( i=0; i<vecLen+2; ++i )
2091 {
2092 pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
2093 pseudoParticle[7].SetMass( tempV[i]->GetMass() );
2094 pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
2095 pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
2096 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
2097 }
2098 }
2099 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2100 }
2101 else
2102 {
2103 theoreticalKinetic -=
2104 ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
2105 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
2106 }
2107 G4double simulatedKinetic =
2108 currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
2109 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
2110 //
2111 // make sure that kinetic energies are correct
2112 // the backward nucleon cluster is not produced within proper kinematics!!!
2113 //
2114
2115 if( simulatedKinetic != 0.0 )
2116 {
2117 wgt = (theoreticalKinetic)/simulatedKinetic;
2118 currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
2119 pp = currentParticle.GetTotalMomentum()/MeV;
2120 pp1 = currentParticle.GetMomentum().mag()/MeV;
2121 if( pp1 < 0.001*MeV )
2122 {
2123 G4double costheta2 = 2.*G4UniformRand() - 1.;
2124 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2125 G4double phi2 = twopi*G4UniformRand();
2126 currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2127 pp*sintheta2*std::sin(phi2)*MeV,
2128 pp*costheta2*MeV ) ;
2129 }
2130 else
2131 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
2132
2133 targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
2134 pp = targetParticle.GetTotalMomentum()/MeV;
2135 pp1 = targetParticle.GetMomentum().mag()/MeV;
2136 if( pp1 < 0.001*MeV )
2137 {
2138 G4double costheta2 = 2.*G4UniformRand() - 1.;
2139 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2140 G4double phi2 = twopi*G4UniformRand();
2141 targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2142 pp*sintheta2*std::sin(phi2)*MeV,
2143 pp*costheta2*MeV ) ;
2144 }
2145 else
2146 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
2147
2148 for( i=0; i<vecLen; ++i )
2149 {
2150 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
2151 pp = vec[i]->GetTotalMomentum()/MeV;
2152 pp1 = vec[i]->GetMomentum().mag()/MeV;
2153 if( pp1 < 0.001 )
2154 {
2155 G4double costheta2 = 2.*G4UniformRand() - 1.;
2156 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2157 G4double phi2 = twopi*G4UniformRand();
2158 vec[i]->SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2159 pp*sintheta2*std::sin(phi2)*MeV,
2160 pp*costheta2*MeV ) ;
2161 }
2162 else
2163 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
2164 }
2165 }
2166 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2167
2168 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
2169 modifiedOriginal, originalIncident, targetNucleus,
2170 currentParticle, targetParticle, vec, vecLen );
2171 //
2172 // add black track particles
2173 // the total number of particles produced is restricted to 198
2174 // this may have influence on very high energies
2175 //
2176 if( atomicWeight >= 1.5 )
2177 {
2178 // npnb is number of proton/neutron black track particles
2179 // ndta is the number of deuterons, tritons, and alphas produced
2180 // epnb is the kinetic energy available for proton/neutron black track
2181 // particles
2182 // edta is the kinetic energy available for deuteron/triton/alpha
2183 // particles
2184
2185 G4int npnb = 0;
2186 G4int ndta = 0;
2187
2188 G4double epnb, edta;
2189 if (veryForward) {
2190 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
2191 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
2192 } else {
2193 epnb = targetNucleus.GetPNBlackTrackEnergy();
2194 edta = targetNucleus.GetDTABlackTrackEnergy();
2195 }
2196
2197 const G4double pnCutOff = 0.001; // GeV
2198 const G4double dtaCutOff = 0.001; // GeV
2199 const G4double kineticMinimum = 1.e-6;
2200 const G4double kineticFactor = -0.005;
2201
2202 G4double sprob = 0.0; // sprob = probability of self-absorption in
2203 // heavy molecules
2204 const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
2205 if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
2206
2207 if( epnb >= pnCutOff )
2208 {
2209 npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
2210 if( numberofFinalStateNucleons + npnb > atomicWeight )
2211 npnb = G4int(atomicWeight - numberofFinalStateNucleons);
2212 npnb = std::min( npnb, 127-vecLen );
2213 }
2214 if( edta >= dtaCutOff )
2215 {
2216 ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
2217 ndta = std::min( ndta, 127-vecLen );
2218 }
2219 if (npnb == 0 && ndta == 0) npnb = 1;
2220
2221 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2222
2223 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
2224 kineticFactor, modifiedOriginal,
2225 PinNucleus, NinNucleus, targetNucleus,
2226 vec, vecLen );
2227 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2228 }
2229 //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
2230 // MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2231 //
2232 // calculate time delay for nuclear reactions
2233 //
2234 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
2235 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
2236 else
2237 currentParticle.SetTOF( 1.0 );
2238
2239 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2240 return true;
2241 }

Referenced by G4InelasticInteraction::CalculateMomenta().

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