PhysicsLists/src/CMSEmStandardPhysicsXS.cc

0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "G4SystemOfUnits.hh"
0005 #include "G4ParticleDefinition.hh"
0006 #include "G4LossTableManager.hh"
0007 #include "G4EmParameters.hh"
0008 #include "G4EmBuilder.hh"
0009
0010 #include "G4ComptonScattering.hh"
0011 #include "G4GammaConversion.hh"
0012 #include "G4PhotoElectricEffect.hh"
0013 #include "G4RayleighScattering.hh"
0014 #include "G4PEEffectFluoModel.hh"
0015 #include "G4KleinNishinaModel.hh"
0016 #include "G4LowEPComptonModel.hh"
0017 #include "G4BetheHeitler5DModel.hh"
0018 #include "G4LivermorePhotoElectricModel.hh"
0019
0020 #include "G4eMultipleScattering.hh"
0021 #include "G4hMultipleScattering.hh"
0022 #include "G4MscStepLimitType.hh"
0023 #include "G4UrbanMscModel.hh"
0024 #include "G4GoudsmitSaundersonMscModel.hh"
0025 #include "G4DummyModel.hh"
0026 #include "G4WentzelVIModel.hh"
0027 #include "G4CoulombScattering.hh"
0028 #include "G4eCoulombScatteringModel.hh"
0029
0030 #include "G4eIonisation.hh"
0031 #include "G4eBremsstrahlung.hh"
0032 #include "G4Generator2BS.hh"
0033 #include "G4SeltzerBergerModel.hh"
0034 #include "G4ePairProduction.hh"
0035 #include "G4UniversalFluctuation.hh"
0036
0037 #include "G4eplusAnnihilation.hh"
0038
0039 #include "G4hIonisation.hh"
0040 #include "G4ionIonisation.hh"
0041
0042 #include "G4Gamma.hh"
0043 #include "G4Electron.hh"
0044 #include "G4Positron.hh"
0045 #include "G4GenericIon.hh"
0046
0047 #include "G4PhysicsListHelper.hh"
0048 #include "G4BuilderType.hh"
0049 #include "G4GammaGeneralProcess.hh"
0050
0051 #include "G4Version.hh"
0052 #if G4VERSION_NUMBER >= 1110
0053 #include "G4ProcessManager.hh"
0054 #include "G4TransportationWithMsc.hh"
0055 #endif
0056
0057 #include "G4RegionStore.hh"
0058 #include "G4Region.hh"
0059 #include "G4GammaGeneralProcess.hh"
0060
0061 #include "G4SystemOfUnits.hh"
0062
0063 CMSEmStandardPhysicsXS::CMSEmStandardPhysicsXS(G4int ver, const edm::ParameterSet& p)
0064     : G4VPhysicsConstructor("CMSEmStandard_emn") {
0065   SetVerboseLevel(ver);
0066   // EM parameters specific for this EM physics configuration
0067   G4EmParameters* param = G4EmParameters::Instance();
0068   param->SetDefaults();
0069   param->SetVerbose(ver);
0070   param->SetApplyCuts(true);
0071   param->SetMinEnergy(100 * CLHEP::eV);
0072   param->SetNumberOfBinsPerDecade(20);
0073   param->SetStepFunction(0.8, 1 * CLHEP::mm);
0074   param->SetMscRangeFactor(0.2);
0075   param->SetMscStepLimitType(fMinimal);
0076   param->SetFluo(true);
0077   param->SetUseMottCorrection(true);  // use Mott-correction for e-/e+ msc gs
0078   SetPhysicsType(bElectromagnetic);
0079   fRangeFactor = p.getParameter<double>("G4MscRangeFactor");
0080   fGeomFactor = p.getParameter<double>("G4MscGeomFactor");
0081   fSafetyFactor = p.getParameter<double>("G4MscSafetyFactor");
0082   fLambdaLimit = p.getParameter<double>("G4MscLambdaLimit") * CLHEP::mm;
0083   std::string msc = p.getParameter<std::string>("G4MscStepLimit");
0084   fStepLimitType = fUseSafety;
0085   if (msc == "UseSafetyPlus") {
0086     fStepLimitType = fUseSafetyPlus;
0087   }
0088   if (msc == "Minimal") {
0089     fStepLimitType = fMinimal;
0090   }
0091   double tcut = p.getParameter<double>("G4TrackingCut") * CLHEP::MeV;
0092   param->SetLowestElectronEnergy(tcut);
0093   param->SetLowestMuHadEnergy(tcut);
0094 }
0095
0096 void CMSEmStandardPhysicsXS::ConstructParticle() {
0097   // minimal set of particles for EM physics
0098   G4EmBuilder::ConstructMinimalEmSet();
0099 }
0100
0101 void CMSEmStandardPhysicsXS::ConstructProcess() {
0102   if (verboseLevel > 0) {
0103     edm::LogVerbatim("PhysicsList") << "### " << GetPhysicsName() << " Construct Processes";
0104   }
0105
0106   // This EM builder takes default models of Geant4 10 EMV.
0107   // Multiple scattering by Urban for all particles
0108   // except e+e- below 100 MeV for which the Urban model is used
0109   G4EmBuilder::PrepareEMPhysics();
0110
0111   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0112
0113   // processes used by several particles
0114   G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc");
0115   G4NuclearStopping* pnuc(nullptr);
0116
0117   // high energy limit for e+- scattering models
0118   G4double highEnergyLimit = G4EmParameters::Instance()->MscEnergyLimit();
0119   const G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
0120   const G4Region* bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
0121
0122   // Add gamma EM Processes
0123   G4ParticleDefinition* particle = G4Gamma::Gamma();
0124
0125   // Photoelectric
0126   G4PhotoElectricEffect* pe = new G4PhotoElectricEffect();
0127   G4VEmModel* theLivermorePEModel = new G4LivermorePhotoElectricModel();
0128   pe->SetEmModel(theLivermorePEModel);
0129
0130   // Compton scattering
0131   G4ComptonScattering* cs = new G4ComptonScattering;
0132   cs->SetEmModel(new G4KleinNishinaModel());
0133
0134   // Gamma conversion
0135   G4GammaConversion* gc = new G4GammaConversion();
0136   G4VEmModel* conv = new G4BetheHeitler5DModel();
0137   gc->SetEmModel(conv);
0138
0139   if (G4EmParameters::Instance()->GeneralProcessActive()) {
0140     G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
0141     sp->AddEmProcess(pe);
0142     sp->AddEmProcess(cs);
0143     sp->AddEmProcess(gc);
0144     sp->AddEmProcess(new G4RayleighScattering());
0145     G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
0146     ph->RegisterProcess(sp, particle);
0147   } else {
0148     ph->RegisterProcess(pe, particle);
0149     ph->RegisterProcess(cs, particle);
0150     ph->RegisterProcess(gc, particle);
0151     ph->RegisterProcess(new G4RayleighScattering(), particle);
0152   }
0153
0154   // e-
0155   particle = G4Electron::Electron();
0156
0157   // multiple scattering
0158   G4UrbanMscModel* msc1 = new G4UrbanMscModel();
0159   G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0160   msc1->SetHighEnergyLimit(highEnergyLimit);
0161   msc2->SetLowEnergyLimit(highEnergyLimit);
0162
0163   // msc for HCAL using the Urban model
0164   G4UrbanMscModel* msc4 = nullptr;
0165   if (nullptr != aRegion) {
0166     msc4 = new G4UrbanMscModel();
0167     msc4->SetHighEnergyLimit(highEnergyLimit);
0168     msc4->SetRangeFactor(fRangeFactor);
0169     msc4->SetGeomFactor(fGeomFactor);
0170     msc4->SetSafetyFactor(fSafetyFactor);
0171     msc4->SetLambdaLimit(fLambdaLimit);
0172     msc4->SetStepLimitType(fStepLimitType);
0173     msc4->SetLocked(true);
0174   }
0175
0176   // msc GS with Mott-correction
0177   G4GoudsmitSaundersonMscModel* msc3 = nullptr;
0178   if (nullptr != bRegion) {
0179     msc3 = new G4GoudsmitSaundersonMscModel();
0180     msc3->SetHighEnergyLimit(highEnergyLimit);
0181     msc3->SetRangeFactor(0.08);
0182     msc3->SetSkin(3.);
0183     msc3->SetStepLimitType(fUseSafetyPlus);
0184     msc3->SetLocked(true);
0185   }
0186
0187 #if G4VERSION_NUMBER >= 1110
0188   G4TransportationWithMscType transportationWithMsc = G4EmParameters::Instance()->TransportationWithMsc();
0189   if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0190     G4ProcessManager* procManager = particle->GetProcessManager();
0191     // Remove default G4Transportation and replace with G4TransportationWithMsc.
0192     G4VProcess* removed = procManager->RemoveProcess(0);
0193     if (removed->GetProcessName() != "Transportation") {
0194       G4Exception("CMSEmStandardPhysics::ConstructProcess",
0195                   "em0050",
0196                   FatalException,
0197                   "replaced process is not G4Transportation!");
0198     }
0199     G4TransportationWithMsc* transportWithMsc =
0200         new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0201     if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0202       transportWithMsc->SetMultipleSteps(true);
0203     }
0204     transportWithMsc->AddMscModel(msc1);
0205     transportWithMsc->AddMscModel(msc2);
0206     if (msc4 != nullptr) {
0207       transportWithMsc->AddMscModel(msc4, -1, aRegion);
0208     }
0209     if (msc3 != nullptr) {
0210       transportWithMsc->AddMscModel(msc3, -1, bRegion);
0211     }
0212     procManager->AddProcess(transportWithMsc, -1, 0, 0);
0213   } else
0214 #endif
0215   {
0216     // Register as a separate process.
0217     G4eMultipleScattering* msc = new G4eMultipleScattering;
0218     msc->SetEmModel(msc1);
0219     msc->SetEmModel(msc2);
0220     if (msc4 != nullptr) {
0221       msc->AddEmModel(-1, msc4, aRegion);
0222     }
0223     if (msc3 != nullptr) {
0224       msc->AddEmModel(-1, msc3, bRegion);
0225     }
0226     ph->RegisterProcess(msc, particle);
0227   }
0228
0229   // single scattering
0230   G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0231   G4CoulombScattering* ss = new G4CoulombScattering();
0232   ss->SetEmModel(ssm);
0233   ss->SetMinKinEnergy(highEnergyLimit);
0234   ssm->SetLowEnergyLimit(highEnergyLimit);
0235   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0236
0237   // ionisation
0238   G4eIonisation* eioni = new G4eIonisation();
0239
0240   // bremsstrahlung
0241   G4eBremsstrahlung* brem = new G4eBremsstrahlung();
0242   G4SeltzerBergerModel* br1 = new G4SeltzerBergerModel();
0243   G4eBremsstrahlungRelModel* br2 = new G4eBremsstrahlungRelModel();
0244   br1->SetAngularDistribution(new G4Generator2BS());
0245   br2->SetAngularDistribution(new G4Generator2BS());
0246   brem->SetEmModel(br1);
0247   brem->SetEmModel(br2);
0248   br1->SetHighEnergyLimit(CLHEP::GeV);
0249
0250   G4ePairProduction* ee = new G4ePairProduction();
0251
0252   // register processes
0253   ph->RegisterProcess(eioni, particle);
0254   ph->RegisterProcess(brem, particle);
0255   ph->RegisterProcess(ee, particle);
0256   ph->RegisterProcess(ss, particle);
0257
0258   // e+
0259   particle = G4Positron::Positron();
0260
0261   // multiple scattering
0262   msc1 = new G4UrbanMscModel();
0263   msc2 = new G4WentzelVIModel();
0264   msc1->SetHighEnergyLimit(highEnergyLimit);
0265   msc2->SetLowEnergyLimit(highEnergyLimit);
0266
0267   // msc for HCAL using the Urban model
0268   if (nullptr != aRegion) {
0269     msc4 = new G4UrbanMscModel();
0270     msc4->SetHighEnergyLimit(highEnergyLimit);
0271     msc4->SetRangeFactor(fRangeFactor);
0272     msc4->SetGeomFactor(fGeomFactor);
0273     msc4->SetSafetyFactor(fSafetyFactor);
0274     msc4->SetLambdaLimit(fLambdaLimit);
0275     msc4->SetStepLimitType(fStepLimitType);
0276     msc4->SetLocked(true);
0277   }
0278
0279   // msc GS with Mott-correction
0280   if (nullptr != bRegion) {
0281     msc3 = new G4GoudsmitSaundersonMscModel();
0282     msc3->SetHighEnergyLimit(highEnergyLimit);
0283     msc3->SetRangeFactor(0.08);
0284     msc3->SetSkin(3.);
0285     msc3->SetStepLimitType(fUseSafetyPlus);
0286     msc3->SetLocked(true);
0287   }
0288
0289 #if G4VERSION_NUMBER >= 1110
0290   if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0291     G4ProcessManager* procManager = particle->GetProcessManager();
0292     // Remove default G4Transportation and replace with G4TransportationWithMsc.
0293     G4VProcess* removed = procManager->RemoveProcess(0);
0294     if (removed->GetProcessName() != "Transportation") {
0295       G4Exception("CMSEmStandardPhysics::ConstructProcess",
0296                   "em0050",
0297                   FatalException,
0298                   "replaced process is not G4Transportation!");
0299     }
0300     G4TransportationWithMsc* transportWithMsc =
0301         new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0302     if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0303       transportWithMsc->SetMultipleSteps(true);
0304     }
0305     transportWithMsc->AddMscModel(msc1);
0306     transportWithMsc->AddMscModel(msc2);
0307     if (msc4 != nullptr) {
0308       transportWithMsc->AddMscModel(msc4, -1, aRegion);
0309     }
0310     if (msc3 != nullptr) {
0311       transportWithMsc->AddMscModel(msc3, -1, bRegion);
0312     }
0313     procManager->AddProcess(transportWithMsc, -1, 0, 0);
0314   } else
0315 #endif
0316   {
0317     // Register as a separate process.
0318     G4eMultipleScattering* msc = new G4eMultipleScattering;
0319     msc->SetEmModel(msc1);
0320     msc->SetEmModel(msc2);
0321     if (msc4 != nullptr) {
0322       msc->AddEmModel(-1, msc4, aRegion);
0323     }
0324     if (msc3 != nullptr) {
0325       msc->AddEmModel(-1, msc3, bRegion);
0326     }
0327     ph->RegisterProcess(msc, particle);
0328   }
0329
0330   // single scattering
0331   ssm = new G4eCoulombScatteringModel();
0332   ss = new G4CoulombScattering();
0333   ss->SetEmModel(ssm);
0334   ss->SetMinKinEnergy(highEnergyLimit);
0335   ssm->SetLowEnergyLimit(highEnergyLimit);
0336   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0337
0338   // ionisation
0339   eioni = new G4eIonisation();
0340
0341   // bremsstrahlung
0342   brem = new G4eBremsstrahlung();
0343   br1 = new G4SeltzerBergerModel();
0344   br2 = new G4eBremsstrahlungRelModel();
0345   br1->SetAngularDistribution(new G4Generator2BS());
0346   br2->SetAngularDistribution(new G4Generator2BS());
0347   brem->SetEmModel(br1);
0348   brem->SetEmModel(br2);
0349   br1->SetHighEnergyLimit(CLHEP::GeV);
0350
0351   // register processes
0352   ph->RegisterProcess(eioni, particle);
0353   ph->RegisterProcess(brem, particle);
0354   ph->RegisterProcess(ee, particle);
0355   ph->RegisterProcess(new G4eplusAnnihilation(), particle);
0356   ph->RegisterProcess(ss, particle);
0357
0358   // generic ion
0359   particle = G4GenericIon::GenericIon();
0360   G4ionIonisation* ionIoni = new G4ionIonisation();
0361   ph->RegisterProcess(hmsc, particle);
0362   ph->RegisterProcess(ionIoni, particle);
0363
0364   // muons, hadrons, ions
0365   G4EmBuilder::ConstructCharged(hmsc, pnuc);
0366 }