Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4AntiNuclElastic Class Reference

#include <G4AntiNuclElastic.hh>

+ Inheritance diagram for G4AntiNuclElastic:

Public Member Functions

 G4AntiNuclElastic ()
 
virtual ~G4AntiNuclElastic ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double SampleThetaCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double SampleThetaLab (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double CalculateParticleBeta (const G4ParticleDefinition *particle, G4double momentum)
 
G4double CalculateZommerfeld (G4double beta, G4double Z1, G4double Z2)
 
G4double CalculateAm (G4double momentum, G4double n, G4double Z)
 
G4double DampFactor (G4double z)
 
G4double BesselJzero (G4double z)
 
G4double BesselJone (G4double z)
 
G4double BesselOneByArg (G4double z)
 
G4double GetcosTeta1 (G4double plab, G4int A)
 
G4ComponentAntiNuclNuclearXS * GetComponentCrossSection ()
 
- Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 
virtual ~G4HadronElastic ()
 
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
void SetLowestEnergyLimit (G4double value)
 
G4double LowestEnergyLimit () const
 
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual void Description () const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
const G4HadronicInteraction * GetMyPointer () const
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4String & GetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
G4bool operator== (const G4HadronicInteraction &right) const
 
G4bool operator!= (const G4HadronicInteraction &right) const
 
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 50 of file G4AntiNuclElastic.hh.

Constructor & Destructor Documentation

◆ G4AntiNuclElastic()

G4AntiNuclElastic::G4AntiNuclElastic ( )

Definition at line 55 of file G4AntiNuclElastic.cc.

56 : G4HadronElastic("AntiAElastic")
57{
58 //V.Ivanchenko commented out
59 //SetMinEnergy( 0.1*GeV );
60 //SetMaxEnergy( 10.*TeV );
61
62
63 theAProton = G4AntiProton::AntiProton();
64 theANeutron = G4AntiNeutron::AntiNeutron();
65 theADeuteron = G4AntiDeuteron::AntiDeuteron();
66 theATriton = G4AntiTriton::AntiTriton();
67 theAAlpha = G4AntiAlpha::AntiAlpha();
68 theAHe3 = G4AntiHe3::AntiHe3();
69
70 theProton = G4Proton::Proton();
71 theNeutron = G4Neutron::Neutron();
72 theDeuteron = G4Deuteron::Deuteron();
73 theAlpha = G4Alpha::Alpha();
74
75
76 cs = new G4ComponentAntiNuclNuclearXS();
77 fParticle = 0;
78 fWaveVector = 0.;
79 fBeta = 0.;
80 fZommerfeld = 0.;
81 fAm = 0.;
82 fTetaCMS = 0.;
83 fRa = 0.;
84 fRef = 0.;
85 fceff = 0.;
86 fptot = 0.;
87 fTmax = 0.;
88 fThetaLab = 0.;
89}
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4AntiAlpha::AntiAlpha
static G4AntiAlpha * AntiAlpha()
Definition: G4AntiAlpha.cc:89
G4AntiDeuteron::AntiDeuteron
static G4AntiDeuteron * AntiDeuteron()
Definition: G4AntiDeuteron.cc:94
G4AntiHe3::AntiHe3
static G4AntiHe3 * AntiHe3()
Definition: G4AntiHe3.cc:94
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4AntiTriton::AntiTriton
static G4AntiTriton * AntiTriton()
Definition: G4AntiTriton.cc:94
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93

◆ ~G4AntiNuclElastic()

G4AntiNuclElastic::~G4AntiNuclElastic ( )
virtual

Definition at line 92 of file G4AntiNuclElastic.cc.

93{
94 delete cs;
95}

Member Function Documentation

◆ BesselJone()

G4double G4AntiNuclElastic::BesselJone ( G4double  z)

Definition at line 583 of file G4AntiNuclElastic.cc.

584{
585 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
586
587 modvalue = std::fabs(value);
588
589 if ( modvalue < 8.0 )
590 {
591 value2 = value*value;
592 fact1 = value*(72362614232.0 + value2*(-7895059235.0
593 + value2*( 242396853.1
594 + value2*(-2972611.439
595 + value2*( 15704.48260
596 + value2*(-30.16036606 ) ) ) ) ) );
597
598 fact2 = 144725228442.0 + value2*(2300535178.0
599 + value2*(18583304.74
600 + value2*(99447.43394
601 + value2*(376.9991397
602 + value2*1.0 ) ) ) );
603 bessel = fact1/fact2;
604 }
605 else
606 {
607 arg = 8.0/modvalue;
608 value2 = arg*arg;
609
610 shift = modvalue - 2.356194491;
611
612 fact1 = 1.0 + value2*( 0.183105e-2
613 + value2*(-0.3516396496e-4
614 + value2*(0.2457520174e-5
615 + value2*(-0.240337019e-6 ) ) ) );
616
617 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
618 + value2*( 0.8449199096e-5
619 + value2*(-0.88228987e-6
620 + value2*0.105787412e-6 ) ) );
621
622 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
623 if (value < 0.0) bessel = -bessel;
624 }
625 return bessel;
626}
G4double
double G4double
Definition: G4Types.hh:64

Referenced by BesselOneByArg().

◆ BesselJzero()

G4double G4AntiNuclElastic::BesselJzero ( G4double  z)

Definition at line 532 of file G4AntiNuclElastic.cc.

533{
534 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
535
536 modvalue = std::fabs(value);
537
538 if ( value < 8.0 && value > -8.0 )
539 {
540 value2 = value*value;
541
542 fact1 = 57568490574.0 + value2*(-13362590354.0
543 + value2*( 651619640.7
544 + value2*(-11214424.18
545 + value2*( 77392.33017
546 + value2*(-184.9052456 ) ) ) ) );
547
548 fact2 = 57568490411.0 + value2*( 1029532985.0
549 + value2*( 9494680.718
550 + value2*(59272.64853
551 + value2*(267.8532712
552 + value2*1.0 ) ) ) );
553
554 bessel = fact1/fact2;
555 }
556 else
557 {
558 arg = 8.0/modvalue;
559
560 value2 = arg*arg;
561
562 shift = modvalue-0.785398164;
563
564 fact1 = 1.0 + value2*(-0.1098628627e-2
565 + value2*(0.2734510407e-4
566 + value2*(-0.2073370639e-5
567 + value2*0.2093887211e-6 ) ) );
568 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
569 + value2*(-0.6911147651e-5
570 + value2*(0.7621095161e-6
571 - value2*0.934945152e-7 ) ) );
572
573 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
574 }
575 return bessel;
576}

Referenced by SampleInvariantT().

◆ BesselOneByArg()

G4double G4AntiNuclElastic::BesselOneByArg ( G4double  z)

Definition at line 630 of file G4AntiNuclElastic.cc.

631{
632 G4double x2, result;
633
634 if( std::fabs(x) < 0.01 )
635 {
636 x *= 0.5;
637 x2 = x*x;
638 result = (2.- x2 + x2*x2/6.)/4.;
639 }
640 else
641 {
642 result = BesselJone(x)/x;
643 }
644 return result;
645}
G4AntiNuclElastic::BesselJone
G4double BesselJone(G4double z)
Definition: G4AntiNuclElastic.cc:583

Referenced by SampleInvariantT().

◆ CalculateAm()

G4double G4AntiNuclElastic::CalculateAm ( G4double  momentum,
G4double  n,
G4double  Z 
)

Definition at line 516 of file G4AntiNuclElastic.cc.

517{
518 G4double k = momentum/hbarc;
519 G4double ch = 1.13 + 3.76*n*n;
520 G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius;
521 G4double zn2 = zn*zn;
522 fAm = ch/zn2;
523
524 return fAm;
525}
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:50
G4Pow::A13
G4double A13(G4double A)
Definition: G4Pow.hh:115

Referenced by SampleInvariantT().

◆ CalculateParticleBeta()

G4double G4AntiNuclElastic::CalculateParticleBeta ( const G4ParticleDefinition *  particle,
G4double  momentum 
)

Definition at line 493 of file G4AntiNuclElastic.cc.

495{
496 G4double mass = particle->GetPDGMass();
497 G4double a = momentum/mass;
498 fBeta = a/std::sqrt(1+a*a);
499
500 return fBeta;
501}

Referenced by SampleInvariantT().

◆ CalculateZommerfeld()

G4double G4AntiNuclElastic::CalculateZommerfeld ( G4double  beta,
G4double  Z1,
G4double  Z2 
)

Definition at line 507 of file G4AntiNuclElastic.cc.

508{
509 fZommerfeld = fine_structure_const*Z1*Z2/beta;
510
511 return fZommerfeld;
512}

Referenced by SampleInvariantT().

◆ DampFactor()

G4double G4AntiNuclElastic::DampFactor ( G4double  z)

Definition at line 473 of file G4AntiNuclElastic.cc.

474{
475 G4double df;
476 G4double f3 = 6.; // first factorials
477
478 if( std::fabs(x) < 0.01 )
479 {
480 df=1./(1.+x*x/f3);
481 }
482 else
483 {
484 df = x/std::sinh(x);
485 }
486 return df;
487}

Referenced by SampleInvariantT().

◆ GetComponentCrossSection()

G4ComponentAntiNuclNuclearXS * G4AntiNuclElastic::GetComponentCrossSection ( )
inline

Definition at line 126 of file G4AntiNuclElastic.hh.

127{
128 return cs;
129}

Referenced by G4HadronDElasticPhysics::ConstructProcess(), G4HadronElasticPhysics::ConstructProcess(), G4HadronHElasticPhysics::ConstructProcess(), and G4HadronQElasticPhysics::ConstructProcess().

◆ GetcosTeta1()

G4double G4AntiNuclElastic::GetcosTeta1 ( G4double  plab,
G4int  A 
)

Definition at line 649 of file G4AntiNuclElastic.cc.

650{
651
652// G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
653 G4double p0 = 1.*hbarc/fermi;
654//G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab);
655 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab);
656//////////////////
657 if(cteta1 < -1.) cteta1 = -1.0;
658 return cteta1;
659}
G4Pow::Z23
G4double Z23(G4int Z)
Definition: G4Pow.hh:134

Referenced by SampleInvariantT().

◆ SampleInvariantT()

G4double G4AntiNuclElastic::SampleInvariantT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 99 of file G4AntiNuclElastic.cc.

101{
102 G4double T;
103 G4double Mproj = particle->GetPDGMass();
104 G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(Plab*Plab+Mproj*Mproj));
105 G4double ctet1 = GetcosTeta1(Plab, A);
106
107 G4double energy=Pproj.e()-Mproj;
108
109 const G4ParticleDefinition* theParticle = particle;
110
111 G4ParticleDefinition * theDef = 0;
112
113 if(Z == 1 && A == 1) theDef = theProton;
114 else if (Z == 1 && A == 2) theDef = theDeuteron;
115 else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
116 else if (Z == 2 && A == 3) theDef = G4He3::He3();
117 else if (Z == 2 && A == 4) theDef = theAlpha;
118
119
120 G4double TargMass =G4NucleiProperties::GetNuclearMass(A,Z);
121
122 //transform to CMS
123
124 G4LorentzVector lv(0.0,0.0,0.0,TargMass);
125 lv += Pproj;
126 G4double S = lv.mag2()/GeV/GeV;
127
128 G4ThreeVector bst = lv.boostVector();
129 Pproj.boost(-bst);
130
131 G4ThreeVector p1 = Pproj.vect();
132 G4double ptot = p1.mag();
133
134 fbst = bst;
135 fptot= ptot;
136 fTmax = 4.0*ptot*ptot;
137
138 if(Plab/std::abs(particle->GetBaryonNumber()) < 100.*MeV) // Uzhi 24 Nov. 2011
139 {return fTmax*G4UniformRand();} // Uzhi 24 Nov. 2011
140
141 G4double Z1 = particle->GetPDGCharge();
142 G4double Z2 = Z;
143
144 G4double beta = CalculateParticleBeta(particle, ptot);
145 G4double n = CalculateZommerfeld( beta, Z1, Z2 );
146 G4double Am = CalculateAm( ptot, n, Z2 );
147 fWaveVector = ptot; // /hbarc;
148
149 G4LorentzVector Fproj(0.,0.,0.,0.);
150 G4double XsCoulomb = sqr(n/fWaveVector)*pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1);
151 XsCoulomb=XsCoulomb*0.38938e+6;
152
153
154 G4double XsElastHad =cs->GetElasticElementCrossSection(particle, energy, Z, (G4double)A);
155 G4double XstotalHad =cs->GetTotalElementCrossSection(particle, energy, Z, (G4double)A);
156
157
158 XsElastHad/=millibarn; XstotalHad/=millibarn;
159
160
161 G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHad);
162
163// G4cout<<" XselastHadron " << XsElastHad << " XsCol "<< XsCoulomb <<G4endl;
164// G4cout <<" XsTotal" << XstotalHad <<G4endl;
165// G4cout<<"XsInel"<< XstotalHad-XsElastHad<<G4endl;
166
167 if(G4UniformRand() < CoulombProb)
168 { // Simulation of Coulomb scattering
169
170 G4double phi = twopi * G4UniformRand();
171 G4double Ksi = G4UniformRand();
172
173 G4double par1 = 2.*(1.+Am)/(1.+ctet1);
174
175// ////sample ThetaCMS in Coulomb part
176
177 G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi);
178
179 G4double PtZ=ptot*cosThetaCMS;
180 Fproj.setPz(PtZ);
181 G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS);
182 G4double PtX= PtProjCMS * std::cos(phi);
183 G4double PtY= PtProjCMS * std::sin(phi);
184 Fproj.setPx(PtX);
185 Fproj.setPy(PtY);
186 Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj));
187 T = -(Pproj-Fproj).mag2();
188 } else
189
190 {
191///////Simulation of strong interaction scattering////////////////////////////
192
193// G4double Qmax = 2.*ptot*197.33; // in fm^-1
194 G4double Qmax = 2.*3.0*197.33; // in fm^-1
195 G4double Amag = 70*70; // A1 in Magora funct:A1*exp(-q*A2)
196 G4double SlopeMag = 2.*3.0; // A2 in Magora funct:A1*exp(-q*A2)
197
198 G4double sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy);
199
200 fRa = 1.113*G4Pow::GetInstance()->Z13(A) -
201 0.227/G4Pow::GetInstance()->Z13(A);
202 if(A == 3) fRa=1.81;
203 if(A == 4) fRa=1.37;
204
205 if((A>=12.) && (A<27) ) fRa=fRa*0.85;
206 if((A>=27.) && (A<48) ) fRa=fRa*0.90;
207 if((A>=48.) && (A<65) ) fRa=fRa*0.95;
208
209 G4double Ref2 = 0;
210 G4double ceff2 =0;
211 G4double rho = 0;
212 if ((theParticle == theAProton) || (theParticle == theANeutron))
213 {
214 if(theDef == theProton)
215 {
216// G4double Mp2=sqr(theDef->GetPDGMass()/GeV );
217
218// change 30 October
219
220 if(Plab < 610.)
221 { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.-
222 13.634*Plab/1000.*Plab/1000.*Plab/1000. ;}
223 if((Plab < 5500.)&&(Plab >= 610.) )
224 { rho = 0.22; }
225 if((Plab >= 5500.)&&(Plab < 12300.) )
226 { rho = -0.32; }
227 if( Plab >= 12300.)
228 { rho = 0.135-2.26/(std::sqrt(S)) ;}
229
230 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(S-4.*0.88))+0.04*std::log(S) ;
231 ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)+1.5) ;
232
233/*
234 Ref2=0.8/std::sqrt(std::sqrt(S-4.*Mp2)) + 0.55;
235 if(S>1000.) Ref2=0.62+0.02*std::log(S) ;
236 ceff2 = 0.035/(sqr(S-4.3)+0.4) + 0.085 * std::log(S) ;
237 if(S>1000.) ceff2 = 0.005 * std::log(S) + 0.29;
238*/
239
240 Ref2=Ref2*Ref2;
241 ceff2 = ceff2*ceff2;
242
243 SlopeMag = 0.5; // Uzhi
244 Amag= 1.; // Uzhi
245 }
246
247 if(Z>2)
248 { Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa;
249 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*std::exp(-0.03*sig_pbarp);
250 }
251 if( (Z==2)&&(A==4) )
252 { Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp;
253 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp);
254 }
255 if( (Z==1)&&(A==3) )
256 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
257 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
258 }
259 if( (Z==2)&&(A==3) )
260 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
261 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
262 }
263 if( (Z==1)&&(A==2) )
264 {
265 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp;
266 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp);
267 }
268 }
269
270 if (theParticle == theADeuteron)
271 {
272 sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy/2.);
273 Ref2 = XstotalHad/10./2./pi ;
274 if(Z>2)
275 {
276 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 * std::exp(-0.03*sig_pbarp);
277 }
278 if(theDef == theProton)
279 {
280 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp);
281 }
282 if(theDef == theDeuteron)
283 {
284 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 * std::exp(-0.03*sig_pbarp);
285 }
286 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
287 {
288 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp);
289 }
290 if(theDef == theAlpha)
291 {
292 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp);
293 }
294 }
295
296 if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
297 {
298 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
299 Ref2 = XstotalHad/10./2./pi ;
300 if(Z>2)
301 {
302 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*std::exp(-0.03*sig_pbarp);
303 }
304 if(theDef == theProton)
305 {
306 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
307 }
308 if(theDef == theDeuteron)
309 {
310 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp);
311 }
312 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
313 {
314 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 * std::exp(-0.02*sig_pbarp);
315 }
316 if(theDef == theAlpha)
317 {
318 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp);
319 }
320 }
321
322
323 if (theParticle == theAAlpha)
324 {
325 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
326 Ref2 = XstotalHad/10./2./pi ;
327 if(Z>2)
328 {
329 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 * std::exp(-0.03*sig_pbarp);
330 }
331 if(theDef == theProton)
332 {
333 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp);
334 }
335 if(theDef == theDeuteron)
336 {
337 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp);
338 }
339 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
340 {
341 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp);
342 }
343 if(theDef == theAlpha)
344 {
345 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 * std::exp(-0.03*sig_pbarp);
346 }
347 }
348
349 fRef=std::sqrt(Ref2);
350 fceff = std::sqrt(ceff2);
351// G4cout<<" Ref "<<fRef<<" c_eff "<<fceff<< " rho "<< rho<<G4endl;
352
353
354 G4double Q = 0.0 ;
355 G4double BracFunct;
356 do
357 {
358 Q = -std::log(1.-(1.- std::exp(-SlopeMag * Qmax))* G4UniformRand() )/SlopeMag;
359 G4double x = fRef * Q;
360 BracFunct = ( ( sqr(BesselOneByArg(x))+sqr(rho/2. * BesselJzero(x)) )
361* sqr(DampFactor(pi*fceff*Q))) /(Amag*std::exp(-SlopeMag*Q));
362
363 BracFunct = BracFunct * Q * sqr(sqr(fRef));
364 }
365 while (G4UniformRand()>BracFunct);
366
367 T= sqr(Q);
368 T*=3.893913e+4; // fm -> MeV^2
369 }
370
371 G4double cosTet=1.0-T/(2.*ptot*ptot);
372 if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov.
373 if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov.
374 fTetaCMS=std::acos(cosTet);
375
376 return T;
377}
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::mag
double mag() const
G4AntiNuclElastic::BesselOneByArg
G4double BesselOneByArg(G4double z)
Definition: G4AntiNuclElastic.cc:630
G4AntiNuclElastic::CalculateParticleBeta
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
Definition: G4AntiNuclElastic.cc:493
G4AntiNuclElastic::CalculateZommerfeld
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Definition: G4AntiNuclElastic.cc:507
G4AntiNuclElastic::DampFactor
G4double DampFactor(G4double z)
Definition: G4AntiNuclElastic.cc:473
G4AntiNuclElastic::GetcosTeta1
G4double GetcosTeta1(G4double plab, G4int A)
Definition: G4AntiNuclElastic.cc:649
G4AntiNuclElastic::BesselJzero
G4double BesselJzero(G4double z)
Definition: G4AntiNuclElastic.cc:532
G4AntiNuclElastic::CalculateAm
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Definition: G4AntiNuclElastic.cc:516
G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentAntiNuclNuclearXS.cc:96
G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentAntiNuclNuclearXS.cc:272
G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc
G4double GetAntiHadronNucleonTotCrSc(const G4ParticleDefinition *aParticle, G4double kinEnergy)
Definition: G4ComponentAntiNuclNuclearXS.cc:295
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:53
G4Pow::Z13
G4double Z13(G4int Z)
Definition: G4Pow.hh:110
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70
sqr
T sqr(const T &x)
Definition: templates.hh:145

Referenced by SampleThetaCMS(), and SampleThetaLab().

◆ SampleThetaCMS()

G4double G4AntiNuclElastic::SampleThetaCMS ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 381 of file G4AntiNuclElastic.cc.

383{
384 G4double T;
385 T = SampleInvariantT( p, plab, Z, A);
386
387 // NaN finder
388 if(!(T < 0.0 || T >= 0.0))
389 {
390 if (verboseLevel > 0)
391 {
392 G4cout << "G4DiffuseElastic:WARNING: A = " << A
393 << " mom(GeV)= " << plab/GeV
394 << " S-wave will be sampled"
395 << G4endl;
396 }
397 T = G4UniformRand()*fTmax;
398
399 }
400
401 if(fptot > 0.) // Uzhi 24 Nov. 2011
402 {
403 G4double cosTet=1.0-T/(2.*fptot*fptot);
404 if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov.
405 if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov.
406 fTetaCMS=std::acos(cosTet);
407 return fTetaCMS;
408 } else // Uzhi 24 Nov. 2011
409 { // Uzhi 24 Nov. 2011
410 return 2.*G4UniformRand()-1.; // Uzhi 24 Nov. 2011
411 } // Uzhi 24 Nov. 2011
412}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4AntiNuclElastic::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4AntiNuclElastic.cc:99

◆ SampleThetaLab()

G4double G4AntiNuclElastic::SampleThetaLab ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 417 of file G4AntiNuclElastic.cc.

419{
420 G4double T;
421 T = SampleInvariantT( p, plab, Z, A);
422
423 // NaN finder
424 if(!(T < 0.0 || T >= 0.0))
425 {
426 if (verboseLevel > 0)
427 {
428 G4cout << "G4DiffuseElastic:WARNING: A = " << A
429 << " mom(GeV)= " << plab/GeV
430 << " S-wave will be sampled"
431 << G4endl;
432 }
433 T = G4UniformRand()*fTmax;
434 }
435
436 G4double phi = G4UniformRand()*twopi;
437
438 G4double cost(1.);
439 if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;} // Uzhi 24 Nov. 2011
440
441 G4double sint;
442 if( cost >= 1.0 )
443 {
444 cost = 1.0;
445 sint = 0.0;
446 }
447 else if( cost <= -1.0)
448 {
449 cost = -1.0;
450 sint = 0.0;
451 }
452 else
453 {
454 sint = std::sqrt((1.0-cost)*(1.0+cost));
455 }
456
457 G4double m1 = p->GetPDGMass();
458 G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost);
459 v *= fptot;
460 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1));
461
462 nlv.boost(fbst);
463
464 G4ThreeVector np = nlv.vect();
465 G4double theta = np.theta();
466 fThetaLab = theta;
467
468 return theta;
469}
CLHEP::Hep3Vector::theta
double theta() const
The documentation for this class was generated from the following files: