| 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 /** \file G4INCLCrossSectionsMultiPionsAndReso
39 * \brief Multipion and mesonic Resonances cro
40 *
41 * \date 4th February 2014
42 * \author Jean-Christophe David
43 */
44
45 #include "G4INCLCrossSectionsMultiPionsAndReso
46 #include "G4INCLKinematicsUtils.hh"
47 #include "G4INCLParticleTable.hh"
48 // #include <cassert>
49
50 namespace G4INCL {
51
52 template<G4int N>
53 struct BystrickyEvaluator {
54 static G4double eval(const G4double pL
55 const G4double pMeV = pLab*1E3;
56 const G4double ekin=std::sqrt(Part
57 const G4double xrat=ekin*oneOverTh
58 const G4double x=std::log(xrat);
59 return HornerEvaluator<N>::eval(x,
60 }
61 };
62
63 const G4int CrossSectionsMultiPionsAndReso
64 const G4int CrossSectionsMultiPionsAndReso
65
66 const G4double CrossSectionsMultiPionsAndR
67 const G4double CrossSectionsMultiPionsAndR
68 const G4double CrossSectionsMultiPionsAndR
69 const G4double CrossSectionsMultiPionsAndR
70 const G4double CrossSectionsMultiPionsAndR
71 const G4double CrossSectionsMultiPionsAndR
72 const G4double CrossSectionsMultiPionsAndR
73 const G4double CrossSectionsMultiPionsAndR
74 const G4double CrossSectionsMultiPionsAndR
75 const G4double CrossSectionsMultiPionsAndR
76
77 CrossSectionsMultiPionsAndResonances::Cros
78 s11pzHC(-2.228000000000294018,8.7560000000
79 s01ppHC(2.0570000000126518344,-6.029000000
80 s01pzHC(0.18030000000000441851,7.870099999
81 s11pmHC(0.20590000000000031866,3.345099999
82 s12pmHC(-0.77235999999999901328,4.26265999
83 s12ppHC(-0.75724999999999975664,2.09343999
84 s12zzHC(-0.89599999999996965072,7.88299999
85 s02pzHC(-1.0579999999999967036,11.11399999
86 s02pmHC(2.4009000000012553286,-7.768000000
87 s12mzHC(-0.21858699999999976269,1.91489999
88 {
89 }
90
91 G4double CrossSectionsMultiPionsAndResonan
92 G4double inelastic;
93 if(p1->isNucleon() && p2->isNucleon())
94 return CrossSectionsMultiPions::NN
95 } else if((p1->isNucleon() && p2->isDe
96 (p1->isDelta() && p2->isNucl
97 inelastic = CrossSectionsMultiPion
98 } else if((p1->isNucleon() && p2->isPi
99 (p1->isPion() && p2->isNucle
100 return CrossSectionsMultiPions::pi
101 } else if((p1->isNucleon() && p2->isEt
102 (p1->isEta() && p2->isNucleo
103 inelastic = etaNToPiN(p1,p2) + eta
104 } else if((p1->isNucleon() && p2->isOm
105 (p1->isOmega() && p2->isNucl
106 inelastic = omegaNInelastic(p1,p2)
107 } else if((p1->isNucleon() && p2->isEt
108 (p1->isEtaPrime() && p2->isN
109 inelastic = etaPrimeNToPiN(p1,p2);
110 } else {
111 inelastic = 0.;
112 }
113
114 return inelastic + elastic(p1, p2);
115 }
116
117
118 G4double CrossSectionsMultiPionsAndResonan
119 if((p1->isNucleon()||p1->isDelta()) &&
120 return CrossSectionsMultiPions::el
121 }
122 else if ((p1->isNucleon() && p2->isPio
123 return CrossSectionsMultiPions::el
124 }
125 else if ((p1->isNucleon() && p2->isEta
126 return etaNElastic(p1, p2);
127 }
128 else if ((p1->isNucleon() && p2->isOme
129 return omegaNElastic(p1, p2);
130 }
131 else {
132 return 0.0;
133 }
134 }
135
136
137 G4double CrossSectionsMultiPionsAndResonan
138 //
139 // pion-Nucleon producing xpi pion
140 //
141 // assert(xpi>1 && xpi<=nMaxPiPiN);
142 // assert((particle1->isNucleon() && particle2
143
144 const G4double oldXS2Pi=CrossSectionsM
145 const G4double oldXS3Pi=CrossSectionsM
146 const G4double oldXS4Pi=CrossSectionsM
147 const G4double xsEta=piNToEtaN(particl
148 const G4double xsOmega=piNToOmegaN(par
149 G4double newXS2Pi=0.;
150 G4double newXS3Pi=0.;
151 G4double newXS4Pi=0.;
152
153 if (xpi == 2) {
154 if (oldXS4Pi != 0.)
155 newXS2Pi=oldXS2Pi;
156 else if (oldXS3Pi != 0.) {
157 newXS3Pi=oldXS3Pi-xsEta-xsOmeg
158 if (newXS3Pi < 1.e-09)
159 newXS2Pi=oldXS2Pi-(xsEta+x
160 else
161 newXS2Pi=oldXS2Pi;
162 }
163 else {
164 newXS2Pi=oldXS2Pi-xsEta-xsOmeg
165 if (newXS2Pi < 1.e-09)
166 newXS2Pi=0.;
167 }
168 return newXS2Pi;
169 }
170 else if (xpi == 3) {
171 if (oldXS4Pi != 0.) {
172 newXS4Pi=oldXS4Pi-xsEta-xsOmeg
173 if (newXS4Pi < 1.e-09)
174 newXS3Pi=oldXS3Pi-(xsEta+x
175 else
176 newXS3Pi=oldXS3Pi;
177 }
178 else {
179 newXS3Pi=oldXS3Pi-xsEta-xsOmeg
180 if (newXS3Pi < 1.e-09)
181 newXS3Pi=0.;
182 }
183 return newXS3Pi;
184 }
185 else if (xpi == 4) {
186 newXS4Pi=oldXS4Pi-xsEta-xsOmega;
187 if (newXS4Pi < 1.e-09)
188 newXS4Pi=0.;
189 return newXS4Pi;
190 }
191 else // should never reach this point
192 return 0.;
193 }
194
195 G4double CrossSectionsMultiPionsAndResonan
196 //
197 // Pion-Nucleon producing Eta cros
198 //
199 // assert((particle1->isNucleon() && particle2
200
201 G4double sigma;
202 sigma=piMinuspToEtaN(particle1,particl
203
204 const G4int isoin = ParticleTable::get
205
206 if (isoin == -1) {
207 if ((particle1->getType()) == Prot
208 else return 0.5 * sigma;
209 }
210 else if (isoin == 1) {
211 if ((particle1->getType()) == Neut
212 else return 0.5 * sigma;
213 }
214 else return 0. ; // should never retur
215
216 // return sigma;
217 }
218
219 G4double CrossSectionsMultiPionsAndResonan
220 //
221 // Pion-Nucleon producing Omega cr
222 //
223 // assert((particle1->isNucleon() && particle2
224
225 G4double sigma;
226 sigma=piMinuspToOmegaN(particle1,parti
227
228 const G4int isoin = ParticleTable::get
229
230 if (isoin == -1) {
231 if ((particle1->getType()) == Prot
232 else return 0.5 * sigma;
233 }
234 else if (isoin == 1) {
235 if ((particle1->getType()) == Neut
236 else return 0.5 * sigma;
237 }
238 else return 0. ; // should never retur
239
240 // return sigma;
241 }
242
243 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT
244 G4double CrossSectionsMultiPionsAndResonan
245 #else
246 G4double CrossSectionsMultiPionsAndResonan
247 #endif
248 //
249 // Pion-Nucleon producing EtaPrime
250 //
251 // assert((particle1->isNucleon() && particle2
252
253 return 0.;
254 }
255
256 G4double CrossSectionsMultiPionsAndResonan
257 //
258 // Eta-Nucleon producing Pion cros
259 //
260 // assert((particle1->isNucleon() && particle2
261
262 const Particle *eta;
263 const Particle *nucleon;
264
265 if (particle1->isEta()) {
266 eta = particle1;
267 nucleon = particle2;
268 }
269 else {
270 eta = particle2;
271 nucleon = particle1;
272 }
273
274 const G4double pLab = KinematicsUtils:
275 G4double sigma=0.;
276
277 if (pLab <= 574.)
278 sigma= 1.511147E-13*std::pow(pLab,6)-
279 else if (pLab <= 850.)
280 sigma= -8.00018E-14*std::pow(pLab,6)+
281 else if (pLab <= 1300.)
282 sigma= 6.56364E-09*std::pow(pLab,3)-
283 else {
284 G4double ECM=KinematicsUtils::tota
285 G4double massPiZero=ParticleTable:
286 G4double massPiMinus=ParticleTable
287 G4double massProton=ParticleTable:
288 G4double masseta;
289 G4double massnucleon;
290 G4double pCM_eta;
291 masseta=eta->getMass();
292 massnucleon=nucleon->getMass();
293 pCM_eta=KinematicsUtils::momentumI
294 G4double pCM_PiZero=KinematicsUtil
295 G4double pCM_PiMinus=KinematicsUti
296 sigma = (piMinuspToEtaN(ECM)/2.) *
297 }
298 if (sigma < 0.) sigma=0.;
299
300 return sigma;
301 }
302
303 G4double CrossSectionsMultiPionsAndResonan
304 //
305 // Eta-Nucleon producing Two Pions
306 //
307 // assert((particle1->isNucleon() && particle2
308
309 G4double sigma=0.;
310
311 const Particle *eta;
312 const Particle *nucleon;
313
314 if (particle1->isEta()) {
315 eta = particle1;
316 nucleon = particle2;
317 }
318 else {
319 eta = particle2;
320 nucleon = particle1;
321 }
322
323 const G4double pLab = KinematicsUtils:
324
325 if (pLab < 450.)
326 sigma = 2.01854221E-13*std::pow(pL
327 else if (pLab < 600.)
328 sigma = 2.01854221E-13*std::pow(45
329 else if (pLab <= 1300.)
330 sigma = -6.32793049e-16*std::pow(p
331 1.37055547e-05*std::pow(pLab,3) -
332 else
333 sigma = etaNToPiN(particle1,partic
334
335 if (sigma < 0.) sigma = 0.;
336 return sigma; // Parameterization from
337 }
338
339
340 G4double CrossSectionsMultiPionsAndResonan
341 //
342 // Eta-Nucleon elastic cross secti
343 //
344 // assert((particle1->isNucleon() && particle2
345
346 G4double sigma=0.;
347
348 const Particle *eta;
349 const Particle *nucleon;
350
351 if (particle1->isEta()) {
352 eta = particle1;
353 nucleon = particle2;
354 }
355 else {
356 eta = particle2;
357 nucleon = particle1;
358 }
359
360 const G4double pLab = KinematicsUtils:
361
362 if (pLab < 700.)
363 sigma = 3.6838e-15*std::pow(pLab,6
364 4.3222e-06*std::pow(pLab,3) + 7.91
365 else if (pLab < 1400.)
366 sigma = 3.562630e-16*std::pow(pLab
367 9.667078e-06*std::pow(pLab,3) + 7.
368 else if (pLab < 2025.)
369 sigma = -1.041950E-03*pLab + 2.110
370 else
371 sigma=0.;
372
373 if (sigma < 0.) sigma = 0.;
374 return sigma; // Parameterization
375 }
376
377 G4double CrossSectionsMultiPionsAndResonan
378 //
379 // Omega-Nucleon inelastic cross s
380 //
381 // assert((particle1->isNucleon() && particle2
382
383 G4double sigma=0.;
384
385 const Particle *omega;
386 const Particle *nucleon;
387
388 if (particle1->isOmega()) {
389 omega = particle1;
390 nucleon = particle2;
391 }
392 else {
393 omega = particle2;
394 nucleon = particle1;
395 }
396
397 const G4double pLab = KinematicsUtils:
398
399 sigma = 20. + 4.0/pLab; // Eq.(24)
400
401 return sigma;
402 }
403
404
405 G4double CrossSectionsMultiPionsAndResonan
406 //
407 // Omega-Nucleon elastic cross sec
408 //
409 // assert((particle1->isNucleon() && particle2
410
411 G4double sigma=0.;
412
413 const Particle *omega;
414 const Particle *nucleon;
415
416 if (particle1->isOmega()) {
417 omega = particle1;
418 nucleon = particle2;
419 }
420 else {
421 omega = particle2;
422 nucleon = particle1;
423 }
424
425 const G4double pLab = KinematicsUtils:
426
427 sigma = 5.4 + 10.*std::exp(-0.6*pLab);
428
429 return sigma;
430 }
431
432
433 G4double CrossSectionsMultiPionsAndResonan
434 //
435 // Omega-Nucleon producing Pion cr
436 //
437 // assert((particle1->isNucleon() && particle2
438
439 G4double ECM=KinematicsUtils::totalEne
440
441 G4double massPiZero=ParticleTable::get
442 G4double massPiMinus=ParticleTable::ge
443 G4double massProton=ParticleTable::get
444
445 G4double massomega;
446 G4double massnucleon;
447 G4double pCM_omega;
448 G4double pLab_omega;
449
450 G4double sigma=0.;
451
452 if (particle1->isOmega()) {
453 massomega=particle1->getMass();
454 massnucleon=particle2->getMass();
455 }
456 else {
457 massomega=particle2->getMass();
458 massnucleon=particle1->getMass();
459 }
460 pCM_omega=KinematicsUtils::momentumInC
461 pLab_omega=KinematicsUtils::momentumIn
462
463 G4double pCM_PiZero=KinematicsUtils::m
464 G4double pCM_PiMinus=KinematicsUtils::
465
466 sigma = (piMinuspToOmegaN(ECM)/2.) * s
467 + piMinuspToOmegaN(ECM) * std::pow((pC
468
469 if (sigma > omegaNInelastic(particle1,
470 // if (sigma > omegaNInelastic(particle
471 sigma = omegaNInelastic(particle1,
472 }
473
474 return sigma;
475 }
476
477
478 G4double CrossSectionsMultiPionsAndResonan
479 //
480 // Omega-Nucleon producing 2 PionS
481 //
482 // assert((particle1->isNucleon() && particle2
483
484 G4double sigma=0.;
485
486 sigma = omegaNInelastic(particle1,part
487
488 return sigma;
489 }
490
491
492 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT
493 G4double CrossSectionsMultiPionsAndResonance
494 #else
495 G4double CrossSectionsMultiPionsAndResonance
496 #endif
497 //
498 // EtaPrime-Nucleon producing Pion
499 //
500 // assert((particle1->isNucleon() && particle2
501
502 return 0.;
503 }
504
505 G4double CrossSectionsMultiPionsAndResonan
506 //
507 // Pion-Nucleon producing Eta cros
508 //
509 // assert((particle1->isNucleon() && particle2
510
511 G4double masspion;
512 G4double massnucleon;
513 if (particle1->isPion()) {
514 masspion=particle1->getMass();
515 massnucleon=particle2->getMass();
516 }
517 else {
518 masspion=particle2->getMass();
519 massnucleon=particle1->getMass();
520 }
521
522 G4double ECM=KinematicsUtils::totalEne
523 G4double plab=KinematicsUtils::momentu
524
525 G4double sigma;
526
527 // new parameterization (JCD) - end of februar
528 if (ECM < 1486.5) sigma=0.;
529 else
530 {
531 if (ECM < 1535.)
532 {
533 sigma = -0.0000003689197974814
534 }
535 else if (ECM < 1670.)
536 {
537 sigma = -0.0000000337986446*st
538 }
539 else if (ECM < 1714.)
540 {
541 sigma = 0.000003737765*std::p
542 }
543 else sigma=1.47*std::pow(plab, -1.
544 }
545 //
546
547 return sigma;
548 }
549
550 G4double CrossSectionsMultiPionsAndResonan
551 //
552 // Pion-Nucleon producing Eta cros
553 //
554
555 const G4double masspion = ParticleTabl
556 const G4double massnucleon = ParticleT
557
558 G4double plab=KinematicsUtils::momentu
559
560 G4double sigma;
561
562 // new parameterization (JCD) - end of februar
563 if (ECM < 1486.5) sigma=0.;
564 else
565 {
566 if (ECM < 1535.)
567 {
568 sigma = -0.0000003689197974814
569 }
570 else if (ECM < 1670.)
571 {
572 sigma = -0.0000000337986446*st
573 }
574 else if (ECM < 1714.)
575 {
576 sigma = 0.000003737765*std::p
577 }
578 else sigma=1.47*std::pow(plab, -1.
579 }
580
581 return sigma;
582 }
583
584 G4double CrossSectionsMultiPionsAndResonan
585 //
586 // Pion-Nucleon producing Omega cr
587 //
588 // assert((particle1->isNucleon() && particle2
589 //jcd to be removed
590 // return 0.;
591 //jcd
592
593 // G4double param=1.095 ; // Deneye (Th
594 G4double param=1.0903 ; // JCD (thresh
595
596 G4double masspion;
597 G4double massnucleon;
598 if (particle1->isPion()) {
599 masspion=particle1->getMass();
600 massnucleon=particle2->getMass();
601 }
602 else {
603 masspion=particle2->getMass();
604 massnucleon=particle1->getMass();
605 }
606 G4double ECM=KinematicsUtils::totalEne
607 G4double plab=KinematicsUtils::momentu
608
609 G4double sigma;
610 if (plab < param) sigma=0.;
611 else sigma=13.76*(plab-param)/(std::po
612
613 return sigma;
614 }
615 G4double CrossSectionsMultiPionsAndResonan
616 //
617 // Pion-Nucleon producing Omega cr
618 //
619 //jcd to be removed
620 // return 0.;
621 //jcd
622
623 // G4double param=1.095 ; // Deneye (Th
624 G4double param=1.0903 ; // JCD (thresh
625
626 const G4double masspion = ParticleTabl
627 const G4double massnucleon = ParticleT
628
629 G4double plab=KinematicsUtils::momentu
630
631 G4double sigma;
632 if (plab < param) sigma=0.;
633 else sigma=13.76*(plab-param)/(std::po
634
635 return sigma;
636 }
637
638 G4double CrossSectionsMultiPionsAndResonan
639
640 const G4double Ecm=0.001*ener;
641 G4double sNNEta; // pp->pp+eta(+X)
642 G4double sNNEta1; // np->np+eta(+X)
643 G4double sNNEta2; // np->d+eta (d will
644 G4double x=Ecm*Ecm/5.88;
645
646 //jcd
647 if (Ecm >= 3.05) {
648 sNNEta = 2.5*std::pow((x-1.),1.47)
649 }
650 else if(Ecm >= 2.6) {
651 sNNEta = -327.29*Ecm*Ecm*Ecm + 287
652 if (sNNEta <= NNToNNEtaExcluIso(en
653 }
654 else {
655 sNNEta = NNToNNEtaExcluIso(ener, 2
656 }
657 //jcd
658 if (sNNEta < 1.e-9) sNNEta = 0.;
659
660 if (iso != 0) {
661 return sNNEta/1000.; // parameteri
662 }
663
664 if(Ecm >= 6.25) {
665 sNNEta1 = sNNEta;
666 }
667 else if (Ecm >= 2.6) {
668 sNNEta1 = sNNEta*std::exp(-(-5.531
669 }
670 else if (Ecm >= 2.525) { // = exclusiv
671 sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ec
672 }
673 else { // = exclusive pn
674 sNNEta1 = 17570.217219*Ecm*Ecm - 8
675 }
676
677 sNNEta2 = -10220.89518466*Ecm*Ecm+5122
678 if (sNNEta2 < 0.) sNNEta2=0.;
679
680 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
681
682 G4double Mn=ParticleTable::getRealMass
683 G4double Mp=ParticleTable::getRealMass
684 G4double Meta=ParticleTable::getRealMa
685 if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta
686
687 return sNNEta/1000.; // parameterizati
688 }
689
690
691 G4double CrossSectionsMultiPionsAndResonan
692
693 const G4double ener=KinematicsUtils::t
694 const G4int iso=ParticleTable::getIsos
695
696 if (iso != 0) {
697 return NNToNNEtaIso(ener, iso);
698 }
699 else {
700 return 0.5*(NNToNNEtaIso(ener, 0)+
701 }
702 }
703
704 G4double CrossSectionsMultiPionsAndResonan
705
706 const G4double Ecm=0.001*ener;
707 G4double sNNEta; // pp->pp+eta
708 G4double sNNEta1; // np->np+eta
709 G4double sNNEta2; // np->d+eta (d wil
710
711 if(Ecm>=3.875) { // By hand (JCD)
712 sNNEta = -13.008*Ecm*Ecm + 84.531*
713 }
714 else if(Ecm>=2.725) { // By hand (JCD)
715 sNNEta = -913.2809*std::pow(Ecm,5)
716 }
717 else if(Ecm>=2.575) { // By hand (JCD)
718 sNNEta = -2640.3*Ecm*Ecm + 14692*E
719 }
720 else {
721 sNNEta = -147043.497285*std::pow(E
722 }
723
724 G4double Mn=ParticleTable::getRealMass
725 G4double Mp=ParticleTable::getRealMass
726 G4double Meta=ParticleTable::getRealMa
727 G4double Thr0=0.;
728 if (iso > 0) {
729 Thr0=2.*Mp+Meta;
730 }
731 else if (iso < 0) {
732 Thr0=2.*Mn+Meta;
733 }
734 else {
735 Thr0=Mn+Mp+Meta;
736 }
737
738 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNE
739
740 if (iso != 0) {
741 return sNNEta/1000.; // parameteri
742 }
743
744 if(Ecm>=3.9) {
745 sNNEta1 = sNNEta;
746 }
747 else if (Ecm >= 3.5) {
748 sNNEta1 = -1916.2*Ecm*Ecm*Ecm + 21
749 }
750 else if (Ecm >= 2.525) {
751 sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ec
752 }
753 else {
754 sNNEta1 = 17570.217219*Ecm*Ecm - 8
755 }
756
757 sNNEta2 = -10220.89518466*Ecm*Ecm+5122
758 if (sNNEta2 < 0.) sNNEta2=0.;
759
760 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
761 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNE
762
763 return sNNEta/1000.; // parameterizati
764
765 }
766
767 G4double CrossSectionsMultiPionsAndResonan
768
769 const G4double ener=KinematicsUtils::t
770 const G4int iso=ParticleTable::getIsos
771
772 if (iso != 0) {
773 return NNToNNEtaExcluIso(ener, iso
774 }
775 else {
776 return 0.5*(NNToNNEtaExcluIso(ener
777 }
778 }
779
780
781 G4double CrossSectionsMultiPionsAndResonan
782
783 const G4double Ecm=0.001*ener;
784 G4double sNNOmega; // pp->pp+eta(+X)
785 G4double sNNOmega1; // np->np+eta(+X)
786 G4double x=Ecm*Ecm/7.06;
787
788 if(Ecm>4.0) {
789 sNNOmega = 2.5*std::pow(x-1, 1.47)
790 }
791 else if(Ecm>2.802) { // 2802 MeV -> th
792 sNNOmega = (568.5254*Ecm*Ecm - 269
793 if (sNNOmega <= NNToNNOmegaExcluIso(en
794 }
795 else {
796 sNNOmega = NNToNNOmegaExcluIso(ene
797 }
798
799 if (sNNOmega < 1.e-9) sNNOmega = 0.;
800
801 if (iso != 0) {
802 return sNNOmega;
803 }
804
805 sNNOmega1 = 3.*sNNOmega; // 3.0: ratio
806
807 sNNOmega = 2.*sNNOmega1-sNNOmega;
808
809 if (sNNOmega < 1.e-9) sNNOmega = 0.;
810
811 return sNNOmega;
812 }
813
814
815 G4double CrossSectionsMultiPionsAndResonan
816
817 const G4double ener=KinematicsUtils::t
818 const G4int iso=ParticleTable::getIsos
819 //jcd to be removed
820 // return 0.;
821 //jcd
822 if (iso != 0) {
823 return NNToNNOmegaIso(ener, iso);
824 }
825 else {
826 return 0.5*(NNToNNOmegaIso(ener, 0
827 }
828 }
829
830 G4double CrossSectionsMultiPionsAndResonan
831
832 const G4double Ecm=0.001*ener;
833 G4double sNNOmega; // pp->pp+eta
834 G4double sNNOmega1; // np->np+eta
835 G4double sthroot=std::sqrt(7.06);
836
837 if(Ecm>=3.0744) { // By hand (JCD)
838 sNNOmega = 330.*(Ecm-sthroot)/(1.0
839 }
840 else if(Ecm>=2.65854){
841 sNNOmega = -1208.09757*std::pow(Ec
842 }
843 else {
844 sNNOmega = 0. ;
845 }
846
847 G4double Mn=ParticleTable::getRealMass
848 G4double Mp=ParticleTable::getRealMass
849 G4double Momega=ParticleTable::getReal
850 G4double Thr0=0.;
851 if (iso > 0) {
852 Thr0=2.*Mp+Momega;
853 }
854 else if (iso < 0) {
855 Thr0=2.*Mn+Momega;
856 }
857 else {
858 Thr0=Mn+Mp+Momega;
859 }
860
861 if (sNNOmega < 1.e-9 || Ecm < Thr0) sN
862
863 if (iso != 0) {
864 return sNNOmega/1000.; // paramete
865 }
866
867 sNNOmega1 = 3*sNNOmega; // 3.0: ratio
868
869 sNNOmega = 2*sNNOmega1-sNNOmega;
870 if (sNNOmega < 1.e-9 || Ecm < Thr0) sN
871
872 return sNNOmega/1000.; // parameteriza
873 }
874
875 G4double CrossSectionsMultiPionsAndResonan
876 //jcd to be removed
877 // return 0.;
878 //jcd
879
880 const G4double ener=KinematicsUtils::t
881 const G4int iso=ParticleTable::getIsos
882
883 if (iso != 0) {
884 return NNToNNOmegaExcluIso(ener, i
885 }
886 else {
887 return 0.5*(NNToNNOmegaExcluIso(en
888 }
889 }
890
891
892 G4double CrossSectionsMultiPionsAndResonan
893 //
894 // Nucleon-Nucleon producing xpi p
895 //
896 // assert(xpi>0 && xpi<=nMaxPiNN);
897 // assert(particle1->isNucleon() && particle2-
898
899 G4double oldXS1Pi=CrossSectionsMultiPi
900 G4double oldXS2Pi=CrossSectionsMultiPi
901 G4double oldXS3Pi=CrossSectionsMultiPi
902 G4double oldXS4Pi=CrossSectionsMultiPi
903 G4double xsEtaOmega=NNToNNEta(particle
904 G4double newXS1Pi=0.;
905 G4double newXS2Pi=0.;
906 G4double newXS3Pi=0.;
907 G4double newXS4Pi=0.;
908
909 if (xpi == 1) {
910 if (oldXS4Pi != 0. || oldXS3Pi !=
911 newXS1Pi=oldXS1Pi;
912 else if (oldXS2Pi != 0.) {
913 newXS2Pi=oldXS2Pi-xsEtaOmega;
914 if (newXS2Pi < 0.)
915 newXS1Pi=oldXS1Pi-(xsEtaOm
916 else
917 newXS1Pi=oldXS1Pi;
918 }
919 else
920 newXS1Pi=oldXS1Pi-xsEtaOmega;
921 return newXS1Pi;
922 }
923 else if (xpi == 2) {
924 if (oldXS4Pi != 0.)
925 newXS2Pi=oldXS2Pi;
926 else if (oldXS3Pi != 0.) {
927 newXS3Pi=oldXS3Pi-xsEtaOmega;
928 if (newXS3Pi < 0.)
929 newXS2Pi=oldXS2Pi-(xsEtaOm
930 else
931 newXS2Pi=oldXS2Pi;
932 }
933 else {
934 newXS2Pi=oldXS2Pi-xsEtaOmega;
935 if (newXS2Pi < 0.)
936 newXS2Pi=0.;
937 }
938 return newXS2Pi;
939 }
940 else if (xpi == 3) {
941 if (oldXS4Pi != 0.) {
942 newXS4Pi=oldXS4Pi-xsEtaOmega;
943 if (newXS4Pi < 0.)
944 newXS3Pi=oldXS3Pi-(xsEtaOm
945 else
946 newXS3Pi=oldXS3Pi;
947 }
948 else {
949 newXS3Pi=oldXS3Pi-xsEtaOmega;
950 if (newXS3Pi < 0.)
951 newXS3Pi=0.;
952 }
953 return newXS3Pi;
954 }
955 else if (xpi == 4) {
956 newXS4Pi=oldXS4Pi-xsEtaOmega;
957 if (newXS4Pi < 0.)
958 newXS4Pi=0.;
959 return newXS4Pi;
960 }
961
962 else // should never reach this point
963 return 0.;
964 }
965
966
967 G4double CrossSectionsMultiPionsAndResonan
968 // Cross section for nucleon-nucleon p
969
970 const G4int iso=ParticleTable::getIsos
971 if (iso!=0)
972 return 0.;
973
974 const G4double ener=KinematicsUtils::t
975 if (ener < 2018.563) return 0.;
976
977 const G4double xsiso2=CrossSectionsMul
978 const G4double xsiso0=CrossSectionsMul
979
980 return 0.25*(CrossSectionsMultiPions::
981 }
982
983 G4double CrossSectionsMultiPionsAndResonan
984 const G4double ener=KinematicsUtils::t
985 if (ener < 2018.563) return 0.;
986 const G4int iso=ParticleTable::getIsos
987
988 const G4double xsiso2=CrossSectionsMul
989 if (iso != 0)
990 return CrossSectionsMultiPions::NN
991 else {
992 const G4double xsiso0=CrossSection
993 return 0.5*(CrossSectionsMultiPion
994 }
995 }
996
997 G4double CrossSectionsMultiPionsAndResonan
998 //
999 // Nucleon-Nucleon producing one e
1000 //
1001 const G4double ener=KinematicsUtils::
1002 if (ener < 2018.563) return 0.;
1003 const G4int iso=ParticleTable::getIso
1004
1005
1006 const G4double xsiso2=CrossSectionsMu
1007 if (iso != 0) {
1008 return CrossSectionsMultiPions::N
1009 }
1010 else {
1011 const G4double xsiso0=CrossSectio
1012 return 0.5*(CrossSectionsMultiPio
1013 }
1014 }
1015
1016 G4double CrossSectionsMultiPionsAndResona
1017 //
1018 // Nucleon-Nucleon producing one
1019 //
1020
1021 const G4double ener=KinematicsUtils::
1022 if (ener < 2018.563) return 0.;
1023 const G4int iso=ParticleTable::getIso
1024
1025
1026 const G4double xsiso2=CrossSectionsMu
1027 const G4double xs1pi2=CrossSectionsMu
1028 const G4double xs2pi2=CrossSectionsMu
1029 if (iso != 0)
1030 return CrossSectionsMultiPions::N
1031 else {
1032 const G4double xsiso0=CrossSectio
1033 const G4double xs1pi0=CrossSectio
1034 const G4double xs2pi0=CrossSectio
1035 return 0.5*(CrossSectionsMultiPio
1036 }
1037 }
1038
1039 G4double CrossSectionsMultiPionsAndResona
1040 //
1041 // Nucleon-Nucleon producing one
1042 //
1043
1044 const G4double ener=KinematicsUtils::
1045 if (ener < 2018.563) return 0.;
1046 const G4double s = ener*ener;
1047 const G4int i = ParticleTable::getIso
1048 G4double xsinelas;
1049 if (i!=0)
1050 xsinelas=CrossSectionsMultiPions:
1051 else
1052 xsinelas=0.5*(CrossSectionsMultiP
1053 if (xsinelas <= 1.e-9) return 0.;
1054 G4double ratio=(NNToNNEta(particle1,
1055 if(s<6.25E6)
1056 return 0.;
1057 const G4double sigma = NNToNNEta(part
1058 return ((sigma>1.e-9) ? sigma : 0.);
1059 }
1060
1061 G4double CrossSectionsMultiPionsAndResona
1062 //
1063 // Nucleon-Nucleon producing one
1064 //
1065 // assert(xpi>0 && xpi<=nMaxPiNN);
1066 // assert(particle1->isNucleon() && particle2
1067
1068 const G4double ener=KinematicsUtils::
1069 if (ener < 2018.563) return 0.;
1070 const G4int i = ParticleTable::getIso
1071 G4double xsinelas;
1072 if (i!=0)
1073 xsinelas=CrossSectionsMultiPions:
1074 else
1075 xsinelas=0.5*(CrossSectionsMultiP
1076 if (xsinelas <= 1.e-9) return 0.;
1077 G4double ratio=(NNToNNEta(particle1,
1078
1079 if (xpi == 1)
1080 return NNToNNEtaOnePi(particle1,
1081 else if (xpi == 2)
1082 return NNToNNEtaTwoPi(particle1,
1083 else if (xpi == 3)
1084 return NNToNNEtaThreePi(particle1
1085 else if (xpi == 4)
1086 return NNToNNEtaFourPi(particle1,
1087 else // should never reach this point
1088 return 0.;
1089 }
1090
1091
1092 G4double CrossSectionsMultiPionsAndResona
1093 // assert(p1->isNucleon() && p2->isNucleon())
1094 const G4int i = ParticleTable::getIso
1095 const G4double ener=KinematicsUtils::
1096 if (ener < 2018.563) return 0.;
1097 G4double xsinelas;
1098 if (i!=0)
1099 xsinelas=CrossSectionsMultiPions:
1100 else
1101 xsinelas=0.5*(CrossSectionsMultiP
1102 if (xsinelas <= 1.e-9) return 0.;
1103 G4double ratio=(NNToNNEta(p1, p2)-NNT
1104 G4double sigma = NNToNNEtaOnePiOrDelt
1105 if(i==0)
1106 sigma *= 0.5;
1107 return sigma;
1108 }
1109
1110
1111 G4double CrossSectionsMultiPionsAndResona
1112 // Cross section for nucleon-nucleon
1113
1114 const G4int iso=ParticleTable::getIso
1115 if (iso!=0)
1116 return 0.;
1117
1118 const G4double ener=KinematicsUtils::
1119 if (ener < 2018.563) return 0.;
1120
1121 const G4double xsiso2=CrossSectionsMu
1122 const G4double xsiso0=CrossSectionsMu
1123
1124 return 0.25*(CrossSectionsMultiPions:
1125 }
1126
1127 G4double CrossSectionsMultiPionsAndResona
1128 const G4double ener=KinematicsUtils::
1129 if (ener < 2018.563) return 0.;
1130 const G4int iso=ParticleTable::getIso
1131
1132 const G4double xsiso2=CrossSectionsMu
1133 if (iso != 0)
1134 return CrossSectionsMultiPions::N
1135 else {
1136 const G4double xsiso0=CrossSectio
1137 return 0.5*(CrossSectionsMultiPio
1138 }
1139 }
1140
1141 G4double CrossSectionsMultiPionsAndResona
1142 //
1143 // Nucleon-Nucleon producing one
1144 //
1145 const G4double ener=KinematicsUtils::
1146 if (ener < 2018.563) return 0.;
1147 const G4int iso=ParticleTable::getIso
1148
1149 const G4double xsiso2=CrossSectionsMu
1150 if (iso != 0) {
1151 return CrossSectionsMultiPions::N
1152 }
1153 else {
1154 const G4double xsiso0=CrossSectio
1155 return 0.5*(CrossSectionsMultiPio
1156 }
1157 }
1158
1159 G4double CrossSectionsMultiPionsAndResona
1160 //
1161 // Nucleon-Nucleon producing one
1162 //
1163
1164 const G4double ener=KinematicsUtils::
1165 if (ener < 2018.563) return 0.;
1166 const G4int iso=ParticleTable::getIso
1167
1168
1169 const G4double xsiso2=CrossSectionsMu
1170 const G4double xs1pi2=CrossSectionsMu
1171 const G4double xs2pi2=CrossSectionsMu
1172 if (iso != 0)
1173 return CrossSectionsMultiPions::N
1174 else {
1175 const G4double xsiso0=CrossSectio
1176 const G4double xs1pi0=CrossSectio
1177 const G4double xs2pi0=CrossSectio
1178 return 0.5*(CrossSectionsMultiPio
1179 }
1180 }
1181
1182 G4double CrossSectionsMultiPionsAndResona
1183 //
1184 // Nucleon-Nucleon producing one
1185 //
1186 //jcd to be removed
1187 // return 0.;
1188 //jcd
1189
1190 const G4double ener=KinematicsUtils::
1191 if (ener < 2018.563) return 0.;
1192 const G4double s = ener*ener;
1193 const G4int i = ParticleTable::getIso
1194 G4double xsinelas;
1195 if (i!=0)
1196 xsinelas=CrossSectionsMultiPions:
1197 else
1198 xsinelas=0.5*(CrossSectionsMultiP
1199 if (xsinelas <= 1.e-9) return 0.;
1200 G4double ratio=(NNToNNOmega(particle1
1201 if(s<6.25E6)
1202 return 0.;
1203 const G4double sigma = NNToNNOmega(pa
1204 return ((sigma>1.e-9) ? sigma : 0.);
1205 }
1206
1207 G4double CrossSectionsMultiPionsAndResona
1208 //
1209 // Nucleon-Nucleon producing one
1210 //
1211 // assert(xpi>0 && xpi<=nMaxPiNN);
1212 // assert(particle1->isNucleon() && particle2
1213 //jcd to be removed
1214 // return 0.;
1215 //jcd
1216
1217 const G4double ener=KinematicsUtils::
1218 if (ener < 2018.563) return 0.;
1219 const G4int i = ParticleTable::getIso
1220 G4double xsinelas;
1221 if (i!=0)
1222 xsinelas=CrossSectionsMultiPions:
1223 else
1224 xsinelas=0.5*(CrossSectionsMultiP
1225 if (xsinelas <= 1.e-9) return 0.;
1226 G4double ratio=(NNToNNOmega(particle1
1227
1228 if (xpi == 1)
1229 return NNToNNOmegaOnePi(particle1
1230 else if (xpi == 2)
1231 return NNToNNOmegaTwoPi(particle1
1232 else if (xpi == 3)
1233 return NNToNNOmegaThreePi(particl
1234 else if (xpi == 4)
1235 return NNToNNOmegaFourPi(particle
1236 else // should never reach this point
1237 return 0.;
1238 }
1239
1240
1241 G4double CrossSectionsMultiPionsAndResona
1242 // assert(p1->isNucleon() && p2->isNucleon())
1243 //jcd to be removed
1244 // return 0.;
1245 //jcd
1246 const G4int i = ParticleTable::getIso
1247 const G4double ener=KinematicsUtils::
1248 if (ener < 2018.563) return 0.;
1249 G4double xsinelas;
1250 if (i!=0)
1251 xsinelas=CrossSectionsMultiPions:
1252 else
1253 xsinelas=0.5*(CrossSectionsMultiP
1254 if (xsinelas <= 1.e-9) return 0.;
1255 G4double ratio=(NNToNNOmega(p1, p2)-N
1256 G4double sigma = NNToNNOmegaOnePiOrDe
1257 if(i==0)
1258 sigma *= 0.5;
1259 return sigma;
1260 }
1261
1262
1263 } // namespace G4INCL
1264
1265