Geant4_10
G4ComponentGGHadronNucleusXsc.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 // author: V. Grichine
27 //
28 // 25.04.12 V. Grichine - first implementation
29 
30 #include "G4ComponentGGHadronNucleusXsc.hh"
31 
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4ParticleTable.hh"
35 #include "G4IonTable.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4DynamicParticle.hh"
38 #include "G4HadronNucleonXsc.hh"
39 
40 
42 //
43 
44 G4ComponentGGHadronNucleusXsc::G4ComponentGGHadronNucleusXsc()
45  : G4VComponentCrossSection("Glauber-Gribov"),
46 // fUpperLimit(100000*GeV),
47  fLowerLimit(10.*MeV),// fLowerLimit(3*GeV),
48  fRadiusConst(1.08*fermi), // 1.1, 1.3 ?
49  fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
50  fDiffractionXsc(0.0)
51 // , fHadronNucleonXsc(0.0)
52 {
53  theGamma = G4Gamma::Gamma();
54  theProton = G4Proton::Proton();
55  theNeutron = G4Neutron::Neutron();
56  theAProton = G4AntiProton::AntiProton();
57  theANeutron = G4AntiNeutron::AntiNeutron();
58  thePiPlus = G4PionPlus::PionPlus();
59  thePiMinus = G4PionMinus::PionMinus();
60  thePiZero = G4PionZero::PionZero();
61  theKPlus = G4KaonPlus::KaonPlus();
62  theKMinus = G4KaonMinus::KaonMinus();
63  theK0S = G4KaonZeroShort::KaonZeroShort();
64  theK0L = G4KaonZeroLong::KaonZeroLong();
65  theL = G4Lambda::Lambda();
66  theAntiL = G4AntiLambda::AntiLambda();
67  theSPlus = G4SigmaPlus::SigmaPlus();
68  theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
69  theSMinus = G4SigmaMinus::SigmaMinus();
70  theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
71  theS0 = G4SigmaZero::SigmaZero();
72  theAS0 = G4AntiSigmaZero::AntiSigmaZero();
73  theXiMinus = G4XiMinus::XiMinus();
74  theXi0 = G4XiZero::XiZero();
75  theAXiMinus = G4AntiXiMinus::AntiXiMinus();
76  theAXi0 = G4AntiXiZero::AntiXiZero();
77  theOmega = G4OmegaMinus::OmegaMinus();
78  theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
79  theD = G4Deuteron::Deuteron();
80  theT = G4Triton::Triton();
81  theA = G4Alpha::Alpha();
82  theHe3 = G4He3::He3();
83 
84  hnXsc = new G4HadronNucleonXsc();
85 }
86 
88 //
89 //
90 
91 G4ComponentGGHadronNucleusXsc::~G4ComponentGGHadronNucleusXsc()
92 {
93  if (hnXsc) delete hnXsc;
94 }
95 
97 
98 G4double G4ComponentGGHadronNucleusXsc::GetTotalIsotopeCrossSection(const G4ParticleDefinition* aParticle,
99  G4double kinEnergy,
100  G4int Z, G4int A)
101 {
102  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
103  kinEnergy);
104  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
105  delete aDP;
106 
107  return fTotalXsc;
108 }
109 
111 
112 G4double G4ComponentGGHadronNucleusXsc::GetTotalElementCrossSection(const G4ParticleDefinition* aParticle,
113  G4double kinEnergy,
114  G4int Z, G4double A)
115 {
116  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
117  kinEnergy);
118  fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
119  delete aDP;
120 
121  return fTotalXsc;
122 }
123 
125 
126 G4double G4ComponentGGHadronNucleusXsc::GetInelasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
127  G4double kinEnergy,
128  G4int Z, G4int A)
129 {
130  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
131  kinEnergy);
132  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
133  delete aDP;
134 
135  return fInelasticXsc;
136 }
137 
139 
140 G4double G4ComponentGGHadronNucleusXsc::GetProductionIsotopeCrossSection(const G4ParticleDefinition* aParticle,
141  G4double kinEnergy,
142  G4int Z, G4int A)
143 {
144  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
145  kinEnergy);
146  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
147  delete aDP;
148 
149  return fProductionXsc;
150 }
151 
153 
154 G4double G4ComponentGGHadronNucleusXsc::GetInelasticElementCrossSection(const G4ParticleDefinition* aParticle,
155  G4double kinEnergy,
156  G4int Z, G4double A)
157 {
158  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
159  kinEnergy);
160  fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
161  delete aDP;
162 
163  return fInelasticXsc;
164 }
165 
167 
168 G4double G4ComponentGGHadronNucleusXsc::GetProductionElementCrossSection(const G4ParticleDefinition* aParticle,
169  G4double kinEnergy,
170  G4int Z, G4double A)
171 {
172  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
173  kinEnergy);
174  fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
175  delete aDP;
176 
177  return fProductionXsc;
178 }
179 
181 
182 G4double G4ComponentGGHadronNucleusXsc::GetElasticElementCrossSection(const G4ParticleDefinition* aParticle,
183  G4double kinEnergy,
184  G4int Z, G4double A)
185 {
186  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
187  kinEnergy);
188  fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
189  delete aDP;
190 
191  return fElasticXsc;
192 }
193 
195 
196 G4double G4ComponentGGHadronNucleusXsc::GetElasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
197  G4double kinEnergy,
198  G4int Z, G4int A)
199 {
200  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
201  kinEnergy);
202  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
203  delete aDP;
204 
205  return fElasticXsc;
206 }
207 
209 
210 G4double G4ComponentGGHadronNucleusXsc::ComputeQuasiElasticRatio(const G4ParticleDefinition* aParticle,
211  G4double kinEnergy,
212  G4int Z, G4int A)
213 {
214  G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
215  kinEnergy);
216  fTotalXsc = GetIsoCrossSection(aDP, Z, A);
217  delete aDP;
218  G4double ratio = 0.;
219 
220  if(fInelasticXsc > 0.)
221  {
222  ratio = (fInelasticXsc - fProductionXsc)/fInelasticXsc;
223  if(ratio < 0.) ratio = 0.;
224  }
225  return ratio;
226 }
227 
228 
229 
230 
232 
233 G4bool
234 G4ComponentGGHadronNucleusXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
235  G4int Z, G4int /*A*/,
236  const G4Element*,
237  const G4Material*)
238 {
239  G4bool applicable = false;
240  // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
241  G4double kineticEnergy = aDP->GetKineticEnergy();
242 
243  const G4ParticleDefinition* theParticle = aDP->GetDefinition();
244 
245  if ( ( kineticEnergy >= fLowerLimit &&
246  Z > 1 && // >= He
247  ( theParticle == theAProton ||
248  theParticle == theGamma ||
249  theParticle == theKPlus ||
250  theParticle == theKMinus ||
251  theParticle == theK0L ||
252  theParticle == theK0S ||
253  theParticle == theSMinus ||
254  theParticle == theProton ||
255  theParticle == theNeutron ||
256  theParticle == thePiPlus ||
257  theParticle == thePiMinus ) ) ) applicable = true;
258 
259  return applicable;
260 }
261 
263 //
264 // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
265 // Glauber model with Gribov correction calculated in the dipole approximation on
266 // light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model.
267 // [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
268 
269 G4double
270 G4ComponentGGHadronNucleusXsc::GetIsoCrossSection(const G4DynamicParticle* aParticle,
271  G4int Z, G4int A,
272  const G4Isotope*,
273  const G4Element*,
274  const G4Material*)
275 {
276  G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
277  G4double hpInXsc(0.), hnInXsc(0.);
278  G4double R = GetNucleusRadius(A);
279 
280  G4int N = A - Z; // number of neutrons
281  if (N < 0) N = 0;
282 
283  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
284 
285  if( theParticle == theProton ||
286  theParticle == theNeutron ||
287  theParticle == thePiPlus ||
288  theParticle == thePiMinus )
289  {
290  // sigma = GetHadronNucleonXscNS(aParticle, A, Z);
291 
292  sigma = Z*hnXsc->GetHadronNucleonXscNS(aParticle, theProton);
293 
294  hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
295 
296  sigma += N*hnXsc->GetHadronNucleonXscNS(aParticle, theNeutron);
297 
298  hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
299 
300  cofInelastic = 2.4;
301  cofTotal = 2.0;
302  }
303  else if( theParticle == theKPlus ||
304  theParticle == theKMinus ||
305  theParticle == theK0S ||
306  theParticle == theK0L )
307  {
308  // sigma = GetKaonNucleonXscVector(aParticle, A, Z);
309 
310  sigma = Z*hnXsc->GetKaonNucleonXscGG(aParticle, theProton);
311 
312  hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
313 
314  sigma += N*hnXsc->GetKaonNucleonXscGG(aParticle, theNeutron);
315 
316  hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
317 
318  cofInelastic = 2.2;
319  cofTotal = 2.0;
320  R = 1.3*fermi;
321  R *= std::pow(G4double(A), 0.3333);
322  }
323  else
324  {
325  sigma = GetHadronNucleonXscNS(aParticle, A, Z);
326  cofInelastic = 2.2;
327  cofTotal = 2.0;
328  }
329  // cofInelastic = 2.0;
330 
331  if( A > 1 )
332  {
333  nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
334  ratio = sigma/nucleusSquare;
335 
336  xsection = nucleusSquare*std::log( 1. + ratio );
337 
338  xsection *= GetParticleBarCorTot(theParticle, Z);
339 
340  fTotalXsc = xsection;
341 
342 
343 
344  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
345 
346  fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);
347 
348  fElasticXsc = fTotalXsc - fInelasticXsc;
349 
350  if(fElasticXsc < 0.) fElasticXsc = 0.;
351 
352  G4double difratio = ratio/(1.+ratio);
353 
354  fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
355 
356 
357  // sigma = GetHNinelasticXsc(aParticle, A, Z);
358 
359  sigma = Z*hpInXsc + N*hnInXsc;
360 
361  ratio = sigma/nucleusSquare;
362 
363  fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
364 
365  fProductionXsc *= GetParticleBarCorIn(theParticle, Z);
366 
367  if (fElasticXsc < 0.) fElasticXsc = 0.;
368  }
369  else // H
370  {
371  fTotalXsc = sigma;
372  xsection = sigma;
373 
374  if ( theParticle != theAProton )
375  {
376  sigma = GetHNinelasticXsc(aParticle, A, Z);
377  fInelasticXsc = sigma;
378  fElasticXsc = fTotalXsc - fInelasticXsc;
379  }
380  else
381  {
382  fElasticXsc = fTotalXsc - fInelasticXsc;
383  }
384  if (fElasticXsc < 0.) fElasticXsc = 0.;
385 
386  }
387  return xsection;
388 }
389 
391 //
392 // Return single-diffraction/inelastic cross-section ratio
393 
394 G4double G4ComponentGGHadronNucleusXsc::
395 GetRatioSD(const G4DynamicParticle* aParticle, G4int A, G4int Z)
396 {
397  G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
398  G4double R = GetNucleusRadius(A);
399 
400  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
401 
402  if( theParticle == theProton ||
403  theParticle == theNeutron ||
404  theParticle == thePiPlus ||
405  theParticle == thePiMinus )
406  {
407  sigma = GetHadronNucleonXscNS(aParticle, A, Z);
408  cofInelastic = 2.4;
409  cofTotal = 2.0;
410  }
411  else
412  {
413  sigma = GetHadronNucleonXscNS(aParticle, A, Z);
414  cofInelastic = 2.2;
415  cofTotal = 2.0;
416  }
417  nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
418  ratio = sigma/nucleusSquare;
419 
420  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
421 
422  G4double difratio = ratio/(1.+ratio);
423 
424  fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
425 
426  if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
427  else ratio = 0.;
428 
429  return ratio;
430 }
431 
433 //
434 // Return suasi-elastic/inelastic cross-section ratio
435 
436 G4double G4ComponentGGHadronNucleusXsc::
437 GetRatioQE(const G4DynamicParticle* aParticle, G4int A, G4int Z)
438 {
439  G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
440  G4double R = GetNucleusRadius(A);
441 
442  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
443 
444  if( theParticle == theProton ||
445  theParticle == theNeutron ||
446  theParticle == thePiPlus ||
447  theParticle == thePiMinus )
448  {
449  sigma = GetHadronNucleonXscNS(aParticle, A, Z);
450  cofInelastic = 2.4;
451  cofTotal = 2.0;
452  }
453  else
454  {
455  sigma = GetHadronNucleonXscNS(aParticle, A, Z);
456  cofInelastic = 2.2;
457  cofTotal = 2.0;
458  }
459  nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
460  ratio = sigma/nucleusSquare;
461 
462  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
463 
464  sigma = GetHNinelasticXsc(aParticle, A, Z);
465  ratio = sigma/nucleusSquare;
466 
467  fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
468 
469  if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
470  else ratio = 0.;
471  if ( ratio < 0. ) ratio = 0.;
472 
473  return ratio;
474 }
475 
477 //
478 // Returns hadron-nucleon Xsc according to differnt parametrisations:
479 // [2] E. Levin, hep-ph/9710546
480 // [3] U. Dersch, et al, hep-ex/9910052
481 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
482 
483 G4double
484 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
485  const G4Element* anElement)
486 {
487  G4int At = G4lrint(anElement->GetN()); // number of nucleons
488  G4int Zt = G4lrint(anElement->GetZ()); // number of protons
489 
490  return GetHadronNucleonXsc(aParticle, At, Zt);
491 }
492 
494 //
495 // Returns hadron-nucleon Xsc according to differnt parametrisations:
496 // [2] E. Levin, hep-ph/9710546
497 // [3] U. Dersch, et al, hep-ex/9910052
498 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
499 
500 G4double
501 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
502  G4int At, G4int /*Zt*/)
503 {
504  G4double xsection;
505 
506  //G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
507 
508  G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
509 
510  G4double proj_mass = aParticle->GetMass();
511  G4double proj_momentum = aParticle->GetMomentum().mag();
512  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
513 
514  sMand /= GeV*GeV; // in GeV for parametrisation
515  proj_momentum /= GeV;
516 
517  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
518 
519  G4double aa = At;
520 
521  if(theParticle == theGamma)
522  {
523  xsection = aa*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
524  }
525  else if(theParticle == theNeutron) // as proton ???
526  {
527  xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
528  }
529  else if(theParticle == theProton)
530  {
531  xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
532  // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
533  // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
534  }
535  else if(theParticle == theAProton)
536  {
537  xsection = aa*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
538  }
539  else if(theParticle == thePiPlus)
540  {
541  xsection = aa*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
542  }
543  else if(theParticle == thePiMinus)
544  {
545  // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
546  xsection = aa*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
547  }
548  else if(theParticle == theKPlus)
549  {
550  xsection = aa*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
551  }
552  else if(theParticle == theKMinus)
553  {
554  xsection = aa*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
555  }
556  else // as proton ???
557  {
558  xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
559  }
560  xsection *= millibarn;
561  return xsection;
562 }
563 
564 
566 //
567 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
568 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
569 
570 G4double
571 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
572  const G4Element* anElement)
573 {
574  G4int At = G4lrint(anElement->GetN()); // number of nucleons
575  G4int Zt = G4lrint(anElement->GetZ()); // number of protons
576 
577  return GetHadronNucleonXscPDG(aParticle, At, Zt);
578 }
579 
580 
581 
582 
584 //
585 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
586 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
587 // At = number of nucleons, Zt = number of protons
588 
589 G4double
590 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
591  G4int At, G4int Zt)
592 {
593  G4double xsection;
594 
595  G4int Nt = At-Zt; // number of neutrons
596  if (Nt < 0) Nt = 0;
597 
598  G4double zz = Zt;
599  G4double aa = At;
600  G4double nn = Nt;
601 
602  G4double targ_mass = G4ParticleTable::GetParticleTable()->
603  GetIonTable()->GetIonMass(Zt, At);
604 
605  targ_mass = 0.939*GeV; // ~mean neutron and proton ???
606 
607  G4double proj_mass = aParticle->GetMass();
608  G4double proj_momentum = aParticle->GetMomentum().mag();
609 
610  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
611 
612  sMand /= GeV*GeV; // in GeV for parametrisation
613 
614  // General PDG fit constants
615 
616  G4double s0 = 5.38*5.38; // in Gev^2
617  G4double eta1 = 0.458;
618  G4double eta2 = 0.458;
619  G4double B = 0.308;
620 
621 
622  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
623 
624 
625  if(theParticle == theNeutron) // proton-neutron fit
626  {
627  xsection = zz*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
628  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
629  xsection += nn*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
630  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
631  }
632  else if(theParticle == theProton)
633  {
634 
635  xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
636  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
637 
638  xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
639  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
640  }
641  else if(theParticle == theAProton)
642  {
643  xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
644  + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
645 
646  xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
647  + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
648  }
649  else if(theParticle == thePiPlus)
650  {
651  xsection = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
652  + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
653  }
654  else if(theParticle == thePiMinus)
655  {
656  xsection = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
657  + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
658  }
659  else if(theParticle == theKPlus || theParticle == theK0L )
660  {
661  xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
662  + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
663 
664  xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
665  + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
666  }
667  else if(theParticle == theKMinus || theParticle == theK0S )
668  {
669  xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
670  + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
671 
672  xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
673  + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
674  }
675  else if(theParticle == theSMinus)
676  {
677  xsection = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
678  - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
679  }
680  else if(theParticle == theGamma) // modify later on
681  {
682  xsection = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
683  + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
684 
685  }
686  else // as proton ???
687  {
688  xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
689  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
690 
691  xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
692  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
693  }
694  xsection *= millibarn; // parametrised in mb
695  return xsection;
696 }
697 
698 
700 //
701 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
702 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
703 
704 G4double
705 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
706  const G4Element* anElement)
707 {
708  G4int At = G4lrint(anElement->GetN()); // number of nucleons
709  G4int Zt = G4lrint(anElement->GetZ()); // number of protons
710 
711  return GetHadronNucleonXscNS(aParticle, At, Zt);
712 }
713 
714 
715 
716 
718 //
719 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
720 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
721 
722 G4double
723 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
724  G4int At, G4int Zt)
725 {
726  G4double xsection(0);
727  // G4double Delta; DHW 19 May 2011: variable set but not used
728  G4double A0, B0;
729  G4double hpXscv(0);
730  G4double hnXscv(0);
731 
732  G4int Nt = At-Zt; // number of neutrons
733  if (Nt < 0) Nt = 0;
734 
735  G4double aa = At;
736  G4double zz = Zt;
737  G4double nn = Nt;
738 
739  G4double targ_mass = G4ParticleTable::GetParticleTable()->
740  GetIonTable()->GetIonMass(Zt, At);
741 
742  targ_mass = 0.939*GeV; // ~mean neutron and proton ???
743 
744  G4double proj_mass = aParticle->GetMass();
745  G4double proj_energy = aParticle->GetTotalEnergy();
746  G4double proj_momentum = aParticle->GetMomentum().mag();
747 
748  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
749 
750  sMand /= GeV*GeV; // in GeV for parametrisation
751  proj_momentum /= GeV;
752  proj_energy /= GeV;
753  proj_mass /= GeV;
754 
755  // General PDG fit constants
756 
757  G4double s0 = 5.38*5.38; // in Gev^2
758  G4double eta1 = 0.458;
759  G4double eta2 = 0.458;
760  G4double B = 0.308;
761 
762 
763  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
764 
765 
766  if(theParticle == theNeutron)
767  {
768  if( proj_momentum >= 373.)
769  {
770  return GetHadronNucleonXscPDG(aParticle,At,Zt);
771  }
772  else if( proj_momentum >= 10.)
773  // if( proj_momentum >= 2.)
774  {
775  // Delta = 1.; // DHW 19 May 2011: variable set but not used
776  // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
777 
778  if(proj_momentum >= 10.)
779  {
780  B0 = 7.5;
781  A0 = 100. - B0*std::log(3.0e7);
782 
783  xsection = A0 + B0*std::log(proj_energy) - 11
784  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
785  0.93827*0.93827,-0.165); // mb
786  }
787  xsection *= zz + nn;
788  }
789  else
790  {
791  // nn to be pp
792 
793  if( proj_momentum < 0.73 )
794  {
795  hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
796  }
797  else if( proj_momentum < 1.05 )
798  {
799  hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
800  (std::log(proj_momentum/0.73));
801  }
802  else // if( proj_momentum < 10. )
803  {
804  hnXscv = 39.0+
805  75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
806  }
807  // pn to be np
808 
809  if( proj_momentum < 0.8 )
810  {
811  hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
812  }
813  else if( proj_momentum < 1.4 )
814  {
815  hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
816  }
817  else // if( proj_momentum < 10. )
818  {
819  hpXscv = 33.3+
820  20.8*(std::pow(proj_momentum,2.0)-1.35)/
821  (std::pow(proj_momentum,2.50)+0.95);
822  }
823  xsection = hpXscv*zz + hnXscv*nn;
824  }
825  }
826  else if(theParticle == theProton)
827  {
828  if( proj_momentum >= 373.)
829  {
830  return GetHadronNucleonXscPDG(aParticle,At,Zt);
831  }
832  else if( proj_momentum >= 10.)
833  // if( proj_momentum >= 2.)
834  {
835  // Delta = 1.; DHW 19 May 2011: variable set but not used
836  // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
837 
838  if(proj_momentum >= 10.)
839  {
840  B0 = 7.5;
841  A0 = 100. - B0*std::log(3.0e7);
842 
843  xsection = A0 + B0*std::log(proj_energy) - 11
844  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
845  0.93827*0.93827,-0.165); // mb
846  }
847  xsection *= zz + nn;
848  }
849  else
850  {
851  // pp
852 
853  if( proj_momentum < 0.73 )
854  {
855  hpXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
856  }
857  else if( proj_momentum < 1.05 )
858  {
859  hpXscv = 23 + 40*(std::log(proj_momentum/0.73))*
860  (std::log(proj_momentum/0.73));
861  }
862  else // if( proj_momentum < 10. )
863  {
864  hpXscv = 39.0+
865  75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
866  }
867  // pn to be np
868 
869  if( proj_momentum < 0.8 )
870  {
871  hnXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
872  }
873  else if( proj_momentum < 1.4 )
874  {
875  hnXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
876  }
877  else // if( proj_momentum < 10. )
878  {
879  hnXscv = 33.3+
880  20.8*(std::pow(proj_momentum,2.0)-1.35)/
881  (std::pow(proj_momentum,2.50)+0.95);
882  }
883  xsection = hpXscv*zz + hnXscv*nn;
884  // xsection = hpXscv*(Zt + Nt);
885  // xsection = hnXscv*(Zt + Nt);
886  }
887  // xsection *= 0.95;
888  }
889  else if( theParticle == theAProton )
890  {
891  // xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
892  // + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
893 
894  // xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
895  // + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
896 
897  G4double logP = std::log(proj_momentum);
898 
899  if( proj_momentum <= 1.0 )
900  {
901  xsection = zz*(65.55 + 53.84/(proj_momentum+1.e-6) );
902  }
903  else
904  {
905  xsection = zz*( 41.1 + 77.2*std::pow( proj_momentum, -0.68)
906  + 0.293*logP*logP - 1.82*logP );
907  }
908  if ( nn > 0.)
909  {
910  xsection += nn*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
911  }
912  else // H
913  {
914  fInelasticXsc = 38.0 + 38.0*std::pow( proj_momentum, -0.96)
915  - 0.169*logP*logP;
916  fInelasticXsc *= millibarn;
917  }
918  }
919  else if( theParticle == thePiPlus )
920  {
921  if(proj_momentum < 0.4)
922  {
923  G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
924  hpXscv = Ex3+20.0;
925  }
926  else if( proj_momentum < 1.15 )
927  {
928  G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
929  hpXscv = Ex4+14.0;
930  }
931  else if(proj_momentum < 3.5)
932  {
933  G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
934  G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
935  hpXscv = Ex1+Ex2+27.5;
936  }
937  else // if(proj_momentum > 3.5) // mb
938  {
939  hpXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
940  }
941  // pi+n = pi-p??
942 
943  if(proj_momentum < 0.37)
944  {
945  hnXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
946  }
947  else if(proj_momentum<0.65)
948  {
949  hnXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
950  }
951  else if(proj_momentum<1.3)
952  {
953  hnXscv = 36.1+
954  10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
955  24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
956  }
957  else if(proj_momentum<3.0)
958  {
959  hnXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
960  3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
961  1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
962  }
963  else // mb
964  {
965  hnXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
966  }
967  xsection = hpXscv*zz + hnXscv*nn;
968  }
969  else if(theParticle == thePiMinus)
970  {
971  // pi-n = pi+p??
972 
973  if(proj_momentum < 0.4)
974  {
975  G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
976  hnXscv = Ex3+20.0;
977  }
978  else if(proj_momentum < 1.15)
979  {
980  G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
981  hnXscv = Ex4+14.0;
982  }
983  else if(proj_momentum < 3.5)
984  {
985  G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
986  G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
987  hnXscv = Ex1+Ex2+27.5;
988  }
989  else // if(proj_momentum > 3.5) // mb
990  {
991  hnXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
992  }
993  // pi-p
994 
995  if(proj_momentum < 0.37)
996  {
997  hpXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
998  }
999  else if(proj_momentum<0.65)
1000  {
1001  hpXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
1002  }
1003  else if(proj_momentum<1.3)
1004  {
1005  hpXscv = 36.1+
1006  10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
1007  24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
1008  }
1009  else if(proj_momentum<3.0)
1010  {
1011  hpXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
1012  3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
1013  1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
1014  }
1015  else // mb
1016  {
1017  hpXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
1018  }
1019  xsection = hpXscv*zz + hnXscv*nn;
1020  }
1021  else if(theParticle == theKPlus)
1022  {
1023  xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
1024  + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
1025 
1026  xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
1027  + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
1028  }
1029  else if(theParticle == theKMinus)
1030  {
1031  xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
1032  + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
1033 
1034  xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
1035  + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
1036  }
1037  else if(theParticle == theSMinus)
1038  {
1039  xsection = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
1040  - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
1041  }
1042  else if(theParticle == theGamma) // modify later on
1043  {
1044  xsection = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
1045  + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
1046 
1047  }
1048  else // as proton ???
1049  {
1050  xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
1051  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
1052 
1053  xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
1054  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
1055  }
1056  xsection *= millibarn; // parametrised in mb
1057  return xsection;
1058 }
1059 
1060 G4double
1061 G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector(const G4DynamicParticle* aParticle,
1062  G4int At, G4int Zt)
1063 {
1064  G4double Tkin, logTkin, xsc, xscP, xscN;
1065  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1066 
1067  G4int Nt = At-Zt; // number of neutrons
1068  if (Nt < 0) Nt = 0;
1069 
1070  Tkin = aParticle->GetKineticEnergy(); // Tkin in MeV
1071 
1072  if( Tkin > 70*GeV ) return GetHadronNucleonXscPDG(aParticle,At,Zt);
1073 
1074  logTkin = std::log(Tkin); // Tkin in MeV!!!
1075 
1076  if( theParticle == theKPlus )
1077  {
1078  xscP = hnXsc->GetKpProtonTotXscVector(logTkin);
1079  xscN = hnXsc->GetKpNeutronTotXscVector(logTkin);
1080  }
1081  else if( theParticle == theKMinus )
1082  {
1083  xscP = hnXsc->GetKmProtonTotXscVector(logTkin);
1084  xscN = hnXsc->GetKmNeutronTotXscVector(logTkin);
1085  }
1086  else // K-zero as half of K+ and K-
1087  {
1088  xscP = (hnXsc->GetKpProtonTotXscVector(logTkin)+hnXsc->GetKmProtonTotXscVector(logTkin))*0.5;
1089  xscN = (hnXsc->GetKpNeutronTotXscVector(logTkin)+hnXsc->GetKmNeutronTotXscVector(logTkin))*0.5;
1090  }
1091  xsc = xscP*Zt + xscN*Nt;
1092  return xsc;
1093 }
1095 //
1096 // Returns hadron-nucleon inelastic cross-section based on proper parametrisation
1097 
1098 G4double
1099 G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
1100  const G4Element* anElement)
1101 {
1102  G4int At = G4lrint(anElement->GetN()); // number of nucleons
1103  G4int Zt = G4lrint(anElement->GetZ()); // number of protons
1104 
1105  return GetHNinelasticXsc(aParticle, At, Zt);
1106 }
1107 
1109 //
1110 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
1111 
1112 G4double
1113 G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
1114  G4int At, G4int Zt)
1115 {
1116  G4ParticleDefinition* hadron = aParticle->GetDefinition();
1117  G4double sumInelastic;
1118  G4int Nt = At - Zt;
1119  if(Nt < 0) Nt = 0;
1120 
1121  if( hadron == theKPlus )
1122  {
1123  sumInelastic = GetHNinelasticXscVU(aParticle, At, Zt);
1124  }
1125  else
1126  {
1127  //sumInelastic = Zt*GetHadronNucleonXscMK(aParticle, theProton);
1128  // sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);
1129  sumInelastic = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
1130  sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);
1131  }
1132  return sumInelastic;
1133 }
1134 
1135 
1137 //
1138 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
1139 
1140 G4double
1141 G4ComponentGGHadronNucleusXsc::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
1142  G4int At, G4int Zt)
1143 {
1144  G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
1145  G4int absPDGcode = std::abs(PDGcode);
1146 
1147  G4double Elab = aParticle->GetTotalEnergy();
1148  // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
1149  G4double Plab = aParticle->GetMomentum().mag();
1150  // std::sqrt(Elab * Elab - 0.88);
1151 
1152  Elab /= GeV;
1153  Plab /= GeV;
1154 
1155  G4double LogPlab = std::log( Plab );
1156  G4double sqrLogPlab = LogPlab * LogPlab;
1157 
1158  //G4cout<<"Plab = "<<Plab<<G4endl;
1159 
1160  G4double NumberOfTargetProtons = G4double(Zt);
1161  G4double NumberOfTargetNucleons = G4double(At);
1162  G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
1163 
1164  if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;
1165 
1166  G4double Xtotal, Xelastic, Xinelastic;
1167 
1168  if( absPDGcode > 1000 ) //------Projectile is baryon --------
1169  {
1170  G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
1171  0.522*sqrLogPlab - 4.51*LogPlab;
1172 
1173  G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
1174  0.513*sqrLogPlab - 4.27*LogPlab;
1175 
1176  G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
1177  0.169*sqrLogPlab - 1.85*LogPlab;
1178 
1179  G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
1180  0.169*sqrLogPlab - 1.85*LogPlab;
1181 
1182  Xtotal = (NumberOfTargetProtons * XtotPP +
1183  NumberOfTargetNeutrons * XtotPN);
1184 
1185  Xelastic = (NumberOfTargetProtons * XelPP +
1186  NumberOfTargetNeutrons * XelPN);
1187  }
1188  else if( PDGcode == 211 ) //------Projectile is PionPlus -------
1189  {
1190  G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) +
1191  0.19 *sqrLogPlab - 0.0 *LogPlab;
1192 
1193  G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) +
1194  0.456*sqrLogPlab - 4.03*LogPlab;
1195 
1196  G4double XelPiP = 0.0 + 11.4*std::pow(Plab,-0.40) +
1197  0.079*sqrLogPlab - 0.0 *LogPlab;
1198 
1199  G4double XelPiN = 1.76 + 11.2*std::pow(Plab,-0.64) +
1200  0.043*sqrLogPlab - 0.0 *LogPlab;
1201 
1202  Xtotal = ( NumberOfTargetProtons * XtotPiP +
1203  NumberOfTargetNeutrons * XtotPiN );
1204 
1205  Xelastic = ( NumberOfTargetProtons * XelPiP +
1206  NumberOfTargetNeutrons * XelPiN );
1207  }
1208  else if( PDGcode == -211 ) //------Projectile is PionMinus -------
1209  {
1210  G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) +
1211  0.456*sqrLogPlab - 4.03*LogPlab;
1212 
1213  G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) +
1214  0.19 *sqrLogPlab - 0.0 *LogPlab;
1215 
1216  G4double XelPiP = 1.76 + 11.2*std::pow(Plab,-0.64) +
1217  0.043*sqrLogPlab - 0.0 *LogPlab;
1218 
1219  G4double XelPiN = 0.0 + 11.4*std::pow(Plab,-0.40) +
1220  0.079*sqrLogPlab - 0.0 *LogPlab;
1221 
1222  Xtotal = ( NumberOfTargetProtons * XtotPiP +
1223  NumberOfTargetNeutrons * XtotPiN );
1224 
1225  Xelastic = ( NumberOfTargetProtons * XelPiP +
1226  NumberOfTargetNeutrons * XelPiN );
1227  }
1228  else if( PDGcode == 111 ) //------Projectile is PionZero -------
1229  {
1230  G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) +
1231  0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
1232  33.0 + 14.0 *std::pow(Plab,-1.36) +
1233  0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
1234 
1235  G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) +
1236  0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
1237  16.4 + 19.3 *std::pow(Plab,-0.42) +
1238  0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1239 
1240  G4double XelPiP =( 0.0 + 11.4*std::pow(Plab,-0.40) +
1241  0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
1242  1.76 + 11.2*std::pow(Plab,-0.64) +
1243  0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1244 
1245  G4double XelPiN =( 1.76 + 11.2*std::pow(Plab,-0.64) +
1246  0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
1247  0.0 + 11.4*std::pow(Plab,-0.40) +
1248  0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
1249 
1250  Xtotal = ( NumberOfTargetProtons * XtotPiP +
1251  NumberOfTargetNeutrons * XtotPiN );
1252 
1253  Xelastic = ( NumberOfTargetProtons * XelPiP +
1254  NumberOfTargetNeutrons * XelPiN );
1255  }
1256  else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
1257  {
1258  G4double XtotKP = 18.1 + 0. *std::pow(Plab, 0. ) +
1259  0.26 *sqrLogPlab - 1.0 *LogPlab;
1260  G4double XtotKN = 18.7 + 0. *std::pow(Plab, 0. ) +
1261  0.21 *sqrLogPlab - 0.89*LogPlab;
1262 
1263  G4double XelKP = 5.0 + 8.1*std::pow(Plab,-1.8 ) +
1264  0.16 *sqrLogPlab - 1.3 *LogPlab;
1265 
1266  G4double XelKN = 7.3 + 0. *std::pow(Plab,-0. ) +
1267  0.29 *sqrLogPlab - 2.4 *LogPlab;
1268 
1269  Xtotal = ( NumberOfTargetProtons * XtotKP +
1270  NumberOfTargetNeutrons * XtotKN );
1271 
1272  Xelastic = ( NumberOfTargetProtons * XelKP +
1273  NumberOfTargetNeutrons * XelKN );
1274  }
1275  else if( PDGcode ==-321 ) //------Projectile is KaonMinus ------
1276  {
1277  G4double XtotKP = 32.1 + 0. *std::pow(Plab, 0. ) +
1278  0.66 *sqrLogPlab - 5.6 *LogPlab;
1279  G4double XtotKN = 25.2 + 0. *std::pow(Plab, 0. ) +
1280  0.38 *sqrLogPlab - 2.9 *LogPlab;
1281 
1282  G4double XelKP = 7.3 + 0. *std::pow(Plab,-0. ) +
1283  0.29 *sqrLogPlab - 2.4 *LogPlab;
1284 
1285  G4double XelKN = 5.0 + 8.1*std::pow(Plab,-1.8 ) +
1286  0.16 *sqrLogPlab - 1.3 *LogPlab;
1287 
1288  Xtotal = ( NumberOfTargetProtons * XtotKP +
1289  NumberOfTargetNeutrons * XtotKN );
1290 
1291  Xelastic = ( NumberOfTargetProtons * XelKP +
1292  NumberOfTargetNeutrons * XelKN );
1293  }
1294  else if( PDGcode == 311 ) //------Projectile is KaonZero ------
1295  {
1296  G4double XtotKP = ( 18.1 + 0. *std::pow(Plab, 0. ) +
1297  0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
1298  32.1 + 0. *std::pow(Plab, 0. ) +
1299  0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
1300 
1301  G4double XtotKN = ( 18.7 + 0. *std::pow(Plab, 0. ) +
1302  0.21 *sqrLogPlab - 0.89*LogPlab + //K+
1303  25.2 + 0. *std::pow(Plab, 0. ) +
1304  0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
1305 
1306  G4double XelKP = ( 5.0 + 8.1*std::pow(Plab,-1.8 )
1307  + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
1308  7.3 + 0. *std::pow(Plab,-0. ) +
1309  0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
1310 
1311  G4double XelKN = ( 7.3 + 0. *std::pow(Plab,-0. ) +
1312  0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
1313  5.0 + 8.1*std::pow(Plab,-1.8 ) +
1314  0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
1315 
1316  Xtotal = ( NumberOfTargetProtons * XtotKP +
1317  NumberOfTargetNeutrons * XtotKN );
1318 
1319  Xelastic = ( NumberOfTargetProtons * XelKP +
1320  NumberOfTargetNeutrons * XelKN );
1321  }
1322  else //------Projectile is undefined, Nucleon assumed
1323  {
1324  G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
1325  0.522*sqrLogPlab - 4.51*LogPlab;
1326 
1327  G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
1328  0.513*sqrLogPlab - 4.27*LogPlab;
1329 
1330  G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
1331  0.169*sqrLogPlab - 1.85*LogPlab;
1332  G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
1333  0.169*sqrLogPlab - 1.85*LogPlab;
1334 
1335  Xtotal = ( NumberOfTargetProtons * XtotPP +
1336  NumberOfTargetNeutrons * XtotPN );
1337 
1338  Xelastic = ( NumberOfTargetProtons * XelPP +
1339  NumberOfTargetNeutrons * XelPN );
1340  }
1341  Xinelastic = Xtotal - Xelastic;
1342 
1343  if( Xinelastic < 0.) Xinelastic = 0.;
1344 
1345  return Xinelastic*= millibarn;
1346 }
1347 
1349 //
1350 //
1351 
1352 G4double
1353 G4ComponentGGHadronNucleusXsc::GetNucleusRadius(const G4DynamicParticle* ,
1354  const G4Element* anElement)
1355 {
1356  G4int At = G4lrint(anElement->GetN());
1357  G4double oneThird = 1.0/3.0;
1358  G4double cubicrAt = std::pow(G4double(At), oneThird);
1359 
1360  G4double R; // = fRadiusConst*cubicrAt;
1361  /*
1362  G4double tmp = std::pow( cubicrAt-1., 3.);
1363  tmp += At;
1364  tmp *= 0.5;
1365 
1366  if (At > 20.) // 20.
1367  {
1368  R = fRadiusConst*std::pow (tmp, oneThird);
1369  }
1370  else
1371  {
1372  R = fRadiusConst*cubicrAt;
1373  }
1374  */
1375 
1376  R = fRadiusConst*cubicrAt;
1377 
1378  G4double meanA = 21.;
1379 
1380  G4double tauA1 = 40.;
1381  G4double tauA2 = 10.;
1382  G4double tauA3 = 5.;
1383 
1384  G4double a1 = 0.85;
1385  G4double b1 = 1. - a1;
1386 
1387  G4double b2 = 0.3;
1388  G4double b3 = 4.;
1389 
1390  if (At > 20) // 20.
1391  {
1392  R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
1393  }
1394  else if (At > 3)
1395  {
1396  R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
1397  }
1398  else
1399  {
1400  R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
1401  }
1402  return R;
1403 
1404 }
1406 //
1407 //
1408 
1409 G4double
1410 G4ComponentGGHadronNucleusXsc::GetNucleusRadius(G4int At)
1411 {
1412  G4double oneThird = 1.0/3.0;
1413  G4double cubicrAt = std::pow(G4double(At), oneThird);
1414 
1415  G4double R; // = fRadiusConst*cubicrAt;
1416 
1417  /*
1418  G4double tmp = std::pow( cubicrAt-1., 3.);
1419  tmp += At;
1420  tmp *= 0.5;
1421 
1422  if (At > 20.)
1423  {
1424  R = fRadiusConst*std::pow (tmp, oneThird);
1425  }
1426  else
1427  {
1428  R = fRadiusConst*cubicrAt;
1429  }
1430  */
1431 
1432  R = fRadiusConst*cubicrAt;
1433 
1434  G4double meanA = 20.;
1435  G4double tauA = 20.;
1436 
1437  if (At > 20) // 20.
1438  {
1439  R *= ( 0.8 + 0.2*std::exp( -(G4double(At) - meanA)/tauA) );
1440  }
1441  else
1442  {
1443  R *= ( 1.0 + 0.1*( 1. - std::exp( (G4double(At) - meanA)/tauA) ) );
1444  }
1445 
1446  return R;
1447 }
1448 
1450 //
1451 //
1452 
1453 G4double G4ComponentGGHadronNucleusXsc::CalculateEcmValue( const G4double mp ,
1454  const G4double mt ,
1455  const G4double Plab )
1456 {
1457  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1458  G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
1459  // G4double Pcm = Plab * mt / Ecm;
1460  // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
1461 
1462  return Ecm ; // KEcm;
1463 }
1464 
1466 //
1467 //
1468 
1469 G4double G4ComponentGGHadronNucleusXsc::CalcMandelstamS( const G4double mp ,
1470  const G4double mt ,
1471  const G4double Plab )
1472 {
1473  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1474  G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1475 
1476  return sMand;
1477 }
1478 
1480 //
1481 //
1482 
1483 void G4ComponentGGHadronNucleusXsc::CrossSectionDescription(std::ostream& outFile) const
1484 {
1485  outFile << "G4ComponentGGHadronNucleusXsc calculates total, inelastic and\n"
1486  << "elastic cross sections for hadron-nucleus cross sections using\n"
1487  << "the Glauber model with Gribov corrections. It is valid for all\n"
1488  << "targets except hydrogen, and for incident p, pbar, n, sigma-,\n"
1489  << "pi+, pi-, K+, K- and gammas with energies above 3 GeV. This is\n"
1490  << "a cross section component which is to be used to build a cross\n"
1491  << "data set.\n";
1492 }
1493 
1494 
1496 //
1497 // Correction arrays for GG <-> Bar changea at ~ 90 GeV
1498 
1499 const G4double G4ComponentGGHadronNucleusXsc::fNeutronBarCorrectionTot[93] = {
1500 
1501  1.0, 1.0, 1.42517e+00, // 1.118517e+00,
1502 1.082002e+00, 1.116171e+00, 1.078747e+00, 1.061315e+00,
1503 1.058205e+00, 1.082663e+00, 1.068500e+00, 1.076912e+00, 1.083475e+00, 1.079117e+00,
1504 1.071856e+00, 1.071990e+00, 1.073774e+00, 1.079356e+00, 1.081314e+00, 1.082056e+00,
1505 1.090772e+00, 1.096776e+00, 1.095828e+00, 1.097678e+00, 1.099157e+00, 1.103677e+00,
1506 1.105132e+00, 1.109806e+00, 1.110816e+00, 1.117378e+00, 1.115165e+00, 1.115710e+00,
1507 1.111855e+00, 1.110482e+00, 1.110112e+00, 1.106676e+00, 1.108706e+00, 1.105549e+00,
1508 1.106318e+00, 1.106242e+00, 1.107672e+00, 1.107342e+00, 1.108119e+00, 1.106655e+00,
1509 1.102588e+00, 1.096657e+00, 1.092920e+00, 1.086629e+00, 1.083592e+00, 1.076030e+00,
1510 1.083777e+00, 1.089460e+00, 1.086545e+00, 1.079924e+00, 1.082218e+00, 1.077798e+00,
1511 1.077062e+00, 1.072825e+00, 1.072241e+00, 1.072104e+00, 1.072490e+00, 1.069829e+00,
1512 1.070398e+00, 1.065458e+00, 1.064968e+00, 1.060524e+00, 1.060048e+00, 1.057620e+00,
1513 1.056428e+00, 1.055366e+00, 1.055017e+00, 1.052304e+00, 1.051767e+00, 1.049728e+00,
1514 1.048745e+00, 1.047399e+00, 1.045876e+00, 1.042972e+00, 1.041824e+00, 1.039993e+00,
1515 1.039021e+00, 1.036627e+00, 1.034176e+00, 1.032526e+00, 1.033633e+00, 1.036107e+00,
1516 1.037803e+00, 1.031266e+00, 1.032991e+00, 1.033284e+00, 1.035015e+00, 1.033945e+00,
1517 1.037075e+00, 1.034721e+00
1518 
1519 };
1520 
1521 const G4double G4ComponentGGHadronNucleusXsc::fNeutronBarCorrectionIn[93] = {
1522 
1523 1.0, 1.0, 1.167421e+00, 1.156250e+00, 1.205364e+00, 1.154225e+00, 1.120391e+00, // 6
1524 1.124632e+00, 1.129460e+00, 1.107863e+00, 1.102152e+00, 1.104593e+00, 1.100285e+00, // 12
1525 1.098450e+00, 1.092677e+00, 1.101124e+00, 1.106461e+00, 1.115049e+00, 1.123903e+00, // 18
1526 1.126661e+00, 1.131259e+00, 1.133949e+00, 1.134185e+00, 1.133767e+00, 1.132813e+00, // 24
1527 1.131515e+00, 1.144338e+00, // 1.130338e+00,
1528 1.134171e+00, 1.139206e+00, 1.148474e+00, // 1.141474e+00,
1529 1.142189e+00,
1530 1.140725e+00, 1.140100e+00, 1.139848e+00, 1.137674e+00, 1.138645e+00, 1.136339e+00,
1531 1.136439e+00, 1.135946e+00, 1.136431e+00, 1.135702e+00, 1.135703e+00, 1.134113e+00,
1532 1.131935e+00, 1.128381e+00, 1.126373e+00, 1.122453e+00, 1.120908e+00, 1.115953e+00,
1533 1.115947e+00, 1.114426e+00, 1.111749e+00, 1.106207e+00, 1.107494e+00, 1.103622e+00,
1534 1.102576e+00, 1.098816e+00, 1.097889e+00, 1.097306e+00, 1.097130e+00, 1.094578e+00,
1535 1.094552e+00, 1.090222e+00, 1.089358e+00, 1.085409e+00, 1.084560e+00, 1.082182e+00,
1536 1.080773e+00, 1.079464e+00, 1.078724e+00, 1.076121e+00, 1.075235e+00, 1.073159e+00,
1537 1.071920e+00, 1.070395e+00, 1.069503e+00, 1.067525e+00, 1.066919e+00, 1.065779e+00,
1538 1.065319e+00, 1.063730e+00, 1.062092e+00, 1.061085e+00, 1.059908e+00, 1.059815e+00,
1539 1.059109e+00, 1.051920e+00, 1.051258e+00, 1.049473e+00, 1.048823e+00, 1.045984e+00,
1540 1.046435e+00, 1.042614e+00
1541 
1542 };
1543 
1544 const G4double G4ComponentGGHadronNucleusXsc::fProtonBarCorrectionTot[93] = {
1545 
1546 1.0, 1.0,
1547 1.118515e+00, 1.082000e+00, 1.116169e+00, 1.078745e+00, 1.061313e+00, 1.058203e+00,
1548 1.082661e+00, 1.068498e+00, 1.076910e+00, 1.083474e+00, 1.079115e+00, 1.071854e+00,
1549 1.071988e+00, 1.073772e+00, 1.079355e+00, 1.081312e+00, 1.082054e+00, 1.090770e+00,
1550 1.096774e+00, 1.095827e+00, 1.097677e+00, 1.099156e+00, 1.103676e+00, 1.105130e+00,
1551 1.109805e+00, 1.110814e+00, 1.117377e+00, 1.115163e+00, 1.115708e+00, 1.111853e+00,
1552 1.110480e+00, 1.110111e+00, 1.106674e+00, 1.108705e+00, 1.105548e+00, 1.106317e+00,
1553 1.106241e+00, 1.107671e+00, 1.107341e+00, 1.108118e+00, 1.106654e+00, 1.102586e+00,
1554 1.096655e+00, 1.092918e+00, 1.086628e+00, 1.083590e+00, 1.076028e+00, 1.083776e+00,
1555 1.089458e+00, 1.086543e+00, 1.079923e+00, 1.082216e+00, 1.077797e+00, 1.077061e+00,
1556 1.072824e+00, 1.072239e+00, 1.072103e+00, 1.072488e+00, 1.069828e+00, 1.070396e+00,
1557 1.065456e+00, 1.064966e+00, 1.060523e+00, 1.060047e+00, 1.057618e+00, 1.056427e+00,
1558 1.055365e+00, 1.055016e+00, 1.052303e+00, 1.051766e+00, 1.049727e+00, 1.048743e+00,
1559 1.047397e+00, 1.045875e+00, 1.042971e+00, 1.041823e+00, 1.039992e+00, 1.039019e+00,
1560 1.036626e+00, 1.034175e+00, 1.032525e+00, 1.033632e+00, 1.036106e+00, 1.037802e+00,
1561 1.031265e+00, 1.032990e+00, 1.033283e+00, 1.035014e+00, 1.033944e+00, 1.037074e+00,
1562 1.034720e+00
1563 
1564 };
1565 
1566 const G4double G4ComponentGGHadronNucleusXsc::fProtonBarCorrectionIn[93] = {
1567 
1568 1.0, 1.0,
1569 1.147419e+00, // 1.167419e+00,
1570 1.156248e+00, 1.205362e+00, 1.154224e+00, 1.120390e+00, 1.124630e+00, // 7
1571 1.129459e+00, 1.107861e+00, 1.102151e+00, 1.104591e+00, 1.100284e+00, 1.098449e+00, // 13
1572 1.092675e+00, 1.101122e+00, 1.106460e+00, 1.115048e+00, 1.123902e+00, 1.126659e+00, // 19
1573 1.131258e+00, 1.133948e+00, 1.134183e+00, 1.133766e+00, 1.132812e+00, 1.131514e+00, // 25
1574 1.140337e+00, // 1.130337e+00,
1575 
1576 1.134170e+00, 1.139205e+00, 1.151472e+00, // 1.141472e+00,
1577 1.142188e+00, 1.140724e+00,
1578 1.140099e+00, 1.139847e+00, 1.137672e+00, 1.138644e+00, 1.136338e+00, 1.136438e+00,
1579 1.135945e+00, 1.136429e+00, 1.135701e+00, 1.135702e+00, 1.134112e+00, 1.131934e+00,
1580 1.128380e+00, 1.126371e+00, 1.122452e+00, 1.120907e+00, 1.115952e+00, 1.115946e+00,
1581 1.114425e+00, 1.111748e+00, 1.106205e+00, 1.107493e+00, 1.103621e+00, 1.102575e+00,
1582 1.098815e+00, 1.097888e+00, 1.097305e+00, 1.097129e+00, 1.094577e+00, 1.094551e+00,
1583 1.090221e+00, 1.089357e+00, 1.085408e+00, 1.084559e+00, 1.082181e+00, 1.080772e+00,
1584 1.079463e+00, 1.078723e+00, 1.076120e+00, 1.075234e+00, 1.073158e+00, 1.071919e+00,
1585 1.070394e+00, 1.069502e+00, 1.067524e+00, 1.066918e+00, 1.065778e+00, 1.065318e+00,
1586 1.063729e+00, 1.062091e+00, 1.061084e+00, 1.059907e+00, 1.059814e+00, 1.059108e+00,
1587 1.051919e+00, 1.051257e+00, 1.049472e+00, 1.048822e+00, 1.045983e+00, 1.046434e+00,
1588 1.042613e+00
1589 
1590 };
1591 
1592 
1593 const G4double G4ComponentGGHadronNucleusXsc::fPionPlusBarCorrectionTot[93] = {
1594 
1595 1.0, 1.0,
1596 1.075927e+00, 1.074407e+00, 1.126098e+00, 1.100127e+00, 1.089742e+00, 1.083536e+00,
1597 1.089988e+00, 1.103566e+00, 1.096922e+00, 1.126573e+00, 1.132734e+00, 1.136512e+00,
1598 1.136629e+00, 1.133086e+00, 1.132428e+00, 1.129299e+00, 1.125622e+00, 1.126992e+00,
1599 1.127840e+00, 1.162670e+00, 1.160392e+00, 1.157864e+00, 1.157227e+00, 1.154627e+00,
1600 1.192555e+00, 1.197243e+00, 1.197911e+00, 1.200326e+00, 1.220053e+00, 1.215019e+00,
1601 1.211703e+00, 1.209080e+00, 1.204248e+00, 1.203328e+00, 1.198671e+00, 1.196840e+00,
1602 1.194392e+00, 1.193037e+00, 1.190408e+00, 1.188583e+00, 1.206127e+00, 1.210028e+00,
1603 1.206434e+00, 1.204456e+00, 1.200547e+00, 1.199058e+00, 1.200174e+00, 1.200276e+00,
1604 1.198912e+00, 1.213048e+00, 1.207160e+00, 1.208020e+00, 1.203814e+00, 1.202380e+00,
1605 1.198306e+00, 1.197002e+00, 1.196027e+00, 1.195449e+00, 1.192563e+00, 1.192135e+00,
1606 1.187556e+00, 1.186308e+00, 1.182124e+00, 1.180900e+00, 1.178224e+00, 1.176471e+00,
1607 1.174811e+00, 1.173702e+00, 1.170827e+00, 1.169581e+00, 1.167205e+00, 1.165626e+00,
1608 1.180244e+00, 1.177626e+00, 1.175121e+00, 1.173903e+00, 1.172192e+00, 1.171128e+00,
1609 1.168997e+00, 1.166826e+00, 1.164130e+00, 1.165412e+00, 1.165504e+00, 1.165020e+00,
1610 1.158462e+00, 1.158014e+00, 1.156519e+00, 1.156081e+00, 1.153602e+00, 1.154190e+00,
1611 1.152974e+00
1612 
1613 };
1614 
1615 const G4double G4ComponentGGHadronNucleusXsc::fPionPlusBarCorrectionIn[93] = {
1616 
1617 1.0, 1.0,
1618 1.140246e+00, 1.097872e+00, 1.104301e+00, 1.068722e+00, 1.056495e+00, 1.062622e+00, // 7
1619 1.047987e+00, 1.037032e+00, 1.035686e+00, 1.042870e+00, 1.052222e+00, 1.075100e+00, // 13
1620 1.084480e+00, 1.078286e+00, 1.081488e+00, 1.089713e+00, 1.099105e+00, 1.098003e+00, // 19
1621 1.102175e+00, 1.117707e+00, 1.121734e+00, 1.125229e+00, 1.126457e+00, 1.128905e+00, // 25
1622 1.163312e+00, 1.126263e+00, 1.126459e+00, 1.135191e+00, 1.116986e+00, 1.117184e+00, // 31
1623 1.117037e+00, 1.116777e+00, 1.115858e+00, 1.115745e+00, 1.114489e+00, 1.113993e+00, // 37
1624 1.113226e+00, 1.112818e+00, 1.111890e+00, 1.111238e+00, 1.111209e+00, 1.111775e+00, // 43
1625 1.110256e+00, 1.109414e+00, 1.107647e+00, 1.106980e+00, 1.106096e+00, 1.107331e+00, // 49
1626 1.107849e+00, 1.106407e+00, 1.103426e+00, 1.103896e+00, 1.101756e+00, 1.101031e+00, // 55
1627 1.098915e+00, 1.098260e+00, 1.097768e+00, 1.097487e+00, 1.095964e+00, 1.095773e+00, // 61
1628 1.093348e+00, 1.092687e+00, 1.090465e+00, 1.089821e+00, 1.088394e+00, 1.087462e+00, // 67
1629 1.086571e+00, 1.085997e+00, 1.084451e+00, 1.083798e+00, 1.082513e+00, 1.081670e+00, // 73
1630 1.080735e+00, 1.075659e+00, 1.074341e+00, 1.073689e+00, 1.072787e+00, 1.072237e+00, // 79
1631 1.071107e+00, 1.069955e+00, 1.074856e+00, 1.065873e+00, 1.065938e+00, 1.065694e+00,
1632 1.062192e+00, 1.061967e+00, 1.061180e+00, 1.060960e+00, 1.059646e+00, 1.059975e+00,
1633 1.059658e+00
1634 
1635 };
1636 
1637 
1638 const G4double G4ComponentGGHadronNucleusXsc::fPionMinusBarCorrectionTot[93] = {
1639 
1640 1.0, 1.0,
1641 1.3956e+00, 1.077959e+00, 1.129145e+00, 1.102088e+00, 1.089765e+00, 1.083542e+00, // 7
1642 1.089995e+00, 1.104895e+00, 1.097154e+00, 1.127663e+00, 1.133063e+00, 1.137425e+00, // 13
1643 1.136724e+00, 1.133859e+00, 1.132498e+00, 1.130276e+00, 1.127896e+00, 1.127656e+00, // 19
1644 1.127905e+00, 1.164210e+00, 1.162259e+00, 1.160075e+00, 1.158978e+00, 1.156649e+00, // 25
1645 1.194157e+00, 1.199177e+00, 1.198983e+00, 1.202325e+00, 1.221967e+00, 1.217548e+00,
1646 1.214389e+00, 1.211760e+00, 1.207335e+00, 1.206081e+00, 1.201766e+00, 1.199779e+00,
1647 1.197283e+00, 1.195706e+00, 1.193071e+00, 1.191115e+00, 1.208838e+00, 1.212681e+00,
1648 1.209235e+00, 1.207163e+00, 1.203451e+00, 1.201807e+00, 1.203283e+00, 1.203388e+00,
1649 1.202244e+00, 1.216509e+00, 1.211066e+00, 1.211504e+00, 1.207539e+00, 1.205991e+00,
1650 1.202143e+00, 1.200724e+00, 1.199595e+00, 1.198815e+00, 1.196025e+00, 1.195390e+00,
1651 1.191137e+00, 1.189791e+00, 1.185888e+00, 1.184575e+00, 1.181996e+00, 1.180229e+00,
1652 1.178545e+00, 1.177355e+00, 1.174616e+00, 1.173312e+00, 1.171016e+00, 1.169424e+00,
1653 1.184120e+00, 1.181478e+00, 1.179085e+00, 1.177817e+00, 1.176124e+00, 1.175003e+00,
1654 1.172947e+00, 1.170858e+00, 1.168170e+00, 1.169397e+00, 1.169304e+00, 1.168706e+00,
1655 1.162774e+00, 1.162217e+00, 1.160740e+00, 1.160196e+00, 1.157857e+00, 1.158220e+00,
1656 1.157267e+00
1657 };
1658 
1659 
1660 const G4double G4ComponentGGHadronNucleusXsc::fPionMinusBarCorrectionIn[93] = {
1661 
1662 1.0, 1.0,
1663 1.463e+00, 1.100898e+00, 1.106773e+00, 1.070289e+00, 1.040514e+00, 1.062628e+00, // 7
1664 1.047992e+00, 1.038041e+00, 1.035862e+00, 1.043679e+00, 1.052466e+00, 1.065780e+00, // 13
1665 1.070551e+00, 1.078869e+00, 1.081541e+00, 1.090455e+00, 1.100847e+00, 1.098511e+00, // 19
1666 1.102226e+00, 1.118865e+00, 1.123143e+00, 1.126904e+00, 1.127785e+00, 1.130444e+00, // 25
1667 1.148502e+00, 1.127678e+00, 1.127244e+00, 1.123634e+00, 1.118347e+00, 1.118988e+00,
1668 1.118957e+00, 1.118696e+00, 1.118074e+00, 1.117722e+00, 1.116717e+00, 1.116111e+00,
1669 1.115311e+00, 1.114745e+00, 1.113814e+00, 1.113069e+00, 1.113141e+00, 1.113660e+00,
1670 1.112249e+00, 1.111343e+00, 1.109718e+00, 1.108942e+00, 1.108310e+00, 1.109549e+00,
1671 1.110227e+00, 1.108846e+00, 1.106183e+00, 1.106354e+00, 1.104388e+00, 1.103583e+00,
1672 1.101632e+00, 1.100896e+00, 1.100296e+00, 1.099873e+00, 1.098420e+00, 1.098082e+00,
1673 1.095892e+00, 1.095162e+00, 1.093144e+00, 1.092438e+00, 1.091083e+00, 1.090142e+00,
1674 1.089236e+00, 1.088604e+00, 1.087159e+00, 1.086465e+00, 1.085239e+00, 1.084388e+00,
1675 1.083473e+00, 1.078373e+00, 1.077136e+00, 1.076450e+00, 1.075561e+00, 1.074973e+00,
1676 1.073898e+00, 1.072806e+00, 1.067706e+00, 1.068684e+00, 1.068618e+00, 1.068294e+00,
1677 1.065241e+00, 1.064939e+00, 1.064166e+00, 1.063872e+00, 1.062659e+00, 1.062828e+00,
1678 1.062699e+00
1679 
1680 };
1681 
1682 
1683 //
1684 //
G4ComponentGGHadronNucleusXsc::GetRatioQE
G4double GetRatioQE(const G4DynamicParticle *, G4int At, G4int Zt)
Definition: G4ComponentGGHadronNucleusXsc.cc:437
G4ComponentGGHadronNucleusXsc::GetNucleusRadius
G4double GetNucleusRadius(const G4DynamicParticle *, const G4Element *)
Definition: G4ComponentGGHadronNucleusXsc.cc:1353
G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG
G4double GetHadronNucleonXscPDG(const G4DynamicParticle *, const G4Element *)
Definition: G4ComponentGGHadronNucleusXsc.cc:571
G4HadronNucleonXsc::GetKmNeutronTotXscVector
G4double GetKmNeutronTotXscVector(G4double logEnergy)
Definition: G4HadronNucleonXsc.hh:99
G4ComponentGGHadronNucleusXsc::GetProductionElementCrossSection
virtual G4double GetProductionElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentGGHadronNucleusXsc.cc:168
G4HadronNucleonXsc::GetKpProtonTotXscVector
G4double GetKpProtonTotXscVector(G4double logEnergy)
Definition: G4HadronNucleonXsc.hh:96
G4AntiOmegaMinus::AntiOmegaMinus
static G4AntiOmegaMinus * AntiOmegaMinus()
Definition: G4AntiOmegaMinus.cc:110
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetTotalEnergy
G4double GetTotalEnergy() const
outFile
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4Element::GetN
G4double GetN() const
Definition: G4Element.hh:134
G4ComponentGGHadronNucleusXsc::CalcMandelstamS
G4double CalcMandelstamS(const G4double, const G4double, const G4double)
Definition: G4ComponentGGHadronNucleusXsc.cc:1469
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
G4DynamicParticle.hh
G4ComponentGGHadronNucleusXsc::GetRatioSD
G4double GetRatioSD(const G4DynamicParticle *, G4int At, G4int Zt)
Definition: G4ComponentGGHadronNucleusXsc.cc:395
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4AntiSigmaPlus::AntiSigmaPlus
static G4AntiSigmaPlus * AntiSigmaPlus()
Definition: G4AntiSigmaPlus.cc:108
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4ComponentGGHadronNucleusXsc::GetIsoCrossSection
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
Definition: G4ComponentGGHadronNucleusXsc.cc:270
G4IonTable.hh
G4SigmaZero::SigmaZero
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99
G4ComponentGGHadronNucleusXsc::GetHNinelasticXscVU
G4double GetHNinelasticXscVU(const G4DynamicParticle *, G4int At, G4int Zt)
Definition: G4ComponentGGHadronNucleusXsc.cc:1141
G4KaonMinus::KaonMinus
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4HadronNucleonXsc::GetHadronNucleonXscNS
G4double GetHadronNucleonXscNS(const G4DynamicParticle *, const G4ParticleDefinition *)
Definition: G4HadronNucleonXsc.cc:365
G4ComponentGGHadronNucleusXsc::GetInelasticIsotopeCrossSection
virtual G4double GetInelasticIsotopeCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
Definition: G4ComponentGGHadronNucleusXsc.cc:126
G4AntiSigmaMinus::AntiSigmaMinus
static G4AntiSigmaMinus * AntiSigmaMinus()
Definition: G4AntiSigmaMinus.cc:106
G4XiZero::XiZero
static G4XiZero * XiZero()
Definition: G4XiZero.cc:106
G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS
G4double GetHadronNucleonXscNS(const G4DynamicParticle *, const G4Element *)
Definition: G4ComponentGGHadronNucleusXsc.cc:705
G4ParticleTable.hh
G4HadronNucleonXsc::GetKaonNucleonXscGG
G4double GetKaonNucleonXscGG(const G4DynamicParticle *, const G4ParticleDefinition *)
Definition: G4HadronNucleonXsc.cc:1069
Z
Float_t Z
Definition: plot.C:39
G4KaonZeroShort::KaonZeroShort
static G4KaonZeroShort * KaonZeroShort()
Definition: G4KaonZeroShort.cc:104
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4DynamicParticle::GetMass
G4double GetMass() const
G4bool
bool G4bool
Definition: G4Types.hh:79
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
G4ParticleDefinition.hh
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4ComponentGGHadronNucleusXsc::GetInelasticElementCrossSection
virtual G4double GetInelasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentGGHadronNucleusXsc.cc:154
G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc
G4double GetHadronNucleonXsc(const G4DynamicParticle *, const G4Element *)
Definition: G4ComponentGGHadronNucleusXsc.cc:484
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc
G4double GetHNinelasticXsc(const G4DynamicParticle *, const G4Element *)
Definition: G4ComponentGGHadronNucleusXsc.cc:1099
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4PhysicalConstants.hh
G4ComponentGGHadronNucleusXsc::GetParticleBarCorIn
G4double GetParticleBarCorIn(const G4ParticleDefinition *theParticle, G4int Z)
Definition: G4ComponentGGHadronNucleusXsc.hh:264
G4SigmaMinus::SigmaMinus
static G4SigmaMinus * SigmaMinus()
Definition: G4SigmaMinus.cc:106
zz
Double_t zz
Definition: macro.C:12
G4ComponentGGHadronNucleusXsc::CrossSectionDescription
virtual void CrossSectionDescription(std::ostream &) const
Definition: G4ComponentGGHadronNucleusXsc.cc:1483
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:95
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:99
G4HadronNucleonXsc::GetKmProtonTotXscVector
G4double GetKmProtonTotXscVector(G4double logEnergy)
Definition: G4HadronNucleonXsc.hh:98
G4ComponentGGHadronNucleusXsc::GetTotalIsotopeCrossSection
virtual G4double GetTotalIsotopeCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
Definition: G4ComponentGGHadronNucleusXsc.cc:98
G4ComponentGGHadronNucleusXsc::GetElasticIsotopeCrossSection
virtual G4double GetElasticIsotopeCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
Definition: G4ComponentGGHadronNucleusXsc.cc:196
G4HadronNucleonXsc.hh
G4AntiXiZero::AntiXiZero
static G4AntiXiZero * AntiXiZero()
Definition: G4AntiXiZero.cc:106
N
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
G4ComponentGGHadronNucleusXsc.hh
G4ComponentGGHadronNucleusXsc::CalculateEcmValue
G4double CalculateEcmValue(const G4double, const G4double, const G4double)
Definition: G4ComponentGGHadronNucleusXsc.cc:1453
G4ComponentGGHadronNucleusXsc::IsIsoApplicable
G4bool IsIsoApplicable(const G4DynamicParticle *aDP, G4int Z, G4int A, const G4Element *elm=0, const G4Material *mat=0)
Definition: G4ComponentGGHadronNucleusXsc.cc:234
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
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
G4ComponentGGHadronNucleusXsc::GetTotalElementCrossSection
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentGGHadronNucleusXsc.cc:112
G4SystemOfUnits.hh
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4ComponentGGHadronNucleusXsc::GetElasticElementCrossSection
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentGGHadronNucleusXsc.cc:182
G4ComponentGGHadronNucleusXsc::GetKaonNucleonXscVector
G4double GetKaonNucleonXscVector(const G4DynamicParticle *, G4int At, G4int Zt)
Definition: G4ComponentGGHadronNucleusXsc.cc:1061
G4ParticleMomentum
G4ThreeVector G4ParticleMomentum
Definition: G4ParticleMomentum.hh:50
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
CLHEP::Hep3Vector::mag
double mag() const
G4INCL::Math::oneThird
const G4double oneThird
Definition: G4INCLGlobals.hh:78
G4ComponentGGHadronNucleusXsc::ComputeQuasiElasticRatio
virtual G4double ComputeQuasiElasticRatio(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
Definition: G4ComponentGGHadronNucleusXsc.cc:210
G4ComponentGGHadronNucleusXsc::GetParticleBarCorTot
G4double GetParticleBarCorTot(const G4ParticleDefinition *theParticle, G4int Z)
Definition: G4ComponentGGHadronNucleusXsc.hh:244
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
G4HadronNucleonXsc::GetInelasticHadronNucleonXsc
G4double GetInelasticHadronNucleonXsc()
Definition: G4HadronNucleonXsc.hh:92
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4HadronNucleonXsc::GetKpNeutronTotXscVector
G4double GetKpNeutronTotXscVector(G4double logEnergy)
Definition: G4HadronNucleonXsc.hh:97
G4ComponentGGHadronNucleusXsc::GetProductionIsotopeCrossSection
virtual G4double GetProductionIsotopeCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
Definition: G4ComponentGGHadronNucleusXsc.cc:140
Generated on Sat Dec 14 2013 14:34:41 for Geant4_10 by   doxygen 1.8.5