| 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. (?)
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. (?)
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 return sigma;
300 return sigma; << 300 }
301 } <<
302 301
303 G4double CrossSectionsMultiPionsAndResonan << 302 G4double CrossSectionsMultiPionsAndResonances::etaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
304 // << 303 //
305 // Eta-Nucleon producing Two Pions << 304 // Eta-Nucleon producing Two Pions cross sections
306 // << 305 //
307 // assert((particle1->isNucleon() && particle2 306 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
308 << 307
309 G4double sigma=0.; << 308 G4double sigma=0.;
310 << 309
311 const Particle *eta; << 310 const Particle *eta;
312 const Particle *nucleon; << 311 const Particle *nucleon;
313 << 312
314 if (particle1->isEta()) { << 313 if (particle1->isEta()) {
315 eta = particle1; << 314 eta = particle1;
316 nucleon = particle2; << 315 nucleon = particle2;
317 } << 316 }
318 else { << 317 else {
319 eta = particle2; << 318 eta = particle2;
320 nucleon = particle1; << 319 nucleon = particle1;
321 } << 320 }
322 << 321
323 const G4double pLab = KinematicsUtils: << 322 const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
324 << 323
325 if (pLab < 450.) << 324 if (pLab < 450.)
326 sigma = 2.01854221E-13*std::pow(pL << 325 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.) << 326 else if (pLab < 600.)
328 sigma = 2.01854221E-13*std::pow(45 << 327 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.) << 328 else if (pLab <= 1300.)
330 sigma = -6.32793049e-16*std::pow(p << 329 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) - << 330 1.37055547e-05*std::pow(pLab,3) - 1.01830486e-02*std::pow(pLab,2) + 3.93492126*pLab - 609.447145;
332 else << 331 else
333 sigma = etaNToPiN(particle1,partic << 332 sigma = etaNToPiN(particle1,particle2);
334 << 333
335 if (sigma < 0.) sigma = 0.; << 334 if (sigma < 0.) sigma = 0.;
336 return sigma; // Parameterization from << 335 return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209] - eta p --> "pi+pi0 n" + "pi0 pi0 p" total XS
337 } << 336 }
338 << 337
339 << 338
340 G4double CrossSectionsMultiPionsAndResonan << 339 G4double CrossSectionsMultiPionsAndResonances::etaNElastic(Particle const * const particle1, Particle const * const particle2) {
341 // << 340 //
342 // Eta-Nucleon elastic cross secti << 341 // Eta-Nucleon elastic cross sections
343 // << 342 //
344 // assert((particle1->isNucleon() && particle2 343 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
345 << 344
346 G4double sigma=0.; << 345 G4double sigma=0.;
347 << 346
348 const Particle *eta; << 347 const Particle *eta;
349 const Particle *nucleon; << 348 const Particle *nucleon;
350 << 349
351 if (particle1->isEta()) { << 350 if (particle1->isEta()) {
352 eta = particle1; << 351 eta = particle1;
353 nucleon = particle2; << 352 nucleon = particle2;
354 } << 353 }
355 else { << 354 else {
356 eta = particle2; << 355 eta = particle2;
357 nucleon = particle1; << 356 nucleon = particle1;
358 } << 357 }
359 << 358
360 const G4double pLab = KinematicsUtils: << 359 const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
361 << 360
362 if (pLab < 700.) << 361 if (pLab < 700.)
363 sigma = 3.6838e-15*std::pow(pLab,6 << 362 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 << 363 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.) << 364 else if (pLab < 1400.)
366 sigma = 3.562630e-16*std::pow(pLab << 365 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. << 366 9.667078e-06*std::pow(pLab,3) + 7.894845e-03*std::pow(pLab,2) - 3.409200*pLab + 609.8501;
368 else if (pLab < 2025.) << 367 else if (pLab < 2025.)
369 sigma = -1.041950E-03*pLab + 2.110 << 368 sigma = -1.041950E-03*pLab + 2.110529E+00;
370 else << 369 else
371 sigma=0.; << 370 sigma=0.;
372 << 371
373 if (sigma < 0.) sigma = 0.; << 372 if (sigma < 0.) sigma = 0.;
374 return sigma; // Parameterization << 373 return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209]
375 } << 374 }
376 << 375
377 G4double CrossSectionsMultiPionsAndResonan << 376 G4double CrossSectionsMultiPionsAndResonances::omegaNInelastic(Particle const * const particle1, Particle const * const particle2) {
378 // << 377 //
379 // Omega-Nucleon inelastic cross s << 378 // Omega-Nucleon inelastic cross sections
380 // << 379 //
381 // assert((particle1->isNucleon() && particle2 380 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
382 << 381
383 G4double sigma=0.; << 382 G4double sigma=0.;
384 << 383
385 const Particle *omega; << 384 const Particle *omega;
386 const Particle *nucleon; << 385 const Particle *nucleon;
387 << 386
388 if (particle1->isOmega()) { << 387 if (particle1->isOmega()) {
389 omega = particle1; << 388 omega = particle1;
390 nucleon = particle2; << 389 nucleon = particle2;
391 } << 390 }
392 else { << 391 else {
393 omega = particle2; << 392 omega = particle2;
394 nucleon = particle1; << 393 nucleon = particle1;
395 } << 394 }
396 << 395
397 const G4double pLab = KinematicsUtils: << 396 const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
398 << 397
399 sigma = 20. + 4.0/pLab; // Eq.(24) << 398 sigma = 20. + 4.0/pLab; // Eq.(24) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
400 << 399
401 return sigma; << 400 return sigma;
402 } << 401 }
403 402
404 403
405 G4double CrossSectionsMultiPionsAndResonan << 404 G4double CrossSectionsMultiPionsAndResonances::omegaNElastic(Particle const * const particle1, Particle const * const particle2) {
406 // << 405 //
407 // Omega-Nucleon elastic cross sec << 406 // Omega-Nucleon elastic cross sections
408 // << 407 //
409 // assert((particle1->isNucleon() && particle2 408 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
410 << 409
411 G4double sigma=0.; << 410 G4double sigma=0.;
412 << 411
413 const Particle *omega; << 412 const Particle *omega;
414 const Particle *nucleon; << 413 const Particle *nucleon;
415 << 414
416 if (particle1->isOmega()) { << 415 if (particle1->isOmega()) {
417 omega = particle1; << 416 omega = particle1;
418 nucleon = particle2; << 417 nucleon = particle2;
419 } << 418 }
420 else { << 419 else {
421 omega = particle2; << 420 omega = particle2;
422 nucleon = particle1; << 421 nucleon = particle1;
423 } << 422 }
424 << 423
425 const G4double pLab = KinematicsUtils: << 424 const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
426 << 425
427 sigma = 5.4 + 10.*std::exp(-0.6*pLab); << 426 sigma = 5.4 + 10.*std::exp(-0.6*pLab); // Eq.(21) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
428 << 427
429 return sigma; << 428 return sigma;
430 } << 429 }
431 430
432 << 431
433 G4double CrossSectionsMultiPionsAndResonan << 432 G4double CrossSectionsMultiPionsAndResonances::omegaNToPiN(Particle const * const particle1, Particle const * const particle2) {
434 // << 433 //
435 // Omega-Nucleon producing Pion cr << 434 // Omega-Nucleon producing Pion cross sections
436 // << 435 //
437 // assert((particle1->isNucleon() && particle2 436 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
438 << 437
439 G4double ECM=KinematicsUtils::totalEne << 438 G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
440 << 439
441 G4double massPiZero=ParticleTable::get << 440 G4double massPiZero=ParticleTable::getINCLMass(PiZero);
442 G4double massPiMinus=ParticleTable::ge << 441 G4double massPiMinus=ParticleTable::getINCLMass(PiMinus);
443 G4double massProton=ParticleTable::get << 442 G4double massProton=ParticleTable::getINCLMass(Proton);
444 << 443
445 G4double massomega; << 444 G4double massomega;
446 G4double massnucleon; << 445 G4double massnucleon;
447 G4double pCM_omega; << 446 G4double pCM_omega;
448 G4double pLab_omega; << 447 G4double pLab_omega;
449 << 448
450 G4double sigma=0.; << 449 G4double sigma=0.;
451 << 450
452 if (particle1->isOmega()) { << 451 if (particle1->isOmega()) {
453 massomega=particle1->getMass(); << 452 massomega=particle1->getMass();
454 massnucleon=particle2->getMass(); << 453 massnucleon=particle2->getMass();
455 } << 454 }
456 else { << 455 else {
457 massomega=particle2->getMass(); << 456 massomega=particle2->getMass();
458 massnucleon=particle1->getMass(); << 457 massnucleon=particle1->getMass();
459 } << 458 }
460 pCM_omega=KinematicsUtils::momentumInC << 459 pCM_omega=KinematicsUtils::momentumInCM(ECM, massomega, massnucleon);
461 pLab_omega=KinematicsUtils::momentumIn << 460 pLab_omega=KinematicsUtils::momentumInLab(ECM*ECM, massomega, massnucleon);
462 << 461
463 G4double pCM_PiZero=KinematicsUtils::m << 462 G4double pCM_PiZero=KinematicsUtils::momentumInCM(ECM, massPiZero, massProton);
464 G4double pCM_PiMinus=KinematicsUtils:: << 463 G4double pCM_PiMinus=KinematicsUtils::momentumInCM(ECM, massPiMinus, massProton); // = pCM_PiPlus (because massPiMinus = massPiPlus)
465 << 464
466 sigma = (piMinuspToOmegaN(ECM)/2.) * s << 465 sigma = (piMinuspToOmegaN(ECM)/2.) * std::pow((pCM_PiZero/pCM_omega), 2)
467 + piMinuspToOmegaN(ECM) * std::pow((pC << 466 + piMinuspToOmegaN(ECM) * std::pow((pCM_PiMinus/pCM_omega), 2);
468 467
469 if (sigma > omegaNInelastic(particle1, << 468 if (sigma > omegaNInelastic(particle1, particle2) || (pLab_omega < 200.)) {
470 // if (sigma > omegaNInelastic(particle << 469 // if (sigma > omegaNInelastic(particle1, particle2)) {
471 sigma = omegaNInelastic(particle1, << 470 sigma = omegaNInelastic(particle1, particle2);
472 } << 471 }
473 472
474 return sigma; << 473 return sigma;
475 } << 474 }
476 << 475
477 476
478 G4double CrossSectionsMultiPionsAndResonan << 477 G4double CrossSectionsMultiPionsAndResonances::omegaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
479 // << 478 //
480 // Omega-Nucleon producing 2 PionS << 479 // Omega-Nucleon producing 2 PionS cross sections
481 // << 480 //
482 // assert((particle1->isNucleon() && particle2 481 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
483 << 482
484 G4double sigma=0.; << 483 G4double sigma=0.;
485 << 484
486 sigma = omegaNInelastic(particle1,part << 485 sigma = omegaNInelastic(particle1,particle2) - omegaNToPiN(particle1,particle2) ;
487 << 486
488 return sigma; << 487 return sigma;
489 } << 488 }
490 489
491 << 490
492 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT 491 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
493 G4double CrossSectionsMultiPionsAndResonance 492 G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const /*particle1*/, Particle const * const /*particle2*/) {
494 #else 493 #else
495 G4double CrossSectionsMultiPionsAndResonance 494 G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const particle1, Particle const * const particle2) {
496 #endif 495 #endif
497 // << 496 //
498 // EtaPrime-Nucleon producing Pion << 497 // EtaPrime-Nucleon producing Pion cross sections
499 // << 498 //
500 // assert((particle1->isNucleon() && particle2 499 // assert((particle1->isNucleon() && particle2->isEtaPrime()) || (particle1->isEtaPrime() && particle2->isNucleon()));
501 << 500
502 return 0.; << 501 return 0.;
503 } << 502 }
504 << 503
505 G4double CrossSectionsMultiPionsAndResonan << 504 G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(Particle const * const particle1, Particle const * const particle2) {
506 // << 505 //
507 // Pion-Nucleon producing Eta cros << 506 // Pion-Nucleon producing Eta cross sections
508 // << 507 //
509 // assert((particle1->isNucleon() && particle2 508 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
510 << 509
511 G4double masspion; << 510 G4double masspion;
512 G4double massnucleon; << 511 G4double massnucleon;
513 if (particle1->isPion()) { << 512 if (particle1->isPion()) {
514 masspion=particle1->getMass(); << 513 masspion=particle1->getMass();
515 massnucleon=particle2->getMass(); << 514 massnucleon=particle2->getMass();
516 } << 515 }
517 else { << 516 else {
518 masspion=particle2->getMass(); << 517 masspion=particle2->getMass();
519 massnucleon=particle1->getMass(); << 518 massnucleon=particle1->getMass();
520 } << 519 }
521 << 520
522 G4double ECM=KinematicsUtils::totalEne << 521 G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
523 G4double plab=KinematicsUtils::momentu << 522 G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
524 << 523
525 G4double sigma; << 524 G4double sigma;
526 << 525
527 // new parameterization (JCD) - end of februar 526 // new parameterization (JCD) - end of february 2016
528 if (ECM < 1486.5) sigma=0.; << 527 if (ECM < 1486.5) sigma=0.;
529 else << 528 else
530 { << 529 {
531 if (ECM < 1535.) << 530 if (ECM < 1535.)
532 { << 531 {
533 sigma = -0.0000003689197974814 << 532 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 } << 533 }
535 else if (ECM < 1670.) << 534 else if (ECM < 1670.)
536 { << 535 {
537 sigma = -0.0000000337986446*st << 536 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 } << 537 }
539 else if (ECM < 1714.) << 538 else if (ECM < 1714.)
540 { << 539 {
541 sigma = 0.000003737765*std::p << 540 sigma = 0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
542 } << 541 }
543 else sigma=1.47*std::pow(plab, -1. << 542 else sigma=1.47*std::pow(plab, -1.68);
544 } << 543 }
545 // << 544 //
546 << 545
547 return sigma; << 546 return sigma;
548 } << 547 }
549 << 548
550 G4double CrossSectionsMultiPionsAndResonan << 549 G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(const G4double ECM) {
551 // << 550 //
552 // Pion-Nucleon producing Eta cros << 551 // Pion-Nucleon producing Eta cross sections
553 // << 552 //
554 << 553
555 const G4double masspion = ParticleTabl << 554 const G4double masspion = ParticleTable::getRealMass(PiMinus);
556 const G4double massnucleon = ParticleT << 555 const G4double massnucleon = ParticleTable::getRealMass(Proton);
557 556
558 G4double plab=KinematicsUtils::momentu << 557 G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
559 << 558
560 G4double sigma; << 559 G4double sigma;
561 << 560
562 // new parameterization (JCD) - end of februar 561 // new parameterization (JCD) - end of february 2016
563 if (ECM < 1486.5) sigma=0.; << 562 if (ECM < 1486.5) sigma=0.;
564 else << 563 else
565 { << 564 {
566 if (ECM < 1535.) << 565 if (ECM < 1535.)
567 { << 566 {
568 sigma = -0.0000003689197974814 << 567 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 } << 568 }
570 else if (ECM < 1670.) << 569 else if (ECM < 1670.)
571 { << 570 {
572 sigma = -0.0000000337986446*st << 571 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 } << 572 }
574 else if (ECM < 1714.) << 573 else if (ECM < 1714.)
575 { << 574 {
576 sigma = 0.000003737765*std::p << 575 sigma = 0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
577 } << 576 }
578 else sigma=1.47*std::pow(plab, -1. << 577 else sigma=1.47*std::pow(plab, -1.68);
579 } << 578 }
580 << 579 //
581 return sigma; << 580
582 } << 581 return sigma;
583 << 582 }
584 G4double CrossSectionsMultiPionsAndResonan << 583
585 // << 584 G4double CrossSectionsMultiPionsAndResonances::piMinuspToOmegaN(Particle const * const particle1, Particle const * const particle2) {
586 // Pion-Nucleon producing Omega cr << 585 //
587 // << 586 // Pion-Nucleon producing Omega cross sections
>> 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
>> 640 const G4double Ecm=0.001*ener;
>> 641 G4double sNNEta; // pp->pp+eta(+X)
>> 642 G4double sNNEta1; // np->np+eta(+X)
>> 643 G4double sNNEta2; // np->d+eta (d wil be considered as np - How far is this right?)
>> 644 G4double x=Ecm*Ecm/5.88;
>> 645
>> 646 if (Ecm > 5.5) {
>> 647 // if (Ecm > 4) {
>> 648 sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.;
>> 649 }
>> 650 else if(Ecm>2.70555) {
>> 651 sNNEta = 167.68*Ecm*Ecm-528.3*Ecm+442.29; // Mimicking the reduction seen in omega by HADES compared to Sibirtsev (5.5)
>> 652 // sNNEta = 710.1*Ecm*Ecm-3719.7*Ecm+5106.2; // Mimicking the reduction seen in omega by HADES compared to Sibirtsev (4)
>> 653 // sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.*std::exp(-std::pow(Ecm-2.70555,0.2)*0.3020904); //exp to reduce from 1 to 0.7 (2.07555 --> 5)
>> 654 if (sNNEta <= NNToNNEtaExcluIso(ener, 2)*1000.) sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
>> 655 }
>> 656 /* if(Ecm>2.70555) {
>> 657 sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.;
>> 658 // sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.*std::exp(-std::pow(Ecm-2.70555,0.2)*0.3020904); //exp to reduce from 1 to 0.7 (2.07555 --> 5)
>> 659 if (sNNEta <= NNToNNEtaExcluIso(ener, 2)*1000.) sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
>> 660 }*/
>> 661 else {
>> 662 sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
>> 663 }
>> 664
>> 665 if (sNNEta < 1.e-9) sNNEta = 0.;
>> 666
>> 667 if (iso != 0) {
>> 668 return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
>> 669 }
>> 670
>> 671 if(Ecm>=5.5) {
>> 672 // if(Ecm>=4.) {
>> 673 // if(Ecm>=5.) {
>> 674 // sNNEta1 = 3*sNNEta;
>> 675 sNNEta1 = sNNEta;
>> 676 }
>> 677 else if (Ecm>=2.70555) {
>> 678 // sNNEta1 = 2349.400647*Ecm - 4794.192041; // with factor 3
>> 679 // sNNEta1 = 329.2183*Ecm + 671.5123; // with factor 1
>> 680 // sNNEta1 = 26.19086*Ecm + 1491.368; // with factor 1 and pp reduced
>> 681 // sNNEta1 = sNNEta*(-4.23*Ecm+17.92); // going from a factor 6.5 (low) to 1 (high) (4)
>> 682 sNNEta1 = sNNEta*(-1.968187*Ecm+11.82503); // going from a factor 6.5 (low) to 1 (high) (5.5)
>> 683 }
>> 684 else {
>> 685 sNNEta1 = 6.5*sNNEta; // 6.5: ratio pn/pp
>> 686 }
>> 687 // sNNEta1 = 6.5*sNNEta; // 6.5: ratio pn/pp
>> 688
>> 689 sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
>> 690 if (sNNEta2 < 0.) sNNEta2=0.;
>> 691
>> 692 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
>> 693
>> 694 G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
>> 695 G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
>> 696 G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
>> 697 if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta) sNNEta = 0.;
639 698
640 const G4double Ecm=0.001*ener; << 699 return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
641 G4double sNNEta; // pp->pp+eta(+X) << 700 }
642 G4double sNNEta1; // np->np+eta(+X) <<
643 G4double sNNEta2; // np->d+eta (d will <<
644 G4double x=Ecm*Ecm/5.88; <<
645 <<
646 //jcd <<
647 if (Ecm >= 3.05) { <<
648 sNNEta = 2.5*std::pow((x-1.),1.47) <<
649 } <<
650 else if(Ecm >= 2.6) { <<
651 sNNEta = -327.29*Ecm*Ecm*Ecm + 287 <<
652 if (sNNEta <= NNToNNEtaExcluIso(en <<
653 } <<
654 else { <<
655 sNNEta = NNToNNEtaExcluIso(ener, 2 <<
656 } <<
657 //jcd <<
658 if (sNNEta < 1.e-9) sNNEta = 0.; <<
659 <<
660 if (iso != 0) { <<
661 return sNNEta/1000.; // parameteri <<
662 } <<
663 <<
664 if(Ecm >= 6.25) { <<
665 sNNEta1 = sNNEta; <<
666 } <<
667 else if (Ecm >= 2.6) { <<
668 sNNEta1 = sNNEta*std::exp(-(-5.531 <<
669 } <<
670 else if (Ecm >= 2.525) { // = exclusiv <<
671 sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ec <<
672 } <<
673 else { // = exclusive pn <<
674 sNNEta1 = 17570.217219*Ecm*Ecm - 8 <<
675 } <<
676 <<
677 sNNEta2 = -10220.89518466*Ecm*Ecm+5122 <<
678 if (sNNEta2 < 0.) sNNEta2=0.; <<
679 <<
680 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta; <<
681 <<
682 G4double Mn=ParticleTable::getRealMass <<
683 G4double Mp=ParticleTable::getRealMass <<
684 G4double Meta=ParticleTable::getRealMa <<
685 if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta <<
686 701
687 return sNNEta/1000.; // parameterizati << 702
>> 703 G4double CrossSectionsMultiPionsAndResonances::NNToNNEta(Particle const * const particle1, Particle const * const particle2) {
>> 704
>> 705 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
>> 706 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
>> 707
>> 708 if (iso != 0) {
>> 709 return NNToNNEtaIso(ener, iso);
>> 710 }
>> 711 else {
>> 712 return 0.5*(NNToNNEtaIso(ener, 0)+NNToNNEtaIso(ener, 2));
>> 713 }
>> 714 }
>> 715
>> 716 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExcluIso(const G4double ener, const G4int iso) {
>> 717
>> 718 const G4double Ecm=0.001*ener;
>> 719 G4double sNNEta; // pp->pp+eta
>> 720 G4double sNNEta1; // np->np+eta
>> 721 G4double sNNEta2; // np->d+eta (d wil be considered as np - How far is this right?)
>> 722
>> 723 if(Ecm>=2.70555) { // By hand (JCD)
>> 724 sNNEta = -70.1*Ecm + 430.23;
>> 725 }
>> 726 else {
>> 727 sNNEta = -147043.497285*std::pow(Ecm,4) + 1487222.5438123*std::pow(Ecm,3) - 5634399.900744*std::pow(Ecm,2) + 9477290.199378*Ecm - 5972174.353438;
>> 728 }
>> 729
>> 730 G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
>> 731 G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
>> 732 G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
>> 733 G4double Thr0=0.;
>> 734 if (iso > 0) {
>> 735 Thr0=2.*Mp+Meta;
>> 736 }
>> 737 else if (iso < 0) {
>> 738 Thr0=2.*Mn+Meta;
>> 739 }
>> 740 else {
>> 741 Thr0=Mn+Mp+Meta;
>> 742 }
>> 743
>> 744 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.; // Thr0: Ecm threshold
>> 745
>> 746 if (iso != 0) {
>> 747 return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
>> 748 }
>> 749
>> 750 if(Ecm>=5.5) {
>> 751 // if(Ecm>=4.) {
>> 752 // if(Ecm>=5.) {
>> 753 // sNNEta1 = 3*sNNEta;
>> 754 sNNEta1 = sNNEta;
>> 755 }
>> 756 else if (Ecm>=2.70555) {
>> 757 // sNNEta1 = -577.2757*Ecm + 3125.5683;
>> 758 // sNNEta1 = -646.775*Ecm + 3313.605;
>> 759 // sNNEta1 = -646.775*Ecm + 3313.605;
>> 760 // sNNEta1 = sNNEta*(-4.23*Ecm+17.92); // going from a factor 6.5 (low) to 1 (high) (4.)
>> 761 sNNEta1 = sNNEta*(-1.968187*Ecm+11.82503); // going from a factor 6.5 (low) to 1 (high) (5.5)
688 } 762 }
689 << 763 else {
>> 764 sNNEta1 = 6.5*sNNEta; // 6.5: ratio pn/pp
>> 765 }
>> 766 // sNNEta1 = 6.5*sNNEta; // 6.5: ratio pn/pp
690 767
691 G4double CrossSectionsMultiPionsAndResonan << 768 sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
692 << 769 if (sNNEta2 < 0.) sNNEta2=0.;
693 const G4double ener=KinematicsUtils::t << 770
694 const G4int iso=ParticleTable::getIsos << 771 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
>> 772 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.; // Thr0: Ecm threshold
>> 773
>> 774 return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
>> 775 }
>> 776
>> 777 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExclu(Particle const * const particle1, Particle const * const particle2) {
>> 778
>> 779 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
>> 780 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
>> 781
>> 782 if (iso != 0) {
>> 783 return NNToNNEtaExcluIso(ener, iso);
>> 784 }
>> 785 else {
>> 786 return 0.5*(NNToNNEtaExcluIso(ener, 0)+NNToNNEtaExcluIso(ener, 2));
>> 787 }
>> 788 }
695 789
696 if (iso != 0) { << 790
697 return NNToNNEtaIso(ener, iso); << 791 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaIso(const G4double ener, const G4int iso) {
698 } << 792
699 else { << 793 const G4double Ecm=0.001*ener;
700 return 0.5*(NNToNNEtaIso(ener, 0)+ << 794 G4double sNNOmega; // pp->pp+eta(+X)
701 } << 795 G4double sNNOmega1; // np->np+eta(+X)
>> 796 G4double x=Ecm*Ecm/7.06;
>> 797
>> 798 if(Ecm>4.0) {
>> 799 sNNOmega = 2.5*std::pow(x-1, 1.47)*std::pow(x, -1.11) ;
702 } 800 }
703 << 801 else if(Ecm>2.802) { // 2802 MeV -> threshold to open inclusive (based on multipion threshold and omega mass)
704 G4double CrossSectionsMultiPionsAndResonan << 802 sNNOmega = (568.5254*Ecm*Ecm - 2694.045*Ecm + 3106.247)/1000.;
705 << 803 if (sNNOmega <= NNToNNOmegaExcluIso(ener, 2)) sNNOmega = NNToNNOmegaExcluIso(ener, 2);
706 const G4double Ecm=0.001*ener; <<
707 G4double sNNEta; // pp->pp+eta <<
708 G4double sNNEta1; // np->np+eta <<
709 G4double sNNEta2; // np->d+eta (d wil <<
710 <<
711 if(Ecm>=3.875) { // By hand (JCD) <<
712 sNNEta = -13.008*Ecm*Ecm + 84.531* <<
713 } <<
714 else if(Ecm>=2.725) { // By hand (JCD) <<
715 sNNEta = -913.2809*std::pow(Ecm,5) <<
716 } <<
717 else if(Ecm>=2.575) { // By hand (JCD) <<
718 sNNEta = -2640.3*Ecm*Ecm + 14692*E <<
719 } <<
720 else { <<
721 sNNEta = -147043.497285*std::pow(E <<
722 } <<
723 <<
724 G4double Mn=ParticleTable::getRealMass <<
725 G4double Mp=ParticleTable::getRealMass <<
726 G4double Meta=ParticleTable::getRealMa <<
727 G4double Thr0=0.; <<
728 if (iso > 0) { <<
729 Thr0=2.*Mp+Meta; <<
730 } <<
731 else if (iso < 0) { <<
732 Thr0=2.*Mn+Meta; <<
733 } <<
734 else { <<
735 Thr0=Mn+Mp+Meta; <<
736 } <<
737 <<
738 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNE <<
739 <<
740 if (iso != 0) { <<
741 return sNNEta/1000.; // parameteri <<
742 } <<
743 <<
744 if(Ecm>=3.9) { <<
745 sNNEta1 = sNNEta; <<
746 } <<
747 else if (Ecm >= 3.5) { <<
748 sNNEta1 = -1916.2*Ecm*Ecm*Ecm + 21 <<
749 } <<
750 else if (Ecm >= 2.525) { <<
751 sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ec <<
752 } <<
753 else { <<
754 sNNEta1 = 17570.217219*Ecm*Ecm - 8 <<
755 } <<
756 <<
757 sNNEta2 = -10220.89518466*Ecm*Ecm+5122 <<
758 if (sNNEta2 < 0.) sNNEta2=0.; <<
759 <<
760 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta; <<
761 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNE <<
762 <<
763 return sNNEta/1000.; // parameterizati <<
764 <<
765 } 804 }
766 << 805 else {
767 G4double CrossSectionsMultiPionsAndResonan << 806 sNNOmega = NNToNNOmegaExcluIso(ener, 2);
768 <<
769 const G4double ener=KinematicsUtils::t <<
770 const G4int iso=ParticleTable::getIsos <<
771 <<
772 if (iso != 0) { <<
773 return NNToNNEtaExcluIso(ener, iso <<
774 } <<
775 else { <<
776 return 0.5*(NNToNNEtaExcluIso(ener <<
777 } <<
778 } 807 }
779 << 808
780 << 809 if (sNNOmega < 1.e-9) sNNOmega = 0.;
781 G4double CrossSectionsMultiPionsAndResonan << 810
782 << 811 if (iso != 0) {
783 const G4double Ecm=0.001*ener; << 812 return sNNOmega;
784 G4double sNNOmega; // pp->pp+eta(+X) <<
785 G4double sNNOmega1; // np->np+eta(+X) <<
786 G4double x=Ecm*Ecm/7.06; <<
787 <<
788 if(Ecm>4.0) { <<
789 sNNOmega = 2.5*std::pow(x-1, 1.47) <<
790 } <<
791 else if(Ecm>2.802) { // 2802 MeV -> th <<
792 sNNOmega = (568.5254*Ecm*Ecm - 269 <<
793 if (sNNOmega <= NNToNNOmegaExcluIso(en <<
794 } <<
795 else { <<
796 sNNOmega = NNToNNOmegaExcluIso(ene <<
797 } <<
798 <<
799 if (sNNOmega < 1.e-9) sNNOmega = 0.; <<
800 <<
801 if (iso != 0) { <<
802 return sNNOmega; <<
803 } <<
804 <<
805 sNNOmega1 = 3.*sNNOmega; // 3.0: ratio <<
806 <<
807 sNNOmega = 2.*sNNOmega1-sNNOmega; <<
808 <<
809 if (sNNOmega < 1.e-9) sNNOmega = 0.; <<
810 <<
811 return sNNOmega; <<
812 } 813 }
813 << 814
814 << 815 sNNOmega1 = 3.*sNNOmega; // 3.0: ratio pn/pp (5 from theory; 2 from experiments)
815 G4double CrossSectionsMultiPionsAndResonan << 816
816 << 817 sNNOmega = 2.*sNNOmega1-sNNOmega;
817 const G4double ener=KinematicsUtils::t << 818
818 const G4int iso=ParticleTable::getIsos << 819 if (sNNOmega < 1.e-9) sNNOmega = 0.;
>> 820
>> 821 return sNNOmega;
>> 822 }
>> 823
>> 824
>> 825 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmega(Particle const * const particle1, Particle const * const particle2) {
>> 826
>> 827 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
>> 828 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
819 //jcd to be removed 829 //jcd to be removed
820 // return 0.; << 830 // return 0.;
821 //jcd 831 //jcd
822 if (iso != 0) { << 832 if (iso != 0) {
823 return NNToNNOmegaIso(ener, iso); << 833 return NNToNNOmegaIso(ener, iso);
824 } <<
825 else { <<
826 return 0.5*(NNToNNOmegaIso(ener, 0 <<
827 } <<
828 } 834 }
829 << 835 else {
830 G4double CrossSectionsMultiPionsAndResonan << 836 return 0.5*(NNToNNOmegaIso(ener, 0)+NNToNNOmegaIso(ener, 2));
831 <<
832 const G4double Ecm=0.001*ener; <<
833 G4double sNNOmega; // pp->pp+eta <<
834 G4double sNNOmega1; // np->np+eta <<
835 G4double sthroot=std::sqrt(7.06); <<
836 <<
837 if(Ecm>=3.0744) { // By hand (JCD) <<
838 sNNOmega = 330.*(Ecm-sthroot)/(1.0 <<
839 } <<
840 else if(Ecm>=2.65854){ <<
841 sNNOmega = -1208.09757*std::pow(Ec <<
842 } <<
843 else { <<
844 sNNOmega = 0. ; <<
845 } <<
846 <<
847 G4double Mn=ParticleTable::getRealMass <<
848 G4double Mp=ParticleTable::getRealMass <<
849 G4double Momega=ParticleTable::getReal <<
850 G4double Thr0=0.; <<
851 if (iso > 0) { <<
852 Thr0=2.*Mp+Momega; <<
853 } <<
854 else if (iso < 0) { <<
855 Thr0=2.*Mn+Momega; <<
856 } <<
857 else { <<
858 Thr0=Mn+Mp+Momega; <<
859 } <<
860 <<
861 if (sNNOmega < 1.e-9 || Ecm < Thr0) sN <<
862 <<
863 if (iso != 0) { <<
864 return sNNOmega/1000.; // paramete <<
865 } <<
866 <<
867 sNNOmega1 = 3*sNNOmega; // 3.0: ratio <<
868 <<
869 sNNOmega = 2*sNNOmega1-sNNOmega; <<
870 if (sNNOmega < 1.e-9 || Ecm < Thr0) sN <<
871 <<
872 return sNNOmega/1000.; // parameteriza <<
873 } 837 }
874 << 838 }
875 G4double CrossSectionsMultiPionsAndResonan << 839
>> 840 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExcluIso(const G4double ener, const G4int iso) {
>> 841
>> 842 const G4double Ecm=0.001*ener;
>> 843 G4double sNNOmega; // pp->pp+eta
>> 844 G4double sNNOmega1; // np->np+eta
>> 845 G4double sthroot=std::sqrt(7.06);
>> 846
>> 847 if(Ecm>=3.0744) { // By hand (JCD)
>> 848 sNNOmega = 330.*(Ecm-sthroot)/(1.05+std::pow((Ecm-sthroot),2));
>> 849 }
>> 850 else if(Ecm>=2.65854){
>> 851 sNNOmega = -1208.09757*std::pow(Ecm,3) + 10773.3322*std::pow(Ecm,2) - 31661.0223*Ecm + 30728.7241 ;
>> 852 }
>> 853 else {
>> 854 sNNOmega = 0. ;
>> 855 }
>> 856
>> 857 G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
>> 858 G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
>> 859 G4double Momega=ParticleTable::getRealMass(Omega)/1000.;
>> 860 G4double Thr0=0.;
>> 861 if (iso > 0) {
>> 862 Thr0=2.*Mp+Momega;
>> 863 }
>> 864 else if (iso < 0) {
>> 865 Thr0=2.*Mn+Momega;
>> 866 }
>> 867 else {
>> 868 Thr0=Mn+Mp+Momega;
>> 869 }
>> 870
>> 871 if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.; // Thr0: Ecm threshold
>> 872
>> 873 if (iso != 0) {
>> 874 return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
>> 875 }
>> 876
>> 877 sNNOmega1 = 3*sNNOmega; // 3.0: ratio pn/pp
>> 878
>> 879 sNNOmega = 2*sNNOmega1-sNNOmega;
>> 880 if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;
>> 881
>> 882 return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
>> 883 }
>> 884
>> 885 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExclu(Particle const * const particle1, Particle const * const particle2) {
876 //jcd to be removed 886 //jcd to be removed
877 // return 0.; << 887 // return 0.;
878 //jcd 888 //jcd
879 << 889
880 const G4double ener=KinematicsUtils::t << 890 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
881 const G4int iso=ParticleTable::getIsos << 891 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
882 << 892
883 if (iso != 0) { << 893 if (iso != 0) {
884 return NNToNNOmegaExcluIso(ener, i << 894 return NNToNNOmegaExcluIso(ener, iso);
885 } << 895 }
886 else { << 896 else {
887 return 0.5*(NNToNNOmegaExcluIso(en << 897 return 0.5*(NNToNNOmegaExcluIso(ener, 0)+NNToNNOmegaExcluIso(ener, 2));
888 } << 898 }
889 } << 899 }
890 << 900
891 << 901
892 G4double CrossSectionsMultiPionsAndResonan << 902 G4double CrossSectionsMultiPionsAndResonances::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
893 // << 903 //
894 // Nucleon-Nucleon producing xpi p << 904 // Nucleon-Nucleon producing xpi pions cross sections
895 // << 905 //
896 // assert(xpi>0 && xpi<=nMaxPiNN); 906 // assert(xpi>0 && xpi<=nMaxPiNN);
897 // assert(particle1->isNucleon() && particle2- 907 // assert(particle1->isNucleon() && particle2->isNucleon());
898 << 908
899 G4double oldXS1Pi=CrossSectionsMultiPi << 909 G4double oldXS1Pi=CrossSectionsMultiPions::NNToxPiNN(1,particle1, particle2);
900 G4double oldXS2Pi=CrossSectionsMultiPi << 910 G4double oldXS2Pi=CrossSectionsMultiPions::NNToxPiNN(2,particle1, particle2);
901 G4double oldXS3Pi=CrossSectionsMultiPi << 911 G4double oldXS3Pi=CrossSectionsMultiPions::NNToxPiNN(3,particle1, particle2);
902 G4double oldXS4Pi=CrossSectionsMultiPi << 912 G4double oldXS4Pi=CrossSectionsMultiPions::NNToxPiNN(4,particle1, particle2);
903 G4double xsEtaOmega=NNToNNEta(particle << 913 G4double xsEtaOmega=NNToNNEta(particle1, particle2)+NNToNNOmega(particle1, particle2);
904 G4double newXS1Pi=0.; << 914 G4double newXS1Pi=0.;
905 G4double newXS2Pi=0.; << 915 G4double newXS2Pi=0.;
906 G4double newXS3Pi=0.; << 916 G4double newXS3Pi=0.;
907 G4double newXS4Pi=0.; << 917 G4double newXS4Pi=0.;
908 << 918
909 if (xpi == 1) { << 919 if (xpi == 1) {
910 if (oldXS4Pi != 0. || oldXS3Pi != << 920 if (oldXS4Pi != 0. || oldXS3Pi != 0.)
911 newXS1Pi=oldXS1Pi; << 921 newXS1Pi=oldXS1Pi;
912 else if (oldXS2Pi != 0.) { << 922 else if (oldXS2Pi != 0.) {
913 newXS2Pi=oldXS2Pi-xsEtaOmega; << 923 newXS2Pi=oldXS2Pi-xsEtaOmega;
914 if (newXS2Pi < 0.) << 924 if (newXS2Pi < 0.)
915 newXS1Pi=oldXS1Pi-(xsEtaOm << 925 newXS1Pi=oldXS1Pi-(xsEtaOmega-oldXS2Pi);
916 else << 926 else
917 newXS1Pi=oldXS1Pi; << 927 newXS1Pi=oldXS1Pi;
918 } << 928 }
919 else << 929 else
920 newXS1Pi=oldXS1Pi-xsEtaOmega; << 930 newXS1Pi=oldXS1Pi-xsEtaOmega;
921 return newXS1Pi; << 931 return newXS1Pi;
922 } << 932 }
923 else if (xpi == 2) { << 933 else if (xpi == 2) {
924 if (oldXS4Pi != 0.) << 934 if (oldXS4Pi != 0.)
925 newXS2Pi=oldXS2Pi; << 935 newXS2Pi=oldXS2Pi;
926 else if (oldXS3Pi != 0.) { << 936 else if (oldXS3Pi != 0.) {
927 newXS3Pi=oldXS3Pi-xsEtaOmega; << 937 newXS3Pi=oldXS3Pi-xsEtaOmega;
928 if (newXS3Pi < 0.) << 938 if (newXS3Pi < 0.)
929 newXS2Pi=oldXS2Pi-(xsEtaOm << 939 newXS2Pi=oldXS2Pi-(xsEtaOmega-oldXS3Pi);
930 else << 940 else
931 newXS2Pi=oldXS2Pi; << 941 newXS2Pi=oldXS2Pi;
932 } << 942 }
933 else { << 943 else {
934 newXS2Pi=oldXS2Pi-xsEtaOmega; << 944 newXS2Pi=oldXS2Pi-xsEtaOmega;
935 if (newXS2Pi < 0.) << 945 if (newXS2Pi < 0.)
936 newXS2Pi=0.; << 946 newXS2Pi=0.;
937 } << 947 }
938 return newXS2Pi; << 948 return newXS2Pi;
939 } << 949 }
940 else if (xpi == 3) { << 950 else if (xpi == 3) {
941 if (oldXS4Pi != 0.) { << 951 if (oldXS4Pi != 0.) {
942 newXS4Pi=oldXS4Pi-xsEtaOmega; << 952 newXS4Pi=oldXS4Pi-xsEtaOmega;
943 if (newXS4Pi < 0.) << 953 if (newXS4Pi < 0.)
944 newXS3Pi=oldXS3Pi-(xsEtaOm << 954 newXS3Pi=oldXS3Pi-(xsEtaOmega-oldXS4Pi);
945 else << 955 else
946 newXS3Pi=oldXS3Pi; << 956 newXS3Pi=oldXS3Pi;
947 } << 957 }
948 else { << 958 else {
949 newXS3Pi=oldXS3Pi-xsEtaOmega; << 959 newXS3Pi=oldXS3Pi-xsEtaOmega;
950 if (newXS3Pi < 0.) << 960 if (newXS3Pi < 0.)
951 newXS3Pi=0.; << 961 newXS3Pi=0.;
952 } << 962 }
953 return newXS3Pi; << 963 return newXS3Pi;
954 } << 964 }
955 else if (xpi == 4) { << 965 else if (xpi == 4) {
956 newXS4Pi=oldXS4Pi-xsEtaOmega; << 966 newXS4Pi=oldXS4Pi-xsEtaOmega;
957 if (newXS4Pi < 0.) << 967 if (newXS4Pi < 0.)
958 newXS4Pi=0.; << 968 newXS4Pi=0.;
959 return newXS4Pi; << 969 return newXS4Pi;
960 } << 970 }
961 << 971
962 else // should never reach this point << 972 else // should never reach this point
963 return 0.; << 973 return 0.;
964 } << 974 }
965 << 975
966 << 976
967 G4double CrossSectionsMultiPionsAndResonan << 977 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePi(Particle const * const particle1, Particle const * const particle2) {
968 // Cross section for nucleon-nucleon p << 978 // Cross section for nucleon-nucleon producing one eta and one pion
969 << 979
970 const G4int iso=ParticleTable::getIsos << 980 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
971 if (iso!=0) << 981 if (iso!=0)
972 return 0.; << 982 return 0.;
973 << 983
974 const G4double ener=KinematicsUtils::t << 984 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 686.987; // 686.987 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.; << 985 if (ener < 2018.563) return 0.;
976 << 986
977 const G4double xsiso2=CrossSectionsMul << 987 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
978 const G4double xsiso0=CrossSectionsMul << 988 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
979 << 989
980 return 0.25*(CrossSectionsMultiPions:: << 990 return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
981 } << 991 }
982 << 992
983 G4double CrossSectionsMultiPionsAndResonan << 993 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
984 const G4double ener=KinematicsUtils::t << 994 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 686.987; // 686.987 MeV translation to open pion production in NNEta
985 if (ener < 2018.563) return 0.; << 995 if (ener < 2018.563) return 0.;
986 const G4int iso=ParticleTable::getIsos << 996 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
987 << 997
988 const G4double xsiso2=CrossSectionsMul << 998 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
989 if (iso != 0) << 999 if (iso != 0)
990 return CrossSectionsMultiPions::NN << 1000 return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
991 else { << 1001 else {
992 const G4double xsiso0=CrossSection << 1002 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
993 return 0.5*(CrossSectionsMultiPion << 1003 return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
994 } << 1004 }
995 } << 1005 }
996 << 1006
997 G4double CrossSectionsMultiPionsAndResonan << 1007 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaTwoPi(Particle const * const particle1, Particle const * const particle2) {
998 // << 1008 //
999 // Nucleon-Nucleon producing one e << 1009 // Nucleon-Nucleon producing one eta and two pions
1000 // << 1010 //
1001 const G4double ener=KinematicsUtils:: << 1011 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 686.987; // 686.987 MeV translation to open pion production in NNEta
1002 if (ener < 2018.563) return 0.; << 1012 if (ener < 2018.563) return 0.;
1003 const G4int iso=ParticleTable::getIso << 1013 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1004 << 1014
1005 << 1015
1006 const G4double xsiso2=CrossSectionsMu << 1016 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1007 if (iso != 0) { << 1017 if (iso != 0) {
1008 return CrossSectionsMultiPions::N << 1018 return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1009 } << 1019 }
1010 else { << 1020 else {
1011 const G4double xsiso0=CrossSectio << 1021 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1012 return 0.5*(CrossSectionsMultiPio << 1022 return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
1013 } << 1023 }
1014 } << 1024 }
1015 << 1025
1016 G4double CrossSectionsMultiPionsAndResona << 1026 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaThreePi(Particle const * const particle1, Particle const * const particle2) {
1017 // << 1027 //
1018 // Nucleon-Nucleon producing one << 1028 // Nucleon-Nucleon producing one eta and three pions
1019 // << 1029 //
1020 << 1030
1021 const G4double ener=KinematicsUtils:: << 1031 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 686.987; // 686.987 MeV translation to open pion production in NNEta
1022 if (ener < 2018.563) return 0.; << 1032 if (ener < 2018.563) return 0.;
1023 const G4int iso=ParticleTable::getIso << 1033 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1024 << 1034
1025 << 1035
1026 const G4double xsiso2=CrossSectionsMu << 1036 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1027 const G4double xs1pi2=CrossSectionsMu << 1037 const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
1028 const G4double xs2pi2=CrossSectionsMu << 1038 const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1029 if (iso != 0) << 1039 if (iso != 0)
1030 return CrossSectionsMultiPions::N << 1040 return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
1031 else { << 1041 else {
1032 const G4double xsiso0=CrossSectio << 1042 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1033 const G4double xs1pi0=CrossSectio << 1043 const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
1034 const G4double xs2pi0=CrossSectio << 1044 const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
1035 return 0.5*(CrossSectionsMultiPio << 1045 return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
1036 } << 1046 }
1037 } << 1047 }
1038 << 1048
1039 G4double CrossSectionsMultiPionsAndResona << 1049 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaFourPi(Particle const * const particle1, Particle const * const particle2) {
1040 // << 1050 //
1041 // Nucleon-Nucleon producing one << 1051 // Nucleon-Nucleon producing one eta and four pions
1042 // << 1052 //
1043 << 1053
1044 const G4double ener=KinematicsUtils:: << 1054 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 686.987; // 686.987 MeV translation to open pion production in NNEta
1045 if (ener < 2018.563) return 0.; << 1055 if (ener < 2018.563) return 0.;
1046 const G4double s = ener*ener; << 1056 const G4double s = ener*ener;
1047 const G4int i = ParticleTable::getIso << 1057 const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1048 G4double xsinelas; << 1058 G4double xsinelas;
1049 if (i!=0) << 1059 if (i!=0)
1050 xsinelas=CrossSectionsMultiPions: << 1060 xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1051 else << 1061 else
1052 xsinelas=0.5*(CrossSectionsMultiP << 1062 xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1053 if (xsinelas <= 1.e-9) return 0.; << 1063 if (xsinelas <= 1.e-9) return 0.;
1054 G4double ratio=(NNToNNEta(particle1, << 1064 G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;
1055 if(s<6.25E6) << 1065 if(s<6.25E6)
1056 return 0.; << 1066 return 0.;
1057 const G4double sigma = NNToNNEta(part << 1067 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.); << 1068 return ((sigma>1.e-9) ? sigma : 0.);
1059 } << 1069 }
1060 << 1070
1061 G4double CrossSectionsMultiPionsAndResona << 1071 G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
1062 // << 1072 //
1063 // Nucleon-Nucleon producing one << 1073 // Nucleon-Nucleon producing one eta and xpi pions
1064 // << 1074 //
1065 // assert(xpi>0 && xpi<=nMaxPiNN); 1075 // assert(xpi>0 && xpi<=nMaxPiNN);
1066 // assert(particle1->isNucleon() && particle2 1076 // assert(particle1->isNucleon() && particle2->isNucleon());
>> 1077
>> 1078 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 686.987; // 686.987 MeV translation to open pion production in NNEta
>> 1079 if (ener < 2018.563) return 0.;
>> 1080 const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
>> 1081 G4double xsinelas;
>> 1082 if (i!=0)
>> 1083 xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
>> 1084 else
>> 1085 xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
>> 1086 if (xsinelas <= 1.e-9) return 0.;
>> 1087 G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;
>> 1088
>> 1089 if (xpi == 1)
>> 1090 return NNToNNEtaOnePi(particle1, particle2)*ratio;
>> 1091 else if (xpi == 2)
>> 1092 return NNToNNEtaTwoPi(particle1, particle2)*ratio;
>> 1093 else if (xpi == 3)
>> 1094 return NNToNNEtaThreePi(particle1, particle2)*ratio;
>> 1095 else if (xpi == 4)
>> 1096 return NNToNNEtaFourPi(particle1, particle2);
>> 1097 else // should never reach this point
>> 1098 return 0.;
>> 1099 }
1067 1100
1068 const G4double ener=KinematicsUtils:: << 1101
1069 if (ener < 2018.563) return 0.; << 1102 G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaEta(Particle const * const p1, Particle const * const p2) {
1070 const G4int i = ParticleTable::getIso <<
1071 G4double xsinelas; <<
1072 if (i!=0) <<
1073 xsinelas=CrossSectionsMultiPions: <<
1074 else <<
1075 xsinelas=0.5*(CrossSectionsMultiP <<
1076 if (xsinelas <= 1.e-9) return 0.; <<
1077 G4double ratio=(NNToNNEta(particle1, <<
1078 <<
1079 if (xpi == 1) <<
1080 return NNToNNEtaOnePi(particle1, <<
1081 else if (xpi == 2) <<
1082 return NNToNNEtaTwoPi(particle1, <<
1083 else if (xpi == 3) <<
1084 return NNToNNEtaThreePi(particle1 <<
1085 else if (xpi == 4) <<
1086 return NNToNNEtaFourPi(particle1, <<
1087 else // should never reach this point <<
1088 return 0.; <<
1089 } <<
1090 <<
1091 <<
1092 G4double CrossSectionsMultiPionsAndResona <<
1093 // assert(p1->isNucleon() && p2->isNucleon()) 1103 // assert(p1->isNucleon() && p2->isNucleon());
1094 const G4int i = ParticleTable::getIso << 1104 const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1095 const G4double ener=KinematicsUtils:: << 1105 const G4double ener=KinematicsUtils::totalEnergyInCM(p1, p2) - 686.987; // 686.987 MeV translation to open pion production in NNEta
1096 if (ener < 2018.563) return 0.; << 1106 if (ener < 2018.563) return 0.;
1097 G4double xsinelas; << 1107 G4double xsinelas;
1098 if (i!=0) << 1108 if (i!=0)
1099 xsinelas=CrossSectionsMultiPions: << 1109 xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1100 else << 1110 else
1101 xsinelas=0.5*(CrossSectionsMultiP << 1111 xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1102 if (xsinelas <= 1.e-9) return 0.; << 1112 if (xsinelas <= 1.e-9) return 0.;
1103 G4double ratio=(NNToNNEta(p1, p2)-NNT << 1113 G4double ratio=(NNToNNEta(p1, p2)-NNToNNEtaExclu(p1, p2))/xsinelas;
1104 G4double sigma = NNToNNEtaOnePiOrDelt << 1114 G4double sigma = NNToNNEtaOnePiOrDelta(p1, p2)*ratio;
1105 if(i==0) << 1115 if(i==0)
1106 sigma *= 0.5; << 1116 sigma *= 0.5;
1107 return sigma; << 1117 return sigma;
1108 } << 1118 }
1109 1119
1110 << 1120
1111 G4double CrossSectionsMultiPionsAndResona << 1121 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePi(Particle const * const particle1, Particle const * const particle2) {
1112 // Cross section for nucleon-nucleon << 1122 // Cross section for nucleon-nucleon producing one omega and one pion
1113 << 1123
1114 const G4int iso=ParticleTable::getIso << 1124 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1115 if (iso!=0) << 1125 if (iso!=0)
1116 return 0.; << 1126 return 0.;
1117 << 1127
1118 const G4double ener=KinematicsUtils:: << 1128 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.; << 1129 if (ener < 2018.563) return 0.;
1120 << 1130
1121 const G4double xsiso2=CrossSectionsMu << 1131 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1122 const G4double xsiso0=CrossSectionsMu << 1132 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1123 << 1133
1124 return 0.25*(CrossSectionsMultiPions: << 1134 return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
1125 } << 1135 }
1126 << 1136
1127 G4double CrossSectionsMultiPionsAndResona << 1137 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
1128 const G4double ener=KinematicsUtils:: << 1138 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.; << 1139 if (ener < 2018.563) return 0.;
1130 const G4int iso=ParticleTable::getIso << 1140 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1131 << 1141
1132 const G4double xsiso2=CrossSectionsMu << 1142 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1133 if (iso != 0) << 1143 if (iso != 0)
1134 return CrossSectionsMultiPions::N << 1144 return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
1135 else { << 1145 else {
1136 const G4double xsiso0=CrossSectio << 1146 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1137 return 0.5*(CrossSectionsMultiPio << 1147 return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
1138 } << 1148 }
1139 } << 1149 }
1140 << 1150
1141 G4double CrossSectionsMultiPionsAndResona << 1151 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaTwoPi(Particle const * const particle1, Particle const * const particle2) {
1142 // << 1152 //
1143 // Nucleon-Nucleon producing one << 1153 // Nucleon-Nucleon producing one omega and two pions
1144 // << 1154 //
1145 const G4double ener=KinematicsUtils:: << 1155 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.; << 1156 if (ener < 2018.563) return 0.;
1147 const G4int iso=ParticleTable::getIso << 1157 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1148 << 1158
1149 const G4double xsiso2=CrossSectionsMu << 1159
1150 if (iso != 0) { << 1160 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1151 return CrossSectionsMultiPions::N << 1161 if (iso != 0) {
1152 } << 1162 return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1153 else { << 1163 }
1154 const G4double xsiso0=CrossSectio << 1164 else {
1155 return 0.5*(CrossSectionsMultiPio << 1165 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1156 } << 1166 return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
1157 } << 1167 }
1158 << 1168 }
1159 G4double CrossSectionsMultiPionsAndResona << 1169
1160 // << 1170 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaThreePi(Particle const * const particle1, Particle const * const particle2) {
1161 // Nucleon-Nucleon producing one << 1171 //
1162 // << 1172 // Nucleon-Nucleon producing one omega and three pions
1163 << 1173 //
1164 const G4double ener=KinematicsUtils:: << 1174
1165 if (ener < 2018.563) return 0.; << 1175 const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1166 const G4int iso=ParticleTable::getIso << 1176 if (ener < 2018.563) return 0.;
1167 << 1177 const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1168 << 1178
1169 const G4double xsiso2=CrossSectionsMu << 1179
1170 const G4double xs1pi2=CrossSectionsMu << 1180 const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1171 const G4double xs2pi2=CrossSectionsMu << 1181 const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
1172 if (iso != 0) << 1182 const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1173 return CrossSectionsMultiPions::N << 1183 if (iso != 0)
1174 else { << 1184 return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
1175 const G4double xsiso0=CrossSectio << 1185 else {
1176 const G4double xs1pi0=CrossSectio << 1186 const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1177 const G4double xs2pi0=CrossSectio << 1187 const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
1178 return 0.5*(CrossSectionsMultiPio << 1188 const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
1179 } << 1189 return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
1180 } << 1190 }
1181 << 1191 }
1182 G4double CrossSectionsMultiPionsAndResona << 1192
1183 // << 1193 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaFourPi(Particle const * const particle1, Particle const * const particle2) {
1184 // Nucleon-Nucleon producing one << 1194 //
1185 // << 1195 // Nucleon-Nucleon producing one omega and four pions
>> 1196 //
1186 //jcd to be removed 1197 //jcd to be removed
1187 // return 0.; << 1198 // return 0.;
1188 //jcd 1199 //jcd
1189 << 1200
1190 const G4double ener=KinematicsUtils:: << 1201 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.; << 1202 if (ener < 2018.563) return 0.;
1192 const G4double s = ener*ener; << 1203 const G4double s = ener*ener;
1193 const G4int i = ParticleTable::getIso << 1204 const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1194 G4double xsinelas; << 1205 G4double xsinelas;
1195 if (i!=0) << 1206 if (i!=0)
1196 xsinelas=CrossSectionsMultiPions: << 1207 xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1197 else << 1208 else
1198 xsinelas=0.5*(CrossSectionsMultiP << 1209 xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1199 if (xsinelas <= 1.e-9) return 0.; << 1210 if (xsinelas <= 1.e-9) return 0.;
1200 G4double ratio=(NNToNNOmega(particle1 << 1211 G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
1201 if(s<6.25E6) << 1212 if(s<6.25E6)
1202 return 0.; << 1213 return 0.;
1203 const G4double sigma = NNToNNOmega(pa << 1214 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.); << 1215 return ((sigma>1.e-9) ? sigma : 0.);
1205 } << 1216 }
1206 << 1217
1207 G4double CrossSectionsMultiPionsAndResona << 1218 G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
1208 // << 1219 //
1209 // Nucleon-Nucleon producing one << 1220 // Nucleon-Nucleon producing one omega and xpi pions
1210 // << 1221 //
1211 // assert(xpi>0 && xpi<=nMaxPiNN); 1222 // assert(xpi>0 && xpi<=nMaxPiNN);
1212 // assert(particle1->isNucleon() && particle2 1223 // assert(particle1->isNucleon() && particle2->isNucleon());
1213 //jcd to be removed 1224 //jcd to be removed
1214 // return 0.; << 1225 // return 0.;
1215 //jcd 1226 //jcd
1216 << 1227
1217 const G4double ener=KinematicsUtils:: << 1228 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.; << 1229 if (ener < 2018.563) return 0.;
1219 const G4int i = ParticleTable::getIso << 1230 const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1220 G4double xsinelas; << 1231 G4double xsinelas;
1221 if (i!=0) << 1232 if (i!=0)
1222 xsinelas=CrossSectionsMultiPions: << 1233 xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1223 else << 1234 else
1224 xsinelas=0.5*(CrossSectionsMultiP << 1235 xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1225 if (xsinelas <= 1.e-9) return 0.; << 1236 if (xsinelas <= 1.e-9) return 0.;
1226 G4double ratio=(NNToNNOmega(particle1 << 1237 G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
1227 << 1238
1228 if (xpi == 1) << 1239 if (xpi == 1)
1229 return NNToNNOmegaOnePi(particle1 << 1240 return NNToNNOmegaOnePi(particle1, particle2)*ratio;
1230 else if (xpi == 2) << 1241 else if (xpi == 2)
1231 return NNToNNOmegaTwoPi(particle1 << 1242 return NNToNNOmegaTwoPi(particle1, particle2)*ratio;
1232 else if (xpi == 3) << 1243 else if (xpi == 3)
1233 return NNToNNOmegaThreePi(particl << 1244 return NNToNNOmegaThreePi(particle1, particle2)*ratio;
1234 else if (xpi == 4) << 1245 else if (xpi == 4)
1235 return NNToNNOmegaFourPi(particle << 1246 return NNToNNOmegaFourPi(particle1, particle2);
1236 else // should never reach this point << 1247 else // should never reach this point
1237 return 0.; << 1248 return 0.;
1238 } << 1249 }
1239 << 1250
1240 << 1251
1241 G4double CrossSectionsMultiPionsAndResona << 1252 G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaOmega(Particle const * const p1, Particle const * const p2) {
1242 // assert(p1->isNucleon() && p2->isNucleon()) 1253 // assert(p1->isNucleon() && p2->isNucleon());
1243 //jcd to be removed 1254 //jcd to be removed
1244 // return 0.; << 1255 // return 0.;
1245 //jcd 1256 //jcd
1246 const G4int i = ParticleTable::getIso << 1257 const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1247 const G4double ener=KinematicsUtils:: << 1258 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.; << 1259 if (ener < 2018.563) return 0.;
1249 G4double xsinelas; << 1260 G4double xsinelas;
1250 if (i!=0) << 1261 if (i!=0)
1251 xsinelas=CrossSectionsMultiPions: << 1262 xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1252 else << 1263 else
1253 xsinelas=0.5*(CrossSectionsMultiP << 1264 xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1254 if (xsinelas <= 1.e-9) return 0.; << 1265 if (xsinelas <= 1.e-9) return 0.;
1255 G4double ratio=(NNToNNOmega(p1, p2)-N << 1266 G4double ratio=(NNToNNOmega(p1, p2)-NNToNNOmegaExclu(p1, p2))/xsinelas;
1256 G4double sigma = NNToNNOmegaOnePiOrDe << 1267 G4double sigma = NNToNNOmegaOnePiOrDelta(p1, p2)*ratio;
1257 if(i==0) << 1268 if(i==0)
1258 sigma *= 0.5; << 1269 sigma *= 0.5;
1259 return sigma; << 1270 return sigma;
1260 } << 1271 }
1261 << 1272
1262 << 1273
1263 } // namespace G4INCL 1274 } // namespace G4INCL
1264 1275
1265 1276