| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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 "G4INCLCrossSectionsINCL46.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLParticleTable.hh"
41 #include "G4INCLLogger.hh"
42 // #include <cassert>
43
44 namespace G4INCL {
45
46 /* G4double elasticNNHighEnergy(const G4dou
47 return 77.0/(momentum + 1.5);
48 }
49
50 G4double elasticProtonNeutron(const G4doub
51 if(momentum < 0.450) {
52 const G4double alp = std::log(momentum
53 return 6.3555*std::exp(-3.2481*alp-0.3
54 } else if(momentum >= 0.450 && momentum
55 return (33.0 + 196.0 * std::sqrt(std::
56 } else if(momentum > 2.0) {
57 return elasticNNHighEnergy(momentum);
58 } else {
59 return 31.0/std::sqrt(momentum);
60 }
61 }
62
63 G4double elasticProtonProtonOrNeutronNeutr
64 {
65 if(momentum < 0.440) {
66 return 34.0*std::pow(momentum/0.4, -2.
67 } else if(momentum < 0.8 && momentum >=
68 return (23.5 + 1000.0*std::pow(momentu
69 } else if(momentum < 2.0) {
70 return (1250.0/(50.0 + momentum) - 4.0
71 } else {
72 return elasticNNHighEnergy(momentum);
73 }
74 }
75
76 G4double elasticNN(Particle const * const
77 G4double momentum = 0.0;
78 momentum = 0.001 * KinematicsUtils::mome
79 if((p1->getType() == Proton && p2->getTy
80 (p1->getType() == Neutron && p2->getT
81 return elasticProtonProtonOrNeutronNeu
82 } else if((p1->getType() == Proton && p2
83 (p1->getType() == Neutron && p
84 return elasticProtonNeutron(momentum);
85 } else {
86 INCL_ERROR("CrossSectionsINCL46::elast
87 << p1->print() << p2->print() <<
88 }
89 return 0.0;
90 }:*/
91
92 G4double CrossSectionsINCL46::elasticNNLeg
93
94
95 G4int i = ParticleTable::getIsospin(pa
96 + ParticleTable::getIsospin(part2->get
97
98 /* The NN cross section is parametrise
99 * of one of the nucleons. For NDelta
100 * assumption is that the cross sectio
101 * total CM energy*. Thus, we calculat
102 * we convert this value to the lab mo
103 * an NN collision*.
104 */
105
106 const G4double s = KinematicsUtils::sq
107 G4double plab = 0.001*KinematicsUtils:
108 if(plab > 2.) { // NN, Delta-Nucleon a
109 return 77./(plab+1.5);
110 }
111 else if (part1->isNucleon() && part2->
112 if (i == 0) { // pn
113 if (plab < 0.450) {
114 G4double alp=std::log(plab
115 return 6.3555*std::exp(-3.
116 }
117 else if (plab < 0.800) {
118 return (33.+196.*std::sqrt
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
127 }
128 else if (plab < 0.800) {
129 return (23.5+1000.*std::po
130 }
131 else {
132 return (1250./(50.+plab)-4
133 }
134 }
135 }
136 else { // Delta-Nucleon and Delta-De
137 if (plab < 0.440) {
138 return 34.*std::pow(plab/0.4,
139 }
140 else if (plab < 0.800) {
141 return (23.5+1000.*std::pow(pl
142 }
143 else {
144 return (1250./(50.+plab)-4.*st
145 }
146 }
147 }
148
149 G4double CrossSectionsINCL46::deltaProductio
150 G4double xs = 0.0;
151 // assert(isospin==-2 || isospin==0 || isospin
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)
161 } else if(pLab >= 1500.0 && pLab < 2000.
162 xs = (41.0 + 60.0*(momentumGeV - 0.9)*
163 } else if(pLab < 1500.0) {
164 xs = (23.5 + 24.6/(1.0 + std::exp(-10.
165 -1250.0/(momentumGeV +50.0)+4.0*
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.
171 xs = (24.2 + 8.9*momentumGeV - 31.1/st
172 } else if(pLab < 1000.0) {
173 xs = (33.0 + 196.0*std::sqrt(std::pow(
174 -31.1/std::sqrt(momentumGeV));
175 }
176 }
177
178 if(xs < 0.0) return 0.0;
179 else return xs;
180 }
181
182 G4double CrossSectionsINCL46::spnPiPlusPHE(c
183 // HE and LE pi- p and pi+ n
184 if(x <= 1750.0) {
185 return -2.33730e-06*std::pow(x, 3)+1.138
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.9169
189 +1.39907e+01*x-9360.76;
190 } else {
191 return -3.18087*std::log(x)+52.9784;
192 }
193 }
194
195 G4double CrossSectionsINCL46::spnPiMinusPHE(
196 // HE pi- p and pi+ n
197 if(x <= 1475.0) {
198 return 0.00120683*(x-1372.52)*(x-1372.52
199 } else if(x > 1475.0 && x <= 1565.0) {
200 return 1.15873e-05*x*x+49965.6/((x-1519.
201 } else if(x > 1565.0 && x <= 2400.0) {
202 return 34.0248+43262.2/((x-1681.65)*(x-1
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
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(Par
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
231 // CONST AT L
232 // COMMON/BL8/RATHR,RAMASS
233 // integer f17
234 // RATHR and RAMASS are always 0.0!!!
235
236 G4double x = KinematicsUtils::totalEnergyI
237 if(x>10000.) return 0.0; // no cross secti
238
239 G4int ipit3 = 0;
240 G4int ind2t3 = 0;
241 G4double ramass = 0.0;
242
243 if(particle1->isPion()) {
244 ipit3 = ParticleTable::getIsospin(partic
245 } else if(particle2->isPion()) {
246 ipit3 = ParticleTable::getIsospin(partic
247 }
248
249 if(particle1->isNucleon()) {
250 ind2t3 = ParticleTable::getIsospin(parti
251 } else if(particle2->isNucleon()) {
252 ind2t3 = ParticleTable::getIsospin(parti
253 }
254
255 G4double y=x*x;
256 G4double q2=(y-1076.0*1076.0)*(y-800.0*800
257 if (q2 <= 0.) {
258 return 0.0;
259 }
260 G4double q3 = std::pow(std::sqrt(q2),3);
261 G4double f3 = q3/(q3 + 5832000.); // 58320
262 G4double spnResult = 326.5/(std::pow((x-12
263 spnResult = spnResult*(1.0-5.0*ramass/1215
264 G4double cg = 4.0 + G4double(ind2t3)*G4dou
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) || (ind2t
274 spnResult=spnPiPlusPHE(x);
275 else if((ind2t3 == 1 && ipit3 == -2) ||
276 spnResult=spnPiMinusPHE(x);
277 else if(ipit3 == 0) spnResult = (spnPiPl
278 else {
279 INCL_ERROR("Unknown configuration!" <<
280 }
281 }
282
283 return spnResult;
284 }
285
286 G4double CrossSectionsINCL46::NDeltaToNN(Par
287 const G4int isospin = ParticleTable::getIs
288 if(isospin==4 || isospin==-4) return 0.0;
289
290 G4double s = KinematicsUtils::squareTotalE
291 G4double Ecm = std::sqrt(s);
292 G4int deltaIsospin;
293 G4double deltaMass;
294 if(p1->isDelta()) {
295 deltaIsospin = ParticleTable::getIsospin
296 deltaMass = p1->getMass();
297 } else {
298 deltaIsospin = ParticleTable::getIsospin
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
311 const G4double x = (s - 4.*ParticleTable::
312 (s - std::pow(ParticleTable::effectiveNu
313 const G4double y = s/(s - std::pow(deltaMa
314 /* Concerning the way we calculate the lab
315 * in CrossSections::elasticNNLegacy().
316 */
317 const G4double pLab = KinematicsUtils::mom
318 G4double result = 0.5 * x * y * deltaProdu
319 result *= 3.*(32.0 + isospin * isospin * (
320 result /= 1.0 + 0.25 * isospin * isospin;
321 return result;
322 }
323
324 G4double CrossSectionsINCL46::NNToNDelta(Par
325 // assert(p1->isNucleon() && p2->isNucleon());
326 const G4double sqrts = KinematicsUtils::to
327 if(sqrts < ParticleTable::effectivePionMas
328 return 0.0;
329 } else {
330 const G4double pLab = KinematicsUtils::m
331 const G4int isospin = ParticleTable::get
332 return deltaProduction(isospin, pLab);
333 }
334 }
335
336 G4double CrossSectionsINCL46::elastic(Partic
337 // if(!p1->isPion() && !p2->isPion())
338 if((p1->isNucleon()||p1->isDelta()) && (p2
339 // return elasticNN(p1, p2); // New i
340 return elasticNNLegacy(p1, p2); // Trans
341 else
342 return 0.0; // No pion-nucleon elastic s
343 }
344
345 G4double CrossSectionsINCL46::calculateNNAng
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-
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
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 poi
365 }
366
367
368 G4double CrossSectionsINCL46::NNToxPiNN(co
369 return 0.;
370 }
371
372 G4double CrossSectionsINCL46::piNToxPiN(co
373 return 0.;
374 }
375
376 G4double CrossSectionsINCL46::piNToEtaN(Pa
377 //
378 // Pion-Nucleon producing Eta cross se
379 //
380 return 0.;
381 }
382
383 G4double CrossSectionsINCL46::piNToOmegaN(
384 //
385 // Pion-Nucleon producing Omega cross
386 //
387 return 0.;
388 }
389
390 G4double CrossSectionsINCL46::piNToEtaPrim
391 //
392 // Pion-Nucleon producing EtaPrime cro
393 //
394 return 0.;
395 }
396
397 G4double CrossSectionsINCL46::etaNToPiN(Pa
398 //
399 // Eta-Nucleon producing Pion cross se
400 //
401 return 0.;
402 }
403
404
405 G4double CrossSectionsINCL46::etaNToPiPiN
406 //
407 // Eta-Nucleon producing Two Pions cro
408 //
409 return 0.;
410 }
411
412 G4double CrossSectionsINCL46::omegaNToPiN(
413 //
414 // Omega-Nucleon producing Pion cross
415 //
416 return 0.;
417 }
418
419 G4double CrossSectionsINCL46::omegaNToPiPi
420 //
421 // Omega-Nucleon producing Two Pions c
422 //
423 return 0.;
424 }
425
426 G4double CrossSectionsINCL46::etaPrimeNToP
427 //
428 // EtaPrime-Nucleon producing Pion cro
429 //
430 return 0.;
431 }
432
433 G4double CrossSectionsINCL46::NNToNNEta(Pa
434 //
435 // Nucleon-Nucleon producing Eta cross
436 //
437 return 0.;
438 }
439
440 G4double CrossSectionsINCL46::NNToNNEtaEx
441 //
442 // Nucleon-Nucleon producing Eta cross
443 //
444 return 0.;
445 }
446
447 G4double CrossSectionsINCL46::NNToNNEtaxPi
448 return 0.;
449 }
450
451 G4double CrossSectionsINCL46::NNToNDeltaE
452 //
453 // Nucleon-Nucleon producing N-Delta-E
454 //
455 return 0.;
456 }
457
458 G4double CrossSectionsINCL46::NNToNNOmega(
459 //
460 // Nucleon-Nucleon producing Omega cro
461 //
462 return 0.;
463 }
464
465 G4double CrossSectionsINCL46::NNToNNOmegaE
466 //
467 // Nucleon-Nucleon producing Omega cro
468 //
469 return 0.;
470 }
471
472 G4double CrossSectionsINCL46::NNToNNOmegax
473 return 0.;
474 }
475
476 G4double CrossSectionsINCL46::NNToNDeltaOm
477 //
478 // Nucleon-Nucleon producing N-Delta-Ome
479 //
480 return 0.;
481 }
482
483
484
485 G4double CrossSectionsINCL46::NYelastic(Pa
486 //
487 // Hyperon-Nucleon elastic cross
488 //
489 return 0.;
490 }
491
492 G4double CrossSectionsINCL46::NKelastic(Pa
493 //
494 // Kaon-Nucleon elastic cross sec
495 //
496 return 0.;
497 }
498
499 G4double CrossSectionsINCL46::NKbelastic(P
500 //
501 // antiKaon-Nucleon elastic cross
502 //
503 return 0.;
504 }
505
506 G4double CrossSectionsINCL46::NNToNLK(Partic
507 //
508 // Nucleon-Nucleon producing N-La
509 //
510 return 0.;
511 }
512
513 G4double CrossSectionsINCL46::NNToNSK(Part
514 //
515 // Nucleon-Nucleon producing N-Si
516 //
517 return 0.;
518 }
519
520 G4double CrossSectionsINCL46::NNToNLKpi(Pa
521 //
522 // Nucleon-Nucleon producing N-La
523 //
524 return 0.;
525 }
526
527 G4double CrossSectionsINCL46::NNToNSKpi(Pa
528 //
529 // Nucleon-Nucleon producing N-Si
530 //
531 return 0.;
532 }
533
534 G4double CrossSectionsINCL46::NNToNLK2pi(P
535 //
536 // Nucleon-Nucleon producing N-Lam
537 //
538 return 0.;
539 }
540
541 G4double CrossSectionsINCL46::NNToNSK2pi(P
542 //
543 // Nucleon-Nucleon producing N-Si
544 //
545 return 0.;
546 }
547
548 G4double CrossSectionsINCL46::NNToNNKKb(Pa
549 //
550 // Nucleon-Nucleon producing Nucl
551 //
552 return 0.;
553 }
554
555 G4double CrossSectionsINCL46::NNToMissingS
556 //
557 // Nucleon-Nucleon missing strang
558 //
559 return 0.;
560 }
561
562 G4double CrossSectionsINCL46::NDeltaToNLK(
563 // Nucleon-Delta producing Nucleon Lam
564 return 0;
565 }
566 G4double CrossSectionsINCL46::NDeltaToNSK(
567 // Nucleon-Delta producing Nucleon Sig
568 return 0;
569 }
570 G4double CrossSectionsINCL46::NDeltaToDelt
571 // Nucleon-Delta producing Delta Lambd
572 return 0;
573 }
574 G4double CrossSectionsINCL46::NDeltaToDelt
575 // Nucleon-Delta producing Delta Sigma
576 return 0;
577 }
578
579 G4double CrossSectionsINCL46::NDeltaToNNKK
580 // Nucleon-Delta producing Nucleon-Nuc
581 return 0;
582 }
583
584 G4double CrossSectionsINCL46::NpiToLK(Part
585 //
586 // Pion-Nucleon producing Lambda-
587 //
588 return 0.;
589 }
590
591 G4double CrossSectionsINCL46::NpiToSK(Part
592 //
593 // Pion-Nucleon producing Sigma-K
594 //
595 return 0.;
596 }
597 G4double CrossSectionsINCL46::p_pimToSmKp(
598 return 0.;
599 }
600 G4double CrossSectionsINCL46::p_pimToSzKz(
601 return 0.;
602 }
603 G4double CrossSectionsINCL46::p_pizToSzKp(
604 return 0.;
605 }
606
607 G4double CrossSectionsINCL46::NpiToLKpi(Pa
608 //
609 // Pion-Nucleon producing Lambda-
610 //
611 return 0.;
612 }
613
614 G4double CrossSectionsINCL46::NpiToSKpi(Pa
615 //
616 // Pion-Nucleon producing Sigma-K
617 //
618 return 0.;
619 }
620
621 G4double CrossSectionsINCL46::NpiToLK2pi(P
622 //
623 // Pion-Nucleon producing Lambda-
624 //
625 return 0.;
626 }
627
628 G4double CrossSectionsINCL46::NpiToSK2pi(P
629 //
630 // Pion-Nucleon producing Lambda-
631 //
632 return 0.;
633 }
634
635 G4double CrossSectionsINCL46::NpiToNKKb(Pa
636 //
637 // Pion-Nucleon producing Nucleon
638 //
639 return 0.;
640 }
641
642 G4double CrossSectionsINCL46::NpiToMissing
643 //
644 // Pion-Nucleon missing strangene
645 //
646 return 0.;
647 }
648
649 G4double CrossSectionsINCL46::NLToNS(Parti
650 //
651 // Nucleon-Hyperon multiplet chan
652 //
653 return 0.;
654 }
655
656 G4double CrossSectionsINCL46::NSToNL(Parti
657 //
658 // Nucleon-Sigma quasi-elastic cr
659 //
660 return 0.;
661 }
662
663 G4double CrossSectionsINCL46::NSToNS(Parti
664 //
665 // Nucleon-Sigma quasi-elastic cr
666 //
667 return 0.;
668 }
669
670 G4double CrossSectionsINCL46::NKToNK(Parti
671 //
672 // Nucleon-Kaon quasi-elastic cro
673 //
674 return 0.;
675 }
676
677 G4double CrossSectionsINCL46::NKToNKpi(Par
678 //
679 // Nucleon-Kaon producing Nucleon
680 //
681 return 0.;
682 }
683
684 G4double CrossSectionsINCL46::NKToNK2pi(Pa
685 //
686 // Nucleon-Kaon producing Nucleon
687 //
688 return 0.;
689 }
690
691 G4double CrossSectionsINCL46::NKbToNKb(Par
692 //
693 // Nucleon-antiKaon quasi-elastic
694 //
695 return 0.;
696 }
697
698 G4double CrossSectionsINCL46::NKbToSpi(Par
699 //
700 // Nucleon-antiKaon producing Sig
701 //
702 return 0.;
703 }
704
705 G4double CrossSectionsINCL46::NKbToLpi(Par
706 //
707 // Nucleon-antiKaon producing Lam
708 //
709 return 0.;
710 }
711
712 G4double CrossSectionsINCL46::NKbToS2pi(Pa
713 //
714 // Nucleon-antiKaon producing Sig
715 //
716 return 0.;
717 }
718
719 G4double CrossSectionsINCL46::NKbToL2pi(Pa
720 //
721 // Nucleon-antiKaon producing Lam
722 //
723 return 0.;
724 }
725
726 G4double CrossSectionsINCL46::NKbToNKbpi(P
727 //
728 // Nucleon-antiKaon producing Nuc
729 //
730 return 0.;
731 }
732
733 G4double CrossSectionsINCL46::NKbToNKb2pi(
734 //
735 // Nucleon-antiKaon producing Nuc
736 //
737 return 0.;
738 }
739
740 G4double CrossSectionsINCL46::NNbarElastic
741 //
742 // Nucleon-AntiNucleon to Nucleon-
743 //
744 return 0.;
745 }
746
747 G4double CrossSectionsINCL46::NNbarCEX(Par
748 //
749 // Nucleon-AntiNucleon charge exch
750 //
751 return 0.;
752 }
753
754 G4double CrossSectionsINCL46::NNbarToLLbar
755 //
756 // Nucleon-AntiNucleon to Lambda-A
757 //
758 return 0.;
759 }
760
761 G4double CrossSectionsINCL46::NNbarToNNbar
762 //
763 // Nucleon-AntiNucleon to Nucleon-
764 //
765 return 0.;
766 }
767
768 G4double CrossSectionsINCL46::NNbarToNNbar
769 //
770 // Nucleon-AntiNucleon to Nucleon-
771 //
772 return 0.;
773 }
774
775 G4double CrossSectionsINCL46::NNbarToNNbar
776 //
777 // Nucleon-AntiNucleon to Nucleon-
778 //
779 return 0.;
780 }
781
782 G4double CrossSectionsINCL46::NNbarToAnnih
783 //
784 // Nucleon-AntiNucleon total annih
785 //
786 return 0.;
787 }
788 } // namespace G4INCL
789
790