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