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

Last change on this file since 1353 was 1350, checked in by garnier, 15 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.