Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-15 01:08:02

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