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 // $Id: G4EmBiasingManager.cc 76765 2013-11-15 09:38:15Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class file 30 // GEANT4 Class file 30 // 31 // 31 // 32 // 32 // File name: G4EmBiasingManager 33 // File name: G4EmBiasingManager 33 // 34 // 34 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko 35 // 36 // 36 // Creation date: 28.07.2011 37 // Creation date: 28.07.2011 37 // 38 // 38 // Modifications: 39 // Modifications: 39 // 40 // 40 // 31-05-12 D. Sawkey put back in high energy 41 // 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 42 // 30-05-12 D. Sawkey brem split gammas are unique; do weight tests for 42 // brem, russian roulette 43 // brem, russian roulette 43 // ------------------------------------------- 44 // ------------------------------------------------------------------- 44 // 45 // 45 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 48 48 #include "G4EmBiasingManager.hh" 49 #include "G4EmBiasingManager.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.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" << 59 #include "G4VEmModel.hh" 58 #include "G4VEmModel.hh" 60 #include "G4LossTableManager.hh" 59 #include "G4LossTableManager.hh" 61 #include "G4ParticleChangeForLoss.hh" 60 #include "G4ParticleChangeForLoss.hh" 62 #include "G4ParticleChangeForGamma.hh" 61 #include "G4ParticleChangeForGamma.hh" 63 #include "G4EmParameters.hh" << 64 62 65 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 66 64 67 G4EmBiasingManager::G4EmBiasingManager() << 65 G4EmBiasingManager::G4EmBiasingManager() 68 : fDirectionalSplittingTarget(0.0,0.0,0.0) << 66 : nForcedRegions(0),nSecBiasedRegions(0),eIonisation(0), >> 67 currentStepLimit(0.0),startTracking(true) 69 { 68 { 70 fSafetyMin = 1.e-6*mm; 69 fSafetyMin = 1.e-6*mm; 71 theElectron = G4Electron::Electron(); 70 theElectron = G4Electron::Electron(); 72 theGamma = G4Gamma::Gamma(); << 73 } 71 } 74 72 75 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 76 74 77 G4EmBiasingManager::~G4EmBiasingManager() = de << 75 G4EmBiasingManager::~G4EmBiasingManager() >> 76 {} 78 77 79 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 80 79 81 void G4EmBiasingManager::Initialise(const G4Pa 80 void G4EmBiasingManager::Initialise(const G4ParticleDefinition& part, 82 const G4St << 81 const G4String& procName, G4int verbose) 83 { 82 { 84 //G4cout << "G4EmBiasingManager::Initialise 83 //G4cout << "G4EmBiasingManager::Initialise for " 85 // << part.GetParticleName() << 84 // << part.GetParticleName() 86 // << " and " << procName << G4endl; << 85 // << " and " << procName << G4endl; 87 const G4ProductionCutsTable* theCoupleTable= 86 const G4ProductionCutsTable* theCoupleTable= 88 G4ProductionCutsTable::GetProductionCutsTa 87 G4ProductionCutsTable::GetProductionCutsTable(); 89 G4int numOfCouples = (G4int)theCoupleTable-> << 88 size_t numOfCouples = theCoupleTable->GetTableSize(); 90 89 91 if(0 < nForcedRegions) { idxForcedCouple.res 90 if(0 < nForcedRegions) { idxForcedCouple.resize(numOfCouples, -1); } 92 if(0 < nSecBiasedRegions) { idxSecBiasedCoup 91 if(0 < nSecBiasedRegions) { idxSecBiasedCouple.resize(numOfCouples, -1); } 93 92 94 // Deexcitation 93 // Deexcitation 95 for (G4int j=0; j<numOfCouples; ++j) { << 94 for (size_t j=0; j<numOfCouples; ++j) { 96 const G4MaterialCutsCouple* couple = 95 const G4MaterialCutsCouple* couple = 97 theCoupleTable->GetMaterialCutsCouple(j) 96 theCoupleTable->GetMaterialCutsCouple(j); 98 const G4ProductionCuts* pcuts = couple->Ge 97 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); 99 if(0 < nForcedRegions) { 98 if(0 < nForcedRegions) { 100 for(G4int i=0; i<nForcedRegions; ++i) { 99 for(G4int i=0; i<nForcedRegions; ++i) { 101 if(forcedRegions[i]) { << 100 if(forcedRegions[i]) { 102 if(pcuts == forcedRegions[i]->GetPro << 101 if(pcuts == forcedRegions[i]->GetProductionCuts()) { 103 idxForcedCouple[j] = i; << 102 idxForcedCouple[j] = i; 104 break; << 103 break; 105 } << 104 } 106 } << 105 } 107 } 106 } 108 } 107 } 109 if(0 < nSecBiasedRegions) { 108 if(0 < nSecBiasedRegions) { 110 for(G4int i=0; i<nSecBiasedRegions; ++i) 109 for(G4int i=0; i<nSecBiasedRegions; ++i) { 111 if(secBiasedRegions[i]) { << 110 if(secBiasedRegions[i]) { 112 if(pcuts == secBiasedRegions[i]->Get << 111 if(pcuts == secBiasedRegions[i]->GetProductionCuts()) { 113 idxSecBiasedCouple[j] = i; << 112 idxSecBiasedCouple[j] = i; 114 break; << 113 break; 115 } << 114 } 116 } << 115 } 117 } 116 } 118 } 117 } 119 } 118 } 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) { 119 if (nForcedRegions > 0 && 0 < verbose) { 129 G4cout << " Forced Interaction is activate 120 G4cout << " Forced Interaction is activated for " 130 << part.GetParticleName() << " and << 121 << part.GetParticleName() << " and " 131 << procName << 122 << procName 132 << " inside G4Regions: " << G4endl; << 123 << " inside G4Regions: " << G4endl; 133 for (G4int i=0; i<nForcedRegions; ++i) { 124 for (G4int i=0; i<nForcedRegions; ++i) { 134 const G4Region* r = forcedRegions[i]; 125 const G4Region* r = forcedRegions[i]; 135 if(r) { G4cout << " " << r->Ge 126 if(r) { G4cout << " " << r->GetName() << G4endl; } 136 } 127 } 137 } 128 } 138 if (nSecBiasedRegions > 0 && 0 < verbose) { 129 if (nSecBiasedRegions > 0 && 0 < verbose) { 139 G4cout << " Secondary biasing is activated 130 G4cout << " Secondary biasing is activated for " 140 << part.GetParticleName() << " and << 131 << part.GetParticleName() << " and " 141 << procName << 132 << procName 142 << " inside G4Regions: " << G4endl; << 133 << " inside G4Regions: " << G4endl; 143 for (G4int i=0; i<nSecBiasedRegions; ++i) 134 for (G4int i=0; i<nSecBiasedRegions; ++i) { 144 const G4Region* r = secBiasedRegions[i]; 135 const G4Region* r = secBiasedRegions[i]; 145 if(r) { 136 if(r) { 146 G4cout << " " << r->GetName( << 137 G4cout << " " << r->GetName() 147 << " BiasingWeight= " << secBi << 138 << " BiasingWeight= " << secBiasedWeight[i] << G4endl; 148 } 139 } 149 } 140 } 150 if (fDirectionalSplitting) { << 151 G4cout << " Directional splitting ac << 152 << fDirectionalSplittingTarget/cm << 153 << " cm; radius: " << 154 << fDirectionalSplittingRadius/cm << 155 << "cm." << G4endl; << 156 } << 157 } 141 } 158 } 142 } 159 143 160 //....oooOO0OOooo........oooOO0OOooo........oo 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 161 145 162 void G4EmBiasingManager::ActivateForcedInterac 146 void G4EmBiasingManager::ActivateForcedInteraction(G4double val, 163 << 147 const G4String& rname) 164 { 148 { 165 G4RegionStore* regionStore = G4RegionStore:: 149 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 166 G4String name = rname; 150 G4String name = rname; 167 if(name == "" || name == "world" || name == 151 if(name == "" || name == "world" || name == "World") { 168 name = "DefaultRegionForTheWorld"; 152 name = "DefaultRegionForTheWorld"; 169 } 153 } 170 const G4Region* reg = regionStore->GetRegion 154 const G4Region* reg = regionStore->GetRegion(name, false); 171 if(!reg) { 155 if(!reg) { 172 G4cout << "### G4EmBiasingManager::ForcedI 156 G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: " 173 << " G4Region <" << 157 << " G4Region <" 174 << rname << "> is unknown" << G4end << 158 << rname << "> is unknown" << G4endl; 175 return; 159 return; 176 } 160 } 177 161 178 // the region is in the list 162 // the region is in the list 179 if (0 < nForcedRegions) { 163 if (0 < nForcedRegions) { 180 for (G4int i=0; i<nForcedRegions; ++i) { 164 for (G4int i=0; i<nForcedRegions; ++i) { 181 if (reg == forcedRegions[i]) { 165 if (reg == forcedRegions[i]) { 182 lengthForRegion[i] = val; << 166 lengthForRegion[i] = val; 183 return; 167 return; 184 } 168 } 185 } 169 } 186 } 170 } 187 if(val < 0.0) { 171 if(val < 0.0) { 188 G4cout << "### G4EmBiasingManager::ForcedI 172 G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: " 189 << val << " < 0.0, so no activation << 173 << val << " < 0.0, so no activation for the G4Region <" 190 << rname << ">" << G4endl; << 174 << rname << ">" << G4endl; 191 return; 175 return; 192 } 176 } 193 177 194 // new region 178 // new region 195 forcedRegions.push_back(reg); 179 forcedRegions.push_back(reg); 196 lengthForRegion.push_back(val); 180 lengthForRegion.push_back(val); 197 ++nForcedRegions; 181 ++nForcedRegions; 198 } 182 } 199 183 200 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 201 185 202 void 186 void 203 G4EmBiasingManager::ActivateSecondaryBiasing(c 187 G4EmBiasingManager::ActivateSecondaryBiasing(const G4String& rname, 204 G << 188 G4double factor, 205 G << 189 G4double energyLimit) 206 { 190 { 207 //G4cout << "G4EmBiasingManager::ActivateSec 191 //G4cout << "G4EmBiasingManager::ActivateSecondaryBiasing: " 208 // << rname << " F= " << factor << " << 192 // << rname << " F= " << factor << " E(MeV)= " << energyLimit/MeV 209 // << G4endl; << 193 // << G4endl; 210 G4RegionStore* regionStore = G4RegionStore:: 194 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 211 G4String name = rname; 195 G4String name = rname; 212 if(name == "" || name == "world" || name == 196 if(name == "" || name == "world" || name == "World") { 213 name = "DefaultRegionForTheWorld"; 197 name = "DefaultRegionForTheWorld"; 214 } 198 } 215 const G4Region* reg = regionStore->GetRegion 199 const G4Region* reg = regionStore->GetRegion(name, false); 216 if(!reg) { 200 if(!reg) { 217 G4cout << "### G4EmBiasingManager::Activat 201 G4cout << "### G4EmBiasingManager::ActivateBremsstrahlungSplitting " 218 << "WARNING: G4Region <" << 202 << "WARNING: G4Region <" 219 << rname << "> is unknown" << G4end << 203 << rname << "> is unknown" << G4endl; 220 return; 204 return; 221 } 205 } 222 206 223 // Range cut 207 // Range cut 224 G4int nsplit = 0; 208 G4int nsplit = 0; 225 G4double w = factor; 209 G4double w = factor; 226 210 227 // splitting 211 // splitting 228 if(factor >= 1.0) { 212 if(factor >= 1.0) { 229 nsplit = G4lrint(factor); 213 nsplit = G4lrint(factor); 230 w = 1.0/G4double(nsplit); 214 w = 1.0/G4double(nsplit); 231 215 232 // Russian roulette 216 // Russian roulette 233 } else if(0.0 < factor) { 217 } else if(0.0 < factor) { 234 nsplit = 1; 218 nsplit = 1; 235 w = 1.0/factor; 219 w = 1.0/factor; 236 } 220 } 237 221 238 // the region is in the list - overwrite par 222 // the region is in the list - overwrite parameters 239 if (0 < nSecBiasedRegions) { 223 if (0 < nSecBiasedRegions) { 240 for (G4int i=0; i<nSecBiasedRegions; ++i) 224 for (G4int i=0; i<nSecBiasedRegions; ++i) { 241 if (reg == secBiasedRegions[i]) { 225 if (reg == secBiasedRegions[i]) { 242 secBiasedWeight[i] = w; << 226 secBiasedWeight[i] = w; 243 nBremSplitting[i] = nsplit; 227 nBremSplitting[i] = nsplit; 244 secBiasedEnegryLimit[i] = energyLimit; << 228 secBiasedEnegryLimit[i] = energyLimit; 245 return; 229 return; 246 } 230 } 247 } 231 } 248 } 232 } 249 /* 233 /* 250 G4cout << "### G4EmBiasingManager::Activat 234 G4cout << "### G4EmBiasingManager::ActivateSecondaryBiasing: " 251 << " nsplit= " << nsplit << " for t << 235 << " nsplit= " << nsplit << " for the G4Region <" 252 << rname << ">" << G4endl; << 236 << rname << ">" << G4endl; 253 */ 237 */ 254 238 255 // new region 239 // new region 256 secBiasedRegions.push_back(reg); 240 secBiasedRegions.push_back(reg); 257 secBiasedWeight.push_back(w); 241 secBiasedWeight.push_back(w); 258 nBremSplitting.push_back(nsplit); 242 nBremSplitting.push_back(nsplit); 259 secBiasedEnegryLimit.push_back(energyLimit); 243 secBiasedEnegryLimit.push_back(energyLimit); 260 ++nSecBiasedRegions; 244 ++nSecBiasedRegions; 261 //G4cout << "nSecBiasedRegions= " << nSecBia 245 //G4cout << "nSecBiasedRegions= " << nSecBiasedRegions << G4endl; 262 } 246 } 263 247 264 //....oooOO0OOooo........oooOO0OOooo........oo 248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 265 249 266 G4double G4EmBiasingManager::GetStepLimit(G4in 250 G4double G4EmBiasingManager::GetStepLimit(G4int coupleIdx, 267 G4do << 251 G4double previousStep) 268 { 252 { 269 if(startTracking) { 253 if(startTracking) { 270 startTracking = false; 254 startTracking = false; 271 G4int i = idxForcedCouple[coupleIdx]; 255 G4int i = idxForcedCouple[coupleIdx]; 272 if(i < 0) { 256 if(i < 0) { 273 currentStepLimit = DBL_MAX; 257 currentStepLimit = DBL_MAX; 274 } else { 258 } else { 275 currentStepLimit = lengthForRegion[i]; 259 currentStepLimit = lengthForRegion[i]; 276 if(currentStepLimit > 0.0) { currentStep 260 if(currentStepLimit > 0.0) { currentStepLimit *= G4UniformRand(); } 277 } 261 } 278 } else { 262 } else { 279 currentStepLimit -= previousStep; 263 currentStepLimit -= previousStep; 280 } 264 } 281 if(currentStepLimit < 0.0) { currentStepLimi 265 if(currentStepLimit < 0.0) { currentStepLimit = 0.0; } 282 return currentStepLimit; 266 return currentStepLimit; 283 } 267 } 284 268 285 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 286 270 287 G4double 271 G4double 288 G4EmBiasingManager::ApplySecondaryBiasing( 272 G4EmBiasingManager::ApplySecondaryBiasing( 289 std::vector<G4DynamicParti << 273 std::vector<G4DynamicParticle*>& vd, 290 const G4Track& track, << 274 const G4Track& track, 291 G4VEmModel* currentModel, << 275 G4VEmModel* currentModel, 292 G4ParticleChangeForLoss* p << 276 G4ParticleChangeForLoss* pPartChange, 293 G4double& eloss, << 277 G4double& eloss, 294 G4int coupleIdx, << 278 G4int coupleIdx, 295 G4double tcut, << 279 G4double tcut, 296 G4double safety) << 280 G4double safety) 297 { 281 { 298 G4int index = idxSecBiasedCouple[coupleIdx]; 282 G4int index = idxSecBiasedCouple[coupleIdx]; 299 G4double weight = 1.; << 283 G4double weight = 1.0; 300 if(0 <= index) { 284 if(0 <= index) { 301 std::size_t n = vd.size(); << 285 size_t n = vd.size(); 302 286 303 // the check cannot be applied per seconda 287 // the check cannot be applied per secondary particle 304 // because weight correction is common, so 288 // because weight correction is common, so the first 305 // secondary is checked 289 // secondary is checked 306 if((0 < n && vd[0]->GetKineticEnergy() < s << 290 if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) { 307 || fDirectionalSplitting) { << 308 291 309 G4int nsplit = nBremSplitting[index]; 292 G4int nsplit = nBremSplitting[index]; 310 293 311 // Range cut 294 // Range cut 312 if(0 == nsplit) { 295 if(0 == nsplit) { 313 if(safety > fSafetyMin) { ApplyRangeCu << 296 if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); } 314 297 315 // Russian Roulette << 298 // Russian Roulette 316 } else if(1 == nsplit) { << 299 } if(1 == nsplit) { 317 weight = ApplyRussianRoulette(vd, inde << 300 weight = ApplyRussianRoulette(vd, index); 318 301 319 // Splitting << 302 // Splitting 320 } else { 303 } else { 321 if (fDirectionalSplitting) { << 304 G4double tmpEnergy = pPartChange->GetProposedKineticEnergy(); 322 weight = ApplyDirectionalSplitting(v << 305 G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection(); 323 } else { << 306 324 G4double tmpEnergy = pPartChange->Ge << 307 weight = ApplySplitting(vd, track, currentModel, index, tcut); 325 G4ThreeVector tmpMomDir = pPartChang << 308 326 << 309 pPartChange->SetProposedKineticEnergy(tmpEnergy); 327 weight = ApplySplitting(vd, track, c << 310 pPartChange->ProposeMomentumDirection(tmpMomDir); 328 << 329 pPartChange->SetProposedKineticEnerg << 330 pPartChange->ProposeMomentumDirectio << 331 } << 332 } 311 } 333 } 312 } 334 } 313 } 335 return weight; 314 return weight; 336 } 315 } 337 316 338 //....oooOO0OOooo........oooOO0OOooo........oo 317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 339 318 340 G4double 319 G4double 341 G4EmBiasingManager::ApplySecondaryBiasing( 320 G4EmBiasingManager::ApplySecondaryBiasing( 342 std::vector<G4DynamicParticl 321 std::vector<G4DynamicParticle*>& vd, 343 const G4Track& track, << 322 const G4Track& track, 344 G4VEmModel* currentModel, << 323 G4VEmModel* currentModel, 345 G4ParticleChangeForGamma* pP << 324 G4ParticleChangeForGamma* pPartChange, 346 G4double& eloss, << 325 G4double& eloss, 347 G4int coupleIdx, << 326 G4int coupleIdx, 348 G4double tcut, << 327 G4double tcut, 349 G4double safety) << 328 G4double safety) 350 { 329 { 351 G4int index = idxSecBiasedCouple[coupleIdx]; 330 G4int index = idxSecBiasedCouple[coupleIdx]; 352 G4double weight = 1.; << 331 G4double weight = 1.0; 353 if(0 <= index) { 332 if(0 <= index) { 354 std::size_t n = vd.size(); << 333 size_t n = vd.size(); 355 334 356 // the check cannot be applied per seconda 335 // the check cannot be applied per secondary particle 357 // because weight correction is common, so 336 // because weight correction is common, so the first 358 // secondary is checked 337 // secondary is checked 359 if((0 < n && vd[0]->GetKineticEnergy() < s << 338 if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) { 360 || fDirectionalSplitting) { << 361 339 362 G4int nsplit = nBremSplitting[index]; 340 G4int nsplit = nBremSplitting[index]; 363 341 364 // Range cut 342 // Range cut 365 if(0 == nsplit) { 343 if(0 == nsplit) { 366 if(safety > fSafetyMin) { ApplyRangeCu << 344 if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); } 367 345 368 // Russian Roulette << 346 // Russian Roulette 369 } else if(1 == nsplit) { 347 } else if(1 == nsplit) { 370 weight = ApplyRussianRoulette(vd, inde << 348 weight = ApplyRussianRoulette(vd, index); 371 349 372 // Splitting << 350 // Splitting 373 } else { 351 } else { 374 if (fDirectionalSplitting) { << 352 G4double tmpEnergy = pPartChange->GetProposedKineticEnergy(); 375 weight = ApplyDirectionalSplitting(v << 353 G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection(); 376 index, tcu << 354 377 } else { << 355 weight = ApplySplitting(vd, track, currentModel, index, tcut); 378 G4double tmpEnergy = pPartChange->Ge << 356 379 G4ThreeVector tmpMomDir = pPartChang << 357 pPartChange->SetProposedKineticEnergy(tmpEnergy); 380 << 358 pPartChange->ProposeMomentumDirection(tmpMomDir); 381 weight = ApplySplitting(vd, track, c << 382 << 383 pPartChange->SetProposedKineticEnerg << 384 pPartChange->ProposeMomentumDirectio << 385 } << 386 } 359 } 387 } 360 } 388 } 361 } 389 return weight; 362 return weight; 390 } 363 } 391 364 392 //....oooOO0OOooo........oooOO0OOooo........oo 365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 393 366 394 G4double 367 G4double 395 G4EmBiasingManager::ApplySecondaryBiasing(std: 368 G4EmBiasingManager::ApplySecondaryBiasing(std::vector<G4Track*>& track, 396 G4in << 369 G4int coupleIdx) 397 { 370 { 398 G4int index = idxSecBiasedCouple[coupleIdx]; 371 G4int index = idxSecBiasedCouple[coupleIdx]; 399 G4double weight = 1.; << 372 G4double weight = 1.0; 400 if(0 <= index) { 373 if(0 <= index) { 401 std::size_t n = track.size(); << 374 size_t n = track.size(); 402 375 403 // the check cannot be applied per seconda 376 // the check cannot be applied per secondary particle 404 // because weight correction is common, so 377 // because weight correction is common, so the first 405 // secondary is checked 378 // secondary is checked 406 if(0 < n && track[0]->GetKineticEnergy() < 379 if(0 < n && track[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) { 407 380 408 G4int nsplit = nBremSplitting[index]; 381 G4int nsplit = nBremSplitting[index]; 409 382 410 // Russian Roulette only << 383 // Russian Roulette only 411 if(1 == nsplit) { 384 if(1 == nsplit) { 412 weight = secBiasedWeight[index]; << 385 weight = secBiasedWeight[index]; 413 for(std::size_t k=0; k<n; ++k) { << 386 for(size_t k=0; k<n; ++k) { 414 if(G4UniformRand()*weight > 1.0) { << 387 if(G4UniformRand()*weight > 1.0) { 415 const G4Track* t = track[k]; << 388 const G4Track* t = track[k]; 416 delete t; << 389 delete t; 417 track[k] = nullptr; << 390 track[k] = 0; 418 } << 391 } 419 } << 392 } 420 } 393 } 421 } 394 } 422 } 395 } 423 return weight; 396 return weight; 424 } 397 } 425 398 426 //....oooOO0OOooo........oooOO0OOooo........oo 399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 427 400 428 void 401 void 429 G4EmBiasingManager::ApplyRangeCut(std::vector< 402 G4EmBiasingManager::ApplyRangeCut(std::vector<G4DynamicParticle*>& vd, 430 const G4Trac << 403 const G4Track& track, 431 G4double& el << 404 G4double& eloss, G4double safety) 432 { 405 { 433 std::size_t n = vd.size(); << 406 size_t n = vd.size(); 434 if(!eIonisation) { 407 if(!eIonisation) { 435 eIonisation = 408 eIonisation = 436 G4LossTableManager::Instance()->GetEnerg 409 G4LossTableManager::Instance()->GetEnergyLossProcess(theElectron); 437 } 410 } 438 if(eIonisation) { 411 if(eIonisation) { 439 for(std::size_t k=0; k<n; ++k) { << 412 for(size_t k=0; k<n; ++k) { 440 const G4DynamicParticle* dp = vd[k]; 413 const G4DynamicParticle* dp = vd[k]; 441 if(dp->GetDefinition() == theElectron) { 414 if(dp->GetDefinition() == theElectron) { 442 G4double e = dp->GetKineticEnergy(); << 415 G4double e = dp->GetKineticEnergy(); 443 if(eIonisation->GetRange(e, track.GetM << 416 if(eIonisation->GetRangeForLoss(e, track.GetMaterialCutsCouple()) 444 eloss += e; << 417 < safety) { 445 delete dp; << 418 eloss += e; 446 vd[k] = nullptr; << 419 delete dp; 447 } << 420 vd[k] = 0; >> 421 } 448 } 422 } 449 } 423 } 450 } 424 } 451 } 425 } 452 426 453 //....oooOO0OOooo........oooOO0OOooo........oo 427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 454 428 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 429 G4double 470 G4EmBiasingManager::ApplySplitting(std::vector 430 G4EmBiasingManager::ApplySplitting(std::vector<G4DynamicParticle*>& vd, 471 const G4Tra << 431 const G4Track& track, 472 G4VEmModel* << 432 G4VEmModel* currentModel, 473 G4int index << 433 G4int index, 474 G4double tc << 434 G4double tcut) 475 { 435 { 476 // method is applied only if 1 secondary cre 436 // method is applied only if 1 secondary created PostStep 477 // in the case of many secondaries there is << 437 // in the case of many secodndaries there is a contrudition 478 G4double weight = 1.; << 438 G4double weight = 1.0; 479 std::size_t n = vd.size(); << 439 size_t n = vd.size(); 480 G4double w = secBiasedWeight[index]; 440 G4double w = secBiasedWeight[index]; 481 441 482 if(1 != n || 1.0 <= w) { return weight; } 442 if(1 != n || 1.0 <= w) { return weight; } 483 443 484 G4double trackWeight = track.GetWeight(); 444 G4double trackWeight = track.GetWeight(); 485 const G4DynamicParticle* dynParticle = track 445 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); 486 446 487 G4int nsplit = nBremSplitting[index]; 447 G4int nsplit = nBremSplitting[index]; 488 448 489 // double splitting is suppressed << 449 // double splitting is supressed 490 if(1 < nsplit && trackWeight>w) { 450 if(1 < nsplit && trackWeight>w) { 491 451 492 weight = w; 452 weight = w; >> 453 // start from 1, because already one secondary created 493 if(nsplit > (G4int)tmpSecondaries.size()) 454 if(nsplit > (G4int)tmpSecondaries.size()) { 494 tmpSecondaries.reserve(nsplit); 455 tmpSecondaries.reserve(nsplit); 495 } 456 } 496 const G4MaterialCutsCouple* couple = track 457 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 497 // start from 1, because already one secon << 458 for(G4int k=1; k<nsplit; ++k) { 498 for(G4int k=1; k<nsplit; ++k) { << 499 tmpSecondaries.clear(); 459 tmpSecondaries.clear(); 500 currentModel->SampleSecondaries(&tmpSeco 460 currentModel->SampleSecondaries(&tmpSecondaries, couple, dynParticle, 501 tcut); << 461 tcut); 502 for (std::size_t kk=0; kk<tmpSecondaries << 462 for (size_t kk=0; kk<tmpSecondaries.size(); ++kk) { 503 vd.push_back(tmpSecondaries[kk]); << 463 vd.push_back(tmpSecondaries[kk]); 504 } 464 } 505 } 465 } 506 } 466 } 507 return weight; 467 return weight; 508 } 468 } 509 469 510 //....oooOO0OOooo........oooOO0OOooo........oo 470 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 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 471