Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/scoring/src/G4EnergySplitter.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/scoring/src/G4EnergySplitter.cc (Version 11.3.0) and /processes/scoring/src/G4EnergySplitter.cc (Version 11.1.2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 #include "G4EnergySplitter.hh"                     26 #include "G4EnergySplitter.hh"
 27                                                <<  27 #include "G4VSolid.hh"
 28 #include "G4EmCalculator.hh"                   <<  28 #include "G4UnitsTable.hh"
                                                   >>  29 #include "G4RegularNavigationHelper.hh"
 29 #include "G4EnergyLossForExtrapolator.hh"          30 #include "G4EnergyLossForExtrapolator.hh"
 30 #include "G4PVParameterised.hh"                <<  31 #include "G4EmCalculator.hh"
 31 #include "G4PhysicalVolumeStore.hh"                32 #include "G4PhysicalVolumeStore.hh"
 32 #include "G4RegularNavigationHelper.hh"        << 
 33 #include "G4Step.hh"                               33 #include "G4Step.hh"
 34 #include "G4UnitsTable.hh"                     <<  34 #include "G4PVParameterised.hh"
 35 #include "G4VSolid.hh"                         << 
 36                                                    35 
 37 //////////////////////////////////////////////     36 ////////////////////////////////////////////////////////////////////////////////
 38 // (Description)                                   37 // (Description)
 39 //                                                 38 //
 40 // Created:                                    <<  39 // Created: 
 41 //                                             <<  40 // 
 42 //////////////////////////////////////////////     41 ///////////////////////////////////////////////////////////////////////////////
 43                                                    42 
 44 G4EnergySplitter::G4EnergySplitter()               43 G4EnergySplitter::G4EnergySplitter()
 45 {                                                  44 {
 46   theElossExt = new G4EnergyLossForExtrapolato <<  45    theElossExt = new G4EnergyLossForExtrapolator(0);
 47   thePhantomParam = nullptr;                   <<  46    thePhantomParam = 0;
 48   theNIterations = 2;                          <<  47    theNIterations = 2;
 49 }                                                  48 }
 50                                                    49 
 51 G4EnergySplitter::~G4EnergySplitter()              50 G4EnergySplitter::~G4EnergySplitter()
 52 {                                                  51 {
 53   delete theElossExt;                          <<  52   delete theElossExt;
 54 }                                                  53 }
 55                                                    54 
 56 G4int G4EnergySplitter::SplitEnergyInVolumes(c <<  55 G4int G4EnergySplitter::SplitEnergyInVolumes(const G4Step* aStep )
 57 {                                                  56 {
 58   theEnergies.clear();                             57   theEnergies.clear();
 59                                                    58 
 60   if (aStep == nullptr) return false;  // it i << 
 61                                                << 
 62   G4double edep = aStep->GetTotalEnergyDeposit     59   G4double edep = aStep->GetTotalEnergyDeposit();
 63                                                <<  60  
 64   if( edep == 0. ) {                           << 
 65     return 0;                                  << 
 66   }                                            << 
 67                                                << 
 68 #ifdef VERBOSE_ENERSPLIT                           61 #ifdef VERBOSE_ENERSPLIT
 69   G4bool verbose = 1;                              62   G4bool verbose = 1;
 70   if (verbose)                                 <<  63   if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit() 
 71     G4cout << "G4EnergySplitter::SplitEnergyIn <<  64            << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size() << G4endl;
 72            << " Nsteps " << G4RegularNavigatio <<  65 #endif    
 73            << G4endl;                          <<  66   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 0 ||
 74 #endif                                         <<  67       aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0)  { // we are only counting dose deposit
 75   if (G4RegularNavigationHelper::Instance()->G << 
 76       || aStep->GetTrack()->GetDefinition()->G << 
 77   {  // we are only counting dose deposit      << 
 78     return (G4int)theEnergies.size();              68     return (G4int)theEnergies.size();
 79   }                                                69   }
 80   if (G4RegularNavigationHelper::Instance()->G <<  70   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
 81     theEnergies.push_back(edep);                   71     theEnergies.push_back(edep);
 82     return (G4int)theEnergies.size();              72     return (G4int)theEnergies.size();
 83   }                                                73   }
 84                                                    74 
 85   if (thePhantomParam == nullptr) GetPhantomPa <<  75   if( !thePhantomParam ) GetPhantomParam(TRUE);
 86                                                <<  76   
 87   //----- Distribute energy deposited in voxel <<  77   if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
 88   std::vector<std::pair<G4int, G4double>> rnsl <<  78   
 89     G4RegularNavigationHelper::Instance()->Get <<  79   //----- Distribute energy deposited in voxels 
                                                   >>  80   std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::Instance()->GetStepLengths(); 
 90                                                    81 
 91   const G4ParticleDefinition* part = aStep->Ge     82   const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
 92   G4double kinEnergyPreOrig = aStep->GetPreSte     83   G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
 93   G4double kinEnergyPre = kinEnergyPreOrig;        84   G4double kinEnergyPre = kinEnergyPreOrig;
 94                                                <<  85   
 95   G4double stepLength = aStep->GetStepLength()     86   G4double stepLength = aStep->GetStepLength();
 96   G4double slSum = 0.;                             87   G4double slSum = 0.;
 97   unsigned int ii;                                 88   unsigned int ii;
 98   for (ii = 0; ii < rnsl.size(); ++ii) {       <<  89   for( ii = 0; ii < rnsl.size(); ++ii ){
 99     G4double sl = rnsl[ii].second;                 90     G4double sl = rnsl[ii].second;
100     slSum += sl;                                   91     slSum += sl;
101 #ifdef VERBOSE_ENERSPLIT                           92 #ifdef VERBOSE_ENERSPLIT
102     if (verbose)                               <<  93     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
103       G4cout << "G4EnergySplitter::SplitEnergy << 
104              << sl << G4endl;                  << 
105 #endif                                             94 #endif
106   }                                                95   }
107                                                <<  96   
108 #ifdef VERBOSE_ENERSPLIT                           97 #ifdef VERBOSE_ENERSPLIT
109   if (verbose)                                 <<  98   if( verbose )
110     G4cout << "G4EnergySplitter RN:  step leng <<  99     G4cout << "G4EnergySplitter RN:  step length geom TOTAL " << slSum 
111            << stepLength << " ratio " << stepL << 100      << " true TOTAL " << stepLength 
112            << aStep->GetPreStepPoint()->GetKin << 101      << " ratio " << stepLength/slSum 
113            << aStep->GetPreStepPoint()->GetMat << 102      << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy() 
114            << rnsl.size() << G4endl;           << 103      << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName() 
115 #endif                                         << 104      << " Number of geom steps " << rnsl.size() << G4endl;
116   //----- No iterations to correct elost and m << 105 #endif
117   // geometrical step length in each voxel     << 106   //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
118   if (theNIterations == 0) {                   << 107   if( theNIterations == 0 ) { 
119     for (ii = 0; ii < rnsl.size(); ++ii) {     << 108     for( ii = 0; ii < rnsl.size(); ++ii ){
120       G4double sl = rnsl[ii].second;              109       G4double sl = rnsl[ii].second;
121       G4double edepStep = edep * sl / slSum;   << 110       G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
122                                                << 
123 #ifdef VERBOSE_ENERSPLIT                          111 #ifdef VERBOSE_ENERSPLIT
124       if (verbose)                             << 112       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii 
125         G4cout << "G4EnergySplitter::SplitEner << 113         << " edep " << edepStep << G4endl;
126 #endif                                            114 #endif
127                                                << 115   
128       theEnergies.push_back(edepStep);            116       theEnergies.push_back(edepStep);
                                                   >> 117   
129     }                                             118     }
130   }                                            << 119   } else { //  1 or more iterations demanded
131   else {  //  1 or more iterations demanded    << 120       
132                                                << 
133 #ifdef VERBOSE_ENERSPLIT                          121 #ifdef VERBOSE_ENERSPLIT
134     // print corrected energy at iteration 0   << 122     // print corrected energy at iteration 0 
135     if (verbose) {                             << 123     if(verbose)  {
136       G4double slSum = 0.;                        124       G4double slSum = 0.;
137       for (ii = 0; ii < rnsl.size(); ++ii) {   << 125       for( ii = 0; ii < rnsl.size(); ++ii ){
138         G4double sl = rnsl[ii].second;         << 126   G4double sl = rnsl[ii].second;
139         slSum += sl;                           << 127   slSum += sl;
140       }                                           128       }
141       for (ii = 0; ii < rnsl.size(); ii++) {   << 129       for( ii = 0; ii < rnsl.size(); ii++ ){
142         G4cout << "G4EnergySplitter::SplitEner << 130   G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
143                << " RN: iter0 corrected energy << 131     << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum  
                                                   >> 132     << G4endl;
144       }                                           133       }
145     }                                             134     }
146 #endif                                            135 #endif
147                                                   136 
148     G4double slRatio = stepLength / slSum;     << 137     G4double slRatio = stepLength/slSum;
149 #ifdef VERBOSE_ENERSPLIT                          138 #ifdef VERBOSE_ENERSPLIT
150     if (verbose)                               << 139     if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes  RN: iter 0, step ratio " << slRatio << G4endl;
151       G4cout << "G4EnergySplitter::SplitEnergy << 
152              << G4endl;                        << 
153 #endif                                            140 #endif
154                                                << 141       
155     //--- energy at each interaction              142     //--- energy at each interaction
156     G4EmCalculator emcalc;                        143     G4EmCalculator emcalc;
157     G4double totalELost = 0.;                     144     G4double totalELost = 0.;
158     std::vector<G4double> stepLengths;            145     std::vector<G4double> stepLengths;
159     for (G4int iiter = 1; iiter <= theNIterati << 146     for( G4int iiter = 1; iiter <= theNIterations; ++iiter ) {
160       //--- iter1: distribute true step length << 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
161       // a constant so that the sum gives the  << 148       if( iiter == 1 ) {
162       if (iiter == 1) {                        << 149   for( ii = 0; ii < rnsl.size(); ++ii ){
163         for (ii = 0; ii < rnsl.size(); ++ii) { << 150     G4double sl = rnsl[ii].second;
164           G4double sl = rnsl[ii].second;       << 151     stepLengths.push_back( sl * slRatio );
165           stepLengths.push_back(sl * slRatio); << 152 #ifdef VERBOSE_ENERSPLIT
166 #ifdef VERBOSE_ENERSPLIT                       << 153     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
167           if (verbose)                         << 154 #endif
168             G4cout << "G4EnergySplitter::Split << 155   }
169                    << " corrected step length  << 156   
170 #endif                                         << 157   for( ii = 0; ii < rnsl.size(); ++ii ){
171         }                                      << 158     const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
172                                                << 159     G4double dEdx = 0.;
173         for (ii = 0; ii < rnsl.size(); ++ii) { << 160     if( kinEnergyPre > 0. ) {  //t check this 
174           const G4Material* mate = thePhantomP << 161       dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
175           G4double dEdx = 0.;                  << 162     }
176           if (kinEnergyPre > 0.) {  // t check << 163     G4double elost = stepLengths[ii] * dEdx;
177             dEdx = emcalc.GetDEDX(kinEnergyPre << 164     
178           }                                    << 165 #ifdef VERBOSE_ENERSPLIT
179           G4double elost = stepLengths[ii] * d << 166     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost "  << elost 
180                                                << 167             << " energy at interaction " << kinEnergyPre 
181 #ifdef VERBOSE_ENERSPLIT                       << 168             << " = stepLength " << stepLengths[ii] 
182           if (verbose)                         << 169             << " * dEdx " << dEdx << G4endl;
183             G4cout << "G4EnergySplitter::Split << 170 #endif
184                    << elost << " energy at int << 171     kinEnergyPre -= elost;
185                    << stepLengths[ii] << " * d << 172     theEnergies.push_back( elost );
186 #endif                                         << 173     totalELost += elost;
187           kinEnergyPre -= elost;               << 174   }
188           theEnergies.push_back(elost);        << 175   
189           totalELost += elost;                 << 176       } else{
190         }                                      << 177   //------ 2nd and other iterations
191       }                                        << 178   //----- Get step lengths corrected by changing geom2true correction
192       else {                                   << 179   //-- Get ratios for each energy 
193         //------ 2nd and other iterations      << 180   slSum = 0.;
194         //----- Get step lengths corrected by  << 181   kinEnergyPre = kinEnergyPreOrig;
195         //-- Get ratios for each energy        << 182   for( ii = 0; ii < rnsl.size(); ++ii ){
196         slSum = 0.;                            << 183     const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
197         kinEnergyPre = kinEnergyPreOrig;       << 184     stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
198         for (ii = 0; ii < rnsl.size(); ++ii) { << 185     kinEnergyPre -= theEnergies[ii];
199           const G4Material* mate = thePhantomP << 186     
200           stepLengths[ii] = theElossExt->TrueS << 187 #ifdef VERBOSE_ENERSPLIT
201           kinEnergyPre -= theEnergies[ii];     << 188     if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii 
202                                                << 189            << " RN: iter" << iiter << " step length geom " << stepLengths[ii] 
203 #ifdef VERBOSE_ENERSPLIT                       << 190            << " geom2true " << rnsl[ii].second / stepLengths[ii]  << G4endl;
204           if (verbose)                         << 191 #endif
205             G4cout << "G4EnergySplitter::Split << 192     
206                    << " step length geom " <<  << 193     slSum += stepLengths[ii];
207                    << rnsl[ii].second / stepLe << 194   }
208 #endif                                         << 195   
209                                                << 196   //Correct step lengths so that they sum the total step length
210           slSum += stepLengths[ii];            << 197   G4double slratio = aStep->GetStepLength()/slSum;
211         }                                      << 198 #ifdef VERBOSE_ENERSPLIT
212                                                << 199   if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
213         // Correct step lengths so that they s << 200 #endif
214         G4double slratio = aStep->GetStepLengt << 201   for( ii = 0; ii < rnsl.size(); ++ii ){
215 #ifdef VERBOSE_ENERSPLIT                       << 202     stepLengths[ii] *= slratio;
216         if (verbose)                           << 203 #ifdef VERBOSE_ENERSPLIT
217           G4cout << "G4EnergySplitter::SplitEn << 204     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
218                  << " step ratio " << slRatio  << 205 #endif
219 #endif                                         << 206   }
220         for (ii = 0; ii < rnsl.size(); ++ii) { << 207   
221           stepLengths[ii] *= slratio;          << 208   //---- Recalculate energy lost with this new step lengths
222 #ifdef VERBOSE_ENERSPLIT                       << 
223           if (verbose)                         << 
224             G4cout << "G4EnergySplitter::Split << 
225                    << " corrected step length  << 
226 #endif                                         << 
227         }                                      << 
228                                                << 
229         //---- Recalculate energy lost with th << 
230         kinEnergyPre = aStep->GetPreStepPoint(    209         kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
231         totalELost = 0.;                       << 210   totalELost = 0.;
232         for (ii = 0; ii < rnsl.size(); ++ii) { << 211   for( ii = 0; ii < rnsl.size(); ++ii ){
233           const G4Material* mate = thePhantomP << 212     const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
234           G4double dEdx = 0.;                  << 213     G4double dEdx = 0.;
235           if (kinEnergyPre > 0.) {             << 214     if( kinEnergyPre > 0. ) {
236             dEdx = emcalc.GetDEDX(kinEnergyPre << 215       dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
237           }                                    << 216     }
238           G4double elost = stepLengths[ii] * d << 217     G4double elost = stepLengths[ii] * dEdx;
239 #ifdef VERBOSE_ENERSPLIT                       << 218 #ifdef VERBOSE_ENERSPLIT
240           if (verbose)                         << 219     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost 
241             G4cout << "G4EnergySplitter::Split << 220             << " energy at interaction " << kinEnergyPre 
242                    << " energy lost " << elost << 221             << " = stepLength " << stepLengths[ii] 
243                    << " = stepLength " << step << 222             << " * dEdx " << dEdx << G4endl;
244 #endif                                         << 223 #endif
245           kinEnergyPre -= elost;               << 224     kinEnergyPre -= elost;
246           theEnergies[ii] = elost;             << 225     theEnergies[ii] = elost;
247           totalELost += elost;                 << 226     totalELost += elost;
248         }                                      << 227   }
                                                   >> 228   
249       }                                           229       }
250                                                << 230       
251       // correct energies so that they reprodu << 231       //correct energies so that they reproduce the real step energy lost
252       G4double enerRatio = (edep / totalELost) << 232       G4double enerRatio = (edep/totalELost);
253                                                << 233       
254 #ifdef VERBOSE_ENERSPLIT                       << 234 #ifdef VERBOSE_ENERSPLIT
255       if (verbose)                             << 235       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
256         G4cout << "G4EnergySplitter::SplitEner << 236 #endif
257                << " energy ratio " << enerRati << 237   
258 #endif                                         << 238 #ifdef VERBOSE_ENERSPLIT
259                                                << 239       G4double elostTot = 0.; 
260 #ifdef VERBOSE_ENERSPLIT                       << 240 #endif
261       G4double elostTot = 0.;                  << 241       for( ii = 0; ii < theEnergies.size(); ++ii ){
262 #endif                                         << 242   theEnergies[ii] *= enerRatio;
263       for (ii = 0; ii < theEnergies.size(); ++ << 243 #ifdef VERBOSE_ENERSPLIT
264         theEnergies[ii] *= enerRatio;          << 244   elostTot += theEnergies[ii];
265 #ifdef VERBOSE_ENERSPLIT                       << 245   if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii] 
266         elostTot += theEnergies[ii];           << 246           << " orig elost " << theEnergies[ii]/enerRatio 
267         if (verbose)                           << 247           << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
268           G4cout << "G4EnergySplitter::SplitEn << 248           << " energy after interaction " << kinEnergyPreOrig-elostTot
269                  << " corrected energy lost "  << 249           << G4endl;
270                  << theEnergies[ii] / enerRati << 
271                  << kinEnergyPreOrig - elostTo << 
272                  << kinEnergyPreOrig - elostTo << 
273 #endif                                            250 #endif
274       }                                           251       }
275     }                                             252     }
                                                   >> 253     
276   }                                               254   }
277                                                << 255   
278   return (G4int)theEnergies.size();               256   return (G4int)theEnergies.size();
279 }                                                 257 }
280                                                   258 
                                                   >> 259 
281 //--------------------------------------------    260 //-----------------------------------------------------------------------
282 void G4EnergySplitter::GetPhantomParam(G4bool     261 void G4EnergySplitter::GetPhantomParam(G4bool mustExist)
283 {                                                 262 {
284   G4PhysicalVolumeStore* pvs = G4PhysicalVolum    263   G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
285   for (const auto pv : *pvs) {                 << 264   for( auto cite = pvs->cbegin(); cite != pvs->cend(); ++cite ) {
286     if (IsPhantomVolume(pv)) {                 << 265     //    G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
287       const auto pvparam = static_cast<const G << 266     if( IsPhantomVolume( *cite ) ) {
                                                   >> 267       const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
288       G4VPVParameterisation* param = pvparam->    268       G4VPVParameterisation* param = pvparam->GetParameterisation();
                                                   >> 269       //    if( static_cast<const G4PhantomParameterisation*>(param) ){
                                                   >> 270       //    if( static_cast<const G4PhantomParameterisation*>(param) ){
                                                   >> 271       //      G4cout << "G4PhantomParameterisation volume found  " << (*cite)->GetName() << G4endl;
289       thePhantomParam = static_cast<G4PhantomP    272       thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
290     }                                             273     }
291   }                                               274   }
292                                                << 275   
293   if ((thePhantomParam == nullptr) && mustExis << 276   if( !thePhantomParam && mustExist )
294     G4Exception("G4EnergySplitter::GetPhantomP << 277     G4Exception("G4EnergySplitter::GetPhantomParam",
                                                   >> 278                 "PhantomParamError", FatalException,
295                 "No G4PhantomParameterisation     279                 "No G4PhantomParameterisation found !");
296 }                                                 280 }
297                                                   281 
                                                   >> 282 
298 //--------------------------------------------    283 //-----------------------------------------------------------------------
299 G4bool G4EnergySplitter::IsPhantomVolume(G4VPh << 284 G4bool G4EnergySplitter::IsPhantomVolume( G4VPhysicalVolume* pv )
300 {                                                 285 {
301   EAxis axis;                                     286   EAxis axis;
302   G4int nReplicas;                                287   G4int nReplicas;
303   G4double width, offset;                      << 288   G4double width,offset;
304   G4bool consuming;                               289   G4bool consuming;
305   pv->GetReplicationData(axis, nReplicas, widt << 290   pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
306   EVolume type = (consuming) ? kReplica : kPar    291   EVolume type = (consuming) ? kReplica : kParameterised;
307   return type == kParameterised && pv->GetRegu << 292   if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {  
308 }                                              << 293     return TRUE;
                                                   >> 294   } else {
                                                   >> 295     return FALSE;
                                                   >> 296   }
                                                   >> 297 
                                                   >> 298 } 
309                                                   299 
310 //--------------------------------------------    300 //-----------------------------------------------------------------------
311 void G4EnergySplitter::GetLastVoxelID(G4int& v << 301 void G4EnergySplitter::GetLastVoxelID( G4int& voxelID)
312 {                                              << 302 { 
313   voxelID = (*(G4RegularNavigationHelper::Inst    303   voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().cbegin())).first;
314 }                                                 304 }
315                                                   305 
316 //--------------------------------------------    306 //-----------------------------------------------------------------------
317 void G4EnergySplitter::GetFirstVoxelID(G4int&  << 307 void G4EnergySplitter::GetFirstVoxelID( G4int& voxelID)
318 {                                                 308 {
319   voxelID = (*(G4RegularNavigationHelper::Inst << 309   voxelID =  (*(G4RegularNavigationHelper::Instance()->GetStepLengths().crbegin())).first;
320 }                                                 310 }
321                                                   311 
322 //--------------------------------------------    312 //-----------------------------------------------------------------------
323 void G4EnergySplitter::GetVoxelID(G4int stepNo << 313 void G4EnergySplitter::GetVoxelID( G4int stepNo, G4int& voxelID )
324 {                                                 314 {
325   if (stepNo < 0 || stepNo >= G4int(G4RegularN << 315   if( stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()) ) {
326   {                                            << 316   G4Exception("G4EnergySplitter::GetVoxelID",
327     G4Exception("G4EnergySplitter::GetVoxelID" << 317         "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
328                 "Invalid stepNo, smaller than  << 318         FatalErrorInArgument,
329                 FatalErrorInArgument,          << 319         G4String("stepNo = " + G4UIcommand::ConvertToString(stepNo) + ", number of voxels = " + G4UIcommand::ConvertToString(G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size())) ).c_str());
330                 G4String("stepNo = " + G4UIcom << 
331                          + ", number of voxels << 
332                          + G4UIcommand::Conver << 
333                            G4int(G4RegularNavi << 
334                   .c_str());                   << 
335   }                                               320   }
336   auto ite = G4RegularNavigationHelper::Instan << 321   std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().cbegin();
337   advance(ite, stepNo);                        << 322   advance( ite, stepNo );
338   voxelID = (*ite).first;                         323   voxelID = (*ite).first;
                                                   >> 324 
339 }                                                 325 }
340                                                   326 
                                                   >> 327 
341 //--------------------------------------------    328 //-----------------------------------------------------------------------
342 void G4EnergySplitter::GetStepLength(G4int ste << 329 void G4EnergySplitter::GetStepLength( G4int stepNo, G4double& stepLength )
343 {                                                 330 {
344   auto ite = G4RegularNavigationHelper::Instan << 331   std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().cbegin();
345   advance(ite, stepNo);                        << 332   advance( ite, stepNo );
346   stepLength = (*ite).second;                     333   stepLength = (*ite).second;
347 }                                                 334 }
348                                                   335