ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HadronNucleonXsc.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4HadronNucleonXsc.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // 14.03.07 V. Grichine - first implementation
27 //
28 // 04.09.18 V. Ivantchenko Major revision of interfaces and implementation
29 // 30.09.18 V. Grichine hyperon-nucleon xsc first implementation
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 "G4LambdacPlus.hh"
46 #include "G4AntiLambdacPlus.hh"
47 #include "G4AntiXibZero.hh"
48 #include "G4OmegacZero.hh"
49 #include "G4SigmacZero.hh"
50 #include "G4AntiLambdab.hh"
51 #include "G4AntiSigmabMinus.hh"
52 #include "G4AntiXicPlus.hh"
53 #include "G4AntiLambdacPlus.hh"
54 #include "G4AntiSigmabPlus.hh"
55 #include "G4AntiXicZero.hh"
56 #include "G4AntiSigmabZero.hh"
57 #include "G4XibMinus.hh"
58 #include "G4AntiSigmacPlus.hh"
59 #include "G4XibZero.hh"
60 #include "G4AntiOmegabMinus.hh"
61 #include "G4AntiSigmacPlusPlus.hh"
62 
63 #include "G4Lambdab.hh"
64 #include "G4SigmabMinus.hh"
65 #include "G4XicPlus.hh"
66 #include "G4AntiOmegacZero.hh"
67 #include "G4AntiSigmacZero.hh"
68 #include "G4LambdacPlus.hh"
69 #include "G4SigmabPlus.hh"
70 #include "G4XicZero.hh"
71 #include "G4SigmabZero.hh"
72 #include "G4SigmacPlus.hh"
73 #include "G4AntiXibMinus.hh"
74 #include "G4OmegabMinus.hh"
75 #include "G4SigmacPlusPlus.hh"
76 
77 #include "G4BMesonZero.hh"
78 #include "G4AntiBMesonZero.hh"
79 #include "G4DMesonZero.hh"
80 #include "G4AntiDMesonZero.hh"
81 #include "G4BsMesonZero.hh"
82 #include "G4AntiBsMesonZero.hh"
83 #include "G4BcMesonPlus.hh"
84 #include "G4BcMesonMinus.hh"
85 #include "G4DsMesonPlus.hh"
86 #include "G4DsMesonMinus.hh"
87 #include "G4Eta.hh"
88 #include "G4EtaPrime.hh"
89 #include "G4Etac.hh"
90 
91 #include "G4BMesonPlus.hh"
92 #include "G4BMesonMinus.hh"
93 
94 #include "G4DMesonPlus.hh"
95 #include "G4DMesonMinus.hh"
96 
97 #include "G4JPsi.hh"
98 #include "G4Upsilon.hh"
99 
100 #include "G4Lambda.hh"
101 #include "G4AntiLambda.hh"
102 #include "G4SigmaPlus.hh"
103 #include "G4AntiSigmaPlus.hh"
104 #include "G4SigmaMinus.hh"
105 #include "G4AntiSigmaMinus.hh"
106 #include "G4SigmaZero.hh"
107 #include "G4AntiSigmaZero.hh"
108 #include "G4XiMinus.hh"
109 #include "G4XiZero.hh"
110 #include "G4AntiXiMinus.hh"
111 #include "G4AntiXiZero.hh"
112 #include "G4OmegaMinus.hh"
113 #include "G4AntiOmegaMinus.hh"
114 
115 static const G4double invGeV = 1.0/CLHEP::GeV;
116 static const G4double invGeV2 = 1.0/(CLHEP::GeV*CLHEP::GeV);
117 // PDG fit constants
118 static const G4double minLogP = 3.5; // min of (lnP-minLogP)^2
119 static const G4double cofLogE = .0557; // elastic (lnP-minLogP)^2
120 static const G4double cofLogT = .3; // total (lnP-minLogP)^2
121 static const G4double pMin = .1; // fast LE calculation
122 static const G4double pMax = 1000.; // fast HE calculation
123 static const G4double ekinmin = 0.1*CLHEP::MeV; // protection against zero ekin
124 
125 G4HadronNucleonXsc::G4HadronNucleonXsc()
126  : fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0)
127 {
128  fHypTotXscCof = 0.88; // for transformation pp(pn) to hyperon-nucleon
129 
130  theGamma = G4Gamma::Gamma();
131  theProton = G4Proton::Proton();
132  theNeutron = G4Neutron::Neutron();
133  theAProton = G4AntiProton::AntiProton();
134  theANeutron = G4AntiNeutron::AntiNeutron();
135  thePiPlus = G4PionPlus::PionPlus();
136  thePiMinus = G4PionMinus::PionMinus();
137  thePiZero = G4PionZero::PionZero();
138  theD = G4Deuteron::Deuteron();
139  theT = G4Triton::Triton();
140  theA = G4Alpha::Alpha();
141  theHe3 = G4He3::He3();
142  // strange
143  theKPlus = G4KaonPlus::KaonPlus();
144  theKMinus = G4KaonMinus::KaonMinus();
145  theK0S = G4KaonZeroShort::KaonZeroShort();
146  theK0L = G4KaonZeroLong::KaonZeroLong();
147  theL = G4Lambda::Lambda();
148  theAntiL = G4AntiLambda::AntiLambda();
149  theSPlus = G4SigmaPlus::SigmaPlus();
150  theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
151  theSMinus = G4SigmaMinus::SigmaMinus();
152  theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
153  theS0 = G4SigmaZero::SigmaZero();
154  theAS0 = G4AntiSigmaZero::AntiSigmaZero();
155  theXiMinus = G4XiMinus::XiMinus();
156  theXi0 = G4XiZero::XiZero();
157  theAXiMinus = G4AntiXiMinus::AntiXiMinus();
158  theAXi0 = G4AntiXiZero::AntiXiZero();
159  theOmega = G4OmegaMinus::OmegaMinus();
160  theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
161  // c- and b- hyperons
162  theLambdaCPlus = G4LambdacPlus::LambdacPlus();
163  theALambdaCPlus = G4AntiLambdacPlus::AntiLambdacPlus();
164  theOmegaC0 = G4OmegacZero::OmegacZero();
165  theAOmegaC0 = G4AntiOmegacZero::AntiOmegacZero();
166  theSigmaCPlus = G4SigmacPlus::SigmacPlus();
167  theASigmaCPlus = G4AntiSigmacPlus::AntiSigmacPlus();
168  theSigmacPP = G4SigmacPlusPlus::SigmacPlusPlus();
169  theASigmacPP = G4AntiSigmacPlusPlus::AntiSigmacPlusPlus();
170  theSigmaC0 = G4SigmacZero::SigmacZero();
171  theASigmaC0 = G4AntiSigmacZero::AntiSigmacZero();
172  theXiCPlus = G4XicPlus::XicPlus();
173  theAXiCPlus = G4AntiXicPlus::AntiXicPlus();
174  theXiC0 = G4XicZero::XicZero();
175  theAXiC0 = G4AntiXicZero::AntiXicZero();
176  theLambdaB = G4Lambdab::Lambdab();
177  theALambdaB = G4AntiLambdab::AntiLambdab();
178  theOmegaBMinus = G4OmegabMinus::OmegabMinus();
179  theAOmegaBMinus = G4AntiOmegabMinus::AntiOmegabMinus();
180  theSigmaBMinus = G4SigmabMinus::SigmabMinus();
181  theASigmaBMinus = G4AntiSigmabMinus::AntiSigmabMinus();
182  theSigmaBPlus = G4SigmabPlus::SigmabPlus();
183  theASigmaBPlus = G4AntiSigmabPlus::AntiSigmabPlus();
184  theSigmaB0 = G4SigmabZero::SigmabZero();
185  theASigmaB0 = G4AntiSigmabZero::AntiSigmabZero();
186  theXiBMinus = G4XibMinus::XibMinus();
187  theAXiBMinus = G4AntiXibMinus::AntiXibMinus();
188  theXiB0 = G4XibZero::XibZero();
189  theAXiB0 = G4AntiXibZero::AntiXibZero();
190  //(s-) c- and b-mesons
191  theBMeson0 = G4BMesonZero::BMesonZero();
192  theABMeson0 = G4AntiBMesonZero::AntiBMesonZero();
193  theDMeson0 = G4DMesonZero::DMesonZero();
194  theADMeson0 = G4AntiDMesonZero::AntiDMesonZero();
195  theBsMeson0 = G4BsMesonZero::BsMesonZero();
196  theABsMeson0 = G4AntiBsMesonZero::AntiBsMesonZero();
197  theBcMesonPlus = G4BcMesonPlus::BcMesonPlus();
198  theBcMesonMinus = G4BcMesonMinus::BcMesonMinus();
199  theDsMesonPlus = G4DsMesonPlus::DsMesonPlus();
200  theDsMesonMinus = G4DsMesonMinus::DsMesonMinus();
201  theDMesonPlus = G4DMesonPlus::DMesonPlus();
202  theDMesonMinus = G4DMesonMinus::DMesonMinus();
203  theBMesonPlus = G4BMesonPlus::BMesonPlus();
204  theBMesonMinus = G4BMesonMinus::BMesonMinus();
205  theEta = G4Eta::Eta();
206  theEtaPrime = G4EtaPrime::EtaPrime();
207  theEtaC = G4Etac::Etac();
208  theJPsi = G4JPsi::JPsi();
209  theUpsilon = G4Upsilon::Upsilon();
210 
211  g4calc = G4Pow::GetInstance();
212 }
213 
214 G4HadronNucleonXsc::~G4HadronNucleonXsc()
215 {}
216 
217 void G4HadronNucleonXsc::CrossSectionDescription(std::ostream& outFile) const
218 {
219  outFile << "G4HadronNucleonXsc calculates the total, inelastic and elastic\n"
220  << "hadron-nucleon cross sections using several different\n"
221  << "parameterizations within the Glauber-Gribov approach. It is\n"
222  << "valid for all incident gammas and long-lived hadrons at\n"
223  << "energies above 30 keV. This is a cross section component which\n"
224  << "is to be used to build a cross section data set.\n";
225 }
226 
227 
228 
229 G4double G4HadronNucleonXsc::HadronNucleonXsc( const G4ParticleDefinition* theParticle,
230  const G4ParticleDefinition* nucleon, G4double ekin)
231 {
232  G4double xsc(0.);
233  G4int pdg = std::abs( theParticle->GetPDGEncoding() );
234 
235  if ( pdg == 2212 || pdg == 2112 || pdg == 211 ) // p, n, pi+-
236  {
237  xsc = HadronNucleonXscNS( theParticle, nucleon, ekin);
238  }
239  else if ( pdg == 321 || pdg == 310 || pdg == 130 ) // K+-, K0, Ks
240  {
241  xsc = KaonNucleonXscNS( theParticle, nucleon, ekin);
242  }
243  else if ( pdg == 3122 || pdg == 3222 || pdg == 3112 || pdg == 3212 || pdg == 3322 || pdg == 3312 || pdg == 3324 ||
244 
245  pdg == 4122 || pdg == 4332 || pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 4112 || pdg == 4232 || pdg == 4132 ||
246 
247  pdg == 5122 || pdg == 5332 || pdg == 5122 || pdg == 5112 || pdg == 5222 || pdg == 5212 || pdg == 5132 || pdg == 5232
248 
249  ) // heavy s-,c-,b-hyperons
250  {
251  xsc = HyperonNucleonXscNS( theParticle, nucleon, ekin);
252  }
253  else if ( pdg == 511 || pdg == 421 || pdg == 531 || pdg == 541 || pdg == 431 || pdg == 411 || pdg == 521 ||
254 
255  pdg == 221 || pdg == 331 || pdg == 441 || pdg == 443 || pdg == 543
256 
257  ) // s-,c-,b-mesons
258  {
259  xsc = SCBMesonNucleonXscNS( theParticle, nucleon, ekin);
260  }
261  else
262  {
263  xsc = HadronNucleonXscNS( theParticle, nucleon, ekin);
264  }
265  return xsc;
266 }
267 
268 
270 //
271 // Returns hadron-nucleon Xsc according to PDG parametrisation (2017):
272 // http://pdg.lbl.gov/2017/reviews/hadronicrpp.pdf
273 
274 G4double G4HadronNucleonXsc::HadronNucleonXscPDG(
275  const G4ParticleDefinition* theParticle,
276  const G4ParticleDefinition* nucleon, G4double ekin)
277 {
278  static const G4double M = 2.1206; // in Gev
279  static const G4double eta1 = 0.4473;
280  static const G4double eta2 = 0.5486;
281  static const G4double H = 0.272;
282 
283  G4double mass1 = theParticle->GetPDGMass();
284  if(theParticle == theGamma) { mass1 = 770.; }
285  G4double mass2 = nucleon->GetPDGMass();
286 
287  G4double sMand = CalcMandelstamS(ekin, mass1, mass2)*invGeV2;
288  G4double x = (mass1 + mass2)*invGeV + M;
289  G4double blog = G4Log(sMand/(x*x));
290 
291  G4double P(0.0), R1(0.0), R2(0.0), del(1.0);
292 
293  G4bool proton = (nucleon == theProton);
294  G4bool neutron = (nucleon == theNeutron);
295 
296  if(theParticle == theNeutron)
297  {
298  if ( proton )
299  {
300  P = 34.71;
301  R1 = 12.52;
302  R2 = -6.66;
303  }
304  else
305  {
306  P = 34.41;
307  R1 = 13.07;
308  R2 = -7.394;
309  }
310  }
311  else if(theParticle == theProton)
312  {
313  if ( neutron )
314  {
315  P = 34.71;
316  R1 = 12.52;
317  R2 = -6.66;
318  }
319  else
320  {
321  P = 34.41;
322  R1 = 13.07;
323  R2 = -7.394;
324  }
325  }
326  else if(theParticle == theAProton)
327  {
328  if ( neutron )
329  {
330  P = 34.71;
331  R1 = 12.52;
332  R2 = 6.66;
333  }
334  else
335  {
336  P = 34.41;
337  R1 = 13.07;
338  R2 = 7.394;
339  }
340  }
341  else if(theParticle == theANeutron)
342  {
343  if ( proton )
344  {
345  P = 34.71;
346  R1 = 12.52;
347  R2 = 6.66;
348  }
349  else
350  {
351  P = 34.41;
352  R1 = 13.07;
353  R2 = 7.394;
354  }
355  }
356  else if(theParticle == thePiPlus)
357  {
358  P = 18.75;
359  R1 = 9.56;
360  R2 = -1.767;
361  }
362  else if(theParticle == thePiMinus)
363  {
364  P = 18.75;
365  R1 = 9.56;
366  R2 = 1.767;
367  }
368  else if(theParticle == theKPlus)
369  {
370  if ( proton )
371  {
372  P = 16.36;
373  R1 = 4.29;
374  R2 = -3.408;
375  }
376  else
377  {
378  P = 16.31;
379  R1 = 3.7;
380  R2 = -1.826;
381  }
382  }
383  else if(theParticle == theKMinus)
384  {
385  if ( proton )
386  {
387  P = 16.36;
388  R1 = 4.29;
389  R2 = 3.408;
390  }
391  else
392  {
393  P = 16.31;
394  R1 = 3.7;
395  R2 = 1.826;
396  }
397  }
398  else if(theParticle == theK0S || theParticle == theK0L)
399  {
400  P = 16.36;
401  R1 = 2.5;
402  R2 = 0.;
403  }
404  else if(theParticle == theSMinus)
405  {
406  P = 34.7;
407  R1 = -46.;
408  R2 = 48.;
409  }
410  else if(theParticle == theGamma) // modify later on
411  {
412  del= 0.003063;
413  P = 34.71*del;
414  R1 = (neutron) ? 0.0231 : 0.0139;
415  R2 = 0.;
416 
417  }
418  else // as proton ???
419  {
420  if ( neutron )
421  {
422  P = 34.71;
423  R1 = 12.52;
424  R2 = -6.66;
425  }
426  else
427  {
428  P = 34.41;
429  R1 = 13.07;
430  R2 = -7.394;
431  }
432  }
433  fTotalXsc = CLHEP::millibarn*
434  (del*(H*blog*blog + P) + R1*G4Exp(-eta1*blog) + R2*G4Exp(-eta2*blog));
435  fInelasticXsc = 0.75*fTotalXsc;
436  fElasticXsc = fTotalXsc - fInelasticXsc;
437 
438  if( proton && theParticle->GetPDGCharge() > 0. && ekin < 100*MeV)
439  {
440  G4double cB = CoulombBarrier(theParticle, nucleon, ekin);
441  fTotalXsc *= cB;
442  fElasticXsc *= cB;
443  fInelasticXsc *= cB;
444  }
445  /*
446  G4cout << "## HadronNucleonXscPDG: ekin(MeV)= " << ekin
447  << " tot= " << fTotalXsc/CLHEP::millibarn
448  << " inel= " << fInelasticXsc/CLHEP::millibarn
449  << " el= " << fElasticXsc/CLHEP::millibarn
450  << G4endl;
451  */
452  return fTotalXsc;
453 }
454 
456 //
457 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
458 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
459 
460 G4double G4HadronNucleonXsc::HadronNucleonXscNS(
461  const G4ParticleDefinition* theParticle,
462  const G4ParticleDefinition* nucleon, G4double ekin0)
463 {
464  const G4double ekin = std::max(ekin0, ekinmin);
465  /*
466  G4cout<< "HadronNucleonXscNS: Ekin(GeV)= " << ekin/GeV << " "
467  << theParticle->GetParticleName() << " + "
468  << nucleon->GetParticleName() << G4endl;
469  */
470  if(theParticle == theAProton || theParticle == theANeutron) {
471  return HadronNucleonXscPDG(theParticle, nucleon, ekin);
472  }
473 
474  G4double pM = theParticle->GetPDGMass();
475  G4double tM = nucleon->GetPDGMass();
476  G4double pE = ekin + pM;
477  G4double pLab = std::sqrt(ekin*(ekin + 2*pM));
478 
479  G4double sMand = CalcMandelstamS(ekin, pM, tM)*invGeV2;
480 
481  pLab *= invGeV;
482  pE *= invGeV;
483 
484  if(pLab >= 10.) {
485  fTotalXsc = HadronNucleonXscPDG(theParticle, nucleon, ekin)/CLHEP::millibarn;
486  } else { fTotalXsc = 0.0; }
487  fElasticXsc = 0.0;
488  //G4cout << "Stot(mb)= " << fTotalXsc << " pLab= " << pLab
489  // << " Smand= " << sMand <<G4endl;
490  G4double logP = G4Log(pLab);
491 
492  G4bool proton = (nucleon == theProton);
493  G4bool neutron = (nucleon == theNeutron);
494 
495  if( theParticle == theNeutron)
496  {
497  if( pLab >= 373.)
498  {
499  fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4Log(sMand/400.)*1.65))
500  + 9.19*G4Exp(-G4Log(sMand)*0.458);
501  }
502  else if( pLab >= 100.)
503  {
504  fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G4Log(sMand/28.9))*1.1)
505  + 9.19*G4Exp(-G4Log(sMand)*0.458);
506  }
507  else if( pLab >= 10.)
508  {
509  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
510  }
511  else // pLab < 10 GeV/c
512  {
513  if( neutron ) // nn to be pp
514  {
515  G4double x = G4Log(pLab/0.73);
516  if( pLab < 0.4 )
517  {
518  fTotalXsc = 23 + 50*std::sqrt(g4calc->powN(-x, 7));
519  fElasticXsc = fTotalXsc;
520  }
521  else if( pLab < 0.73 )
522  {
523  fTotalXsc = 23 + 50*std::sqrt(g4calc->powN(-x, 7));
524  fElasticXsc = fTotalXsc;
525  }
526  else if( pLab < 1.05 )
527  {
528  fTotalXsc = 23 + 40*x*x;
529  fElasticXsc = 23 + 20*x*x;
530  }
531  else // 1.05 - 10 GeV/c
532  {
533  fTotalXsc = 39.0+75*(pLab - 1.2)/(g4calc->powN(pLab,3) + 0.15);
534  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
535  }
536  }
537  if( proton ) // pn to be np
538  {
539  if( pLab < 0.02 )
540  {
541  fTotalXsc = 4100+30*G4Exp(G4Log(G4Log(1.3/pLab))*3.6);
542  fElasticXsc = fTotalXsc;
543  }
544  else if( pLab < 0.8 )
545  {
546  fTotalXsc = 33+30*g4calc->powN(G4Log(pLab/1.3),4);
547  fElasticXsc = fTotalXsc;
548  }
549  else if( pLab < 1.4 )
550  {
551  fTotalXsc = 33+30*g4calc->powN(G4Log(pLab/0.95),2);
552  G4double x = G4Log(0.511/pLab);
553  fElasticXsc = 6 + 52/( x*x + 1.6 );
554  }
555  else // 1.4 < pLab < 10. )
556  {
557  fTotalXsc = 33.3 + 20.8*(pLab*pLab - 1.35)
558  /(std::sqrt(g4calc->powN(pLab,5)) + 0.95);
559  fElasticXsc = 6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
560  }
561  }
562  }
563  }
565  else if(theParticle == theProton)
566  {
567  if( pLab >= 373.) // pdg due to TOTEM data
568  {
569  fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4Log(sMand/400.))*1.65)
570  + 9.19*G4Exp(-G4Log(sMand)*0.458);
571  }
572  else if( pLab >= 100.)
573  {
574  fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G4Log(sMand/28.9))*1.1)
575  + 9.19*G4Exp(-G4Log(sMand)*0.458);
576  }
577  else if( pLab >= 10.)
578  {
579  fElasticXsc = 6. + 20./( (logP-0.182)*(logP-0.182) + 1.0 );
580  }
581  else
582  {
583  // pp
584  if( proton )
585  {
586  if( pLab < 0.73 )
587  {
588  fTotalXsc = 23 + 50*std::sqrt(g4calc->powN(G4Log(0.73/pLab),7));
589  fElasticXsc = fTotalXsc;
590  }
591  else if( pLab < 1.05 )
592  {
593  G4double x = G4Log(pLab/0.73);
594  fTotalXsc = 23 + 40*x*x;
595  fElasticXsc = 23 + 20*x*x;
596  }
597  else // 1.05 - 10 GeV/c
598  {
599  fTotalXsc = 39.0+75*(pLab - 1.2)/(g4calc->powN(pLab,3) + 0.15);
600  fElasticXsc = 6. + 20./( (logP-0.182)*(logP-0.182) + 1.0 );
601  }
602  }
603  else if( neutron ) // pn to be np
604  {
605  if( pLab < 0.02 )
606  {
607  fTotalXsc = 4100+30*G4Exp(G4Log(G4Log(1.3/pLab))*3.6);
608  fElasticXsc = fTotalXsc;
609  }
610  else if( pLab < 0.8 )
611  {
612  fTotalXsc = 33+30*g4calc->powN(G4Log(pLab/1.3),4);
613  fElasticXsc = fTotalXsc;
614  }
615  else if( pLab < 1.4 )
616  {
617  G4double x1 = G4Log(pLab/0.95);
618  G4double x2 = G4Log(0.511/pLab);
619  fTotalXsc = 33+30*x1*x1;
620  fElasticXsc = 6 + 52/(x2*x2 + 1.6);
621  }
622  else // 1.4 < pLab < 10. )
623  {
624  fTotalXsc = 33.3 + 20.8*(pLab*pLab - 1.35)
625  /(std::sqrt(g4calc->powN(pLab,5)) + 0.95);
626  fElasticXsc = 6. + 20./( (logP-0.182)*(logP-0.182) + 1.0 );
627  }
628  }
629  }
630  }
631  // pi+ p; pi- n
632  else if((theParticle == thePiPlus && proton) ||
633  (theParticle == thePiMinus && neutron))
634  {
635  if( pLab < 0.28 )
636  {
637  fTotalXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
638  fElasticXsc = fTotalXsc;
639  }
640  else if( pLab < 0.68 )
641  {
642  fTotalXsc = 14./((logP + 1.273)*(logP + 1.273) + 0.07);
643  fElasticXsc = fTotalXsc;
644  }
645  else if( pLab < 0.85 )
646  {
647  G4double x = G4Log(pLab/0.77);
648  fTotalXsc = 88.*x*x + 14.9;
649  fElasticXsc = fTotalXsc*G4Exp(-3.*(pLab - 0.68));
650  }
651  else if( pLab < 1.15 )
652  {
653  G4double x = G4Log(pLab/0.77);
654  fTotalXsc = 88.*x*x + 14.9;
655  fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
656  }
657  else if( pLab < 1.4) // ns original
658  {
659  G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
660  G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
661  fTotalXsc = Ex1 + Ex2 + 27.5;
662  fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
663  }
664  else if( pLab < 2.0 ) // ns original
665  {
666  G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
667  G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
668  fTotalXsc = Ex1 + Ex2 + 27.5;
669  fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);
670  }
671  else if( pLab < 3.5 ) // ns original
672  {
673  G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
674  G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
675  fTotalXsc = Ex1 + Ex2 + 27.5;
676  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
677  }
678  else if( pLab < 10. ) // my
679  {
680  fTotalXsc = 10.6 + 2.*G4Log(pE) + 25*G4Exp(-G4Log(pE)*0.43 );
681  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
682  }
683  else // pLab > 10 // my
684  {
685  fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);
686  }
687  }
688  // pi+ n; pi- p
689  else if((theParticle == thePiPlus && neutron) ||
690  (theParticle == thePiMinus && proton))
691  {
692  if( pLab < 0.28 )
693  {
694  fTotalXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
695  fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
696  }
697  else if( pLab < 0.395676 ) // first peak
698  {
699  fTotalXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
700  fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
701  }
702  else if( pLab < 0.5 )
703  {
704  G4double y = G4Log(pLab/0.48);
705  fTotalXsc = 26 + 110*y*y;
706  fElasticXsc = 0.37*fTotalXsc;
707  }
708  else if( pLab < 0.65 )
709  {
710  G4double x = G4Log(pLab/0.48);
711  fTotalXsc = 26. + 110.*x*x;
712  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
713  }
714  else if( pLab < 0.72 )
715  {
716  fTotalXsc = 36.1 + 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
717  24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
718  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
719  }
720  else if( pLab < 0.88 )
721  {
722  fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
723  24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
724  fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
725  }
726  else if( pLab < 1.03 )
727  {
728  fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
729  24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
730  fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
731  }
732  else if( pLab < 1.15 )
733  {
734  fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
735  24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
736  fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
737  }
738  else if( pLab < 1.3 )
739  {
740  fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
741  24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
742  fElasticXsc = 3. + 13./pLab;
743  }
744  else if( pLab < 10.) // < 3.0) // ns original
745  {
746  fTotalXsc = 36.1 + 0.079-4.313*logP+
747  3*G4Exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
748  1.5*G4Exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
749  fElasticXsc = 3. + 13./pLab;
750  }
751  else // mb
752  {
753  fElasticXsc = 3. + 13./pLab;
754  }
755  }
756  else if( (theParticle == theKMinus) && proton ) // K-p
757  {
758  if( pLab < pMin)
759  {
760  G4double psp = pLab*std::sqrt(pLab);
761  fElasticXsc = 5.2/psp;
762  fTotalXsc = 14./psp;
763  }
764  else if( pLab > pMax )
765  {
766  G4double ld = logP - minLogP;
767  G4double ld2 = ld*ld;
768  fElasticXsc = cofLogE*ld2 + 2.23;
769  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
770  }
771  else
772  {
773  G4double ld = logP - minLogP;
774  G4double ld2 = ld*ld;
775  G4double sp = std::sqrt(pLab);
776  G4double psp = pLab*sp;
777  G4double p2 = pLab*pLab;
778  G4double p4 = p2*p2;
779 
780  G4double lh = pLab - 1.01;
781  G4double hd = lh*lh + .011;
782 
783  G4double lm = pLab - .39;
784  G4double md = lm*lm + .000356;
785  G4double lh1 = pLab - 0.78;
786  G4double hd1 = lh1*lh1 + .00166;
787  G4double lh2 = pLab - 1.63;
788  G4double hd2 = lh2*lh2 + .007;
789 
790  // small peaks were added but commented out now
791  fElasticXsc = 5.2/psp + (1.1*cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4)
792  + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd;
793 
794  fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4)
795  + .006/md + 0.01/hd1+ 0.02/hd2 + .20/hd ;
796  }
797  }
798  else if( (theParticle == theKMinus) && neutron ) // K-n
799  {
800  if( pLab > pMax )
801  {
802  G4double ld = logP - minLogP;
803  G4double ld2 = ld*ld;
804  fElasticXsc = cofLogE*ld2 + 2.23;
805  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
806  }
807  else
808  {
809  G4double lh = pLab - 0.98;
810  G4double hd = lh*lh + .021;
811  G4double sqrLogPlab = logP*logP;
812 
813  fElasticXsc = 5.0 + 8.1*G4Exp(-logP*1.8) + 0.16*sqrLogPlab
814  - 1.3*logP + .15/hd;
815  fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2.9*logP + 0.30/hd;
816  }
817  }
818  else if( (theParticle == theKPlus) && proton ) // K+p
819  {
820  // VI: modified low-energy part
821  if( pLab < 0.631 )
822  {
823  fElasticXsc = fTotalXsc = 12.03;
824  }
825  else if( pLab > pMax )
826  {
827  G4double ld = logP - minLogP;
828  G4double ld2 = ld*ld;
829  fElasticXsc = cofLogE*ld2 + 2.23;
830  fTotalXsc = cofLogT*ld2 + 19.2;
831  }
832  else
833  {
834  G4double ld = logP - minLogP;
835  G4double ld2 = ld*ld;
836  G4double lr = pLab - .38;
837  G4double LE = .7/(lr*lr + .076);
838  G4double sp = std::sqrt(pLab);
839  G4double p2 = pLab*pLab;
840  G4double p4 = p2*p2;
841  // VI: tuned elastic
842  fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4)
843  + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652);
844  fTotalXsc = LE + (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4)
845  + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392);
846  }
847  }
848  else if( (theParticle == theKPlus) && neutron) // K+n
849  {
850  if( pLab < pMin )
851  {
852  G4double lm = pLab - 0.94;
853  G4double md = lm*lm + .392;
854  fElasticXsc = 2./md;
855  fTotalXsc = 4.6/md;
856  }
857  else if( pLab > pMax )
858  {
859  G4double ld = logP - minLogP;
860  G4double ld2 = ld*ld;
861  fElasticXsc = cofLogE*ld2 + 2.23;
862  fTotalXsc = cofLogT*ld2 + 19.2;
863  }
864  else
865  {
866  G4double ld = logP - minLogP;
867  G4double ld2 = ld*ld;
868  G4double sp = std::sqrt(pLab);
869  G4double p2 = pLab*pLab;
870  G4double p4 = p2*p2;
871  G4double lm = pLab - 0.94;
872  G4double md = lm*lm + .392;
873  fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
874  fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 4.6/md;
875  }
876  }
877  fTotalXsc *= CLHEP::millibarn;
878  fElasticXsc *= CLHEP::millibarn;
879  fElasticXsc = std::min(fElasticXsc, fTotalXsc);
880 
881  if( proton && theParticle->GetPDGCharge() > 0. && ekin < 100*MeV)
882  {
883  G4double cB = G4NuclearRadii::CoulombFactor(theParticle, nucleon, ekin);
884  fTotalXsc *= cB;
885  fElasticXsc *= cB;
886  }
887  fInelasticXsc = std::max(fTotalXsc - fElasticXsc,0.0);
888  /*
889  G4cout<< "HNXsc: Ekin(GeV)= " << ekin/GeV << "; tot(mb)= " << fTotalXsc/millibarn
890  <<"; el(mb)= " <<fElasticXsc/millibarn
891  <<"; inel(mb)= " <<fInelasticXsc/millibarn<<" "
892  << theParticle->GetParticleName() << " + "
893  << nucleon->GetParticleName() << G4endl;
894  */
895  return fTotalXsc;
896 }
897 
899 //
900 // Returns kaon-nucleon cross-section based on smoothed NS
901 // tuned for the Glauber-Gribov hadron model for Z>1
902 
903 G4double G4HadronNucleonXsc::KaonNucleonXscGG(
904  const G4ParticleDefinition* theParticle,
905  const G4ParticleDefinition* nucleon, G4double ekin)
906 {
907  fTotalXsc = fElasticXsc = fInelasticXsc = 0.0;
908  if(theParticle == theKMinus || theParticle == theKPlus) {
909  KaonNucleonXscVG(theParticle, nucleon, ekin);
910 
911  } else if(theParticle == theK0S || theParticle == theK0L) {
912  G4double stot = KaonNucleonXscVG(theKMinus, nucleon, ekin);
913  G4double sel = fElasticXsc;
914  G4double sinel = fInelasticXsc;
915  stot += KaonNucleonXscVG(theKPlus, nucleon, ekin);
916  sel += fElasticXsc;
917  sinel += fInelasticXsc;
918  fTotalXsc = stot*0.5;
919  fElasticXsc = sel*0.5;
920  fInelasticXsc = sinel*0.5;
921  }
922  return fTotalXsc;
923 }
924 
926 //
927 // Returns kaon-nucleon cross-section using NS
928 
929 G4double G4HadronNucleonXsc::KaonNucleonXscNS(
930  const G4ParticleDefinition* theParticle,
931  const G4ParticleDefinition* nucleon, G4double ekin)
932 {
933  fTotalXsc = fElasticXsc = fInelasticXsc = 0.0;
934  if(theParticle == theKMinus || theParticle == theKPlus) {
935  HadronNucleonXscNS(theParticle, nucleon, ekin);
936 
937  } else if(theParticle == theK0S || theParticle == theK0L) {
938  G4double stot = HadronNucleonXscNS(theKMinus, nucleon, ekin);
939  G4double sel = fElasticXsc;
940  G4double sinel = fInelasticXsc;
941  stot += HadronNucleonXscNS(theKPlus, nucleon, ekin);
942  sel += fElasticXsc;
943  sinel += fInelasticXsc;
944  fTotalXsc = stot*0.5;
945  fElasticXsc = sel*0.5;
946  fInelasticXsc = sinel*0.5;
947  }
948  return fTotalXsc;
949 }
950 
952 //
953 // Returns kaon-nucleon cross-section with smoothed NS parameterisation
954 
955 G4double G4HadronNucleonXsc::KaonNucleonXscVG(
956  const G4ParticleDefinition* theParticle,
957  const G4ParticleDefinition* nucleon, G4double ekin)
958 {
959  G4double pM = theParticle->GetPDGMass();
960  G4double pLab = std::sqrt(ekin*(ekin + 2*pM));
961 
962  pLab *= invGeV;
963  G4double logP = G4Log(pLab);
964 
965  fTotalXsc = 0.0;
966 
967  G4bool proton = (nucleon == theProton);
968  G4bool neutron = (nucleon == theNeutron);
969 
970  if( (theParticle == theKMinus) && proton ) // K-p
971  {
972  if( pLab < pMin)
973  {
974  G4double psp = pLab*std::sqrt(pLab);
975  fElasticXsc = 5.2/psp;
976  fTotalXsc = 14./psp;
977  }
978  else if( pLab > pMax )
979  {
980  G4double ld = logP - minLogP;
981  G4double ld2 = ld*ld;
982  fElasticXsc = cofLogE*ld2 + 2.23;
983  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
984  }
985  else
986  {
987  G4double ld = logP - minLogP;
988  G4double ld2 = ld*ld;
989  G4double sp = std::sqrt(pLab);
990  G4double psp = pLab*sp;
991  G4double p2 = pLab*pLab;
992  G4double p4 = p2*p2;
993 
994  G4double lh = pLab - 1.01;
995  G4double hd = lh*lh + .011;
996  fElasticXsc = 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .15/hd;
997  fTotalXsc = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4) + .60/hd;
998  }
999  }
1000  else if( (theParticle == theKMinus) && neutron ) // K-n
1001  {
1002  if( pLab > pMax )
1003  {
1004  G4double ld = logP - minLogP;
1005  G4double ld2 = ld*ld;
1006  fElasticXsc = cofLogE*ld2 + 2.23;
1007  fTotalXsc = 1.1*cofLogT*ld2 + 19.7;
1008  }
1009  else
1010  {
1011  G4double lh = pLab - 0.98;
1012  G4double hd = lh*lh + .045; // vg version
1013  G4double sqrLogPlab = logP*logP;
1014 
1015  fElasticXsc = 5.0 + 8.1*G4Exp(-logP*1.8) + 0.16*sqrLogPlab
1016  - 1.3*logP + .15/hd;
1017  fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2.9*logP + 0.60/hd; // vg version
1018  }
1019  }
1020  else if( (theParticle == theKPlus) && proton ) // K+p
1021  {
1022  // VI: modified low-energy part
1023  if( pLab < 0.631 )
1024  {
1025  fElasticXsc = fTotalXsc = 12.03;
1026  }
1027  else if( pLab > pMax )
1028  {
1029  G4double ld = logP - minLogP;
1030  G4double ld2 = ld*ld;
1031  fElasticXsc = cofLogE*ld2 + 2.23;
1032  fTotalXsc = cofLogT*ld2 + 19.2;
1033  }
1034  else
1035  {
1036  G4double ld = logP - minLogP;
1037  G4double ld2 = ld*ld;
1038  G4double lr = pLab - .38;
1039  G4double LE = .7/(lr*lr + .076);
1040  G4double sp = std::sqrt(pLab);
1041  G4double p2 = pLab*pLab;
1042  G4double p4 = p2*p2;
1043  // VI: tuned elastic
1044  fElasticXsc = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4)
1045  + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652);
1046  fTotalXsc = LE + (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4)
1047  + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392);
1048  }
1049  }
1050  else if( (theParticle == theKPlus) && neutron) // K+n
1051  {
1052  if( pLab < pMin )
1053  {
1054  G4double lm = pLab - 0.94;
1055  G4double md = lm*lm + .392;
1056  fElasticXsc = 2./md;
1057  fTotalXsc = 4.6/md;
1058  }
1059  else if( pLab > pMax )
1060  {
1061  G4double ld = logP - minLogP;
1062  G4double ld2 = ld*ld;
1063  fElasticXsc = cofLogE*ld2 + 2.23;
1064  fTotalXsc = cofLogT*ld2 + 19.2;
1065  }
1066  else
1067  {
1068  G4double ld = logP - minLogP;
1069  G4double ld2 = ld*ld;
1070  G4double sp = std::sqrt(pLab);
1071  G4double p2 = pLab*pLab;
1072  G4double p4 = p2*p2;
1073  G4double lm = pLab - 0.94;
1074  G4double md = lm*lm + .392;
1075  fElasticXsc = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
1076  fTotalXsc = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 4.6/md;
1077  }
1078  }
1079 
1080  fTotalXsc *= CLHEP::millibarn;
1081  fElasticXsc *= CLHEP::millibarn;
1082 
1083  if( proton && theParticle->GetPDGCharge() > 0. )
1084  {
1085  G4double cB = G4NuclearRadii::CoulombFactor(theParticle, nucleon, ekin);
1086  fTotalXsc *= cB;
1087  fElasticXsc *= cB;
1088  }
1089  fElasticXsc = std::min(fElasticXsc, fTotalXsc);
1090  fInelasticXsc = std::max(fTotalXsc - fElasticXsc,0.0);
1091  /*
1092  G4cout << "HNXscVG: E= " << ekin << " " << theParticle->GetParticleName()
1093  << " P: " << proton << " xtot(b)= " << fTotalXsc/barn
1094  << " xel(b)= " << fElasticXsc/barn << " xinel(b)= " << fInelasticXsc/barn
1095  << G4endl;
1096  */
1097  return fTotalXsc;
1098 }
1099 
1101 //
1102 // Returns hyperon-nucleon cross-section using NS x-section for protons
1103 
1104 G4double G4HadronNucleonXsc::HyperonNucleonXscNS(
1105  const G4ParticleDefinition* theParticle,
1106  const G4ParticleDefinition* nucleon, G4double ekin)
1107 {
1108  G4double coeff = 1.0;
1109 
1110  static const G4double lBarCof1S = 0.88;
1111  static const G4double lBarCof2S = 0.76;
1112  static const G4double lBarCof3S = 0.64;
1113  static const G4double lBarCof1C = 0.784378;
1114  static const G4double lBarCofSC = 0.664378;
1115  static const G4double lBarCof2SC = 0.544378;
1116  static const G4double lBarCof1B = 0.740659;
1117  static const G4double lBarCofSB = 0.620659;
1118  static const G4double lBarCof2SB = 0.500659;
1119 
1120  if( theParticle == theL || theParticle == theSPlus ||
1121  theParticle == theSMinus || theParticle == theS0 ||
1122  theParticle == theAntiL || theParticle == theASPlus ||
1123  theParticle == theASMinus || theParticle == theAS0 )
1124  {
1125  coeff = lBarCof1S;
1126 
1127  } else if( theParticle == theXiMinus || theParticle == theXi0 ||
1128  theParticle == theAXiMinus || theParticle == theAXi0 )
1129  {
1130  coeff = lBarCof2S;
1131  }
1132  else if( theParticle == theOmega || theParticle == theAOmega)
1133  {
1134  coeff = lBarCof3S;
1135  }
1136  else if( theParticle == theLambdaCPlus || theParticle == theALambdaCPlus ||
1137  theParticle == theSigmaCPlus || theParticle == theASigmaCPlus ||
1138  theParticle == theSigmacPP || theParticle == theASigmacPP ||
1139  theParticle == theSigmaC0 || theParticle == theASigmaC0
1140  )
1141  {
1142  coeff = lBarCof1C;
1143  }
1144  else if( theParticle == theOmegaC0 || theParticle == theAOmegaC0 )
1145  {
1146  coeff = lBarCof2SC;
1147  }
1148  else if( theParticle == theXiCPlus || theParticle == theXiC0 ||
1149  theParticle == theAXiCPlus || theParticle == theAXiC0)
1150  {
1151  coeff = lBarCofSC;
1152  }
1153  else if( theParticle == theLambdaB || theParticle == theALambdaB ||
1154  theParticle == theSigmaBPlus || theParticle == theASigmaBPlus ||
1155  theParticle == theSigmaBMinus || theParticle == theASigmaBMinus ||
1156  theParticle == theSigmaB0 || theParticle == theASigmaB0
1157  )
1158  {
1159  coeff = lBarCof1B;
1160  }
1161  else if( theParticle == theOmegaBMinus || theParticle == theAOmegaBMinus)
1162  {
1163  coeff = lBarCof2SB;
1164  }
1165  else if( theParticle == theXiBMinus || theParticle == theXiB0 ||
1166  theParticle == theAXiBMinus || theParticle == theAXiB0)
1167  {
1168  coeff = lBarCofSB;
1169  }
1170  fTotalXsc = coeff*HadronNucleonXscNS( theProton, nucleon, ekin);
1171  fInelasticXsc *= coeff;
1172  fElasticXsc *= coeff;
1173 
1174  return fTotalXsc;
1175 }
1176 
1178 //
1179 // Returns hyperon-nucleon cross-section using NS x-section for protons
1180 
1181 G4double G4HadronNucleonXsc::SCBMesonNucleonXscNS(
1182  const G4ParticleDefinition* theParticle,
1183  const G4ParticleDefinition* nucleon, G4double ekin )
1184 {
1185  G4double coeff(1.0);
1186  // static const G4double lMesCof1S = 0.82; // Kp/piP
1187  static const G4double llMesCof1C = 0.676568;
1188  static const G4double llMesCof1B = 0.610989;
1189  static const G4double llMesCof2C = 0.353135;
1190  static const G4double llMesCof2B = 0.221978;
1191  static const G4double llMesCofSC = 0.496568;
1192  static const G4double llMesCofSB = 0.430989;
1193  static const G4double llMesCofCB = 0.287557;
1194  static const G4double llMesCofEtaP = 0.88;
1195  static const G4double llMesCofEta = 0.76;
1196 
1197  if( theParticle == theBMeson0 || theParticle == theABMeson0 ||
1198  theParticle == theBMesonPlus || theParticle == theBMesonMinus )
1199  {
1200  coeff = llMesCof1B;
1201  }
1202  else if(theParticle == theDMeson0 || theParticle == theADMeson0 ||
1203  theParticle == theDMesonPlus || theParticle == theDMesonMinus )
1204  {
1205  coeff = llMesCof1C;
1206  }
1207  else if(theParticle == theBsMeson0 || theParticle == theABsMeson0 )
1208  {
1209  coeff = llMesCofSB;
1210  }
1211  else if(theParticle == theBcMesonPlus || theParticle == theBcMesonMinus )
1212  {
1213  coeff = llMesCofCB;
1214  }
1215  else if(theParticle == theDsMesonPlus || theParticle == theDsMesonMinus )
1216  {
1217  coeff = llMesCofSC;
1218  }
1219  else if(theParticle == theBMesonPlus || theParticle == theBMesonMinus )
1220  {
1221  coeff = llMesCof1B;
1222  }
1223  else if(theParticle == theDMesonPlus || theParticle == theDMesonMinus )
1224  {
1225  coeff = llMesCof1C;
1226  }
1227  else if(theParticle == theEtaC || theParticle == theJPsi )
1228  {
1229  coeff = llMesCof2C;
1230  }
1231  else if(theParticle == theUpsilon )
1232  {
1233  coeff = llMesCof2B;
1234  }
1235  else if(theParticle == theEta )
1236  {
1237  coeff = llMesCofEta;
1238  }
1239  else if(theParticle == theEtaPrime )
1240  {
1241  coeff = llMesCofEtaP;
1242  }
1243  fTotalXsc = coeff*HadronNucleonXscNS( thePiPlus, nucleon, ekin);
1244  fElasticXsc *= coeff;
1245  fInelasticXsc *= coeff;
1246  return fTotalXsc;
1247 }
1249 //
1250 // Returns hadron-nucleon cross-section based on V. Uzjinsky parametrisation of
1251 // data from G4FTFCrossSection class
1252 
1253 G4double G4HadronNucleonXsc::HadronNucleonXscVU(
1254  const G4ParticleDefinition* theParticle,
1255  const G4ParticleDefinition* nucleon, G4double ekin)
1256 {
1257  G4int PDGcode = theParticle->GetPDGEncoding();
1258  G4int absPDGcode = std::abs(PDGcode);
1259  G4double mass = theParticle->GetPDGMass();
1260  G4double Elab = ekin + mass;
1261  // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
1262  G4double Plab = std::sqrt(ekin*(ekin + 2.*mass));
1263 
1264  Elab *= invGeV;
1265  Plab *= invGeV;
1266 
1267  G4double logPlab = G4Log( Plab );
1268  G4double sqrLogPlab = logPlab * logPlab;
1269 
1270  G4bool proton = (nucleon == theProton);
1271  G4bool neutron = (nucleon == theNeutron);
1272 
1273  if( absPDGcode > 1000) //------Projectile is baryon -
1274  {
1275  if(proton)
1276  {
1277  fTotalXsc = 48.0 + 0.522*sqrLogPlab - 4.51*logPlab;
1278  fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1279  }
1280  if(neutron)
1281  {
1282  fTotalXsc = 47.3 + 0.513*sqrLogPlab - 4.27*logPlab;
1283  fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1284  }
1285  }
1286  else if( PDGcode == 211) //------Projectile is PionPlus ----
1287  {
1288  if(proton)
1289  {
1290  fTotalXsc = 16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab;
1291  fElasticXsc = 0.0 + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab;
1292  }
1293  if(neutron)
1294  {
1295  fTotalXsc = 33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab;
1296  fElasticXsc = 1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab;
1297  }
1298  }
1299  else if( PDGcode == -211) //------Projectile is PionMinus ----
1300  {
1301  if(proton)
1302  {
1303  fTotalXsc = 33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab;
1304  fElasticXsc = 1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab;
1305  }
1306  if(neutron)
1307  {
1308  fTotalXsc = 16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab;
1309  fElasticXsc = 0.0 + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab;
1310  }
1311  }
1312  else if( PDGcode == 111) //------Projectile is PionZero --
1313  {
1314  if(proton)
1315  {
1316  fTotalXsc = (16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab + //Pi+
1317  33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab)/2; //Pi-
1318 
1319  fElasticXsc = (0.0 + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab + //Pi+
1320  1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab)/2;//Pi-
1321 
1322  }
1323  if(neutron)
1324  {
1325  fTotalXsc = (33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab + //Pi+
1326  16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab)/2; //Pi-
1327  fElasticXsc = (1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab + //Pi+
1328  0.0 + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab)/2;//Pi-
1329  }
1330  }
1331  else if( PDGcode == 321 ) //------Projectile is KaonPlus --
1332  {
1333  if(proton)
1334  {
1335  fTotalXsc = 18.1 + 0.26 *sqrLogPlab - 1.0 *logPlab;
1336  fElasticXsc = 5.0 + 8.1*G4Exp(-logPlab*1.8 ) + 0.16 *sqrLogPlab - 1.3 *logPlab;
1337  }
1338  if(neutron)
1339  {
1340  fTotalXsc = 18.7 + 0.21 *sqrLogPlab - 0.89*logPlab;
1341  fElasticXsc = 7.3 + 0.29 *sqrLogPlab - 2.4 *logPlab;
1342  }
1343  }
1344  else if( PDGcode ==-321 ) //------Projectile is KaonMinus ----
1345  {
1346  if(proton)
1347  {
1348  fTotalXsc = 32.1 + 0.66*sqrLogPlab - 5.6*logPlab;
1349  fElasticXsc = 7.3 + 0.29*sqrLogPlab - 2.4*logPlab;
1350  }
1351  if(neutron)
1352  {
1353  fTotalXsc = 25.2 + 0.38*sqrLogPlab - 2.9*logPlab;
1354  fElasticXsc = 5.0 + 8.1*G4Exp(-logPlab*1.8 ) + 0.16*sqrLogPlab - 1.3*logPlab;
1355  }
1356  }
1357  else if( PDGcode == 311 ) //------Projectile is KaonZero -----
1358  {
1359  if(proton)
1360  {
1361  fTotalXsc = ( 18.1 + 0.26 *sqrLogPlab - 1.0 *logPlab + //K+
1362  32.1 + 0.66 *sqrLogPlab - 5.6 *logPlab)/2; //K-
1363  fElasticXsc = ( 5.0 + 8.1*G4Exp(-logPlab*1.8 ) + 0.16 *sqrLogPlab - 1.3 *logPlab + //K+
1364  7.3 + 0.29 *sqrLogPlab - 2.4 *logPlab)/2; //K-
1365  }
1366  if(neutron)
1367  {
1368  fTotalXsc = ( 18.7 + 0.21 *sqrLogPlab - 0.89*logPlab + //K+
1369  25.2 + 0.38 *sqrLogPlab - 2.9 *logPlab)/2; //K-
1370  fElasticXsc = ( 7.3 + 0.29 *sqrLogPlab - 2.4 *logPlab + //K+
1371  5.0 + 8.1*G4Exp(-logPlab*1.8 ) + 0.16 *sqrLogPlab - 1.3 *logPlab)/2; //K-
1372  }
1373  }
1374  else //------Projectile is undefined, Nucleon assumed
1375  {
1376  if(proton)
1377  {
1378  fTotalXsc = 48.0 + 0.522*sqrLogPlab - 4.51*logPlab;
1379  fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1380  }
1381  if(neutron)
1382  {
1383  fTotalXsc = 47.3 + 0.513*sqrLogPlab - 4.27*logPlab;
1384  fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1385  }
1386  }
1387 
1388  fTotalXsc *= CLHEP::millibarn;
1389  fElasticXsc *= CLHEP::millibarn;
1390  fElasticXsc = std::min(fElasticXsc, fTotalXsc);
1391  fInelasticXsc = fTotalXsc - fElasticXsc;
1392 
1393  return fTotalXsc;
1394 }
1395 
1397 //
1398 // Returns hadron-nucleon Xsc according to differnt parametrisations:
1399 // [2] E. Levin, hep-ph/9710546
1400 // [3] U. Dersch, et al, hep-ex/9910052
1401 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
1402 
1403 G4double G4HadronNucleonXsc::HadronNucleonXscEL(
1404  const G4ParticleDefinition* theParticle,
1405  const G4ParticleDefinition*, G4double ekin)
1406 {
1407  G4double xsection(0.);
1408  static const G4double targ_mass =
1409  0.5*(CLHEP::proton_mass_c2 + CLHEP::neutron_mass_c2);
1410  G4double sMand =
1411  CalcMandelstamS(ekin, theParticle->GetPDGMass(), targ_mass)*invGeV2;
1412 
1413  G4double x1 = G4Exp(G4Log(sMand)*0.0808);
1414  G4double x2 = G4Exp(G4Log(-sMand)*0.4525);
1415 
1416  if(theParticle == theGamma)
1417  {
1418  xsection = 0.0677*x1 + 0.129*x2;
1419  }
1420  else if(theParticle == theNeutron)
1421  {
1422  xsection = 21.70*x1 + 56.08*x2;
1423  }
1424  else if(theParticle == theProton)
1425  {
1426  xsection = 21.70*x1 + 56.08*x2;
1427  }
1428  else if(theParticle == theAProton)
1429  {
1430  xsection = 21.70*x1 + 98.39*x2;
1431  }
1432  else if(theParticle == thePiPlus)
1433  {
1434  xsection = 13.63*x1 + 27.56*x2;
1435  }
1436  else if(theParticle == thePiMinus)
1437  {
1438  xsection = 13.63*x1 + 36.02*x2;
1439  }
1440  else if(theParticle == theKPlus)
1441  {
1442  xsection = 11.82*x1 + 8.15*x2;
1443  }
1444  else if(theParticle == theKMinus)
1445  {
1446  xsection = 11.82*x1 + 26.36*x2;
1447  }
1448  else if(theParticle == theK0S || theParticle == theK0L)
1449  {
1450  xsection = 11.82*x1 + 17.25*x2;
1451  }
1452  else
1453  {
1454  xsection = 21.70*x1 + 56.08*x2;
1455  }
1456  fTotalXsc = xsection*CLHEP::millibarn;
1457  fInelasticXsc = 0.83*fTotalXsc;
1458  fElasticXsc = fTotalXsc - fInelasticXsc;
1459  return fTotalXsc;
1460 }
1461 
1463 
1464 G4double
1465 G4HadronNucleonXsc::CoulombBarrier(const G4ParticleDefinition* theParticle,
1466  const G4ParticleDefinition* nucleon,
1467  G4double ekin)
1468 {
1469  G4double tR = 0.895*CLHEP::fermi;
1470  G4double pR = 0.5*CLHEP::fermi;
1471 
1472  if ( theParticle == theProton ) pR = 0.895*fermi;
1473  else if( theParticle == thePiPlus ) pR = 0.663*fermi;
1474  else if( theParticle == theKPlus ) pR = 0.340*fermi;
1475 
1476  G4double pZ = theParticle->GetPDGCharge();
1477  G4double tZ = nucleon->GetPDGCharge();
1478 
1479  G4double pM = theParticle->GetPDGMass();
1480  G4double tM = nucleon->GetPDGMass();
1481 
1482  G4double pElab = ekin + pM;
1483 
1484  G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
1485 
1486  G4double totTcm = totEcm - pM -tM;
1487 
1488  G4double bC = fine_structure_const*hbarc*pZ*tZ/(2.*(pR + tR));
1489 
1490  //G4cout<<"##CB: ekin = "<<ekin<<"; pElab= " << pElab
1491  // <<"; bC = "<<bC<<"; Tcm = "<< totTcm <<G4endl;
1492 
1493  G4double ratio = (totTcm > bC) ? 1. - bC/totTcm : 0.0;
1494 
1495  // G4cout <<"ratio = "<<ratio<<G4endl;
1496  return ratio;
1497 }
1498 
1500