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

#include <G4HadronElastic.hh>

Inheritance diagram for G4HadronElastic:

Public Member Functions

 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 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 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

Protected Attributes

G4double pLocalTmax
G4int secID
Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
G4int verboseLevel
G4double theMinEnergy
G4double theMaxEnergy
G4bool isBlocked

Additional Inherited Members

Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
G4bool IsBlocked () const
void Block ()

Detailed Description

Definition at line 48 of file G4HadronElastic.hh.

Constructor & Destructor Documentation

◆ G4HadronElastic()

G4HadronElastic::G4HadronElastic ( const G4String & name = "hElasticLHEP")
explicit

Definition at line 49 of file G4HadronElastic.cc.

50 : G4HadronicInteraction(name), secID(-1)
51{
52 SetMinEnergy( 0.0*GeV );
54 lowestEnergyLimit= 1.e-6*eV;
55 pLocalTmax = 0.0;
56 nwarn = 0;
57
58 theProton = G4Proton::Proton();
59 theNeutron = G4Neutron::Neutron();
60 theDeuteron = G4Deuteron::Deuteron();
61 theAlpha = G4Alpha::Alpha();
62
63 secID = G4PhysicsModelCatalog::GetModelID( "model_" + name );
64}
static G4Alpha * Alpha()
Definition G4Alpha.cc:83
static G4Deuteron * Deuteron()
Definition G4Deuteron.cc:90
void SetMinEnergy(G4double anEnergy)
G4HadronicInteraction(const G4String &modelName="HadronicModel")
void SetMaxEnergy(const G4double anEnergy)
static G4HadronicParameters * Instance()
static G4Neutron * Neutron()
Definition G4Neutron.cc:101
static G4int GetModelID(const G4int modelIndex)
static G4Proton * Proton()
Definition G4Proton.cc:90

Referenced by G4AntiNuclElastic::G4AntiNuclElastic(), G4ChipsElasticModel::G4ChipsElasticModel(), G4DiffuseElastic::G4DiffuseElastic(), G4DiffuseElasticV2::G4DiffuseElasticV2(), G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE(), G4hhElastic::G4hhElastic(), G4hhElastic::G4hhElastic(), G4hhElastic::G4hhElastic(), G4LEHadronProtonElastic::G4LEHadronProtonElastic(), G4LEnp::G4LEnp(), G4LEpp::G4LEpp(), G4LowEHadronElastic::G4LowEHadronElastic(), G4NeutrinoElectronNcModel::G4NeutrinoElectronNcModel(), G4NeutronElectronElModel::G4NeutronElectronElModel(), and G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic().

◆ ~G4HadronElastic()

G4HadronElastic::~G4HadronElastic ( )
override

Definition at line 66 of file G4HadronElastic.cc.

67{}

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4HadronElastic::ApplyYourself ( const G4HadProjectile & aTrack,
G4Nucleus & targetNucleus )
overridevirtual

Reimplemented from G4HadronicInteraction.

Reimplemented in G4LEHadronProtonElastic, G4LEnp, G4LEpp, G4NeutrinoElectronNcModel, and G4NeutronElectronElModel.

Definition at line 81 of file G4HadronElastic.cc.

83{
84 theParticleChange.Clear();
85
86 const G4HadProjectile* aParticle = &aTrack;
87 G4double ekin = aParticle->GetKineticEnergy();
88
89 // no scattering below the limit
90 if(ekin <= lowestEnergyLimit) {
91 theParticleChange.SetEnergyChange(ekin);
92 theParticleChange.SetMomentumChange(0.,0.,1.);
93 return &theParticleChange;
94 }
95
96 G4int A = targetNucleus.GetA_asInt();
97 G4int Z = targetNucleus.GetZ_asInt();
98
99 // Scattered particle referred to axis of incident particle
100 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
101 G4double m1 = theParticle->GetPDGMass();
102 G4double plab = std::sqrt(ekin*(ekin + 2.0*m1));
103
104 if (verboseLevel>1) {
105 G4cout << "G4HadronElastic: "
106 << aParticle->GetDefinition()->GetParticleName()
107 << " Plab(GeV/c)= " << plab/GeV
108 << " Ekin(MeV) = " << ekin/MeV
109 << " scattered off Z= " << Z
110 << " A= " << A
111 << G4endl;
112 }
113
115 G4double e1 = m1 + ekin;
116 G4LorentzVector lv(0.0,0.0,plab,e1+mass2);
117 G4ThreeVector bst = lv.boostVector();
118 G4double momentumCMS = plab*mass2/std::sqrt(m1*m1 + mass2*mass2 + 2.*mass2*e1);
119
120 pLocalTmax = 4.0*momentumCMS*momentumCMS;
121
122 // Sampling in CM system
123 G4double t = SampleInvariantT(theParticle, plab, Z, A);
124
125 if(t < 0.0 || t > pLocalTmax) {
126 // For the very rare cases where cos(theta) is greater than 1 or smaller than -1,
127 // print some debugging information via a "JustWarning" exception, and resample
128 // using the default algorithm
129#ifdef G4VERBOSE
130 if(nwarn < 2) {
132 ed << GetModelName() << " wrong sampling t= " << t << " tmax= " << pLocalTmax
133 << " for " << aParticle->GetDefinition()->GetParticleName()
134 << " ekin=" << ekin << " MeV"
135 << " off (Z,A)=(" << Z << "," << A << ") - will be resampled" << G4endl;
136 G4Exception( "G4HadronElastic::ApplyYourself", "hadEla001", JustWarning, ed);
137 ++nwarn;
138 }
139#endif
140 t = G4HadronElastic::SampleInvariantT(theParticle, plab, Z, A);
141 }
142
143 G4double phi = G4UniformRand()*CLHEP::twopi;
144 G4double cost = 1. - 2.0*t/pLocalTmax;
145
146 // if cos(theta) negative, there is a numerical problem
147 // instead of making scattering backward, make in this case
148 // no scattering
149 if (std::abs(cost) > 1.0) { cost = 1.0; }
150
151 G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
152
153 if (verboseLevel>1) {
154 G4cout << " t= " << t << " tmax(GeV^2)= " << pLocalTmax/(GeV*GeV)
155 << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost
156 << " sin(t)=" << sint << G4endl;
157 }
158 G4LorentzVector nlv1(momentumCMS*sint*std::cos(phi),
159 momentumCMS*sint*std::sin(phi),
160 momentumCMS*cost,
161 std::sqrt(momentumCMS*momentumCMS + m1*m1));
162
163 nlv1.boost(bst);
164
165 G4double eFinal = nlv1.e() - m1;
166 if (verboseLevel > 1) {
167 G4cout <<"G4HadronElastic: m= " << m1 << " Efin(MeV)= " << eFinal
168 << " 4-M Final: " << nlv1
169 << G4endl;
170 }
171
172 if(eFinal <= 0.0) {
173 theParticleChange.SetMomentumChange(0.0,0.0,1.0);
174 theParticleChange.SetEnergyChange(0.0);
175 } else {
176 theParticleChange.SetMomentumChange(nlv1.vect().unit());
177 theParticleChange.SetEnergyChange(eFinal);
178 }
179 lv -= nlv1;
180 G4double erec = std::max(lv.e() - mass2, 0.0);
181 if (verboseLevel > 1) {
182 G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
183 << " 4-mom: " << lv
184 << G4endl;
185 }
186
187 // the recoil is created if kinetic energy above the threshold
188 if(erec > GetRecoilEnergyThreshold()) {
189 G4ParticleDefinition * theDef = nullptr;
190 if(Z == 1 && A == 1) { theDef = theProton; }
191 else if (Z == 1 && A == 2) { theDef = theDeuteron; }
192 else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
193 else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
194 else if (Z == 2 && A == 4) { theDef = theAlpha; }
195 else {
196 theDef =
198 }
199 G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv.vect().unit(), erec);
200 theParticleChange.AddSecondary(aSec, secID);
201 } else {
202 theParticleChange.SetLocalEnergyDeposit(erec);
203 }
204
205 return &theParticleChange;
206}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition G4Types.hh:83
int G4int
Definition G4Types.hh:85
const G4double A[17]
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition Randomize.hh:52
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
G4double GetRecoilEnergyThreshold() const
const G4String & GetModelName() const
static G4He3 * He3()
Definition G4He3.cc:90
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetA_asInt() const
Definition G4Nucleus.hh:78
G4int GetZ_asInt() const
Definition G4Nucleus.hh:84
const G4String & GetParticleName() const
G4IonTable * GetIonTable() const
static G4ParticleTable * GetParticleTable()
static G4Triton * Triton()
Definition G4Triton.cc:90

◆ ComputeMomentumCMS()

G4double G4HadronElastic::ComputeMomentumCMS ( const G4ParticleDefinition * p,
G4double plab,
G4int Z,
G4int A )
inline

Definition at line 102 of file G4HadronElastic.hh.

104{
105 G4double m1 = p->GetPDGMass();
106 G4double m12= m1*m1;
108 return plab*mass2/std::sqrt(m12 + mass2*mass2 + 2.*mass2*std::sqrt(m12 + plab*plab));
109}

◆ GetSlopeCof()

G4double G4HadronElastic::GetSlopeCof ( const G4int pdg)

Definition at line 279 of file G4HadronElastic.cc.

280{
281 // The input parameter "pdg" should be the absolute value of the PDG code
282 // (i.e. the same value for a particle and its antiparticle).
283
284 G4double coeff = 1.0;
285
286 // heavy barions
287
288 static const G4double lBarCof1S = 0.88;
289 static const G4double lBarCof2S = 0.76;
290 static const G4double lBarCof3S = 0.64;
291 static const G4double lBarCof1C = 0.784378;
292 static const G4double lBarCofSC = 0.664378;
293 static const G4double lBarCof2SC = 0.544378;
294 static const G4double lBarCof1B = 0.740659;
295 static const G4double lBarCofSB = 0.620659;
296 static const G4double lBarCof2SB = 0.500659;
297
298 if( pdg == 3122 || pdg == 3222 || pdg == 3112 || pdg == 3212 )
299 {
300 coeff = lBarCof1S; // Lambda, Sigma+, Sigma-, Sigma0
301
302 } else if( pdg == 3322 || pdg == 3312 )
303 {
304 coeff = lBarCof2S; // Xi-, Xi0
305 }
306 else if( pdg == 3324)
307 {
308 coeff = lBarCof3S; // Omega
309 }
310 else if( pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 4112 )
311 {
312 coeff = lBarCof1C; // LambdaC+, SigmaC+, SigmaC++, SigmaC0
313 }
314 else if( pdg == 4332 )
315 {
316 coeff = lBarCof2SC; // OmegaC
317 }
318 else if( pdg == 4232 || pdg == 4132 )
319 {
320 coeff = lBarCofSC; // XiC+, XiC0
321 }
322 else if( pdg == 5122 || pdg == 5222 || pdg == 5112 || pdg == 5212 )
323 {
324 coeff = lBarCof1B; // LambdaB, SigmaB+, SigmaB-, SigmaB0
325 }
326 else if( pdg == 5332 )
327 {
328 coeff = lBarCof2SB; // OmegaB-
329 }
330 else if( pdg == 5132 || pdg == 5232 ) // XiB-, XiB0
331 {
332 coeff = lBarCofSB;
333 }
334 // heavy mesons Kaons?
335 static const G4double lMesCof1S = 0.82; // Kp/piP kaons?
336 static const G4double llMesCof1C = 0.676568;
337 static const G4double llMesCof1B = 0.610989;
338 static const G4double llMesCof2C = 0.353135;
339 static const G4double llMesCof2B = 0.221978;
340 static const G4double llMesCofSC = 0.496568;
341 static const G4double llMesCofSB = 0.430989;
342 static const G4double llMesCofCB = 0.287557;
343 static const G4double llMesCofEtaP = 0.88;
344 static const G4double llMesCofEta = 0.76;
345
346 if( pdg == 321 || pdg == 311 || pdg == 310 )
347 {
348 coeff = lMesCof1S; //K+-0
349 }
350 else if( pdg == 511 || pdg == 521 )
351 {
352 coeff = llMesCof1B; // BMeson0, BMeson+
353 }
354 else if(pdg == 421 || pdg == 411 )
355 {
356 coeff = llMesCof1C; // DMeson+, DMeson0
357 }
358 else if( pdg == 531 )
359 {
360 coeff = llMesCofSB; // BSMeson0
361 }
362 else if( pdg == 541 )
363 {
364 coeff = llMesCofCB; // BCMeson+-
365 }
366 else if(pdg == 431 )
367 {
368 coeff = llMesCofSC; // DSMeson+-
369 }
370 else if(pdg == 441 || pdg == 443 )
371 {
372 coeff = llMesCof2C; // Etac, JPsi
373 }
374 else if(pdg == 553 )
375 {
376 coeff = llMesCof2B; // Upsilon
377 }
378 else if(pdg == 221 )
379 {
380 coeff = llMesCofEta; // Eta
381 }
382 else if(pdg == 331 )
383 {
384 coeff = llMesCofEtaP; // Eta'
385 }
386 return coeff;
387}

◆ LowestEnergyLimit()

G4double G4HadronElastic::LowestEnergyLimit ( ) const
inline

Definition at line 96 of file G4HadronElastic.hh.

97{
98 return lowestEnergyLimit;
99}

Referenced by G4NeutrinoElectronNcModel::ApplyYourself(), and G4NeutronElectronElModel::ApplyYourself().

◆ ModelDescription()

void G4HadronElastic::ModelDescription ( std::ostream & outFile) const
overridevirtual

Reimplemented from G4HadronicInteraction.

Reimplemented in G4NeutrinoElectronNcModel, and G4NeutronElectronElModel.

Definition at line 70 of file G4HadronElastic.cc.

71{
72 outFile << "G4HadronElastic is the base class for all hadron-nucleus\n"
73 << "elastic scattering models except HP.\n"
74 << "By default it uses the Gheisha two-exponential momentum\n"
75 << "transfer parameterization. The model is fully relativistic\n"
76 << "as opposed to the original Gheisha model which was not.\n"
77 << "This model may be used for all long-lived hadrons at all\n"
78 << "incident energies but fit the data only for relativistic scattering.\n";
79}

◆ SampleInvariantT()

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

Reimplemented from G4HadronicInteraction.

Reimplemented in G4hhElastic, G4LEHadronProtonElastic, G4LEnp, G4LEpp, G4LowEHadronElastic, and G4NuclNuclDiffuseElastic.

Definition at line 210 of file G4HadronElastic.cc.

212{
213 const G4double plabLowLimit = 400.0*CLHEP::MeV;
214 const G4double GeV2 = CLHEP::GeV*CLHEP::GeV;
215 const G4double z07in13 = std::pow(0.7, 0.3333333333);
216 const G4double numLimit = 18.;
217
218 G4int pdg = std::abs(part->GetPDGEncoding());
219 G4double tmax = pLocalTmax/GeV2;
220
221 G4double aa, bb, cc, dd;
222 G4Pow* g4pow = G4Pow::GetInstance();
223 if (A <= 62) {
224 if (pdg == 211){ //Pions
225 if(mom >= plabLowLimit){ //High energy
226 bb = 14.5*g4pow->Z23(A);/*14.5*/
227 dd = 10.;
228 cc = 0.075*g4pow->Z13(A)/dd;//1.4
229 //aa = g4pow->powZ(A, 1.93)/bb;//1.63
230 aa = (A*A)/bb;//1.63
231 } else { //Low energy
232 bb = 29.*z07in13*z07in13*g4pow->Z23(A);
233 dd = 15.;
234 cc = 0.04*g4pow->Z13(A)/dd;//1.4
235 aa = g4pow->powZ(A, 1.63)/bb;//1.63
236 }
237 } else { //Other particles
238 bb = 14.5*g4pow->Z23(A);
239 dd = 20.;
240 aa = (A*A)/bb;//1.63
241 cc = 1.4*g4pow->Z13(A)/dd;
242 }
243 //===========================
244 } else { //(A>62)
245 if (pdg == 211) {
246 if(mom >= plabLowLimit){ //high
247 bb = 60.*z07in13*g4pow->Z13(A);//60
248 dd = 30.;
249 aa = 0.5*(A*A)/bb;//1.33
250 cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4 --- 2: 0.4
251 } else { //low
252 bb = 120.*z07in13*g4pow->Z13(A);//60
253 dd = 30.;
254 aa = 2.*g4pow->powZ(A,1.33)/bb;
255 cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4 --- 2: 0.4
256 }
257 } else {
258 bb = 60.*g4pow->Z13(A);
259 dd = 25.;
260 aa = g4pow->powZ(A,1.33)/bb;//1.33
261 cc = 0.2*g4pow->powZ(A,0.4)/dd;//1:0.4 --- 2: 0.4
262 }
263 }
264 G4double q1 = 1.0 - G4Exp(-std::min(bb*tmax, numLimit));
265 G4double q2 = 1.0 - G4Exp(-std::min(dd*tmax, numLimit));
266 G4double s1 = q1*aa;
267 G4double s2 = q2*cc;
268 if ((s1 + s2)*G4UniformRand() < s2) {
269 q1 = q2;
270 bb = dd;
271 }
272 return -GeV2*G4Log(1.0 - G4UniformRand()*q1)/bb;
273}
const G4double GeV2
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
G4double G4Log(G4double x)
Definition G4Log.hh:169
static G4Pow * GetInstance()
Definition G4Pow.cc:41
G4double powZ(G4int Z, G4double y) const
Definition G4Pow.hh:225
G4double Z13(G4int Z) const
Definition G4Pow.hh:123
G4double Z23(G4int Z) const
Definition G4Pow.hh:125

Referenced by ApplyYourself(), G4ChipsElasticModel::SampleInvariantT(), G4ElasticHadrNucleusHE::SampleInvariantT(), and G4LowEHadronElastic::SampleInvariantT().

◆ SetLowestEnergyLimit()

void G4HadronElastic::SetLowestEnergyLimit ( G4double value)
inline

Definition at line 91 of file G4HadronElastic.hh.

92{
93 lowestEnergyLimit = value;
94}

Referenced by G4NeutrinoElectronNcModel::G4NeutrinoElectronNcModel(), and G4NeutronElectronElModel::G4NeutronElectronElModel().

Member Data Documentation

◆ pLocalTmax

◆ secID


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