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