| 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 (size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
559 if (tmpSecondaries[kk]->GetParticleDef 567 if (tmpSecondaries[kk]->GetParticleDefinition() == theGamma) {
560 if (CheckDirection(pos, tmpSecondari 568 if (CheckDirection(pos, tmpSecondaries[kk]->GetMomentumDirection())){
561 vd.push_back(tmpSecondaries[kk]); 569 vd.push_back(tmpSecondaries[kk]);
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(tmpSecondaries[kk]);
565 fDirectionalSplittingWeights.push_ 573 fDirectionalSplittingWeights.push_back(1./weight);
566 } else { 574 } else {
567 delete tmpSecondaries[kk]; 575 delete tmpSecondaries[kk];
568 tmpSecondaries[kk] = nullptr; 576 tmpSecondaries[kk] = nullptr;
569 } 577 }
570 } else if (k==0) { // keep charged 2ry 578 } else if (k==0) { // keep charged 2ry from first splitting
571 vd.push_back(tmpSecondaries[kk]); 579 vd.push_back(tmpSecondaries[kk]);
572 fDirectionalSplittingWeights.push_ba 580 fDirectionalSplittingWeights.push_back(1./weight);
573 } else { 581 } else {
574 delete tmpSecondaries[kk]; 582 delete tmpSecondaries[kk];
575 tmpSecondaries[kk] = nullptr; 583 tmpSecondaries[kk] = nullptr;
576 } 584 }
577 } 585 }
578 586
579 // primary 587 // primary
580 G4double en = partChange->GetProposedKin 588 G4double en = partChange->GetProposedKineticEnergy();
581 if (en>0.) { // don't add if kinetic ene 589 if (en>0.) { // don't add if kinetic energy = 0
582 G4ThreeVector momdir = partChange->Get 590 G4ThreeVector momdir = partChange->GetProposedMomentumDirection();
583 if (CheckDirection(pos,momdir)) { 591 if (CheckDirection(pos,momdir)) {
584 // keep only one primary; others are 592 // keep only one primary; others are secondaries
585 if (!foundPrimaryParticle) { 593 if (!foundPrimaryParticle) {
586 primaryEnergy = en; 594 primaryEnergy = en;
587 primaryMomdir = momdir; 595 primaryMomdir = momdir;
588 foundPrimaryParticle = true; 596 foundPrimaryParticle = true;
589 primaryWeight = weight; 597 primaryWeight = weight;
590 } else { 598 } else {
591 auto dp = new G4DynamicParticle(th << 599 G4DynamicParticle* dp = new G4DynamicParticle(theGamma,
592 partChange->GetPropo << 600 partChange->GetProposedMomentumDirection(),
593 partChange->GetPropo << 601 partChange->GetProposedKineticEnergy());
594 vd.push_back(dp); 602 vd.push_back(dp);
595 fDirectionalSplittingWeights.push_ 603 fDirectionalSplittingWeights.push_back(1.);
596 } 604 }
597 } else if (G4UniformRand()<w) { // not 605 } else if (G4UniformRand()<w) { // not going to target. play RR.
598 if (!foundPrimaryParticle) { 606 if (!foundPrimaryParticle) {
599 foundPrimaryParticle = true; 607 foundPrimaryParticle = true;
600 primaryEnergy = en; 608 primaryEnergy = en;
601 primaryMomdir = momdir; 609 primaryMomdir = momdir;
602 primaryWeight = 1.; 610 primaryWeight = 1.;
603 } else { 611 } else {
604 auto dp = new G4DynamicParticle(th << 612 G4DynamicParticle* dp = new G4DynamicParticle(theGamma,
605 partChange->GetPropo << 613 partChange->GetProposedMomentumDirection(),
606 partChange->GetPropo << 614 partChange->GetProposedKineticEnergy());
607 vd.push_back(dp); 615 vd.push_back(dp);
608 fDirectionalSplittingWeights.push_ 616 fDirectionalSplittingWeights.push_back(1./weight);
609 } 617 }
610 } 618 }
611 } 619 }
612 } // end of loop over nsplit 620 } // end of loop over nsplit
613 621
614 partChange->ProposeWeight(primaryWeight); 622 partChange->ProposeWeight(primaryWeight);
615 partChange->SetProposedKineticEnergy(prima 623 partChange->SetProposedKineticEnergy(primaryEnergy);
616 partChange->ProposeMomentumDirection(prima 624 partChange->ProposeMomentumDirection(primaryMomdir);
617 } else { 625 } else {
618 for (std::size_t i = 0; i < vd.size(); ++i << 626 for (size_t i = 0; i < vd.size(); ++i) {
619 fDirectionalSplittingWeights.push_back(1 627 fDirectionalSplittingWeights.push_back(1.);
620 } 628 }
621 } 629 }
622 630
623 return weight; 631 return weight;
624 } 632 }
625 633
626 //....oooOO0OOooo........oooOO0OOooo........oo 634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
627 635
628 G4double G4EmBiasingManager::GetWeight(G4int i 636 G4double G4EmBiasingManager::GetWeight(G4int i)
629 { 637 {
630 // normally return 1. If a directionally spl 638 // normally return 1. If a directionally split particle survives RR,
631 // return 1./(splitting factor) 639 // return 1./(splitting factor)
632 if (fDirectionalSplittingWeights.size() >= ( 640 if (fDirectionalSplittingWeights.size() >= (unsigned int)(i+1) ) {
633 G4double w = fDirectionalSplittingWeights[ 641 G4double w = fDirectionalSplittingWeights[i];
634 fDirectionalSplittingWeights[i] = 1.; // e 642 fDirectionalSplittingWeights[i] = 1.; // ensure it's not used again
635 return w; 643 return w;
636 } else { 644 } else {
637 return 1.; 645 return 1.;
638 } 646 }
639 } 647 }
640 648
641 G4double 649 G4double
642 G4EmBiasingManager::ApplyDirectionalSplitting( 650 G4EmBiasingManager::ApplyDirectionalSplitting(
643 std::vector< 651 std::vector<G4DynamicParticle*>& vd,
644 const G4Trac 652 const G4Track& track,
645 G4VEmModel* 653 G4VEmModel* currentModel,
646 G4int index, 654 G4int index,
647 G4double tcu 655 G4double tcut)
648 { 656 {
649 // primary is not a gamma. Do nothing with p 657 // primary is not a gamma. Do nothing with primary
650 658
651 G4double weight = 1.0; 659 G4double weight = 1.0;
652 G4double w = secBiasedWeight[index]; 660 G4double w = secBiasedWeight[index];
653 661
654 fDirectionalSplittingWeights.clear(); 662 fDirectionalSplittingWeights.clear();
655 if(1.0 <= w) { 663 if(1.0 <= w) {
656 fDirectionalSplittingWeights.push_back(wei 664 fDirectionalSplittingWeights.push_back(weight);
657 return weight; 665 return weight;
658 } 666 }
659 667
660 G4double trackWeight = track.GetWeight(); 668 G4double trackWeight = track.GetWeight();
661 G4int nsplit = nBremSplitting[index]; 669 G4int nsplit = nBremSplitting[index];
662 670
663 // double splitting is suppressed 671 // double splitting is suppressed
664 if(1 < nsplit && trackWeight>w) { 672 if(1 < nsplit && trackWeight>w) {
665 673
666 weight = w; 674 weight = w;
667 const G4ThreeVector pos = track.GetPositio 675 const G4ThreeVector pos = track.GetPosition();
668 676
669 tmpSecondaries = vd; 677 tmpSecondaries = vd;
670 vd.clear(); 678 vd.clear();
671 vd.reserve(nsplit); 679 vd.reserve(nsplit);
672 for (G4int k=0; k<nsplit; ++k) { 680 for (G4int k=0; k<nsplit; ++k) {
673 if (k>0) { 681 if (k>0) {
674 tmpSecondaries.clear(); 682 tmpSecondaries.clear();
675 currentModel->SampleSecondaries(&tmpSe 683 currentModel->SampleSecondaries(&tmpSecondaries,
676 track. 684 track.GetMaterialCutsCouple(),
677 track. 685 track.GetDynamicParticle(), tcut);
678 } 686 }
679 //for (auto sec : tmpSecondaries) { 687 //for (auto sec : tmpSecondaries) {
680 for (std::size_t kk=0; kk < tmpSecondari << 688 for (size_t kk=0; kk < tmpSecondaries.size(); ++kk) {
681 if (CheckDirection(pos, tmpSecondaries 689 if (CheckDirection(pos, tmpSecondaries[kk]->GetMomentumDirection())) {
682 vd.push_back(tmpSecondaries[kk]); 690 vd.push_back(tmpSecondaries[kk]);
683 fDirectionalSplittingWeights.push_ba 691 fDirectionalSplittingWeights.push_back(1.);
684 } else if (G4UniformRand()<w) { 692 } else if (G4UniformRand()<w) {
685 vd.push_back(tmpSecondaries[kk]); 693 vd.push_back(tmpSecondaries[kk]);
686 fDirectionalSplittingWeights.push_ba 694 fDirectionalSplittingWeights.push_back(1./weight);
687 } else { 695 } else {
688 delete tmpSecondaries[kk]; 696 delete tmpSecondaries[kk];
689 tmpSecondaries[kk] = nullptr; 697 tmpSecondaries[kk] = nullptr;
690 } 698 }
691 } 699 }
692 } // end of loop over nsplit 700 } // end of loop over nsplit
693 } else { // no splitting was done; still nee 701 } else { // no splitting was done; still need weights
694 for (std::size_t i = 0; i < vd.size(); ++i << 702 for (size_t i = 0; i < vd.size(); ++i) {
695 fDirectionalSplittingWeights.push_back(1 703 fDirectionalSplittingWeights.push_back(1.0);
696 } 704 }
697 } 705 }
698 return weight; 706 return weight;
699 } 707 }
700 708