Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•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//
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
52G4EnergyLossTablesHelper *G4EnergyLossTables::t = nullptr;
53G4EnergyLossTablesHelper *G4EnergyLossTables::null_loss = nullptr;
54G4ParticleDefinition* G4EnergyLossTables::lastParticle = nullptr;
55G4double G4EnergyLossTables::QQPositron = 1.0; // e_squared
56G4double G4EnergyLossTables::Chargesquare ;
57G4int G4EnergyLossTables::oldIndex = -1 ;
58G4double G4EnergyLossTables::rmin = 0. ;
59G4double G4EnergyLossTables::rmax = 0. ;
60G4double G4EnergyLossTables::Thigh = 0. ;
61G4int G4EnergyLossTables::let_counter = 0;
62G4int G4EnergyLossTables::let_max_num_warnings = 100;
63G4bool G4EnergyLossTables::first_loss = true;
64
65G4EnergyLossTables::helper_map *G4EnergyLossTables::dict = nullptr;
66
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68
69G4EnergyLossTablesHelper::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
92G4EnergyLossTablesHelper::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
104void 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(
130 nullptr, nullptr, nullptr, nullptr, nullptr, 0.0, 0.0, 0.0, 0);
131 first_loss = false;
132 }
133}
134
135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136
137const 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 nullptr;
143 return (*it).second.theDEDXTable;
144}
145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147
148const 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 nullptr;
154 return (*it).second.theRangeTable;
155}
156
157//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158
159const 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 nullptr;
165 return (*it).second.theInverseRangeTable;
166}
167
168//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
169
170const 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 nullptr;
176 return (*it).second.theLabTimeTable;
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180
181const 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 nullptr;
187 return (*it).second.theProperTimeTable;
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191
192G4EnergyLossTablesHelper 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
207G4double 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 = (G4int)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
258G4double 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 = (G4int)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
308G4double 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 = (G4int)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
392G4double 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 = (G4int)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
442G4double 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 = (G4int)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
526G4double 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 = (G4int)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
581G4double 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 = (G4int)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 = (G4int)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 = (G4int)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
760G4double 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
813G4double 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
869G4double 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
940G4double 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
988G4double 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
1045void 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
1066void
1067G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition*,
1068 const G4String& /*q*/)
1069{
1070}
1071
1072//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4EnergyLossTables::GetRange
static G4double GetRange(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:526
G4EnergyLossTables::GetLabTimeTable
static const G4PhysicsTable * GetLabTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:170
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:104
G4EnergyLossTables::GetDEDX
static G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:207
G4EnergyLossTables::GetDeltaProperTime
static G4double GetDeltaProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:442
G4EnergyLossTables::GetLabTime
static G4double GetLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:258
G4EnergyLossTables::GetPreciseDEDX
static G4double GetPreciseDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:650
G4EnergyLossTables::GetPreciseEnergyFromRange
static G4double GetPreciseEnergyFromRange(const G4ParticleDefinition *aParticle, G4double range, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:581
G4EnergyLossTables::GetDEDXTable
static const G4PhysicsTable * GetDEDXTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:137
G4EnergyLossTables::GetRangeTable
static const G4PhysicsTable * GetRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:148
G4EnergyLossTables::GetInverseRangeTable
static const G4PhysicsTable * GetInverseRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:159
G4EnergyLossTables::GetDeltaLabTime
static G4double GetDeltaLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:308
G4EnergyLossTables::GetPreciseRangeFromEnergy
static G4double GetPreciseRangeFromEnergy(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:701
G4EnergyLossTables::GetProperTime
static G4double GetProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:392
G4EnergyLossTables::GetProperTimeTable
static const G4PhysicsTable * GetProperTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:181
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4LossTableManager::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *aParticle, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.hh:352
G4LossTableManager::GetRange
G4double GetRange(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.hh:341
G4LossTableManager::GetDEDX
G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.hh:307
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:255
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:113
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62