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 9.6.p1)


  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;                          << 
 54 }                                              << 
 55                                                    52 
 56 G4int G4EnergySplitter::SplitEnergyInVolumes(c <<  53 G4int G4EnergySplitter::SplitEnergyInVolumes(const G4Step* aStep )
 57 {                                                  54 {
 58   theEnergies.clear();                             55   theEnergies.clear();
 59                                                    56 
 60   if (aStep == nullptr) return false;  // it i << 
 61                                                << 
 62   G4double edep = aStep->GetTotalEnergyDeposit     57   G4double edep = aStep->GetTotalEnergyDeposit();
 63                                                <<  58  
 64   if( edep == 0. ) {                           << 
 65     return 0;                                  << 
 66   }                                            << 
 67                                                << 
 68 #ifdef VERBOSE_ENERSPLIT                           59 #ifdef VERBOSE_ENERSPLIT
 69   G4bool verbose = 1;                              60   G4bool verbose = 1;
 70   if (verbose)                                 <<  61   if( verbose ) G4cout << "G4EnergySplitter::SplitEnergyInVolumes totalEdepo " << aStep->GetTotalEnergyDeposit() 
 71     G4cout << "G4EnergySplitter::SplitEnergyIn <<  62            << " Nsteps " << G4RegularNavigationHelper::Instance()->GetStepLengths().size() << G4endl;
 72            << " Nsteps " << G4RegularNavigatio <<  63 #endif    
 73            << G4endl;                          <<  64   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 0 ||
 74 #endif                                         <<  65       aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0)  { // we are only counting dose deposit
 75   if (G4RegularNavigationHelper::Instance()->G <<  66     return theEnergies.size();
 76       || aStep->GetTrack()->GetDefinition()->G << 
 77   {  // we are only counting dose deposit      << 
 78     return (G4int)theEnergies.size();          << 
 79   }                                                67   }
 80   if (G4RegularNavigationHelper::Instance()->G <<  68   if( G4RegularNavigationHelper::Instance()->GetStepLengths().size() == 1 ) {
 81     theEnergies.push_back(edep);                   69     theEnergies.push_back(edep);
 82     return (G4int)theEnergies.size();          <<  70     return theEnergies.size();
 83   }                                                71   }
 84                                                    72 
 85   if (thePhantomParam == nullptr) GetPhantomPa <<  73   if( !thePhantomParam ) GetPhantomParam(TRUE);
 86                                                <<  74   
 87   //----- Distribute energy deposited in voxel <<  75   if( aStep == 0 ) return FALSE; // it is 0 when called by GmScoringMgr after last event
 88   std::vector<std::pair<G4int, G4double>> rnsl <<  76   
 89     G4RegularNavigationHelper::Instance()->Get <<  77   //----- Distribute energy deposited in voxels 
                                                   >>  78   std::vector< std::pair<G4int,G4double> > rnsl = G4RegularNavigationHelper::Instance()->GetStepLengths(); 
 90                                                    79 
 91   const G4ParticleDefinition* part = aStep->Ge     80   const G4ParticleDefinition* part = aStep->GetTrack()->GetDefinition();
 92   G4double kinEnergyPreOrig = aStep->GetPreSte     81   G4double kinEnergyPreOrig = aStep->GetPreStepPoint()->GetKineticEnergy();
 93   G4double kinEnergyPre = kinEnergyPreOrig;        82   G4double kinEnergyPre = kinEnergyPreOrig;
 94                                                <<  83   
 95   G4double stepLength = aStep->GetStepLength()     84   G4double stepLength = aStep->GetStepLength();
 96   G4double slSum = 0.;                             85   G4double slSum = 0.;
 97   unsigned int ii;                                 86   unsigned int ii;
 98   for (ii = 0; ii < rnsl.size(); ++ii) {       <<  87   for( ii = 0; ii < rnsl.size(); ii++ ){
 99     G4double sl = rnsl[ii].second;                 88     G4double sl = rnsl[ii].second;
100     slSum += sl;                                   89     slSum += sl;
101 #ifdef VERBOSE_ENERSPLIT                           90 #ifdef VERBOSE_ENERSPLIT
102     if (verbose)                               <<  91     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 step length geom " << sl << G4endl;
103       G4cout << "G4EnergySplitter::SplitEnergy << 
104              << sl << G4endl;                  << 
105 #endif                                             92 #endif
106   }                                                93   }
107                                                <<  94   
108 #ifdef VERBOSE_ENERSPLIT                           95 #ifdef VERBOSE_ENERSPLIT
109   if (verbose)                                 <<  96   if( verbose )
110     G4cout << "G4EnergySplitter RN:  step leng <<  97     G4cout << "G4EnergySplitter RN:  step length geom TOTAL " << slSum 
111            << stepLength << " ratio " << stepL <<  98      << " true TOTAL " << stepLength 
112            << aStep->GetPreStepPoint()->GetKin <<  99      << " ratio " << stepLength/slSum 
113            << aStep->GetPreStepPoint()->GetMat << 100      << " Energy " << aStep->GetPreStepPoint()->GetKineticEnergy() 
114            << rnsl.size() << G4endl;           << 101      << " Material " << aStep->GetPreStepPoint()->GetMaterial()->GetName() 
115 #endif                                         << 102      << " Number of geom steps " << rnsl.size() << G4endl;
116   //----- No iterations to correct elost and m << 103 #endif
117   // geometrical step length in each voxel     << 104   //----- No iterations to correct elost and msc => distribute energy deposited according to geometrical step length in each voxel
118   if (theNIterations == 0) {                   << 105   if( theNIterations == 0 ) { 
119     for (ii = 0; ii < rnsl.size(); ++ii) {     << 106     for( ii = 0; ii < rnsl.size(); ii++ ){
120       G4double sl = rnsl[ii].second;              107       G4double sl = rnsl[ii].second;
121       G4double edepStep = edep * sl / slSum;   << 108       G4double edepStep = edep * sl/slSum; //divide edep along steps, proportional to step length
122                                                << 
123 #ifdef VERBOSE_ENERSPLIT                          109 #ifdef VERBOSE_ENERSPLIT
124       if (verbose)                             << 110       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii 
125         G4cout << "G4EnergySplitter::SplitEner << 111         << " edep " << edepStep << G4endl;
126 #endif                                            112 #endif
127                                                << 113   
128       theEnergies.push_back(edepStep);            114       theEnergies.push_back(edepStep);
                                                   >> 115   
129     }                                             116     }
130   }                                            << 117   } else { //  1 or more iterations demanded
131   else {  //  1 or more iterations demanded    << 118       
132                                                << 
133 #ifdef VERBOSE_ENERSPLIT                          119 #ifdef VERBOSE_ENERSPLIT
134     // print corrected energy at iteration 0   << 120     // print corrected energy at iteration 0 
135     if (verbose) {                             << 121     if(verbose)  {
136       G4double slSum = 0.;                        122       G4double slSum = 0.;
137       for (ii = 0; ii < rnsl.size(); ++ii) {   << 123       for( ii = 0; ii < rnsl.size(); ii++ ){
138         G4double sl = rnsl[ii].second;         << 124   G4double sl = rnsl[ii].second;
139         slSum += sl;                           << 125   slSum += sl;
140       }                                           126       }
141       for (ii = 0; ii < rnsl.size(); ii++) {   << 127       for( ii = 0; ii < rnsl.size(); ii++ ){
142         G4cout << "G4EnergySplitter::SplitEner << 128   G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii
143                << " RN: iter0 corrected energy << 129     << " RN: iter0 corrected energy lost " << edep*rnsl[ii].second/slSum  
                                                   >> 130     << G4endl;
144       }                                           131       }
145     }                                             132     }
146 #endif                                            133 #endif
147                                                   134 
148     G4double slRatio = stepLength / slSum;     << 135     G4double slRatio = stepLength/slSum;
149 #ifdef VERBOSE_ENERSPLIT                          136 #ifdef VERBOSE_ENERSPLIT
150     if (verbose)                               << 137     if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes  RN: iter 0, step ratio " << slRatio << G4endl;
151       G4cout << "G4EnergySplitter::SplitEnergy << 
152              << G4endl;                        << 
153 #endif                                            138 #endif
154                                                << 139       
155     //--- energy at each interaction              140     //--- energy at each interaction
156     G4EmCalculator emcalc;                        141     G4EmCalculator emcalc;
157     G4double totalELost = 0.;                     142     G4double totalELost = 0.;
158     std::vector<G4double> stepLengths;            143     std::vector<G4double> stepLengths;
159     for (G4int iiter = 1; iiter <= theNIterati << 144     for( int iiter = 1; iiter <= theNIterations; iiter++ ) {
160       //--- iter1: distribute true step length << 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
161       // a constant so that the sum gives the  << 146       if( iiter == 1 ) {
162       if (iiter == 1) {                        << 147   for( ii = 0; ii < rnsl.size(); ii++ ){
163         for (ii = 0; ii < rnsl.size(); ++ii) { << 148     G4double sl = rnsl[ii].second;
164           G4double sl = rnsl[ii].second;       << 149     stepLengths.push_back( sl * slRatio );
165           stepLengths.push_back(sl * slRatio); << 150 #ifdef VERBOSE_ENERSPLIT
166 #ifdef VERBOSE_ENERSPLIT                       << 151     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << sl*slRatio << G4endl;
167           if (verbose)                         << 152 #endif
168             G4cout << "G4EnergySplitter::Split << 153   }
169                    << " corrected step length  << 154   
170 #endif                                         << 155   for( ii = 0; ii < rnsl.size(); ii++ ){
171         }                                      << 156     const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
172                                                << 157     G4double dEdx = 0.;
173         for (ii = 0; ii < rnsl.size(); ++ii) { << 158     if( kinEnergyPre > 0. ) {  //t check this 
174           const G4Material* mate = thePhantomP << 159       dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
175           G4double dEdx = 0.;                  << 160     }
176           if (kinEnergyPre > 0.) {  // t check << 161     G4double elost = stepLengths[ii] * dEdx;
177             dEdx = emcalc.GetDEDX(kinEnergyPre << 162     
178           }                                    << 163 #ifdef VERBOSE_ENERSPLIT
179           G4double elost = stepLengths[ii] * d << 164     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter1 energy lost "  << elost 
180                                                << 165             << " energy at interaction " << kinEnergyPre 
181 #ifdef VERBOSE_ENERSPLIT                       << 166             << " = stepLength " << stepLengths[ii] 
182           if (verbose)                         << 167             << " * dEdx " << dEdx << G4endl;
183             G4cout << "G4EnergySplitter::Split << 168 #endif
184                    << elost << " energy at int << 169     kinEnergyPre -= elost;
185                    << stepLengths[ii] << " * d << 170     theEnergies.push_back( elost );
186 #endif                                         << 171     totalELost += elost;
187           kinEnergyPre -= elost;               << 172   }
188           theEnergies.push_back(elost);        << 173   
189           totalELost += elost;                 << 174       } else{
190         }                                      << 175   //------ 2nd and other iterations
191       }                                        << 176   //----- Get step lengths corrected by changing geom2true correction
192       else {                                   << 177   //-- Get ratios for each energy 
193         //------ 2nd and other iterations      << 178   slSum = 0.;
194         //----- Get step lengths corrected by  << 179   kinEnergyPre = kinEnergyPreOrig;
195         //-- Get ratios for each energy        << 180   for( ii = 0; ii < rnsl.size(); ii++ ){
196         slSum = 0.;                            << 181     const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
197         kinEnergyPre = kinEnergyPreOrig;       << 182     stepLengths[ii] = theElossExt->TrueStepLength( kinEnergyPre, rnsl[ii].second , mate, part );
198         for (ii = 0; ii < rnsl.size(); ++ii) { << 183     kinEnergyPre -= theEnergies[ii];
199           const G4Material* mate = thePhantomP << 184     
200           stepLengths[ii] = theElossExt->TrueS << 185 #ifdef VERBOSE_ENERSPLIT
201           kinEnergyPre -= theEnergies[ii];     << 186     if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii 
202                                                << 187            << " RN: iter" << iiter << " step length geom " << stepLengths[ii] 
203 #ifdef VERBOSE_ENERSPLIT                       << 188            << " geom2true " << rnsl[ii].second / stepLengths[ii]  << G4endl;
204           if (verbose)                         << 189 #endif
205             G4cout << "G4EnergySplitter::Split << 190     
206                    << " step length geom " <<  << 191     slSum += stepLengths[ii];
207                    << rnsl[ii].second / stepLe << 192   }
208 #endif                                         << 193   
209                                                << 194   //Correct step lengths so that they sum the total step length
210           slSum += stepLengths[ii];            << 195   G4double slratio = aStep->GetStepLength()/slSum;
211         }                                      << 196 #ifdef VERBOSE_ENERSPLIT
212                                                << 197   if(verbose) G4cout << "G4EnergySplitter::SplitEnergyInVolumes" << ii << " RN: iter" << iiter << " step ratio " << slRatio << G4endl;
213         // Correct step lengths so that they s << 198 #endif
214         G4double slratio = aStep->GetStepLengt << 199   for( ii = 0; ii < rnsl.size(); ii++ ){
215 #ifdef VERBOSE_ENERSPLIT                       << 200     stepLengths[ii] *= slratio;
216         if (verbose)                           << 201 #ifdef VERBOSE_ENERSPLIT
217           G4cout << "G4EnergySplitter::SplitEn << 202     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " corrected step length " << stepLengths[ii] << G4endl;
218                  << " step ratio " << slRatio  << 203 #endif
219 #endif                                         << 204   }
220         for (ii = 0; ii < rnsl.size(); ++ii) { << 205   
221           stepLengths[ii] *= slratio;          << 206   //---- 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(    207         kinEnergyPre = aStep->GetPreStepPoint()->GetKineticEnergy();
231         totalELost = 0.;                       << 208   totalELost = 0.;
232         for (ii = 0; ii < rnsl.size(); ++ii) { << 209   for( ii = 0; ii < rnsl.size(); ii++ ){
233           const G4Material* mate = thePhantomP << 210     const G4Material* mate = thePhantomParam->GetMaterial( rnsl[ii].first );
234           G4double dEdx = 0.;                  << 211     G4double dEdx = 0.;
235           if (kinEnergyPre > 0.) {             << 212     if( kinEnergyPre > 0. ) {
236             dEdx = emcalc.GetDEDX(kinEnergyPre << 213       dEdx = emcalc.GetDEDX(kinEnergyPre, part, mate);
237           }                                    << 214     }
238           G4double elost = stepLengths[ii] * d << 215     G4double elost = stepLengths[ii] * dEdx;
239 #ifdef VERBOSE_ENERSPLIT                       << 216 #ifdef VERBOSE_ENERSPLIT
240           if (verbose)                         << 217     if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy lost " << elost 
241             G4cout << "G4EnergySplitter::Split << 218             << " energy at interaction " << kinEnergyPre 
242                    << " energy lost " << elost << 219             << " = stepLength " << stepLengths[ii] 
243                    << " = stepLength " << step << 220             << " * dEdx " << dEdx << G4endl;
244 #endif                                         << 221 #endif
245           kinEnergyPre -= elost;               << 222     kinEnergyPre -= elost;
246           theEnergies[ii] = elost;             << 223     theEnergies[ii] = elost;
247           totalELost += elost;                 << 224     totalELost += elost;
248         }                                      << 225   }
                                                   >> 226   
249       }                                           227       }
250                                                << 228       
251       // correct energies so that they reprodu << 229       //correct energies so that they reproduce the real step energy lost
252       G4double enerRatio = (edep / totalELost) << 230       G4double enerRatio = (edep/totalELost);
253                                                << 231       
254 #ifdef VERBOSE_ENERSPLIT                       << 232 #ifdef VERBOSE_ENERSPLIT
255       if (verbose)                             << 233       if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes"<< ii << " RN: iter" << iiter << " energy ratio " << enerRatio << G4endl;
256         G4cout << "G4EnergySplitter::SplitEner << 234 #endif
257                << " energy ratio " << enerRati << 235   
258 #endif                                         << 236 #ifdef VERBOSE_ENERSPLIT
259                                                << 237       G4double elostTot = 0.; 
260 #ifdef VERBOSE_ENERSPLIT                       << 238 #endif
261       G4double elostTot = 0.;                  << 239       for( ii = 0; ii < theEnergies.size(); ii++ ){
262 #endif                                         << 240   theEnergies[ii] *= enerRatio;
263       for (ii = 0; ii < theEnergies.size(); ++ << 241 #ifdef VERBOSE_ENERSPLIT
264         theEnergies[ii] *= enerRatio;          << 242   elostTot += theEnergies[ii];
265 #ifdef VERBOSE_ENERSPLIT                       << 243   if(verbose) G4cout  << "G4EnergySplitter::SplitEnergyInVolumes "<< ii << " RN: iter" << iiter << " corrected energy lost " << theEnergies[ii] 
266         elostTot += theEnergies[ii];           << 244           << " orig elost " << theEnergies[ii]/enerRatio 
267         if (verbose)                           << 245           << " energy before interaction " << kinEnergyPreOrig-elostTot+theEnergies[ii]
268           G4cout << "G4EnergySplitter::SplitEn << 246           << " energy after interaction " << kinEnergyPreOrig-elostTot
269                  << " corrected energy lost "  << 247           << G4endl;
270                  << theEnergies[ii] / enerRati << 
271                  << kinEnergyPreOrig - elostTo << 
272                  << kinEnergyPreOrig - elostTo << 
273 #endif                                            248 #endif
274       }                                           249       }
275     }                                             250     }
                                                   >> 251     
276   }                                               252   }
277                                                << 253   
278   return (G4int)theEnergies.size();            << 254   return theEnergies.size();
279 }                                                 255 }
280                                                   256 
                                                   >> 257 
281 //--------------------------------------------    258 //-----------------------------------------------------------------------
282 void G4EnergySplitter::GetPhantomParam(G4bool     259 void G4EnergySplitter::GetPhantomParam(G4bool mustExist)
283 {                                                 260 {
284   G4PhysicalVolumeStore* pvs = G4PhysicalVolum    261   G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
285   for (const auto pv : *pvs) {                 << 262   std::vector<G4VPhysicalVolume*>::iterator cite;
286     if (IsPhantomVolume(pv)) {                 << 263   for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
287       const auto pvparam = static_cast<const G << 264     //    G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
                                                   >> 265     if( IsPhantomVolume( *cite ) ) {
                                                   >> 266       const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
288       G4VPVParameterisation* param = pvparam->    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;
289       thePhantomParam = static_cast<G4PhantomP    271       thePhantomParam = static_cast<G4PhantomParameterisation*>(param);
290     }                                             272     }
291   }                                               273   }
292                                                << 274   
293   if ((thePhantomParam == nullptr) && mustExis << 275   if( !thePhantomParam && mustExist )
294     G4Exception("G4EnergySplitter::GetPhantomP << 276     G4Exception("G4EnergySplitter::GetPhantomParam",
                                                   >> 277                 "PhantomParamError", FatalException,
295                 "No G4PhantomParameterisation     278                 "No G4PhantomParameterisation found !");
296 }                                                 279 }
297                                                   280 
                                                   >> 281 
298 //--------------------------------------------    282 //-----------------------------------------------------------------------
299 G4bool G4EnergySplitter::IsPhantomVolume(G4VPh << 283 G4bool G4EnergySplitter::IsPhantomVolume( G4VPhysicalVolume* pv )
300 {                                                 284 {
301   EAxis axis;                                     285   EAxis axis;
302   G4int nReplicas;                                286   G4int nReplicas;
303   G4double width, offset;                      << 287   G4double width,offset;
304   G4bool consuming;                               288   G4bool consuming;
305   pv->GetReplicationData(axis, nReplicas, widt << 289   pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
306   EVolume type = (consuming) ? kReplica : kPar    290   EVolume type = (consuming) ? kReplica : kParameterised;
307   return type == kParameterised && pv->GetRegu << 291   if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {  
308 }                                              << 292     return TRUE;
                                                   >> 293   } else {
                                                   >> 294     return FALSE;
                                                   >> 295   }
                                                   >> 296 
                                                   >> 297 } 
309                                                   298 
310 //--------------------------------------------    299 //-----------------------------------------------------------------------
311 void G4EnergySplitter::GetLastVoxelID(G4int& v << 300 void G4EnergySplitter::GetLastVoxelID( G4int& voxelID)
312 {                                              << 301 { 
313   voxelID = (*(G4RegularNavigationHelper::Inst << 302   voxelID = (*(G4RegularNavigationHelper::Instance()->GetStepLengths().begin())).first;
314 }                                                 303 }
315                                                   304 
316 //--------------------------------------------    305 //-----------------------------------------------------------------------
317 void G4EnergySplitter::GetFirstVoxelID(G4int&  << 306 void G4EnergySplitter::GetFirstVoxelID( G4int& voxelID)
318 {                                                 307 {
319   voxelID = (*(G4RegularNavigationHelper::Inst << 308   voxelID =  (*(G4RegularNavigationHelper::Instance()->GetStepLengths().rbegin())).first;
320 }                                                 309 }
321                                                   310 
322 //--------------------------------------------    311 //-----------------------------------------------------------------------
323 void G4EnergySplitter::GetVoxelID(G4int stepNo << 312 void G4EnergySplitter::GetVoxelID( G4int stepNo, G4int& voxelID )
324 {                                                 313 {
325   if (stepNo < 0 || stepNo >= G4int(G4RegularN << 314   if( stepNo < 0 || stepNo >= G4int(G4RegularNavigationHelper::Instance()->GetStepLengths().size()) ) {
326   {                                            << 315   G4Exception("G4EnergySplitter::GetVoxelID",
327     G4Exception("G4EnergySplitter::GetVoxelID" << 316         "Invalid stepNo, smaller than 0 or bigger or equal to number of voxels traversed",
328                 "Invalid stepNo, smaller than  << 317         FatalErrorInArgument,
329                 FatalErrorInArgument,          << 318         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   }                                               319   }
336   auto ite = G4RegularNavigationHelper::Instan << 320   std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
337   advance(ite, stepNo);                        << 321   advance( ite, stepNo );
338   voxelID = (*ite).first;                         322   voxelID = (*ite).first;
                                                   >> 323 
339 }                                                 324 }
340                                                   325 
                                                   >> 326 
341 //--------------------------------------------    327 //-----------------------------------------------------------------------
342 void G4EnergySplitter::GetStepLength(G4int ste << 328 void G4EnergySplitter::GetStepLength( G4int stepNo, G4double& stepLength )
343 {                                                 329 {
344   auto ite = G4RegularNavigationHelper::Instan << 330   std::vector< std::pair<G4int,G4double> >::const_iterator ite = G4RegularNavigationHelper::Instance()->GetStepLengths().begin();
345   advance(ite, stepNo);                        << 331   advance( ite, stepNo );
346   stepLength = (*ite).second;                     332   stepLength = (*ite).second;
347 }                                                 333 }
348                                                   334