Application/src/ParametrisedEMPhysics.cc

0001 //
0002 // Joanna Weng 08.2005
0003 // Physics process for Gflash parameterisation
0004 // modified by Soon Yung Jun, Dongwook Jang
0005 // V.Ivanchenko rename the class, cleanup, and move
0006 //              to SimG4Core/Application - 2012/08/14
0007
0008 #include "SimG4Core/Application/interface/ParametrisedEMPhysics.h"
0009 #include "SimG4Core/Application/interface/GFlashEMShowerModel.h"
0010 #include "SimG4Core/Application/interface/GFlashHadronShowerModel.h"
0011 #include "SimG4Core/Application/interface/LowEnergyFastSimModel.h"
0012 #include "SimG4Core/Application/interface/ElectronLimiter.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "G4FastSimulationManagerProcess.hh"
0016 #include "G4ProcessManager.hh"
0017 #include "G4HadronicParameters.hh"
0018 #include "G4EmBuilder.hh"
0019
0020 #include "G4RegionStore.hh"
0021 #include "G4Electron.hh"
0022 #include "G4Positron.hh"
0023 #include "G4MuonMinus.hh"
0024 #include "G4MuonPlus.hh"
0025 #include "G4PionMinus.hh"
0026 #include "G4PionPlus.hh"
0027 #include "G4KaonMinus.hh"
0028 #include "G4KaonPlus.hh"
0029 #include "G4Proton.hh"
0030 #include "G4AntiProton.hh"
0031
0032 #include "G4EmParameters.hh"
0033 #include "G4LossTableManager.hh"
0034 #include "G4PhysicsListHelper.hh"
0035 #include "G4ProcessManager.hh"
0036 #include <CLHEP/Units/SystemOfUnits.h>
0037 #include "G4Transportation.hh"
0038 #include "G4UAtomicDeexcitation.hh"
0039 #include <memory>
0040
0041 #include <string>
0042 #include <vector>
0043
0044 const G4int NREG = 7;
0045 const G4String rname[NREG] = {"EcalRegion",
0046                               "HcalRegion",
0047                               "HGcalRegion",
0048                               "MuonIron",
0049                               "PreshowerRegion",
0050                               "CastorRegion",
0051                               "DefaultRegionForTheWorld"};
0052
0053 struct ParametrisedEMPhysics::TLSmod {
0054   std::unique_ptr<GFlashEMShowerModel> theEcalEMShowerModel;
0055   std::unique_ptr<LowEnergyFastSimModel> theLowEnergyFastSimModel;
0056   std::unique_ptr<GFlashEMShowerModel> theHcalEMShowerModel;
0057   std::unique_ptr<GFlashHadronShowerModel> theEcalHadShowerModel;
0058   std::unique_ptr<GFlashHadronShowerModel> theHcalHadShowerModel;
0059   std::unique_ptr<G4FastSimulationManagerProcess> theFastSimulationManagerProcess;
0060 };
0061
0062 G4ThreadLocal ParametrisedEMPhysics::TLSmod* ParametrisedEMPhysics::m_tpmod = nullptr;
0063
0064 ParametrisedEMPhysics::ParametrisedEMPhysics(const std::string& name, const edm::ParameterSet& p)
0065     : G4VPhysicsConstructor(name), theParSet(p) {
0066   // EM parameters common for any EM physics configuration
0067   G4EmParameters* param = G4EmParameters::Instance();
0068   G4int verb = theParSet.getUntrackedParameter<int>("Verbosity", 0);
0069   param->SetVerbose(verb);
0070
0071   G4double bremth = theParSet.getParameter<double>("G4BremsstrahlungThreshold") * CLHEP::GeV;
0072   param->SetBremsstrahlungTh(bremth);
0073   G4double mubrth = theParSet.getParameter<double>("G4MuonBremsstrahlungThreshold") * CLHEP::GeV;
0074   param->SetMuHadBremsstrahlungTh(mubrth);
0075
0076   bool genp = theParSet.getParameter<bool>("G4GammaGeneralProcess");
0077   param->SetGeneralProcessActive(genp);
0078
0079   bool pe = p.getParameter<bool>("PhotoeffectBelowKShell");
0080   param->SetPhotoeffectBelowKShell(pe);
0081
0082   int type = p.getParameter<int>("G4TransportWithMSC");
0083   G4TransportationWithMscType trtype = G4TransportationWithMscType::fDisabled;
0084   if (type == 1) {
0085     trtype = G4TransportationWithMscType::fEnabled;
0086   } else if (type == 2) {
0087     trtype = G4TransportationWithMscType::fMultipleSteps;
0088   }
0089   param->SetTransportationWithMsc(trtype);
0090
0091   bool mudat = theParSet.getParameter<bool>("ReadMuonData");
0092   param->SetRetrieveMuDataFromFile(mudat);
0093
0094   bool fluo = theParSet.getParameter<bool>("FlagFluo");
0095   param->SetFluo(fluo);
0096
0097   bool modifyT = theParSet.getParameter<bool>("ModifyTransportation");
0098   double th1 = theParSet.getUntrackedParameter<double>("ThresholdWarningEnergy") * CLHEP::MeV;
0099   double th2 = theParSet.getUntrackedParameter<double>("ThresholdImportantEnergy") * CLHEP::MeV;
0100   int nt = theParSet.getUntrackedParameter<int>("ThresholdTrials");
0101
0102   edm::LogVerbatim("SimG4CoreApplication")
0103       << "### ParametrisedEMPhysics parameters:"
0104       << "\n verbosity= " << verb << "\n fluoFlag: " << fluo << "\n fluo below K-shell: " << pe
0105       << "\n transportation with msc: " << type << "\n modifyTransport: " << modifyT << "  Ntrials= " << nt
0106       << "\n ThWarning(MeV)= " << th1 / CLHEP::MeV << "\n ThException(MeV)= " << th2 / CLHEP::MeV
0107       << "\n read muon data: " << mudat << "\n bremsstrahlung threshold Eth(GeV)= " << bremth / CLHEP::GeV;
0108
0109   // Russian roulette and tracking cut for e+-
0110   double energyLim = theParSet.getParameter<double>("RusRoElectronEnergyLimit") * CLHEP::MeV;
0111   if (energyLim > 0.0) {
0112     G4double rrfact[NREG] = {1.0};
0113
0114     rrfact[0] = theParSet.getParameter<double>("RusRoEcalElectron");
0115     rrfact[1] = theParSet.getParameter<double>("RusRoHcalElectron");
0116     rrfact[2] = theParSet.getParameter<double>("RusRoMuonIronElectron");
0117     rrfact[3] = theParSet.getParameter<double>("RusRoPreShowerElectron");
0118     rrfact[4] = theParSet.getParameter<double>("RusRoCastorElectron");
0119     rrfact[5] = theParSet.getParameter<double>("RusRoWorldElectron");
0120     for (int i = 0; i < NREG; ++i) {
0121       if (rrfact[i] < 1.0) {
0122         param->ActivateSecondaryBiasing("eIoni", rname[i], rrfact[i], energyLim);
0123         param->ActivateSecondaryBiasing("hIoni", rname[i], rrfact[i], energyLim);
0124         edm::LogVerbatim("SimG4CoreApplication")
0125             << "ParametrisedEMPhysics: Russian Roulette"
0126             << " for e- Prob= " << rrfact[i] << " Elimit(MeV)= " << energyLim / CLHEP::MeV << " inside " << rname[i];
0127       }
0128     }
0129   }
0130 }
0131
0132 ParametrisedEMPhysics::~ParametrisedEMPhysics() {
0133   delete m_tpmod;
0134   m_tpmod = nullptr;
0135 }
0136
0137 void ParametrisedEMPhysics::ConstructParticle() {
0138   // minimal set of particles for EM physics
0139   G4EmBuilder::ConstructMinimalEmSet();
0140 }
0141
0142 void ParametrisedEMPhysics::ConstructProcess() {
0143   edm::LogVerbatim("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess() started";
0144
0145   // GFlash part
0146   bool gem = theParSet.getParameter<bool>("GflashEcal");
0147   bool lowEnergyGem = theParSet.getParameter<bool>("LowEnergyGflashEcal");
0148   bool ghad = theParSet.getParameter<bool>("GflashHcal");
0149   bool gemHad = theParSet.getParameter<bool>("GflashEcalHad");
0150   bool ghadHad = theParSet.getParameter<bool>("GflashHcalHad");
0151
0152   if (gem || ghad || lowEnergyGem || gemHad || ghadHad) {
0153     if (!m_tpmod) {
0154       m_tpmod = new TLSmod;
0155     }
0156     edm::LogVerbatim("SimG4CoreApplication")
0157         << "ParametrisedEMPhysics: GFlash Construct for e+-: " << gem << "  " << ghad << " " << lowEnergyGem
0158         << " for hadrons: " << gemHad << "  " << ghadHad;
0159
0160     m_tpmod->theFastSimulationManagerProcess = std::make_unique<G4FastSimulationManagerProcess>();
0161
0162     if (gem || ghad) {
0163       G4Electron::Electron()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0164       G4Positron::Positron()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0165     } else if (lowEnergyGem) {
0166       G4Electron::Electron()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0167       G4Positron::Positron()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0168     }
0169
0170     if (gemHad || ghadHad) {
0171       G4Proton::Proton()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0172       G4AntiProton::AntiProton()->GetProcessManager()->AddDiscreteProcess(
0173           m_tpmod->theFastSimulationManagerProcess.get());
0174       G4PionPlus::PionPlus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0175       G4PionMinus::PionMinus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0176       G4KaonPlus::KaonPlus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0177       G4KaonMinus::KaonMinus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
0178     }
0179
0180     if (gem || gemHad || lowEnergyGem) {
0181       G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("EcalRegion", false);
0182
0183       if (!aRegion) {
0184         edm::LogWarning("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess: "
0185                                                 << "EcalRegion is not defined, GFlash will not be enabled for ECAL!";
0186
0187       } else {
0188         if (gem) {
0189           //Electromagnetic Shower Model for ECAL
0190           m_tpmod->theEcalEMShowerModel =
0191               std::make_unique<GFlashEMShowerModel>("GflashEcalEMShowerModel", aRegion, theParSet);
0192         } else if (lowEnergyGem) {
0193           //Low energy electromagnetic Shower Model for ECAL
0194           m_tpmod->theLowEnergyFastSimModel =
0195               std::make_unique<LowEnergyFastSimModel>("LowEnergyFastSimModel", aRegion, theParSet);
0196         }
0197
0198         if (gemHad) {
0199           //Electromagnetic Shower Model for ECAL
0200           m_tpmod->theEcalHadShowerModel =
0201               std::make_unique<GFlashHadronShowerModel>("GflashEcalHadShowerModel", aRegion, theParSet);
0202         }
0203       }
0204     }
0205     if (ghad || ghadHad) {
0206       G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
0207       if (!aRegion) {
0208         edm::LogWarning("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess: "
0209                                                 << "HcalRegion is not defined, GFlash will not be enabled for HCAL!";
0210
0211       } else {
0212         if (ghad) {
0213           //Electromagnetic Shower Model for HCAL
0214           m_tpmod->theHcalEMShowerModel =
0215               std::make_unique<GFlashEMShowerModel>("GflashHcalEMShowerModel", aRegion, theParSet);
0216         }
0217         if (ghadHad) {
0218           //Electromagnetic Shower Model for ECAL
0219           m_tpmod->theHcalHadShowerModel =
0220               std::make_unique<GFlashHadronShowerModel>("GflashHcalHadShowerModel", aRegion, theParSet);
0221         }
0222       }
0223     }
0224   }
0225
0226   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0227
0228   // Step limiters for e+-
0229   bool eLimiter = theParSet.getParameter<bool>("ElectronStepLimit");
0230   bool rLimiter = theParSet.getParameter<bool>("ElectronRangeTest");
0231   bool pLimiter = theParSet.getParameter<bool>("PositronStepLimit");
0232   // Step limiters for hadrons
0233   bool pTCut = theParSet.getParameter<bool>("ProtonRegionLimit");
0234   bool piTCut = theParSet.getParameter<bool>("PionRegionLimit");
0235
0236   std::vector<std::string> regnames = theParSet.getParameter<std::vector<std::string> >("LimitsPerRegion");
0237   std::vector<double> limitsE = theParSet.getParameter<std::vector<double> >("EnergyLimitsE");
0238   std::vector<double> limitsH = theParSet.getParameter<std::vector<double> >("EnergyLimitsH");
0239   std::vector<double> facE = theParSet.getParameter<std::vector<double> >("EnergyFactorsE");
0240   std::vector<double> rmsE = theParSet.getParameter<std::vector<double> >("EnergyRMSE");
0241   int nlimits = regnames.size();
0242   int nlimitsH = 0;
0243   std::vector<const G4Region*> reg;
0244   std::vector<G4double> rlimE;
0245   std::vector<G4double> rlimH;
0246   std::vector<G4double> factE;
0247   std::vector<G4double> rmsvE;
0248   if (0 < nlimits) {
0249     G4RegionStore* store = G4RegionStore::GetInstance();
0250     for (int i = 0; i < nlimits; ++i) {
0251       // apply limiter for whole CMS
0252       if (regnames[i] == "all") {
0253         reg.clear();
0254         rlimE.clear();
0255         rlimH.clear();
0256         factE.clear();
0257         rmsvE.clear();
0258         reg.emplace_back(nullptr);
0259         rlimE.emplace_back(limitsE[i] * CLHEP::MeV);
0260         rlimH.emplace_back(limitsH[i] * CLHEP::MeV);
0261         factE.emplace_back(facE[i]);
0262         rmsvE.emplace_back(rmsE[i]);
0263         nlimitsH = (limitsH[i] > 0) ? 1 : 0;
0264         break;
0265       }
0266       const G4Region* r = store->GetRegion(regnames[i], false);
0267       // apply for concrete G4Region
0268       if (r && (limitsE[i] > 0.0 || limitsH[i] > 0.0)) {
0269         reg.emplace_back(r);
0270         rlimE.emplace_back(limitsE[i] * CLHEP::MeV);
0271         rlimH.emplace_back(limitsH[i] * CLHEP::MeV);
0272         factE.emplace_back(facE[i]);
0273         rmsvE.emplace_back(rmsE[i]);
0274         if (limitsH[i] > 0) {
0275           ++nlimitsH;
0276         }
0277       }
0278     }
0279     nlimits = reg.size();
0280   }
0281
0282   if (eLimiter || rLimiter || 0 < nlimits) {
0283     ElectronLimiter* elim = new ElectronLimiter(theParSet, G4Electron::Electron());
0284     elim->SetRangeCheckFlag(rLimiter);
0285     elim->SetFieldCheckFlag(eLimiter);
0286     elim->SetTrackingCutPerRegion(reg, rlimE, factE, rmsvE);
0287     ph->RegisterProcess(elim, G4Electron::Electron());
0288   }
0289
0290   if (pLimiter || 0 < nlimits) {
0291     ElectronLimiter* plim = new ElectronLimiter(theParSet, G4Positron::Positron());
0292     plim->SetFieldCheckFlag(pLimiter);
0293     plim->SetTrackingCutPerRegion(reg, rlimE, factE, rmsvE);
0294     ph->RegisterProcess(plim, G4Positron::Positron());
0295   }
0296   if (0 < nlimits && 0 < nlimitsH) {
0297     if (pTCut) {
0298       ElectronLimiter* plim = new ElectronLimiter(theParSet, G4Proton::Proton());
0299       plim->SetFieldCheckFlag(pLimiter);
0300       plim->SetTrackingCutPerRegion(reg, rlimH, factE, rmsvE);
0301       ph->RegisterProcess(plim, G4Proton::Proton());
0302     }
0303     if (piTCut) {
0304       ElectronLimiter* plim = new ElectronLimiter(theParSet, G4PionPlus::PionPlus());
0305       plim->SetFieldCheckFlag(pLimiter);
0306       plim->SetTrackingCutPerRegion(reg, rlimH, factE, rmsvE);
0307       ph->RegisterProcess(plim, G4PionPlus::PionPlus());
0308       plim = new ElectronLimiter(theParSet, G4PionMinus::PionMinus());
0309       plim->SetFieldCheckFlag(pLimiter);
0310       plim->SetTrackingCutPerRegion(reg, rlimH, factE, rmsvE);
0311       ph->RegisterProcess(plim, G4PionMinus::PionMinus());
0312     }
0313   }
0314   // enable fluorescence
0315   bool fluo = theParSet.getParameter<bool>("FlagFluo");
0316   if (fluo && !G4LossTableManager::Instance()->AtomDeexcitation()) {
0317     G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
0318     G4LossTableManager::Instance()->SetAtomDeexcitation(de);
0319   }
0320   // change parameters of transportation
0321   bool modifyT = theParSet.getParameter<bool>("ModifyTransportation");
0322   if (modifyT) {
0323     double th1 = theParSet.getUntrackedParameter<double>("ThresholdWarningEnergy") * CLHEP::MeV;
0324     double th2 = theParSet.getUntrackedParameter<double>("ThresholdImportantEnergy") * CLHEP::MeV;
0325     int nt = theParSet.getUntrackedParameter<int>("ThresholdTrials");
0326     ModifyTransportation(G4Electron::Electron(), nt, th1, th2);
0327   }
0328   edm::LogVerbatim("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess() is done";
0329 }
0330
0331 void ParametrisedEMPhysics::ModifyTransportation(const G4ParticleDefinition* part, int ntry, double th1, double th2) {
0332   G4ProcessManager* man = part->GetProcessManager();
0333   G4Transportation* trans = (G4Transportation*)((*(man->GetProcessList()))[0]);
0334   if (trans) {
0335     trans->SetThresholdWarningEnergy(th1);
0336     trans->SetThresholdImportantEnergy(th2);
0337     trans->SetThresholdTrials(ntry);
0338     edm::LogVerbatim("SimG4CoreApplication")
0339         << "ParametrisedEMPhysics: printout level changed for " << part->GetParticleName();
0340   }
0341 }