Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EnergyLossForExtrapolator.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4EnergyLossForExtrapolator.cc 97392 2016-06-02 10:10:32Z gcosmo $
27 //
28 //---------------------------------------------------------------------------
29 //
30 // ClassName: G4EnergyLossForExtrapolator
31 //
32 // Description: This class provide calculation of energy loss, fluctuation,
33 // and msc angle
34 //
35 // Author: 09.12.04 V.Ivanchenko
36 //
37 // Modification:
38 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
39 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
40 // 21-03-06 Add verbosity defined in the constructor and Initialisation
41 // start only when first public method is called (V.Ivanchenko)
42 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
43 // 12-05-06 SEt linLossLimit=0.001 (VI)
44 //
45 //----------------------------------------------------------------------------
46 //
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 
51 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4Material.hh"
55 #include "G4MaterialCutsCouple.hh"
56 #include "G4Electron.hh"
57 #include "G4Positron.hh"
58 #include "G4Proton.hh"
59 #include "G4MuonPlus.hh"
60 #include "G4MuonMinus.hh"
61 #include "G4ParticleTable.hh"
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
65 G4TablesForExtrapolator* G4EnergyLossForExtrapolator::tables = nullptr;
66 
68  : maxEnergyTransfer(DBL_MAX), verbose(verb)
69 {
70  currentParticle = nullptr;
71  currentMaterial = nullptr;
72 
73  linLossLimit = 0.01;
74  emin = 1.*MeV;
75  emax = 10.*TeV;
76  nbins = 70;
77 
78  nmat = index = 0;
79 
80  mass = charge2 = electronDensity = radLength = bg2 = beta2
81  = kineticEnergy = tmax = 0.0;
82  gam = 1.0;
83 
84  idxDedxElectron = idxDedxPositron = idxDedxMuon = idxDedxProton
85  = idxRangeElectron = idxRangePositron = idxRangeMuon = idxRangeProton
86  = idxInvRangeElectron = idxInvRangePositron = idxInvRangeMuon
87  = idxInvRangeProton = idxMscElectron = 0;
88 
89  electron = positron = proton = muonPlus = muonMinus = nullptr;
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 
94 G4double
96  G4double stepLength,
97  const G4Material* mat,
98  const G4ParticleDefinition* part)
99 {
100  if(0 == nmat) { Initialisation(); }
101  G4double kinEnergyFinal = kinEnergy;
102  if(SetupKinematics(part, mat, kinEnergy)) {
103  G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
104  G4double r = ComputeRange(kinEnergy,part);
105  if(r <= step) {
106  kinEnergyFinal = 0.0;
107  } else if(step < linLossLimit*r) {
108  kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
109  } else {
110  G4double r1 = r - step;
111  kinEnergyFinal = ComputeEnergy(r1,part);
112  }
113  }
114  return kinEnergyFinal;
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 
119 G4double
121  G4double stepLength,
122  const G4Material* mat,
123  const G4ParticleDefinition* part)
124 {
125  //G4cout << "G4EnergyLossForExtrapolator::EnergyBeforeStep" << G4endl;
126  if(0 == nmat) { Initialisation(); }
127  G4double kinEnergyFinal = kinEnergy;
128 
129  if(SetupKinematics(part, mat, kinEnergy)) {
130  G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
131  G4double r = ComputeRange(kinEnergy,part);
132 
133  if(step < linLossLimit*r) {
134  kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
135  } else {
136  G4double r1 = r + step;
137  kinEnergyFinal = ComputeEnergy(r1,part);
138  }
139  }
140  return kinEnergyFinal;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144 
145 G4double
147  G4double stepLength,
148  const G4Material* mat,
149  const G4ParticleDefinition* part)
150 {
151  G4double res = stepLength;
152  if(0 == nmat) { Initialisation(); }
153  //G4cout << "## G4EnergyLossForExtrapolator::TrueStepLength L= " << res
154  // << " " << part->GetParticleName() << G4endl;
155  if(SetupKinematics(part, mat, kinEnergy)) {
156  if(part == electron || part == positron) {
157  G4double x = stepLength*
158  ComputeValue(kinEnergy, GetPhysicsTable(fMscElectron), idxMscElectron);
159  //G4cout << " x= " << x << G4endl;
160  if(x < 0.2) { res *= (1.0 + 0.5*x + x*x/3.0); }
161  else if(x < 0.9999) { res = -G4Log(1.0 - x)*stepLength/x; }
162  else { res = ComputeRange(kinEnergy,part); }
163  } else {
164  res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
165  }
166  }
167  return res;
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 
172 G4bool
173 G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
174  const G4Material* mat,
175  G4double kinEnergy)
176 {
177  if(0 == nmat) { Initialisation(); }
178  if(!part || !mat || kinEnergy < keV) { return false; }
179  G4bool flag = false;
180  if(part != currentParticle) {
181  flag = true;
182  currentParticle = part;
183  mass = part->GetPDGMass();
184  G4double q = part->GetPDGCharge()/eplus;
185  charge2 = q*q;
186  }
187  if(mat != currentMaterial) {
188  G4int i = mat->GetIndex();
189  if(i >= nmat) {
190  G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= "
191  << i << " is out of table - NO extrapolation" << G4endl;
192  } else {
193  flag = true;
194  currentMaterial = mat;
195  electronDensity = mat->GetElectronDensity();
196  radLength = mat->GetRadlen();
197  index = i;
198  }
199  }
200  if(flag || kinEnergy != kineticEnergy) {
201  kineticEnergy = kinEnergy;
202  G4double tau = kinEnergy/mass;
203 
204  gam = tau + 1.0;
205  bg2 = tau * (tau + 2.0);
206  beta2 = bg2/(gam*gam);
207  tmax = kinEnergy;
208  if(part == electron) tmax *= 0.5;
209  else if(part != positron) {
210  G4double r = electron_mass_c2/mass;
211  tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
212  }
213  if(tmax > maxEnergyTransfer) { tmax = maxEnergyTransfer; }
214  }
215  return true;
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
219 
220 const G4ParticleDefinition*
221 G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
222 {
223  const G4ParticleDefinition* p = nullptr;
224  if(name != currentParticleName) {
226  if(!p) {
227  G4cout << "### G4EnergyLossForExtrapolator WARNING: "
228  << "FindParticle() fails to find "
229  << name << G4endl;
230  }
231  } else {
232  p = currentParticle;
233  }
234  return p;
235 }
236 
237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
238 
239 G4double
241  const G4ParticleDefinition* part)
242 {
243  G4double x = 0.0;
244  if(part == electron) {
245  x = ComputeValue(ekin, GetPhysicsTable(fDedxElectron), idxDedxElectron);
246  } else if(part == positron) {
247  x = ComputeValue(ekin, GetPhysicsTable(fDedxPositron), idxDedxPositron);
248  } else if(part == muonPlus || part == muonMinus) {
249  x = ComputeValue(ekin, GetPhysicsTable(fDedxMuon), idxDedxMuon);
250  } else {
251  G4double e = ekin*proton_mass_c2/mass;
252  x = ComputeValue(e, GetPhysicsTable(fDedxProton), idxDedxProton)*charge2;
253  }
254  return x;
255 }
256 
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
258 
259 G4double
261  const G4ParticleDefinition* part)
262 {
263  G4double x = 0.0;
264  if(part == electron) {
265  x = ComputeValue(ekin, GetPhysicsTable(fRangeElectron), idxRangeElectron);
266  } else if(part == positron) {
267  x = ComputeValue(ekin, GetPhysicsTable(fRangePositron), idxRangePositron);
268  } else if(part == muonPlus || part == muonMinus) {
269  x = ComputeValue(ekin, GetPhysicsTable(fRangeMuon), idxRangeMuon);
270  } else {
271  G4double massratio = proton_mass_c2/mass;
272  G4double e = ekin*massratio;
273  x = ComputeValue(e, GetPhysicsTable(fRangeProton), idxRangeProton)
274  /(charge2*massratio);
275  }
276  return x;
277 }
278 
279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
280 
281 G4double
283  const G4ParticleDefinition* part)
284 {
285  G4double x = 0.0;
286  if(part == electron) {
287  x = ComputeValue(range, GetPhysicsTable(fInvRangeElectron),
288  idxInvRangeElectron);
289  } else if(part == positron) {
290  x = ComputeValue(range, GetPhysicsTable(fInvRangePositron),
291  idxInvRangePositron);
292  } else if(part == muonPlus || part == muonMinus) {
293  x = ComputeValue(range, GetPhysicsTable(fInvRangeMuon), idxInvRangeMuon);
294  } else {
295  G4double massratio = proton_mass_c2/mass;
296  G4double r = range*massratio*charge2;
297  x = ComputeValue(r, GetPhysicsTable(fInvRangeProton),
298  idxInvRangeProton)/massratio;
299  }
300  return x;
301 }
302 
303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
304 
305 void G4EnergyLossForExtrapolator::Initialisation()
306 {
307  if(verbose>1) {
308  G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
309  }
310  currentParticle = nullptr;
311  currentMaterial = nullptr;
312  kineticEnergy = 0.0;
313 
314  electron = G4Electron::Electron();
315  positron = G4Positron::Positron();
316  proton = G4Proton::Proton();
317  muonPlus = G4MuonPlus::MuonPlus();
318  muonMinus= G4MuonMinus::MuonMinus();
319 
321  currentParticleName = "";
322  BuildTables();
323 }
324 
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
326 
327 #include "G4AutoLock.hh"
328 namespace { G4Mutex G4EnergyLossForExtrapolatorMutex = G4MUTEX_INITIALIZER; }
329 
330 void G4EnergyLossForExtrapolator::BuildTables()
331 {
332  G4AutoLock l(&G4EnergyLossForExtrapolatorMutex);
333 
334  if(!tables) {
335  if(verbose > 0) {
336  G4cout << "### G4EnergyLossForExtrapolator::BuildTables for "
337  << nmat << " materials Nbins= " << nbins
338  << " Emin(MeV)= " << emin << " Emax(MeV)= " << emax << G4endl;
339  }
340  tables = new G4TablesForExtrapolator(verbose, nbins, emin, emax);
341  }
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 
347 {
348  G4AutoLock l(&G4EnergyLossForExtrapolatorMutex);
349  delete tables;
350  tables = nullptr;
351 }
352 
353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
354 
const XML_Char * name
Definition: expat.h:151
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
size_t GetIndex() const
Definition: G4Material.hh:262
const char * p
Definition: xmltok.h:285
tuple x
Definition: test.py:50
G4double TrueStepLength(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *part)
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
static constexpr double TeV
Definition: G4SIunits.hh:218
G4double ComputeTrueStep(const G4Material *, const G4ParticleDefinition *part, G4double kinEnergy, G4double stepLength)
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition const G4Material *G4double range
G4double GetElectronDensity() const
Definition: G4Material.hh:217
bool G4bool
Definition: G4Types.hh:79
static constexpr double eplus
Definition: G4SIunits.hh:199
static G4Proton * Proton()
Definition: G4Proton.cc:93
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
G4double ComputeDEDX(G4double kinEnergy, const G4ParticleDefinition *)
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
G4double ComputeRange(G4double kinEnergy, const G4ParticleDefinition *)
G4double GetRadlen() const
Definition: G4Material.hh:220
G4double EnergyBeforeStep(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *)
G4double EnergyAfterStep(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4int G4Mutex
Definition: G4Threading.hh:173
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4double ComputeEnergy(G4double range, const G4ParticleDefinition *)
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
static constexpr double keV
Definition: G4SIunits.hh:216