| 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