PhysicsLists/src/CMSEmStandardPhysics.cc

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