PhysicsLists/src/CMSEmStandardPhysics.cc

0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysics.h"
0002 #include "SimG4Core/Physics/interface/CMSG4TrackInterface.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 #include "G4HepEmTrackingManager.hh"
0048 #include "G4HepEmConfig.hh"
0049
0050 CMSEmStandardPhysics::CMSEmStandardPhysics(G4int ver, const edm::ParameterSet& p)
0051     : G4VPhysicsConstructor("CMSEmStandard_emm") {
0052   SetVerboseLevel(ver);
0053   // EM parameters specific for this EM physics configuration
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   fG4HepEmActive = p.getParameter<bool>("G4HepEmActive");
0079   std::string type = p.getParameter<std::string>("type");
0080   if (type == "SimG4Core/Physics/FTFP_BERT_EMH") {
0081     int id = CMSG4TrackInterface::instance()->getThreadID();
0082     edm::LogVerbatim("PhysicsList") << "EMM -> EMH: Forcing usage of G4HepEm; threadID=" << id;
0083     fG4HepEmActive = true;
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     int id = CMSG4TrackInterface::instance()->getThreadID();
0095     edm::LogVerbatim("PhysicsList") << "### " << GetPhysicsName() << " Construct EM Processes; EMH=" << fG4HepEmActive
0096                                     << " threadID=" << id;
0097   }
0098
0099   // This EM builder takes default models of Geant4 10 EMV.
0100   // Multiple scattering by WentzelVI for all particles except:
0101   //   a) e+e- below 100 MeV for which the Urban model is used
0102   //   b) ions for which Urban model is used
0103   G4EmBuilder::PrepareEMPhysics();
0104
0105   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
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   if (!fG4HepEmActive) {
0118     // Add gamma EM Processes
0119     G4ParticleDefinition* particle = G4Gamma::Gamma();
0120
0121     G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
0122
0123     if (param->GeneralProcessActive()) {
0124       G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
0125       sp->AddEmProcess(pee);
0126       sp->AddEmProcess(new G4ComptonScattering());
0127       sp->AddEmProcess(new G4GammaConversion());
0128       G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
0129       ph->RegisterProcess(sp, particle);
0130
0131     } else {
0132       ph->RegisterProcess(pee, particle);
0133       ph->RegisterProcess(new G4ComptonScattering(), particle);
0134       ph->RegisterProcess(new G4GammaConversion(), particle);
0135     }
0136
0137     // e-
0138     particle = G4Electron::Electron();
0139
0140     G4UrbanMscModel* msc1 = new G4UrbanMscModel();
0141     G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0142     msc1->SetHighEnergyLimit(highEnergyLimit);
0143     msc2->SetLowEnergyLimit(highEnergyLimit);
0144
0145     // e-/e+ msc for HCAL and HGCAL using the Urban model
0146     G4UrbanMscModel* msc3 = nullptr;
0147     if (nullptr != aRegion || nullptr != bRegion) {
0148       msc3 = new G4UrbanMscModel();
0149       msc3->SetHighEnergyLimit(highEnergyLimit);
0150       msc3->SetRangeFactor(fRangeFactor);
0151       msc3->SetGeomFactor(fGeomFactor);
0152       msc3->SetSafetyFactor(fSafetyFactor);
0153       msc3->SetLambdaLimit(fLambdaLimit);
0154       msc3->SetStepLimitType(fStepLimitType);
0155       msc3->SetLocked(true);
0156     }
0157
0158     G4TransportationWithMscType transportationWithMsc = param->TransportationWithMsc();
0159     if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0160       // Remove default G4Transportation and replace with G4TransportationWithMsc.
0161       G4ProcessManager* procManager = particle->GetProcessManager();
0162       G4VProcess* removed = procManager->RemoveProcess(0);
0163       if (removed->GetProcessName() != "Transportation") {
0164         G4Exception("CMSEmStandardPhysics::ConstructProcess",
0165                     "em0050",
0166                     FatalException,
0167                     "replaced process is not G4Transportation!");
0168       }
0169       G4TransportationWithMsc* transportWithMsc =
0170           new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0171       if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0172         transportWithMsc->SetMultipleSteps(true);
0173       }
0174       transportWithMsc->AddMscModel(msc1);
0175       transportWithMsc->AddMscModel(msc2);
0176       if (nullptr != aRegion) {
0177         transportWithMsc->AddMscModel(msc3, -1, aRegion);
0178       }
0179       if (nullptr != bRegion) {
0180         transportWithMsc->AddMscModel(msc3, -1, bRegion);
0181       }
0182       procManager->AddProcess(transportWithMsc, -1, 0, 0);
0183     } else {
0184       // Multiple scattering is registered as a separate process
0185       G4eMultipleScattering* msc = new G4eMultipleScattering;
0186       msc->SetEmModel(msc1);
0187       msc->SetEmModel(msc2);
0188       if (nullptr != aRegion) {
0189         msc->AddEmModel(-1, msc3, aRegion);
0190       }
0191       if (nullptr != bRegion) {
0192         msc->AddEmModel(-1, msc3, bRegion);
0193       }
0194       ph->RegisterProcess(msc, particle);
0195     }
0196
0197     // single scattering
0198     G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0199     G4CoulombScattering* ss = new G4CoulombScattering();
0200     ss->SetEmModel(ssm);
0201     ss->SetMinKinEnergy(highEnergyLimit);
0202     ssm->SetLowEnergyLimit(highEnergyLimit);
0203     ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0204
0205     ph->RegisterProcess(new G4eIonisation(), particle);
0206     ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0207     ph->RegisterProcess(ss, particle);
0208
0209     // e+
0210     particle = G4Positron::Positron();
0211
0212     msc1 = new G4UrbanMscModel();
0213     msc2 = new G4WentzelVIModel();
0214     msc1->SetHighEnergyLimit(highEnergyLimit);
0215     msc2->SetLowEnergyLimit(highEnergyLimit);
0216
0217     // e-/e+ msc for HCAL and HGCAL using the Urban model
0218     if (nullptr != aRegion || nullptr != bRegion) {
0219       msc3 = new G4UrbanMscModel();
0220       msc3->SetHighEnergyLimit(highEnergyLimit);
0221       msc3->SetRangeFactor(fRangeFactor);
0222       msc3->SetGeomFactor(fGeomFactor);
0223       msc3->SetSafetyFactor(fSafetyFactor);
0224       msc3->SetLambdaLimit(fLambdaLimit);
0225       msc3->SetStepLimitType(fStepLimitType);
0226       msc3->SetLocked(true);
0227     }
0228
0229     if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0230       G4ProcessManager* procManager = particle->GetProcessManager();
0231       // Remove default G4Transportation and replace with G4TransportationWithMsc.
0232       G4VProcess* removed = procManager->RemoveProcess(0);
0233       if (removed->GetProcessName() != "Transportation") {
0234         G4Exception("CMSEmStandardPhysics::ConstructProcess",
0235                     "em0050",
0236                     FatalException,
0237                     "replaced process is not G4Transportation!");
0238       }
0239       G4TransportationWithMsc* transportWithMsc =
0240           new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0241       if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0242         transportWithMsc->SetMultipleSteps(true);
0243       }
0244       transportWithMsc->AddMscModel(msc1);
0245       transportWithMsc->AddMscModel(msc2);
0246       if (nullptr != aRegion) {
0247         transportWithMsc->AddMscModel(msc3, -1, aRegion);
0248       }
0249       if (nullptr != bRegion) {
0250         transportWithMsc->AddMscModel(msc3, -1, bRegion);
0251       }
0252       procManager->AddProcess(transportWithMsc, -1, 0, 0);
0253     } else {
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(new G4eIonisation(), particle);
0276     ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0277     ph->RegisterProcess(new G4eplusAnnihilation(), particle);
0278     ph->RegisterProcess(ss, particle);
0279
0280   } else {
0281     // G4HepEm is Active
0282     if (verboseLevel > 0) {
0283       edm::LogVerbatim("PhysicsList") << "G4HepEm is active, registering G4HepEmTrackingManager";
0284     }
0285     auto* hepEmTM = new G4HepEmTrackingManager(verboseLevel);
0286
0287     // configuration
0288     G4HepEmConfig* config = hepEmTM->GetConfig();
0289     // First set global configuration parameters:
0290     // ------------------------------------------
0291     // The default MSC `RangeFactor`, `SafetyFactor`, `StepLimitType` parameters
0292     // as well as  `SetApplyCuts` and `SetLowestElectronEnergy` are taken from
0293     // `G4EmParameters` so we set only the step function parameters here.
0294
0295     config->SetEnergyLossStepLimitFunctionParameters(0.8, 1.0 * CLHEP::mm);
0296
0297     // Then set special configuration for some regions:
0298     // ------------------------------------------------
0299     if (nullptr != aRegion) {
0300       // HCal region
0301       const G4String& rname = aRegion->GetName();
0302       config->SetMinimalMSCStepLimit(fStepLimitType == fMinimal, rname);
0303       config->SetMSCRangeFactor(fRangeFactor, rname);
0304       config->SetMSCSafetyFactor(fSafetyFactor, rname);
0305     }
0306
0307     if (nullptr != bRegion) {
0308       // HGCal region
0309       const G4String& rname = bRegion->GetName();
0310       config->SetMinimalMSCStepLimit(fStepLimitType == fMinimal, rname);
0311       config->SetMSCRangeFactor(fRangeFactor, rname);
0312       config->SetMSCSafetyFactor(fSafetyFactor, rname);
0313
0314       config->SetWoodcockTrackingRegion(rname);
0315       config->SetWDTEnergyLimit(0.5 * CLHEP::MeV);
0316     }
0317
0318     G4Electron::Electron()->SetTrackingManager(hepEmTM);
0319     G4Positron::Positron()->SetTrackingManager(hepEmTM);
0320     G4Gamma::Gamma()->SetTrackingManager(hepEmTM);
0321   }
0322
0323   // generic ion
0324   G4ParticleDefinition* particle = G4GenericIon::GenericIon();
0325   G4ionIonisation* ionIoni = new G4ionIonisation();
0326   ph->RegisterProcess(hmsc, particle);
0327   ph->RegisterProcess(ionIoni, particle);
0328
0329   // muons, hadrons ions
0330   G4EmBuilder::ConstructCharged(hmsc, pnuc);
0331 }