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