ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCrossSectionsINCL46.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCrossSectionsINCL46.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 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLParticleTable.hh"
41 #include "G4INCLLogger.hh"
42 // #include <cassert>
43 
44 namespace G4INCL {
45 
46 /* G4double elasticNNHighEnergy(const G4double momentum) {
47  return 77.0/(momentum + 1.5);
48  }
49 
50  G4double elasticProtonNeutron(const G4double momentum) {
51  if(momentum < 0.450) {
52  const G4double alp = std::log(momentum);
53  return 6.3555*std::exp(-3.2481*alp-0.377*alp*alp);
54  } else if(momentum >= 0.450 && momentum < 0.8) {
55  return (33.0 + 196.0 * std::sqrt(std::pow(std::abs(momentum - 0.95), 5)));
56  } else if(momentum > 2.0) {
57  return elasticNNHighEnergy(momentum);
58  } else {
59  return 31.0/std::sqrt(momentum);
60  }
61  }
62 
63  G4double elasticProtonProtonOrNeutronNeutron(const G4double momentum)
64  {
65  if(momentum < 0.440) {
66  return 34.0*std::pow(momentum/0.4, -2.104);
67  } else if(momentum < 0.8 && momentum >= 0.440) {
68  return (23.5 + 1000.0*std::pow(momentum-0.7, 4));
69  } else if(momentum < 2.0) {
70  return (1250.0/(50.0 + momentum) - 4.0*std::pow(momentum-1.3, 2));
71  } else {
72  return elasticNNHighEnergy(momentum);
73  }
74  }
75 
76  G4double elasticNN(Particle const * const p1, Particle const * const p2) {
77  G4double momentum = 0.0;
78  momentum = 0.001 * KinematicsUtils::momentumInLab(p1, p2);
79  if((p1->getType() == Proton && p2->getType() == Proton) ||
80  (p1->getType() == Neutron && p2->getType() == Neutron)) {
81  return elasticProtonProtonOrNeutronNeutron(momentum);
82  } else if((p1->getType() == Proton && p2->getType() == Neutron) ||
83  (p1->getType() == Neutron && p2->getType() == Proton)) {
84  return elasticProtonNeutron(momentum);
85  } else {
86  INCL_ERROR("CrossSectionsINCL46::elasticNN: Bad input!" << '\n'
87  << p1->print() << p2->print() << '\n');
88  }
89  return 0.0;
90  }:*/
91 
92  G4double CrossSectionsINCL46::elasticNNLegacy(Particle const * const part1, Particle const * const part2) {
93 
94 
97 
98  /* The NN cross section is parametrised as a function of the lab momentum
99  * of one of the nucleons. For NDelta or DeltaDelta, the physical
100  * assumption is that the cross section is the same as NN *for the same
101  * total CM energy*. Thus, we calculate s from the particles involved, and
102  * we convert this value to the lab momentum of a nucleon *as if this were
103  * an NN collision*.
104  */
105 
106  const G4double s = KinematicsUtils::squareTotalEnergyInCM(part1, part2);
108  if(plab > 2.) { // NN, Delta-Nucleon and Delta-Delta for plab > 2.0 GeV
109  return 77./(plab+1.5);
110  }
111  else if (part1->isNucleon() && part2->isNucleon()){ // NN
112  if (i == 0) { // pn
113  if (plab < 0.450) {
114  G4double alp=std::log(plab);
115  return 6.3555*std::exp(-3.2481*alp-0.377*alp*alp);
116  }
117  else if (plab < 0.800) {
118  return (33.+196.*std::sqrt(std::pow(std::abs(plab-0.95),5)));
119  }
120  else {
121  return 31./std::sqrt(plab);
122  }
123  }
124  else { // nn and pp
125  if (plab < 0.440) {
126  return 34.*std::pow(plab/0.4, (-2.104));
127  }
128  else if (plab < 0.800) {
129  return (23.5+1000.*std::pow(plab-0.7, 4));
130  }
131  else {
132  return (1250./(50.+plab)-4.*std::pow(plab-1.3, 2));
133  }
134  }
135  }
136  else { // Delta-Nucleon and Delta-Delta
137  if (plab < 0.440) {
138  return 34.*std::pow(plab/0.4, (-2.104));
139  }
140  else if (plab < 0.800) {
141  return (23.5+1000.*std::pow(plab-0.7, 4));
142  }
143  else {
144  return (1250./(50.+plab)-4.*std::pow(plab-1.3, 2));
145  }
146  }
147  }
148 
150  G4double xs = 0.0;
151 // assert(isospin==-2 || isospin==0 || isospin==2);
152 
153  const G4double momentumGeV = 0.001 * pLab;
154  if(pLab < 800.0) {
155  return 0.0;
156  }
157 
158  if(isospin==2 || isospin==-2) { // pp, nn
159  if(pLab >= 2000.0) {
160  xs = (41.0 + (60.0*momentumGeV - 54.0)*std::exp(-1.2*momentumGeV) - 77.0/(momentumGeV + 1.5));
161  } else if(pLab >= 1500.0 && pLab < 2000.0) {
162  xs = (41.0 + 60.0*(momentumGeV - 0.9)*std::exp(-1.2*momentumGeV) - 1250.0/(momentumGeV+50.0)+ 4.0*std::pow(momentumGeV - 1.3, 2));
163  } else if(pLab < 1500.0) {
164  xs = (23.5 + 24.6/(1.0 + std::exp(-10.0*momentumGeV + 12.0))
165  -1250.0/(momentumGeV +50.0)+4.0*std::pow(momentumGeV - 1.3,2));
166  }
167  } else if(isospin==0) { // pn
168  if(pLab >= 2000.0) {
169  xs = (42.0 - 77.0/(momentumGeV + 1.5));
170  } else if(pLab >= 1000.0 && pLab < 2000.0) {
171  xs = (24.2 + 8.9*momentumGeV - 31.1/std::sqrt(momentumGeV));
172  } else if(pLab < 1000.0) {
173  xs = (33.0 + 196.0*std::sqrt(std::pow(std::abs(momentumGeV - 0.95),5))
174  -31.1/std::sqrt(momentumGeV));
175  }
176  }
177 
178  if(xs < 0.0) return 0.0;
179  else return xs;
180  }
181 
183  // HE and LE pi- p and pi+ n
184  if(x <= 1750.0) {
185  return -2.33730e-06*std::pow(x, 3)+1.13819e-02*std::pow(x,2)
186  -1.83993e+01*x+9893.4;
187  } else if(x > 1750.0 && x <= 2175.0) {
188  return 1.13531e-06*std::pow(x, 3)-6.91694e-03*std::pow(x, 2)
189  +1.39907e+01*x-9360.76;
190  } else {
191  return -3.18087*std::log(x)+52.9784;
192  }
193  }
194 
196  // HE pi- p and pi+ n
197  if(x <= 1475.0) {
198  return 0.00120683*(x-1372.52)*(x-1372.52)+26.2058;
199  } else if(x > 1475.0 && x <= 1565.0) {
200  return 1.15873e-05*x*x+49965.6/((x-1519.59)*(x-1519.59)+2372.55);
201  } else if(x > 1565.0 && x <= 2400.0) {
202  return 34.0248+43262.2/((x-1681.65)*(x-1681.65)+1689.35);
203  } else if(x > 2400.0 && x <= 7500.0) {
204  return 3.3e-7*(x-7500.0)*(x-7500.0)+24.5;
205  } else {
206  return 24.5;
207  }
208  }
209 
210  G4double CrossSectionsINCL46::total(Particle const * const p1, Particle const * const p2) {
211  G4double inelastic = 0.0;
212  if(p1->isNucleon() && p2->isNucleon()) {
213  inelastic = NNToNDelta(p1, p2);
214  } else if((p1->isNucleon() && p2->isDelta()) ||
215  (p1->isDelta() && p2->isNucleon())) {
216  inelastic = NDeltaToNN(p1, p2);
217  } else if((p1->isNucleon() && p2->isPion()) ||
218  (p1->isPion() && p2->isNucleon())) {
219  inelastic = piNToDelta(p1, p2);
220  } else {
221  inelastic = 0.0;
222  }
223 
224  return inelastic + elastic(p1, p2);
225  }
226 
227  G4double CrossSectionsINCL46::piNToDelta(Particle const * const particle1, Particle const * const particle2) {
228  // FUNCTION SPN(X,IND2T3,IPIT3,f17)
229  // SIGMA(PI+ + P) IN THE (3,3) REGION
230  // NEW FIT BY J.VANDERMEULEN + FIT BY Th AOUST ABOVE (3,3) RES
231  // CONST AT LOW AND VERY HIGH ENERGY
232  // COMMON/BL8/RATHR,RAMASS REL21800
233  // integer f17
234  // RATHR and RAMASS are always 0.0!!!
235 
236  G4double x = KinematicsUtils::totalEnergyInCM(particle1, particle2);
237  if(x>10000.) return 0.0; // no cross section above this value
238 
239  G4int ipit3 = 0;
240  G4int ind2t3 = 0;
241  G4double ramass = 0.0;
242 
243  if(particle1->isPion()) {
244  ipit3 = ParticleTable::getIsospin(particle1->getType());
245  } else if(particle2->isPion()) {
246  ipit3 = ParticleTable::getIsospin(particle2->getType());
247  }
248 
249  if(particle1->isNucleon()) {
250  ind2t3 = ParticleTable::getIsospin(particle1->getType());
251  } else if(particle2->isNucleon()) {
252  ind2t3 = ParticleTable::getIsospin(particle2->getType());
253  }
254 
255  G4double y=x*x;
256  G4double q2=(y-1076.0*1076.0)*(y-800.0*800.0)/y/4.0;
257  if (q2 <= 0.) {
258  return 0.0;
259  }
260  G4double q3 = std::pow(std::sqrt(q2),3);
261  G4double f3 = q3/(q3 + 5832000.); // 5832000 = 180^3
262  G4double spnResult = 326.5/(std::pow((x-1215.0-ramass)*2.0/(110.0-ramass), 2)+1.0);
263  spnResult = spnResult*(1.0-5.0*ramass/1215.0);
264  G4double cg = 4.0 + G4double(ind2t3)*G4double(ipit3);
265  spnResult = spnResult*f3*cg/6.0;
266 
267  if(x < 1200.0 && spnResult < 5.0) {
268  spnResult = 5.0;
269  }
270 
271  // HE pi+ p and pi- n
272  if(x > 1290.0) {
273  if((ind2t3 == 1 && ipit3 == 2) || (ind2t3 == -1 && ipit3 == -2))
274  spnResult=spnPiPlusPHE(x);
275  else if((ind2t3 == 1 && ipit3 == -2) || (ind2t3 == -1 && ipit3 == 2))
276  spnResult=spnPiMinusPHE(x);
277  else if(ipit3 == 0) spnResult = (spnPiPlusPHE(x) + spnPiMinusPHE(x))/2.0; // (spnpipphe(x)+spnpimphe(x))/2.0
278  else {
279  INCL_ERROR("Unknown configuration!" << '\n');
280  }
281  }
282 
283  return spnResult;
284  }
285 
286  G4double CrossSectionsINCL46::NDeltaToNN(Particle const * const p1, Particle const * const p2) {
288  if(isospin==4 || isospin==-4) return 0.0;
289 
291  G4double Ecm = std::sqrt(s);
292  G4int deltaIsospin;
293  G4double deltaMass;
294  if(p1->isDelta()) {
295  deltaIsospin = ParticleTable::getIsospin(p1->getType());
296  deltaMass = p1->getMass();
297  } else {
298  deltaIsospin = ParticleTable::getIsospin(p2->getType());
299  deltaMass = p2->getMass();
300  }
301 
302  if(Ecm <= 938.3 + deltaMass) {
303  return 0.0;
304  }
305 
306  if(Ecm < 938.3 + deltaMass + 2.0) {
307  Ecm = 938.3 + deltaMass + 2.0;
308  s = Ecm*Ecm;
309  }
310 
312  (s - std::pow(ParticleTable::effectiveNucleonMass + deltaMass, 2));
313  const G4double y = s/(s - std::pow(deltaMass - ParticleTable::effectiveNucleonMass, 2));
314  /* Concerning the way we calculate the lab momentum, see the considerations
315  * in CrossSections::elasticNNLegacy().
316  */
318  G4double result = 0.5 * x * y * deltaProduction(isospin, pLab);
319  result *= 3.*(32.0 + isospin * isospin * (deltaIsospin * deltaIsospin - 5))/64.0;
320  result /= 1.0 + 0.25 * isospin * isospin;
321  return result;
322  }
323 
324  G4double CrossSectionsINCL46::NNToNDelta(Particle const * const p1, Particle const * const p2) {
325 // assert(p1->isNucleon() && p2->isNucleon());
326  const G4double sqrts = KinematicsUtils::totalEnergyInCM(p1,p2);
327  if(sqrts < ParticleTable::effectivePionMass + 2*ParticleTable::effectiveNucleonMass + 50.) { // approximately yields INCL4.6's hard-coded threshold in collis, 2065 MeV
328  return 0.0;
329  } else {
330  const G4double pLab = KinematicsUtils::momentumInLab(p1,p2);
332  return deltaProduction(isospin, pLab);
333  }
334  }
335 
336  G4double CrossSectionsINCL46::elastic(Particle const * const p1, Particle const * const p2) {
337 // if(!p1->isPion() && !p2->isPion())
338  if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta()))
339  // return elasticNN(p1, p2); // New implementation
340  return elasticNNLegacy(p1, p2); // Translated from INCL4.6 FORTRAN
341  else
342  return 0.0; // No pion-nucleon elastic scattering
343  }
344 
346  G4double x = 0.001 * pl; // Change to GeV
347  if(iso != 0) {
348  if(pl <= 2000.0) {
349  x = std::pow(x, 8);
350  return 5.5e-6 * x/(7.7 + x);
351  } else {
352  return (5.34 + 0.67*(x - 2.0)) * 1.0e-6;
353  }
354  } else {
355  if(pl < 800.0) {
356  G4double b = (7.16 - 1.63*x) * 1.0e-6;
357  return b/(1.0 + std::exp(-(x - 0.45)/0.05));
358  } else if(pl < 1100.0) {
359  return (9.87 - 4.88 * x) * 1.0e-6;
360  } else {
361  return (3.68 + 0.76*x) * 1.0e-6;
362  }
363  }
364  return 0.0; // Should never reach this point
365  }
366 
367 
368  G4double CrossSectionsINCL46::NNToxPiNN(const G4int, Particle const * const, Particle const * const) {
369  return 0.;
370  }
371 
372  G4double CrossSectionsINCL46::piNToxPiN(const G4int, Particle const * const, Particle const * const) {
373  return 0.;
374  }
375 
377  //
378  // Pion-Nucleon producing Eta cross sections
379  //
380  return 0.;
381  }
382 
384  //
385  // Pion-Nucleon producing Omega cross sections
386  //
387  return 0.;
388  }
389 
391  //
392  // Pion-Nucleon producing EtaPrime cross sections
393  //
394  return 0.;
395  }
396 
398  //
399  // Eta-Nucleon producing Pion cross sections
400  //
401  return 0.;
402  }
403 
404 
406  //
407  // Eta-Nucleon producing Two Pions cross sections
408  //
409  return 0.;
410  }
411 
413  //
414  // Omega-Nucleon producing Pion cross sections
415  //
416  return 0.;
417  }
418 
420  //
421  // Omega-Nucleon producing Two Pions cross sections
422  //
423  return 0.;
424  }
425 
427  //
428  // EtaPrime-Nucleon producing Pion cross sections
429  //
430  return 0.;
431  }
432 
434  //
435  // Nucleon-Nucleon producing Eta cross sections
436  //
437  return 0.;
438  }
439 
441  //
442  // Nucleon-Nucleon producing Eta cross sections
443  //
444  return 0.;
445  }
446 
447  G4double CrossSectionsINCL46::NNToNNEtaxPi(const G4int, Particle const * const, Particle const * const) {
448  return 0.;
449  }
450 
452  //
453  // Nucleon-Nucleon producing N-Delta-Eta cross sections
454  //
455  return 0.;
456  }
457 
459  //
460  // Nucleon-Nucleon producing Omega cross sections
461  //
462  return 0.;
463  }
464 
466  //
467  // Nucleon-Nucleon producing Omega cross sections
468  //
469  return 0.;
470  }
471 
472  G4double CrossSectionsINCL46::NNToNNOmegaxPi(const G4int, Particle const * const, Particle const * const) {
473  return 0.;
474  }
475 
477  //
478  // Nucleon-Nucleon producing N-Delta-Omega cross sections
479  //
480  return 0.;
481  }
482 
483 
484 
485  G4double CrossSectionsINCL46::NYelastic(Particle const * const , Particle const * const ) {
486  //
487  // Hyperon-Nucleon elastic cross sections
488  //
489  return 0.;
490  }
491 
492  G4double CrossSectionsINCL46::NKelastic(Particle const * const , Particle const * const ) {
493  //
494  // Kaon-Nucleon elastic cross sections
495  //
496  return 0.;
497  }
498 
499  G4double CrossSectionsINCL46::NKbelastic(Particle const * const , Particle const * const ) {
500  //
501  // antiKaon-Nucleon elastic cross sections
502  //
503  return 0.;
504  }
505 
506  G4double CrossSectionsINCL46::NNToNLK(Particle const * const, Particle const * const) {
507  //
508  // Nucleon-Nucleon producing N-Lambda-Kaon cross sections
509  //
510  return 0.;
511  }
512 
513  G4double CrossSectionsINCL46::NNToNSK(Particle const * const, Particle const * const) {
514  //
515  // Nucleon-Nucleon producing N-Sigma-Kaon cross sections
516  //
517  return 0.;
518  }
519 
521  //
522  // Nucleon-Nucleon producing N-Lambda-Kaon-pion cross sections
523  //
524  return 0.;
525  }
526 
528  //
529  // Nucleon-Nucleon producing N-Sigma-Kaon-pion cross sections
530  //
531  return 0.;
532  }
533 
535  //
536  // Nucleon-Nucleon producing N-Lambda-Kaon-2pion cross sections
537  //
538  return 0.;
539  }
540 
542  //
543  // Nucleon-Nucleon producing N-Sigma-Kaon-2pion cross sections
544  //
545  return 0.;
546  }
547 
549  //
550  // Nucleon-Nucleon producing Nucleon-Nucleon-Kaon-antiKaon cross sections
551  //
552  return 0.;
553  }
554 
556  //
557  // Nucleon-Nucleon missing strangeness production cross sections
558  //
559  return 0.;
560  }
561 
563  // Nucleon-Delta producing Nucleon Lambda Kaon cross section
564  return 0;
565  }
567  // Nucleon-Delta producing Nucleon Sigma Kaon cross section
568  return 0;
569  }
571  // Nucleon-Delta producing Delta Lambda Kaon cross section
572  return 0;
573  }
575  // Nucleon-Delta producing Delta Sigma Kaon cross section
576  return 0;
577  }
578 
580  // Nucleon-Delta producing Nucleon-Nucleon Kaon antiKaon cross section
581  return 0;
582  }
583 
584  G4double CrossSectionsINCL46::NpiToLK(Particle const * const, Particle const * const) {
585  //
586  // Pion-Nucleon producing Lambda-Kaon cross sections
587  //
588  return 0.;
589  }
590 
591  G4double CrossSectionsINCL46::NpiToSK(Particle const * const, Particle const * const) {
592  //
593  // Pion-Nucleon producing Sigma-Kaon cross sections
594  //
595  return 0.;
596  }
598  return 0.;
599  }
601  return 0.;
602  }
604  return 0.;
605  }
606 
608  //
609  // Pion-Nucleon producing Lambda-Kaon-pion cross sections
610  //
611  return 0.;
612  }
613 
615  //
616  // Pion-Nucleon producing Sigma-Kaon-pion cross sections
617  //
618  return 0.;
619  }
620 
622  //
623  // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
624  //
625  return 0.;
626  }
627 
629  //
630  // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
631  //
632  return 0.;
633  }
634 
636  //
637  // Pion-Nucleon producing Nucleon-Kaon-antiKaon cross sections
638  //
639  return 0.;
640  }
641 
643  //
644  // Pion-Nucleon missing strangeness production cross sections
645  //
646  return 0.;
647  }
648 
649  G4double CrossSectionsINCL46::NLToNS(Particle const * const, Particle const * const) {
650  //
651  // Nucleon-Hyperon multiplet changing cross sections
652  //
653  return 0.;
654  }
655 
656  G4double CrossSectionsINCL46::NSToNL(Particle const * const, Particle const * const) {
657  //
658  // Nucleon-Sigma quasi-elastic cross sections
659  //
660  return 0.;
661  }
662 
663  G4double CrossSectionsINCL46::NSToNS(Particle const * const, Particle const * const) {
664  //
665  // Nucleon-Sigma quasi-elastic cross sections
666  //
667  return 0.;
668  }
669 
670  G4double CrossSectionsINCL46::NKToNK(Particle const * const, Particle const * const) {
671  //
672  // Nucleon-Kaon quasi-elastic cross sections
673  //
674  return 0.;
675  }
676 
677  G4double CrossSectionsINCL46::NKToNKpi(Particle const * const, Particle const * const) {
678  //
679  // Nucleon-Kaon producing Nucleon-Kaon-pion cross sections
680  //
681  return 0.;
682  }
683 
685  //
686  // Nucleon-Kaon producing Nucleon-Kaon-2pion cross sections
687  //
688  return 0.;
689  }
690 
691  G4double CrossSectionsINCL46::NKbToNKb(Particle const * const, Particle const * const) {
692  //
693  // Nucleon-antiKaon quasi-elastic cross sections
694  //
695  return 0.;
696  }
697 
698  G4double CrossSectionsINCL46::NKbToSpi(Particle const * const, Particle const * const) {
699  //
700  // Nucleon-antiKaon producing Sigma-pion cross sections
701  //
702  return 0.;
703  }
704 
705  G4double CrossSectionsINCL46::NKbToLpi(Particle const * const, Particle const * const) {
706  //
707  // Nucleon-antiKaon producing Lambda-pion cross sections
708  //
709  return 0.;
710  }
711 
713  //
714  // Nucleon-antiKaon producing Sigma-2pion cross sections
715  //
716  return 0.;
717  }
718 
720  //
721  // Nucleon-antiKaon producing Lambda-2pion cross sections
722  //
723  return 0.;
724  }
725 
727  //
728  // Nucleon-antiKaon producing Nucleon-antiKaon-pion cross sections
729  //
730  return 0.;
731  }
732 
734  //
735  // Nucleon-antiKaon producing Nucleon-antiKaon-2pion cross sections
736  //
737  return 0.;
738  }
739 
740 } // namespace G4INCL
741