| 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 // 14.03.07 V. Grichine - first implementation 26 // 14.03.07 V. Grichine - first implementation
27 // 27 //
28 // 04.09.18 V. Ivantchenko Major revision of i <<
29 // 30.09.18 V. Grichine hyperon-nucleon xsc fi <<
30 28
31 #include "G4HadronNucleonXsc.hh" 29 #include "G4HadronNucleonXsc.hh"
32 #include "G4Element.hh" << 30
33 #include "G4Proton.hh" <<
34 #include "G4Nucleus.hh" <<
35 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
37 #include "G4ParticleTable.hh" 33 #include "G4ParticleTable.hh"
38 #include "G4IonTable.hh" 34 #include "G4IonTable.hh"
>> 35 #include "G4ParticleDefinition.hh"
39 #include "G4HadTmpUtil.hh" 36 #include "G4HadTmpUtil.hh"
40 #include "G4Log.hh" 37 #include "G4Log.hh"
41 #include "G4Exp.hh" 38 #include "G4Exp.hh"
42 #include "G4Pow.hh" 39 #include "G4Pow.hh"
43 #include "G4NuclearRadii.hh" <<
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 <<
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 :
>> 43 // fUpperLimit( 10000 * GeV ),
>> 44 fLowerLimit( 0.03 * MeV ),
>> 45 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0)
>> 46 // , fHadronNucleonXsc(0.0)
67 { 47 {
68 // basic hadrons << 48 theGamma = G4Gamma::Gamma();
69 theProton = G4Proton::Proton(); 49 theProton = G4Proton::Proton();
70 theNeutron = G4Neutron::Neutron(); 50 theNeutron = G4Neutron::Neutron();
>> 51 theAProton = G4AntiProton::AntiProton();
>> 52 theANeutron = G4AntiNeutron::AntiNeutron();
71 thePiPlus = G4PionPlus::PionPlus(); 53 thePiPlus = G4PionPlus::PionPlus();
72 << 54 thePiMinus = G4PionMinus::PionMinus();
73 // basic strange mesons << 55 thePiZero = G4PionZero::PionZero();
74 theKPlus = G4KaonPlus::KaonPlus(); 56 theKPlus = G4KaonPlus::KaonPlus();
75 theKMinus = G4KaonMinus::KaonMinus(); 57 theKMinus = G4KaonMinus::KaonMinus();
76 theK0S = G4KaonZeroShort::KaonZeroShort 58 theK0S = G4KaonZeroShort::KaonZeroShort();
77 theK0L = G4KaonZeroLong::KaonZeroLong() 59 theK0L = G4KaonZeroLong::KaonZeroLong();
>> 60 theL = G4Lambda::Lambda();
>> 61 theAntiL = G4AntiLambda::AntiLambda();
>> 62 theSPlus = G4SigmaPlus::SigmaPlus();
>> 63 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
>> 64 theSMinus = G4SigmaMinus::SigmaMinus();
>> 65 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
>> 66 theS0 = G4SigmaZero::SigmaZero();
>> 67 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
>> 68 theXiMinus = G4XiMinus::XiMinus();
>> 69 theXi0 = G4XiZero::XiZero();
>> 70 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
>> 71 theAXi0 = G4AntiXiZero::AntiXiZero();
>> 72 theOmega = G4OmegaMinus::OmegaMinus();
>> 73 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
>> 74 theD = G4Deuteron::Deuteron();
>> 75 theT = G4Triton::Triton();
>> 76 theA = G4Alpha::Alpha();
>> 77 theHe3 = G4He3::He3();
78 78
79 g4calc = G4Pow::GetInstance(); << 79 // InitialiseKaonNucleonTotXsc();
80 } 80 }
81 81
>> 82
>> 83 G4HadronNucleonXsc::~G4HadronNucleonXsc()
>> 84 {}
>> 85
82 void G4HadronNucleonXsc::CrossSectionDescripti 86 void G4HadronNucleonXsc::CrossSectionDescription(std::ostream& outFile) const
83 { 87 {
84 outFile << "G4HadronNucleonXsc calculates th 88 outFile << "G4HadronNucleonXsc calculates the total, inelastic and elastic\n"
85 << "hadron-nucleon cross sections us 89 << "hadron-nucleon cross sections using several different\n"
86 << "parameterizations within the Gla 90 << "parameterizations within the Glauber-Gribov approach. It is\n"
87 << "valid for all incident gammas an 91 << "valid for all incident gammas and long-lived hadrons at\n"
88 << "energies above 30 keV. This is 92 << "energies above 30 keV. This is a cross section component which\n"
89 << "is to be used to build a cross s 93 << "is to be used to build a cross section data set.\n";
90 } 94 }
91 95
92 G4double G4HadronNucleonXsc::HadronNucleonXsc( << 96 G4bool
93 const G4Particl << 97 G4HadronNucleonXsc::IsApplicable(const G4DynamicParticle* aDP,
>> 98 const G4Element* anElement)
94 { 99 {
95 G4double xsc(0.); << 100 G4int Z = G4lrint(anElement->GetZ());
96 G4int pdg = std::abs(theParticle->GetPDGEnco << 101 G4int A = G4lrint(anElement->GetN());
>> 102 return IsIsoApplicable(aDP, Z, A);
>> 103 }
>> 104
>> 105 ////////////////////////////////////////////////////////////////////////////////////////
>> 106 //
>> 107
>> 108 G4bool
>> 109 G4HadronNucleonXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
>> 110 G4int Z, G4int)
>> 111 {
>> 112 G4bool applicable = false;
>> 113 // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
>> 114 G4double kineticEnergy = aDP->GetKineticEnergy();
>> 115
>> 116 const G4ParticleDefinition* theParticle = aDP->GetDefinition();
>> 117
>> 118 if ( ( kineticEnergy >= fLowerLimit &&
>> 119 Z > 1 && // >= He
>> 120 ( theParticle == theAProton ||
>> 121 theParticle == theGamma ||
>> 122 theParticle == theKPlus ||
>> 123 theParticle == theKMinus ||
>> 124 theParticle == theSMinus) ) ||
>> 125
>> 126 ( kineticEnergy >= 0.1*fLowerLimit &&
>> 127 Z > 1 && // >= He
>> 128 ( theParticle == theProton ||
>> 129 theParticle == theNeutron ||
>> 130 theParticle == thePiPlus ||
>> 131 theParticle == thePiMinus ) ) ) applicable = true;
>> 132
>> 133 return applicable;
>> 134 }
>> 135
>> 136
>> 137 /////////////////////////////////////////////////////////////////////////////////////
>> 138 //
>> 139 // Returns hadron-nucleon Xsc according to differnt parametrisations:
>> 140 // [2] E. Levin, hep-ph/9710546
>> 141 // [3] U. Dersch, et al, hep-ex/9910052
>> 142 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
>> 143
>> 144 G4double
>> 145 G4HadronNucleonXsc::GetHadronNucleonXscEL(const G4DynamicParticle* aParticle,
>> 146 const G4ParticleDefinition* nucleon )
>> 147 {
>> 148 G4double xsection;
>> 149
>> 150
>> 151 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
>> 152
>> 153 G4double proj_mass = aParticle->GetMass();
>> 154 G4double proj_momentum = aParticle->GetMomentum().mag();
>> 155 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
>> 156
>> 157 sMand /= GeV*GeV; // in GeV for parametrisation
>> 158 proj_momentum /= GeV;
>> 159
>> 160 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 161
>> 162 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
97 163
98 // p, n, pi+-, pbar, nbar << 164
99 if ( pdg == 2212 || pdg == 2112 || pdg == 21 << 165 if(theParticle == theGamma && pORn )
100 xsc = HadronNucleonXscNS(theParticle, nucl <<
101 } <<
102 else if ( pdg == 22 ) // gamma <<
103 { 166 {
104 xsc = HadronNucleonXscPDG(theParticle, nuc << 167 xsection = (0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
105 } << 168 }
106 else if ( pdg == 321 || pdg == 310 || pdg == << 169 else if(theParticle == theNeutron && pORn ) // as proton ???
>> 170 {
>> 171 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
>> 172 }
>> 173 else if(theParticle == theProton && pORn )
>> 174 {
>> 175 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
>> 176
>> 177 // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*G4Log(sMand)*G4Log(sMand) );
>> 178 // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
>> 179 }
>> 180 else if(theParticle == theAProton && pORn )
>> 181 {
>> 182 xsection = ( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
>> 183 }
>> 184 else if(theParticle == thePiPlus && pORn )
>> 185 {
>> 186 xsection = (13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
>> 187 }
>> 188 else if(theParticle == thePiMinus && pORn )
>> 189 {
>> 190 // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*G4Log(sMand)*G4Log(sMand) );
>> 191 xsection = (13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
>> 192 }
>> 193 else if(theParticle == theKPlus && pORn )
>> 194 {
>> 195 xsection = (11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
>> 196 }
>> 197 else if(theParticle == theKMinus && pORn )
107 { 198 {
108 xsc = KaonNucleonXscNS(theParticle, nucleo << 199 xsection = (11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
109 } 200 }
110 else if (pdg > 3000) << 201 else // as proton ???
111 { 202 {
112 if (pdg == 3122 || pdg == 3222 || pdg == 3 << 203 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 } 204 }
136 return xsc; << 205 xsection *= millibarn;
>> 206
>> 207 fTotalXsc = xsection;
>> 208 fInelasticXsc = 0.83*xsection;
>> 209 fElasticXsc = fTotalXsc - fInelasticXsc;
>> 210 if (fElasticXsc < 0.)fElasticXsc = 0.;
>> 211
>> 212 return xsection;
137 } 213 }
138 214
139 ////////////////////////////////////////////// << 215
>> 216 /////////////////////////////////////////////////////////////////////////////////////
140 // 217 //
141 // Returns hadron-nucleon Xsc according to PDG << 218 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
142 // http://pdg.lbl.gov/2017/reviews/hadronicrpp << 219 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
>> 220 // At = number of nucleons, Zt = number of protons
143 221
144 G4double G4HadronNucleonXsc::HadronNucleonXscP << 222 G4double
145 const G4ParticleD << 223 G4HadronNucleonXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
146 const G4ParticleDefinition* nucleon << 224 const G4ParticleDefinition* nucleon )
147 { 225 {
148 static const G4double M = 2.1206; // in G << 226 G4double xsection(0);
149 static const G4double eta1 = 0.4473; << 227 G4int Zt=1, Nt=1, At=1;
150 static const G4double eta2 = 0.5486; << 228
151 static const G4double H = 0.272; << 229 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
>> 230
>> 231 G4double proj_mass = aParticle->GetMass();
>> 232 G4double proj_momentum = aParticle->GetMomentum().mag();
152 233
153 G4int pdg = theParticle->GetPDGEncoding(); << 234 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
154 235
155 G4double mass1 = (pdg == 22) ? 770. : thePar << 236 sMand /= GeV*GeV; // in GeV for parametrisation
156 G4double mass2 = nucleon->GetPDGMass(); <<
157 237
158 G4double sMand = CalcMandelstamS(ekin, mass1 << 238 // General PDG fit constants
159 G4double x = (mass1 + mass2)*invGeV + M; <<
160 G4double blog = G4Log(sMand/(x*x)); <<
161 239
162 G4double P(0.0), R1(0.0), R2(0.0), del(1.0); << 240 G4double s0 = 5.38*5.38; // in Gev^2
>> 241 G4double eta1 = 0.458;
>> 242 G4double eta2 = 0.458;
>> 243 G4double B = 0.308;
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 if(theParticle == theNeutron) // proton-neutron fit
168 { 252 {
169 if ( proton ) 253 if ( proton )
170 { 254 {
171 P = 34.71; << 255 xsection = Zt*( 35.80 + B*std::pow(G4Log(sMand/s0),2.)
172 R1 = 12.52; << 256 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));// on p
173 R2 = -6.66; <<
174 } 257 }
175 else << 258 if ( neutron )
176 { 259 {
177 P = 34.41; << 260 xsection = Nt*( 35.45 + B*std::pow(G4Log(sMand/s0),2.)
178 R1 = 13.07; << 261 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // on n pp for nn
179 R2 = -7.394; <<
180 } 262 }
181 } 263 }
182 else if(theParticle == theProton) 264 else if(theParticle == theProton)
183 { 265 {
184 if ( neutron ) << 266 if ( proton )
185 { << 267 {
186 P = 34.71; << 268 xsection = Zt*( 35.45 + B*std::pow(G4Log(sMand/s0),2.)
187 R1 = 12.52; << 269 + 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 } 270 }
196 } <<
197 // pbar <<
198 else if(pdg == -2212) <<
199 { <<
200 if ( neutron ) 271 if ( neutron )
201 { 272 {
202 P = 34.71; << 273 xsection = Nt*( 35.80 + B*std::pow(G4Log(sMand/s0),2.)
203 R1 = 12.52; << 274 + 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 } 275 }
212 } 276 }
213 // nbar << 277 else if(theParticle == theAProton)
214 else if(pdg == -2112) <<
215 { 278 {
216 if ( proton ) 279 if ( proton )
217 { << 280 {
218 P = 34.71; << 281 xsection = Zt*( 35.45 + B*std::pow(G4Log(sMand/s0),2.)
219 R1 = 12.52; << 282 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
220 R2 = 6.66; <<
221 } 283 }
222 else << 284 if ( neutron )
223 { 285 {
224 P = 34.41; << 286 xsection = Nt*( 35.80 + B*std::pow(G4Log(sMand/s0),2.)
225 R1 = 13.07; << 287 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
226 R2 = 7.394; <<
227 } 288 }
228 } 289 }
229 // pi+ << 290 else if(theParticle == thePiPlus && pORn )
230 else if(pdg == 211) <<
231 { 291 {
232 P = 18.75; << 292 xsection = At*( 20.86 + B*std::pow(G4Log(sMand/s0),2.)
233 R1 = 9.56; << 293 + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
234 R2 = -1.767; <<
235 } 294 }
236 // pi- << 295 else if(theParticle == thePiMinus && pORn )
237 else if(pdg == -211) <<
238 { 296 {
239 P = 18.75; << 297 xsection = At*( 20.86 + B*std::pow(G4Log(sMand/s0),2.)
240 R1 = 9.56; << 298 + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
241 R2 = 1.767; <<
242 } 299 }
243 else if(theParticle == theKPlus) 300 else if(theParticle == theKPlus)
244 { 301 {
245 if ( proton ) 302 if ( proton )
246 { << 303 {
247 P = 16.36; << 304 xsection = Zt*( 17.91 + B*std::pow(G4Log(sMand/s0),2.)
248 R1 = 4.29; << 305 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
249 R2 = -3.408; <<
250 } 306 }
251 else << 307 if ( neutron )
252 { 308 {
253 P = 16.31; << 309 xsection = Nt*( 17.87 + B*std::pow(G4Log(sMand/s0),2.)
254 R1 = 3.7; << 310 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
255 R2 = -1.826; <<
256 } 311 }
257 } 312 }
258 else if(theParticle == theKMinus) 313 else if(theParticle == theKMinus)
259 { 314 {
260 if ( proton ) 315 if ( proton )
261 { << 316 {
262 P = 16.36; << 317 xsection = Zt*( 17.91 + B*std::pow(G4Log(sMand/s0),2.)
263 R1 = 4.29; << 318 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
264 R2 = 3.408; <<
265 } 319 }
266 else << 320 if ( neutron )
267 { 321 {
268 P = 16.31; << 322 xsection = Nt*( 17.87 + B*std::pow(G4Log(sMand/s0),2.)
269 R1 = 3.7; << 323 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2) );
270 R2 = 1.826; <<
271 } 324 }
272 } 325 }
273 else if(theParticle == theK0S || theParticle << 326 else if(theParticle == theSMinus && pORn )
274 { 327 {
275 P = 16.36; << 328 xsection = At*( 35.20 + B*std::pow(G4Log(sMand/s0),2.)
276 R1 = 2.5; << 329 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2) );
277 R2 = 0.; <<
278 } <<
279 // sigma- <<
280 else if(pdg == 3112) <<
281 { <<
282 P = 34.7; <<
283 R1 = -46.; <<
284 R2 = 48.; <<
285 } 330 }
286 // gamma << 331 else if(theParticle == theGamma && pORn ) // modify later on
287 else if(pdg == 22) // modify later on <<
288 { 332 {
289 del= 0.003063; << 333 xsection = At*( 0.0 + B*std::pow(G4Log(sMand/s0),2.)
290 P = 34.71*del; << 334 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2) );
291 R1 = (neutron) ? 0.0231 : 0.0139; << 335
292 R2 = 0.; <<
293 } 336 }
294 else // as proton ??? 337 else // as proton ???
295 { 338 {
296 if ( neutron ) << 339 if ( proton )
297 { << 340 {
298 P = 34.71; << 341 xsection = Zt*( 35.45 + B*std::pow(G4Log(sMand/s0),2.)
299 R1 = 12.52; << 342 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2) );
300 R2 = -6.66; <<
301 } 343 }
302 else << 344 if ( neutron )
303 { 345 {
304 P = 34.41; << 346 xsection = Nt*( 35.80 + B*std::pow(G4Log(sMand/s0),2.)
305 R1 = 13.07; << 347 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
306 R2 = -7.394; <<
307 } 348 }
308 } 349 }
309 fTotalXsc = CLHEP::millibarn* << 350 xsection *= millibarn; // parametrised in mb
310 (del*(H*blog*blog + P) + R1*G4Exp(-eta1*bl << 351
311 fInelasticXsc = 0.75*fTotalXsc; << 352 fTotalXsc = xsection;
>> 353 fInelasticXsc = 0.75*xsection;
312 fElasticXsc = fTotalXsc - fInelasticXsc; 354 fElasticXsc = fTotalXsc - fInelasticXsc;
>> 355 if (fElasticXsc < 0.) fElasticXsc = 0.;
313 356
314 if( proton && theParticle->GetPDGCharge() > << 357 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 } 358 }
330 359
331 ////////////////////////////////////////////// << 360
>> 361 /////////////////////////////////////////////////////////////////////////////////////
332 // 362 //
333 // Returns hadron-nucleon cross-section based 363 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
334 // data from mainly http://wwwppds.ihep.su:800 364 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
335 365
336 G4double G4HadronNucleonXsc::HadronNucleonXscN << 366 G4double
337 const G4ParticleD << 367 G4HadronNucleonXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
338 const G4ParticleDefinition* nucleon << 368 const G4ParticleDefinition* nucleon )
339 { << 369 {
340 const G4double ekin = std::max(ekin0, ekinmi << 370 G4double xsection(0);
341 G4int pdg = theParticle->GetPDGEncoding(); << 371
342 /* << 372 G4double A0, B0;
343 G4cout<< "HadronNucleonXscNS: Ekin(GeV)= " < << 373 G4double hpXsc(0);
344 << theParticle->GetParticleName() << " << 374 G4double hnXsc(0);
345 << nucleon->GetParticleName() << G4end << 375
346 */ << 376
347 if(pdg == -2212 || pdg == -2112) { << 377 G4double tM = 0.939*GeV; // ~mean neutron and proton ???
348 return HadronNucleonXscPDG(theParticle, nu << 378
349 } << 379 G4double pM = aParticle->GetMass();
350 << 380 G4double pE = aParticle->GetTotalEnergy();
351 G4double pM = theParticle->GetPDGMass(); << 381 G4double pLab = aParticle->GetMomentum().mag();
352 G4double tM = nucleon->GetPDGMass(); << 382
353 G4double pE = ekin + pM; << 383 G4double sMand = CalcMandelstamS ( pM , tM , pLab );
354 G4double pLab = std::sqrt(ekin*(ekin + 2*pM) << 384
355 << 385 sMand /= GeV*GeV; // in GeV for parametrisation
356 G4double sMand = CalcMandelstamS(ekin, pM, t << 386 pLab /= GeV;
357 << 387 pE /= GeV;
358 pLab *= invGeV; << 388 pM /= GeV;
359 pE *= invGeV; << 389
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); 390 G4double logP = G4Log(pLab);
368 391
>> 392
>> 393 // General PDG fit constants
>> 394
>> 395 G4double s0 = 5.38*5.38; // in Gev^2
>> 396 G4double eta1 = 0.458;
>> 397 G4double eta2 = 0.458;
>> 398 G4double B = 0.308;
>> 399 G4double minLogP = 3.5; // min of (lnP-minLogP)^2
>> 400 G4double cofLogE = .0557; // elastic (lnP-minLogP)^2
>> 401 G4double cofLogT = .3; // total (lnP-minLogP)^2
>> 402 G4double pMin = .1; // fast LE calculation
>> 403 G4double pMax = 1000.; // fast HE calculation
>> 404
>> 405
>> 406 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 407
>> 408 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
369 G4bool proton = (nucleon == theProton); 409 G4bool proton = (nucleon == theProton);
370 G4bool neutron = (nucleon == theNeutron); 410 G4bool neutron = (nucleon == theNeutron);
371 411
372 if( theParticle == theNeutron) << 412 if( theParticle == theNeutron && pORn )
373 { 413 {
374 if( pLab >= 373.) 414 if( pLab >= 373.)
375 { 415 {
376 fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4 << 416 xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
377 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 417
>> 418 fElasticXsc = 6.5 + 0.308*std::pow(G4Log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
>> 419
>> 420 fTotalXsc = xsection;
>> 421
378 } 422 }
379 else if( pLab >= 100.) 423 else if( pLab >= 100.)
380 { 424 {
381 fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G << 425 B0 = 7.5;
382 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 426 A0 = 100. - B0*G4Log(3.0e7);
>> 427
>> 428 xsection = A0 + B0*G4Log(pE) - 11
>> 429 // + 103*std::pow(2*0.93827*pE + pM*pM+0.93827*0.93827,-0.165); // mb
>> 430 + 103*std::pow(sMand,-0.165); // mb
>> 431
>> 432 fElasticXsc = 5.53 + 0.308*std::pow(G4Log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
>> 433
>> 434 fTotalXsc = xsection;
383 } 435 }
384 else if( pLab >= 10.) 436 else if( pLab >= 10.)
385 { 437 {
>> 438 B0 = 7.5;
>> 439 A0 = 100. - B0*G4Log(3.0e7);
>> 440
>> 441 xsection = A0 + B0*G4Log(pE) - 11
>> 442 + 103*std::pow(2*0.93827*pE + pM*pM+
>> 443 0.93827*0.93827,-0.165); // mb
>> 444 fTotalXsc = xsection;
386 fElasticXsc = 6 + 20/( (logP-0.182)*(lo 445 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
387 } 446 }
388 else // pLab < 10 GeV/c 447 else // pLab < 10 GeV/c
389 { 448 {
390 if( neutron ) // nn to be pp 449 if( neutron ) // nn to be pp
391 { 450 {
392 G4double x = G4Log(pLab/0.73); <<
393 if( pLab < 0.4 ) 451 if( pLab < 0.4 )
394 { 452 {
395 fTotalXsc = 23 + 50*std::sqrt(g4calc << 453 hnXsc = 23 + 50*( std::pow( G4Log(0.73/pLab), 3.5 ) );
396 fElasticXsc = fTotalXsc; << 454 fElasticXsc = hnXsc;
397 } 455 }
398 else if( pLab < 0.73 ) 456 else if( pLab < 0.73 )
399 { 457 {
400 fTotalXsc = 23 + 50*std::sqrt(g4calc << 458 hnXsc = 23 + 50*( std::pow( G4Log(0.73/pLab), 3.5 ) );
401 fElasticXsc = fTotalXsc; << 459 fElasticXsc = hnXsc;
402 } 460 }
403 else if( pLab < 1.05 ) 461 else if( pLab < 1.05 )
404 { 462 {
405 fTotalXsc = 23 + 40*x*x; << 463 hnXsc = 23 + 40*(G4Log(pLab/0.73))*
406 fElasticXsc = 23 + 20*x*x; << 464 (G4Log(pLab/0.73));
>> 465 fElasticXsc = 23 + 20*(G4Log(pLab/0.73))*
>> 466 (G4Log(pLab/0.73));
407 } 467 }
408 else // 1.05 - 10 GeV/c 468 else // 1.05 - 10 GeV/c
409 { 469 {
410 fTotalXsc = 39.0+75*(pLab - 1.2)/(g4 << 470 hnXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
>> 471
411 fElasticXsc = 6 + 20/( (logP-0.182) 472 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
412 } 473 }
>> 474 fTotalXsc = hnXsc;
413 } 475 }
414 if( proton ) // pn to be np 476 if( proton ) // pn to be np
415 { 477 {
416 if( pLab < 0.02 ) 478 if( pLab < 0.02 )
417 { 479 {
418 fTotalXsc = 4100+30*G4Exp(G4Log(G4Lo << 480 hpXsc = 4100+30*std::pow(G4Log(1.3/pLab),3.6); // was as pLab < 0.8
419 fElasticXsc = fTotalXsc; << 481 fElasticXsc = hpXsc;
420 } 482 }
421 else if( pLab < 0.8 ) 483 else if( pLab < 0.8 )
422 { 484 {
423 fTotalXsc = 33+30*g4calc->powN(G4Log << 485 hpXsc = 33+30*std::pow(G4Log(pLab/1.3),4.0);
424 fElasticXsc = fTotalXsc; << 486 fElasticXsc = hpXsc;
425 } 487 }
>> 488 else if( pLab < 1.05 )
>> 489 {
>> 490 hpXsc = 33+30*std::pow(G4Log(pLab/0.95),2.0);
>> 491 fElasticXsc = 6 + 52/( G4Log(0.511/pLab)*G4Log(0.511/pLab) + 1.6 );
>> 492 }
426 else if( pLab < 1.4 ) 493 else if( pLab < 1.4 )
427 { 494 {
428 fTotalXsc = 33+30*g4calc->powN(G4Log << 495 hpXsc = 33+30*std::pow(G4Log(pLab/0.95),2.0);
429 G4double x = G4Log(0.511/pLab); << 496 fElasticXsc = 6 + 52/( G4Log(0.511/pLab)*G4Log(0.511/pLab) + 1.6 );
430 fElasticXsc = 6 + 52/( x*x + 1.6 ); <<
431 } 497 }
432 else // 1.4 < pLab < 10. ) 498 else // 1.4 < pLab < 10. )
433 { 499 {
434 fTotalXsc = 33.3 + 20.8*(pLab*pLab - << 500 hpXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
435 /(std::sqrt(g4calc->powN(pLab,5)) + 0.95 << 501
436 fElasticXsc = 6 + 20/( (logP-0.182) 502 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
437 } 503 }
>> 504 fTotalXsc = hpXsc;
438 } 505 }
439 } 506 }
440 } << 507 }
441 ////// proton ////////////////////////////// << 508 else if( theParticle == theProton && pORn ) ////// proton //////////////////////////////////////////////
442 else if(theParticle == theProton) <<
443 { 509 {
444 if( pLab >= 373.) // pdg due to TOTEM data 510 if( pLab >= 373.) // pdg due to TOTEM data
445 { 511 {
446 fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4 << 512 xsection = GetHadronNucleonXscPDG(aParticle, nucleon)/millibarn;
447 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 513
>> 514 fElasticXsc = 6.5 + 0.308*std::pow(G4Log(sMand/400.),1.65) + 9.19*std::pow(sMand,-0.458);
>> 515
>> 516 fTotalXsc = xsection;
448 } 517 }
449 else if( pLab >= 100.) 518 else if( pLab >= 100.)
450 { 519 {
451 fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G << 520 B0 = 7.5;
452 + 9.19*G4Exp(-G4Log(sMand)*0.458); << 521 A0 = 100. - B0*G4Log(3.0e7);
>> 522
>> 523 xsection = A0 + B0*G4Log(pE) - 11 + 103*std::pow(sMand,-0.165); // mb
>> 524
>> 525 fElasticXsc = 5.53 + 0.308*std::pow(G4Log(sMand/28.9),1.1) + 9.19*std::pow(sMand,-0.458);
>> 526
>> 527 fTotalXsc = xsection;
453 } 528 }
454 else if( pLab >= 10.) 529 else if( pLab >= 10.)
455 { 530 {
456 fElasticXsc = 6. + 20./( (logP-0.182)*(l << 531 B0 = 7.5;
>> 532 A0 = 100. - B0*G4Log(3.0e7);
>> 533
>> 534 xsection = A0 + B0*G4Log(pE) - 11 + 103*std::pow(sMand,-0.165); // mb
>> 535
>> 536 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
>> 537
>> 538 fTotalXsc = xsection;
457 } 539 }
458 else 540 else
459 { 541 {
460 // pp 542 // pp
>> 543
461 if( proton ) 544 if( proton )
462 { 545 {
463 if( pLab < 0.73 ) << 546 if( pLab < 0.4 )
>> 547 {
>> 548 hpXsc = 23 + 50*( std::pow( G4Log(0.73/pLab), 3.5 ) );
>> 549 fElasticXsc = hpXsc;
>> 550 }
>> 551 else if( pLab < 0.73 )
464 { 552 {
465 fTotalXsc = 23 + 50*std::sqrt(g4calc << 553 hpXsc = 23 + 50*( std::pow( G4Log(0.73/pLab), 3.5 ) );
466 fElasticXsc = fTotalXsc; << 554 fElasticXsc = hpXsc;
467 } 555 }
468 else if( pLab < 1.05 ) 556 else if( pLab < 1.05 )
469 { 557 {
470 G4double x = G4Log(pLab/0.73); << 558 hpXsc = 23 + 40*(G4Log(pLab/0.73))*
471 fTotalXsc = 23 + 40*x*x; << 559 (G4Log(pLab/0.73));
472 fElasticXsc = 23 + 20*x*x; << 560 fElasticXsc = 23 + 20*(G4Log(pLab/0.73))*
>> 561 (G4Log(pLab/0.73));
473 } 562 }
474 else // 1.05 - 10 GeV/c 563 else // 1.05 - 10 GeV/c
475 { 564 {
476 fTotalXsc = 39.0+75*(pLab - 1.2)/(g4 << 565 hpXsc = 39.0+75*(pLab - 1.2)/(std::pow(pLab,3.0) + 0.15);
477 fElasticXsc = 6. + 20./( (logP-0.182 << 566
>> 567 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
478 } 568 }
>> 569 fTotalXsc = hpXsc;
479 } 570 }
480 else if( neutron ) // pn to be np << 571 if( neutron ) // pn to be np
481 { 572 {
482 if( pLab < 0.02 ) 573 if( pLab < 0.02 )
483 { 574 {
484 fTotalXsc = 4100+30*G4Exp(G4Log(G4Lo << 575 hnXsc = 4100+30*std::pow(G4Log(1.3/pLab),3.6); // was as pLab < 0.8
485 fElasticXsc = fTotalXsc; << 576 fElasticXsc = hnXsc;
486 } 577 }
487 else if( pLab < 0.8 ) 578 else if( pLab < 0.8 )
488 { 579 {
489 fTotalXsc = 33+30*g4calc->powN(G4Log << 580 hnXsc = 33+30*std::pow(G4Log(pLab/1.3),4.0);
490 fElasticXsc = fTotalXsc; << 581 fElasticXsc = hnXsc;
491 } 582 }
>> 583 else if( pLab < 1.05 )
>> 584 {
>> 585 hnXsc = 33+30*std::pow(G4Log(pLab/0.95),2.0);
>> 586 fElasticXsc = 6 + 52/( G4Log(0.511/pLab)*G4Log(0.511/pLab) + 1.6 );
>> 587 }
492 else if( pLab < 1.4 ) 588 else if( pLab < 1.4 )
493 { 589 {
494 G4double x1 = G4Log(pLab/0.95); << 590 hnXsc = 33+30*std::pow(G4Log(pLab/0.95),2.0);
495 G4double x2 = G4Log(0.511/pLab); << 591 fElasticXsc = 6 + 52/( G4Log(0.511/pLab)*G4Log(0.511/pLab) + 1.6 );
496 fTotalXsc = 33+30*x1*x1; <<
497 fElasticXsc = 6 + 52/(x2*x2 + 1.6); <<
498 } 592 }
499 else // 1.4 < pLab < 10. ) 593 else // 1.4 < pLab < 10. )
500 { 594 {
501 fTotalXsc = 33.3 + 20.8*(pLab*pLab - << 595 hnXsc = 33.3 + 20.8*(std::pow(pLab,2.0) - 1.35)/(std::pow(pLab,2.50) + 0.95);
502 /(std::sqrt(g4calc->powN(pLab,5)) + 0.95 << 596
503 fElasticXsc = 6. + 20./( (logP-0.182 << 597 fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
504 } 598 }
>> 599 fTotalXsc = hnXsc;
505 } 600 }
506 } 601 }
507 } << 602 }
508 // pi+ p; pi- n << 603 else if( theParticle == theAProton && pORn ) /////////////////// p_bar ///////////////////////////
509 else if((pdg == 211 && proton) || (pdg == -2 <<
510 { 604 {
511 if( pLab < 0.28 ) << 605 if( proton )
512 { 606 {
513 fTotalXsc = 10./((logP + 1.273)*(logP + << 607 xsection = 35.45 + B*std::pow(G4Log(sMand/s0),2.)
514 fElasticXsc = fTotalXsc; << 608 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2);
515 } 609 }
516 else if( pLab < 0.68 ) << 610 if( neutron ) // ???
517 { 611 {
518 fTotalXsc = 14./((logP + 1.273)*(logP + << 612 xsection = 35.80 + B*std::pow(G4Log(sMand/s0),2.)
519 fElasticXsc = fTotalXsc; << 613 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2);
520 } 614 }
521 else if( pLab < 0.85 ) << 615 fTotalXsc = xsection;
522 { << 616 }
523 G4double x = G4Log(pLab/0.77); << 617 else if( theParticle == thePiPlus && pORn ) // pi+ /////////////////////////////////////////////
524 fTotalXsc = 88.*x*x + 14.9; <<
525 fElasticXsc = fTotalXsc*G4Exp(-3.*(pLab <<
526 } <<
527 else if( pLab < 1.15 ) <<
528 { <<
529 G4double x = G4Log(pLab/0.77); <<
530 fTotalXsc = 88.*x*x + 14.9; <<
531 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( <<
532 } <<
533 else if( pLab < 1.4) // ns original <<
534 { <<
535 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p <<
536 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL <<
537 fTotalXsc = Ex1 + Ex2 + 27.5; <<
538 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( <<
539 } <<
540 else if( pLab < 2.0 ) // ns original <<
541 { <<
542 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p <<
543 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL <<
544 fTotalXsc = Ex1 + Ex2 + 27.5; <<
545 fElasticXsc = 3.0 + 1.36/( (logP - 0.336 <<
546 } <<
547 else if( pLab < 3.5 ) // ns original <<
548 { <<
549 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p <<
550 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL <<
551 fTotalXsc = Ex1 + Ex2 + 27.5; <<
552 fElasticXsc = 3.0 + 6.20/( (logP - 0.336 <<
553 } <<
554 else if( pLab < 10. ) // my <<
555 { <<
556 fTotalXsc = 10.6 + 2.*G4Log(pE) + 25*G4E <<
557 fElasticXsc = 3.0 + 6.20/( (logP - 0.336 <<
558 } <<
559 else // pLab > 10 // my <<
560 { <<
561 fElasticXsc = 3.0 + 6.20/( (logP - 0.336 <<
562 } <<
563 } <<
564 // pi+ n; pi- p <<
565 else if((pdg == 211 && neutron) || (pdg == - <<
566 { 618 {
567 if( pLab < 0.28 ) << 619 if( proton ) // pi+ p
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 { <<
574 fTotalXsc = 0.648/((pLab - 0.28)*(pLab - <<
575 fElasticXsc = 0.257/((pLab - 0.28)*(pLab <<
576 } <<
577 else if( pLab < 0.5 ) <<
578 { <<
579 G4double y = G4Log(pLab/0.48); <<
580 fTotalXsc = 26 + 110*y*y; <<
581 fElasticXsc = 0.37*fTotalXsc; <<
582 } <<
583 else if( pLab < 0.65 ) <<
584 { <<
585 G4double x = G4Log(pLab/0.48); <<
586 fTotalXsc = 26. + 110.*x*x; <<
587 fElasticXsc = 0.95/((pLab - 0.72)*(pLab <<
588 } <<
589 else if( pLab < 0.72 ) <<
590 { <<
591 fTotalXsc = 36.1 + 10*G4Exp(-(pLab-0.72) <<
592 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. <<
593 fElasticXsc = 0.95/((pLab - 0.72)*(pLab <<
594 } <<
595 else if( pLab < 0.88 ) <<
596 { <<
597 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 <<
598 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. <<
599 fElasticXsc = 0.95/((pLab - 0.72)*(pLab <<
600 } <<
601 else if( pLab < 1.03 ) <<
602 { 620 {
603 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 << 621 if( pLab < 0.28 )
604 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 622 {
605 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*( << 623 hpXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
606 } << 624 fElasticXsc = hpXsc;
607 else if( pLab < 1.15 ) << 625 }
608 { << 626 else if( pLab < 0.4 )
609 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 << 627 {
610 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 628 hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
611 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*( << 629 fElasticXsc = hpXsc;
612 } << 630 }
613 else if( pLab < 1.3 ) << 631 else if( pLab < 0.68 )
>> 632 {
>> 633 hpXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
>> 634 fElasticXsc = hpXsc;
>> 635 }
>> 636 else if( pLab < 0.85 )
>> 637 {
>> 638 G4double Ex4 = 88*(G4Log(pLab/0.77))*(G4Log(pLab/0.77));
>> 639 hpXsc = Ex4 + 14.9;
>> 640 fElasticXsc = hpXsc*G4Exp(-3.*(pLab - 0.68));
>> 641 }
>> 642 else if( pLab < 1.15 )
>> 643 {
>> 644 G4double Ex4 = 88*(G4Log(pLab/0.77))*(G4Log(pLab/0.77));
>> 645 hpXsc = Ex4 + 14.9;
>> 646
>> 647 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
>> 648 }
>> 649 else if( pLab < 1.4) // ns original
>> 650 {
>> 651 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
>> 652 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
>> 653 hpXsc = Ex1 + Ex2 + 27.5;
>> 654 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
>> 655 }
>> 656 else if( pLab < 2.0 ) // ns original
>> 657 {
>> 658 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
>> 659 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
>> 660 hpXsc = Ex1 + Ex2 + 27.5;
>> 661 fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
>> 662 }
>> 663 else if( pLab < 3.5 ) // ns original
>> 664 {
>> 665 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
>> 666 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
>> 667 hpXsc = Ex1 + Ex2 + 27.5;
>> 668 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
>> 669 }
>> 670 else if( pLab < 200. ) // my
>> 671 {
>> 672 hpXsc = 10.6 + 2.*G4Log(pE) + 25*std::pow(pE, -0.43 ); // ns original
>> 673 // hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
>> 674 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
>> 675 }
>> 676 else // pLab > 100 // my
>> 677 {
>> 678 hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
>> 679 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
>> 680 }
>> 681 fTotalXsc = hpXsc;
>> 682 }
>> 683 if( neutron ) // pi+ n = pi- p??
614 { 684 {
615 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72 << 685 if( pLab < 0.28 )
616 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0. << 686 {
617 fElasticXsc = 3. + 13./pLab; << 687 hnXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
>> 688 fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
>> 689 }
>> 690 else if( pLab < 0.395676 ) // first peak
>> 691 {
>> 692 hnXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
>> 693 fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
>> 694 }
>> 695 else if( pLab < 0.5 )
>> 696 {
>> 697 hnXsc = 26 + 110*(G4Log(pLab/0.48))*(G4Log(pLab/0.48));
>> 698 fElasticXsc = 0.37*hnXsc;
>> 699 }
>> 700 else if( pLab < 0.65 )
>> 701 {
>> 702 hnXsc = 26 + 110*(G4Log(pLab/0.48))*(G4Log(pLab/0.48));
>> 703 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
>> 704 }
>> 705 else if( pLab < 0.72 )
>> 706 {
>> 707 hnXsc = 36.1 + 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 708 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 709 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
>> 710 }
>> 711 else if( pLab < 0.88 )
>> 712 {
>> 713 hnXsc = 36.1 + 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 714 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 715 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
>> 716 }
>> 717 else if( pLab < 1.03 )
>> 718 {
>> 719 hnXsc = 36.1 + 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 720 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 721 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
>> 722 }
>> 723 else if( pLab < 1.15 )
>> 724 {
>> 725 hnXsc = 36.1 + 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 726 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 727 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
>> 728 }
>> 729 else if( pLab < 1.3 )
>> 730 {
>> 731 hnXsc = 36.1 + 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 732 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 733 fElasticXsc = 3. + 13./pLab;
>> 734 }
>> 735 else if( pLab < 2.6 ) // < 3.0) // ns original
>> 736 {
>> 737 hnXsc = 36.1 + 0.079-4.313*G4Log(pLab)+
>> 738 3*G4Exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
>> 739 1.5*G4Exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
>> 740 fElasticXsc = 3. + 13./pLab;
>> 741 }
>> 742 else if( pLab < 20. ) // < 3.0) // ns original
>> 743 {
>> 744 hnXsc = 36.1 + 0.079 - 4.313*G4Log(pLab)+
>> 745 3*G4Exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
>> 746 1.5*G4Exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
>> 747 fElasticXsc = 3. + 13./pLab;
>> 748 }
>> 749 else // mb
>> 750 {
>> 751 hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
>> 752 fElasticXsc = 3. + 13./pLab;
>> 753 }
>> 754 fTotalXsc = hnXsc;
618 } 755 }
619 else if( pLab < 10.) // < 3.0) // ns origi << 756 }
>> 757 else if( theParticle == thePiMinus && pORn ) /// pi- ////////////////////////////////////////////
>> 758 {
>> 759 if( neutron ) // pi- n = pi+ p??
620 { 760 {
621 fTotalXsc = 36.1 + 0.079-4.313*logP+ << 761 if( pLab < 0.28 )
622 3*G4Exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+ << 762 {
623 1.5*G4Exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12); << 763 hnXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
624 fElasticXsc = 3. + 13./pLab; << 764 fElasticXsc = hnXsc;
>> 765 }
>> 766 else if( pLab < 0.4 )
>> 767 {
>> 768 hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
>> 769 fElasticXsc = hnXsc;
>> 770 }
>> 771 else if( pLab < 0.68 )
>> 772 {
>> 773 hnXsc = 14./( (logP + 1.273)*(logP + 1.273) + 0.07);
>> 774 fElasticXsc = hnXsc;
>> 775 }
>> 776 else if( pLab < 0.85 )
>> 777 {
>> 778 G4double Ex4 = 88*(G4Log(pLab/0.77))*(G4Log(pLab/0.77));
>> 779 hnXsc = Ex4 + 14.9;
>> 780 fElasticXsc = hnXsc*G4Exp(-3.*(pLab - 0.68));
>> 781 }
>> 782 else if( pLab < 1.15 )
>> 783 {
>> 784 G4double Ex4 = 88*(G4Log(pLab/0.77))*(G4Log(pLab/0.77));
>> 785 hnXsc = Ex4 + 14.9;
>> 786
>> 787 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
>> 788 }
>> 789 else if( pLab < 1.4) // ns original
>> 790 {
>> 791 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
>> 792 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
>> 793 hnXsc = Ex1 + Ex2 + 27.5;
>> 794 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
>> 795 }
>> 796 else if( pLab < 2.0 ) // ns original
>> 797 {
>> 798 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
>> 799 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
>> 800 hnXsc = Ex1 + Ex2 + 27.5;
>> 801 fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
>> 802 }
>> 803 else if( pLab < 3.5 ) // ns original
>> 804 {
>> 805 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
>> 806 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
>> 807 hnXsc = Ex1 + Ex2 + 27.5;
>> 808 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
>> 809 }
>> 810 else if( pLab < 200. ) // my
>> 811 {
>> 812 hnXsc = 10.6 + 2.*G4Log(pE) + 25*std::pow(pE, -0.43 ); // ns original
>> 813 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
>> 814 }
>> 815 else // pLab > 100 // my
>> 816 {
>> 817 hnXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
>> 818 fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
>> 819 }
>> 820 fTotalXsc = hnXsc;
625 } 821 }
626 else // mb << 822 if( proton ) // pi- p
627 { 823 {
628 fElasticXsc = 3. + 13./pLab; << 824 if( pLab < 0.28 )
>> 825 {
>> 826 hpXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
>> 827 fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
>> 828 }
>> 829 else if( pLab < 0.395676 ) // first peak
>> 830 {
>> 831 hpXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
>> 832 fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
>> 833 }
>> 834 else if( pLab < 0.5 )
>> 835 {
>> 836 hpXsc = 26 + 110*(G4Log(pLab/0.48))*(G4Log(pLab/0.48));
>> 837 fElasticXsc = 0.37*hpXsc;
>> 838 }
>> 839 else if( pLab < 0.65 )
>> 840 {
>> 841 hpXsc = 26 + 110*(G4Log(pLab/0.48))*(G4Log(pLab/0.48));
>> 842 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
>> 843 }
>> 844 else if( pLab < 0.72 )
>> 845 {
>> 846 hpXsc = 36.1+
>> 847 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 848 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 849 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
>> 850 }
>> 851 else if( pLab < 0.88 )
>> 852 {
>> 853 hpXsc = 36.1+
>> 854 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 855 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 856 fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
>> 857 }
>> 858 else if( pLab < 1.03 )
>> 859 {
>> 860 hpXsc = 36.1+
>> 861 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 862 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 863 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
>> 864 }
>> 865 else if( pLab < 1.15 )
>> 866 {
>> 867 hpXsc = 36.1+
>> 868 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 869 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 870 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
>> 871 }
>> 872 else if( pLab < 1.3 )
>> 873 {
>> 874 hpXsc = 36.1+
>> 875 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
>> 876 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
>> 877 fElasticXsc = 3. + 13./pLab;
>> 878 }
>> 879 else if( pLab < 2.6 ) // < 3.0) // ns original
>> 880 {
>> 881 hpXsc = 36.1+0.079-4.313*G4Log(pLab)+
>> 882 3*G4Exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
>> 883 1.5*G4Exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
>> 884 fElasticXsc = 3. +13./pLab; // *G4Log(pLab*6.79);
>> 885 }
>> 886 else // mb
>> 887 {
>> 888 hpXsc = GetHadronNucleonXscPDG(aParticle, nucleon )/millibarn;
>> 889 fElasticXsc = 3. + 13./pLab;
>> 890 }
>> 891 fTotalXsc = hpXsc;
629 } 892 }
630 } << 893 }
631 else if( (theParticle == theKMinus) && proto << 894 else if( (theParticle == theKMinus || theParticle == theK0S) && proton ) // Kmp/K0p /////////////////////////////////
632 { 895 {
633 if( pLab < pMin) 896 if( pLab < pMin)
634 { 897 {
635 G4double psp = pLab*std::sqrt(pLab); 898 G4double psp = pLab*std::sqrt(pLab);
636 fElasticXsc = 5.2/psp; 899 fElasticXsc = 5.2/psp;
637 fTotalXsc = 14./psp; 900 fTotalXsc = 14./psp;
638 } 901 }
639 else if( pLab > pMax ) 902 else if( pLab > pMax )
640 { 903 {
641 G4double ld = logP - minLogP; << 904 G4double ld = G4Log(pLab) - minLogP;
642 G4double ld2 = ld*ld; 905 G4double ld2 = ld*ld;
643 fElasticXsc = cofLogE*ld2 + 2.23; << 906 fElasticXsc = cofLogE*ld2 + 2.23;
644 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 907 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
645 } 908 }
646 else 909 else
647 { 910 {
648 G4double ld = logP - minLogP; << 911 G4double ld = G4Log(pLab) - minLogP;
649 G4double ld2 = ld*ld; 912 G4double ld2 = ld*ld;
650 G4double sp = std::sqrt(pLab); 913 G4double sp = std::sqrt(pLab);
651 G4double psp = pLab*sp; 914 G4double psp = pLab*sp;
652 G4double p2 = pLab*pLab; 915 G4double p2 = pLab*pLab;
653 G4double p4 = p2*p2; 916 G4double p4 = p2*p2;
654 <<
655 G4double lh = pLab - 1.01; <<
656 G4double hd = lh*lh + .011; <<
657 <<
658 G4double lm = pLab - .39; 917 G4double lm = pLab - .39;
659 G4double md = lm*lm + .000356; 918 G4double md = lm*lm + .000356;
>> 919
660 G4double lh1 = pLab - 0.78; 920 G4double lh1 = pLab - 0.78;
661 G4double hd1 = lh1*lh1 + .00166; 921 G4double hd1 = lh1*lh1 + .00166;
>> 922
>> 923 G4double lh = pLab - 1.01;
>> 924 G4double hd = lh*lh + .011;
>> 925
662 G4double lh2 = pLab - 1.63; 926 G4double lh2 = pLab - 1.63;
663 G4double hd2 = lh2*lh2 + .007; 927 G4double hd2 = lh2*lh2 + .007;
664 928
665 // small peaks were added but commented <<
666 fElasticXsc = 5.2/psp + (1.1*cofLogE*ld 929 fElasticXsc = 5.2/psp + (1.1*cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4)
667 + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd; << 930 + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd; // small peaks were added
668 931
669 fTotalXsc = 14./psp + (1.1*cofLogT*ld2 << 932 fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
670 + .006/md + 0.01/hd1+ 0.02/hd2 + .20/ << 933 + .006/md + 0.01/hd1+ 0.02/hd2 + .20/hd ;
671 } 934 }
672 } 935 }
673 else if( (theParticle == theKMinus) && neutr << 936 else if( (theParticle == theKMinus || theParticle == theK0S) && neutron ) // Kmn/K0n //////////////////////////////
674 { 937 {
675 if( pLab > pMax ) 938 if( pLab > pMax )
676 { 939 {
677 G4double ld = logP - minLogP; << 940 G4double ld = G4Log(pLab) - minLogP;
678 G4double ld2 = ld*ld; 941 G4double ld2 = ld*ld;
679 fElasticXsc = cofLogE*ld2 + 2.23; << 942 fElasticXsc = cofLogE*ld2 + 2.23;
680 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 943 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
681 } 944 }
682 else 945 else
683 { << 946 {
>> 947
684 G4double lh = pLab - 0.98; 948 G4double lh = pLab - 0.98;
685 G4double hd = lh*lh + .021; << 949 G4double hd = lh*lh + .021;
686 G4double sqrLogPlab = logP*logP; <<
687 950
688 fElasticXsc = 5.0 + 8.1*G4Exp(-logP*1.8 << 951 G4double LogPlab = G4Log( pLab );
689 - 1.3*logP + .15/hd; << 952 G4double sqrLogPlab = LogPlab * LogPlab;
690 fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2. << 953
>> 954 fElasticXsc = // 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md
>> 955 5.0 + 8.1*std::pow(pLab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab + .15/hd;
>> 956 fTotalXsc = // 14./psp +
>> 957 // (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
>> 958 // WP 25.2 + 0. *std::pow(pLab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab
>> 959 25.2 + 0.38*sqrLogPlab - 2.9*LogPlab
>> 960 // + .006/md + 0.01/hd1+ 0.02/hd2
>> 961 + 0.30/hd ;
691 } 962 }
692 } 963 }
693 else if( (theParticle == theKPlus) && proton << 964 else if( (theParticle == theKPlus || theParticle == theK0L) && proton ) // Kpp/aKp ////////////////////////
694 { 965 {
695 // VI: modified low-energy part << 966 if( pLab < pMin )
696 if( pLab < 0.631 ) <<
697 { 967 {
698 fElasticXsc = fTotalXsc = 12.03; << 968 G4double lr = pLab - .38;
>> 969 G4double lm = pLab - 1.;
>> 970 G4double md = lm*lm + .392;
>> 971 fElasticXsc = .7/(lr*lr + .076) + 2./md;
>> 972 fTotalXsc = .7/(lr*lr + .076) + 2.6/md;
699 } 973 }
700 else if( pLab > pMax ) 974 else if( pLab > pMax )
701 { 975 {
702 G4double ld = logP - minLogP; << 976 G4double ld = G4Log(pLab) - minLogP;
703 G4double ld2 = ld*ld; 977 G4double ld2 = ld*ld;
704 fElasticXsc = cofLogE*ld2 + 2.23; << 978 fElasticXsc = cofLogE*ld2 + 2.23;
705 fTotalXsc = cofLogT*ld2 + 19.2; << 979 fTotalXsc = cofLogT*ld2 + 19.2;
706 } 980 }
707 else 981 else
708 { 982 {
709 G4double ld = logP - minLogP; << 983 G4double ld = G4Log(pLab) - minLogP;
710 G4double ld2 = ld*ld; 984 G4double ld2 = ld*ld;
711 G4double lr = pLab - .38; 985 G4double lr = pLab - .38;
712 G4double LE = .7/(lr*lr + .076); 986 G4double LE = .7/(lr*lr + .076);
713 G4double sp = std::sqrt(pLab); 987 G4double sp = std::sqrt(pLab);
714 G4double p2 = pLab*pLab; 988 G4double p2 = pLab*pLab;
715 G4double p4 = p2*p2; 989 G4double p4 = p2*p2;
716 // VI: tuned elastic << 990 G4double lm = pLab - 1.;
717 fElasticXsc = LE + (cofLogE*ld2 + 2.23) << 991 G4double md = lm*lm + .392;
718 + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652); << 992 fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
719 fTotalXsc = LE + (cofLogT*ld2 + 19.5) << 993 fTotalXsc = LE + (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 2.6/md;
720 + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392); <<
721 } 994 }
722 } 995 }
723 else if( (theParticle == theKPlus) && neutr << 996 else if( (theParticle == theKPlus || theParticle == theK0L) && neutron ) // Kpn/aKn ///////////////////////
724 { 997 {
725 if( pLab < pMin ) 998 if( pLab < pMin )
726 { 999 {
727 G4double lm = pLab - 0.94; 1000 G4double lm = pLab - 0.94;
728 G4double md = lm*lm + .392; 1001 G4double md = lm*lm + .392;
729 fElasticXsc = 2./md; 1002 fElasticXsc = 2./md;
730 fTotalXsc = 4.6/md; 1003 fTotalXsc = 4.6/md;
731 } 1004 }
732 else if( pLab > pMax ) 1005 else if( pLab > pMax )
733 { 1006 {
734 G4double ld = logP - minLogP; << 1007 G4double ld = G4Log(pLab) - minLogP;
735 G4double ld2 = ld*ld; 1008 G4double ld2 = ld*ld;
736 fElasticXsc = cofLogE*ld2 + 2.23; << 1009 fElasticXsc = cofLogE*ld2 + 2.23;
737 fTotalXsc = cofLogT*ld2 + 19.2; << 1010 fTotalXsc = cofLogT*ld2 + 19.2;
738 } 1011 }
739 else 1012 else
740 { 1013 {
741 G4double ld = logP - minLogP; << 1014 G4double ld = G4Log(pLab) - minLogP;
742 G4double ld2 = ld*ld; 1015 G4double ld2 = ld*ld;
743 G4double sp = std::sqrt(pLab); 1016 G4double sp = std::sqrt(pLab);
744 G4double p2 = pLab*pLab; 1017 G4double p2 = pLab*pLab;
745 G4double p4 = p2*p2; 1018 G4double p4 = p2*p2;
746 G4double lm = pLab - 0.94; 1019 G4double lm = pLab - 0.94;
747 G4double md = lm*lm + .392; << 1020 G4double md = lm*lm + .392;
748 fElasticXsc = (cofLogE*ld2 + 2.23)/(1. 1021 fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
749 fTotalXsc = (cofLogT*ld2 + 19.5)/(1. 1022 fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 4.6/md;
750 } 1023 }
751 } 1024 }
752 fTotalXsc *= CLHEP::millibarn; << 1025 else if( theParticle == theSMinus && pORn )
753 fElasticXsc *= CLHEP::millibarn; << 1026 {
754 fElasticXsc = std::min(fElasticXsc, fTotalX << 1027 xsection = 35.20 + B*std::pow(G4Log(sMand/s0),2.)
>> 1028 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2);
>> 1029 }
>> 1030 else if( theParticle == theGamma && pORn ) // modify later on
>> 1031 {
>> 1032 xsection = 0.0 + B*std::pow(G4Log(sMand/s0),2.)
>> 1033 + 0.032*std::pow(sMand,-eta1); // WP - 0.0*std::pow(sMand,-eta2);
>> 1034 fTotalXsc = xsection;
>> 1035 }
>> 1036 else // other then p,n,pi+,pi-,K+,K- as proton ???
>> 1037 {
>> 1038 if( proton )
>> 1039 {
>> 1040 xsection = 35.45 + B*std::pow(G4Log(sMand/s0),2.)
>> 1041 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2);
>> 1042 }
>> 1043 if( neutron )
>> 1044 {
>> 1045 xsection += 35.80 + B*std::pow(G4Log(sMand/s0),2.)
>> 1046 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2);
>> 1047 }
>> 1048 fTotalXsc = xsection;
>> 1049 }
>> 1050 fTotalXsc *= millibarn; // parametrised in mb
>> 1051 fElasticXsc *= millibarn; // parametrised in mb
755 1052
756 if( proton && theParticle->GetPDGCharge() > << 1053 if( proton && aParticle->GetDefinition()->GetPDGCharge() > 0. )
757 { 1054 {
758 G4double cB = G4NuclearRadii::CoulombFacto << 1055 G4double cB = GetCoulombBarrier(aParticle, nucleon);
759 fTotalXsc *= cB; 1056 fTotalXsc *= cB;
760 fElasticXsc *= cB; 1057 fElasticXsc *= cB;
761 } 1058 }
762 fInelasticXsc = std::max(fTotalXsc - fElasti << 1059 fInelasticXsc = fTotalXsc - fElasticXsc;
763 /* << 1060 if( fInelasticXsc < 0. ) fInelasticXsc = 0.;
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 } <<
772 <<
773 ////////////////////////////////////////////// <<
774 // <<
775 // Returns kaon-nucleon cross-section based on <<
776 // tuned for the Glauber-Gribov hadron model f <<
777 <<
778 G4double G4HadronNucleonXsc::KaonNucleonXscGG( <<
779 const G4ParticleD <<
780 const G4ParticleDefinition* nucleon <<
781 { <<
782 fTotalXsc = fElasticXsc = fInelasticXsc = 0. <<
783 if(theParticle == theKMinus || theParticle = <<
784 KaonNucleonXscVG(theParticle, nucleon, eki <<
785 <<
786 } else if(theParticle == theK0S || thePartic <<
787 G4double stot = KaonNucleonXscVG(theKMinu <<
788 G4double sel = fElasticXsc; <<
789 G4double sinel = fInelasticXsc; <<
790 stot += KaonNucleonXscVG(theKPlus, nucleo <<
791 sel += fElasticXsc; <<
792 sinel += fInelasticXsc; <<
793 fTotalXsc = stot*0.5; <<
794 fElasticXsc = sel*0.5; <<
795 fInelasticXsc = sinel*0.5; <<
796 } <<
797 return fTotalXsc; <<
798 } <<
799 1061
800 ////////////////////////////////////////////// << 1062 // G4cout<<fTotalXsc/millibarn<<"; "<<fElasticXsc/millibarn<<"; "<<fInelasticXsc/millibarn<<G4endl;
801 // <<
802 // Returns kaon-nucleon cross-section using NS <<
803 1063
804 G4double G4HadronNucleonXsc::KaonNucleonXscNS( <<
805 const G4ParticleD <<
806 const G4ParticleDefinition* nucleon <<
807 { <<
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 } <<
837 return fTotalXsc; 1064 return fTotalXsc;
838 } 1065 }
839 1066
840 ////////////////////////////////////////////// << 1067 /////////////////////////////////////////////////////////////////////////////////////
841 // 1068 //
842 // Returns kaon-nucleon cross-section with smo << 1069 // Returns kaon-nucleon cross-section based on smoothed NS for GG model
843 1070
844 G4double G4HadronNucleonXsc::KaonNucleonXscVG( << 1071 G4double
845 const G4ParticleDefinition* t << 1072 G4HadronNucleonXsc::GetKaonNucleonXscGG(const G4DynamicParticle* aParticle,
846 const G4ParticleDefinition* nucleon, G4do << 1073 const G4ParticleDefinition* nucleon )
847 { 1074 {
848 G4double pM = theParticle->GetPDGMass(); << 1075 G4double pLab = aParticle->GetMomentum().mag();
849 G4double pLab = std::sqrt(ekin*(ekin + 2*pM) <<
850 1076
851 pLab *= invGeV; << 1077 pLab /= GeV;
852 G4double logP = G4Log(pLab); << 1078 G4double LogPlab = G4Log( pLab );
>> 1079 G4double sqrLogPlab = LogPlab * LogPlab;
>> 1080
>> 1081 G4double minLogP = 3.5; // min of (lnP-minLogP)^2
>> 1082 G4double cofLogE = .0557; // elastic (lnP-minLogP)^2
>> 1083 G4double cofLogT = .3; // total (lnP-minLogP)^2
>> 1084 G4double pMin = .1; // fast LE calculation
>> 1085 G4double pMax = 1000.; // fast HE calculation
853 1086
854 fTotalXsc = 0.0; << 1087 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
855 1088
856 G4bool proton = (nucleon == theProton); 1089 G4bool proton = (nucleon == theProton);
857 G4bool neutron = (nucleon == theNeutron); 1090 G4bool neutron = (nucleon == theNeutron);
858 1091
859 if( (theParticle == theKMinus) && proton ) << 1092 if( (theParticle == theKMinus || theParticle == theK0S) && proton ) // (K-,K0)on p ////////////////////////////
860 { 1093 {
>> 1094
861 if( pLab < pMin) 1095 if( pLab < pMin)
862 { 1096 {
863 G4double psp = pLab*std::sqrt(pLab); 1097 G4double psp = pLab*std::sqrt(pLab);
864 fElasticXsc = 5.2/psp; 1098 fElasticXsc = 5.2/psp;
865 fTotalXsc = 14./psp; 1099 fTotalXsc = 14./psp;
866 } 1100 }
867 else if( pLab > pMax ) 1101 else if( pLab > pMax )
868 { 1102 {
869 G4double ld = logP - minLogP; << 1103 G4double ld = G4Log(pLab) - minLogP;
870 G4double ld2 = ld*ld; 1104 G4double ld2 = ld*ld;
871 fElasticXsc = cofLogE*ld2 + 2.23; << 1105 fElasticXsc = cofLogE*ld2 + 2.23;
872 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 1106 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
873 } 1107 }
874 else 1108 else
875 { 1109 {
876 G4double ld = logP - minLogP; << 1110 G4double ld = G4Log(pLab) - minLogP;
877 G4double ld2 = ld*ld; 1111 G4double ld2 = ld*ld;
878 G4double sp = std::sqrt(pLab); 1112 G4double sp = std::sqrt(pLab);
879 G4double psp = pLab*sp; 1113 G4double psp = pLab*sp;
880 G4double p2 = pLab*pLab; 1114 G4double p2 = pLab*pLab;
881 G4double p4 = p2*p2; 1115 G4double p4 = p2*p2;
>> 1116
>> 1117 G4double lh = pLab - 0.98;
>> 1118 G4double hd = lh*lh + .045;
882 1119
883 G4double lh = pLab - 1.01; << 1120
884 G4double hd = lh*lh + .011; << 1121 fElasticXsc = 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) // + .004/md
885 fElasticXsc = 5.2/psp + (cofLogE*ld2 + << 1122 + .15/hd;
886 fTotalXsc = 14./psp + (1.1*cofLogT*ld << 1123 fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
>> 1124 // + .006/md + 0.01/hd1 + 0.02/hd2
>> 1125 + .60/hd;
887 } 1126 }
888 } 1127 }
889 else if( (theParticle == theKMinus) && neutr << 1128 else if( (theParticle == theKMinus || theParticle == theK0S) && neutron ) // Kmn/K0n /////////////////////////////
890 { 1129 {
891 if( pLab > pMax ) 1130 if( pLab > pMax )
892 { 1131 {
893 G4double ld = logP - minLogP; << 1132 G4double ld = G4Log(pLab) - minLogP;
894 G4double ld2 = ld*ld; 1133 G4double ld2 = ld*ld;
895 fElasticXsc = cofLogE*ld2 + 2.23; << 1134 fElasticXsc = cofLogE*ld2 + 2.23;
896 fTotalXsc = 1.1*cofLogT*ld2 + 19.7; << 1135 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
897 } 1136 }
898 else 1137 else
899 { << 1138 {
>> 1139
900 G4double lh = pLab - 0.98; 1140 G4double lh = pLab - 0.98;
901 G4double hd = lh*lh + .045; // vg ve << 1141 G4double hd = lh*lh + .045;
902 G4double sqrLogPlab = logP*logP; <<
903 1142
904 fElasticXsc = 5.0 + 8.1*G4Exp(-logP*1.8 << 1143 fElasticXsc = // 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md
905 - 1.3*logP + .15/hd; << 1144 5.0 + 8.1*std::pow(pLab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab + .15/hd;
906 fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2. << 1145 fTotalXsc = // 14./psp +
>> 1146 // (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
>> 1147 // WP 25.2 + 0. *std::pow(pLab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab
>> 1148 25.2 + 0.38*sqrLogPlab - 2.9*LogPlab
>> 1149 // + .006/md + 0.01/hd1+ 0.02/hd2
>> 1150 + 0.60/hd ;
907 } 1151 }
908 } 1152 }
909 else if( (theParticle == theKPlus) && proton << 1153 else if( (theParticle == theKPlus || theParticle == theK0L) && proton ) // Kpp/aKp //////////////////////
910 { 1154 {
911 // VI: modified low-energy part << 1155 if( pLab < pMin )
912 if( pLab < 0.631 ) <<
913 { 1156 {
914 fElasticXsc = fTotalXsc = 12.03; << 1157 G4double lr = pLab - .38;
>> 1158 G4double lm = pLab - 1.;
>> 1159 G4double md = lm*lm + .392;
>> 1160 fElasticXsc = .7/(lr*lr + .076) + 2./md;
>> 1161 fTotalXsc = // .7/(lr*lr + .076) +
>> 1162 2.6/md;
915 } 1163 }
916 else if( pLab > pMax ) 1164 else if( pLab > pMax )
917 { 1165 {
918 G4double ld = logP - minLogP; << 1166 G4double ld = G4Log(pLab) - minLogP;
919 G4double ld2 = ld*ld; 1167 G4double ld2 = ld*ld;
920 fElasticXsc = cofLogE*ld2 + 2.23; << 1168 fElasticXsc = cofLogE*ld2 + 2.23;
921 fTotalXsc = cofLogT*ld2 + 19.2; << 1169 fTotalXsc = cofLogT*ld2 + 19.2;
922 } 1170 }
923 else 1171 else
924 { 1172 {
925 G4double ld = logP - minLogP; << 1173 G4double ld = G4Log(pLab) - minLogP;
926 G4double ld2 = ld*ld; 1174 G4double ld2 = ld*ld;
927 G4double lr = pLab - .38; 1175 G4double lr = pLab - .38;
928 G4double LE = .7/(lr*lr + .076); 1176 G4double LE = .7/(lr*lr + .076);
929 G4double sp = std::sqrt(pLab); 1177 G4double sp = std::sqrt(pLab);
930 G4double p2 = pLab*pLab; 1178 G4double p2 = pLab*pLab;
931 G4double p4 = p2*p2; 1179 G4double p4 = p2*p2;
932 // VI: tuned elastic << 1180 G4double lm = pLab - 0.8;
933 fElasticXsc = LE + (cofLogE*ld2 + 2.23) << 1181 G4double md = lm*lm + .652;
934 + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652); << 1182 fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
935 fTotalXsc = LE + (cofLogT*ld2 + 19.5) << 1183 fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 7.6/md; // + LE;
936 + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392); <<
937 } 1184 }
938 } 1185 }
939 else if( (theParticle == theKPlus) && neutr << 1186 else if( (theParticle == theKPlus || theParticle == theK0L) && neutron ) // Kpn/aKn //////////////////////////////////
940 { 1187 {
941 if( pLab < pMin ) 1188 if( pLab < pMin )
942 { 1189 {
943 G4double lm = pLab - 0.94; 1190 G4double lm = pLab - 0.94;
944 G4double md = lm*lm + .392; 1191 G4double md = lm*lm + .392;
945 fElasticXsc = 2./md; 1192 fElasticXsc = 2./md;
946 fTotalXsc = 4.6/md; 1193 fTotalXsc = 4.6/md;
947 } 1194 }
948 else if( pLab > pMax ) 1195 else if( pLab > pMax )
949 { 1196 {
950 G4double ld = logP - minLogP; << 1197 G4double ld = G4Log(pLab) - minLogP;
951 G4double ld2 = ld*ld; 1198 G4double ld2 = ld*ld;
952 fElasticXsc = cofLogE*ld2 + 2.23; << 1199 fElasticXsc = cofLogE*ld2 + 2.23;
953 fTotalXsc = cofLogT*ld2 + 19.2; << 1200 fTotalXsc = cofLogT*ld2 + 19.2;
954 } 1201 }
955 else 1202 else
956 { 1203 {
957 G4double ld = logP - minLogP; << 1204 G4double ld = G4Log(pLab) - minLogP;
958 G4double ld2 = ld*ld; 1205 G4double ld2 = ld*ld;
959 G4double sp = std::sqrt(pLab); 1206 G4double sp = std::sqrt(pLab);
960 G4double p2 = pLab*pLab; 1207 G4double p2 = pLab*pLab;
961 G4double p4 = p2*p2; 1208 G4double p4 = p2*p2;
962 G4double lm = pLab - 0.94; << 1209 G4double lm = pLab - 0.8;
963 G4double md = lm*lm + .392; << 1210 G4double md = lm*lm + .652;
964 fElasticXsc = (cofLogE*ld2 + 2.23)/(1. 1211 fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
965 fTotalXsc = (cofLogT*ld2 + 19.5)/(1. << 1212 fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 7.6/md;
966 } 1213 }
967 } 1214 }
>> 1215 fTotalXsc *= millibarn; // parametrised in mb
>> 1216 fElasticXsc *= millibarn; // parametrised in mb
968 1217
969 fTotalXsc *= CLHEP::millibarn; << 1218 if( proton && aParticle->GetDefinition()->GetPDGCharge() > 0. )
970 fElasticXsc *= CLHEP::millibarn; <<
971 <<
972 if( proton && theParticle->GetPDGCharge() > <<
973 { 1219 {
974 G4double cB = G4NuclearRadii::CoulombFacto << 1220 G4double cB = GetCoulombBarrier(aParticle, nucleon);
975 fTotalXsc *= cB; 1221 fTotalXsc *= cB;
976 fElasticXsc *= cB; 1222 fElasticXsc *= cB;
977 } 1223 }
978 fElasticXsc = std::min(fElasticXsc, fTotalXs << 1224 fInelasticXsc = fTotalXsc - fElasticXsc;
979 fInelasticXsc = std::max(fTotalXsc - fElasti << 1225 if( fInelasticXsc < 0. ) fInelasticXsc = 0.;
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 1226
989 ////////////////////////////////////////////// << 1227 // G4cout<<fTotalXsc/millibarn<<"; "<<fElasticXsc/millibarn<<"; "<<fInelasticXsc/millibarn<<G4endl;
990 // <<
991 // Returns hyperon-nucleon cross-section using <<
992 1228
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; 1229 return fTotalXsc;
1050 } 1230 }
1051 1231
1052 ///////////////////////////////////////////// << 1232 /////////////////////////////////////////////////////////////////////////////////////
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 { <<
1101 coeff = 0.76; <<
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 // 1233 //
1115 // Returns hadron-nucleon cross-section based 1234 // Returns hadron-nucleon cross-section based on V. Uzjinsky parametrisation of
1116 // data from G4FTFCrossSection class 1235 // data from G4FTFCrossSection class
1117 1236
1118 G4double G4HadronNucleonXsc::HadronNucleonXsc << 1237 G4double
1119 const G4Particle << 1238 G4HadronNucleonXsc::GetHadronNucleonXscVU(const G4DynamicParticle* aParticle,
1120 const G4ParticleDefinition* nucleo << 1239 const G4ParticleDefinition* nucleon )
1121 { 1240 {
1122 G4int PDGcode = theParticle->GetPDGEncoding << 1241 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
1123 G4int absPDGcode = std::abs(PDGcode); 1242 G4int absPDGcode = std::abs(PDGcode);
1124 G4double mass = theParticle->GetPDGMass(); << 1243 G4double Elab = aParticle->GetTotalEnergy();
1125 G4double Plab = std::sqrt(ekin*(ekin + 2.*m << 1244 // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
>> 1245 G4double Plab = aParticle->GetMomentum().mag();
>> 1246 // std::sqrt(Elab * Elab - 0.88);
>> 1247
>> 1248 Elab /= GeV;
>> 1249 Plab /= GeV;
1126 1250
1127 G4double logPlab = G4Log( Plab ); << 1251 G4double LogPlab = G4Log( Plab );
1128 G4double sqrLogPlab = logPlab * logPlab; << 1252 G4double sqrLogPlab = LogPlab * LogPlab;
1129 1253
>> 1254 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
1130 G4bool proton = (nucleon == theProton); 1255 G4bool proton = (nucleon == theProton);
1131 G4bool neutron = (nucleon == theNeutron); 1256 G4bool neutron = (nucleon == theNeutron);
>> 1257
1132 1258
1133 if( absPDGcode > 1000) //------Projectile << 1259 if( absPDGcode > 1000 && pORn ) //------Projectile is baryon -
1134 { 1260 {
1135 if(proton) 1261 if(proton)
1136 { 1262 {
1137 fTotalXsc = 48.0 + 0.522*sqrLogPlab - << 1263 // WP fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
1138 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1264 fTotalXsc = 48.0 + 0.522*sqrLogPlab - 4.51*LogPlab;
>> 1265 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1139 } 1266 }
1140 if(neutron) 1267 if(neutron)
1141 { 1268 {
1142 fTotalXsc = 47.3 + 0.513*sqrLogPlab << 1269 // WP fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
1143 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1270 fTotalXsc = 47.3 + 0.513*sqrLogPlab - 4.27*LogPlab;
>> 1271 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1144 } 1272 }
1145 } 1273 }
1146 else if( PDGcode == 211) //------Projecti << 1274 else if( PDGcode == 211 && pORn ) //------Projectile is PionPlus ----
1147 { 1275 {
1148 if(proton) 1276 if(proton)
1149 { 1277 {
1150 fTotalXsc = 16.4 + 19.3 *G4Exp(-logPla << 1278 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
1151 fElasticXsc = 0.0 + 11.4*G4Exp(-logPla << 1279 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
1152 } 1280 }
1153 if(neutron) 1281 if(neutron)
1154 { 1282 {
1155 fTotalXsc = 33.0 + 14.0 *G4Exp(-logP << 1283 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
1156 fElasticXsc = 1.76 + 11.2*G4Exp(-logPla << 1284 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
1157 } 1285 }
1158 } 1286 }
1159 else if( PDGcode == -211) //------Projecti << 1287 else if( PDGcode == -211 && pORn ) //------Projectile is PionMinus ----
1160 { 1288 {
1161 if(proton) 1289 if(proton)
1162 { 1290 {
1163 fTotalXsc = 33.0 + 14.0 *G4Exp(-logPl << 1291 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
1164 fElasticXsc = 1.76 + 11.2*G4Exp(-logPla << 1292 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
1165 } 1293 }
1166 if(neutron) 1294 if(neutron)
1167 { 1295 {
1168 fTotalXsc = 16.4 + 19.3 *G4Exp(-logPl << 1296 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
1169 fElasticXsc = 0.0 + 11.4*G4Exp(-logPla << 1297 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
1170 } 1298 }
1171 } 1299 }
1172 else if( PDGcode == 111) //------Projecti << 1300 else if( PDGcode == 111 && pORn ) //------Projectile is PionZero --
1173 { 1301 {
1174 if(proton) 1302 if(proton)
1175 { 1303 {
1176 fTotalXsc = (16.4 + 19.3 *G4Exp(-logPla << 1304 fTotalXsc = (16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
1177 33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0. << 1305 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
1178 1306
1179 fElasticXsc = (0.0 + 11.4*G4Exp(-logPla << 1307 fElasticXsc = ( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
1180 1.76 + 11.2*G4Exp(-logPlab*0.64) + 0 << 1308 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1181 1309
1182 } 1310 }
1183 if(neutron) 1311 if(neutron)
1184 { 1312 {
1185 fTotalXsc = (33.0 + 14.0 *G4Exp(-logPla << 1313 fTotalXsc = (33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
1186 16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0. << 1314 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1187 fElasticXsc = (1.76 + 11.2*G4Exp(-logPl << 1315 fElasticXsc = ( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
1188 0.0 + 11.4*G4Exp(-logPlab*0.40) + 0 << 1316 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1189 } 1317 }
1190 } 1318 }
1191 else if( PDGcode == 321 ) //------Projec << 1319 else if( PDGcode == 321 && pORn ) //------Projectile is KaonPlus --
1192 { 1320 {
1193 if(proton) 1321 if(proton)
1194 { 1322 {
1195 fTotalXsc = 18.1 + 0.26 *sqrLogPlab << 1323 // WP fTotalXsc = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
1196 fElasticXsc = 5.0 + 8.1*G4Exp(-logPla << 1324 fTotalXsc = 18.1 + 0.26 *sqrLogPlab - 1.0 *LogPlab;
>> 1325 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
1197 } 1326 }
1198 if(neutron) 1327 if(neutron)
1199 { 1328 {
1200 fTotalXsc = 18.7 + 0.21 *sqrLogPlab << 1329 // WP fTotalXsc = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
1201 fElasticXsc = 7.3 + 0.29 *sqrLogPlab << 1330 // WP fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
>> 1331 fTotalXsc = 18.7 + 0.21 *sqrLogPlab - 0.89*LogPlab;
>> 1332 fElasticXsc = 7.3 + 0.29 *sqrLogPlab - 2.4 *LogPlab;
1202 } 1333 }
1203 } 1334 }
1204 else if( PDGcode ==-321 ) //------Projecti << 1335 else if( PDGcode ==-321 && pORn ) //------Projectile is KaonMinus ----
1205 { 1336 {
1206 if(proton) 1337 if(proton)
1207 { 1338 {
1208 fTotalXsc = 32.1 + 0.66*sqrLogPlab - << 1339 // WP fTotalXsc = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66*sqrLogPlab - 5.6*LogPlab;
1209 fElasticXsc = 7.3 + 0.29*sqrLogPlab - << 1340 // WP fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29*sqrLogPlab - 2.4*LogPlab;
>> 1341 fTotalXsc = 32.1 + 0.66*sqrLogPlab - 5.6*LogPlab;
>> 1342 fElasticXsc = 7.3 + 0.29*sqrLogPlab - 2.4*LogPlab;
1210 } 1343 }
1211 if(neutron) 1344 if(neutron)
1212 { 1345 {
1213 fTotalXsc = 25.2 + 0.38*sqrLogPlab - << 1346 // WP fTotalXsc = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab;
1214 fElasticXsc = 5.0 + 8.1*G4Exp(-logPla << 1347 fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2.9*LogPlab;
>> 1348 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab;
1215 } 1349 }
1216 } 1350 }
1217 else if( PDGcode == 311 ) //------Projecti << 1351 else if( PDGcode == 311 && pORn ) //------Projectile is KaonZero -----
1218 { 1352 {
1219 if(proton) 1353 if(proton)
1220 { 1354 {
1221 fTotalXsc = ( 18.1 + 0.26 *sqrLogPlab << 1355 // WP fTotalXsc = ( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
1222 32.1 + 0.66 *sqrLogPlab << 1356 // WP 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
1223 fElasticXsc = ( 5.0 + 8.1*G4Exp(-logP << 1357 fTotalXsc = ( 18.1 + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
1224 7.3 + 0.29 *sqrLogPl << 1358 32.1 + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
>> 1359 fElasticXsc = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
>> 1360 7.3 + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
>> 1361 // WP 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
1225 } 1362 }
1226 if(neutron) 1363 if(neutron)
1227 { 1364 {
1228 fTotalXsc = ( 18.7 + 0.21 *sqrLogPlab << 1365 // WP fTotalXsc = ( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
1229 25.2 + 0.38 *sqrLogPlab << 1366 // WP 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
1230 fElasticXsc = ( 7.3 + 0.29 *sqrLogPlab << 1367 fTotalXsc = ( 18.7 + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
1231 5.0 + 8.1*G4Exp(-logPl << 1368 25.2 + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
>> 1369 // WP fElasticXsc = ( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
>> 1370 fElasticXsc = ( 7.3 + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
>> 1371 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
1232 } 1372 }
1233 } 1373 }
1234 else //------Projectile is undefined, Nucl 1374 else //------Projectile is undefined, Nucleon assumed
1235 { 1375 {
1236 if(proton) 1376 if(proton)
1237 { 1377 {
1238 fTotalXsc = 48.0 + 0.522*sqrLogPlab - << 1378 // WP fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
1239 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1379 fTotalXsc = 48.0 + 0.522*sqrLogPlab - 4.51*LogPlab;
>> 1380 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1240 } 1381 }
1241 if(neutron) 1382 if(neutron)
1242 { 1383 {
1243 fTotalXsc = 47.3 + 0.513*sqrLogPlab - << 1384 // WP fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
1244 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla << 1385 fTotalXsc = 47.3 + 0.513*sqrLogPlab - 4.27*LogPlab;
>> 1386 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
1245 } 1387 }
1246 } 1388 }
1247 << 1389 fTotalXsc *= millibarn;
1248 fTotalXsc *= CLHEP::millibarn; << 1390 fElasticXsc *= millibarn;
1249 fElasticXsc *= CLHEP::millibarn; << 1391 fInelasticXsc = fTotalXsc - fElasticXsc;
1250 fElasticXsc = std::min(fElasticXsc, fTotalX << 1392 if (fInelasticXsc < 0.) fInelasticXsc = 0.;
1251 fInelasticXsc = fTotalXsc - fElasticXsc; <<
1252 1393
1253 return fTotalXsc; 1394 return fTotalXsc;
1254 } 1395 }
1255 1396
1256 ///////////////////////////////////////////// << 1397 ////////////////////////////////////////////////////////////////////////////////////
>> 1398 //
1257 // 1399 //
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 1400
1263 G4double G4HadronNucleonXsc::HadronNucleonXsc << 1401 G4double G4HadronNucleonXsc::CalculateEcmValue( const G4double mp ,
1264 const G4Particle << 1402 const G4double mt ,
1265 const G4ParticleDefinition*, G4dou << 1403 const G4double Plab )
1266 { 1404 {
1267 G4int pdg = theParticle->GetPDGEncoding(); << 1405 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1268 G4double xsection(0.); << 1406 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
1269 static const G4double targ_mass = << 1407 // G4double Pcm = Plab * mt / Ecm;
1270 0.5*(CLHEP::proton_mass_c2 + CLHEP::neutr << 1408 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
1271 G4double sMand = <<
1272 CalcMandelstamS(ekin, theParticle->GetPDG <<
1273 1409
1274 G4double x1 = G4Exp(G4Log(sMand)*0.0808); << 1410 return Ecm ; // KEcm;
1275 G4double x2 = G4Exp(G4Log(-sMand)*0.4525); << 1411 }
1276 1412
1277 if(pdg == 22) << 1413
1278 { << 1414 ////////////////////////////////////////////////////////////////////////////////////
1279 xsection = 0.0677*x1 + 0.129*x2; << 1415 //
1280 } << 1416 //
1281 else if(theParticle == theNeutron) << 1417
1282 { << 1418 G4double G4HadronNucleonXsc::CalcMandelstamS( const G4double mp ,
1283 xsection = 21.70*x1 + 56.08*x2; << 1419 const G4double mt ,
1284 } << 1420 const G4double Plab )
1285 else if(theParticle == theProton) << 1421 {
1286 { << 1422 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1287 xsection = 21.70*x1 + 56.08*x2; << 1423 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1288 } << 1424
1289 // pbar << 1425 return sMand;
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; <<
1323 } 1426 }
1324 1427
>> 1428
1325 ///////////////////////////////////////////// 1429 ///////////////////////////////////////////////////////////////////////////////
>> 1430 //
>> 1431 //
1326 1432
1327 G4double << 1433 G4double G4HadronNucleonXsc::GetCoulombBarrier(const G4DynamicParticle* aParticle,
1328 G4HadronNucleonXsc::CoulombBarrier(const G4Pa << 1434 const G4ParticleDefinition* nucleon )
1329 const G4ParticleDefinition* nucleo <<
1330 G4double ekin) <<
1331 { 1435 {
1332 G4double tR = 0.895*CLHEP::fermi; << 1436 G4double ratio;
1333 G4double pR = 0.5*CLHEP::fermi; << 1437
>> 1438 G4double tR = 0.895*fermi, pR;
1334 1439
1335 if ( theParticle == theProton ) pR = 0. << 1440 if ( aParticle->GetDefinition() == theProton ) pR = 0.895*fermi;
1336 else if( theParticle == thePiPlus ) pR = 0. << 1441 else if( aParticle->GetDefinition() == thePiPlus ) pR = 0.663*fermi;
1337 else if( theParticle == theKPlus ) pR = 0. << 1442 else if( aParticle->GetDefinition() == theKPlus ) pR = 0.340*fermi;
>> 1443 else pR = 0.500*fermi;
1338 1444
1339 G4double pZ = theParticle->GetPDGCharge(); << 1445 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
1340 G4double tZ = nucleon->GetPDGCharge(); 1446 G4double tZ = nucleon->GetPDGCharge();
1341 1447
1342 G4double pM = theParticle->GetPDGMass(); << 1448 G4double pTkin = aParticle->GetKineticEnergy();
1343 G4double tM = nucleon->GetPDGMass(); << 1449
>> 1450 G4double pM = aParticle->GetDefinition()->GetPDGMass();
>> 1451 G4double tM = nucleon->GetPDGMass();
1344 1452
1345 G4double pElab = ekin + pM; << 1453 G4double pElab = pTkin + pM;
1346 1454
1347 G4double totEcm = std::sqrt(pM*pM + tM*tM 1455 G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
1348 1456
1349 G4double totTcm = totEcm - pM -tM; 1457 G4double totTcm = totEcm - pM -tM;
1350 1458
1351 G4double bC = fine_structure_const*hbarc*pZ << 1459 G4double bC = fine_structure_const*hbarc*pZ*tZ;
>> 1460 bC /= pR + tR;
>> 1461 bC /= 2.; // 4., 2. parametrisation cof ??? vmg
>> 1462
>> 1463 // G4cout<<"pTkin = "<<pTkin/GeV<<"; pPlab = "
>> 1464 // <<pPlab/GeV<<"; bC = "<<bC/GeV<<"; pTcm = "<<pTcm/GeV<<G4endl;
1352 1465
1353 //G4cout<<"##CB: ekin = "<<ekin<<"; pElab= << 1466 if( totTcm <= bC ) ratio = 0.;
1354 // <<"; bC = "<<bC<<"; Tcm = "<< totTcm << << 1467 else ratio = 1. - bC/totTcm;
1355 1468
1356 G4double ratio = (totTcm > bC) ? 1. - bC/to << 1469 // if(ratio < DBL_MIN) ratio = DBL_MIN;
>> 1470 if( ratio < 0.) ratio = 0.;
1357 1471
1358 // G4cout <<"ratio = "<<ratio<<G4endl; 1472 // G4cout <<"ratio = "<<ratio<<G4endl;
1359 return ratio; 1473 return ratio;
1360 } 1474 }
1361 1475
1362 ///////////////////////////////////////////// << 1476
>> 1477
>> 1478 /*
>> 1479
>> 1480 ////////////////////////////////////////////////////////////////////////////////////
>> 1481 //
>> 1482 // Initialaise K(p,m)-(p,n) total cross section vectors
>> 1483
>> 1484
>> 1485 void G4HadronNucleonXsc::InitialiseKaonNucleonTotXsc()
>> 1486 {
>> 1487 G4int i = 0, iMax;
>> 1488 G4double tmpxsc[106];
>> 1489
>> 1490 // Kp-proton tot xsc
>> 1491
>> 1492 iMax = 66;
>> 1493 fKpProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKpProtonTotTkin[0], fKpProtonTotTkin[iMax-1]);
>> 1494 fKpProtonTotXscVector.SetSpline(true);
>> 1495
>> 1496 for( i = 0; i < iMax; i++)
>> 1497 {
>> 1498 tmpxsc[i] = 0.;
>> 1499
>> 1500 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpProtonTotXsc[i];
>> 1501 else tmpxsc[i] = 0.5*(fKpProtonTotXsc[i-1]+fKpProtonTotXsc[i+1]);
>> 1502
>> 1503 fKpProtonTotXscVector.PutValues(size_t(i), fKpProtonTotTkin[i], tmpxsc[i]*millibarn);
>> 1504 }
>> 1505
>> 1506 // Kp-neutron tot xsc
>> 1507
>> 1508 iMax = 75;
>> 1509 fKpNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKpNeutronTotTkin[0], fKpNeutronTotTkin[iMax-1]);
>> 1510 fKpNeutronTotXscVector.SetSpline(true);
>> 1511
>> 1512 for( i = 0; i < iMax; i++)
>> 1513 {
>> 1514 tmpxsc[i] = 0.;
>> 1515 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKpNeutronTotXsc[i];
>> 1516 else tmpxsc[i] = 0.5*(fKpNeutronTotXsc[i-1]+fKpNeutronTotXsc[i+1]);
>> 1517
>> 1518 fKpNeutronTotXscVector.PutValues(size_t(i), fKpNeutronTotTkin[i], tmpxsc[i]*millibarn);
>> 1519 }
>> 1520
>> 1521 // Km-proton tot xsc
>> 1522
>> 1523 iMax = 106;
>> 1524 fKmProtonTotXscVector = G4LPhysicsFreeVector(iMax, fKmProtonTotTkin[0], fKmProtonTotTkin[iMax-1]);
>> 1525 fKmProtonTotXscVector.SetSpline(true);
>> 1526
>> 1527 for( i = 0; i < iMax; i++)
>> 1528 {
>> 1529 tmpxsc[i] = 0.;
>> 1530
>> 1531 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmProtonTotXsc[i];
>> 1532 else tmpxsc[i] = 0.5*(fKmProtonTotXsc[i-1]+fKmProtonTotXsc[i+1]);
>> 1533
>> 1534 fKmProtonTotXscVector.PutValues(size_t(i), fKmProtonTotTkin[i], tmpxsc[i]*millibarn);
>> 1535 }
>> 1536
>> 1537 // Km-neutron tot xsc
>> 1538
>> 1539 iMax = 68;
>> 1540 fKmNeutronTotXscVector = G4LPhysicsFreeVector(iMax, fKmNeutronTotTkin[0], fKmNeutronTotTkin[iMax-1]);
>> 1541 fKmNeutronTotXscVector.SetSpline(true);
>> 1542
>> 1543 for( i = 0; i < iMax; i++)
>> 1544 {
>> 1545 tmpxsc[i] = 0.;
>> 1546 if( i == 0 || i == iMax-1 ) tmpxsc[i] = fKmNeutronTotXsc[i];
>> 1547 else tmpxsc[i] = 0.5*(fKmNeutronTotXsc[i-1]+fKmNeutronTotXsc[i+1]);
>> 1548
>> 1549 fKmNeutronTotXscVector.PutValues(size_t(i), fKmNeutronTotTkin[i], tmpxsc[i]*millibarn);
>> 1550 }
>> 1551 }
>> 1552
>> 1553 ///////////////////////////////////////////////////////
>> 1554 //
>> 1555 // K-nucleon tot xsc (mb) fit data, G4Log(Tkin(MeV))
>> 1556
>> 1557 const G4double G4HadronNucleonXsc::fKpProtonTotXsc[66] = {
>> 1558 0.000000e+00, 1.592400e-01, 3.184700e-01, 7.961800e-01, 1.433120e+00, 2.070060e+00,
>> 1559 2.866240e+00, 3.582800e+00, 4.378980e+00, 5.015920e+00, 5.573250e+00, 6.449040e+00,
>> 1560 7.404460e+00, 8.200640e+00, 8.837580e+00, 9.713380e+00, 1.027070e+01, 1.090764e+01,
>> 1561 1.130573e+01, 1.170382e+01, 1.242038e+01, 1.281847e+01, 1.321656e+01, 1.337580e+01,
>> 1562 1.345541e+01, 1.329618e+01, 1.265924e+01, 1.242038e+01, 1.250000e+01, 1.305732e+01,
>> 1563 1.369427e+01, 1.425159e+01, 1.544586e+01, 1.648089e+01, 1.751592e+01, 1.791401e+01,
>> 1564 1.791401e+01, 1.775478e+01, 1.751592e+01, 1.735669e+01, 1.719745e+01, 1.711783e+01,
>> 1565 1.703822e+01, 1.695860e+01, 1.695860e+01, 1.695860e+01, 1.695860e+01, 1.687898e+01,
>> 1566 1.687898e+01, 1.703822e+01, 1.719745e+01, 1.735669e+01, 1.751592e+01, 1.767516e+01,
>> 1567 1.783439e+01, 1.799363e+01, 1.815287e+01, 1.839172e+01, 1.855095e+01, 1.871019e+01,
>> 1568 1.886943e+01, 1.918790e+01, 1.942675e+01, 1.966561e+01, 2.006369e+01, 2.054140e+01
>> 1569 }; // 66
>> 1570
>> 1571
>> 1572 const G4double G4HadronNucleonXsc::fKpProtonTotTkin[66] = {
>> 1573 2.100000e+00, 2.180770e+00, 2.261540e+00, 2.396150e+00, 2.476920e+00, 2.557690e+00,
>> 1574 2.557690e+00, 2.584620e+00, 2.638460e+00, 2.665380e+00, 2.719230e+00, 2.746150e+00,
>> 1575 2.800000e+00, 2.853850e+00, 2.934620e+00, 3.042310e+00, 3.150000e+00, 3.311540e+00,
>> 1576 3.446150e+00, 3.607690e+00, 3.930770e+00, 4.226920e+00, 4.361540e+00, 4.846150e+00,
>> 1577 4.980770e+00, 5.088460e+00, 5.465380e+00, 5.653850e+00, 5.950000e+00, 6.084620e+00,
>> 1578 6.246150e+00, 6.300000e+00, 6.380770e+00, 6.515380e+00, 6.730770e+00, 6.838460e+00,
>> 1579 7.000000e+00, 7.161540e+00, 7.323080e+00, 7.457690e+00, 7.619230e+00, 7.780770e+00,
>> 1580 7.915380e+00, 8.130770e+00, 8.265380e+00, 8.453850e+00, 8.642310e+00, 8.803850e+00,
>> 1581 9.019230e+00, 9.234620e+00, 9.530770e+00, 9.773080e+00, 1.001538e+01, 1.017692e+01,
>> 1582 1.033846e+01, 1.058077e+01, 1.082308e+01, 1.098462e+01, 1.114615e+01, 1.138846e+01,
>> 1583 1.160385e+01, 1.173846e+01, 1.192692e+01, 1.216923e+01, 1.238461e+01, 1.257308e+01
>> 1584 }; // 66
>> 1585
>> 1586 const G4double G4HadronNucleonXsc::fKpNeutronTotXsc[75] = {
>> 1587 3.980900e-01, 3.184700e-01, 3.184700e-01, 3.980900e-01, 3.980900e-01, 3.980900e-01,
>> 1588 3.980900e-01, 3.980900e-01, 3.980900e-01, 4.777100e-01, 3.980900e-01, 3.980900e-01,
>> 1589 4.777100e-01, 5.573200e-01, 1.035030e+00, 1.512740e+00, 2.149680e+00, 2.786620e+00,
>> 1590 3.503180e+00, 4.219750e+00, 5.015920e+00, 5.652870e+00, 6.289810e+00, 7.245220e+00,
>> 1591 8.121020e+00, 8.837580e+00, 9.633760e+00, 1.042994e+01, 1.114650e+01, 1.194268e+01,
>> 1592 1.265924e+01, 1.329618e+01, 1.393312e+01, 1.449045e+01, 1.496815e+01, 1.552548e+01,
>> 1593 1.592357e+01, 1.664013e+01, 1.727707e+01, 1.783439e+01, 1.831210e+01, 1.902866e+01,
>> 1594 1.902866e+01, 1.878981e+01, 1.847134e+01, 1.831210e+01, 1.807325e+01, 1.791401e+01,
>> 1595 1.783439e+01, 1.767516e+01, 1.759554e+01, 1.743631e+01, 1.743631e+01, 1.751592e+01,
>> 1596 1.743631e+01, 1.735669e+01, 1.751592e+01, 1.759554e+01, 1.767516e+01, 1.783439e+01,
>> 1597 1.783439e+01, 1.791401e+01, 1.815287e+01, 1.823248e+01, 1.847134e+01, 1.878981e+01,
>> 1598 1.894905e+01, 1.902866e+01, 1.934713e+01, 1.966561e+01, 1.990446e+01, 2.014331e+01,
>> 1599 2.030255e+01, 2.046178e+01, 2.085987e+01
>> 1600 }; // 75
>> 1601
>> 1602 const G4double G4HadronNucleonXsc::fKpNeutronTotTkin[75] = {
>> 1603 2.692000e-02, 1.615400e-01, 2.961500e-01, 4.576900e-01, 6.461500e-01, 7.538500e-01,
>> 1604 8.884600e-01, 1.103850e+00, 1.211540e+00, 1.400000e+00, 1.561540e+00, 1.776920e+00,
>> 1605 1.992310e+00, 2.126920e+00, 2.342310e+00, 2.423080e+00, 2.557690e+00, 2.692310e+00,
>> 1606 2.800000e+00, 2.988460e+00, 3.203850e+00, 3.365380e+00, 3.500000e+00, 3.688460e+00,
>> 1607 3.850000e+00, 4.011540e+00, 4.173080e+00, 4.415380e+00, 4.630770e+00, 4.873080e+00,
>> 1608 5.061540e+00, 5.276920e+00, 5.492310e+00, 5.707690e+00, 5.896150e+00, 6.030770e+00,
>> 1609 6.138460e+00, 6.219230e+00, 6.273080e+00, 6.326920e+00, 6.407690e+00, 6.650000e+00,
>> 1610 6.784620e+00, 7.026920e+00, 7.242310e+00, 7.350000e+00, 7.484620e+00, 7.619230e+00,
>> 1611 7.807690e+00, 7.915380e+00, 8.050000e+00, 8.211540e+00, 8.453850e+00, 8.588460e+00,
>> 1612 8.830770e+00, 9.073080e+00, 9.288460e+00, 9.476920e+00, 9.665380e+00, 9.826920e+00,
>> 1613 1.004231e+01, 1.031154e+01, 1.052692e+01, 1.071538e+01, 1.095769e+01, 1.120000e+01,
>> 1614 1.138846e+01, 1.155000e+01, 1.176538e+01, 1.190000e+01, 1.214231e+01, 1.222308e+01,
>> 1615 1.238461e+01, 1.246538e+01, 1.265385e+01
>> 1616 }; // 75
>> 1617
>> 1618 const G4double G4HadronNucleonXsc::fKmProtonTotXsc[106] = {
>> 1619 1.136585e+02, 9.749129e+01, 9.275262e+01, 8.885017e+01, 8.334146e+01, 7.955401e+01,
>> 1620 7.504530e+01, 7.153658e+01, 6.858537e+01, 6.674913e+01, 6.525784e+01, 6.448781e+01,
>> 1621 6.360279e+01, 6.255401e+01, 6.127526e+01, 6.032404e+01, 5.997910e+01, 5.443554e+01,
>> 1622 5.376307e+01, 5.236934e+01, 5.113937e+01, 5.090941e+01, 4.967944e+01, 4.844948e+01,
>> 1623 4.705575e+01, 4.638327e+01, 4.571080e+01, 4.475958e+01, 4.397213e+01, 4.257840e+01,
>> 1624 4.102090e+01, 4.090592e+01, 3.906969e+01, 3.839721e+01, 3.756097e+01, 3.644599e+01,
>> 1625 3.560976e+01, 3.533101e+01, 3.533101e+01, 3.644599e+01, 3.811847e+01, 3.839721e+01,
>> 1626 3.979094e+01, 4.090592e+01, 4.257840e+01, 4.341463e+01, 4.425087e+01, 4.564460e+01,
>> 1627 4.759582e+01, 4.703833e+01, 4.843206e+01, 4.787457e+01, 4.452962e+01, 4.202090e+01,
>> 1628 4.034843e+01, 3.839721e+01, 3.616725e+01, 3.365854e+01, 3.170732e+01, 3.087108e+01,
>> 1629 3.170732e+01, 3.254355e+01, 3.310104e+01, 3.254355e+01, 3.142857e+01, 3.059233e+01,
>> 1630 2.947735e+01, 2.891986e+01, 2.836237e+01, 2.752613e+01, 2.696864e+01, 2.641115e+01,
>> 1631 2.501742e+01, 2.473868e+01, 2.418118e+01, 2.362369e+01, 2.334495e+01, 2.278746e+01,
>> 1632 2.250871e+01, 2.222997e+01, 2.167247e+01, 2.139373e+01, 2.139373e+01, 2.139373e+01,
>> 1633 2.111498e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01,
>> 1634 2.083624e+01, 2.055749e+01, 2.055749e+01, 2.055749e+01, 2.027875e+01, 2.000000e+01,
>> 1635 2.055749e+01, 2.027875e+01, 2.083624e+01, 2.083624e+01, 2.055749e+01, 2.083624e+01,
>> 1636 2.083624e+01, 2.083624e+01, 2.139373e+01, 2.139373e+01
>> 1637 }; // 106
>> 1638
>> 1639 const G4double G4HadronNucleonXsc::fKmProtonTotTkin[106] = {
>> 1640 4.017980e+00, 4.125840e+00, 4.179780e+00, 4.251690e+00, 4.287640e+00, 4.341570e+00,
>> 1641 4.395510e+00, 4.467420e+00, 4.503370e+00, 4.575280e+00, 4.683150e+00, 4.737080e+00,
>> 1642 4.773030e+00, 4.826970e+00, 4.880900e+00, 4.916850e+00, 4.952810e+00, 4.988760e+00,
>> 1643 4.988761e+00, 5.006740e+00, 5.006741e+00, 5.042700e+00, 5.078650e+00, 5.114610e+00,
>> 1644 5.132580e+00, 5.150560e+00, 5.186520e+00, 5.204490e+00, 5.276400e+00, 5.348310e+00,
>> 1645 5.366290e+00, 5.384270e+00, 5.456180e+00, 5.564040e+00, 5.600000e+00, 5.671910e+00,
>> 1646 5.743820e+00, 5.833710e+00, 5.905620e+00, 5.977530e+00, 6.085390e+00, 6.085390e+00,
>> 1647 6.157300e+00, 6.175280e+00, 6.211240e+00, 6.229210e+00, 6.247190e+00, 6.337080e+00,
>> 1648 6.391010e+00, 6.516850e+00, 6.462920e+00, 6.498880e+00, 6.570790e+00, 6.606740e+00,
>> 1649 6.660670e+00, 6.678650e+00, 6.696630e+00, 6.732580e+00, 6.804490e+00, 6.876400e+00,
>> 1650 6.948310e+00, 7.020220e+00, 7.074160e+00, 7.182020e+00, 7.235960e+00, 7.289890e+00,
>> 1651 7.397750e+00, 7.523600e+00, 7.631460e+00, 7.757300e+00, 7.901120e+00, 8.062920e+00,
>> 1652 8.260670e+00, 8.386520e+00, 8.530340e+00, 8.674160e+00, 8.817980e+00, 8.943820e+00,
>> 1653 9.087640e+00, 9.267420e+00, 9.429210e+00, 9.573030e+00, 9.698880e+00, 9.896630e+00,
>> 1654 1.002247e+01, 1.016629e+01, 1.031011e+01, 1.048989e+01, 1.063371e+01, 1.077753e+01,
>> 1655 1.095730e+01, 1.108315e+01, 1.120899e+01, 1.135281e+01, 1.149663e+01, 1.162247e+01,
>> 1656 1.174831e+01, 1.187416e+01, 1.200000e+01, 1.212584e+01, 1.221573e+01, 1.234157e+01,
>> 1657 1.239551e+01, 1.250337e+01, 1.261124e+01, 1.273708e+01
>> 1658 }; // 106
>> 1659
>> 1660 const G4double G4HadronNucleonXsc::fKmNeutronTotXsc[68] = {
>> 1661 2.621810e+01, 2.741123e+01, 2.868413e+01, 2.963889e+01, 3.067343e+01, 3.178759e+01,
>> 1662 3.282148e+01, 3.417466e+01, 3.536778e+01, 3.552620e+01, 3.544576e+01, 3.496756e+01,
>> 1663 3.433030e+01, 3.401166e+01, 3.313537e+01, 3.257772e+01, 3.178105e+01, 3.138264e+01,
>> 1664 3.074553e+01, 2.970952e+01, 2.891301e+01, 2.827542e+01, 2.787700e+01, 2.715978e+01,
>> 1665 2.660181e+01, 2.612394e+01, 2.564574e+01, 2.516721e+01, 2.421098e+01, 2.365235e+01,
>> 1666 2.317366e+01, 2.261437e+01, 2.237389e+01, 2.205427e+01, 2.181395e+01, 2.165357e+01,
>> 1667 2.149335e+01, 2.133297e+01, 2.109232e+01, 2.093128e+01, 2.069030e+01, 2.052992e+01,
>> 1668 2.028927e+01, 2.012824e+01, 1.996737e+01, 1.996590e+01, 1.988530e+01, 1.964432e+01,
>> 1669 1.948361e+01, 1.940236e+01, 1.940040e+01, 1.931882e+01, 1.947593e+01, 1.947429e+01,
>> 1670 1.939320e+01, 1.939157e+01, 1.946922e+01, 1.962715e+01, 1.970481e+01, 1.970301e+01,
>> 1671 1.993958e+01, 2.009669e+01, 2.025380e+01, 2.033178e+01, 2.049003e+01, 2.064747e+01,
>> 1672 2.080540e+01, 2.096333e+01
>> 1673 }; // 68
>> 1674
>> 1675 const G4double G4HadronNucleonXsc::fKmNeutronTotTkin[68] = {
>> 1676 5.708500e+00, 5.809560e+00, 5.896270e+00, 5.954120e+00, 5.997630e+00, 6.041160e+00,
>> 1677 6.142160e+00, 6.171410e+00, 6.272470e+00, 6.344390e+00, 6.416230e+00, 6.459180e+00,
>> 1678 6.487690e+00, 6.501940e+00, 6.544740e+00, 6.573280e+00, 6.616110e+00, 6.644710e+00,
>> 1679 6.658840e+00, 6.744700e+00, 6.773150e+00, 6.830410e+00, 6.859010e+00, 6.916240e+00,
>> 1680 6.973530e+00, 6.987730e+00, 7.030670e+00, 7.102360e+00, 7.173880e+00, 7.288660e+00,
>> 1681 7.374720e+00, 7.547000e+00, 7.690650e+00, 7.791150e+00, 7.920420e+00, 8.020980e+00,
>> 1682 8.107160e+00, 8.207720e+00, 8.365740e+00, 8.523790e+00, 8.710560e+00, 8.811110e+00,
>> 1683 8.969140e+00, 9.127190e+00, 9.270860e+00, 9.400230e+00, 9.486440e+00, 9.673210e+00,
>> 1684 9.802510e+00, 9.946220e+00, 1.011870e+01, 1.029116e+01, 1.047808e+01, 1.062181e+01,
>> 1685 1.075114e+01, 1.089488e+01, 1.106739e+01, 1.118244e+01, 1.135496e+01, 1.151307e+01,
>> 1686 1.171439e+01, 1.190130e+01, 1.208822e+01, 1.223199e+01, 1.231829e+01, 1.247646e+01,
>> 1687 1.259150e+01, 1.270655e+01
>> 1688 }; // 68
>> 1689
>> 1690
>> 1691 */
>> 1692
>> 1693
>> 1694 //
>> 1695 //
>> 1696 ///////////////////////////////////////////////////////////////////////////////////////
1363 1697