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 // 26 // 14.03.07 V. Grichine - first implementation 27 // 14.03.07 V. Grichine - first implementation 27 // 28 // 28 // 04.09.18 V. Ivantchenko Major revision of i << 29 // 30.09.18 V. Grichine hyperon-nucleon xsc fi << 30 29 31 #include "G4HadronNucleonXsc.hh" 30 #include "G4HadronNucleonXsc.hh" 32 #include "G4Element.hh" << 31 33 #include "G4Proton.hh" << 34 #include "G4Nucleus.hh" << 35 #include "G4PhysicalConstants.hh" << 36 #include "G4SystemOfUnits.hh" << 37 #include "G4ParticleTable.hh" 32 #include "G4ParticleTable.hh" 38 #include "G4IonTable.hh" 33 #include "G4IonTable.hh" 39 #include "G4HadTmpUtil.hh" << 34 #include "G4ParticleDefinition.hh" 40 #include "G4Log.hh" << 35 41 #include "G4Exp.hh" << 36 ////////////////////////////////////////////////////////////////////////////////////// 42 #include "G4Pow.hh" << 37 // 43 #include "G4NuclearRadii.hh" << 38 // 44 << 45 #include "G4Neutron.hh" << 46 #include "G4PionPlus.hh" << 47 #include "G4KaonPlus.hh" << 48 #include "G4KaonMinus.hh" << 49 #include "G4KaonZeroShort.hh" << 50 #include "G4KaonZeroLong.hh" << 51 39 52 namespace << 53 { << 54 const G4double invGeV = 1.0/CLHEP::GeV; << 55 const G4double invGeV2 = 1.0/(CLHEP::GeV*CLH << 56 // PDG fit constants << 57 const G4double minLogP = 3.5; // min of ( << 58 const G4double cofLogE = .0557; // elastic << 59 const G4double cofLogT = .3; // total (l << 60 const G4double pMin = .1; // fast LE << 61 const G4double pMax = 1000.; // fast HE << 62 const G4double ekinmin = 0.1*CLHEP::MeV; / << 63 const G4double ekinmaxQB = 100*CLHEP::MeV; / << 64 } << 65 40 66 G4HadronNucleonXsc::G4HadronNucleonXsc() 41 G4HadronNucleonXsc::G4HadronNucleonXsc() >> 42 : fUpperLimit( 10000 * GeV ), >> 43 fLowerLimit( 0.03 * MeV ) 67 { 44 { 68 // basic hadrons << 45 theGamma = G4Gamma::Gamma(); 69 theProton = G4Proton::Proton(); 46 theProton = G4Proton::Proton(); 70 theNeutron = G4Neutron::Neutron(); 47 theNeutron = G4Neutron::Neutron(); >> 48 theAProton = G4AntiProton::AntiProton(); >> 49 theANeutron = G4AntiNeutron::AntiNeutron(); 71 thePiPlus = G4PionPlus::PionPlus(); 50 thePiPlus = G4PionPlus::PionPlus(); 72 << 51 thePiMinus = G4PionMinus::PionMinus(); 73 // basic strange mesons << 52 thePiZero = G4PionZero::PionZero(); 74 theKPlus = G4KaonPlus::KaonPlus(); 53 theKPlus = G4KaonPlus::KaonPlus(); 75 theKMinus = G4KaonMinus::KaonMinus(); 54 theKMinus = G4KaonMinus::KaonMinus(); 76 theK0S = G4KaonZeroShort::KaonZeroShort 55 theK0S = G4KaonZeroShort::KaonZeroShort(); 77 theK0L = G4KaonZeroLong::KaonZeroLong() 56 theK0L = G4KaonZeroLong::KaonZeroLong(); >> 57 theL = G4Lambda::Lambda(); >> 58 theAntiL = G4AntiLambda::AntiLambda(); >> 59 theSPlus = G4SigmaPlus::SigmaPlus(); >> 60 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus(); >> 61 theSMinus = G4SigmaMinus::SigmaMinus(); >> 62 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus(); >> 63 theS0 = G4SigmaZero::SigmaZero(); >> 64 theAS0 = G4AntiSigmaZero::AntiSigmaZero(); >> 65 theXiMinus = G4XiMinus::XiMinus(); >> 66 theXi0 = G4XiZero::XiZero(); >> 67 theAXiMinus = G4AntiXiMinus::AntiXiMinus(); >> 68 theAXi0 = G4AntiXiZero::AntiXiZero(); >> 69 theOmega = G4OmegaMinus::OmegaMinus(); >> 70 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus(); >> 71 theD = G4Deuteron::Deuteron(); >> 72 theT = G4Triton::Triton(); >> 73 theA = G4Alpha::Alpha(); >> 74 theHe3 = G4He3::He3(); >> 75 } 78 76 79 g4calc = G4Pow::GetInstance(); << 77 /////////////////////////////////////////////////////////////////////////////////////// >> 78 // >> 79 // >> 80 >> 81 G4HadronNucleonXsc::~G4HadronNucleonXsc() >> 82 { 80 } 83 } 81 84 82 void G4HadronNucleonXsc::CrossSectionDescripti << 85 >> 86 //////////////////////////////////////////////////////////////////////////////////////// >> 87 // >> 88 // >> 89 >> 90 >> 91 G4bool >> 92 G4HadronNucleonXsc::IsApplicable(const G4DynamicParticle* aDP, >> 93 const G4Element* anElement) >> 94 { >> 95 return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN()); >> 96 } >> 97 >> 98 //////////////////////////////////////////////////////////////////////////////////////// >> 99 // >> 100 // >> 101 >> 102 G4bool >> 103 G4HadronNucleonXsc::IsZAApplicable(const G4DynamicParticle* aDP, >> 104 G4double Z, G4double) 83 { 105 { 84 outFile << "G4HadronNucleonXsc calculates th << 106 G4bool applicable = false; 85 << "hadron-nucleon cross sections us << 107 // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber(); 86 << "parameterizations within the Gla << 108 G4double kineticEnergy = aDP->GetKineticEnergy(); 87 << "valid for all incident gammas an << 109 88 << "energies above 30 keV. This is << 110 const G4ParticleDefinition* theParticle = aDP->GetDefinition(); 89 << "is to be used to build a cross s << 111 >> 112 if ( ( kineticEnergy >= fLowerLimit && >> 113 Z > 1.5 && // >= He >> 114 ( theParticle == theAProton || >> 115 theParticle == theGamma || >> 116 theParticle == theKPlus || >> 117 theParticle == theKMinus || >> 118 theParticle == theSMinus) ) || >> 119 >> 120 ( kineticEnergy >= 0.1*fLowerLimit && >> 121 Z > 1.5 && // >= He >> 122 ( theParticle == theProton || >> 123 theParticle == theNeutron || >> 124 theParticle == thePiPlus || >> 125 theParticle == thePiMinus ) ) ) applicable = true; >> 126 >> 127 return applicable; 90 } 128 } 91 129 92 G4double G4HadronNucleonXsc::HadronNucleonXsc( << 130 93 const G4Particl << 131 >> 132 >> 133 ///////////////////////////////////////////////////////////////////////////////////// >> 134 // >> 135 // Returns hadron-nucleon Xsc according to differnt parametrisations: >> 136 // [2] E. Levin, hep-ph/9710546 >> 137 // [3] U. Dersch, et al, hep-ex/9910052 >> 138 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725 >> 139 >> 140 G4double >> 141 G4HadronNucleonXsc::GetHadronNucleonXscEL(const G4DynamicParticle* aParticle, >> 142 const G4ParticleDefinition* nucleon ) 94 { 143 { 95 G4double xsc(0.); << 144 G4double xsection; 96 G4int pdg = std::abs(theParticle->GetPDGEnco << 145 >> 146 >> 147 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ??? >> 148 >> 149 G4double proj_mass = aParticle->GetMass(); >> 150 G4double proj_momentum = aParticle->GetMomentum().mag(); >> 151 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum ); >> 152 >> 153 sMand /= GeV*GeV; // in GeV for parametrisation >> 154 proj_momentum /= GeV; >> 155 >> 156 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); >> 157 >> 158 G4bool pORn = (nucleon == theProton || nucleon == theNeutron ); 97 159 98 // p, n, pi+-, pbar, nbar << 160 99 if ( pdg == 2212 || pdg == 2112 || pdg == 21 << 161 if(theParticle == theGamma && pORn ) 100 xsc = HadronNucleonXscNS(theParticle, nucl << 101 } << 102 else if ( pdg == 22 ) // gamma << 103 { 162 { 104 xsc = HadronNucleonXscPDG(theParticle, nuc << 163 xsection = (0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525)); 105 } << 164 } 106 else if ( pdg == 321 || pdg == 310 || pdg == << 165 else if(theParticle == theNeutron && pORn ) // as proton ??? >> 166 { >> 167 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525)); >> 168 } >> 169 else if(theParticle == theProton && pORn ) >> 170 { >> 171 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525)); >> 172 >> 173 // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) ); >> 174 // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) ); >> 175 } >> 176 else if(theParticle == theAProton && pORn ) >> 177 { >> 178 xsection = ( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525)); >> 179 } >> 180 else if(theParticle == thePiPlus && pORn ) 107 { 181 { 108 xsc = KaonNucleonXscNS(theParticle, nucleo << 182 xsection = (13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525)); >> 183 } >> 184 else if(theParticle == thePiMinus && pORn ) >> 185 { >> 186 // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) ); >> 187 xsection = (13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525)); >> 188 } >> 189 else if(theParticle == theKPlus && pORn ) >> 190 { >> 191 xsection = (11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525)); >> 192 } >> 193 else if(theParticle == theKMinus && pORn ) >> 194 { >> 195 xsection = (11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525)); 109 } 196 } 110 else if (pdg > 3000) << 197 else // as proton ??? 111 { 198 { 112 if (pdg == 3122 || pdg == 3222 || pdg == 3 << 199 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525)); 113 pdg == 4122 || pdg == 4332 || pdg == 4122 || << 114 pdg == 5122 || pdg == 5332 || pdg == 5 << 115 ) // heavy s-,c-,b-hyperons << 116 { << 117 xsc = HyperonNucleonXscNS(theParticle, n << 118 } << 119 else << 120 { << 121 xsc = HadronNucleonXscPDG(theParticle, n << 122 } << 123 } else if (pdg > 220) { << 124 if (pdg == 511 || pdg == 421 || pdg == 531 << 125 pdg == 221 || pdg == 331 || pdg == 441 || pd << 126 { << 127 xsc = SCBMesonNucleonXscNS(theParticle, << 128 } << 129 else << 130 { << 131 xsc = HadronNucleonXscPDG(theParticle, n << 132 } << 133 } else { << 134 xsc = HadronNucleonXscPDG(theParticle, nuc << 135 } 200 } 136 return xsc; << 201 xsection *= millibarn; >> 202 >> 203 fTotalXsc = xsection; >> 204 fInelasticXsc = 0.83*xsection; >> 205 fElasticXsc = fTotalXsc - fInelasticXsc; >> 206 if (fElasticXsc < 0.)fElasticXsc = 0.; >> 207 >> 208 return xsection; 137 } 209 } 138 210 139 ////////////////////////////////////////////// << 211 >> 212 >> 213 >> 214 >> 215 ///////////////////////////////////////////////////////////////////////////////////// 140 // 216 // 141 // Returns hadron-nucleon Xsc according to PDG << 217 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005): 142 // http://pdg.lbl.gov/2017/reviews/hadronicrpp << 218 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf >> 219 // At = number of nucleons, Zt = number of protons 143 220 144 G4double G4HadronNucleonXsc::HadronNucleonXscP << 221 G4double 145 const G4ParticleD << 222 G4HadronNucleonXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle, 146 const G4ParticleDefinition* nucleon << 223 const G4ParticleDefinition* nucleon ) 147 { 224 { 148 static const G4double M = 2.1206; // in G << 225 G4double xsection(0); 149 static const G4double eta1 = 0.4473; << 226 G4int Zt=1, Nt=1, At=1; 150 static const G4double eta2 = 0.5486; << 227 151 static const G4double H = 0.272; << 228 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ??? >> 229 >> 230 G4double proj_mass = aParticle->GetMass(); >> 231 G4double proj_momentum = aParticle->GetMomentum().mag(); >> 232 >> 233 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum ); 152 234 153 G4int pdg = theParticle->GetPDGEncoding(); << 235 sMand /= GeV*GeV; // in GeV for parametrisation 154 236 155 G4double mass1 = (pdg == 22) ? 770. : thePar << 237 // General PDG fit constants 156 G4double mass2 = nucleon->GetPDGMass(); << 157 238 158 G4double sMand = CalcMandelstamS(ekin, mass1 << 239 G4double s0 = 5.38*5.38; // in Gev^2 159 G4double x = (mass1 + mass2)*invGeV + M; << 240 G4double eta1 = 0.458; 160 G4double blog = G4Log(sMand/(x*x)); << 241 G4double eta2 = 0.458; >> 242 G4double B = 0.308; 161 243 162 G4double P(0.0), R1(0.0), R2(0.0), del(1.0); << 163 244 164 G4bool proton = (nucleon == theProton); << 245 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); >> 246 >> 247 G4bool pORn = (nucleon == theProton || nucleon == theNeutron ); >> 248 G4bool proton = (nucleon == theProton); 165 G4bool neutron = (nucleon == theNeutron); 249 G4bool neutron = (nucleon == theNeutron); 166 250 167 if(theParticle == theNeutron) << 251 >> 252 if(theParticle == theNeutron) // proton-neutron fit 168 { 253 { 169 if ( proton ) 254 if ( proton ) 170 { 255 { 171 P = 34.71; << 256 xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 172 R1 = 12.52; << 257 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));// on p 173 R2 = -6.66; << 174 } 258 } 175 else << 259 if ( neutron ) 176 { 260 { 177 P = 34.41; << 261 xsection = Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 178 R1 = 13.07; << 262 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // on n pp for nn 179 R2 = -7.394; << 180 } 263 } 181 } 264 } 182 else if(theParticle == theProton) 265 else if(theParticle == theProton) 183 { 266 { 184 if ( neutron ) << 267 if ( proton ) 185 { << 268 { 186 P = 34.71; << 269 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 187 R1 = 12.52; << 270 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); 188 R2 = -6.66; << 189 } << 190 else << 191 { << 192 P = 34.41; << 193 R1 = 13.07; << 194 R2 = -7.394; << 195 } 271 } 196 } << 197 // pbar << 198 else if(pdg == -2212) << 199 { << 200 if ( neutron ) 272 if ( neutron ) 201 { 273 { 202 P = 34.71; << 274 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 203 R1 = 12.52; << 275 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2)); 204 R2 = 6.66; << 205 } << 206 else << 207 { << 208 P = 34.41; << 209 R1 = 13.07; << 210 R2 = 7.394; << 211 } 276 } 212 } 277 } 213 // nbar << 278 else if(theParticle == theAProton) 214 else if(pdg == -2112) << 215 { 279 { 216 if ( proton ) 280 if ( proton ) 217 { << 281 { 218 P = 34.71; << 282 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 219 R1 = 12.52; << 283 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2)); 220 R2 = 6.66; << 221 } 284 } 222 else << 285 if ( neutron ) 223 { 286 { 224 P = 34.41; << 287 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 225 R1 = 13.07; << 288 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2)); 226 R2 = 7.394; << 227 } 289 } 228 } 290 } 229 // pi+ << 291 else if(theParticle == thePiPlus && pORn ) 230 else if(pdg == 211) << 231 { 292 { 232 P = 18.75; << 293 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 233 R1 = 9.56; << 294 + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2)); 234 R2 = -1.767; << 235 } 295 } 236 // pi- << 296 else if(theParticle == thePiMinus && pORn ) 237 else if(pdg == -211) << 238 { 297 { 239 P = 18.75; << 298 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.) 240 R1 = 9.56; << 299 + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2)); 241 R2 = 1.767; << 242 } 300 } 243 else if(theParticle == theKPlus) 301 else if(theParticle == theKPlus) 244 { 302 { 245 if ( proton ) 303 if ( proton ) 246 { << 304 { 247 P = 16.36; << 305 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 248 R1 = 4.29; << 306 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2)); 249 R2 = -3.408; << 250 } 307 } 251 else << 308 if ( neutron ) 252 { 309 { 253 P = 16.31; << 310 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 254 R1 = 3.7; << 311 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2)); 255 R2 = -1.826; << 256 } 312 } 257 } 313 } 258 else if(theParticle == theKMinus) 314 else if(theParticle == theKMinus) 259 { 315 { 260 if ( proton ) 316 if ( proton ) 261 { << 317 { 262 P = 16.36; << 318 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 263 R1 = 4.29; << 319 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2)); 264 R2 = 3.408; << 265 } 320 } 266 else << 321 if ( neutron ) 267 { 322 { 268 P = 16.31; << 323 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 269 R1 = 3.7; << 324 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2) ); 270 R2 = 1.826; << 271 } 325 } 272 } 326 } 273 else if(theParticle == theK0S || theParticle << 327 else if(theParticle == theSMinus && pORn ) 274 { << 275 P = 16.36; << 276 R1 = 2.5; << 277 R2 = 0.; << 278 } << 279 // sigma- << 280 else if(pdg == 3112) << 281 { 328 { 282 P = 34.7; << 329 xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.) 283 R1 = -46.; << 330 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2) ); 284 R2 = 48.; << 285 } 331 } 286 // gamma << 332 else if(theParticle == theGamma && pORn ) // modify later on 287 else if(pdg == 22) // modify later on << 288 { 333 { 289 del= 0.003063; << 334 xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.) 290 P = 34.71*del; << 335 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2) ); 291 R1 = (neutron) ? 0.0231 : 0.0139; << 336 292 R2 = 0.; << 293 } 337 } 294 else // as proton ??? 338 else // as proton ??? 295 { 339 { 296 if ( neutron ) << 340 if ( proton ) 297 { << 341 { 298 P = 34.71; << 342 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 299 R1 = 12.52; << 343 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2) ); 300 R2 = -6.66; << 301 } 344 } 302 else << 345 if ( neutron ) 303 { 346 { 304 P = 34.41; << 347 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 305 R1 = 13.07; << 348 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2)); 306 R2 = -7.394; << 307 } 349 } 308 } 350 } 309 fTotalXsc = CLHEP::millibarn* << 351 xsection *= millibarn; // parametrised in mb 310 (del*(H*blog*blog + P) + R1*G4Exp(-eta1*bl << 352 311 fInelasticXsc = 0.75*fTotalXsc; << 353 fTotalXsc = xsection; >> 354 fInelasticXsc = 0.83*xsection; 312 fElasticXsc = fTotalXsc - fInelasticXsc; 355 fElasticXsc = fTotalXsc - fInelasticXsc; >> 356 if (fElasticXsc < 0.)fElasticXsc = 0.; 313 357 314 if( proton && theParticle->GetPDGCharge() > << 358 return xsection; 315 { << 316 G4double cB = CoulombBarrier(theParticle, << 317 fTotalXsc *= cB; << 318 fElasticXsc *= cB; << 319 fInelasticXsc *= cB; << 320 } << 321 /* << 322 G4cout << "## HadronNucleonXscPDG: ekin(MeV) << 323 << " tot= " << fTotalXsc/CLHEP::millibarn << 324 << " inel= " << fInelasticXsc/CLHEP::millib << 325 << " el= " << fElasticXsc/CLHEP::millibarn << 326 << G4endl; << 327 */ << 328 return fTotalXsc; << 329 } 359 } 330 360 331 ////////////////////////////////////////////// << 361 >> 362 >> 363 ///////////////////////////////////////////////////////////////////////////////////// 332 // 364 // 333 // Returns hadron-nucleon cross-section based 365 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of 334 // data from mainly http://wwwppds.ihep.su:800 366 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database 335 367 336 G4double G4HadronNucleonXsc::HadronNucleonXscN << 368 G4double 337 const G4ParticleD << 369 G4HadronNucleonXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle, 338 const G4ParticleDefinition* nucleon << 370 const G4ParticleDefinition* nucleon ) 339 { << 371 { 340 const G4double ekin = std::max(ekin0, ekinmi << 372 G4double xsection(0), Delta, A0, B0; 341 G4int pdg = theParticle->GetPDGEncoding(); << 373 G4int Zt=1, Nt=1, At=1; 342 /* << 374 G4double hpXsc(0); 343 G4cout<< "HadronNucleonXscNS: Ekin(GeV)= " < << 375 G4double hnXsc(0); 344 << theParticle->GetParticleName() << " << 376 345 << nucleon->GetParticleName() << G4end << 346 */ << 347 if(pdg == -2212 || pdg == -2112) { << 348 return HadronNucleonXscPDG(theParticle, nu << 349 } << 350 << 351 G4double pM = theParticle->GetPDGMass(); << 352 G4double tM = nucleon->GetPDGMass(); << 353 G4double pE = ekin + pM; << 354 G4double pLab = std::sqrt(ekin*(ekin + 2*pM) << 355 << 356 G4double sMand = CalcMandelstamS(ekin, pM, t << 357 << 358 pLab *= invGeV; << 359 pE *= invGeV; << 360 << 361 if(pLab >= 10.) { << 362 fTotalXsc = HadronNucleonXscPDG(theParticl << 363 } else { fTotalXsc = 0.0; } << 364 fElasticXsc = 0.0; << 365 //G4cout << "Stot(mb)= " << fTotalXsc << " p << 366 // << " Smand= " << sMand <<G4endl; << 367 G4double logP = G4Log(pLab); << 368 377 >> 378 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ??? >> 379 >> 380 G4double proj_mass = aParticle->GetMass(); >> 381 G4double proj_energy = aParticle->GetTotalEnergy(); >> 382 G4double proj_momentum = aParticle->GetMomentum().mag(); >> 383 >> 384 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum ); >> 385 >> 386 sMand /= GeV*GeV; // in GeV for parametrisation >> 387 proj_momentum /= GeV; >> 388 proj_energy /= GeV; >> 389 proj_mass /= GeV; >> 390 >> 391 // General PDG fit constants >> 392 >> 393 G4double s0 = 5.38*5.38; // in Gev^2 >> 394 G4double eta1 = 0.458; >> 395 G4double eta2 = 0.458; >> 396 G4double B = 0.308; >> 397 >> 398 >> 399 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); >> 400 >> 401 G4bool pORn = (nucleon == theProton || nucleon == theNeutron ); 369 G4bool proton = (nucleon == theProton); 402 G4bool proton = (nucleon == theProton); 370 G4bool neutron = (nucleon == theNeutron); 403 G4bool neutron = (nucleon == theNeutron); 371 404 372 if( theParticle == theNeutron) << 405 if( theParticle == theNeutron && pORn) 373 { 406 { 374 if( pLab >= 373.) << 407 if( proj_momentum >= 10.) >> 408 // if( proj_momentum >= 2.) 375 { 409 { 376 fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4 << 410 Delta = 1.; 377 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 411 378 } << 412 if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy; 379 else if( pLab >= 100.) << 413 380 { << 414 if(proj_momentum >= 10.) 381 fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G << 415 { 382 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 416 B0 = 7.5; 383 } << 417 A0 = 100. - B0*std::log(3.0e7); 384 else if( pLab >= 10.) << 418 385 { << 419 xsection = A0 + B0*std::log(proj_energy) - 11 386 fElasticXsc = 6 + 20/( (logP-0.182)*(lo << 420 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+ >> 421 0.93827*0.93827,-0.165); // mb >> 422 } >> 423 fTotalXsc = xsection; 387 } 424 } 388 else // pLab < 10 GeV/c << 425 else 389 { 426 { 390 if( neutron ) // nn to be pp << 427 // nn to be pp >> 428 >> 429 if(neutron) 391 { 430 { 392 G4double x = G4Log(pLab/0.73); << 431 if( proj_momentum < 0.73 ) 393 if( pLab < 0.4 ) << 394 { << 395 fTotalXsc = 23 + 50*std::sqrt(g4calc << 396 fElasticXsc = fTotalXsc; << 397 } << 398 else if( pLab < 0.73 ) << 399 { 432 { 400 fTotalXsc = 23 + 50*std::sqrt(g4calc << 433 hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) ); 401 fElasticXsc = fTotalXsc; << 402 } 434 } 403 else if( pLab < 1.05 ) << 435 else if( proj_momentum < 1.05 ) 404 { 436 { 405 fTotalXsc = 23 + 40*x*x; << 437 hnXsc = 23 + 40*(std::log(proj_momentum/0.73))* 406 fElasticXsc = 23 + 20*x*x; << 438 (std::log(proj_momentum/0.73)); 407 } 439 } 408 else // 1.05 - 10 GeV/c << 440 else // if( proj_momentum < 10. ) 409 { 441 { 410 fTotalXsc = 39.0+75*(pLab - 1.2)/(g4 << 442 hnXsc = 39.0+ 411 fElasticXsc = 6 + 20/( (logP-0.182) << 443 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15); 412 } 444 } >> 445 fTotalXsc = hnXsc; 413 } 446 } 414 if( proton ) // pn to be np << 447 // pn to be np >> 448 >> 449 if(proton) 415 { 450 { 416 if( pLab < 0.02 ) << 451 if( proj_momentum < 0.8 ) 417 { << 418 fTotalXsc = 4100+30*G4Exp(G4Log(G4Lo << 419 fElasticXsc = fTotalXsc; << 420 } << 421 else if( pLab < 0.8 ) << 422 { 452 { 423 fTotalXsc = 33+30*g4calc->powN(G4Log << 453 hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0); 424 fElasticXsc = fTotalXsc; << 425 } 454 } 426 else if( pLab < 1.4 ) << 455 else if( proj_momentum < 1.4 ) 427 { 456 { 428 fTotalXsc = 33+30*g4calc->powN(G4Log << 457 hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0); 429 G4double x = G4Log(0.511/pLab); << 430 fElasticXsc = 6 + 52/( x*x + 1.6 ); << 431 } 458 } 432 else // 1.4 < pLab < 10. ) << 459 else // if( proj_momentum < 10. ) 433 { 460 { 434 fTotalXsc = 33.3 + 20.8*(pLab*pLab - << 461 hpXsc = 33.3+ 435 /(std::sqrt(g4calc->powN(pLab,5)) + 0.95 << 462 20.8*(std::pow(proj_momentum,2.0)-1.35)/ 436 fElasticXsc = 6 + 20/( (logP-0.182) << 463 (std::pow(proj_momentum,2.50)+0.95); 437 } 464 } >> 465 fTotalXsc = hpXsc; 438 } 466 } >> 467 // xsection = hpXsc*Zt + hnXsc*Nt; 439 } 468 } 440 } << 469 } 441 ////// proton ////////////////////////////// << 470 else if(theParticle == theProton && pORn) 442 else if(theParticle == theProton) << 443 { 471 { 444 if( pLab >= 373.) // pdg due to TOTEM data << 472 if( proj_momentum >= 10.) 445 { << 473 // if( proj_momentum >= 2.) 446 fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4 << 447 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 448 } << 449 else if( pLab >= 100.) << 450 { 474 { 451 fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G << 475 Delta = 1.; 452 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 476 453 } << 477 if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy; 454 else if( pLab >= 10.) << 478 455 { << 479 if(proj_momentum >= 10.) 456 fElasticXsc = 6. + 20./( (logP-0.182)*(l << 480 { >> 481 B0 = 7.5; >> 482 A0 = 100. - B0*std::log(3.0e7); >> 483 >> 484 xsection = A0 + B0*std::log(proj_energy) - 11 >> 485 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+ >> 486 0.93827*0.93827,-0.165); // mb >> 487 } >> 488 fTotalXsc = xsection; 457 } 489 } 458 else 490 else 459 { 491 { 460 // pp 492 // pp 461 if( proton ) << 493 >> 494 if(proton) 462 { 495 { 463 if( pLab < 0.73 ) << 496 if( proj_momentum < 0.73 ) 464 { 497 { 465 fTotalXsc = 23 + 50*std::sqrt(g4calc << 498 hpXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) ); 466 fElasticXsc = fTotalXsc; << 467 } 499 } 468 else if( pLab < 1.05 ) << 500 else if( proj_momentum < 1.05 ) 469 { 501 { 470 G4double x = G4Log(pLab/0.73); << 502 hpXsc = 23 + 40*(std::log(proj_momentum/0.73))* 471 fTotalXsc = 23 + 40*x*x; << 503 (std::log(proj_momentum/0.73)); 472 fElasticXsc = 23 + 20*x*x; << 473 } 504 } 474 else // 1.05 - 10 GeV/c << 505 else // if( proj_momentum < 10. ) 475 { 506 { 476 fTotalXsc = 39.0+75*(pLab - 1.2)/(g4 << 507 hpXsc = 39.0+ 477 fElasticXsc = 6. + 20./( (logP-0.182 << 508 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15); 478 } 509 } >> 510 fTotalXsc = hpXsc; 479 } 511 } 480 else if( neutron ) // pn to be np << 512 // pn to be np >> 513 >> 514 if(neutron) 481 { 515 { 482 if( pLab < 0.02 ) << 516 if( proj_momentum < 0.8 ) 483 { << 484 fTotalXsc = 4100+30*G4Exp(G4Log(G4Lo << 485 fElasticXsc = fTotalXsc; << 486 } << 487 else if( pLab < 0.8 ) << 488 { 517 { 489 fTotalXsc = 33+30*g4calc->powN(G4Log << 518 hnXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0); 490 fElasticXsc = fTotalXsc; << 491 } 519 } 492 else if( pLab < 1.4 ) << 520 else if( proj_momentum < 1.4 ) 493 { 521 { 494 G4double x1 = G4Log(pLab/0.95); << 522 hnXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0); 495 G4double x2 = G4Log(0.511/pLab); << 496 fTotalXsc = 33+30*x1*x1; << 497 fElasticXsc = 6 + 52/(x2*x2 + 1.6); << 498 } 523 } 499 else // 1.4 < pLab < 10. ) << 524 else // if( proj_momentum < 10. ) 500 { 525 { 501 fTotalXsc = 33.3 + 20.8*(pLab*pLab - << 526 hnXsc = 33.3+ 502 /(std::sqrt(g4calc->powN(pLab,5)) + 0.95 << 527 20.8*(std::pow(proj_momentum,2.0)-1.35)/ 503 fElasticXsc = 6. + 20./( (logP-0.182 << 528 (std::pow(proj_momentum,2.50)+0.95); 504 } 529 } >> 530 fTotalXsc = hnXsc; 505 } 531 } >> 532 // xsection = hpXsc*Zt + hnXsc*Nt; >> 533 // xsection = hpXsc*(Zt + Nt); >> 534 // xsection = hnXsc*(Zt + Nt); 506 } 535 } 507 } << 536 // xsection *= 0.95; 508 // pi+ p; pi- n << 537 } 509 else if((pdg == 211 && proton) || (pdg == -2 << 538 else if(theParticle == theAProton && pORn) 510 { 539 { 511 if( pLab < 0.28 ) << 540 if(proton) 512 { << 513 fTotalXsc = 10./((logP + 1.273)*(logP + << 514 fElasticXsc = fTotalXsc; << 515 } << 516 else if( pLab < 0.68 ) << 517 { << 518 fTotalXsc = 14./((logP + 1.273)*(logP + << 519 fElasticXsc = fTotalXsc; << 520 } << 521 else if( pLab < 0.85 ) << 522 { << 523 G4double x = G4Log(pLab/0.77); << 524 fTotalXsc = 88.*x*x + 14.9; << 525 fElasticXsc = fTotalXsc*G4Exp(-3.*(pLab << 526 } << 527 else if( pLab < 1.15 ) << 528 { 541 { 529 G4double x = G4Log(pLab/0.77); << 542 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 530 fTotalXsc = 88.*x*x + 14.9; << 543 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2)); 531 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( << 532 } 544 } 533 else if( pLab < 1.4) // ns original << 545 if(proton) 534 { 546 { 535 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p << 547 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 536 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL << 548 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2)); 537 fTotalXsc = Ex1 + Ex2 + 27.5; << 538 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( << 539 } 549 } 540 else if( pLab < 2.0 ) // ns original << 550 fTotalXsc = xsection; >> 551 } >> 552 else if(theParticle == thePiPlus && pORn) >> 553 { >> 554 if(proton) 541 { 555 { 542 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p << 556 if(proj_momentum < 0.4) 543 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL << 557 { 544 fTotalXsc = Ex1 + Ex2 + 27.5; << 558 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085); 545 fElasticXsc = 3.0 + 1.36/( (logP - 0.336 << 559 hpXsc = Ex3+20.0; 546 } << 560 } 547 else if( pLab < 3.5 ) // ns original << 561 else if(proj_momentum < 1.15) >> 562 { >> 563 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75)); >> 564 hpXsc = Ex4+14.0; >> 565 } >> 566 else if(proj_momentum < 3.5) >> 567 { >> 568 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55); >> 569 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225); >> 570 hpXsc = Ex1+Ex2+27.5; >> 571 } >> 572 else // if(proj_momentum > 3.5) // mb >> 573 { >> 574 hpXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43); >> 575 } >> 576 fTotalXsc = hpXsc; >> 577 } >> 578 >> 579 // pi+n = pi-p?? >> 580 >> 581 if(neutron) 548 { 582 { 549 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p << 583 if(proj_momentum < 0.37) 550 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL << 584 { 551 fTotalXsc = Ex1 + Ex2 + 27.5; << 585 hnXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07); 552 fElasticXsc = 3.0 + 6.20/( (logP - 0.336 << 586 } >> 587 else if(proj_momentum<0.65) >> 588 { >> 589 hnXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48)); >> 590 } >> 591 else if(proj_momentum<1.3) >> 592 { >> 593 hnXsc = 36.1+ >> 594 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+ >> 595 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075); >> 596 } >> 597 else if(proj_momentum<3.0) >> 598 { >> 599 hnXsc = 36.1+0.079-4.313*std::log(proj_momentum)+ >> 600 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+ >> 601 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12); >> 602 } >> 603 else // mb >> 604 { >> 605 hnXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); >> 606 } >> 607 fTotalXsc = hnXsc; 553 } 608 } 554 else if( pLab < 10. ) // my << 609 // xsection = hpXsc*Zt + hnXsc*Nt; >> 610 } >> 611 else if(theParticle == thePiMinus && pORn) >> 612 { >> 613 // pi-n = pi+p?? >> 614 >> 615 if(neutron) 555 { 616 { 556 fTotalXsc = 10.6 + 2.*G4Log(pE) + 25*G4E << 617 if(proj_momentum < 0.4) 557 fElasticXsc = 3.0 + 6.20/( (logP - 0.336 << 618 { >> 619 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085); >> 620 hnXsc = Ex3+20.0; >> 621 } >> 622 else if(proj_momentum < 1.15) >> 623 { >> 624 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75)); >> 625 hnXsc = Ex4+14.0; >> 626 } >> 627 else if(proj_momentum < 3.5) >> 628 { >> 629 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55); >> 630 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225); >> 631 hnXsc = Ex1+Ex2+27.5; >> 632 } >> 633 else // if(proj_momentum > 3.5) // mb >> 634 { >> 635 hnXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43); >> 636 } >> 637 fTotalXsc = hnXsc; 558 } 638 } 559 else // pLab > 10 // my << 639 // pi-p >> 640 >> 641 if(proton) 560 { 642 { 561 fElasticXsc = 3.0 + 6.20/( (logP - 0.336 << 643 if(proj_momentum < 0.37) >> 644 { >> 645 hpXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07); >> 646 } >> 647 else if(proj_momentum<0.65) >> 648 { >> 649 hpXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48)); >> 650 } >> 651 else if(proj_momentum<1.3) >> 652 { >> 653 hpXsc = 36.1+ >> 654 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+ >> 655 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075); >> 656 } >> 657 else if(proj_momentum<3.0) >> 658 { >> 659 hpXsc = 36.1+0.079-4.313*std::log(proj_momentum)+ >> 660 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+ >> 661 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12); >> 662 } >> 663 else // mb >> 664 { >> 665 hpXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); >> 666 } >> 667 fTotalXsc = hpXsc; 562 } 668 } 563 } << 669 // xsection = hpXsc*Zt + hnXsc*Nt; 564 // pi+ n; pi- p << 670 } 565 else if((pdg == 211 && neutron) || (pdg == - << 671 else if(theParticle == theKPlus && pORn) 566 { 672 { 567 if( pLab < 0.28 ) << 673 if(proton) 568 { << 569 fTotalXsc = 0.288/((pLab - 0.28)*(pLab - << 570 fElasticXsc = 1.8/((logP + 1.273)*(logP << 571 } << 572 else if( pLab < 0.395676 ) // first peak << 573 { 674 { 574 fTotalXsc = 0.648/((pLab - 0.28)*(pLab - << 675 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 575 fElasticXsc = 0.257/((pLab - 0.28)*(pLab << 676 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2)); 576 } 677 } 577 else if( pLab < 0.5 ) << 678 if(neutron) 578 { 679 { 579 G4double y = G4Log(pLab/0.48); << 680 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 580 fTotalXsc = 26 + 110*y*y; << 681 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2)); 581 fElasticXsc = 0.37*fTotalXsc; << 582 } 682 } 583 else if( pLab < 0.65 ) << 683 fTotalXsc = xsection; >> 684 } >> 685 else if(theParticle == theKMinus && pORn) >> 686 { >> 687 if(proton) 584 { 688 { 585 G4double x = G4Log(pLab/0.48); << 689 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) 586 fTotalXsc = 26. + 110.*x*x; << 690 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2)); 587 fElasticXsc = 0.95/((pLab - 0.72)*(pLab << 588 } 691 } 589 else if( pLab < 0.72 ) << 692 if(neutron) 590 { 693 { 591 fTotalXsc = 36.1 + 10*G4Exp(-(pLab-0.72) << 694 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) 592 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 695 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2)); 593 fElasticXsc = 0.95/((pLab - 0.72)*(pLab << 594 } 696 } 595 else if( pLab < 0.88 ) << 697 fTotalXsc = xsection; >> 698 } >> 699 else if(theParticle == theSMinus && pORn) >> 700 { >> 701 xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.) >> 702 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2)); >> 703 } >> 704 else if(theParticle == theGamma && pORn) // modify later on >> 705 { >> 706 xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.) >> 707 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2)); >> 708 fTotalXsc = xsection; >> 709 } >> 710 else // as proton ??? >> 711 { >> 712 if(proton) 596 { 713 { 597 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 << 714 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 598 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 715 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); 599 fElasticXsc = 0.95/((pLab - 0.72)*(pLab << 600 } 716 } 601 else if( pLab < 1.03 ) << 717 if(neutron) 602 { 718 { 603 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 << 719 xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 604 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 720 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2)); 605 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*( << 606 } 721 } 607 else if( pLab < 1.15 ) << 722 fTotalXsc = xsection; >> 723 } >> 724 fTotalXsc *= millibarn; // parametrised in mb >> 725 // xsection *= millibarn; // parametrised in mb >> 726 >> 727 fInelasticXsc = 0.83*fTotalXsc; >> 728 fElasticXsc = fTotalXsc - fInelasticXsc; >> 729 if (fElasticXsc < 0.)fElasticXsc = 0.; >> 730 >> 731 return fTotalXsc; >> 732 } >> 733 >> 734 ///////////////////////////////////////////////////////////////////////////////////// >> 735 // >> 736 // Returns hadron-nucleon cross-section based on V. Uzjinsky parametrisation of >> 737 // data from G4FTFCrossSection class >> 738 >> 739 G4double >> 740 G4HadronNucleonXsc::GetHadronNucleonXscVU(const G4DynamicParticle* aParticle, >> 741 const G4ParticleDefinition* nucleon ) >> 742 { >> 743 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding(); >> 744 G4int absPDGcode = std::abs(PDGcode); >> 745 G4double Elab = aParticle->GetTotalEnergy(); >> 746 // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV; >> 747 G4double Plab = aParticle->GetMomentum().mag(); >> 748 // std::sqrt(Elab * Elab - 0.88); >> 749 >> 750 Elab /= GeV; >> 751 Plab /= GeV; >> 752 >> 753 G4double LogPlab = std::log( Plab ); >> 754 G4double sqrLogPlab = LogPlab * LogPlab; >> 755 >> 756 G4bool pORn = (nucleon == theProton || nucleon == theNeutron ); >> 757 G4bool proton = (nucleon == theProton); >> 758 G4bool neutron = (nucleon == theNeutron); >> 759 >> 760 >> 761 >> 762 >> 763 if( absPDGcode > 1000 && pORn ) //------Projectile is baryon - >> 764 { >> 765 if(proton) 608 { 766 { 609 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 << 767 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab; 610 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 768 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab; 611 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*( << 612 } 769 } 613 else if( pLab < 1.3 ) << 770 if(neutron) 614 { << 771 { 615 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 << 772 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab; 616 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 773 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab; 617 fElasticXsc = 3. + 13./pLab; << 618 } 774 } 619 else if( pLab < 10.) // < 3.0) // ns origi << 775 } >> 776 else if( PDGcode == 211 && pORn ) //------Projectile is PionPlus ---- >> 777 { >> 778 if(proton) 620 { 779 { 621 fTotalXsc = 36.1 + 0.079-4.313*logP+ << 780 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab; 622 3*G4Exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+ << 781 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab; 623 1.5*G4Exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12); << 624 fElasticXsc = 3. + 13./pLab; << 625 } 782 } 626 else // mb << 783 if(neutron) 627 { << 784 { 628 fElasticXsc = 3. + 13./pLab; << 785 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab; >> 786 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab; 629 } 787 } 630 } 788 } 631 else if( (theParticle == theKMinus) && proto << 789 else if( PDGcode == -211 && pORn ) //------Projectile is PionMinus ---- 632 { 790 { 633 if( pLab < pMin) << 791 if(proton) 634 { 792 { 635 G4double psp = pLab*std::sqrt(pLab); << 793 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab; 636 fElasticXsc = 5.2/psp; << 794 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab; 637 fTotalXsc = 14./psp; << 638 } 795 } 639 else if( pLab > pMax ) << 796 if(neutron) 640 { << 797 { 641 G4double ld = logP - minLogP; << 798 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab; 642 G4double ld2 = ld*ld; << 799 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab; 643 fElasticXsc = cofLogE*ld2 + 2.23; << 644 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 645 } 800 } 646 else << 801 } >> 802 else if( PDGcode == 111 && pORn ) //------Projectile is PionZero -- >> 803 { >> 804 if(proton) 647 { 805 { 648 G4double ld = logP - minLogP; << 806 fTotalXsc = (16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+ 649 G4double ld2 = ld*ld; << 807 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi- 650 G4double sp = std::sqrt(pLab); << 651 G4double psp = pLab*sp; << 652 G4double p2 = pLab*pLab; << 653 G4double p4 = p2*p2; << 654 << 655 G4double lh = pLab - 1.01; << 656 G4double hd = lh*lh + .011; << 657 808 658 G4double lm = pLab - .39; << 809 fElasticXsc = ( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+ 659 G4double md = lm*lm + .000356; << 810 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi- 660 G4double lh1 = pLab - 0.78; << 661 G4double hd1 = lh1*lh1 + .00166; << 662 G4double lh2 = pLab - 1.63; << 663 G4double hd2 = lh2*lh2 + .007; << 664 811 665 // small peaks were added but commented << 812 } 666 fElasticXsc = 5.2/psp + (1.1*cofLogE*ld << 813 if(neutron) 667 + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd; << 814 { 668 << 815 fTotalXsc = (33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+ 669 fTotalXsc = 14./psp + (1.1*cofLogT*ld2 << 816 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi- 670 + .006/md + 0.01/hd1+ 0.02/hd2 + .20/ << 817 fElasticXsc = ( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+ >> 818 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi- 671 } 819 } 672 } 820 } 673 else if( (theParticle == theKMinus) && neutr << 821 else if( PDGcode == 321 && pORn ) //------Projectile is KaonPlus -- 674 { 822 { 675 if( pLab > pMax ) << 823 if(proton) 676 { 824 { 677 G4double ld = logP - minLogP; << 825 fTotalXsc = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab; 678 G4double ld2 = ld*ld; << 826 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab; 679 fElasticXsc = cofLogE*ld2 + 2.23; << 680 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 681 } 827 } 682 else << 828 if(neutron) 683 { << 829 { 684 G4double lh = pLab - 0.98; << 830 fTotalXsc = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab; 685 G4double hd = lh*lh + .021; << 831 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab; 686 G4double sqrLogPlab = logP*logP; << 687 << 688 fElasticXsc = 5.0 + 8.1*G4Exp(-logP*1.8 << 689 - 1.3*logP + .15/hd; << 690 fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2. << 691 } 832 } 692 } 833 } 693 else if( (theParticle == theKPlus) && proton << 834 else if( PDGcode ==-321 && pORn ) //------Projectile is KaonMinus ---- 694 { 835 { 695 // VI: modified low-energy part << 836 if(proton) 696 if( pLab < 0.631 ) << 697 { 837 { 698 fElasticXsc = fTotalXsc = 12.03; << 838 fTotalXsc = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66*sqrLogPlab - 5.6*LogPlab; >> 839 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29*sqrLogPlab - 2.4*LogPlab; 699 } 840 } 700 else if( pLab > pMax ) << 841 if(neutron) 701 { << 842 { 702 G4double ld = logP - minLogP; << 843 fTotalXsc = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab; 703 G4double ld2 = ld*ld; << 844 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab; 704 fElasticXsc = cofLogE*ld2 + 2.23; << 705 fTotalXsc = cofLogT*ld2 + 19.2; << 706 } 845 } 707 else << 846 } >> 847 else if( PDGcode == 311 && pORn ) //------Projectile is KaonZero ----- >> 848 { >> 849 if(proton) 708 { 850 { 709 G4double ld = logP - minLogP; << 851 fTotalXsc = ( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+ 710 G4double ld2 = ld*ld; << 852 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K- 711 G4double lr = pLab - .38; << 853 fElasticXsc = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+ 712 G4double LE = .7/(lr*lr + .076); << 854 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K- 713 G4double sp = std::sqrt(pLab); << 855 } 714 G4double p2 = pLab*pLab; << 856 if(neutron) 715 G4double p4 = p2*p2; << 857 { 716 // VI: tuned elastic << 858 fTotalXsc = ( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+ 717 fElasticXsc = LE + (cofLogE*ld2 + 2.23) << 859 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K- 718 + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652); << 860 fElasticXsc = ( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+ 719 fTotalXsc = LE + (cofLogT*ld2 + 19.5) << 861 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K- 720 + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392); << 721 } 862 } 722 } 863 } 723 else if( (theParticle == theKPlus) && neutr << 864 else //------Projectile is undefined, Nucleon assumed 724 { 865 { 725 if( pLab < pMin ) << 866 if(proton) 726 { 867 { 727 G4double lm = pLab - 0.94; << 868 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab; 728 G4double md = lm*lm + .392; << 869 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab; 729 fElasticXsc = 2./md; << 730 fTotalXsc = 4.6/md; << 731 } 870 } 732 else if( pLab > pMax ) << 871 if(neutron) 733 { << 872 { 734 G4double ld = logP - minLogP; << 873 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab; 735 G4double ld2 = ld*ld; << 874 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab; 736 fElasticXsc = cofLogE*ld2 + 2.23; << 737 fTotalXsc = cofLogT*ld2 + 19.2; << 738 } 875 } 739 else << 876 } 740 { << 877 fTotalXsc *= millibarn; 741 G4double ld = logP - minLogP; << 878 fElasticXsc *= millibarn; 742 G4double ld2 = ld*ld; << 879 fInelasticXsc = fTotalXsc - fElasticXsc; 743 G4double sp = std::sqrt(pLab); << 880 if (fInelasticXsc < 0.) fInelasticXsc = 0.; 744 G4double p2 = pLab*pLab; << 881 745 G4double p4 = p2*p2; << 882 return fTotalXsc; 746 G4double lm = pLab - 0.94; << 747 G4double md = lm*lm + .392; << 748 fElasticXsc = (cofLogE*ld2 + 2.23)/(1. << 749 fTotalXsc = (cofLogT*ld2 + 19.5)/(1. << 750 } << 751 } << 752 fTotalXsc *= CLHEP::millibarn; << 753 fElasticXsc *= CLHEP::millibarn; << 754 fElasticXsc = std::min(fElasticXsc, fTotalX << 755 << 756 if( proton && theParticle->GetPDGCharge() > << 757 { << 758 G4double cB = G4NuclearRadii::CoulombFacto << 759 fTotalXsc *= cB; << 760 fElasticXsc *= cB; << 761 } << 762 fInelasticXsc = std::max(fTotalXsc - fElasti << 763 /* << 764 G4cout<< "HNXsc: Ekin(GeV)= " << ekin/GeV << << 765 <<"; el(mb)= " <<fElasticXsc/millibarn << 766 <<"; inel(mb)= " <<fInelasticXsc/milli << 767 << theParticle->GetParticleName() << " << 768 << nucleon->GetParticleName() << G4end << 769 */ << 770 return fTotalXsc; << 771 } 883 } 772 884 773 ////////////////////////////////////////////// << 885 ///////////////////////////////////////////////////////////////////////////////////// 774 // 886 // 775 // Returns kaon-nucleon cross-section based on << 887 // Returns hadron-nucleon cross-section based on Mikhail Kossov parametrisation of 776 // tuned for the Glauber-Gribov hadron model f << 888 // data from G4QuasiFreeRatios class 777 889 778 G4double G4HadronNucleonXsc::KaonNucleonXscGG( << 890 G4double 779 const G4ParticleD << 891 G4HadronNucleonXsc::GetHadronNucleonXscMK(const G4DynamicParticle* aParticle, 780 const G4ParticleDefinition* nucleon << 892 const G4ParticleDefinition* nucleon ) 781 { 893 { 782 fTotalXsc = fElasticXsc = fInelasticXsc = 0. << 894 G4int I = -1; 783 if(theParticle == theKMinus || theParticle = << 895 G4int PDG = aParticle->GetDefinition()->GetPDGEncoding(); 784 KaonNucleonXscVG(theParticle, nucleon, eki << 896 785 << 897 // G4int absPDG = std::abs(PDG); 786 } else if(theParticle == theK0S || thePartic << 898 787 G4double stot = KaonNucleonXscVG(theKMinu << 899 G4double p = aParticle->GetMomentum().mag()/GeV; 788 G4double sel = fElasticXsc; << 900 789 G4double sinel = fInelasticXsc; << 901 G4bool F = false; 790 stot += KaonNucleonXscVG(theKPlus, nucleo << 902 if(nucleon == theProton) F = true; 791 sel += fElasticXsc; << 903 else if(nucleon == theNeutron) F = false; 792 sinel += fInelasticXsc; << 904 else 793 fTotalXsc = stot*0.5; << 905 { 794 fElasticXsc = sel*0.5; << 906 G4cout << "nucleon is not proton or neutron, return xsc for proton" << G4endl; 795 fInelasticXsc = sinel*0.5; << 907 F = true; 796 } 908 } 797 return fTotalXsc; << 798 } << 799 909 800 ////////////////////////////////////////////// << 910 G4bool kfl = true; // Flag of K0/aK0 oscillation 801 // << 911 G4bool kf = false; 802 // Returns kaon-nucleon cross-section using NS << 803 912 804 G4double G4HadronNucleonXsc::KaonNucleonXscNS( << 913 if( PDG == 130 || PDG == 310 ) 805 const G4ParticleD << 914 { 806 const G4ParticleDefinition* nucleon << 915 kf = true; 807 { << 916 if( G4UniformRand() > .5 ) kfl = false; 808 fTotalXsc = fElasticXsc = fInelasticXsc = 0. << 809 if(theParticle == theKMinus || theParticle = << 810 HadronNucleonXscNS(theParticle, nucleon, e << 811 << 812 } else if(theParticle == theK0S || thePartic << 813 G4double fact = 0.5; << 814 G4double stot = 0.0; << 815 G4double sel = 0.0; << 816 G4double sinel= 0.0; << 817 if(ekin > ekinmaxQB) { << 818 stot = HadronNucleonXscNS(theKMinus, nu << 819 sel = fElasticXsc; << 820 sinel = fInelasticXsc; << 821 stot += HadronNucleonXscNS(theKPlus, nu << 822 sel += fElasticXsc; << 823 sinel += fInelasticXsc; << 824 } else { << 825 fact *= std::sqrt(ekinmaxQB/std::max(eki << 826 stot = HadronNucleonXscNS(theKMinus, nu << 827 sel = fElasticXsc; << 828 sinel = fInelasticXsc; << 829 stot += HadronNucleonXscNS(theKPlus, nu << 830 sel += fElasticXsc; << 831 sinel += fInelasticXsc; << 832 } << 833 fTotalXsc = stot*fact; << 834 fElasticXsc = sel*fact; << 835 fInelasticXsc = sinel*fact; << 836 } 917 } 837 return fTotalXsc; << 918 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) I = 0; // pp/nn 838 } << 919 else if( (PDG == 2112 && F) || (PDG == 2212 && !F) ) I = 1; // np/pn 839 920 840 ////////////////////////////////////////////// << 921 else if( (PDG == -211 && F) || (PDG == 211 && !F) ) I = 2; // pimp/pipn 841 // << 922 else if( (PDG == 211 && F) || (PDG ==-211 && !F) ) I = 3; // pipp/pimn 842 // Returns kaon-nucleon cross-section with smo << 843 923 844 G4double G4HadronNucleonXsc::KaonNucleonXscVG( << 924 else if( PDG == -321 || PDG == -311 || ( kf && !kfl ) ) I = 4; // KmN/K0N 845 const G4ParticleDefinition* t << 925 else if( PDG == 321 || PDG == 311 || ( kf && kfl ) ) I = 5; // KpN/aK0N 846 const G4ParticleDefinition* nucleon, G4do << 847 { << 848 G4double pM = theParticle->GetPDGMass(); << 849 G4double pLab = std::sqrt(ekin*(ekin + 2*pM) << 850 926 851 pLab *= invGeV; << 927 else if( PDG > 3000 && PDG < 3335) I = 6; // @@ for all hyperons - take Lambda 852 G4double logP = G4Log(pLab); << 928 else if( PDG < -2000 && PDG > -3335) I = 7; // @@ for all anti-baryons - anti-p/anti-n >> 929 else >> 930 { >> 931 G4cout<<"MK PDG = "<<PDG >> 932 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; >> 933 G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash"); >> 934 } 853 935 854 fTotalXsc = 0.0; << 936 // Each parameter set can have not more than nPoints = 128 parameters 855 937 856 G4bool proton = (nucleon == theProton); << 938 static const G4double lmi = 3.5; // min of (lnP-lmi)^2 parabola 857 G4bool neutron = (nucleon == theNeutron); << 939 static const G4double pbe = .0557; // elastic (lnP-lmi)^2 parabola coefficient >> 940 static const G4double pbt = .3; // total (lnP-lmi)^2 parabola coefficient >> 941 static const G4double pmi = .1; // Below that fast LE calculation is made >> 942 static const G4double pma = 1000.; // Above that fast HE calculation is made >> 943 >> 944 if( p <= 0.) >> 945 { >> 946 G4cout<<" p = "<<p<<" is zero or negative"<<G4endl; >> 947 >> 948 fElasticXsc = 0.; >> 949 fInelasticXsc = 0.; >> 950 fTotalXsc = 0.; 858 951 859 if( (theParticle == theKMinus) && proton ) << 952 return fTotalXsc; >> 953 } >> 954 if (!I) // pp/nn 860 { 955 { 861 if( pLab < pMin) << 956 if( p < pmi ) 862 { 957 { 863 G4double psp = pLab*std::sqrt(pLab); << 958 G4double p2 = p*p; 864 fElasticXsc = 5.2/psp; << 959 fElasticXsc = 1./(.00012 + p2*.2); 865 fTotalXsc = 14./psp; << 960 fTotalXsc = fElasticXsc; 866 } 961 } 867 else if( pLab > pMax ) << 962 else if(p>pma) 868 { 963 { 869 G4double ld = logP - minLogP; << 964 G4double lp = std::log(p)-lmi; 870 G4double ld2 = ld*ld; << 965 G4double lp2 = lp*lp; 871 fElasticXsc = cofLogE*ld2 + 2.23; << 966 fElasticXsc = pbe*lp2 + 6.72; 872 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 967 fTotalXsc = pbt*lp2 + 38.2; 873 } 968 } 874 else 969 else 875 { 970 { 876 G4double ld = logP - minLogP; << 971 G4double p2 = p*p; 877 G4double ld2 = ld*ld; << 972 G4double LE = 1./( .00012 + p2*.2); 878 G4double sp = std::sqrt(pLab); << 973 G4double lp = std::log(p) - lmi; 879 G4double psp = pLab*sp; << 974 G4double lp2 = lp*lp; 880 G4double p2 = pLab*pLab; << 975 G4double rp2 = 1./p2; 881 G4double p4 = p2*p2; << 976 fElasticXsc = LE + ( pbe*lp2 + 6.72+32.6/p)/( 1. + rp2/p); 882 << 977 fTotalXsc = LE + ( pbt*lp2 + 38.2+52.7*rp2)/( 1. + 2.72*rp2*rp2); 883 G4double lh = pLab - 1.01; << 884 G4double hd = lh*lh + .011; << 885 fElasticXsc = 5.2/psp + (cofLogE*ld2 + << 886 fTotalXsc = 14./psp + (1.1*cofLogT*ld << 887 } 978 } 888 } 979 } 889 else if( (theParticle == theKMinus) && neutr << 980 else if( I==1 ) // np/pn 890 { 981 { 891 if( pLab > pMax ) << 982 if( p < pmi ) 892 { 983 { 893 G4double ld = logP - minLogP; << 984 G4double p2 = p*p; 894 G4double ld2 = ld*ld; << 985 fElasticXsc = 1./( .00012 + p2*( .051 + .1*p2)); 895 fElasticXsc = cofLogE*ld2 + 2.23; << 986 fTotalXsc = fElasticXsc; 896 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 987 } >> 988 else if( p > pma ) >> 989 { >> 990 G4double lp = std::log(p) - lmi; >> 991 G4double lp2 = lp*lp; >> 992 fElasticXsc = pbe*lp2 + 6.72; >> 993 fTotalXsc = pbt*lp2 + 38.2; 897 } 994 } 898 else 995 else 899 { << 996 { 900 G4double lh = pLab - 0.98; << 997 G4double p2 = p*p; 901 G4double hd = lh*lh + .045; // vg ve << 998 G4double LE = 1./( .00012 + p2*( .051 + .1*p2 ) ); 902 G4double sqrLogPlab = logP*logP; << 999 G4double lp = std::log(p) - lmi; 903 << 1000 G4double lp2 = lp*lp; 904 fElasticXsc = 5.0 + 8.1*G4Exp(-logP*1.8 << 1001 G4double rp2 = 1./p2; 905 - 1.3*logP + .15/hd; << 1002 fElasticXsc = LE + (pbe*lp2 + 6.72 + 30./p)/( 1. + .49*rp2/p); 906 fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2. << 1003 fTotalXsc = LE + (pbt*lp2 + 38.2)/( 1. + .54*rp2*rp2); 907 } 1004 } 908 } 1005 } 909 else if( (theParticle == theKPlus) && proton << 1006 else if( I == 2 ) // pimp/pipn 910 { 1007 { 911 // VI: modified low-energy part << 1008 G4double lp = std::log(p); 912 if( pLab < 0.631 ) << 1009 >> 1010 if(p<pmi) 913 { 1011 { 914 fElasticXsc = fTotalXsc = 12.03; << 1012 G4double lr = lp + 1.27; >> 1013 fElasticXsc = 1.53/( lr*lr + .0676); >> 1014 fTotalXsc = fElasticXsc*3; 915 } 1015 } 916 else if( pLab > pMax ) << 1016 else if( p > pma ) 917 { 1017 { 918 G4double ld = logP - minLogP; << 1018 G4double ld = lp - lmi; 919 G4double ld2 = ld*ld; 1019 G4double ld2 = ld*ld; 920 fElasticXsc = cofLogE*ld2 + 2.23; << 1020 G4double sp = std::sqrt(p); 921 fTotalXsc = cofLogT*ld2 + 19.2; << 1021 fElasticXsc = pbe*ld2 + 2.4 + 7./sp; >> 1022 fTotalXsc = pbt*ld2 + 22.3 + 12./sp; 922 } 1023 } 923 else 1024 else 924 { 1025 { 925 G4double ld = logP - minLogP; << 1026 G4double lr = lp + 1.27; >> 1027 G4double LE = 1.53/( lr*lr + .0676); >> 1028 G4double ld = lp - lmi; 926 G4double ld2 = ld*ld; 1029 G4double ld2 = ld*ld; 927 G4double lr = pLab - .38; << 1030 G4double p2 = p*p; 928 G4double LE = .7/(lr*lr + .076); << 929 G4double sp = std::sqrt(pLab); << 930 G4double p2 = pLab*pLab; << 931 G4double p4 = p2*p2; 1031 G4double p4 = p2*p2; 932 // VI: tuned elastic << 1032 G4double sp = std::sqrt(p); 933 fElasticXsc = LE + (cofLogE*ld2 + 2.23) << 1033 G4double lm = lp + .36; 934 + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652); << 1034 G4double md = lm*lm + .04; 935 fTotalXsc = LE + (cofLogT*ld2 + 19.5) << 1035 G4double lh = lp - .017; 936 + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392); << 1036 G4double hd = lh*lh + .0025; >> 1037 fElasticXsc = LE + (pbe*ld2 + 2.4 + 7./sp)/( 1. + .7/p4) + .6/md + .05/hd; >> 1038 fTotalXsc = LE*3 + (pbt*ld2 + 22.3 + 12./sp)/(1. + .4/p4) + 1./md + .06/hd; 937 } 1039 } 938 } 1040 } 939 else if( (theParticle == theKPlus) && neutr << 1041 else if( I == 3 ) // pipp/pimn 940 { 1042 { 941 if( pLab < pMin ) << 1043 G4double lp = std::log(p); >> 1044 >> 1045 if( p < pmi ) 942 { 1046 { 943 G4double lm = pLab - 0.94; << 1047 G4double lr = lp + 1.27; 944 G4double md = lm*lm + .392; << 1048 G4double lr2 = lr*lr; 945 fElasticXsc = 2./md; << 1049 fElasticXsc = 13./( lr2 + lr2*lr2 + .0676); 946 fTotalXsc = 4.6/md; << 1050 fTotalXsc = fElasticXsc; 947 } 1051 } 948 else if( pLab > pMax ) << 1052 else if( p > pma ) 949 { 1053 { 950 G4double ld = logP - minLogP; << 1054 G4double ld = lp - lmi; 951 G4double ld2 = ld*ld; 1055 G4double ld2 = ld*ld; 952 fElasticXsc = cofLogE*ld2 + 2.23; << 1056 G4double sp = std::sqrt(p); 953 fTotalXsc = cofLogT*ld2 + 19.2; << 1057 fElasticXsc = pbe*ld2 + 2.4 + 6./sp; >> 1058 fTotalXsc = pbt*ld2 + 22.3 + 5./sp; 954 } 1059 } 955 else 1060 else 956 { 1061 { 957 G4double ld = logP - minLogP; << 1062 G4double lr = lp + 1.27; >> 1063 G4double lr2 = lr*lr; >> 1064 G4double LE = 13./(lr2 + lr2*lr2 + .0676); >> 1065 G4double ld = lp - lmi; 958 G4double ld2 = ld*ld; 1066 G4double ld2 = ld*ld; 959 G4double sp = std::sqrt(pLab); << 1067 G4double p2 = p*p; 960 G4double p2 = pLab*pLab; << 961 G4double p4 = p2*p2; 1068 G4double p4 = p2*p2; 962 G4double lm = pLab - 0.94; << 1069 G4double sp = std::sqrt(p); 963 G4double md = lm*lm + .392; << 1070 G4double lm = lp - .32; 964 fElasticXsc = (cofLogE*ld2 + 2.23)/(1. << 1071 G4double md = lm*lm + .0576; 965 fTotalXsc = (cofLogT*ld2 + 19.5)/(1. << 1072 fElasticXsc = LE + (pbe*ld2 + 2.4 + 6./sp)/(1. + 3./p4) + .7/md; >> 1073 fTotalXsc = LE + (pbt*ld2 + 22.3 + 5./sp)/(1. + 1./p4) + .8/md; 966 } 1074 } 967 } 1075 } 968 << 1076 else if( I == 4 ) // Kmp/Kmn/K0p/K0n 969 fTotalXsc *= CLHEP::millibarn; << 970 fElasticXsc *= CLHEP::millibarn; << 971 << 972 if( proton && theParticle->GetPDGCharge() > << 973 { << 974 G4double cB = G4NuclearRadii::CoulombFacto << 975 fTotalXsc *= cB; << 976 fElasticXsc *= cB; << 977 } << 978 fElasticXsc = std::min(fElasticXsc, fTotalXs << 979 fInelasticXsc = std::max(fTotalXsc - fElasti << 980 /* << 981 G4cout << "HNXscVG: E= " << ekin << " " << t << 982 << " P: " << proton << " xtot(b)= " << fTot << 983 << " xel(b)= " << fElasticXsc/barn << " xi << 984 << G4endl; << 985 */ << 986 return fTotalXsc; << 987 } << 988 << 989 ////////////////////////////////////////////// << 990 // << 991 // Returns hyperon-nucleon cross-section using << 992 << 993 G4double G4HadronNucleonXsc::HyperonNucleonXsc << 994 const G4ParticleDefinition* thePartic << 995 const G4ParticleDefinition* nucleon, G4doub << 996 { << 997 G4double coeff = 1.0; << 998 G4int pdg = std::abs(theParticle->GetPDGEnco << 999 << 1000 // lambda, sigma+-0 and anti-hyperons << 1001 if( pdg == 3122 || pdg == 3112 || pdg == 32 << 1002 { << 1003 coeff = 0.88; << 1004 } << 1005 // Xi hyperons and anti-hyperons << 1006 else if( pdg == 3312 || pdg == 3322 ) << 1007 { << 1008 coeff = 0.76; << 1009 } << 1010 // omega, anti_omega << 1011 else if( pdg == 3334 ) << 1012 { << 1013 coeff = 0.64; << 1014 } << 1015 // lambdaC, sigmaC+-0 and anti-hyperonsC << 1016 else if( pdg == 4122 || pdg == 4112 || pdg << 1017 { << 1018 coeff = 0.784378; << 1019 } << 1020 // omegaC0, anti_omegaC0 << 1021 else if( pdg == 4332 ) << 1022 { << 1023 coeff = 0.544378; << 1024 } << 1025 // XiC+0 and anti-hyperonC << 1026 else if( pdg == 4132 || pdg == 4232 ) << 1027 { << 1028 coeff = 0.664378; << 1029 } << 1030 // lambdaB, sigmaB+-0 and anti-hyperonsB << 1031 else if( pdg == 5122 || pdg == 5112 || pdg << 1032 { << 1033 coeff = 0.740659; << 1034 } << 1035 // omegaB0, anti_omegaB0 << 1036 else if( pdg == 5332 ) << 1037 { << 1038 coeff = 0.500659; << 1039 } << 1040 // XiB+0 and anti-hyperonB << 1041 else if( pdg == 5132 || pdg == 5232 ) << 1042 { << 1043 coeff = 0.620659; << 1044 } << 1045 fTotalXsc = coeff*HadronNucleonXscNS( thePr << 1046 fInelasticXsc *= coeff; << 1047 fElasticXsc *= coeff; << 1048 << 1049 return fTotalXsc; << 1050 } << 1051 << 1052 ///////////////////////////////////////////// << 1053 // << 1054 // Returns hyperon-nucleon cross-section usin << 1055 << 1056 G4double G4HadronNucleonXsc::SCBMesonNucleonX << 1057 const G4ParticleDefinition* theParti << 1058 const G4ParticleDefinition* nucleon, << 1059 { << 1060 G4double coeff(1.0); << 1061 G4int pdg = std::abs(theParticle->GetPDGEnc << 1062 << 1063 // B+-0 anti << 1064 if( pdg == 511 || pdg == 521 ) << 1065 { << 1066 coeff = 0.610989; << 1067 } << 1068 // D+-0 anti << 1069 else if( pdg == 411 || pdg == 421 ) << 1070 { << 1071 coeff = 0.676568; << 1072 } << 1073 // Bs, antiBs << 1074 else if( pdg == 531 ) << 1075 { << 1076 coeff = 0.430989; << 1077 } << 1078 // Bc+- << 1079 else if( pdg == 541 ) << 1080 { << 1081 coeff = 0.287557; << 1082 } << 1083 // Ds+- << 1084 else if( pdg == 431 ) << 1085 { << 1086 coeff = 0.496568; << 1087 } << 1088 // etaC, J/Psi << 1089 else if( pdg == 441 || pdg == 443 ) << 1090 { << 1091 coeff = 0.353135; << 1092 } << 1093 // Upsilon << 1094 else if( pdg == 553 ) << 1095 { << 1096 coeff = 0.221978; << 1097 } << 1098 // eta << 1099 else if( pdg == 221 ) << 1100 { 1077 { 1101 coeff = 0.76; << 1078 if( p < pmi) 1102 } << 1103 // eta' << 1104 else if( pdg == 331 ) << 1105 { << 1106 coeff = 0.88; << 1107 } << 1108 fTotalXsc = coeff*HadronNucleonXscNS(thePiP << 1109 fElasticXsc *= coeff; << 1110 fInelasticXsc *= coeff; << 1111 return fTotalXsc; << 1112 } << 1113 ///////////////////////////////////////////// << 1114 // << 1115 // Returns hadron-nucleon cross-section based << 1116 // data from G4FTFCrossSection class << 1117 << 1118 G4double G4HadronNucleonXsc::HadronNucleonXsc << 1119 const G4Particle << 1120 const G4ParticleDefinition* nucleo << 1121 { << 1122 G4int PDGcode = theParticle->GetPDGEncoding << 1123 G4int absPDGcode = std::abs(PDGcode); << 1124 G4double mass = theParticle->GetPDGMass(); << 1125 G4double Plab = std::sqrt(ekin*(ekin + 2.*m << 1126 << 1127 G4double logPlab = G4Log( Plab ); << 1128 G4double sqrLogPlab = logPlab * logPlab; << 1129 << 1130 G4bool proton = (nucleon == theProton); << 1131 G4bool neutron = (nucleon == theNeutron); << 1132 << 1133 if( absPDGcode > 1000) //------Projectile << 1134 { << 1135 if(proton) << 1136 { 1079 { 1137 fTotalXsc = 48.0 + 0.522*sqrLogPlab - << 1080 G4double psp = p*std::sqrt(p); 1138 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1081 fElasticXsc = 5.2/psp; 1139 } << 1082 fTotalXsc = 14./psp; 1140 if(neutron) << 1141 { << 1142 fTotalXsc = 47.3 + 0.513*sqrLogPlab << 1143 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1144 } 1083 } 1145 } << 1084 else if( p > pma ) 1146 else if( PDGcode == 211) //------Projecti << 1147 { << 1148 if(proton) << 1149 { 1085 { 1150 fTotalXsc = 16.4 + 19.3 *G4Exp(-logPla << 1086 G4double ld = std::log(p) - lmi; 1151 fElasticXsc = 0.0 + 11.4*G4Exp(-logPla << 1087 G4double ld2 = ld*ld; >> 1088 fElasticXsc = pbe*ld2 + 2.23; >> 1089 fTotalXsc = pbt*ld2 + 19.5; 1152 } 1090 } 1153 if(neutron) << 1091 else 1154 { << 1092 { 1155 fTotalXsc = 33.0 + 14.0 *G4Exp(-logP << 1093 G4double ld = std::log(p) - lmi; 1156 fElasticXsc = 1.76 + 11.2*G4Exp(-logPla << 1094 G4double ld2 = ld*ld; >> 1095 G4double sp = std::sqrt(p); >> 1096 G4double psp = p*sp; >> 1097 G4double p2 = p*p; >> 1098 G4double p4 = p2*p2; >> 1099 G4double lm = p - .39; >> 1100 G4double md = lm*lm + .000156; >> 1101 G4double lh = p - 1.; >> 1102 G4double hd = lh*lh + .0156; >> 1103 fElasticXsc = 5.2/psp + (pbe*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md + .15/hd; >> 1104 fTotalXsc = 14./psp + (pbt*ld2 + 19.5)/(1. - .21/sp + .52/p4) + .006/md + .30/hd; 1157 } 1105 } 1158 } 1106 } 1159 else if( PDGcode == -211) //------Projecti << 1107 else if( I == 5 ) // Kpp/Kpn/aKp/aKn 1160 { 1108 { 1161 if(proton) << 1109 if( p < pmi ) 1162 { 1110 { 1163 fTotalXsc = 33.0 + 14.0 *G4Exp(-logPl << 1111 G4double lr = p - .38; 1164 fElasticXsc = 1.76 + 11.2*G4Exp(-logPla << 1112 G4double lm = p - 1.; 1165 } << 1113 G4double md = lm*lm + .372; 1166 if(neutron) << 1114 fElasticXsc = .7/(lr*lr + .0676) + 2./md; 1167 { << 1115 fTotalXsc = fElasticXsc + .6/md; 1168 fTotalXsc = 16.4 + 19.3 *G4Exp(-logPl << 1169 fElasticXsc = 0.0 + 11.4*G4Exp(-logPla << 1170 } 1116 } 1171 } << 1117 else if( p > pma ) 1172 else if( PDGcode == 111) //------Projecti << 1173 { << 1174 if(proton) << 1175 { 1118 { 1176 fTotalXsc = (16.4 + 19.3 *G4Exp(-logPla << 1119 G4double ld = std::log(p) - lmi; 1177 33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0. << 1120 G4double ld2 = ld*ld; 1178 << 1121 fElasticXsc = pbe*ld2 + 2.23; 1179 fElasticXsc = (0.0 + 11.4*G4Exp(-logPla << 1122 fTotalXsc = pbt*ld2 + 19.5; 1180 1.76 + 11.2*G4Exp(-logPlab*0.64) + 0 << 1181 << 1182 } 1123 } 1183 if(neutron) << 1124 else 1184 { << 1125 { 1185 fTotalXsc = (33.0 + 14.0 *G4Exp(-logPla << 1126 G4double ld = std::log(p) - lmi; 1186 16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0. << 1127 G4double ld2 = ld*ld; 1187 fElasticXsc = (1.76 + 11.2*G4Exp(-logPl << 1128 G4double lr = p - .38; 1188 0.0 + 11.4*G4Exp(-logPlab*0.40) + 0 << 1129 G4double LE = .7/(lr*lr + .0676); >> 1130 G4double sp = std::sqrt(p); >> 1131 G4double p2 = p*p; >> 1132 G4double p4 = p2*p2; >> 1133 G4double lm = p - 1.; >> 1134 G4double md = lm*lm + .372; >> 1135 fElasticXsc = LE + (pbe*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md; >> 1136 fTotalXsc = LE + (pbt*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 2.6/md; 1189 } 1137 } 1190 } 1138 } 1191 else if( PDGcode == 321 ) //------Projec << 1139 else if( I == 6 ) // hyperon-N 1192 { 1140 { 1193 if(proton) << 1141 if( p < pmi ) 1194 { 1142 { 1195 fTotalXsc = 18.1 + 0.26 *sqrLogPlab << 1143 G4double p2 = p*p; 1196 fElasticXsc = 5.0 + 8.1*G4Exp(-logPla << 1144 fElasticXsc = 1./(.002 + p2*(.12 + p2)); 1197 } << 1145 fTotalXsc = fElasticXsc; 1198 if(neutron) << 1199 { << 1200 fTotalXsc = 18.7 + 0.21 *sqrLogPlab << 1201 fElasticXsc = 7.3 + 0.29 *sqrLogPlab << 1202 } 1146 } 1203 } << 1147 else if( p > pma ) 1204 else if( PDGcode ==-321 ) //------Projecti << 1205 { << 1206 if(proton) << 1207 { 1148 { 1208 fTotalXsc = 32.1 + 0.66*sqrLogPlab - << 1149 G4double lp = std::log(p) - lmi; 1209 fElasticXsc = 7.3 + 0.29*sqrLogPlab - << 1150 G4double lp2 = lp*lp; >> 1151 G4double sp = std::sqrt(p); >> 1152 fElasticXsc = (pbe*lp2 + 6.72)/(1. + 2./sp); >> 1153 fTotalXsc = (pbt*lp2 + 38.2 + 900./sp)/(1. + 27./sp); 1210 } 1154 } 1211 if(neutron) << 1155 else 1212 { << 1156 { 1213 fTotalXsc = 25.2 + 0.38*sqrLogPlab - << 1157 G4double p2 = p*p; 1214 fElasticXsc = 5.0 + 8.1*G4Exp(-logPla << 1158 G4double LE = 1./(.002 + p2*(.12 + p2)); >> 1159 G4double lp = std::log(p) - lmi; >> 1160 G4double lp2 = lp*lp; >> 1161 G4double p4 = p2*p2; >> 1162 G4double sp = std::sqrt(p); >> 1163 fElasticXsc = LE + (pbe*lp2 + 6.72 + 99./p2)/(1. + 2./sp + 2./p4); >> 1164 fTotalXsc = LE + (pbt*lp2 + 38.2 + 900./sp)/(1. + 27./sp + 3./p4); 1215 } 1165 } 1216 } 1166 } 1217 else if( PDGcode == 311 ) //------Projecti << 1167 else if( I == 7 ) // antibaryon-N 1218 { 1168 { 1219 if(proton) << 1169 if( p > pma ) 1220 { 1170 { 1221 fTotalXsc = ( 18.1 + 0.26 *sqrLogPlab << 1171 G4double lp = std::log(p) - lmi; 1222 32.1 + 0.66 *sqrLogPlab << 1172 G4double lp2 = lp*lp; 1223 fElasticXsc = ( 5.0 + 8.1*G4Exp(-logP << 1173 fElasticXsc = pbe*lp2 + 6.72; 1224 7.3 + 0.29 *sqrLogPl << 1174 fTotalXsc = pbt*lp2 + 38.2; 1225 } 1175 } 1226 if(neutron) << 1176 else 1227 { << 1177 { 1228 fTotalXsc = ( 18.7 + 0.21 *sqrLogPlab << 1178 G4double ye = std::pow(p, 1.25); 1229 25.2 + 0.38 *sqrLogPlab << 1179 G4double yt = std::pow(p, .35); 1230 fElasticXsc = ( 7.3 + 0.29 *sqrLogPlab << 1180 G4double lp = std::log(p) - lmi; 1231 5.0 + 8.1*G4Exp(-logPl << 1181 G4double lp2 = lp*lp; >> 1182 fElasticXsc = 80./(ye + 1.) + pbe*lp2 + 6.72; >> 1183 fTotalXsc = (80./yt + .3)/yt +pbt*lp2 + 38.2; 1232 } 1184 } 1233 } 1185 } 1234 else //------Projectile is undefined, Nucl << 1186 else 1235 { 1187 { 1236 if(proton) << 1188 G4cout<<"PDG incoding = "<<I<<" is not defined (0-7)"<<G4endl; 1237 { << 1189 1238 fTotalXsc = 48.0 + 0.522*sqrLogPlab - << 1239 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1240 } << 1241 if(neutron) << 1242 { << 1243 fTotalXsc = 47.3 + 0.513*sqrLogPlab - << 1244 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1245 } << 1246 } 1190 } >> 1191 if( fElasticXsc > fTotalXsc ) fElasticXsc = fTotalXsc; 1247 1192 1248 fTotalXsc *= CLHEP::millibarn; << 1193 fTotalXsc *= millibarn; 1249 fElasticXsc *= CLHEP::millibarn; << 1194 fElasticXsc *= millibarn; 1250 fElasticXsc = std::min(fElasticXsc, fTotalX << 1195 fInelasticXsc = fTotalXsc - fElasticXsc; 1251 fInelasticXsc = fTotalXsc - fElasticXsc; << 1196 if (fInelasticXsc < 0.) fInelasticXsc = 0.; 1252 1197 1253 return fTotalXsc; << 1254 } << 1255 << 1256 ///////////////////////////////////////////// << 1257 // << 1258 // Returns hadron-nucleon Xsc according to di << 1259 // [2] E. Levin, hep-ph/9710546 << 1260 // [3] U. Dersch, et al, hep-ex/9910052 << 1261 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 ( << 1262 << 1263 G4double G4HadronNucleonXsc::HadronNucleonXsc << 1264 const G4Particle << 1265 const G4ParticleDefinition*, G4dou << 1266 { << 1267 G4int pdg = theParticle->GetPDGEncoding(); << 1268 G4double xsection(0.); << 1269 static const G4double targ_mass = << 1270 0.5*(CLHEP::proton_mass_c2 + CLHEP::neutr << 1271 G4double sMand = << 1272 CalcMandelstamS(ekin, theParticle->GetPDG << 1273 << 1274 G4double x1 = G4Exp(G4Log(sMand)*0.0808); << 1275 G4double x2 = G4Exp(G4Log(-sMand)*0.4525); << 1276 1198 1277 if(pdg == 22) << 1278 { << 1279 xsection = 0.0677*x1 + 0.129*x2; << 1280 } << 1281 else if(theParticle == theNeutron) << 1282 { << 1283 xsection = 21.70*x1 + 56.08*x2; << 1284 } << 1285 else if(theParticle == theProton) << 1286 { << 1287 xsection = 21.70*x1 + 56.08*x2; << 1288 } << 1289 // pbar << 1290 else if(pdg == -2212) << 1291 { << 1292 xsection = 21.70*x1 + 98.39*x2; << 1293 } << 1294 else if(theParticle == thePiPlus) << 1295 { << 1296 xsection = 13.63*x1 + 27.56*x2; << 1297 } << 1298 // pi- << 1299 else if(pdg == -211) << 1300 { << 1301 xsection = 13.63*x1 + 36.02*x2; << 1302 } << 1303 else if(theParticle == theKPlus) << 1304 { << 1305 xsection = 11.82*x1 + 8.15*x2; << 1306 } << 1307 else if(theParticle == theKMinus) << 1308 { << 1309 xsection = 11.82*x1 + 26.36*x2; << 1310 } << 1311 else if(theParticle == theK0S || theParticl << 1312 { << 1313 xsection = 11.82*x1 + 17.25*x2; << 1314 } << 1315 else << 1316 { << 1317 xsection = 21.70*x1 + 56.08*x2; << 1318 } << 1319 fTotalXsc = xsection*CLHEP::millibarn; << 1320 fInelasticXsc = 0.83*fTotalXsc; << 1321 fElasticXsc = fTotalXsc - fInelasticXsc; << 1322 return fTotalXsc; 1199 return fTotalXsc; 1323 } 1200 } 1324 1201 1325 ///////////////////////////////////////////// << 1326 << 1327 G4double << 1328 G4HadronNucleonXsc::CoulombBarrier(const G4Pa << 1329 const G4ParticleDefinition* nucleo << 1330 G4double ekin) << 1331 { << 1332 G4double tR = 0.895*CLHEP::fermi; << 1333 G4double pR = 0.5*CLHEP::fermi; << 1334 1202 1335 if ( theParticle == theProton ) pR = 0. << 1336 else if( theParticle == thePiPlus ) pR = 0. << 1337 else if( theParticle == theKPlus ) pR = 0. << 1338 1203 1339 G4double pZ = theParticle->GetPDGCharge(); << 1340 G4double tZ = nucleon->GetPDGCharge(); << 1341 1204 1342 G4double pM = theParticle->GetPDGMass(); << 1343 G4double tM = nucleon->GetPDGMass(); << 1344 1205 1345 G4double pElab = ekin + pM; << 1206 //////////////////////////////////////////////////////////////////////////////////// >> 1207 // >> 1208 // 1346 1209 1347 G4double totEcm = std::sqrt(pM*pM + tM*tM << 1210 G4double G4HadronNucleonXsc::CalculateEcmValue( const G4double mp , >> 1211 const G4double mt , >> 1212 const G4double Plab ) >> 1213 { >> 1214 G4double Elab = std::sqrt ( mp * mp + Plab * Plab ); >> 1215 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt ); >> 1216 // G4double Pcm = Plab * mt / Ecm; >> 1217 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp; 1348 1218 1349 G4double totTcm = totEcm - pM -tM; << 1219 return Ecm ; // KEcm; >> 1220 } 1350 1221 1351 G4double bC = fine_structure_const*hbarc*pZ << 1352 1222 1353 //G4cout<<"##CB: ekin = "<<ekin<<"; pElab= << 1223 //////////////////////////////////////////////////////////////////////////////////// 1354 // <<"; bC = "<<bC<<"; Tcm = "<< totTcm << << 1224 // >> 1225 // 1355 1226 1356 G4double ratio = (totTcm > bC) ? 1. - bC/to << 1227 G4double G4HadronNucleonXsc::CalcMandelstamS( const G4double mp , >> 1228 const G4double mt , >> 1229 const G4double Plab ) >> 1230 { >> 1231 G4double Elab = std::sqrt ( mp * mp + Plab * Plab ); >> 1232 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ; 1357 1233 1358 // G4cout <<"ratio = "<<ratio<<G4endl; << 1234 return sMand; 1359 return ratio; << 1360 } 1235 } 1361 1236 1362 ///////////////////////////////////////////// << 1237 >> 1238 // >> 1239 // >> 1240 /////////////////////////////////////////////////////////////////////////////////////// 1363 1241