Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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
53G4EnergyLossTablesHelper G4EnergyLossTables::t ;
54G4EnergyLossTablesHelper G4EnergyLossTables::null_loss ;
55const G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0;
56G4double G4EnergyLossTables::QQPositron = eplus*eplus ;
57G4double G4EnergyLossTables::Chargesquare ;
58G4int G4EnergyLossTables::oldIndex = -1 ;
59G4double G4EnergyLossTables::rmin = 0. ;
60G4double G4EnergyLossTables::rmax = 0. ;
61G4double G4EnergyLossTables::Thigh = 0. ;
62G4int G4EnergyLossTables::let_counter = 0;
63G4int G4EnergyLossTables::let_max_num_warnings = 100;
64G4bool G4EnergyLossTables::first_loss = true;
65
66G4EnergyLossTables::helper_map G4EnergyLossTables::dict;
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70G4EnergyLossTablesHelper::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
93G4EnergyLossTablesHelper::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
103void 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
131const 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
141const 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
151const 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
161const 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
171const 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
181G4EnergyLossTablesHelper 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
196G4double 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
245G4double 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
293G4double 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
375G4double 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
423G4double 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
505G4double 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
558G4double 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
731G4double 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
782G4double 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
836G4double 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
905G4double 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
952G4double 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
1007void 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
1032void
1033G4EnergyLossTables::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......
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4EnergyLossTables::GetRange
static G4double GetRange(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:505
G4EnergyLossTables::GetLabTimeTable
static const G4PhysicsTable * GetLabTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:161
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:103
G4EnergyLossTables::GetDEDX
static G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:196
G4EnergyLossTables::GetDeltaProperTime
static G4double GetDeltaProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:423
G4EnergyLossTables::GetLabTime
static G4double GetLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:245
G4EnergyLossTables::GetPreciseDEDX
static G4double GetPreciseDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:625
G4EnergyLossTables::GetPreciseEnergyFromRange
static G4double GetPreciseEnergyFromRange(const G4ParticleDefinition *aParticle, G4double range, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:558
G4EnergyLossTables::GetDEDXTable
static const G4PhysicsTable * GetDEDXTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:131
G4EnergyLossTables::GetRangeTable
static const G4PhysicsTable * GetRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:141
G4EnergyLossTables::GetInverseRangeTable
static const G4PhysicsTable * GetInverseRangeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:151
G4EnergyLossTables::GetDeltaLabTime
static G4double GetDeltaLabTime(const G4ParticleDefinition *aParticle, G4double KineticEnergyStart, G4double KineticEnergyEnd, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:293
G4EnergyLossTables::GetPreciseRangeFromEnergy
static G4double GetPreciseRangeFromEnergy(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:674
G4EnergyLossTables::GetProperTime
static G4double GetProperTime(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
Definition: G4EnergyLossTables.cc:375
G4EnergyLossTables::GetProperTimeTable
static const G4PhysicsTable * GetProperTimeTable(const G4ParticleDefinition *p)
Definition: G4EnergyLossTables.cc:171
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:107
G4LossTableManager::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *aParticle, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.cc:1212
G4LossTableManager::GetRange
G4double GetRange(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.cc:1198
G4LossTableManager::GetDEDX
G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.cc:1145
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:261
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:119
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83