Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4EmBiasingManager.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/electromagnetic/utils/src/G4EmBiasingManager.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4EmBiasingManager.cc (Version 10.5)


  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 //                                                 26 //
 27 // -------------------------------------------     27 // -------------------------------------------------------------------
 28 //                                                 28 //
 29 // GEANT4 Class file                               29 // GEANT4 Class file
 30 //                                                 30 //
 31 //                                                 31 //
 32 // File name:     G4EmBiasingManager               32 // File name:     G4EmBiasingManager
 33 //                                                 33 //
 34 // Author:        Vladimir Ivanchenko              34 // Author:        Vladimir Ivanchenko 
 35 //                                                 35 //
 36 // Creation date: 28.07.2011                       36 // Creation date: 28.07.2011
 37 //                                                 37 //
 38 // Modifications:                                  38 // Modifications:
 39 //                                                 39 //
 40 // 31-05-12 D. Sawkey put back in high energy      40 // 31-05-12 D. Sawkey put back in high energy limit for brem, russian roulette 
 41 // 30-05-12 D. Sawkey  brem split gammas are u     41 // 30-05-12 D. Sawkey  brem split gammas are unique; do weight tests for 
 42 //          brem, russian roulette                 42 //          brem, russian roulette
 43 // -------------------------------------------     43 // -------------------------------------------------------------------
 44 //                                                 44 //
 45 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 46 //....oooOO0OOooo........oooOO0OOooo........oo     46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 47                                                    47 
 48 #include "G4EmBiasingManager.hh"                   48 #include "G4EmBiasingManager.hh"
 49 #include "G4SystemOfUnits.hh"                      49 #include "G4SystemOfUnits.hh"
 50 #include "G4PhysicalConstants.hh"                  50 #include "G4PhysicalConstants.hh"
 51 #include "G4MaterialCutsCouple.hh"                 51 #include "G4MaterialCutsCouple.hh"
 52 #include "G4ProductionCutsTable.hh"                52 #include "G4ProductionCutsTable.hh"
 53 #include "G4ProductionCuts.hh"                     53 #include "G4ProductionCuts.hh"
 54 #include "G4Region.hh"                             54 #include "G4Region.hh"
 55 #include "G4RegionStore.hh"                        55 #include "G4RegionStore.hh"
 56 #include "G4Track.hh"                              56 #include "G4Track.hh"
 57 #include "G4Electron.hh"                           57 #include "G4Electron.hh"
 58 #include "G4Gamma.hh"                              58 #include "G4Gamma.hh"
 59 #include "G4VEmModel.hh"                           59 #include "G4VEmModel.hh"
 60 #include "G4LossTableManager.hh"                   60 #include "G4LossTableManager.hh"
 61 #include "G4ParticleChangeForLoss.hh"              61 #include "G4ParticleChangeForLoss.hh"
 62 #include "G4ParticleChangeForGamma.hh"             62 #include "G4ParticleChangeForGamma.hh"
 63 #include "G4EmParameters.hh"                       63 #include "G4EmParameters.hh"
 64                                                    64 
 65 //....oooOO0OOooo........oooOO0OOooo........oo     65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 66                                                    66 
 67 G4EmBiasingManager::G4EmBiasingManager()       <<  67 G4EmBiasingManager::G4EmBiasingManager() 
 68   : fDirectionalSplittingTarget(0.0,0.0,0.0)   <<  68   : nForcedRegions(0),nSecBiasedRegions(0),eIonisation(nullptr),
                                                   >>  69     currentStepLimit(0.0),startTracking(true)
 69 {                                                  70 {
 70   fSafetyMin = 1.e-6*mm;                           71   fSafetyMin = 1.e-6*mm;
 71   theElectron = G4Electron::Electron();            72   theElectron = G4Electron::Electron();
 72   theGamma    = G4Gamma::Gamma();                  73   theGamma    = G4Gamma::Gamma();
                                                   >>  74 
                                                   >>  75   fDirectionalSplitting = false;
                                                   >>  76   fDirectionalSplittingRadius = 0.;
                                                   >>  77   fDirectionalSplittingTarget = G4ThreeVector(0.,0.,0.);
                                                   >>  78   fDirectionalSplittingWeights.clear();
                                                   >>  79   fWeight = 1.;
 73 }                                                  80 }
 74                                                    81 
 75 //....oooOO0OOooo........oooOO0OOooo........oo     82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 76                                                    83 
 77 G4EmBiasingManager::~G4EmBiasingManager() = de <<  84 G4EmBiasingManager::~G4EmBiasingManager()
                                                   >>  85 {}
 78                                                    86 
 79 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 80                                                    88 
 81 void G4EmBiasingManager::Initialise(const G4Pa     89 void G4EmBiasingManager::Initialise(const G4ParticleDefinition& part,
 82                                     const G4St     90                                     const G4String& procName, G4int verbose)
 83 {                                                  91 {
 84   //G4cout << "G4EmBiasingManager::Initialise      92   //G4cout << "G4EmBiasingManager::Initialise for "
 85   //         << part.GetParticleName()             93   //         << part.GetParticleName()
 86   //         << " and " << procName << G4endl;     94   //         << " and " << procName << G4endl;
 87   const G4ProductionCutsTable* theCoupleTable=     95   const G4ProductionCutsTable* theCoupleTable=
 88     G4ProductionCutsTable::GetProductionCutsTa     96     G4ProductionCutsTable::GetProductionCutsTable();
 89   G4int numOfCouples = (G4int)theCoupleTable-> <<  97   size_t numOfCouples = theCoupleTable->GetTableSize();
 90                                                    98 
 91   if(0 < nForcedRegions) { idxForcedCouple.res     99   if(0 < nForcedRegions) { idxForcedCouple.resize(numOfCouples, -1); }
 92   if(0 < nSecBiasedRegions) { idxSecBiasedCoup    100   if(0 < nSecBiasedRegions) { idxSecBiasedCouple.resize(numOfCouples, -1); }
 93                                                   101 
 94   // Deexcitation                                 102   // Deexcitation
 95   for (G4int j=0; j<numOfCouples; ++j) {       << 103   for (size_t j=0; j<numOfCouples; ++j) {
 96     const G4MaterialCutsCouple* couple =          104     const G4MaterialCutsCouple* couple =
 97       theCoupleTable->GetMaterialCutsCouple(j)    105       theCoupleTable->GetMaterialCutsCouple(j);
 98     const G4ProductionCuts* pcuts = couple->Ge    106     const G4ProductionCuts* pcuts = couple->GetProductionCuts();
 99     if(0 <  nForcedRegions) {                     107     if(0 <  nForcedRegions) {
100       for(G4int i=0; i<nForcedRegions; ++i) {     108       for(G4int i=0; i<nForcedRegions; ++i) {
101         if(forcedRegions[i]) {                    109         if(forcedRegions[i]) {
102           if(pcuts == forcedRegions[i]->GetPro    110           if(pcuts == forcedRegions[i]->GetProductionCuts()) { 
103             idxForcedCouple[j] = i;               111             idxForcedCouple[j] = i;
104             break;                                112             break; 
105           }                                       113           }
106         }                                         114         }
107       }                                           115       }
108     }                                             116     }
109     if(0 < nSecBiasedRegions) {                   117     if(0 < nSecBiasedRegions) { 
110       for(G4int i=0; i<nSecBiasedRegions; ++i)    118       for(G4int i=0; i<nSecBiasedRegions; ++i) {
111         if(secBiasedRegions[i]) {                 119         if(secBiasedRegions[i]) {
112           if(pcuts == secBiasedRegions[i]->Get    120           if(pcuts == secBiasedRegions[i]->GetProductionCuts()) { 
113             idxSecBiasedCouple[j] = i;            121             idxSecBiasedCouple[j] = i;
114             break;                                122             break; 
115           }                                       123           }
116         }                                         124         }
117       }                                           125       }
118     }                                             126     }
119   }                                               127   }
120                                                   128 
121   G4EmParameters* param = G4EmParameters::Inst    129   G4EmParameters* param = G4EmParameters::Instance();
122   SetDirectionalSplitting(param->GetDirectiona    130   SetDirectionalSplitting(param->GetDirectionalSplitting());
123   if (fDirectionalSplitting) {                    131   if (fDirectionalSplitting) {
124     SetDirectionalSplittingTarget(param->GetDi    132     SetDirectionalSplittingTarget(param->GetDirectionalSplittingTarget());
125     SetDirectionalSplittingRadius(param->GetDi    133     SetDirectionalSplittingRadius(param->GetDirectionalSplittingRadius());
126   }                                               134   }
127                                                   135 
128   if (nForcedRegions > 0 && 0 < verbose) {        136   if (nForcedRegions > 0 && 0 < verbose) {
129     G4cout << " Forced Interaction is activate    137     G4cout << " Forced Interaction is activated for "
130            << part.GetParticleName() << " and     138            << part.GetParticleName() << " and " 
131            << procName                            139            << procName 
132            << " inside G4Regions: " << G4endl;    140            << " inside G4Regions: " << G4endl;
133     for (G4int i=0; i<nForcedRegions; ++i) {      141     for (G4int i=0; i<nForcedRegions; ++i) {
134       const G4Region* r = forcedRegions[i];       142       const G4Region* r = forcedRegions[i];
135       if(r) { G4cout << "           " << r->Ge    143       if(r) { G4cout << "           " << r->GetName() << G4endl; }
136     }                                             144     }
137   }                                               145   }
138   if (nSecBiasedRegions > 0 && 0 < verbose) {     146   if (nSecBiasedRegions > 0 && 0 < verbose) {
139     G4cout << " Secondary biasing is activated    147     G4cout << " Secondary biasing is activated for " 
140            << part.GetParticleName() << " and     148            << part.GetParticleName() << " and " 
141            << procName                            149            << procName 
142            << " inside G4Regions: " << G4endl;    150            << " inside G4Regions: " << G4endl;
143     for (G4int i=0; i<nSecBiasedRegions; ++i)     151     for (G4int i=0; i<nSecBiasedRegions; ++i) {
144       const G4Region* r = secBiasedRegions[i];    152       const G4Region* r = secBiasedRegions[i];
145       if(r) {                                     153       if(r) { 
146         G4cout << "           " << r->GetName(    154         G4cout << "           " << r->GetName() 
147                << "  BiasingWeight= " << secBi    155                << "  BiasingWeight= " << secBiasedWeight[i] << G4endl; 
148       }                                           156       }
149     }                                             157     }
150     if (fDirectionalSplitting) {                  158     if (fDirectionalSplitting) {
151       G4cout << "     Directional splitting ac    159       G4cout << "     Directional splitting activated, with target position: "
152              << fDirectionalSplittingTarget/cm    160              << fDirectionalSplittingTarget/cm
153              << " cm; radius: "                   161              << " cm; radius: "
154              << fDirectionalSplittingRadius/cm    162              << fDirectionalSplittingRadius/cm
155              << "cm." << G4endl;                  163              << "cm." << G4endl;
156     }                                             164     }
157   }                                               165   }
158 }                                                 166 }
159                                                   167 
160 //....oooOO0OOooo........oooOO0OOooo........oo    168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161                                                   169 
162 void G4EmBiasingManager::ActivateForcedInterac    170 void G4EmBiasingManager::ActivateForcedInteraction(G4double val, 
163                                                   171                                                    const G4String& rname)
164 {                                                 172 {
165   G4RegionStore* regionStore = G4RegionStore::    173   G4RegionStore* regionStore = G4RegionStore::GetInstance();
166   G4String name = rname;                          174   G4String name = rname;
167   if(name == "" || name == "world" || name ==     175   if(name == "" || name == "world" || name == "World") {
168     name = "DefaultRegionForTheWorld";            176     name = "DefaultRegionForTheWorld";
169   }                                               177   }
170   const G4Region* reg = regionStore->GetRegion    178   const G4Region* reg = regionStore->GetRegion(name, false);
171   if(!reg) {                                      179   if(!reg) { 
172     G4cout << "### G4EmBiasingManager::ForcedI    180     G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
173            << " G4Region <"                       181            << " G4Region <"
174            << rname << "> is unknown" << G4end    182            << rname << "> is unknown" << G4endl;
175     return;                                       183     return; 
176   }                                               184   }
177                                                   185 
178   // the region is in the list                    186   // the region is in the list
179   if (0 < nForcedRegions) {                       187   if (0 < nForcedRegions) {
180     for (G4int i=0; i<nForcedRegions; ++i) {      188     for (G4int i=0; i<nForcedRegions; ++i) {
181       if (reg == forcedRegions[i]) {              189       if (reg == forcedRegions[i]) {
182         lengthForRegion[i] = val;                 190         lengthForRegion[i] = val; 
183         return;                                   191         return;
184       }                                           192       }
185     }                                             193     }
186   }                                               194   }
187   if(val < 0.0) {                                 195   if(val < 0.0) { 
188     G4cout << "### G4EmBiasingManager::ForcedI    196     G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
189            << val << " < 0.0, so no activation    197            << val << " < 0.0, so no activation for the G4Region <"
190            << rname << ">" << G4endl;             198            << rname << ">" << G4endl;
191     return;                                       199     return; 
192   }                                               200   }
193                                                   201 
194   // new region                                   202   // new region 
195   forcedRegions.push_back(reg);                   203   forcedRegions.push_back(reg);
196   lengthForRegion.push_back(val);                 204   lengthForRegion.push_back(val);
197   ++nForcedRegions;                               205   ++nForcedRegions;
198 }                                                 206 }
199                                                   207 
200 //....oooOO0OOooo........oooOO0OOooo........oo    208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201                                                   209 
202 void                                              210 void 
203 G4EmBiasingManager::ActivateSecondaryBiasing(c    211 G4EmBiasingManager::ActivateSecondaryBiasing(const G4String& rname, 
204                                              G    212                                              G4double factor,
205                                              G    213                                              G4double energyLimit)
206 {                                                 214 {
207   //G4cout << "G4EmBiasingManager::ActivateSec    215   //G4cout << "G4EmBiasingManager::ActivateSecondaryBiasing: "
208   //         << rname << " F= " << factor << "    216   //         << rname << " F= " << factor << " E(MeV)= " << energyLimit/MeV
209   //         << G4endl;                           217   //         << G4endl; 
210   G4RegionStore* regionStore = G4RegionStore::    218   G4RegionStore* regionStore = G4RegionStore::GetInstance();
211   G4String name = rname;                          219   G4String name = rname;
212   if(name == "" || name == "world" || name ==     220   if(name == "" || name == "world" || name == "World") {
213     name = "DefaultRegionForTheWorld";            221     name = "DefaultRegionForTheWorld";
214   }                                               222   }
215   const G4Region* reg = regionStore->GetRegion    223   const G4Region* reg = regionStore->GetRegion(name, false);
216   if(!reg) {                                      224   if(!reg) { 
217     G4cout << "### G4EmBiasingManager::Activat    225     G4cout << "### G4EmBiasingManager::ActivateBremsstrahlungSplitting "
218            << "WARNING: G4Region <"               226            << "WARNING: G4Region <"
219            << rname << "> is unknown" << G4end    227            << rname << "> is unknown" << G4endl;
220     return;                                       228     return; 
221   }                                               229   }
222                                                   230 
223   // Range cut                                    231   // Range cut
224   G4int nsplit = 0;                               232   G4int nsplit = 0;
225   G4double w = factor;                            233   G4double w = factor;
226                                                   234 
227   // splitting                                    235   // splitting
228   if(factor >= 1.0) {                             236   if(factor >= 1.0) {
229     nsplit = G4lrint(factor);                     237     nsplit = G4lrint(factor);
230     w = 1.0/G4double(nsplit);                     238     w = 1.0/G4double(nsplit);
231                                                   239 
232     // Russian roulette                           240     // Russian roulette 
233   } else if(0.0 < factor) {                       241   } else if(0.0 < factor) { 
234     nsplit = 1;                                   242     nsplit = 1;
235     w = 1.0/factor;                               243     w = 1.0/factor;
236   }                                               244   }
237                                                   245 
238   // the region is in the list - overwrite par    246   // the region is in the list - overwrite parameters
239   if (0 < nSecBiasedRegions) {                    247   if (0 < nSecBiasedRegions) {
240     for (G4int i=0; i<nSecBiasedRegions; ++i)     248     for (G4int i=0; i<nSecBiasedRegions; ++i) {
241       if (reg == secBiasedRegions[i]) {           249       if (reg == secBiasedRegions[i]) {
242         secBiasedWeight[i] = w;                   250         secBiasedWeight[i] = w;
243         nBremSplitting[i]  = nsplit;              251         nBremSplitting[i]  = nsplit; 
244         secBiasedEnegryLimit[i] = energyLimit;    252         secBiasedEnegryLimit[i] = energyLimit;
245         return;                                   253         return;
246       }                                           254       }
247     }                                             255     }
248   }                                               256   }
249   /*                                              257   /*
250     G4cout << "### G4EmBiasingManager::Activat    258     G4cout << "### G4EmBiasingManager::ActivateSecondaryBiasing: "
251            << " nsplit= " << nsplit << " for t    259            << " nsplit= " << nsplit << " for the G4Region <"
252            << rname << ">" << G4endl;             260            << rname << ">" << G4endl; 
253   */                                              261   */
254                                                   262 
255   // new region                                   263   // new region 
256   secBiasedRegions.push_back(reg);                264   secBiasedRegions.push_back(reg);
257   secBiasedWeight.push_back(w);                   265   secBiasedWeight.push_back(w);
258   nBremSplitting.push_back(nsplit);               266   nBremSplitting.push_back(nsplit);
259   secBiasedEnegryLimit.push_back(energyLimit);    267   secBiasedEnegryLimit.push_back(energyLimit);
260   ++nSecBiasedRegions;                            268   ++nSecBiasedRegions;
261   //G4cout << "nSecBiasedRegions= " << nSecBia    269   //G4cout << "nSecBiasedRegions= " << nSecBiasedRegions << G4endl;
262 }                                                 270 }
263                                                   271 
264 //....oooOO0OOooo........oooOO0OOooo........oo    272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265                                                   273 
266 G4double G4EmBiasingManager::GetStepLimit(G4in    274 G4double G4EmBiasingManager::GetStepLimit(G4int coupleIdx, 
267                                           G4do    275                                           G4double previousStep)
268 {                                                 276 {
269   if(startTracking) {                             277   if(startTracking) {
270     startTracking = false;                        278     startTracking = false;
271     G4int i = idxForcedCouple[coupleIdx];         279     G4int i = idxForcedCouple[coupleIdx];
272     if(i < 0) {                                   280     if(i < 0) {
273       currentStepLimit = DBL_MAX;                 281       currentStepLimit = DBL_MAX;
274     } else {                                      282     } else {
275       currentStepLimit = lengthForRegion[i];      283       currentStepLimit = lengthForRegion[i];
276       if(currentStepLimit > 0.0) { currentStep    284       if(currentStepLimit > 0.0) { currentStepLimit *= G4UniformRand(); }
277     }                                             285     }
278   } else {                                        286   } else {
279     currentStepLimit -= previousStep;             287     currentStepLimit -= previousStep;
280   }                                               288   }
281   if(currentStepLimit < 0.0) { currentStepLimi    289   if(currentStepLimit < 0.0) { currentStepLimit = 0.0; }
282   return currentStepLimit;                        290   return currentStepLimit;
283 }                                                 291 }
284                                                   292 
285 //....oooOO0OOooo........oooOO0OOooo........oo    293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286                                                   294 
287 G4double                                          295 G4double 
288 G4EmBiasingManager::ApplySecondaryBiasing(        296 G4EmBiasingManager::ApplySecondaryBiasing(
289                     std::vector<G4DynamicParti    297                     std::vector<G4DynamicParticle*>& vd,
290                     const G4Track& track,         298                     const G4Track& track,
291                     G4VEmModel* currentModel,     299                     G4VEmModel* currentModel,
292                     G4ParticleChangeForLoss* p    300                     G4ParticleChangeForLoss* pPartChange,
293                     G4double& eloss,              301                     G4double& eloss,  
294                     G4int coupleIdx,              302                     G4int coupleIdx,
295                     G4double tcut,                303                     G4double tcut, 
296                     G4double safety)              304                     G4double safety)
297 {                                                 305 {
298   G4int index = idxSecBiasedCouple[coupleIdx];    306   G4int index = idxSecBiasedCouple[coupleIdx];
299   G4double weight = 1.;                        << 307   fWeight = 1.0;
300   if(0 <= index) {                                308   if(0 <= index) {
301     std::size_t n = vd.size();                 << 309     size_t n = vd.size();
302                                                   310 
303     // the check cannot be applied per seconda    311     // the check cannot be applied per secondary particle
304     // because weight correction is common, so    312     // because weight correction is common, so the first
305     // secondary is checked                       313     // secondary is checked
306     if((0 < n && vd[0]->GetKineticEnergy() < s    314     if((0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index])
307           || fDirectionalSplitting) {             315           || fDirectionalSplitting) {
308                                                   316 
309       G4int nsplit = nBremSplitting[index];       317       G4int nsplit = nBremSplitting[index];
310                                                   318 
311       // Range cut                                319       // Range cut
312       if(0 == nsplit) {                           320       if(0 == nsplit) { 
313         if(safety > fSafetyMin) { ApplyRangeCu    321         if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
314                                                   322 
315         // Russian Roulette                       323         // Russian Roulette
316       } else if(1 == nsplit) {                    324       } else if(1 == nsplit) { 
317         weight = ApplyRussianRoulette(vd, inde << 325         fWeight = ApplyRussianRoulette(vd, index);
318                                                   326 
319         // Splitting                              327         // Splitting
320       } else {                                    328       } else {
321         if (fDirectionalSplitting) {              329         if (fDirectionalSplitting) {
322           weight = ApplyDirectionalSplitting(v << 330           ApplyDirectionalSplitting(vd, track, currentModel, index, tcut);
                                                   >> 331           fWeight = 1.;
323         } else {                                  332         } else {
324           G4double tmpEnergy = pPartChange->Ge    333           G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
325           G4ThreeVector tmpMomDir = pPartChang    334           G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
326                                                   335 
327           weight = ApplySplitting(vd, track, c << 336           fWeight = ApplySplitting(vd, track, currentModel, index, tcut);
328                                                   337 
329           pPartChange->SetProposedKineticEnerg    338           pPartChange->SetProposedKineticEnergy(tmpEnergy);
330           pPartChange->ProposeMomentumDirectio    339           pPartChange->ProposeMomentumDirection(tmpMomDir);
331         }                                         340         }
332       }                                           341       }
333     }                                             342     }
334   }                                               343   }
335   return weight;                               << 344   return fWeight;
336 }                                                 345 }
337                                                   346 
338 //....oooOO0OOooo........oooOO0OOooo........oo    347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339                                                   348 
340 G4double                                          349 G4double 
341 G4EmBiasingManager::ApplySecondaryBiasing(        350 G4EmBiasingManager::ApplySecondaryBiasing(
342                   std::vector<G4DynamicParticl    351                   std::vector<G4DynamicParticle*>& vd,
343                   const G4Track& track,           352                   const G4Track& track,
344                   G4VEmModel* currentModel,       353                   G4VEmModel* currentModel, 
345                   G4ParticleChangeForGamma* pP    354                   G4ParticleChangeForGamma* pPartChange,
346                   G4double& eloss,                355                   G4double& eloss,  
347                   G4int coupleIdx,                356                   G4int coupleIdx,
348                   G4double tcut,                  357                   G4double tcut, 
349                   G4double safety)                358                   G4double safety)
350 {                                                 359 {
351   G4int index = idxSecBiasedCouple[coupleIdx];    360   G4int index = idxSecBiasedCouple[coupleIdx];
352   G4double weight = 1.;                        << 361   fWeight = 1.0;
353   if(0 <= index) {                                362   if(0 <= index) {
354     std::size_t n = vd.size();                 << 363     size_t n = vd.size();
355                                                   364 
356     // the check cannot be applied per seconda    365     // the check cannot be applied per secondary particle
357     // because weight correction is common, so    366     // because weight correction is common, so the first
358     // secondary is checked                       367     // secondary is checked
359     if((0 < n && vd[0]->GetKineticEnergy() < s    368     if((0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index])
360           || fDirectionalSplitting) {             369           || fDirectionalSplitting) {
361                                                   370 
362       G4int nsplit = nBremSplitting[index];       371       G4int nsplit = nBremSplitting[index];
363                                                   372 
364       // Range cut                                373       // Range cut
365       if(0 == nsplit) {                           374       if(0 == nsplit) { 
366         if(safety > fSafetyMin) { ApplyRangeCu    375         if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
367                                                   376 
368         // Russian Roulette                       377         // Russian Roulette
369       } else if(1 == nsplit) {                    378       } else if(1 == nsplit) { 
370         weight = ApplyRussianRoulette(vd, inde << 379         fWeight = ApplyRussianRoulette(vd, index);
371                                                   380 
372         // Splitting                              381         // Splitting
373       } else {                                    382       } else {
374         if (fDirectionalSplitting) {              383         if (fDirectionalSplitting) {
375           weight = ApplyDirectionalSplitting(v << 384           ApplyDirectionalSplitting(vd, track, currentModel,
376                                     index, tcu    385                                     index, tcut, pPartChange);
                                                   >> 386           fWeight = 1.;
377         } else {                                  387         } else {
378           G4double tmpEnergy = pPartChange->Ge    388           G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
379           G4ThreeVector tmpMomDir = pPartChang    389           G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
380                                                   390 
381           weight = ApplySplitting(vd, track, c << 391           fWeight = ApplySplitting(vd, track, currentModel, index, tcut);
382                                                   392 
383           pPartChange->SetProposedKineticEnerg    393           pPartChange->SetProposedKineticEnergy(tmpEnergy);
384           pPartChange->ProposeMomentumDirectio    394           pPartChange->ProposeMomentumDirection(tmpMomDir);
385         }                                         395         }
386       }                                           396       }
387     }                                             397     }
388   }                                               398   }
389   return weight;                               << 399   return fWeight;
390 }                                                 400 }
391                                                   401 
392 //....oooOO0OOooo........oooOO0OOooo........oo    402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
393                                                   403 
394 G4double                                          404 G4double 
395 G4EmBiasingManager::ApplySecondaryBiasing(std:    405 G4EmBiasingManager::ApplySecondaryBiasing(std::vector<G4Track*>& track,
396                                           G4in    406                                           G4int coupleIdx)
397 {                                                 407 {
398   G4int index = idxSecBiasedCouple[coupleIdx];    408   G4int index = idxSecBiasedCouple[coupleIdx];
399   G4double weight = 1.;                        << 409   fWeight = 1.0;
400   if(0 <= index) {                                410   if(0 <= index) {
401     std::size_t n = track.size();              << 411     size_t n = track.size();
402                                                   412 
403     // the check cannot be applied per seconda    413     // the check cannot be applied per secondary particle
404     // because weight correction is common, so    414     // because weight correction is common, so the first
405     // secondary is checked                       415     // secondary is checked
406     if(0 < n && track[0]->GetKineticEnergy() <    416     if(0 < n && track[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
407                                                   417 
408       G4int nsplit = nBremSplitting[index];       418       G4int nsplit = nBremSplitting[index];
409                                                   419 
410         // Russian Roulette only                  420         // Russian Roulette only
411       if(1 == nsplit) {                           421       if(1 == nsplit) { 
412         weight = secBiasedWeight[index];       << 422         fWeight = secBiasedWeight[index];
413         for(std::size_t k=0; k<n; ++k) {       << 423         for(size_t k=0; k<n; ++k) {
414           if(G4UniformRand()*weight > 1.0) {   << 424           if(G4UniformRand()*fWeight > 1.0) {
415             const G4Track* t = track[k];          425             const G4Track* t = track[k];
416             delete t;                             426             delete t;
417             track[k] = nullptr;                << 427             track[k] = 0;
418           }                                       428           }
419         }                                         429         }
420       }                                           430       }
421     }                                             431     }
422   }                                               432   }
423   return weight;                               << 433   return fWeight;
424 }                                                 434 }
425                                                   435 
426 //....oooOO0OOooo........oooOO0OOooo........oo    436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
427                                                   437 
428 void                                              438 void
429 G4EmBiasingManager::ApplyRangeCut(std::vector<    439 G4EmBiasingManager::ApplyRangeCut(std::vector<G4DynamicParticle*>& vd,
430                                   const G4Trac    440                                   const G4Track& track,
431                                   G4double& el    441                                   G4double& eloss, G4double safety)
432 {                                                 442 {
433   std::size_t n = vd.size();                   << 443   size_t n = vd.size();
434   if(!eIonisation) {                              444   if(!eIonisation) { 
435     eIonisation =                                 445     eIonisation = 
436       G4LossTableManager::Instance()->GetEnerg    446       G4LossTableManager::Instance()->GetEnergyLossProcess(theElectron);
437   }                                               447   }
438   if(eIonisation) {                               448   if(eIonisation) { 
439     for(std::size_t k=0; k<n; ++k) {           << 449     for(size_t k=0; k<n; ++k) {
440       const G4DynamicParticle* dp = vd[k];        450       const G4DynamicParticle* dp = vd[k];
441       if(dp->GetDefinition() == theElectron) {    451       if(dp->GetDefinition() == theElectron) {
442         G4double e = dp->GetKineticEnergy();      452         G4double e = dp->GetKineticEnergy();
443         if(eIonisation->GetRange(e, track.GetM << 453         if(eIonisation->GetRangeForLoss(e, track.GetMaterialCutsCouple()) 
                                                   >> 454            < safety) {
444           eloss += e;                             455           eloss += e;
445           delete dp;                              456           delete dp;
446           vd[k] = nullptr;                     << 457           vd[k] = 0;
447         }                                         458         }
448       }                                           459       }
449     }                                             460     }
450   }                                               461   }
451 }                                                 462 }
452                                                   463 
453 //....oooOO0OOooo........oooOO0OOooo........oo    464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
454                                                   465 
455 G4bool G4EmBiasingManager::CheckDirection(G4Th    466 G4bool G4EmBiasingManager::CheckDirection(G4ThreeVector pos,
456                                           G4Th    467                                           G4ThreeVector momdir) const
457 {                                                 468 {
458   G4ThreeVector delta = fDirectionalSplittingT    469   G4ThreeVector delta = fDirectionalSplittingTarget - pos;
459   G4double angle = momdir.angle(delta);           470   G4double angle = momdir.angle(delta);
460   G4double dist = delta.cross(momdir).mag();      471   G4double dist = delta.cross(momdir).mag();
461   if (dist <= fDirectionalSplittingRadius && a    472   if (dist <= fDirectionalSplittingRadius && angle < halfpi) {
462     return true;                                  473     return true;
463   }                                               474   }
464   return false;                                   475   return false;
465 }                                                 476 }
466                                                   477 
467 //....oooOO0OOooo........oooOO0OOooo........oo    478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
468                                                   479 
469 G4double                                          480 G4double
470 G4EmBiasingManager::ApplySplitting(std::vector    481 G4EmBiasingManager::ApplySplitting(std::vector<G4DynamicParticle*>& vd,
471                                    const G4Tra    482                                    const G4Track& track,
472                                    G4VEmModel*    483                                    G4VEmModel* currentModel, 
473                                    G4int index    484                                    G4int index,
474                                    G4double tc    485                                    G4double tcut)
475 {                                                 486 {
476   // method is applied only if 1 secondary cre    487   // method is applied only if 1 secondary created PostStep 
477   // in the case of many secondaries there is     488   // in the case of many secondaries there is a contradiction
478   G4double weight = 1.;                        << 489   fWeight = 1.0;
479   std::size_t n = vd.size();                   << 490   size_t n = vd.size();
480   G4double w = secBiasedWeight[index];            491   G4double w = secBiasedWeight[index];
481                                                   492 
482   if(1 != n || 1.0 <= w) { return weight; }    << 493   if(1 != n || 1.0 <= w) { return fWeight; }
483                                                   494 
484   G4double trackWeight = track.GetWeight();       495   G4double trackWeight = track.GetWeight();
485   const G4DynamicParticle* dynParticle = track    496   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
486                                                   497 
487   G4int nsplit = nBremSplitting[index];           498   G4int nsplit = nBremSplitting[index];
488                                                   499 
489   // double splitting is suppressed               500   // double splitting is suppressed 
490   if(1 < nsplit && trackWeight>w) {               501   if(1 < nsplit && trackWeight>w) {
491                                                   502 
492     weight = w;                                << 503     fWeight = w;
493     if(nsplit > (G4int)tmpSecondaries.size())     504     if(nsplit > (G4int)tmpSecondaries.size()) { 
494       tmpSecondaries.reserve(nsplit);             505       tmpSecondaries.reserve(nsplit);
495     }                                             506     }
496     const G4MaterialCutsCouple* couple = track    507     const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
497     // start from 1, because already one secon    508     // start from 1, because already one secondary created
498     for(G4int k=1; k<nsplit; ++k) {               509     for(G4int k=1; k<nsplit; ++k) {
499       tmpSecondaries.clear();                     510       tmpSecondaries.clear();
500       currentModel->SampleSecondaries(&tmpSeco    511       currentModel->SampleSecondaries(&tmpSecondaries, couple, dynParticle, 
501                                       tcut);      512                                       tcut);
502       for (std::size_t kk=0; kk<tmpSecondaries << 513       for (size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
503         vd.push_back(tmpSecondaries[kk]);         514         vd.push_back(tmpSecondaries[kk]);
504       }                                           515       }
505     }                                             516     }
506   }                                               517   }
507   return weight;                               << 518   return fWeight;
508 }                                                 519 }
509                                                   520 
510 //....oooOO0OOooo........oooOO0OOooo........oo    521 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
511                                                   522 
512 G4double                                       << 523 void
513 G4EmBiasingManager::ApplyDirectionalSplitting(    524 G4EmBiasingManager::ApplyDirectionalSplitting(
514                                    std::vector    525                                    std::vector<G4DynamicParticle*>& vd,
515                                    const G4Tra    526                                    const G4Track& track,
516                                    G4VEmModel*    527                                    G4VEmModel* currentModel,
517                                    G4int index    528                                    G4int index,
518                                    G4double tc    529                                    G4double tcut,
519                                    G4ParticleC    530                                    G4ParticleChangeForGamma* partChange)
520 {                                                 531 {
521   // primary is gamma. do splitting/RR as appr    532   // primary is gamma. do splitting/RR as appropriate
522   // method applied for any number of secondar    533   // method applied for any number of secondaries
523                                                   534 
524   G4double weight = 1.0;                          535   G4double weight = 1.0;
525   G4double w = secBiasedWeight[index];            536   G4double w = secBiasedWeight[index];
526                                                   537 
527   fDirectionalSplittingWeights.clear();           538   fDirectionalSplittingWeights.clear();
528   if(1.0 <= w) {                                  539   if(1.0 <= w) {
529     fDirectionalSplittingWeights.push_back(wei << 540     fDirectionalSplittingWeights.push_back(fWeight);
530     return weight;                             << 541     return;
531   }                                               542   }
532                                                   543 
533   G4double trackWeight = track.GetWeight();       544   G4double trackWeight = track.GetWeight();
534   G4int nsplit = nBremSplitting[index];           545   G4int nsplit = nBremSplitting[index];
535                                                   546 
536   // double splitting is suppressed               547   // double splitting is suppressed
537   if(1 < nsplit && trackWeight>w) {               548   if(1 < nsplit && trackWeight>w) {
538                                                   549 
539     weight = w;                                   550     weight = w;
540     const G4ThreeVector pos = track.GetPositio    551     const G4ThreeVector pos = track.GetPosition();
541                                                   552 
542     G4bool foundPrimaryParticle = false;          553     G4bool foundPrimaryParticle = false;
543     G4double primaryEnergy = 0.;                  554     G4double primaryEnergy = 0.;
544     G4ThreeVector primaryMomdir(0.,0.,0.);        555     G4ThreeVector primaryMomdir(0.,0.,0.);
545     G4double primaryWeight = trackWeight;         556     G4double primaryWeight = trackWeight;
546                                                   557 
547     tmpSecondaries = vd;                          558     tmpSecondaries = vd;
548     vd.clear();                                   559     vd.clear();
549     vd.reserve(nsplit);                           560     vd.reserve(nsplit);
550     for (G4int k=0; k<nsplit; ++k) {              561     for (G4int k=0; k<nsplit; ++k) {
551       if (k>0) {  // for k==0, SampleSecondari    562       if (k>0) {  // for k==0, SampleSecondaries has already been called
552         tmpSecondaries.clear();                   563         tmpSecondaries.clear();
553         // SampleSecondaries modifies primary     564         // SampleSecondaries modifies primary info stored in partChange
554         currentModel->SampleSecondaries(&tmpSe    565         currentModel->SampleSecondaries(&tmpSecondaries,
555                                         track.    566                                         track.GetMaterialCutsCouple(),
556                                         track.    567                                         track.GetDynamicParticle(), tcut);
557       }                                           568       }
558       for (std::size_t kk=0; kk<tmpSecondaries << 569       for (auto sec : tmpSecondaries) {
559         if (tmpSecondaries[kk]->GetParticleDef << 570         if (sec->GetParticleDefinition() == theGamma) {
560           if (CheckDirection(pos, tmpSecondari << 571           if (CheckDirection(pos, sec->GetMomentumDirection())) {
561             vd.push_back(tmpSecondaries[kk]);  << 572             vd.push_back(sec);
562             fDirectionalSplittingWeights.push_ << 573             fDirectionalSplittingWeights.push_back(weight);
563           } else if (G4UniformRand() < w) {       574           } else if (G4UniformRand() < w) {
564             vd.push_back(tmpSecondaries[kk]);  << 575             vd.push_back(sec);
565             fDirectionalSplittingWeights.push_ << 576             fDirectionalSplittingWeights.push_back(1.0);
566           } else {                             << 
567             delete tmpSecondaries[kk];         << 
568             tmpSecondaries[kk] = nullptr;      << 
569           }                                       577           }
570         } else if (k==0) { // keep charged 2ry << 578         } else if (k==0) { // not gamma
571           vd.push_back(tmpSecondaries[kk]);    << 579           vd.push_back(sec);
572           fDirectionalSplittingWeights.push_ba << 580           fDirectionalSplittingWeights.push_back(1.);
573         } else {                               << 
574           delete tmpSecondaries[kk];           << 
575           tmpSecondaries[kk] = nullptr;        << 
576         }                                         581         }
577       }                                           582       }
578                                                   583 
579       // primary                                  584       // primary
580       G4double en = partChange->GetProposedKin    585       G4double en = partChange->GetProposedKineticEnergy();
581       if (en>0.) { // don't add if kinetic ene    586       if (en>0.) { // don't add if kinetic energy = 0
582         G4ThreeVector momdir = partChange->Get    587         G4ThreeVector momdir = partChange->GetProposedMomentumDirection();
583         if (CheckDirection(pos,momdir)) {         588         if (CheckDirection(pos,momdir)) {
584           // keep only one primary; others are    589           // keep only one primary; others are secondaries
585           if (!foundPrimaryParticle) {            590           if (!foundPrimaryParticle) {
586             primaryEnergy = en;                   591             primaryEnergy = en;
587             primaryMomdir = momdir;               592             primaryMomdir = momdir;
588             foundPrimaryParticle = true;          593             foundPrimaryParticle = true;
589             primaryWeight = weight;               594             primaryWeight = weight;
590           } else {                                595           } else {
591             auto dp = new G4DynamicParticle(th << 596             G4DynamicParticle* dp = new G4DynamicParticle(theGamma,
592                           partChange->GetPropo << 597                                     partChange->GetProposedMomentumDirection(),
593                           partChange->GetPropo << 598                                     partChange->GetProposedKineticEnergy());
594             vd.push_back(dp);                     599             vd.push_back(dp);
595             fDirectionalSplittingWeights.push_ << 600             fDirectionalSplittingWeights.push_back(weight);
596           }                                       601           }
597         } else if (G4UniformRand()<w) { // not    602         } else if (G4UniformRand()<w) { // not going to target. play RR.
598           if (!foundPrimaryParticle) {            603           if (!foundPrimaryParticle) {
599             foundPrimaryParticle = true;          604             foundPrimaryParticle = true;
600             primaryEnergy = en;                   605             primaryEnergy = en;
601             primaryMomdir = momdir;               606             primaryMomdir = momdir;
602             primaryWeight = 1.;                   607             primaryWeight = 1.;
603           } else {                                608           } else {
604             auto dp = new G4DynamicParticle(th << 609             G4DynamicParticle* dp = new G4DynamicParticle(theGamma,
605                           partChange->GetPropo << 610                                     partChange->GetProposedMomentumDirection(),
606                           partChange->GetPropo << 611                                     partChange->GetProposedKineticEnergy());
607             vd.push_back(dp);                     612             vd.push_back(dp);
608             fDirectionalSplittingWeights.push_ << 613             fDirectionalSplittingWeights.push_back(1.0);
609           }                                       614           }
610         }                                         615         }
611       }                                           616       }
612     }  // end of loop over nsplit                 617     }  // end of loop over nsplit
613                                                   618 
614     partChange->ProposeWeight(primaryWeight);     619     partChange->ProposeWeight(primaryWeight);
615     partChange->SetProposedKineticEnergy(prima    620     partChange->SetProposedKineticEnergy(primaryEnergy);
616     partChange->ProposeMomentumDirection(prima    621     partChange->ProposeMomentumDirection(primaryMomdir);
617   } else {                                        622   } else {
618     for (std::size_t i = 0; i < vd.size(); ++i << 623     for (size_t i = 0; i < vd.size(); ++i) {
619       fDirectionalSplittingWeights.push_back(1 << 624       fDirectionalSplittingWeights.push_back(trackWeight);
620     }                                             625     }
621   }                                               626   }
622                                                   627 
623   return weight;                               << 628   return;
624 }                                                 629 }
625                                                   630 
626 //....oooOO0OOooo........oooOO0OOooo........oo    631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
627                                                   632 
628 G4double G4EmBiasingManager::GetWeight(G4int i    633 G4double G4EmBiasingManager::GetWeight(G4int i)
629 {                                                 634 {
630   // normally return 1. If a directionally spl << 
631   //  return 1./(splitting factor)             << 
632   if (fDirectionalSplittingWeights.size() >= (    635   if (fDirectionalSplittingWeights.size() >= (unsigned int)(i+1) ) {
633     G4double w = fDirectionalSplittingWeights[ << 636     return fDirectionalSplittingWeights[i];
634     fDirectionalSplittingWeights[i] = 1.; // e << 637   }
635     return w;                                  << 638   else {
636   } else {                                     << 639     return fWeight; 
637     return 1.;                                 << 
638   }                                               640   }
639 }                                                 641 }
640                                                   642 
641 G4double                                       << 643 void
642 G4EmBiasingManager::ApplyDirectionalSplitting(    644 G4EmBiasingManager::ApplyDirectionalSplitting(
643                                   std::vector<    645                                   std::vector<G4DynamicParticle*>& vd,
644                                   const G4Trac    646                                   const G4Track& track,
645                                   G4VEmModel*     647                                   G4VEmModel* currentModel,
646                                   G4int index,    648                                   G4int index,
647                                   G4double tcu    649                                   G4double tcut)
648 {                                                 650 {
649   // primary is not a gamma. Do nothing with p << 651   // Do nothing with primary
650                                                   652 
651   G4double weight = 1.0;                          653   G4double weight = 1.0;
652   G4double w = secBiasedWeight[index];            654   G4double w = secBiasedWeight[index];
653                                                   655 
654   fDirectionalSplittingWeights.clear();           656   fDirectionalSplittingWeights.clear();
655   if(1.0 <= w) {                                  657   if(1.0 <= w) {
656     fDirectionalSplittingWeights.push_back(wei    658     fDirectionalSplittingWeights.push_back(weight);
657     return weight;                             << 659     return;
658   }                                               660   }
659                                                   661 
660   G4double trackWeight = track.GetWeight();       662   G4double trackWeight = track.GetWeight();
661   G4int nsplit = nBremSplitting[index];           663   G4int nsplit = nBremSplitting[index];
662                                                   664 
663   // double splitting is suppressed               665   // double splitting is suppressed
664   if(1 < nsplit && trackWeight>w) {               666   if(1 < nsplit && trackWeight>w) {
665                                                   667 
666     weight = w;                                   668     weight = w;
667     const G4ThreeVector pos = track.GetPositio    669     const G4ThreeVector pos = track.GetPosition();
668                                                   670 
669     tmpSecondaries = vd;                          671     tmpSecondaries = vd;
670     vd.clear();                                   672     vd.clear();
671     vd.reserve(nsplit);                           673     vd.reserve(nsplit);
672     for (G4int k=0; k<nsplit; ++k) {              674     for (G4int k=0; k<nsplit; ++k) {
673       if (k>0) {                                  675       if (k>0) {
674         tmpSecondaries.clear();                   676         tmpSecondaries.clear();
675         currentModel->SampleSecondaries(&tmpSe    677         currentModel->SampleSecondaries(&tmpSecondaries,
676                                         track.    678                                         track.GetMaterialCutsCouple(),
677                                         track.    679                                         track.GetDynamicParticle(), tcut);
678       }                                           680       }
679       //for (auto sec : tmpSecondaries) {      << 681       for (auto sec : tmpSecondaries) {
680       for (std::size_t kk=0; kk < tmpSecondari << 682         if (CheckDirection(pos, sec->GetMomentumDirection())) {
681         if (CheckDirection(pos, tmpSecondaries << 683           vd.push_back(sec);
682           vd.push_back(tmpSecondaries[kk]);    << 684           fDirectionalSplittingWeights.push_back(weight);
683           fDirectionalSplittingWeights.push_ba << 
684         } else if (G4UniformRand()<w) {           685         } else if (G4UniformRand()<w) {
685           vd.push_back(tmpSecondaries[kk]);    << 686           vd.push_back(sec);
686           fDirectionalSplittingWeights.push_ba << 687           fDirectionalSplittingWeights.push_back(1.0);
687         } else {                               << 
688           delete tmpSecondaries[kk];           << 
689           tmpSecondaries[kk] = nullptr;        << 
690         }                                         688         }
691       }                                           689       }
692     }  // end of loop over nsplit                 690     }  // end of loop over nsplit
693   } else { // no splitting was done; still nee    691   } else { // no splitting was done; still need weights
694     for (std::size_t i = 0; i < vd.size(); ++i << 692     for (size_t i = 0; i < vd.size(); ++i) {
695       fDirectionalSplittingWeights.push_back(1 << 693       fDirectionalSplittingWeights.push_back(trackWeight);
696     }                                             694     }
697   }                                               695   }
698   return weight;                               << 696   return;
699 }                                                 697 }
700                                                   698