Geant4_10
G4EnergyLossTables.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 //
27 // $Id: G4EnergyLossTables.cc 77290 2013-11-22 10:57:50Z gcosmo $
28 //
29 // -------------------------------------------------------------------
30 // first version created by P.Urban , 06/04/1998
31 // modifications + "precise" functions added by L.Urban , 27/05/98
32 // modifications , TOF functions , 26/10/98, L.Urban
33 // cache mechanism in order to gain time, 11/02/99, L.Urban
34 // bug fixed , 12/04/99 , L.Urban
35 // 10.11.99: moved from RWT hash dictionary to STL map, G.Barrand, M.Maire
36 // 27.09.01 L.Urban , bug fixed (negative energy deposit)
37 // 26.10.01 all static functions moved from .icc files (mma)
38 // 15.01.03 Add interfaces required for "cut per region" (V.Ivanchenko)
39 // 12.03.03 Add warnings to obsolete interfaces (V.Ivanchenko)
40 // 10.04.03 Add call to G4LossTableManager is particle is not registered (V.Ivanchenko)
41 //
42 // -------------------------------------------------------------------
43 
44 #include "G4EnergyLossTables.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4MaterialCutsCouple.hh"
47 #include "G4RegionStore.hh"
48 #include "G4LossTableManager.hh"
49 
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 /*
53 G4ThreadLocal G4EnergyLossTablesHelper *G4EnergyLossTables::t = 0;
54 G4ThreadLocal G4EnergyLossTablesHelper *G4EnergyLossTables::null_loss = 0;
55 G4ThreadLocal G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0;
56 G4ThreadLocal G4double G4EnergyLossTables::QQPositron = 1.0; // e_squared
57 G4ThreadLocal G4double G4EnergyLossTables::Chargesquare ;
58 G4ThreadLocal G4int G4EnergyLossTables::oldIndex = -1 ;
59 G4ThreadLocal G4double G4EnergyLossTables::rmin = 0. ;
60 G4ThreadLocal G4double G4EnergyLossTables::rmax = 0. ;
61 G4ThreadLocal G4double G4EnergyLossTables::Thigh = 0. ;
62 G4ThreadLocal G4int G4EnergyLossTables::let_counter = 0;
63 G4ThreadLocal G4int G4EnergyLossTables::let_max_num_warnings = 100;
64 G4ThreadLocal G4bool G4EnergyLossTables::first_loss = true;
65 
66 G4ThreadLocal G4EnergyLossTables::helper_map *G4EnergyLossTables::dict = 0;
67 */
68 G4EnergyLossTablesHelper *G4EnergyLossTables::t = 0;
69 G4EnergyLossTablesHelper *G4EnergyLossTables::null_loss = 0;
70 G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0;
71 G4double G4EnergyLossTables::QQPositron = 1.0; // e_squared
72 G4double G4EnergyLossTables::Chargesquare ;
73 G4int G4EnergyLossTables::oldIndex = -1 ;
74 G4double G4EnergyLossTables::rmin = 0. ;
75 G4double G4EnergyLossTables::rmax = 0. ;
76 G4double G4EnergyLossTables::Thigh = 0. ;
77 G4int G4EnergyLossTables::let_counter = 0;
78 G4int G4EnergyLossTables::let_max_num_warnings = 100;
79 G4bool G4EnergyLossTables::first_loss = true;
80 
81 G4EnergyLossTables::helper_map *G4EnergyLossTables::dict = 0;
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
85 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper(
86  const G4PhysicsTable* aDEDXTable,
87  const G4PhysicsTable* aRangeTable,
88  const G4PhysicsTable* anInverseRangeTable,
89  const G4PhysicsTable* aLabTimeTable,
90  const G4PhysicsTable* aProperTimeTable,
91  G4double aLowestKineticEnergy,
92  G4double aHighestKineticEnergy,
93  G4double aMassRatio,
94  G4int aNumberOfBins)
95  :
96  theDEDXTable(aDEDXTable), theRangeTable(aRangeTable),
97  theInverseRangeTable(anInverseRangeTable),
98  theLabTimeTable(aLabTimeTable),
99  theProperTimeTable(aProperTimeTable),
100  theLowestKineticEnergy(aLowestKineticEnergy),
101  theHighestKineticEnergy(aHighestKineticEnergy),
102  theMassRatio(aMassRatio),
103  theNumberOfBins(aNumberOfBins)
104 { }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 
108 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper()
109 {
110  theLowestKineticEnergy = 0.0;
111  theHighestKineticEnergy= 0.0;
112  theMassRatio = 0.0;
113  theNumberOfBins = 0;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 void G4EnergyLossTables::Register(
119  const G4ParticleDefinition* p,
120  const G4PhysicsTable* tDEDX,
121  const G4PhysicsTable* tRange,
122  const G4PhysicsTable* tInverseRange,
123  const G4PhysicsTable* tLabTime,
124  const G4PhysicsTable* tProperTime,
125  G4double lowestKineticEnergy,
126  G4double highestKineticEnergy,
127  G4double massRatio,
128  G4int NumberOfBins)
129 {
130  if (!dict) dict = new G4EnergyLossTables::helper_map;
131  if (!null_loss) null_loss = new G4EnergyLossTablesHelper;
132  if (!t) t = new G4EnergyLossTablesHelper;
133 
134  (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, tRange,tInverseRange,
135  tLabTime,tProperTime,lowestKineticEnergy,
136  highestKineticEnergy, massRatio,NumberOfBins);
137 
138  *t = GetTables(p) ; // important for cache !!!!!
139  lastParticle = (G4ParticleDefinition*) p ;
140  Chargesquare = (p->GetPDGCharge())*(p->GetPDGCharge())/
141  QQPositron ;
142  if (first_loss ) {
143  *null_loss = G4EnergyLossTablesHelper(0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0);
144  first_loss = false;
145  }
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
150 const G4PhysicsTable* G4EnergyLossTables::GetDEDXTable(
151  const G4ParticleDefinition* p)
152 {
153  if (!dict) dict = new G4EnergyLossTables::helper_map;
154  helper_map::iterator it;
155  if((it=dict->find(p))==dict->end()) return 0;
156  return (*it).second.theDEDXTable;
157 }
158 
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
160 
161 const G4PhysicsTable* G4EnergyLossTables::GetRangeTable(
162  const G4ParticleDefinition* p)
163 {
164  if (!dict) dict = new G4EnergyLossTables::helper_map;
165  helper_map::iterator it;
166  if((it=dict->find(p))==dict->end()) return 0;
167  return (*it).second.theRangeTable;
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171 
172 const G4PhysicsTable* G4EnergyLossTables::GetInverseRangeTable(
173  const G4ParticleDefinition* p)
174 {
175  if (!dict) dict = new G4EnergyLossTables::helper_map;
176  helper_map::iterator it;
177  if((it=dict->find(p))==dict->end()) return 0;
178  return (*it).second.theInverseRangeTable;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
183 const G4PhysicsTable* G4EnergyLossTables::GetLabTimeTable(
184  const G4ParticleDefinition* p)
185 {
186  if (!dict) dict = new G4EnergyLossTables::helper_map;
187  helper_map::iterator it;
188  if((it=dict->find(p))==dict->end()) return 0;
189  return (*it).second.theLabTimeTable;
190 }
191 
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193 
194 const G4PhysicsTable* G4EnergyLossTables::GetProperTimeTable(
195  const G4ParticleDefinition* p)
196 {
197  if (!dict) dict = new G4EnergyLossTables::helper_map;
198  helper_map::iterator it;
199  if((it=dict->find(p))==dict->end()) return 0;
200  return (*it).second.theProperTimeTable;
201 }
202 
203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
204 
205 G4EnergyLossTablesHelper G4EnergyLossTables::GetTables(
206  const G4ParticleDefinition* p)
207 {
208  if (!dict) dict = new G4EnergyLossTables::helper_map;
209  if (!null_loss) null_loss = new G4EnergyLossTablesHelper;
210 
211  helper_map::iterator it;
212  if ((it=dict->find(p))==dict->end()) {
213  return *null_loss;
214  }
215  return (*it).second;
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 
220 G4double G4EnergyLossTables::GetDEDX(
221  const G4ParticleDefinition *aParticle,
222  G4double KineticEnergy,
223  const G4Material *aMaterial)
224 {
225  if (!t) t = new G4EnergyLossTablesHelper;
226 
227  CPRWarning();
228  if(aParticle != (const G4ParticleDefinition*) lastParticle)
229  {
230  *t= GetTables(aParticle);
231  lastParticle = (G4ParticleDefinition*) aParticle ;
232  Chargesquare = (aParticle->GetPDGCharge())*
233  (aParticle->GetPDGCharge())/
234  QQPositron ;
235  oldIndex = -1 ;
236  }
237  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
238  if (!dEdxTable) {
239  ParticleHaveNoLoss(aParticle,"dEdx");
240  return 0.0;
241  }
242 
243  G4int materialIndex = aMaterial->GetIndex();
244  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
245  G4double dEdx;
246  G4bool isOut;
247 
248  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
249 
250  dEdx =(*dEdxTable)(materialIndex)->GetValue(
251  t->theLowestKineticEnergy,isOut)
252  *std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy);
253 
254  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
255 
256  dEdx = (*dEdxTable)(materialIndex)->GetValue(
257  t->theHighestKineticEnergy,isOut);
258 
259  } else {
260 
261  dEdx = (*dEdxTable)(materialIndex)->GetValue(
262  scaledKineticEnergy,isOut);
263 
264  }
265 
266  return dEdx*Chargesquare;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 
271 G4double G4EnergyLossTables::GetLabTime(
272  const G4ParticleDefinition *aParticle,
273  G4double KineticEnergy,
274  const G4Material *aMaterial)
275 {
276  if (!t) t = new G4EnergyLossTablesHelper;
277 
278  CPRWarning();
279  if(aParticle != (const G4ParticleDefinition*) lastParticle)
280  {
281  *t= GetTables(aParticle);
282  lastParticle = (G4ParticleDefinition*) aParticle ;
283  oldIndex = -1 ;
284  }
285  const G4PhysicsTable* labtimeTable= t->theLabTimeTable;
286  if (!labtimeTable) {
287  ParticleHaveNoLoss(aParticle,"LabTime");
288  return 0.0;
289  }
290 
291  const G4double parlowen=0.4 , ppar=0.5-parlowen ;
292  G4int materialIndex = aMaterial->GetIndex();
293  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
294  G4double time;
295  G4bool isOut;
296 
297  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
298 
299  time = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))*
300  (*labtimeTable)(materialIndex)->GetValue(
301  t->theLowestKineticEnergy,isOut);
302 
303 
304  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
305 
306  time = (*labtimeTable)(materialIndex)->GetValue(
307  t->theHighestKineticEnergy,isOut);
308 
309  } else {
310 
311  time = (*labtimeTable)(materialIndex)->GetValue(
312  scaledKineticEnergy,isOut);
313 
314  }
315 
316  return time/t->theMassRatio ;
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320 
321 G4double G4EnergyLossTables::GetDeltaLabTime(
322  const G4ParticleDefinition *aParticle,
323  G4double KineticEnergyStart,
324  G4double KineticEnergyEnd,
325  const G4Material *aMaterial)
326 {
327  if (!t) t = new G4EnergyLossTablesHelper;
328 
329  CPRWarning();
330  if(aParticle != (const G4ParticleDefinition*) lastParticle)
331  {
332  *t= GetTables(aParticle);
333  lastParticle = (G4ParticleDefinition*) aParticle ;
334  oldIndex = -1 ;
335  }
336  const G4PhysicsTable* labtimeTable= t->theLabTimeTable;
337  if (!labtimeTable) {
338  ParticleHaveNoLoss(aParticle,"LabTime");
339  return 0.0;
340  }
341 
342  const G4double parlowen=0.4 , ppar=0.5-parlowen ;
343  const G4double dToverT = 0.05 , facT = 1. -dToverT ;
344  G4double timestart,timeend,deltatime,dTT;
345  G4bool isOut;
346 
347  G4int materialIndex = aMaterial->GetIndex();
348  G4double scaledKineticEnergy = KineticEnergyStart*t->theMassRatio;
349 
350  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
351 
352  timestart = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))*
353  (*labtimeTable)(materialIndex)->GetValue(
354  t->theLowestKineticEnergy,isOut);
355 
356 
357  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
358 
359  timestart = (*labtimeTable)(materialIndex)->GetValue(
360  t->theHighestKineticEnergy,isOut);
361 
362  } else {
363 
364  timestart = (*labtimeTable)(materialIndex)->GetValue(
365  scaledKineticEnergy,isOut);
366 
367  }
368 
369  dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ;
370 
371  if( dTT < dToverT )
372  scaledKineticEnergy = facT*KineticEnergyStart*t->theMassRatio;
373  else
374  scaledKineticEnergy = KineticEnergyEnd*t->theMassRatio;
375 
376  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
377 
378  timeend = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))*
379  (*labtimeTable)(materialIndex)->GetValue(
380  t->theLowestKineticEnergy,isOut);
381 
382 
383  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
384 
385  timeend = (*labtimeTable)(materialIndex)->GetValue(
386  t->theHighestKineticEnergy,isOut);
387 
388  } else {
389 
390  timeend = (*labtimeTable)(materialIndex)->GetValue(
391  scaledKineticEnergy,isOut);
392 
393  }
394 
395  deltatime = timestart - timeend ;
396 
397  if( dTT < dToverT )
398  deltatime *= dTT/dToverT;
399 
400  return deltatime/t->theMassRatio ;
401 }
402 
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 
405 G4double G4EnergyLossTables::GetProperTime(
406  const G4ParticleDefinition *aParticle,
407  G4double KineticEnergy,
408  const G4Material *aMaterial)
409 {
410  if (!t) t = new G4EnergyLossTablesHelper;
411 
412  CPRWarning();
413  if(aParticle != (const G4ParticleDefinition*) lastParticle)
414  {
415  *t= GetTables(aParticle);
416  lastParticle = (G4ParticleDefinition*) aParticle ;
417  oldIndex = -1 ;
418  }
419  const G4PhysicsTable* propertimeTable= t->theProperTimeTable;
420  if (!propertimeTable) {
421  ParticleHaveNoLoss(aParticle,"ProperTime");
422  return 0.0;
423  }
424 
425  const G4double parlowen=0.4 , ppar=0.5-parlowen ;
426  G4int materialIndex = aMaterial->GetIndex();
427  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
428  G4double time;
429  G4bool isOut;
430 
431  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
432 
433  time = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))*
434  (*propertimeTable)(materialIndex)->GetValue(
435  t->theLowestKineticEnergy,isOut);
436 
437 
438  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
439 
440  time = (*propertimeTable)(materialIndex)->GetValue(
441  t->theHighestKineticEnergy,isOut);
442 
443  } else {
444 
445  time = (*propertimeTable)(materialIndex)->GetValue(
446  scaledKineticEnergy,isOut);
447 
448  }
449 
450  return time/t->theMassRatio ;
451 }
452 
453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
454 
455 G4double G4EnergyLossTables::GetDeltaProperTime(
456  const G4ParticleDefinition *aParticle,
457  G4double KineticEnergyStart,
458  G4double KineticEnergyEnd,
459  const G4Material *aMaterial)
460 {
461  if (!t) t = new G4EnergyLossTablesHelper;
462 
463  CPRWarning();
464  if(aParticle != (const G4ParticleDefinition*) lastParticle)
465  {
466  *t= GetTables(aParticle);
467  lastParticle = (G4ParticleDefinition*) aParticle ;
468  oldIndex = -1 ;
469  }
470  const G4PhysicsTable* propertimeTable= t->theProperTimeTable;
471  if (!propertimeTable) {
472  ParticleHaveNoLoss(aParticle,"ProperTime");
473  return 0.0;
474  }
475 
476  const G4double parlowen=0.4 , ppar=0.5-parlowen ;
477  const G4double dToverT = 0.05 , facT = 1. -dToverT ;
478  G4double timestart,timeend,deltatime,dTT;
479  G4bool isOut;
480 
481  G4int materialIndex = aMaterial->GetIndex();
482  G4double scaledKineticEnergy = KineticEnergyStart*t->theMassRatio;
483 
484  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
485 
486  timestart = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))*
487  (*propertimeTable)(materialIndex)->GetValue(
488  t->theLowestKineticEnergy,isOut);
489 
490 
491  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
492 
493  timestart = (*propertimeTable)(materialIndex)->GetValue(
494  t->theHighestKineticEnergy,isOut);
495 
496  } else {
497 
498  timestart = (*propertimeTable)(materialIndex)->GetValue(
499  scaledKineticEnergy,isOut);
500 
501  }
502 
503  dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ;
504 
505  if( dTT < dToverT )
506  scaledKineticEnergy = facT*KineticEnergyStart*t->theMassRatio;
507  else
508  scaledKineticEnergy = KineticEnergyEnd*t->theMassRatio;
509 
510  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
511 
512  timeend = std::exp(ppar*std::log(scaledKineticEnergy/t->theLowestKineticEnergy))*
513  (*propertimeTable)(materialIndex)->GetValue(
514  t->theLowestKineticEnergy,isOut);
515 
516 
517  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
518 
519  timeend = (*propertimeTable)(materialIndex)->GetValue(
520  t->theHighestKineticEnergy,isOut);
521 
522  } else {
523 
524  timeend = (*propertimeTable)(materialIndex)->GetValue(
525  scaledKineticEnergy,isOut);
526 
527  }
528 
529  deltatime = timestart - timeend ;
530 
531  if( dTT < dToverT )
532  deltatime *= dTT/dToverT ;
533 
534  return deltatime/t->theMassRatio ;
535 }
536 
537 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
538 
539 G4double G4EnergyLossTables::GetRange(
540  const G4ParticleDefinition *aParticle,
541  G4double KineticEnergy,
542  const G4Material *aMaterial)
543 {
544  if (!t) t = new G4EnergyLossTablesHelper;
545 
546  CPRWarning();
547  if(aParticle != (const G4ParticleDefinition*) lastParticle)
548  {
549  *t= GetTables(aParticle);
550  lastParticle = (G4ParticleDefinition*) aParticle ;
551  Chargesquare = (aParticle->GetPDGCharge())*
552  (aParticle->GetPDGCharge())/
553  QQPositron ;
554  oldIndex = -1 ;
555  }
556  const G4PhysicsTable* rangeTable= t->theRangeTable;
557  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
558  if (!rangeTable) {
559  ParticleHaveNoLoss(aParticle,"Range");
560  return 0.0;
561  }
562 
563  G4int materialIndex = aMaterial->GetIndex();
564  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
565  G4double Range;
566  G4bool isOut;
567 
568  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
569 
570  Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)*
571  (*rangeTable)(materialIndex)->GetValue(
572  t->theLowestKineticEnergy,isOut);
573 
574  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
575 
576  Range = (*rangeTable)(materialIndex)->GetValue(
577  t->theHighestKineticEnergy,isOut)+
578  (scaledKineticEnergy-t->theHighestKineticEnergy)/
579  (*dEdxTable)(materialIndex)->GetValue(
580  t->theHighestKineticEnergy,isOut);
581 
582  } else {
583 
584  Range = (*rangeTable)(materialIndex)->GetValue(
585  scaledKineticEnergy,isOut);
586 
587  }
588 
589  return Range/(Chargesquare*t->theMassRatio);
590 }
591 
592 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
593 
594 G4double G4EnergyLossTables::GetPreciseEnergyFromRange(
595  const G4ParticleDefinition *aParticle,
596  G4double range,
597  const G4Material *aMaterial)
598 // it returns the value of the kinetic energy for a given range
599 {
600  if (!t) t = new G4EnergyLossTablesHelper;
601 
602  CPRWarning();
603  if( aParticle != (const G4ParticleDefinition*) lastParticle)
604  {
605  *t= GetTables(aParticle);
606  lastParticle = (G4ParticleDefinition*) aParticle;
607  Chargesquare = (aParticle->GetPDGCharge())*
608  (aParticle->GetPDGCharge())/
609  QQPositron ;
610  oldIndex = -1 ;
611  }
612  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
613  const G4PhysicsTable* inverseRangeTable= t->theInverseRangeTable;
614  if (!inverseRangeTable) {
615  ParticleHaveNoLoss(aParticle,"InverseRange");
616  return 0.0;
617  }
618 
619  G4double scaledrange,scaledKineticEnergy ;
620  G4bool isOut ;
621 
622  G4int materialIndex = aMaterial->GetIndex() ;
623 
624  if(materialIndex != oldIndex)
625  {
626  oldIndex = materialIndex ;
627  rmin = (*inverseRangeTable)(materialIndex)->
628  GetLowEdgeEnergy(0) ;
629  rmax = (*inverseRangeTable)(materialIndex)->
630  GetLowEdgeEnergy(t->theNumberOfBins-2) ;
631  Thigh = (*inverseRangeTable)(materialIndex)->
632  GetValue(rmax,isOut) ;
633  }
634 
635  scaledrange = range*Chargesquare*t->theMassRatio ;
636 
637  if(scaledrange < rmin)
638  {
639  scaledKineticEnergy = t->theLowestKineticEnergy*
640  scaledrange*scaledrange/(rmin*rmin) ;
641  }
642  else
643  {
644  if(scaledrange < rmax)
645  {
646  scaledKineticEnergy = (*inverseRangeTable)(materialIndex)->
647  GetValue( scaledrange,isOut) ;
648  }
649  else
650  {
651  scaledKineticEnergy = Thigh +
652  (scaledrange-rmax)*
653  (*dEdxTable)(materialIndex)->
654  GetValue(Thigh,isOut) ;
655  }
656  }
657 
658  return scaledKineticEnergy/t->theMassRatio ;
659 }
660 
661 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
662 
663  G4double G4EnergyLossTables::GetPreciseDEDX(
664  const G4ParticleDefinition *aParticle,
665  G4double KineticEnergy,
666  const G4Material *aMaterial)
667 {
668  if (!t) t = new G4EnergyLossTablesHelper;
669 
670  CPRWarning();
671  if( aParticle != (const G4ParticleDefinition*) lastParticle)
672  {
673  *t= GetTables(aParticle);
674  lastParticle = (G4ParticleDefinition*) aParticle;
675  Chargesquare = (aParticle->GetPDGCharge())*
676  (aParticle->GetPDGCharge())/
677  QQPositron ;
678  oldIndex = -1 ;
679  }
680  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
681  if (!dEdxTable) {
682  ParticleHaveNoLoss(aParticle,"dEdx");
683  return 0.0;
684  }
685 
686  G4int materialIndex = aMaterial->GetIndex();
687  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
688  G4double dEdx;
689  G4bool isOut;
690 
691  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
692 
693  dEdx = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)
694  *(*dEdxTable)(materialIndex)->GetValue(
695  t->theLowestKineticEnergy,isOut);
696 
697  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
698 
699  dEdx = (*dEdxTable)(materialIndex)->GetValue(
700  t->theHighestKineticEnergy,isOut);
701 
702  } else {
703 
704  dEdx = (*dEdxTable)(materialIndex)->GetValue(
705  scaledKineticEnergy,isOut) ;
706 
707  }
708 
709  return dEdx*Chargesquare;
710 }
711 
712 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
713 
714  G4double G4EnergyLossTables::GetPreciseRangeFromEnergy(
715  const G4ParticleDefinition *aParticle,
716  G4double KineticEnergy,
717  const G4Material *aMaterial)
718 {
719  if (!t) t = new G4EnergyLossTablesHelper;
720 
721  CPRWarning();
722  if( aParticle != (const G4ParticleDefinition*) lastParticle)
723  {
724  *t= GetTables(aParticle);
725  lastParticle = (G4ParticleDefinition*) aParticle;
726  Chargesquare = (aParticle->GetPDGCharge())*
727  (aParticle->GetPDGCharge())/
728  QQPositron ;
729  oldIndex = -1 ;
730  }
731  const G4PhysicsTable* rangeTable= t->theRangeTable;
732  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
733  if (!rangeTable) {
734  ParticleHaveNoLoss(aParticle,"Range");
735  return 0.0;
736  }
737  G4int materialIndex = aMaterial->GetIndex();
738 
739  G4double Thighr = t->theHighestKineticEnergy*t->theLowestKineticEnergy/
740  (*rangeTable)(materialIndex)->
741  GetLowEdgeEnergy(1) ;
742 
743  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
744  G4double Range;
745  G4bool isOut;
746 
747  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
748 
749  Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)*
750  (*rangeTable)(materialIndex)->GetValue(
751  t->theLowestKineticEnergy,isOut);
752 
753  } else if (scaledKineticEnergy>Thighr) {
754 
755  Range = (*rangeTable)(materialIndex)->GetValue(
756  Thighr,isOut)+
757  (scaledKineticEnergy-Thighr)/
758  (*dEdxTable)(materialIndex)->GetValue(
759  Thighr,isOut);
760 
761  } else {
762 
763  Range = (*rangeTable)(materialIndex)->GetValue(
764  scaledKineticEnergy,isOut) ;
765 
766  }
767 
768  return Range/(Chargesquare*t->theMassRatio);
769 }
770 
771 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
772 
773 G4double G4EnergyLossTables::GetDEDX(
774  const G4ParticleDefinition *aParticle,
775  G4double KineticEnergy,
776  const G4MaterialCutsCouple *couple,
777  G4bool check)
778 {
779  if (!t) t = new G4EnergyLossTablesHelper;
780 
781  if(aParticle != (const G4ParticleDefinition*) lastParticle)
782  {
783  *t= GetTables(aParticle);
784  lastParticle = (G4ParticleDefinition*) aParticle ;
785  Chargesquare = (aParticle->GetPDGCharge())*
786  (aParticle->GetPDGCharge())/
787  QQPositron ;
788  oldIndex = -1 ;
789  }
790  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
791 
792  if (!dEdxTable ) {
793  if (check) return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple);
794  else ParticleHaveNoLoss(aParticle, "dEdx");
795  return 0.0;
796  }
797 
798  G4int materialIndex = couple->GetIndex();
799  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
800  G4double dEdx;
801  G4bool isOut;
802 
803  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
804 
805  dEdx =(*dEdxTable)(materialIndex)->GetValue(
806  t->theLowestKineticEnergy,isOut)
807  *std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy);
808 
809  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
810 
811  dEdx = (*dEdxTable)(materialIndex)->GetValue(
812  t->theHighestKineticEnergy,isOut);
813 
814  } else {
815 
816  dEdx = (*dEdxTable)(materialIndex)->GetValue(
817  scaledKineticEnergy,isOut);
818 
819  }
820 
821  return dEdx*Chargesquare;
822 }
823 
824 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
825 
826 G4double G4EnergyLossTables::GetRange(
827  const G4ParticleDefinition *aParticle,
828  G4double KineticEnergy,
829  const G4MaterialCutsCouple *couple,
830  G4bool check)
831 {
832  if (!t) t = new G4EnergyLossTablesHelper;
833 
834  if(aParticle != (const G4ParticleDefinition*) lastParticle)
835  {
836  *t= GetTables(aParticle);
837  lastParticle = (G4ParticleDefinition*) aParticle ;
838  Chargesquare = (aParticle->GetPDGCharge())*
839  (aParticle->GetPDGCharge())/
840  QQPositron ;
841  oldIndex = -1 ;
842  }
843  const G4PhysicsTable* rangeTable= t->theRangeTable;
844  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
845  if (!rangeTable) {
846  if(check) return G4LossTableManager::Instance()->GetRange(aParticle,KineticEnergy,couple);
847  else return DBL_MAX;
848  //ParticleHaveNoLoss(aParticle,"Range");
849  }
850 
851  G4int materialIndex = couple->GetIndex();
852  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
853  G4double Range;
854  G4bool isOut;
855 
856  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
857 
858  Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)*
859  (*rangeTable)(materialIndex)->GetValue(
860  t->theLowestKineticEnergy,isOut);
861 
862  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
863 
864  Range = (*rangeTable)(materialIndex)->GetValue(
865  t->theHighestKineticEnergy,isOut)+
866  (scaledKineticEnergy-t->theHighestKineticEnergy)/
867  (*dEdxTable)(materialIndex)->GetValue(
868  t->theHighestKineticEnergy,isOut);
869 
870  } else {
871 
872  Range = (*rangeTable)(materialIndex)->GetValue(
873  scaledKineticEnergy,isOut);
874 
875  }
876 
877  return Range/(Chargesquare*t->theMassRatio);
878 }
879 
880 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
881 
882 G4double G4EnergyLossTables::GetPreciseEnergyFromRange(
883  const G4ParticleDefinition *aParticle,
884  G4double range,
885  const G4MaterialCutsCouple *couple,
886  G4bool check)
887 // it returns the value of the kinetic energy for a given range
888 {
889  if (!t) t = new G4EnergyLossTablesHelper;
890 
891  if( aParticle != (const G4ParticleDefinition*) lastParticle)
892  {
893  *t= GetTables(aParticle);
894  lastParticle = (G4ParticleDefinition*) aParticle;
895  Chargesquare = (aParticle->GetPDGCharge())*
896  (aParticle->GetPDGCharge())/
897  QQPositron ;
898  oldIndex = -1 ;
899  }
900  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
901  const G4PhysicsTable* inverseRangeTable= t->theInverseRangeTable;
902 
903  if (!inverseRangeTable) {
904  if(check) return G4LossTableManager::Instance()->GetEnergy(aParticle,range,couple);
905  else return DBL_MAX;
906  // else ParticleHaveNoLoss(aParticle,"InverseRange");
907  }
908 
909  G4double scaledrange,scaledKineticEnergy ;
910  G4bool isOut ;
911 
912  G4int materialIndex = couple->GetIndex() ;
913 
914  if(materialIndex != oldIndex)
915  {
916  oldIndex = materialIndex ;
917  rmin = (*inverseRangeTable)(materialIndex)->
918  GetLowEdgeEnergy(0) ;
919  rmax = (*inverseRangeTable)(materialIndex)->
920  GetLowEdgeEnergy(t->theNumberOfBins-2) ;
921  Thigh = (*inverseRangeTable)(materialIndex)->
922  GetValue(rmax,isOut) ;
923  }
924 
925  scaledrange = range*Chargesquare*t->theMassRatio ;
926 
927  if(scaledrange < rmin)
928  {
929  scaledKineticEnergy = t->theLowestKineticEnergy*
930  scaledrange*scaledrange/(rmin*rmin) ;
931  }
932  else
933  {
934  if(scaledrange < rmax)
935  {
936  scaledKineticEnergy = (*inverseRangeTable)(materialIndex)->
937  GetValue( scaledrange,isOut) ;
938  }
939  else
940  {
941  scaledKineticEnergy = Thigh +
942  (scaledrange-rmax)*
943  (*dEdxTable)(materialIndex)->
944  GetValue(Thigh,isOut) ;
945  }
946  }
947 
948  return scaledKineticEnergy/t->theMassRatio ;
949 }
950 
951 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
952 
953 G4double G4EnergyLossTables::GetPreciseDEDX(
954  const G4ParticleDefinition *aParticle,
955  G4double KineticEnergy,
956  const G4MaterialCutsCouple *couple)
957 {
958  if (!t) t = new G4EnergyLossTablesHelper;
959 
960  if( aParticle != (const G4ParticleDefinition*) lastParticle)
961  {
962  *t= GetTables(aParticle);
963  lastParticle = (G4ParticleDefinition*) aParticle;
964  Chargesquare = (aParticle->GetPDGCharge())*
965  (aParticle->GetPDGCharge())/
966  QQPositron ;
967  oldIndex = -1 ;
968  }
969  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
970  if ( !dEdxTable )
971  return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple);
972 
973  G4int materialIndex = couple->GetIndex();
974  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
975  G4double dEdx;
976  G4bool isOut;
977 
978  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
979 
980  dEdx = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)
981  *(*dEdxTable)(materialIndex)->GetValue(
982  t->theLowestKineticEnergy,isOut);
983 
984  } else if (scaledKineticEnergy>t->theHighestKineticEnergy) {
985 
986  dEdx = (*dEdxTable)(materialIndex)->GetValue(
987  t->theHighestKineticEnergy,isOut);
988 
989  } else {
990 
991  dEdx = (*dEdxTable)(materialIndex)->GetValue(
992  scaledKineticEnergy,isOut) ;
993 
994  }
995 
996  return dEdx*Chargesquare;
997 }
998 
999 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1000 
1001 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy(
1002  const G4ParticleDefinition *aParticle,
1003  G4double KineticEnergy,
1004  const G4MaterialCutsCouple *couple)
1005 {
1006  if (!t) t = new G4EnergyLossTablesHelper;
1007 
1008  if( aParticle != (const G4ParticleDefinition*) lastParticle)
1009  {
1010  *t= GetTables(aParticle);
1011  lastParticle = (G4ParticleDefinition*) aParticle;
1012  Chargesquare = (aParticle->GetPDGCharge())*
1013  (aParticle->GetPDGCharge())/
1014  QQPositron ;
1015  oldIndex = -1 ;
1016  }
1017  const G4PhysicsTable* rangeTable= t->theRangeTable;
1018  const G4PhysicsTable* dEdxTable= t->theDEDXTable;
1019  if ( !dEdxTable || !rangeTable)
1020  return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple);
1021 
1022  G4int materialIndex = couple->GetIndex();
1023 
1024  G4double Thighr = t->theHighestKineticEnergy*t->theLowestKineticEnergy/
1025  (*rangeTable)(materialIndex)->
1026  GetLowEdgeEnergy(1) ;
1027 
1028  G4double scaledKineticEnergy = KineticEnergy*t->theMassRatio;
1029  G4double Range;
1030  G4bool isOut;
1031 
1032  if (scaledKineticEnergy<t->theLowestKineticEnergy) {
1033 
1034  Range = std::sqrt(scaledKineticEnergy/t->theLowestKineticEnergy)*
1035  (*rangeTable)(materialIndex)->GetValue(
1036  t->theLowestKineticEnergy,isOut);
1037 
1038  } else if (scaledKineticEnergy>Thighr) {
1039 
1040  Range = (*rangeTable)(materialIndex)->GetValue(
1041  Thighr,isOut)+
1042  (scaledKineticEnergy-Thighr)/
1043  (*dEdxTable)(materialIndex)->GetValue(
1044  Thighr,isOut);
1045 
1046  } else {
1047 
1048  Range = (*rangeTable)(materialIndex)->GetValue(
1049  scaledKineticEnergy,isOut) ;
1050 
1051  }
1052 
1053  return Range/(Chargesquare*t->theMassRatio);
1054 }
1055 
1056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1057 
1058 void G4EnergyLossTables::CPRWarning()
1059 {
1060  if (let_counter < let_max_num_warnings) {
1061 
1062  G4cout << G4endl;
1063  G4cout << "##### G4EnergyLossTable WARNING: The obsolete interface is used!" << G4endl;
1064  G4cout << "##### RESULTS ARE NOT GARANTEED!" << G4endl;
1065  G4cout << "##### Please, substitute G4Material by G4MaterialCutsCouple" << G4endl;
1066  G4cout << "##### Obsolete interface will be removed soon" << G4endl;
1067  G4cout << G4endl;
1068  let_counter++;
1069 
1070  } else if (let_counter == let_max_num_warnings) {
1071 
1072  G4cout << "##### G4EnergyLossTable WARNING closed" << G4endl;
1073  let_counter++;
1074  }
1075 }
1076 
1077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1078 
1079 void
1080 G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition*,
1081  const G4String& /*q*/)
1082 {
1083 }
1084 
1085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4EnergyLossTables::GetDeltaProperTime
static G4double GetDeltaProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:455
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:111
G4LossTableManager::GetDEDX
G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.cc:1312
G4EnergyLossTables::GetDEDXTable
static const G4PhysicsTable * GetDEDXTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:150
G4MaterialCutsCouple.hh
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:260
p
const char * p
Definition: xmltok.h:285
G4LossTableManager.hh
G4EnergyLossTables::GetDeltaLabTime
static G4double GetDeltaLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:321
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
pyG4EmCalculator::range
const G4ParticleDefinition const G4Material *G4double range
Definition: pyG4EmCalculator.cc:185
G4RegionStore.hh
G4bool
bool G4bool
Definition: G4Types.hh:79
G4EnergyLossTables::GetLabTime
static G4double GetLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:271
G4EnergyLossTables::GetPreciseRangeFromEnergy
static G4double GetPreciseRangeFromEnergy(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:714
G4LossTableManager::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *aParticle, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.cc:1373
G4EnergyLossTables::GetLabTimeTable
static const G4PhysicsTable * GetLabTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:183
G4EnergyLossTables::GetRangeTable
static const G4PhysicsTable * GetRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:161
G4endl
#define G4endl
Definition: G4ios.hh:61
G4EnergyLossTables::GetDEDX
static G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:220
G4EnergyLossTables::GetProperTime
static G4double GetProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:405
G4double
double G4double
Definition: G4Types.hh:76
G4EnergyLossTables::Register
static void Register(const G4ParticleDefinition *p, const G4PhysicsTable *tDEDX, const G4PhysicsTable *tRange, const G4PhysicsTable *tInverseRange, const G4PhysicsTable *tLabTime, const G4PhysicsTable *tProperTime, G4double lowestKineticEnergy, G4double highestKineticEnergy, G4double massRatio, G4int NumberOfBins)
Definition: G4EnergyLossTables.cc:118
G4SystemOfUnits.hh
G4EnergyLossTables::GetInverseRangeTable
static const G4PhysicsTable * GetInverseRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:172
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4EnergyLossTables::GetPreciseEnergyFromRange
static G4double GetPreciseEnergyFromRange(const G4ParticleDefinition *aParticle, G4double range, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:594
G4EnergyLossTables::GetPreciseDEDX
static G4double GetPreciseDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:663
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4LossTableManager::GetRange
G4double GetRange(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.cc:1361
G4EnergyLossTables.hh
G4EnergyLossTables::GetProperTimeTable
static const G4PhysicsTable * GetProperTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:194
G4EnergyLossTables::GetRange
static G4double GetRange(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:539
Generated on Sat Dec 14 2013 14:34:38 for Geant4_10 by   doxygen 1.8.5