| 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 // Implementation of a custom tracking manager
27 // the same processes as defined in G4EmStanda
28 //
29 // Original author: Jonas Hahnfeld, 2021
30
31 #include "EmStandardPhysicsTrackingManager.hh"
32
33 #include "TrackingManagerHelper.hh"
34
35 #include "G4ComptonScattering.hh"
36 #include "G4CoulombScattering.hh"
37 #include "G4Electron.hh"
38 #include "G4EmParameters.hh"
39 #include "G4Gamma.hh"
40 #include "G4GammaConversion.hh"
41 #include "G4KleinNishinaModel.hh"
42 #include "G4LivermorePhotoElectricModel.hh"
43 #include "G4LivermorePolarizedRayleighModel.hh
44 #include "G4PhotoElectricAngularGeneratorPolar
45 #include "G4PhotoElectricEffect.hh"
46 #include "G4Positron.hh"
47 #include "G4RayleighScattering.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4UrbanMscModel.hh"
50 #include "G4WentzelVIModel.hh"
51 #include "G4eBremsstrahlung.hh"
52 #include "G4eCoulombScatteringModel.hh"
53 #include "G4eIonisation.hh"
54 #include "G4eMultipleScattering.hh"
55 #include "G4eplusAnnihilation.hh"
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 EmStandardPhysicsTrackingManager* EmStandardPh
60 nullptr;
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63
64 EmStandardPhysicsTrackingManager::EmStandardPh
65 {
66 G4EmParameters* param = G4EmParameters::Inst
67 G4double highEnergyLimit = param->MscEnergyL
68 G4bool polar = param->EnablePolarisation();
69
70 // e-
71 {
72 G4eMultipleScattering* msc = new G4eMultip
73 G4UrbanMscModel* msc1 = new G4UrbanMscMode
74 G4WentzelVIModel* msc2 = new G4WentzelVIMo
75 msc1->SetHighEnergyLimit(highEnergyLimit);
76 msc2->SetLowEnergyLimit(highEnergyLimit);
77 msc->SetEmModel(msc1);
78 msc->SetEmModel(msc2);
79 fElectronProcs.msc = msc;
80
81 fElectronProcs.ioni = new G4eIonisation;
82 fElectronProcs.brems = new G4eBremsstrahlu
83
84 G4CoulombScattering* ss = new G4CoulombSca
85 G4eCoulombScatteringModel* ssm = new G4eCo
86 ssm->SetLowEnergyLimit(highEnergyLimit);
87 ssm->SetActivationLowEnergyLimit(highEnerg
88 ss->SetEmModel(ssm);
89 ss->SetMinKinEnergy(highEnergyLimit);
90 fElectronProcs.ss = ss;
91 }
92
93 // e+
94 {
95 G4eMultipleScattering* msc = new G4eMultip
96 G4UrbanMscModel* msc1 = new G4UrbanMscMode
97 G4WentzelVIModel* msc2 = new G4WentzelVIMo
98 msc1->SetHighEnergyLimit(highEnergyLimit);
99 msc2->SetLowEnergyLimit(highEnergyLimit);
100 msc->SetEmModel(msc1);
101 msc->SetEmModel(msc2);
102 fPositronProcs.msc = msc;
103
104 fPositronProcs.ioni = new G4eIonisation;
105 fPositronProcs.brems = new G4eBremsstrahlu
106 fPositronProcs.annihilation = new G4eplusA
107
108 G4CoulombScattering* ss = new G4CoulombSca
109 G4eCoulombScatteringModel* ssm = new G4eCo
110 ssm->SetLowEnergyLimit(highEnergyLimit);
111 ssm->SetActivationLowEnergyLimit(highEnerg
112 ss->SetEmModel(ssm);
113 ss->SetMinKinEnergy(highEnergyLimit);
114 fPositronProcs.ss = ss;
115 }
116
117 {
118 G4PhotoElectricEffect* pe = new G4PhotoEle
119 G4VEmModel* peModel = new G4LivermorePhoto
120 if (polar) {
121 peModel->SetAngularDistribution(new G4Ph
122 }
123 pe->SetEmModel(peModel);
124 fGammaProcs.pe = pe;
125
126 G4ComptonScattering* cs = new G4ComptonSca
127 if (polar) {
128 cs->SetEmModel(new G4KleinNishinaModel);
129 }
130 fGammaProcs.compton = cs;
131
132 fGammaProcs.conversion = new G4GammaConver
133
134 G4RayleighScattering* rl = new G4RayleighS
135 if (polar) {
136 rl->SetEmModel(new G4LivermorePolarizedR
137 }
138 fGammaProcs.rayleigh = rl;
139 }
140
141 if (fMasterTrackingManager == nullptr) {
142 fMasterTrackingManager = this;
143 }
144 else {
145 fElectronProcs.msc->SetMasterProcess(fMast
146 fElectronProcs.ss->SetMasterProcess(fMaste
147 fElectronProcs.ioni->SetMasterProcess(fMas
148 fElectronProcs.brems->SetMasterProcess(fMa
149
150 fPositronProcs.msc->SetMasterProcess(fMast
151 fPositronProcs.ss->SetMasterProcess(fMaste
152 fPositronProcs.ioni->SetMasterProcess(fMas
153 fPositronProcs.brems->SetMasterProcess(fMa
154 fPositronProcs.annihilation->SetMasterProc
155 fMasterTrackingManager->fPositronProcs.a
156
157 fGammaProcs.pe->SetMasterProcess(fMasterTr
158 fGammaProcs.compton->SetMasterProcess(fMas
159 fGammaProcs.conversion->SetMasterProcess(f
160 fGammaProcs.rayleigh->SetMasterProcess(fMa
161 }
162 }
163
164 //....oooOO0OOooo........oooOO0OOooo........oo
165
166 EmStandardPhysicsTrackingManager::~EmStandardP
167 {
168 if (fMasterTrackingManager == this) {
169 fMasterTrackingManager = nullptr;
170 }
171 }
172
173 //....oooOO0OOooo........oooOO0OOooo........oo
174
175 void EmStandardPhysicsTrackingManager::BuildPh
176 {
177 if (&part == G4Electron::Definition()) {
178 fElectronProcs.msc->BuildPhysicsTable(part
179 fElectronProcs.ioni->BuildPhysicsTable(par
180 fElectronProcs.brems->BuildPhysicsTable(pa
181 fElectronProcs.ss->BuildPhysicsTable(part)
182 }
183 else if (&part == G4Positron::Definition())
184 fPositronProcs.msc->BuildPhysicsTable(part
185 fPositronProcs.ioni->BuildPhysicsTable(par
186 fPositronProcs.brems->BuildPhysicsTable(pa
187 fPositronProcs.annihilation->BuildPhysicsT
188 fPositronProcs.ss->BuildPhysicsTable(part)
189 }
190 else if (&part == G4Gamma::Definition()) {
191 fGammaProcs.pe->BuildPhysicsTable(part);
192 fGammaProcs.compton->BuildPhysicsTable(par
193 fGammaProcs.conversion->BuildPhysicsTable(
194 fGammaProcs.rayleigh->BuildPhysicsTable(pa
195 }
196 }
197
198 //....oooOO0OOooo........oooOO0OOooo........oo
199
200 void EmStandardPhysicsTrackingManager::Prepare
201 {
202 if (&part == G4Electron::Definition()) {
203 fElectronProcs.msc->PreparePhysicsTable(pa
204 fElectronProcs.ioni->PreparePhysicsTable(p
205 fElectronProcs.brems->PreparePhysicsTable(
206 fElectronProcs.ss->PreparePhysicsTable(par
207 }
208 else if (&part == G4Positron::Definition())
209 fPositronProcs.msc->PreparePhysicsTable(pa
210 fPositronProcs.ioni->PreparePhysicsTable(p
211 fPositronProcs.brems->PreparePhysicsTable(
212 fPositronProcs.annihilation->PreparePhysic
213 fPositronProcs.ss->PreparePhysicsTable(par
214 }
215 else if (&part == G4Gamma::Definition()) {
216 fGammaProcs.pe->PreparePhysicsTable(part);
217 fGammaProcs.compton->PreparePhysicsTable(p
218 fGammaProcs.conversion->PreparePhysicsTabl
219 fGammaProcs.rayleigh->PreparePhysicsTable(
220 }
221 }
222
223 //....oooOO0OOooo........oooOO0OOooo........oo
224
225 void EmStandardPhysicsTrackingManager::TrackEl
226 {
227 class ElectronPhysics final : public Trackin
228 {
229 public:
230 ElectronPhysics(EmStandardPhysicsTrackin
231
232 void StartTracking(G4Track* aTrack) over
233 {
234 auto& electronProcs = fMgr.fElectronPr
235
236 electronProcs.msc->StartTracking(aTrac
237 electronProcs.ioni->StartTracking(aTra
238 electronProcs.brems->StartTracking(aTr
239 electronProcs.ss->StartTracking(aTrack
240
241 fPreviousStepLength = 0;
242 }
243 void EndTracking() override
244 {
245 auto& electronProcs = fMgr.fElectronPr
246
247 electronProcs.msc->EndTracking();
248 electronProcs.ioni->EndTracking();
249 electronProcs.brems->EndTracking();
250 electronProcs.ss->EndTracking();
251 }
252
253 G4double GetPhysicalInteractionLength(co
254 {
255 auto& electronProcs = fMgr.fElectronPr
256 G4double physIntLength, proposedSafety
257 G4ForceCondition condition;
258 G4GPILSelection selection;
259
260 fProposedStep = DBL_MAX;
261 fSelected = -1;
262
263 physIntLength = electronProcs.ss->Post
264 if (physIntLength < fProposedStep) {
265 fProposedStep = physIntLength;
266 fSelected = 0;
267 }
268
269 physIntLength = electronProcs.brems->P
270 if (physIntLength < fProposedStep) {
271 fProposedStep = physIntLength;
272 fSelected = 1;
273 }
274
275 physIntLength = electronProcs.ioni->Po
276 if (physIntLength < fProposedStep) {
277 fProposedStep = physIntLength;
278 fSelected = 2;
279 }
280
281 physIntLength = electronProcs.ioni->Al
282
283 if (physIntLength < fProposedStep) {
284 fProposedStep = physIntLength;
285 fSelected = -1;
286 }
287
288 physIntLength = electronProcs.msc->Alo
289
290 if (physIntLength < fProposedStep) {
291 fProposedStep = physIntLength;
292 // Check if MSC actually wants to wi
293 // step size.
294 if (selection == CandidateForSelecti
295 fSelected = -1;
296 }
297 }
298
299 return fProposedStep;
300 }
301
302 void AlongStepDoIt(G4Track& track, G4Ste
303 {
304 if (step.GetStepLength() == fProposedS
305 step.GetPostStepPoint()->SetStepStat
306 }
307 else {
308 // Remember that the step was limite
309 fSelected = -1;
310 }
311 auto& electronProcs = fMgr.fElectronPr
312 G4VParticleChange* particleChange;
313
314 particleChange = electronProcs.msc->Al
315 particleChange->UpdateStepForAlongStep
316 track.SetTrackStatus(particleChange->G
317 particleChange->Clear();
318
319 particleChange = electronProcs.ioni->A
320 particleChange->UpdateStepForAlongStep
321 track.SetTrackStatus(particleChange->G
322 particleChange->Clear();
323
324 fPreviousStepLength = step.GetStepLeng
325 }
326
327 void PostStepDoIt(G4Track& track, G4Step
328 {
329 if (fSelected < 0) {
330 return;
331 }
332 step.GetPostStepPoint()->SetStepStatus
333
334 auto& electronProcs = fMgr.fElectronPr
335 G4VProcess* process = nullptr;
336 G4VParticleChange* particleChange = nu
337
338 switch (fSelected) {
339 case 0:
340 process = electronProcs.ss;
341 particleChange = electronProcs.ss-
342 break;
343 case 1:
344 process = electronProcs.brems;
345 particleChange = electronProcs.bre
346 break;
347 case 2:
348 process = electronProcs.ioni;
349 particleChange = electronProcs.ion
350 break;
351 }
352
353 particleChange->UpdateStepForPostStep(
354 step.UpdateTrack();
355
356 G4int numSecondaries = particleChange-
357 for (G4int i = 0; i < numSecondaries;
358 G4Track* secondary = particleChange-
359 secondary->SetParentID(track.GetTrac
360 secondary->SetCreatorProcess(process
361 secondaries.push_back(secondary);
362 }
363
364 track.SetTrackStatus(particleChange->G
365 particleChange->Clear();
366 }
367
368 private:
369 EmStandardPhysicsTrackingManager& fMgr;
370 G4double fPreviousStepLength;
371 G4double fProposedStep;
372 G4int fSelected;
373 };
374
375 ElectronPhysics physics(*this);
376 TrackingManagerHelper::TrackChargedParticle(
377 }
378
379 //....oooOO0OOooo........oooOO0OOooo........oo
380
381 void EmStandardPhysicsTrackingManager::TrackPo
382 {
383 class PositronPhysics final : public Trackin
384 {
385 public:
386 PositronPhysics(EmStandardPhysicsTrackin
387
388 void StartTracking(G4Track* aTrack) over
389 {
390 auto& positronProcs = fMgr.fPositronPr
391
392 positronProcs.msc->StartTracking(aTrac
393 positronProcs.ioni->StartTracking(aTra
394 positronProcs.brems->StartTracking(aTr
395 positronProcs.annihilation->StartTrack
396 positronProcs.ss->StartTracking(aTrack
397
398 fPreviousStepLength = 0;
399 }
400 void EndTracking() override
401 {
402 auto& positronProcs = fMgr.fPositronPr
403
404 positronProcs.msc->EndTracking();
405 positronProcs.ioni->EndTracking();
406 positronProcs.brems->EndTracking();
407 positronProcs.annihilation->EndTrackin
408 positronProcs.ss->EndTracking();
409 }
410
411 G4double GetPhysicalInteractionLength(co
412 {
413 auto& positronProcs = fMgr.fPositronPr
414 G4double physIntLength, proposedSafety
415 G4ForceCondition condition;
416 G4GPILSelection selection;
417
418 fProposedStep = DBL_MAX;
419 fSelected = -1;
420
421 physIntLength = positronProcs.ss->Post
422 if (physIntLength < fProposedStep) {
423 fProposedStep = physIntLength;
424 fSelected = 0;
425 }
426
427 physIntLength =
428 positronProcs.annihilation->PostStep
429 if (physIntLength < fProposedStep) {
430 fProposedStep = physIntLength;
431 fSelected = 1;
432 }
433
434 physIntLength = positronProcs.brems->P
435 if (physIntLength < fProposedStep) {
436 fProposedStep = physIntLength;
437 fSelected = 2;
438 }
439
440 physIntLength = positronProcs.ioni->Po
441 if (physIntLength < fProposedStep) {
442 fProposedStep = physIntLength;
443 fSelected = 3;
444 }
445
446 physIntLength = positronProcs.ioni->Al
447
448 if (physIntLength < fProposedStep) {
449 fProposedStep = physIntLength;
450 fSelected = -1;
451 }
452
453 physIntLength = positronProcs.msc->Alo
454
455 if (physIntLength < fProposedStep) {
456 fProposedStep = physIntLength;
457 // Check if MSC actually wants to wi
458 // step size.
459 if (selection == CandidateForSelecti
460 fSelected = -1;
461 }
462 }
463
464 return fProposedStep;
465 }
466
467 void AlongStepDoIt(G4Track& track, G4Ste
468 {
469 if (step.GetStepLength() == fProposedS
470 step.GetPostStepPoint()->SetStepStat
471 }
472 else {
473 // Remember that the step was limite
474 fSelected = -1;
475 }
476 auto& positronProcs = fMgr.fPositronPr
477 G4VParticleChange* particleChange;
478
479 particleChange = positronProcs.msc->Al
480 particleChange->UpdateStepForAlongStep
481 track.SetTrackStatus(particleChange->G
482 particleChange->Clear();
483
484 particleChange = positronProcs.ioni->A
485 particleChange->UpdateStepForAlongStep
486 track.SetTrackStatus(particleChange->G
487 particleChange->Clear();
488
489 fPreviousStepLength = step.GetStepLeng
490 }
491
492 void PostStepDoIt(G4Track& track, G4Step
493 {
494 if (fSelected < 0) {
495 return;
496 }
497 step.GetPostStepPoint()->SetStepStatus
498
499 auto& positronProcs = fMgr.fPositronPr
500 G4VProcess* process = nullptr;
501 G4VParticleChange* particleChange = nu
502
503 switch (fSelected) {
504 case 0:
505 process = positronProcs.ss;
506 particleChange = positronProcs.ss-
507 break;
508 case 1:
509 process = positronProcs.annihilati
510 particleChange = positronProcs.ann
511 break;
512 case 2:
513 process = positronProcs.brems;
514 particleChange = positronProcs.bre
515 break;
516 case 3:
517 process = positronProcs.ioni;
518 particleChange = positronProcs.ion
519 break;
520 }
521
522 particleChange->UpdateStepForPostStep(
523 step.UpdateTrack();
524
525 G4int numSecondaries = particleChange-
526 for (G4int i = 0; i < numSecondaries;
527 G4Track* secondary = particleChange-
528 secondary->SetParentID(track.GetTrac
529 secondary->SetCreatorProcess(process
530 secondaries.push_back(secondary);
531 }
532
533 track.SetTrackStatus(particleChange->G
534 particleChange->Clear();
535 }
536
537 G4bool HasAtRestProcesses() override { r
538
539 void AtRestDoIt(G4Track& track, G4Step&
540 {
541 auto& positronProcs = fMgr.fPositronPr
542 // Annihilate the positron at rest.
543 G4VParticleChange* particleChange = po
544 particleChange->UpdateStepForAtRest(&s
545 step.UpdateTrack();
546
547 G4int numSecondaries = particleChange-
548 for (G4int i = 0; i < numSecondaries;
549 G4Track* secondary = particleChange-
550 secondary->SetParentID(track.GetTrac
551 secondary->SetCreatorProcess(positro
552 secondaries.push_back(secondary);
553 }
554
555 track.SetTrackStatus(particleChange->G
556 particleChange->Clear();
557 }
558
559 private:
560 EmStandardPhysicsTrackingManager& fMgr;
561 G4double fPreviousStepLength;
562 G4double fProposedStep;
563 G4int fSelected;
564 };
565
566 PositronPhysics physics(*this);
567 TrackingManagerHelper::TrackChargedParticle(
568 }
569
570 //....oooOO0OOooo........oooOO0OOooo........oo
571
572 void EmStandardPhysicsTrackingManager::TrackGa
573 {
574 class GammaPhysics final : public TrackingMa
575 {
576 public:
577 GammaPhysics(EmStandardPhysicsTrackingMa
578
579 void StartTracking(G4Track* aTrack) over
580 {
581 auto& gammaProcs = fMgr.fGammaProcs;
582
583 gammaProcs.pe->StartTracking(aTrack);
584 gammaProcs.compton->StartTracking(aTra
585 gammaProcs.conversion->StartTracking(a
586 gammaProcs.rayleigh->StartTracking(aTr
587
588 fPreviousStepLength = 0;
589 }
590 void EndTracking() override
591 {
592 auto& gammaProcs = fMgr.fGammaProcs;
593
594 gammaProcs.pe->EndTracking();
595 gammaProcs.compton->EndTracking();
596 gammaProcs.conversion->EndTracking();
597 gammaProcs.rayleigh->EndTracking();
598 }
599
600 G4double GetPhysicalInteractionLength(co
601 {
602 auto& gammaProcs = fMgr.fGammaProcs;
603 G4double physIntLength;
604 G4ForceCondition condition;
605
606 fProposedStep = DBL_MAX;
607 fSelected = -1;
608
609 physIntLength = gammaProcs.rayleigh->P
610 if (physIntLength < fProposedStep) {
611 fProposedStep = physIntLength;
612 fSelected = 0;
613 }
614
615 physIntLength = gammaProcs.conversion-
616 if (physIntLength < fProposedStep) {
617 fProposedStep = physIntLength;
618 fSelected = 1;
619 }
620
621 physIntLength = gammaProcs.compton->Po
622 if (physIntLength < fProposedStep) {
623 fProposedStep = physIntLength;
624 fSelected = 2;
625 }
626
627 physIntLength = gammaProcs.pe->PostSte
628 if (physIntLength < fProposedStep) {
629 fProposedStep = physIntLength;
630 fSelected = 3;
631 }
632
633 return fProposedStep;
634 }
635
636 void AlongStepDoIt(G4Track&, G4Step& ste
637 {
638 if (step.GetStepLength() == fProposedS
639 step.GetPostStepPoint()->SetStepStat
640 }
641 else {
642 // Remember that the step was limite
643 fSelected = -1;
644 }
645 fPreviousStepLength = step.GetStepLeng
646 }
647
648 void PostStepDoIt(G4Track& track, G4Step
649 {
650 if (fSelected < 0) {
651 return;
652 }
653 step.GetPostStepPoint()->SetStepStatus
654
655 auto& gammaProcs = fMgr.fGammaProcs;
656 G4VProcess* process = nullptr;
657 G4VParticleChange* particleChange = nu
658
659 switch (fSelected) {
660 case 0:
661 process = gammaProcs.rayleigh;
662 particleChange = gammaProcs.raylei
663 break;
664 case 1:
665 process = gammaProcs.conversion;
666 particleChange = gammaProcs.conver
667 break;
668 case 2:
669 process = gammaProcs.compton;
670 particleChange = gammaProcs.compto
671 break;
672 case 3:
673 process = gammaProcs.pe;
674 particleChange = gammaProcs.pe->Po
675 break;
676 }
677
678 particleChange->UpdateStepForPostStep(
679 step.UpdateTrack();
680
681 G4int numSecondaries = particleChange-
682 for (G4int i = 0; i < numSecondaries;
683 G4Track* secondary = particleChange-
684 secondary->SetParentID(track.GetTrac
685 secondary->SetCreatorProcess(process
686 secondaries.push_back(secondary);
687 }
688
689 track.SetTrackStatus(particleChange->G
690 particleChange->Clear();
691 }
692
693 private:
694 EmStandardPhysicsTrackingManager& fMgr;
695 G4double fPreviousStepLength;
696 G4double fProposedStep;
697 G4int fSelected;
698 };
699
700 GammaPhysics physics(*this);
701 TrackingManagerHelper::TrackNeutralParticle(
702 }
703
704 //....oooOO0OOooo........oooOO0OOooo........oo
705
706 void EmStandardPhysicsTrackingManager::HandOve
707 {
708 const G4ParticleDefinition* part = aTrack->G
709
710 if (part == G4Electron::Definition()) {
711 TrackElectron(aTrack);
712 }
713 else if (part == G4Positron::Definition()) {
714 TrackPositron(aTrack);
715 }
716 else if (part == G4Gamma::Definition()) {
717 TrackGamma(aTrack);
718 }
719
720 aTrack->SetTrackStatus(fStopAndKill);
721 delete aTrack;
722 }
723
724 //....oooOO0OOooo........oooOO0OOooo........oo
725