Geant4 11.1.1
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)
 
G4ComponentAntiNuclNuclearXSGetComponentCrossSection ()
 
- Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 
 ~G4HadronElastic () override
 
G4HadFinalStateApplyYourself (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 G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
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 G4StringGetModelName () 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, G4doubleGetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4doubleGetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 
 G4HadronicInteraction (const G4HadronicInteraction &right)=delete
 
const G4HadronicInteractionoperator= (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
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}
static G4Alpha * Alpha()
Definition: G4Alpha.cc:88
static G4AntiAlpha * AntiAlpha()
Definition: G4AntiAlpha.cc:88
static G4AntiDeuteron * AntiDeuteron()
static G4AntiHe3 * AntiHe3()
Definition: G4AntiHe3.cc:93
static G4AntiNeutron * AntiNeutron()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:92
static G4AntiTriton * AntiTriton()
Definition: G4AntiTriton.cc:93
G4VComponentCrossSection * GetComponentCrossSection(const G4String &name)
static G4CrossSectionDataSetRegistry * Instance()
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
static G4Proton * Proton()
Definition: G4Proton.cc:92

◆ ~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 570 of file G4AntiNuclElastic.cc.

571{
572 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
573
574 modvalue = std::fabs(value);
575
576 if ( modvalue < 8.0 )
577 {
578 value2 = value*value;
579 fact1 = value*(72362614232.0 + value2*(-7895059235.0
580 + value2*( 242396853.1
581 + value2*(-2972611.439
582 + value2*( 15704.48260
583 + value2*(-30.16036606 ) ) ) ) ) );
584
585 fact2 = 144725228442.0 + value2*(2300535178.0
586 + value2*(18583304.74
587 + value2*(99447.43394
588 + value2*(376.9991397
589 + value2*1.0 ) ) ) );
590 bessel = fact1/fact2;
591 }
592 else
593 {
594 arg = 8.0/modvalue;
595 value2 = arg*arg;
596
597 shift = modvalue - 2.356194491;
598
599 fact1 = 1.0 + value2*( 0.183105e-2
600 + value2*(-0.3516396496e-4
601 + value2*(0.2457520174e-5
602 + value2*(-0.240337019e-6 ) ) ) );
603
604 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
605 + value2*( 0.8449199096e-5
606 + value2*(-0.88228987e-6
607 + value2*0.105787412e-6 ) ) );
608
609 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
610 if (value < 0.0) bessel = -bessel;
611 }
612 return bessel;
613}
double G4double
Definition: G4Types.hh:83

Referenced by BesselOneByArg().

◆ BesselJzero()

G4double G4AntiNuclElastic::BesselJzero ( G4double  z)

Definition at line 519 of file G4AntiNuclElastic.cc.

520{
521 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
522
523 modvalue = std::fabs(value);
524
525 if ( value < 8.0 && value > -8.0 )
526 {
527 value2 = value*value;
528
529 fact1 = 57568490574.0 + value2*(-13362590354.0
530 + value2*( 651619640.7
531 + value2*(-11214424.18
532 + value2*( 77392.33017
533 + value2*(-184.9052456 ) ) ) ) );
534
535 fact2 = 57568490411.0 + value2*( 1029532985.0
536 + value2*( 9494680.718
537 + value2*(59272.64853
538 + value2*(267.8532712
539 + value2*1.0 ) ) ) );
540
541 bessel = fact1/fact2;
542 }
543 else
544 {
545 arg = 8.0/modvalue;
546
547 value2 = arg*arg;
548
549 shift = modvalue-0.785398164;
550
551 fact1 = 1.0 + value2*(-0.1098628627e-2
552 + value2*(0.2734510407e-4
553 + value2*(-0.2073370639e-5
554 + value2*0.2093887211e-6 ) ) );
555 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
556 + value2*(-0.6911147651e-5
557 + value2*(0.7621095161e-6
558 - value2*0.934945152e-7 ) ) );
559
560 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
561 }
562 return bessel;
563}

Referenced by SampleInvariantT().

◆ BesselOneByArg()

G4double G4AntiNuclElastic::BesselOneByArg ( G4double  z)

Definition at line 617 of file G4AntiNuclElastic.cc.

618{
619 G4double x2, result;
620
621 if( std::fabs(x) < 0.01 )
622 {
623 x *= 0.5;
624 x2 = x*x;
625 result = (2.- x2 + x2*x2/6.)/4.;
626 }
627 else
628 {
629 result = BesselJone(x)/x;
630 }
631 return result;
632}
G4double BesselJone(G4double z)

Referenced by SampleInvariantT().

◆ CalculateAm()

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

Definition at line 503 of file G4AntiNuclElastic.cc.

504{
505 G4double k = momentum/hbarc;
506 G4double ch = 1.13 + 3.76*n*n;
507 G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius;
508 G4double zn2 = zn*zn;
509 fAm = ch/zn2;
510
511 return fAm;
512}
const G4int Z[17]
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double A13(G4double A) const
Definition: G4Pow.cc:116

Referenced by SampleInvariantT().

◆ CalculateParticleBeta()

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

Definition at line 480 of file G4AntiNuclElastic.cc.

482{
483 G4double mass = particle->GetPDGMass();
484 G4double a = momentum/mass;
485 fBeta = a/std::sqrt(1+a*a);
486
487 return fBeta;
488}

Referenced by SampleInvariantT().

◆ CalculateZommerfeld()

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

Definition at line 494 of file G4AntiNuclElastic.cc.

495{
496 fZommerfeld = fine_structure_const*Z1*Z2/beta;
497
498 return fZommerfeld;
499}

Referenced by SampleInvariantT().

◆ DampFactor()

G4double G4AntiNuclElastic::DampFactor ( G4double  z)

Definition at line 460 of file G4AntiNuclElastic.cc.

461{
462 G4double df;
463 G4double f3 = 6.; // first factorials
464
465 if( std::fabs(x) < 0.01 )
466 {
467 df=1./(1.+x*x/f3);
468 }
469 else
470 {
471 df = x/std::sinh(x);
472 }
473 return df;
474}

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 636 of file G4AntiNuclElastic.cc.

637{
638
639// G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
640 G4double p0 = 1.*hbarc/fermi;
641//G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab);
642 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab);
643//////////////////
644 if(cteta1 < -1.) cteta1 = -1.0;
645 return cteta1;
646}
const G4double A[17]
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
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}
G4double S(G4double temp)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double G4Log(G4double x)
Definition: G4Log.hh:227
int G4int
Definition: G4Types.hh:85
#define G4UniformRand()
Definition: Randomize.hh:52
double mag() const
G4double BesselOneByArg(G4double z)
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double DampFactor(G4double z)
G4double GetcosTeta1(G4double plab, G4int A)
G4double BesselJzero(G4double z)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
G4double GetAntiHadronNucleonTotCrSc(const G4ParticleDefinition *aParticle, G4double kinEnergy)
static G4He3 * He3()
Definition: G4He3.cc:93
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double Z13(G4int Z) const
Definition: G4Pow.hh:123
static G4Triton * Triton()
Definition: G4Triton.cc:93
G4double energy(const ThreeVector &p, const G4double m)
const G4double pi
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;
372 T = SampleInvariantT( p, plab, Z, A);
373
374 // NaN finder
375 if(!(T < 0.0 || T >= 0.0))
376 {
377 if (verboseLevel > 0)
378 {
379 G4cout << "G4DiffuseElastic:WARNING: A = " << A
380 << " mom(GeV)= " << plab/GeV
381 << " S-wave will be sampled"
382 << G4endl;
383 }
384 T = G4UniformRand()*fTmax;
385
386 }
387
388 if(fptot > 0.)
389 {
390 G4double cosTet=1.0-T/(2.*fptot*fptot);
391 if(cosTet > 1.0 ) cosTet= 1.;
392 if(cosTet < -1.0 ) cosTet=-1.;
393 fTetaCMS=std::acos(cosTet);
394 return fTetaCMS;
395 } else
396 {
397 return 2.*G4UniformRand()-1.;
398 }
399}
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override

◆ SampleThetaLab()

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

Definition at line 404 of file G4AntiNuclElastic.cc.

406{
407 G4double T;
408 T = SampleInvariantT( p, plab, Z, A);
409
410 // NaN finder
411 if(!(T < 0.0 || T >= 0.0))
412 {
413 if (verboseLevel > 0)
414 {
415 G4cout << "G4DiffuseElastic:WARNING: A = " << A
416 << " mom(GeV)= " << plab/GeV
417 << " S-wave will be sampled"
418 << G4endl;
419 }
420 T = G4UniformRand()*fTmax;
421 }
422
423 G4double phi = G4UniformRand()*twopi;
424
425 G4double cost(1.);
426 if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;}
427
428 G4double sint;
429 if( cost >= 1.0 )
430 {
431 cost = 1.0;
432 sint = 0.0;
433 }
434 else if( cost <= -1.0)
435 {
436 cost = -1.0;
437 sint = 0.0;
438 }
439 else
440 {
441 sint = std::sqrt((1.0-cost)*(1.0+cost));
442 }
443
444 G4double m1 = p->GetPDGMass();
445 G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost);
446 v *= fptot;
447 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1));
448
449 nlv.boost(fbst);
450
451 G4ThreeVector np = nlv.vect();
452 G4double theta = np.theta();
453 fThetaLab = theta;
454
455 return theta;
456}
double theta() const

The documentation for this class was generated from the following files: