G4EnergySplitter.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 #include "G4EnergySplitter.hh"
00027 #include "G4VSolid.hh"
00028 #include "G4UnitsTable.hh"
00029 #include "G4RegularNavigationHelper.hh"
00030 #include "G4EnergyLossForExtrapolator.hh"
00031 #include "G4EmCalculator.hh"
00032 #include "G4PhysicalVolumeStore.hh"
00033 #include "G4Step.hh"
00034 #include "G4PVParameterised.hh"
00035 
00037 // (Description)
00038 //
00039 // Created: 
00040 // 
00042 
00043 G4EnergySplitter::G4EnergySplitter()
00044 {
00045    theElossExt = new G4EnergyLossForExtrapolator(0);
00046    thePhantomParam = 0;
00047    theNIterations = 2;
00048 }
00049 
00050 G4EnergySplitter::~G4EnergySplitter()
00051 {;}
00052 
00053 G4int G4EnergySplitter::SplitEnergyInVolumes(const G4Step* aStep )
00054 {
00055   theEnergies.clear();
00056 
00057   G4double edep = aStep->GetTotalEnergyDeposit();
00058  
00059 #ifdef VERBOSE_ENERSPLIT
00060   G4bool verbose = 1;
00061   if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit() 
00062                        << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size() << G4endl;
00063 #endif    
00064   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 0 ||
00065       aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0)  { // we are only counting dose deposit
00066     return theEnergies.size();
00067   }
00068   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
00069     theEnergies.push_back(edep);
00070     return theEnergies.size();
00071   }
00072 
00073   if( !thePhantomParam ) GetPhantomParam(TRUE);
00074   
00075   if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
00076   
00077   //----- Distribute energy deposited in voxels 
00078   std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::Instance()->GetStepLengths(); 
00079 
00080   const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
00081   G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
00082   G4double kinEnergyPre = kinEnergyPreOrig;
00083   
00084   G4double stepLength = aStep->GetStepLength();
00085   G4double slSum = 0.;
00086   unsigned int ii;
00087   for( ii = 0; ii < rnsl.size(); ii++ ){
00088     G4double sl = rnsl[ii].second;
00089     slSum += sl;
00090 #ifdef VERBOSE_ENERSPLIT
00091     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
00092 #endif
00093   }
00094   
00095 #ifdef VERBOSE_ENERSPLIT
00096   if( verbose )
00097     G4cout << "G4EnergySplitter RN:  step length geom TOTAL " << slSum 
00098            << " true TOTAL " << stepLength 
00099            << " ratio " << stepLength/slSum 
00100            << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy() 
00101            << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName() 
00102            << " Number of geom steps " << rnsl.size() << G4endl;
00103 #endif
00104   //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
00105   if( theNIterations == 0 ) { 
00106     for( ii = 0; ii < rnsl.size(); ii++ ){
00107       G4double sl = rnsl[ii].second;
00108       G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
00109 #ifdef VERBOSE_ENERSPLIT
00110       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii 
00111                           << " edep " << edepStep << G4endl;
00112 #endif
00113         
00114       theEnergies.push_back(edepStep);
00115         
00116     }
00117   } else { //  1 or more iterations demanded
00118       
00119 #ifdef VERBOSE_ENERSPLIT
00120     // print corrected energy at iteration 0 
00121     if(verbose)  {
00122       G4double slSum = 0.;
00123       for( ii = 0; ii < rnsl.size(); ii++ ){
00124         G4double sl = rnsl[ii].second;
00125         slSum += sl;
00126       }
00127       for( ii = 0; ii < rnsl.size(); ii++ ){
00128         G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
00129                 << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum  
00130                 << G4endl;
00131       }
00132     }
00133 #endif
00134 
00135     G4double slRatio = stepLength/slSum;
00136 #ifdef VERBOSE_ENERSPLIT
00137     if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes  RN: iter 0, step ratio " << slRatio << G4endl;
00138 #endif
00139       
00140     //--- energy at each interaction
00141     G4EmCalculator emcalc;
00142     G4double totalELost = 0.;
00143     std::vector<G4double> stepLengths;
00144     for( int iiter = 1; iiter <= theNIterations; iiter++ ) {
00145       //--- iter1: distribute true step length in each voxel: geom SL in each voxel is multiplied by a constant so that the sum gives the total true step length
00146       if( iiter == 1 ) {
00147         for( ii = 0; ii < rnsl.size(); ii++ ){
00148           G4double sl = rnsl[ii].second;
00149           stepLengths.push_back( sl * slRatio );
00150 #ifdef VERBOSE_ENERSPLIT
00151           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
00152 #endif
00153         }
00154         
00155         for( ii = 0; ii < rnsl.size(); ii++ ){
00156           const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
00157           G4double dEdx = 0.;
00158           if( kinEnergyPre > 0. ) {  //t check this 
00159             dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
00160           }
00161           G4double elost = stepLengths[ii] * dEdx;
00162           
00163 #ifdef VERBOSE_ENERSPLIT
00164           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost "  << elost 
00165                               << " energy at interaction " << kinEnergyPre 
00166                               << " = stepLength " << stepLengths[ii] 
00167                               << " * dEdx " << dEdx << G4endl;
00168 #endif
00169           kinEnergyPre -= elost;
00170           theEnergies.push_back( elost );
00171           totalELost += elost;
00172         }
00173         
00174       } else{
00175         //------ 2nd and other iterations
00176         //----- Get step lengths corrected by changing geom2true correction
00177         //-- Get ratios for each energy 
00178         slSum = 0.;
00179         kinEnergyPre = kinEnergyPreOrig;
00180         for( ii = 0; ii < rnsl.size(); ii++ ){
00181           const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
00182           stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
00183           kinEnergyPre -= theEnergies[ii];
00184           
00185 #ifdef VERBOSE_ENERSPLIT
00186           if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii 
00187                              << " RN: iter" << iiter << " step length geom " << stepLengths[ii] 
00188                              << " geom2true " << rnsl[ii].second / stepLengths[ii]  << G4endl;
00189 #endif
00190           
00191           slSum += stepLengths[ii];
00192         }
00193         
00194         //Correct step lengths so that they sum the total step length
00195         G4double slratio = aStep->GetStepLength()/slSum;
00196 #ifdef VERBOSE_ENERSPLIT
00197         if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
00198 #endif
00199         for( ii = 0; ii < rnsl.size(); ii++ ){
00200           stepLengths[ii] *= slratio;
00201 #ifdef VERBOSE_ENERSPLIT
00202           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
00203 #endif
00204         }
00205         
00206         //---- Recalculate energy lost with this new step lengths
00207         kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
00208         totalELost = 0.;
00209         for( ii = 0; ii < rnsl.size(); ii++ ){
00210           const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
00211           G4double dEdx = 0.;
00212           if( kinEnergyPre > 0. ) {
00213             dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
00214           }
00215           G4double elost = stepLengths[ii] * dEdx;
00216 #ifdef VERBOSE_ENERSPLIT
00217           if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost 
00218                               << " energy at interaction " << kinEnergyPre 
00219                               << " = stepLength " << stepLengths[ii] 
00220                               << " * dEdx " << dEdx << G4endl;
00221 #endif
00222           kinEnergyPre -= elost;
00223           theEnergies[ii] = elost;
00224           totalELost += elost;
00225         }
00226         
00227       }
00228       
00229       //correct energies so that they reproduce the real step energy lost
00230       G4double enerRatio = (edep/totalELost);
00231       
00232 #ifdef VERBOSE_ENERSPLIT
00233       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
00234 #endif
00235         
00236 #ifdef VERBOSE_ENERSPLIT
00237       G4double elostTot = 0.; 
00238 #endif
00239       for( ii = 0; ii < theEnergies.size(); ii++ ){
00240         theEnergies[ii] *= enerRatio;
00241 #ifdef VERBOSE_ENERSPLIT
00242         elostTot += theEnergies[ii];
00243         if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii] 
00244                             << " orig elost " << theEnergies[ii]/enerRatio 
00245                             << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
00246                             << " energy after interaction " << kinEnergyPreOrig-elostTot
00247                             << G4endl;
00248 #endif
00249       }
00250     }
00251     
00252   }
00253   
00254   return theEnergies.size();
00255 }
00256 
00257 
00258 //-----------------------------------------------------------------------
00259 void G4EnergySplitter::GetPhantomParam(G4bool mustExist)
00260 {
00261   G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
00262   std::vector<G4VPhysicalVolume*>::iterator cite;
00263   for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
00264     //    G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
00265     if( IsPhantomVolume( *cite ) ) {
00266       const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
00267       G4VPVParameterisation* param = pvparam->GetParameterisation();
00268       //    if( static_cast<const G4PhantomParameterisation*>(param) ){
00269       //    if( static_cast<const G4PhantomParameterisation*>(param) ){
00270       //      G4cout << "G4PhantomParameterisation volume found  " << (*cite)->GetName() << G4endl;
00271       thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
00272     }
00273   }
00274   
00275   if( !thePhantomParam && mustExist )
00276     G4Exception("G4EnergySplitter::GetPhantomParam",
00277                 "PhantomParamError", FatalException,
00278                 "No G4PhantomParameterisation found !");
00279 }
00280 
00281 
00282 //-----------------------------------------------------------------------
00283 G4bool G4EnergySplitter::IsPhantomVolume( G4VPhysicalVolume* pv )
00284 {
00285   EAxis axis;
00286   G4int nReplicas;
00287   G4double width,offset;
00288   G4bool consuming;
00289   pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
00290   EVolume type = (consuming) ? kReplica : kParameterised;
00291   if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {  
00292     return TRUE;
00293   } else {
00294     return FALSE;
00295   }
00296 
00297 } 
00298 
00299 //-----------------------------------------------------------------------
00300 void G4EnergySplitter::GetLastVoxelID( G4int& voxelID)
00301 {       
00302   voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().begin())).first;
00303 }
00304 
00305 //-----------------------------------------------------------------------
00306 void G4EnergySplitter::GetFirstVoxelID( G4int& voxelID)
00307 {
00308   voxelID =  (*(G4RegularNavigationHelper::Instance()->GetStepLengths().rbegin())).first;
00309 }
00310 
00311 //-----------------------------------------------------------------------
00312 void G4EnergySplitter::GetVoxelID( G4int stepNo, G4int& voxelID )
00313 {
00314   if( stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()) ) {
00315   G4Exception("G4EnergySplitter::GetVoxelID",
00316               "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
00317               FatalErrorInArgument,
00318               G4String("stepNo = " + G4UIcommand::ConvertToString(stepNo) + ", number of voxels = " + G4UIcommand::ConvertToString(G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size())) ).c_str());
00319   }
00320   std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
00321   advance( ite, stepNo );
00322   voxelID = (*ite).first;
00323 
00324 }
00325 
00326 
00327 //-----------------------------------------------------------------------
00328 void G4EnergySplitter::GetStepLength( G4int stepNo, G4double& stepLength )
00329 {
00330   std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
00331   advance( ite, stepNo );
00332   stepLength = (*ite).second;
00333 }

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