ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EnergyLossForExtrapolator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EnergyLossForExtrapolator.cc
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 //
27 //---------------------------------------------------------------------------
28 //
29 // ClassName: G4EnergyLossForExtrapolator
30 //
31 // Description: This class provide calculation of energy loss, fluctuation,
32 // and msc angle
33 //
34 // Author: 09.12.04 V.Ivanchenko
35 //
36 // Modification:
37 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
38 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
39 // 21-03-06 Add verbosity defined in the constructor and Initialisation
40 // start only when first public method is called (V.Ivanchenko)
41 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
42 // 12-05-06 SEt linLossLimit=0.001 (VI)
43 //
44 //----------------------------------------------------------------------------
45 //
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 #include "G4EnergyLossForExtrapolator.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4ParticleDefinition.hh"
53 #include "G4Material.hh"
54 #include "G4MaterialCutsCouple.hh"
55 #include "G4Electron.hh"
56 #include "G4Positron.hh"
57 #include "G4Proton.hh"
58 #include "G4MuonPlus.hh"
59 #include "G4MuonMinus.hh"
60 #include "G4ParticleTable.hh"
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
64 #ifdef G4MULTITHREADED
65 G4Mutex G4EnergyLossForExtrapolator::extrapolatorMutex = G4MUTEX_INITIALIZER;
66 #endif
67 
68 G4TablesForExtrapolator* G4EnergyLossForExtrapolator::tables = nullptr;
69 
70 G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
71  : maxEnergyTransfer(DBL_MAX), verbose(verb)
72 {
73  currentParticle = nullptr;
74  currentMaterial = nullptr;
75 
76  linLossLimit = 0.01;
77  emin = 1.*MeV;
78  emax = 10.*TeV;
79  nbins = 70;
80 
81  nmat = index = 0;
82 
83  mass = charge2 = electronDensity = radLength = bg2 = beta2
84  = kineticEnergy = tmax = 0.0;
85  gam = 1.0;
86 
87  idxDedxElectron = idxDedxPositron = idxDedxMuon = idxDedxProton
88  = idxRangeElectron = idxRangePositron = idxRangeMuon = idxRangeProton
89  = idxInvRangeElectron = idxInvRangePositron = idxInvRangeMuon
90  = idxInvRangeProton = idxMscElectron = 0;
91 
92  electron = positron = proton = muonPlus = muonMinus = nullptr;
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 
97 G4EnergyLossForExtrapolator::~G4EnergyLossForExtrapolator()
98 {
99  if(tables) {
100 #ifdef G4MULTITHREADED
101  G4MUTEXLOCK(&extrapolatorMutex);
102  if (tables) {
103 #endif
104  delete tables;
105  tables = nullptr;
106 #ifdef G4MULTITHREADED
107  }
108  G4MUTEXUNLOCK(&extrapolatorMutex);
109 #endif
110  }
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114 
115 G4double
116 G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
117  G4double stepLength,
118  const G4Material* mat,
119  const G4ParticleDefinition* part)
120 {
121  if(0 == nmat) { Initialisation(); }
122  G4double kinEnergyFinal = kinEnergy;
123  if(SetupKinematics(part, mat, kinEnergy)) {
124  G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
125  G4double r = ComputeRange(kinEnergy,part);
126  if(r <= step) {
127  kinEnergyFinal = 0.0;
128  } else if(step < linLossLimit*r) {
129  kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
130  } else {
131  G4double r1 = r - step;
132  kinEnergyFinal = ComputeEnergy(r1,part);
133  }
134  }
135  return kinEnergyFinal;
136 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 
140 G4double
141 G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
142  G4double stepLength,
143  const G4Material* mat,
144  const G4ParticleDefinition* part)
145 {
146  //G4cout << "G4EnergyLossForExtrapolator::EnergyBeforeStep" << G4endl;
147  if(0 == nmat) { Initialisation(); }
148  G4double kinEnergyFinal = kinEnergy;
149 
150  if(SetupKinematics(part, mat, kinEnergy)) {
151  G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
152  G4double r = ComputeRange(kinEnergy,part);
153 
154  if(step < linLossLimit*r) {
155  kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
156  } else {
157  G4double r1 = r + step;
158  kinEnergyFinal = ComputeEnergy(r1,part);
159  }
160  }
161  return kinEnergyFinal;
162 }
163 
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
165 
166 G4double
167 G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
168  G4double stepLength,
169  const G4Material* mat,
170  const G4ParticleDefinition* part)
171 {
172  G4double res = stepLength;
173  if(0 == nmat) { Initialisation(); }
174  //G4cout << "## G4EnergyLossForExtrapolator::TrueStepLength L= " << res
175  // << " " << part->GetParticleName() << G4endl;
176  if(SetupKinematics(part, mat, kinEnergy)) {
177  if(part == electron || part == positron) {
178  G4double x = stepLength*
179  ComputeValue(kinEnergy, GetPhysicsTable(fMscElectron), idxMscElectron);
180  //G4cout << " x= " << x << G4endl;
181  if(x < 0.2) { res *= (1.0 + 0.5*x + x*x/3.0); }
182  else if(x < 0.9999) { res = -G4Log(1.0 - x)*stepLength/x; }
183  else { res = ComputeRange(kinEnergy,part); }
184  } else {
185  res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
186  }
187  }
188  return res;
189 }
190 
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
192 
193 G4bool
194 G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
195  const G4Material* mat,
196  G4double kinEnergy)
197 {
198  if(0 == nmat) { Initialisation(); }
199  if(!part || !mat || kinEnergy < keV) { return false; }
200  G4bool flag = false;
201  if(part != currentParticle) {
202  flag = true;
203  currentParticle = part;
204  mass = part->GetPDGMass();
205  G4double q = part->GetPDGCharge()/eplus;
206  charge2 = q*q;
207  }
208  if(mat != currentMaterial) {
209  G4int i = mat->GetIndex();
210  if(i >= nmat) {
211  G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= "
212  << i << " is out of table - NO extrapolation" << G4endl;
213  } else {
214  flag = true;
215  currentMaterial = mat;
216  electronDensity = mat->GetElectronDensity();
217  radLength = mat->GetRadlen();
218  index = i;
219  }
220  }
221  if(flag || kinEnergy != kineticEnergy) {
222  kineticEnergy = kinEnergy;
223  G4double tau = kinEnergy/mass;
224 
225  gam = tau + 1.0;
226  bg2 = tau * (tau + 2.0);
227  beta2 = bg2/(gam*gam);
228  tmax = kinEnergy;
229  if(part == electron) tmax *= 0.5;
230  else if(part != positron) {
231  G4double r = electron_mass_c2/mass;
232  tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
233  }
234  if(tmax > maxEnergyTransfer) { tmax = maxEnergyTransfer; }
235  }
236  return true;
237 }
238 
239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
240 
241 const G4ParticleDefinition*
242 G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
243 {
244  const G4ParticleDefinition* p = nullptr;
245  if(name != currentParticleName) {
246  p = G4ParticleTable::GetParticleTable()->FindParticle(name);
247  if(!p) {
248  G4cout << "### G4EnergyLossForExtrapolator WARNING: "
249  << "FindParticle() fails to find "
250  << name << G4endl;
251  }
252  } else {
253  p = currentParticle;
254  }
255  return p;
256 }
257 
258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
259 
260 G4double
261 G4EnergyLossForExtrapolator::ComputeDEDX(G4double ekin,
262  const G4ParticleDefinition* part)
263 {
264  G4double x = 0.0;
265  if(part == electron) {
266  x = ComputeValue(ekin, GetPhysicsTable(fDedxElectron), idxDedxElectron);
267  } else if(part == positron) {
268  x = ComputeValue(ekin, GetPhysicsTable(fDedxPositron), idxDedxPositron);
269  } else if(part == muonPlus || part == muonMinus) {
270  x = ComputeValue(ekin, GetPhysicsTable(fDedxMuon), idxDedxMuon);
271  } else {
272  G4double e = ekin*proton_mass_c2/mass;
273  x = ComputeValue(e, GetPhysicsTable(fDedxProton), idxDedxProton)*charge2;
274  }
275  return x;
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
280 G4double
281 G4EnergyLossForExtrapolator::ComputeRange(G4double ekin,
282  const G4ParticleDefinition* part)
283 {
284  G4double x = 0.0;
285  if(part == electron) {
286  x = ComputeValue(ekin, GetPhysicsTable(fRangeElectron), idxRangeElectron);
287  } else if(part == positron) {
288  x = ComputeValue(ekin, GetPhysicsTable(fRangePositron), idxRangePositron);
289  } else if(part == muonPlus || part == muonMinus) {
290  x = ComputeValue(ekin, GetPhysicsTable(fRangeMuon), idxRangeMuon);
291  } else {
292  G4double massratio = proton_mass_c2/mass;
293  G4double e = ekin*massratio;
294  x = ComputeValue(e, GetPhysicsTable(fRangeProton), idxRangeProton)
295  /(charge2*massratio);
296  }
297  return x;
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
301 
302 G4double
303 G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
304  const G4ParticleDefinition* part)
305 {
306  G4double x = 0.0;
307  if(part == electron) {
308  x = ComputeValue(range, GetPhysicsTable(fInvRangeElectron),
309  idxInvRangeElectron);
310  } else if(part == positron) {
311  x = ComputeValue(range, GetPhysicsTable(fInvRangePositron),
312  idxInvRangePositron);
313  } else if(part == muonPlus || part == muonMinus) {
314  x = ComputeValue(range, GetPhysicsTable(fInvRangeMuon), idxInvRangeMuon);
315  } else {
316  G4double massratio = proton_mass_c2/mass;
317  G4double r = range*massratio*charge2;
318  x = ComputeValue(r, GetPhysicsTable(fInvRangeProton),
319  idxInvRangeProton)/massratio;
320  }
321  return x;
322 }
323 
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325 
326 void G4EnergyLossForExtrapolator::Initialisation()
327 {
328  if(verbose>1) {
329  G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
330  }
331  currentParticle = nullptr;
332  currentMaterial = nullptr;
333  kineticEnergy = 0.0;
334 
335  electron = G4Electron::Electron();
336  positron = G4Positron::Positron();
337  proton = G4Proton::Proton();
338  muonPlus = G4MuonPlus::MuonPlus();
339  muonMinus= G4MuonMinus::MuonMinus();
340 
341  currentParticleName = "";
342  BuildTables();
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346 
347 void G4EnergyLossForExtrapolator::BuildTables()
348 {
349 #ifdef G4MULTITHREADED
350  G4MUTEXLOCK(&extrapolatorMutex);
351 #endif
352  if(verbose > 0) {
353  G4cout << "### G4EnergyLossForExtrapolator::BuildTables for "
354  << G4Material::GetNumberOfMaterials() << " materials Nbins= "
355  << nbins << " Emin(MeV)= " << emin << " Emax(MeV)= " << emax
356  << G4endl;
357  }
358  G4int newmat = G4Material::GetNumberOfMaterials();
359  if(!tables) {
360  tables = new G4TablesForExtrapolator(verbose, nbins, emin, emax);
361  } else if(nmat != newmat) {
362  tables->Initialisation();
363  }
364  nmat = newmat;
365 #ifdef G4MULTITHREADED
366  G4MUTEXUNLOCK(&extrapolatorMutex);
367 #endif
368 }
369 
370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
371