PhysicsLists/src/CMSEmStandardPhysicsLPM.cc

0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsLPM.h"
0002 #include "SimG4Core/PhysicsLists/interface/EmParticleList.h"
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "G4EmParameters.hh"
0007 #include "G4ParticleTable.hh"
0008
0009 #include "G4ParticleDefinition.hh"
0010 #include "G4LossTableManager.hh"
0011
0012 #include "G4ComptonScattering.hh"
0013 #include "G4GammaConversion.hh"
0014 #include "G4PhotoElectricEffect.hh"
0015
0016 #include "G4hMultipleScattering.hh"
0017 #include "G4eMultipleScattering.hh"
0018 #include "G4MuMultipleScattering.hh"
0019 #include "G4CoulombScattering.hh"
0020 #include "G4eCoulombScatteringModel.hh"
0021 #include "G4WentzelVIModel.hh"
0022 #include "G4UrbanMscModel.hh"
0023
0024 #include "G4eIonisation.hh"
0025 #include "G4eBremsstrahlung.hh"
0026 #include "G4eplusAnnihilation.hh"
0027 #include "G4UAtomicDeexcitation.hh"
0028
0029 #include "G4MuIonisation.hh"
0030 #include "G4MuBremsstrahlung.hh"
0031 #include "G4MuPairProduction.hh"
0032
0033 #include "G4hIonisation.hh"
0034 #include "G4ionIonisation.hh"
0035 #include "G4hBremsstrahlung.hh"
0036 #include "G4hPairProduction.hh"
0037
0038 #include "G4Gamma.hh"
0039 #include "G4Electron.hh"
0040 #include "G4Positron.hh"
0041 #include "G4MuonPlus.hh"
0042 #include "G4MuonMinus.hh"
0043 #include "G4PionPlus.hh"
0044 #include "G4PionMinus.hh"
0045 #include "G4KaonPlus.hh"
0046 #include "G4KaonMinus.hh"
0047 #include "G4Proton.hh"
0048 #include "G4AntiProton.hh"
0049 #include "G4GenericIon.hh"
0050
0051 #include "G4PhysicsListHelper.hh"
0052 #include "G4BuilderType.hh"
0053 #include "G4RegionStore.hh"
0054 #include "G4Region.hh"
0055 #include "G4GammaGeneralProcess.hh"
0056 #include "G4EmBuilder.hh"
0057
0058 #include <CLHEP/Units/SystemOfUnits.h>
0059
0060 CMSEmStandardPhysicsLPM::CMSEmStandardPhysicsLPM(G4int ver, const edm::ParameterSet& p)
0061     : G4VPhysicsConstructor("CMSEmStandard_emm") {
0062   SetVerboseLevel(ver);
0063   // EM parameters specific for this EM physics configuration
0064   G4EmParameters* param = G4EmParameters::Instance();
0065   param->SetDefaults();
0066   param->SetVerbose(ver);
0067   param->SetApplyCuts(true);
0068   param->SetStepFunction(0.8, 1 * CLHEP::mm);
0069   param->SetMscRangeFactor(0.2);
0070   param->SetMscStepLimitType(fMinimal);
0071   SetPhysicsType(bElectromagnetic);
0072   double tcut = p.getParameter<double>("G4TrackingCut") * CLHEP::MeV;
0073   param->SetLowestElectronEnergy(tcut);
0074   param->SetLowestMuHadEnergy(tcut);
0075 }
0076
0077 void CMSEmStandardPhysicsLPM::ConstructParticle() {
0078   // minimal set of particles for EM physics
0079   G4EmBuilder::ConstructMinimalEmSet();
0080 }
0081
0082 void CMSEmStandardPhysicsLPM::ConstructProcess() {
0083   if (verboseLevel > 1) {
0084     edm::LogVerbatim("PhysicsList") << "### " << GetPhysicsName() << " Construct Processes";
0085   }
0086
0087   // This EM builder takes default models of Geant4 10 EMV.
0088   // Multiple scattering by Urban for all particles
0089   // except e+e- below 100 MeV for which the Urban model is used
0090
0091   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0092   G4LossTableManager* man = G4LossTableManager::Instance();
0093
0094   // muon & hadron bremsstrahlung and pair production
0095   G4MuBremsstrahlung* mub = nullptr;
0096   G4MuPairProduction* mup = nullptr;
0097   G4hBremsstrahlung* pib = nullptr;
0098   G4hPairProduction* pip = nullptr;
0099   G4hBremsstrahlung* kb = nullptr;
0100   G4hPairProduction* kp = nullptr;
0101   G4hBremsstrahlung* pb = nullptr;
0102   G4hPairProduction* pp = nullptr;
0103
0104   // muon & hadron multiple scattering
0105   G4MuMultipleScattering* mumsc = nullptr;
0106   G4hMultipleScattering* pimsc = nullptr;
0107   G4hMultipleScattering* kmsc = nullptr;
0108   G4hMultipleScattering* pmsc = nullptr;
0109   G4hMultipleScattering* hmsc = nullptr;
0110
0111   // muon and hadron single scattering
0112   G4CoulombScattering* muss = nullptr;
0113   G4CoulombScattering* piss = nullptr;
0114   G4CoulombScattering* kss = nullptr;
0115
0116   // high energy limit for e+- scattering models and bremsstrahlung
0117   G4double highEnergyLimit = 100 * CLHEP::MeV;
0118
0119   G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
0120   G4Region* bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
0121   if (verboseLevel > 1) {
0122     edm::LogVerbatim("PhysicsList") << "CMSEmStandardPhysicsLPM: HcalRegion " << aRegion << "; HGCalRegion " << bRegion;
0123   }
0124   G4ParticleTable* table = G4ParticleTable::GetParticleTable();
0125   EmParticleList emList;
0126   for (const auto& particleName : emList.PartNames()) {
0127     G4ParticleDefinition* particle = table->FindParticle(particleName);
0128
0129     if (particleName == "gamma") {
0130       G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
0131
0132       if (G4EmParameters::Instance()->GeneralProcessActive()) {
0133         G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
0134         sp->AddEmProcess(pee);
0135         sp->AddEmProcess(new G4ComptonScattering());
0136         sp->AddEmProcess(new G4GammaConversion());
0137         man->SetGammaGeneralProcess(sp);
0138         ph->RegisterProcess(sp, particle);
0139
0140       } else {
0141         ph->RegisterProcess(pee, particle);
0142         ph->RegisterProcess(new G4ComptonScattering(), particle);
0143         ph->RegisterProcess(new G4GammaConversion(), particle);
0144       }
0145
0146     } else if (particleName == "e-") {
0147       G4eIonisation* eioni = new G4eIonisation();
0148
0149       G4eMultipleScattering* msc = new G4eMultipleScattering;
0150       G4UrbanMscModel* msc1 = new G4UrbanMscModel();
0151       G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0152       msc1->SetHighEnergyLimit(highEnergyLimit);
0153       msc2->SetLowEnergyLimit(highEnergyLimit);
0154       msc->SetEmModel(msc1);
0155       msc->SetEmModel(msc2);
0156
0157       // e-/e+ msc for HCAL and HGCAL using the Urban model
0158       if (nullptr != aRegion || nullptr != bRegion) {
0159         G4UrbanMscModel* msc3 = new G4UrbanMscModel();
0160         msc3->SetHighEnergyLimit(highEnergyLimit);
0161         msc3->SetLocked(true);
0162
0163         if (nullptr != aRegion) {
0164           msc->AddEmModel(-1, msc3, aRegion);
0165         }
0166         if (nullptr != bRegion) {
0167           msc->AddEmModel(-1, msc3, bRegion);
0168         }
0169       }
0170
0171       G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0172       G4CoulombScattering* ss = new G4CoulombScattering();
0173       ss->SetEmModel(ssm);
0174       ss->SetMinKinEnergy(highEnergyLimit);
0175       ssm->SetLowEnergyLimit(highEnergyLimit);
0176       ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0177
0178       ph->RegisterProcess(msc, particle);
0179       ph->RegisterProcess(eioni, particle);
0180       ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0181       ph->RegisterProcess(ss, particle);
0182
0183     } else if (particleName == "e+") {
0184       G4eIonisation* eioni = new G4eIonisation();
0185
0186       G4eMultipleScattering* msc = new G4eMultipleScattering;
0187       G4UrbanMscModel* msc1 = new G4UrbanMscModel();
0188       G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0189       msc1->SetHighEnergyLimit(highEnergyLimit);
0190       msc2->SetLowEnergyLimit(highEnergyLimit);
0191       msc->SetEmModel(msc1);
0192       msc->SetEmModel(msc2);
0193
0194       // e-/e+ msc for HCAL and HGCAL using the Urban model
0195       if (nullptr != aRegion || nullptr != bRegion) {
0196         G4UrbanMscModel* msc3 = new G4UrbanMscModel();
0197         msc3->SetHighEnergyLimit(highEnergyLimit);
0198         msc3->SetLocked(true);
0199
0200         if (nullptr != aRegion) {
0201           msc->AddEmModel(-1, msc3, aRegion);
0202         }
0203         if (nullptr != bRegion) {
0204           msc->AddEmModel(-1, msc3, bRegion);
0205         }
0206       }
0207
0208       G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0209       G4CoulombScattering* ss = new G4CoulombScattering();
0210       ss->SetEmModel(ssm);
0211       ss->SetMinKinEnergy(highEnergyLimit);
0212       ssm->SetLowEnergyLimit(highEnergyLimit);
0213       ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0214
0215       ph->RegisterProcess(msc, particle);
0216       ph->RegisterProcess(eioni, particle);
0217       ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0218       ph->RegisterProcess(new G4eplusAnnihilation(), particle);
0219       ph->RegisterProcess(ss, particle);
0220
0221     } else if (particleName == "mu+" || particleName == "mu-") {
0222       if (nullptr == mub) {
0223         mub = new G4MuBremsstrahlung();
0224         mup = new G4MuPairProduction();
0225         mumsc = new G4MuMultipleScattering();
0226         mumsc->SetEmModel(new G4WentzelVIModel());
0227         muss = new G4CoulombScattering();
0228       }
0229       ph->RegisterProcess(mumsc, particle);
0230       ph->RegisterProcess(new G4MuIonisation(), particle);
0231       ph->RegisterProcess(mub, particle);
0232       ph->RegisterProcess(mup, particle);
0233       ph->RegisterProcess(muss, particle);
0234
0235     } else if (particleName == "alpha" || particleName == "He3") {
0236       ph->RegisterProcess(new G4hMultipleScattering(), particle);
0237       ph->RegisterProcess(new G4ionIonisation(), particle);
0238
0239     } else if (particleName == "GenericIon") {
0240       if (nullptr == hmsc) {
0241         hmsc = new G4hMultipleScattering("ionmsc");
0242       }
0243       ph->RegisterProcess(hmsc, particle);
0244       ph->RegisterProcess(new G4ionIonisation(), particle);
0245
0246     } else if (particleName == "pi+" || particleName == "pi-") {
0247       if (nullptr == pib) {
0248         pib = new G4hBremsstrahlung();
0249         pip = new G4hPairProduction();
0250         pimsc = new G4hMultipleScattering();
0251         pimsc->SetEmModel(new G4WentzelVIModel());
0252         piss = new G4CoulombScattering();
0253       }
0254       ph->RegisterProcess(pimsc, particle);
0255       ph->RegisterProcess(new G4hIonisation(), particle);
0256       ph->RegisterProcess(pib, particle);
0257       ph->RegisterProcess(pip, particle);
0258       ph->RegisterProcess(piss, particle);
0259
0260     } else if (particleName == "kaon+" || particleName == "kaon-") {
0261       if (nullptr == kb) {
0262         kb = new G4hBremsstrahlung();
0263         kp = new G4hPairProduction();
0264         kmsc = new G4hMultipleScattering();
0265         kmsc->SetEmModel(new G4WentzelVIModel());
0266         kss = new G4CoulombScattering();
0267       }
0268       ph->RegisterProcess(kmsc, particle);
0269       ph->RegisterProcess(new G4hIonisation(), particle);
0270       ph->RegisterProcess(kb, particle);
0271       ph->RegisterProcess(kp, particle);
0272       ph->RegisterProcess(kss, particle);
0273
0274     } else if (particleName == "proton" || particleName == "anti_proton") {
0275       if (nullptr == pb) {
0276         pb = new G4hBremsstrahlung();
0277         pp = new G4hPairProduction();
0278       }
0279       pmsc = new G4hMultipleScattering();
0280       pmsc->SetEmModel(new G4WentzelVIModel());
0281
0282       ph->RegisterProcess(pmsc, particle);
0283       ph->RegisterProcess(new G4hIonisation(), particle);
0284       ph->RegisterProcess(pb, particle);
0285       ph->RegisterProcess(pp, particle);
0286       ph->RegisterProcess(new G4CoulombScattering(), particle);
0287
0288     } else if (particle->GetPDGCharge() != 0.0) {
0289       if (nullptr == hmsc) {
0290         hmsc = new G4hMultipleScattering("ionmsc");
0291       }
0292       ph->RegisterProcess(hmsc, particle);
0293       ph->RegisterProcess(new G4hIonisation(), particle);
0294     }
0295   }
0296   edm::LogVerbatim("PhysicsList") << "CMSEmStandardPhysicsLPM: EM physics is instantiated";
0297 }