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