ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmBiasingManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EmBiasingManager.cc
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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4EmBiasingManager
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 28.07.2011
37 //
38 // Modifications:
39 //
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 unique; do weight tests for
42 // brem, russian roulette
43 // -------------------------------------------------------------------
44 //
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
48 #include "G4EmBiasingManager.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4PhysicalConstants.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 "G4Gamma.hh"
59 #include "G4VEmModel.hh"
60 #include "G4LossTableManager.hh"
63 #include "G4EmParameters.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68  : nForcedRegions(0),nSecBiasedRegions(0),eIonisation(nullptr),
69  currentStepLimit(0.0),startTracking(true)
70 {
71  fSafetyMin = 1.e-6*mm;
74 
75  fDirectionalSplitting = false;
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 
84 {}
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 
89  const G4String& procName, G4int verbose)
90 {
91  //G4cout << "G4EmBiasingManager::Initialise for "
92  // << part.GetParticleName()
93  // << " and " << procName << G4endl;
94  const G4ProductionCutsTable* theCoupleTable=
96  size_t numOfCouples = theCoupleTable->GetTableSize();
97 
98  if(0 < nForcedRegions) { idxForcedCouple.resize(numOfCouples, -1); }
99  if(0 < nSecBiasedRegions) { idxSecBiasedCouple.resize(numOfCouples, -1); }
100 
101  // Deexcitation
102  for (size_t j=0; j<numOfCouples; ++j) {
103  const G4MaterialCutsCouple* couple =
104  theCoupleTable->GetMaterialCutsCouple(j);
105  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
106  if(0 < nForcedRegions) {
107  for(G4int i=0; i<nForcedRegions; ++i) {
108  if(forcedRegions[i]) {
109  if(pcuts == forcedRegions[i]->GetProductionCuts()) {
110  idxForcedCouple[j] = i;
111  break;
112  }
113  }
114  }
115  }
116  if(0 < nSecBiasedRegions) {
117  for(G4int i=0; i<nSecBiasedRegions; ++i) {
118  if(secBiasedRegions[i]) {
119  if(pcuts == secBiasedRegions[i]->GetProductionCuts()) {
120  idxSecBiasedCouple[j] = i;
121  break;
122  }
123  }
124  }
125  }
126  }
127 
130  if (fDirectionalSplitting) {
133  }
134 
135  if (nForcedRegions > 0 && 0 < verbose) {
136  G4cout << " Forced Interaction is activated for "
137  << part.GetParticleName() << " and "
138  << procName
139  << " inside G4Regions: " << G4endl;
140  for (G4int i=0; i<nForcedRegions; ++i) {
141  const G4Region* r = forcedRegions[i];
142  if(r) { G4cout << " " << r->GetName() << G4endl; }
143  }
144  }
145  if (nSecBiasedRegions > 0 && 0 < verbose) {
146  G4cout << " Secondary biasing is activated for "
147  << part.GetParticleName() << " and "
148  << procName
149  << " inside G4Regions: " << G4endl;
150  for (G4int i=0; i<nSecBiasedRegions; ++i) {
151  const G4Region* r = secBiasedRegions[i];
152  if(r) {
153  G4cout << " " << r->GetName()
154  << " BiasingWeight= " << secBiasedWeight[i] << G4endl;
155  }
156  }
157  if (fDirectionalSplitting) {
158  G4cout << " Directional splitting activated, with target position: "
160  << " cm; radius: "
162  << "cm." << G4endl;
163  }
164  }
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168 
170  const G4String& rname)
171 {
172  G4RegionStore* regionStore = G4RegionStore::GetInstance();
173  G4String name = rname;
174  if(name == "" || name == "world" || name == "World") {
175  name = "DefaultRegionForTheWorld";
176  }
177  const G4Region* reg = regionStore->GetRegion(name, false);
178  if(!reg) {
179  G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
180  << " G4Region <"
181  << rname << "> is unknown" << G4endl;
182  return;
183  }
184 
185  // the region is in the list
186  if (0 < nForcedRegions) {
187  for (G4int i=0; i<nForcedRegions; ++i) {
188  if (reg == forcedRegions[i]) {
189  lengthForRegion[i] = val;
190  return;
191  }
192  }
193  }
194  if(val < 0.0) {
195  G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
196  << val << " < 0.0, so no activation for the G4Region <"
197  << rname << ">" << G4endl;
198  return;
199  }
200 
201  // new region
202  forcedRegions.push_back(reg);
203  lengthForRegion.push_back(val);
204  ++nForcedRegions;
205 }
206 
207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208 
209 void
211  G4double factor,
212  G4double energyLimit)
213 {
214  //G4cout << "G4EmBiasingManager::ActivateSecondaryBiasing: "
215  // << rname << " F= " << factor << " E(MeV)= " << energyLimit/MeV
216  // << G4endl;
217  G4RegionStore* regionStore = G4RegionStore::GetInstance();
218  G4String name = rname;
219  if(name == "" || name == "world" || name == "World") {
220  name = "DefaultRegionForTheWorld";
221  }
222  const G4Region* reg = regionStore->GetRegion(name, false);
223  if(!reg) {
224  G4cout << "### G4EmBiasingManager::ActivateBremsstrahlungSplitting "
225  << "WARNING: G4Region <"
226  << rname << "> is unknown" << G4endl;
227  return;
228  }
229 
230  // Range cut
231  G4int nsplit = 0;
232  G4double w = factor;
233 
234  // splitting
235  if(factor >= 1.0) {
236  nsplit = G4lrint(factor);
237  w = 1.0/G4double(nsplit);
238 
239  // Russian roulette
240  } else if(0.0 < factor) {
241  nsplit = 1;
242  w = 1.0/factor;
243  }
244 
245  // the region is in the list - overwrite parameters
246  if (0 < nSecBiasedRegions) {
247  for (G4int i=0; i<nSecBiasedRegions; ++i) {
248  if (reg == secBiasedRegions[i]) {
249  secBiasedWeight[i] = w;
250  nBremSplitting[i] = nsplit;
251  secBiasedEnegryLimit[i] = energyLimit;
252  return;
253  }
254  }
255  }
256  /*
257  G4cout << "### G4EmBiasingManager::ActivateSecondaryBiasing: "
258  << " nsplit= " << nsplit << " for the G4Region <"
259  << rname << ">" << G4endl;
260  */
261 
262  // new region
263  secBiasedRegions.push_back(reg);
264  secBiasedWeight.push_back(w);
265  nBremSplitting.push_back(nsplit);
266  secBiasedEnegryLimit.push_back(energyLimit);
268  //G4cout << "nSecBiasedRegions= " << nSecBiasedRegions << G4endl;
269 }
270 
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
272 
274  G4double previousStep)
275 {
276  if(startTracking) {
277  startTracking = false;
278  G4int i = idxForcedCouple[coupleIdx];
279  if(i < 0) {
281  } else {
284  }
285  } else {
286  currentStepLimit -= previousStep;
287  }
288  if(currentStepLimit < 0.0) { currentStepLimit = 0.0; }
289  return currentStepLimit;
290 }
291 
292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
293 
294 G4double
296  std::vector<G4DynamicParticle*>& vd,
297  const G4Track& track,
298  G4VEmModel* currentModel,
299  G4ParticleChangeForLoss* pPartChange,
300  G4double& eloss,
301  G4int coupleIdx,
302  G4double tcut,
303  G4double safety)
304 {
305  G4int index = idxSecBiasedCouple[coupleIdx];
306  G4double weight = 1.;
307  if(0 <= index) {
308  size_t n = vd.size();
309 
310  // the check cannot be applied per secondary particle
311  // because weight correction is common, so the first
312  // secondary is checked
313  if((0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index])
315 
316  G4int nsplit = nBremSplitting[index];
317 
318  // Range cut
319  if(0 == nsplit) {
320  if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
321 
322  // Russian Roulette
323  } else if(1 == nsplit) {
324  weight = ApplyRussianRoulette(vd, index);
325 
326  // Splitting
327  } else {
328  if (fDirectionalSplitting) {
329  weight = ApplyDirectionalSplitting(vd, track, currentModel, index, tcut);
330  } else {
331  G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
332  G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
333 
334  weight = ApplySplitting(vd, track, currentModel, index, tcut);
335 
336  pPartChange->SetProposedKineticEnergy(tmpEnergy);
337  pPartChange->ProposeMomentumDirection(tmpMomDir);
338  }
339  }
340  }
341  }
342  return weight;
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346 
347 G4double
349  std::vector<G4DynamicParticle*>& vd,
350  const G4Track& track,
351  G4VEmModel* currentModel,
352  G4ParticleChangeForGamma* pPartChange,
353  G4double& eloss,
354  G4int coupleIdx,
355  G4double tcut,
356  G4double safety)
357 {
358  G4int index = idxSecBiasedCouple[coupleIdx];
359  G4double weight = 1.;
360  if(0 <= index) {
361  size_t n = vd.size();
362 
363  // the check cannot be applied per secondary particle
364  // because weight correction is common, so the first
365  // secondary is checked
366  if((0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index])
368 
369  G4int nsplit = nBremSplitting[index];
370 
371  // Range cut
372  if(0 == nsplit) {
373  if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
374 
375  // Russian Roulette
376  } else if(1 == nsplit) {
377  weight = ApplyRussianRoulette(vd, index);
378 
379  // Splitting
380  } else {
381  if (fDirectionalSplitting) {
382  weight = ApplyDirectionalSplitting(vd, track, currentModel,
383  index, tcut, pPartChange);
384  } else {
385  G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
386  G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
387 
388  weight = ApplySplitting(vd, track, currentModel, index, tcut);
389 
390  pPartChange->SetProposedKineticEnergy(tmpEnergy);
391  pPartChange->ProposeMomentumDirection(tmpMomDir);
392  }
393  }
394  }
395  }
396  return weight;
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 
401 G4double
403  G4int coupleIdx)
404 {
405  G4int index = idxSecBiasedCouple[coupleIdx];
406  G4double weight = 1.;
407  if(0 <= index) {
408  size_t n = track.size();
409 
410  // the check cannot be applied per secondary particle
411  // because weight correction is common, so the first
412  // secondary is checked
413  if(0 < n && track[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
414 
415  G4int nsplit = nBremSplitting[index];
416 
417  // Russian Roulette only
418  if(1 == nsplit) {
419  weight = secBiasedWeight[index];
420  for(size_t k=0; k<n; ++k) {
421  if(G4UniformRand()*weight > 1.0) {
422  const G4Track* t = track[k];
423  delete t;
424  track[k] = 0;
425  }
426  }
427  }
428  }
429  }
430  return weight;
431 }
432 
433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
434 
435 void
436 G4EmBiasingManager::ApplyRangeCut(std::vector<G4DynamicParticle*>& vd,
437  const G4Track& track,
438  G4double& eloss, G4double safety)
439 {
440  size_t n = vd.size();
441  if(!eIonisation) {
442  eIonisation =
444  }
445  if(eIonisation) {
446  for(size_t k=0; k<n; ++k) {
447  const G4DynamicParticle* dp = vd[k];
448  if(dp->GetDefinition() == theElectron) {
449  G4double e = dp->GetKineticEnergy();
451  < safety) {
452  eloss += e;
453  delete dp;
454  vd[k] = 0;
455  }
456  }
457  }
458  }
459 }
460 
461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
462 
464  G4ThreeVector momdir) const
465 {
467  G4double angle = momdir.angle(delta);
468  G4double dist = delta.cross(momdir).mag();
469  if (dist <= fDirectionalSplittingRadius && angle < halfpi) {
470  return true;
471  }
472  return false;
473 }
474 
475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476 
477 G4double
478 G4EmBiasingManager::ApplySplitting(std::vector<G4DynamicParticle*>& vd,
479  const G4Track& track,
480  G4VEmModel* currentModel,
481  G4int index,
482  G4double tcut)
483 {
484  // method is applied only if 1 secondary created PostStep
485  // in the case of many secondaries there is a contradiction
486  G4double weight = 1.;
487  size_t n = vd.size();
488  G4double w = secBiasedWeight[index];
489 
490  if(1 != n || 1.0 <= w) { return weight; }
491 
492  G4double trackWeight = track.GetWeight();
493  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
494 
495  G4int nsplit = nBremSplitting[index];
496 
497  // double splitting is suppressed
498  if(1 < nsplit && trackWeight>w) {
499 
500  weight = w;
501  if(nsplit > (G4int)tmpSecondaries.size()) {
502  tmpSecondaries.reserve(nsplit);
503  }
504  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
505  // start from 1, because already one secondary created
506  for(G4int k=1; k<nsplit; ++k) {
507  tmpSecondaries.clear();
508  currentModel->SampleSecondaries(&tmpSecondaries, couple, dynParticle,
509  tcut);
510  for (size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
511  vd.push_back(tmpSecondaries[kk]);
512  }
513  }
514  }
515  return weight;
516 }
517 
518 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
519 
520 G4double
522  std::vector<G4DynamicParticle*>& vd,
523  const G4Track& track,
524  G4VEmModel* currentModel,
525  G4int index,
526  G4double tcut,
527  G4ParticleChangeForGamma* partChange)
528 {
529  // primary is gamma. do splitting/RR as appropriate
530  // method applied for any number of secondaries
531 
532  G4double weight = 1.0;
533  G4double w = secBiasedWeight[index];
534 
536  if(1.0 <= w) {
537  fDirectionalSplittingWeights.push_back(weight);
538  return weight;
539  }
540 
541  G4double trackWeight = track.GetWeight();
542  G4int nsplit = nBremSplitting[index];
543 
544  // double splitting is suppressed
545  if(1 < nsplit && trackWeight>w) {
546 
547  weight = w;
548  const G4ThreeVector pos = track.GetPosition();
549 
550  G4bool foundPrimaryParticle = false;
551  G4double primaryEnergy = 0.;
552  G4ThreeVector primaryMomdir(0.,0.,0.);
553  G4double primaryWeight = trackWeight;
554 
555  tmpSecondaries = vd;
556  vd.clear();
557  vd.reserve(nsplit);
558  for (G4int k=0; k<nsplit; ++k) {
559  if (k>0) { // for k==0, SampleSecondaries has already been called
560  tmpSecondaries.clear();
561  // SampleSecondaries modifies primary info stored in partChange
562  currentModel->SampleSecondaries(&tmpSecondaries,
563  track.GetMaterialCutsCouple(),
564  track.GetDynamicParticle(), tcut);
565  }
566  for (size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
567  if (tmpSecondaries[kk]->GetParticleDefinition() == theGamma) {
568  if (CheckDirection(pos, tmpSecondaries[kk]->GetMomentumDirection())){
569  vd.push_back(tmpSecondaries[kk]);
570  fDirectionalSplittingWeights.push_back(1.);
571  } else if (G4UniformRand() < w) {
572  vd.push_back(tmpSecondaries[kk]);
573  fDirectionalSplittingWeights.push_back(1./weight);
574  } else {
575  delete tmpSecondaries[kk];
576  tmpSecondaries[kk] = nullptr;
577  }
578  } else if (k==0) { // keep charged 2ry from first splitting
579  vd.push_back(tmpSecondaries[kk]);
580  fDirectionalSplittingWeights.push_back(1./weight);
581  } else {
582  delete tmpSecondaries[kk];
583  tmpSecondaries[kk] = nullptr;
584  }
585  }
586 
587  // primary
588  G4double en = partChange->GetProposedKineticEnergy();
589  if (en>0.) { // don't add if kinetic energy = 0
590  G4ThreeVector momdir = partChange->GetProposedMomentumDirection();
591  if (CheckDirection(pos,momdir)) {
592  // keep only one primary; others are secondaries
593  if (!foundPrimaryParticle) {
594  primaryEnergy = en;
595  primaryMomdir = momdir;
596  foundPrimaryParticle = true;
597  primaryWeight = weight;
598  } else {
600  partChange->GetProposedMomentumDirection(),
601  partChange->GetProposedKineticEnergy());
602  vd.push_back(dp);
603  fDirectionalSplittingWeights.push_back(1.);
604  }
605  } else if (G4UniformRand()<w) { // not going to target. play RR.
606  if (!foundPrimaryParticle) {
607  foundPrimaryParticle = true;
608  primaryEnergy = en;
609  primaryMomdir = momdir;
610  primaryWeight = 1.;
611  } else {
613  partChange->GetProposedMomentumDirection(),
614  partChange->GetProposedKineticEnergy());
615  vd.push_back(dp);
616  fDirectionalSplittingWeights.push_back(1./weight);
617  }
618  }
619  }
620  } // end of loop over nsplit
621 
622  partChange->ProposeWeight(primaryWeight);
623  partChange->SetProposedKineticEnergy(primaryEnergy);
624  partChange->ProposeMomentumDirection(primaryMomdir);
625  } else {
626  for (size_t i = 0; i < vd.size(); ++i) {
627  fDirectionalSplittingWeights.push_back(1.);
628  }
629  }
630 
631  return weight;
632 }
633 
634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
635 
637 {
638  // normally return 1. If a directionally split particle survives RR,
639  // return 1./(splitting factor)
640  if (fDirectionalSplittingWeights.size() >= (unsigned int)(i+1) ) {
642  fDirectionalSplittingWeights[i] = 1.; // ensure it's not used again
643  return w;
644  } else {
645  return 1.;
646  }
647 }
648 
649 G4double
651  std::vector<G4DynamicParticle*>& vd,
652  const G4Track& track,
653  G4VEmModel* currentModel,
654  G4int index,
655  G4double tcut)
656 {
657  // primary is not a gamma. Do nothing with primary
658 
659  G4double weight = 1.0;
660  G4double w = secBiasedWeight[index];
661 
663  if(1.0 <= w) {
664  fDirectionalSplittingWeights.push_back(weight);
665  return weight;
666  }
667 
668  G4double trackWeight = track.GetWeight();
669  G4int nsplit = nBremSplitting[index];
670 
671  // double splitting is suppressed
672  if(1 < nsplit && trackWeight>w) {
673 
674  weight = w;
675  const G4ThreeVector pos = track.GetPosition();
676 
677  tmpSecondaries = vd;
678  vd.clear();
679  vd.reserve(nsplit);
680  for (G4int k=0; k<nsplit; ++k) {
681  if (k>0) {
682  tmpSecondaries.clear();
683  currentModel->SampleSecondaries(&tmpSecondaries,
684  track.GetMaterialCutsCouple(),
685  track.GetDynamicParticle(), tcut);
686  }
687  //for (auto sec : tmpSecondaries) {
688  for (size_t kk=0; kk < tmpSecondaries.size(); ++kk) {
689  if (CheckDirection(pos, tmpSecondaries[kk]->GetMomentumDirection())) {
690  vd.push_back(tmpSecondaries[kk]);
691  fDirectionalSplittingWeights.push_back(1.);
692  } else if (G4UniformRand()<w) {
693  vd.push_back(tmpSecondaries[kk]);
694  fDirectionalSplittingWeights.push_back(1./weight);
695  } else {
696  delete tmpSecondaries[kk];
697  tmpSecondaries[kk] = nullptr;
698  }
699  }
700  } // end of loop over nsplit
701  } else { // no splitting was done; still need weights
702  for (size_t i = 0; i < vd.size(); ++i) {
703  fDirectionalSplittingWeights.push_back(1.0);
704  }
705  }
706  return weight;
707 }