Geant4-11
•All Data Structures 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 = 0;
53G4EnergyLossTablesHelper *G4EnergyLossTables::null_loss = 0;
54G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0;
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 = 0;
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(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
136const 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
147const 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
158const 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
169const 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
180const 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
191G4EnergyLossTablesHelper 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
206G4double 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
257G4double 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
307G4double 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
391G4double 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
441G4double 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
525G4double 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
580G4double 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
759G4double 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
812G4double 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
868G4double 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
939G4double 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
987G4double 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
1044void 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
1065void
1066G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition*,
1067 const G4String& /*q*/)
1068{
1069}
1070
1071//....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
G4EnergyLossTablesHelper::theProperTimeTable
const G4PhysicsTable * theProperTimeTable
Definition: G4EnergyLossTables.hh:92
G4EnergyLossTablesHelper::theInverseRangeTable
const G4PhysicsTable * theInverseRangeTable
Definition: G4EnergyLossTables.hh:90
G4EnergyLossTablesHelper::theDEDXTable
const G4PhysicsTable * theDEDXTable
Definition: G4EnergyLossTables.hh:88
G4EnergyLossTablesHelper::theRangeTable
const G4PhysicsTable * theRangeTable
Definition: G4EnergyLossTables.hh:89
G4EnergyLossTablesHelper::theLabTimeTable
const G4PhysicsTable * theLabTimeTable
Definition: G4EnergyLossTables.hh:91
G4EnergyLossTables::helper_map
std::map< K, G4EnergyLossTablesHelper, std::less< K > > helper_map
Definition: G4EnergyLossTables.hh:226
G4EnergyLossTables::GetRange
static G4double GetRange(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:525
G4EnergyLossTables::GetLabTimeTable
static const G4PhysicsTable * GetLabTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:169
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::GetTables
static G4EnergyLossTablesHelper GetTables(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:191
G4EnergyLossTables::GetDEDX
static G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:206
G4EnergyLossTables::GetDeltaProperTime
static G4double GetDeltaProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:441
G4EnergyLossTables::GetLabTime
static G4double GetLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:257
G4EnergyLossTables::QQPositron
static G4double QQPositron
Definition: G4EnergyLossTables.hh:234
G4EnergyLossTables::GetPreciseDEDX
static G4double GetPreciseDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:649
G4EnergyLossTables::GetPreciseEnergyFromRange
static G4double GetPreciseEnergyFromRange(const G4ParticleDefinition *aParticle, G4double range, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:580
G4EnergyLossTables::GetDEDXTable
static const G4PhysicsTable * GetDEDXTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:136
G4EnergyLossTables::GetRangeTable
static const G4PhysicsTable * GetRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:147
G4EnergyLossTables::GetInverseRangeTable
static const G4PhysicsTable * GetInverseRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:158
G4EnergyLossTables::let_max_num_warnings
static G4int let_max_num_warnings
Definition: G4EnergyLossTables.hh:239
G4EnergyLossTables::dict
static helper_map * dict
Definition: G4EnergyLossTables.hh:227
G4EnergyLossTables::GetDeltaLabTime
static G4double GetDeltaLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:307
G4EnergyLossTables::null_loss
static G4EnergyLossTablesHelper * null_loss
Definition: G4EnergyLossTables.hh:232
G4EnergyLossTables::GetPreciseRangeFromEnergy
static G4double GetPreciseRangeFromEnergy(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:700
G4EnergyLossTables::GetProperTime
static G4double GetProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:391
G4EnergyLossTables::ParticleHaveNoLoss
static void ParticleHaveNoLoss(const G4ParticleDefinition *aParticle, const G4String &)
Definition: G4EnergyLossTables.cc:1066
G4EnergyLossTables::GetProperTimeTable
static const G4PhysicsTable * GetProperTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:180
G4EnergyLossTables::lastParticle
static G4ParticleDefinition * lastParticle
Definition: G4EnergyLossTables.hh:233
G4EnergyLossTables::Chargesquare
static G4double Chargesquare
Definition: G4EnergyLossTables.hh:235
G4EnergyLossTables::t
static G4EnergyLossTablesHelper * t
Definition: G4EnergyLossTables.hh:231
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:256
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:113
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62
Generated by doxygen 1.9.3