PhysicsLists/src/CMSEmStandardPhysicsXS.cc

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