| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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
41 // 30-05-12 D. Sawkey brem split gammas are u
42 // brem, russian roulette
43 // -------------------------------------------
44 //
45 //....oooOO0OOooo........oooOO0OOooo........oo
46 //....oooOO0OOooo........oooOO0OOooo........oo
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"
61 #include "G4ParticleChangeForLoss.hh"
62 #include "G4ParticleChangeForGamma.hh"
63 #include "G4EmParameters.hh"
64
65 //....oooOO0OOooo........oooOO0OOooo........oo
66
67 G4EmBiasingManager::G4EmBiasingManager()
68 : fDirectionalSplittingTarget(0.0,0.0,0.0)
69 {
70 fSafetyMin = 1.e-6*mm;
71 theElectron = G4Electron::Electron();
72 theGamma = G4Gamma::Gamma();
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oo
76
77 G4EmBiasingManager::~G4EmBiasingManager() = de
78
79 //....oooOO0OOooo........oooOO0OOooo........oo
80
81 void G4EmBiasingManager::Initialise(const G4Pa
82 const G4St
83 {
84 //G4cout << "G4EmBiasingManager::Initialise
85 // << part.GetParticleName()
86 // << " and " << procName << G4endl;
87 const G4ProductionCutsTable* theCoupleTable=
88 G4ProductionCutsTable::GetProductionCutsTa
89 G4int numOfCouples = (G4int)theCoupleTable->
90
91 if(0 < nForcedRegions) { idxForcedCouple.res
92 if(0 < nSecBiasedRegions) { idxSecBiasedCoup
93
94 // Deexcitation
95 for (G4int j=0; j<numOfCouples; ++j) {
96 const G4MaterialCutsCouple* couple =
97 theCoupleTable->GetMaterialCutsCouple(j)
98 const G4ProductionCuts* pcuts = couple->Ge
99 if(0 < nForcedRegions) {
100 for(G4int i=0; i<nForcedRegions; ++i) {
101 if(forcedRegions[i]) {
102 if(pcuts == forcedRegions[i]->GetPro
103 idxForcedCouple[j] = i;
104 break;
105 }
106 }
107 }
108 }
109 if(0 < nSecBiasedRegions) {
110 for(G4int i=0; i<nSecBiasedRegions; ++i)
111 if(secBiasedRegions[i]) {
112 if(pcuts == secBiasedRegions[i]->Get
113 idxSecBiasedCouple[j] = i;
114 break;
115 }
116 }
117 }
118 }
119 }
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) {
129 G4cout << " Forced Interaction is activate
130 << part.GetParticleName() << " and
131 << procName
132 << " inside G4Regions: " << G4endl;
133 for (G4int i=0; i<nForcedRegions; ++i) {
134 const G4Region* r = forcedRegions[i];
135 if(r) { G4cout << " " << r->Ge
136 }
137 }
138 if (nSecBiasedRegions > 0 && 0 < verbose) {
139 G4cout << " Secondary biasing is activated
140 << part.GetParticleName() << " and
141 << procName
142 << " inside G4Regions: " << G4endl;
143 for (G4int i=0; i<nSecBiasedRegions; ++i)
144 const G4Region* r = secBiasedRegions[i];
145 if(r) {
146 G4cout << " " << r->GetName(
147 << " BiasingWeight= " << secBi
148 }
149 }
150 if (fDirectionalSplitting) {
151 G4cout << " Directional splitting ac
152 << fDirectionalSplittingTarget/cm
153 << " cm; radius: "
154 << fDirectionalSplittingRadius/cm
155 << "cm." << G4endl;
156 }
157 }
158 }
159
160 //....oooOO0OOooo........oooOO0OOooo........oo
161
162 void G4EmBiasingManager::ActivateForcedInterac
163
164 {
165 G4RegionStore* regionStore = G4RegionStore::
166 G4String name = rname;
167 if(name == "" || name == "world" || name ==
168 name = "DefaultRegionForTheWorld";
169 }
170 const G4Region* reg = regionStore->GetRegion
171 if(!reg) {
172 G4cout << "### G4EmBiasingManager::ForcedI
173 << " G4Region <"
174 << rname << "> is unknown" << G4end
175 return;
176 }
177
178 // the region is in the list
179 if (0 < nForcedRegions) {
180 for (G4int i=0; i<nForcedRegions; ++i) {
181 if (reg == forcedRegions[i]) {
182 lengthForRegion[i] = val;
183 return;
184 }
185 }
186 }
187 if(val < 0.0) {
188 G4cout << "### G4EmBiasingManager::ForcedI
189 << val << " < 0.0, so no activation
190 << rname << ">" << G4endl;
191 return;
192 }
193
194 // new region
195 forcedRegions.push_back(reg);
196 lengthForRegion.push_back(val);
197 ++nForcedRegions;
198 }
199
200 //....oooOO0OOooo........oooOO0OOooo........oo
201
202 void
203 G4EmBiasingManager::ActivateSecondaryBiasing(c
204 G
205 G
206 {
207 //G4cout << "G4EmBiasingManager::ActivateSec
208 // << rname << " F= " << factor << "
209 // << G4endl;
210 G4RegionStore* regionStore = G4RegionStore::
211 G4String name = rname;
212 if(name == "" || name == "world" || name ==
213 name = "DefaultRegionForTheWorld";
214 }
215 const G4Region* reg = regionStore->GetRegion
216 if(!reg) {
217 G4cout << "### G4EmBiasingManager::Activat
218 << "WARNING: G4Region <"
219 << rname << "> is unknown" << G4end
220 return;
221 }
222
223 // Range cut
224 G4int nsplit = 0;
225 G4double w = factor;
226
227 // splitting
228 if(factor >= 1.0) {
229 nsplit = G4lrint(factor);
230 w = 1.0/G4double(nsplit);
231
232 // Russian roulette
233 } else if(0.0 < factor) {
234 nsplit = 1;
235 w = 1.0/factor;
236 }
237
238 // the region is in the list - overwrite par
239 if (0 < nSecBiasedRegions) {
240 for (G4int i=0; i<nSecBiasedRegions; ++i)
241 if (reg == secBiasedRegions[i]) {
242 secBiasedWeight[i] = w;
243 nBremSplitting[i] = nsplit;
244 secBiasedEnegryLimit[i] = energyLimit;
245 return;
246 }
247 }
248 }
249 /*
250 G4cout << "### G4EmBiasingManager::Activat
251 << " nsplit= " << nsplit << " for t
252 << rname << ">" << G4endl;
253 */
254
255 // new region
256 secBiasedRegions.push_back(reg);
257 secBiasedWeight.push_back(w);
258 nBremSplitting.push_back(nsplit);
259 secBiasedEnegryLimit.push_back(energyLimit);
260 ++nSecBiasedRegions;
261 //G4cout << "nSecBiasedRegions= " << nSecBia
262 }
263
264 //....oooOO0OOooo........oooOO0OOooo........oo
265
266 G4double G4EmBiasingManager::GetStepLimit(G4in
267 G4do
268 {
269 if(startTracking) {
270 startTracking = false;
271 G4int i = idxForcedCouple[coupleIdx];
272 if(i < 0) {
273 currentStepLimit = DBL_MAX;
274 } else {
275 currentStepLimit = lengthForRegion[i];
276 if(currentStepLimit > 0.0) { currentStep
277 }
278 } else {
279 currentStepLimit -= previousStep;
280 }
281 if(currentStepLimit < 0.0) { currentStepLimi
282 return currentStepLimit;
283 }
284
285 //....oooOO0OOooo........oooOO0OOooo........oo
286
287 G4double
288 G4EmBiasingManager::ApplySecondaryBiasing(
289 std::vector<G4DynamicParti
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 {
398 G4int index = idxSecBiasedCouple[coupleIdx];
399 G4double weight = 1.;
400 if(0 <= index) {
401 std::size_t n = track.size();
402
403 // the check cannot be applied per seconda
404 // because weight correction is common, so
405 // secondary is checked
406 if(0 < n && track[0]->GetKineticEnergy() <
407
408 G4int nsplit = nBremSplitting[index];
409
410 // Russian Roulette only
411 if(1 == nsplit) {
412 weight = secBiasedWeight[index];
413 for(std::size_t k=0; k<n; ++k) {
414 if(G4UniformRand()*weight > 1.0) {
415 const G4Track* t = track[k];
416 delete t;
417 track[k] = nullptr;
418 }
419 }
420 }
421 }
422 }
423 return weight;
424 }
425
426 //....oooOO0OOooo........oooOO0OOooo........oo
427
428 void
429 G4EmBiasingManager::ApplyRangeCut(std::vector<
430 const G4Trac
431 G4double& el
432 {
433 std::size_t n = vd.size();
434 if(!eIonisation) {
435 eIonisation =
436 G4LossTableManager::Instance()->GetEnerg
437 }
438 if(eIonisation) {
439 for(std::size_t k=0; k<n; ++k) {
440 const G4DynamicParticle* dp = vd[k];
441 if(dp->GetDefinition() == theElectron) {
442 G4double e = dp->GetKineticEnergy();
443 if(eIonisation->GetRange(e, track.GetM
444 eloss += e;
445 delete dp;
446 vd[k] = nullptr;
447 }
448 }
449 }
450 }
451 }
452
453 //....oooOO0OOooo........oooOO0OOooo........oo
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