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