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