Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
G4LEnp Class Reference

#include <G4LEnp.hh>

Inheritance diagram for G4LEnp:
Inheritance graph
[legend]
Collaboration diagram for G4LEnp:
Collaboration graph
[legend]

Public Member Functions

 G4LEnp ()
 
 ~G4LEnp ()
 
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
- Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 
virtual ~G4HadronElastic ()
 
void SetLowestEnergyLimit (G4double value)
 
G4double LowestEnergyLimit () const
 
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual void ModelDescription (std::ostream &) const
 
- 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 ()
 

Additional Inherited Members

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

Detailed Description

Definition at line 57 of file G4LEnp.hh.

Constructor & Destructor Documentation

G4LEnp::G4LEnp ( )

Definition at line 45 of file G4LEnp.cc.

45  :
46  G4HadronElastic("G4LEnp") // G4HadronicInteraction("G4LEnp")
47 {
48  // theParticleChange.SetNumberOfSecondaries(1);
49 
50  // SetMinEnergy(10.*MeV);
51  // SetMaxEnergy(1200.*MeV);
52  SetMinEnergy(0.);
53  SetMaxEnergy(5.*GeV);
54 }
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition: G4HadronElastic.cc:48
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103

Here is the call graph for this function:

G4LEnp::~G4LEnp ( )

Definition at line 56 of file G4LEnp.cc.

57 {
58  theParticleChange.Clear();
59 }

Here is the call graph for this function:

Member Function Documentation

G4HadFinalState * G4LEnp::ApplyYourself ( const G4HadProjectile &  aTrack,
G4Nucleus &  targetNucleus 
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 62 of file G4LEnp.cc.

63 {
64  theParticleChange.Clear();
65  const G4HadProjectile* aParticle = &aTrack;
66 
67  G4double P = aParticle->GetTotalMomentum();
68  G4double Px = aParticle->Get4Momentum().x();
69  G4double Py = aParticle->Get4Momentum().y();
70  G4double Pz = aParticle->Get4Momentum().z();
71  G4double ek = aParticle->GetKineticEnergy();
72  G4ThreeVector theInitial = aParticle->Get4Momentum().vect();
73 
74  if (verboseLevel > 1) {
75  G4double E = aParticle->GetTotalEnergy();
76  G4double E0 = aParticle->GetDefinition()->GetPDGMass();
77  G4double Q = aParticle->GetDefinition()->GetPDGCharge();
78  G4int A = targetNucleus.GetA_asInt();
79  G4int Z = targetNucleus.GetZ_asInt();
80  G4cout << "G4LEnp:ApplyYourself: incident particle: "
81  << aParticle->GetDefinition()->GetParticleName() << G4endl;
82  G4cout << "P = " << P/GeV << " GeV/c"
83  << ", Px = " << Px/GeV << " GeV/c"
84  << ", Py = " << Py/GeV << " GeV/c"
85  << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
86  G4cout << "E = " << E/GeV << " GeV"
87  << ", kinetic energy = " << ek/GeV << " GeV"
88  << ", mass = " << E0/GeV << " GeV"
89  << ", charge = " << Q << G4endl;
90  G4cout << "G4LEnp:ApplyYourself: material:" << G4endl;
91  G4cout << "A = " << A
92  << ", Z = " << Z
93  << ", atomic mass "
94  << G4Proton::Proton()->GetPDGMass()/GeV << "GeV"
95  << G4endl;
96  //
97  // GHEISHA ADD operation to get total energy, mass, charge
98  //
99  E += proton_mass_c2;
100  G4double E02 = E*E - P*P;
101  E0 = std::sqrt(std::abs(E02));
102  if (E02 < 0)E0 *= -1;
103  Q += Z;
104  G4cout << "G4LEnp:ApplyYourself: total:" << G4endl;
105  G4cout << "E = " << E/GeV << " GeV"
106  << ", mass = " << E0/GeV << " GeV"
107  << ", charge = " << Q << G4endl;
108  }
109 
110  // Find energy bin
111 
112  G4int je1 = 0;
113  G4int je2 = NENERGY - 1;
114  ek = ek/GeV;
115  do {
116  G4int midBin = (je1 + je2)/2;
117  if (ek < elab[midBin])
118  je2 = midBin;
119  else
120  je1 = midBin;
121  } while (je2 - je1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
122  G4double delab = elab[je2] - elab[je1];
123 
124  // Sample the angle
125 
126  G4double sample = G4UniformRand();
127  G4int ke1 = 0;
128  G4int ke2 = NANGLE - 1;
129  G4double dsig = sig[je2][0] - sig[je1][0];
130  G4double rc = dsig/delab;
131  G4double b = sig[je1][0] - rc*elab[je1];
132  G4double sigint1 = rc*ek + b;
133  G4double sigint2 = 0.;
134 
135  if (verboseLevel > 1) {
136  G4cout << "sample=" << sample << G4endl
137  << ke1 << " " << ke2 << " "
138  << sigint1 << " " << sigint2 << G4endl;
139  }
140  do {
141  G4int midBin = (ke1 + ke2)/2;
142  dsig = sig[je2][midBin] - sig[je1][midBin];
143  rc = dsig/delab;
144  b = sig[je1][midBin] - rc*elab[je1];
145  G4double sigint = rc*ek + b;
146  if (sample < sigint) {
147  ke2 = midBin;
148  sigint2 = sigint;
149  }
150  else {
151  ke1 = midBin;
152  sigint1 = sigint;
153  }
154  if (verboseLevel > 1) {
155  G4cout << ke1 << " " << ke2 << " "
156  << sigint1 << " " << sigint2 << G4endl;
157  }
158  } while (ke2 - ke1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
159 
160  dsig = sigint2 - sigint1;
161  rc = 1./dsig;
162  b = ke1 - rc*sigint1;
163  G4double kint = rc*sample + b;
164  G4double theta = (0.5 + kint)*pi/180.;
165 
166  if (verboseLevel > 1) {
167  G4cout << " energy bin " << je1 << " energy=" << elab[je1] << G4endl;
168  G4cout << " angle bin " << kint << " angle=" << theta/degree << G4endl;
169  }
170 
171  // Get the target particle
172 
173  G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
174 
175  G4double E1 = aParticle->GetTotalEnergy();
176  G4double M1 = aParticle->GetDefinition()->GetPDGMass();
177  G4double E2 = targetParticle->GetTotalEnergy();
178  G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
179  G4double totalEnergy = E1 + E2;
180  G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
181 
182  // Transform into centre of mass system
183 
184  G4double px = (M2/pseudoMass)*Px;
185  G4double py = (M2/pseudoMass)*Py;
186  G4double pz = (M2/pseudoMass)*Pz;
187  G4double p = std::sqrt(px*px + py*py + pz*pz);
188 
189  if (verboseLevel > 1) {
190  G4cout << " E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
191  G4cout << " E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
192  G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
193  << pz/GeV << " " << p/GeV << G4endl;
194  }
195 
196  // First scatter w.r.t. Z axis
197  G4double phi = G4UniformRand()*twopi;
198  G4double pxnew = p*std::sin(theta)*std::cos(phi);
199  G4double pynew = p*std::sin(theta)*std::sin(phi);
200  G4double pznew = p*std::cos(theta);
201 
202  // Rotate according to the direction of the incident particle
203  if (px*px + py*py > 0) {
204  G4double cost, sint, ph, cosp, sinp;
205  cost = pz/p;
206  sint = (std::sqrt(std::fabs((1-cost)*(1+cost))) + std::sqrt(px*px+py*py)/p)/2;
207  py < 0 ? ph = 3*halfpi : ph = halfpi;
208  if (std::abs(px) > 0.000001*GeV) ph = std::atan2(py,px);
209  cosp = std::cos(ph);
210  sinp = std::sin(ph);
211  px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
212  py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
213  pz = (-sint*pxnew + cost*pznew);
214  }
215  else {
216  px = pxnew;
217  py = pynew;
218  pz = pznew;
219  }
220 
221  if (verboseLevel > 1) {
222  G4cout << " AFTER SCATTER..." << G4endl;
223  G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
224  << pz/GeV << " " << p/GeV << G4endl;
225  }
226 
227  // Transform to lab system
228 
229  G4double E1pM2 = E1 + M2;
230  G4double betaCM = P/E1pM2;
231  G4double betaCMx = Px/E1pM2;
232  G4double betaCMy = Py/E1pM2;
233  G4double betaCMz = Pz/E1pM2;
234  G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
235 
236  if (verboseLevel > 1) {
237  G4cout << " betaCM " << betaCMx << " " << betaCMy << " "
238  << betaCMz << " " << betaCM << G4endl;
239  G4cout << " gammaCM " << gammaCM << G4endl;
240  }
241 
242  // Now following GLOREN...
243 
244  G4double BETA[5], PA[5], PB[5];
245  BETA[1] = -betaCMx;
246  BETA[2] = -betaCMy;
247  BETA[3] = -betaCMz;
248  BETA[4] = gammaCM;
249 
250  //The incident particle...
251 
252  PA[1] = px;
253  PA[2] = py;
254  PA[3] = pz;
255  PA[4] = std::sqrt(M1*M1 + p*p);
256 
257  G4double BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
258  G4double BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
259 
260  PB[1] = PA[1] + BPGAM * BETA[1];
261  PB[2] = PA[2] + BPGAM * BETA[2];
262  PB[3] = PA[3] + BPGAM * BETA[3];
263  PB[4] = (PA[4] - BETPA) * BETA[4];
264 
265  G4DynamicParticle* newP = new G4DynamicParticle;
266  newP->SetDefinition(aParticle->GetDefinition());
267  newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
268 
269  //The target particle...
270 
271  PA[1] = -px;
272  PA[2] = -py;
273  PA[3] = -pz;
274  PA[4] = std::sqrt(M2*M2 + p*p);
275 
276  BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
277  BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
278 
279  PB[1] = PA[1] + BPGAM * BETA[1];
280  PB[2] = PA[2] + BPGAM * BETA[2];
281  PB[3] = PA[3] + BPGAM * BETA[3];
282  PB[4] = (PA[4] - BETPA) * BETA[4];
283 
284  targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
285 
286  if (verboseLevel > 1) {
287  G4cout << " particle 1 momentum in LAB "
288  << newP->GetMomentum()*(1./GeV)
289  << " " << newP->GetTotalMomentum()/GeV << G4endl;
290  G4cout << " particle 2 momentum in LAB "
291  << targetParticle->GetMomentum()*(1./GeV)
292  << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
293  G4cout << " TOTAL momentum in LAB "
294  << (newP->GetMomentum()+targetParticle->GetMomentum())*(1./GeV)
295  << " "
296  << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
297  << G4endl;
298  }
299 
300  theParticleChange.SetMomentumChange(newP->GetMomentumDirection());
301  theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
302  delete newP;
303  theParticleChange.AddSecondary(targetParticle);
304 
305  return &theParticleChange;
306 }
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4DynamicParticle::SetMomentum
void SetMomentum(const G4ThreeVector &momentum)
Definition: G4DynamicParticle.cc:371
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4DynamicParticle::GetTotalEnergy
G4double GetTotalEnergy() const
p
const char * p
Definition: xmltok.h:285
Q
static double Q[]
Definition: nf_gammaFunctions.cc:91
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
G4Nucleus::ReturnTargetParticle
G4DynamicParticle * ReturnTargetParticle() const
Definition: G4Nucleus.cc:241
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
P
static double P[]
Definition: nf_gammaFunctions.cc:89
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
test.b
tuple b
Definition: test.py:12
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
degree
static constexpr double degree
Definition: G4SIunits.hh:144
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
ek
G4double ek
Definition: G4ParticleHPJENDLHEData.cc:259
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:86
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:39
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4HadFinalState::AddSecondary
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
Definition: G4HadFinalState.hh:61
halfpi
static constexpr double halfpi
Definition: G4SIunits.hh:77
G4double
double G4double
Definition: G4Types.hh:76
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:305
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.hh:56
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4HadProjectile::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4HadProjectile.hh:101
G4HadProjectile::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4HadProjectile.hh:96

Here is the call graph for this function:

G4double G4LEnp::SampleInvariantT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 312 of file G4LEnp.cc.

314 {
315  G4double nMass = p->GetPDGMass(); // 939.565346*MeV;
316  G4double ek = std::sqrt(plab*plab+nMass*nMass) - nMass;
317 
318  // Find energy bin
319 
320  G4int je1 = 0;
321  G4int je2 = NENERGY - 1;
322  ek = ek/GeV;
323 
324  do
325  {
326  G4int midBin = (je1 + je2)/2;
327  if (ek < elab[midBin])
328  je2 = midBin;
329  else
330  je1 = midBin;
331  } while (je2 - je1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
332 
333  G4double delab = elab[je2] - elab[je1];
334 
335  // Sample the angle
336 
337  G4double sample = G4UniformRand();
338  G4int ke1 = 0;
339  G4int ke2 = NANGLE - 1;
340  G4double dsig = sig[je2][0] - sig[je1][0];
341  G4double rc = dsig/delab;
342  G4double b = sig[je1][0] - rc*elab[je1];
343  G4double sigint1 = rc*ek + b;
344  G4double sigint2 = 0.;
345 
346  do
347  {
348  G4int midBin = (ke1 + ke2)/2;
349  dsig = sig[je2][midBin] - sig[je1][midBin];
350  rc = dsig/delab;
351  b = sig[je1][midBin] - rc*elab[je1];
352  G4double sigint = rc*ek + b;
353 
354  if (sample < sigint)
355  {
356  ke2 = midBin;
357  sigint2 = sigint;
358  }
359  else
360  {
361  ke1 = midBin;
362  sigint1 = sigint;
363  }
364  } while (ke2 - ke1 > 1); /* Loop checking, 10.08.2015, A.Ribon */
365 
366  dsig = sigint2 - sigint1;
367  rc = 1./dsig;
368  b = ke1 - rc*sigint1;
369 
370  G4double kint = rc*sample + b;
371  G4double theta = (0.5 + kint)*pi/180.;
372  G4double t = 0.5*plab*plab*(1-std::cos(theta));
373 
374  return t;
375 }
G4int
int G4int
Definition: G4Types.hh:78
test.b
tuple b
Definition: test.py:12
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
ek
G4double ek
Definition: G4ParticleHPJENDLHEData.cc:259
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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