source: trunk/source/processes/scoring/src/G4EnergySplitter.cc

Last change on this file was 1350, checked in by garnier, 14 years ago

update to last version 4.9.4
File size: 11.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26#include "G4EnergySplitter.hh"
27#include "G4VSolid.hh"
28#include "G4UnitsTable.hh"
29#include "G4RegularNavigationHelper.hh"
30#include "G4EnergyLossForExtrapolator.hh"
31#include "G4EmCalculator.hh"
32#include "G4PhysicalVolumeStore.hh"
33#include "G4Step.hh"
34#include "G4PVParameterised.hh"
35
36////////////////////////////////////////////////////////////////////////////////
37// (Description)
38//
39// Created:
40//
41///////////////////////////////////////////////////////////////////////////////
42
43G4EnergySplitter::G4EnergySplitter()
44{
45   theElossExt = new G4EnergyLossForExtrapolator(0);
46   thePhantomParam = 0;
47   theNIterations = 2;
48}
49
50G4EnergySplitter::~G4EnergySplitter()
51{;}
52
53G4int G4EnergySplitter::SplitEnergyInVolumes(const G4Step* aStep )
54{
55  theEnergies.clear();
56
57  G4double edep = aStep->GetTotalEnergyDeposit();
58 
59#ifdef VERBOSE_ENERSPLIT
60  G4bool verbose = 1;
61  if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit() 
62                       << " Nsteps " << G4RegularNavigationHelper::theStepLengths.size() << G4endl;
63#endif   
64  if( G4RegularNavigationHelper::theStepLengths.size() == 0 ||
65      aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0)  { // we are only counting dose deposit
66    return theEnergies.size();
67  }
68  if( G4RegularNavigationHelper::theStepLengths.size() == 1 ) {
69    theEnergies.push_back(edep);
70    return theEnergies.size();
71  }
72
73  if( !thePhantomParam ) GetPhantomParam(TRUE);
74 
75  if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
76 
77  //----- Distribute energy deposited in voxels
78  std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::theStepLengths; 
79
80  const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
81  G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
82  G4double kinEnergyPre = kinEnergyPreOrig;
83 
84  G4double stepLength = aStep->GetStepLength();
85  G4double slSum = 0.;
86  unsigned int ii;
87  for( ii = 0; ii < rnsl.size(); ii++ ){
88    G4double sl = rnsl[ii].second;
89    slSum += sl;
90#ifdef VERBOSE_ENERSPLIT
91    if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
92#endif
93  }
94 
95#ifdef VERBOSE_ENERSPLIT
96  if( verbose )
97    G4cout << "G4EnergySplitter RN:  step length geom TOTAL " << slSum
98           << " true TOTAL " << stepLength
99           << " ratio " << stepLength/slSum
100           << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy() 
101           << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName() 
102           << " Number of geom steps " << rnsl.size() << G4endl;
103#endif
104  //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
105  if( theNIterations == 0 ) { 
106    for( unsigned int ii = 0; ii < rnsl.size(); ii++ ){
107      G4double sl = G4RegularNavigationHelper::theStepLengths[ii].second;
108      G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
109#ifdef VERBOSE_ENERSPLIT
110      if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii
111                          << " edep " << edepStep << G4endl;
112#endif
113       
114      theEnergies.push_back(edepStep);
115       
116    }
117  } else { //  1 or more iterations demanded
118     
119#ifdef VERBOSE_ENERSPLIT
120    // print corrected energy at iteration 0
121    if(verbose)  {
122      G4double slSum = 0.;
123      for( ii = 0; ii < rnsl.size(); ii++ ){
124        G4double sl = rnsl[ii].second;
125        slSum += sl;
126      }
127      for( ii = 0; ii < rnsl.size(); ii++ ){
128        G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
129                << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum 
130                << G4endl;
131      }
132    }
133#endif
134
135    G4double slRatio = stepLength/slSum;
136#ifdef VERBOSE_ENERSPLIT
137    if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes  RN: iter 0, step ratio " << slRatio << G4endl;
138#endif
139     
140    //--- energy at each interaction
141    G4EmCalculator emcalc;
142    G4double totalELost = 0.;
143    std::vector<G4double> stepLengths;
144    for( int iiter = 1; iiter <= theNIterations; iiter++ ) {
145      //--- 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
146      if( iiter == 1 ) {
147        for( ii = 0; ii < rnsl.size(); ii++ ){
148          G4double sl = rnsl[ii].second;
149          stepLengths.push_back( sl * slRatio );
150#ifdef VERBOSE_ENERSPLIT
151          if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
152#endif
153        }
154       
155        for( ii = 0; ii < rnsl.size(); ii++ ){
156          const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
157          G4double dEdx = 0.;
158          if( kinEnergyPre > 0. ) {  //t check this
159            dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
160          }
161          G4double elost = stepLengths[ii] * dEdx;
162         
163#ifdef VERBOSE_ENERSPLIT
164          if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost "  << elost
165                              << " energy at interaction " << kinEnergyPre
166                              << " = stepLength " << stepLengths[ii] 
167                              << " * dEdx " << dEdx << G4endl;
168#endif
169          kinEnergyPre -= elost;
170          theEnergies.push_back( elost );
171          totalELost += elost;
172        }
173       
174      } else{
175        //------ 2nd and other iterations
176        //----- Get step lengths corrected by changing geom2true correction
177        //-- Get ratios for each energy
178        slSum = 0.;
179        kinEnergyPre = kinEnergyPreOrig;
180        for( ii = 0; ii < rnsl.size(); ii++ ){
181          const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
182          stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
183          kinEnergyPre -= theEnergies[ii];
184         
185#ifdef VERBOSE_ENERSPLIT
186          if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii
187                             << " RN: iter" << iiter << " step length geom " << stepLengths[ii] 
188                             << " geom2true " << rnsl[ii].second / stepLengths[ii]  << G4endl;
189#endif
190         
191          slSum += stepLengths[ii];
192        }
193       
194        //Correct step lengths so that they sum the total step length
195        G4double slratio = aStep->GetStepLength()/slSum;
196#ifdef VERBOSE_ENERSPLIT
197        if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
198#endif
199        for( ii = 0; ii < rnsl.size(); ii++ ){
200          stepLengths[ii] *= slratio;
201#ifdef VERBOSE_ENERSPLIT
202          if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
203#endif
204        }
205       
206        //---- Recalculate energy lost with this new step lengths
207        G4double kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
208        totalELost = 0.;
209        for( ii = 0; ii < rnsl.size(); ii++ ){
210          const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
211          G4double dEdx = 0.;
212          if( kinEnergyPre > 0. ) {
213            dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
214          }
215          G4double elost = stepLengths[ii] * dEdx;
216#ifdef VERBOSE_ENERSPLIT
217          if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost
218                              << " energy at interaction " << kinEnergyPre
219                              << " = stepLength " << stepLengths[ii] 
220                              << " * dEdx " << dEdx << G4endl;
221#endif
222          kinEnergyPre -= elost;
223          theEnergies[ii] = elost;
224          totalELost += elost;
225        }
226       
227      }
228     
229      //correct energies so that they reproduce the real step energy lost
230      G4double enerRatio = (edep/totalELost);
231     
232#ifdef VERBOSE_ENERSPLIT
233      if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
234#endif
235       
236#ifdef VERBOSE_ENERSPLIT
237      G4double elostTot = 0.; 
238#endif
239      for( ii = 0; ii < theEnergies.size(); ii++ ){
240        theEnergies[ii] *= enerRatio;
241#ifdef VERBOSE_ENERSPLIT
242        elostTot += theEnergies[ii];
243        if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii] 
244                            << " orig elost " << theEnergies[ii]/enerRatio
245                            << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
246                            << " energy after interaction " << kinEnergyPreOrig-elostTot
247                            << G4endl;
248#endif
249      }
250    }
251   
252  }
253 
254  return theEnergies.size();
255}
256
257
258//-----------------------------------------------------------------------
259void G4EnergySplitter::GetPhantomParam(G4bool mustExist)
260{
261  G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
262  std::vector<G4VPhysicalVolume*>::iterator cite;
263  for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
264    //    G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
265    if( IsPhantomVolume( *cite ) ) {
266      const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
267      G4VPVParameterisation* param = pvparam->GetParameterisation();
268      //    if( static_cast<const G4PhantomParameterisation*>(param) ){
269      //    if( static_cast<const G4PhantomParameterisation*>(param) ){
270      //      G4cout << "G4PhantomParameterisation volume found  " << (*cite)->GetName() << G4endl;
271      thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
272    }
273  }
274 
275  if( !thePhantomParam && mustExist ) G4Exception("GmRegularParamUtils::GetPhantomParam:  No G4PhantomParameterisation found ");
276 
277 
278}
279
280
281//-----------------------------------------------------------------------
282G4bool G4EnergySplitter::IsPhantomVolume( G4VPhysicalVolume* pv )
283{
284  EAxis axis;
285  G4int nReplicas;
286  G4double width,offset;
287  G4bool consuming;
288  pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
289  EVolume type = (consuming) ? kReplica : kParameterised;
290  if( type == kParameterised && pv->GetRegularStructureId() == 1 ) { 
291    return TRUE;
292  } else {
293    return FALSE;
294  }
295
296} 
297
Note: See TracBrowser for help on using the repository browser.