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