| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // -------------------------------------------------------------------
29 // first version created by P.Urban , 06/04/1998
30 // modifications + "precise" functions added by L.Urban , 27/05/98
31 // modifications , TOF functions , 26/10/98, L.Urban
32 // cache mechanism in order to gain time, 11/02/99, L.Urban
33 // bug fixed , 12/04/99 , L.Urban
34 // 10.11.99: moved from RWT hash dictionary to STL map, G.Barrand, M.Maire
35 // 27.09.01 L.Urban , bug fixed (negative energy deposit)
36 // 26.10.01 all static functions moved from .icc files (mma)
37 // 15.01.03 Add interfaces required for "cut per region" (V.Ivanchenko)
38 // 12.03.03 Add warnings to obsolete interfaces (V.Ivanchenko)
39 // 10.04.03 Add call to G4LossTableManager is particle is not registered (V.Ivanchenko)
40 //
41 // -------------------------------------------------------------------
42
43 #include "G4EnergyLossTables.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4MaterialCutsCouple.hh"
46 #include "G4RegionStore.hh"
47 #include "G4LossTableManager.hh"
48
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52 G4EnergyLossTablesHelper *G4EnergyLossTables::t = nullptr;
53 G4EnergyLossTablesHelper *G4EnergyLossTables::null_loss = nullptr;
54 G4ParticleDefinition* G4EnergyLossTables::lastParticle = nullptr;
55 G4double G4EnergyLossTables::QQPositron = 1.0; // e_squared
56 G4double G4EnergyLossTables::Chargesquare ;
57 G4int G4EnergyLossTables::oldIndex = -1 ;
58 G4double G4EnergyLossTables::rmin = 0. ;
59 G4double G4EnergyLossTables::rmax = 0. ;
60 G4double G4EnergyLossTables::Thigh = 0. ;
61 G4int G4EnergyLossTables::let_counter = 0;
62 G4int G4EnergyLossTables::let_max_num_warnings = 100;
63 G4bool G4EnergyLossTables::first_loss = true;
64
65 G4EnergyLossTables::helper_map *G4EnergyLossTables::dict = nullptr;
66
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68
69 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper(
70 const G4PhysicsTable* aDEDXTable,
71 const G4PhysicsTable* aRangeTable,
72 const G4PhysicsTable* anInverseRangeTable,
73 const G4PhysicsTable* aLabTimeTable,
74 const G4PhysicsTable* aProperTimeTable,
75 G4double aLowestKineticEnergy,
76 G4double aHighestKineticEnergy,
77 G4double aMassRatio,
78 G4int aNumberOfBins)
79 :
80 theDEDXTable(aDEDXTable), theRangeTable(aRangeTable),
81 theInverseRangeTable(anInverseRangeTable),
82 theLabTimeTable(aLabTimeTable),
83 theProperTimeTable(aProperTimeTable),
84 theLowestKineticEnergy(aLowestKineticEnergy),
85 theHighestKineticEnergy(aHighestKineticEnergy),
86 theMassRatio(aMassRatio),
87 theNumberOfBins(aNumberOfBins)
88 { }
89
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91
92 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper()
93 {
94 theLowestKineticEnergy = 0.0;
95 theHighestKineticEnergy= 0.0;
96 theMassRatio = 0.0;
97 theNumberOfBins = 0;
98 theDEDXTable = theRangeTable = theInverseRangeTable = theLabTimeTable
99 = theProperTimeTable = nullptr;
100 }
101
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103
104 void G4EnergyLossTables::Register(
105 const G4ParticleDefinition* p,
106 const G4PhysicsTable* tDEDX,
107 const G4PhysicsTable* tRange,
108 const G4PhysicsTable* tInverseRange,
109 const G4PhysicsTable* tLabTime,
110 const G4PhysicsTable* tProperTime,
111 G4double lowestKineticEnergy,
112 G4double highestKineticEnergy,
113 G4double massRatio,
114 G4int NumberOfBins)
115 {
116 if (!dict) dict = new G4EnergyLossTables::helper_map;
117 if (!null_loss) null_loss = new G4EnergyLossTablesHelper;
118 if (!t) t = new G4EnergyLossTablesHelper;
119
120 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, tRange,tInverseRange,
121 tLabTime,tProperTime,lowestKineticEnergy,
122 highestKineticEnergy, massRatio,NumberOfBins);
123
124 *t = GetTables(p) ; // important for cache !!!!!
125 lastParticle = (G4ParticleDefinition*) p ;
126 Chargesquare = (p->GetPDGCharge())*(p->GetPDGCharge())/
127 QQPositron ;
128 if (first_loss ) {
129 *null_loss = G4EnergyLossTablesHelper(
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
137 const 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
148 const 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
159 const 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
170 const 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
181 const 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
192 G4EnergyLossTablesHelper 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
207 G4double 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
258 G4double 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
308 G4double 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
392 G4double 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
442 G4double 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
526 G4double 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
581 G4double 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
760 G4double 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
813 G4double 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
869 G4double 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
940 G4double 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
988 G4double 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
1045 void 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
1066 void
1067 G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition*,
1068 const G4String& /*q*/)
1069 {
1070 }
1071
1072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1073