Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4AntiNuclElastic Class Reference

#include <G4AntiNuclElastic.hh>

Inheritance diagram for G4AntiNuclElastic:

Public Member Functions

 G4AntiNuclElastic ()
 ~G4AntiNuclElastic () override
G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
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")
 ~G4HadronElastic () override
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) override
G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
G4double GetSlopeCof (const G4int pdg)
void SetLowestEnergyLimit (G4double value)
G4double LowestEnergyLimit () const
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
void ModelDescription (std::ostream &) const override
Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
virtual ~G4HadronicInteraction ()
virtual G4bool IsApplicable (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
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)
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
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels () const
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels () const
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
virtual void InitialiseModel ()
 G4HadronicInteraction (const G4HadronicInteraction &right)=delete
const G4HadronicInteraction & operator= (const G4HadronicInteraction &right)=delete
G4bool operator== (const G4HadronicInteraction &right) const =delete
G4bool operator!= (const G4HadronicInteraction &right) const =delete

Additional Inherited Members

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

Detailed Description

Definition at line 46 of file G4AntiNuclElastic.hh.

Constructor & Destructor Documentation

◆ G4AntiNuclElastic()

G4AntiNuclElastic::G4AntiNuclElastic ( )
explicit

Definition at line 56 of file G4AntiNuclElastic.cc.

57 : G4HadronElastic("AntiAElastic")
58{
59 //V.Ivanchenko commented out
60 //SetMinEnergy( 0.1*GeV );
61 //SetMaxEnergy( 10.*TeV );
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 G4CrossSectionDataSetRegistry* reg = G4CrossSectionDataSetRegistry::Instance();
76 cs = static_cast<G4ComponentAntiNuclNuclearXS*>(reg->GetComponentCrossSection("AntiAGlauber"));
77 if(!cs) { cs = new G4ComponentAntiNuclNuclearXS(); }
78
79 fParticle = 0;
80 fWaveVector = 0.;
81 fBeta = 0.;
82 fZommerfeld = 0.;
83 fAm = 0.;
84 fTetaCMS = 0.;
85 fRa = 0.;
86 fRef = 0.;
87 fceff = 0.;
88 fptot = 0.;
89 fTmax = 0.;
90 fThetaLab = 0.;
91}
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition G4Alpha.cc:83
G4AntiAlpha::AntiAlpha
static G4AntiAlpha * AntiAlpha()
Definition G4AntiAlpha.cc:84
G4AntiDeuteron::AntiDeuteron
static G4AntiDeuteron * AntiDeuteron()
Definition G4AntiDeuteron.cc:90
G4AntiHe3::AntiHe3
static G4AntiHe3 * AntiHe3()
Definition G4AntiHe3.cc:90
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition G4AntiNeutron.cc:98
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition G4AntiProton.cc:90
G4AntiTriton::AntiTriton
static G4AntiTriton * AntiTriton()
Definition G4AntiTriton.cc:90
G4CrossSectionDataSetRegistry::GetComponentCrossSection
G4VComponentCrossSection * GetComponentCrossSection(const G4String &name)
Definition G4CrossSectionDataSetRegistry.cc:185
G4CrossSectionDataSetRegistry::Instance
static G4CrossSectionDataSetRegistry * Instance()
Definition G4CrossSectionDataSetRegistry.cc:70
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition G4Deuteron.cc:90
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition G4HadronElastic.cc:49
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition G4Neutron.cc:101
G4Proton::Proton
static G4Proton * Proton()
Definition G4Proton.cc:90

◆ ~G4AntiNuclElastic()

G4AntiNuclElastic::~G4AntiNuclElastic ( )
override

Definition at line 94 of file G4AntiNuclElastic.cc.

95{}

Member Function Documentation

◆ BesselJone()

G4double G4AntiNuclElastic::BesselJone ( G4double z)

Definition at line 542 of file G4AntiNuclElastic.cc.

543{
544 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
545
546 modvalue = std::fabs(value);
547
548 if ( modvalue < 8.0 )
549 {
550 value2 = value*value;
551 fact1 = value*(72362614232.0 + value2*(-7895059235.0
552 + value2*( 242396853.1
553 + value2*(-2972611.439
554 + value2*( 15704.48260
555 + value2*(-30.16036606 ) ) ) ) ) );
556
557 fact2 = 144725228442.0 + value2*(2300535178.0
558 + value2*(18583304.74
559 + value2*(99447.43394
560 + value2*(376.9991397
561 + value2*1.0 ) ) ) );
562 bessel = fact1/fact2;
563 }
564 else
565 {
566 arg = 8.0/modvalue;
567 value2 = arg*arg;
568
569 shift = modvalue - 2.356194491;
570
571 fact1 = 1.0 + value2*( 0.183105e-2
572 + value2*(-0.3516396496e-4
573 + value2*(0.2457520174e-5
574 + value2*(-0.240337019e-6 ) ) ) );
575
576 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
577 + value2*( 0.8449199096e-5
578 + value2*(-0.88228987e-6
579 + value2*0.105787412e-6 ) ) );
580
581 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
582 if (value < 0.0) bessel = -bessel;
583 }
584 return bessel;
585}
G4double
double G4double
Definition G4Types.hh:83

Referenced by BesselOneByArg().

◆ BesselJzero()

G4double G4AntiNuclElastic::BesselJzero ( G4double z)

Definition at line 491 of file G4AntiNuclElastic.cc.

492{
493 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
494
495 modvalue = std::fabs(value);
496
497 if ( value < 8.0 && value > -8.0 )
498 {
499 value2 = value*value;
500
501 fact1 = 57568490574.0 + value2*(-13362590354.0
502 + value2*( 651619640.7
503 + value2*(-11214424.18
504 + value2*( 77392.33017
505 + value2*(-184.9052456 ) ) ) ) );
506
507 fact2 = 57568490411.0 + value2*( 1029532985.0
508 + value2*( 9494680.718
509 + value2*(59272.64853
510 + value2*(267.8532712
511 + value2*1.0 ) ) ) );
512
513 bessel = fact1/fact2;
514 }
515 else
516 {
517 arg = 8.0/modvalue;
518
519 value2 = arg*arg;
520
521 shift = modvalue-0.785398164;
522
523 fact1 = 1.0 + value2*(-0.1098628627e-2
524 + value2*(0.2734510407e-4
525 + value2*(-0.2073370639e-5
526 + value2*0.2093887211e-6 ) ) );
527 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
528 + value2*(-0.6911147651e-5
529 + value2*(0.7621095161e-6
530 - value2*0.934945152e-7 ) ) );
531
532 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
533 }
534 return bessel;
535}

Referenced by SampleInvariantT().

◆ BesselOneByArg()

G4double G4AntiNuclElastic::BesselOneByArg ( G4double z)

Definition at line 589 of file G4AntiNuclElastic.cc.

590{
591 G4double x2, result;
592
593 if( std::fabs(x) < 0.01 )
594 {
595 x *= 0.5;
596 x2 = x*x;
597 result = (2.- x2 + x2*x2/6.)/4.;
598 }
599 else
600 {
601 result = BesselJone(x)/x;
602 }
603 return result;
604}
G4AntiNuclElastic::BesselJone
G4double BesselJone(G4double z)
Definition G4AntiNuclElastic.cc:542

Referenced by SampleInvariantT().

◆ CalculateAm()

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

Definition at line 475 of file G4AntiNuclElastic.cc.

476{
477 G4double k = momentum/hbarc;
478 G4double ch = 1.13 + 3.76*n*n;
479 G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius;
480 G4double zn2 = zn*zn;
481 fAm = ch/zn2;
482
483 return fAm;
484}
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition G4Pow.cc:41
G4Pow::A13
G4double A13(G4double A) const
Definition G4Pow.cc:116

Referenced by SampleInvariantT().

◆ CalculateParticleBeta()

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

Definition at line 452 of file G4AntiNuclElastic.cc.

454{
455 G4double mass = particle->GetPDGMass();
456 G4double a = momentum/mass;
457 fBeta = a/std::sqrt(1+a*a);
458
459 return fBeta;
460}

Referenced by SampleInvariantT().

◆ CalculateZommerfeld()

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

Definition at line 466 of file G4AntiNuclElastic.cc.

467{
468 fZommerfeld = fine_structure_const*Z1*Z2/beta;
469
470 return fZommerfeld;
471}

Referenced by SampleInvariantT().

◆ DampFactor()

G4double G4AntiNuclElastic::DampFactor ( G4double z)

Definition at line 432 of file G4AntiNuclElastic.cc.

433{
434 G4double df;
435 G4double f3 = 6.; // first factorials
436
437 if( std::fabs(x) < 0.01 )
438 {
439 df=1./(1.+x*x/f3);
440 }
441 else
442 {
443 df = x/std::sinh(x);
444 }
445 return df;
446}

Referenced by SampleInvariantT().

◆ GetComponentCrossSection()

G4ComponentAntiNuclNuclearXS * G4AntiNuclElastic::GetComponentCrossSection ( )
inline

Definition at line 121 of file G4AntiNuclElastic.hh.

122{
123 return cs;
124}

Referenced by LBE::ConstructHad().

◆ GetcosTeta1()

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

Definition at line 608 of file G4AntiNuclElastic.cc.

609{
610
611// G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
612 G4double p0 = 1.*hbarc/fermi;
613//G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab);
614 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab);
615//////////////////
616 if(cteta1 < -1.) cteta1 = -1.0;
617 return cteta1;
618}
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4Pow::Z23
G4double Z23(G4int Z) const
Definition G4Pow.hh:125

Referenced by SampleInvariantT().

◆ SampleInvariantT()

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

Reimplemented from G4HadronicInteraction.

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 * theTargetDef = 0;
112
113 if (Z == 1 && A == 1) theTargetDef = theProton;
114 else if (Z == 1 && A == 2) theTargetDef = theDeuteron;
115 else if (Z == 1 && A == 3) theTargetDef = G4Triton::Triton();
116 else if (Z == 2 && A == 3) theTargetDef = G4He3::He3();
117 else if (Z == 2 && A == 4) theTargetDef = 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; // In (MeV/c)^2
137
138 if(Plab < (std::abs(particle->GetBaryonNumber())*100)*MeV)
139 {return fTmax*G4UniformRand();}
140
141 // Calculation of NN collision properties
142 G4double PlabPerN = Plab/std::abs(theParticle->GetBaryonNumber());
143 G4double NucleonMass = 0.5*( theProton->GetPDGMass() + theNeutron->GetPDGMass() );
144 G4double PrNucleonMass(0.); // Projectile average nucleon mass
145 if( std::abs(theParticle->GetBaryonNumber()) == 1 ) { PrNucleonMass = theParticle->GetPDGMass(); }
146 else { PrNucleonMass = NucleonMass; }
147 G4double energyPerN = std::sqrt( sqr(PlabPerN) + sqr(PrNucleonMass));
148 energyPerN -= PrNucleonMass;
149 //---
150
151 G4double Z1 = particle->GetPDGCharge();
152 G4double Z2 = Z;
153
154 G4double beta = CalculateParticleBeta(particle, ptot);
155 G4double n = CalculateZommerfeld( beta, Z1, Z2 );
156 G4double Am = CalculateAm( ptot, n, Z2 );
157 fWaveVector = ptot; // /hbarc;
158
159 G4LorentzVector Fproj(0.,0.,0.,0.);
160 const G4double mevToBarn = 0.38938e+6;
161 G4double XsCoulomb = mevToBarn*sqr(n/fWaveVector)*pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1);
162
163 G4double XsElastHadronic =cs->GetElasticElementCrossSection(particle, energy, Z, (G4double)A);
164 G4double XsTotalHadronic =cs->GetTotalElementCrossSection(particle, energy, Z, (G4double)A);
165
166 XsElastHadronic/=millibarn; XsTotalHadronic/=millibarn;
167
168 G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHadronic);
169
170 if(G4UniformRand() < CoulombProb)
171 { // Simulation of Coulomb scattering
172
173 G4double phi = twopi * G4UniformRand();
174 G4double Ksi = G4UniformRand();
175
176 G4double par1 = 2.*(1.+Am)/(1.+ctet1);
177
178 // ////sample ThetaCMS in Coulomb part
179
180 G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi);
181
182 G4double PtZ=ptot*cosThetaCMS;
183 Fproj.setPz(PtZ);
184 G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS);
185 G4double PtX= PtProjCMS * std::cos(phi);
186 G4double PtY= PtProjCMS * std::sin(phi);
187 Fproj.setPx(PtX);
188 Fproj.setPy(PtY);
189 Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj));
190 T = -(Pproj-Fproj).mag2();
191 }
192 else
193 {
194 // Simulation of strong interaction scattering
195
196 G4double Qmax = 2.*ptot/197.33; // in fm^-1
197
198 G4double Amag = 1.0; // A1 in Majorant funct:A1*exp(-q*A2)
199 G4double SlopeMag = 0.5; // A2 in Majorant funct:A1*exp(-q*A2)
200
201 G4double sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(theAProton,energyPerN); //mb
202
203 fRa = 1.113*G4Pow::GetInstance()->Z13(A) -
204 0.227/G4Pow::GetInstance()->Z13(A);
205 if(A == 3) fRa=1.81;
206 if(A == 4) fRa=1.37;
207
208 if((A>=12.) && (A<27) ) fRa=fRa*0.85;
209 if((A>=27.) && (A<48) ) fRa=fRa*0.90;
210 if((A>=48.) && (A<65) ) fRa=fRa*0.95;
211
212 G4double Ref2 = XsTotalHadronic/10./2./pi; // in fm^2
213 G4double ceff2 = 0.0;
214 G4double rho = 0.0;
215
216 if ((theParticle == theAProton) || (theParticle == theANeutron))
217 {
218 if(theTargetDef == theProton)
219 {
220 // Determination of the real part of Pbar+N amplitude
221 if(Plab < 610.)
222 { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.-
223 13.634*Plab/1000.*Plab/1000.*Plab/1000. ;}
224 if((Plab < 5500.)&&(Plab >= 610.) )
225 { rho = 0.22; }
226 if((Plab >= 5500.)&&(Plab < 12300.) )
227 { rho = -0.32; }
228 if( Plab >= 12300.)
229 { rho = 0.135-2.26/(std::sqrt(S)) ;}
230 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(S-4.*0.88))+0.04*G4Log(S) ;
231 ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)+1.5) ;
232 Ref2 =Ref2*Ref2;
233 ceff2 = ceff2*ceff2;
234 }
235
236 if( (Z==1)&&(A==2) )
237 {
238 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp;
239 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*G4Exp(-0.03*sig_pbarp);
240 }
241 if( (Z==1)&&(A==3) )
242 {
243 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
244 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*G4Exp(-0.03*sig_pbarp);
245 }
246 if( (Z==2)&&(A==3) )
247 {
248 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
249 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*G4Exp(-0.03*sig_pbarp);
250 }
251 if( (Z==2)&&(A==4) )
252 {
253 Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp;
254 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*G4Exp(-0.03*sig_pbarp);
255 }
256 if(Z>2)
257 {
258 Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa;
259 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*G4Exp(-0.03*sig_pbarp);
260 }
261 } // End of if ((theParticle == theAProton) || (theParticle == theANeutron))
262
263 if (theParticle == theADeuteron)
264 {
265 if(theTargetDef == theProton)
266 {
267 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*G4Exp(-0.03*sig_pbarp);
268 }
269 if(theTargetDef == theDeuteron)
270 {
271 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 * G4Exp(-0.03*sig_pbarp);
272 }
273 if( (theTargetDef == G4Triton::Triton()) || (theTargetDef == G4He3::He3() ) )
274 {
275 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * G4Exp(-0.02*sig_pbarp);
276 }
277 if(theTargetDef == theAlpha)
278 {
279 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * G4Exp(-0.02*sig_pbarp);
280 }
281 if(Z>2)
282 {
283 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 * G4Exp(-0.03*sig_pbarp);
284 }
285 }
286
287 if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
288 {
289 if(theTargetDef == theProton)
290 {
291 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*G4Exp(-0.03*sig_pbarp);
292 }
293 if(theTargetDef == theDeuteron)
294 {
295 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * G4Exp(-0.02*sig_pbarp);
296 }
297 if( (theTargetDef == G4Triton::Triton()) || (theTargetDef == G4He3::He3() ) )
298 {
299 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 * G4Exp(-0.02*sig_pbarp);
300 }
301 if(theTargetDef == theAlpha)
302 {
303 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * G4Exp(-0.03*sig_pbarp);
304 }
305 if(Z>2)
306 {
307 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*G4Exp(-0.03*sig_pbarp);
308 }
309 }
310
311 if ( (theParticle == theAAlpha) || (ceff2 == 0.0) )
312 {
313 if(theTargetDef == theProton)
314 {
315 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*G4Exp(-0.03*sig_pbarp);
316 }
317 if(theTargetDef == theDeuteron)
318 {
319 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * G4Exp(-0.02*sig_pbarp);
320 }
321 if( (theTargetDef == G4Triton::Triton()) || (theTargetDef == G4He3::He3() ) )
322 {
323 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * G4Exp(-0.03*sig_pbarp);
324 }
325 if(theTargetDef == theAlpha)
326 {
327 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 * G4Exp(-0.03*sig_pbarp);
328 }
329 if(Z>2)
330 {
331 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 * G4Exp(-0.03*sig_pbarp);
332 }
333 }
334
335 fRef=std::sqrt(Ref2);
336 fceff = std::sqrt(ceff2);
337
338 G4double Q = 0.0 ;
339 G4double BracFunct;
340
341 const G4int maxNumberOfLoops = 10000;
342 G4int loopCounter = 0;
343 do
344 {
345 Q = -G4Log(1.-(1.- G4Exp(-SlopeMag * Qmax))* G4UniformRand() )/SlopeMag;
346 G4double x = fRef * Q;
347 BracFunct = ( ( sqr(BesselOneByArg(x))+sqr(rho/2. * BesselJzero(x)) )
348 * sqr(DampFactor(pi*fceff*Q))) /(Amag*G4Exp(-SlopeMag*Q));
349 BracFunct = BracFunct * Q;
350 }
351 while ( (G4UniformRand()>BracFunct) &&
352 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
353 if ( loopCounter >= maxNumberOfLoops ) {
354 fTetaCMS = 0.0;
355 return 0.0;
356 }
357
358 T= sqr(Q);
359 T*=3.893913e+4; // fm^(-2) -> MeV^2
360
361 } // End of simulation of strong interaction scattering
362
363 return T;
364}
S
G4double S(G4double temp)
Definition G4DNAElectronHoleRecombination.cc:79
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:169
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition G4LorentzVector.hh:33
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4int
int G4int
Definition G4Types.hh:85
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::mag
double mag() const
G4AntiNuclElastic::BesselOneByArg
G4double BesselOneByArg(G4double z)
Definition G4AntiNuclElastic.cc:589
G4AntiNuclElastic::CalculateParticleBeta
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
Definition G4AntiNuclElastic.cc:452
G4AntiNuclElastic::CalculateZommerfeld
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Definition G4AntiNuclElastic.cc:466
G4AntiNuclElastic::DampFactor
G4double DampFactor(G4double z)
Definition G4AntiNuclElastic.cc:432
G4AntiNuclElastic::GetcosTeta1
G4double GetcosTeta1(G4double plab, G4int A)
Definition G4AntiNuclElastic.cc:608
G4AntiNuclElastic::BesselJzero
G4double BesselJzero(G4double z)
Definition G4AntiNuclElastic.cc:491
G4AntiNuclElastic::CalculateAm
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Definition G4AntiNuclElastic.cc:475
G4He3::He3
static G4He3 * He3()
Definition G4He3.cc:90
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition G4NucleiProperties.cc:49
G4Pow::Z13
G4double Z13(G4int Z) const
Definition G4Pow.hh:123
G4Triton::Triton
static G4Triton * Triton()
Definition G4Triton.cc:90
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185
G4INCL::Math::pi
const G4double pi
Definition G4INCLGlobals.hh:69
sqr
T sqr(const T &x)
Definition templates.hh:128

Referenced by SampleThetaCMS(), and SampleThetaLab().

◆ SampleThetaCMS()

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

Definition at line 368 of file G4AntiNuclElastic.cc.

370{
371 G4double T = SampleInvariantT( p, plab, Z, A);
372 if (T <= 0.0 || fTmax <= 0.0) { return 1.0; }
373
374 // NaN finder substituted by simple check
375 if (T > fTmax)
376 {
377 if (verboseLevel > 0)
378 {
379 G4cout << "G4AntiNuclElastic::SampleThetaCMS WARNING: A = " << A
380 << " mom(GeV)=" << plab/GeV << " t(GeV2)=" << T/(GeV*GeV)
381 << " > Tmax(GeV2)=" << fTmax/(GeV*GeV) << " S-wave will be sampled"
382 << G4endl;
383 }
384 T = G4UniformRand()*fTmax;
385 }
386 G4double cosTet = 1.0 - 2*T/fTmax;
387 return cosTet;
388}
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4AntiNuclElastic::SampleInvariantT
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
Definition G4AntiNuclElastic.cc:99

◆ SampleThetaLab()

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

Definition at line 393 of file G4AntiNuclElastic.cc.

395{
396 G4double T = SampleInvariantT( p, plab, Z, A);
397 if (T <= 0.0 || fTmax <= 0.0) { return 1.0; }
398
399 // NaN finder substituted by simple check
400 if (T > fTmax)
401 {
402 if (verboseLevel > 0)
403 {
404 G4cout << "G4AntiNuclElastic::SampleThetaLab WARNING: A = " << A
405 << " mom(GeV)=" << plab/GeV << " t(GeV2)=" << T/(GeV*GeV)
406 << " > Tmax(GeV2)=" << fTmax/(GeV*GeV) << " S-wave will be sampled"
407 << G4endl;
408 }
409 T = G4UniformRand()*fTmax;
410 }
411
412 G4double cost = 1.0 - 2*T/fTmax;
413 G4double phi = G4UniformRand()*twopi;
414 G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
415
416 G4double m1 = p->GetPDGMass();
417 G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost);
418 v *= fptot;
419 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1));
420
421 nlv.boost(fbst);
422
423 G4ThreeVector np = nlv.vect();
424 G4double theta = np.theta();
425 fThetaLab = theta;
426
427 return theta;
428}
CLHEP::Hep3Vector::theta
double theta() const
The documentation for this class was generated from the following files: