Geant4  10.02.p01
G4RPGInelastic.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4RPGInelastic.cc 94214 2015-11-09 08:18:05Z gcosmo $
27 //
28 
29 #include "G4RPGInelastic.hh"
30 #include "G4Log.hh"
31 #include "Randomize.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
36 #include "G4RPGTwoBody.hh"
37 
38 
40  : G4HadronicInteraction(modelName)
41 {
42  cache = 0.0;
61 
62  G4cout << " **************************************************** " << G4endl;
63  G4cout << " * The RPG model is currently under development and * " << G4endl;
64  G4cout << " * should not be used. * " << G4endl;
65  G4cout << " **************************************************** " << G4endl;
66 }
67 
68 
70  G4int n, G4double b, G4double c)
71 {
72  const G4double expxu = 82.; // upper bound for arg. of exp
73  const G4double expxl = -expxu; // lower bound for arg. of exp
74  G4double npf = 0.0;
75  G4double nmf = 0.0;
76  G4double nzf = 0.0;
77  G4int i;
78  for( i=2; i<=np; i++ )npf += G4Log((double)i);
79  for( i=2; i<=nneg; i++ )nmf += G4Log((double)i);
80  for( i=2; i<=nz; i++ )nzf += G4Log((double)i);
81  G4double r;
82  r = std::min( expxu, std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
83  return G4Exp(r);
84 }
85 
86 
88 {
89  G4int j = std::min(n,10);
90  G4int result = 1;
91  if (j <= 1) return result;
92  for (G4int i = 2; i <= j; ++i) result *= i;
93  return result;
94 }
95 
96 
98  const G4ReactionProduct &currentParticle,
99  const G4ReactionProduct &targetParticle,
100  G4ReactionProduct &leadParticle )
101 {
102  // The following was in GenerateXandPt and TwoCluster.
103  // Add a parameter to the GenerateXandPt function telling it about the
104  // strange particle.
105  //
106  // Assumes that the original particle was a strange particle
107  //
108  G4bool lead = false;
109  if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
110  (currentParticle.GetDefinition() != G4Proton::Proton()) &&
111  (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
112  {
113  lead = true;
114  leadParticle = currentParticle; // set lead to the incident particle
115  }
116  else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
117  (targetParticle.GetDefinition() != G4Proton::Proton()) &&
118  (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
119  {
120  lead = true;
121  leadParticle = targetParticle; // set lead to the target particle
122  }
123  return lead;
124 }
125 
126 
127  void G4RPGInelastic::SetUpPions(const G4int np, const G4int nneg,
128  const G4int nz,
130  G4int &vecLen)
131  {
132  if( np+nneg+nz == 0 )return;
133  G4int i;
135  for( i=0; i<np; ++i )
136  {
137  p = new G4ReactionProduct;
139  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
140  vec.SetElement( vecLen++, p );
141  }
142  for( i=np; i<np+nneg; ++i )
143  {
144  p = new G4ReactionProduct;
146  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
147  vec.SetElement( vecLen++, p );
148  }
149  for( i=np+nneg; i<np+nneg+nz; ++i )
150  {
151  p = new G4ReactionProduct;
153  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
154  vec.SetElement( vecLen++, p );
155  }
156  }
157 
158 
160  const G4double energy, // MeV, <0 means annihilation channels
161  G4double &n,
162  G4double &anpn )
163  {
164  const G4double expxu = 82.; // upper bound for arg. of exp
165  const G4double expxl = -expxu; // lower bound for arg. of exp
166  const G4int numSec = 60;
167  //
168  // the only difference between the calculation for annihilation channels
169  // and normal is the starting value, iBegin, for the loop below
170  //
171  G4int iBegin = 1;
172  G4double en = energy;
173  if( energy < 0.0 )
174  {
175  iBegin = 2;
176  en *= -1.0;
177  }
178  //
179  // number of total particles vs. centre of mass Energy - 2*proton mass
180  //
181  G4double aleab = G4Log(en/GeV);
182  n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
183  n -= 2.0;
184  //
185  // normalization constant for kno-distribution
186  //
187  anpn = 0.0;
188  G4double test, temp;
189  for( G4int i=iBegin; i<=numSec; ++i )
190  {
191  temp = pi*i/(2.0*n*n);
192  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
193  if( temp < 1.0 )
194  {
195  if( test >= 1.0e-10 )anpn += temp*test;
196  }
197  else
198  anpn += temp*test;
199  }
200  }
201 
202 void
204  G4int& vecLen,
205  const G4HadProjectile* originalIncident,
206  const G4DynamicParticle* originalTarget,
207  G4ReactionProduct& modifiedOriginal,
208  G4Nucleus& targetNucleus,
209  G4ReactionProduct& currentParticle,
210  G4ReactionProduct& targetParticle,
211  G4bool& incidentHasChanged,
212  G4bool& targetHasChanged,
213  G4bool quasiElastic)
214 {
215  cache = 0;
216  what = originalIncident->Get4Momentum().vect();
217 
218  G4ReactionProduct leadingStrangeParticle;
219 
220  // strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
221  // incidentHasChanged, originalTarget,
222  // targetParticle, targetHasChanged,
223  // targetNucleus, currentParticle,
224  // vec, vecLen,
225  // false, leadingStrangeParticle);
226 
227  if( quasiElastic )
228  {
229  twoBody.ReactionStage(originalIncident, modifiedOriginal,
230  incidentHasChanged, originalTarget,
231  targetParticle, targetHasChanged,
232  targetNucleus, currentParticle,
233  vec, vecLen,
234  false, leadingStrangeParticle);
235  return;
236  }
237 
238  G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
239  targetParticle,
240  leadingStrangeParticle );
241  //
242  // Note: the number of secondaries can be reduced in GenerateXandPt
243  // and TwoCluster
244  //
245  G4bool finishedGenXPt = false;
246  G4bool annihilation = false;
247  if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
248  currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
249  {
250  // original was an anti-particle and annihilation has taken place
251  annihilation = true;
252  G4double ekcor = 1.0;
253  G4double ek = originalIncident->GetKineticEnergy();
254  G4double ekOrg = ek;
255 
256  const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
257  if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
258  const G4double atomicWeight = targetNucleus.GetA_asInt();
259  ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
260  G4double tkin = targetNucleus.Cinema(ek);
261  ek += tkin;
262  ekOrg += tkin;
263  // modifiedOriginal.SetKineticEnergy( ekOrg );
264  //
265  // evaporation -- re-calculate black track energies
266  // this was Done already just before the cascade
267  //
268  tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
269  ekOrg -= tkin;
270  ekOrg = std::max( 0.0001*GeV, ekOrg );
271  modifiedOriginal.SetKineticEnergy( ekOrg );
272  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
273  G4double et = ekOrg + amas;
274  G4double p = std::sqrt( std::abs(et*et-amas*amas) );
275  G4double pp = modifiedOriginal.GetMomentum().mag();
276  if( pp > 0.0 )
277  {
278  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
279  modifiedOriginal.SetMomentum( momentum * (p/pp) );
280  }
281  if( ekOrg <= 0.0001 )
282  {
283  modifiedOriginal.SetKineticEnergy( 0.0 );
284  modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
285  }
286  }
287 
288  // twsup gives percentage of time two-cluster model is called
289 
290  const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
291  G4double rand1 = G4UniformRand();
292  G4double rand2 = G4UniformRand();
293 
294  // Cache current, target, and secondaries
295  G4ReactionProduct saveCurrent = currentParticle;
296  G4ReactionProduct saveTarget = targetParticle;
297  std::vector<G4ReactionProduct> savevec;
298  for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
299 
300  // Call fragmentation code if
301  // 1) there is annihilation, or
302  // 2) there are more than 5 secondaries, or
303  // 3) incident KE is > 1 GeV AND
304  // ( incident is a kaon AND rand < 0.5 OR twsup )
305  //
306 
307  if( annihilation || vecLen > 5 ||
308  ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
309 
310  (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
311  originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
312  originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
313  originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
314  rand1 < 0.5)
315  || rand2 > twsup[vecLen]) ) )
316 
317  finishedGenXPt =
318  fragmentation.ReactionStage(originalIncident, modifiedOriginal,
319  incidentHasChanged, originalTarget,
320  targetParticle, targetHasChanged,
321  targetNucleus, currentParticle,
322  vec, vecLen,
323  leadFlag, leadingStrangeParticle);
324 
325  if (finishedGenXPt) return;
326 
327  G4bool finishedTwoClu = false;
328 
329  if (modifiedOriginal.GetTotalMomentum() < 1.0) {
330  for (G4int i = 0; i < vecLen; i++) delete vec[i];
331  vecLen = 0;
332 
333  } else {
334  // Occaisionally, GenerateXandPt will fail in the annihilation channel.
335  // Restore current, target and secondaries to pre-GenerateXandPt state
336  // before trying annihilation in TwoCluster
337 
338  if (!finishedGenXPt && annihilation) {
339  currentParticle = saveCurrent;
340  targetParticle = saveTarget;
341  for (G4int i = 0; i < vecLen; i++) delete vec[i];
342  vecLen = 0;
343  vec.Initialize( 0 );
344  for (G4int i = 0; i < G4int(savevec.size()); i++) {
346  *p = savevec[i];
347  vec.SetElement( vecLen++, p );
348  }
349  }
350 
351  // Big violations of energy conservation in this method - don't use
352  //
353  // pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
354  // incidentHasChanged, originalTarget,
355  // targetParticle, targetHasChanged,
356  // targetNucleus, currentParticle,
357  // vec, vecLen,
358  // false, leadingStrangeParticle);
359 
360  try
361  {
362  finishedTwoClu =
363  twoCluster.ReactionStage(originalIncident, modifiedOriginal,
364  incidentHasChanged, originalTarget,
365  targetParticle, targetHasChanged,
366  targetNucleus, currentParticle,
367  vec, vecLen,
368  leadFlag, leadingStrangeParticle);
369  }
370  catch(G4HadReentrentException aC)
371  {
372  aC.Report(G4cout);
373  throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
374  }
375  }
376 
377  if (finishedTwoClu) return;
378 
379  twoBody.ReactionStage(originalIncident, modifiedOriginal,
380  incidentHasChanged, originalTarget,
381  targetParticle, targetHasChanged,
382  targetNucleus, currentParticle,
383  vec, vecLen,
384  false, leadingStrangeParticle);
385 }
386 
387 /*
388  void G4RPGInelastic::
389  Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen)
390  {
391  G4double rotation = 2.*pi*G4UniformRand();
392  cache = rotation;
393  G4int i;
394  for( i=0; i<vecLen; ++i )
395  {
396  G4ThreeVector momentum = vec[i]->GetMomentum();
397  momentum = momentum.rotate(rotation, what);
398  vec[i]->SetMomentum(momentum);
399  }
400  }
401 */
402 
403 void
405  G4int& vecLen,
406  G4ReactionProduct& currentParticle,
407  G4ReactionProduct& targetParticle,
408  G4bool& incidentHasChanged )
409 {
413  G4int i;
414 
415  if (currentParticle.GetDefinition() == particleDef[k0]) {
416  if (G4UniformRand() < 0.5) {
417  currentParticle.SetDefinitionAndUpdateE(aKaonZL);
418  incidentHasChanged = true;
419  } else {
420  currentParticle.SetDefinitionAndUpdateE(aKaonZS);
421  }
422  } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
423  if (G4UniformRand() < 0.5) {
424  currentParticle.SetDefinitionAndUpdateE(aKaonZL);
425  } else {
426  currentParticle.SetDefinitionAndUpdateE(aKaonZS);
427  incidentHasChanged = true;
428  }
429  }
430 
431  if (targetParticle.GetDefinition() == particleDef[k0] ||
432  targetParticle.GetDefinition() == particleDef[k0b] ) {
433  if (G4UniformRand() < 0.5) {
434  targetParticle.SetDefinitionAndUpdateE(aKaonZL);
435  } else {
436  targetParticle.SetDefinitionAndUpdateE(aKaonZS);
437  }
438  }
439 
440  for (i = 0; i < vecLen; ++i) {
441  if (vec[i]->GetDefinition() == particleDef[k0] ||
442  vec[i]->GetDefinition() == particleDef[k0b] ) {
443  if (G4UniformRand() < 0.5) {
444  vec[i]->SetDefinitionAndUpdateE(aKaonZL);
445  } else {
446  vec[i]->SetDefinitionAndUpdateE(aKaonZS);
447  }
448  }
449  }
450 
451  if (incidentHasChanged) {
453  p0->SetDefinition(currentParticle.GetDefinition() );
454  p0->SetMomentum(currentParticle.GetMomentum() );
458 
459  } else {
460  G4double p = currentParticle.GetMomentum().mag()/MeV;
461  G4ThreeVector mom = currentParticle.GetMomentum();
462  if (p > DBL_MIN)
463  theParticleChange.SetMomentumChange(mom.x()/p, mom.y()/p, mom.z()/p );
464  else
465  theParticleChange.SetMomentumChange(0.0, 0.0, 1.0);
466 
467  G4double aE = currentParticle.GetKineticEnergy();
468  if (std::fabs(aE)<.1*eV) aE=.1*eV;
470  }
471 
472  if (targetParticle.GetMass() > 0.0) // Tgt particle can be eliminated in TwoBody
473  {
474  G4ThreeVector momentum = targetParticle.GetMomentum();
475  momentum = momentum.rotate(cache, what);
476  G4double targKE = targetParticle.GetKineticEnergy();
477  G4ThreeVector dir(0.0, 0.0, 1.0);
478  if (targKE < DBL_MIN)
479  targKE = DBL_MIN;
480  else
481  dir = momentum/momentum.mag();
482 
483  G4DynamicParticle* p1 =
484  new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
485 
487  }
488 
490  for (i = 0; i < vecLen; ++i) {
491  G4double secKE = vec[i]->GetKineticEnergy();
492  G4ThreeVector momentum = vec[i]->GetMomentum();
493  G4ThreeVector dir(0.0, 0.0, 1.0);
494  if (secKE < DBL_MIN)
495  secKE = DBL_MIN;
496  else
497  dir = momentum/momentum.mag();
498 
499  p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
501  delete vec[i];
502  }
503 }
504 
505 
506 std::pair<G4int, G4double>
508 {
509  G4int index = 29;
510  G4double fraction = 0.0;
511 
512  for (G4int i = 1; i < 30; i++) {
513  if (e < energyScale[i]) {
514  index = i-1;
515  fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
516  break;
517  }
518  }
519  return std::pair<G4int, G4double>(index, fraction);
520 }
521 
522 
523 G4int
524 G4RPGInelastic::sampleFlat(std::vector<G4double> sigma) const
525 {
526  G4int i;
527  G4double sum(0.);
528  for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];
529 
530  G4double fsum = sum*G4UniformRand();
531  G4double partialSum = 0.0;
532  G4int channel = 0;
533 
534  for (i = 0; i < G4int(sigma.size()); i++) {
535  partialSum += sigma[i];
536  if (fsum < partialSum) {
537  channel = i;
538  break;
539  }
540  }
541 
542  return channel;
543 }
544 
545 
547  G4int &vecLen,
548  G4ReactionProduct &currentParticle,
549  G4ReactionProduct &targetParticle,
551 {
552  const G4ParticleDefinition* projDef = currentParticle.GetDefinition();
553  const G4ParticleDefinition* targDef = targetParticle.GetDefinition();
554  G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
555  G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
556  G4double strangenessSum = projDef->GetQuarkContent(3) -
557  projDef->GetAntiQuarkContent(3) +
558  targDef->GetQuarkContent(3) -
559  targDef->GetAntiQuarkContent(3);
560 
561  const G4ParticleDefinition* secDef = 0;
562  for (G4int i = 0; i < vecLen; i++) {
563  secDef = vec[i]->GetDefinition();
564  chargeSum += secDef->GetPDGCharge();
565  baryonSum += secDef->GetBaryonNumber();
566  strangenessSum += secDef->GetQuarkContent(3)
567  - secDef->GetAntiQuarkContent(3);
568  }
569 
570  G4bool OK = true;
571  if (chargeSum != Q) {
572  G4cout << " Charge not conserved " << G4endl;
573  OK = false;
574  }
575  if (baryonSum != B) {
576  G4cout << " Baryon number not conserved " << G4endl;
577  OK = false;
578  }
579  if (strangenessSum != S) {
580  G4cout << " Strangeness not conserved " << G4endl;
581  OK = false;
582  }
583 
584  if (!OK) {
585  G4cout << " projectile: " << projDef->GetParticleName()
586  << " target: " << targDef->GetParticleName() << G4endl;
587  for (G4int i = 0; i < vecLen; i++) {
588  secDef = vec[i]->GetDefinition();
589  G4cout << secDef->GetParticleName() << " " ;
590  }
591  G4cout << G4endl;
592  }
593 
594 }
595 
596 
598  0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
599  0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
600  2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
601 
602 
603 /* end of file */
static G4KaonZero * KaonZero()
Definition: G4KaonZero.cc:104
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
static const double MeV
Definition: G4SIunits.hh:211
G4bool MarkLeadingStrangeParticle(const G4ReactionProduct &currentParticle, const G4ReactionProduct &targetParticle, G4ReactionProduct &leadParticle)
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
double S(double temp)
CLHEP::Hep3Vector G4ThreeVector
void SetUpChange(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
static G4OmegaMinus * OmegaMinus()
static G4KaonZeroLong * KaonZeroLong()
void SetSide(const G4int sid)
G4int GetAntiQuarkContent(G4int flavor) const
void CalculateMomenta(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4DynamicParticle *originalTarget, G4ReactionProduct &modifiedOriginal, G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool quasiElastic)
static double Q[]
G4int sampleFlat(std::vector< G4double > sigma) const
G4ParticleDefinition * GetDefinition() const
G4ThreeVector what
double B(double temperature)
void Initialize(G4int items)
Definition: G4FastVector.hh:63
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:102
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
void SetStatusChange(G4HadFinalStateStatus aS)
G4double Pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
static G4XiZero * XiZero()
Definition: G4XiZero.cc:106
const G4ParticleDefinition * GetDefinition() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
G4ParticleDefinition * particleDef[18]
static G4KaonZeroShort * KaonZeroShort()
static const G4double energyScale[30]
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
static G4XiMinus * XiMinus()
Definition: G4XiMinus.cc:106
G4int GetQuarkContent(G4int flavor) const
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static const double GeV
Definition: G4SIunits.hh:214
void SetDefinitionAndUpdateE(const G4ParticleDefinition *aParticleDefinition)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4int n
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
Definition: G4PionZero.cc:108
static G4SigmaMinus * SigmaMinus()
G4double GetKineticEnergy() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const double eV
Definition: G4SIunits.hh:212
void SetEnergyChange(G4double anEnergy)
G4double AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
Definition: G4Nucleus.cc:336
static G4AntiKaonZero * AntiKaonZero()
G4double GetPDGMass() const
G4RPGTwoCluster twoCluster
static const double pi
Definition: G4SIunits.hh:74
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
Definition: G4RPGTwoBody.cc:45
G4RPGInelastic(const G4String &modelName="RPGInelastic")
G4double Cinema(G4double kineticEnergy)
Definition: G4Nucleus.cc:381
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
#define DBL_MIN
Definition: templates.hh:75
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4ThreeVector GetMomentum() const
G4int Factorial(G4int n)
#define G4endl
Definition: G4ios.hh:61
void CheckQnums(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4double Q, G4double B, G4double S)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
void GetNormalizationConstant(const G4double availableEnergy, G4double &n, G4double &anpn)
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4double GetPDGCharge() const
void SetUpPions(const G4int np, const G4int nm, const G4int nz, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
G4RPGFragmentation fragmentation
std::pair< G4int, G4double > interpolateEnergy(G4double ke) const
void SetMomentumChange(const G4ThreeVector &aV)
void Report(std::ostream &aS)
G4double GetMass() const
static G4AntiNeutron * AntiNeutron()
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
G4RPGTwoBody twoBody