PhysicsLists/src/CMSEmStandardPhysics.cc

0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysics.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "G4SystemOfUnits.hh"
0005 #include "G4ParticleDefinition.hh"
0006 #include "G4EmParameters.hh"
0007 #include "G4EmBuilder.hh"
0008
0009 #include "G4ComptonScattering.hh"
0010 #include "G4GammaConversion.hh"
0011 #include "G4PhotoElectricEffect.hh"
0012 #include "G4LivermorePhotoElectricModel.hh"
0013
0014 #include "G4MscStepLimitType.hh"
0015
0016 #include "G4eMultipleScattering.hh"
0017 #include "G4hMultipleScattering.hh"
0018 #include "G4eCoulombScatteringModel.hh"
0019 #include "G4CoulombScattering.hh"
0020 #include "G4WentzelVIModel.hh"
0021 #include "G4UrbanMscModel.hh"
0022
0023 #include "G4eIonisation.hh"
0024 #include "G4eBremsstrahlung.hh"
0025 #include "G4eplusAnnihilation.hh"
0026
0027 #include "G4hIonisation.hh"
0028 #include "G4ionIonisation.hh"
0029
0030 #include "G4ParticleTable.hh"
0031 #include "G4Gamma.hh"
0032 #include "G4Electron.hh"
0033 #include "G4Positron.hh"
0034 #include "G4GenericIon.hh"
0035
0036 #include "G4PhysicsListHelper.hh"
0037 #include "G4BuilderType.hh"
0038 #include "G4GammaGeneralProcess.hh"
0039 #include "G4LossTableManager.hh"
0040
0041 #include "G4Version.hh"
0042 #if G4VERSION_NUMBER >= 1110
0043 #include "G4ProcessManager.hh"
0044 #include "G4TransportationWithMsc.hh"
0045 #endif
0046
0047 #include "G4RegionStore.hh"
0048 #include "G4Region.hh"
0049 #include <string>
0050
0051 CMSEmStandardPhysics::CMSEmStandardPhysics(G4int ver, const edm::ParameterSet& p)
0052     : G4VPhysicsConstructor("CMSEmStandard_emm") {
0053   SetVerboseLevel(ver);
0054   G4EmParameters* param = G4EmParameters::Instance();
0055   param->SetDefaults();
0056   param->SetVerbose(ver);
0057   param->SetApplyCuts(true);
0058   param->SetStepFunction(0.8, 1 * CLHEP::mm);
0059   param->SetMscRangeFactor(0.2);
0060   param->SetMscStepLimitType(fMinimal);
0061   param->SetFluo(false);
0062   SetPhysicsType(bElectromagnetic);
0063   fRangeFactor = p.getParameter<double>("G4MscRangeFactor");
0064   fGeomFactor = p.getParameter<double>("G4MscGeomFactor");
0065   fSafetyFactor = p.getParameter<double>("G4MscSafetyFactor");
0066   fLambdaLimit = p.getParameter<double>("G4MscLambdaLimit") * CLHEP::mm;
0067   std::string msc = p.getParameter<std::string>("G4MscStepLimit");
0068   fStepLimitType = fUseSafety;
0069   if (msc == "UseSafetyPlus") {
0070     fStepLimitType = fUseSafetyPlus;
0071   }
0072   if (msc == "Minimal") {
0073     fStepLimitType = fMinimal;
0074   }
0075   double tcut = p.getParameter<double>("G4TrackingCut") * CLHEP::MeV;
0076   param->SetLowestElectronEnergy(tcut);
0077   param->SetLowestMuHadEnergy(tcut);
0078 }
0079
0080 CMSEmStandardPhysics::~CMSEmStandardPhysics() {}
0081
0082 void CMSEmStandardPhysics::ConstructParticle() {
0083   // minimal set of particles for EM physics
0084   G4EmBuilder::ConstructMinimalEmSet();
0085 }
0086
0087 void CMSEmStandardPhysics::ConstructProcess() {
0088   if (verboseLevel > 0) {
0089     edm::LogVerbatim("PhysicsList") << "### " << GetPhysicsName() << " Construct EM Processes";
0090   }
0091
0092   // This EM builder takes default models of Geant4 10 EMV.
0093   // Multiple scattering by WentzelVI for all particles except:
0094   //   a) e+e- below 100 MeV for which the Urban model is used
0095   //   b) ions for which Urban model is used
0096   G4EmBuilder::PrepareEMPhysics();
0097
0098   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0099   // processes used by several particles
0100   G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc");
0101   G4NuclearStopping* pnuc(nullptr);
0102
0103   // high energy limit for e+- scattering models
0104   G4double highEnergyLimit = G4EmParameters::Instance()->MscEnergyLimit();
0105
0106   const G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
0107   const G4Region* bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
0108
0109   // Add gamma EM Processes
0110   G4ParticleDefinition* particle = G4Gamma::Gamma();
0111
0112   G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
0113
0114   if (G4EmParameters::Instance()->GeneralProcessActive()) {
0115     G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
0116     sp->AddEmProcess(pee);
0117     sp->AddEmProcess(new G4ComptonScattering());
0118     sp->AddEmProcess(new G4GammaConversion());
0119     G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
0120     ph->RegisterProcess(sp, particle);
0121
0122   } else {
0123     ph->RegisterProcess(pee, particle);
0124     ph->RegisterProcess(new G4ComptonScattering(), particle);
0125     ph->RegisterProcess(new G4GammaConversion(), particle);
0126   }
0127
0128   // e-
0129   particle = G4Electron::Electron();
0130
0131   G4eIonisation* eioni = new G4eIonisation();
0132
0133   G4UrbanMscModel* msc1 = new G4UrbanMscModel();
0134   G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0135   msc1->SetHighEnergyLimit(highEnergyLimit);
0136   msc2->SetLowEnergyLimit(highEnergyLimit);
0137
0138   // e-/e+ msc for HCAL and HGCAL using the Urban model
0139   G4UrbanMscModel* msc3 = nullptr;
0140   if (nullptr != aRegion || nullptr != bRegion) {
0141     msc3 = new G4UrbanMscModel();
0142     msc3->SetHighEnergyLimit(highEnergyLimit);
0143     msc3->SetRangeFactor(fRangeFactor);
0144     msc3->SetGeomFactor(fGeomFactor);
0145     msc3->SetSafetyFactor(fSafetyFactor);
0146     msc3->SetLambdaLimit(fLambdaLimit);
0147     msc3->SetStepLimitType(fStepLimitType);
0148     msc3->SetLocked(true);
0149   }
0150
0151 #if G4VERSION_NUMBER >= 1110
0152   G4TransportationWithMscType transportationWithMsc = G4EmParameters::Instance()->TransportationWithMsc();
0153   if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0154     G4ProcessManager* procManager = particle->GetProcessManager();
0155     // Remove default G4Transportation and replace with G4TransportationWithMsc.
0156     G4VProcess* removed = procManager->RemoveProcess(0);
0157     if (removed->GetProcessName() != "Transportation") {
0158       G4Exception("CMSEmStandardPhysics::ConstructProcess",
0159                   "em0050",
0160                   FatalException,
0161                   "replaced process is not G4Transportation!");
0162     }
0163     G4TransportationWithMsc* transportWithMsc =
0164         new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0165     if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0166       transportWithMsc->SetMultipleSteps(true);
0167     }
0168     transportWithMsc->AddMscModel(msc1);
0169     transportWithMsc->AddMscModel(msc2);
0170     if (nullptr != aRegion) {
0171       transportWithMsc->AddMscModel(msc3, -1, aRegion);
0172     }
0173     if (nullptr != bRegion) {
0174       transportWithMsc->AddMscModel(msc3, -1, bRegion);
0175     }
0176     procManager->AddProcess(transportWithMsc, -1, 0, 0);
0177   } else
0178 #endif
0179   {
0180     // Register as a separate process.
0181     G4eMultipleScattering* msc = new G4eMultipleScattering;
0182     msc->SetEmModel(msc1);
0183     msc->SetEmModel(msc2);
0184     if (nullptr != aRegion) {
0185       msc->AddEmModel(-1, msc3, aRegion);
0186     }
0187     if (nullptr != bRegion) {
0188       msc->AddEmModel(-1, msc3, bRegion);
0189     }
0190     ph->RegisterProcess(msc, particle);
0191   }
0192
0193   // single scattering
0194   G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0195   G4CoulombScattering* ss = new G4CoulombScattering();
0196   ss->SetEmModel(ssm);
0197   ss->SetMinKinEnergy(highEnergyLimit);
0198   ssm->SetLowEnergyLimit(highEnergyLimit);
0199   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0200
0201   ph->RegisterProcess(eioni, particle);
0202   ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0203   ph->RegisterProcess(ss, particle);
0204
0205   // e+
0206   particle = G4Positron::Positron();
0207   eioni = new G4eIonisation();
0208
0209   msc1 = new G4UrbanMscModel();
0210   msc2 = new G4WentzelVIModel();
0211   msc1->SetHighEnergyLimit(highEnergyLimit);
0212   msc2->SetLowEnergyLimit(highEnergyLimit);
0213
0214   // e-/e+ msc for HCAL and HGCAL using the Urban model
0215   if (nullptr != aRegion || nullptr != bRegion) {
0216     msc3 = new G4UrbanMscModel();
0217     msc3->SetHighEnergyLimit(highEnergyLimit);
0218     msc3->SetRangeFactor(fRangeFactor);
0219     msc3->SetGeomFactor(fGeomFactor);
0220     msc3->SetSafetyFactor(fSafetyFactor);
0221     msc3->SetLambdaLimit(fLambdaLimit);
0222     msc3->SetStepLimitType(fStepLimitType);
0223     msc3->SetLocked(true);
0224   }
0225
0226 #if G4VERSION_NUMBER >= 1110
0227   if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0228     G4ProcessManager* procManager = particle->GetProcessManager();
0229     // Remove default G4Transportation and replace with G4TransportationWithMsc.
0230     G4VProcess* removed = procManager->RemoveProcess(0);
0231     if (removed->GetProcessName() != "Transportation") {
0232       G4Exception("CMSEmStandardPhysics::ConstructProcess",
0233                   "em0050",
0234                   FatalException,
0235                   "replaced process is not G4Transportation!");
0236     }
0237     G4TransportationWithMsc* transportWithMsc =
0238         new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0239     if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0240       transportWithMsc->SetMultipleSteps(true);
0241     }
0242     transportWithMsc->AddMscModel(msc1);
0243     transportWithMsc->AddMscModel(msc2);
0244     if (nullptr != aRegion) {
0245       transportWithMsc->AddMscModel(msc3, -1, aRegion);
0246     }
0247     if (nullptr != bRegion) {
0248       transportWithMsc->AddMscModel(msc3, -1, bRegion);
0249     }
0250     procManager->AddProcess(transportWithMsc, -1, 0, 0);
0251   } else
0252 #endif
0253   {
0254     // Register as a separate process.
0255     G4eMultipleScattering* msc = new G4eMultipleScattering;
0256     msc->SetEmModel(msc1);
0257     msc->SetEmModel(msc2);
0258     if (nullptr != aRegion) {
0259       msc->AddEmModel(-1, msc3, aRegion);
0260     }
0261     if (nullptr != bRegion) {
0262       msc->AddEmModel(-1, msc3, bRegion);
0263     }
0264     ph->RegisterProcess(msc, particle);
0265   }
0266
0267   // single scattering
0268   ssm = new G4eCoulombScatteringModel();
0269   ss = new G4CoulombScattering();
0270   ss->SetEmModel(ssm);
0271   ss->SetMinKinEnergy(highEnergyLimit);
0272   ssm->SetLowEnergyLimit(highEnergyLimit);
0273   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0274
0275   ph->RegisterProcess(eioni, particle);
0276   ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0277   ph->RegisterProcess(new G4eplusAnnihilation(), particle);
0278   ph->RegisterProcess(ss, particle);
0279
0280   // generic ion
0281   particle = G4GenericIon::GenericIon();
0282   G4ionIonisation* ionIoni = new G4ionIonisation();
0283   ph->RegisterProcess(hmsc, particle);
0284   ph->RegisterProcess(ionIoni, particle);
0285
0286   // muons, hadrons ions
0287   G4EmBuilder::ConstructCharged(hmsc, pnuc);
0288 }