G4VRangeToEnergyConverter.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 // --------------------------------------------------------------
00031 //      GEANT 4 class implementation file/  History:
00032 //    5 Oct. 2002, H.Kuirashige : Structure created based on object model
00033 // --------------------------------------------------------------
00034 
00035 #include "G4VRangeToEnergyConverter.hh"
00036 #include "G4ParticleTable.hh"
00037 #include "G4Material.hh"
00038 #include "G4MaterialTable.hh"
00039 #include "G4PhysicsLogVector.hh"
00040 
00041 #include "G4ios.hh"
00042 #include "G4SystemOfUnits.hh"
00043 
00044 // energy range
00045 G4double  G4VRangeToEnergyConverter::LowestEnergy = 0.99e-3*MeV;
00046 G4double  G4VRangeToEnergyConverter::HighestEnergy = 100.0e6*MeV;
00047 
00048 // max energy cut
00049 G4double  G4VRangeToEnergyConverter::MaxEnergyCut = 10.0*GeV;
00050 
00051 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter():
00052   theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(300),
00053   verboseLevel(1)
00054 {
00055   fMaxEnergyCut = 0.;
00056 }
00057 
00058 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right) :  theParticle(right.theParticle), theLossTable(0), TotBin(right.TotBin)
00059 {
00060   fMaxEnergyCut = 0.;
00061   *this = right;
00062 }
00063 
00064 G4VRangeToEnergyConverter & G4VRangeToEnergyConverter::operator=(const G4VRangeToEnergyConverter &right)
00065 {
00066   if (this == &right) return *this;
00067   if (theLossTable) {
00068     theLossTable->clearAndDestroy();
00069     delete theLossTable;
00070     theLossTable=0;
00071  }
00072 
00073   NumberOfElements = right.NumberOfElements;
00074   theParticle = right.theParticle;
00075   verboseLevel = right.verboseLevel;
00076   
00077   // create the loss table
00078   theLossTable = new G4LossTable();
00079   theLossTable->reserve(G4Element::GetNumberOfElements());  
00080   // fill the loss table
00081   for (size_t j=0; j<size_t(NumberOfElements); j++){
00082     G4LossVector* aVector= new
00083             G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
00084     for (size_t i=0; i<size_t(TotBin); i++) {
00085       G4double Value = (*((*right.theLossTable)[j]))[i];
00086       aVector->PutValue(i,Value);
00087     }
00088     theLossTable->insert(aVector);
00089   }
00090 
00091   // clean up range vector store
00092   for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
00093     delete fRangeVectorStore.at(idx);
00094   }
00095   fRangeVectorStore.clear();
00096 
00097   // copy range vector store
00098   for (size_t j=0; j<((right.fRangeVectorStore).size()); j++){
00099     G4RangeVector* vector = (right.fRangeVectorStore).at(j);
00100     G4RangeVector* rangeVector = 0; 
00101     if (vector !=0 ) {
00102       rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
00103       fMaxEnergyCut = MaxEnergyCut;   
00104       for (size_t i=0; i<size_t(TotBin); i++) {
00105         G4double Value = (*vector)[i];
00106         rangeVector->PutValue(i,Value);
00107       }
00108     }
00109     fRangeVectorStore.push_back(rangeVector);
00110   }
00111   return *this;
00112 }
00113 
00114 
00115 G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
00116 { 
00117   Reset(); 
00118 }
00119 
00120 G4int G4VRangeToEnergyConverter::operator==(const G4VRangeToEnergyConverter &right) const
00121 {
00122   return this == &right;
00123 }
00124 
00125 G4int G4VRangeToEnergyConverter::operator!=(const G4VRangeToEnergyConverter &right) const
00126 {
00127   return this != &right;
00128 }
00129 
00130 
00131 // **********************************************************************
00132 // ************************* Convert  ***********************************
00133 // **********************************************************************
00134 G4double G4VRangeToEnergyConverter::Convert(G4double rangeCut, 
00135                                             const G4Material* material) 
00136 {
00137 #ifdef G4VERBOSE
00138     if (GetVerboseLevel()>3) {
00139       G4cout << "G4VRangeToEnergyConverter::Convert() ";
00140       G4cout << "Convert for " << material->GetName() 
00141              << " with Range Cut " << rangeCut/mm << "[mm]" << G4endl;
00142     }
00143 #endif
00144 
00145   G4double theKineticEnergyCuts = 0.;
00146 
00147   if (fMaxEnergyCut != MaxEnergyCut) {
00148     fMaxEnergyCut = MaxEnergyCut;      
00149     // clear loss table and renge vectors
00150     Reset();
00151   }
00152  
00153   // Build the energy loss table
00154   BuildLossTable();
00155   
00156   // Build range vector for every material, convert cut into energy-cut,
00157   // fill theKineticEnergyCuts and delete the range vector
00158   G4double tune = 0.025*mm*g/cm3 ,lowen = 30.*keV ; 
00159 
00160   // check density
00161   G4double density = material->GetDensity() ;
00162   if(density <= 0.) {
00163  #ifdef G4VERBOSE
00164     if (GetVerboseLevel()>0) {
00165       G4cout << "G4VRangeToEnergyConverter::Convert() ";
00166       G4cout << material->GetName() << "has zero density "
00167              << "( " << density << ")" << G4endl;
00168     }
00169 #endif
00170     return 0.;
00171   }
00172  
00173    // initialize RangeVectorStore
00174   const G4MaterialTable* table = G4Material::GetMaterialTable();
00175   G4int ext_size = table->size() - fRangeVectorStore.size();
00176   for (int i=0; i<ext_size; i++) fRangeVectorStore.push_back(0);
00177   
00178   // Build Range Vector
00179   G4int idx = material->GetIndex(); 
00180   G4RangeVector* rangeVector = fRangeVectorStore.at(idx);
00181   if (rangeVector == 0) {
00182     rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin); 
00183     BuildRangeVector(material, rangeVector);
00184     fRangeVectorStore.at(idx) = rangeVector;
00185   }
00186 
00187   // Convert Range Cut ro Kinetic Energy Cut 
00188   theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
00189   
00190   if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
00191       && (theKineticEnergyCuts < lowen) ) {
00192     //  corr. should be switched on smoothly   
00193     theKineticEnergyCuts /= (1.+(1.-theKineticEnergyCuts/lowen)*
00194                              tune/(rangeCut*density)); 
00195   }
00196   
00197   if(theKineticEnergyCuts < LowestEnergy) {
00198     theKineticEnergyCuts = LowestEnergy ;
00199   } else if(theKineticEnergyCuts > MaxEnergyCut) {
00200     theKineticEnergyCuts = MaxEnergyCut;
00201   }
00202   
00203   return theKineticEnergyCuts;
00204 }
00205 
00206 // **********************************************************************
00207 // ************************ SetEnergyRange  *****************************
00208 // **********************************************************************
00209 void G4VRangeToEnergyConverter::SetEnergyRange(G4double lowedge, 
00210                                                G4double highedge)
00211 {
00212   // check LowestEnergy/ HighestEnergy 
00213   if ( (lowedge<0.0)||(highedge<=lowedge) ){
00214 #ifdef G4VERBOSE
00215     G4cerr << "Error in G4VRangeToEnergyConverter::SetEnergyRange";
00216     G4cerr << " :  illegal energy range" << "(" << lowedge/GeV;
00217     G4cerr << "," << highedge/GeV << ") [GeV]" << G4endl;
00218 #endif
00219     G4Exception( "G4VRangeToEnergyConverter::SetEnergyRange()",
00220                  "ProcCuts101",
00221                  JustWarning, "Illegal energy range ");
00222   } else {
00223     LowestEnergy = lowedge;
00224     HighestEnergy = highedge;
00225   }
00226 }
00227 
00228 
00229 G4double G4VRangeToEnergyConverter::GetLowEdgeEnergy()
00230 {
00231   return LowestEnergy;
00232 }
00233     
00234 
00235 G4double G4VRangeToEnergyConverter::GetHighEdgeEnergy()
00236 {
00237   return HighestEnergy;
00238 }
00239 
00240 // **********************************************************************
00241 // ******************* Get/SetMaxEnergyCut  *****************************
00242 // **********************************************************************
00243 G4double G4VRangeToEnergyConverter::GetMaxEnergyCut()
00244 {
00245   return MaxEnergyCut;
00246 }
00247 
00248 void G4VRangeToEnergyConverter::SetMaxEnergyCut(G4double value)
00249 {
00250   MaxEnergyCut = value;
00251 }
00252 
00253 // **********************************************************************
00254 // ************************ Reset  **************************************
00255 // **********************************************************************
00256 void G4VRangeToEnergyConverter::Reset()
00257 {
00258   // delete loss table
00259   if (theLossTable) {  
00260     theLossTable->clearAndDestroy();
00261     delete theLossTable;
00262   }
00263   theLossTable=0;
00264   NumberOfElements=0;
00265   
00266   //clear RangeVectorStore
00267   for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
00268     delete fRangeVectorStore.at(idx);
00269   }
00270   fRangeVectorStore.clear();
00271 } 
00272 
00273 
00274 // **********************************************************************
00275 // ************************ BuildLossTable ******************************
00276 // **********************************************************************
00277 //   create Energy Loss Table for charged particles 
00278 //   (cross section tabel for neutral )
00279 void G4VRangeToEnergyConverter::BuildLossTable()
00280 {
00281   if (size_t(NumberOfElements) == G4Element::GetNumberOfElements()) return;
00282   
00283   // clear Loss table and Range vectors
00284   Reset();
00285 
00286   //  Build dE/dx tables for elements
00287   NumberOfElements = G4Element::GetNumberOfElements();
00288   theLossTable = new G4LossTable();
00289   theLossTable->reserve(G4Element::GetNumberOfElements());
00290 #ifdef G4VERBOSE
00291   if (GetVerboseLevel()>3) {
00292     G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
00293     G4cout << "Create theLossTable[" << theLossTable << "]";
00294     G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
00295   }
00296 #endif
00297   
00298   
00299   // fill the loss table
00300   for (size_t j=0; j<size_t(NumberOfElements); j++){
00301     G4double Value;
00302     G4LossVector* aVector= 0;
00303     aVector= new G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
00304     for (size_t i=0; i<size_t(TotBin); i++) {
00305       Value = ComputeLoss(  (*G4Element::GetElementTable())[j]->GetZ(),
00306                             aVector->GetLowEdgeEnergy(i)
00307                             );
00308       aVector->PutValue(i,Value);
00309     }
00310     theLossTable->insert(aVector);
00311   }
00312 }
00313 
00314 // **********************************************************************
00315 // ************************ BuildRangeVector ****************************
00316 // **********************************************************************
00317 void G4VRangeToEnergyConverter::BuildRangeVector(const G4Material* aMaterial,
00318                                              G4PhysicsLogVector* rangeVector)
00319 {
00320   //  create range vector for a material
00321   const G4ElementVector* elementVector = aMaterial->GetElementVector();
00322   const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
00323   G4int NumEl = aMaterial->GetNumberOfElements();
00324 
00325   // calculate parameters of the low energy part first
00326   size_t i;
00327   std::vector<G4double> lossV;
00328   for ( size_t ib=0; ib<size_t(TotBin); ib++) {
00329     G4double loss=0.;
00330     for (i=0; i<size_t(NumEl); i++) {
00331       G4int IndEl = (*elementVector)[i]->GetIndex();
00332       loss += atomicNumDensityVector[i]*
00333                 (*((*theLossTable)[IndEl]))[ib];
00334     }
00335     lossV.push_back(loss);
00336   }
00337    
00338   // Integrate with Simpson formula with logarithmic binning
00339   G4double ltt = std::log(MaxEnergyCut/LowestEnergy);
00340   G4double dltau = ltt/TotBin;
00341 
00342   G4double s0 = 0.;
00343   G4double Value;
00344   for ( i=0; i<size_t(TotBin); i++) {
00345     G4double t = rangeVector->GetLowEdgeEnergy(i);
00346     G4double q = t/lossV[i];
00347     if (i==0) s0 += 0.5*q;
00348     else s0 += q;
00349     
00350     if (i==0) {
00351        Value = (s0 + 0.5*q)*dltau ;
00352     } else {
00353       Value = (s0 - 0.5*q)*dltau ;
00354     }
00355     rangeVector->PutValue(i,Value);
00356   }
00357 } 
00358 
00359 // **********************************************************************
00360 // ****************** ConvertCutToKineticEnergy *************************
00361 // **********************************************************************
00362 G4double G4VRangeToEnergyConverter::ConvertCutToKineticEnergy(
00363                                     G4RangeVector* rangeVector,
00364                                     G4double       theCutInLength, 
00365 #ifdef G4VERBOSE
00366                                     size_t         materialIndex
00367 #else
00368                                     size_t
00369 #endif
00370                                                               ) const
00371 {
00372   const G4double epsilon=0.01;
00373 
00374   //  find max. range and the corresponding energy (rmax,Tmax)
00375   G4double rmax= -1.e10*mm;
00376 
00377   G4double T1 = LowestEnergy;
00378   G4double r1 =(*rangeVector)[0] ;
00379 
00380   G4double T2 = MaxEnergyCut;
00381 
00382   // check theCutInLength < r1 
00383   if ( theCutInLength <= r1 ) {  return T1; }
00384 
00385   // scan range vector to find nearest bin 
00386   // ( suppose that r(Ti) > r(Tj) if Ti >Tj )
00387   for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
00388     G4double T=rangeVector->GetLowEdgeEnergy(ibin);
00389     G4double r=(*rangeVector)[ibin];
00390     if ( r>rmax )   rmax=r;
00391     if (r <theCutInLength ) {
00392       T1 = T;
00393       r1 = r;
00394     } else if (r >theCutInLength ) {
00395       T2 = T;
00396       break;
00397     }
00398   }
00399 
00400   // check cut in length is smaller than range max
00401   if ( theCutInLength >= rmax )  {
00402 #ifdef G4VERBOSE
00403     if (GetVerboseLevel()>2) {
00404       G4cout << "G4VRangeToEnergyConverter::ConvertCutToKineticEnergy ";
00405       G4cout << "  for " << theParticle->GetParticleName() << G4endl;
00406       G4cout << "The cut in range [" << theCutInLength/mm << " (mm)]  ";
00407       G4cout << " is too big  " ;
00408       G4cout << " for material  idx=" << materialIndex <<G4endl; 
00409     }
00410 #endif
00411     return  MaxEnergyCut;
00412   }
00413   
00414   // convert range to energy
00415   G4double T3 = std::sqrt(T1*T2);
00416   G4double r3 = rangeVector->Value(T3);
00417   while ( std::fabs(1.-r3/theCutInLength)>epsilon ) {
00418     if ( theCutInLength <= r3 ) {
00419       T2 = T3;
00420     } else {
00421       T1 = T3;
00422     }
00423     T3 = std::sqrt(T1*T2);
00424     r3 = rangeVector->Value(T3);
00425   }
00426 
00427   return T3;
00428 }
00429 

Generated on Mon May 27 17:50:20 2013 for Geant4 by  doxygen 1.4.7