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