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

#include <G4LEHadronProtonElastic.hh>

Inheritance diagram for G4LEHadronProtonElastic:

Public Member Functions

 G4LEHadronProtonElastic ()
 ~G4LEHadronProtonElastic () override
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) override
G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
G4double RandCosThetaDipPen ()
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 40 of file G4LEHadronProtonElastic.hh.

Constructor & Destructor Documentation

◆ G4LEHadronProtonElastic()

G4LEHadronProtonElastic::G4LEHadronProtonElastic ( )

Definition at line 44 of file G4LEHadronProtonElastic.cc.

44 :
45 G4HadronElastic("G4LEHadronProtonElastic")
46{
47 secID = G4PhysicsModelCatalog::GetModelID( "model_" + GetModelName() );
48 SetMinEnergy(0.);
49 SetMaxEnergy(20.*MeV);
50}
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

◆ ~G4LEHadronProtonElastic()

G4LEHadronProtonElastic::~G4LEHadronProtonElastic ( )
override

Definition at line 52 of file G4LEHadronProtonElastic.cc.

53{
54 theParticleChange.Clear();
55}

Member Function Documentation

◆ ApplyYourself()

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

Reimplemented from G4HadronElastic.

Definition at line 58 of file G4LEHadronProtonElastic.cc.

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

◆ RandCosThetaDipPen()

G4double G4LEHadronProtonElastic::RandCosThetaDipPen ( )

Definition at line 278 of file G4LEHadronProtonElastic.cc.

279{
280 G4double x, cosTheta, signX, modX, power = 1./3.;
281
282 if( G4UniformRand() > 0.25)
283 {
284 cosTheta = 2.*G4UniformRand()-1.;
285 }
286 else
287 {
288 x = 2.*G4UniformRand()-1.;
289
290 if ( x < 0. )
291 {
292 modX = -x;
293 signX = -1.;
294 }
295 else
296 {
297 modX = x;
298 signX = 1.;
299 }
300 cosTheta = signX*G4Pow::GetInstance()->powA(modX,power);
301 }
302 return cosTheta;
303}
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition G4Pow.cc:41
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition G4Pow.hh:230

Referenced by SampleInvariantT().

◆ SampleInvariantT()

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

Reimplemented from G4HadronElastic.

Definition at line 255 of file G4LEHadronProtonElastic.cc.

257{
258 G4double hMass = p->GetPDGMass();
259 G4double pCMS = 0.5*plab;
260 // pCMS *= 50;
261 G4double hEcms = std::sqrt(pCMS*pCMS+hMass*hMass);
262 // G4double gamma = hEcms/hMass;
263 // gamma *= 15;
264 G4double beta = pCMS/hEcms; // std::sqrt(1-1./gamma/gamma); //
265 // beta /= 0.8; // 0.95; // 1.0; // 1.1 // 0.5*pi; // pi; twopi;
266 G4double cosDipole = RandCosThetaDipPen();
267 G4double cosTheta = cosDipole + beta;
268 cosTheta /= 1. + cosDipole*beta;
269 G4double t = 2.*pCMS*pCMS*(1.-cosTheta);
270 return t;
271
272}
The documentation for this class was generated from the following files: