CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
HCalSD.cc
Go to the documentation of this file.
1 // File: HCalSD.cc
3 // Description: Sensitive Detector class for calorimeters
5 
13 
18 
19 #include "G4LogicalVolumeStore.hh"
20 #include "G4LogicalVolume.hh"
21 #include "G4Step.hh"
22 #include "G4Track.hh"
23 
24 #include "G4SystemOfUnits.hh"
25 #include "G4PhysicalConstants.hh"
26 #include "Randomize.hh"
27 
28 #include "DD4hep/Filter.h"
29 
30 #include <iostream>
31 #include <fstream>
32 #include <iomanip>
33 #include <sstream>
34 
35 //#define EDM_ML_DEBUG
36 //#define plotDebug
37 //#define printDebug
38 
39 #ifdef plotDebug
40 #include <TH1F.h>
41 #endif
42 
44  const HcalDDDSimConstants* hcns,
45  const HcalDDDRecConstants* hcnr,
46  const HcalSimulationConstants* hscs,
47  const HBHEDarkening* hbd,
48  const HBHEDarkening* hed,
49  const SensitiveDetectorCatalog& clg,
50  edm::ParameterSet const& p,
51  const SimTrackManager* manager)
52  : CaloSD(name,
53  clg,
54  p,
55  manager,
56  (float)(p.getParameter<edm::ParameterSet>("HCalSD").getParameter<double>("TimeSliceUnit")),
57  p.getParameter<edm::ParameterSet>("HCalSD").getParameter<bool>("IgnoreTrackID")),
58  hcalConstants_(hcns),
59  hcalSimConstants_(hscs),
60  m_HBDarkening(hbd),
61  m_HEDarkening(hed),
62  isHF(false),
63  weight_(1.0),
64  depth_(1) {
65  numberingFromDDD.reset(nullptr);
66  numberingScheme.reset(nullptr);
67  showerLibrary.reset(nullptr);
68  hfshower.reset(nullptr);
69  showerParam.reset(nullptr);
70  showerPMT.reset(nullptr);
71  showerBundle.reset(nullptr);
72  m_HFDarkening.reset(nullptr);
73  m_HcalTestNS.reset(nullptr);
74 
75  //static SimpleConfigurable<double> bk1(0.013, "HCalSD:BirkC1");
76  //static SimpleConfigurable<double> bk2(0.0568,"HCalSD:BirkC2");
77  //static SimpleConfigurable<double> bk3(1.75, "HCalSD:BirkC3");
78  // Values from NIM 80 (1970) 239-244: as implemented in Geant3
79 
80  bool dd4hep = p.getParameter<bool>("g4GeometryDD4hepSource");
81  edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("HCalSD");
82  useBirk = m_HC.getParameter<bool>("UseBirkLaw");
83  double bunit = (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
84  birk1 = m_HC.getParameter<double>("BirkC1") * bunit;
85  birk2 = m_HC.getParameter<double>("BirkC2");
86  birk3 = m_HC.getParameter<double>("BirkC3");
87  useShowerLibrary = m_HC.getParameter<bool>("UseShowerLibrary");
88  useParam = m_HC.getParameter<bool>("UseParametrize");
89  testNumber = m_HC.getParameter<bool>("TestNumberingScheme");
90  neutralDensity = m_HC.getParameter<bool>("doNeutralDensityFilter");
91  usePMTHit = m_HC.getParameter<bool>("UsePMTHits");
92  betaThr = m_HC.getParameter<double>("BetaThreshold");
93  eminHitHB = m_HC.getParameter<double>("EminHitHB") * MeV;
94  eminHitHE = m_HC.getParameter<double>("EminHitHE") * MeV;
95  eminHitHO = m_HC.getParameter<double>("EminHitHO") * MeV;
96  eminHitHF = m_HC.getParameter<double>("EminHitHF") * MeV;
97  useFibreBundle = m_HC.getParameter<bool>("UseFibreBundleHits");
98  deliveredLumi = m_HC.getParameter<double>("DelivLuminosity");
99  agingFlagHB = m_HC.getParameter<bool>("HBDarkening");
100  agingFlagHE = m_HC.getParameter<bool>("HEDarkening");
101  bool agingFlagHF = m_HC.getParameter<bool>("HFDarkening");
102  useHF = m_HC.getUntrackedParameter<bool>("UseHF", true);
103  bool forTBHC = m_HC.getUntrackedParameter<bool>("ForTBHCAL", false);
104  bool forTBH2 = m_HC.getUntrackedParameter<bool>("ForTBH2", false);
105  useLayerWt = m_HC.getUntrackedParameter<bool>("UseLayerWt", false);
106  std::string file = m_HC.getUntrackedParameter<std::string>("WtFile", "None");
107  testNS_ = m_HC.getUntrackedParameter<bool>("TestNS", false);
108  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
109  applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
110  bool dumpGeom = m_HC.getUntrackedParameter<bool>("DumpGeometry", false);
111 
112 #ifdef EDM_ML_DEBUG
113  edm::LogVerbatim("HcalSim") << "***************************************************"
114  << "\n"
115  << "* Constructing a HCalSD with name " << name << "\n"
116  << "\n"
117  << "***************************************************";
118 #endif
119  edm::LogVerbatim("HcalSim") << "HCalSD:: Use of HF code is set to " << useHF
120  << "\n Use of shower parametrization set to " << useParam
121  << "\n Use of shower library is set to " << useShowerLibrary
122  << "\n Use PMT Hit is set to " << usePMTHit << " with beta Threshold " << betaThr
123  << "\n USe of FibreBundle Hit set to " << useFibreBundle
124  << "\n Use of Birks law is set to " << useBirk
125  << " with three constants kB = " << birk1 / bunit << ", C1 = " << birk2
126  << ", C2 = " << birk3;
127  edm::LogVerbatim("HcalSim") << "HCalSD:: Suppression Flag " << suppressHeavy << " protons below " << kmaxProton
128  << " MeV,"
129  << " neutrons below " << kmaxNeutron << " MeV and"
130  << " ions below " << kmaxIon << " MeV\n"
131  << " Threshold for storing hits in HB: " << eminHitHB << " HE: " << eminHitHE
132  << " HO: " << eminHitHO << " HF: " << eminHitHF << "\n"
133  << "Delivered luminosity for Darkening " << deliveredLumi << " Flag (HE) " << agingFlagHE
134  << " Flag (HB) " << agingFlagHB << " Flag (HF) " << agingFlagHF << "\n"
135  << "Application of Fiducial Cut " << applyFidCut
136  << "Flag for test number|neutral density filter " << testNumber << " " << neutralDensity;
137 
138  if (forTBHC) {
139  useHF = false;
140  matNames.emplace_back("Scintillator");
141  } else {
143  }
144 
146  if (testNumber || forTBH2) {
147  scheme = dynamic_cast<HcalNumberingScheme*>(new HcalTestNumberingScheme(forTBH2));
148  } else {
149  scheme = new HcalNumberingScheme();
150  }
152 
153  // always call getFromLibrary() method to identify HF region
154  setParameterized(true);
155 
156  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
157  // std::vector<G4LogicalVolume *>::const_iterator lvcite;
158  const G4LogicalVolume* lv;
159  std::string attribute, value;
160 
161  if (useHF) {
162  if (useParam) {
163  showerParam = std::make_unique<HFShowerParam>(name, hcalConstants_, hcalSimConstants_->hcalsimpar(), p);
164  } else {
165  if (useShowerLibrary) {
166  showerLibrary = std::make_unique<HFShowerLibrary>(name, hcalConstants_, hcalSimConstants_->hcalsimpar(), p);
167  }
168  hfshower = std::make_unique<HFShower>(name, hcalConstants_, hcalSimConstants_->hcalsimpar(), p, 0);
169  }
170 
171  // HF volume names
173  const std::vector<int>& temp = hcalSimConstants_->hcalsimpar()->hfLevels_;
174 #ifdef EDM_ML_DEBUG
175  std::stringstream ss0;
176  ss0 << "HCalSD: Names to be tested for Volume = HF has " << hfNames.size() << " elements";
177 #endif
178  int addlevel = dd4hep ? 1 : 0;
179  for (unsigned int i = 0; i < hfNames.size(); ++i) {
180  G4String namv(static_cast<std::string>(dd4hep::dd::noNamespace(hfNames[i])));
181  lv = nullptr;
182  for (auto lvol : *lvs) {
183  if (dd4hep::dd::noNamespace(lvol->GetName()) == namv) {
184  lv = lvol;
185  break;
186  }
187  }
188  hfLV.emplace_back(lv);
189  hfLevels.emplace_back(temp[i] + addlevel);
190 #ifdef EDM_ML_DEBUG
191  ss0 << "\n HF[" << i << "] = " << namv << " LV " << lv << " at level " << (temp[i] + addlevel);
192 #endif
193  }
194 #ifdef EDM_ML_DEBUG
195  edm::LogVerbatim("HcalSim") << ss0.str();
196 #endif
197  // HF Fibre volume names
200  const std::vector<std::string>& pmtNames = hcalSimConstants_->hcalsimpar()->hfPMTNames_;
201  fillLogVolumeVector("HFPMT", pmtNames, pmtLV);
202  const std::vector<std::string>& straightNames = hcalSimConstants_->hcalsimpar()->hfFibreStraightNames_;
203  fillLogVolumeVector("HFFibreBundleStraight", straightNames, fibre1LV);
204  const std::vector<std::string>& conicalNames = hcalSimConstants_->hcalsimpar()->hfFibreConicalNames_;
205  fillLogVolumeVector("HFFibreBundleConical", conicalNames, fibre2LV);
206  }
207 
208  //Material list for HB/HE/HO sensitive detectors
209  const G4MaterialTable* matTab = G4Material::GetMaterialTable();
210  std::vector<G4Material*>::const_iterator matite;
211  for (auto const& namx : matNames) {
212  const G4Material* mat = nullptr;
213  for (matite = matTab->begin(); matite != matTab->end(); ++matite) {
214  if (static_cast<std::string>(dd4hep::dd::noNamespace((*matite)->GetName())) == namx) {
215  mat = (*matite);
216  break;
217  }
218  }
219  materials.emplace_back(mat);
220  }
221 #ifdef EDM_ML_DEBUG
222  std::stringstream ss1;
223  for (unsigned int i = 0; i < matNames.size(); ++i) {
224  if (i / 10 * 10 == i) {
225  ss1 << "\n";
226  }
227  ss1 << " " << matNames[i];
228  }
229  edm::LogVerbatim("HcalSim") << "HCalSD: Material names for HCAL: " << ss1.str();
230 #endif
231  if (useLayerWt) {
233  }
234  numberingFromDDD = std::make_unique<HcalNumberingFromDDD>(hcalConstants_);
235 
236  //Special Geometry parameters
238 #ifdef EDM_ML_DEBUG
239  std::stringstream ss2;
240  for (unsigned int ig = 0; ig < gpar.size(); ig++) {
241  ss2 << "\n gpar[" << ig << "] = " << gpar[ig] / cm << " cm";
242  }
243  edm::LogVerbatim("HcalSim") << "Maximum depth for HF " << hcalConstants_->getMaxDepth(2) << gpar.size()
244  << " gpar (cm)" << ss2.str();
245 #endif
246 
247  //Test Hcal Numbering Scheme
248  if (testNS_)
249  m_HcalTestNS = std::make_unique<HcalTestNS>(hcnr);
250 
251  for (int i = 0; i < 9; ++i) {
252  hit_[i] = time_[i] = dist_[i] = nullptr;
253  }
254  hzvem = hzvhad = nullptr;
255 
256  if (agingFlagHF) {
257  m_HFDarkening = std::make_unique<HFDarkening>(m_HC.getParameter<edm::ParameterSet>("HFDarkeningParameterBlock"));
258  }
259 #ifdef plotDebug
261 
262  if (tfile.isAvailable()) {
263  static const char* const labels[] = {"HB",
264  "HE",
265  "HO",
266  "HF Absorber",
267  "HF PMT",
268  "HF Absorber Long",
269  "HF Absorber Short",
270  "HF PMT Long",
271  "HF PMT Short"};
272  TFileDirectory hcDir = tfile->mkdir("ProfileFromHCalSD");
273  char name[20], title[60];
274  for (int i = 0; i < 9; ++i) {
275  sprintf(title, "Hit energy in %s", labels[i]);
276  sprintf(name, "HCalSDHit%d", i);
277  hit_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
278  sprintf(title, "Energy (MeV)");
279  hit_[i]->GetXaxis()->SetTitle(title);
280  hit_[i]->GetYaxis()->SetTitle("Hits");
281  sprintf(title, "Time of the hit in %s", labels[i]);
282  sprintf(name, "HCalSDTime%d", i);
283  time_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
284  sprintf(title, "Time (ns)");
285  time_[i]->GetXaxis()->SetTitle(title);
286  time_[i]->GetYaxis()->SetTitle("Hits");
287  sprintf(title, "Longitudinal profile in %s", labels[i]);
288  sprintf(name, "HCalSDDist%d", i);
289  dist_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
290  sprintf(title, "Distance (mm)");
291  dist_[i]->GetXaxis()->SetTitle(title);
292  dist_[i]->GetYaxis()->SetTitle("Hits");
293  }
294  if (useHF && (!useParam)) {
295  hzvem = hcDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part)", 330, 0.0, 1650.0);
296  hzvem->GetXaxis()->SetTitle("Longitudinal Profile (EM Part)");
297  hzvhad = hcDir.make<TH1F>("hzvhad", "Longitudinal Profile (Had Part)", 330, 0.0, 1650.0);
298  hzvhad->GetXaxis()->SetTitle("Longitudinal Profile (Hadronic Part)");
299  }
300  }
301 #endif
302  if (dumpGeom) {
303  std::unique_ptr<HcalNumberingFromDDD> hcn = std::make_unique<HcalNumberingFromDDD>(hcalConstants_);
304  const auto& lvNames = clg.logicalNames(name);
305  HcalDumpGeometry geom(lvNames, hcn.get(), testNumber, false);
306  geom.update();
307  }
308 }
309 
311  const std::vector<std::string>& lvnames,
312  std::vector<const G4LogicalVolume*>& lvvec) {
313  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
314  const G4LogicalVolume* lv;
315  std::stringstream ss3;
316  ss3 << "HCalSD: " << lvnames.size() << " names to be tested for Volume <" << value << ">:";
317  for (unsigned int i = 0; i < lvnames.size(); ++i) {
318  G4String namv(static_cast<std::string>(dd4hep::dd::noNamespace(lvnames[i])));
319  lv = nullptr;
320  for (auto lvol : *lvs) {
321  if (dd4hep::dd::noNamespace(lvol->GetName()) == namv) {
322  lv = lvol;
323  break;
324  }
325  }
326  lvvec.emplace_back(lv);
327  if (i / 10 * 10 == i) {
328  ss3 << "\n";
329  }
330  ss3 << " " << namv;
331  }
332  edm::LogVerbatim("HcalSim") << ss3.str();
333 }
334 
335 bool HCalSD::getFromLibrary(const G4Step* aStep) {
336  auto const track = aStep->GetTrack();
337  depth_ = (aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(0)) % 10;
338  weight_ = 1.0;
339  bool kill(false);
340  isHF = isItHF(aStep);
341 #ifdef EDM_ML_DEBUG
342  edm::LogVerbatim("HcalSim") << "GetFromLibrary: isHF " << isHF << " darken " << (m_HFDarkening != nullptr)
343  << " useParam " << useParam << " useShowerLibrary " << useShowerLibrary << " Muon? "
344  << G4TrackToParticleID::isMuon(track) << " electron? "
345  << G4TrackToParticleID::isGammaElectronPositron(track) << " Stable Hadron? "
347 #endif
348  if (isHF) {
349  if (m_HFDarkening) {
350  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
351  const double invcm = 1. / CLHEP::cm;
352  double r = hitPoint.perp() * invcm;
353  double z = std::abs(hitPoint.z()) * invcm;
354  double dose_acquired = 0.;
356  unsigned int hfZLayer = (unsigned int)((z - HFDarkening::lowZLimit) / 5);
357  if (hfZLayer >= HFDarkening::upperZLimit)
358  hfZLayer = (HFDarkening::upperZLimit - 1);
359  float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
360  for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
361  dose_acquired = m_HFDarkening->dose(i, r);
362  weight_ *= m_HFDarkening->degradation(normalized_lumi * dose_acquired);
363  }
364  }
365 #ifdef EDM_ML_DEBUG
366  edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary: HFLumiDarkening at "
367  << "r= " << r << ", z= " << z << " Dose= " << dose_acquired << " weight= " << weight_;
368 #endif
369  }
370 
371  if (useParam) {
372  getFromParam(aStep, kill);
373 #ifdef EDM_ML_DEBUG
374  G4String nameVolume = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
375  edm::LogVerbatim("HcalSim") << "HCalSD: " << getNumberOfHits() << " hits from parametrization in " << nameVolume
376  << " for Track " << track->GetTrackID() << " ("
377  << track->GetDefinition()->GetParticleName() << ")";
378 #endif
381 #ifdef EDM_ML_DEBUG
382  auto nameVolume = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
383  edm::LogVerbatim("HcalSim") << "HCalSD: Starts shower library from " << nameVolume << " for Track "
384  << track->GetTrackID() << " (" << track->GetDefinition()->GetParticleName() << ")";
385 
386 #endif
387  getFromHFLibrary(aStep, kill);
388  }
389  }
390  }
391 #ifdef EDM_ML_DEBUG
392  edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary ID= " << track->GetTrackID() << " ("
393  << track->GetDefinition()->GetParticleName() << ") kill= " << kill
394  << " weight= " << weight_ << " depth= " << depth_ << " isHF: " << isHF;
395 #endif
396  return kill;
397 }
398 
399 double HCalSD::getEnergyDeposit(const G4Step* aStep) {
400  double destep(0.0);
401  auto const lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
402  auto const theTrack = aStep->GetTrack();
403 
404  if (isHF) {
405  if (useShowerLibrary && G4TrackToParticleID::isMuon(theTrack) && isItFibre(lv)) {
406 #ifdef EDM_ML_DEBUG
407  edm::LogVerbatim("HcalSim") << "HCalSD: Hit at Fibre in LV " << lv->GetName() << " for track "
408  << aStep->GetTrack()->GetTrackID() << " ("
409  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
410 #endif
411  hitForFibre(aStep);
412  }
413  return destep;
414  }
415 
416  if (isItPMT(lv)) {
417  if (usePMTHit && showerPMT) {
418  getHitPMT(aStep);
419  }
420 #ifdef EDM_ML_DEBUG
421  edm::LogVerbatim("HcalSim") << "HCalSD: Hit from PMT parametrization in LV " << lv->GetName() << " for Track "
422  << aStep->GetTrack()->GetTrackID() << " ("
423  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
424 #endif
425  return destep;
426 
427  } else if (isItStraightBundle(lv)) {
428  if (useFibreBundle && showerBundle) {
429  getHitFibreBundle(aStep, false);
430  }
431 #ifdef EDM_ML_DEBUG
432  edm::LogVerbatim("HcalSim") << "HCalSD: Hit from straight FibreBundle in LV: " << lv->GetName() << " for track "
433  << aStep->GetTrack()->GetTrackID() << " ("
434  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
435 #endif
436  return destep;
437 
438  } else if (isItConicalBundle(lv)) {
439  if (useFibreBundle && showerBundle) {
440  getHitFibreBundle(aStep, true);
441  }
442 #ifdef EDM_ML_DEBUG
443  edm::LogVerbatim("HcalSim") << "HCalSD: Hit from conical FibreBundle PV: " << lv->GetName() << " for track "
444  << aStep->GetTrack()->GetTrackID() << " ("
445  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
446 #endif
447  return destep;
448  }
449 
450  // normal hit
451  destep = aStep->GetTotalEnergyDeposit();
452 
453  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
454  uint32_t detid = setDetUnitId(aStep);
455  int det(0), ieta(0), phi(0), z(0), lay, depth(-1);
456  if (testNumber) {
457  HcalTestNumbering::unpackHcalIndex(detid, det, z, depth, ieta, phi, lay);
458  if (z == 0) {
459  z = -1;
460  }
461  } else {
462  HcalDetId hcid(detid);
463  det = hcid.subdetId();
464  ieta = hcid.ietaAbs();
465  phi = hcid.iphi();
466  z = hcid.zside();
467  }
468  lay = (touch->GetReplicaNumber(0) / 10) % 100 + 1;
469 #ifdef EDM_ML_DEBUG
470  edm::LogVerbatim("HcalSim") << "HCalSD: det: " << det << " ieta: " << ieta << " iphi: " << phi << " zside " << z
471  << " lay: " << lay - 2;
472 #endif
473  if (depth_ == 0 && (det == 1 || det == 2) && ((!testNumber) || neutralDensity))
475  if (useLayerWt) {
476  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
477  weight_ = layerWeight(det + 2, hitPoint, depth_, lay);
478  }
479 
480  if (agingFlagHB && m_HBDarkening && det == 1) {
481  double dweight = m_HBDarkening->degradation(deliveredLumi, ieta, lay);
482  weight_ *= dweight;
483 #ifdef EDM_ML_DEBUG
484  edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi << " coefficient = " << dweight
485  << " Weight= " << weight_;
486 #endif
487  }
488 
489  if (agingFlagHE && m_HEDarkening && det == 2) {
490  double dweight = m_HEDarkening->degradation(deliveredLumi, ieta, lay);
491  weight_ *= dweight;
492 #ifdef EDM_ML_DEBUG
493  edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi << " coefficient = " << dweight
494  << " Weight= " << weight_;
495 #endif
496  }
497 
498  if (suppressHeavy) {
499  TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
500  if (trkInfo) {
501  int pdg = theTrack->GetDefinition()->GetPDGEncoding();
502  if (!(trkInfo->isPrimary())) { // Only secondary particles
503  double ke = theTrack->GetKineticEnergy();
504  if (pdg / 1000000000 == 1 && (pdg / 10000) % 100 > 0 && (pdg / 10) % 100 > 0 && ke < kmaxIon)
505  weight_ = 0;
506  if ((pdg == 2212) && (ke < kmaxProton))
507  weight_ = 0;
508  if ((pdg == 2112) && (ke < kmaxNeutron))
509  weight_ = 0;
510  }
511  }
512  }
513  double wt0(1.0);
514  if (useBirk) {
515  const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
516  if (isItScintillator(mat))
517  wt0 = getAttenuation(aStep, birk1, birk2, birk3);
518  }
519  weight_ *= wt0;
520  double wt1 = getResponseWt(theTrack);
521  double wt2 = theTrack->GetWeight();
522  double edep = weight_ * wt1 * destep;
523  if (wt2 > 0.0) {
524  edep *= wt2;
525  }
526 #ifdef EDM_ML_DEBUG
527  edm::LogVerbatim("HcalSim") << "HCalSD: edep= " << edep << " Det: " << det + 2 << " depth= " << depth_
528  << " weight= " << weight_ << " wt0= " << wt0 << " wt1= " << wt1 << " wt2= " << wt2;
529 #endif
530  return edep;
531 }
532 
533 uint32_t HCalSD::setDetUnitId(const G4Step* aStep) {
534  auto const prePoint = aStep->GetPreStepPoint();
535  auto const touch = prePoint->GetTouchable();
536  const G4ThreeVector& hitPoint = prePoint->GetPosition();
537 
538  int depth = (touch->GetReplicaNumber(0)) % 10 + 1;
539  int lay = (touch->GetReplicaNumber(0) / 10) % 100 + 1;
540  int det = (touch->GetReplicaNumber(1)) / 1000;
541 
542  uint32_t idx = setDetUnitId(det, hitPoint, depth, lay);
543 #ifdef EDM_ML_DEBUG
544  if (depth == 1) {
545  edm::LogVerbatim("HcalSim") << "HCalSD: Check for " << det << ":" << depth << ":" << lay << " ID " << std::hex
546  << idx << std::dec;
547  int det0, z0, depth0, eta0, phi0, lay0(-1);
548  if (testNumber) {
549  HcalTestNumbering::unpackHcalIndex(idx, det0, z0, depth0, eta0, phi0, lay0);
550  } else {
551  HcalDetId hcid0(idx);
552  det0 = hcid0.subdetId();
553  eta0 = hcid0.ietaAbs();
554  phi0 = hcid0.iphi();
555  z0 = hcid0.zside();
556  depth0 = hcid0.depth();
557  }
558  edm::LogVerbatim("HcalSim") << "HCalSD: det|z|depth|eta|phi|lay " << det0 << ":" << z0 << ":" << depth0 << ":"
559  << eta0 << ":" << phi0 << ":" << lay0;
560  printVolume(touch);
561  }
562 #endif
563  return idx;
564 }
565 
567  if (scheme != nullptr) {
568  edm::LogVerbatim("HcalSim") << "HCalSD: updates numbering scheme for " << GetName();
569  numberingScheme.reset(scheme);
570  }
571 }
572 
573 void HCalSD::update(const BeginOfJob* job) {}
574 
576 
577 bool HCalSD::filterHit(CaloG4Hit* aHit, double time) {
578  double threshold = 0;
579  DetId theId(aHit->getUnitID());
580  switch (theId.subdetId()) {
581  case HcalBarrel:
583  break;
584  case HcalEndcap:
586  break;
587  case HcalOuter:
589  break;
590  case HcalForward:
592  break;
593  default:
594  break;
595  }
596  return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > threshold));
597 }
598 
599 uint32_t HCalSD::setDetUnitId(int det, const G4ThreeVector& pos, int depth, int lay = 1) {
600  uint32_t id = 0;
601  if (det == 0) {
602 #ifdef printDebug
603  double eta = std::abs(pos.eta());
604 #endif
605  if (std::abs(pos.z()) > maxZ_) {
606  det = 5;
607 #ifdef printDebug
608  if (eta < 2.868)
609  ++detNull_[2];
610 #endif
611  } else if (!(hcalConstants_->isHE())) {
612  det = 3;
613 #ifdef printDebug
614  ++detNull_[0];
615 #endif
616  } else {
617  double minR = minRoff_ + slopeHE_ * std::abs(pos.z());
618  double maxR = maxRoff_ + slopeHE_ * std::abs(pos.z());
619  det = ((pos.perp() > minR) && (pos.perp() < maxR)) ? 4 : 3;
620 #ifdef printDebug
621  ++detNull_[det - 3];
622 #endif
623  }
624 #ifdef printDEBUG
625  edm::LogVerbatim("HcalSim") << "Position " << pos.perp() << ":" << std::abs(pos.z()) << " Limits "
626  << !(hcalConstants_->isHE()) << ":" << maxZ_ << " det " << det;
627  } else {
628  ++detNull_[3];
629 #endif
630  }
631 
632  if (numberingFromDDD.get()) {
633  //get the ID's as eta, phi, depth, ... indices
635  numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
636  id = setDetUnitId(tmp);
637  }
638  return id;
639 }
640 
642  modifyDepth(tmp);
643  uint32_t id = (numberingScheme.get()) ? numberingScheme->getUnitID(tmp) : 0;
644  if ((!testNumber) && m_HcalTestNS.get()) {
645  bool ok = m_HcalTestNS->compare(tmp, id);
646  if (!ok)
647  edm::LogWarning("HcalSim") << "Det ID from HCalSD " << HcalDetId(id) << " " << std::hex << id << std::dec
648  << " does not match one from relabller for " << tmp.subdet << ":" << tmp.etaR << ":"
649  << tmp.phi << ":" << tmp.phis << ":" << tmp.depth << ":" << tmp.lay << std::endl;
650  }
651  return id;
652 }
653 
654 bool HCalSD::isItHF(const G4Step* aStep) {
655  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
656  int levels = (touch->GetHistoryDepth()) + 1;
657  for (unsigned int it = 0; it < hfNames.size(); ++it) {
658  if (levels >= hfLevels[it]) {
659  const G4LogicalVolume* lv = touch->GetVolume(levels - hfLevels[it])->GetLogicalVolume();
660  if (lv == hfLV[it])
661  return true;
662  }
663  }
664  return false;
665 }
666 
667 bool HCalSD::isItHF(const G4String& name) {
668  for (const auto& nam : hfNames)
669  if (name == static_cast<G4String>(nam)) {
670  return true;
671  }
672  return false;
673 }
674 
675 bool HCalSD::isItFibre(const G4LogicalVolume* lv) {
676  for (auto lvol : fibreLV)
677  if (lv == lvol) {
678  return true;
679  }
680  return false;
681 }
682 
683 bool HCalSD::isItFibre(const G4String& name) {
684  for (const auto& nam : fibreNames)
685  if (name == static_cast<G4String>(nam)) {
686  return true;
687  }
688  return false;
689 }
690 
691 bool HCalSD::isItPMT(const G4LogicalVolume* lv) {
692  for (auto lvol : pmtLV)
693  if (lv == lvol) {
694  return true;
695  }
696  return false;
697 }
698 
699 bool HCalSD::isItStraightBundle(const G4LogicalVolume* lv) {
700  for (auto lvol : fibre1LV)
701  if (lv == lvol) {
702  return true;
703  }
704  return false;
705 }
706 
707 bool HCalSD::isItConicalBundle(const G4LogicalVolume* lv) {
708  for (auto lvol : fibre2LV)
709  if (lv == lvol) {
710  return true;
711  }
712  return false;
713 }
714 
715 bool HCalSD::isItScintillator(const G4Material* mat) {
716  for (auto amat : materials)
717  if (amat == mat) {
718  return true;
719  }
720  return false;
721 }
722 
723 bool HCalSD::isItinFidVolume(const G4ThreeVector& hitPoint) {
724  bool flag = true;
725  if (applyFidCut) {
726  int npmt = HFFibreFiducial::PMTNumber(hitPoint);
727 #ifdef EDM_ML_DEBUG
728  edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume:#PMT= " << npmt << " for hit point " << hitPoint;
729 #endif
730  if (npmt <= 0)
731  flag = false;
732  }
733 #ifdef EDM_ML_DEBUG
734  edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume: point " << hitPoint << " return flag " << flag;
735 #endif
736  return flag;
737 }
738 
739 void HCalSD::getFromHFLibrary(const G4Step* aStep, bool& isKilled) {
740  std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, isKilled, weight_, false);
741  if (!isKilled || hits.empty()) {
742  return;
743  }
744 
745  int primaryID = setTrackID(aStep);
746 
747  // Reset entry point for new primary
748  resetForNewPrimary(aStep);
749 
750  auto const theTrack = aStep->GetTrack();
751  int det = 5;
752 
754  edepositEM = 1. * GeV;
755  edepositHAD = 0.;
756  } else {
757  edepositEM = 0.;
758  edepositHAD = 1. * GeV;
759  }
760 #ifdef EDM_ML_DEBUG
761  edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
762  << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
763  << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV << " GeV";
764 #endif
765  for (unsigned int i = 0; i < hits.size(); ++i) {
766  G4ThreeVector hitPoint = hits[i].position;
767  if (isItinFidVolume(hitPoint)) {
768  int depth = hits[i].depth;
769  double time = hits[i].time;
770  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
771  currentID.setID(unitID, time, primaryID, 0);
772 #ifdef plotDebug
773  plotProfile(aStep, hitPoint, 1.0 * GeV, time, depth);
774  bool emType = G4TrackToParticleID::isGammaElectronPositron(theTrack->GetDefinition()->GetPDGEncoding());
775  plotHF(hitPoint, emType);
776 #endif
777  processHit(aStep);
778  }
779  }
780 }
781 
782 void HCalSD::hitForFibre(const G4Step* aStep) { // if not ParamShower
783 
784  std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight_);
785  if (hits.empty()) {
786  return;
787  }
788 
789  auto const theTrack = aStep->GetTrack();
790  int primaryID = setTrackID(aStep);
791  int det = 5;
792 
794  edepositEM = 1. * GeV;
795  edepositHAD = 0.;
796  } else {
797  edepositEM = 0.;
798  edepositHAD = 1. * GeV;
799  }
800 
801 #ifdef EDM_ML_DEBUG
802  edm::LogVerbatim("HcalSim") << "HCalSD::hitForFibre " << hits.size() << " hits for " << GetName() << " of "
803  << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
804  << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV << " GeV in detector type " << det;
805 #endif
806 
807  for (unsigned int i = 0; i < hits.size(); ++i) {
808  G4ThreeVector hitPoint = hits[i].position;
809  if (isItinFidVolume(hitPoint)) {
810  int depth = hits[i].depth;
811  double time = hits[i].time;
812  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
813  currentID.setID(unitID, time, primaryID, 0);
814 #ifdef plotDebug
815  plotProfile(aStep, hitPoint, edepositEM, time, depth);
816  bool emType = (edepositEM > 0.) ? true : false;
817  plotHF(hitPoint, emType);
818 #endif
819  processHit(aStep);
820  }
821  }
822 }
823 
824 void HCalSD::getFromParam(const G4Step* aStep, bool& isKilled) {
825  std::vector<HFShowerParam::Hit> hits = showerParam->getHits(aStep, weight_, isKilled);
826  if (!isKilled || hits.empty()) {
827  return;
828  }
829 
830  int primaryID = setTrackID(aStep);
831  int det = 5;
832 
833 #ifdef EDM_ML_DEBUG
834  edm::LogVerbatim("HcalSim") << "HCalSD::getFromParam " << hits.size() << " hits for " << GetName() << " of "
835  << primaryID << " with " << aStep->GetTrack()->GetDefinition()->GetParticleName()
836  << " of " << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV
837  << " GeV in detector type " << det;
838 #endif
839  for (unsigned int i = 0; i < hits.size(); ++i) {
840  G4ThreeVector hitPoint = hits[i].position;
841  int depth = hits[i].depth;
842  double time = hits[i].time;
843  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
844  currentID.setID(unitID, time, primaryID, 0);
845  edepositEM = hits[i].edep * GeV;
846  edepositHAD = 0.;
847 #ifdef plotDebug
848  plotProfile(aStep, hitPoint, edepositEM, time, depth);
849 #endif
850  processHit(aStep);
851  }
852 }
853 
854 void HCalSD::getHitPMT(const G4Step* aStep) {
855  auto const preStepPoint = aStep->GetPreStepPoint();
856  auto const theTrack = aStep->GetTrack();
857  double edep = showerPMT->getHits(aStep);
858 
859  if (edep >= 0.) {
860  edep *= GeV;
861  double etrack = preStepPoint->GetKineticEnergy();
862  int primaryID = 0;
863  if (etrack >= energyCut) {
864  primaryID = theTrack->GetTrackID();
865  } else {
866  primaryID = theTrack->GetParentID();
867  if (primaryID == 0)
868  primaryID = theTrack->GetTrackID();
869  }
870  // Reset entry point for new primary
871  resetForNewPrimary(aStep);
872  //
873  int det = static_cast<int>(HcalForward);
874  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
875  double rr = (hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y());
876  double phi = (rr == 0. ? 0. : atan2(hitPoint.y(), hitPoint.x()));
877  double etaR = showerPMT->getRadius();
878  int depth = 1;
879  if (etaR < 0) {
880  depth = 2;
881  etaR = -etaR;
882  }
883  if (hitPoint.z() < 0)
884  etaR = -etaR;
885 #ifdef EDM_ML_DEBUG
886  edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR " << etaR << " phi " << phi / deg
887  << " depth " << depth;
888 #endif
889  double time = (aStep->GetPostStepPoint()->GetGlobalTime());
890  uint32_t unitID = 0;
891  if (numberingFromDDD) {
892  HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, etaR, phi, depth, 1);
893  unitID = setDetUnitId(tmp);
894  }
895  currentID.setID(unitID, time, primaryID, 1);
896 
897  edepositHAD = aStep->GetTotalEnergyDeposit();
898  edepositEM = -edepositHAD + edep;
899 #ifdef plotDebug
900  plotProfile(aStep, hitPoint, edep, time, depth);
901 #endif
902 #ifdef EDM_ML_DEBUG
903  double beta = preStepPoint->GetBeta();
904  edm::LogVerbatim("HcalSim") << "HCalSD::getHitPMT 1 hit for " << GetName() << " of " << primaryID << " with "
905  << theTrack->GetDefinition()->GetParticleName() << " of "
906  << preStepPoint->GetKineticEnergy() / GeV << " GeV with velocity " << beta << " UnitID "
907  << std::hex << unitID << std::dec;
908 #endif
909  processHit(aStep);
910  }
911 }
912 
913 void HCalSD::getHitFibreBundle(const G4Step* aStep, bool type) {
914  auto const preStepPoint = aStep->GetPreStepPoint();
915  auto const theTrack = aStep->GetTrack();
916  double edep = showerBundle->getHits(aStep, type);
917 
918  if (edep >= 0.0) {
919  edep *= GeV;
920  double etrack = preStepPoint->GetKineticEnergy();
921  int primaryID = 0;
922  if (etrack >= energyCut) {
923  primaryID = theTrack->GetTrackID();
924  } else {
925  primaryID = theTrack->GetParentID();
926  if (primaryID == 0)
927  primaryID = theTrack->GetTrackID();
928  }
929  // Reset entry point for new primary
930  resetForNewPrimary(aStep);
931  //
932  int det = static_cast<int>(HcalForward);
933  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
934  double rr = hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y();
935  double phi = rr == 0. ? 0. : atan2(hitPoint.y(), hitPoint.x());
936  double etaR = showerBundle->getRadius();
937  int depth = 1;
938  if (etaR < 0.) {
939  depth = 2;
940  etaR = -etaR;
941  }
942  if (hitPoint.z() < 0.)
943  etaR = -etaR;
944 #ifdef EDM_ML_DEBUG
945  edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR " << etaR << " phi " << phi / deg
946  << " depth " << depth;
947 #endif
948  double time = (aStep->GetPostStepPoint()->GetGlobalTime());
949  uint32_t unitID = 0;
950  if (numberingFromDDD) {
951  HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, etaR, phi, depth, 1);
952  unitID = setDetUnitId(tmp);
953  }
954  if (type)
955  currentID.setID(unitID, time, primaryID, 3);
956  else
957  currentID.setID(unitID, time, primaryID, 2);
958 
959  edepositHAD = aStep->GetTotalEnergyDeposit();
960  edepositEM = -edepositHAD + edep;
961 #ifdef plotDebug
962  plotProfile(aStep, hitPoint, edep, time, depth);
963 #endif
964 #ifdef EDM_ML_DEBUG
965  double beta = preStepPoint->GetBeta();
966  edm::LogVerbatim("HcalSim") << "HCalSD::getHitFibreBundle 1 hit for " << GetName() << " of " << primaryID
967  << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
968  << preStepPoint->GetKineticEnergy() / GeV << " GeV with velocity " << beta << " UnitID "
969  << std::hex << unitID << std::dec;
970 #endif
971  processHit(aStep);
972  } // non-zero energy deposit
973 }
974 
976  std::ifstream infile;
977  int entry = 0;
978  infile.open(fName.c_str(), std::ios::in);
979  if (infile) {
980  int det, zside, etaR, phi, lay;
981  double wt;
982  while (infile >> det >> zside >> etaR >> phi >> lay >> wt) {
983  uint32_t id = HcalTestNumbering::packHcalIndex(det, zside, 1, etaR, phi, lay);
984  layerWeights.insert(std::pair<uint32_t, double>(id, wt));
985  ++entry;
986 #ifdef EDM_ML_DEBUG
987  edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile:Entry " << entry << " ID " << std::hex << id
988  << std::dec << " (" << det << "/" << zside << "/1/" << etaR << "/" << phi << "/"
989  << lay << ") Weight " << wt;
990 #endif
991  }
992  infile.close();
993  }
994  edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile: reads " << entry << " weights from " << fName;
995  if (entry <= 0)
996  useLayerWt = false;
997 }
998 
999 double HCalSD::layerWeight(int det, const G4ThreeVector& pos, int depth, int lay) {
1000  double wt = 1.;
1001  if (numberingFromDDD) {
1002  //get the ID's as eta, phi, depth, ... indices
1004  numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
1005  modifyDepth(tmp);
1006  uint32_t id = HcalTestNumbering::packHcalIndex(tmp.subdet, tmp.zside, 1, tmp.etaR, tmp.phis, tmp.lay);
1007  std::map<uint32_t, double>::const_iterator ite = layerWeights.find(id);
1008  if (ite != layerWeights.end())
1009  wt = ite->second;
1010 #ifdef EDM_ML_DEBUG
1011  edm::LogVerbatim("HcalSim") << "HCalSD::layerWeight: ID " << std::hex << id << std::dec << " (" << tmp.subdet << "/"
1012  << tmp.zside << "/1/" << tmp.etaR << "/" << tmp.phis << "/" << tmp.lay << ") Weight "
1013  << wt;
1014 #endif
1015  }
1016  return wt;
1017 }
1018 
1019 void HCalSD::plotProfile(const G4Step* aStep, const G4ThreeVector& global, double edep, double time, int id) {
1020  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
1021  static const unsigned int names = 10;
1022  static const G4String modName[names] = {
1023  "HEModule", "HVQF", "HBModule", "MBAT", "MBBT", "MBBTC", "MBBT_R1P", "MBBT_R1M", "MBBT_R1PX", "MBBT_R1MX"};
1024  G4ThreeVector local;
1025  bool found = false;
1026  double depth = -2000;
1027  int idx = 4;
1028  for (int n = 0; n < touch->GetHistoryDepth(); ++n) {
1029  G4String name(static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(n)->GetName())));
1030 #ifdef EDM_ML_DEBUG
1031  edm::LogVerbatim("HcalSim") << "plotProfile Depth " << n << " Name " << name;
1032 #endif
1033  for (unsigned int ii = 0; ii < names; ++ii) {
1034  if (name == modName[ii]) {
1035  found = true;
1036  int dn = touch->GetHistoryDepth() - n;
1037  local = touch->GetHistory()->GetTransform(dn).TransformPoint(global);
1038  if (ii == 0) {
1039  depth = local.z() - 4006.5;
1040  idx = 1;
1041  } else if (ii == 1) {
1042  depth = local.z() + 825.0;
1043  idx = 3;
1044  } else if (ii == 2) {
1045  depth = local.x() - 1775.;
1046  idx = 0;
1047  } else {
1048  depth = local.y() + 15.;
1049  idx = 2;
1050  }
1051  break;
1052  }
1053  }
1054  if (found)
1055  break;
1056  }
1057  if (!found)
1058  depth = std::abs(global.z()) - 11500;
1059 #ifdef EDM_ML_DEBUG
1060  edm::LogVerbatim("HcalSim") << "plotProfile Found " << found << " Global " << global << " Local " << local
1061  << " depth " << depth << " ID " << id << " EDEP " << edep << " Time " << time;
1062 #endif
1063  if (hit_[idx] != nullptr)
1064  hit_[idx]->Fill(edep);
1065  if (time_[idx] != nullptr)
1066  time_[idx]->Fill(time, edep);
1067  if (dist_[idx] != nullptr)
1068  dist_[idx]->Fill(depth, edep);
1069  int jd = 2 * idx + id - 7;
1070  if (jd >= 0 && jd < 4) {
1071  jd += 5;
1072  if (hit_[jd] != nullptr)
1073  hit_[jd]->Fill(edep);
1074  if (time_[jd] != nullptr)
1075  time_[jd]->Fill(time, edep);
1076  if (dist_[jd] != nullptr)
1077  dist_[jd]->Fill(depth, edep);
1078  }
1079 }
1080 
1081 void HCalSD::plotHF(const G4ThreeVector& hitPoint, bool emType) {
1082  double zv = std::abs(hitPoint.z()) - gpar[4];
1083  if (emType) {
1084  if (hzvem != nullptr)
1085  hzvem->Fill(zv);
1086  } else {
1087  if (hzvhad != nullptr)
1088  hzvhad->Fill(zv);
1089  }
1090 }
1091 
1093  if (id.subdet == 4) {
1094  int ieta = (id.zside == 0) ? -id.etaR : id.etaR;
1095  if (hcalConstants_->maxHFDepth(ieta, id.phis) > 2) {
1096  if (id.depth <= 2) {
1097  if (G4UniformRand() > 0.5)
1098  id.depth += 2;
1099  }
1100  }
1101  } else if ((id.subdet == 1 || id.subdet == 2) && testNumber) {
1102  id.depth = (depth_ == 0) ? 1 : 2;
1103  }
1104 }
1105 
1107 #ifdef printDebug
1108  detNull_ = {0, 0, 0, 0};
1109 #endif
1110 }
1111 
1113 #ifdef printDebug
1114  int sum = detNull_[0] + detNull_[1] + detNull_[2];
1115  if (sum > 0)
1116  edm::LogVerbatim("HcalSim") << "NullDets " << detNull_[0] << " " << detNull_[1] << " " << detNull_[2] << " "
1117  << detNull_[3] << " " << (static_cast<float>(sum) / (sum + detNull_[3]));
1118 #endif
1119 }
1120 
1121 void HCalSD::printVolume(const G4VTouchable* touch) const {
1122  if (touch) {
1123 #ifdef EDM_ML_DEBUG
1124  int level = ((touch->GetHistoryDepth()) + 1);
1125  edm::LogVerbatim("CaloSimX") << "HCalSD::printVolume with " << level << " levels";
1126  static const std::string unknown("Unknown");
1127  //Get name and copy numbers
1128  for (int ii = 0; ii < level; ii++) {
1129  int i = level - ii - 1;
1130  G4VPhysicalVolume* pv = touch->GetVolume(i);
1131  G4String name = (pv != nullptr) ? pv->GetName() : unknown;
1132  G4int copyno = touch->GetReplicaNumber(i);
1133  edm::LogVerbatim("HcalSim") << "[" << ii << "] " << name << ":" << copyno;
1134  }
1135 #endif
1136  }
1137 }
float edepositEM
Definition: CaloSD.h:140
double energyCut
Definition: CaloSD.h:144
Log< level::Info, true > LogVerbatim
std::vector< std::string > matNames
Definition: HCalSD.h:121
void readWeightFromFile(const std::string &)
Definition: HCalSD.cc:975
bool isPrimary() const
std::vector< const G4LogicalVolume * > hfLV
Definition: HCalSD.h:123
bool useParam
Definition: HCalSD.h:112
double eminHitHE
Definition: HCalSD.h:113
TH1F * time_[9]
Definition: HCalSD.h:125
double kmaxNeutron
Definition: CaloSD.h:149
void hitForFibre(const G4Step *step)
Definition: HCalSD.cc:782
void initEvent(const BeginOfEvent *) override
Definition: HCalSD.cc:1106
static bool isMuon(int pdgCode)
constexpr int zside() const
get the z-side of the cell (1/-1)
Definition: HcalDetId.h:141
bool useLayerWt
Definition: HCalSD.h:109
Definition: CaloSD.h:40
std::vector< std::string > fibreNames
Definition: HCalSD.h:120
static constexpr double slopeHE_
Definition: HCalSD.h:106
std::unique_ptr< HFShowerParam > showerParam
Definition: HCalSD.h:92
double weight_
Definition: HCalSD.h:115
std::vector< double > gpar
Definition: HCalSD.h:117
bool useFibreBundle
Definition: HCalSD.h:109
constexpr int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.h:148
const HcalDDDSimConstants * hcalConstants_
Definition: HCalSD.h:96
double betaThr
Definition: HCalSD.h:111
std::vector< std::string > hfNames_
double deliveredLumi
Definition: HCalSD.h:114
int maxHFDepth(const int &ieta, const int &iphi) const
float degradation(float intlumi, int ieta, int lay) const
float *__restrict__ zv
std::vector< std::string > hfFibreNames_
void endEvent() override
Definition: HCalSD.cc:1112
double eminHitHB
Definition: HCalSD.h:113
bool useShowerLibrary
Definition: HCalSD.h:112
std::vector< std::string > hfNames
Definition: HCalSD.h:119
void plotProfile(const G4Step *step, const G4ThreeVector &pos, double edep, double time, int id)
Definition: HCalSD.cc:1019
double birk2
Definition: HCalSD.h:111
bool usePMTHit
Definition: HCalSD.h:109
double getLayer0Wt(const int &det, const int &phi, const int &zside) const
std::vector< const G4LogicalVolume * > fibre2LV
Definition: HCalSD.h:123
int zside(DetId const &)
static const unsigned int numberOfZLayers
Definition: HFDarkening.h:24
void modifyDepth(HcalNumberingFromDDD::HcalID &id)
Definition: HCalSD.cc:1092
void update(const BeginOfJob *) override
This routine will be called when the appropriate signal arrives.
Definition: HCalSD.cc:573
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:566
void processHit(const G4Step *step)
Definition: CaloSD.h:113
bool agingFlagHE
Definition: HCalSD.h:108
double birk1
Definition: HCalSD.h:111
double kmaxProton
Definition: CaloSD.h:149
const std::vector< std::string_view > logicalNames(const std::string &readoutName) const
std::unique_ptr< HFShowerPMT > showerPMT
Definition: HCalSD.h:93
const std::vector< double > & getGparHF() const
const std::string names[nVars_]
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
TH1F * hzvhad
Definition: HCalSD.h:125
bool isItConicalBundle(const G4LogicalVolume *)
Definition: HCalSD.cc:707
HCalSD(const std::string &, const HcalDDDSimConstants *, const HcalDDDRecConstants *, const HcalSimulationConstants *, const HBHEDarkening *, const HBHEDarkening *, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)
Definition: HCalSD.cc:43
double eminHitHF
Definition: HCalSD.h:113
const HcalSimulationParameters * hcalsimpar() const
T * make(const Args &...args) const
make new ROOT object
static constexpr double minRoff_
Definition: HCalSD.h:104
static uint32_t packHcalIndex(int det, int z, int depth, int eta, int phi, int lay)
std::unique_ptr< HcalNumberingFromDDD > numberingFromDDD
Definition: HCalSD.h:88
TH1F * hit_[9]
Definition: HCalSD.h:125
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > XYZVectorD
spatial vector with cartesian internal representation
Definition: Vector3D.h:8
#define EDM_ML_DEBUG
Definition: MPUnpacker.cc:1
bool agingFlagHB
Definition: HCalSD.h:108
bool testNumber
Definition: HCalSD.h:110
double kmaxIon
Definition: CaloSD.h:149
bool suppressHeavy
Definition: CaloSD.h:148
std::vector< std::string > hfFibreConicalNames_
bool useBirk
Definition: HCalSD.h:109
float edepositHAD
Definition: CaloSD.h:140
std::unique_ptr< HFDarkening > m_HFDarkening
Definition: HCalSD.h:100
double birk3
Definition: HCalSD.h:111
void resetForNewPrimary(const G4Step *)
Definition: CaloSD.cc:652
Definition: tfile.py:1
const HBHEDarkening * m_HEDarkening
Definition: HCalSD.h:99
bool applyFidCut
Definition: HCalSD.h:112
TH1F * dist_[9]
Definition: HCalSD.h:125
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void plotHF(const G4ThreeVector &pos, bool emType)
Definition: HCalSD.cc:1081
std::unique_ptr< HcalTestNS > m_HcalTestNS
Definition: HCalSD.h:101
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
int getNumberOfHits()
Definition: CaloSD.cc:461
Definition: value.py:1
std::vector< int > hfLevels
Definition: HCalSD.h:118
void getFromParam(const G4Step *step, bool &isKilled)
Definition: HCalSD.cc:824
void getFromHFLibrary(const G4Step *step, bool &isKilled)
Definition: HCalSD.cc:739
uint32_t setDetUnitId(const G4Step *step) override
Definition: HCalSD.cc:533
bool isItStraightBundle(const G4LogicalVolume *)
Definition: HCalSD.cc:699
bool testNS_
Definition: HCalSD.h:110
ii
Definition: cuy.py:589
bool getFromLibrary(const G4Step *) override
Definition: HCalSD.cc:335
double eminHitHO
Definition: HCalSD.h:113
double tmaxHit
Definition: CaloSD.h:144
static bool isStableHadronIon(const G4Track *)
void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth=0)
Definition: CaloHitID.cc:41
TH1F * hzvem
Definition: HCalSD.h:125
std::vector< const G4LogicalVolume * > pmtLV
Definition: HCalSD.h:123
int PMTNumber(const G4ThreeVector &pe_effect)
Definition: DetId.h:17
static const unsigned int lowZLimit
Definition: HFDarkening.h:27
bool filterHit(CaloG4Hit *, double) override
Definition: HCalSD.cc:577
CaloHitID currentID
Definition: CaloSD.h:142
static const unsigned int upperZLimit
Definition: HFDarkening.h:28
bool neutralDensity
Definition: HCalSD.h:110
bool isHF
Definition: HCalSD.h:107
void fillLogVolumeVector(const std::string &, const std::vector< std::string > &, std::vector< const G4LogicalVolume *> &)
Definition: HCalSD.cc:310
static constexpr double maxZ_
Definition: HCalSD.h:103
void getHitFibreBundle(const G4Step *step, bool type)
Definition: HCalSD.cc:913
std::vector< std::string > hcalMaterialNames_
double getEnergyDeposit(const G4Step *) override
Definition: HCalSD.cc:399
double getAttenuation(const G4Step *aStep, double birk1, double birk2, double birk3) const
Definition: CaloSD.cc:665
double getEnergyDeposit() const
Definition: CaloG4Hit.h:79
std::unique_ptr< HFShowerFibreBundle > showerBundle
Definition: HCalSD.h:94
std::vector< std::string > hfPMTNames_
bool isItPMT(const G4LogicalVolume *)
Definition: HCalSD.cc:691
bool isItinFidVolume(const G4ThreeVector &)
Definition: HCalSD.cc:723
uint32_t getUnitID() const
Definition: CaloG4Hit.h:66
bool isHF(int etabin, int depth)
void printVolume(const G4VTouchable *touch) const
Definition: HCalSD.cc:1121
void getHitPMT(const G4Step *step)
Definition: HCalSD.cc:854
std::vector< int > detNull_
Definition: HCalSD.h:126
std::unique_ptr< HFShowerLibrary > showerLibrary
Definition: HCalSD.h:90
double layerWeight(int, const G4ThreeVector &, int, int)
Definition: HCalSD.cc:999
const HcalSimulationConstants * hcalSimConstants_
Definition: HCalSD.h:97
int depth_
Definition: HCalSD.h:116
HLT enums.
virtual int setTrackID(const G4Step *)
Definition: CaloSD.cc:825
std::map< uint32_t, double > layerWeights
Definition: HCalSD.h:124
bool isItHF(const G4Step *)
Definition: HCalSD.cc:654
static bool isGammaElectronPositron(int pdgCode)
bool isItScintillator(const G4Material *)
Definition: HCalSD.cc:715
std::vector< std::string > hfFibreStraightNames_
std::unique_ptr< HcalNumberingScheme > numberingScheme
Definition: HCalSD.h:89
std::unique_ptr< HFShower > hfshower
Definition: HCalSD.h:91
std::vector< const G4LogicalVolume * > fibreLV
Definition: HCalSD.h:123
Log< level::Warning, false > LogWarning
double getResponseWt(const G4Track *)
Definition: CaloSD.cc:860
std::vector< const G4LogicalVolume * > fibre1LV
Definition: HCalSD.h:123
tmp
align.sh
Definition: createJobs.py:716
std::vector< const G4Material * > materials
Definition: HCalSD.h:122
int getMaxDepth(const int &type) const
void initRun() override
Definition: HCalSD.cc:575
const HBHEDarkening * m_HBDarkening
Definition: HCalSD.h:98
constexpr int iphi() const
get the cell iphi
Definition: HcalDetId.h:157
static constexpr double maxRoff_
Definition: HCalSD.h:105
bool useHF
Definition: HCalSD.h:112
bool isItFibre(const G4LogicalVolume *)
Definition: HCalSD.cc:675
void setParameterized(bool val)
Definition: CaloSD.h:110
constexpr int depth() const
get the tower depth
Definition: HcalDetId.h:164