| 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 G4INCLCrossSectionsStrangeness.cc 38 /** \file G4INCLCrossSectionsStrangeness.cc
39 * \brief Multipion, mesonic Resonances and st 39 * \brief Multipion, mesonic Resonances and strange cross sections
40 * 40 *
41 * \date 1st March 2016 41 * \date 1st March 2016
42 * \author Jason Hirtz 42 * \author Jason Hirtz
43 */ 43 */
44 44
45 #include "G4INCLCrossSectionsStrangeness.hh" 45 #include "G4INCLCrossSectionsStrangeness.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 CrossSectionsStrangeness::nMax 63 const G4int CrossSectionsStrangeness::nMaxPiNN = 4;
64 const G4int CrossSectionsStrangeness::nMax 64 const G4int CrossSectionsStrangeness::nMaxPiPiN = 4;
65 65
66 CrossSectionsStrangeness::CrossSectionsStr 66 CrossSectionsStrangeness::CrossSectionsStrangeness() :
67 s11pzHC(-2.228000000000294018,8.7560000000 67 s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
68 s01ppHC(2.0570000000126518344,-6.029000000 68 s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
69 s01pzHC(0.18030000000000441851,7.870099999 69 s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
70 s11pmHC(0.20590000000000031866,3.345099999 70 s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
71 s12pmHC(-0.77235999999999901328,4.26265999 71 s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
72 s12ppHC(-0.75724999999999975664,2.09343999 72 s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
73 s12zzHC(-0.89599999999996965072,7.88299999 73 s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
74 s02pzHC(-1.0579999999999967036,11.11399999 74 s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
75 s02pmHC(2.4009000000012553286,-7.768000000 75 s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
76 s12mzHC(-0.21858699999999976269,1.91489999 76 s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
77 { 77 {
78 } 78 }
79 79
80 /// \brief redefining previous cross section 80 /// \brief redefining previous cross sections
81 81
82 G4double CrossSectionsStrangeness::total(P 82 G4double CrossSectionsStrangeness::total(Particle const * const p1, Particle const * const p2) {
83 G4double inelastic; 83 G4double inelastic;
84 if(p1->isNucleon() && p2->isNucleon()) 84 if(p1->isNucleon() && p2->isNucleon()) {
85 return CrossSectionsMultiPions::NN 85 return CrossSectionsMultiPions::NNTot(p1, p2);
86 } else if((p1->isNucleon() && p2->isDe 86 } else if((p1->isNucleon() && p2->isDelta()) ||
87 (p1->isDelta() && p2->isNucl 87 (p1->isDelta() && p2->isNucleon())) {
88 inelastic = CrossSectionsMultiPion 88 inelastic = CrossSectionsMultiPions::NDeltaToNN(p1, p2) + NDeltaToNLK(p1, p2) + NDeltaToNSK(p1, p2) + NDeltaToDeltaLK(p1, p2) + NDeltaToDeltaSK(p1, p2) + NDeltaToNNKKb(p1, p2);
89 } else if((p1->isNucleon() && p2->isPi 89 } else if((p1->isNucleon() && p2->isPion()) ||
90 (p1->isPion() && p2->isNucle 90 (p1->isPion() && p2->isNucleon())) {
91 return CrossSectionsMultiPions::pi 91 return CrossSectionsMultiPions::piNTot(p1,p2);
92 } else if((p1->isNucleon() && p2->isEt 92 } else if((p1->isNucleon() && p2->isEta()) ||
93 (p1->isEta() && p2->isNucleo 93 (p1->isEta() && p2->isNucleon())) {
94 inelastic = CrossSectionsMultiPion 94 inelastic = CrossSectionsMultiPionsAndResonances::etaNToPiN(p1,p2) + CrossSectionsMultiPionsAndResonances::etaNToPiPiN(p1,p2);
95 } else if((p1->isNucleon() && p2->isOm 95 } else if((p1->isNucleon() && p2->isOmega()) ||
96 (p1->isOmega() && p2->isNucl 96 (p1->isOmega() && p2->isNucleon())) {
97 inelastic = CrossSectionsMultiPion 97 inelastic = CrossSectionsMultiPionsAndResonances::omegaNInelastic(p1,p2);
98 } else if((p1->isNucleon() && p2->isEt 98 } else if((p1->isNucleon() && p2->isEtaPrime()) ||
99 (p1->isEtaPrime() && p2->isN 99 (p1->isEtaPrime() && p2->isNucleon())) {
100 inelastic = CrossSectionsMultiPion 100 inelastic = CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(p1,p2);
101 } else if((p1->isNucleon() && p2->isLa 101 } else if((p1->isNucleon() && p2->isLambda()) ||
102 (p1->isLambda() && p2->isNuc 102 (p1->isLambda() && p2->isNucleon())) {
103 inelastic = NLToNS(p1,p2); 103 inelastic = NLToNS(p1,p2);
104 } else if((p1->isNucleon() && p2->isSi 104 } else if((p1->isNucleon() && p2->isSigma()) ||
105 (p1->isSigma() && p2->isNucl 105 (p1->isSigma() && p2->isNucleon())) {
106 inelastic = NSToNL(p1,p2) + NSToNS 106 inelastic = NSToNL(p1,p2) + NSToNS(p1,p2);
107 } else if((p1->isNucleon() && p2->isKa 107 } else if((p1->isNucleon() && p2->isKaon()) ||
108 (p1->isKaon() && p2->isNucle 108 (p1->isKaon() && p2->isNucleon())) {
109 inelastic = NKToNK(p1,p2) + NKToNK 109 inelastic = NKToNK(p1,p2) + NKToNKpi(p1,p2) + NKToNK2pi(p1,p2);
110 } else if((p1->isNucleon() && p2->isAn 110 } else if((p1->isNucleon() && p2->isAntiKaon()) ||
111 (p1->isAntiKaon() && p2->isN 111 (p1->isAntiKaon() && p2->isNucleon())) {
112 inelastic = NKbToLpi(p1,p2) + NKbT 112 inelastic = NKbToLpi(p1,p2) + NKbToSpi(p1,p2) + NKbToL2pi(p1,p2) + NKbToS2pi(p1,p2) + NKbToNKb(p1,p2) + NKbToNKbpi(p1,p2) + NKbToNKb2pi(p1,p2);
113 } else { 113 } else {
114 inelastic = 0.; 114 inelastic = 0.;
115 } 115 }
116 return inelastic + elastic(p1, p2); 116 return inelastic + elastic(p1, p2);
117 } 117 }
118 118
119 G4double CrossSectionsStrangeness::elastic 119 G4double CrossSectionsStrangeness::elastic(Particle const * const p1, Particle const * const p2) {
120 if((p1->isNucleon()||p1->isDelta()) && 120 if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){ // N-N, N-Delta, Delta-Delta
121 return CrossSectionsMultiPions::el 121 return CrossSectionsMultiPions::elastic(p1, p2);
122 } 122 }
123 else if ((p1->isNucleon() && p2->isPio 123 else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
124 return CrossSectionsMultiPions::el 124 return CrossSectionsMultiPions::elastic(p1, p2);
125 } 125 }
126 else if ((p1->isNucleon() && p2->isEta 126 else if ((p1->isNucleon() && p2->isEta()) || (p2->isNucleon() && p1->isEta())){
127 return CrossSectionsMultiPionsAndR 127 return CrossSectionsMultiPionsAndResonances::etaNElastic(p1, p2);
128 } 128 }
129 else if ((p1->isNucleon() && p2->isHyp 129 else if ((p1->isNucleon() && p2->isHyperon()) || (p2->isNucleon() && p1->isHyperon())){
130 return NYelastic(p1, p2); 130 return NYelastic(p1, p2);
131 } 131 }
132 else if ((p1->isNucleon() && p2->isKao 132 else if ((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon())){
133 return NKelastic(p1, p2); 133 return NKelastic(p1, p2);
134 } 134 }
135 else if ((p1->isNucleon() && p2->isAnt 135 else if ((p1->isNucleon() && p2->isAntiKaon()) || (p2->isNucleon() && p1->isAntiKaon())){
136 return NKbelastic(p1, p2); 136 return NKbelastic(p1, p2);
137 } 137 }
138 else { 138 else {
139 return 0.0; 139 return 0.0;
140 } 140 }
141 } 141 }
142 142
143 G4double CrossSectionsStrangeness::piNToxP 143 G4double CrossSectionsStrangeness::piNToxPiN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
144 // 144 //
145 // pion-Nucleon producing xpi pion 145 // pion-Nucleon producing xpi pions cross sections (corrected due to new particles)
146 // 146 //
147 // assert(xpi>1 && xpi<=nMaxPiPiN); 147 // assert(xpi>1 && xpi<=nMaxPiPiN);
148 // assert((particle1->isNucleon() && particle2 148 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
149 149
150 const G4double oldXS2Pi=CrossSecti 150 const G4double oldXS2Pi=CrossSectionsMultiPions::piNToxPiN(2,particle1, particle2);
151 const G4double oldXS3Pi=CrossSectionsM 151 const G4double oldXS3Pi=CrossSectionsMultiPions::piNToxPiN(3,particle1, particle2);
152 const G4double oldXS4Pi=CrossSectionsM 152 const G4double oldXS4Pi=CrossSectionsMultiPions::piNToxPiN(4,particle1, particle2);
153 const G4double xsEta=CrossSectionsMult 153 const G4double xsEta=CrossSectionsMultiPionsAndResonances::piNToEtaN(particle1, particle2);
154 const G4double xsOmega=CrossSectionsMu 154 const G4double xsOmega=CrossSectionsMultiPionsAndResonances::piNToOmegaN(particle1, particle2);
155 const G4double xs1=NpiToLK(particle2, 155 const G4double xs1=NpiToLK(particle2, particle1);
156 const G4double xs2=NpiToSK(particle1, 156 const G4double xs2=NpiToSK(particle1, particle2);
157 const G4double xs3=NpiToLKpi(particle1 157 const G4double xs3=NpiToLKpi(particle1, particle2);
158 const G4double xs4=NpiToSKpi(particle1 158 const G4double xs4=NpiToSKpi(particle1, particle2);
159 const G4double xs5=NpiToLK2pi(particle 159 const G4double xs5=NpiToLK2pi(particle1, particle2);
160 const G4double xs6=NpiToSK2pi(particle 160 const G4double xs6=NpiToSK2pi(particle1, particle2);
161 const G4double xs7=NpiToNKKb(particle1 161 const G4double xs7=NpiToNKKb(particle1, particle2);
162 const G4double xs8=NpiToMissingStrange 162 const G4double xs8=NpiToMissingStrangeness(particle1, particle2);
163 const G4double xs0 = xs1 + xs2 + xs3 + 163 const G4double xs0 = xs1 + xs2 + xs3 + xs4 + xs5 + xs6 + xs7 +xs8;
164 G4double newXS2Pi=0.; 164 G4double newXS2Pi=0.;
165 G4double newXS3Pi=0.; 165 G4double newXS3Pi=0.;
166 G4double newXS4Pi=0.; 166 G4double newXS4Pi=0.;
167 167
168 if (xpi == 2) { 168 if (xpi == 2) {
169 if (oldXS4Pi != 0.) 169 if (oldXS4Pi != 0.)
170 newXS2Pi=oldXS2Pi; 170 newXS2Pi=oldXS2Pi;
171 else if (oldXS3Pi != 0.) { 171 else if (oldXS3Pi != 0.) {
172 newXS3Pi=oldXS3Pi-xsEta-xsOmeg 172 newXS3Pi=oldXS3Pi-xsEta-xsOmega-xs0;
173 if (newXS3Pi < 1.e-09) 173 if (newXS3Pi < 1.e-09)
174 newXS2Pi=oldXS2Pi-(xsEta+x 174 newXS2Pi=oldXS2Pi-(xsEta+xsOmega+xs0-oldXS3Pi);
175 else 175 else
176 newXS2Pi=oldXS2Pi; 176 newXS2Pi=oldXS2Pi;
177 } 177 }
178 else { 178 else {
179 newXS2Pi=oldXS2Pi-xsEta-xsOmeg 179 newXS2Pi=oldXS2Pi-xsEta-xsOmega-xs0;
180 if (newXS2Pi < 1.e-09 && newXS 180 if (newXS2Pi < 1.e-09 && newXS2Pi!=0.){
181 newXS2Pi=0.; 181 newXS2Pi=0.;
182 } 182 }
183 } 183 }
184 return newXS2Pi; 184 return newXS2Pi;
185 } 185 }
186 else if (xpi == 3) { 186 else if (xpi == 3) {
187 if (oldXS4Pi != 0.) { 187 if (oldXS4Pi != 0.) {
188 newXS4Pi=oldXS4Pi-xsEta-xsOmeg 188 newXS4Pi=oldXS4Pi-xsEta-xsOmega-xs0;
189 if (newXS4Pi < 1.e-09) 189 if (newXS4Pi < 1.e-09)
190 newXS3Pi=oldXS3Pi-(xsEta+x 190 newXS3Pi=oldXS3Pi-(xsEta+xsOmega+xs0-oldXS4Pi);
191 else 191 else
192 newXS3Pi=oldXS3Pi; 192 newXS3Pi=oldXS3Pi;
193 } 193 }
194 else { 194 else {
195 newXS3Pi=oldXS3Pi-xsEta-xsOmeg 195 newXS3Pi=oldXS3Pi-xsEta-xsOmega-xs0;
196 if (newXS3Pi < 1.e-09){ 196 if (newXS3Pi < 1.e-09){
197 newXS3Pi=0.; 197 newXS3Pi=0.;
198 } 198 }
199 } 199 }
200 return newXS3Pi; 200 return newXS3Pi;
201 } 201 }
202 else if (xpi == 4) { 202 else if (xpi == 4) {
203 newXS4Pi=oldXS4Pi-xsEta-xsOmega-xs 203 newXS4Pi=oldXS4Pi-xsEta-xsOmega-xs0;
204 if (newXS4Pi < 1.e-09){ 204 if (newXS4Pi < 1.e-09){
205 newXS4Pi=0.; 205 newXS4Pi=0.;
206 } 206 }
207 return newXS4Pi; 207 return newXS4Pi;
208 } 208 }
209 else // should never reach this point 209 else // should never reach this point
210 return 0.; 210 return 0.;
211 } 211 }
212 212
213 G4double CrossSectionsStrangeness::NNToxPi 213 G4double CrossSectionsStrangeness::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
214 // 214 //
215 // Nucleon-Nucleon producing xpi p 215 // Nucleon-Nucleon producing xpi pions cross sections corrected
216 // 216 //
217 // assert(xpi>0 && xpi<=nMaxPiNN); 217 // assert(xpi>0 && xpi<=nMaxPiNN);
218 // assert(particle1->isNucleon() && particle2- 218 // assert(particle1->isNucleon() && particle2->isNucleon());
219 219
220 G4double oldXS1Pi=CrossSectionsMultiPi 220 G4double oldXS1Pi=CrossSectionsMultiPions::NNToxPiNN(1,particle1, particle2);
221 G4double oldXS2Pi=CrossSectionsMultiPi 221 G4double oldXS2Pi=CrossSectionsMultiPions::NNToxPiNN(2,particle1, particle2);
222 G4double oldXS3Pi=CrossSectionsMultiPi 222 G4double oldXS3Pi=CrossSectionsMultiPions::NNToxPiNN(3,particle1, particle2);
223 G4double oldXS4Pi=CrossSectionsMultiPi 223 G4double oldXS4Pi=CrossSectionsMultiPions::NNToxPiNN(4,particle1, particle2);
224 G4double xsEta=CrossSectionsMultiPions 224 G4double xsEta=CrossSectionsMultiPionsAndResonances::NNToNNEta(particle1, particle2)+CrossSectionsMultiPionsAndResonances::NNToNNOmega(particle1, particle2);;
225 225
226 const G4double xs1=NNToNLK(particle1, 226 const G4double xs1=NNToNLK(particle1, particle2);
227 const G4double xs2=NNToNSK(particle1, 227 const G4double xs2=NNToNSK(particle1, particle2);
228 const G4double xs3=NNToNLKpi(particle1 228 const G4double xs3=NNToNLKpi(particle1, particle2);
229 const G4double xs4=NNToNSKpi(particle1 229 const G4double xs4=NNToNSKpi(particle1, particle2);
230 const G4double xs5=NNToNLK2pi(particle 230 const G4double xs5=NNToNLK2pi(particle1, particle2);
231 const G4double xs6=NNToNSK2pi(particle 231 const G4double xs6=NNToNSK2pi(particle1, particle2);
232 const G4double xs7=NNToNNKKb(particle1 232 const G4double xs7=NNToNNKKb(particle1, particle2);
233 const G4double xs8=NNToMissingStrangen 233 const G4double xs8=NNToMissingStrangeness(particle1, particle2);
234 const G4double xs0 = xs1 + xs2 + xs3 + 234 const G4double xs0 = xs1 + xs2 + xs3 + xs4 + xs5 + xs6 + xs7 + xs8;
235 G4double newXS1Pi=0.; 235 G4double newXS1Pi=0.;
236 G4double newXS2Pi=0.; 236 G4double newXS2Pi=0.;
237 G4double newXS3Pi=0.; 237 G4double newXS3Pi=0.;
238 G4double newXS4Pi=0.; 238 G4double newXS4Pi=0.;
239 239
240 if (xpi == 1) { 240 if (xpi == 1) {
241 if (oldXS4Pi != 0. || oldXS3Pi != 241 if (oldXS4Pi != 0. || oldXS3Pi != 0.)
242 newXS1Pi=oldXS1Pi; 242 newXS1Pi=oldXS1Pi;
243 else if (oldXS2Pi != 0.) { 243 else if (oldXS2Pi != 0.) {
244 newXS2Pi=oldXS2Pi-xsEta-xs0; 244 newXS2Pi=oldXS2Pi-xsEta-xs0;
245 if (newXS2Pi < 0.) 245 if (newXS2Pi < 0.)
246 newXS1Pi=oldXS1Pi-(xsEta+x 246 newXS1Pi=oldXS1Pi-(xsEta+xs0-oldXS2Pi);
247 else 247 else
248 newXS1Pi=oldXS1Pi; 248 newXS1Pi=oldXS1Pi;
249 } 249 }
250 else 250 else
251 newXS1Pi=oldXS1Pi-xsEta-xs0; 251 newXS1Pi=oldXS1Pi-xsEta-xs0;
252 return newXS1Pi; 252 return newXS1Pi;
253 } 253 }
254 else if (xpi == 2) { 254 else if (xpi == 2) {
255 if (oldXS4Pi != 0.){ 255 if (oldXS4Pi != 0.){
256 newXS2Pi=oldXS2Pi; 256 newXS2Pi=oldXS2Pi;
257 } 257 }
258 else if (oldXS3Pi != 0.) { 258 else if (oldXS3Pi != 0.) {
259 newXS3Pi=oldXS3Pi-xsEta-xs0; 259 newXS3Pi=oldXS3Pi-xsEta-xs0;
260 if (newXS3Pi < 0.) 260 if (newXS3Pi < 0.)
261 newXS2Pi = oldXS2Pi-(xsEta 261 newXS2Pi = oldXS2Pi-(xsEta+xs0-oldXS3Pi);
262 else 262 else
263 newXS2Pi = oldXS2Pi; 263 newXS2Pi = oldXS2Pi;
264 } 264 }
265 else { 265 else {
266 newXS2Pi = oldXS2Pi-xsEta-xs0; 266 newXS2Pi = oldXS2Pi-xsEta-xs0;
267 if (newXS2Pi < 0.) 267 if (newXS2Pi < 0.)
268 newXS2Pi = 0.; 268 newXS2Pi = 0.;
269 } 269 }
270 return newXS2Pi; 270 return newXS2Pi;
271 } 271 }
272 else if (xpi == 3) { 272 else if (xpi == 3) {
273 if (oldXS4Pi != 0.) { 273 if (oldXS4Pi != 0.) {
274 newXS4Pi=oldXS4Pi-xsEta-xs0; 274 newXS4Pi=oldXS4Pi-xsEta-xs0;
275 if (newXS4Pi < 0.) 275 if (newXS4Pi < 0.)
276 newXS3Pi=oldXS3Pi-(xsEta+xs0- 276 newXS3Pi=oldXS3Pi-(xsEta+xs0-oldXS4Pi);
277 else 277 else
278 newXS3Pi=oldXS3Pi; 278 newXS3Pi=oldXS3Pi;
279 } 279 }
280 else { 280 else {
281 newXS3Pi=oldXS3Pi-xsEta-xs0; 281 newXS3Pi=oldXS3Pi-xsEta-xs0;
282 if (newXS3Pi < 0.) 282 if (newXS3Pi < 0.)
283 newXS3Pi=0.; 283 newXS3Pi=0.;
284 } 284 }
285 return newXS3Pi; 285 return newXS3Pi;
286 } 286 }
287 else if (xpi == 4) { 287 else if (xpi == 4) {
288 newXS4Pi=oldXS4Pi-xsEta-xs0; 288 newXS4Pi=oldXS4Pi-xsEta-xs0;
289 if (newXS4Pi < 0.) 289 if (newXS4Pi < 0.)
290 newXS4Pi=0.; 290 newXS4Pi=0.;
291 return newXS4Pi; 291 return newXS4Pi;
292 } 292 }
293 293
294 else // should never reach this point 294 else // should never reach this point
295 return 0.; 295 return 0.;
296 } 296 }
297 297
298 /// \brief elastic cross sections 298 /// \brief elastic cross sections
299 299
300 G4double CrossSectionsStrangeness::NYelast 300 G4double CrossSectionsStrangeness::NYelastic(Particle const * const p1, Particle const * const p2) {
301 // 301 //
302 // Hyperon-Nucleon elastic cross 302 // Hyperon-Nucleon elastic cross sections
303 // 303 //
304 // assert((p1->isNucleon() && p2->isHyperon()) 304 // assert((p1->isNucleon() && p2->isHyperon()) || (p1->isHyperon() && p2->isNucleon()));
305 305
306 G4double sigma = 0.; 306 G4double sigma = 0.;
307 307
308 const Particle *hyperon; 308 const Particle *hyperon;
309 const Particle *nucleon; 309 const Particle *nucleon;
310 310
311 if (p1->isHyperon()) { 311 if (p1->isHyperon()) {
312 hyperon = p1; 312 hyperon = p1;
313 nucleon = p2; 313 nucleon = p2;
314 } 314 }
315 else { 315 else {
316 hyperon = p2; 316 hyperon = p2;
317 nucleon = p1; 317 nucleon = p1;
318 } 318 }
319 319
320 const G4double pLab = KinematicsUtils: 320 const G4double pLab = KinematicsUtils::momentumInLab(hyperon, nucleon); // MeV
321 321
322 if (pLab < 145.) 322 if (pLab < 145.)
323 sigma = 200.; 323 sigma = 200.;
324 else if (pLab < 425.) 324 else if (pLab < 425.)
325 sigma = 869.*std::exp(-pLab/100.); 325 sigma = 869.*std::exp(-pLab/100.);
326 else if (pLab < 30000.) 326 else if (pLab < 30000.)
327 sigma = 12.8*std::exp(-6.2e-5*pLab 327 sigma = 12.8*std::exp(-6.2e-5*pLab);
328 else 328 else
329 sigma=0.; 329 sigma=0.;
330 330
331 if (sigma < 0.) sigma = 0.; // should 331 if (sigma < 0.) sigma = 0.; // should never happen
332 return sigma; 332 return sigma;
333 } 333 }
334 334
335 G4double CrossSectionsStrangeness::NKelast 335 G4double CrossSectionsStrangeness::NKelastic(Particle const * const p1, Particle const * const p2) {
336 // 336 //
337 // Kaon-Nucleon elastic cross sec 337 // Kaon-Nucleon elastic cross sections
338 // 338 //
339 // assert((p1->isNucleon() && p2->isKaon()) || 339 // assert((p1->isNucleon() && p2->isKaon()) || (p1->isKaon() && p2->isNucleon()));
340 340
341 G4double sigma=0.; 341 G4double sigma=0.;
342 342
343 const Particle *kaon; 343 const Particle *kaon;
344 const Particle *nucleon; 344 const Particle *nucleon;
345 345
346 if (p1->isKaon()) { 346 if (p1->isKaon()) {
347 kaon = p1; 347 kaon = p1;
348 nucleon = p2; 348 nucleon = p2;
349 } 349 }
350 else { 350 else {
351 kaon = p2; 351 kaon = p2;
352 nucleon = p1; 352 nucleon = p1;
353 } 353 }
354 354
355 const G4double pLab = KinematicsUtils: 355 const G4double pLab = KinematicsUtils::momentumInLab(kaon, nucleon); // MeV
356 356
357 if (pLab < 935.) 357 if (pLab < 935.)
358 sigma = 12.; 358 sigma = 12.;
359 else if (pLab < 2080.) 359 else if (pLab < 2080.)
360 sigma = 17.4-3.*std::exp(6.3e-4*pL 360 sigma = 17.4-3.*std::exp(6.3e-4*pLab);
361 else if (pLab < 5500.) 361 else if (pLab < 5500.)
362 sigma = 832.*std::pow(pLab,-0.64); 362 sigma = 832.*std::pow(pLab,-0.64);
363 else if (pLab < 30000.) 363 else if (pLab < 30000.)
364 sigma = 3.36; 364 sigma = 3.36;
365 else 365 else
366 sigma=0.; 366 sigma=0.;
367 367
368 if (sigma < 0.) sigma = 0.; // should 368 if (sigma < 0.) sigma = 0.; // should never happen
369 return sigma; 369 return sigma;
370 } 370 }
371 371
372 G4double CrossSectionsStrangeness::NKbelas 372 G4double CrossSectionsStrangeness::NKbelastic(Particle const * const p1, Particle const * const p2) {
373 // 373 //
374 // antiKaon-Nucleon elastic cross 374 // antiKaon-Nucleon elastic cross sections
375 // 375 //
376 // assert((p1->isNucleon() && p2->isAntiKaon() 376 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
377 377
378 G4double sigma=0.; 378 G4double sigma=0.;
379 379
380 const Particle *antikaon; 380 const Particle *antikaon;
381 const Particle *nucleon; 381 const Particle *nucleon;
382 382
383 if (p1->isAntiKaon()) { 383 if (p1->isAntiKaon()) {
384 antikaon = p1; 384 antikaon = p1;
385 nucleon = p2; 385 nucleon = p2;
386 } 386 }
387 else { 387 else {
388 antikaon = p2; 388 antikaon = p2;
389 nucleon = p1; 389 nucleon = p1;
390 } 390 }
391 391
392 const G4double pLab = 0.001*Kinematics 392 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
393 393
394 if(pLab > 1E-6) // sigma = 287.823 [mb 394 if(pLab > 1E-6) // sigma = 287.823 [mb] -> rise very slowly, not cut needed
395 sigma = 6.132*std::pow(pLab,-0.2437 395 sigma = 6.132*std::pow(pLab,-0.2437)+12.98*std::exp(-std::pow(pLab-0.9902,2)/0.05558)+2.928*std::exp(-std::pow(pLab-1.649,2)/0.772)+564.3*std::exp(-std::pow(pLab+0.9901,2)/0.5995);
396 396
397 if (sigma < 0.) sigma = 0.; // should 397 if (sigma < 0.) sigma = 0.; // should never happen
398 return sigma; 398 return sigma;
399 } 399 }
400 400
401 /// \brief NN to strange cross sections 401 /// \brief NN to strange cross sections
402 402
403 G4double CrossSectionsStrangeness::NNToNLK 403 G4double CrossSectionsStrangeness::NNToNLK(Particle const * const p1, Particle const * const p2) {
404 // 404 //
405 // Nucleon-Nucleon producing N-La 405 // Nucleon-Nucleon producing N-Lambda-Kaon cross sections
406 // 406 //
407 // ratio 407 // ratio
408 // p p (1) p n (1) 408 // p p (1) p n (1)
409 // 409 //
410 // p p -> p L K+ (1) 410 // p p -> p L K+ (1)
411 // p n -> p L K0 (1/2) 411 // p n -> p L K0 (1/2)
412 // p n -> n L K+ (1/2) 412 // p n -> n L K+ (1/2)
413 // assert(p1->isNucleon() && p2->isNucleon()); 413 // assert(p1->isNucleon() && p2->isNucleon());
414 414
415 const Particle *particle1; 415 const Particle *particle1;
416 const Particle *particle2; 416 const Particle *particle2;
417 417
418 if(p2->getType() == Proton && p1->getT 418 if(p2->getType() == Proton && p1->getType() == Neutron){
419 particle1 = p2; 419 particle1 = p2;
420 particle2 = p1; 420 particle2 = p1;
421 } 421 }
422 else{ 422 else{
423 particle1 = p1; 423 particle1 = p1;
424 particle2 = p2; 424 particle2 = p2;
425 } 425 }
426 426
427 G4double sigma = 0.; 427 G4double sigma = 0.;
428 428
429 const G4double pLab = 0.001 * Kinemati 429 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1, particle2); // GeV
430 430
431 if(particle2->getType() == Proton){ 431 if(particle2->getType() == Proton){
432 if(pLab < 2.3393) return 0.; 432 if(pLab < 2.3393) return 0.;
433 else if (pLab < 30.) sigma = 1.118 433 else if (pLab < 30.) sigma = 1.11875*std::pow((pLab-2.3393),1.0951)/std::pow((pLab+2.3393),2.0958); // pLab = 30 Gev -> exess of energie = 5 GeV
434 else return 0.; 434 else return 0.;
435 } 435 }
436 else{ 436 else{
437 if(pLab < 2.3508) return 0.; 437 if(pLab < 2.3508) return 0.;
438 else if (pLab < 30.) sigma = 1.118 438 else if (pLab < 30.) sigma = 1.11875*std::pow((pLab-2.3508),1.0951)/std::pow((pLab+2.3508),2.0958);
439 else return 0.; 439 else return 0.;
440 } 440 }
441 441
442 return sigma; 442 return sigma;
443 } 443 }
444 444
445 G4double CrossSectionsStrangeness::NNToNSK 445 G4double CrossSectionsStrangeness::NNToNSK(Particle const * const p1, Particle const * const p2) {
446 // 446 //
447 // Nucleon-Nucleon producing N-Si 447 // Nucleon-Nucleon producing N-Sigma-Kaon cross sections
448 // 448 //
449 // Meson symmetry 449 // Meson symmetry
450 // pp->pS+K0 (1/4) 450 // pp->pS+K0 (1/4)
451 // pp->pS0K+ (1/8) // HEM 451 // pp->pS0K+ (1/8) // HEM
452 // pp->pS0K+ (1/4) // Data 452 // pp->pS0K+ (1/4) // Data
453 // pp->nS+K+ (1) 453 // pp->nS+K+ (1)
454 // 454 //
455 // pn->nS+K0 (1/4) 455 // pn->nS+K0 (1/4)
456 // pn->pS-K+ (1/4) 456 // pn->pS-K+ (1/4)
457 // pn->nS0K+ (5/8) 457 // pn->nS0K+ (5/8)
458 // pn->pS0K0 (5/8) 458 // pn->pS0K0 (5/8)
459 // 459 //
460 // assert(p1->isNucleon() && p2->isNucleon()); 460 // assert(p1->isNucleon() && p2->isNucleon());
461 461
462 const Particle *particle1; 462 const Particle *particle1;
463 const Particle *particle2; 463 const Particle *particle2;
464 464
465 if(p2->getType() == Proton && p1->getT 465 if(p2->getType() == Proton && p1->getType() == Neutron){
466 particle1 = p2; 466 particle1 = p2;
467 particle2 = p1; 467 particle2 = p1;
468 } 468 }
469 else{ 469 else{
470 particle1 = p1; 470 particle1 = p1;
471 particle2 = p2; 471 particle2 = p2;
472 } 472 }
473 473
474 G4double sigma = 0.; 474 G4double sigma = 0.;
475 475
476 const G4double pLab = 0.001 * Kinemati 476 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1, particle2); // GeV
477 477
478 if(pLab < 2.593) 478 if(pLab < 2.593)
479 return 0.; 479 return 0.;
480 480
481 if(p2->getType() == p1->getType()) 481 if(p2->getType() == p1->getType())
482 // sigma = 1.375*2*6.38*std::pow(pL 482 // sigma = 1.375*2*6.38*std::pow(pLab-2.57,2.1)/std::pow(pLab,4.162);
483 sigma = 1.5*6.38*std::pow(pLab-2.5 483 sigma = 1.5*6.38*std::pow(pLab-2.593,2.1)/std::pow(pLab,4.162);
484 else 484 else
485 sigma = 1.75*6.38*std::pow(pLab-2. 485 sigma = 1.75*6.38*std::pow(pLab-2.593,2.1)/std::pow(pLab,4.162);
486 486
487 487
488 return sigma; 488 return sigma;
489 } 489 }
490 490
491 G4double CrossSectionsStrangeness::NNToNLK 491 G4double CrossSectionsStrangeness::NNToNLKpi(Particle const * const p1, Particle const * const p2) {
492 // 492 //
493 // Nucleon-Nucleon producing N-La 493 // Nucleon-Nucleon producing N-Lambda-Kaon-pion cross sections
494 // 494 //
495 // ratio (pure NN -> DLK) 495 // ratio (pure NN -> DLK)
496 // pp (12) pn (8) 496 // pp (12) pn (8)
497 // 497 //
498 // pp -> p pi+ L K0 (9)(3) 498 // pp -> p pi+ L K0 (9)(3)
499 // pp -> p pi0 L K+ (2)(1*2/3) 499 // pp -> p pi0 L K+ (2)(1*2/3)
500 // pp -> n pi+ L K+ (1)(1*1/3) 500 // pp -> n pi+ L K+ (1)(1*1/3)
501 // 501 //
502 // pn -> p pi- L K+ (2)(2*1/3) 502 // pn -> p pi- L K+ (2)(2*1/3)
503 // pn -> n pi0 L K+ (4)(2*2/3) 503 // pn -> n pi0 L K+ (4)(2*2/3)
504 // pn -> p pi0 L K0 (4) 504 // pn -> p pi0 L K0 (4)
505 // pn -> n pi+ L K0 (2) 505 // pn -> n pi+ L K0 (2)
506 506
507 // assert(p1->isNucleon() && p2->isNucleon()); 507 // assert(p1->isNucleon() && p2->isNucleon());
508 508
509 G4double sigma = 0.; 509 G4double sigma = 0.;
510 G4double ratio = 0.; 510 G4double ratio = 0.;
511 G4double ratio1 = 0.; 511 G4double ratio1 = 0.;
512 G4double ratio2 = 0.; 512 G4double ratio2 = 0.;
513 const G4double ener = KinematicsUtils: 513 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 540.;
514 if( ener < p1->getMass() + p2->getMass 514 if( ener < p1->getMass() + p2->getMass())
515 return 0; 515 return 0;
516 const G4int iso = ParticleTable::getIs 516 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
517 517
518 const G4double xsiso2 = CrossSectionsM 518 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
519 if (iso != 0){ 519 if (iso != 0){
520 ratio1 = CrossSectionsMultiPions:: 520 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
521 ratio2 = CrossSectionsMultiPions:: 521 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
522 } 522 }
523 else { 523 else {
524 const G4double xsiso0 = CrossSecti 524 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
525 ratio1 = 0.5*(CrossSectionsMultiPi 525 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
526 ratio2 = 0.5*(CrossSectionsMultiPi 526 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
527 } 527 }
528 528
529 if( ratio1 == 0 || ratio2 == 0) 529 if( ratio1 == 0 || ratio2 == 0)
530 return 0.; 530 return 0.;
531 531
532 ratio = ratio2/ratio1; 532 ratio = ratio2/ratio1;
533 533
534 sigma = ratio * NNToNLK(p1,p2) * 3; 534 sigma = ratio * NNToNLK(p1,p2) * 3;
535 535
536 /* const G4double pLab = 0.001 * Kinema 536 /* const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1, p2); // GeV
537 if(pLab <= 2.77) return 0.; 537 if(pLab <= 2.77) return 0.;
538 sigma = 0.4 * std::pow(pLab-2.77,1.603 538 sigma = 0.4 * std::pow(pLab-2.77,1.603)/std::pow(pLab,1.492);*/
539 539
540 return sigma; 540 return sigma;
541 } 541 }
542 542
543 G4double CrossSectionsStrangeness::NNToNSK 543 G4double CrossSectionsStrangeness::NNToNSKpi(Particle const * const p1, Particle const * const p2) {
544 // 544 //
545 // Nucleon-Nucleon producing N-Sigma-K 545 // Nucleon-Nucleon producing N-Sigma-Kaon-pion cross sections
546 // 546 //
547 // ratio (pure NN -> DSK) 547 // ratio (pure NN -> DSK)
548 // pp (36) pn (36) 548 // pp (36) pn (36)
549 // 549 //
550 // pp -> p pi+ S- K+ (9) 550 // pp -> p pi+ S- K+ (9)
551 // pp -> p pi+ S0 K0 (9) 551 // pp -> p pi+ S0 K0 (9)
552 // pp -> p pi0 S+ K0 (4) 552 // pp -> p pi0 S+ K0 (4)
553 // pp -> n pi+ S+ K0 (2) 553 // pp -> n pi+ S+ K0 (2)
554 // pp -> p pi0 S0 K+ (4) 554 // pp -> p pi0 S0 K+ (4)
555 // pp -> n pi+ S0 K+ (2) 555 // pp -> n pi+ S0 K+ (2)
556 // pp -> p pi- S+ K+ (2) 556 // pp -> p pi- S+ K+ (2)
557 // pp -> n pi0 S+ K+ (4) 557 // pp -> n pi0 S+ K+ (4)
558 558
559 // pn -> p pi0 S- K+ (4) 559 // pn -> p pi0 S- K+ (4)
560 // pn -> n pi+ S- K+ (2) 560 // pn -> n pi+ S- K+ (2)
561 // pn -> p pi0 S0 K0 (2) 561 // pn -> p pi0 S0 K0 (2)
562 // pn -> n pi+ S0 K0 (1) 562 // pn -> n pi+ S0 K0 (1)
563 // pn -> p pi+ S- K0 (9) 563 // pn -> p pi+ S- K0 (9)
564 564
565 // assert(p1->isNucleon() && p2->isNucleon()); 565 // assert(p1->isNucleon() && p2->isNucleon());
566 566
567 G4double sigma = 0.; 567 G4double sigma = 0.;
568 G4double ratio = 0.; 568 G4double ratio = 0.;
569 G4double ratio1 = 0.; 569 G4double ratio1 = 0.;
570 G4double ratio2 = 0.; 570 G4double ratio2 = 0.;
571 const G4double ener = KinematicsUtils: 571 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 620.;
572 if( ener < p1->getMass() + p2->getMass 572 if( ener < p1->getMass() + p2->getMass())
573 return 0; 573 return 0;
574 const G4int iso = ParticleTable::getIs 574 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
575 575
576 const G4double xsiso2 = CrossSectionsM 576 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
577 if (iso != 0){ 577 if (iso != 0){
578 ratio1 = CrossSectionsMultiPions:: 578 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
579 ratio2 = CrossSectionsMultiPions:: 579 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
580 } 580 }
581 else { 581 else {
582 const G4double xsiso0 = CrossSecti 582 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
583 ratio1 = 0.5*(CrossSectionsMultiPi 583 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
584 ratio2 = 0.5*(CrossSectionsMultiPi 584 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
585 } 585 }
586 586
587 if( ratio1 == 0 || ratio2 == 0) 587 if( ratio1 == 0 || ratio2 == 0)
588 return 0.; 588 return 0.;
589 589
590 ratio = ratio2/ratio1; 590 ratio = ratio2/ratio1;
591 591
592 sigma = ratio * NNToNSK(p1,p2) * 3; 592 sigma = ratio * NNToNSK(p1,p2) * 3;
593 593
594 return sigma; 594 return sigma;
595 } 595 }
596 596
597 G4double CrossSectionsStrangeness::NNToNLK 597 G4double CrossSectionsStrangeness::NNToNLK2pi(Particle const * const p1, Particle const * const p2) {
598 // 598 //
599 // Nucleon-Nucleon producing N-La 599 // Nucleon-Nucleon producing N-Lambda-Kaon-pion cross sections
600 // 600 //
601 // assert(p1->isNucleon() && p2->isNucleon()); 601 // assert(p1->isNucleon() && p2->isNucleon());
602 602
603 G4double sigma = 0.; 603 G4double sigma = 0.;
604 G4double ratio = 0.; 604 G4double ratio = 0.;
605 G4double ratio1 = 0.; 605 G4double ratio1 = 0.;
606 G4double ratio2 = 0.; 606 G4double ratio2 = 0.;
607 const G4double ener = KinematicsUtils: 607 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 675.;
608 if( ener < p1->getMass() + p2->getMass 608 if( ener < p1->getMass() + p2->getMass())
609 return 0; 609 return 0;
610 const G4int iso = ParticleTable::getIs 610 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
611 611
612 const G4double xsiso2 = CrossSectionsM 612 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
613 if (iso != 0){ 613 if (iso != 0){
614 ratio1 = CrossSectionsMultiPions:: 614 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
615 ratio2 = CrossSectionsMultiPions:: 615 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
616 } 616 }
617 else { 617 else {
618 const G4double xsiso0 = CrossSecti 618 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
619 ratio1 = 0.5*(CrossSectionsMultiPi 619 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
620 ratio2 = 0.5*(CrossSectionsMultiPi 620 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
621 } 621 }
622 622
623 if( ratio1 == 0 || ratio2 == 0) 623 if( ratio1 == 0 || ratio2 == 0)
624 return 0.; 624 return 0.;
625 625
626 ratio = ratio2/ratio1; 626 ratio = ratio2/ratio1;
627 627
628 sigma = ratio * NNToNLKpi(p1,p2); 628 sigma = ratio * NNToNLKpi(p1,p2);
629 629
630 return sigma; 630 return sigma;
631 } 631 }
632 632
633 G4double CrossSectionsStrangeness::NNToNSK 633 G4double CrossSectionsStrangeness::NNToNSK2pi(Particle const * const p1, Particle const * const p2) {
634 // 634 //
635 // Nucleon-Nucleon producing N-Si 635 // Nucleon-Nucleon producing N-Sigma-Kaon-pion cross sections
636 // 636 //
637 // assert(p1->isNucleon() && p2->isNucleon()); 637 // assert(p1->isNucleon() && p2->isNucleon());
638 638
639 G4double sigma = 0.; 639 G4double sigma = 0.;
640 G4double ratio = 0.; 640 G4double ratio = 0.;
641 G4double ratio1 = 0.; 641 G4double ratio1 = 0.;
642 G4double ratio2 = 0.; 642 G4double ratio2 = 0.;
643 const G4double ener = KinematicsUtils: 643 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 755.;
644 if( ener < p1->getMass() + p2->getMass 644 if( ener < p1->getMass() + p2->getMass())
645 return 0; 645 return 0;
646 const G4int iso = ParticleTable::getIs 646 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
647 647
648 const G4double xsiso2 = CrossSectionsM 648 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
649 if (iso != 0){ 649 if (iso != 0){
650 ratio1 = CrossSectionsMultiPions:: 650 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
651 ratio2 = CrossSectionsMultiPions:: 651 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
652 } 652 }
653 else { 653 else {
654 const G4double xsiso0 = CrossSecti 654 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
655 ratio1 = 0.5*(CrossSectionsMultiPi 655 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
656 ratio2 = 0.5*(CrossSectionsMultiPi 656 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
657 } 657 }
658 658
659 if( ratio1 == 0 || ratio2 == 0) 659 if( ratio1 == 0 || ratio2 == 0)
660 return 0.; 660 return 0.;
661 661
662 ratio = ratio2/ratio1; 662 ratio = ratio2/ratio1;
663 663
664 sigma = ratio * NNToNSKpi(p1,p2); 664 sigma = ratio * NNToNSKpi(p1,p2);
665 665
666 return sigma; 666 return sigma;
667 } 667 }
668 668
669 G4double CrossSectionsStrangeness::NNToNNK 669 G4double CrossSectionsStrangeness::NNToNNKKb(Particle const * const p1, Particle const * const p2) {
670 // 670 //
671 // Nucleon-Nucleon producing Nucl 671 // Nucleon-Nucleon producing Nucleon-Nucleon-Kaon-antiKaon cross sections
672 // 672 //
673 // Channel strongly resonant; fit from 673 // Channel strongly resonant; fit from Sibirtesev - Z. Phys. A 358, 101-106 (1997) (eq.21)
674 // ratio 674 // ratio
675 // pp (6) pn (13)*2 675 // pp (6) pn (13)*2
676 // pp -> pp K+ K- (1) 676 // pp -> pp K+ K- (1)
677 // pp -> pp K0 K0 (1) 677 // pp -> pp K0 K0 (1)
678 // pp -> pn K+ K0 (4) 678 // pp -> pn K+ K0 (4)
679 // pn -> pp K0 K- (4) 679 // pn -> pp K0 K- (4)
680 // pn -> pn K+ K- (9) 680 // pn -> pn K+ K- (9)
681 // 681 //
682 682
683 // assert(p1->isNucleon() && p2->isNucleon()); 683 // assert(p1->isNucleon() && p2->isNucleon());
684 684
685 G4double sigma = 0.; 685 G4double sigma = 0.;
686 const G4int iso = ParticleTable::getIs 686 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
687 const G4double ener = 0.001*Kinematics 687 const G4double ener = 0.001*KinematicsUtils::totalEnergyInCM(p1, p2); // GeV
688 688
689 if(ener < 2.872) 689 if(ener < 2.872)
690 return 0.; 690 return 0.;
691 691
692 if(iso == 0) 692 if(iso == 0)
693 sigma = 26 * 5./19. * 0.3 *std::po 693 sigma = 26 * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
694 else 694 else
695 sigma = 6 * 5./19. * 0.3 *std::pow 695 sigma = 6 * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
696 696
697 return sigma; 697 return sigma;
698 } 698 }
699 699
700 G4double CrossSectionsStrangeness::NNToMis 700 G4double CrossSectionsStrangeness::NNToMissingStrangeness(Particle const * const p1, Particle const * const p2) {
701 // 701 //
702 // Nucleon-Nucleon missing strang 702 // Nucleon-Nucleon missing strangeness production cross sections
703 // 703 //
704 // assert(p1->isNucleon() && p2->isNucleon()); 704 // assert(p1->isNucleon() && p2->isNucleon());
705 705
706 G4double sigma = 0.; 706 G4double sigma = 0.;
707 707
708 const G4double pLab = 0.001 * Kinemati 708 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
709 const G4int iso = ParticleTable::getIs 709 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
710 710
711 if(pLab < 6.) return 0.; 711 if(pLab < 6.) return 0.;
712 712
713 if(iso == 0){ 713 if(iso == 0){
714 if(pLab < 30.) sigma = 10.15*std:: 714 if(pLab < 30.) sigma = 10.15*std::pow((pLab - 6.),2.157)/std::pow(pLab,2.333);
715 else return 0.; 715 else return 0.;
716 } 716 }
717 else{ 717 else{
718 if(pLab < 30.) sigma = 8.12*std::p 718 if(pLab < 30.) sigma = 8.12*std::pow((pLab - 6.),2.157)/std::pow(pLab,2.333);
719 else return 0.; 719 else return 0.;
720 } 720 }
721 return sigma; 721 return sigma;
722 } 722 }
723 723
724 /** \brief NDelta to strange cross sections 724 /** \brief NDelta to strange cross sections
725 * 725 *
726 * No experimental data 726 * No experimental data
727 * Parametrization from Phys.Rev.C 59 1 (369 727 * Parametrization from Phys.Rev.C 59 1 (369) (1999)
728 * 728 *
729 * Correction are applied on the isospin sym 729 * Correction are applied on the isospin symetry provided in the paper
730 */ 730 */
731 731
732 G4double CrossSectionsStrangeness::NDeltaT 732 G4double CrossSectionsStrangeness::NDeltaToNLK(Particle const * const p1, Particle const * const p2) {
733 // Nucleon-Delta producing Nucleon Lam 733 // Nucleon-Delta producing Nucleon Lambda Kaon cross section
734 // 734 //
735 // XS from K. Tsushima, A. Sibirtsev, 735 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
736 // 736 //
737 // ratio <<
738 // D++ n -> p L K+ (3) 737 // D++ n -> p L K+ (3)
739 // 738 //
740 // D+ p -> p L K+ (1) 739 // D+ p -> p L K+ (1)
741 // 740 //
742 // D+ n -> p L K0 (1) 741 // D+ n -> p L K0 (1)
743 // D+ n -> n L K+ (1) 742 // D+ n -> n L K+ (1)
744 << 743 //return 0.;
745 G4double a = 4.169; <<
746 G4double b = 2.227; <<
747 G4double c = 2.511; <<
748 G4double n_channel = 4.; // number of <<
749 <<
750 // assert((p1->isNucleon() && p2->isResonance( 744 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
751 745
752 const G4int iso = ParticleTable::getIs 746 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
753 if(std::abs(iso) == 4) return 0.; 747 if(std::abs(iso) == 4) return 0.;
754 748
755 G4double sigma = 0.; 749 G4double sigma = 0.;
756 750
757 const G4double s = KinematicsUtils::sq 751 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // MeV^2
758 const G4double s0 = 6.511E6; // MeV^2 752 const G4double s0 = 6.511E6; // MeV^2
759 753
760 if(s <= s0) return 0.; 754 if(s <= s0) return 0.;
761 755
762 sigma = n_channel*a*std::pow(s/s0-1,b) << 756 sigma = 4.*4.169*std::pow(s/s0-1,2.227)*std::pow(s0/s,2.511);
763 757
764 //const G4double pLab = sdt::sqrt(s*s/ 758 //const G4double pLab = sdt::sqrt(s*s/(4*ParticleTable::effectiveNucleonMass2)-s)*0.001;
765 //sigma = 3*1.11875*std::pow((pLab-2.3 759 //sigma = 3*1.11875*std::pow((pLab-2.3508),1.0951)/std::pow((pLab+2.3508),2.0958); // NDelta sim to NN
766 760
767 if(iso == 0){ // D+ n 761 if(iso == 0){ // D+ n
768 sigma *= 2./6.; 762 sigma *= 2./6.;
769 } 763 }
770 else if (ParticleTable::getIsospin(p1- 764 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType())){// D+ p
771 sigma *= 1./6.; 765 sigma *= 1./6.;
772 } 766 }
773 else{// D++ n 767 else{// D++ n
774 sigma *= 3./6.; 768 sigma *= 3./6.;
775 } 769 }
776 return sigma; 770 return sigma;
777 } 771 }
778 G4double CrossSectionsStrangeness::NDeltaT 772 G4double CrossSectionsStrangeness::NDeltaToNSK(Particle const * const p1, Particle const * const p2) {
779 // Nucleon-Delta producing Nucleon Sig 773 // Nucleon-Delta producing Nucleon Sigma Kaon cross section
780 // 774 //
781 // XS from K. Tsushima, A. Sibirtsev, 775 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369 ( X 1.25 (124/99) for isospin consideration)
782 // 776 //
783 // ratio <<
784 // D++ p -> p S+ K+ (6) 777 // D++ p -> p S+ K+ (6)
785 // 778 //
786 // D++ n -> p S+ K0 (3) **** 779 // D++ n -> p S+ K0 (3) ****
787 // D++ n -> p S0 K+ (3) 780 // D++ n -> p S0 K+ (3)
788 // D++ n -> n S+ K+ (3) 781 // D++ n -> n S+ K+ (3)
789 // 782 //
790 // D+ p -> p S+ K0 (2) 783 // D+ p -> p S+ K0 (2)
791 // D+ p -> p S0 K+ (2) 784 // D+ p -> p S0 K+ (2)
792 // D+ p -> n S+ K+ (3) 785 // D+ p -> n S+ K+ (3)
793 // 786 //
794 // D+ n -> p S0 K0 (3) 787 // D+ n -> p S0 K0 (3)
795 // D+ n -> p S- K+ (2) 788 // D+ n -> p S- K+ (2)
796 // D+ n -> n S+ K0 (2) 789 // D+ n -> n S+ K0 (2)
797 // D+ n -> n S0 K+ (2) 790 // D+ n -> n S0 K+ (2)
798 791
799 G4double a = 39.54; <<
800 G4double b = 2.799; <<
801 G4double c = 6.303; <<
802 G4double n_channel = 11.; <<
803 <<
804 // assert((p1->isNucleon() && p2->isResonance( 792 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
805 793
806 G4double sigma = 0.; 794 G4double sigma = 0.;
807 795
808 const G4double s = KinematicsUtils::sq 796 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
809 const G4double s0 = 6.935E6; // Mev^2 797 const G4double s0 = 6.935E6; // Mev^2
810 const G4int iso = ParticleTable::getIs 798 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
811 799
812 if(s <= s0) 800 if(s <= s0)
813 return 0.; 801 return 0.;
814 802
815 sigma = n_channel*a*std::pow(s/s0-1,b) << 803 sigma = 11.*39.54*std::pow(s/s0-1,2.799)*std::pow(s0/s,6.303);
816 804
817 //const G4double pLab = sdt::sqrt(s*s/ 805 //const G4double pLab = sdt::sqrt(s*s/(4*ParticleTable::effectiveNucleonMass2)-s)*0.001;
818 //sigma = 22./12./2. * 4.75*6.38*std:: 806 //sigma = 22./12./2. * 4.75*6.38*std::pow(pLab-2.593,2.1)/std::pow(pLab,4.162); // NDelta sim to NN
819 807
820 if(iso == 0)// D+ n 808 if(iso == 0)// D+ n
821 sigma *= 9./31.; 809 sigma *= 9./31.;
822 else if (ParticleTable::getIsospin(p1- 810 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
823 sigma *= 7./31.; 811 sigma *= 7./31.;
824 else if (std::abs(iso) == 2)// D++ n 812 else if (std::abs(iso) == 2)// D++ n
825 sigma *= 9./31.; 813 sigma *= 9./31.;
826 else // D++ p 814 else // D++ p
827 sigma *= 6./31.; 815 sigma *= 6./31.;
828 816
829 return sigma; 817 return sigma;
830 } 818 }
831 G4double CrossSectionsStrangeness::NDeltaT 819 G4double CrossSectionsStrangeness::NDeltaToDeltaLK(Particle const * const p1, Particle const * const p2) {
832 // Nucleon-Delta producing Delta Lambd 820 // Nucleon-Delta producing Delta Lambda Kaon cross section
833 // 821 //
834 // XS from K. Tsushima, A. Sibirtsev, 822 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
835 // 823 //
836 // ratio 824 // ratio
837 // D++ p -> L K+ D++ (4) 825 // D++ p -> L K+ D++ (4)
838 // 826 //
839 // D++ n -> L K+ D+ (3) 827 // D++ n -> L K+ D+ (3)
840 // D++ n -> L K0 D++ (4) 828 // D++ n -> L K0 D++ (4)
841 // 829 //
842 // D+ p -> L K0 D++ (3) 830 // D+ p -> L K0 D++ (3)
843 // D+ p -> L K+ D+ (2) 831 // D+ p -> L K+ D+ (2)
844 // 832 //
845 // D+ n -> L K+ D0 (4) 833 // D+ n -> L K+ D0 (4)
846 // D+ n -> L K0 D+ (2) 834 // D+ n -> L K0 D+ (2)
847 << 835 //return 0.;
848 G4double a = 2.679; <<
849 G4double b = 2.280; <<
850 G4double c = 5.086; <<
851 G4double n_channel = 7.; <<
852 836
853 // assert((p1->isNucleon() && p2->isResonance( 837 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
854 838
855 const G4double s = KinematicsUtils::sq 839 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
856 const G4double s0 = 8.096E6; // Mev^2 840 const G4double s0 = 8.096E6; // Mev^2
857 const G4int iso = ParticleTable::getIs 841 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
858 842
859 if(s <= s0) 843 if(s <= s0)
860 return 0.; 844 return 0.;
861 845
862 G4double sigma = n_channel*a*std::pow( << 846 G4double sigma = 7.*2.679*std::pow(s/s0-1,2.280)*std::pow(s0/s,5.086);
863 847
864 if(iso == 0)// D+ n 848 if(iso == 0)// D+ n
865 sigma *= 6./22.; 849 sigma *= 6./22.;
866 else if (ParticleTable::getIsospin(p1- 850 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
867 sigma *= 5./22.; 851 sigma *= 5./22.;
868 else if (std::abs(iso) == 2)// D++ n 852 else if (std::abs(iso) == 2)// D++ n
869 sigma *= 7./22.; 853 sigma *= 7./22.;
870 else // D++ p 854 else // D++ p
871 sigma *= 4./22.; 855 sigma *= 4./22.;
872 856
873 return sigma; 857 return sigma;
874 } 858 }
875 G4double CrossSectionsStrangeness::NDeltaT 859 G4double CrossSectionsStrangeness::NDeltaToDeltaSK(Particle const * const p1, Particle const * const p2) {
876 // Nucleon-Delta producing Delta Sigma 860 // Nucleon-Delta producing Delta Sigma Kaon cross section
877 // 861 //
878 // XS from K. Tsushima, A. Sibirtsev, 862 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
879 // 863 //
880 // D++ p (9) 864 // D++ p (9)
881 // D++ n (15) 865 // D++ n (15)
882 // D+ p (11) 866 // D+ p (11)
883 // D+ n (13) 867 // D+ n (13)
884 // 868 //
885 // ratio 869 // ratio
886 // D++ p -> S+ K+ D+ (a) (2) 870 // D++ p -> S+ K+ D+ (a) (2)
887 // D++ p -> S0 K+ D++ (b) (1) 871 // D++ p -> S0 K+ D++ (b) (1)
888 // D++ p -> S+ K0 D++ (c) (6) 872 // D++ p -> S+ K0 D++ (c) (6)
889 // 873 //
890 // D++ n -> S+ K+ D0 *(d)* (2) 874 // D++ n -> S+ K+ D0 *(d)* (2)
891 // D++ n -> S0 K+ D+ (e) (4) 875 // D++ n -> S0 K+ D+ (e) (4)
892 // D++ n -> S- K+ D++ (f) (6)(c)* 876 // D++ n -> S- K+ D++ (f) (6)(c)*
893 // D++ n -> S+ K0 D+ (a)* (2) 877 // D++ n -> S+ K0 D+ (a)* (2)
894 // D++ n -> S0 K0 D++ (b)* (1)* 878 // D++ n -> S0 K0 D++ (b)* (1)*
895 // 879 //
896 // D+ p -> S+ K+ D0 (i) (2)* 880 // D+ p -> S+ K+ D0 (i) (2)*
897 // D+ p -> S0 K+ D+ (j) (1) 881 // D+ p -> S0 K+ D+ (j) (1)
898 // D+ p -> S- K+ D++ (k) (2)(g=a)* 882 // D+ p -> S- K+ D++ (k) (2)(g=a)*
899 // D+ p -> S+ K0 D+ (l) (2) 883 // D+ p -> S+ K0 D+ (l) (2)
900 // D+ p -> S0 K0 D++ (m) (4)(e)* 884 // D+ p -> S0 K0 D++ (m) (4)(e)*
901 // 885 //
902 // D+ n -> S+ K+ D- *(d)* (2) 886 // D+ n -> S+ K+ D- *(d)* (2)
903 // D+ n -> S0 K+ D0 (o) (4) 887 // D+ n -> S0 K+ D0 (o) (4)
904 // D+ n -> S- K+ D+ (p) (2)* 888 // D+ n -> S- K+ D+ (p) (2)*
905 // D+ n -> S+ K0 D0 (i)* (2)* 889 // D+ n -> S+ K0 D0 (i)* (2)*
906 // D+ n -> S0 K0 D+ (j)* (1)* 890 // D+ n -> S0 K0 D+ (j)* (1)*
907 // D+ n -> S- K0 D++ (k)* (2)* 891 // D+ n -> S- K0 D++ (k)* (2)*
908 << 892 //return 0.;
909 G4double a = 8.407; <<
910 G4double b = 2.743; <<
911 G4double c = 21.18; <<
912 G4double n_channel = 19.; <<
913 893
914 // assert((p1->isNucleon() && p2->isResonance( 894 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
915 895
916 const G4double s = KinematicsUtils::sq 896 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
917 const G4double s0 = 8.568E6; // Mev^2 897 const G4double s0 = 8.568E6; // Mev^2
918 const G4int iso = ParticleTable::getIs 898 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
919 899
920 if(s <= s0) 900 if(s <= s0)
921 return 0.; 901 return 0.;
922 902
923 G4double sigma = n_channel*a*std::pow( << 903 G4double sigma = 19.*21.18*std::pow(s/s0-1,2.743)*std::pow(s0/s,8.407);
924 904
925 if(iso == 0)// D+ n 905 if(iso == 0)// D+ n
926 sigma *= 13./48.; 906 sigma *= 13./48.;
927 else if (ParticleTable::getIsospin(p1- 907 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
928 sigma *= 11./48.; 908 sigma *= 11./48.;
929 else if (std::abs(iso) == 2)// D++ n 909 else if (std::abs(iso) == 2)// D++ n
930 sigma *= 15./48.; 910 sigma *= 15./48.;
931 else // D++ p 911 else // D++ p
932 sigma *= 9./48.; 912 sigma *= 9./48.;
933 913
934 return sigma; 914 return sigma;
935 } 915 }
936 916
937 G4double CrossSectionsStrangeness::NDeltaT 917 G4double CrossSectionsStrangeness::NDeltaToNNKKb(Particle const * const p1, Particle const * const p2) {
938 // Nucleon-Delta producing Nucleon-Nuc 918 // Nucleon-Delta producing Nucleon-Nucleon Kaon antiKaon cross section
939 // 919 //
940 // Total = sigma(NN->NNKKb)*10 920 // Total = sigma(NN->NNKKb)*10
941 // 921 //
942 // D++ p (6) 922 // D++ p (6)
943 // D++ n (9) 923 // D++ n (9)
944 // D+ p (7) 924 // D+ p (7)
945 // D+ n (8) 925 // D+ n (8)
946 // 926 //
947 // ratio 927 // ratio
948 // D++ p -> p p K+ K0b (6) 928 // D++ p -> p p K+ K0b (6)
949 // 929 //
950 // D++ n -> p p K+ K- (3) 930 // D++ n -> p p K+ K- (3)
951 // D++ n -> p p K0 K0b (3) 931 // D++ n -> p p K0 K0b (3)
952 // D++ n -> p n K+ K0b (3) 932 // D++ n -> p n K+ K0b (3)
953 // 933 //
954 // D+ p -> p p K+ K- (3) 934 // D+ p -> p p K+ K- (3)
955 // D+ p -> p p K0 K0b (1) 935 // D+ p -> p p K0 K0b (1)
956 // D+ p -> p n K+ K0b (3) 936 // D+ p -> p n K+ K0b (3)
957 // 937 //
958 // D+ n -> p p K0 K- (2) 938 // D+ n -> p p K0 K- (2)
959 // D+ n -> p n K+ K- (1) 939 // D+ n -> p n K+ K- (1)
960 // D+ n -> p n K0 K0b (3) 940 // D+ n -> p n K0 K0b (3)
961 // D+ n -> n n K+ K0b (2) 941 // D+ n -> n n K+ K0b (2)
962 // 942 //
963 943
964 // assert((p1->isNucleon() && p2->isResonance( 944 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
965 945
966 G4double sigma = 0.; 946 G4double sigma = 0.;
967 const G4int iso = ParticleTable::getIs 947 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
968 const G4double ener = 0.001*Kinematics 948 const G4double ener = 0.001*KinematicsUtils::totalEnergyInCM(p1, p2); // GeV
969 949
970 if(ener <= 2.872) 950 if(ener <= 2.872)
971 return 0.; 951 return 0.;
972 952
973 if(iso == 0)// D+ n 953 if(iso == 0)// D+ n
974 sigma = 8* 22./60. * 3. *std::pow( << 954 sigma = 8* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
975 else if (ParticleTable::getIsospin(p1- 955 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
976 sigma = 7* 22./60. * 3. *std::pow( << 956 sigma = 7* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
977 else if (std::abs(iso) == 2)// D++ n 957 else if (std::abs(iso) == 2)// D++ n
978 sigma = 9* 22./60. * 3. *std::pow( << 958 sigma = 9* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
979 else // D++ p 959 else // D++ p
980 sigma = 6* 22./60. * 3. *std::pow( << 960 sigma = 6* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
981 961
982 return sigma; 962 return sigma;
983 } 963 }
984 964
985 /// \brief piN to strange cross sections 965 /// \brief piN to strange cross sections
986 966
987 G4double CrossSectionsStrangeness::NpiToLK 967 G4double CrossSectionsStrangeness::NpiToLK(Particle const * const p1, Particle const * const p2) {
988 // 968 //
989 // Pion-Nucleon producing Lambda- 969 // Pion-Nucleon producing Lambda-Kaon cross sections
990 // 970 //
991 // ratio 971 // ratio
992 // p pi0 -> L K+ (1/2) 972 // p pi0 -> L K+ (1/2)
993 // p pi- -> L K0 (1) 973 // p pi- -> L K0 (1)
994 974
995 // assert((p1->isPion() && p2->isNucleon()) || 975 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
996 976
997 const Particle *pion; 977 const Particle *pion;
998 const Particle *nucleon; 978 const Particle *nucleon;
999 const G4int iso=ParticleTable::getIsos 979 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1000 if(iso == 3 || iso == -3) 980 if(iso == 3 || iso == -3)
1001 return 0.; 981 return 0.;
1002 982
1003 if(p1->isPion()){ 983 if(p1->isPion()){
1004 pion = p1; 984 pion = p1;
1005 nucleon = p2; 985 nucleon = p2;
1006 } 986 }
1007 else{ 987 else{
1008 nucleon = p1; 988 nucleon = p1;
1009 pion = p2; 989 pion = p2;
1010 } 990 }
1011 G4double sigma = 0.; 991 G4double sigma = 0.;
1012 992
1013 if(pion->getType() == PiZero) 993 if(pion->getType() == PiZero)
1014 sigma = 0.5 * p_pimToLK0(pion,nuc 994 sigma = 0.5 * p_pimToLK0(pion,nucleon);
1015 else 995 else
1016 sigma = p_pimToLK0(pion,nucleon); 996 sigma = p_pimToLK0(pion,nucleon);
1017 return sigma; 997 return sigma;
1018 } 998 }
1019 999
1020 G4double CrossSectionsStrangeness::p_pimT 1000 G4double CrossSectionsStrangeness::p_pimToLK0(Particle const * const p1, Particle const * const p2) {
1021 1001
1022 G4double sigma = 0.; 1002 G4double sigma = 0.;
1023 const G4double pLab = 0.001 * Kinemat 1003 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1024 1004
1025 if(pLab < 0.911) 1005 if(pLab < 0.911)
1026 return 0.; 1006 return 0.;
1027 1007
1028 sigma = 0.3936*std::pow(pLab,-1.357)- 1008 sigma = 0.3936*std::pow(pLab,-1.357)-6.052*std::exp(-std::pow(pLab-0.7154,2)/0.02026)-0.16*std::exp(-std::pow(pLab-0.9684,2)/0.001432)+0.489*std::exp(-std::pow(pLab-0.8886,2)/0.08378);
1029 if(sigma < 0.) return 0; 1009 if(sigma < 0.) return 0;
1030 return sigma; 1010 return sigma;
1031 } 1011 }
1032 1012
1033 G4double CrossSectionsStrangeness::NpiToS 1013 G4double CrossSectionsStrangeness::NpiToSK(Particle const * const p1, Particle const * const p2) {
1034 // 1014 //
1035 // Pion-Nucleon producing Sigma- 1015 // Pion-Nucleon producing Sigma-Kaon cross sections
1036 // 1016 //
1037 // ratio 1017 // ratio
1038 // p pi+ (5/3) p pi0 (11/6) p p 1018 // p pi+ (5/3) p pi0 (11/6) p pi- (2)
1039 // 1019 //
1040 // p pi+ -> S+ K+ (10) 1020 // p pi+ -> S+ K+ (10)
1041 // p pi0 -> S+ K0 (6)* 1021 // p pi0 -> S+ K0 (6)*
1042 // p pi0 -> S0 K+ (5) 1022 // p pi0 -> S0 K+ (5)
1043 // p pi- -> S0 K0 (6)* 1023 // p pi- -> S0 K0 (6)*
1044 // p pi- -> S- K+ (6) 1024 // p pi- -> S- K+ (6)
1045 1025
1046 // assert((p1->isPion() && p2->isNucleon()) | 1026 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1047 1027
1048 const Particle *pion; 1028 const Particle *pion;
1049 const Particle *nucleon; 1029 const Particle *nucleon;
1050 const G4int iso=ParticleTable::getIso 1030 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1051 1031
1052 if(p1->isPion()){ 1032 if(p1->isPion()){
1053 pion = p1; 1033 pion = p1;
1054 nucleon = p2; 1034 nucleon = p2;
1055 } 1035 }
1056 else{ 1036 else{
1057 nucleon = p1; 1037 nucleon = p1;
1058 pion = p2; 1038 pion = p2;
1059 } 1039 }
1060 G4double sigma = 0.; 1040 G4double sigma = 0.;
1061 1041
1062 if(iso == 3 || iso == -3) 1042 if(iso == 3 || iso == -3)
1063 sigma = p_pipToSpKp(pion,nucleon) 1043 sigma = p_pipToSpKp(pion,nucleon);
1064 else if(pion->getType() == PiZero) 1044 else if(pion->getType() == PiZero)
1065 sigma = p_pizToSzKp(pion,nucleon) 1045 sigma = p_pizToSzKp(pion,nucleon)+p_pimToSzKz(pion,nucleon);
1066 else if(iso == 1 || iso == -1) 1046 else if(iso == 1 || iso == -1)
1067 sigma = p_pimToSzKz(pion,nucleon) 1047 sigma = p_pimToSzKz(pion,nucleon)+p_pimToSmKp(pion,nucleon);
1068 else // should never append 1048 else // should never append
1069 sigma = 0.; 1049 sigma = 0.;
1070 1050
1071 return sigma; 1051 return sigma;
1072 } 1052 }
1073 G4double CrossSectionsStrangeness::p_pimT 1053 G4double CrossSectionsStrangeness::p_pimToSmKp(Particle const * const p1, Particle const * const p2) {
1074 1054
1075 G4double sigma = 0.; 1055 G4double sigma = 0.;
1076 const G4double pLab = 0.001 * Kinemat 1056 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1077 1057
1078 if(pLab < 1.0356) 1058 if(pLab < 1.0356)
1079 return 0.; 1059 return 0.;
1080 1060
1081 sigma = 4.352*std::pow(pLab-1.0356,1. 1061 sigma = 4.352*std::pow(pLab-1.0356,1.006)/(std::pow(pLab+1.0356,0.0978)*std::pow(pLab,5.375));
1082 1062
1083 if(sigma < 0.) // should never append 1063 if(sigma < 0.) // should never append
1084 return 0; 1064 return 0;
1085 1065
1086 return sigma; 1066 return sigma;
1087 } 1067 }
1088 G4double CrossSectionsStrangeness::p_pipT 1068 G4double CrossSectionsStrangeness::p_pipToSpKp(Particle const * const p1, Particle const * const p2) {
1089 1069
1090 G4double sigma = 0.; 1070 G4double sigma = 0.;
1091 const G4double pLab = 0.001 * Kinemat 1071 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1092 1072
1093 if(pLab < 1.0428) 1073 if(pLab < 1.0428)
1094 return 0.; 1074 return 0.;
1095 1075
1096 sigma = 0.001897*std::pow(pLab-1.0428 1076 sigma = 0.001897*std::pow(pLab-1.0428,2.869)/(std::pow(pLab+1.0428,-16.68)*std::pow(pLab,19.1));
1097 1077
1098 if(sigma < 0.) // should never append 1078 if(sigma < 0.) // should never append
1099 return 0; 1079 return 0;
1100 1080
1101 return sigma; 1081 return sigma;
1102 } 1082 }
1103 G4double CrossSectionsStrangeness::p_pizT 1083 G4double CrossSectionsStrangeness::p_pizToSzKp(Particle const * const p1, Particle const * const p2) {
1104 1084
1105 G4double sigma = 0.; 1085 G4double sigma = 0.;
1106 const G4double pLab = 0.001 * Kinemat 1086 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1107 1087
1108 if(pLab < 1.0356) 1088 if(pLab < 1.0356)
1109 return 0.; 1089 return 0.;
1110 1090
1111 sigma = 3.624*std::pow(pLab-1.0356,1. 1091 sigma = 3.624*std::pow(pLab-1.0356,1.4)/std::pow(pLab,5.14);
1112 1092
1113 if(sigma < 0.) // should never append 1093 if(sigma < 0.) // should never append
1114 return 0; 1094 return 0;
1115 1095
1116 return sigma; 1096 return sigma;
1117 } 1097 }
1118 G4double CrossSectionsStrangeness::p_pimT 1098 G4double CrossSectionsStrangeness::p_pimToSzKz(Particle const * const p1, Particle const * const p2) {
1119 1099
1120 G4double sigma = 0.; 1100 G4double sigma = 0.;
1121 const G4double pLab = 0.001 * Kinemat 1101 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1122 1102
1123 if((p1->getType() == PiZero && pLab < 1103 if((p1->getType() == PiZero && pLab < 1.0356) || (pLab < 1.034))
1124 return 0.; 1104 return 0.;
1125 1105
1126 sigma = 0.3474*std::pow(pLab-1.034,0. 1106 sigma = 0.3474*std::pow(pLab-1.034,0.07678)/std::pow(pLab,1.627);
1127 1107
1128 if(sigma < 0.) // should never append 1108 if(sigma < 0.) // should never append
1129 return 0; 1109 return 0;
1130 1110
1131 return sigma; 1111 return sigma;
1132 } 1112 }
1133 1113
1134 G4double CrossSectionsStrangeness::NpiToL 1114 G4double CrossSectionsStrangeness::NpiToLKpi(Particle const * const p1, Particle const * const p2) {
1135 // 1115 //
1136 // Pion-Nucleon producing Lambda 1116 // Pion-Nucleon producing Lambda-Kaon-pion cross sections
1137 // 1117 //
1138 // ratio 1118 // ratio
1139 // p pi+ (1) p pi0 (3/2) p 1119 // p pi+ (1) p pi0 (3/2) p pi- (2)
1140 // 1120 //
1141 // p pi0 -> L K+ pi0 (1/2) 1121 // p pi0 -> L K+ pi0 (1/2)
1142 // all the others (1) 1122 // all the others (1)
1143 // 1123 //
1144 // assert((p1->isPion() && p2->isNucleon()) | 1124 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1145 1125
1146 G4double sigma=0.; 1126 G4double sigma=0.;
1147 const Particle *pion; 1127 const Particle *pion;
1148 const Particle *nucleon; 1128 const Particle *nucleon;
1149 1129
1150 if(p1->isPion()){ 1130 if(p1->isPion()){
1151 pion = p1; 1131 pion = p1;
1152 nucleon = p2; 1132 nucleon = p2;
1153 } 1133 }
1154 else{ 1134 else{
1155 nucleon = p1; 1135 nucleon = p1;
1156 pion = p2; 1136 pion = p2;
1157 } 1137 }
1158 const G4int iso=ParticleTable::getIso 1138 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1159 const G4double pLab = 0.001*Kinematic 1139 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1160 1140
1161 if(pLab < 1.147) 1141 if(pLab < 1.147)
1162 return 0.; 1142 return 0.;
1163 1143
1164 if(iso == 3 || iso == -3) 1144 if(iso == 3 || iso == -3)
1165 sigma = 146.2*std::pow(pLab-1.147 1145 sigma = 146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1166 else if(pion->getType() == PiZero) 1146 else if(pion->getType() == PiZero)
1167 sigma = 1.5*146.2*std::pow(pLab-1 1147 sigma = 1.5*146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1168 else 1148 else
1169 sigma = 2*146.2*std::pow(pLab-1.1 1149 sigma = 2*146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1170 1150
1171 return sigma; 1151 return sigma;
1172 } 1152 }
1173 1153
1174 G4double CrossSectionsStrangeness::NpiToS 1154 G4double CrossSectionsStrangeness::NpiToSKpi(Particle const * const p1, Particle const * const p2) {
1175 // 1155 //
1176 // Pion-Nucleon producing Sigma- 1156 // Pion-Nucleon producing Sigma-Kaon-pion cross sections
1177 // 1157 //
1178 //ratio 1158 //ratio
1179 // p pi+ (2.25) p pi0 (2.625) p 1159 // p pi+ (2.25) p pi0 (2.625) p pi-(3)
1180 // 1160 //
1181 // p pi+ -> S+ pi+ K0 (5/4) 1161 // p pi+ -> S+ pi+ K0 (5/4)
1182 // p pi+ -> S+ pi0 K+ (3/4) 1162 // p pi+ -> S+ pi0 K+ (3/4)
1183 // p pi+ -> S0 pi+ K+ (1/4) 1163 // p pi+ -> S0 pi+ K+ (1/4)
1184 // p pi0 -> S+ pi0 K0 (1/2) 1164 // p pi0 -> S+ pi0 K0 (1/2)
1185 // p pi0 -> S+ pi- K+ (1/2) 1165 // p pi0 -> S+ pi- K+ (1/2)
1186 // p pi0 -> S0 pi+ K0 (3/4) 1166 // p pi0 -> S0 pi+ K0 (3/4)
1187 // p pi0 -> S0 pi0 K+ (3/8) 1167 // p pi0 -> S0 pi0 K+ (3/8)
1188 // p pi0 -> S- pi+ K+ (1/2) 1168 // p pi0 -> S- pi+ K+ (1/2)
1189 // p pi- -> S+ pi- K0 (3/8) 1169 // p pi- -> S+ pi- K0 (3/8)
1190 // p pi- -> S0 pi0 K0 (5/8) 1170 // p pi- -> S0 pi0 K0 (5/8)
1191 // p pi- -> S0 pi- K+ (5/8) 1171 // p pi- -> S0 pi- K+ (5/8)
1192 // p pi- -> S- pi+ K0 (1) 1172 // p pi- -> S- pi+ K0 (1)
1193 // p pi- -> S- pi0 K+ (3/8) 1173 // p pi- -> S- pi0 K+ (3/8)
1194 1174
1195 // assert((p1->isPion() && p2->isNucleon()) | 1175 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1196 1176
1197 G4double sigma=0.; 1177 G4double sigma=0.;
1198 const Particle *pion; 1178 const Particle *pion;
1199 const Particle *nucleon; 1179 const Particle *nucleon;
1200 1180
1201 if(p1->isPion()){ 1181 if(p1->isPion()){
1202 pion = p1; 1182 pion = p1;
1203 nucleon = p2; 1183 nucleon = p2;
1204 } 1184 }
1205 else{ 1185 else{
1206 nucleon = p1; 1186 nucleon = p1;
1207 pion = p2; 1187 pion = p2;
1208 } 1188 }
1209 const G4int iso=ParticleTable::getIso 1189 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1210 const G4double pLab = 0.001*Kinematic 1190 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1211 1191
1212 if(pLab <= 1.3041) 1192 if(pLab <= 1.3041)
1213 return 0.; 1193 return 0.;
1214 1194
1215 if(iso == 3 || iso == -3) 1195 if(iso == 3 || iso == -3)
1216 sigma = 2.25*8.139*std::pow(pLab- 1196 sigma = 2.25*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1217 else if(pion->getType() == PiZero) 1197 else if(pion->getType() == PiZero)
1218 sigma = 2.625*8.139*std::pow(pLab 1198 sigma = 2.625*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1219 else 1199 else
1220 sigma = 3.*8.139*std::pow(pLab-1. 1200 sigma = 3.*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1221 1201
1222 return sigma; 1202 return sigma;
1223 } 1203 }
1224 1204
1225 G4double CrossSectionsStrangeness::NpiToL 1205 G4double CrossSectionsStrangeness::NpiToLK2pi(Particle const * const p1, Particle const * const p2) {
1226 // 1206 //
1227 // Pion-Nucleon producing Lambda 1207 // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1228 // 1208 //
1229 // p pi+ (2) p pi0 (1.75) p pi- 1209 // p pi+ (2) p pi0 (1.75) p pi- (2.5)
1230 // 1210 //
1231 // p pi+ -> L K+ pi+ pi0 (1) 1211 // p pi+ -> L K+ pi+ pi0 (1)
1232 // p pi+ -> L K0 pi+ pi+ (1) 1212 // p pi+ -> L K0 pi+ pi+ (1)
1233 // p pi0 -> L K+ pi0 pi0 (1/4) 1213 // p pi0 -> L K+ pi0 pi0 (1/4)
1234 // p pi0 -> L K+ pi+ pi- (1) 1214 // p pi0 -> L K+ pi+ pi- (1)
1235 // p pi0 -> L K0 pi+ pi0 (1/2) 1215 // p pi0 -> L K0 pi+ pi0 (1/2)
1236 // p pi- -> L K+ pi0 pi- (1) 1216 // p pi- -> L K+ pi0 pi- (1)
1237 // p pi- -> L K0 pi+ pi- (1) 1217 // p pi- -> L K0 pi+ pi- (1)
1238 // p pi- -> L K0 pi0 pi0 (1/2) 1218 // p pi- -> L K0 pi0 pi0 (1/2)
1239 1219
1240 // assert((p1->isPion() && p2->isNucleon()) | 1220 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1241 1221
1242 G4double sigma=0.; 1222 G4double sigma=0.;
1243 const Particle *pion; 1223 const Particle *pion;
1244 const Particle *nucleon; 1224 const Particle *nucleon;
1245 1225
1246 if(p1->isPion()){ 1226 if(p1->isPion()){
1247 pion = p1; 1227 pion = p1;
1248 nucleon = p2; 1228 nucleon = p2;
1249 } 1229 }
1250 else{ 1230 else{
1251 nucleon = p1; 1231 nucleon = p1;
1252 pion = p2; 1232 pion = p2;
1253 } 1233 }
1254 const G4int iso=ParticleTable::getIso 1234 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1255 const G4double pLab = 0.001*Kinematic 1235 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1256 1236
1257 if(pLab <= 1.4162) 1237 if(pLab <= 1.4162)
1258 return 0.; 1238 return 0.;
1259 1239
1260 if(iso == 3 || iso == -3) 1240 if(iso == 3 || iso == -3)
1261 sigma = 2*18.77*std::pow(pLab-1.4 1241 sigma = 2*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1262 else if(pion->getType() == PiZero) 1242 else if(pion->getType() == PiZero)
1263 sigma = 1.75*18.77*std::pow(pLab- 1243 sigma = 1.75*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1264 else 1244 else
1265 sigma = 2.5*18.77*std::pow(pLab-1 1245 sigma = 2.5*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1266 1246
1267 return sigma; 1247 return sigma;
1268 } 1248 }
1269 1249
1270 G4double CrossSectionsStrangeness::NpiToS 1250 G4double CrossSectionsStrangeness::NpiToSK2pi(Particle const * const p1, Particle const * const p2) {
1271 // 1251 //
1272 // Pion-Nucleon producing Lambda 1252 // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1273 // 1253 //
1274 // ratio 1254 // ratio
1275 // p pi+ (3.25) p pi0 (3.5) p p 1255 // p pi+ (3.25) p pi0 (3.5) p pi- (3.75)
1276 // 1256 //
1277 // p pi+ -> S+ K+ pi+ pi- (1) 1257 // p pi+ -> S+ K+ pi+ pi- (1)
1278 // p pi+ -> S+ K+ pi0 pi0 (1/4) 1258 // p pi+ -> S+ K+ pi0 pi0 (1/4)
1279 // p pi+ -> S0 K+ pi+ pi0 (1/2) 1259 // p pi+ -> S0 K+ pi+ pi0 (1/2)
1280 // p pi+ -> S- K+ pi+ pi+ (1/4) 1260 // p pi+ -> S- K+ pi+ pi+ (1/4)
1281 // p pi+ -> S+ K0 pi+ pi0 (1) 1261 // p pi+ -> S+ K0 pi+ pi0 (1)
1282 // p pi+ -> S0 K0 pi+ pi+ (1/4) 1262 // p pi+ -> S0 K0 pi+ pi+ (1/4)
1283 // 1263 //
1284 // p pi0 -> S+ K+ pi0 pi- (1/2) 1264 // p pi0 -> S+ K+ pi0 pi- (1/2)
1285 // p pi0 -> S0 K+ pi+ pi- (1/2) 1265 // p pi0 -> S0 K+ pi+ pi- (1/2)
1286 // p pi0 -> S0 K+ pi0 pi0 (1/4) 1266 // p pi0 -> S0 K+ pi0 pi0 (1/4)
1287 // p pi0 -> S- K+ pi+ pi0 (1/4) 1267 // p pi0 -> S- K+ pi+ pi0 (1/4)
1288 // p pi0 -> S+ K0 pi+ pi- (1) 1268 // p pi0 -> S+ K0 pi+ pi- (1)
1289 // p pi0 -> S+ K0 pi0 pi0 (1/4) 1269 // p pi0 -> S+ K0 pi0 pi0 (1/4)
1290 // p pi0 -> S0 K0 pi+ pi0 (1/4) 1270 // p pi0 -> S0 K0 pi+ pi0 (1/4)
1291 // p pi0 -> S- K0 pi+ pi+ (1/2) 1271 // p pi0 -> S- K0 pi+ pi+ (1/2)
1292 // 1272 //
1293 // p pi- -> S+ K+ pi- pi- (1/4) 1273 // p pi- -> S+ K+ pi- pi- (1/4)
1294 // p pi- -> S0 K+ pi0 pi- (1/2) 1274 // p pi- -> S0 K+ pi0 pi- (1/2)
1295 // p pi- -> S- K+ pi+ pi- (1/4) 1275 // p pi- -> S- K+ pi+ pi- (1/4)
1296 // p pi- -> S- K+ pi0 pi0 (1/4) 1276 // p pi- -> S- K+ pi0 pi0 (1/4)
1297 // p pi- -> S+ K0 pi0 pi- (1/2) 1277 // p pi- -> S+ K0 pi0 pi- (1/2)
1298 // p pi- -> S0 K0 pi+ pi- (1) 1278 // p pi- -> S0 K0 pi+ pi- (1)
1299 // p pi- -> S0 K0 pi0 pi0 (1/2) 1279 // p pi- -> S0 K0 pi0 pi0 (1/2)
1300 // p pi- -> S- K0 pi+ pi0 (1/2) 1280 // p pi- -> S- K0 pi+ pi0 (1/2)
1301 1281
1302 // assert((p1->isPion() && p2->isNucleon()) | 1282 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1303 1283
1304 G4double sigma=0.; 1284 G4double sigma=0.;
1305 const Particle *pion; 1285 const Particle *pion;
1306 const Particle *nucleon; 1286 const Particle *nucleon;
1307 1287
1308 if(p1->isPion()){ 1288 if(p1->isPion()){
1309 pion = p1; 1289 pion = p1;
1310 nucleon = p2; 1290 nucleon = p2;
1311 } 1291 }
1312 else{ 1292 else{
1313 nucleon = p1; 1293 nucleon = p1;
1314 pion = p2; 1294 pion = p2;
1315 } 1295 }
1316 const G4int iso=ParticleTable::getIso 1296 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1317 const G4double pLab = 0.001*Kinematic 1297 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1318 1298
1319 if(pLab <= 1.5851) 1299 if(pLab <= 1.5851)
1320 return 0.; 1300 return 0.;
1321 1301
1322 if(iso == 3 || iso == -3) 1302 if(iso == 3 || iso == -3)
1323 sigma = 3.25*137.6*std::pow(pLab- 1303 sigma = 3.25*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1324 else if(pion->getType() == PiZero) 1304 else if(pion->getType() == PiZero)
1325 sigma = 3.5*137.6*std::pow(pLab-1 1305 sigma = 3.5*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1326 else 1306 else
1327 sigma = 3.75*137.6*std::pow(pLab- 1307 sigma = 3.75*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1328 1308
1329 return sigma; 1309 return sigma;
1330 } 1310 }
1331 1311
1332 G4double CrossSectionsStrangeness::NpiToN 1312 G4double CrossSectionsStrangeness::NpiToNKKb(Particle const * const p1, Particle const * const p2) {
1333 // 1313 //
1334 // Pion-Nucleon producing Nucleo 1314 // Pion-Nucleon producing Nucleon-Kaon-antiKaon cross sections
1335 // 1315 //
1336 // ratio 1316 // ratio
1337 // p pi+ (1/2) p pi0 (3/2) p pi 1317 // p pi+ (1/2) p pi0 (3/2) p pi- (5/2)
1338 // 1318 //
1339 // p pi+ -> p K+ K0b (1/2) 1319 // p pi+ -> p K+ K0b (1/2)
1340 // p pi0 -> p K0 K0b (1/4) 1320 // p pi0 -> p K0 K0b (1/4)
1341 // p pi0 -> p K+ K- (1/4) 1321 // p pi0 -> p K+ K- (1/4)
1342 // p pi0 -> n K+ K0b (1) 1322 // p pi0 -> n K+ K0b (1)
1343 // p pi- -> p K0 K- (1/2) 1323 // p pi- -> p K0 K- (1/2)
1344 // p pi- -> n K+ K- (1) 1324 // p pi- -> n K+ K- (1)
1345 // p pi- -> n K0 K0b (1) 1325 // p pi- -> n K0 K0b (1)
1346 1326
1347 // assert((p1->isPion() && p2->isNucleon()) | 1327 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1348 1328
1349 const Particle *particle1; 1329 const Particle *particle1;
1350 const Particle *particle2; 1330 const Particle *particle2;
1351 1331
1352 if(p1->isPion()){ 1332 if(p1->isPion()){
1353 particle1 = p1; 1333 particle1 = p1;
1354 particle2 = p2; 1334 particle2 = p2;
1355 } 1335 }
1356 else{ 1336 else{
1357 particle1 = p2; 1337 particle1 = p2;
1358 particle2 = p1; 1338 particle2 = p1;
1359 } 1339 }
1360 1340
1361 G4double sigma = 0.; 1341 G4double sigma = 0.;
1362 1342
1363 const G4double pLab = 0.001 * Kinemat 1343 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1364 1344
1365 if(particle1->getType() == PiZero){ 1345 if(particle1->getType() == PiZero){
1366 if(pLab < 1.5066) return 0.; 1346 if(pLab < 1.5066) return 0.;
1367 else if(pLab < 30.) sigma = 3./2. 1347 else if(pLab < 30.) sigma = 3./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1368 else return 0.; 1348 else return 0.;
1369 } 1349 }
1370 else if((particle1->getType() == PiPl 1350 else if((particle1->getType() == PiPlus && particle2->getType() == Neutron) || (particle1->getType() == PiMinus && particle2->getType() == Proton)){
1371 if(pLab < 1.5066) return 0.; 1351 if(pLab < 1.5066) return 0.;
1372 else if(pLab < 30.) sigma = 5./2. 1352 else if(pLab < 30.) sigma = 5./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1373 else return 0.; 1353 else return 0.;
1374 } 1354 }
1375 else{ 1355 else{
1376 if(pLab < 1.5066) return 0.; 1356 if(pLab < 1.5066) return 0.;
1377 else if(pLab < 30.) sigma = 1./2. 1357 else if(pLab < 30.) sigma = 1./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1378 else return 0.; 1358 else return 0.;
1379 } 1359 }
1380 return sigma; 1360 return sigma;
1381 } 1361 }
1382 1362
1383 G4double CrossSectionsStrangeness::NpiToM 1363 G4double CrossSectionsStrangeness::NpiToMissingStrangeness(Particle const * const p1, Particle const * const p2) {
1384 // 1364 //
1385 // Pion-Nucleon missing strangen 1365 // Pion-Nucleon missing strangeness production cross sections
1386 // 1366 //
1387 // assert((p1->isPion() && p2->isNucleon()) | 1367 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1388 1368
1389 const Particle *pion; 1369 const Particle *pion;
1390 const Particle *nucleon; 1370 const Particle *nucleon;
1391 1371
1392 if(p1->isPion()){ 1372 if(p1->isPion()){
1393 pion = p1; 1373 pion = p1;
1394 nucleon = p2; 1374 nucleon = p2;
1395 } 1375 }
1396 else{ 1376 else{
1397 pion = p2; 1377 pion = p2;
1398 nucleon = p1; 1378 nucleon = p1;
1399 } 1379 }
1400 1380
1401 G4double sigma = 0.; 1381 G4double sigma = 0.;
1402 1382
1403 const G4double pLab = 0.001 * Kinemat 1383 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(pion,nucleon); // GeV
1404 if(pLab < 2.2) return 0.; 1384 if(pLab < 2.2) return 0.;
1405 1385
1406 if(pion->getType() == PiZero){ 1386 if(pion->getType() == PiZero){
1407 if(pLab < 30.) sigma = 4.4755*std 1387 if(pLab < 30.) sigma = 4.4755*std::pow((pLab - 2.2),1.927)/std::pow(pLab,1.89343);
1408 else return 0.; 1388 else return 0.;
1409 } 1389 }
1410 else if((pion->getType() == PiPlus && 1390 else if((pion->getType() == PiPlus && nucleon->getType() == Neutron) || (pion->getType() == PiMinus && nucleon->getType() == Proton)){
1411 if(pLab < 30.) sigma = 5.1*std::p 1391 if(pLab < 30.) sigma = 5.1*std::pow((pLab - 2.2),1.854)/std::pow(pLab,1.904);
1412 else return 0.; 1392 else return 0.;
1413 } 1393 }
1414 else{ 1394 else{
1415 if(pLab < 30.) sigma = 3.851*std: 1395 if(pLab < 30.) sigma = 3.851*std::pow((pLab - 2.2),2)/std::pow(pLab,1.88286);
1416 else return 0.; 1396 else return 0.;
1417 } 1397 }
1418 return sigma; 1398 return sigma;
1419 } 1399 }
1420 1400
1421 /// \brief NY cross sections 1401 /// \brief NY cross sections
1422 1402
1423 G4double CrossSectionsStrangeness::NLToNS 1403 G4double CrossSectionsStrangeness::NLToNS(Particle const * const p1, Particle const * const p2) {
1424 // 1404 //
1425 // Nucleon-Lambda producing Nucl 1405 // Nucleon-Lambda producing Nucleon-Sigma cross sections
1426 // 1406 //
1427 // ratio 1407 // ratio
1428 // p L -> p S0 (1/2) 1408 // p L -> p S0 (1/2)
1429 // p L -> n S+ (1) 1409 // p L -> n S+ (1)
1430 1410
1431 1411
1432 // assert((p1->isLambda() && p2->isNucleon()) 1412 // assert((p1->isLambda() && p2->isNucleon()) || (p2->isLambda() && p1->isNucleon()));
1433 1413
1434 G4double sigma = 0.; 1414 G4double sigma = 0.;
1435 1415
1436 const Particle *particle1; 1416 const Particle *particle1;
1437 const Particle *particle2; 1417 const Particle *particle2;
1438 1418
1439 if(p1->isLambda()){ 1419 if(p1->isLambda()){
1440 particle1 = p1; 1420 particle1 = p1;
1441 particle2 = p2; 1421 particle2 = p2;
1442 } 1422 }
1443 else{ 1423 else{
1444 particle1 = p2; 1424 particle1 = p2;
1445 particle2 = p1; 1425 particle2 = p1;
1446 } 1426 }
1447 1427
1448 const G4double pLab = 0.001 * Kinemat 1428 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2);
1449 1429
1450 if(pLab < 0.664) 1430 if(pLab < 0.664)
1451 return 0.; 1431 return 0.;
1452 1432
1453 sigma = 3 * 8.74*std::pow((pLab-0.664 1433 sigma = 3 * 8.74*std::pow((pLab-0.664),0.438)/std::pow(pLab,2.717); // 3 * L p -> S0 p
1454 1434
1455 return sigma; 1435 return sigma;
1456 } 1436 }
1457 1437
1458 G4double CrossSectionsStrangeness::NSToNL 1438 G4double CrossSectionsStrangeness::NSToNL(Particle const * const p1, Particle const * const p2) {
1459 // 1439 //
1460 // Nucleon-Lambda producing Nucl 1440 // Nucleon-Lambda producing Nucleon-Sigma cross sections
1461 // 1441 //
1462 // ratio 1442 // ratio
1463 // p S0 -> p L (1/2) 1443 // p S0 -> p L (1/2)
1464 // p S- -> n L (1) 1444 // p S- -> n L (1)
1465 1445
1466 // assert((p1->isSigma() && p2->isNucleon()) 1446 // assert((p1->isSigma() && p2->isNucleon()) || (p2->isSigma() && p1->isNucleon()));
1467 1447
1468 const G4int iso=ParticleTable::getIso 1448 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1469 if(iso == 3 || iso == -3) 1449 if(iso == 3 || iso == -3)
1470 return 0.; 1450 return 0.;
1471 1451
1472 G4double sigma; 1452 G4double sigma;
1473 const Particle *particle1; 1453 const Particle *particle1;
1474 const Particle *particle2; 1454 const Particle *particle2;
1475 1455
1476 if(p1->isSigma()){ 1456 if(p1->isSigma()){
1477 particle1 = p1; 1457 particle1 = p1;
1478 particle2 = p2; 1458 particle2 = p2;
1479 } 1459 }
1480 else{ 1460 else{
1481 particle1 = p2; 1461 particle1 = p2;
1482 particle2 = p1; 1462 particle2 = p1;
1483 } 1463 }
1484 const G4double pLab = 0.001 * Kinemat 1464 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1485 1465
1486 if(particle1->getType() == SigmaZero) 1466 if(particle1->getType() == SigmaZero){
1487 if(pLab < 0.1) return 100.; // cu 1467 if(pLab < 0.1) return 100.; // cut-max
1488 sigma = 8.23*std::pow(pLab,-1.087 1468 sigma = 8.23*std::pow(pLab,-1.087);
1489 } 1469 }
1490 else{ 1470 else{
1491 if(pLab < 0.1) return 200.; // cu 1471 if(pLab < 0.1) return 200.; // cut-max
1492 sigma = 16.46*std::pow(pLab,-1.08 1472 sigma = 16.46*std::pow(pLab,-1.087);
1493 } 1473 }
1494 return sigma; 1474 return sigma;
1495 } 1475 }
1496 1476
1497 G4double CrossSectionsStrangeness::NSToNS 1477 G4double CrossSectionsStrangeness::NSToNS(Particle const * const p1, Particle const * const p2) {
1498 1478
1499 // assert((p1->isSigma() && p2->isNucleon()) 1479 // assert((p1->isSigma() && p2->isNucleon()) || (p2->isSigma() && p1->isNucleon()));
1500 1480
1501 const G4int iso=ParticleTable::getIso 1481 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1502 if(iso == 3 || iso == -3) 1482 if(iso == 3 || iso == -3)
1503 return 0.; // only quasi-elastic 1483 return 0.; // only quasi-elastic here
1504 1484
1505 const Particle *particle1; 1485 const Particle *particle1;
1506 const Particle *particle2; 1486 const Particle *particle2;
1507 1487
1508 if(p1->isSigma()){ 1488 if(p1->isSigma()){
1509 particle1 = p1; 1489 particle1 = p1;
1510 particle2 = p2; 1490 particle2 = p2;
1511 } 1491 }
1512 else{ 1492 else{
1513 particle1 = p2; 1493 particle1 = p2;
1514 particle2 = p1; 1494 particle2 = p1;
1515 } 1495 }
1516 const G4double pLab = 0.001 * Kinemat 1496 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1517 1497
1518 if(particle2->getType() == Neutron && 1498 if(particle2->getType() == Neutron && pLab < 0.162) return 0.;
1519 else if(pLab < 0.1035) return 200.; / 1499 else if(pLab < 0.1035) return 200.; // cut-max
1520 1500
1521 return 13.79*std::pow(pLab,-1.181); 1501 return 13.79*std::pow(pLab,-1.181);
1522 } 1502 }
1523 1503
1524 /// \brief NK cross sections 1504 /// \brief NK cross sections
1525 1505
1526 G4double CrossSectionsStrangeness::NKToNK 1506 G4double CrossSectionsStrangeness::NKToNK(Particle const * const p1, Particle const * const p2) {
1527 // 1507 //
1528 // Nucleon-Kaon quasi-elastic cr 1508 // Nucleon-Kaon quasi-elastic cross sections
1529 // 1509 //
1530 // assert((p1->isNucleon() && p2->isKaon()) | 1510 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1531 1511
1532 const G4int iso=ParticleTable::getIso 1512 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1533 1513
1534 if(iso != 0) 1514 if(iso != 0)
1535 return 0.; 1515 return 0.;
1536 1516
1537 const Particle *particle1; 1517 const Particle *particle1;
1538 const Particle *particle2; 1518 const Particle *particle2;
1539 1519
1540 if(p1->isKaon()){ 1520 if(p1->isKaon()){
1541 particle1 = p1; 1521 particle1 = p1;
1542 particle2 = p2; 1522 particle2 = p2;
1543 } 1523 }
1544 else{ 1524 else{
1545 particle1 = p2; 1525 particle1 = p2;
1546 particle2 = p1; 1526 particle2 = p1;
1547 } 1527 }
1548 1528
1549 G4double sigma = 0.; 1529 G4double sigma = 0.;
1550 G4double pLab = 0.001 * KinematicsUti 1530 G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1551 1531
1552 if(particle1->getType() == Proton) 1532 if(particle1->getType() == Proton)
1553 pLab += 2*0.0774; 1533 pLab += 2*0.0774;
1554 1534
1555 if(pLab <= 0.0774) 1535 if(pLab <= 0.0774)
1556 return 0.; 1536 return 0.;
1557 1537
1558 sigma = 12.84*std::pow((pLab-0.0774), 1538 sigma = 12.84*std::pow((pLab-0.0774),18.19)/std::pow((pLab),20.41);
1559 1539
1560 return sigma; 1540 return sigma;
1561 } 1541 }
1562 1542
1563 G4double CrossSectionsStrangeness::NKToNK 1543 G4double CrossSectionsStrangeness::NKToNKpi(Particle const * const p1, Particle const * const p2) {
1564 // 1544 //
1565 // Nucleon-Kaon producing Nucleo 1545 // Nucleon-Kaon producing Nucleon-Kaon-pion cross sections
1566 // 1546 //
1567 // Ratio determined by meson symmetry 1547 // Ratio determined by meson symmetry using only "resonante" diagram (with Delta or K*)
1568 // 1548 //
1569 // ratio: K+ p (5) K0 p (5.545) 1549 // ratio: K+ p (5) K0 p (5.545)
1570 // 1550 //
1571 // K+ p -> p K+ pi0 1.2 1551 // K+ p -> p K+ pi0 1.2
1572 // K+ p -> p K0 pi+ 3 1552 // K+ p -> p K0 pi+ 3
1573 // K+ p -> n K+ pi+ 0.8 1553 // K+ p -> n K+ pi+ 0.8
1574 // K0 p -> p K+ pi- 1 1554 // K0 p -> p K+ pi- 1
1575 // K0 p -> p K0 pi0 0.845 1555 // K0 p -> p K0 pi0 0.845
1576 // K0 p -> n K+ pi0 1.47 1556 // K0 p -> n K+ pi0 1.47
1577 // K0 p -> n K0 pi+ 2.23 1557 // K0 p -> n K0 pi+ 2.23
1578 // assert((p1->isNucleon() && p2->isKaon()) | 1558 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1579 1559
1580 const G4int iso=ParticleTable::getIso 1560 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1581 1561
1582 const Particle *particle1; 1562 const Particle *particle1;
1583 const Particle *particle2; 1563 const Particle *particle2;
1584 1564
1585 if(p1->isKaon()){ 1565 if(p1->isKaon()){
1586 particle1 = p1; 1566 particle1 = p1;
1587 particle2 = p2; 1567 particle2 = p2;
1588 } 1568 }
1589 else{ 1569 else{
1590 particle1 = p2; 1570 particle1 = p2;
1591 particle2 = p1; 1571 particle2 = p1;
1592 } 1572 }
1593 1573
1594 G4double sigma = 0.; 1574 G4double sigma = 0.;
1595 const G4double pLab = 0.001 * Kinemat 1575 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1596 1576
1597 if(pLab <= 0.53) 1577 if(pLab <= 0.53)
1598 return 0.; 1578 return 0.;
1599 1579
1600 if(iso == 0) 1580 if(iso == 0)
1601 sigma = 5.55*116.8*std::pow(pLab- 1581 sigma = 5.55*116.8*std::pow(pLab-0.53,6.874)/std::pow(pLab,10.11);
1602 else 1582 else
1603 sigma = 5.*116.8*std::pow(pLab-0. 1583 sigma = 5.*116.8*std::pow(pLab-0.53,6.874)/std::pow(pLab,10.11);;
1604 return sigma; 1584 return sigma;
1605 } 1585 }
1606 1586
1607 G4double CrossSectionsStrangeness::NKToNK 1587 G4double CrossSectionsStrangeness::NKToNK2pi(Particle const * const p1, Particle const * const p2) {
1608 // 1588 //
1609 // Nucleon-Kaon producing Nucleo 1589 // Nucleon-Kaon producing Nucleon-Kaon-2pion cross sections
1610 // 1590 //
1611 // p K+ (2.875) p K0 (3.125) 1591 // p K+ (2.875) p K0 (3.125)
1612 // 1592 //
1613 // p K+ -> p K+ pi+ pi- (1) 1593 // p K+ -> p K+ pi+ pi- (1)
1614 // p K+ -> p K+ pi0 pi0 (1/8) 1594 // p K+ -> p K+ pi0 pi0 (1/8)
1615 // p K+ -> p K0 pi+ pi0 (1) 1595 // p K+ -> p K0 pi+ pi0 (1)
1616 // p K+ -> n K+ pi+ pi0 (1/2) 1596 // p K+ -> n K+ pi+ pi0 (1/2)
1617 // p K+ -> n K0 pi+ pi+ (1/4) 1597 // p K+ -> n K0 pi+ pi+ (1/4)
1618 // p K0 -> p K+ pi0 pi- (1) 1598 // p K0 -> p K+ pi0 pi- (1)
1619 // p K0 -> p K0 pi+ pi- (1) 1599 // p K0 -> p K0 pi+ pi- (1)
1620 // p K0 -> p K0 pi0 pi0 (1/8) 1600 // p K0 -> p K0 pi0 pi0 (1/8)
1621 // p K0 -> n K+ pi+ pi- (1/4) 1601 // p K0 -> n K+ pi+ pi- (1/4)
1622 // p K0 -> n K+ pi0 pi0 (1/4) 1602 // p K0 -> n K+ pi0 pi0 (1/4)
1623 // p K0 -> n K0 pi+ pi0 (1/2) 1603 // p K0 -> n K0 pi+ pi0 (1/2)
1624 1604
1625 // assert((p1->isNucleon() && p2->isKaon()) | 1605 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1626 1606
1627 const G4int iso=ParticleTable::getIso 1607 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1628 1608
1629 const Particle *particle1; 1609 const Particle *particle1;
1630 const Particle *particle2; 1610 const Particle *particle2;
1631 1611
1632 if(p1->isKaon()){ 1612 if(p1->isKaon()){
1633 particle1 = p1; 1613 particle1 = p1;
1634 particle2 = p2; 1614 particle2 = p2;
1635 } 1615 }
1636 else{ 1616 else{
1637 particle1 = p2; 1617 particle1 = p2;
1638 particle2 = p1; 1618 particle2 = p1;
1639 } 1619 }
1640 1620
1641 G4double sigma = 0.; 1621 G4double sigma = 0.;
1642 const G4double pLab = 0.001 * Kinemat 1622 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1643 1623
1644 if(pLab < 0.812) 1624 if(pLab < 0.812)
1645 sigma = 0.; 1625 sigma = 0.;
1646 else if(pLab < 1.744) 1626 else if(pLab < 1.744)
1647 sigma = 26.41*std::pow(pLab-0.812 1627 sigma = 26.41*std::pow(pLab-0.812,7.138)/std::pow(pLab,5.337);
1648 else if(pLab < 3.728) 1628 else if(pLab < 3.728)
1649 sigma = 1572.*std::pow(pLab-0.812 1629 sigma = 1572.*std::pow(pLab-0.812,9.069)/std::pow(pLab,12.44);
1650 else 1630 else
1651 sigma = 60.23*std::pow(pLab-0.812 1631 sigma = 60.23*std::pow(pLab-0.812,5.084)/std::pow(pLab,6.72);
1652 1632
1653 if(iso == 0) 1633 if(iso == 0)
1654 sigma *= 3.125; 1634 sigma *= 3.125;
1655 else 1635 else
1656 sigma *= 2.875; 1636 sigma *= 2.875;
1657 1637
1658 return sigma; 1638 return sigma;
1659 } 1639 }
1660 1640
1661 /// \brief NKb cross sections 1641 /// \brief NKb cross sections
1662 1642
1663 G4double CrossSectionsStrangeness::NKbToN 1643 G4double CrossSectionsStrangeness::NKbToNKb(Particle const * const p1, Particle const * const p2) {
1664 // 1644 //
1665 // Nucleon-antiKaon quasi-elasti 1645 // Nucleon-antiKaon quasi-elastic cross sections
1666 // 1646 //
1667 // assert((p1->isNucleon() && p2->isAntiKaon( 1647 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1668 1648
1669 G4double sigma=0.; 1649 G4double sigma=0.;
1670 const G4int iso=ParticleTable::getIso 1650 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1671 1651
1672 const Particle *antikaon; 1652 const Particle *antikaon;
1673 const Particle *nucleon; 1653 const Particle *nucleon;
1674 1654
1675 if (p1->isAntiKaon()) { 1655 if (p1->isAntiKaon()) {
1676 antikaon = p1; 1656 antikaon = p1;
1677 nucleon = p2; 1657 nucleon = p2;
1678 } 1658 }
1679 else { 1659 else {
1680 antikaon = p2; 1660 antikaon = p2;
1681 nucleon = p1; 1661 nucleon = p1;
1682 } 1662 }
1683 1663
1684 const G4double pLab = 0.001*Kinematic 1664 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1685 1665
1686 if(iso != 0) // K0b p and K- n -> for 1666 if(iso != 0) // K0b p and K- n -> forbidden: quasi-elastic diffusion only
1687 return 0; 1667 return 0;
1688 else if(nucleon->getType() == Proton) 1668 else if(nucleon->getType() == Proton){ // K- p -> K0b n
1689 if(pLab < 0.08921) 1669 if(pLab < 0.08921)
1690 return 0.; 1670 return 0.;
1691 else if(pLab < 0.2) 1671 else if(pLab < 0.2)
1692 sigma = 0.4977*std::pow(pLab 1672 sigma = 0.4977*std::pow(pLab - 0.08921,0.5581)/std::pow(pLab,2.704);
1693 else if(pLab < 0.73) 1673 else if(pLab < 0.73)
1694 sigma = 2.*std::pow(pLab,-1.2 1674 sigma = 2.*std::pow(pLab,-1.2) + 6.493*std::exp(-0.5*std::pow((pLab-0.3962)/0.02,2));
1695 else if(pLab < 1.38) 1675 else if(pLab < 1.38)
1696 sigma = 2.3*std::pow(pLab,-0. 1676 sigma = 2.3*std::pow(pLab,-0.9) + 1.1*std::exp(-0.5*std::pow((pLab-0.82)/0.04,2)) + 5.*std::exp(-0.5*std::pow((pLab-1.04)/0.1,2));
1697 else if(pLab < 30.) 1677 else if(pLab < 30.)
1698 sigma = 2.5*std::pow(pLab,-1. 1678 sigma = 2.5*std::pow(pLab,-1.68) + 0.7*std::exp(-0.5*std::pow((pLab-1.6)/0.2,2)) + 0.2*std::exp(-0.5*std::pow((pLab-2.3)/0.2,2));
1699 else sigma = 0.; 1679 else sigma = 0.;
1700 } 1680 }
1701 else{ // K0b n -> K- p (same as K- p 1681 else{ // K0b n -> K- p (same as K- p but without threshold)
1702 if(pLab < 0.1) 1682 if(pLab < 0.1)
1703 sigma = 30.; 1683 sigma = 30.;
1704 else if(pLab < 0.73) 1684 else if(pLab < 0.73)
1705 sigma = 2.*std::pow(pLab,-1.2 1685 sigma = 2.*std::pow(pLab,-1.2) + 6.493*std::exp(-0.5*std::pow((pLab-0.3962)/0.02,2));
1706 else if(pLab < 1.38) 1686 else if(pLab < 1.38)
1707 sigma = 2.3*std::pow(pLab,-0. 1687 sigma = 2.3*std::pow(pLab,-0.9) + 1.1*std::exp(-0.5*std::pow((pLab-0.82)/0.04,2)) + 5.*std::exp(-0.5*std::pow((pLab-1.04)/0.1,2));
1708 else if(pLab < 30.) 1688 else if(pLab < 30.)
1709 sigma = 2.5*std::pow(pLab,-1. 1689 sigma = 2.5*std::pow(pLab,-1.68) + 0.7*std::exp(-0.5*std::pow((pLab-1.6)/0.2,2)) + 0.2*std::exp(-0.5*std::pow((pLab-2.3)/0.2,2));
1710 else sigma = 0.; 1690 else sigma = 0.;
1711 } 1691 }
1712 return sigma; 1692 return sigma;
1713 } 1693 }
1714 1694
1715 G4double CrossSectionsStrangeness::NKbToS 1695 G4double CrossSectionsStrangeness::NKbToSpi(Particle const * const p1, Particle const * const p2) {
1716 // 1696 //
1717 // Nucleon-antiKaon producing Si 1697 // Nucleon-antiKaon producing Sigma-pion cross sections
1718 // 1698 //
1719 // ratio 1699 // ratio
1720 // p K0b (4/3) p K- (13/6) 1700 // p K0b (4/3) p K- (13/6)
1721 // 1701 //
1722 // p K0b -> S+ pi0 (2/3) 1702 // p K0b -> S+ pi0 (2/3)
1723 // p K0b -> S0 pi+ (2/3) 1703 // p K0b -> S0 pi+ (2/3)
1724 // p K- -> S+ pi- (1) 1704 // p K- -> S+ pi- (1)
1725 // p K- -> S0 pi0 (1/2) 1705 // p K- -> S0 pi0 (1/2)
1726 // p K- -> S- pi+ (2/3) 1706 // p K- -> S- pi+ (2/3)
1727 1707
1728 // assert((p1->isNucleon() && p2->isAntiKaon( 1708 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1729 1709
1730 G4double sigma=0.; 1710 G4double sigma=0.;
1731 const G4int iso=ParticleTable::getIso 1711 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1732 1712
1733 const Particle *antikaon; 1713 const Particle *antikaon;
1734 const Particle *nucleon; 1714 const Particle *nucleon;
1735 1715
1736 if (p1->isAntiKaon()) { 1716 if (p1->isAntiKaon()) {
1737 antikaon = p1; 1717 antikaon = p1;
1738 nucleon = p2; 1718 nucleon = p2;
1739 } 1719 }
1740 else { 1720 else {
1741 antikaon = p2; 1721 antikaon = p2;
1742 nucleon = p1; 1722 nucleon = p1;
1743 } 1723 }
1744 1724
1745 const G4double pLab = 0.001*Kinematic 1725 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1746 1726
1747 if(iso == 0){ 1727 if(iso == 0){
1748 if(pLab < 0.1) 1728 if(pLab < 0.1)
1749 return 152.0; // 70.166*13/6 1729 return 152.0; // 70.166*13/6
1750 else 1730 else
1751 sigma = 13./6.*(1.4*std::pow( 1731 sigma = 13./6.*(1.4*std::pow(pLab,-1.7)+1.88*std::exp(-std::pow(pLab-0.747,2)/0.005)+8*std::exp(-std::pow(pLab-0.4,2)/0.002)+0.8*std::exp(-std::pow(pLab-1.07,2)/0.01));
1752 } 1732 }
1753 else{ 1733 else{
1754 if(pLab < 0.1) 1734 if(pLab < 0.1)
1755 return 93.555; // 70.166*4/3 1735 return 93.555; // 70.166*4/3
1756 else 1736 else
1757 sigma = 4./3.*(1.4*std::pow(p 1737 sigma = 4./3.*(1.4*std::pow(pLab,-1.7)+1.88*std::exp(-std::pow(pLab-0.747,2)/0.005)+8*std::exp(-std::pow(pLab-0.4,2)/0.002)+0.8*std::exp(-std::pow(pLab-1.07,2)/0.01));
1758 //sigma = 4./3.*(1.4*std::pow 1738 //sigma = 4./3.*(1.4*std::pow(pLab,-1.7)+1.88*std::exp(-std::pow(pLab-0.747,2)/0.005)+0.8*std::exp(-std::pow(pLab-1.07,2)/0.01));
1759 } 1739 }
1760 1740
1761 return sigma; 1741 return sigma;
1762 } 1742 }
1763 1743
1764 G4double CrossSectionsStrangeness::NKbToL 1744 G4double CrossSectionsStrangeness::NKbToLpi(Particle const * const p1, Particle const * const p2) {
1765 // 1745 //
1766 // Nucleon-antiKaon producing La 1746 // Nucleon-antiKaon producing Lambda-pion cross sections
1767 // 1747 //
1768 // ratio 1748 // ratio
1769 // p K0b (1) p K- (1/2) 1749 // p K0b (1) p K- (1/2)
1770 // 1750 //
1771 // p K- -> L pi0 (1/2) 1751 // p K- -> L pi0 (1/2)
1772 // p K0b -> L pi+ (1) 1752 // p K0b -> L pi+ (1)
1773 // assert((p1->isNucleon() && p2->isAntiKaon( 1753 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1774 1754
1775 G4double sigma = 0.; 1755 G4double sigma = 0.;
1776 const G4int iso=ParticleTable::getIso 1756 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1777 1757
1778 const Particle *antikaon; 1758 const Particle *antikaon;
1779 const Particle *nucleon; 1759 const Particle *nucleon;
1780 1760
1781 if (p1->isAntiKaon()) { 1761 if (p1->isAntiKaon()) {
1782 antikaon = p1; 1762 antikaon = p1;
1783 nucleon = p2; 1763 nucleon = p2;
1784 } 1764 }
1785 else { 1765 else {
1786 antikaon = p2; 1766 antikaon = p2;
1787 nucleon = p1; 1767 nucleon = p1;
1788 } 1768 }
1789 if(iso == 0) 1769 if(iso == 0)
1790 sigma = p_kmToL_pz(antikaon,nucle 1770 sigma = p_kmToL_pz(antikaon,nucleon);
1791 else 1771 else
1792 sigma = 2*p_kmToL_pz(antikaon,nuc 1772 sigma = 2*p_kmToL_pz(antikaon,nucleon);
1793 1773
1794 return sigma; 1774 return sigma;
1795 } 1775 }
1796 G4double CrossSectionsStrangeness::p_kmTo 1776 G4double CrossSectionsStrangeness::p_kmToL_pz(Particle const * const p1, Particle const * const p2) {
1797 const G4double pLab = 0.001*Kinematic 1777 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(p1, p2); // GeV
1798 G4double sigma = 0.; 1778 G4double sigma = 0.;
1799 if(pLab < 0.086636) 1779 if(pLab < 0.086636)
1800 sigma = 40.24; 1780 sigma = 40.24;
1801 else if(pLab < 0.5) 1781 else if(pLab < 0.5)
1802 sigma = 0.97*std::pow(pLab,-1.523 1782 sigma = 0.97*std::pow(pLab,-1.523);
1803 else if(pLab < 2.) 1783 else if(pLab < 2.)
1804 sigma = 1.23*std::pow(pLab,-1.467 1784 sigma = 1.23*std::pow(pLab,-1.467)+0.872*std::exp(-std::pow(pLab-0.749,2)/0.0045)+2.337*std::exp(-std::pow(pLab-0.957,2)/0.017)+0.476*std::exp(-std::pow(pLab-1.434,2)/0.136);
1805 else if(pLab < 30.) 1785 else if(pLab < 30.)
1806 sigma = 3.*std::pow(pLab,-2.57); 1786 sigma = 3.*std::pow(pLab,-2.57);
1807 else 1787 else
1808 sigma = 0.; 1788 sigma = 0.;
1809 1789
1810 return sigma; 1790 return sigma;
1811 } 1791 }
1812 1792
1813 G4double CrossSectionsStrangeness::NKbToS 1793 G4double CrossSectionsStrangeness::NKbToS2pi(Particle const * const p1, Particle const * const p2) {
1814 // 1794 //
1815 // Nucleon-antiKaon producing Si 1795 // Nucleon-antiKaon producing Sigma-2pion cross sections
1816 // 1796 //
1817 // ratio 1797 // ratio
1818 // p K0b (29/12) p K- (59/24) 1798 // p K0b (29/12) p K- (59/24)
1819 // 1799 //
1820 // p K0b -> S+ pi+ pi- (2/3) 1800 // p K0b -> S+ pi+ pi- (2/3)
1821 // p K0b -> S+ pi0 pi0 (1/4) 1801 // p K0b -> S+ pi0 pi0 (1/4)
1822 // p K0b -> S0 pi+ pi0 (5/6) 1802 // p K0b -> S0 pi+ pi0 (5/6)
1823 // p K0b -> S- pi+ pi+ (2/3) 1803 // p K0b -> S- pi+ pi+ (2/3)
1824 // p K- -> S+ pi0 pi- (1) 1804 // p K- -> S+ pi0 pi- (1)
1825 // p K- -> S0 pi+ pi- (2/3) 1805 // p K- -> S0 pi+ pi- (2/3)
1826 // p K- -> S0 pi0 pi0 (1/8) 1806 // p K- -> S0 pi0 pi0 (1/8)
1827 // p K- -> S- pi+ pi0 (2/3) 1807 // p K- -> S- pi+ pi0 (2/3)
1828 1808
1829 // assert((p1->isNucleon() && p2->isAntiKaon( 1809 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1830 1810
1831 G4double sigma=0.; 1811 G4double sigma=0.;
1832 const G4int iso=ParticleTable::getIso 1812 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1833 1813
1834 const Particle *antikaon; 1814 const Particle *antikaon;
1835 const Particle *nucleon; 1815 const Particle *nucleon;
1836 1816
1837 if (p1->isAntiKaon()) { 1817 if (p1->isAntiKaon()) {
1838 antikaon = p1; 1818 antikaon = p1;
1839 nucleon = p2; 1819 nucleon = p2;
1840 } 1820 }
1841 else { 1821 else {
1842 antikaon = p2; 1822 antikaon = p2;
1843 nucleon = p1; 1823 nucleon = p1;
1844 } 1824 }
1845 1825
1846 const G4double pLab = 0.001*Kinematic 1826 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1847 1827
1848 if(pLab < 0.260) 1828 if(pLab < 0.260)
1849 return 0.; 1829 return 0.;
1850 1830
1851 if(iso == 0) 1831 if(iso == 0)
1852 sigma = 29./12.*3./2.*(49.96*std: 1832 sigma = 29./12.*3./2.*(49.96*std::pow(pLab-0.260,6.398)/std::pow(pLab+0.260,9.732)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));
1853 else 1833 else
1854 sigma = 54./24.*3./2.*(49.96*std: 1834 sigma = 54./24.*3./2.*(49.96*std::pow(pLab-0.260,6.398)/std::pow(pLab+0.260,9.732)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));
1855 1835
1856 /* 1836 /*
1857 if(iso == 0) 1837 if(iso == 0)
1858 sigma = 29./12.*(85.46*std::pow(p 1838 sigma = 29./12.*(85.46*std::pow(pLab-0.226,8.118)/std::pow(pLab+0.226,11.69)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));
1859 else 1839 else
1860 sigma = 54./24.*(85.46*std::pow(p 1840 sigma = 54./24.*(85.46*std::pow(pLab-0.226,8.118)/std::pow(pLab+0.226,11.69)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));*/
1861 1841
1862 return sigma; 1842 return sigma;
1863 } 1843 }
1864 1844
1865 G4double CrossSectionsStrangeness::NKbToL 1845 G4double CrossSectionsStrangeness::NKbToL2pi(Particle const * const p1, Particle const * const p2) {
1866 // 1846 //
1867 // Nucleon-antiKaon producing La 1847 // Nucleon-antiKaon producing Lambda-2pion cross sections
1868 // 1848 //
1869 // ratio 1849 // ratio
1870 // p K0b -> L pi+ pi0 (1) 1850 // p K0b -> L pi+ pi0 (1)
1871 // p K- -> L pi+ pi- (1) 1851 // p K- -> L pi+ pi- (1)
1872 // p K- -> L pi0 pi0 (1/4) 1852 // p K- -> L pi0 pi0 (1/4)
1873 1853
1874 // assert((p1->isNucleon() && p2->isAntiKaon( 1854 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1875 1855
1876 G4double sigma = 0.; 1856 G4double sigma = 0.;
1877 const G4int iso=ParticleTable::getIso 1857 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1878 1858
1879 const Particle *antikaon; 1859 const Particle *antikaon;
1880 const Particle *nucleon; 1860 const Particle *nucleon;
1881 1861
1882 if (p1->isAntiKaon()) { 1862 if (p1->isAntiKaon()) {
1883 antikaon = p1; 1863 antikaon = p1;
1884 nucleon = p2; 1864 nucleon = p2;
1885 } 1865 }
1886 else { 1866 else {
1887 antikaon = p2; 1867 antikaon = p2;
1888 nucleon = p1; 1868 nucleon = p1;
1889 } 1869 }
1890 1870
1891 if(iso == 0) 1871 if(iso == 0)
1892 sigma = 1.25*p_kmToL_pp_pm(antika 1872 sigma = 1.25*p_kmToL_pp_pm(antikaon,nucleon);
1893 else 1873 else
1894 sigma = p_kmToL_pp_pm(antikaon,nu 1874 sigma = p_kmToL_pp_pm(antikaon,nucleon);
1895 1875
1896 return sigma; 1876 return sigma;
1897 } 1877 }
1898 G4double CrossSectionsStrangeness::p_kmTo 1878 G4double CrossSectionsStrangeness::p_kmToL_pp_pm(Particle const * const p1, Particle const * const p2) {
1899 const G4double pLab = 0.001*Kinematic 1879 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(p1, p2); // GeV
1900 G4double sigma = 0.; 1880 G4double sigma = 0.;
1901 if(pLab < 0.97) 1881 if(pLab < 0.97)
1902 sigma = 6364.*std::pow(pLab,6.07) 1882 sigma = 6364.*std::pow(pLab,6.07)/std::pow(pLab+1.,10.58)+2.158*std::exp(-std::pow((pLab-0.395)/.01984,2)/2.);
1903 else if(pLab < 30) 1883 else if(pLab < 30)
1904 sigma = 46.3*std::pow(pLab,0.62)/ 1884 sigma = 46.3*std::pow(pLab,0.62)/std::pow(pLab+1.,3.565);
1905 else 1885 else
1906 sigma = 0.; 1886 sigma = 0.;
1907 1887
1908 return sigma; 1888 return sigma;
1909 } 1889 }
1910 1890
1911 G4double CrossSectionsStrangeness::NKbToN 1891 G4double CrossSectionsStrangeness::NKbToNKbpi(Particle const * const p1, Particle const * const p2) {
1912 // 1892 //
1913 // Nucleon-antiKaon producing Nu 1893 // Nucleon-antiKaon producing Nucleon-antiKaon-pion cross sections
1914 // 1894 //
1915 // ratio 1895 // ratio
1916 // p K- (28) p K0b (20) << 1896 // p K0b (2) p K- (4)
1917 // <<
1918 // p K- -> p K- pi0 (6)* <<
1919 // p K- -> p K0b pi- (7)* <<
1920 // p K- -> n K- pi+ (9)* <<
1921 // p K- -> n K0b pi0 (6) <<
1922 // p K0b -> p K0b pi0 (4) <<
1923 // p K0b -> p K- pi+ (10)* <<
1924 // p K0b -> n K0b pi+ (6)* <<
1925 // 1897 //
>> 1898 // p K0b -> p K0b pi0 (1/2)
>> 1899 // p K0b -> p K- pi+ (1)*
>> 1900 // p K0b -> n K0b pi+ (1/2)*
>> 1901 // p K- -> p K- pi0 (1/2)*
>> 1902 // p K- -> p K0b pi- (2/3)*
>> 1903 // p K- -> n K- pi+ (3/4)*
>> 1904 // p K- -> n K0b pi0 (2)
1926 1905
>> 1906 //
1927 // assert((p1->isNucleon() && p2->isAntiKaon( 1907 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1928 1908
1929 G4double sigma=0.; 1909 G4double sigma=0.;
1930 const G4int iso=ParticleTable::getIso 1910 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1931 1911
1932 const Particle *antikaon; 1912 const Particle *antikaon;
1933 const Particle *nucleon; 1913 const Particle *nucleon;
1934 1914
1935 if (p1->isAntiKaon()) { 1915 if (p1->isAntiKaon()) {
1936 antikaon = p1; 1916 antikaon = p1;
1937 nucleon = p2; 1917 nucleon = p2;
1938 } 1918 }
1939 else { 1919 else {
1940 antikaon = p2; 1920 antikaon = p2;
1941 nucleon = p1; 1921 nucleon = p1;
1942 } 1922 }
1943 1923
1944 const G4double pLab = 0.001*Kinematic 1924 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1945 1925
1946 if(pLab < 0.526) 1926 if(pLab < 0.526)
1947 return 0.; 1927 return 0.;
1948 1928
1949 if(iso == 0) 1929 if(iso == 0)
1950 sigma = 28. * 10.13*std::pow(pLab << 1930 sigma = 2. * 101.3*std::pow(pLab-0.526,5.846)/std::pow(pLab,8.343);
1951 else 1931 else
1952 sigma = 20. * 10.13*std::pow(pLab << 1932 sigma = 4. * 101.3*std::pow(pLab-0.526,5.846)/std::pow(pLab,8.343);
1953 1933
1954 return sigma; 1934 return sigma;
1955 } 1935 }
1956 1936
1957 G4double CrossSectionsStrangeness::NKbToN 1937 G4double CrossSectionsStrangeness::NKbToNKb2pi(Particle const * const p1, Particle const * const p2) {
1958 // 1938 //
1959 // Nucleon-antiKaon producing Nu 1939 // Nucleon-antiKaon producing Nucleon-antiKaon-2pion cross sections
1960 // 1940 //
1961 // ratio 1941 // ratio
1962 // p K0b (4.25) p K- (4.75) 1942 // p K0b (4.25) p K- (4.75)
1963 // 1943 //
1964 // p K0b -> p K0b pi+ pi- (1) 1944 // p K0b -> p K0b pi+ pi- (1)
1965 // p K0b -> p K0b pi0 pi0 (1/4) 1945 // p K0b -> p K0b pi0 pi0 (1/4)
1966 // p K0b -> p K- pi+ pi0 (1) 1946 // p K0b -> p K- pi+ pi0 (1)
1967 // p K0b -> n K0b pi+ pi0 (1) 1947 // p K0b -> n K0b pi+ pi0 (1)
1968 // p K0b -> n K- pi+ pi+ (1) 1948 // p K0b -> n K- pi+ pi+ (1)
1969 // p K- -> p K0b pi0 pi- (1) 1949 // p K- -> p K0b pi0 pi- (1)
1970 // p K- -> p K- pi+ pi- (1) 1950 // p K- -> p K- pi+ pi- (1)
1971 // p K- -> p K- pi0 pi0 (1/4) 1951 // p K- -> p K- pi0 pi0 (1/4)
1972 // p K- -> n K0b pi+ pi- (1) 1952 // p K- -> n K0b pi+ pi- (1)
1973 // p K- -> n K0b pi0 pi0 (1/2) 1953 // p K- -> n K0b pi0 pi0 (1/2)
1974 // p K- -> n K- pi+ pi0 (1) 1954 // p K- -> n K- pi+ pi0 (1)
1975 1955
1976 // assert((p1->isNucleon() && p2->isAntiKaon( 1956 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1977 1957
1978 G4double sigma=0.; 1958 G4double sigma=0.;
1979 const G4int iso=ParticleTable::getIso 1959 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1980 1960
1981 const Particle *antikaon; 1961 const Particle *antikaon;
1982 const Particle *nucleon; 1962 const Particle *nucleon;
1983 1963
1984 if (p1->isAntiKaon()) { 1964 if (p1->isAntiKaon()) {
1985 antikaon = p1; 1965 antikaon = p1;
1986 nucleon = p2; 1966 nucleon = p2;
1987 } 1967 }
1988 else { 1968 else {
1989 antikaon = p2; 1969 antikaon = p2;
1990 nucleon = p1; 1970 nucleon = p1;
1991 } 1971 }
1992 1972
1993 const G4double pLab = 0.001*Kinematic 1973 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1994 1974
1995 if(pLab < 0.85) 1975 if(pLab < 0.85)
1996 return 0.; 1976 return 0.;
1997 1977
1998 if(iso == 0) 1978 if(iso == 0)
1999 sigma = 4.75 * 26.8*std::pow(pLab 1979 sigma = 4.75 * 26.8*std::pow(pLab-0.85,4.9)/std::pow(pLab,6.34);
2000 else 1980 else
2001 sigma = 4.25 * 26.8*std::pow(pLab 1981 sigma = 4.25 * 26.8*std::pow(pLab-0.85,4.9)/std::pow(pLab,6.34);
2002 1982
2003 return sigma; 1983 return sigma;
2004 } 1984 }
2005 1985
2006 1986
2007 } // namespace G4INCL 1987 } // namespace G4INCL
2008 1988
2009 1989