| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 #include "G4INCLCrossSectionsMultiPions.hh" 38 #include "G4INCLCrossSectionsMultiPions.hh"
39 #include "G4INCLKinematicsUtils.hh" 39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLParticleTable.hh" 40 #include "G4INCLParticleTable.hh"
41 #include "G4INCLLogger.hh" 41 #include "G4INCLLogger.hh"
42 // #include <cassert> 42 // #include <cassert>
43 43
44 namespace G4INCL { 44 namespace G4INCL {
45 45
46 template<G4int N> 46 template<G4int N>
47 struct BystrickyEvaluator { 47 struct BystrickyEvaluator {
48 static G4double eval(const G4double pLab 48 static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients<N> const &coeffs) {
49 const G4double pMeV = pLab*1E3; 49 const G4double pMeV = pLab*1E3;
50 const G4double ekin=std::sqrt(Particle 50 const G4double ekin=std::sqrt(ParticleTable::effectiveNucleonMass2+pMeV*pMeV)-ParticleTable::effectiveNucleonMass;
51 const G4double xrat=ekin*oneOverThresh 51 const G4double xrat=ekin*oneOverThreshold;
52 const G4double x=std::log(xrat); 52 const G4double x=std::log(xrat);
53 return HornerEvaluator<N>::eval(x, coe 53 return HornerEvaluator<N>::eval(x, coeffs) * x * std::exp(-0.5*x);
54 } 54 }
55 }; 55 };
56 56
57 const G4int CrossSectionsMultiPions::nMaxPiN <<
58 const G4int CrossSectionsMultiPions::nMaxPiP <<
59 <<
60 const G4double CrossSectionsMultiPions::s11p 57 const G4double CrossSectionsMultiPions::s11pzOOT = 0.0035761542037692665889;
61 const G4double CrossSectionsMultiPions::s01p 58 const G4double CrossSectionsMultiPions::s01ppOOT = 0.003421025623481919853;
62 const G4double CrossSectionsMultiPions::s01p 59 const G4double CrossSectionsMultiPions::s01pzOOT = 0.0035739814152966403123;
63 const G4double CrossSectionsMultiPions::s11p 60 const G4double CrossSectionsMultiPions::s11pmOOT = 0.0034855350296270480281;
64 const G4double CrossSectionsMultiPions::s12p 61 const G4double CrossSectionsMultiPions::s12pmOOT = 0.0016672224074691565119;
65 const G4double CrossSectionsMultiPions::s12p 62 const G4double CrossSectionsMultiPions::s12ppOOT = 0.0016507643038726931312;
66 const G4double CrossSectionsMultiPions::s12z 63 const G4double CrossSectionsMultiPions::s12zzOOT = 0.0011111111111111111111;
67 const G4double CrossSectionsMultiPions::s02p 64 const G4double CrossSectionsMultiPions::s02pzOOT = 0.00125;
68 const G4double CrossSectionsMultiPions::s02p 65 const G4double CrossSectionsMultiPions::s02pmOOT = 0.0016661112962345883443;
69 const G4double CrossSectionsMultiPions::s12m 66 const G4double CrossSectionsMultiPions::s12mzOOT = 0.0017047391749062392793;
70 67
71 CrossSectionsMultiPions::CrossSectionsMultiP 68 CrossSectionsMultiPions::CrossSectionsMultiPions() :
72 s11pzHC(-2.228000000000294018,8.7560000000 69 s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
73 s01ppHC(2.0570000000126518344,-6.029000000 70 s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
74 s01pzHC(0.18030000000000441851,7.870099999 71 s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
75 s11pmHC(0.20590000000000031866,3.345099999 72 s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
76 s12pmHC(-0.77235999999999901328,4.26265999 73 s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
77 s12ppHC(-0.75724999999999975664,2.09343999 74 s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
78 s12zzHC(-0.89599999999996965072,7.88299999 75 s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
79 s02pzHC(-1.0579999999999967036,11.11399999 76 s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
80 s02pmHC(2.4009000000012553286,-7.768000000 77 s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
81 s12mzHC(-0.21858699999999976269,1.91489999 78 s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
82 { << 79 {
83 } << 80 }
84 81
85 G4double CrossSectionsMultiPions::NNElastic( 82 G4double CrossSectionsMultiPions::NNElastic(Particle const * const part1, Particle const * const part2) {
86 83
87 /* The NN cross section is parametrised as 84 /* The NN cross section is parametrised as a function of the lab momentum
88 * of one of the nucleons. For NDelta or D 85 * of one of the nucleons. For NDelta or DeltaDelta, the physical
89 * assumption is that the cross section is 86 * assumption is that the cross section is the same as NN *for the same
90 * total CM energy*. Thus, we calculate s 87 * total CM energy*. Thus, we calculate s from the particles involved, and
91 * we convert this value to the lab moment 88 * we convert this value to the lab momentum of a nucleon *as if this were
92 * an NN collision*. 89 * an NN collision*.
93 */ 90 */
94 const G4double s = KinematicsUtils::square 91 const G4double s = KinematicsUtils::squareTotalEnergyInCM(part1, part2);
95 92
96 if(part1->isNucleon() && part2->isNucleon( 93 if(part1->isNucleon() && part2->isNucleon()) { // NN
97 const G4int i = ParticleTable::getIsospi 94 const G4int i = ParticleTable::getIsospin(part1->getType())
98 + ParticleTable::getIsospin(part2->get 95 + ParticleTable::getIsospin(part2->getType());
99 return NNElasticFixed(s, i); 96 return NNElasticFixed(s, i);
100 } 97 }
101 else { // Nucleon-Delta and Delta-Delta 98 else { // Nucleon-Delta and Delta-Delta
102 const G4double plab = 0.001*KinematicsUt 99 const G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
103 if (plab < 0.440) { 100 if (plab < 0.440) {
104 return 34.*std::pow(plab/0.4, (-2.104) 101 return 34.*std::pow(plab/0.4, (-2.104));
105 } 102 }
106 else if (plab < 0.800) { 103 else if (plab < 0.800) {
107 return 23.5+1000.*std::pow(plab-0.7, 4 104 return 23.5+1000.*std::pow(plab-0.7, 4);
108 } 105 }
109 else if (plab <= 2.0) { 106 else if (plab <= 2.0) {
110 return 1250./(50.+plab)-4.*std::pow(pl 107 return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
111 } 108 }
112 else { 109 else {
113 return 77./(plab+1.5); 110 return 77./(plab+1.5);
114 } 111 }
115 } 112 }
116 } 113 }
117 114
118 G4double CrossSectionsMultiPions::NNElasti 115 G4double CrossSectionsMultiPions::NNElasticFixed(const G4double s, const G4int i) {
119 116
120 /* From NNElastic, with isospin fixed an 117 /* From NNElastic, with isospin fixed and for NN only.
121 */ 118 */
122 119
123 G4double plab = 0.001*KinematicsUtils::m 120 G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
124 G4double sigma = 0.; <<
125 121
126 if (i == 0) { // pn 122 if (i == 0) { // pn
127 if (plab < 0.446) { 123 if (plab < 0.446) {
128 G4double alp=std::log(plab); 124 G4double alp=std::log(plab);
129 sigma = 6.3555*std::exp(-3.2481*alp- << 125 return 6.3555*std::exp(-3.2481*alp-0.377*alp*alp);
130 } 126 }
131 else if (plab < 0.851) { 127 else if (plab < 0.851) {
132 sigma = 33.+196.*std::pow(std::fabs( << 128 return 33.+196.*std::pow(std::fabs(plab-0.95),2.5);
133 } 129 }
134 else if (plab <= 2.0) { 130 else if (plab <= 2.0) {
135 sigma = 31./std::sqrt(plab); << 131 return 31./std::sqrt(plab);
136 } 132 }
137 else { 133 else {
138 sigma = 77./(plab+1.5); << 134 return 77./(plab+1.5);
139 } 135 }
140 //if(plab < 0.9 && plab > 0.802) sigma <<
141 //if(plab < 1.4 && plab > 1.31) sigma <<
142 return sigma; <<
143 } 136 }
144 else { // pp and nn 137 else { // pp and nn
145 if (plab < 0.440) { 138 if (plab < 0.440) {
146 return 34.*std::pow(plab/0.4, (-2.10 139 return 34.*std::pow(plab/0.4, (-2.104));
147 } 140 }
148 else if (plab < 0.8067) { 141 else if (plab < 0.8067) {
149 return 23.5+1000.*std::pow(plab-0.7, 142 return 23.5+1000.*std::pow(plab-0.7, 4);
150 } 143 }
151 else if (plab <= 2.0) { 144 else if (plab <= 2.0) {
152 return 1250./(50.+plab)-4.*std::pow( 145 return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
153 } 146 }
154 else if (plab <= 3.0956) { 147 else if (plab <= 3.0956) {
155 return 77./(plab+1.5); 148 return 77./(plab+1.5);
156 } 149 }
157 else { 150 else {
158 G4double alp=std::log(plab); 151 G4double alp=std::log(plab);
159 return 11.2+25.5*std::pow(plab, -1.1 152 return 11.2+25.5*std::pow(plab, -1.12)+0.151*std::pow(alp, 2)-1.62*alp;
160 } 153 }
161 } 154 }
162 } 155 }
163 156
164 G4double CrossSectionsMultiPions::NNTot(Pa 157 G4double CrossSectionsMultiPions::NNTot(Particle const * const part1, Particle const * const part2) {
165 158
166 G4int i = ParticleTable::getIsospin(pa 159 G4int i = ParticleTable::getIsospin(part1->getType())
167 + ParticleTable::getIsospin(part2->get 160 + ParticleTable::getIsospin(part2->getType());
168 161
169 if(part1->isNucleon() && part2->isNucl 162 if(part1->isNucleon() && part2->isNucleon()) { // NN
170 const G4double s = KinematicsUtils:: 163 const G4double s = KinematicsUtils::squareTotalEnergyInCM(part1, part2);
171 return NNTotFixed(s, i); 164 return NNTotFixed(s, i);
172 } 165 }
173 else if (part1->isDelta() && part2->is 166 else if (part1->isDelta() && part2->isDelta()) { // Delta-Delta
174 return elastic(part1, part2); 167 return elastic(part1, part2);
175 } 168 }
176 else { // Nucleon-Delta 169 else { // Nucleon-Delta
177 return NDeltaToNN(part1, part2) + 170 return NDeltaToNN(part1, part2) + elastic(part1, part2);
178 } 171 }
179 } 172 }
180 173
181 G4double CrossSectionsMultiPions::NNTotFix 174 G4double CrossSectionsMultiPions::NNTotFixed(const G4double s, const G4int i) {
182 175
183 /* From NNTot, with isospin fixed and fo 176 /* From NNTot, with isospin fixed and for NN only.
184 */ 177 */
185 178
186 G4double plab = 0.001*KinematicsUtils::m 179 G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
187 180
188 if (i == 0) { // pn 181 if (i == 0) { // pn
189 if (plab < 0.446) { 182 if (plab < 0.446) {
190 G4double alp=std::log(plab); 183 G4double alp=std::log(plab);
191 return 6.3555*std::exp(-3.2481*alp-0 184 return 6.3555*std::exp(-3.2481*alp-0.377*std::pow(alp, 2));
192 } 185 }
193 else if (plab < 1.0) { 186 else if (plab < 1.0) {
194 return 33.+196.*std::sqrt(std::pow(s 187 return 33.+196.*std::sqrt(std::pow(std::fabs(plab-0.95),5));
195 } 188 }
196 else if (plab < 1.924) { 189 else if (plab < 1.924) {
197 return 24.2+8.9*plab; 190 return 24.2+8.9*plab;
198 } 191 }
199 else { 192 else {
200 G4double alp=std::log(plab); 193 G4double alp=std::log(plab);
201 return 48.9-33.7*std::pow(plab, -3.0 194 return 48.9-33.7*std::pow(plab, -3.08)+0.619*std::pow(alp, 2)-5.12*alp;
202 } 195 }
203 } 196 }
204 else { // pp and nn 197 else { // pp and nn
205 if (plab < 0.440) { 198 if (plab < 0.440) {
206 return 34.*std::pow(plab/0.4, (-2.10 199 return 34.*std::pow(plab/0.4, (-2.104));
207 } 200 }
208 else if (plab < 0.8734) { 201 else if (plab < 0.8734) {
209 return 23.5+1000.*std::pow(plab-0.7, 202 return 23.5+1000.*std::pow(plab-0.7, 4);
210 } 203 }
211 else if (plab < 1.5) { 204 else if (plab < 1.5) {
212 return 23.5+24.6/(1.+std::exp(-10.*( 205 return 23.5+24.6/(1.+std::exp(-10.*(plab-1.2)));
213 } 206 }
214 else if (plab < 3.0044) { 207 else if (plab < 3.0044) {
215 return 41.+60.*(plab-0.9)*std::exp(- 208 return 41.+60.*(plab-0.9)*std::exp(-1.2*plab);
216 } 209 }
217 else { 210 else {
218 G4double alp=std::log(plab); 211 G4double alp=std::log(plab);
219 return 45.6+219.*std::pow(plab, -4.2 212 return 45.6+219.*std::pow(plab, -4.23)+0.41*std::pow(alp, 2)-3.41*alp;
220 } 213 }
221 } 214 }
222 } 215 }
223 216
224 G4double CrossSectionsMultiPions::NNInelas 217 G4double CrossSectionsMultiPions::NNInelasticIso(const G4double ener, const G4int iso) {
225 218
226 const G4double s = ener*ener; 219 const G4double s = ener*ener;
227 G4double sincl; 220 G4double sincl;
228 221
229 if (iso != 0) { 222 if (iso != 0) {
230 if(s>=4074595.287720512986) { // plab> 223 if(s>=4074595.287720512986) { // plab>800 MeV/c
231 sincl = NNTotFixed(s, 2)-NNElasticFi 224 sincl = NNTotFixed(s, 2)-NNElasticFixed(s, 2);
232 } 225 }
233 else { 226 else {
234 sincl = 0. ; 227 sincl = 0. ;
235 } 228 }
236 } else { 229 } else {
237 if(s>=4074595.287720512986) { // plab> 230 if(s>=4074595.287720512986) { // plab>800 MeV/c
238 sincl = 2*(NNTotFixed(s, 0)-NNElasti 231 sincl = 2*(NNTotFixed(s, 0)-NNElasticFixed(s, 0))-(NNTotFixed(s, 2)-NNElasticFixed(s, 2));
239 } 232 }
240 else { 233 else {
241 return 0. ; 234 return 0. ;
242 } 235 }
243 } 236 }
244 if (sincl < 0.) sincl = 0.; 237 if (sincl < 0.) sincl = 0.;
245 return sincl; 238 return sincl;
246 } 239 }
247 240
248 G4double CrossSectionsMultiPions::NNOnePiO 241 G4double CrossSectionsMultiPions::NNOnePiOrDelta(const G4double ener, const G4int iso, const G4double xsiso) {
249 242
250 /* Article J. Physique 48 (1987)1901-1 243 /* Article J. Physique 48 (1987)1901-1924 "Energy dependence of
251 nucleon-cucleon inelastic total cross 244 nucleon-cucleon inelastic total cross-sections."
252 J. Bystricky, P. La France, F. Lehar, 245 J. Bystricky, P. La France, F. Lehar, F. Perrot, T. Siemiarczuk & P. Winternitz
253 S11PZ= section pp->pp pi0 246 S11PZ= section pp->pp pi0
254 S01PP= section pp->pn pi+ 247 S01PP= section pp->pn pi+
255 S01PZ= section pn->pn pi0 248 S01PZ= section pn->pn pi0
256 S11PM= section pn->pp pi- 249 S11PM= section pn->pp pi-
257 S= X-Section, 1st number : 1 if pp an 250 S= X-Section, 1st number : 1 if pp and 0 if pn
258 2nd number = number of pions, PP= pi+ 251 2nd number = number of pions, PP= pi+; PZ= pi0 ; PM= pi-
259 */ 252 */
260 253
261 const G4double s = ener*ener; 254 const G4double s = ener*ener;
262 G4double plab = 0.001*KinematicsUtils: 255 G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
263 256
264 G4double snnpit1=0.; 257 G4double snnpit1=0.;
265 G4double snnpit=0.; 258 G4double snnpit=0.;
266 G4double s11pz=0.; 259 G4double s11pz=0.;
267 G4double s01pp=0.; 260 G4double s01pp=0.;
268 G4double s01pz=0.; 261 G4double s01pz=0.;
269 G4double s11pm=0.; 262 G4double s11pm=0.;
270 263
271 if ((iso != 0) && (plab < 2.1989)) { 264 if ((iso != 0) && (plab < 2.1989)) {
272 snnpit = xsiso - NNTwoPi(ener, iso 265 snnpit = xsiso - NNTwoPi(ener, iso, xsiso);
273 if (snnpit < 1.e-8) snnpit=0.; 266 if (snnpit < 1.e-8) snnpit=0.;
274 return snnpit; 267 return snnpit;
275 } 268 }
276 else if ((iso == 0) && (plab < 1.7369) 269 else if ((iso == 0) && (plab < 1.7369)) {
277 snnpit = xsiso; 270 snnpit = xsiso;
278 if (snnpit < 1.e-8) snnpit=0.; 271 if (snnpit < 1.e-8) snnpit=0.;
279 return snnpit; 272 return snnpit;
280 } 273 }
281 274
282 //s11pz 275 //s11pz
283 if (plab > 18.) { 276 if (plab > 18.) {
284 s11pz=55.185/std::pow((0.1412*plab 277 s11pz=55.185/std::pow((0.1412*plab+5),2);
285 } 278 }
286 else if (plab > 13.9) { 279 else if (plab > 13.9) {
287 G4double alp=std::log(plab); 280 G4double alp=std::log(plab);
288 s11pz=6.67-13.3*std::pow(plab, -6. 281 s11pz=6.67-13.3*std::pow(plab, -6.18)+0.456*alp*alp-3.29*alp;
289 } 282 }
290 else if (plab >= 0.7765) { 283 else if (plab >= 0.7765) {
291 const G4double b=BystrickyEvaluato 284 const G4double b=BystrickyEvaluator<7>::eval(plab,s11pzOOT,s11pzHC);
292 s11pz=b*b; 285 s11pz=b*b;
293 } 286 }
294 //s01pp 287 //s01pp
295 if (plab >= 0.79624) { 288 if (plab >= 0.79624) {
296 const G4double b=BystrickyEvaluato 289 const G4double b=BystrickyEvaluator<8>::eval(plab,s01ppOOT,s01ppHC);
297 s01pp=b*b; 290 s01pp=b*b;
298 } 291 }
299 292
300 // channel T=1 293 // channel T=1
301 snnpit1=s11pz+s01pp; 294 snnpit1=s11pz+s01pp;
302 if (snnpit1 < 1.e-8) snnpit1=0.; 295 if (snnpit1 < 1.e-8) snnpit1=0.;
303 if (iso != 0) { 296 if (iso != 0) {
304 return snnpit1; 297 return snnpit1;
305 } 298 }
306 299
307 //s01pz 300 //s01pz
308 if (plab > 4.5) { 301 if (plab > 4.5) {
309 s01pz=15289.4/std::pow((11.573*pla 302 s01pz=15289.4/std::pow((11.573*plab+5),2);
310 } 303 }
311 else if (plab >= 0.777) { 304 else if (plab >= 0.777) {
312 const G4double b=BystrickyEvaluato 305 const G4double b=BystrickyEvaluator<4>::eval(plab,s01pzOOT,s01pzHC);
313 s01pz=b*b; 306 s01pz=b*b;
314 } 307 }
315 //s11pm 308 //s11pm
316 if (plab > 14.) { 309 if (plab > 14.) {
317 s11pm=46.68/std::pow((0.2231*plab+ 310 s11pm=46.68/std::pow((0.2231*plab+5),2);
318 } 311 }
319 else if (plab >= 0.788) { 312 else if (plab >= 0.788) {
320 const G4double b=BystrickyEvaluato 313 const G4double b=BystrickyEvaluator<4>::eval(plab,s11pmOOT,s11pmHC);
321 s11pm=b*b; 314 s11pm=b*b;
322 } 315 }
323 316
324 // channel T=0 317 // channel T=0
325 // snnpit=s01pz+2*s11pm-snnpit1; //modi << 318 snnpit=2*(s01pz+2*s11pm)-snnpit1;
326 snnpit = 2*(s01pz+2*s11pm)-snnpit1; <<
327 if (snnpit < 1.e-8) snnpit=0.; 319 if (snnpit < 1.e-8) snnpit=0.;
328 return snnpit; 320 return snnpit;
329 } 321 }
330 322
331 G4double CrossSectionsMultiPions::NNTwoPi( 323 G4double CrossSectionsMultiPions::NNTwoPi(const G4double ener, const G4int iso, const G4double xsiso) {
332 324
333 /* Article J. Physique 48 (1987)1901-1 325 /* Article J. Physique 48 (1987)1901-1924 "Energy dependence of nucleon-cucleon inelastic total cross-sections."
334 J. Bystricky, P. La France, F. Leha 326 J. Bystricky, P. La France, F. Lehar, F. Perrot, T. Siemiarczuk & P. Winternitz
335 S12PM : pp -> pp Pi+ Pi- 327 S12PM : pp -> pp Pi+ Pi-
336 S12ZZ : pp -> pp Pi0 Pi0 328 S12ZZ : pp -> pp Pi0 Pi0
337 S12PP : pp -> nn Pi+ Pi+ 329 S12PP : pp -> nn Pi+ Pi+
338 S02PZ : pp -> pn Pi+ Pi0 330 S02PZ : pp -> pn Pi+ Pi0
339 S02PM : pn -> pn Pi+ Pi- 331 S02PM : pn -> pn Pi+ Pi-
340 S12MZ : pn -> pp Pi- Pi0 332 S12MZ : pn -> pp Pi- Pi0
341 */ 333 */
342 334
343 const G4double s = ener*ener; 335 const G4double s = ener*ener;
344 G4double plab = 0.001*KinematicsUtils: 336 G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
345 337
346 G4double snn2pit=0.; 338 G4double snn2pit=0.;
347 G4double s12pm=0.; 339 G4double s12pm=0.;
348 G4double s12pp=0.; 340 G4double s12pp=0.;
349 G4double s12zz=0.; 341 G4double s12zz=0.;
350 G4double s02pz=0.; 342 G4double s02pz=0.;
351 G4double s02pm=0.; 343 G4double s02pm=0.;
352 G4double s12mz=0.; 344 G4double s12mz=0.;
353 345
354 if (iso==0 && plab<3.33) { << 346 if (iso==0 && plab<3.3) {
355 snn2pit = xsiso - NNOnePiOrDelta(e 347 snn2pit = xsiso - NNOnePiOrDelta(ener, iso, xsiso);
356 if (snn2pit < 1.e-8) snn2pit=0.; 348 if (snn2pit < 1.e-8) snn2pit=0.;
357 return snn2pit; 349 return snn2pit;
358 } 350 }
359 351
360 if (iso != 0) { 352 if (iso != 0) {
361 //s12pm 353 //s12pm
362 if (plab > 15.) { 354 if (plab > 15.) {
363 s12pm=25.977/plab; 355 s12pm=25.977/plab;
364 } 356 }
365 else if (plab >= 1.3817) { 357 else if (plab >= 1.3817) {
366 const G4double b=BystrickyEvaluato 358 const G4double b=BystrickyEvaluator<5>::eval(plab,s12pmOOT,s12pmHC);
367 s12pm=b*b; 359 s12pm=b*b;
368 } 360 }
369 //s12pp 361 //s12pp
370 if (plab > 10.) { 362 if (plab > 10.) {
371 s12pp=141.505/std::pow((-0.1016*pl 363 s12pp=141.505/std::pow((-0.1016*plab-7),2);
372 } 364 }
373 else if (plab >= 1.5739) { 365 else if (plab >= 1.5739) {
374 const G4double b=BystrickyEvaluato 366 const G4double b=BystrickyEvaluator<3>::eval(plab,s12ppOOT,s12ppHC);
375 s12pp=b*b; 367 s12pp=b*b;
376 } 368 }
377 } 369 }
378 //s12zz 370 //s12zz
379 if (plab > 4.) { 371 if (plab > 4.) {
380 s12zz=97.355/std::pow((1.1579*plab 372 s12zz=97.355/std::pow((1.1579*plab+5),2);
381 } 373 }
382 else if (plab >= 1.72207) { 374 else if (plab >= 1.72207) {
383 const G4double b=BystrickyEvaluato 375 const G4double b=BystrickyEvaluator<4>::eval(plab,s12zzOOT,s12zzHC);
384 s12zz=b*b; 376 s12zz=b*b;
385 } 377 }
386 //s02pz 378 //s02pz
387 if (plab > 4.5) { 379 if (plab > 4.5) {
388 s02pz=178.082/std::pow((0.2014*pla 380 s02pz=178.082/std::pow((0.2014*plab+5),2);
389 } 381 }
390 else if (plab >= 1.5656) { 382 else if (plab >= 1.5656) {
391 const G4double b=BystrickyEvaluato 383 const G4double b=BystrickyEvaluator<4>::eval(plab,s02pzOOT,s02pzHC);
392 s02pz=b*b; 384 s02pz=b*b;
393 } 385 }
394 386
395 // channel T=1 387 // channel T=1
396 if (iso != 0) { 388 if (iso != 0) {
397 snn2pit=s12pm+s12pp+s12zz+s02pz; 389 snn2pit=s12pm+s12pp+s12zz+s02pz;
398 if (snn2pit < 1.e-8) snn2pit=0.; 390 if (snn2pit < 1.e-8) snn2pit=0.;
399 return snn2pit; 391 return snn2pit;
400 } 392 }
401 393
402 //s02pm 394 //s02pm
403 if (plab > 5.) { 395 if (plab > 5.) {
404 s02pm=135.826/std::pow(plab,2); 396 s02pm=135.826/std::pow(plab,2);
405 } 397 }
406 else if (plab >= 1.21925) { 398 else if (plab >= 1.21925) {
407 const G4double b=BystrickyEvaluato 399 const G4double b=BystrickyEvaluator<6>::eval(plab,s02pmOOT,s02pmHC);
408 s02pm=b*b; 400 s02pm=b*b;
409 } 401 }
410 //s12mz 402 //s12mz
411 if (plab >= 1.29269) { 403 if (plab >= 1.29269) {
412 const G4double b=BystrickyEvaluato 404 const G4double b=BystrickyEvaluator<4>::eval(plab,s12mzOOT,s12mzHC);
413 s12mz=b*b; 405 s12mz=b*b;
414 } 406 }
415 407
416 // channel T=0 408 // channel T=0
417 // snn2pit=3*(0.5*s02pm+0.5*s12mz-0.5*s <<
418 snn2pit=3*(s02pm+0.5*s12mz-0.5*s02pz-s 409 snn2pit=3*(s02pm+0.5*s12mz-0.5*s02pz-s12zz);
419 if (snn2pit < 1.e-8) snn2pit=0.; 410 if (snn2pit < 1.e-8) snn2pit=0.;
420 return snn2pit; 411 return snn2pit;
421 } 412 }
422 413
423 G4double CrossSectionsMultiPions::NNThreeP 414 G4double CrossSectionsMultiPions::NNThreePi(const G4double ener, const G4int iso, const G4double xsiso, const G4double xs1pi, const G4double xs2pi) {
424 415
425 const G4double s = ener*ener; 416 const G4double s = ener*ener;
426 G4double plab = 0.001*KinematicsUtils: 417 G4double plab = 0.001*KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
427 418
428 G4double snn3pit=0.; 419 G4double snn3pit=0.;
429 420
430 if (iso == 0) { 421 if (iso == 0) {
431 // channel T=0 422 // channel T=0
432 if (plab > 7.2355) { 423 if (plab > 7.2355) {
433 return 46.72/std::pow((plab - 424 return 46.72/std::pow((plab - 5.8821),2);
434 } 425 }
435 else { 426 else {
436 snn3pit=xsiso-xs1pi-xs2pi; 427 snn3pit=xsiso-xs1pi-xs2pi;
437 if (snn3pit < 1.e-8) snn3pit=0 428 if (snn3pit < 1.e-8) snn3pit=0.;
438 return snn3pit; 429 return snn3pit;
439 } 430 }
440 } 431 }
441 else { 432 else {
442 // channel T=1 433 // channel T=1
443 if (plab > 7.206) { 434 if (plab > 7.206) {
444 return 5592.92/std::pow((plab+ 435 return 5592.92/std::pow((plab+14.9764),2);
445 } 436 }
446 else if (plab > 2.1989){ 437 else if (plab > 2.1989){
447 snn3pit=xsiso-xs1pi-xs2pi; 438 snn3pit=xsiso-xs1pi-xs2pi;
448 if (snn3pit < 1.e-8) snn3pit=0 439 if (snn3pit < 1.e-8) snn3pit=0.;
449 return snn3pit; 440 return snn3pit;
450 } 441 }
451 else return snn3pit; 442 else return snn3pit;
452 } 443 }
453 } 444 }
454 445
455 G4double CrossSectionsMultiPions::NNOnePi( 446 G4double CrossSectionsMultiPions::NNOnePi(Particle const * const particle1, Particle const * const particle2) {
456 // Cross section for nucleon-nucleon d 447 // Cross section for nucleon-nucleon directly producing one pion
457 448
458 const G4int iso=ParticleTable::getIsos 449 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
459 if (iso!=0) // If pp or nn we choose t << 450 if (iso!=0)
460 return 0.; 451 return 0.;
461 452
462 const G4double ener=KinematicsUtils::t 453 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
463 454
464 const G4double xsiso2=NNInelasticIso(e 455 const G4double xsiso2=NNInelasticIso(ener, 2);
465 const G4double xsiso0=NNInelasticIso(e 456 const G4double xsiso0=NNInelasticIso(ener, 0);
466 return 0.25*(NNOnePiOrDelta(ener, 0, x 457 return 0.25*(NNOnePiOrDelta(ener, 0, xsiso0)+ NNOnePiOrDelta(ener, 2, xsiso2));
467 } 458 }
468 459
469 G4double CrossSectionsMultiPions::NNOnePiO 460 G4double CrossSectionsMultiPions::NNOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
470 // Cross section for nucleon-nucleon d 461 // Cross section for nucleon-nucleon directly producing one pion or producing a nucleon-delta pair
471 const G4double ener=KinematicsUtils::t 462 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
472 const G4int iso=ParticleTable::getIsos 463 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
473 464
474 const G4double xsiso2=NNInelasticIso(e 465 const G4double xsiso2=NNInelasticIso(ener, 2);
475 if (iso != 0) 466 if (iso != 0)
476 return NNOnePiOrDelta(ener, iso, xsi 467 return NNOnePiOrDelta(ener, iso, xsiso2);
477 else { 468 else {
478 const G4double xsiso0=NNInelasticIso 469 const G4double xsiso0=NNInelasticIso(ener, 0);
479 return 0.5*(NNOnePiOrDelta(ener, 0, 470 return 0.5*(NNOnePiOrDelta(ener, 0, xsiso0)+ NNOnePiOrDelta(ener, 2, xsiso2));
480 } 471 }
481 } 472 }
482 473
483 G4double CrossSectionsMultiPions::NNTwoPi( 474 G4double CrossSectionsMultiPions::NNTwoPi(Particle const * const particle1, Particle const * const particle2) {
484 // 475 //
485 // Nucleon-Nucleon producing one p 476 // Nucleon-Nucleon producing one pion cross sections
486 // 477 //
487 const G4double ener=KinematicsUtils::t 478 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
488 const G4int iso=ParticleTable::getIsos 479 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
489 480
490 481
491 const G4double xsiso2=NNInelasticIso(e 482 const G4double xsiso2=NNInelasticIso(ener, 2);
492 if (iso != 0) { 483 if (iso != 0) {
493 return NNTwoPi(ener, 2, xsiso2); 484 return NNTwoPi(ener, 2, xsiso2);
494 } 485 }
495 else { 486 else {
496 const G4double xsiso0=NNInelasticI 487 const G4double xsiso0=NNInelasticIso(ener, 0);
497 return 0.5*(NNTwoPi(ener, 0, xsiso 488 return 0.5*(NNTwoPi(ener, 0, xsiso0)+ NNTwoPi(ener, 2, xsiso2));
498 } 489 }
499 return 0.0; // Should never reach this 490 return 0.0; // Should never reach this point
500 } 491 }
501 492
502 G4double CrossSectionsMultiPions::NNThreeP 493 G4double CrossSectionsMultiPions::NNThreePi(Particle const * const particle1, Particle const * const particle2) {
503 // 494 //
504 // Nucleon-Nucleon producing one p 495 // Nucleon-Nucleon producing one pion cross sections
505 // 496 //
506 497
507 const G4double ener=KinematicsUtils::t 498 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
508 const G4int iso=ParticleTable::getIsos 499 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
509 500
510 501
511 const G4double xsiso2=NNInelasticIso(e 502 const G4double xsiso2=NNInelasticIso(ener, 2);
512 const G4double xs1pi2=NNOnePiOrDelta(e 503 const G4double xs1pi2=NNOnePiOrDelta(ener, 2, xsiso2);
513 const G4double xs2pi2=NNTwoPi(ener, 2, 504 const G4double xs2pi2=NNTwoPi(ener, 2, xsiso2);
514 if (iso != 0) 505 if (iso != 0)
515 return NNThreePi(ener, 2, xsiso2, xs 506 return NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
516 else { 507 else {
517 const G4double xsiso0=NNInelasticIso 508 const G4double xsiso0=NNInelasticIso(ener, 0);
518 const G4double xs1pi0=NNOnePiOrDelta 509 const G4double xs1pi0=NNOnePiOrDelta(ener, 0, xsiso0);
519 const G4double xs2pi0=NNTwoPi(ener, 510 const G4double xs2pi0=NNTwoPi(ener, 0, xsiso0);
520 return 0.5*(NNThreePi(ener, 0, xsiso 511 return 0.5*(NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
521 } 512 }
522 } 513 }
523 514
524 G4double CrossSectionsMultiPions::NNFourPi 515 G4double CrossSectionsMultiPions::NNFourPi(Particle const * const particle1, Particle const * const particle2) {
525 const G4double s = KinematicsUtils::squa 516 const G4double s = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
526 if(s<6.25E6) 517 if(s<6.25E6)
527 return 0.; 518 return 0.;
528 const G4double sigma = NNTot(particle1, 519 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.); << 520 return ((sigma>0.) ? sigma : 0.);
530 } 521 }
531 522
532 G4double CrossSectionsMultiPions::NNToxPiN 523 G4double CrossSectionsMultiPions::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
533 // 524 //
534 // Nucleon-Nucleon producing xpi pio 525 // Nucleon-Nucleon producing xpi pions cross sections
535 // 526 //
536 // assert(xpi>0 && xpi<=nMaxPiNN); << 527 // assert(xpi>0 && xpi<5);
537 // assert(particle1->isNucleon() && particle2- 528 // assert(particle1->isNucleon() && particle2->isNucleon());
538 529
539 if (xpi == 1) 530 if (xpi == 1)
540 return NNOnePi(particle1, particle2); 531 return NNOnePi(particle1, particle2);
541 else if (xpi == 2) 532 else if (xpi == 2)
542 return NNTwoPi(particle1, particle2); 533 return NNTwoPi(particle1, particle2);
543 else if (xpi == 3) 534 else if (xpi == 3)
544 return NNThreePi(particle1, particle2) 535 return NNThreePi(particle1, particle2);
545 else if (xpi == 4) << 536 else // if (xpi == 4)
546 return NNFourPi(particle1, particle2); 537 return NNFourPi(particle1, particle2);
547 else // should never reach this point <<
548 return 0.; <<
549 } 538 }
550 539
551 540
552 G4double CrossSectionsMultiPions::spnPiPlusP 541 G4double CrossSectionsMultiPions::spnPiPlusPHE(const G4double x) {
553 // HE and LE pi- p and pi+ n 542 // HE and LE pi- p and pi+ n
554 G4double ramass = 0.0; 543 G4double ramass = 0.0;
555 544
556 if(x <= 1306.78) { << 545 if(x <= 1306.0) {
557 G4double y = x*x; 546 G4double y = x*x;
558 G4double q2; 547 G4double q2;
559 q2=(y-std::pow(1076.0, 2))*(y-std::pow( 548 q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
560 if (q2 > 0.) { 549 if (q2 > 0.) {
561 G4double q3=std::pow(q2, 3./2.); 550 G4double q3=std::pow(q2, 3./2.);
562 G4double f3=q3/(q3+std::pow(180.0, 3 551 G4double f3=q3/(q3+std::pow(180.0, 3));
563 G4double sdel; 552 G4double sdel;
564 sdel=326.5/(std::pow((x-1215.0-ramass)*2.0 553 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); 554 return sdel*f3*(1.0-5.0*ramass/1215.0);
566 } 555 }
567 else { 556 else {
568 return 0; 557 return 0;
569 } 558 }
570 } 559 }
571 if(x <= 1754.0) { 560 if(x <= 1754.0) {
572 return -2.33730e-06*std::pow(x, 3)+1.138 561 return -2.33730e-06*std::pow(x, 3)+1.13819e-02*std::pow(x,2)
573 -1.83993e+01*x+9893.4; 562 -1.83993e+01*x+9893.4;
574 } else if (x <= 2150.0) { 563 } else if (x <= 2150.0) {
575 return 1.13531e-06*std::pow(x, 3)-6.9169 564 return 1.13531e-06*std::pow(x, 3)-6.91694e-03*std::pow(x, 2)
576 +1.39907e+01*x-9360.76; 565 +1.39907e+01*x-9360.76;
577 } else { 566 } else {
578 return -3.18087*std::log(x)+52.9784; 567 return -3.18087*std::log(x)+52.9784;
579 } 568 }
580 } 569 }
581 570
582 G4double CrossSectionsMultiPions::spnPiMinus 571 G4double CrossSectionsMultiPions::spnPiMinusPHE(const G4double x) {
583 // HE pi- p and pi+ n 572 // HE pi- p and pi+ n
584 G4double ramass = 0.0; 573 G4double ramass = 0.0;
585 574
586 if(x <= 1275.8) { 575 if(x <= 1275.8) {
587 G4double y = x*x; 576 G4double y = x*x;
588 G4double q2; 577 G4double q2;
589 q2=(y-std::pow(1076.0, 2))*(y-std::pow( 578 q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
590 if (q2 > 0.) { 579 if (q2 > 0.) {
591 G4double q3=std::pow(q2, 3./2.); 580 G4double q3=std::pow(q2, 3./2.);
592 G4double f3=q3/(q3+std::pow(180.0, 3 581 G4double f3=q3/(q3+std::pow(180.0, 3));
593 G4double sdel; 582 G4double sdel;
594 sdel=326.5/(std::pow((x-1215.0-ramass)*2.0 583 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.; 584 return sdel*f3*(1.0-5.0*ramass/1215.0)/3.;
596 } 585 }
597 else { 586 else {
598 return 0; 587 return 0;
599 } 588 }
600 } 589 }
601 if(x <= 1495.0) { 590 if(x <= 1495.0) {
602 return 0.00120683*(x-1372.52)*(x-1372.52 591 return 0.00120683*(x-1372.52)*(x-1372.52)+26.2058;
603 } else if(x <= 1578.0) { 592 } else if(x <= 1578.0) {
604 return 1.15873e-05*x*x+49965.6/((x-1519. 593 return 1.15873e-05*x*x+49965.6/((x-1519.59)*(x-1519.59)+2372.55);
605 } else if(x <= 2028.4) { 594 } else if(x <= 2028.4) {
606 return 34.0248+43262.2/((x-1681.65)*(x-1 595 return 34.0248+43262.2/((x-1681.65)*(x-1681.65)+1689.35);
607 } else if(x <= 7500.0) { 596 } else if(x <= 7500.0) {
608 return 3.3e-7*(x-7500.0)*(x-7500.0)+24.5 597 return 3.3e-7*(x-7500.0)*(x-7500.0)+24.5;
609 } else { 598 } else {
610 return 24.5; 599 return 24.5;
611 } 600 }
612 } 601 }
613 602
614 G4double CrossSectionsMultiPions::total(Part 603 G4double CrossSectionsMultiPions::total(Particle const * const p1, Particle const * const p2) {
615 G4double inelastic; 604 G4double inelastic;
616 if(p1->isNucleon() && p2->isNucleon()) { 605 if(p1->isNucleon() && p2->isNucleon()) {
617 return NNTot(p1, p2); 606 return NNTot(p1, p2);
618 } else if((p1->isNucleon() && p2->isDelta( 607 } else if((p1->isNucleon() && p2->isDelta()) ||
619 (p1->isDelta() && p2->isNucleon( 608 (p1->isDelta() && p2->isNucleon())) {
620 inelastic = NDeltaToNN(p1, p2); 609 inelastic = NDeltaToNN(p1, p2);
621 } else if((p1->isNucleon() && p2->isPion() 610 } else if((p1->isNucleon() && p2->isPion()) ||
622 (p1->isPion() && p2->isNucleon() 611 (p1->isPion() && p2->isNucleon())) {
623 return piNTot(p1,p2); 612 return piNTot(p1,p2);
624 } else { 613 } else {
625 inelastic = 0.; 614 inelastic = 0.;
626 } 615 }
627 616
628 return inelastic + elastic(p1, p2); 617 return inelastic + elastic(p1, p2);
629 } 618 }
630 619
631 620
632 G4double CrossSectionsMultiPions::piNIne(Par << 621 G4double CrossSectionsMultiPions::piNIne(Particle const * const particle1, Particle const * const particle2) {
633 // piN inelastic cross section (Delta << 622 // piN inelastic cross section (Delta excluded)
634 << 623
635 const Particle *pion; << 624 const Particle *pion;
636 const Particle *nucleon; << 625 const Particle *nucleon;
637 if(particle1->isNucleon()) { << 626 if(particle1->isNucleon()) {
638 nucleon = particle1; << 627 nucleon = particle1;
639 pion = particle2; << 628 pion = particle2;
640 } else { << 629 } else {
641 pion = particle1; << 630 pion = particle1;
642 nucleon = particle2; << 631 nucleon = particle2;
643 } << 632 }
644 // assert(pion->isPion()); 633 // assert(pion->isPion());
645 << 634
646 const G4double pLab = KinematicsUtils::mom << 635 const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
647 << 636
648 // these limits correspond to sqrt(s)=1230 << 637 // these limits correspond to sqrt(s)=1230 and 20000 MeV
649 if(pLab>212677. || pLab<296.367) << 638 if(pLab>212677. || pLab<296.367)
650 return 0.0; << 639 return 0.0;
651 << 640
652 const G4int ipit3 = ParticleTable::getIsos << 641 const G4int ipit3 = ParticleTable::getIsospin(pion->getType());
653 const G4int ind2t3 = ParticleTable::getIso << 642 const G4int ind2t3 = ParticleTable::getIsospin(nucleon->getType());
654 const G4int cg = 4 + ind2t3*ipit3; << 643 const G4int cg = 4 + ind2t3*ipit3;
655 // assert(cg==2 || cg==4 || cg==6); 644 // assert(cg==2 || cg==4 || cg==6);
656 << 645
657 // const G4double p1=1e-3*pLab; << 646 const G4double p1=1e-3*pLab;
658 // const G4double p2=std::log(p1); << 647 const G4double p2=std::log(p1);
659 G4double xpipp = 0.0; << 648 G4double xpipp = 0.0;
660 G4double xpimp = 0.0; << 649 G4double xpimp = 0.0;
661 << 650
662 if(cg!=2) { << 651 if(cg!=2) {
663 // x-section pi+ p inelastique : << 652 // x-section pi+ p inelastique :
664 xpipp=piPluspIne(pion,nucleon); << 653 if(p1<=0.75)
665 << 654 xpipp=17.965*std::pow(p1, 5.4606);
666 if(cg==6) // cas pi+ p et pi- n << 655 else
667 return xpipp; << 656 xpipp=24.3-12.3*std::pow(p1, -1.91)+0.324*p2*p2-2.44*p2;
668 } << 657 if(cg==6)
669 << 658 // cas pi+ p et pi- n
670 // x-section pi- p inelastique : << 659 return xpipp;
671 xpimp=piMinuspIne(pion,nucleon); << 660 }
672 << 661
673 if(cg==2) // cas pi- p et pi+ n << 662 // x-section pi- p inelastique :
674 return xpimp; << 663 if(p1 <= 0.4731)
675 else // cas pi0 p et pi0 n << 664 xpimp=0;
676 return 0.5*(xpipp+xpimp); << 665 else
677 } << 666 xpimp=26.6-7.18*std::pow(p1, -1.86)+0.327*p2*p2-2.81*p2;
678 << 667 if(xpimp<0.)
>> 668 xpimp=0;
>> 669
>> 670 if(cg==2) // cas pi- p et pi+ n
>> 671 return xpimp;
>> 672 else // cas pi0 p et pi0 n
>> 673 return 0.5*(xpipp+xpimp);
>> 674 }
>> 675
679 G4double CrossSectionsMultiPions::piNToDelta 676 G4double CrossSectionsMultiPions::piNToDelta(Particle const * const particle1, Particle const * const particle2) {
680 // piN Delta production 677 // piN Delta production
681 678
682 G4double x = KinematicsUtils::totalEnergyI 679 G4double x = KinematicsUtils::totalEnergyInCM(particle1, particle2);
683 if(x>20000.) return 0.0; // no cross secti 680 if(x>20000.) return 0.0; // no cross section above this value
684 681
685 G4int ipit3 = 0; 682 G4int ipit3 = 0;
686 G4int ind2t3 = 0; 683 G4int ind2t3 = 0;
687 const G4double ramass = 0.0; 684 const G4double ramass = 0.0;
688 685
689 if(particle1->isPion()) { 686 if(particle1->isPion()) {
690 ipit3 = ParticleTable::getIsospin(partic 687 ipit3 = ParticleTable::getIsospin(particle1->getType());
691 ind2t3 = ParticleTable::getIsospin(parti 688 ind2t3 = ParticleTable::getIsospin(particle2->getType());
692 } else if(particle2->isPion()) { 689 } else if(particle2->isPion()) {
693 ipit3 = ParticleTable::getIsospin(partic 690 ipit3 = ParticleTable::getIsospin(particle2->getType());
694 ind2t3 = ParticleTable::getIsospin(parti 691 ind2t3 = ParticleTable::getIsospin(particle1->getType());
695 } 692 }
696 693
697 const G4double y=x*x; 694 const G4double y=x*x;
698 const G4double q2=(y-1076.0*1076.0)*(y-800 695 const G4double q2=(y-1076.0*1076.0)*(y-800.0*800.0)/y/4.0;
699 if (q2 <= 0.) { 696 if (q2 <= 0.) {
700 return 0.0; 697 return 0.0;
701 } 698 }
702 const G4double q3 = std::pow(std::sqrt(q2) 699 const G4double q3 = std::pow(std::sqrt(q2),3);
703 const G4double f3 = q3/(q3 + 5832000.); // 700 const G4double f3 = q3/(q3 + 5832000.); // 5832000 = 180^3
704 G4double sdelResult = 326.5/(std::pow((x-1 701 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/12 702 sdelResult = sdelResult*(1.0-5.0*ramass/1215.0);
706 const G4int cg = 4 + ind2t3*ipit3; 703 const G4int cg = 4 + ind2t3*ipit3;
707 sdelResult = sdelResult*f3*cg/6.0; 704 sdelResult = sdelResult*f3*cg/6.0;
708 705
709 return sdelResult; 706 return sdelResult;
710 } 707 }
711 708
712 G4double CrossSectionsMultiPions::piNTot(Par 709 G4double CrossSectionsMultiPions::piNTot(Particle const * const particle1, Particle const * const particle2) {
713 // FUNCTION SPN(X,IND2T3,IPIT3,f17) 710 // FUNCTION SPN(X,IND2T3,IPIT3,f17)
714 // SIGMA(PI+ + P) IN THE (3,3) REGION 711 // SIGMA(PI+ + P) IN THE (3,3) REGION
715 // NEW FIT BY J.VANDERMEULEN + FIT BY Th 712 // NEW FIT BY J.VANDERMEULEN + FIT BY Th AOUST ABOVE (3,3) RES
716 // CONST AT L 713 // CONST AT LOW AND VERY HIGH ENERGY
717 // COMMON/BL8/RATHR,RAMASS 714 // COMMON/BL8/RATHR,RAMASS REL21800
718 // integer f17 715 // integer f17
719 // RATHR and RAMASS are always 0.0!!! 716 // RATHR and RAMASS are always 0.0!!!
720 717
721 G4double x = KinematicsUtils::totalEnergyI 718 G4double x = KinematicsUtils::totalEnergyInCM(particle1, particle2);
722 719
723 G4int ipit3 = 0; 720 G4int ipit3 = 0;
724 G4int ind2t3 = 0; 721 G4int ind2t3 = 0;
725 722
726 if(particle1->isPion()) { 723 if(particle1->isPion()) {
727 ipit3 = ParticleTable::getIsospin(partic 724 ipit3 = ParticleTable::getIsospin(particle1->getType());
728 ind2t3 = ParticleTable::getIsospin(parti 725 ind2t3 = ParticleTable::getIsospin(particle2->getType());
729 } else if(particle2->isPion()) { 726 } else if(particle2->isPion()) {
730 ipit3 = ParticleTable::getIsospin(partic 727 ipit3 = ParticleTable::getIsospin(particle2->getType());
731 ind2t3 = ParticleTable::getIsospin(parti 728 ind2t3 = ParticleTable::getIsospin(particle1->getType());
732 } 729 }
733 730
734 G4double spnResult=0.0; 731 G4double spnResult=0.0;
735 732
736 // HE pi+ p and pi- n 733 // HE pi+ p and pi- n
737 if((ind2t3 == 1 && ipit3 == 2) || (ind2t 734 if((ind2t3 == 1 && ipit3 == 2) || (ind2t3 == -1 && ipit3 == -2))
738 spnResult=spnPiPlusPHE(x); 735 spnResult=spnPiPlusPHE(x);
739 else if((ind2t3 == 1 && ipit3 == -2) || 736 else if((ind2t3 == 1 && ipit3 == -2) || (ind2t3 == -1 && ipit3 == 2))
740 spnResult=spnPiMinusPHE(x); 737 spnResult=spnPiMinusPHE(x);
741 else if(ipit3 == 0) spnResult = (spnPiPl 738 else if(ipit3 == 0) spnResult = (spnPiPlusPHE(x) + spnPiMinusPHE(x))/2.0; // (spnpipphe(x)+spnpimphe(x))/2.0
742 else { 739 else {
743 INCL_ERROR("Unknown configuration!\n" 740 INCL_ERROR("Unknown configuration!\n" << particle1->print() << particle2->print() << '\n');
744 } 741 }
745 742
746 return spnResult; 743 return spnResult;
747 } 744 }
748 745
749 G4double CrossSectionsMultiPions::NDeltaToNN 746 G4double CrossSectionsMultiPions::NDeltaToNN(Particle const * const p1, Particle const * const p2) {
750 const G4int isospin = ParticleTable::getIs 747 const G4int isospin = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
751 if(isospin==4 || isospin==-4) return 0.0; 748 if(isospin==4 || isospin==-4) return 0.0;
752 749
753 G4double s = KinematicsUtils::squareTotalE 750 G4double s = KinematicsUtils::squareTotalEnergyInCM(p1, p2);
754 G4double Ecm = std::sqrt(s); 751 G4double Ecm = std::sqrt(s);
755 G4int deltaIsospin; 752 G4int deltaIsospin;
756 G4double deltaMass; 753 G4double deltaMass;
757 if(p1->isDelta()) { 754 if(p1->isDelta()) {
758 deltaIsospin = ParticleTable::getIsospin 755 deltaIsospin = ParticleTable::getIsospin(p1->getType());
759 deltaMass = p1->getMass(); 756 deltaMass = p1->getMass();
760 } else { 757 } else {
761 deltaIsospin = ParticleTable::getIsospin 758 deltaIsospin = ParticleTable::getIsospin(p2->getType());
762 deltaMass = p2->getMass(); 759 deltaMass = p2->getMass();
763 } 760 }
764 761
765 if(Ecm <= 938.3 + deltaMass) { 762 if(Ecm <= 938.3 + deltaMass) {
766 return 0.0; 763 return 0.0;
767 } 764 }
768 765
769 if(Ecm < 938.3 + deltaMass + 2.0) { 766 if(Ecm < 938.3 + deltaMass + 2.0) {
770 Ecm = 938.3 + deltaMass + 2.0; 767 Ecm = 938.3 + deltaMass + 2.0;
771 s = Ecm*Ecm; 768 s = Ecm*Ecm;
772 } 769 }
773 770
774 const G4double x = (s - 4.*ParticleTable:: 771 const G4double x = (s - 4.*ParticleTable::effectiveNucleonMass2) /
775 (s - std::pow(ParticleTable::effectiveNu 772 (s - std::pow(ParticleTable::effectiveNucleonMass + deltaMass, 2));
776 const G4double y = s/(s - std::pow(deltaMa 773 const G4double y = s/(s - std::pow(deltaMass - ParticleTable::effectiveNucleonMass, 2));
777 /* Concerning the way we calculate the lab 774 /* Concerning the way we calculate the lab momentum, see the considerations
778 * in CrossSections::elasticNNLegacy(). 775 * in CrossSections::elasticNNLegacy().
779 */ 776 */
780 G4double sDelta; 777 G4double sDelta;
781 const G4double xsiso2=NNInelasticIso(Ecm, 778 const G4double xsiso2=NNInelasticIso(Ecm, 2);
782 if (isospin != 0) 779 if (isospin != 0)
783 sDelta = NNOnePiOrDelta(Ecm, isospin, xs 780 sDelta = NNOnePiOrDelta(Ecm, isospin, xsiso2);
784 else { 781 else {
785 const G4double xsiso0=NNInelasticIso(Ecm 782 const G4double xsiso0=NNInelasticIso(Ecm, 0);
786 sDelta = 0.25*(NNOnePiOrDelta(Ecm, 0, xs 783 sDelta = 0.25*(NNOnePiOrDelta(Ecm, 0, xsiso0)+ NNOnePiOrDelta(Ecm, 2, xsiso2));
787 } 784 }
788 G4double result = 0.5 * x * y * sDelta; 785 G4double result = 0.5 * x * y * sDelta;
789 /* modification for pion-induced cascade ( 786 /* modification for pion-induced cascade (see JC and MC LEMAIRE,NPA489(88)781
790 * result=3.*result 787 * result=3.*result
791 * pi absorption increased also for intern 788 * pi absorption increased also for internal pions (7/3/01)
792 */ 789 */
793 result *= 3.*(32.0 + isospin * isospin * ( 790 result *= 3.*(32.0 + isospin * isospin * (deltaIsospin * deltaIsospin - 5))/64.0;
794 result /= 1.0 + 0.25 * (isospin * isospin) 791 result /= 1.0 + 0.25 * (isospin * isospin);
795 return result; 792 return result;
796 } 793 }
797 794
798 G4double CrossSectionsMultiPions::NNToNDelta 795 G4double CrossSectionsMultiPions::NNToNDelta(Particle const * const p1, Particle const * const p2) {
799 // assert(p1->isNucleon() && p2->isNucleon()); 796 // assert(p1->isNucleon() && p2->isNucleon());
800 const G4int isospin = ParticleTable::getIs 797 const G4int isospin = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
801 G4double sigma = NNOnePiOrDelta(p1, p2); 798 G4double sigma = NNOnePiOrDelta(p1, p2);
802 if(isospin==0) 799 if(isospin==0)
803 sigma *= 0.5; 800 sigma *= 0.5;
804 return sigma; 801 return sigma;
805 } 802 }
806 803
807 G4double CrossSectionsMultiPions::elastic(Pa 804 G4double CrossSectionsMultiPions::elastic(Particle const * const p1, Particle const * const p2) {
808 // if(!p1->isPion() && !p2->isPion()){ << 805 if(!p1->isPion() && !p2->isPion()){
809 if((p1->isNucleon()||p1->isDelta()) && (p2 <<
810 return NNElastic(p1, p2); 806 return NNElastic(p1, p2);
811 } 807 }
812 // else if (p1->isNucleon() || p2->isNucleo << 808 else if (p1->isNucleon() || p2->isNucleon()){
813 else if ((p1->isNucleon() && p2->isPion()) | <<
814 G4double pielas = piNTot(p1,p2) - piNIne 809 G4double pielas = piNTot(p1,p2) - piNIne(p1,p2) - piNToDelta(p1,p2);
815 if (pielas < 0.){ 810 if (pielas < 0.){
816 pielas = 0.; 811 pielas = 0.;
817 } 812 }
818 // return piNTot(p1,p2) - piNIne(p1,p2) 813 // return piNTot(p1,p2) - piNIne(p1,p2) - piNToDelta(p1,p2);
819 return pielas; 814 return pielas;
820 } 815 }
821 else { 816 else {
822 return 0.0; 817 return 0.0;
823 } 818 }
824 } 819 }
825 820
826 G4double CrossSectionsMultiPions::calculateN 821 G4double CrossSectionsMultiPions::calculateNNAngularSlope(G4double pl, G4int iso) {
827 G4double x = 0.001 * pl; // Change to GeV 822 G4double x = 0.001 * pl; // Change to GeV
828 if(iso != 0) { 823 if(iso != 0) {
829 if(pl <= 2000.0) { 824 if(pl <= 2000.0) {
830 x = std::pow(x, 8); 825 x = std::pow(x, 8);
831 return 5.5e-6 * x/(7.7 + x); 826 return 5.5e-6 * x/(7.7 + x);
832 } else { 827 } else {
833 return (5.34 + 0.67*(x - 2.0)) * 1.0e- 828 return (5.34 + 0.67*(x - 2.0)) * 1.0e-6;
834 } 829 }
835 } else { 830 } else {
836 if(pl < 800.0) { 831 if(pl < 800.0) {
837 G4double b = (7.16 - 1.63*x) * 1.0e-6; 832 G4double b = (7.16 - 1.63*x) * 1.0e-6;
838 return b/(1.0 + std::exp(-(x - 0.45)/0 833 return b/(1.0 + std::exp(-(x - 0.45)/0.05));
839 } else if(pl < 1100.0) { 834 } else if(pl < 1100.0) {
840 return (9.87 - 4.88 * x) * 1.0e-6; 835 return (9.87 - 4.88 * x) * 1.0e-6;
841 } else { 836 } else {
842 return (3.68 + 0.76*x) * 1.0e-6; 837 return (3.68 + 0.76*x) * 1.0e-6;
843 } 838 }
844 } 839 }
845 return 0.0; // Should never reach this poi 840 return 0.0; // Should never reach this point
846 } 841 }
847 842
848 843
849 G4double CrossSectionsMultiPions::piNToxPi 844 G4double CrossSectionsMultiPions::piNToxPiN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
850 // 845 //
851 // pion-Nucleon producing xpi pion 846 // pion-Nucleon producing xpi pions cross sections
852 // 847 //
853 const Particle *pion; << 848 // assert(xpi>1 && xpi<5);
854 const Particle *nucleon; << 849 if (xpi == 2)
855 if(particle1->isNucleon()) { << 850 return piNOnePi(particle1,particle2);
856 nucleon = particle1; << 851 else if (xpi == 3)
857 pion = particle2; << 852 return piNTwoPi(particle1,particle2);
858 } else { <<
859 pion = particle1; <<
860 nucleon = particle2; <<
861 } <<
862 // assert(xpi>1 && xpi<=nMaxPiPiN); <<
863 // assert((particle1->isNucleon() && particle2 <<
864 const G4double plab = KinematicsUtils: <<
865 if (xpi == 2) { <<
866 G4double OnePi=piNOnePi(particle1,partic <<
867 if (OnePi < 1.e-09) OnePi = 0.; <<
868 return OnePi; <<
869 } <<
870 else if (xpi == 3){ <<
871 G4double TwoPi=piNTwoPi(particle1,partic <<
872 if (TwoPi < 1.e-09) TwoPi = 0.; <<
873 return TwoPi; <<
874 } <<
875 else if (xpi == 4) { 853 else if (xpi == 4) {
876 G4double piNThreePi = piNIne(parti << 854 const G4double piNThreePi = piNIne(particle1,particle2) - piNOnePi(particle1,particle2) - piNTwoPi(particle1,particle2);
877 if (piNThreePi < 1.e-09 || plab < <<
878 return piNThreePi; 855 return piNThreePi;
879 } else // should never reach this poin << 856 }
880 return 0.0; << 857 return 0.0; // Should never reach this point
881 } 858 }
882 859
883 G4double CrossSectionsMultiPions::piNOnePi(P << 860 G4double CrossSectionsMultiPions::piNOnePi(Particle const * const particle1, Particle const * const particle2) {
884 const Particle *pion; << 861 const Particle *pion;
885 const Particle *nucleon; << 862 const Particle *nucleon;
886 if(particle1->isNucleon()) { << 863 if(particle1->isNucleon()) {
887 nucleon = particle1; << 864 nucleon = particle1;
888 pion = particle2; << 865 pion = particle2;
889 } else { << 866 } else {
890 pion = particle1; << 867 pion = particle1;
891 nucleon = particle2; << 868 nucleon = particle2;
892 } << 869 }
893 // assert(pion->isPion()); <<
894 <<
895 const G4double pLab = KinematicsUtils::mom <<
896 <<
897 // this limit corresponds to sqrt(s)=1230 <<
898 if(pLab<296.367) <<
899 return 0.0; <<
900 <<
901 const G4int ipi = ParticleTable::getIsospi <<
902 const G4int ind2 = ParticleTable::getIsosp <<
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, p <<
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 <<
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(partic <<
929 if(s1pin >= elas && pLab < 410.) s1pin = <<
930 if (s1pin > inelastic) <<
931 s1pin = inelastic; <<
932 return s1pin; <<
933 } <<
934 } <<
935 <<
936 G4double CrossSectionsMultiPions::piNTwoPi(P <<
937 // <<
938 // pion-nucleon interaction, producing <<
939 // fit from Landolt-Bornstein multipli <<
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::mom <<
954 const G4double elas = elastic(pion, nucleo <<
955 <<
956 // this limit corresponds to sqrt(s)=1230 <<
957 if(pLab<296.367) <<
958 return 0.0; <<
959 <<
960 const G4int ipi = ParticleTable::getIsospi <<
961 const G4int ind2 = ParticleTable::getIsosp <<
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 <<
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 <<
988 // piPlusP inelastic cross section (D <<
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::mo <<
1002 <<
1003 // these limits correspond to sqrt(s)=123 <<
1004 if(pLab>212677. || pLab<296.367) <<
1005 return 0.0; <<
1006 <<
1007 // const G4int ipit3 = ParticleTable::getI <<
1008 // const G4int ind2t3 = ParticleTable::get <<
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.3 <<
1021 // cas pi+ p et pi- n <<
1022 return xpipp; <<
1023 <<
1024 } <<
1025 <<
1026 G4double CrossSectionsMultiPions::piMinuspI <<
1027 // piMinusp inelastic cross section <<
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::mo <<
1041 <<
1042 // these limits correspond to sqrt(s)=123 <<
1043 if(pLab>212677. || pLab<296.367) <<
1044 return 0.0; <<
1045 <<
1046 // const G4int ipit3 = ParticleTable::getI <<
1047 // const G4int ind2t3 = ParticleTable::get <<
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.3 <<
1060 if(xpimp<0.) <<
1061 xpimp=0; <<
1062 <<
1063 // cas pi- p et pi+ n <<
1064 return xpimp; <<
1065 <<
1066 } <<
1067 <<
1068 G4double CrossSectionsMultiPions::piPluspOn <<
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::mo <<
1081 <<
1082 // this limit corresponds to sqrt(s)=1230 <<
1083 if(pLab<296.367) <<
1084 return 0.0; <<
1085 <<
1086 // const G4int ipi = ParticleTable::getI <<
1087 // const G4int ind2 = ParticleTable::get <<
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 <<
1096 tamp6=piPluspIne(particle1, particle2); <<
1097 else <<
1098 tamp6=0.204+18.2*std::pow(p1, -1.72)+6. <<
1099 <<
1100 // CAS PI+ P ET PI- N <<
1101 return tamp6; <<
1102 <<
1103 } <<
1104 <<
1105 G4double CrossSectionsMultiPions::piMinuspO <<
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::mo <<
1118 <<
1119 // this limit corresponds to sqrt(s)=1230 <<
1120 if(pLab<296.367) <<
1121 return 0.0; <<
1122 <<
1123 // const G4int ipi = ParticleTable::getI <<
1124 // const G4int ind2 = ParticleTable::get <<
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 sqr <<
1133 tamp2=piMinuspIne(particle1, particle2) <<
1134 else <<
1135 tamp2=9.04*std::pow(p1, -1.17)+18.*std: <<
1136 if (tamp2 < 0.0) tamp2=0; <<
1137 <<
1138 // CAS PI- P ET PI+ N <<
1139 return tamp2; <<
1140 } <<
1141 <<
1142 G4double CrossSectionsMultiPions::piPluspTw <<
1143 // <<
1144 // pion-nucleon interaction, producin <<
1145 // fit from Landolt-Bornstein multipl <<
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::mo <<
1160 <<
1161 // this limit corresponds to sqrt(s)=1230 <<
1162 if(pLab<296.367) <<
1163 return 0.0; <<
1164 <<
1165 // const G4int ipi = ParticleTable::getI <<
1166 // const G4int ind2 = ParticleTable::get <<
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( <<
1175 tamp6=piPluspIne(particle1, particle2)- <<
1176 else <<
1177 tamp6=1.59+25.5*std::pow(p1, -1.04); // <<
1178 <<
1179 // CAS PI+ P ET PI- N <<
1180 return tamp6; <<
1181 } <<
1182 <<
1183 G4double CrossSectionsMultiPions::piMinus <<
1184 // <<
1185 // pion-nucleon interaction, producing <<
1186 // fit from Landolt-Bornstein multiplie <<
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()); 870 // assert(pion->isPion());
1199 <<
1200 const G4double pLab = KinematicsUtils::mome <<
1201 <<
1202 // this limit corresponds to sqrt(s)=1230 M <<
1203 if(pLab<296.367) <<
1204 return 0.0; <<
1205 <<
1206 // const G4int ipi = ParticleTable::getIso <<
1207 // const G4int ind2 = ParticleTable::getIs <<
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)= <<
1216 tamp2=piMinuspIne(particle1, particle2)-p <<
1217 else <<
1218 tamp2=2.457794117647+18.066176470588*std: <<
1219 <<
1220 // CAS PI- P ET PI+ N <<
1221 return tamp2; <<
1222 } <<
1223 <<
1224 <<
1225 <<
1226 <<
1227 G4double CrossSectionsMultiPions::piNToEt <<
1228 // <<
1229 // Pion-Nucleon producing Eta cross s <<
1230 // <<
1231 return 0.; <<
1232 } <<
1233 <<
1234 G4double CrossSectionsMultiPions::piNToOm <<
1235 // <<
1236 // Pion-Nucleon producing Omega cross <<
1237 // <<
1238 return 0.; <<
1239 } <<
1240 <<
1241 G4double CrossSectionsMultiPions::piNToEt <<
1242 // <<
1243 // Pion-Nucleon producing EtaPrime cr <<
1244 // <<
1245 return 0.; <<
1246 } <<
1247 <<
1248 G4double CrossSectionsMultiPions::etaNToP <<
1249 // <<
1250 // Eta-Nucleon producing Pion cross s <<
1251 // <<
1252 return 0.; <<
1253 } <<
1254 <<
1255 <<
1256 G4double CrossSectionsMultiPions::etaNTo <<
1257 // <<
1258 // Eta-Nucleon producing Two Pions cr <<
1259 // <<
1260 return 0.; <<
1261 } <<
1262 <<
1263 <<
1264 G4double CrossSectionsMultiPions::omegaNT <<
1265 // <<
1266 // Omega-Nucleon producing Pion cross <<
1267 // <<
1268 return 0.; <<
1269 } <<
1270 <<
1271 G4double CrossSectionsMultiPions::omegaNT <<
1272 // <<
1273 // Omega-Nucleon producing Two Pions <<
1274 // <<
1275 return 0.; <<
1276 } <<
1277 <<
1278 G4double CrossSectionsMultiPions::etaPrim <<
1279 // <<
1280 // EtaPrime-Nucleon producing Pion cr <<
1281 // <<
1282 return 0.; <<
1283 } <<
1284 <<
1285 G4double CrossSectionsMultiPions::NNToNNE <<
1286 // <<
1287 // Nucleon-Nucleon producing Eta cros <<
1288 // <<
1289 return 0.; <<
1290 } <<
1291 <<
1292 G4double CrossSectionsMultiPions::NNToNNEta <<
1293 // <<
1294 // Nucleon-Nucleon producing Eta cros <<
1295 // <<
1296 return 0.; <<
1297 } <<
1298 <<
1299 G4double CrossSectionsMultiPions::NNToNNEta <<
1300 return 0.; <<
1301 } <<
1302 <<
1303 G4double CrossSectionsMultiPions::NNToNDe <<
1304 // <<
1305 // Nucleon-Nucleon producing N-Delta- <<
1306 // <<
1307 return 0.; <<
1308 } <<
1309 <<
1310 G4double CrossSectionsMultiPions::NNToNNO <<
1311 // <<
1312 // Nucleon-Nucleon producing Omega cr <<
1313 // <<
1314 return 0.; <<
1315 } <<
1316 <<
1317 G4double CrossSectionsMultiPions::NNToNNO <<
1318 // <<
1319 // Nucleon-Nucleon producing Omega cr <<
1320 // <<
1321 return 0.; <<
1322 } <<
1323 <<
1324 G4double CrossSectionsMultiPions::NNToNNO <<
1325 return 0.; <<
1326 } <<
1327 <<
1328 G4double CrossSectionsMultiPions::NNToNDe <<
1329 // <<
1330 // Nucleon-Nucleon producing N-Delta-Om <<
1331 // <<
1332 return 0.; <<
1333 } <<
1334 <<
1335 <<
1336 <<
1337 <<
1338 G4double CrossSectionsMultiPions::NYelast <<
1339 // <<
1340 // Hyperon-Nucleon elastic cross <<
1341 // <<
1342 return 0.; <<
1343 } <<
1344 <<
1345 G4double CrossSectionsMultiPions::NKelast <<
1346 // <<
1347 // Kaon-Nucleon elastic cross se <<
1348 // <<
1349 return 0.; <<
1350 } <<
1351 <<
1352 G4double CrossSectionsMultiPions::NKbelas <<
1353 // <<
1354 // antiKaon-Nucleon elastic cros <<
1355 // <<
1356 return 0.; <<
1357 } <<
1358 <<
1359 <<
1360 G4double CrossSectionsMultiPions::NNToNLK(P <<
1361 // <<
1362 // Nucleon-Nucleon producing N-L <<
1363 // <<
1364 return 0.; <<
1365 } <<
1366 <<
1367 G4double CrossSectionsMultiPions::NNToNSK <<
1368 // <<
1369 // Nucleon-Nucleon producing N-S <<
1370 // <<
1371 return 0.; <<
1372 } <<
1373 <<
1374 G4double CrossSectionsMultiPions::NNToNLK <<
1375 // <<
1376 // Nucleon-Nucleon producing N-L <<
1377 // <<
1378 return 0.; <<
1379 } <<
1380 <<
1381 G4double CrossSectionsMultiPions::NNToNSK <<
1382 // <<
1383 // Nucleon-Nucleon producing N-S <<
1384 // <<
1385 return 0.; <<
1386 } <<
1387 <<
1388 G4double CrossSectionsMultiPions::NNToNLK <<
1389 // <<
1390 // Nucleon-Nucleon producing N-La <<
1391 // <<
1392 return 0.; <<
1393 } <<
1394 <<
1395 G4double CrossSectionsMultiPions::NNToNSK <<
1396 // <<
1397 // Nucleon-Nucleon producing N-S <<
1398 // <<
1399 return 0.; <<
1400 } <<
1401 <<
1402 G4double CrossSectionsMultiPions::NNToNNK <<
1403 // <<
1404 // Nucleon-Nucleon producing Nuc <<
1405 // <<
1406 return 0.; <<
1407 } <<
1408 <<
1409 G4double CrossSectionsMultiPions::NNToMis <<
1410 // <<
1411 // Nucleon-Nucleon missing stran <<
1412 // <<
1413 return 0.; <<
1414 } <<
1415 <<
1416 G4double CrossSectionsMultiPions::NDeltaT <<
1417 // Nucleon-Delta producing Nucleon La <<
1418 return 0; <<
1419 } <<
1420 G4double CrossSectionsMultiPions::NDeltaT <<
1421 // Nucleon-Delta producing Nucleon Si <<
1422 return 0; <<
1423 } <<
1424 G4double CrossSectionsMultiPions::NDeltaT <<
1425 // Nucleon-Delta producing Delta Lamb <<
1426 return 0; <<
1427 } <<
1428 G4double CrossSectionsMultiPions::NDeltaT <<
1429 // Nucleon-Delta producing Delta Sigm <<
1430 return 0; <<
1431 } <<
1432 <<
1433 G4double CrossSectionsMultiPions::NDeltaT <<
1434 // Nucleon-Delta producing Nucleon-Nu <<
1435 return 0; <<
1436 } <<
1437 <<
1438 <<
1439 G4double CrossSectionsMultiPions::NpiToLK <<
1440 // <<
1441 // Pion-Nucleon producing Lambda <<
1442 // <<
1443 return 0.; <<
1444 } <<
1445 <<
1446 G4double CrossSectionsMultiPions::NpiToSK <<
1447 // <<
1448 // Pion-Nucleon producing Sigma- <<
1449 // <<
1450 return 0.; <<
1451 } <<
1452 G4double CrossSectionsMultiPions::p_pimTo <<
1453 return 0.; <<
1454 } <<
1455 G4double CrossSectionsMultiPions::p_pimTo <<
1456 return 0.; <<
1457 } <<
1458 G4double CrossSectionsMultiPions::p_pizTo <<
1459 return 0.; <<
1460 } <<
1461 <<
1462 G4double CrossSectionsMultiPions::NpiToLK <<
1463 // <<
1464 // Pion-Nucleon producing Lambda <<
1465 // <<
1466 return 0.; <<
1467 } <<
1468 <<
1469 G4double CrossSectionsMultiPions::NpiToSK <<
1470 // <<
1471 // Pion-Nucleon producing Sigma- <<
1472 // <<
1473 return 0.; <<
1474 } <<
1475 <<
1476 G4double CrossSectionsMultiPions::NpiToLK <<
1477 // <<
1478 // Pion-Nucleon producing Lambda <<
1479 // <<
1480 return 0.; <<
1481 } <<
1482 <<
1483 G4double CrossSectionsMultiPions::NpiToSK <<
1484 // <<
1485 // Pion-Nucleon producing Lambda <<
1486 // <<
1487 return 0.; <<
1488 } <<
1489 <<
1490 G4double CrossSectionsMultiPions::NpiToNK <<
1491 // <<
1492 // Pion-Nucleon producing Nucleo <<
1493 // <<
1494 return 0.; <<
1495 } <<
1496 <<
1497 G4double CrossSectionsMultiPions::NpiToMi <<
1498 // <<
1499 // Pion-Nucleon missing strangen <<
1500 // <<
1501 return 0.; <<
1502 } <<
1503 <<
1504 G4double CrossSectionsMultiPions::NLToNS( <<
1505 // <<
1506 // Nucleon-Hyperon multiplet cha <<
1507 // <<
1508 return 0.; <<
1509 } <<
1510 <<
1511 G4double CrossSectionsMultiPions::NSToNL( <<
1512 // <<
1513 // Nucleon-Sigma quasi-elastic c <<
1514 // <<
1515 return 0.; <<
1516 } <<
1517 <<
1518 G4double CrossSectionsMultiPions::NSToNS( <<
1519 // <<
1520 // Nucleon-Sigma quasi-elastic c <<
1521 // <<
1522 return 0.; <<
1523 } <<
1524 871
1525 G4double CrossSectionsMultiPions::NKToNK( << 872 const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1526 // <<
1527 // Nucleon-Kaon quasi-elastic cr <<
1528 // <<
1529 return 0.; <<
1530 } <<
1531 873
1532 G4double CrossSectionsMultiPions::NKToNKp << 874 // this limit corresponds to sqrt(s)=1230 MeV
1533 // << 875 if(pLab<296.367)
1534 // Nucleon-Kaon producing Nucleo << 876 return 0.0;
1535 // << 877
1536 return 0.; << 878 const G4int ipi = ParticleTable::getIsospin(pion->getType());
1537 } << 879 const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
1538 << 880 const G4int cg = 4 + ind2*ipi;
1539 G4double CrossSectionsMultiPions::NKToNK2 << 881 // assert(cg==2 || cg==4 || cg==6);
1540 // <<
1541 // Nucleon-Kaon producing Nucleo <<
1542 // <<
1543 return 0.; <<
1544 } <<
1545 882
1546 G4double CrossSectionsMultiPions::NKbToNK << 883 const G4double p1=1e-3*pLab;
1547 // << 884 G4double tamp6=0.;
1548 // Nucleon-antiKaon quasi-elasti << 885 G4double tamp2=0.;
1549 // << 886
1550 return 0.; << 887 // X-SECTION PI+ P INELASTIQUE :
1551 } << 888 if(cg != 2) {
>> 889 if(pLab < 1532.52) // corresponds to sqrt(s)=1946 MeV
>> 890 tamp6=piNIne(particle1, particle2);
>> 891 else
>> 892 tamp6=0.204+18.2*std::pow(p1, -1.72)+6.33*std::pow(p1, -1.13);
>> 893 if (cg == 6) // CAS PI+ P ET PI- N
>> 894 return tamp6;
>> 895 }
>> 896
>> 897 // X-SECTION PI- P INELASTIQUE :
>> 898 if (pLab < 1228.06) // corresponds to sqrt(s)=1794 MeV
>> 899 tamp2=piNIne(particle1, particle2);
>> 900 else
>> 901 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.;
>> 902 if (tamp2 < 0.0) tamp2=0;
1552 903
1553 G4double CrossSectionsMultiPions::NKbToSp << 904 if (cg == 2) // CAS PI- P ET PI+ N
1554 // << 905 return tamp2;
1555 // Nucleon-antiKaon producing Si << 906 else {
1556 // << 907 // CAS PI0 P ET PI0 N
1557 return 0.; << 908 G4double s1pin = 0.5*(tamp6+tamp2);
>> 909 const G4double inelastic = piNIne(particle1, particle2);
>> 910 if (s1pin > inelastic)
>> 911 s1pin = inelastic;
>> 912 return s1pin;
>> 913 }
1558 } 914 }
1559 915
1560 G4double CrossSectionsMultiPions::NKbToLp << 916 G4double CrossSectionsMultiPions::piNTwoPi(Particle const * const particle1, Particle const * const particle2) {
1561 // << 917 //
1562 // Nucleon-antiKaon producing La << 918 // pion-nucleon interaction, producing 2 pions
1563 // << 919 // fit from Landolt-Bornstein multiplied by factor determined with evaluation of total xs
1564 return 0.; << 920 //
1565 } <<
1566 <<
1567 G4double CrossSectionsMultiPions::NKbToS2 <<
1568 // <<
1569 // Nucleon-antiKaon producing Si <<
1570 // <<
1571 return 0.; <<
1572 } <<
1573 921
1574 G4double CrossSectionsMultiPions::NKbToL2 << 922 const Particle *pion;
1575 // << 923 const Particle *nucleon;
1576 // Nucleon-antiKaon producing La << 924 if(particle1->isNucleon()) {
1577 // << 925 nucleon = particle1;
1578 return 0.; << 926 pion = particle2;
1579 } << 927 } else {
>> 928 pion = particle1;
>> 929 nucleon = particle2;
>> 930 }
>> 931 // assert(pion->isPion());
1580 932
1581 G4double CrossSectionsMultiPions::NKbToNK << 933 const G4double pLab = KinematicsUtils::momentumInLab(pion, nucleon);
1582 // <<
1583 // Nucleon-antiKaon producing Nu <<
1584 // <<
1585 return 0.; <<
1586 } <<
1587 934
1588 G4double CrossSectionsMultiPions::NKbToNK << 935 // this limit corresponds to sqrt(s)=1230 MeV
1589 // << 936 if(pLab<296.367)
1590 // Nucleon-antiKaon producing Nu << 937 return 0.0;
1591 // << 938
1592 return 0.; << 939 const G4int ipi = ParticleTable::getIsospin(pion->getType());
1593 } << 940 const G4int ind2 = ParticleTable::getIsospin(nucleon->getType());
1594 << 941 const G4int cg = 4 + ind2*ipi;
1595 G4double CrossSectionsMultiPions::NNbarEl << 942 // assert(cg==2 || cg==4 || cg==6);
1596 // <<
1597 // Nucleon-AntiNucleon to Nucleon <<
1598 // <<
1599 return 0.; <<
1600 } <<
1601 943
1602 G4double CrossSectionsMultiPions::NNbarCE << 944 const G4double p1=1e-3*pLab;
1603 // << 945 G4double tamp6=0.;
1604 // Nucleon-AntiNucleon charge exc << 946 G4double tamp2=0.;
1605 // << 947
1606 return 0.; << 948 // X-SECTION PI+ P INELASTIQUE :
1607 } << 949 if(cg!=2) {
>> 950 if(pLab < 2444.7) // corresponds to sqrt(s)=2344 MeV
>> 951 tamp6=piNIne(particle1, particle2)-piNOnePi(particle1, particle2);
>> 952 else
>> 953 tamp6=1.59+25.5*std::pow(p1, -1.04); // tamp6=(0.636+10.2*std::pow(p1, -1.04))*15./6.;
>> 954
>> 955 if(cg==6) // CAS PI+ P ET PI- N
>> 956 return tamp6;
>> 957 }
>> 958
>> 959 // X-SECTION PI- P INELASTIQUE :
>> 960 if(pLab<2083.63) // corresponds to sqrt(s)=2195 MeV
>> 961 tamp2=piNIne(particle1, particle2)-piNOnePi(particle1, particle2);
>> 962 else
>> 963 tamp2=2.457794117647+18.066176470588*std::pow(p1, -0.92); // tamp2=(0.619+4.55*std::pow(p1, -0.92))*135./34.;
1608 964
1609 G4double CrossSectionsMultiPions::NNbarTo << 965 if(cg==2) // CAS PI- P ET PI+ N
1610 // << 966 return tamp2;
1611 // Nucleon-AntiNucleon to Lambda- << 967 else {
1612 // << 968 // CAS PI0 P ET PI0 N
1613 return 0.; << 969 const G4double s2pin=0.5*(tamp6+tamp2);
1614 } << 970 return s2pin;
1615 << 971 }
1616 G4double CrossSectionsMultiPions::NNbarTo <<
1617 // <<
1618 // Nucleon-AntiNucleon to Nucleon <<
1619 // <<
1620 return 0.; <<
1621 } 972 }
1622 973
1623 G4double CrossSectionsMultiPions::NNbarTo <<
1624 // <<
1625 // Nucleon-AntiNucleon to Nucleon <<
1626 // <<
1627 return 0.; <<
1628 } <<
1629 <<
1630 G4double CrossSectionsMultiPions::NNbarTo <<
1631 // <<
1632 // Nucleon-AntiNucleon to Nucleon <<
1633 // <<
1634 return 0.; <<
1635 } <<
1636 <<
1637 G4double CrossSectionsMultiPions::NNbarTo <<
1638 // <<
1639 // Nucleon-AntiNucleon total anni <<
1640 // <<
1641 return 0.; <<
1642 } <<
1643 } // namespace G4INCL 974 } // namespace G4INCL
1644 975
1645 976