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