Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-11 22:29:58

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