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