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