00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "G4EnergyLossTables.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "G4MaterialCutsCouple.hh"
00047 #include "G4RegionStore.hh"
00048 #include "G4LossTableManager.hh"
00049
00050
00051
00052
00053 G4EnergyLossTablesHelper G4EnergyLossTables::t ;
00054 G4EnergyLossTablesHelper G4EnergyLossTables::null_loss ;
00055 const G4ParticleDefinition* G4EnergyLossTables::lastParticle = 0;
00056 G4double G4EnergyLossTables::QQPositron = eplus*eplus ;
00057 G4double G4EnergyLossTables::Chargesquare ;
00058 G4int G4EnergyLossTables::oldIndex = -1 ;
00059 G4double G4EnergyLossTables::rmin = 0. ;
00060 G4double G4EnergyLossTables::rmax = 0. ;
00061 G4double G4EnergyLossTables::Thigh = 0. ;
00062 G4int G4EnergyLossTables::let_counter = 0;
00063 G4int G4EnergyLossTables::let_max_num_warnings = 100;
00064 G4bool G4EnergyLossTables::first_loss = true;
00065
00066 G4EnergyLossTables::helper_map G4EnergyLossTables::dict;
00067
00068
00069
00070 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper(
00071 const G4PhysicsTable* aDEDXTable,
00072 const G4PhysicsTable* aRangeTable,
00073 const G4PhysicsTable* anInverseRangeTable,
00074 const G4PhysicsTable* aLabTimeTable,
00075 const G4PhysicsTable* aProperTimeTable,
00076 G4double aLowestKineticEnergy,
00077 G4double aHighestKineticEnergy,
00078 G4double aMassRatio,
00079 G4int aNumberOfBins)
00080 :
00081 theDEDXTable(aDEDXTable), theRangeTable(aRangeTable),
00082 theInverseRangeTable(anInverseRangeTable),
00083 theLabTimeTable(aLabTimeTable),
00084 theProperTimeTable(aProperTimeTable),
00085 theLowestKineticEnergy(aLowestKineticEnergy),
00086 theHighestKineticEnergy(aHighestKineticEnergy),
00087 theMassRatio(aMassRatio),
00088 theNumberOfBins(aNumberOfBins)
00089 { }
00090
00091
00092
00093 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper()
00094 {
00095 theLowestKineticEnergy = 0.0;
00096 theHighestKineticEnergy= 0.0;
00097 theMassRatio = 0.0;
00098 theNumberOfBins = 0;
00099 }
00100
00101
00102
00103 void G4EnergyLossTables::Register(
00104 const G4ParticleDefinition* p,
00105 const G4PhysicsTable* tDEDX,
00106 const G4PhysicsTable* tRange,
00107 const G4PhysicsTable* tInverseRange,
00108 const G4PhysicsTable* tLabTime,
00109 const G4PhysicsTable* tProperTime,
00110 G4double lowestKineticEnergy,
00111 G4double highestKineticEnergy,
00112 G4double massRatio,
00113 G4int NumberOfBins)
00114 {
00115 dict[p]= G4EnergyLossTablesHelper(tDEDX, tRange,tInverseRange,
00116 tLabTime,tProperTime,lowestKineticEnergy,
00117 highestKineticEnergy, massRatio,NumberOfBins);
00118
00119 t = GetTables(p) ;
00120 lastParticle = p ;
00121 Chargesquare = (p->GetPDGCharge())*(p->GetPDGCharge())/
00122 QQPositron ;
00123 if (first_loss ) {
00124 null_loss = G4EnergyLossTablesHelper(0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0);
00125 first_loss = false;
00126 }
00127 }
00128
00129
00130
00131 const G4PhysicsTable* G4EnergyLossTables::GetDEDXTable(
00132 const G4ParticleDefinition* p)
00133 {
00134 helper_map::iterator it;
00135 if((it=dict.find(p))==dict.end()) return 0;
00136 return (*it).second.theDEDXTable;
00137 }
00138
00139
00140
00141 const G4PhysicsTable* G4EnergyLossTables::GetRangeTable(
00142 const G4ParticleDefinition* p)
00143 {
00144 helper_map::iterator it;
00145 if((it=dict.find(p))==dict.end()) return 0;
00146 return (*it).second.theRangeTable;
00147 }
00148
00149
00150
00151 const G4PhysicsTable* G4EnergyLossTables::GetInverseRangeTable(
00152 const G4ParticleDefinition* p)
00153 {
00154 helper_map::iterator it;
00155 if((it=dict.find(p))==dict.end()) return 0;
00156 return (*it).second.theInverseRangeTable;
00157 }
00158
00159
00160
00161 const G4PhysicsTable* G4EnergyLossTables::GetLabTimeTable(
00162 const G4ParticleDefinition* p)
00163 {
00164 helper_map::iterator it;
00165 if((it=dict.find(p))==dict.end()) return 0;
00166 return (*it).second.theLabTimeTable;
00167 }
00168
00169
00170
00171 const G4PhysicsTable* G4EnergyLossTables::GetProperTimeTable(
00172 const G4ParticleDefinition* p)
00173 {
00174 helper_map::iterator it;
00175 if((it=dict.find(p))==dict.end()) return 0;
00176 return (*it).second.theProperTimeTable;
00177 }
00178
00179
00180
00181 G4EnergyLossTablesHelper G4EnergyLossTables::GetTables(
00182 const G4ParticleDefinition* p)
00183 {
00184 helper_map::iterator it;
00185 if ((it=dict.find(p))==dict.end()) {
00186
00187
00188
00189 return null_loss;
00190 }
00191 return (*it).second;
00192 }
00193
00194
00195
00196 G4double G4EnergyLossTables::GetDEDX(
00197 const G4ParticleDefinition *aParticle,
00198 G4double KineticEnergy,
00199 const G4Material *aMaterial)
00200 {
00201 CPRWarning();
00202 if(aParticle != lastParticle)
00203 {
00204 t= GetTables(aParticle);
00205 lastParticle = aParticle ;
00206 Chargesquare = (aParticle->GetPDGCharge())*
00207 (aParticle->GetPDGCharge())/
00208 QQPositron ;
00209 oldIndex = -1 ;
00210 }
00211 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00212 if (!dEdxTable) {
00213 ParticleHaveNoLoss(aParticle,"dEdx");
00214 return 0.0;
00215 }
00216
00217 G4int materialIndex = aMaterial->GetIndex();
00218 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00219 G4double dEdx;
00220 G4bool isOut;
00221
00222 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00223
00224 dEdx =(*dEdxTable)(materialIndex)->GetValue(
00225 t.theLowestKineticEnergy,isOut)
00226 *std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy);
00227
00228 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00229
00230 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00231 t.theHighestKineticEnergy,isOut);
00232
00233 } else {
00234
00235 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00236 scaledKineticEnergy,isOut);
00237
00238 }
00239
00240 return dEdx*Chargesquare;
00241 }
00242
00243
00244
00245 G4double G4EnergyLossTables::GetLabTime(
00246 const G4ParticleDefinition *aParticle,
00247 G4double KineticEnergy,
00248 const G4Material *aMaterial)
00249 {
00250 CPRWarning();
00251 if(aParticle != lastParticle)
00252 {
00253 t= GetTables(aParticle);
00254 lastParticle = aParticle ;
00255 oldIndex = -1 ;
00256 }
00257 const G4PhysicsTable* labtimeTable= t.theLabTimeTable;
00258 if (!labtimeTable) {
00259 ParticleHaveNoLoss(aParticle,"LabTime");
00260 return 0.0;
00261 }
00262
00263 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
00264 G4int materialIndex = aMaterial->GetIndex();
00265 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00266 G4double time;
00267 G4bool isOut;
00268
00269 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00270
00271 time = std::exp(ppar*std::log(scaledKineticEnergy/t.theLowestKineticEnergy))*
00272 (*labtimeTable)(materialIndex)->GetValue(
00273 t.theLowestKineticEnergy,isOut);
00274
00275
00276 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00277
00278 time = (*labtimeTable)(materialIndex)->GetValue(
00279 t.theHighestKineticEnergy,isOut);
00280
00281 } else {
00282
00283 time = (*labtimeTable)(materialIndex)->GetValue(
00284 scaledKineticEnergy,isOut);
00285
00286 }
00287
00288 return time/t.theMassRatio ;
00289 }
00290
00291
00292
00293 G4double G4EnergyLossTables::GetDeltaLabTime(
00294 const G4ParticleDefinition *aParticle,
00295 G4double KineticEnergyStart,
00296 G4double KineticEnergyEnd,
00297 const G4Material *aMaterial)
00298 {
00299 CPRWarning();
00300 if(aParticle != lastParticle)
00301 {
00302 t= GetTables(aParticle);
00303 lastParticle = aParticle ;
00304 oldIndex = -1 ;
00305 }
00306 const G4PhysicsTable* labtimeTable= t.theLabTimeTable;
00307 if (!labtimeTable) {
00308 ParticleHaveNoLoss(aParticle,"LabTime");
00309 return 0.0;
00310 }
00311
00312 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
00313 const G4double dToverT = 0.05 , facT = 1. -dToverT ;
00314 G4double timestart,timeend,deltatime,dTT;
00315 G4bool isOut;
00316
00317 G4int materialIndex = aMaterial->GetIndex();
00318 G4double scaledKineticEnergy = KineticEnergyStart*t.theMassRatio;
00319
00320 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00321
00322 timestart = std::exp(ppar*std::log(scaledKineticEnergy/t.theLowestKineticEnergy))*
00323 (*labtimeTable)(materialIndex)->GetValue(
00324 t.theLowestKineticEnergy,isOut);
00325
00326
00327 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00328
00329 timestart = (*labtimeTable)(materialIndex)->GetValue(
00330 t.theHighestKineticEnergy,isOut);
00331
00332 } else {
00333
00334 timestart = (*labtimeTable)(materialIndex)->GetValue(
00335 scaledKineticEnergy,isOut);
00336
00337 }
00338
00339 dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ;
00340
00341 if( dTT < dToverT )
00342 scaledKineticEnergy = facT*KineticEnergyStart*t.theMassRatio;
00343 else
00344 scaledKineticEnergy = KineticEnergyEnd*t.theMassRatio;
00345
00346 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00347
00348 timeend = std::exp(ppar*std::log(scaledKineticEnergy/t.theLowestKineticEnergy))*
00349 (*labtimeTable)(materialIndex)->GetValue(
00350 t.theLowestKineticEnergy,isOut);
00351
00352
00353 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00354
00355 timeend = (*labtimeTable)(materialIndex)->GetValue(
00356 t.theHighestKineticEnergy,isOut);
00357
00358 } else {
00359
00360 timeend = (*labtimeTable)(materialIndex)->GetValue(
00361 scaledKineticEnergy,isOut);
00362
00363 }
00364
00365 deltatime = timestart - timeend ;
00366
00367 if( dTT < dToverT )
00368 deltatime *= dTT/dToverT;
00369
00370 return deltatime/t.theMassRatio ;
00371 }
00372
00373
00374
00375 G4double G4EnergyLossTables::GetProperTime(
00376 const G4ParticleDefinition *aParticle,
00377 G4double KineticEnergy,
00378 const G4Material *aMaterial)
00379 {
00380 CPRWarning();
00381 if(aParticle != lastParticle)
00382 {
00383 t= GetTables(aParticle);
00384 lastParticle = aParticle ;
00385 oldIndex = -1 ;
00386 }
00387 const G4PhysicsTable* propertimeTable= t.theProperTimeTable;
00388 if (!propertimeTable) {
00389 ParticleHaveNoLoss(aParticle,"ProperTime");
00390 return 0.0;
00391 }
00392
00393 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
00394 G4int materialIndex = aMaterial->GetIndex();
00395 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00396 G4double time;
00397 G4bool isOut;
00398
00399 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00400
00401 time = std::exp(ppar*std::log(scaledKineticEnergy/t.theLowestKineticEnergy))*
00402 (*propertimeTable)(materialIndex)->GetValue(
00403 t.theLowestKineticEnergy,isOut);
00404
00405
00406 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00407
00408 time = (*propertimeTable)(materialIndex)->GetValue(
00409 t.theHighestKineticEnergy,isOut);
00410
00411 } else {
00412
00413 time = (*propertimeTable)(materialIndex)->GetValue(
00414 scaledKineticEnergy,isOut);
00415
00416 }
00417
00418 return time/t.theMassRatio ;
00419 }
00420
00421
00422
00423 G4double G4EnergyLossTables::GetDeltaProperTime(
00424 const G4ParticleDefinition *aParticle,
00425 G4double KineticEnergyStart,
00426 G4double KineticEnergyEnd,
00427 const G4Material *aMaterial)
00428 {
00429 CPRWarning();
00430 if(aParticle != lastParticle)
00431 {
00432 t= GetTables(aParticle);
00433 lastParticle = aParticle ;
00434 oldIndex = -1 ;
00435 }
00436 const G4PhysicsTable* propertimeTable= t.theProperTimeTable;
00437 if (!propertimeTable) {
00438 ParticleHaveNoLoss(aParticle,"ProperTime");
00439 return 0.0;
00440 }
00441
00442 const G4double parlowen=0.4 , ppar=0.5-parlowen ;
00443 const G4double dToverT = 0.05 , facT = 1. -dToverT ;
00444 G4double timestart,timeend,deltatime,dTT;
00445 G4bool isOut;
00446
00447 G4int materialIndex = aMaterial->GetIndex();
00448 G4double scaledKineticEnergy = KineticEnergyStart*t.theMassRatio;
00449
00450 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00451
00452 timestart = std::exp(ppar*std::log(scaledKineticEnergy/t.theLowestKineticEnergy))*
00453 (*propertimeTable)(materialIndex)->GetValue(
00454 t.theLowestKineticEnergy,isOut);
00455
00456
00457 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00458
00459 timestart = (*propertimeTable)(materialIndex)->GetValue(
00460 t.theHighestKineticEnergy,isOut);
00461
00462 } else {
00463
00464 timestart = (*propertimeTable)(materialIndex)->GetValue(
00465 scaledKineticEnergy,isOut);
00466
00467 }
00468
00469 dTT = (KineticEnergyStart - KineticEnergyEnd)/KineticEnergyStart ;
00470
00471 if( dTT < dToverT )
00472 scaledKineticEnergy = facT*KineticEnergyStart*t.theMassRatio;
00473 else
00474 scaledKineticEnergy = KineticEnergyEnd*t.theMassRatio;
00475
00476 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00477
00478 timeend = std::exp(ppar*std::log(scaledKineticEnergy/t.theLowestKineticEnergy))*
00479 (*propertimeTable)(materialIndex)->GetValue(
00480 t.theLowestKineticEnergy,isOut);
00481
00482
00483 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00484
00485 timeend = (*propertimeTable)(materialIndex)->GetValue(
00486 t.theHighestKineticEnergy,isOut);
00487
00488 } else {
00489
00490 timeend = (*propertimeTable)(materialIndex)->GetValue(
00491 scaledKineticEnergy,isOut);
00492
00493 }
00494
00495 deltatime = timestart - timeend ;
00496
00497 if( dTT < dToverT )
00498 deltatime *= dTT/dToverT ;
00499
00500 return deltatime/t.theMassRatio ;
00501 }
00502
00503
00504
00505 G4double G4EnergyLossTables::GetRange(
00506 const G4ParticleDefinition *aParticle,
00507 G4double KineticEnergy,
00508 const G4Material *aMaterial)
00509 {
00510 CPRWarning();
00511 if(aParticle != lastParticle)
00512 {
00513 t= GetTables(aParticle);
00514 lastParticle = aParticle ;
00515 Chargesquare = (aParticle->GetPDGCharge())*
00516 (aParticle->GetPDGCharge())/
00517 QQPositron ;
00518 oldIndex = -1 ;
00519 }
00520 const G4PhysicsTable* rangeTable= t.theRangeTable;
00521 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00522 if (!rangeTable) {
00523 ParticleHaveNoLoss(aParticle,"Range");
00524 return 0.0;
00525 }
00526
00527 G4int materialIndex = aMaterial->GetIndex();
00528 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00529 G4double Range;
00530 G4bool isOut;
00531
00532 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00533
00534 Range = std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)*
00535 (*rangeTable)(materialIndex)->GetValue(
00536 t.theLowestKineticEnergy,isOut);
00537
00538 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00539
00540 Range = (*rangeTable)(materialIndex)->GetValue(
00541 t.theHighestKineticEnergy,isOut)+
00542 (scaledKineticEnergy-t.theHighestKineticEnergy)/
00543 (*dEdxTable)(materialIndex)->GetValue(
00544 t.theHighestKineticEnergy,isOut);
00545
00546 } else {
00547
00548 Range = (*rangeTable)(materialIndex)->GetValue(
00549 scaledKineticEnergy,isOut);
00550
00551 }
00552
00553 return Range/(Chargesquare*t.theMassRatio);
00554 }
00555
00556
00557
00558 G4double G4EnergyLossTables::GetPreciseEnergyFromRange(
00559 const G4ParticleDefinition *aParticle,
00560 G4double range,
00561 const G4Material *aMaterial)
00562
00563 {
00564 CPRWarning();
00565 if( aParticle != lastParticle)
00566 {
00567 t= GetTables(aParticle);
00568 lastParticle = aParticle;
00569 Chargesquare = (aParticle->GetPDGCharge())*
00570 (aParticle->GetPDGCharge())/
00571 QQPositron ;
00572 oldIndex = -1 ;
00573 }
00574 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00575 const G4PhysicsTable* inverseRangeTable= t.theInverseRangeTable;
00576 if (!inverseRangeTable) {
00577 ParticleHaveNoLoss(aParticle,"InverseRange");
00578 return 0.0;
00579 }
00580
00581 G4double scaledrange,scaledKineticEnergy ;
00582 G4bool isOut ;
00583
00584 G4int materialIndex = aMaterial->GetIndex() ;
00585
00586 if(materialIndex != oldIndex)
00587 {
00588 oldIndex = materialIndex ;
00589 rmin = (*inverseRangeTable)(materialIndex)->
00590 GetLowEdgeEnergy(0) ;
00591 rmax = (*inverseRangeTable)(materialIndex)->
00592 GetLowEdgeEnergy(t.theNumberOfBins-2) ;
00593 Thigh = (*inverseRangeTable)(materialIndex)->
00594 GetValue(rmax,isOut) ;
00595 }
00596
00597 scaledrange = range*Chargesquare*t.theMassRatio ;
00598
00599 if(scaledrange < rmin)
00600 {
00601 scaledKineticEnergy = t.theLowestKineticEnergy*
00602 scaledrange*scaledrange/(rmin*rmin) ;
00603 }
00604 else
00605 {
00606 if(scaledrange < rmax)
00607 {
00608 scaledKineticEnergy = (*inverseRangeTable)(materialIndex)->
00609 GetValue( scaledrange,isOut) ;
00610 }
00611 else
00612 {
00613 scaledKineticEnergy = Thigh +
00614 (scaledrange-rmax)*
00615 (*dEdxTable)(materialIndex)->
00616 GetValue(Thigh,isOut) ;
00617 }
00618 }
00619
00620 return scaledKineticEnergy/t.theMassRatio ;
00621 }
00622
00623
00624
00625 G4double G4EnergyLossTables::GetPreciseDEDX(
00626 const G4ParticleDefinition *aParticle,
00627 G4double KineticEnergy,
00628 const G4Material *aMaterial)
00629 {
00630 CPRWarning();
00631 if( aParticle != lastParticle)
00632 {
00633 t= GetTables(aParticle);
00634 lastParticle = aParticle;
00635 Chargesquare = (aParticle->GetPDGCharge())*
00636 (aParticle->GetPDGCharge())/
00637 QQPositron ;
00638 oldIndex = -1 ;
00639 }
00640 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00641 if (!dEdxTable) {
00642 ParticleHaveNoLoss(aParticle,"dEdx");
00643 return 0.0;
00644 }
00645
00646 G4int materialIndex = aMaterial->GetIndex();
00647 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00648 G4double dEdx;
00649 G4bool isOut;
00650
00651 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00652
00653 dEdx = std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)
00654 *(*dEdxTable)(materialIndex)->GetValue(
00655 t.theLowestKineticEnergy,isOut);
00656
00657 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00658
00659 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00660 t.theHighestKineticEnergy,isOut);
00661
00662 } else {
00663
00664 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00665 scaledKineticEnergy,isOut) ;
00666
00667 }
00668
00669 return dEdx*Chargesquare;
00670 }
00671
00672
00673
00674 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy(
00675 const G4ParticleDefinition *aParticle,
00676 G4double KineticEnergy,
00677 const G4Material *aMaterial)
00678 {
00679 CPRWarning();
00680 if( aParticle != lastParticle)
00681 {
00682 t= GetTables(aParticle);
00683 lastParticle = aParticle;
00684 Chargesquare = (aParticle->GetPDGCharge())*
00685 (aParticle->GetPDGCharge())/
00686 QQPositron ;
00687 oldIndex = -1 ;
00688 }
00689 const G4PhysicsTable* rangeTable= t.theRangeTable;
00690 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00691 if (!rangeTable) {
00692 ParticleHaveNoLoss(aParticle,"Range");
00693 return 0.0;
00694 }
00695 G4int materialIndex = aMaterial->GetIndex();
00696
00697 G4double Thighr = t.theHighestKineticEnergy*t.theLowestKineticEnergy/
00698 (*rangeTable)(materialIndex)->
00699 GetLowEdgeEnergy(1) ;
00700
00701 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00702 G4double Range;
00703 G4bool isOut;
00704
00705 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00706
00707 Range = std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)*
00708 (*rangeTable)(materialIndex)->GetValue(
00709 t.theLowestKineticEnergy,isOut);
00710
00711 } else if (scaledKineticEnergy>Thighr) {
00712
00713 Range = (*rangeTable)(materialIndex)->GetValue(
00714 Thighr,isOut)+
00715 (scaledKineticEnergy-Thighr)/
00716 (*dEdxTable)(materialIndex)->GetValue(
00717 Thighr,isOut);
00718
00719 } else {
00720
00721 Range = (*rangeTable)(materialIndex)->GetValue(
00722 scaledKineticEnergy,isOut) ;
00723
00724 }
00725
00726 return Range/(Chargesquare*t.theMassRatio);
00727 }
00728
00729
00730
00731 G4double G4EnergyLossTables::GetDEDX(
00732 const G4ParticleDefinition *aParticle,
00733 G4double KineticEnergy,
00734 const G4MaterialCutsCouple *couple,
00735 G4bool check)
00736 {
00737 if(aParticle != lastParticle)
00738 {
00739 t= GetTables(aParticle);
00740 lastParticle = aParticle ;
00741 Chargesquare = (aParticle->GetPDGCharge())*
00742 (aParticle->GetPDGCharge())/
00743 QQPositron ;
00744 oldIndex = -1 ;
00745 }
00746 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00747
00748 if (!dEdxTable ) {
00749 if (check) return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple);
00750 else ParticleHaveNoLoss(aParticle, "dEdx");
00751 return 0.0;
00752 }
00753
00754 G4int materialIndex = couple->GetIndex();
00755 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00756 G4double dEdx;
00757 G4bool isOut;
00758
00759 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00760
00761 dEdx =(*dEdxTable)(materialIndex)->GetValue(
00762 t.theLowestKineticEnergy,isOut)
00763 *std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy);
00764
00765 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00766
00767 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00768 t.theHighestKineticEnergy,isOut);
00769
00770 } else {
00771
00772 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00773 scaledKineticEnergy,isOut);
00774
00775 }
00776
00777 return dEdx*Chargesquare;
00778 }
00779
00780
00781
00782 G4double G4EnergyLossTables::GetRange(
00783 const G4ParticleDefinition *aParticle,
00784 G4double KineticEnergy,
00785 const G4MaterialCutsCouple *couple,
00786 G4bool check)
00787 {
00788 if(aParticle != lastParticle)
00789 {
00790 t= GetTables(aParticle);
00791 lastParticle = aParticle ;
00792 Chargesquare = (aParticle->GetPDGCharge())*
00793 (aParticle->GetPDGCharge())/
00794 QQPositron ;
00795 oldIndex = -1 ;
00796 }
00797 const G4PhysicsTable* rangeTable= t.theRangeTable;
00798 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00799 if (!rangeTable) {
00800 if(check) return G4LossTableManager::Instance()->GetRange(aParticle,KineticEnergy,couple);
00801 else return DBL_MAX;
00802
00803 }
00804
00805 G4int materialIndex = couple->GetIndex();
00806 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00807 G4double Range;
00808 G4bool isOut;
00809
00810 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00811
00812 Range = std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)*
00813 (*rangeTable)(materialIndex)->GetValue(
00814 t.theLowestKineticEnergy,isOut);
00815
00816 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00817
00818 Range = (*rangeTable)(materialIndex)->GetValue(
00819 t.theHighestKineticEnergy,isOut)+
00820 (scaledKineticEnergy-t.theHighestKineticEnergy)/
00821 (*dEdxTable)(materialIndex)->GetValue(
00822 t.theHighestKineticEnergy,isOut);
00823
00824 } else {
00825
00826 Range = (*rangeTable)(materialIndex)->GetValue(
00827 scaledKineticEnergy,isOut);
00828
00829 }
00830
00831 return Range/(Chargesquare*t.theMassRatio);
00832 }
00833
00834
00835
00836 G4double G4EnergyLossTables::GetPreciseEnergyFromRange(
00837 const G4ParticleDefinition *aParticle,
00838 G4double range,
00839 const G4MaterialCutsCouple *couple,
00840 G4bool check)
00841
00842 {
00843 if( aParticle != lastParticle)
00844 {
00845 t= GetTables(aParticle);
00846 lastParticle = aParticle;
00847 Chargesquare = (aParticle->GetPDGCharge())*
00848 (aParticle->GetPDGCharge())/
00849 QQPositron ;
00850 oldIndex = -1 ;
00851 }
00852 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00853 const G4PhysicsTable* inverseRangeTable= t.theInverseRangeTable;
00854
00855 if (!inverseRangeTable) {
00856 if(check) return G4LossTableManager::Instance()->GetEnergy(aParticle,range,couple);
00857 else return DBL_MAX;
00858
00859 }
00860
00861 G4double scaledrange,scaledKineticEnergy ;
00862 G4bool isOut ;
00863
00864 G4int materialIndex = couple->GetIndex() ;
00865
00866 if(materialIndex != oldIndex)
00867 {
00868 oldIndex = materialIndex ;
00869 rmin = (*inverseRangeTable)(materialIndex)->
00870 GetLowEdgeEnergy(0) ;
00871 rmax = (*inverseRangeTable)(materialIndex)->
00872 GetLowEdgeEnergy(t.theNumberOfBins-2) ;
00873 Thigh = (*inverseRangeTable)(materialIndex)->
00874 GetValue(rmax,isOut) ;
00875 }
00876
00877 scaledrange = range*Chargesquare*t.theMassRatio ;
00878
00879 if(scaledrange < rmin)
00880 {
00881 scaledKineticEnergy = t.theLowestKineticEnergy*
00882 scaledrange*scaledrange/(rmin*rmin) ;
00883 }
00884 else
00885 {
00886 if(scaledrange < rmax)
00887 {
00888 scaledKineticEnergy = (*inverseRangeTable)(materialIndex)->
00889 GetValue( scaledrange,isOut) ;
00890 }
00891 else
00892 {
00893 scaledKineticEnergy = Thigh +
00894 (scaledrange-rmax)*
00895 (*dEdxTable)(materialIndex)->
00896 GetValue(Thigh,isOut) ;
00897 }
00898 }
00899
00900 return scaledKineticEnergy/t.theMassRatio ;
00901 }
00902
00903
00904
00905 G4double G4EnergyLossTables::GetPreciseDEDX(
00906 const G4ParticleDefinition *aParticle,
00907 G4double KineticEnergy,
00908 const G4MaterialCutsCouple *couple)
00909 {
00910
00911 if( aParticle != lastParticle)
00912 {
00913 t= GetTables(aParticle);
00914 lastParticle = aParticle;
00915 Chargesquare = (aParticle->GetPDGCharge())*
00916 (aParticle->GetPDGCharge())/
00917 QQPositron ;
00918 oldIndex = -1 ;
00919 }
00920 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00921 if ( !dEdxTable )
00922 return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple);
00923
00924 G4int materialIndex = couple->GetIndex();
00925 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00926 G4double dEdx;
00927 G4bool isOut;
00928
00929 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00930
00931 dEdx = std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)
00932 *(*dEdxTable)(materialIndex)->GetValue(
00933 t.theLowestKineticEnergy,isOut);
00934
00935 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) {
00936
00937 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00938 t.theHighestKineticEnergy,isOut);
00939
00940 } else {
00941
00942 dEdx = (*dEdxTable)(materialIndex)->GetValue(
00943 scaledKineticEnergy,isOut) ;
00944
00945 }
00946
00947 return dEdx*Chargesquare;
00948 }
00949
00950
00951
00952 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy(
00953 const G4ParticleDefinition *aParticle,
00954 G4double KineticEnergy,
00955 const G4MaterialCutsCouple *couple)
00956 {
00957 if( aParticle != lastParticle)
00958 {
00959 t= GetTables(aParticle);
00960 lastParticle = aParticle;
00961 Chargesquare = (aParticle->GetPDGCharge())*
00962 (aParticle->GetPDGCharge())/
00963 QQPositron ;
00964 oldIndex = -1 ;
00965 }
00966 const G4PhysicsTable* rangeTable= t.theRangeTable;
00967 const G4PhysicsTable* dEdxTable= t.theDEDXTable;
00968 if ( !dEdxTable || !rangeTable)
00969 return G4LossTableManager::Instance()->GetDEDX(aParticle,KineticEnergy,couple);
00970
00971 G4int materialIndex = couple->GetIndex();
00972
00973 G4double Thighr = t.theHighestKineticEnergy*t.theLowestKineticEnergy/
00974 (*rangeTable)(materialIndex)->
00975 GetLowEdgeEnergy(1) ;
00976
00977 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio;
00978 G4double Range;
00979 G4bool isOut;
00980
00981 if (scaledKineticEnergy<t.theLowestKineticEnergy) {
00982
00983 Range = std::sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)*
00984 (*rangeTable)(materialIndex)->GetValue(
00985 t.theLowestKineticEnergy,isOut);
00986
00987 } else if (scaledKineticEnergy>Thighr) {
00988
00989 Range = (*rangeTable)(materialIndex)->GetValue(
00990 Thighr,isOut)+
00991 (scaledKineticEnergy-Thighr)/
00992 (*dEdxTable)(materialIndex)->GetValue(
00993 Thighr,isOut);
00994
00995 } else {
00996
00997 Range = (*rangeTable)(materialIndex)->GetValue(
00998 scaledKineticEnergy,isOut) ;
00999
01000 }
01001
01002 return Range/(Chargesquare*t.theMassRatio);
01003 }
01004
01005
01006
01007 void G4EnergyLossTables::CPRWarning()
01008 {
01009 if (let_counter < let_max_num_warnings) {
01010
01011 G4cout << G4endl;
01012 G4cout << "##### G4EnergyLossTable WARNING: The obsolete interface is used!" << G4endl;
01013 G4cout << "##### RESULTS ARE NOT GARANTEED!" << G4endl;
01014 G4cout << "##### Please, substitute G4Material by G4MaterialCutsCouple" << G4endl;
01015 G4cout << "##### Obsolete interface will be removed soon" << G4endl;
01016 G4cout << G4endl;
01017 let_counter++;
01018
01019
01020
01021
01022
01023 } else if (let_counter == let_max_num_warnings) {
01024
01025 G4cout << "##### G4EnergyLossTable WARNING closed" << G4endl;
01026 let_counter++;
01027 }
01028 }
01029
01030
01031
01032 void
01033 G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition*,
01034 const G4String& )
01035 {
01036
01037
01038
01039
01040 }
01041
01042