ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EnergySplitter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EnergySplitter.cc
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"
31 #include "G4EmCalculator.hh"
32 #include "G4PhysicalVolumeStore.hh"
33 #include "G4Step.hh"
34 #include "G4PVParameterised.hh"
35 
37 // (Description)
38 //
39 // Created:
40 //
42 
44 {
46  thePhantomParam = 0;
47  theNIterations = 2;
48 }
49 
51 {
52  delete theElossExt;
53 }
54 
56 {
57  theEnergies.clear();
58 
60 
61 #ifdef VERBOSE_ENERSPLIT
62  G4bool verbose = 1;
63  if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit()
64  << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size() << G4endl;
65 #endif
66  if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 0 ||
67  aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0) { // we are only counting dose deposit
68  return theEnergies.size();
69  }
70  if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
71  theEnergies.push_back(edep);
72  return theEnergies.size();
73  }
74 
76 
77  if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
78 
79  //----- Distribute energy deposited in voxels
80  std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::Instance()->GetStepLengths();
81 
82  const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
83  G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
84  G4double kinEnergyPre = kinEnergyPreOrig;
85 
86  G4double stepLength = aStep->GetStepLength();
87  G4double slSum = 0.;
88  unsigned int ii;
89  for( ii = 0; ii < rnsl.size(); ii++ ){
90  G4double sl = rnsl[ii].second;
91  slSum += sl;
92 #ifdef VERBOSE_ENERSPLIT
93  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
94 #endif
95  }
96 
97 #ifdef VERBOSE_ENERSPLIT
98  if( verbose )
99  G4cout << "G4EnergySplitter RN: step length geom TOTAL " << slSum
100  << " true TOTAL " << stepLength
101  << " ratio " << stepLength/slSum
102  << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy()
103  << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName()
104  << " Number of geom steps " << rnsl.size() << G4endl;
105 #endif
106  //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
107  if( theNIterations == 0 ) {
108  for( ii = 0; ii < rnsl.size(); ii++ ){
109  G4double sl = rnsl[ii].second;
110  G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
111 #ifdef VERBOSE_ENERSPLIT
112  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii
113  << " edep " << edepStep << G4endl;
114 #endif
115 
116  theEnergies.push_back(edepStep);
117 
118  }
119  } else { // 1 or more iterations demanded
120 
121 #ifdef VERBOSE_ENERSPLIT
122  // print corrected energy at iteration 0
123  if(verbose) {
124  G4double slSum = 0.;
125  for( ii = 0; ii < rnsl.size(); ii++ ){
126  G4double sl = rnsl[ii].second;
127  slSum += sl;
128  }
129  for( ii = 0; ii < rnsl.size(); ii++ ){
130  G4cout << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
131  << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum
132  << G4endl;
133  }
134  }
135 #endif
136 
137  G4double slRatio = stepLength/slSum;
138 #ifdef VERBOSE_ENERSPLIT
139  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes RN: iter 0, step ratio " << slRatio << G4endl;
140 #endif
141 
142  //--- energy at each interaction
143  G4EmCalculator emcalc;
144  G4double totalELost = 0.;
145  std::vector<G4double> stepLengths;
146  for( int iiter = 1; iiter <= theNIterations; iiter++ ) {
147  //--- 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
148  if( iiter == 1 ) {
149  for( ii = 0; ii < rnsl.size(); ii++ ){
150  G4double sl = rnsl[ii].second;
151  stepLengths.push_back( sl * slRatio );
152 #ifdef VERBOSE_ENERSPLIT
153  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
154 #endif
155  }
156 
157  for( ii = 0; ii < rnsl.size(); ii++ ){
158  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
159  G4double dEdx = 0.;
160  if( kinEnergyPre > 0. ) { //t check this
161  dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
162  }
163  G4double elost = stepLengths[ii] * dEdx;
164 
165 #ifdef VERBOSE_ENERSPLIT
166  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost " << elost
167  << " energy at interaction " << kinEnergyPre
168  << " = stepLength " << stepLengths[ii]
169  << " * dEdx " << dEdx << G4endl;
170 #endif
171  kinEnergyPre -= elost;
172  theEnergies.push_back( elost );
173  totalELost += elost;
174  }
175 
176  } else{
177  //------ 2nd and other iterations
178  //----- Get step lengths corrected by changing geom2true correction
179  //-- Get ratios for each energy
180  slSum = 0.;
181  kinEnergyPre = kinEnergyPreOrig;
182  for( ii = 0; ii < rnsl.size(); ii++ ){
183  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
184  stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
185  kinEnergyPre -= theEnergies[ii];
186 
187 #ifdef VERBOSE_ENERSPLIT
188  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii
189  << " RN: iter" << iiter << " step length geom " << stepLengths[ii]
190  << " geom2true " << rnsl[ii].second / stepLengths[ii] << G4endl;
191 #endif
192 
193  slSum += stepLengths[ii];
194  }
195 
196  //Correct step lengths so that they sum the total step length
197  G4double slratio = aStep->GetStepLength()/slSum;
198 #ifdef VERBOSE_ENERSPLIT
199  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
200 #endif
201  for( ii = 0; ii < rnsl.size(); ii++ ){
202  stepLengths[ii] *= slratio;
203 #ifdef VERBOSE_ENERSPLIT
204  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
205 #endif
206  }
207 
208  //---- Recalculate energy lost with this new step lengths
209  kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
210  totalELost = 0.;
211  for( ii = 0; ii < rnsl.size(); ii++ ){
212  const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
213  G4double dEdx = 0.;
214  if( kinEnergyPre > 0. ) {
215  dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
216  }
217  G4double elost = stepLengths[ii] * dEdx;
218 #ifdef VERBOSE_ENERSPLIT
219  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost
220  << " energy at interaction " << kinEnergyPre
221  << " = stepLength " << stepLengths[ii]
222  << " * dEdx " << dEdx << G4endl;
223 #endif
224  kinEnergyPre -= elost;
225  theEnergies[ii] = elost;
226  totalELost += elost;
227  }
228 
229  }
230 
231  //correct energies so that they reproduce the real step energy lost
232  G4double enerRatio = (edep/totalELost);
233 
234 #ifdef VERBOSE_ENERSPLIT
235  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
236 #endif
237 
238 #ifdef VERBOSE_ENERSPLIT
239  G4double elostTot = 0.;
240 #endif
241  for( ii = 0; ii < theEnergies.size(); ii++ ){
242  theEnergies[ii] *= enerRatio;
243 #ifdef VERBOSE_ENERSPLIT
244  elostTot += theEnergies[ii];
245  if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii]
246  << " orig elost " << theEnergies[ii]/enerRatio
247  << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
248  << " energy after interaction " << kinEnergyPreOrig-elostTot
249  << G4endl;
250 #endif
251  }
252  }
253 
254  }
255 
256  return theEnergies.size();
257 }
258 
259 
260 //-----------------------------------------------------------------------
262 {
264  std::vector<G4VPhysicalVolume*>::iterator cite;
265  for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
266  // G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
267  if( IsPhantomVolume( *cite ) ) {
268  const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
269  G4VPVParameterisation* param = pvparam->GetParameterisation();
270  // if( static_cast<const G4PhantomParameterisation*>(param) ){
271  // if( static_cast<const G4PhantomParameterisation*>(param) ){
272  // G4cout << "G4PhantomParameterisation volume found " << (*cite)->GetName() << G4endl;
273  thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
274  }
275  }
276 
277  if( !thePhantomParam && mustExist )
278  G4Exception("G4EnergySplitter::GetPhantomParam",
279  "PhantomParamError", FatalException,
280  "No G4PhantomParameterisation found !");
281 }
282 
283 
284 //-----------------------------------------------------------------------
286 {
287  EAxis axis;
288  G4int nReplicas;
290  G4bool consuming;
291  pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
292  EVolume type = (consuming) ? kReplica : kParameterised;
293  if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {
294  return TRUE;
295  } else {
296  return FALSE;
297  }
298 
299 }
300 
301 //-----------------------------------------------------------------------
303 {
304  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().begin())).first;
305 }
306 
307 //-----------------------------------------------------------------------
309 {
310  voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().rbegin())).first;
311 }
312 
313 //-----------------------------------------------------------------------
314 void G4EnergySplitter::GetVoxelID( G4int stepNo, G4int& voxelID )
315 {
316  if( stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()) ) {
317  G4Exception("G4EnergySplitter::GetVoxelID",
318  "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
320  G4String("stepNo = " + G4UIcommand::ConvertToString(stepNo) + ", number of voxels = " + G4UIcommand::ConvertToString(G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size())) ).c_str());
321  }
322  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
323  advance( ite, stepNo );
324  voxelID = (*ite).first;
325 
326 }
327 
328 
329 //-----------------------------------------------------------------------
330 void G4EnergySplitter::GetStepLength( G4int stepNo, G4double& stepLength )
331 {
332  std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
333  advance( ite, stepNo );
334  stepLength = (*ite).second;
335 }