ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCrossSectionsMultiPions.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCrossSectionsMultiPions.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 "G4INCLCrossSectionsMultiPions.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLParticleTable.hh"
41 #include "G4INCLLogger.hh"
42 // #include <cassert>
43 
44 namespace G4INCL {
45 
46  template<G4int N>
47  struct BystrickyEvaluator {
48  static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients<N> const &coeffs) {
49  const G4double pMeV = pLab*1E3;
50  const G4double ekin=std::sqrt(ParticleTable::effectiveNucleonMass2+pMeV*pMeV)-ParticleTable::effectiveNucleonMass;
51  const G4double xrat=ekin*oneOverThreshold;
52  const G4double x=std::log(xrat);
53  return HornerEvaluator<N>::eval(x, coeffs) * x * std::exp(-0.5*x);
54  }
55  };
56 
57  const G4int CrossSectionsMultiPions::nMaxPiNN = 4;
58  const G4int CrossSectionsMultiPions::nMaxPiPiN = 4;
59 
60  const G4double CrossSectionsMultiPions::s11pzOOT = 0.0035761542037692665889;
61  const G4double CrossSectionsMultiPions::s01ppOOT = 0.003421025623481919853;
62  const G4double CrossSectionsMultiPions::s01pzOOT = 0.0035739814152966403123;
63  const G4double CrossSectionsMultiPions::s11pmOOT = 0.0034855350296270480281;
64  const G4double CrossSectionsMultiPions::s12pmOOT = 0.0016672224074691565119;
65  const G4double CrossSectionsMultiPions::s12ppOOT = 0.0016507643038726931312;
66  const G4double CrossSectionsMultiPions::s12zzOOT = 0.0011111111111111111111;
67  const G4double CrossSectionsMultiPions::s02pzOOT = 0.00125;
68  const G4double CrossSectionsMultiPions::s02pmOOT = 0.0016661112962345883443;
69  const G4double CrossSectionsMultiPions::s12mzOOT = 0.0017047391749062392793;
70 
71  CrossSectionsMultiPions::CrossSectionsMultiPions() :
72  s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
73  s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
74  s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
75  s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
76  s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
77  s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
78  s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
79  s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
80  s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
81  s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
82  {
83  }
84 
85  G4double CrossSectionsMultiPions::NNElastic(Particle const * const part1, Particle const * const part2) {
86 
87  /* The NN cross section is parametrised as a function of the lab momentum
88  * of one of the nucleons. For NDelta or DeltaDelta, the physical
89  * assumption is that the cross section is the same as NN *for the same
90  * total CM energy*. Thus, we calculate s from the particles involved, and
91  * we convert this value to the lab momentum of a nucleon *as if this were
92  * an NN collision*.
93  */
94  const G4double s = KinematicsUtils::squareTotalEnergyInCM(part1, part2);
95 
96  if(part1->isNucleon() && part2->isNucleon()) { // NN
97  const G4int i = ParticleTable::getIsospin(part1->getType())
98  + ParticleTable::getIsospin(part2->getType());
99  return NNElasticFixed(s, i);
100  }
101  else { // Nucleon-Delta and Delta-Delta
102  const G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
103  if (plab < 0.440) {
104  return 34.*std::pow(plab/0.4, (-2.104));
105  }
106  else if (plab < 0.800) {
107  return 23.5+1000.*std::pow(plab-0.7, 4);
108  }
109  else if (plab <= 2.0) {
110  return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
111  }
112  else {
113  return 77./(plab+1.5);
114  }
115  }
116  }
117 
118  G4double CrossSectionsMultiPions::NNElasticFixed(const G4double s, const G4int i) {
119 
120  /* From NNElastic, with isospin fixed and for NN only.
121  */
122 
123  G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
124  G4double sigma = 0.;
125 
126  if (i == 0) { // pn
127  if (plab < 0.446) {
128  G4double alp=std::log(plab);
129  sigma = 6.3555*std::exp(-3.2481*alp-0.377*alp*alp);
130  }
131  else if (plab < 0.851) {
132  sigma = 33.+196.*std::pow(std::fabs(plab-0.95),2.5);
133  }
134  else if (plab <= 2.0) {
135  sigma = 31./std::sqrt(plab);
136  }
137  else {
138  sigma = 77./(plab+1.5);
139  }
140  //if(plab < 0.9 && plab > 0.802) sigma -= 0.1387*std::exp(-std::pow((plab-0.861),2)/0.0006861); //correction if totalcx-sumcx < 0.1
141  //if(plab < 1.4 && plab > 1.31) sigma -= 0.1088*std::exp(-std::pow((plab-1.35),2)/0.00141); //correction if totalcx-sumcx < 0.1
142  return sigma;
143  }
144  else { // pp and nn
145  if (plab < 0.440) {
146  return 34.*std::pow(plab/0.4, (-2.104));
147  }
148  else if (plab < 0.8067) {
149  return 23.5+1000.*std::pow(plab-0.7, 4);
150  }
151  else if (plab <= 2.0) {
152  return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
153  }
154  else if (plab <= 3.0956) {
155  return 77./(plab+1.5);
156  }
157  else {
158  G4double alp=std::log(plab);
159  return 11.2+25.5*std::pow(plab, -1.12)+0.151*std::pow(alp, 2)-1.62*alp;
160  }
161  }
162  }
163 
164  G4double CrossSectionsMultiPions::NNTot(Particle const * const part1, Particle const * const part2) {
165 
166  G4int i = ParticleTable::getIsospin(part1->getType())
167  + ParticleTable::getIsospin(part2->getType());
168 
169  if(part1->isNucleon() && part2->isNucleon()) { // NN
170  const G4double s = KinematicsUtils::squareTotalEnergyInCM(part1, part2);
171  return NNTotFixed(s, i);
172  }
173  else if (part1->isDelta() && part2->isDelta()) { // Delta-Delta
174  return elastic(part1, part2);
175  }
176  else { // Nucleon-Delta
177  return NDeltaToNN(part1, part2) + elastic(part1, part2);
178  }
179  }
180 
181  G4double CrossSectionsMultiPions::NNTotFixed(const G4double s, const G4int i) {
182 
183  /* From NNTot, with isospin fixed and for NN only.
184  */
185 
186  G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
187 
188  if (i == 0) { // pn
189  if (plab < 0.446) {
190  G4double alp=std::log(plab);
191  return 6.3555*std::exp(-3.2481*alp-0.377*std::pow(alp, 2));
192  }
193  else if (plab < 1.0) {
194  return 33.+196.*std::sqrt(std::pow(std::fabs(plab-0.95),5));
195  }
196  else if (plab < 1.924) {
197  return 24.2+8.9*plab;
198  }
199  else {
200  G4double alp=std::log(plab);
201  return 48.9-33.7*std::pow(plab, -3.08)+0.619*std::pow(alp, 2)-5.12*alp;
202  }
203  }
204  else { // pp and nn
205  if (plab < 0.440) {
206  return 34.*std::pow(plab/0.4, (-2.104));
207  }
208  else if (plab < 0.8734) {
209  return 23.5+1000.*std::pow(plab-0.7, 4);
210  }
211  else if (plab < 1.5) {
212  return 23.5+24.6/(1.+std::exp(-10.*(plab-1.2)));
213  }
214  else if (plab < 3.0044) {
215  return 41.+60.*(plab-0.9)*std::exp(-1.2*plab);
216  }
217  else {
218  G4double alp=std::log(plab);
219  return 45.6+219.*std::pow(plab, -4.23)+0.41*std::pow(alp, 2)-3.41*alp;
220  }
221  }
222  }
223 
224  G4double CrossSectionsMultiPions::NNInelasticIso(const G4double ener, const G4int iso) {
225 
226  const G4double s = ener*ener;
227  G4double sincl;
228 
229  if (iso != 0) {
230  if(s>=4074595.287720512986) { // plab>800 MeV/c
231  sincl = NNTotFixed(s, 2)-NNElasticFixed(s, 2);
232  }
233  else {
234  sincl = 0. ;
235  }
236  } else {
237  if(s>=4074595.287720512986) { // plab>800 MeV/c
238  sincl = 2*(NNTotFixed(s, 0)-NNElasticFixed(s, 0))-(NNTotFixed(s, 2)-NNElasticFixed(s, 2));
239  }
240  else {
241  return 0. ;
242  }
243  }
244  if (sincl < 0.) sincl = 0.;
245  return sincl;
246  }
247 
248  G4double CrossSectionsMultiPions::NNOnePiOrDelta(const G4double ener, const G4int iso, const G4double xsiso) {
249 
250  /* Article J. Physique 48 (1987)1901-1924 "Energy dependence of
251  nucleon-cucleon inelastic total cross-sections."
252  J. Bystricky, P. La France, F. Lehar, F. Perrot, T. Siemiarczuk & P. Winternitz
253  S11PZ= section pp->pp pi0
254  S01PP= section pp->pn pi+
255  S01PZ= section pn->pn pi0
256  S11PM= section pn->pp pi-
257  S= X-Section, 1st number : 1 if pp and 0 if pn
258  2nd number = number of pions, PP= pi+; PZ= pi0 ; PM= pi-
259  */
260 
261  const G4double s = ener*ener;
262  G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
263 
264  G4double snnpit1=0.;
265  G4double snnpit=0.;
266  G4double s11pz=0.;
267  G4double s01pp=0.;
268  G4double s01pz=0.;
269  G4double s11pm=0.;
270 
271  if ((iso != 0) && (plab < 2.1989)) {
272  snnpit = xsiso - NNTwoPi(ener, iso, xsiso);
273  if (snnpit < 1.e-8) snnpit=0.;
274  return snnpit;
275  }
276  else if ((iso == 0) && (plab < 1.7369)) {
277  snnpit = xsiso;
278  if (snnpit < 1.e-8) snnpit=0.;
279  return snnpit;
280  }
281 
282 //s11pz
283  if (plab > 18.) {
284  s11pz=55.185/std::pow((0.1412*plab+5),2);
285  }
286  else if (plab > 13.9) {
287  G4double alp=std::log(plab);
288  s11pz=6.67-13.3*std::pow(plab, -6.18)+0.456*alp*alp-3.29*alp;
289  }
290  else if (plab >= 0.7765) {
291  const G4double b=BystrickyEvaluator<7>::eval(plab,s11pzOOT,s11pzHC);
292  s11pz=b*b;
293  }
294 //s01pp
295  if (plab >= 0.79624) {
296  const G4double b=BystrickyEvaluator<8>::eval(plab,s01ppOOT,s01ppHC);
297  s01pp=b*b;
298  }
299 
300 // channel T=1
301  snnpit1=s11pz+s01pp;
302  if (snnpit1 < 1.e-8) snnpit1=0.;
303  if (iso != 0) {
304  return snnpit1;
305  }
306 
307 //s01pz
308  if (plab > 4.5) {
309  s01pz=15289.4/std::pow((11.573*plab+5),2);
310  }
311  else if (plab >= 0.777) {
312  const G4double b=BystrickyEvaluator<4>::eval(plab,s01pzOOT,s01pzHC);
313  s01pz=b*b;
314  }
315 //s11pm
316  if (plab > 14.) {
317  s11pm=46.68/std::pow((0.2231*plab+5),2);
318  }
319  else if (plab >= 0.788) {
320  const G4double b=BystrickyEvaluator<4>::eval(plab,s11pmOOT,s11pmHC);
321  s11pm=b*b;
322  }
323 
324 // channel T=0
325 // snnpit=s01pz+2*s11pm-snnpit1; //modif 2*(s01pz+2*s11pm)-snnpit1;
326  snnpit = 2*(s01pz+2*s11pm)-snnpit1;
327  if (snnpit < 1.e-8) snnpit=0.;
328  return snnpit;
329  }
330 
331  G4double CrossSectionsMultiPions::NNTwoPi(const G4double ener, const G4int iso, const G4double xsiso) {
332 
333  /* Article J. Physique 48 (1987)1901-1924 "Energy dependence of nucleon-cucleon inelastic total cross-sections."
334  J. Bystricky, P. La France, F. Lehar, F. Perrot, T. Siemiarczuk & P. Winternitz
335  S12PM : pp -> pp Pi+ Pi-
336  S12ZZ : pp -> pp Pi0 Pi0
337  S12PP : pp -> nn Pi+ Pi+
338  S02PZ : pp -> pn Pi+ Pi0
339  S02PM : pn -> pn Pi+ Pi-
340  S12MZ : pn -> pp Pi- Pi0
341  */
342 
343  const G4double s = ener*ener;
344  G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
345 
346  G4double snn2pit=0.;
347  G4double s12pm=0.;
348  G4double s12pp=0.;
349  G4double s12zz=0.;
350  G4double s02pz=0.;
351  G4double s02pm=0.;
352  G4double s12mz=0.;
353 
354  if (iso==0 && plab<3.33) {
355  snn2pit = xsiso - NNOnePiOrDelta(ener, iso, xsiso);
356  if (snn2pit < 1.e-8) snn2pit=0.;
357  return snn2pit;
358  }
359 
360  if (iso != 0) {
361 //s12pm
362  if (plab > 15.) {
363  s12pm=25.977/plab;
364  }
365  else if (plab >= 1.3817) {
366  const G4double b=BystrickyEvaluator<5>::eval(plab,s12pmOOT,s12pmHC);
367  s12pm=b*b;
368  }
369 //s12pp
370  if (plab > 10.) {
371  s12pp=141.505/std::pow((-0.1016*plab-7),2);
372  }
373  else if (plab >= 1.5739) {
374  const G4double b=BystrickyEvaluator<3>::eval(plab,s12ppOOT,s12ppHC);
375  s12pp=b*b;
376  }
377  }
378 //s12zz
379  if (plab > 4.) {
380  s12zz=97.355/std::pow((1.1579*plab+5),2);
381  }
382  else if (plab >= 1.72207) {
383  const G4double b=BystrickyEvaluator<4>::eval(plab,s12zzOOT,s12zzHC);
384  s12zz=b*b;
385  }
386 //s02pz
387  if (plab > 4.5) {
388  s02pz=178.082/std::pow((0.2014*plab+5),2);
389  }
390  else if (plab >= 1.5656) {
391  const G4double b=BystrickyEvaluator<4>::eval(plab,s02pzOOT,s02pzHC);
392  s02pz=b*b;
393  }
394 
395 // channel T=1
396  if (iso != 0) {
397  snn2pit=s12pm+s12pp+s12zz+s02pz;
398  if (snn2pit < 1.e-8) snn2pit=0.;
399  return snn2pit;
400  }
401 
402 //s02pm
403  if (plab > 5.) {
404  s02pm=135.826/std::pow(plab,2);
405  }
406  else if (plab >= 1.21925) {
407  const G4double b=BystrickyEvaluator<6>::eval(plab,s02pmOOT,s02pmHC);
408  s02pm=b*b;
409  }
410 //s12mz
411  if (plab >= 1.29269) {
412  const G4double b=BystrickyEvaluator<4>::eval(plab,s12mzOOT,s12mzHC);
413  s12mz=b*b;
414  }
415 
416 // channel T=0
417 // snn2pit=3*(0.5*s02pm+0.5*s12mz-0.5*s02pz-s12zz); //modif snn2pit=3*(s02pm+0.5*s12mz-0.5*s02pz-s12zz);
418  snn2pit=3*(s02pm+0.5*s12mz-0.5*s02pz-s12zz);
419  if (snn2pit < 1.e-8) snn2pit=0.;
420  return snn2pit;
421  }
422 
423  G4double CrossSectionsMultiPions::NNThreePi(const G4double ener, const G4int iso, const G4double xsiso, const G4double xs1pi, const G4double xs2pi) {
424 
425  const G4double s = ener*ener;
426  G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
427 
428  G4double snn3pit=0.;
429 
430  if (iso == 0) {
431 // channel T=0
432  if (plab > 7.2355) {
433  return 46.72/std::pow((plab - 5.8821),2);
434  }
435  else {
436  snn3pit=xsiso-xs1pi-xs2pi;
437  if (snn3pit < 1.e-8) snn3pit=0.;
438  return snn3pit;
439  }
440  }
441  else {
442 // channel T=1
443  if (plab > 7.206) {
444  return 5592.92/std::pow((plab+14.9764),2);
445  }
446  else if (plab > 2.1989){
447  snn3pit=xsiso-xs1pi-xs2pi;
448  if (snn3pit < 1.e-8) snn3pit=0.;
449  return snn3pit;
450  }
451  else return snn3pit;
452  }
453  }
454 
455  G4double CrossSectionsMultiPions::NNOnePi(Particle const * const particle1, Particle const * const particle2) {
456  // Cross section for nucleon-nucleon directly producing one pion
457 
458  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
459  if (iso!=0) // If pp or nn we choose to always pass by the N-N to N-Delta channel
460  return 0.;
461 
462  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
463 
464  const G4double xsiso2=NNInelasticIso(ener, 2);
465  const G4double xsiso0=NNInelasticIso(ener, 0);
466  return 0.25*(NNOnePiOrDelta(ener, 0, xsiso0)+ NNOnePiOrDelta(ener, 2, xsiso2));
467  }
468 
469  G4double CrossSectionsMultiPions::NNOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
470  // Cross section for nucleon-nucleon directly producing one pion or producing a nucleon-delta pair
471  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
472  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
473 
474  const G4double xsiso2=NNInelasticIso(ener, 2);
475  if (iso != 0)
476  return NNOnePiOrDelta(ener, iso, xsiso2);
477  else {
478  const G4double xsiso0=NNInelasticIso(ener, 0);
479  return 0.5*(NNOnePiOrDelta(ener, 0, xsiso0)+ NNOnePiOrDelta(ener, 2, xsiso2));
480  }
481  }
482 
483  G4double CrossSectionsMultiPions::NNTwoPi(Particle const * const particle1, Particle const * const particle2) {
484  //
485  // Nucleon-Nucleon producing one pion cross sections
486  //
487  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
488  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
489 
490 
491  const G4double xsiso2=NNInelasticIso(ener, 2);
492  if (iso != 0) {
493  return NNTwoPi(ener, 2, xsiso2);
494  }
495  else {
496  const G4double xsiso0=NNInelasticIso(ener, 0);
497  return 0.5*(NNTwoPi(ener, 0, xsiso0)+ NNTwoPi(ener, 2, xsiso2));
498  }
499  return 0.0; // Should never reach this point
500  }
501 
502  G4double CrossSectionsMultiPions::NNThreePi(Particle const * const particle1, Particle const * const particle2) {
503  //
504  // Nucleon-Nucleon producing one pion cross sections
505  //
506 
507  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
508  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
509 
510 
511  const G4double xsiso2=NNInelasticIso(ener, 2);
512  const G4double xs1pi2=NNOnePiOrDelta(ener, 2, xsiso2);
513  const G4double xs2pi2=NNTwoPi(ener, 2, xsiso2);
514  if (iso != 0)
515  return NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
516  else {
517  const G4double xsiso0=NNInelasticIso(ener, 0);
518  const G4double xs1pi0=NNOnePiOrDelta(ener, 0, xsiso0);
519  const G4double xs2pi0=NNTwoPi(ener, 0, xsiso0);
520  return 0.5*(NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
521  }
522  }
523 
524  G4double CrossSectionsMultiPions::NNFourPi(Particle const * const particle1, Particle const * const particle2) {
525  const G4double s = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
526  if(s<6.25E6)
527  return 0.;
528  const G4double sigma = NNTot(particle1, particle2) - NNElastic(particle1, particle2) - NNOnePiOrDelta(particle1, particle2) - NNTwoPi(particle1, particle2) - NNThreePi(particle1, particle2);
529  return ((sigma>1.e-9) ? sigma : 0.);
530  }
531 
532  G4double CrossSectionsMultiPions::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
533  //
534  // Nucleon-Nucleon producing xpi pions cross sections
535  //
536 // assert(xpi>0 && xpi<=nMaxPiNN);
537 // assert(particle1->isNucleon() && particle2->isNucleon());
538 
539  if (xpi == 1)
540  return NNOnePi(particle1, particle2);
541  else if (xpi == 2)
542  return NNTwoPi(particle1, particle2);
543  else if (xpi == 3)
544  return NNThreePi(particle1, particle2);
545  else if (xpi == 4)
546  return NNFourPi(particle1, particle2);
547  else // should never reach this point
548  return 0.;
549  }
550 
551 
552  G4double CrossSectionsMultiPions::spnPiPlusPHE(const G4double x) {
553  // HE and LE pi- p and pi+ n
554  G4double ramass = 0.0;
555 
556  if(x <= 1306.78) {
557  G4double y = x*x;
558  G4double q2;
559  q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
560  if (q2 > 0.) {
561  G4double q3=std::pow(q2, 3./2.);
562  G4double f3=q3/(q3+std::pow(180.0, 3));
563  G4double sdel;
564  sdel=326.5/(std::pow((x-1215.0-ramass)*2.0/110.0,2)+1.0);
565  return sdel*f3*(1.0-5.0*ramass/1215.0);
566  }
567  else {
568  return 0;
569  }
570  }
571  if(x <= 1754.0) {
572  return -2.33730e-06*std::pow(x, 3)+1.13819e-02*std::pow(x,2)
573  -1.83993e+01*x+9893.4;
574  } else if (x <= 2150.0) {
575  return 1.13531e-06*std::pow(x, 3)-6.91694e-03*std::pow(x, 2)
576  +1.39907e+01*x-9360.76;
577  } else {
578  return -3.18087*std::log(x)+52.9784;
579  }
580  }
581 
582  G4double CrossSectionsMultiPions::spnPiMinusPHE(const G4double x) {
583  // HE pi- p and pi+ n
584  G4double ramass = 0.0;
585 
586  if(x <= 1275.8) {
587  G4double y = x*x;
588  G4double q2;
589  q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
590  if (q2 > 0.) {
591  G4double q3=std::pow(q2, 3./2.);
592  G4double f3=q3/(q3+std::pow(180.0, 3));
593  G4double sdel;
594  sdel=326.5/(std::pow((x-1215.0-ramass)*2.0/110.0,2)+1.0);
595  return sdel*f3*(1.0-5.0*ramass/1215.0)/3.;
596  }
597  else {
598  return 0;
599  }
600  }
601  if(x <= 1495.0) {
602  return 0.00120683*(x-1372.52)*(x-1372.52)+26.2058;
603  } else if(x <= 1578.0) {
604  return 1.15873e-05*x*x+49965.6/((x-1519.59)*(x-1519.59)+2372.55);
605  } else if(x <= 2028.4) {
606  return 34.0248+43262.2/((x-1681.65)*(x-1681.65)+1689.35);
607  } else if(x <= 7500.0) {
608  return 3.3e-7*(x-7500.0)*(x-7500.0)+24.5;
609  } else {
610  return 24.5;
611  }
612  }
613 
614  G4double CrossSectionsMultiPions::total(Particle const * const p1, Particle const * const p2) {
615  G4double inelastic;
616  if(p1->isNucleon() && p2->isNucleon()) {
617  return NNTot(p1, p2);
618  } else if((p1->isNucleon() && p2->isDelta()) ||
619  (p1->isDelta() && p2->isNucleon())) {
620  inelastic = NDeltaToNN(p1, p2);
621  } else if((p1->isNucleon() && p2->isPion()) ||
622  (p1->isPion() && p2->isNucleon())) {
623  return piNTot(p1,p2);
624  } else {
625  inelastic = 0.;
626  }
627 
628  return inelastic + elastic(p1, p2);
629  }
630 
631 
632  G4double CrossSectionsMultiPions::piNIne(Particle const * const particle1, Particle const * const particle2) {
633  // piN inelastic cross section (Delta excluded)
634 
635  const Particle *pion;
636  const Particle *nucleon;
637  if(particle1->isNucleon()) {
638  nucleon = particle1;
639  pion = particle2;
640  } else {
641  pion = particle1;
642  nucleon = particle2;
643  }
644 // assert(pion->isPion());
645 
646  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
647 
648  // these limits correspond to sqrt(s)=1230 and 20000 MeV
649  if(pLab>212677. || pLab<296.367)
650  return 0.0;
651 
652  const G4int ipit3 = ParticleTable::getIsospin(pion->getType());
653  const G4int ind2t3 = ParticleTable::getIsospin(nucleon->getType());
654  const G4int cg = 4 + ind2t3*ipit3;
655 // assert(cg==2 || cg==4 || cg==6);
656 
657 // const G4double p1=1e-3*pLab;
658 // const G4double p2=std::log(p1);
659  G4double xpipp = 0.0;
660  G4double xpimp = 0.0;
661 
662  if(cg!=2) {
663  // x-section pi+ p inelastique :
664  xpipp=piPluspIne(pion,nucleon);
665 
666  if(cg==6) // cas pi+ p et pi- n
667  return xpipp;
668  }
669 
670  // x-section pi- p inelastique :
671  xpimp=piMinuspIne(pion,nucleon);
672 
673  if(cg==2) // cas pi- p et pi+ n
674  return xpimp;
675  else // cas pi0 p et pi0 n
676  return 0.5*(xpipp+xpimp);
677  }
678 
679  G4double CrossSectionsMultiPions::piNToDelta(Particle const * const particle1, Particle const * const particle2) {
680  // piN Delta production
681 
682  G4double x = KinematicsUtils::totalEnergyInCM(particle1, particle2);
683  if(x>20000.) return 0.0; // no cross section above this value
684 
685  G4int ipit3 = 0;
686  G4int ind2t3 = 0;
687  const G4double ramass = 0.0;
688 
689  if(particle1->isPion()) {
690  ipit3 = ParticleTable::getIsospin(particle1->getType());
691  ind2t3 = ParticleTable::getIsospin(particle2->getType());
692  } else if(particle2->isPion()) {
693  ipit3 = ParticleTable::getIsospin(particle2->getType());
694  ind2t3 = ParticleTable::getIsospin(particle1->getType());
695  }
696 
697  const G4double y=x*x;
698  const G4double q2=(y-1076.0*1076.0)*(y-800.0*800.0)/y/4.0;
699  if (q2 <= 0.) {
700  return 0.0;
701  }
702  const G4double q3 = std::pow(std::sqrt(q2),3);
703  const G4double f3 = q3/(q3 + 5832000.); // 5832000 = 180^3
704  G4double sdelResult = 326.5/(std::pow((x-1215.0-ramass)*2.0/(110.0-ramass), 2)+1.0);
705  sdelResult = sdelResult*(1.0-5.0*ramass/1215.0);
706  const G4int cg = 4 + ind2t3*ipit3;
707  sdelResult = sdelResult*f3*cg/6.0;
708 
709  return sdelResult;
710  }
711 
712  G4double CrossSectionsMultiPions::piNTot(Particle const * const particle1, Particle const * const particle2) {
713  // FUNCTION SPN(X,IND2T3,IPIT3,f17)
714  // SIGMA(PI+ + P) IN THE (3,3) REGION
715  // NEW FIT BY J.VANDERMEULEN + FIT BY Th AOUST ABOVE (3,3) RES
716  // CONST AT LOW AND VERY HIGH ENERGY
717  // COMMON/BL8/RATHR,RAMASS REL21800
718  // integer f17
719  // RATHR and RAMASS are always 0.0!!!
720 
721  G4double x = KinematicsUtils::totalEnergyInCM(particle1, particle2);
722 
723  G4int ipit3 = 0;
724  G4int ind2t3 = 0;
725 
726  if(particle1->isPion()) {
727  ipit3 = ParticleTable::getIsospin(particle1->getType());
728  ind2t3 = ParticleTable::getIsospin(particle2->getType());
729  } else if(particle2->isPion()) {
730  ipit3 = ParticleTable::getIsospin(particle2->getType());
731  ind2t3 = ParticleTable::getIsospin(particle1->getType());
732  }
733 
734  G4double spnResult=0.0;
735 
736  // HE pi+ p and pi- n
737  if((ind2t3 == 1 && ipit3 == 2) || (ind2t3 == -1 && ipit3 == -2))
738  spnResult=spnPiPlusPHE(x);
739  else if((ind2t3 == 1 && ipit3 == -2) || (ind2t3 == -1 && ipit3 == 2))
740  spnResult=spnPiMinusPHE(x);
741  else if(ipit3 == 0) spnResult = (spnPiPlusPHE(x) + spnPiMinusPHE(x))/2.0; // (spnpipphe(x)+spnpimphe(x))/2.0
742  else {
743  INCL_ERROR("Unknown configuration!\n" << particle1->print() << particle2->print() << '\n');
744  }
745 
746  return spnResult;
747  }
748 
749  G4double CrossSectionsMultiPions::NDeltaToNN(Particle const * const p1, Particle const * const p2) {
750  const G4int isospin = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
751  if(isospin==4 || isospin==-4) return 0.0;
752 
753  G4double s = KinematicsUtils::squareTotalEnergyInCM(p1, p2);
754  G4double Ecm = std::sqrt(s);
755  G4int deltaIsospin;
756  G4double deltaMass;
757  if(p1->isDelta()) {
758  deltaIsospin = ParticleTable::getIsospin(p1->getType());
759  deltaMass = p1->getMass();
760  } else {
761  deltaIsospin = ParticleTable::getIsospin(p2->getType());
762  deltaMass = p2->getMass();
763  }
764 
765  if(Ecm <= 938.3 + deltaMass) {
766  return 0.0;
767  }
768 
769  if(Ecm < 938.3 + deltaMass + 2.0) {
770  Ecm = 938.3 + deltaMass + 2.0;
771  s = Ecm*Ecm;
772  }
773 
774  const G4double x = (s - 4.*ParticleTable::effectiveNucleonMass2) /
775  (s - std::pow(ParticleTable::effectiveNucleonMass + deltaMass, 2));
776  const G4double y = s/(s - std::pow(deltaMass - ParticleTable::effectiveNucleonMass, 2));
777  /* Concerning the way we calculate the lab momentum, see the considerations
778  * in CrossSections::elasticNNLegacy().
779  */
780  G4double sDelta;
781  const G4double xsiso2=NNInelasticIso(Ecm, 2);
782  if (isospin != 0)
783  sDelta = NNOnePiOrDelta(Ecm, isospin, xsiso2);
784  else {
785  const G4double xsiso0=NNInelasticIso(Ecm, 0);
786  sDelta = 0.25*(NNOnePiOrDelta(Ecm, 0, xsiso0)+ NNOnePiOrDelta(Ecm, 2, xsiso2));
787  }
788  G4double result = 0.5 * x * y * sDelta;
789  /* modification for pion-induced cascade (see JC and MC LEMAIRE,NPA489(88)781
790  * result=3.*result
791  * pi absorption increased also for internal pions (7/3/01)
792  */
793  result *= 3.*(32.0 + isospin * isospin * (deltaIsospin * deltaIsospin - 5))/64.0;
794  result /= 1.0 + 0.25 * (isospin * isospin);
795  return result;
796  }
797 
798  G4double CrossSectionsMultiPions::NNToNDelta(Particle const * const p1, Particle const * const p2) {
799 // assert(p1->isNucleon() && p2->isNucleon());
800  const G4int isospin = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
801  G4double sigma = NNOnePiOrDelta(p1, p2);
802  if(isospin==0)
803  sigma *= 0.5;
804  return sigma;
805  }
806 
807  G4double CrossSectionsMultiPions::elastic(Particle const * const p1, Particle const * const p2) {
808 // if(!p1->isPion() && !p2->isPion()){
809  if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){
810  return NNElastic(p1, p2);
811  }
812 // else if (p1->isNucleon() || p2->isNucleon()){
813  else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
814  G4double pielas = piNTot(p1,p2) - piNIne(p1,p2) - piNToDelta(p1,p2);
815  if (pielas < 0.){
816  pielas = 0.;
817  }
818 // return piNTot(p1,p2) - piNIne(p1,p2) - piNToDelta(p1,p2);
819  return pielas;
820  }
821  else {
822  return 0.0;
823  }
824  }
825 
826  G4double CrossSectionsMultiPions::calculateNNAngularSlope(G4double pl, G4int iso) {
827  G4double x = 0.001 * pl; // Change to GeV
828  if(iso != 0) {
829  if(pl <= 2000.0) {
830  x = std::pow(x, 8);
831  return 5.5e-6 * x/(7.7 + x);
832  } else {
833  return (5.34 + 0.67*(x - 2.0)) * 1.0e-6;
834  }
835  } else {
836  if(pl < 800.0) {
837  G4double b = (7.16 - 1.63*x) * 1.0e-6;
838  return b/(1.0 + std::exp(-(x - 0.45)/0.05));
839  } else if(pl < 1100.0) {
840  return (9.87 - 4.88 * x) * 1.0e-6;
841  } else {
842  return (3.68 + 0.76*x) * 1.0e-6;
843  }
844  }
845  return 0.0; // Should never reach this point
846  }
847 
848 
849  G4double CrossSectionsMultiPions::piNToxPiN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
850  //
851  // pion-Nucleon producing xpi pions cross sections
852  //
853  const Particle *pion;
854  const Particle *nucleon;
855  if(particle1->isNucleon()) {
856  nucleon = particle1;
857  pion = particle2;
858  } else {
859  pion = particle1;
860  nucleon = particle2;
861  }
862 // assert(xpi>1 && xpi<=nMaxPiPiN);
863 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
864  const G4double plab = KinematicsUtils::momentumInLab(pion,nucleon);
865  if (xpi == 2) {
866  G4double OnePi=piNOnePi(particle1,particle2);
867  if (OnePi < 1.e-09) OnePi = 0.;
868  return OnePi;
869  }
870  else if (xpi == 3){
871  G4double TwoPi=piNTwoPi(particle1,particle2);
872  if (TwoPi < 1.e-09) TwoPi = 0.;
873  return TwoPi;
874  }
875  else if (xpi == 4) {
876  G4double piNThreePi = piNIne(particle1,particle2) - piNOnePi(particle1,particle2) - piNTwoPi(particle1,particle2);
877  if (piNThreePi < 1.e-09 || plab < 2000.) piNThreePi = 0.;
878  return piNThreePi;
879  } else // should never reach this point
880  return 0.0;
881  }
882 
883  G4double CrossSectionsMultiPions::piNOnePi(Particle const * const particle1, Particle const * const particle2) {
884  const Particle *pion;
885  const Particle *nucleon;
886  if(particle1->isNucleon()) {
887  nucleon = particle1;
888  pion = particle2;
889  } else {
890  pion = particle1;
891  nucleon = particle2;
892  }
893 // assert(pion->isPion());
894 
895  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
896 
897  // this limit corresponds to sqrt(s)=1230 MeV
898  if(pLab<296.367)
899  return 0.0;
900 
901  const G4int ipi = ParticleTable::getIsospin(pion->getType());
902  const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
903  const G4int cg = 4 + ind2*ipi;
904 // assert(cg==2 || cg==4 || cg==6);
905 
906  // const G4double p1=1e-3*pLab;
907  G4double tamp6=0.;
908  G4double tamp2=0.;
909  const G4double elas = elastic(particle1, particle2);
910 
911  // X-SECTION PI+ P INELASTIQUE :
912  if(cg != 2) {
913  tamp6=piPluspOnePi(particle1,particle2);
914  if (cg == 6){ // CAS PI+ P ET PI- N
915  if(tamp6 >= elas && pLab < 410.) tamp6 = elas;
916  return tamp6;
917  }
918  }
919 
920  // X-SECTION PI- P INELASTIQUE :
921  tamp2=piMinuspOnePi(particle1,particle2);
922  if (tamp2 < 0.0) tamp2=0;
923 
924  if (cg == 2) // CAS PI- P ET PI+ N
925  return tamp2;
926  else { // CAS PI0 P ET PI0 N
927  G4double s1pin = 0.5*(tamp6+tamp2);
928  const G4double inelastic = piNIne(particle1, particle2);
929  if(s1pin >= elas && pLab < 410.) s1pin = 0.;
930  if (s1pin > inelastic)
931  s1pin = inelastic;
932  return s1pin;
933  }
934  }
935 
936  G4double CrossSectionsMultiPions::piNTwoPi(Particle const * const particle1, Particle const * const particle2) {
937  //
938  // pion-nucleon interaction, producing 2 pions
939  // fit from Landolt-Bornstein multiplied by factor determined with evaluation of total xs
940  //
941 
942  const Particle *pion;
943  const Particle *nucleon;
944  if(particle1->isNucleon()) {
945  nucleon = particle1;
946  pion = particle2;
947  } else {
948  pion = particle1;
949  nucleon = particle2;
950  }
951 // assert(pion->isPion());
952 
953  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
954  const G4double elas = elastic(pion, nucleon);
955 
956  // this limit corresponds to sqrt(s)=1230 MeV
957  if(pLab<296.367)
958  return 0.0;
959 
960  const G4int ipi = ParticleTable::getIsospin(pion->getType());
961  const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
962  const G4int cg = 4 + ind2*ipi;
963 // assert(cg==2 || cg==4 || cg==6);
964 
965  G4double tamp6=0.;
966  G4double tamp2=0.;
967 
968  // X-SECTION PI+ P INELASTIQUE :
969  if(cg!=2) {
970  tamp6=piPluspTwoPi(particle1,particle2);
971  if(cg==6){ // CAS PI+ P ET PI- N
972  if(tamp6 >= elas && pLab < 410.) tamp6 = 0.;
973  return tamp6;}
974  }
975 
976  // X-SECTION PI- P INELASTIQUE :
977  tamp2=piMinuspTwoPi(particle1,particle2);
978 
979  if(cg==2) // CAS PI- P ET PI+ N
980  return tamp2;
981  else { // CAS PI0 P ET PI0 N
982  const G4double s2pin=0.5*(tamp6+tamp2);
983  return s2pin;
984  }
985  }
986 
987  G4double CrossSectionsMultiPions::piPluspIne(Particle const * const particle1, Particle const * const particle2) {
988  // piPlusP inelastic cross section (Delta excluded)
989 
990  const Particle *pion;
991  const Particle *nucleon;
992  if(particle1->isNucleon()) {
993  nucleon = particle1;
994  pion = particle2;
995  } else {
996  pion = particle1;
997  nucleon = particle2;
998  }
999 // assert(pion->isPion());
1000 
1001  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1002 
1003  // these limits correspond to sqrt(s)=1230 and 20000 MeV
1004  if(pLab>212677. || pLab<296.367)
1005  return 0.0;
1006 
1007 // const G4int ipit3 = ParticleTable::getIsospin(pion->getType());
1008 // const G4int ind2t3 = ParticleTable::getIsospin(nucleon->getType());
1009 // const G4int cg = 4 + ind2t3*ipit3;
1010 // assert(cg==2 || cg==4 || cg==6);
1011 
1012  const G4double p1=1e-3*pLab;
1013  const G4double p2=std::log(p1);
1014  G4double xpipp = 0.0;
1015 
1016  // x-section pi+ p inelastique :
1017  if(p1<=0.75)
1018  xpipp=17.965*std::pow(p1, 5.4606);
1019  else
1020  xpipp=24.3-12.3*std::pow(p1, -1.91)+0.324*p2*p2-2.44*p2;
1021  // cas pi+ p et pi- n
1022  return xpipp;
1023 
1024  }
1025 
1026  G4double CrossSectionsMultiPions::piMinuspIne(Particle const * const particle1, Particle const * const particle2) {
1027  // piMinusp inelastic cross section (Delta excluded)
1028 
1029  const Particle *pion;
1030  const Particle *nucleon;
1031  if(particle1->isNucleon()) {
1032  nucleon = particle1;
1033  pion = particle2;
1034  } else {
1035  pion = particle1;
1036  nucleon = particle2;
1037  }
1038 // assert(pion->isPion());
1039 
1040  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1041 
1042  // these limits correspond to sqrt(s)=1230 and 20000 MeV
1043  if(pLab>212677. || pLab<296.367)
1044  return 0.0;
1045 
1046 // const G4int ipit3 = ParticleTable::getIsospin(pion->getType());
1047 // const G4int ind2t3 = ParticleTable::getIsospin(nucleon->getType());
1048 // const G4int cg = 4 + ind2t3*ipit3;
1049 // assert(cg==2 || cg==4 || cg==6);
1050 
1051  const G4double p1=1e-3*pLab;
1052  const G4double p2=std::log(p1);
1053  G4double xpimp = 0.0;
1054 
1055  // x-section pi- p inelastique :
1056  if(p1 <= 0.4731)
1057  xpimp=0;
1058  else
1059  xpimp=26.6-7.18*std::pow(p1, -1.86)+0.327*p2*p2-2.81*p2;
1060  if(xpimp<0.)
1061  xpimp=0;
1062 
1063  // cas pi- p et pi+ n
1064  return xpimp;
1065 
1066  }
1067 
1068  G4double CrossSectionsMultiPions::piPluspOnePi(Particle const * const particle1, Particle const * const particle2) {
1069  const Particle *pion;
1070  const Particle *nucleon;
1071  if(particle1->isNucleon()) {
1072  nucleon = particle1;
1073  pion = particle2;
1074  } else {
1075  pion = particle1;
1076  nucleon = particle2;
1077  }
1078 // assert(pion->isPion());
1079 
1080  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1081 
1082  // this limit corresponds to sqrt(s)=1230 MeV
1083  if(pLab<296.367)
1084  return 0.0;
1085 
1086  // const G4int ipi = ParticleTable::getIsospin(pion->getType());
1087  // const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
1088  // const G4int cg = 4 + ind2*ipi;
1089  // assert(cg==2 || cg==4 || cg==6);
1090 
1091  const G4double p1=1e-3*pLab;
1092  G4double tamp6=0.;
1093 
1094  // X-SECTION PI+ P INELASTIQUE :
1095  if(pLab < 1532.52) // corresponds to sqrt(s)=1946 MeV
1096  tamp6=piPluspIne(particle1, particle2);
1097  else
1098  tamp6=0.204+18.2*std::pow(p1, -1.72)+6.33*std::pow(p1, -1.13);
1099 
1100  // CAS PI+ P ET PI- N
1101  return tamp6;
1102 
1103  }
1104 
1105  G4double CrossSectionsMultiPions::piMinuspOnePi(Particle const * const particle1, Particle const * const particle2) {
1106  const Particle *pion;
1107  const Particle *nucleon;
1108  if(particle1->isNucleon()) {
1109  nucleon = particle1;
1110  pion = particle2;
1111  } else {
1112  pion = particle1;
1113  nucleon = particle2;
1114  }
1115 // assert(pion->isPion());
1116 
1117  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1118 
1119  // this limit corresponds to sqrt(s)=1230 MeV
1120  if(pLab<296.367)
1121  return 0.0;
1122 
1123  // const G4int ipi = ParticleTable::getIsospin(pion->getType());
1124  // const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
1125  // const G4int cg = 4 + ind2*ipi;
1126  // assert(cg==2 || cg==4 || cg==6);
1127 
1128  const G4double p1=1e-3*pLab;
1129  G4double tamp2=0.;
1130 
1131  // X-SECTION PI- P INELASTIQUE :
1132  if (pLab < 1228.06) // corresponds to sqrt(s)=1794 MeV
1133  tamp2=piMinuspIne(particle1, particle2);
1134  else
1135  tamp2=9.04*std::pow(p1, -1.17)+18.*std::pow(p1, -1.21); // tamp2=9.04*std::pow(p1, -1.17)+(13.5*std::pow(p1, -1.21))*4./3.;
1136  if (tamp2 < 0.0) tamp2=0;
1137 
1138  // CAS PI- P ET PI+ N
1139  return tamp2;
1140  }
1141 
1142  G4double CrossSectionsMultiPions::piPluspTwoPi(Particle const * const particle1, Particle const * const particle2) {
1143  //
1144  // pion-nucleon interaction, producing 2 pions
1145  // fit from Landolt-Bornstein multiplied by factor determined with evaluation of total xs
1146  //
1147 
1148  const Particle *pion;
1149  const Particle *nucleon;
1150  if(particle1->isNucleon()) {
1151  nucleon = particle1;
1152  pion = particle2;
1153  } else {
1154  pion = particle1;
1155  nucleon = particle2;
1156  }
1157 // assert(pion->isPion());
1158 
1159  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1160 
1161  // this limit corresponds to sqrt(s)=1230 MeV
1162  if(pLab<296.367)
1163  return 0.0;
1164 
1165  // const G4int ipi = ParticleTable::getIsospin(pion->getType());
1166  // const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
1167  // const G4int cg = 4 + ind2*ipi;
1168  // assert(cg==2 || cg==4 || cg==6);
1169 
1170  const G4double p1=1e-3*pLab;
1171  G4double tamp6=0.;
1172 
1173  // X-SECTION PI+ P INELASTIQUE :
1174  if(pLab < 2444.7) // corresponds to sqrt(s)=2344 MeV
1175  tamp6=piPluspIne(particle1, particle2)-piPluspOnePi(particle1, particle2);
1176  else
1177  tamp6=1.59+25.5*std::pow(p1, -1.04); // tamp6=(0.636+10.2*std::pow(p1, -1.04))*15./6.;
1178 
1179  // CAS PI+ P ET PI- N
1180  return tamp6;
1181  }
1182 
1183  G4double CrossSectionsMultiPions::piMinuspTwoPi(Particle const * const particle1, Particle const * const particle2) {
1184  //
1185  // pion-nucleon interaction, producing 2 pions
1186  // fit from Landolt-Bornstein multiplied by factor determined with evaluation of total xs
1187  //
1188 
1189  const Particle *pion;
1190  const Particle *nucleon;
1191  if(particle1->isNucleon()) {
1192  nucleon = particle1;
1193  pion = particle2;
1194  } else {
1195  pion = particle1;
1196  nucleon = particle2;
1197  }
1198 // assert(pion->isPion());
1199 
1200  const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1201 
1202  // this limit corresponds to sqrt(s)=1230 MeV
1203  if(pLab<296.367)
1204  return 0.0;
1205 
1206  // const G4int ipi = ParticleTable::getIsospin(pion->getType());
1207  // const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
1208  // const G4int cg = 4 + ind2*ipi;
1209  // assert(cg==2 || cg==4 || cg==6);
1210 
1211  const G4double p1=1e-3*pLab;
1212  G4double tamp2=0.;
1213 
1214  // X-SECTION PI- P INELASTIQUE :
1215  if(pLab<2083.63) // corresponds to sqrt(s)=2195 MeV
1216  tamp2=piMinuspIne(particle1, particle2)-piMinuspOnePi(particle1, particle2);
1217  else
1218  tamp2=2.457794117647+18.066176470588*std::pow(p1, -0.92); // tamp2=(0.619+4.55*std::pow(p1, -0.92))*135./34.;
1219 
1220  // CAS PI- P ET PI+ N
1221  return tamp2;
1222 }
1223 
1224 
1225 
1226 
1227  G4double CrossSectionsMultiPions::piNToEtaN(Particle const * const, Particle const * const) {
1228  //
1229  // Pion-Nucleon producing Eta cross sections
1230  //
1231  return 0.;
1232  }
1233 
1234  G4double CrossSectionsMultiPions::piNToOmegaN(Particle const * const, Particle const * const) {
1235  //
1236  // Pion-Nucleon producing Omega cross sections
1237  //
1238  return 0.;
1239  }
1240 
1241  G4double CrossSectionsMultiPions::piNToEtaPrimeN(Particle const * const, Particle const * const) {
1242  //
1243  // Pion-Nucleon producing EtaPrime cross sections
1244  //
1245  return 0.;
1246  }
1247 
1248  G4double CrossSectionsMultiPions::etaNToPiN(Particle const * const, Particle const * const) {
1249  //
1250  // Eta-Nucleon producing Pion cross sections
1251  //
1252  return 0.;
1253  }
1254 
1255 
1256  G4double CrossSectionsMultiPions::etaNToPiPiN(Particle const * const, Particle const * const) {
1257  //
1258  // Eta-Nucleon producing Two Pions cross sections
1259  //
1260  return 0.;
1261  }
1262 
1263 
1264  G4double CrossSectionsMultiPions::omegaNToPiN(Particle const * const, Particle const * const) {
1265  //
1266  // Omega-Nucleon producing Pion cross sections
1267  //
1268  return 0.;
1269  }
1270 
1271  G4double CrossSectionsMultiPions::omegaNToPiPiN(Particle const * const, Particle const * const) {
1272  //
1273  // Omega-Nucleon producing Two Pions cross sections
1274  //
1275  return 0.;
1276  }
1277 
1278  G4double CrossSectionsMultiPions::etaPrimeNToPiN(Particle const * const, Particle const * const) {
1279  //
1280  // EtaPrime-Nucleon producing Pion cross sections
1281  //
1282  return 0.;
1283  }
1284 
1285  G4double CrossSectionsMultiPions::NNToNNEta(Particle const * const, Particle const * const) {
1286  //
1287  // Nucleon-Nucleon producing Eta cross sections
1288  //
1289  return 0.;
1290  }
1291 
1292  G4double CrossSectionsMultiPions::NNToNNEtaExclu(Particle const * const, Particle const * const) {
1293  //
1294  // Nucleon-Nucleon producing Eta cross sections
1295  //
1296  return 0.;
1297  }
1298 
1299  G4double CrossSectionsMultiPions::NNToNNEtaxPi(const G4int, Particle const * const, Particle const * const) {
1300  return 0.;
1301  }
1302 
1303  G4double CrossSectionsMultiPions::NNToNDeltaEta(Particle const * const, Particle const * const) {
1304  //
1305  // Nucleon-Nucleon producing N-Delta-Eta cross sections
1306  //
1307  return 0.;
1308  }
1309 
1310  G4double CrossSectionsMultiPions::NNToNNOmega(Particle const * const, Particle const * const) {
1311  //
1312  // Nucleon-Nucleon producing Omega cross sections
1313  //
1314  return 0.;
1315  }
1316 
1317  G4double CrossSectionsMultiPions::NNToNNOmegaExclu(Particle const * const, Particle const * const) {
1318  //
1319  // Nucleon-Nucleon producing Omega cross sections
1320  //
1321  return 0.;
1322  }
1323 
1324  G4double CrossSectionsMultiPions::NNToNNOmegaxPi(const G4int, Particle const * const, Particle const * const) {
1325  return 0.;
1326  }
1327 
1328  G4double CrossSectionsMultiPions::NNToNDeltaOmega(Particle const * const, Particle const * const) {
1329  //
1330  // Nucleon-Nucleon producing N-Delta-Omega cross sections
1331  //
1332  return 0.;
1333  }
1334 
1335 
1336 
1337 
1338  G4double CrossSectionsMultiPions::NYelastic(Particle const * const , Particle const * const ) {
1339  //
1340  // Hyperon-Nucleon elastic cross sections
1341  //
1342  return 0.;
1343  }
1344 
1345  G4double CrossSectionsMultiPions::NKelastic(Particle const * const , Particle const * const ) {
1346  //
1347  // Kaon-Nucleon elastic cross sections
1348  //
1349  return 0.;
1350  }
1351 
1352  G4double CrossSectionsMultiPions::NKbelastic(Particle const * const , Particle const * const ) {
1353  //
1354  // antiKaon-Nucleon elastic cross sections
1355  //
1356  return 0.;
1357  }
1358 
1359 
1360  G4double CrossSectionsMultiPions::NNToNLK(Particle const * const, Particle const * const) {
1361  //
1362  // Nucleon-Nucleon producing N-Lambda-Kaon cross sections
1363  //
1364  return 0.;
1365  }
1366 
1367  G4double CrossSectionsMultiPions::NNToNSK(Particle const * const, Particle const * const) {
1368  //
1369  // Nucleon-Nucleon producing N-Sigma-Kaon cross sections
1370  //
1371  return 0.;
1372  }
1373 
1374  G4double CrossSectionsMultiPions::NNToNLKpi(Particle const * const, Particle const * const) {
1375  //
1376  // Nucleon-Nucleon producing N-Lambda-Kaon-pion cross sections
1377  //
1378  return 0.;
1379  }
1380 
1381  G4double CrossSectionsMultiPions::NNToNSKpi(Particle const * const, Particle const * const) {
1382  //
1383  // Nucleon-Nucleon producing N-Sigma-Kaon-pion cross sections
1384  //
1385  return 0.;
1386  }
1387 
1388  G4double CrossSectionsMultiPions::NNToNLK2pi(Particle const * const, Particle const * const) {
1389  //
1390  // Nucleon-Nucleon producing N-Lambda-Kaon-2pion cross sections
1391  //
1392  return 0.;
1393  }
1394 
1395  G4double CrossSectionsMultiPions::NNToNSK2pi(Particle const * const, Particle const * const) {
1396  //
1397  // Nucleon-Nucleon producing N-Sigma-Kaon-2pion cross sections
1398  //
1399  return 0.;
1400  }
1401 
1402  G4double CrossSectionsMultiPions::NNToNNKKb(Particle const * const, Particle const * const) {
1403  //
1404  // Nucleon-Nucleon producing Nucleon-Nucleon-Kaon-antiKaon cross sections
1405  //
1406  return 0.;
1407  }
1408 
1409  G4double CrossSectionsMultiPions::NNToMissingStrangeness(Particle const * const, Particle const * const) {
1410  //
1411  // Nucleon-Nucleon missing strangeness production cross sections
1412  //
1413  return 0.;
1414  }
1415 
1416  G4double CrossSectionsMultiPions::NDeltaToNLK(Particle const * const, Particle const * const) {
1417  // Nucleon-Delta producing Nucleon Lambda Kaon cross section
1418  return 0;
1419  }
1420  G4double CrossSectionsMultiPions::NDeltaToNSK(Particle const * const, Particle const * const) {
1421  // Nucleon-Delta producing Nucleon Sigma Kaon cross section
1422  return 0;
1423  }
1424  G4double CrossSectionsMultiPions::NDeltaToDeltaLK(Particle const * const, Particle const * const) {
1425  // Nucleon-Delta producing Delta Lambda Kaon cross section
1426  return 0;
1427  }
1428  G4double CrossSectionsMultiPions::NDeltaToDeltaSK(Particle const * const, Particle const * const) {
1429  // Nucleon-Delta producing Delta Sigma Kaon cross section
1430  return 0;
1431  }
1432 
1433  G4double CrossSectionsMultiPions::NDeltaToNNKKb(Particle const * const, Particle const * const) {
1434  // Nucleon-Delta producing Nucleon-Nucleon Kaon antiKaon cross section
1435  return 0;
1436  }
1437 
1438 
1439  G4double CrossSectionsMultiPions::NpiToLK(Particle const * const, Particle const * const) {
1440  //
1441  // Pion-Nucleon producing Lambda-Kaon cross sections
1442  //
1443  return 0.;
1444  }
1445 
1446  G4double CrossSectionsMultiPions::NpiToSK(Particle const * const, Particle const * const) {
1447  //
1448  // Pion-Nucleon producing Sigma-Kaon cross sections
1449  //
1450  return 0.;
1451  }
1452  G4double CrossSectionsMultiPions::p_pimToSmKp(Particle const * const, Particle const * const) {
1453  return 0.;
1454  }
1455  G4double CrossSectionsMultiPions::p_pimToSzKz(Particle const * const, Particle const * const) {
1456  return 0.;
1457  }
1458  G4double CrossSectionsMultiPions::p_pizToSzKp(Particle const * const, Particle const * const) {
1459  return 0.;
1460  }
1461 
1462  G4double CrossSectionsMultiPions::NpiToLKpi(Particle const * const, Particle const * const) {
1463  //
1464  // Pion-Nucleon producing Lambda-Kaon-pion cross sections
1465  //
1466  return 0.;
1467  }
1468 
1469  G4double CrossSectionsMultiPions::NpiToSKpi(Particle const * const, Particle const * const) {
1470  //
1471  // Pion-Nucleon producing Sigma-Kaon-pion cross sections
1472  //
1473  return 0.;
1474  }
1475 
1476  G4double CrossSectionsMultiPions::NpiToLK2pi(Particle const * const, Particle const * const) {
1477  //
1478  // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1479  //
1480  return 0.;
1481  }
1482 
1483  G4double CrossSectionsMultiPions::NpiToSK2pi(Particle const * const, Particle const * const) {
1484  //
1485  // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1486  //
1487  return 0.;
1488  }
1489 
1490  G4double CrossSectionsMultiPions::NpiToNKKb(Particle const * const, Particle const * const) {
1491  //
1492  // Pion-Nucleon producing Nucleon-Kaon-antiKaon cross sections
1493  //
1494  return 0.;
1495  }
1496 
1497  G4double CrossSectionsMultiPions::NpiToMissingStrangeness(Particle const * const, Particle const * const) {
1498  //
1499  // Pion-Nucleon missing strangeness production cross sections
1500  //
1501  return 0.;
1502  }
1503 
1504  G4double CrossSectionsMultiPions::NLToNS(Particle const * const, Particle const * const) {
1505  //
1506  // Nucleon-Hyperon multiplet changing cross sections
1507  //
1508  return 0.;
1509  }
1510 
1511  G4double CrossSectionsMultiPions::NSToNL(Particle const * const, Particle const * const) {
1512  //
1513  // Nucleon-Sigma quasi-elastic cross sections
1514  //
1515  return 0.;
1516  }
1517 
1518  G4double CrossSectionsMultiPions::NSToNS(Particle const * const, Particle const * const) {
1519  //
1520  // Nucleon-Sigma quasi-elastic cross sections
1521  //
1522  return 0.;
1523  }
1524 
1525  G4double CrossSectionsMultiPions::NKToNK(Particle const * const, Particle const * const) {
1526  //
1527  // Nucleon-Kaon quasi-elastic cross sections
1528  //
1529  return 0.;
1530  }
1531 
1532  G4double CrossSectionsMultiPions::NKToNKpi(Particle const * const, Particle const * const) {
1533  //
1534  // Nucleon-Kaon producing Nucleon-Kaon-pion cross sections
1535  //
1536  return 0.;
1537  }
1538 
1539  G4double CrossSectionsMultiPions::NKToNK2pi(Particle const * const, Particle const * const) {
1540  //
1541  // Nucleon-Kaon producing Nucleon-Kaon-2pion cross sections
1542  //
1543  return 0.;
1544  }
1545 
1546  G4double CrossSectionsMultiPions::NKbToNKb(Particle const * const, Particle const * const) {
1547  //
1548  // Nucleon-antiKaon quasi-elastic cross sections
1549  //
1550  return 0.;
1551  }
1552 
1553  G4double CrossSectionsMultiPions::NKbToSpi(Particle const * const, Particle const * const) {
1554  //
1555  // Nucleon-antiKaon producing Sigma-pion cross sections
1556  //
1557  return 0.;
1558  }
1559 
1560  G4double CrossSectionsMultiPions::NKbToLpi(Particle const * const, Particle const * const) {
1561  //
1562  // Nucleon-antiKaon producing Lambda-pion cross sections
1563  //
1564  return 0.;
1565  }
1566 
1567  G4double CrossSectionsMultiPions::NKbToS2pi(Particle const * const, Particle const * const) {
1568  //
1569  // Nucleon-antiKaon producing Sigma-2pion cross sections
1570  //
1571  return 0.;
1572  }
1573 
1574  G4double CrossSectionsMultiPions::NKbToL2pi(Particle const * const, Particle const * const) {
1575  //
1576  // Nucleon-antiKaon producing Lambda-2pion cross sections
1577  //
1578  return 0.;
1579  }
1580 
1581  G4double CrossSectionsMultiPions::NKbToNKbpi(Particle const * const, Particle const * const) {
1582  //
1583  // Nucleon-antiKaon producing Nucleon-antiKaon-pion cross sections
1584  //
1585  return 0.;
1586  }
1587 
1588  G4double CrossSectionsMultiPions::NKbToNKb2pi(Particle const * const, Particle const * const) {
1589  //
1590  // Nucleon-antiKaon producing Nucleon-antiKaon-2pion cross sections
1591  //
1592  return 0.;
1593  }
1594 
1595 
1596 
1597 
1598 } // namespace G4INCL
1599