CMS 3D CMS Logo
ParametrisedEMPhysics.cc
Go to the documentation of this file.
1 //
2 // Joanna Weng 08.2005
3 // Physics process for Gflash parameterisation
4 // modified by Soon Yung Jun, Dongwook Jang
5 // V.Ivanchenko rename the class, cleanup, and move
6 // to SimG4Core/Application - 2012/08/14
7 
8 #include "SimG4Core/Application/interface/ParametrisedEMPhysics.h"
9 #include "SimG4Core/Application/interface/GFlashEMShowerModel.h"
10 #include "SimG4Core/Application/interface/GFlashHadronShowerModel.h"
11 #include "SimG4Core/Application/interface/LowEnergyFastSimModel.h"
12 #include "SimG4Core/Application/interface/ElectronLimiter.h"
13 #include "FWCore/MessageLogger/interface/MessageLogger.h"
14 
15 #include "G4FastSimulationManagerProcess.hh"
16 #include "G4ProcessManager.hh"
17 
18 #include "G4LeptonConstructor.hh"
19 #include "G4MesonConstructor.hh"
20 #include "G4BaryonConstructor.hh"
21 #include "G4ShortLivedConstructor.hh"
22 #include "G4IonConstructor.hh"
23 #include "G4RegionStore.hh"
24 #include "G4Electron.hh"
25 #include "G4Positron.hh"
26 #include "G4MuonMinus.hh"
27 #include "G4MuonPlus.hh"
28 #include "G4PionMinus.hh"
29 #include "G4PionPlus.hh"
30 #include "G4KaonMinus.hh"
31 #include "G4KaonPlus.hh"
32 #include "G4Proton.hh"
33 #include "G4AntiProton.hh"
34 
35 #include "G4EmParameters.hh"
36 #include "G4PhysicsListHelper.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4UAtomicDeexcitation.hh"
39 #include "G4LossTableManager.hh"
40 #include "G4ProcessManager.hh"
41 #include "G4Transportation.hh"
42 #include <string>
43 #include <vector>
44 
45 const G4int NREG = 7;
46 const G4String rname[NREG] = {"EcalRegion",
47  "HcalRegion",
48  "HGcalRegion",
49  "MuonIron",
50  "PreshowerRegion",
51  "CastorRegion",
52  "DefaultRegionForTheWorld"};
53 
54 struct ParametrisedEMPhysics::TLSmod {
55  std::unique_ptr<GFlashEMShowerModel> theEcalEMShowerModel;
56  std::unique_ptr<LowEnergyFastSimModel> theLowEnergyFastSimModel;
57  std::unique_ptr<GFlashEMShowerModel> theHcalEMShowerModel;
58  std::unique_ptr<GFlashHadronShowerModel> theEcalHadShowerModel;
59  std::unique_ptr<GFlashHadronShowerModel> theHcalHadShowerModel;
60  std::unique_ptr<G4FastSimulationManagerProcess> theFastSimulationManagerProcess;
61 };
62 
63 G4ThreadLocal ParametrisedEMPhysics::TLSmod* ParametrisedEMPhysics::m_tpmod = nullptr;
64 
65 ParametrisedEMPhysics::ParametrisedEMPhysics(const std::string& name, const edm::ParameterSet& p)
66  : G4VPhysicsConstructor(name), theParSet(p) {
67  // bremsstrahlung threshold and EM verbosity
68  G4EmParameters* param = G4EmParameters::Instance();
69  G4int verb = theParSet.getUntrackedParameter<int>("Verbosity", 0);
70  param->SetVerbose(verb);
71 
72  G4double bremth = theParSet.getParameter<double>("G4BremsstrahlungThreshold") * GeV;
73  param->SetBremsstrahlungTh(bremth);
74 
75  bool fluo = theParSet.getParameter<bool>("FlagFluo");
76  param->SetFluo(fluo);
77 
78  bool modifyT = theParSet.getParameter<bool>("ModifyTransportation");
79  double th1 = theParSet.getUntrackedParameter<double>("ThresholdWarningEnergy") * MeV;
80  double th2 = theParSet.getUntrackedParameter<double>("ThresholdImportantEnergy") * MeV;
81  int nt = theParSet.getUntrackedParameter<int>("ThresholdTrials");
82 
83  edm::LogVerbatim("SimG4CoreApplication")
84  << "ParametrisedEMPhysics::ConstructProcess: bremsstrahlung threshold Eth= " << bremth / GeV << " GeV"
85  << "\n verbosity= " << verb << " fluoFlag: " << fluo
86  << " modifyTransport: " << modifyT << " Ntrials= " << nt
87  << "\n ThWarning(MeV)= " << th1 << " ThException(MeV)= " << th2;
88 
89  // Russian roulette and tracking cut for e+-
90  double energyLim = theParSet.getParameter<double>("RusRoElectronEnergyLimit") * MeV;
91  if (energyLim > 0.0) {
92  G4double rrfact[NREG] = {1.0};
93 
94  rrfact[0] = theParSet.getParameter<double>("RusRoEcalElectron");
95  rrfact[1] = theParSet.getParameter<double>("RusRoHcalElectron");
96  rrfact[2] = theParSet.getParameter<double>("RusRoMuonIronElectron");
97  rrfact[3] = theParSet.getParameter<double>("RusRoPreShowerElectron");
98  rrfact[4] = theParSet.getParameter<double>("RusRoCastorElectron");
99  rrfact[5] = theParSet.getParameter<double>("RusRoWorldElectron");
100  for (int i = 0; i < NREG; ++i) {
101  if (rrfact[i] < 1.0) {
102  param->ActivateSecondaryBiasing("eIoni", rname[i], rrfact[i], energyLim);
103  param->ActivateSecondaryBiasing("hIoni", rname[i], rrfact[i], energyLim);
104  edm::LogVerbatim("SimG4CoreApplication")
105  << "ParametrisedEMPhysics: Russian Roulette"
106  << " for e- Prob= " << rrfact[i] << " Elimit(MeV)= " << energyLim / CLHEP::MeV << " inside " << rname[i];
107  }
108  }
109  }
110 }
111 
112 ParametrisedEMPhysics::~ParametrisedEMPhysics() {
113  if (m_tpmod) {
114  delete m_tpmod;
115  m_tpmod = nullptr;
116  }
117 }
118 
119 void ParametrisedEMPhysics::ConstructParticle() {
120  G4LeptonConstructor pLeptonConstructor;
121  pLeptonConstructor.ConstructParticle();
122 
123  G4MesonConstructor pMesonConstructor;
124  pMesonConstructor.ConstructParticle();
125 
126  G4BaryonConstructor pBaryonConstructor;
127  pBaryonConstructor.ConstructParticle();
128 
129  G4ShortLivedConstructor pShortLivedConstructor;
130  pShortLivedConstructor.ConstructParticle();
131 
132  G4IonConstructor pConstructor;
133  pConstructor.ConstructParticle();
134 }
135 
136 void ParametrisedEMPhysics::ConstructProcess() {
137  edm::LogVerbatim("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess() started";
138 
139  // GFlash part
140  bool gem = theParSet.getParameter<bool>("GflashEcal");
141  bool lowEnergyGem = theParSet.getParameter<bool>("LowEnergyGflashEcal");
142  bool ghad = theParSet.getParameter<bool>("GflashHcal");
143  bool gemHad = theParSet.getParameter<bool>("GflashEcalHad");
144  bool ghadHad = theParSet.getParameter<bool>("GflashHcalHad");
145 
146  if (gem || ghad || lowEnergyGem || gemHad || ghadHad) {
147  if (!m_tpmod) {
148  m_tpmod = new TLSmod;
149  }
150  edm::LogVerbatim("SimG4CoreApplication")
151  << "ParametrisedEMPhysics: GFlash Construct for e+-: " << gem << " " << ghad << " " << lowEnergyGem
152  << " for hadrons: " << gemHad << " " << ghadHad;
153 
154  m_tpmod->theFastSimulationManagerProcess.reset(new G4FastSimulationManagerProcess());
155 
156  if (gem || ghad) {
157  G4Electron::Electron()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
158  G4Positron::Positron()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
159  } else if (lowEnergyGem) {
160  G4Electron::Electron()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
161  }
162 
163  if (gemHad || ghadHad) {
164  G4Proton::Proton()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
165  G4AntiProton::AntiProton()->GetProcessManager()->AddDiscreteProcess(
166  m_tpmod->theFastSimulationManagerProcess.get());
167  G4PionPlus::PionPlus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
168  G4PionMinus::PionMinus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
169  G4KaonPlus::KaonPlus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
170  G4KaonMinus::KaonMinus()->GetProcessManager()->AddDiscreteProcess(m_tpmod->theFastSimulationManagerProcess.get());
171  }
172 
173  if (gem || gemHad || lowEnergyGem) {
174  G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("EcalRegion", false);
175 
176  if (!aRegion) {
177  edm::LogWarning("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess: "
178  << "EcalRegion is not defined, GFlash will not be enabled for ECAL!";
179 
180  } else {
181  if (gem) {
182  //Electromagnetic Shower Model for ECAL
183  m_tpmod->theEcalEMShowerModel.reset(new GFlashEMShowerModel("GflashEcalEMShowerModel", aRegion, theParSet));
184  } else if (lowEnergyGem) {
185  //Low energy electromagnetic Shower Model for ECAL
186  m_tpmod->theLowEnergyFastSimModel.reset(
187  new LowEnergyFastSimModel("LowEnergyFastSimModel", aRegion, theParSet));
188  }
189 
190  if (gemHad) {
191  //Electromagnetic Shower Model for ECAL
192  m_tpmod->theEcalHadShowerModel.reset(
193  new GFlashHadronShowerModel("GflashEcalHadShowerModel", aRegion, theParSet));
194  }
195  }
196  }
197  if (ghad || ghadHad) {
198  G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
199  if (!aRegion) {
200  edm::LogWarning("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess: "
201  << "HcalRegion is not defined, GFlash will not be enabled for HCAL!";
202 
203  } else {
204  if (ghad) {
205  //Electromagnetic Shower Model for HCAL
206  m_tpmod->theHcalEMShowerModel.reset(new GFlashEMShowerModel("GflashHcalEMShowerModel", aRegion, theParSet));
207  }
208  if (ghadHad) {
209  //Electromagnetic Shower Model for ECAL
210  m_tpmod->theHcalHadShowerModel.reset(
211  new GFlashHadronShowerModel("GflashHcalHadShowerModel", aRegion, theParSet));
212  }
213  }
214  }
215  }
216 
217  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
218 
219  // Step limiters for e+-
220  bool eLimiter = theParSet.getParameter<bool>("ElectronStepLimit");
221  bool rLimiter = theParSet.getParameter<bool>("ElectronRangeTest");
222  bool pLimiter = theParSet.getParameter<bool>("PositronStepLimit");
223  // Step limiters for hadrons
224  bool pTCut = theParSet.getParameter<bool>("ProtonRegionLimit");
225  bool piTCut = theParSet.getParameter<bool>("PionRegionLimit");
226 
227  std::vector<std::string> regnames = theParSet.getParameter<std::vector<std::string> >("LimitsPerRegion");
228  std::vector<double> limitsE = theParSet.getParameter<std::vector<double> >("EnergyLimitsE");
229  std::vector<double> limitsH = theParSet.getParameter<std::vector<double> >("EnergyLimitsH");
230  std::vector<double> facE = theParSet.getParameter<std::vector<double> >("EnergyFactorsE");
231  std::vector<double> rmsE = theParSet.getParameter<std::vector<double> >("EnergyRMSE");
232  int nlimits = regnames.size();
233  int nlimitsH = 0;
234  std::vector<const G4Region*> reg;
235  std::vector<G4double> rlimE;
236  std::vector<G4double> rlimH;
237  std::vector<G4double> factE;
238  std::vector<G4double> rmsvE;
239  if (0 < nlimits) {
240  G4RegionStore* store = G4RegionStore::GetInstance();
241  for (int i = 0; i < nlimits; ++i) {
242  // apply limiter for whole CMS
243  if (regnames[i] == "all") {
244  reg.clear();
245  rlimE.clear();
246  rlimH.clear();
247  factE.clear();
248  rmsvE.clear();
249  reg.emplace_back(nullptr);
250  rlimE.emplace_back(limitsE[i] * CLHEP::MeV);
251  rlimH.emplace_back(limitsH[i] * CLHEP::MeV);
252  factE.emplace_back(facE[i]);
253  rmsvE.emplace_back(rmsE[i]);
254  nlimitsH = (limitsH[i] > 0) ? 1 : 0;
255  break;
256  }
257  const G4Region* r = store->GetRegion(regnames[i], false);
258  // apply for concrete G4Region
259  if (r && (limitsE[i] > 0.0 || limitsH[i] > 0.0)) {
260  reg.emplace_back(r);
261  rlimE.emplace_back(limitsE[i] * CLHEP::MeV);
262  rlimH.emplace_back(limitsH[i] * CLHEP::MeV);
263  factE.emplace_back(facE[i]);
264  rmsvE.emplace_back(rmsE[i]);
265  if (limitsH[i] > 0) {
266  ++nlimitsH;
267  }
268  }
269  }
270  nlimits = reg.size();
271  }
272 
273  if (eLimiter || rLimiter || 0 < nlimits) {
274  ElectronLimiter* elim = new ElectronLimiter(theParSet, G4Electron::Electron());
275  elim->SetRangeCheckFlag(rLimiter);
276  elim->SetFieldCheckFlag(eLimiter);
277  elim->SetTrackingCutPerRegion(reg, rlimE, factE, rmsvE);
278  ph->RegisterProcess(elim, G4Electron::Electron());
279  }
280 
281  if (pLimiter || 0 < nlimits) {
282  ElectronLimiter* plim = new ElectronLimiter(theParSet, G4Positron::Positron());
283  plim->SetFieldCheckFlag(pLimiter);
284  plim->SetTrackingCutPerRegion(reg, rlimE, factE, rmsvE);
285  ph->RegisterProcess(plim, G4Positron::Positron());
286  }
287  if (0 < nlimits && 0 < nlimitsH) {
288  if (pTCut) {
289  ElectronLimiter* plim = new ElectronLimiter(theParSet, G4Proton::Proton());
290  plim->SetFieldCheckFlag(pLimiter);
291  plim->SetTrackingCutPerRegion(reg, rlimH, factE, rmsvE);
292  ph->RegisterProcess(plim, G4Proton::Proton());
293  }
294  if (piTCut) {
295  ElectronLimiter* plim = new ElectronLimiter(theParSet, G4PionPlus::PionPlus());
296  plim->SetFieldCheckFlag(pLimiter);
297  plim->SetTrackingCutPerRegion(reg, rlimH, factE, rmsvE);
298  ph->RegisterProcess(plim, G4PionPlus::PionPlus());
299  plim = new ElectronLimiter(theParSet, G4PionMinus::PionMinus());
300  plim->SetFieldCheckFlag(pLimiter);
301  plim->SetTrackingCutPerRegion(reg, rlimH, factE, rmsvE);
302  ph->RegisterProcess(plim, G4PionMinus::PionMinus());
303  }
304  }
305  // enable fluorescence
306  bool fluo = theParSet.getParameter<bool>("FlagFluo");
307  if (fluo && !G4LossTableManager::Instance()->AtomDeexcitation()) {
308  G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
309  G4LossTableManager::Instance()->SetAtomDeexcitation(de);
310  }
311  // change parameters of transportation
312  bool modifyT = theParSet.getParameter<bool>("ModifyTransportation");
313  if (modifyT) {
314  double th1 = theParSet.getUntrackedParameter<double>("ThresholdWarningEnergy") * MeV;
315  double th2 = theParSet.getUntrackedParameter<double>("ThresholdImportantEnergy") * MeV;
316  int nt = theParSet.getUntrackedParameter<int>("ThresholdTrials");
317  ModifyTransportation(G4Electron::Electron(), nt, th1, th2);
318  ModifyTransportation(G4Positron::Positron(), nt, th1, th2);
319  ModifyTransportation(G4Proton::Proton(), nt, th1, th2);
320  }
321  edm::LogVerbatim("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess() is done";
322 }
323 
324 void ParametrisedEMPhysics::ModifyTransportation(const G4ParticleDefinition* part, int ntry, double th1, double th2) {
325  G4ProcessManager* man = part->GetProcessManager();
326  G4Transportation* trans = (G4Transportation*)((*(man->GetProcessList()))[0]);
327  if (trans) {
328  trans->SetThresholdWarningEnergy(th1);
329  trans->SetThresholdImportantEnergy(th2);
330  trans->SetThresholdTrials(ntry);
331  edm::LogVerbatim("SimG4CoreApplication")
332  << "ParametrisedEMPhysics: printout level changed for " << part->GetParticleName();
333  }
334 }
ParametrisedEMPhysics::TLSmod::theFastSimulationManagerProcess
std::unique_ptr< G4FastSimulationManagerProcess > theFastSimulationManagerProcess
Definition: ParametrisedEMPhysics.cc:60
ElectronLimiter
Definition: ElectronLimiter.h:23
mps_fire.i
i
Definition: mps_fire.py:355
MessageLogger.h
ParametrisedEMPhysics::TLSmod::theEcalEMShowerModel
std::unique_ptr< GFlashEMShowerModel > theEcalEMShowerModel
Definition: ParametrisedEMPhysics.cc:55
nt
int nt
Definition: AMPTWrapper.h:42
ElectronLimiter::SetTrackingCutPerRegion
void SetTrackingCutPerRegion(std::vector< const G4Region * > &, std::vector< G4double > &, std::vector< G4double > &, std::vector< G4double > &)
Definition: ElectronLimiter.cc:111
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
ParametrisedEMPhysics::TLSmod::theHcalEMShowerModel
std::unique_ptr< GFlashEMShowerModel > theHcalEMShowerModel
Definition: ParametrisedEMPhysics.cc:57
ParametrisedEMPhysics::TLSmod::theLowEnergyFastSimModel
std::unique_ptr< LowEnergyFastSimModel > theLowEnergyFastSimModel
Definition: ParametrisedEMPhysics.cc:56
LowEnergyFastSimModel.h
MeV
const double MeV
GFlashEMShowerModel
Definition: GFlashEMShowerModel.h:29
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
GFlashEMShowerModel.h
ParametrisedEMPhysics::ParametrisedEMPhysics
ParametrisedEMPhysics(const std::string &name, const edm::ParameterSet &p)
Definition: ParametrisedEMPhysics.cc:65
ParametrisedEMPhysics::TLSmod::theHcalHadShowerModel
std::unique_ptr< GFlashHadronShowerModel > theHcalHadShowerModel
Definition: ParametrisedEMPhysics.cc:59
ParametrisedEMPhysics::ModifyTransportation
void ModifyTransportation(const G4ParticleDefinition *, int ntry, double th1, double th2)
Definition: ParametrisedEMPhysics.cc:324
part
part
Definition: HCALResponse.h:20
ElectronLimiter.h
ParametrisedEMPhysics::TLSmod::theEcalHadShowerModel
std::unique_ptr< GFlashHadronShowerModel > theEcalHadShowerModel
Definition: ParametrisedEMPhysics.cc:58
ElectronLimiter::SetFieldCheckFlag
void SetFieldCheckFlag(G4bool)
Definition: ElectronLimiter.h:70
ParametrisedEMPhysics::TLSmod
Definition: ParametrisedEMPhysics.cc:54
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::LogWarning
Definition: MessageLogger.h:141
edm::ParameterSet
Definition: ParameterSet.h:36
ParametrisedEMPhysics.h
rname
const G4String rname[NREG]
Definition: ParametrisedEMPhysics.cc:46
GeV
const double GeV
Definition: MathUtil.h:16
ParametrisedEMPhysics::theParSet
edm::ParameterSet theParSet
Definition: ParametrisedEMPhysics.h:27
edm::LogVerbatim
Definition: MessageLogger.h:297
gem
Definition: AMC13Event.h:6
NREG
const G4int NREG
Definition: ParametrisedEMPhysics.cc:45
alignCSCRings.r
r
Definition: alignCSCRings.py:93
ElectronLimiter::SetRangeCheckFlag
void SetRangeCheckFlag(G4bool)
Definition: ElectronLimiter.h:68
GFlashHadronShowerModel
Definition: GFlashHadronShowerModel.h:20
nanoDQM_cff.Electron
Electron
Definition: nanoDQM_cff.py:62
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
LowEnergyFastSimModel
Definition: LowEnergyFastSimModel.h:13
ParametrisedEMPhysics::~ParametrisedEMPhysics
~ParametrisedEMPhysics() override
Definition: ParametrisedEMPhysics.cc:112
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
ParametrisedEMPhysics::ConstructParticle
void ConstructParticle() override
Definition: ParametrisedEMPhysics.cc:119
GFlashHadronShowerModel.h
ParametrisedEMPhysics::ConstructProcess
void ConstructProcess() override
Definition: ParametrisedEMPhysics.cc:136
ParametrisedEMPhysics::m_tpmod
static G4ThreadLocal TLSmod * m_tpmod
Definition: ParametrisedEMPhysics.h:29
Generated for CMSSW Reference Manual by doxygen 1.8.16