PhysicsLists/src/CMSEmStandardPhysicsTrackingManager.cc

0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsTrackingManager.h"
0002 #include "TrackingManagerHelper.h"
0003
0004 #include "G4CoulombScattering.hh"
0005 #include "G4UrbanMscModel.hh"
0006 #include "G4WentzelVIModel.hh"
0007 #include "G4eBremsstrahlung.hh"
0008 #include "G4eCoulombScatteringModel.hh"
0009 #include "G4eIonisation.hh"
0010 #include "G4eMultipleScattering.hh"
0011 #include "G4eplusAnnihilation.hh"
0012
0013 #include "G4ElectroVDNuclearModel.hh"
0014 #include "G4ElectronNuclearProcess.hh"
0015 #include "G4PositronNuclearProcess.hh"
0016
0017 #include "G4ComptonScattering.hh"
0018 #include "G4GammaConversion.hh"
0019 #include "G4PhotoElectricEffect.hh"
0020
0021 #include "G4CascadeInterface.hh"
0022 #include "G4CrossSectionDataSetRegistry.hh"
0023 #include "G4ExcitedStringDecay.hh"
0024 #include "G4GammaNuclearXS.hh"
0025 #include "G4GammaParticipants.hh"
0026 #include "G4GeneratorPrecompoundInterface.hh"
0027 #include "G4HadronInelasticProcess.hh"
0028 #include "G4HadronicParameters.hh"
0029 #include "G4LowEGammaNuclearModel.hh"
0030 #include "G4PhotoNuclearCrossSection.hh"
0031 #include "G4QGSMFragmentation.hh"
0032 #include "G4QGSModel.hh"
0033 #include "G4TheoFSGenerator.hh"
0034
0035 #include "G4GammaGeneralProcess.hh"
0036 #include "G4LossTableManager.hh"
0037
0038 #include "G4EmParameters.hh"
0039 #include <CLHEP/Units/SystemOfUnits.h>
0040
0041 #include "G4Electron.hh"
0042 #include "G4Gamma.hh"
0043 #include "G4Positron.hh"
0044
0045 #include "G4RegionStore.hh"
0046 #include "G4Region.hh"
0047 #include <string>
0048
0049 CMSEmStandardPhysicsTrackingManager *CMSEmStandardPhysicsTrackingManager::masterTrackingManager = nullptr;
0050
0051 CMSEmStandardPhysicsTrackingManager::CMSEmStandardPhysicsTrackingManager(const edm::ParameterSet &p) {
0052   fRangeFactor = p.getParameter<double>("G4MscRangeFactor");
0053   fGeomFactor = p.getParameter<double>("G4MscGeomFactor");
0054   fSafetyFactor = p.getParameter<double>("G4MscSafetyFactor");
0055   fLambdaLimit = p.getParameter<double>("G4MscLambdaLimit") * CLHEP::mm;
0056   std::string msc = p.getParameter<std::string>("G4MscStepLimit");
0057   fStepLimitType = fUseSafety;
0058   if (msc == "UseSafetyPlus") {
0059     fStepLimitType = fUseSafetyPlus;
0060   } else if (msc == "Minimal") {
0061     fStepLimitType = fMinimal;
0062   }
0063
0064   G4EmParameters *param = G4EmParameters::Instance();
0065   G4double highEnergyLimit = param->MscEnergyLimit();
0066
0067   const G4Region *aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
0068   const G4Region *bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
0069
0070   // e-
0071   {
0072     G4eMultipleScattering *msc = new G4eMultipleScattering;
0073     G4UrbanMscModel *msc1 = new G4UrbanMscModel;
0074     G4WentzelVIModel *msc2 = new G4WentzelVIModel;
0075     msc1->SetHighEnergyLimit(highEnergyLimit);
0076     msc2->SetLowEnergyLimit(highEnergyLimit);
0077     msc->SetEmModel(msc1);
0078     msc->SetEmModel(msc2);
0079
0080     // e-/e+ msc for HCAL and HGCAL using the Urban model
0081     if (nullptr != aRegion || nullptr != bRegion) {
0082       G4UrbanMscModel *msc3 = new G4UrbanMscModel();
0083       msc3->SetHighEnergyLimit(highEnergyLimit);
0084       msc3->SetRangeFactor(fRangeFactor);
0085       msc3->SetGeomFactor(fGeomFactor);
0086       msc3->SetSafetyFactor(fSafetyFactor);
0087       msc3->SetLambdaLimit(fLambdaLimit);
0088       msc3->SetStepLimitType(fStepLimitType);
0089       msc3->SetLocked(true);
0090       if (nullptr != aRegion) {
0091         msc->AddEmModel(-1, msc3, aRegion);
0092       }
0093       if (nullptr != bRegion) {
0094         msc->AddEmModel(-1, msc3, bRegion);
0095       }
0096     }
0097
0098     electron.msc = msc;
0099
0100     electron.ioni = new G4eIonisation;
0101     electron.brems = new G4eBremsstrahlung;
0102
0103     G4CoulombScattering *ss = new G4CoulombScattering;
0104     G4eCoulombScatteringModel *ssm = new G4eCoulombScatteringModel;
0105     ssm->SetLowEnergyLimit(highEnergyLimit);
0106     ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0107     ss->SetEmModel(ssm);
0108     ss->SetMinKinEnergy(highEnergyLimit);
0109     electron.ss = ss;
0110   }
0111
0112   // e+
0113   {
0114     G4eMultipleScattering *msc = new G4eMultipleScattering;
0115     G4UrbanMscModel *msc1 = new G4UrbanMscModel;
0116     G4WentzelVIModel *msc2 = new G4WentzelVIModel;
0117     msc1->SetHighEnergyLimit(highEnergyLimit);
0118     msc2->SetLowEnergyLimit(highEnergyLimit);
0119     msc->SetEmModel(msc1);
0120     msc->SetEmModel(msc2);
0121
0122     // e-/e+ msc for HCAL and HGCAL using the Urban model
0123     if (nullptr != aRegion || nullptr != bRegion) {
0124       G4UrbanMscModel *msc3 = new G4UrbanMscModel();
0125       msc3->SetHighEnergyLimit(highEnergyLimit);
0126       msc3->SetRangeFactor(fRangeFactor);
0127       msc3->SetGeomFactor(fGeomFactor);
0128       msc3->SetSafetyFactor(fSafetyFactor);
0129       msc3->SetLambdaLimit(fLambdaLimit);
0130       msc3->SetStepLimitType(fStepLimitType);
0131       msc3->SetLocked(true);
0132       if (nullptr != aRegion) {
0133         msc->AddEmModel(-1, msc3, aRegion);
0134       }
0135       if (nullptr != bRegion) {
0136         msc->AddEmModel(-1, msc3, bRegion);
0137       }
0138     }
0139
0140     positron.msc = msc;
0141
0142     positron.ioni = new G4eIonisation;
0143     positron.brems = new G4eBremsstrahlung;
0144     positron.annihilation = new G4eplusAnnihilation;
0145
0146     G4CoulombScattering *ss = new G4CoulombScattering;
0147     G4eCoulombScatteringModel *ssm = new G4eCoulombScatteringModel;
0148     ssm->SetLowEnergyLimit(highEnergyLimit);
0149     ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0150     ss->SetEmModel(ssm);
0151     ss->SetMinKinEnergy(highEnergyLimit);
0152     positron.ss = ss;
0153   }
0154
0155   // gamma
0156   {
0157     gammaProc = new G4GammaGeneralProcess;
0158
0159     gammaProc->AddEmProcess(new G4PhotoElectricEffect);
0160     gammaProc->AddEmProcess(new G4ComptonScattering);
0161     gammaProc->AddEmProcess(new G4GammaConversion);
0162
0163     G4HadronInelasticProcess *nuc = new G4HadronInelasticProcess("photonNuclear", G4Gamma::Definition());
0164     auto xsreg = G4CrossSectionDataSetRegistry::Instance();
0165     G4VCrossSectionDataSet *xs = nullptr;
0166     bool useGammaNuclearXS = true;
0167     if (useGammaNuclearXS) {
0168       xs = xsreg->GetCrossSectionDataSet("GammaNuclearXS");
0169       if (nullptr == xs)
0170         xs = new G4GammaNuclearXS;
0171     } else {
0172       xs = xsreg->GetCrossSectionDataSet("PhotoNuclearXS");
0173       if (nullptr == xs)
0174         xs = new G4PhotoNuclearCrossSection;
0175     }
0176     nuc->AddDataSet(xs);
0177
0178     G4QGSModel<G4GammaParticipants> *theStringModel = new G4QGSModel<G4GammaParticipants>;
0179     G4QGSMFragmentation *theFrag = new G4QGSMFragmentation;
0180     G4ExcitedStringDecay *theStringDecay = new G4ExcitedStringDecay(theFrag);
0181     theStringModel->SetFragmentationModel(theStringDecay);
0182
0183     G4GeneratorPrecompoundInterface *theCascade = new G4GeneratorPrecompoundInterface;
0184
0185     G4TheoFSGenerator *theModel = new G4TheoFSGenerator;
0186     theModel->SetTransport(theCascade);
0187     theModel->SetHighEnergyGenerator(theStringModel);
0188
0189     G4HadronicParameters *param = G4HadronicParameters::Instance();
0190
0191     G4CascadeInterface *cascade = new G4CascadeInterface;
0192
0193     // added low-energy model LEND disabled
0194     G4double gnLowEnergyLimit = 200 * CLHEP::MeV;
0195     if (gnLowEnergyLimit > 0.0) {
0196       G4LowEGammaNuclearModel *lemod = new G4LowEGammaNuclearModel;
0197       lemod->SetMaxEnergy(gnLowEnergyLimit);
0198       nuc->RegisterMe(lemod);
0199       cascade->SetMinEnergy(gnLowEnergyLimit - CLHEP::MeV);
0200     }
0201     cascade->SetMaxEnergy(param->GetMaxEnergyTransitionFTF_Cascade());
0202     nuc->RegisterMe(cascade);
0203     theModel->SetMinEnergy(param->GetMinEnergyTransitionFTF_Cascade());
0204     theModel->SetMaxEnergy(param->GetMaxEnergy());
0205     nuc->RegisterMe(theModel);
0206
0207     gammaProc->AddHadProcess(nuc);
0208     G4LossTableManager::Instance()->SetGammaGeneralProcess(gammaProc);
0209   }
0210
0211   // Create lepto-nuclear processes last, they access cross section data from
0212   // the gamma nuclear process!
0213   G4ElectroVDNuclearModel *eModel = new G4ElectroVDNuclearModel;
0214
0215   {
0216     G4ElectronNuclearProcess *nuc = new G4ElectronNuclearProcess;
0217     nuc->RegisterMe(eModel);
0218     electron.nuc = nuc;
0219   }
0220   {
0221     G4PositronNuclearProcess *nuc = new G4PositronNuclearProcess;
0222     nuc->RegisterMe(eModel);
0223     positron.nuc = nuc;
0224   }
0225
0226   if (masterTrackingManager == nullptr) {
0227     masterTrackingManager = this;
0228   } else {
0229     electron.msc->SetMasterProcess(masterTrackingManager->electron.msc);
0230     electron.ss->SetMasterProcess(masterTrackingManager->electron.ss);
0231     electron.ioni->SetMasterProcess(masterTrackingManager->electron.ioni);
0232     electron.brems->SetMasterProcess(masterTrackingManager->electron.brems);
0233     electron.nuc->SetMasterProcess(masterTrackingManager->electron.nuc);
0234
0235     positron.msc->SetMasterProcess(masterTrackingManager->positron.msc);
0236     positron.ss->SetMasterProcess(masterTrackingManager->positron.ss);
0237     positron.ioni->SetMasterProcess(masterTrackingManager->positron.ioni);
0238     positron.brems->SetMasterProcess(masterTrackingManager->positron.brems);
0239     positron.annihilation->SetMasterProcess(masterTrackingManager->positron.annihilation);
0240     positron.nuc->SetMasterProcess(masterTrackingManager->positron.nuc);
0241
0242     gammaProc->SetMasterProcess(masterTrackingManager->gammaProc);
0243   }
0244 }
0245
0246 CMSEmStandardPhysicsTrackingManager::~CMSEmStandardPhysicsTrackingManager() {
0247   if (masterTrackingManager == this) {
0248     masterTrackingManager = nullptr;
0249   }
0250 }
0251
0252 void CMSEmStandardPhysicsTrackingManager::BuildPhysicsTable(const G4ParticleDefinition &part) {
0253   if (&part == G4Electron::Definition()) {
0254     electron.msc->BuildPhysicsTable(part);
0255     electron.ioni->BuildPhysicsTable(part);
0256     electron.brems->BuildPhysicsTable(part);
0257     electron.ss->BuildPhysicsTable(part);
0258     electron.nuc->BuildPhysicsTable(part);
0259   } else if (&part == G4Positron::Definition()) {
0260     positron.msc->BuildPhysicsTable(part);
0261     positron.ioni->BuildPhysicsTable(part);
0262     positron.brems->BuildPhysicsTable(part);
0263     positron.annihilation->BuildPhysicsTable(part);
0264     positron.ss->BuildPhysicsTable(part);
0265     positron.nuc->BuildPhysicsTable(part);
0266   } else if (&part == G4Gamma::Definition()) {
0267     gammaProc->BuildPhysicsTable(part);
0268   }
0269 }
0270
0271 void CMSEmStandardPhysicsTrackingManager::PreparePhysicsTable(const G4ParticleDefinition &part) {
0272   if (&part == G4Electron::Definition()) {
0273     electron.msc->PreparePhysicsTable(part);
0274     electron.ioni->PreparePhysicsTable(part);
0275     electron.brems->PreparePhysicsTable(part);
0276     electron.ss->PreparePhysicsTable(part);
0277     electron.nuc->PreparePhysicsTable(part);
0278   } else if (&part == G4Positron::Definition()) {
0279     positron.msc->PreparePhysicsTable(part);
0280     positron.ioni->PreparePhysicsTable(part);
0281     positron.brems->PreparePhysicsTable(part);
0282     positron.annihilation->PreparePhysicsTable(part);
0283     positron.ss->PreparePhysicsTable(part);
0284     positron.nuc->PreparePhysicsTable(part);
0285   } else if (&part == G4Gamma::Definition()) {
0286     gammaProc->PreparePhysicsTable(part);
0287   }
0288 }
0289
0290 void CMSEmStandardPhysicsTrackingManager::TrackElectron(G4Track *aTrack) {
0291   class ElectronPhysics final : public TrackingManagerHelper::Physics {
0292   public:
0293     ElectronPhysics(CMSEmStandardPhysicsTrackingManager &mgr) : fMgr(mgr) {}
0294
0295     void StartTracking(G4Track *aTrack) override {
0296       auto &electron = fMgr.electron;
0297
0298       electron.msc->StartTracking(aTrack);
0299       electron.ioni->StartTracking(aTrack);
0300       electron.brems->StartTracking(aTrack);
0301       electron.ss->StartTracking(aTrack);
0302       electron.nuc->StartTracking(aTrack);
0303
0304       fPreviousStepLength = 0;
0305     }
0306     void EndTracking() override {
0307       auto &electron = fMgr.electron;
0308
0309       electron.msc->EndTracking();
0310       electron.ioni->EndTracking();
0311       electron.brems->EndTracking();
0312       electron.ss->EndTracking();
0313       electron.nuc->EndTracking();
0314     }
0315
0316     G4double GetPhysicalInteractionLength(const G4Track &track) override {
0317       auto &electron = fMgr.electron;
0318       G4double physIntLength, proposedSafety = DBL_MAX;
0319       G4ForceCondition condition;
0320       G4GPILSelection selection;
0321
0322       fProposedStep = DBL_MAX;
0323       fSelected = -1;
0324
0325       physIntLength = electron.nuc->PostStepGPIL(track, fPreviousStepLength, &condition);
0326       if (physIntLength < fProposedStep) {
0327         fProposedStep = physIntLength;
0328         fSelected = 3;
0329       }
0330
0331       physIntLength = electron.ss->PostStepGPIL(track, fPreviousStepLength, &condition);
0332       if (physIntLength < fProposedStep) {
0333         fProposedStep = physIntLength;
0334         fSelected = 0;
0335       }
0336
0337       physIntLength = electron.brems->PostStepGPIL(track, fPreviousStepLength, &condition);
0338       if (physIntLength < fProposedStep) {
0339         fProposedStep = physIntLength;
0340         fSelected = 1;
0341       }
0342
0343       physIntLength = electron.ioni->PostStepGPIL(track, fPreviousStepLength, &condition);
0344       if (physIntLength < fProposedStep) {
0345         fProposedStep = physIntLength;
0346         fSelected = 2;
0347       }
0348
0349       physIntLength =
0350           electron.ioni->AlongStepGPIL(track, fPreviousStepLength, fProposedStep, proposedSafety, &selection);
0351       if (physIntLength < fProposedStep) {
0352         fProposedStep = physIntLength;
0353         fSelected = -1;
0354       }
0355
0356       physIntLength =
0357           electron.msc->AlongStepGPIL(track, fPreviousStepLength, fProposedStep, proposedSafety, &selection);
0358       if (physIntLength < fProposedStep) {
0359         fProposedStep = physIntLength;
0360         // Check if MSC actually wants to win, in most cases it only limits the
0361         // step size.
0362         if (selection == CandidateForSelection) {
0363           fSelected = -1;
0364         }
0365       }
0366
0367       return fProposedStep;
0368     }
0369
0370     void AlongStepDoIt(G4Track &track, G4Step &step, G4TrackVector &) override {
0371       if (step.GetStepLength() == fProposedStep) {
0372         step.GetPostStepPoint()->SetStepStatus(fAlongStepDoItProc);
0373       } else {
0374         // Remember that the step was limited by geometry.
0375         fSelected = -1;
0376       }
0377       auto &electron = fMgr.electron;
0378       G4VParticleChange *particleChange;
0379
0380       particleChange = electron.msc->AlongStepDoIt(track, step);
0381       particleChange->UpdateStepForAlongStep(&step);
0382       track.SetTrackStatus(particleChange->GetTrackStatus());
0383       particleChange->Clear();
0384
0385       particleChange = electron.ioni->AlongStepDoIt(track, step);
0386       particleChange->UpdateStepForAlongStep(&step);
0387       track.SetTrackStatus(particleChange->GetTrackStatus());
0388       particleChange->Clear();
0389
0390       fPreviousStepLength = step.GetStepLength();
0391     }
0392
0393     void PostStepDoIt(G4Track &track, G4Step &step, G4TrackVector &secondaries) override {
0394       if (fSelected < 0) {
0395         return;
0396       }
0397       step.GetPostStepPoint()->SetStepStatus(fPostStepDoItProc);
0398
0399       auto &electron = fMgr.electron;
0400       G4VProcess *process = nullptr;
0401       G4VParticleChange *particleChange = nullptr;
0402
0403       switch (fSelected) {
0404         case 3:
0405           process = electron.nuc;
0406           particleChange = electron.nuc->PostStepDoIt(track, step);
0407           break;
0408         case 0:
0409           process = electron.ss;
0410           particleChange = electron.ss->PostStepDoIt(track, step);
0411           break;
0412         case 1:
0413           process = electron.brems;
0414           particleChange = electron.brems->PostStepDoIt(track, step);
0415           break;
0416         case 2:
0417           process = electron.ioni;
0418           particleChange = electron.ioni->PostStepDoIt(track, step);
0419           break;
0420       }
0421
0422       particleChange->UpdateStepForPostStep(&step);
0423       step.UpdateTrack();
0424
0425       int numSecondaries = particleChange->GetNumberOfSecondaries();
0426       for (int i = 0; i < numSecondaries; i++) {
0427         G4Track *secondary = particleChange->GetSecondary(i);
0428         secondary->SetParentID(track.GetTrackID());
0429         secondary->SetCreatorProcess(process);
0430         secondaries.push_back(secondary);
0431       }
0432
0433       track.SetTrackStatus(particleChange->GetTrackStatus());
0434       particleChange->Clear();
0435     }
0436
0437   private:
0438     CMSEmStandardPhysicsTrackingManager &fMgr;
0439     G4double fPreviousStepLength;
0440     G4double fProposedStep;
0441     G4int fSelected;
0442   };
0443
0444   ElectronPhysics physics(*this);
0445   TrackingManagerHelper::TrackChargedParticle(aTrack, physics);
0446 }
0447
0448 void CMSEmStandardPhysicsTrackingManager::TrackPositron(G4Track *aTrack) {
0449   class PositronPhysics final : public TrackingManagerHelper::Physics {
0450   public:
0451     PositronPhysics(CMSEmStandardPhysicsTrackingManager &mgr) : fMgr(mgr) {}
0452
0453     void StartTracking(G4Track *aTrack) override {
0454       auto &positron = fMgr.positron;
0455
0456       positron.msc->StartTracking(aTrack);
0457       positron.ioni->StartTracking(aTrack);
0458       positron.brems->StartTracking(aTrack);
0459       positron.annihilation->StartTracking(aTrack);
0460       positron.ss->StartTracking(aTrack);
0461       positron.nuc->StartTracking(aTrack);
0462
0463       fPreviousStepLength = 0;
0464     }
0465     void EndTracking() override {
0466       auto &positron = fMgr.positron;
0467
0468       positron.msc->EndTracking();
0469       positron.ioni->EndTracking();
0470       positron.brems->EndTracking();
0471       positron.annihilation->EndTracking();
0472       positron.ss->EndTracking();
0473       positron.nuc->EndTracking();
0474     }
0475
0476     G4double GetPhysicalInteractionLength(const G4Track &track) override {
0477       auto &positron = fMgr.positron;
0478       G4double physIntLength, proposedSafety = DBL_MAX;
0479       G4ForceCondition condition;
0480       G4GPILSelection selection;
0481
0482       fProposedStep = DBL_MAX;
0483       fSelected = -1;
0484
0485       physIntLength = positron.nuc->PostStepGPIL(track, fPreviousStepLength, &condition);
0486       if (physIntLength < fProposedStep) {
0487         fProposedStep = physIntLength;
0488         fSelected = 4;
0489       }
0490
0491       physIntLength = positron.ss->PostStepGPIL(track, fPreviousStepLength, &condition);
0492       if (physIntLength < fProposedStep) {
0493         fProposedStep = physIntLength;
0494         fSelected = 0;
0495       }
0496
0497       physIntLength = positron.annihilation->PostStepGPIL(track, fPreviousStepLength, &condition);
0498       if (physIntLength < fProposedStep) {
0499         fProposedStep = physIntLength;
0500         fSelected = 1;
0501       }
0502
0503       physIntLength = positron.brems->PostStepGPIL(track, fPreviousStepLength, &condition);
0504       if (physIntLength < fProposedStep) {
0505         fProposedStep = physIntLength;
0506         fSelected = 2;
0507       }
0508
0509       physIntLength = positron.ioni->PostStepGPIL(track, fPreviousStepLength, &condition);
0510       if (physIntLength < fProposedStep) {
0511         fProposedStep = physIntLength;
0512         fSelected = 3;
0513       }
0514
0515       physIntLength =
0516           positron.ioni->AlongStepGPIL(track, fPreviousStepLength, fProposedStep, proposedSafety, &selection);
0517       if (physIntLength < fProposedStep) {
0518         fProposedStep = physIntLength;
0519         fSelected = -1;
0520       }
0521
0522       physIntLength =
0523           positron.msc->AlongStepGPIL(track, fPreviousStepLength, fProposedStep, proposedSafety, &selection);
0524       if (physIntLength < fProposedStep) {
0525         fProposedStep = physIntLength;
0526         // Check if MSC actually wants to win, in most cases it only limits the
0527         // step size.
0528         if (selection == CandidateForSelection) {
0529           fSelected = -1;
0530         }
0531       }
0532
0533       return fProposedStep;
0534     }
0535
0536     void AlongStepDoIt(G4Track &track, G4Step &step, G4TrackVector &) override {
0537       if (step.GetStepLength() == fProposedStep) {
0538         step.GetPostStepPoint()->SetStepStatus(fAlongStepDoItProc);
0539       } else {
0540         // Remember that the step was limited by geometry.
0541         fSelected = -1;
0542       }
0543       auto &positron = fMgr.positron;
0544       G4VParticleChange *particleChange;
0545
0546       particleChange = positron.msc->AlongStepDoIt(track, step);
0547       particleChange->UpdateStepForAlongStep(&step);
0548       track.SetTrackStatus(particleChange->GetTrackStatus());
0549       particleChange->Clear();
0550
0551       particleChange = positron.ioni->AlongStepDoIt(track, step);
0552       particleChange->UpdateStepForAlongStep(&step);
0553       track.SetTrackStatus(particleChange->GetTrackStatus());
0554       particleChange->Clear();
0555
0556       fPreviousStepLength = step.GetStepLength();
0557     }
0558
0559     void PostStepDoIt(G4Track &track, G4Step &step, G4TrackVector &secondaries) override {
0560       if (fSelected < 0) {
0561         return;
0562       }
0563       step.GetPostStepPoint()->SetStepStatus(fPostStepDoItProc);
0564
0565       auto &positron = fMgr.positron;
0566       G4VProcess *process;
0567       G4VParticleChange *particleChange = nullptr;
0568
0569       switch (fSelected) {
0570         case 4:
0571           process = positron.nuc;
0572           particleChange = positron.nuc->PostStepDoIt(track, step);
0573           break;
0574         case 0:
0575           process = positron.ss;
0576           particleChange = positron.ss->PostStepDoIt(track, step);
0577           break;
0578         case 1:
0579           process = positron.annihilation;
0580           particleChange = positron.annihilation->PostStepDoIt(track, step);
0581           break;
0582         case 2:
0583           process = positron.brems;
0584           particleChange = positron.brems->PostStepDoIt(track, step);
0585           break;
0586         case 3:
0587           process = positron.ioni;
0588           particleChange = positron.ioni->PostStepDoIt(track, step);
0589           break;
0590       }
0591
0592       particleChange->UpdateStepForPostStep(&step);
0593       step.UpdateTrack();
0594
0595       int numSecondaries = particleChange->GetNumberOfSecondaries();
0596       for (int i = 0; i < numSecondaries; i++) {
0597         G4Track *secondary = particleChange->GetSecondary(i);
0598         secondary->SetParentID(track.GetTrackID());
0599         secondary->SetCreatorProcess(process);
0600         secondaries.push_back(secondary);
0601       }
0602
0603       track.SetTrackStatus(particleChange->GetTrackStatus());
0604       particleChange->Clear();
0605     }
0606
0607     G4bool HasAtRestProcesses() override { return true; }
0608
0609     void AtRestDoIt(G4Track &track, G4Step &step, G4TrackVector &secondaries) override {
0610       auto &positron = fMgr.positron;
0611       // Annihilate the positron at rest.
0612       G4VParticleChange *particleChange = positron.annihilation->AtRestDoIt(track, step);
0613       particleChange->UpdateStepForAtRest(&step);
0614       step.UpdateTrack();
0615
0616       int numSecondaries = particleChange->GetNumberOfSecondaries();
0617       for (int i = 0; i < numSecondaries; i++) {
0618         G4Track *secondary = particleChange->GetSecondary(i);
0619         secondary->SetParentID(track.GetTrackID());
0620         secondary->SetCreatorProcess(positron.annihilation);
0621         secondaries.push_back(secondary);
0622       }
0623
0624       track.SetTrackStatus(particleChange->GetTrackStatus());
0625       particleChange->Clear();
0626     }
0627
0628   private:
0629     CMSEmStandardPhysicsTrackingManager &fMgr;
0630     G4double fPreviousStepLength;
0631     G4double fProposedStep;
0632     G4int fSelected;
0633   };
0634
0635   PositronPhysics physics(*this);
0636   TrackingManagerHelper::TrackChargedParticle(aTrack, physics);
0637 }
0638
0639 void CMSEmStandardPhysicsTrackingManager::TrackGamma(G4Track *aTrack) {
0640   class GammaPhysics final : public TrackingManagerHelper::Physics {
0641   public:
0642     GammaPhysics(CMSEmStandardPhysicsTrackingManager &mgr) : fMgr(mgr) {}
0643
0644     void StartTracking(G4Track *aTrack) override {
0645       fMgr.gammaProc->StartTracking(aTrack);
0646
0647       fPreviousStepLength = 0;
0648     }
0649     void EndTracking() override { fMgr.gammaProc->EndTracking(); }
0650
0651     G4double GetPhysicalInteractionLength(const G4Track &track) override {
0652       G4double physIntLength;
0653       G4ForceCondition condition;
0654
0655       fProposedStep = DBL_MAX;
0656       fSelected = -1;
0657
0658       physIntLength = fMgr.gammaProc->PostStepGPIL(track, fPreviousStepLength, &condition);
0659       if (physIntLength < fProposedStep) {
0660         fProposedStep = physIntLength;
0661         fSelected = 0;
0662       }
0663
0664       return fProposedStep;
0665     }
0666
0667     void AlongStepDoIt(G4Track &, G4Step &step, G4TrackVector &) override {
0668       if (step.GetStepLength() == fProposedStep) {
0669         step.GetPostStepPoint()->SetStepStatus(fAlongStepDoItProc);
0670       } else {
0671         // Remember that the step was limited by geometry.
0672         fSelected = -1;
0673       }
0674       fPreviousStepLength = step.GetStepLength();
0675     }
0676
0677     void PostStepDoIt(G4Track &track, G4Step &step, G4TrackVector &secondaries) override {
0678       if (fSelected < 0) {
0679         return;
0680       }
0681       step.GetPostStepPoint()->SetStepStatus(fPostStepDoItProc);
0682
0683       G4VProcess *process = fMgr.gammaProc;
0684       G4VParticleChange *particleChange = fMgr.gammaProc->PostStepDoIt(track, step);
0685
0686       particleChange->UpdateStepForPostStep(&step);
0687       step.UpdateTrack();
0688
0689       int numSecondaries = particleChange->GetNumberOfSecondaries();
0690       for (int i = 0; i < numSecondaries; i++) {
0691         G4Track *secondary = particleChange->GetSecondary(i);
0692         secondary->SetParentID(track.GetTrackID());
0693         secondary->SetCreatorProcess(process);
0694         secondaries.push_back(secondary);
0695       }
0696
0697       track.SetTrackStatus(particleChange->GetTrackStatus());
0698       particleChange->Clear();
0699     }
0700
0701   private:
0702     CMSEmStandardPhysicsTrackingManager &fMgr;
0703     G4double fPreviousStepLength;
0704     G4double fProposedStep;
0705     G4int fSelected;
0706   };
0707
0708   GammaPhysics physics(*this);
0709   TrackingManagerHelper::TrackNeutralParticle(aTrack, physics);
0710 }
0711
0712 void CMSEmStandardPhysicsTrackingManager::HandOverOneTrack(G4Track *aTrack) {
0713   const G4ParticleDefinition *part = aTrack->GetParticleDefinition();
0714
0715   if (part == G4Electron::Definition()) {
0716     TrackElectron(aTrack);
0717   } else if (part == G4Positron::Definition()) {
0718     TrackPositron(aTrack);
0719   } else if (part == G4Gamma::Definition()) {
0720     TrackGamma(aTrack);
0721   }
0722
0723   aTrack->SetTrackStatus(fStopAndKill);
0724   delete aTrack;
0725 }