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