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