ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCrossSections.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCrossSections.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 
38 #include "G4INCLCrossSections.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLParticleTable.hh"
41 #include "G4INCLLogger.hh"
42 #include "G4INCLCrossSectionsINCL46.hh"
43 #include "G4INCLCrossSectionsMultiPions.hh"
44 #include "G4INCLCrossSectionsTruncatedMultiPions.hh"
45 #include "G4INCLCrossSectionsMultiPionsAndResonances.hh"
46 #include "G4INCLCrossSectionsStrangeness.hh"
47 // #include <cassert>
48 
49 namespace G4INCL {
50 
51  namespace {
52  G4ThreadLocal ICrossSections *theCrossSections;
53  }
54 
55  namespace CrossSections {
56  G4double elastic(Particle const * const p1, Particle const * const p2) {
57  return theCrossSections->elastic(p1,p2);
58  }
59 
60  G4double total(Particle const * const p1, Particle const * const p2) {
61  return theCrossSections->total(p1,p2);
62  }
63 
64  G4double NDeltaToNN(Particle const * const p1, Particle const * const p2) {
65  return theCrossSections->NDeltaToNN(p1,p2);
66  }
67 
68  G4double NNToNDelta(Particle const * const p1, Particle const * const p2) {
69  return theCrossSections->NNToNDelta(p1,p2);
70  }
71 
72  G4double NNToxPiNN(const G4int xpi, Particle const * const p1, Particle const * const p2) {
73  return theCrossSections->NNToxPiNN(xpi,p1,p2);
74  }
75 
76  G4double piNToDelta(Particle const * const p1, Particle const * const p2) {
77  return theCrossSections->piNToDelta(p1,p2);
78  }
79 
80  G4double piNToxPiN(const G4int xpi, Particle const * const p1, Particle const * const p2) {
81  return theCrossSections->piNToxPiN(xpi,p1,p2);
82  }
83 
84  G4double piNToEtaN(Particle const * const p1, Particle const * const p2) {
85  return theCrossSections->piNToEtaN(p1,p2);
86  }
87 
88  G4double piNToOmegaN(Particle const * const p1, Particle const * const p2) {
89  return theCrossSections->piNToOmegaN(p1,p2);
90  }
91 
92  G4double piNToEtaPrimeN(Particle const * const p1, Particle const * const p2) {
93  return theCrossSections->piNToEtaPrimeN(p1,p2);
94  }
95 
96  G4double etaNToPiN(Particle const * const p1, Particle const * const p2) {
97  return theCrossSections->etaNToPiN(p1,p2);
98  }
99 
100  G4double etaNToPiPiN(Particle const * const p1, Particle const * const p2) {
101  return theCrossSections->etaNToPiPiN(p1,p2);
102  }
103 
104  G4double omegaNToPiN(Particle const * const p1, Particle const * const p2) {
105  return theCrossSections->omegaNToPiN(p1,p2);
106  }
107 
108  G4double omegaNToPiPiN(Particle const * const p1, Particle const * const p2) {
109  return theCrossSections->omegaNToPiPiN(p1,p2);
110  }
111 
112  G4double etaPrimeNToPiN(Particle const * const p1, Particle const * const p2) {
113  return theCrossSections->etaPrimeNToPiN(p1,p2);
114  }
115 
116  G4double NNToNNEta(Particle const * const p1, Particle const * const p2) {
117  return theCrossSections->NNToNNEta(p1,p2);
118  }
119 
120  G4double NNToNNEtaExclu(Particle const * const p1, Particle const * const p2) {
121  return theCrossSections->NNToNNEtaExclu(p1,p2);
122  }
123 
124  G4double NNToNNEtaxPi(const G4int xpi, Particle const * const p1, Particle const * const p2) {
125  return theCrossSections->NNToNNEtaxPi(xpi,p1,p2);
126  }
127 
128  G4double NNToNDeltaEta(Particle const * const p1, Particle const * const p2) {
129  return theCrossSections->NNToNDeltaEta(p1,p2);
130  }
131 
132 
133  G4double NNToNNOmega(Particle const * const p1, Particle const * const p2) {
134  return theCrossSections->NNToNNOmega(p1,p2);
135  }
136 
137  G4double NNToNNOmegaExclu(Particle const * const p1, Particle const * const p2) {
138  return theCrossSections->NNToNNOmegaExclu(p1,p2);
139  }
140 
141  G4double NNToNNOmegaxPi(const G4int xpi, Particle const * const p1, Particle const * const p2) {
142  return theCrossSections->NNToNNOmegaxPi(xpi,p1,p2);
143  }
144 
145  G4double NNToNDeltaOmega(Particle const * const p1, Particle const * const p2) {
146  return theCrossSections->NNToNDeltaOmega(p1,p2);
147  }
148 
149 
150  G4double NYelastic(Particle const * const p1, Particle const * const p2) {
151  return theCrossSections->NYelastic(p1,p2);
152  }
153 
154  G4double NKbelastic(Particle const * const p1, Particle const * const p2) {
155  return theCrossSections->NKbelastic(p1,p2);
156  }
157 
158  G4double NKelastic(Particle const * const p1, Particle const * const p2) {
159  return theCrossSections->NKelastic(p1,p2);
160  }
161 
162  G4double NNToNLK(Particle const * const p1, Particle const * const p2) {
163  return theCrossSections->NNToNLK(p1,p2);
164  }
165 
166  G4double NNToNSK(Particle const * const p1, Particle const * const p2) {
167  return theCrossSections->NNToNSK(p1,p2);
168  }
169 
170  G4double NNToNLKpi(Particle const * const p1, Particle const * const p2) {
171  return theCrossSections->NNToNLKpi(p1,p2);
172  }
173 
174  G4double NNToNSKpi(Particle const * const p1, Particle const * const p2) {
175  return theCrossSections->NNToNSKpi(p1,p2);
176  }
177 
178  G4double NNToNLK2pi(Particle const * const p1, Particle const * const p2) {
179  return theCrossSections->NNToNLK2pi(p1,p2);
180  }
181 
182  G4double NNToNSK2pi(Particle const * const p1, Particle const * const p2) {
183  return theCrossSections->NNToNSK2pi(p1,p2);
184  }
185 
186  G4double NNToNNKKb(Particle const * const p1, Particle const * const p2) {
187  return theCrossSections->NNToNNKKb(p1,p2);
188  }
189 
190  G4double NNToMissingStrangeness(Particle const * const p1, Particle const * const p2) {
191  return theCrossSections->NNToMissingStrangeness(p1,p2);
192  }
193 
194  G4double NDeltaToNLK(Particle const * const p1, Particle const * const p2) {
195  return theCrossSections->NDeltaToNLK(p1,p2);
196  }
197  G4double NDeltaToNSK(Particle const * const p1, Particle const * const p2) {
198  return theCrossSections->NDeltaToNSK(p1,p2);
199  }
200  G4double NDeltaToDeltaLK(Particle const * const p1, Particle const * const p2) {
201  return theCrossSections->NDeltaToDeltaLK(p1,p2);
202  }
203  G4double NDeltaToDeltaSK(Particle const * const p1, Particle const * const p2) {
204  return theCrossSections->NDeltaToDeltaSK(p1,p2);
205  }
206 
207  G4double NDeltaToNNKKb(Particle const * const p1, Particle const * const p2) {
208  return theCrossSections->NDeltaToNNKKb(p1,p2);
209  }
210 
211  G4double NpiToLK(Particle const * const p1, Particle const * const p2) {
212  return theCrossSections->NpiToLK(p1,p2);
213  }
214 
215  G4double NpiToSK(Particle const * const p1, Particle const * const p2) {
216  return theCrossSections->NpiToSK(p1,p2);
217  }
218 
219  G4double p_pimToSmKp(Particle const * const p1, Particle const * const p2) {
220  return theCrossSections->p_pimToSmKp(p1,p2);
221  }
222 
223  G4double p_pimToSzKz(Particle const * const p1, Particle const * const p2) {
224  return theCrossSections->p_pimToSzKz(p1,p2);
225  }
226 
227  G4double p_pizToSzKp(Particle const * const p1, Particle const * const p2) {
228  return theCrossSections->p_pizToSzKp(p1,p2);
229  }
230 
231  G4double NpiToLKpi(Particle const * const p1, Particle const * const p2) {
232  return theCrossSections->NpiToLKpi(p1,p2);
233  }
234 
235  G4double NpiToSKpi(Particle const * const p1, Particle const * const p2) {
236  return theCrossSections->NpiToSKpi(p1,p2);
237  }
238 
239  G4double NpiToLK2pi(Particle const * const p1, Particle const * const p2) {
240  return theCrossSections->NpiToLK2pi(p1,p2);
241  }
242 
243  G4double NpiToSK2pi(Particle const * const p1, Particle const * const p2) {
244  return theCrossSections->NpiToSK2pi(p1,p2);
245  }
246 
247  G4double NpiToNKKb(Particle const * const p1, Particle const * const p2) {
248  return theCrossSections->NpiToNKKb(p1,p2);
249  }
250 
251  G4double NpiToMissingStrangeness(Particle const * const p1, Particle const * const p2) {
252  return theCrossSections->NpiToMissingStrangeness(p1,p2);
253  }
254 
255  G4double NLToNS(Particle const * const p1, Particle const * const p2) {
256  return theCrossSections->NLToNS(p1,p2);
257  }
258 
259  G4double NSToNL(Particle const * const p1, Particle const * const p2) {
260  return theCrossSections->NSToNL(p1,p2);
261  }
262 
263  G4double NSToNS(Particle const * const p1, Particle const * const p2) {
264  return theCrossSections->NSToNS(p1,p2);
265  }
266 
267  G4double NKToNK(Particle const * const p1, Particle const * const p2) {
268  return theCrossSections->NKToNK(p1,p2);
269  }
270 
271  G4double NKToNKpi(Particle const * const p1, Particle const * const p2) {
272  return theCrossSections->NKToNKpi(p1,p2);
273  }
274 
275  G4double NKToNK2pi(Particle const * const p1, Particle const * const p2) {
276  return theCrossSections->NKToNK2pi(p1,p2);
277  }
278 
279  G4double NKbToNKb(Particle const * const p1, Particle const * const p2) {
280  return theCrossSections->NKbToNKb(p1,p2);
281  }
282 
283  G4double NKbToSpi(Particle const * const p1, Particle const * const p2) {
284  return theCrossSections->NKbToSpi(p1,p2);
285  }
286 
287  G4double NKbToLpi(Particle const * const p1, Particle const * const p2) {
288  return theCrossSections->NKbToLpi(p1,p2);
289  }
290 
291  G4double NKbToS2pi(Particle const * const p1, Particle const * const p2) {
292  return theCrossSections->NKbToS2pi(p1,p2);
293  }
294 
295  G4double NKbToL2pi(Particle const * const p1, Particle const * const p2) {
296  return theCrossSections->NKbToL2pi(p1,p2);
297  }
298 
299  G4double NKbToNKbpi(Particle const * const p1, Particle const * const p2) {
300  return theCrossSections->NKbToNKbpi(p1,p2);
301  }
302 
303  G4double NKbToNKb2pi(Particle const * const p1, Particle const * const p2) {
304  return theCrossSections->NKbToNKb2pi(p1,p2);
305  }
306 
307 
308  G4double calculateNNAngularSlope(G4double energyCM, G4int iso) {
309  return theCrossSections->calculateNNAngularSlope(energyCM, iso);
310  }
311 
312  G4double interactionDistancePiN(const G4double projectileKineticEnergy) {
313  ThreeVector nullVector;
314  ThreeVector unitVector(0., 0., 1.);
315 
316  Particle piPlusProjectile(PiPlus, unitVector, nullVector);
317  piPlusProjectile.setEnergy(piPlusProjectile.getMass()+projectileKineticEnergy);
318  piPlusProjectile.adjustMomentumFromEnergy();
319  Particle piZeroProjectile(PiZero, unitVector, nullVector);
320  piZeroProjectile.setEnergy(piZeroProjectile.getMass()+projectileKineticEnergy);
321  piZeroProjectile.adjustMomentumFromEnergy();
322  Particle piMinusProjectile(PiMinus, unitVector, nullVector);
323  piMinusProjectile.setEnergy(piMinusProjectile.getMass()+projectileKineticEnergy);
324  piMinusProjectile.adjustMomentumFromEnergy();
325 
326  Particle protonTarget(Proton, nullVector, nullVector);
327  Particle neutronTarget(Neutron, nullVector, nullVector);
328  const G4double sigmapipp = total(&piPlusProjectile, &protonTarget);
329  const G4double sigmapipn = total(&piPlusProjectile, &neutronTarget);
330  const G4double sigmapi0p = total(&piZeroProjectile, &protonTarget);
331  const G4double sigmapi0n = total(&piZeroProjectile, &neutronTarget);
332  const G4double sigmapimp = total(&piMinusProjectile, &protonTarget);
333  const G4double sigmapimn = total(&piMinusProjectile, &neutronTarget);
334  /* We compute the interaction distance from the largest of the pi-N cross
335  * sections. Note that this is different from INCL4.6, which just takes the
336  * average of the six, and will in general lead to a different geometrical
337  * cross section.
338  */
339  const G4double largestSigma = std::max(sigmapipp, std::max(sigmapipn, std::max(sigmapi0p, std::max(sigmapi0n, std::max(sigmapimp,sigmapimn)))));
340  const G4double interactionDistance = std::sqrt(largestSigma/Math::tenPi);
341 
342  return interactionDistance;
343  }
344 
345  G4double interactionDistanceNN(const ParticleSpecies &aSpecies, const G4double kineticEnergy) {
346 // assert(aSpecies.theType==Proton || aSpecies.theType==Neutron || aSpecies.theType==Composite);
347 // assert(aSpecies.theA>0);
348  ThreeVector nullVector;
349  ThreeVector unitVector(0.,0.,1.);
350 
351  const G4double kineticEnergyPerNucleon = kineticEnergy / aSpecies.theA;
352 
353  Particle protonProjectile(Proton, unitVector, nullVector);
354  protonProjectile.setEnergy(protonProjectile.getMass()+kineticEnergyPerNucleon);
355  protonProjectile.adjustMomentumFromEnergy();
356  Particle neutronProjectile(Neutron, unitVector, nullVector);
357  neutronProjectile.setEnergy(neutronProjectile.getMass()+kineticEnergyPerNucleon);
358  neutronProjectile.adjustMomentumFromEnergy();
359 
360  Particle protonTarget(Proton, nullVector, nullVector);
361  Particle neutronTarget(Neutron, nullVector, nullVector);
362  const G4double sigmapp = total(&protonProjectile, &protonTarget);
363  const G4double sigmapn = total(&protonProjectile, &neutronTarget);
364  const G4double sigmann = total(&neutronProjectile, &neutronTarget);
365  /* We compute the interaction distance from the largest of the NN cross
366  * sections. Note that this is different from INCL4.6, which just takes the
367  * average of the four, and will in general lead to a different geometrical
368  * cross section.
369  */
370  const G4double largestSigma = std::max(sigmapp, std::max(sigmapn, sigmann));
371  const G4double interactionDistance = std::sqrt(largestSigma/Math::tenPi);
372 
373  return interactionDistance;
374  }
375 
376  G4double interactionDistanceKN(const G4double kineticEnergy) {
377  ThreeVector nullVector;
378  ThreeVector unitVector(0.,0.,1.);
379 
380  Particle kpProjectile(KPlus, unitVector, nullVector);
381  kpProjectile.setEnergy(kpProjectile.getMass()+kineticEnergy);
382  kpProjectile.adjustMomentumFromEnergy();
383  Particle kzProjectile(KZero, unitVector, nullVector);
384  kzProjectile.setEnergy(kzProjectile.getMass()+kineticEnergy);
385  kzProjectile.adjustMomentumFromEnergy();
386 
387  Particle protonTarget(Proton, nullVector, nullVector);
388  Particle neutronTarget(Neutron, nullVector, nullVector);
389  const G4double sigmakpp = total(&kpProjectile, &protonTarget);
390  const G4double sigmakpn = total(&kpProjectile, &neutronTarget);
391  const G4double sigmakzp = total(&kzProjectile, &protonTarget);
392  const G4double sigmakzn = total(&kzProjectile, &neutronTarget);
393 
394  const G4double largestSigma = std::max(sigmakpp, std::max(sigmakpn, std::max(sigmakzp, sigmakzn)));
395  const G4double interactionDistance = std::sqrt(largestSigma/Math::tenPi);
396 
397  return interactionDistance;
398  }
399 
400  G4double interactionDistanceKbarN(const G4double kineticEnergy) {
401  ThreeVector nullVector;
402  ThreeVector unitVector(0.,0.,1.);
403 
404  Particle kmProjectile(KMinus, unitVector, nullVector);
405  kmProjectile.setEnergy(kmProjectile.getMass()+kineticEnergy);
406  kmProjectile.adjustMomentumFromEnergy();
407  Particle kzProjectile(KZeroBar, unitVector, nullVector);
408  kzProjectile.setEnergy(kzProjectile.getMass()+kineticEnergy);
409  kzProjectile.adjustMomentumFromEnergy();
410 
411  Particle protonTarget(Proton, nullVector, nullVector);
412  Particle neutronTarget(Neutron, nullVector, nullVector);
413  const G4double sigmakmp = total(&kmProjectile, &protonTarget);
414  const G4double sigmakmn = total(&kmProjectile, &neutronTarget);
415  const G4double sigmakzp = total(&kzProjectile, &protonTarget);
416  const G4double sigmakzn = total(&kzProjectile, &neutronTarget);
417 
418  const G4double largestSigma = std::max(sigmakmp, std::max(sigmakmn, std::max(sigmakzp, sigmakzn)));
419  const G4double interactionDistance = std::sqrt(largestSigma/Math::tenPi);
420 
421  return interactionDistance;
422  }
423 
424  G4double interactionDistanceYN(const G4double kineticEnergy) {
425  ThreeVector nullVector;
426  ThreeVector unitVector(0.,0.,1.);
427 
428  Particle lProjectile(Lambda, unitVector, nullVector);
429  lProjectile.setEnergy(lProjectile.getMass()+kineticEnergy);
430  lProjectile.adjustMomentumFromEnergy();
431  Particle spProjectile(SigmaPlus, unitVector, nullVector);
432  spProjectile.setEnergy(spProjectile.getMass()+kineticEnergy);
433  spProjectile.adjustMomentumFromEnergy();
434  Particle szProjectile(SigmaZero, unitVector, nullVector);
435  szProjectile.setEnergy(szProjectile.getMass()+kineticEnergy);
436  szProjectile.adjustMomentumFromEnergy();
437  Particle smProjectile(SigmaMinus, unitVector, nullVector);
438  smProjectile.setEnergy(smProjectile.getMass()+kineticEnergy);
439  smProjectile.adjustMomentumFromEnergy();
440 
441  Particle protonTarget(Proton, nullVector, nullVector);
442  Particle neutronTarget(Neutron, nullVector, nullVector);
443  const G4double sigmalp = total(&lProjectile, &protonTarget);
444  const G4double sigmaln = total(&lProjectile, &neutronTarget);
445  const G4double sigmaspp = total(&spProjectile, &protonTarget);
446  const G4double sigmaspn = total(&spProjectile, &neutronTarget);
447  const G4double sigmaszp = total(&szProjectile, &protonTarget);
448  const G4double sigmaszn = total(&szProjectile, &neutronTarget);
449  const G4double sigmasmp = total(&smProjectile, &protonTarget);
450  const G4double sigmasmn = total(&smProjectile, &neutronTarget);
451 
452  const G4double largestSigma = std::max(sigmalp, std::max(sigmaln, std::max(sigmaspp, std::max(sigmaspn, std::max(sigmaszp, std::max(sigmaszn, std::max(sigmasmp, sigmasmn)))))));
453  const G4double interactionDistance = std::sqrt(largestSigma/Math::tenPi);
454 
455  return interactionDistance;
456  }
457 
458  void setCrossSections(ICrossSections *c) {
459  theCrossSections = c;
460  }
461 
462  void deleteCrossSections() {
463  delete theCrossSections;
464  theCrossSections = NULL;
465  }
466 
467  void initialize(Config const * const theConfig) {
468  CrossSectionsType crossSections = theConfig->getCrossSectionsType();
469  if(crossSections == INCL46CrossSections)
470  setCrossSections(new CrossSectionsINCL46);
471  else if(crossSections == MultiPionsCrossSections)
472  setCrossSections(new CrossSectionsMultiPions);
473  else if(crossSections == TruncatedMultiPionsCrossSections) {
474  const G4int nMaxPi = theConfig->getMaxNumberMultipions();
475  if(nMaxPi>0)
476  setCrossSections(new CrossSectionsTruncatedMultiPions(nMaxPi));
477  else {
478  INCL_WARN("Truncated multipion cross sections were requested, but the specified maximum\n"
479  << "number of pions is <=0. Falling back to standard multipion cross-sections.\n");
480  setCrossSections(new CrossSectionsMultiPions);
481  }
482  } else if(crossSections == MultiPionsAndResonancesCrossSections)
483  setCrossSections(new CrossSectionsMultiPionsAndResonances);
484  else if(crossSections == StrangenessCrossSections)
485  setCrossSections(new CrossSectionsStrangeness);
486  }
487  }
488 }