Geant4.10
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
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 ( )
inline

Definition at line 52 of file G4ReactionDynamics.hh.

52 {}
virtual G4ReactionDynamics::~G4ReactionDynamics ( )
inlinevirtual

Definition at line 54 of file G4ReactionDynamics.hh.

54 {}

Member Function Documentation

G4int G4ReactionDynamics::Factorial ( G4int  n)

Definition at line 2743 of file G4ReactionDynamics.cc.

References G4INCL::Math::min().

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:78
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
virtual G4double G4ReactionDynamics::FindInelasticity ( )
inlinevirtual

Definition at line 56 of file G4ReactionDynamics.hh.

57  { return 0.0; }
virtual G4double G4ReactionDynamics::FindTimeDelay ( )
inlinevirtual

Definition at line 59 of file G4ReactionDynamics.hh.

60  { return 0.0; }
G4double G4ReactionDynamics::GenerateNBodyEvent ( const G4double  totalEnergy,
const G4bool  constantCrossSection,
G4FastVector< G4ReactionProduct, GHADLISTSIZE > &  vec,
G4int &  vecLen 
)

Definition at line 2491 of file G4ReactionDynamics.cc.

References test::a, test::b, test::c, energy(), G4cerr, G4endl, G4UniformRand, python.hepunit::GeV, G4INCL::Math::max(), G4INCL::Math::min(), G4InuclParticleNames::s0, G4InuclParticleNames::sm, and python.hepunit::twopi.

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

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  }
G4int
int G4int
Definition: G4Types.hh:78
energy
double precision function energy(A, Z)
Definition: dpm25nuc6.f:4106
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4bool
bool G4bool
Definition: G4Types.hh:79
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4cerr
G4GLOB_DLL std::ostream G4cerr
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.

References FatalException, G4cerr, G4endl, G4Exception(), G4UniformRand, GenerateNBodyEvent(), G4Nucleus::GetA_asInt(), G4Nucleus::GetAnnihilationDTABlackTrackEnergy(), G4Nucleus::GetAnnihilationPNBlackTrackEnergy(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4Nucleus::GetN_asInt(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetSide(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), python.hepunit::GeV, G4FastVector< Type, N >::Initialize(), G4Lambda::Lambda(), G4ReactionProduct::Lorentz(), CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), python.hepunit::MeV, G4INCL::Math::min(), G4Neutron::Neutron(), python.hepunit::pi, G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4InuclParticleNames::pp, G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetNewlyAdded(), G4ReactionProduct::SetSide(), G4ReactionProduct::SetTotalEnergy(), G4ReactionProduct::SetZero(), sqr(), python.hepunit::twopi, test::x, CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

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  }
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4FastVector::SetElement
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:199
CLHEP::Hep3Vector::x
double x() const
G4Nucleus::GetAnnihilationPNBlackTrackEnergy
G4double GetAnnihilationPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:156
G4ReactionProduct::GetSide
G4int GetSide() const
Definition: G4ReactionProduct.hh:162
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:165
G4ReactionProduct::SetSide
void SetSide(const G4int sid)
Definition: G4ReactionProduct.hh:159
G4Nucleus::GetDTABlackTrackEnergy
G4double GetDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:153
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:181
G4FastVector::Initialize
void Initialize(G4int items)
Definition: G4FastVector.hh:63
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4ReactionProduct::SetNewlyAdded
void SetNewlyAdded(const G4bool f)
Definition: G4ReactionProduct.hh:165
G4ReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen)
Definition: G4ReactionDynamics.cc:2491
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:147
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4Nucleus::GetN_asInt
G4int GetN_asInt() const
Definition: G4Nucleus.hh:112
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4Nucleus::GetAnnihilationDTABlackTrackEnergy
G4double GetAnnihilationDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:159
G4bool
bool G4bool
Definition: G4Types.hh:79
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
CLHEP::Hep3Vector::y
double y() const
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
FatalException
Definition: G4ExceptionSeverity.hh:60
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:154
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4endl
#define G4endl
Definition: G4ios.hh:61
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4Nucleus::GetPNBlackTrackEnergy
G4double GetPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:150
G4cerr
G4GLOB_DLL std::ostream G4cerr
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.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4UniformRand, G4Gamma::Gamma(), GenerateNBodyEvent(), G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalMomentum(), python.hepunit::GeV, G4FastVector< Type, N >::Initialize(), G4ReactionProduct::Lorentz(), CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), python.hepunit::MeV, G4INCL::Math::min(), G4Neutron::Neutron(), G4InuclParticleNames::pp, G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTotalEnergy(), G4Triton::Triton(), python.hepunit::twopi, and test::v.

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  }
G4FastVector::SetElement
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:199
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:165
p
const char * p
Definition: xmltok.h:285
G4FastVector::Initialize
void Initialize(G4int items)
Definition: G4FastVector.hh:63
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4ReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen)
Definition: G4ReactionDynamics.cc:2491
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:147
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4bool
bool G4bool
Definition: G4Types.hh:79
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:154
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
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.

References G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4UniformRand, G4ReactionProduct::GetDefinition(), G4DynamicParticle::GetDefinition(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), python.hepunit::GeV, G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4INCL::Math::max(), G4Neutron::Neutron(), G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetMayBeKilled(), G4ReactionProduct::SetSide(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), and G4SigmaZero::SigmaZero().

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  }
G4FastVector::SetElement
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
G4ReactionProduct::SetMayBeKilled
void SetMayBeKilled(const G4bool f)
Definition: G4ReactionProduct.hh:171
G4KaonZeroLong::KaonZeroLong
static G4KaonZeroLong * KaonZeroLong()
Definition: G4KaonZeroLong.cc:115
G4ReactionProduct::SetSide
void SetSide(const G4int sid)
Definition: G4ReactionProduct.hh:159
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4AntiSigmaPlus::AntiSigmaPlus
static G4AntiSigmaPlus * AntiSigmaPlus()
Definition: G4AntiSigmaPlus.cc:108
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::SetDefinitionAndUpdateE
void SetDefinitionAndUpdateE(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:142
G4SigmaZero::SigmaZero
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99
G4KaonMinus::KaonMinus
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4AntiSigmaMinus::AntiSigmaMinus
static G4AntiSigmaMinus * AntiSigmaMinus()
Definition: G4AntiSigmaMinus.cc:106
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4KaonZeroShort::KaonZeroShort
static G4KaonZeroShort * KaonZeroShort()
Definition: G4KaonZeroShort.cc:104
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4SigmaMinus::SigmaMinus
static G4SigmaMinus * SigmaMinus()
Definition: G4SigmaMinus.cc:106
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4AntiLambda::AntiLambda
static G4AntiLambda * AntiLambda()
Definition: G4AntiLambda.cc:108
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4AntiSigmaZero::AntiSigmaZero
static G4AntiSigmaZero * AntiSigmaZero()
Definition: G4AntiSigmaZero.cc:99
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:154
G4SigmaPlus::SigmaPlus
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4double
double G4double
Definition: G4Types.hh:76
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
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.

References G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4UniformRand, G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), python.hepunit::GeV, G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4Proton::Proton(), and G4ReactionProduct::SetDefinitionAndUpdateE().

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  }
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4AntiOmegaMinus::AntiOmegaMinus
static G4AntiOmegaMinus * AntiOmegaMinus()
Definition: G4AntiOmegaMinus.cc:110
G4AntiSigmaPlus::AntiSigmaPlus
static G4AntiSigmaPlus * AntiSigmaPlus()
Definition: G4AntiSigmaPlus.cc:108
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::SetDefinitionAndUpdateE
void SetDefinitionAndUpdateE(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:142
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4AntiSigmaMinus::AntiSigmaMinus
static G4AntiSigmaMinus * AntiSigmaMinus()
Definition: G4AntiSigmaMinus.cc:106
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4bool
bool G4bool
Definition: G4Types.hh:79
G4AntiXiMinus::AntiXiMinus
static G4AntiXiMinus * AntiXiMinus()
Definition: G4AntiXiMinus.cc:105
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4AntiLambda::AntiLambda
static G4AntiLambda * AntiLambda()
Definition: G4AntiLambda.cc:108
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4AntiXiZero::AntiXiZero
static G4AntiXiZero * AntiXiZero()
Definition: G4AntiXiZero.cc:106
G4double
double G4double
Definition: G4Types.hh:76
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
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.

References test::b, G4UniformRand, G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4Nucleus::GetN_asInt(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), python.hepunit::GeV, G4ReactionProduct::Lorentz(), CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), python.hepunit::MeV, G4INCL::Math::min(), G4InuclParticleNames::pp, G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTOF(), G4ReactionProduct::SetTotalEnergy(), and python.hepunit::twopi.

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  }
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:199
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:165
G4Nucleus::GetDTABlackTrackEnergy
G4double GetDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:153
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:181
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:147
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4Nucleus::GetN_asInt
G4int GetN_asInt() const
Definition: G4Nucleus.hh:112
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4ReactionProduct::SetTOF
void SetTOF(const G4double t)
Definition: G4ReactionProduct.hh:153
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4Nucleus::GetPNBlackTrackEnergy
G4double GetPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:150
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.

References test::a, G4UniformRand, GenerateNBodyEvent(), G4Nucleus::GetA_asInt(), G4Nucleus::GetAnnihilationDTABlackTrackEnergy(), G4Nucleus::GetAnnihilationPNBlackTrackEnergy(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4Nucleus::GetN_asInt(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetSide(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), python.hepunit::GeV, G4FastVector< Type, N >::Initialize(), G4Lambda::Lambda(), G4ReactionProduct::Lorentz(), CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), python.hepunit::MeV, G4INCL::Math::min(), G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4InuclParticleNames::pp, G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetNewlyAdded(), G4ReactionProduct::SetSide(), G4ReactionProduct::SetTOF(), G4ReactionProduct::SetTotalEnergy(), and python.hepunit::twopi.

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  }
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4FastVector::SetElement
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:199
G4Nucleus::GetAnnihilationPNBlackTrackEnergy
G4double GetAnnihilationPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:156
G4ReactionProduct::GetSide
G4int GetSide() const
Definition: G4ReactionProduct.hh:162
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:165
G4ReactionProduct::SetSide
void SetSide(const G4int sid)
Definition: G4ReactionProduct.hh:159
G4Nucleus::GetDTABlackTrackEnergy
G4double GetDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:153
G4FastVector::Initialize
void Initialize(G4int items)
Definition: G4FastVector.hh:63
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4ReactionProduct::SetNewlyAdded
void SetNewlyAdded(const G4bool f)
Definition: G4ReactionProduct.hh:165
G4ReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen)
Definition: G4ReactionDynamics.cc:2491
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:147
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4Nucleus::GetN_asInt
G4int GetN_asInt() const
Definition: G4Nucleus.hh:112
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4Nucleus::GetAnnihilationDTABlackTrackEnergy
G4double GetAnnihilationDTABlackTrackEnergy() const
Definition: G4Nucleus.hh:159
G4bool
bool G4bool
Definition: G4Types.hh:79
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:154
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4ReactionProduct::SetTOF
void SetTOF(const G4double t)
Definition: G4ReactionProduct.hh:153
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4Nucleus::GetPNBlackTrackEnergy
G4double GetPNBlackTrackEnergy() const
Definition: G4Nucleus.hh:150
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:57:00 for Geant4.10 by   doxygen 1.8.7