Geant4_10
G4EmBiasingManager.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4EmBiasingManager.cc 76765 2013-11-15 09:38:15Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4EmBiasingManager
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 28.07.2011
38 //
39 // Modifications:
40 //
41 // 31-05-12 D. Sawkey put back in high energy limit for brem, russian roulette
42 // 30-05-12 D. Sawkey brem split gammas are unique; do weight tests for
43 // brem, russian roulette
44 // -------------------------------------------------------------------
45 //
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 #include "G4EmBiasingManager.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4MaterialCutsCouple.hh"
52 #include "G4ProductionCutsTable.hh"
53 #include "G4ProductionCuts.hh"
54 #include "G4Region.hh"
55 #include "G4RegionStore.hh"
56 #include "G4Track.hh"
57 #include "G4Electron.hh"
58 #include "G4VEmModel.hh"
59 #include "G4LossTableManager.hh"
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
66  : nForcedRegions(0),nSecBiasedRegions(0),eIonisation(0),
67  currentStepLimit(0.0),startTracking(true)
68 {
69  fSafetyMin = 1.e-6*mm;
70  theElectron = G4Electron::Electron();
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
76 {}
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 
81  const G4String& procName, G4int verbose)
82 {
83  //G4cout << "G4EmBiasingManager::Initialise for "
84  // << part.GetParticleName()
85  // << " and " << procName << G4endl;
86  const G4ProductionCutsTable* theCoupleTable=
88  size_t numOfCouples = theCoupleTable->GetTableSize();
89 
90  if(0 < nForcedRegions) { idxForcedCouple.resize(numOfCouples, -1); }
91  if(0 < nSecBiasedRegions) { idxSecBiasedCouple.resize(numOfCouples, -1); }
92 
93  // Deexcitation
94  for (size_t j=0; j<numOfCouples; ++j) {
95  const G4MaterialCutsCouple* couple =
96  theCoupleTable->GetMaterialCutsCouple(j);
97  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
98  if(0 < nForcedRegions) {
99  for(G4int i=0; i<nForcedRegions; ++i) {
100  if(forcedRegions[i]) {
101  if(pcuts == forcedRegions[i]->GetProductionCuts()) {
102  idxForcedCouple[j] = i;
103  break;
104  }
105  }
106  }
107  }
108  if(0 < nSecBiasedRegions) {
109  for(G4int i=0; i<nSecBiasedRegions; ++i) {
110  if(secBiasedRegions[i]) {
111  if(pcuts == secBiasedRegions[i]->GetProductionCuts()) {
112  idxSecBiasedCouple[j] = i;
113  break;
114  }
115  }
116  }
117  }
118  }
119  if (nForcedRegions > 0 && 0 < verbose) {
120  G4cout << " Forced Interaction is activated for "
121  << part.GetParticleName() << " and "
122  << procName
123  << " inside G4Regions: " << G4endl;
124  for (G4int i=0; i<nForcedRegions; ++i) {
125  const G4Region* r = forcedRegions[i];
126  if(r) { G4cout << " " << r->GetName() << G4endl; }
127  }
128  }
129  if (nSecBiasedRegions > 0 && 0 < verbose) {
130  G4cout << " Secondary biasing is activated for "
131  << part.GetParticleName() << " and "
132  << procName
133  << " inside G4Regions: " << G4endl;
134  for (G4int i=0; i<nSecBiasedRegions; ++i) {
135  const G4Region* r = secBiasedRegions[i];
136  if(r) {
137  G4cout << " " << r->GetName()
138  << " BiasingWeight= " << secBiasedWeight[i] << G4endl;
139  }
140  }
141  }
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
147  const G4String& rname)
148 {
149  G4RegionStore* regionStore = G4RegionStore::GetInstance();
150  G4String name = rname;
151  if(name == "" || name == "world" || name == "World") {
152  name = "DefaultRegionForTheWorld";
153  }
154  const G4Region* reg = regionStore->GetRegion(name, false);
155  if(!reg) {
156  G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
157  << " G4Region <"
158  << rname << "> is unknown" << G4endl;
159  return;
160  }
161 
162  // the region is in the list
163  if (0 < nForcedRegions) {
164  for (G4int i=0; i<nForcedRegions; ++i) {
165  if (reg == forcedRegions[i]) {
166  lengthForRegion[i] = val;
167  return;
168  }
169  }
170  }
171  if(val < 0.0) {
172  G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
173  << val << " < 0.0, so no activation for the G4Region <"
174  << rname << ">" << G4endl;
175  return;
176  }
177 
178  // new region
179  forcedRegions.push_back(reg);
180  lengthForRegion.push_back(val);
181  ++nForcedRegions;
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 
186 void
188  G4double factor,
189  G4double energyLimit)
190 {
191  //G4cout << "G4EmBiasingManager::ActivateSecondaryBiasing: "
192  // << rname << " F= " << factor << " E(MeV)= " << energyLimit/MeV
193  // << G4endl;
194  G4RegionStore* regionStore = G4RegionStore::GetInstance();
195  G4String name = rname;
196  if(name == "" || name == "world" || name == "World") {
197  name = "DefaultRegionForTheWorld";
198  }
199  const G4Region* reg = regionStore->GetRegion(name, false);
200  if(!reg) {
201  G4cout << "### G4EmBiasingManager::ActivateBremsstrahlungSplitting "
202  << "WARNING: G4Region <"
203  << rname << "> is unknown" << G4endl;
204  return;
205  }
206 
207  // Range cut
208  G4int nsplit = 0;
209  G4double w = factor;
210 
211  // splitting
212  if(factor >= 1.0) {
213  nsplit = G4lrint(factor);
214  w = 1.0/G4double(nsplit);
215 
216  // Russian roulette
217  } else if(0.0 < factor) {
218  nsplit = 1;
219  w = 1.0/factor;
220  }
221 
222  // the region is in the list - overwrite parameters
223  if (0 < nSecBiasedRegions) {
224  for (G4int i=0; i<nSecBiasedRegions; ++i) {
225  if (reg == secBiasedRegions[i]) {
226  secBiasedWeight[i] = w;
227  nBremSplitting[i] = nsplit;
228  secBiasedEnegryLimit[i] = energyLimit;
229  return;
230  }
231  }
232  }
233  /*
234  G4cout << "### G4EmBiasingManager::ActivateSecondaryBiasing: "
235  << " nsplit= " << nsplit << " for the G4Region <"
236  << rname << ">" << G4endl;
237  */
238 
239  // new region
240  secBiasedRegions.push_back(reg);
241  secBiasedWeight.push_back(w);
242  nBremSplitting.push_back(nsplit);
243  secBiasedEnegryLimit.push_back(energyLimit);
244  ++nSecBiasedRegions;
245  //G4cout << "nSecBiasedRegions= " << nSecBiasedRegions << G4endl;
246 }
247 
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249 
251  G4double previousStep)
252 {
253  if(startTracking) {
254  startTracking = false;
255  G4int i = idxForcedCouple[coupleIdx];
256  if(i < 0) {
257  currentStepLimit = DBL_MAX;
258  } else {
259  currentStepLimit = lengthForRegion[i];
260  if(currentStepLimit > 0.0) { currentStepLimit *= G4UniformRand(); }
261  }
262  } else {
263  currentStepLimit -= previousStep;
264  }
265  if(currentStepLimit < 0.0) { currentStepLimit = 0.0; }
266  return currentStepLimit;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270 
271 G4double
273  std::vector<G4DynamicParticle*>& vd,
274  const G4Track& track,
275  G4VEmModel* currentModel,
276  G4ParticleChangeForLoss* pPartChange,
277  G4double& eloss,
278  G4int coupleIdx,
279  G4double tcut,
280  G4double safety)
281 {
282  G4int index = idxSecBiasedCouple[coupleIdx];
283  G4double weight = 1.0;
284  if(0 <= index) {
285  size_t n = vd.size();
286 
287  // the check cannot be applied per secondary particle
288  // because weight correction is common, so the first
289  // secondary is checked
290  if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
291 
292  G4int nsplit = nBremSplitting[index];
293 
294  // Range cut
295  if(0 == nsplit) {
296  if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
297 
298  // Russian Roulette
299  } if(1 == nsplit) {
300  weight = ApplyRussianRoulette(vd, index);
301 
302  // Splitting
303  } else {
304  G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
305  G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
306 
307  weight = ApplySplitting(vd, track, currentModel, index, tcut);
308 
309  pPartChange->SetProposedKineticEnergy(tmpEnergy);
310  pPartChange->ProposeMomentumDirection(tmpMomDir);
311  }
312  }
313  }
314  return weight;
315 }
316 
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
318 
319 G4double
321  std::vector<G4DynamicParticle*>& vd,
322  const G4Track& track,
323  G4VEmModel* currentModel,
324  G4ParticleChangeForGamma* pPartChange,
325  G4double& eloss,
326  G4int coupleIdx,
327  G4double tcut,
328  G4double safety)
329 {
330  G4int index = idxSecBiasedCouple[coupleIdx];
331  G4double weight = 1.0;
332  if(0 <= index) {
333  size_t n = vd.size();
334 
335  // the check cannot be applied per secondary particle
336  // because weight correction is common, so the first
337  // secondary is checked
338  if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
339 
340  G4int nsplit = nBremSplitting[index];
341 
342  // Range cut
343  if(0 == nsplit) {
344  if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
345 
346  // Russian Roulette
347  } else if(1 == nsplit) {
348  weight = ApplyRussianRoulette(vd, index);
349 
350  // Splitting
351  } else {
352  G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
353  G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
354 
355  weight = ApplySplitting(vd, track, currentModel, index, tcut);
356 
357  pPartChange->SetProposedKineticEnergy(tmpEnergy);
358  pPartChange->ProposeMomentumDirection(tmpMomDir);
359  }
360  }
361  }
362  return weight;
363 }
364 
365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
366 
367 G4double
368 G4EmBiasingManager::ApplySecondaryBiasing(std::vector<G4Track*>& track,
369  G4int coupleIdx)
370 {
371  G4int index = idxSecBiasedCouple[coupleIdx];
372  G4double weight = 1.0;
373  if(0 <= index) {
374  size_t n = track.size();
375 
376  // the check cannot be applied per secondary particle
377  // because weight correction is common, so the first
378  // secondary is checked
379  if(0 < n && track[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
380 
381  G4int nsplit = nBremSplitting[index];
382 
383  // Russian Roulette only
384  if(1 == nsplit) {
385  weight = secBiasedWeight[index];
386  for(size_t k=0; k<n; ++k) {
387  if(G4UniformRand()*weight > 1.0) {
388  const G4Track* t = track[k];
389  delete t;
390  track[k] = 0;
391  }
392  }
393  }
394  }
395  }
396  return weight;
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 
401 void
402 G4EmBiasingManager::ApplyRangeCut(std::vector<G4DynamicParticle*>& vd,
403  const G4Track& track,
404  G4double& eloss, G4double safety)
405 {
406  size_t n = vd.size();
407  if(!eIonisation) {
408  eIonisation =
410  }
411  if(eIonisation) {
412  for(size_t k=0; k<n; ++k) {
413  const G4DynamicParticle* dp = vd[k];
414  if(dp->GetDefinition() == theElectron) {
415  G4double e = dp->GetKineticEnergy();
416  if(eIonisation->GetRangeForLoss(e, track.GetMaterialCutsCouple())
417  < safety) {
418  eloss += e;
419  delete dp;
420  vd[k] = 0;
421  }
422  }
423  }
424  }
425 }
426 
427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
428 
429 G4double
430 G4EmBiasingManager::ApplySplitting(std::vector<G4DynamicParticle*>& vd,
431  const G4Track& track,
432  G4VEmModel* currentModel,
433  G4int index,
434  G4double tcut)
435 {
436  // method is applied only if 1 secondary created PostStep
437  // in the case of many secodndaries there is a contrudition
438  G4double weight = 1.0;
439  size_t n = vd.size();
440  G4double w = secBiasedWeight[index];
441 
442  if(1 != n || 1.0 <= w) { return weight; }
443 
444  G4double trackWeight = track.GetWeight();
445  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
446 
447  G4int nsplit = nBremSplitting[index];
448 
449  // double splitting is supressed
450  if(1 < nsplit && trackWeight>w) {
451 
452  weight = w;
453  // start from 1, because already one secondary created
454  if(nsplit > (G4int)tmpSecondaries.size()) {
455  tmpSecondaries.reserve(nsplit);
456  }
457  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
458  for(G4int k=1; k<nsplit; ++k) {
459  tmpSecondaries.clear();
460  currentModel->SampleSecondaries(&tmpSecondaries, couple, dynParticle,
461  tcut);
462  for (size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
463  vd.push_back(tmpSecondaries[kk]);
464  }
465  }
466  }
467  return weight;
468 }
469 
470 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
const G4String & GetName() const
static G4LossTableManager * Instance()
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="")
G4double ApplySecondaryBiasing(std::vector< G4DynamicParticle * > &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForGamma *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
G4double GetKineticEnergy() const
Int_t index
Definition: macro.C:9
const G4DynamicParticle * GetDynamicParticle() const
G4double GetRangeForLoss(G4double &kineticEnergy, const G4MaterialCutsCouple *)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const XML_Char * name
Definition: expat.h:151
const G4ThreeVector & GetProposedMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
double weight
Definition: plottest35.C:25
G4double GetProposedKineticEnergy() const
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void Initialise(const G4ParticleDefinition &part, const G4String &procName, G4int verbose)
static G4RegionStore * GetInstance()
Char_t n[5]
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4double GetStepLimit(G4int coupleIdx, G4double previousStep)
TString part[npart]
Definition: Style.C:32
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ActivateSecondaryBiasing(const G4String &region, G4double factor, G4double energyLimit)
G4VEnergyLossProcess * GetEnergyLossProcess(const G4ParticleDefinition *)
static G4ProductionCutsTable * GetProductionCutsTable()
jump r
Definition: plot.C:36
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
G4double GetProposedKineticEnergy() const
G4double GetWeight() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
const G4ThreeVector & GetProposedMomentumDirection() const