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

#include <G4LEpp.hh>

Inheritance diagram for G4LEpp:

Public Member Functions

 G4LEpp ()
 ~G4LEpp () override
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) override
G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 ~G4HadronElastic () 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 41 of file G4LEpp.hh.

Constructor & Destructor Documentation

◆ G4LEpp()

G4LEpp::G4LEpp ( )
explicit

Definition at line 44 of file G4LEpp.cc.

44 :G4HadronElastic("G4LEpp")
45{
46 secID = G4PhysicsModelCatalog::GetModelID( "model_" + GetModelName() );
47 SetMinEnergy(0.);
48 SetMaxEnergy(5.*GeV);
49}
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition G4HadronElastic.cc:49
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition G4HadronicInteraction.hh:89
G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition G4HadronicInteraction.hh:115
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition G4HadronicInteraction.hh:102
G4PhysicsModelCatalog::GetModelID
static G4int GetModelID(const G4int modelIndex)
Definition G4PhysicsModelCatalog.cc:709

◆ ~G4LEpp()

G4LEpp::~G4LEpp ( )
override

Definition at line 51 of file G4LEpp.cc.

52{}

Member Function Documentation

◆ ApplyYourself()

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

Reimplemented from G4HadronElastic.

Definition at line 55 of file G4LEpp.cc.

56{
57 theParticleChange.Clear();
58 const G4HadProjectile* aParticle = &aTrack;
59
60 G4double P = aParticle->GetTotalMomentum();
61 G4double Px = aParticle->Get4Momentum().x();
62 G4double Py = aParticle->Get4Momentum().y();
63 G4double Pz = aParticle->Get4Momentum().z();
64 G4double E = aParticle->GetTotalEnergy();
65 G4ThreeVector theInitial = aParticle->Get4Momentum().vect().unit();
66
67 if (verboseLevel > 1) {
68 G4double ek = aParticle->GetKineticEnergy();
69 G4double E0 = aParticle->GetDefinition()->GetPDGMass();
70 G4double Q = aParticle->GetDefinition()->GetPDGCharge();
71 G4int A = targetNucleus.GetA_asInt();
72 G4int Z = targetNucleus.GetZ_asInt();
73 G4cout << "G4LEpp:ApplyYourself: incident particle: "
74 << aParticle->GetDefinition()->GetParticleName() << G4endl;
75 G4cout << "P = " << P/GeV << " GeV/c"
76 << ", Px = " << Px/GeV << " GeV/c"
77 << ", Py = " << Py/GeV << " GeV/c"
78 << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
79 G4cout << "E = " << E/GeV << " GeV"
80 << ", kinetic energy = " << ek/GeV << " GeV"
81 << ", mass = " << E0/GeV << " GeV"
82 << ", charge = " << Q << G4endl;
83 G4cout << "G4LEpp:ApplyYourself: material:" << G4endl;
84 G4cout << "A = " << A
85 << ", Z = " << Z
86 << ", atomic mass "
87 << G4Proton::Proton()->GetPDGMass()/GeV << "GeV"
88 << G4endl;
89 //
90 // GHEISHA ADD operation to get total energy, mass, charge
91 //
92 E += proton_mass_c2;
93 G4double E02 = E*E - P*P;
94 E0 = std::sqrt(std::fabs(E02));
95 if (E02 < 0)E0 *= -1;
96 Q += Z;
97 G4cout << "G4LEpp:ApplyYourself: total:" << G4endl;
98 G4cout << "E = " << E/GeV << " GeV"
99 << ", mass = " << E0/GeV << " GeV"
100 << ", charge = " << Q << G4endl;
101 }
102 G4double t = SampleInvariantT(aParticle->GetDefinition(), P, 0, 0);
103 G4double cost = 1.0 - 2*t/(P*P);
104 if(cost > 1.0) { cost = 1.0; }
105 if(cost <-1.0) { cost =-1.0; }
106 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
107 G4double phi = twopi*G4UniformRand();
108 // Get the target particle
109 G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
110
111 G4double E1 = aParticle->GetTotalEnergy();
112 G4double M1 = aParticle->GetDefinition()->GetPDGMass();
113 G4double E2 = targetParticle->GetTotalEnergy();
114 G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
115 G4double totalEnergy = E1 + E2;
116 G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
117
118 // Transform into centre of mass system
119
120 G4double px = (M2/pseudoMass)*Px;
121 G4double py = (M2/pseudoMass)*Py;
122 G4double pz = (M2/pseudoMass)*Pz;
123 G4double p = std::sqrt(px*px + py*py + pz*pz);
124
125 if (verboseLevel > 1) {
126 G4cout << " E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
127 G4cout << " E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
128 G4cout << " particle 1 momentum in CM " << px/GeV
129 << " " << py/GeV << " "
130 << pz/GeV << " " << p/GeV << G4endl;
131 }
132
133 // First scatter w.r.t. Z axis
134 G4double pxnew = p*sint*std::cos(phi);
135 G4double pynew = p*sint*std::sin(phi);
136 G4double pznew = p*cost;
137
138 // Rotate according to the direction of the incident particle
139 if (px*px + py*py > 0) {
140 G4double ph, cosp, sinp;
141 cost = pz/p;
142 sint = (std::sqrt((1-cost)*(1+cost)) + std::sqrt(px*px+py*py)/p)/2;
143 py < 0 ? ph = 3*halfpi : ph = halfpi;
144 if (std::fabs(px) > 0.000001*GeV) ph = std::atan2(py,px);
145 cosp = std::cos(ph);
146 sinp = std::sin(ph);
147 px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
148 py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
149 pz = (-sint*pxnew + cost*pznew);
150 }
151 else {
152 px = pxnew;
153 py = pynew;
154 pz = pznew;
155 }
156
157 if (verboseLevel > 1) {
158 G4cout << " AFTER SCATTER..." << G4endl;
159 G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
160 << pz/GeV << " " << p/GeV << G4endl;
161 }
162
163 // Transform to lab system
164
165 G4double E1pM2 = E1 + M2;
166 G4double betaCM = P/E1pM2;
167 G4double betaCMx = Px/E1pM2;
168 G4double betaCMy = Py/E1pM2;
169 G4double betaCMz = Pz/E1pM2;
170 G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
171
172 if (verboseLevel > 1) {
173 G4cout << " betaCM " << betaCMx << " " << betaCMy << " "
174 << betaCMz << " " << betaCM << G4endl;
175 G4cout << " gammaCM " << gammaCM << G4endl;
176 }
177
178 // Now following GLOREN...
179
180 G4double BETA[5], PA[5], PB[5];
181 BETA[1] = -betaCMx;
182 BETA[2] = -betaCMy;
183 BETA[3] = -betaCMz;
184 BETA[4] = gammaCM;
185
186 //The incident particle...
187
188 PA[1] = px;
189 PA[2] = py;
190 PA[3] = pz;
191 PA[4] = std::sqrt(M1*M1 + p*p);
192
193 G4double BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
194 G4double BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
195
196 PB[1] = PA[1] + BPGAM * BETA[1];
197 PB[2] = PA[2] + BPGAM * BETA[2];
198 PB[3] = PA[3] + BPGAM * BETA[3];
199 PB[4] = (PA[4] - BETPA) * BETA[4];
200
201 G4DynamicParticle* newP = new G4DynamicParticle;
202 newP->SetDefinition(aParticle->GetDefinition());
203 newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
204
205 //The target particle...
206
207 PA[1] = -px;
208 PA[2] = -py;
209 PA[3] = -pz;
210 PA[4] = std::sqrt(M2*M2 + p*p);
211
212 BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
213 BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
214
215 PB[1] = PA[1] + BPGAM * BETA[1];
216 PB[2] = PA[2] + BPGAM * BETA[2];
217 PB[3] = PA[3] + BPGAM * BETA[3];
218 PB[4] = (PA[4] - BETPA) * BETA[4];
219
220 targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
221
222 if (verboseLevel > 1) {
223 G4cout << " particle 1 momentum in LAB "
224 << newP->GetMomentum()/GeV
225 << " " << newP->GetTotalMomentum()/GeV << G4endl;
226 G4cout << " particle 2 momentum in LAB "
227 << targetParticle->GetMomentum()/GeV
228 << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
229 G4cout << " TOTAL momentum in LAB "
230 << (newP->GetMomentum()+targetParticle->GetMomentum())/GeV
231 << " "
232 << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
233 << G4endl;
234 }
235
236 theParticleChange.SetMomentumChange( newP->GetMomentumDirection());
237 theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
238 delete newP;
239
240 // Recoil particle
241 theParticleChange.AddSecondary(targetParticle, secID);
242 return &theParticleChange;
243}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition G4DynamicParticle.cc:279
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetTotalEnergy
G4double GetTotalEnergy() const
G4DynamicParticle::SetMomentum
void SetMomentum(const G4ThreeVector &momentum)
Definition G4DynamicParticle.cc:332
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4HadProjectile::GetTotalMomentum
G4double GetTotalMomentum() const
Definition G4HadProjectile.hh:111
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition G4HadProjectile.hh:86
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition G4HadProjectile.hh:116
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition G4HadProjectile.hh:91
G4HadProjectile::GetTotalEnergy
G4double GetTotalEnergy() const
Definition G4HadProjectile.hh:106
G4LEpp::SampleInvariantT
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
Definition G4LEpp.cc:249
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition G4Nucleus.hh:78
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition G4Nucleus.hh:84
G4Nucleus::ReturnTargetParticle
G4DynamicParticle * ReturnTargetParticle() const
Definition G4Nucleus.cc:352
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96
G4Proton::Proton
static G4Proton * Proton()
Definition G4Proton.cc:90

◆ SampleInvariantT()

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

Reimplemented from G4HadronElastic.

Definition at line 249 of file G4LEpp.cc.

251{
252 G4double nMass = p->GetPDGMass(); // 939.565346*MeV;
253 G4double ek = std::sqrt(plab*plab+nMass*nMass) - nMass;
254
255 // Find energy bin
256
257 G4int je1 = 0;
258 G4int je2 = NENERGY - 1;
259 ek /= GeV;
260
261 do
262 {
263 G4int midBin = (je1 + je2)/2;
264
265 if (ek < elab[midBin]) je2 = midBin;
266 else je1 = midBin;
267 }
268 while (je2 - je1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
269
270 G4double delab = elab[je2] - elab[je1];
271
272 // Sample the angle
273
274 G4double sample = G4UniformRand();
275 G4int ke1 = 0;
276 G4int ke2 = NANGLE - 1;
277 G4double dsig, b, rc;
278
279 dsig = Sig[je2][0] - Sig[je1][0];
280 rc = dsig/delab;
281 b = Sig[je1][0] - rc*elab[je1];
282
283 G4double sigint1 = rc*ek + b;
284 G4double sigint2 = 0.;
285
286 do
287 {
288 G4int midBin = (ke1 + ke2)/2;
289 dsig = Sig[je2][midBin] - Sig[je1][midBin];
290 rc = dsig/delab;
291 b = Sig[je1][midBin] - rc*elab[je1];
292 G4double sigint = rc*ek + b;
293
294 if (sample < sigint)
295 {
296 ke2 = midBin;
297 sigint2 = sigint;
298 }
299 else
300 {
301 ke1 = midBin;
302 sigint1 = sigint;
303 }
304 }
305 while (ke2 - ke1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
306
307 dsig = sigint2 - sigint1;
308 rc = 1./dsig;
309 b = ke1 - rc*sigint1;
310
311 G4double kint = rc*sample + b;
312 G4double theta = (0.5 + kint)*pi/180.;
313 G4double t = 0.5*plab*plab*(1 - std::cos(theta));
314
315 return t;
316}

Referenced by ApplyYourself().

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