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