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