G4EnergyLossTables.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // -------------------------------------------------------------------
00030 // first version created by P.Urban , 06/04/1998
00031 // modifications + "precise" functions added by L.Urban , 27/05/98
00032 // modifications , TOF functions , 26/10/98, L.Urban
00033 // cache mechanism in order to gain time, 11/02/99, L.Urban
00034 // bug fixed , 12/04/99 , L.Urban
00035 // 10.11.99: moved from RWT hash dictionary to STL map, G.Barrand, M.Maire
00036 // 27.09.01 L.Urban , bug fixed (negative energy deposit)
00037 // 26.10.01 all static functions moved from .icc files (mma)
00038 // 15.01.03 Add interfaces required for "cut per region" (V.Ivanchenko)
00039 // 12.03.03 Add warnings to obsolete interfaces (V.Ivanchenko)
00040 // 10.04.03 Add call to G4LossTableManager is particle is not registered (V.Ivanchenko)
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00092 
00093 G4EnergyLossTablesHelper::G4EnergyLossTablesHelper()
00094 { 
00095   theLowestKineticEnergy = 0.0;
00096   theHighestKineticEnergy= 0.0;
00097   theMassRatio = 0.0;
00098   theNumberOfBins = 0;
00099 }
00100 
00101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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) ;    // important for cache !!!!!
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00180 
00181 G4EnergyLossTablesHelper G4EnergyLossTables::GetTables(
00182   const G4ParticleDefinition* p)
00183 {
00184   helper_map::iterator it;
00185   if ((it=dict.find(p))==dict.end()) {
00186 //    G4cout << "Table is not found out for " << p->GetParticleName() << G4endl;
00187 //    G4Exception("G4EnergyLossTables::GetTables: table not found!");
00188 //    exit(1);
00189     return null_loss;
00190   }
00191   return (*it).second;
00192 }
00193 
00194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00557 
00558 G4double G4EnergyLossTables::GetPreciseEnergyFromRange(
00559                                      const G4ParticleDefinition *aParticle,
00560                                            G4double range,
00561                                      const G4Material *aMaterial)
00562 // it returns the value of the kinetic energy for a given range
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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       //ParticleHaveNoLoss(aParticle,"Range");
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00835 
00836 G4double G4EnergyLossTables::GetPreciseEnergyFromRange(
00837                                      const G4ParticleDefinition *aParticle,
00838                                            G4double range,
00839                                      const G4MaterialCutsCouple *couple,
00840                                            G4bool check)
00841 // it returns the value of the kinetic energy for a given range
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     //    else      ParticleHaveNoLoss(aParticle,"InverseRange");
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //  if ((G4RegionStore::GetInstance())->size() > 1) {
01019 //     G4Exception("G4EnergyLossTables:: More than 1 region - table can't be accessed with obsolete interface");
01020 //     exit(1);
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01031 
01032 void 
01033 G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition*, 
01034                                        const G4String& /*q*/)
01035 {
01036   //G4String s = " " + q + " table not found for "
01037   //           + aParticle->GetParticleName() + " !";
01038   //G4Exception("G4EnergyLossTables::ParticleHaveNoLoss", "EM01",
01039   //          FatalException, s);
01040 }
01041 
01042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Generated on Mon May 27 17:48:10 2013 for Geant4 by  doxygen 1.4.7