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