PhysicsLists/src/CMSEmStandardPhysicsTrackingManager.cc

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