ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGInelastic.cc
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 //
27 
28 #include "G4RPGInelastic.hh"
29 #include "G4Log.hh"
30 #include "Randomize.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
35 #include "G4RPGTwoBody.hh"
36 
37 
39  : G4HadronicInteraction(modelName)
40 {
41  cache = 0.0;
60 
61  G4cout << " **************************************************** " << G4endl;
62  G4cout << " * The RPG model is currently under development and * " << G4endl;
63  G4cout << " * should not be used. * " << G4endl;
64  G4cout << " **************************************************** " << G4endl;
65 }
66 
67 
70 {
71  const G4double expxu = 82.; // upper bound for arg. of exp
72  const G4double expxl = -expxu; // lower bound for arg. of exp
73  G4double npf = 0.0;
74  G4double nmf = 0.0;
75  G4double nzf = 0.0;
76  G4int i;
77  for( i=2; i<=np; i++ )npf += G4Log((double)i);
78  for( i=2; i<=nneg; i++ )nmf += G4Log((double)i);
79  for( i=2; i<=nz; i++ )nzf += G4Log((double)i);
80  G4double r;
81  r = std::min( expxu, std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
82  return G4Exp(r);
83 }
84 
85 
87 {
88  G4int j = std::min(n,10);
89  G4int result = 1;
90  if (j <= 1) return result;
91  for (G4int i = 2; i <= j; ++i) result *= i;
92  return result;
93 }
94 
95 
97  const G4ReactionProduct &currentParticle,
98  const G4ReactionProduct &targetParticle,
99  G4ReactionProduct &leadParticle )
100 {
101  // The following was in GenerateXandPt and TwoCluster.
102  // Add a parameter to the GenerateXandPt function telling it about the
103  // strange particle.
104  //
105  // Assumes that the original particle was a strange particle
106  //
107  G4bool lead = false;
108  if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
109  (currentParticle.GetDefinition() != G4Proton::Proton()) &&
110  (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
111  {
112  lead = true;
113  leadParticle = currentParticle; // set lead to the incident particle
114  }
115  else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
116  (targetParticle.GetDefinition() != G4Proton::Proton()) &&
117  (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
118  {
119  lead = true;
120  leadParticle = targetParticle; // set lead to the target particle
121  }
122  return lead;
123 }
124 
125 
126  void G4RPGInelastic::SetUpPions(const G4int np, const G4int nneg,
127  const G4int nz,
129  G4int &vecLen)
130  {
131  if( np+nneg+nz == 0 )return;
132  G4int i;
134  for( i=0; i<np; ++i )
135  {
136  p = new G4ReactionProduct;
138  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
139  vec.SetElement( vecLen++, p );
140  }
141  for( i=np; i<np+nneg; ++i )
142  {
143  p = new G4ReactionProduct;
145  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
146  vec.SetElement( vecLen++, p );
147  }
148  for( i=np+nneg; i<np+nneg+nz; ++i )
149  {
150  p = new G4ReactionProduct;
152  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
153  vec.SetElement( vecLen++, p );
154  }
155  }
156 
157 
159  const G4double energy, // MeV, <0 means annihilation channels
160  G4double &n,
161  G4double &anpn )
162  {
163  const G4double expxu = 82.; // upper bound for arg. of exp
164  const G4double expxl = -expxu; // lower bound for arg. of exp
165  const G4int numSec = 60;
166  //
167  // the only difference between the calculation for annihilation channels
168  // and normal is the starting value, iBegin, for the loop below
169  //
170  G4int iBegin = 1;
171  G4double en = energy;
172  if( energy < 0.0 )
173  {
174  iBegin = 2;
175  en *= -1.0;
176  }
177  //
178  // number of total particles vs. centre of mass Energy - 2*proton mass
179  //
180  G4double aleab = G4Log(en/GeV);
181  n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
182  n -= 2.0;
183  //
184  // normalization constant for kno-distribution
185  //
186  anpn = 0.0;
187  G4double test, temp;
188  for( G4int i=iBegin; i<=numSec; ++i )
189  {
190  temp = pi*i/(2.0*n*n);
191  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
192  if( temp < 1.0 )
193  {
194  if( test >= 1.0e-10 )anpn += temp*test;
195  }
196  else
197  anpn += temp*test;
198  }
199  }
200 
201 void
203  G4int& vecLen,
204  const G4HadProjectile* originalIncident,
205  const G4DynamicParticle* originalTarget,
206  G4ReactionProduct& modifiedOriginal,
207  G4Nucleus& targetNucleus,
208  G4ReactionProduct& currentParticle,
209  G4ReactionProduct& targetParticle,
210  G4bool& incidentHasChanged,
211  G4bool& targetHasChanged,
212  G4bool quasiElastic)
213 {
214  cache = 0;
215  what = originalIncident->Get4Momentum().vect();
216 
217  G4ReactionProduct leadingStrangeParticle;
218 
219  // strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
220  // incidentHasChanged, originalTarget,
221  // targetParticle, targetHasChanged,
222  // targetNucleus, currentParticle,
223  // vec, vecLen,
224  // false, leadingStrangeParticle);
225 
226  if( quasiElastic )
227  {
228  twoBody.ReactionStage(originalIncident, modifiedOriginal,
229  incidentHasChanged, originalTarget,
230  targetParticle, targetHasChanged,
231  targetNucleus, currentParticle,
232  vec, vecLen,
233  false, leadingStrangeParticle);
234  return;
235  }
236 
237  G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
238  targetParticle,
239  leadingStrangeParticle );
240  //
241  // Note: the number of secondaries can be reduced in GenerateXandPt
242  // and TwoCluster
243  //
244  G4bool finishedGenXPt = false;
245  G4bool annihilation = false;
246  if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
247  currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
248  {
249  // original was an anti-particle and annihilation has taken place
250  annihilation = true;
251  G4double ekcor = 1.0;
252  G4double ek = originalIncident->GetKineticEnergy();
253  G4double ekOrg = ek;
254 
255  const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
256  if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
257  const G4double atomicWeight = targetNucleus.GetA_asInt();
258  ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
259  G4double tkin = targetNucleus.Cinema(ek);
260  ek += tkin;
261  ekOrg += tkin;
262  // modifiedOriginal.SetKineticEnergy( ekOrg );
263  //
264  // evaporation -- re-calculate black track energies
265  // this was Done already just before the cascade
266  //
267  tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
268  ekOrg -= tkin;
269  ekOrg = std::max( 0.0001*GeV, ekOrg );
270  modifiedOriginal.SetKineticEnergy( ekOrg );
271  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
272  G4double et = ekOrg + amas;
273  G4double p = std::sqrt( std::abs(et*et-amas*amas) );
274  G4double pp = modifiedOriginal.GetMomentum().mag();
275  if( pp > 0.0 )
276  {
277  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
278  modifiedOriginal.SetMomentum( momentum * (p/pp) );
279  }
280  if( ekOrg <= 0.0001 )
281  {
282  modifiedOriginal.SetKineticEnergy( 0.0 );
283  modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
284  }
285  }
286 
287  // twsup gives percentage of time two-cluster model is called
288 
289  const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
292 
293  // Cache current, target, and secondaries
294  G4ReactionProduct saveCurrent = currentParticle;
295  G4ReactionProduct saveTarget = targetParticle;
296  std::vector<G4ReactionProduct> savevec;
297  for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
298 
299  // Call fragmentation code if
300  // 1) there is annihilation, or
301  // 2) there are more than 5 secondaries, or
302  // 3) incident KE is > 1 GeV AND
303  // ( incident is a kaon AND rand < 0.5 OR twsup )
304  //
305 
306  if( annihilation || vecLen > 5 ||
307  ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
308 
309  (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
310  originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
311  originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
312  originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
313  rand1 < 0.5)
314  || rand2 > twsup[vecLen]) ) )
315 
316  finishedGenXPt =
317  fragmentation.ReactionStage(originalIncident, modifiedOriginal,
318  incidentHasChanged, originalTarget,
319  targetParticle, targetHasChanged,
320  targetNucleus, currentParticle,
321  vec, vecLen,
322  leadFlag, leadingStrangeParticle);
323 
324  if (finishedGenXPt) return;
325 
326  G4bool finishedTwoClu = false;
327 
328  if (modifiedOriginal.GetTotalMomentum() < 1.0) {
329  for (G4int i = 0; i < vecLen; i++) delete vec[i];
330  vecLen = 0;
331 
332  } else {
333  // Occaisionally, GenerateXandPt will fail in the annihilation channel.
334  // Restore current, target and secondaries to pre-GenerateXandPt state
335  // before trying annihilation in TwoCluster
336 
337  if (!finishedGenXPt && annihilation) {
338  currentParticle = saveCurrent;
339  targetParticle = saveTarget;
340  for (G4int i = 0; i < vecLen; i++) delete vec[i];
341  vecLen = 0;
342  vec.Initialize( 0 );
343  for (G4int i = 0; i < G4int(savevec.size()); i++) {
345  *p = savevec[i];
346  vec.SetElement( vecLen++, p );
347  }
348  }
349 
350  // Big violations of energy conservation in this method - don't use
351  //
352  // pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
353  // incidentHasChanged, originalTarget,
354  // targetParticle, targetHasChanged,
355  // targetNucleus, currentParticle,
356  // vec, vecLen,
357  // false, leadingStrangeParticle);
358 
359  try
360  {
361  finishedTwoClu =
362  twoCluster.ReactionStage(originalIncident, modifiedOriginal,
363  incidentHasChanged, originalTarget,
364  targetParticle, targetHasChanged,
365  targetNucleus, currentParticle,
366  vec, vecLen,
367  leadFlag, leadingStrangeParticle);
368  }
369  catch(G4HadReentrentException & aC)
370  {
371  aC.Report(G4cout);
372  throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
373  }
374  }
375 
376  if (finishedTwoClu) return;
377 
378  twoBody.ReactionStage(originalIncident, modifiedOriginal,
379  incidentHasChanged, originalTarget,
380  targetParticle, targetHasChanged,
381  targetNucleus, currentParticle,
382  vec, vecLen,
383  false, leadingStrangeParticle);
384 }
385 
386 /*
387  void G4RPGInelastic::
388  Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen)
389  {
390  G4double rotation = 2.*pi*G4UniformRand();
391  cache = rotation;
392  G4int i;
393  for( i=0; i<vecLen; ++i )
394  {
395  G4ThreeVector momentum = vec[i]->GetMomentum();
396  momentum = momentum.rotate(rotation, what);
397  vec[i]->SetMomentum(momentum);
398  }
399  }
400 */
401 
402 void
404  G4int& vecLen,
405  G4ReactionProduct& currentParticle,
406  G4ReactionProduct& targetParticle,
407  G4bool& incidentHasChanged )
408 {
412  G4int i;
413 
414  if (currentParticle.GetDefinition() == particleDef[k0]) {
415  if (G4UniformRand() < 0.5) {
416  currentParticle.SetDefinitionAndUpdateE(aKaonZL);
417  incidentHasChanged = true;
418  } else {
419  currentParticle.SetDefinitionAndUpdateE(aKaonZS);
420  }
421  } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
422  if (G4UniformRand() < 0.5) {
423  currentParticle.SetDefinitionAndUpdateE(aKaonZL);
424  } else {
425  currentParticle.SetDefinitionAndUpdateE(aKaonZS);
426  incidentHasChanged = true;
427  }
428  }
429 
430  if (targetParticle.GetDefinition() == particleDef[k0] ||
431  targetParticle.GetDefinition() == particleDef[k0b] ) {
432  if (G4UniformRand() < 0.5) {
433  targetParticle.SetDefinitionAndUpdateE(aKaonZL);
434  } else {
435  targetParticle.SetDefinitionAndUpdateE(aKaonZS);
436  }
437  }
438 
439  for (i = 0; i < vecLen; ++i) {
440  if (vec[i]->GetDefinition() == particleDef[k0] ||
441  vec[i]->GetDefinition() == particleDef[k0b] ) {
442  if (G4UniformRand() < 0.5) {
443  vec[i]->SetDefinitionAndUpdateE(aKaonZL);
444  } else {
445  vec[i]->SetDefinitionAndUpdateE(aKaonZS);
446  }
447  }
448  }
449 
450  if (incidentHasChanged) {
452  p0->SetDefinition(currentParticle.GetDefinition() );
453  p0->SetMomentum(currentParticle.GetMomentum() );
457 
458  } else {
459  G4double p = currentParticle.GetMomentum().mag()/MeV;
460  G4ThreeVector mom = currentParticle.GetMomentum();
461  if (p > DBL_MIN)
462  theParticleChange.SetMomentumChange(mom.x()/p, mom.y()/p, mom.z()/p );
463  else
464  theParticleChange.SetMomentumChange(0.0, 0.0, 1.0);
465 
466  G4double aE = currentParticle.GetKineticEnergy();
467  if (std::fabs(aE)<.1*eV) aE=.1*eV;
469  }
470 
471  if (targetParticle.GetMass() > 0.0) // Tgt particle can be eliminated in TwoBody
472  {
473  G4ThreeVector momentum = targetParticle.GetMomentum();
474  momentum = momentum.rotate(cache, what);
475  G4double targKE = targetParticle.GetKineticEnergy();
476  G4ThreeVector dir(0.0, 0.0, 1.0);
477  if (targKE < DBL_MIN)
478  targKE = DBL_MIN;
479  else
480  dir = momentum/momentum.mag();
481 
482  G4DynamicParticle* p1 =
483  new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
484 
486  }
487 
489  for (i = 0; i < vecLen; ++i) {
490  G4double secKE = vec[i]->GetKineticEnergy();
491  G4ThreeVector momentum = vec[i]->GetMomentum();
492  G4ThreeVector dir(0.0, 0.0, 1.0);
493  if (secKE < DBL_MIN)
494  secKE = DBL_MIN;
495  else
496  dir = momentum/momentum.mag();
497 
498  p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
500  delete vec[i];
501  }
502 }
503 
504 
505 std::pair<G4int, G4double>
507 {
508  G4int index = 29;
509  G4double fraction = 0.0;
510 
511  for (G4int i = 1; i < 30; i++) {
512  if (e < energyScale[i]) {
513  index = i-1;
514  fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
515  break;
516  }
517  }
518  return std::pair<G4int, G4double>(index, fraction);
519 }
520 
521 
522 G4int
523 G4RPGInelastic::sampleFlat(std::vector<G4double> sigma) const
524 {
525  G4int i;
526  G4double sum(0.);
527  for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];
528 
529  G4double fsum = sum*G4UniformRand();
530  G4double partialSum = 0.0;
531  G4int channel = 0;
532 
533  for (i = 0; i < G4int(sigma.size()); i++) {
534  partialSum += sigma[i];
535  if (fsum < partialSum) {
536  channel = i;
537  break;
538  }
539  }
540 
541  return channel;
542 }
543 
544 
546  G4int &vecLen,
547  G4ReactionProduct &currentParticle,
548  G4ReactionProduct &targetParticle,
550 {
551  const G4ParticleDefinition* projDef = currentParticle.GetDefinition();
552  const G4ParticleDefinition* targDef = targetParticle.GetDefinition();
553  G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
554  G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
555  G4double strangenessSum = projDef->GetQuarkContent(3) -
556  projDef->GetAntiQuarkContent(3) +
557  targDef->GetQuarkContent(3) -
558  targDef->GetAntiQuarkContent(3);
559 
560  const G4ParticleDefinition* secDef = 0;
561  for (G4int i = 0; i < vecLen; i++) {
562  secDef = vec[i]->GetDefinition();
563  chargeSum += secDef->GetPDGCharge();
564  baryonSum += secDef->GetBaryonNumber();
565  strangenessSum += secDef->GetQuarkContent(3)
566  - secDef->GetAntiQuarkContent(3);
567  }
568 
569  G4bool OK = true;
570  if (chargeSum != Q) {
571  G4cout << " Charge not conserved " << G4endl;
572  OK = false;
573  }
574  if (baryonSum != B) {
575  G4cout << " Baryon number not conserved " << G4endl;
576  OK = false;
577  }
578  if (strangenessSum != S) {
579  G4cout << " Strangeness not conserved " << G4endl;
580  OK = false;
581  }
582 
583  if (!OK) {
584  G4cout << " projectile: " << projDef->GetParticleName()
585  << " target: " << targDef->GetParticleName() << G4endl;
586  for (G4int i = 0; i < vecLen; i++) {
587  secDef = vec[i]->GetDefinition();
588  G4cout << secDef->GetParticleName() << " " ;
589  }
590  G4cout << G4endl;
591  }
592 
593 }
594 
595 
597  0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
598  0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
599  2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
600 
601 
602 /* end of file */