| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // 14.03.07 V. Grichine - first implementation
27 //
28 // 04.09.18 V. Ivantchenko Major revision of i
29 // 30.09.18 V. Grichine hyperon-nucleon xsc fi
30
31 #include "G4HadronNucleonXsc.hh"
32 #include "G4Element.hh"
33 #include "G4Proton.hh"
34 #include "G4Nucleus.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4ParticleTable.hh"
38 #include "G4IonTable.hh"
39 #include "G4HadTmpUtil.hh"
40 #include "G4Log.hh"
41 #include "G4Exp.hh"
42 #include "G4Pow.hh"
43 #include "G4NuclearRadii.hh"
44
45 #include "G4Neutron.hh"
46 #include "G4PionPlus.hh"
47 #include "G4KaonPlus.hh"
48 #include "G4KaonMinus.hh"
49 #include "G4KaonZeroShort.hh"
50 #include "G4KaonZeroLong.hh"
51
52 namespace
53 {
54 const G4double invGeV = 1.0/CLHEP::GeV;
55 const G4double invGeV2 = 1.0/(CLHEP::GeV*CLH
56 // PDG fit constants
57 const G4double minLogP = 3.5; // min of (
58 const G4double cofLogE = .0557; // elastic
59 const G4double cofLogT = .3; // total (l
60 const G4double pMin = .1; // fast LE
61 const G4double pMax = 1000.; // fast HE
62 const G4double ekinmin = 0.1*CLHEP::MeV; /
63 const G4double ekinmaxQB = 100*CLHEP::MeV; /
64 }
65
66 G4HadronNucleonXsc::G4HadronNucleonXsc()
67 {
68 // basic hadrons
69 theProton = G4Proton::Proton();
70 theNeutron = G4Neutron::Neutron();
71 thePiPlus = G4PionPlus::PionPlus();
72
73 // basic strange mesons
74 theKPlus = G4KaonPlus::KaonPlus();
75 theKMinus = G4KaonMinus::KaonMinus();
76 theK0S = G4KaonZeroShort::KaonZeroShort
77 theK0L = G4KaonZeroLong::KaonZeroLong()
78
79 g4calc = G4Pow::GetInstance();
80 }
81
82 void G4HadronNucleonXsc::CrossSectionDescripti
83 {
84 outFile << "G4HadronNucleonXsc calculates th
85 << "hadron-nucleon cross sections us
86 << "parameterizations within the Gla
87 << "valid for all incident gammas an
88 << "energies above 30 keV. This is
89 << "is to be used to build a cross s
90 }
91
92 G4double G4HadronNucleonXsc::HadronNucleonXsc(
93 const G4Particl
94 {
95 G4double xsc(0.);
96 G4int pdg = std::abs(theParticle->GetPDGEnco
97
98 // p, n, pi+-, pbar, nbar
99 if ( pdg == 2212 || pdg == 2112 || pdg == 21
100 xsc = HadronNucleonXscNS(theParticle, nucl
101 }
102 else if ( pdg == 22 ) // gamma
103 {
104 xsc = HadronNucleonXscPDG(theParticle, nuc
105 }
106 else if ( pdg == 321 || pdg == 310 || pdg ==
107 {
108 xsc = KaonNucleonXscNS(theParticle, nucleo
109 }
110 else if (pdg > 3000)
111 {
112 if (pdg == 3122 || pdg == 3222 || pdg == 3
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 }
136 return xsc;
137 }
138
139 //////////////////////////////////////////////
140 //
141 // Returns hadron-nucleon Xsc according to PDG
142 // http://pdg.lbl.gov/2017/reviews/hadronicrpp
143
144 G4double G4HadronNucleonXsc::HadronNucleonXscP
145 const G4ParticleD
146 const G4ParticleDefinition* nucleon
147 {
148 static const G4double M = 2.1206; // in G
149 static const G4double eta1 = 0.4473;
150 static const G4double eta2 = 0.5486;
151 static const G4double H = 0.272;
152
153 G4int pdg = theParticle->GetPDGEncoding();
154
155 G4double mass1 = (pdg == 22) ? 770. : thePar
156 G4double mass2 = nucleon->GetPDGMass();
157
158 G4double sMand = CalcMandelstamS(ekin, mass1
159 G4double x = (mass1 + mass2)*invGeV + M;
160 G4double blog = G4Log(sMand/(x*x));
161
162 G4double P(0.0), R1(0.0), R2(0.0), del(1.0);
163
164 G4bool proton = (nucleon == theProton);
165 G4bool neutron = (nucleon == theNeutron);
166
167 if(theParticle == theNeutron)
168 {
169 if ( proton )
170 {
171 P = 34.71;
172 R1 = 12.52;
173 R2 = -6.66;
174 }
175 else
176 {
177 P = 34.41;
178 R1 = 13.07;
179 R2 = -7.394;
180 }
181 }
182 else if(theParticle == theProton)
183 {
184 if ( neutron )
185 {
186 P = 34.71;
187 R1 = 12.52;
188 R2 = -6.66;
189 }
190 else
191 {
192 P = 34.41;
193 R1 = 13.07;
194 R2 = -7.394;
195 }
196 }
197 // pbar
198 else if(pdg == -2212)
199 {
200 if ( neutron )
201 {
202 P = 34.71;
203 R1 = 12.52;
204 R2 = 6.66;
205 }
206 else
207 {
208 P = 34.41;
209 R1 = 13.07;
210 R2 = 7.394;
211 }
212 }
213 // nbar
214 else if(pdg == -2112)
215 {
216 if ( proton )
217 {
218 P = 34.71;
219 R1 = 12.52;
220 R2 = 6.66;
221 }
222 else
223 {
224 P = 34.41;
225 R1 = 13.07;
226 R2 = 7.394;
227 }
228 }
229 // pi+
230 else if(pdg == 211)
231 {
232 P = 18.75;
233 R1 = 9.56;
234 R2 = -1.767;
235 }
236 // pi-
237 else if(pdg == -211)
238 {
239 P = 18.75;
240 R1 = 9.56;
241 R2 = 1.767;
242 }
243 else if(theParticle == theKPlus)
244 {
245 if ( proton )
246 {
247 P = 16.36;
248 R1 = 4.29;
249 R2 = -3.408;
250 }
251 else
252 {
253 P = 16.31;
254 R1 = 3.7;
255 R2 = -1.826;
256 }
257 }
258 else if(theParticle == theKMinus)
259 {
260 if ( proton )
261 {
262 P = 16.36;
263 R1 = 4.29;
264 R2 = 3.408;
265 }
266 else
267 {
268 P = 16.31;
269 R1 = 3.7;
270 R2 = 1.826;
271 }
272 }
273 else if(theParticle == theK0S || theParticle
274 {
275 P = 16.36;
276 R1 = 2.5;
277 R2 = 0.;
278 }
279 // sigma-
280 else if(pdg == 3112)
281 {
282 P = 34.7;
283 R1 = -46.;
284 R2 = 48.;
285 }
286 // gamma
287 else if(pdg == 22) // modify later on
288 {
289 del= 0.003063;
290 P = 34.71*del;
291 R1 = (neutron) ? 0.0231 : 0.0139;
292 R2 = 0.;
293 }
294 else // as proton ???
295 {
296 if ( neutron )
297 {
298 P = 34.71;
299 R1 = 12.52;
300 R2 = -6.66;
301 }
302 else
303 {
304 P = 34.41;
305 R1 = 13.07;
306 R2 = -7.394;
307 }
308 }
309 fTotalXsc = CLHEP::millibarn*
310 (del*(H*blog*blog + P) + R1*G4Exp(-eta1*bl
311 fInelasticXsc = 0.75*fTotalXsc;
312 fElasticXsc = fTotalXsc - fInelasticXsc;
313
314 if( proton && theParticle->GetPDGCharge() >
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 }
330
331 //////////////////////////////////////////////
332 //
333 // Returns hadron-nucleon cross-section based
334 // data from mainly http://wwwppds.ihep.su:800
335
336 G4double G4HadronNucleonXsc::HadronNucleonXscN
337 const G4ParticleD
338 const G4ParticleDefinition* nucleon
339 {
340 const G4double ekin = std::max(ekin0, ekinmi
341 G4int pdg = theParticle->GetPDGEncoding();
342 /*
343 G4cout<< "HadronNucleonXscNS: Ekin(GeV)= " <
344 << theParticle->GetParticleName() << "
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
369 G4bool proton = (nucleon == theProton);
370 G4bool neutron = (nucleon == theNeutron);
371
372 if( theParticle == theNeutron)
373 {
374 if( pLab >= 373.)
375 {
376 fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4
377 + 9.19*G4Exp(-G4Log(sMand)*0.458);
378 }
379 else if( pLab >= 100.)
380 {
381 fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G
382 + 9.19*G4Exp(-G4Log(sMand)*0.458);
383 }
384 else if( pLab >= 10.)
385 {
386 fElasticXsc = 6 + 20/( (logP-0.182)*(lo
387 }
388 else // pLab < 10 GeV/c
389 {
390 if( neutron ) // nn to be pp
391 {
392 G4double x = G4Log(pLab/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 {
400 fTotalXsc = 23 + 50*std::sqrt(g4calc
401 fElasticXsc = fTotalXsc;
402 }
403 else if( pLab < 1.05 )
404 {
405 fTotalXsc = 23 + 40*x*x;
406 fElasticXsc = 23 + 20*x*x;
407 }
408 else // 1.05 - 10 GeV/c
409 {
410 fTotalXsc = 39.0+75*(pLab - 1.2)/(g4
411 fElasticXsc = 6 + 20/( (logP-0.182)
412 }
413 }
414 if( proton ) // pn to be np
415 {
416 if( pLab < 0.02 )
417 {
418 fTotalXsc = 4100+30*G4Exp(G4Log(G4Lo
419 fElasticXsc = fTotalXsc;
420 }
421 else if( pLab < 0.8 )
422 {
423 fTotalXsc = 33+30*g4calc->powN(G4Log
424 fElasticXsc = fTotalXsc;
425 }
426 else if( pLab < 1.4 )
427 {
428 fTotalXsc = 33+30*g4calc->powN(G4Log
429 G4double x = G4Log(0.511/pLab);
430 fElasticXsc = 6 + 52/( x*x + 1.6 );
431 }
432 else // 1.4 < pLab < 10. )
433 {
434 fTotalXsc = 33.3 + 20.8*(pLab*pLab -
435 /(std::sqrt(g4calc->powN(pLab,5)) + 0.95
436 fElasticXsc = 6 + 20/( (logP-0.182)
437 }
438 }
439 }
440 }
441 ////// proton //////////////////////////////
442 else if(theParticle == theProton)
443 {
444 if( pLab >= 373.) // pdg due to TOTEM data
445 {
446 fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4
447 + 9.19*G4Exp(-G4Log(sMand)*0.458);
448 }
449 else if( pLab >= 100.)
450 {
451 fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G
452 + 9.19*G4Exp(-G4Log(sMand)*0.458);
453 }
454 else if( pLab >= 10.)
455 {
456 fElasticXsc = 6. + 20./( (logP-0.182)*(l
457 }
458 else
459 {
460 // pp
461 if( proton )
462 {
463 if( pLab < 0.73 )
464 {
465 fTotalXsc = 23 + 50*std::sqrt(g4calc
466 fElasticXsc = fTotalXsc;
467 }
468 else if( pLab < 1.05 )
469 {
470 G4double x = G4Log(pLab/0.73);
471 fTotalXsc = 23 + 40*x*x;
472 fElasticXsc = 23 + 20*x*x;
473 }
474 else // 1.05 - 10 GeV/c
475 {
476 fTotalXsc = 39.0+75*(pLab - 1.2)/(g4
477 fElasticXsc = 6. + 20./( (logP-0.182
478 }
479 }
480 else if( neutron ) // pn to be np
481 {
482 if( pLab < 0.02 )
483 {
484 fTotalXsc = 4100+30*G4Exp(G4Log(G4Lo
485 fElasticXsc = fTotalXsc;
486 }
487 else if( pLab < 0.8 )
488 {
489 fTotalXsc = 33+30*g4calc->powN(G4Log
490 fElasticXsc = fTotalXsc;
491 }
492 else if( pLab < 1.4 )
493 {
494 G4double x1 = G4Log(pLab/0.95);
495 G4double x2 = G4Log(0.511/pLab);
496 fTotalXsc = 33+30*x1*x1;
497 fElasticXsc = 6 + 52/(x2*x2 + 1.6);
498 }
499 else // 1.4 < pLab < 10. )
500 {
501 fTotalXsc = 33.3 + 20.8*(pLab*pLab -
502 /(std::sqrt(g4calc->powN(pLab,5)) + 0.95
503 fElasticXsc = 6. + 20./( (logP-0.182
504 }
505 }
506 }
507 }
508 // pi+ p; pi- n
509 else if((pdg == 211 && proton) || (pdg == -2
510 {
511 if( pLab < 0.28 )
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 {
529 G4double x = G4Log(pLab/0.77);
530 fTotalXsc = 88.*x*x + 14.9;
531 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*(
532 }
533 else if( pLab < 1.4) // ns original
534 {
535 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p
536 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL
537 fTotalXsc = Ex1 + Ex2 + 27.5;
538 fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*(
539 }
540 else if( pLab < 2.0 ) // ns original
541 {
542 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p
543 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL
544 fTotalXsc = Ex1 + Ex2 + 27.5;
545 fElasticXsc = 3.0 + 1.36/( (logP - 0.336
546 }
547 else if( pLab < 3.5 ) // ns original
548 {
549 G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(p
550 G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pL
551 fTotalXsc = Ex1 + Ex2 + 27.5;
552 fElasticXsc = 3.0 + 6.20/( (logP - 0.336
553 }
554 else if( pLab < 10. ) // my
555 {
556 fTotalXsc = 10.6 + 2.*G4Log(pE) + 25*G4E
557 fElasticXsc = 3.0 + 6.20/( (logP - 0.336
558 }
559 else // pLab > 10 // my
560 {
561 fElasticXsc = 3.0 + 6.20/( (logP - 0.336
562 }
563 }
564 // pi+ n; pi- p
565 else if((pdg == 211 && neutron) || (pdg == -
566 {
567 if( pLab < 0.28 )
568 {
569 fTotalXsc = 0.288/((pLab - 0.28)*(pLab -
570 fElasticXsc = 1.8/((logP + 1.273)*(logP
571 }
572 else if( pLab < 0.395676 ) // first peak
573 {
574 fTotalXsc = 0.648/((pLab - 0.28)*(pLab -
575 fElasticXsc = 0.257/((pLab - 0.28)*(pLab
576 }
577 else if( pLab < 0.5 )
578 {
579 G4double y = G4Log(pLab/0.48);
580 fTotalXsc = 26 + 110*y*y;
581 fElasticXsc = 0.37*fTotalXsc;
582 }
583 else if( pLab < 0.65 )
584 {
585 G4double x = G4Log(pLab/0.48);
586 fTotalXsc = 26. + 110.*x*x;
587 fElasticXsc = 0.95/((pLab - 0.72)*(pLab
588 }
589 else if( pLab < 0.72 )
590 {
591 fTotalXsc = 36.1 + 10*G4Exp(-(pLab-0.72)
592 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.
593 fElasticXsc = 0.95/((pLab - 0.72)*(pLab
594 }
595 else if( pLab < 0.88 )
596 {
597 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72
598 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.
599 fElasticXsc = 0.95/((pLab - 0.72)*(pLab
600 }
601 else if( pLab < 1.03 )
602 {
603 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72
604 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.
605 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(
606 }
607 else if( pLab < 1.15 )
608 {
609 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72
610 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.
611 fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(
612 }
613 else if( pLab < 1.3 )
614 {
615 fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72
616 24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.
617 fElasticXsc = 3. + 13./pLab;
618 }
619 else if( pLab < 10.) // < 3.0) // ns origi
620 {
621 fTotalXsc = 36.1 + 0.079-4.313*logP+
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);
624 fElasticXsc = 3. + 13./pLab;
625 }
626 else // mb
627 {
628 fElasticXsc = 3. + 13./pLab;
629 }
630 }
631 else if( (theParticle == theKMinus) && proto
632 {
633 if( pLab < pMin)
634 {
635 G4double psp = pLab*std::sqrt(pLab);
636 fElasticXsc = 5.2/psp;
637 fTotalXsc = 14./psp;
638 }
639 else if( pLab > pMax )
640 {
641 G4double ld = logP - minLogP;
642 G4double ld2 = ld*ld;
643 fElasticXsc = cofLogE*ld2 + 2.23;
644 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
645 }
646 else
647 {
648 G4double ld = logP - minLogP;
649 G4double ld2 = ld*ld;
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
658 G4double lm = pLab - .39;
659 G4double md = lm*lm + .000356;
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
665 // small peaks were added but commented
666 fElasticXsc = 5.2/psp + (1.1*cofLogE*ld
667 + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd;
668
669 fTotalXsc = 14./psp + (1.1*cofLogT*ld2
670 + .006/md + 0.01/hd1+ 0.02/hd2 + .20/
671 }
672 }
673 else if( (theParticle == theKMinus) && neutr
674 {
675 if( pLab > pMax )
676 {
677 G4double ld = logP - minLogP;
678 G4double ld2 = ld*ld;
679 fElasticXsc = cofLogE*ld2 + 2.23;
680 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
681 }
682 else
683 {
684 G4double lh = pLab - 0.98;
685 G4double hd = lh*lh + .021;
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 }
692 }
693 else if( (theParticle == theKPlus) && proton
694 {
695 // VI: modified low-energy part
696 if( pLab < 0.631 )
697 {
698 fElasticXsc = fTotalXsc = 12.03;
699 }
700 else if( pLab > pMax )
701 {
702 G4double ld = logP - minLogP;
703 G4double ld2 = ld*ld;
704 fElasticXsc = cofLogE*ld2 + 2.23;
705 fTotalXsc = cofLogT*ld2 + 19.2;
706 }
707 else
708 {
709 G4double ld = logP - minLogP;
710 G4double ld2 = ld*ld;
711 G4double lr = pLab - .38;
712 G4double LE = .7/(lr*lr + .076);
713 G4double sp = std::sqrt(pLab);
714 G4double p2 = pLab*pLab;
715 G4double p4 = p2*p2;
716 // VI: tuned elastic
717 fElasticXsc = LE + (cofLogE*ld2 + 2.23)
718 + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652);
719 fTotalXsc = LE + (cofLogT*ld2 + 19.5)
720 + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392);
721 }
722 }
723 else if( (theParticle == theKPlus) && neutr
724 {
725 if( pLab < pMin )
726 {
727 G4double lm = pLab - 0.94;
728 G4double md = lm*lm + .392;
729 fElasticXsc = 2./md;
730 fTotalXsc = 4.6/md;
731 }
732 else if( pLab > pMax )
733 {
734 G4double ld = logP - minLogP;
735 G4double ld2 = ld*ld;
736 fElasticXsc = cofLogE*ld2 + 2.23;
737 fTotalXsc = cofLogT*ld2 + 19.2;
738 }
739 else
740 {
741 G4double ld = logP - minLogP;
742 G4double ld2 = ld*ld;
743 G4double sp = std::sqrt(pLab);
744 G4double p2 = pLab*pLab;
745 G4double p4 = p2*p2;
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 }
772
773 //////////////////////////////////////////////
774 //
775 // Returns kaon-nucleon cross-section based on
776 // tuned for the Glauber-Gribov hadron model f
777
778 G4double G4HadronNucleonXsc::KaonNucleonXscGG(
779 const G4ParticleD
780 const G4ParticleDefinition* nucleon
781 {
782 fTotalXsc = fElasticXsc = fInelasticXsc = 0.
783 if(theParticle == theKMinus || theParticle =
784 KaonNucleonXscVG(theParticle, nucleon, eki
785
786 } else if(theParticle == theK0S || thePartic
787 G4double stot = KaonNucleonXscVG(theKMinu
788 G4double sel = fElasticXsc;
789 G4double sinel = fInelasticXsc;
790 stot += KaonNucleonXscVG(theKPlus, nucleo
791 sel += fElasticXsc;
792 sinel += fInelasticXsc;
793 fTotalXsc = stot*0.5;
794 fElasticXsc = sel*0.5;
795 fInelasticXsc = sinel*0.5;
796 }
797 return fTotalXsc;
798 }
799
800 //////////////////////////////////////////////
801 //
802 // Returns kaon-nucleon cross-section using NS
803
804 G4double G4HadronNucleonXsc::KaonNucleonXscNS(
805 const G4ParticleD
806 const G4ParticleDefinition* nucleon
807 {
808 fTotalXsc = fElasticXsc = fInelasticXsc = 0.
809 if(theParticle == theKMinus || theParticle =
810 HadronNucleonXscNS(theParticle, nucleon, e
811
812 } else if(theParticle == theK0S || thePartic
813 G4double fact = 0.5;
814 G4double stot = 0.0;
815 G4double sel = 0.0;
816 G4double sinel= 0.0;
817 if(ekin > ekinmaxQB) {
818 stot = HadronNucleonXscNS(theKMinus, nu
819 sel = fElasticXsc;
820 sinel = fInelasticXsc;
821 stot += HadronNucleonXscNS(theKPlus, nu
822 sel += fElasticXsc;
823 sinel += fInelasticXsc;
824 } else {
825 fact *= std::sqrt(ekinmaxQB/std::max(eki
826 stot = HadronNucleonXscNS(theKMinus, nu
827 sel = fElasticXsc;
828 sinel = fInelasticXsc;
829 stot += HadronNucleonXscNS(theKPlus, nu
830 sel += fElasticXsc;
831 sinel += fInelasticXsc;
832 }
833 fTotalXsc = stot*fact;
834 fElasticXsc = sel*fact;
835 fInelasticXsc = sinel*fact;
836 }
837 return fTotalXsc;
838 }
839
840 //////////////////////////////////////////////
841 //
842 // Returns kaon-nucleon cross-section with smo
843
844 G4double G4HadronNucleonXsc::KaonNucleonXscVG(
845 const G4ParticleDefinition* t
846 const G4ParticleDefinition* nucleon, G4do
847 {
848 G4double pM = theParticle->GetPDGMass();
849 G4double pLab = std::sqrt(ekin*(ekin + 2*pM)
850
851 pLab *= invGeV;
852 G4double logP = G4Log(pLab);
853
854 fTotalXsc = 0.0;
855
856 G4bool proton = (nucleon == theProton);
857 G4bool neutron = (nucleon == theNeutron);
858
859 if( (theParticle == theKMinus) && proton )
860 {
861 if( pLab < pMin)
862 {
863 G4double psp = pLab*std::sqrt(pLab);
864 fElasticXsc = 5.2/psp;
865 fTotalXsc = 14./psp;
866 }
867 else if( pLab > pMax )
868 {
869 G4double ld = logP - minLogP;
870 G4double ld2 = ld*ld;
871 fElasticXsc = cofLogE*ld2 + 2.23;
872 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
873 }
874 else
875 {
876 G4double ld = logP - minLogP;
877 G4double ld2 = ld*ld;
878 G4double sp = std::sqrt(pLab);
879 G4double psp = pLab*sp;
880 G4double p2 = pLab*pLab;
881 G4double p4 = p2*p2;
882
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 }
888 }
889 else if( (theParticle == theKMinus) && neutr
890 {
891 if( pLab > pMax )
892 {
893 G4double ld = logP - minLogP;
894 G4double ld2 = ld*ld;
895 fElasticXsc = cofLogE*ld2 + 2.23;
896 fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
897 }
898 else
899 {
900 G4double lh = pLab - 0.98;
901 G4double hd = lh*lh + .045; // vg ve
902 G4double sqrLogPlab = logP*logP;
903
904 fElasticXsc = 5.0 + 8.1*G4Exp(-logP*1.8
905 - 1.3*logP + .15/hd;
906 fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2.
907 }
908 }
909 else if( (theParticle == theKPlus) && proton
910 {
911 // VI: modified low-energy part
912 if( pLab < 0.631 )
913 {
914 fElasticXsc = fTotalXsc = 12.03;
915 }
916 else if( pLab > pMax )
917 {
918 G4double ld = logP - minLogP;
919 G4double ld2 = ld*ld;
920 fElasticXsc = cofLogE*ld2 + 2.23;
921 fTotalXsc = cofLogT*ld2 + 19.2;
922 }
923 else
924 {
925 G4double ld = logP - minLogP;
926 G4double ld2 = ld*ld;
927 G4double lr = pLab - .38;
928 G4double LE = .7/(lr*lr + .076);
929 G4double sp = std::sqrt(pLab);
930 G4double p2 = pLab*pLab;
931 G4double p4 = p2*p2;
932 // VI: tuned elastic
933 fElasticXsc = LE + (cofLogE*ld2 + 2.23)
934 + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652);
935 fTotalXsc = LE + (cofLogT*ld2 + 19.5)
936 + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392);
937 }
938 }
939 else if( (theParticle == theKPlus) && neutr
940 {
941 if( pLab < pMin )
942 {
943 G4double lm = pLab - 0.94;
944 G4double md = lm*lm + .392;
945 fElasticXsc = 2./md;
946 fTotalXsc = 4.6/md;
947 }
948 else if( pLab > pMax )
949 {
950 G4double ld = logP - minLogP;
951 G4double ld2 = ld*ld;
952 fElasticXsc = cofLogE*ld2 + 2.23;
953 fTotalXsc = cofLogT*ld2 + 19.2;
954 }
955 else
956 {
957 G4double ld = logP - minLogP;
958 G4double ld2 = ld*ld;
959 G4double sp = std::sqrt(pLab);
960 G4double p2 = pLab*pLab;
961 G4double p4 = p2*p2;
962 G4double lm = pLab - 0.94;
963 G4double md = lm*lm + .392;
964 fElasticXsc = (cofLogE*ld2 + 2.23)/(1.
965 fTotalXsc = (cofLogT*ld2 + 19.5)/(1.
966 }
967 }
968
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 {
1101 coeff = 0.76;
1102 }
1103 // eta'
1104 else if( pdg == 331 )
1105 {
1106 coeff = 0.88;
1107 }
1108 fTotalXsc = coeff*HadronNucleonXscNS(thePiP
1109 fElasticXsc *= coeff;
1110 fInelasticXsc *= coeff;
1111 return fTotalXsc;
1112 }
1113 /////////////////////////////////////////////
1114 //
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 {
1137 fTotalXsc = 48.0 + 0.522*sqrLogPlab -
1138 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla
1139 }
1140 if(neutron)
1141 {
1142 fTotalXsc = 47.3 + 0.513*sqrLogPlab
1143 fElasticXsc = 11.9 + 26.9*G4Exp(-logPla
1144 }
1145 }
1146 else if( PDGcode == 211) //------Projecti
1147 {
1148 if(proton)
1149 {
1150 fTotalXsc = 16.4 + 19.3 *G4Exp(-logPla
1151 fElasticXsc = 0.0 + 11.4*G4Exp(-logPla
1152 }
1153 if(neutron)
1154 {
1155 fTotalXsc = 33.0 + 14.0 *G4Exp(-logP
1156 fElasticXsc = 1.76 + 11.2*G4Exp(-logPla
1157 }
1158 }
1159 else if( PDGcode == -211) //------Projecti
1160 {
1161 if(proton)
1162 {
1163 fTotalXsc = 33.0 + 14.0 *G4Exp(-logPl
1164 fElasticXsc = 1.76 + 11.2*G4Exp(-logPla
1165 }
1166 if(neutron)
1167 {
1168 fTotalXsc = 16.4 + 19.3 *G4Exp(-logPl
1169 fElasticXsc = 0.0 + 11.4*G4Exp(-logPla
1170 }
1171 }
1172 else if( PDGcode == 111) //------Projecti
1173 {
1174 if(proton)
1175 {
1176 fTotalXsc = (16.4 + 19.3 *G4Exp(-logPla
1177 33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.
1178
1179 fElasticXsc = (0.0 + 11.4*G4Exp(-logPla
1180 1.76 + 11.2*G4Exp(-logPlab*0.64) + 0
1181
1182 }
1183 if(neutron)
1184 {
1185 fTotalXsc = (33.0 + 14.0 *G4Exp(-logPla
1186 16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.
1187 fElasticXsc = (1.76 + 11.2*G4Exp(-logPl
1188 0.0 + 11.4*G4Exp(-logPlab*0.40) + 0
1189 }
1190 }
1191 else if( PDGcode == 321 ) //------Projec
1192 {
1193 if(proton)
1194 {
1195 fTotalXsc = 18.1 + 0.26 *sqrLogPlab
1196 fElasticXsc = 5.0 + 8.1*G4Exp(-logPla
1197 }
1198 if(neutron)
1199 {
1200 fTotalXsc = 18.7 + 0.21 *sqrLogPlab
1201 fElasticXsc = 7.3 + 0.29 *sqrLogPlab
1202 }
1203 }
1204 else if( PDGcode ==-321 ) //------Projecti
1205 {
1206 if(proton)
1207 {
1208 fTotalXsc = 32.1 + 0.66*sqrLogPlab -
1209 fElasticXsc = 7.3 + 0.29*sqrLogPlab -
1210 }
1211 if(neutron)
1212 {
1213 fTotalXsc = 25.2 + 0.38*sqrLogPlab -
1214 fElasticXsc = 5.0 + 8.1*G4Exp(-logPla
1215 }
1216 }
1217 else if( PDGcode == 311 ) //------Projecti
1218 {
1219 if(proton)
1220 {
1221 fTotalXsc = ( 18.1 + 0.26 *sqrLogPlab
1222 32.1 + 0.66 *sqrLogPlab
1223 fElasticXsc = ( 5.0 + 8.1*G4Exp(-logP
1224 7.3 + 0.29 *sqrLogPl
1225 }
1226 if(neutron)
1227 {
1228 fTotalXsc = ( 18.7 + 0.21 *sqrLogPlab
1229 25.2 + 0.38 *sqrLogPlab
1230 fElasticXsc = ( 7.3 + 0.29 *sqrLogPlab
1231 5.0 + 8.1*G4Exp(-logPl
1232 }
1233 }
1234 else //------Projectile is undefined, Nucl
1235 {
1236 if(proton)
1237 {
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 }
1247
1248 fTotalXsc *= CLHEP::millibarn;
1249 fElasticXsc *= CLHEP::millibarn;
1250 fElasticXsc = std::min(fElasticXsc, fTotalX
1251 fInelasticXsc = fTotalXsc - fElasticXsc;
1252
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
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;
1323 }
1324
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
1335 if ( theParticle == theProton ) pR = 0.
1336 else if( theParticle == thePiPlus ) pR = 0.
1337 else if( theParticle == theKPlus ) pR = 0.
1338
1339 G4double pZ = theParticle->GetPDGCharge();
1340 G4double tZ = nucleon->GetPDGCharge();
1341
1342 G4double pM = theParticle->GetPDGMass();
1343 G4double tM = nucleon->GetPDGMass();
1344
1345 G4double pElab = ekin + pM;
1346
1347 G4double totEcm = std::sqrt(pM*pM + tM*tM
1348
1349 G4double totTcm = totEcm - pM -tM;
1350
1351 G4double bC = fine_structure_const*hbarc*pZ
1352
1353 //G4cout<<"##CB: ekin = "<<ekin<<"; pElab=
1354 // <<"; bC = "<<bC<<"; Tcm = "<< totTcm <<
1355
1356 G4double ratio = (totTcm > bC) ? 1. - bC/to
1357
1358 // G4cout <<"ratio = "<<ratio<<G4endl;
1359 return ratio;
1360 }
1361
1362 /////////////////////////////////////////////
1363