CMS 3D CMS Logo

HCalSD.cc
Go to the documentation of this file.
1 // File: HCalSD.cc
3 // Description: Sensitive Detector class for calorimeters
5 
17 
26 
27 #include "G4LogicalVolumeStore.hh"
28 #include "G4LogicalVolume.hh"
29 #include "G4Step.hh"
30 #include "G4Track.hh"
31 #include "G4ParticleTable.hh"
32 #include "G4VProcess.hh"
33 
34 #include "G4SystemOfUnits.hh"
35 #include "G4PhysicalConstants.hh"
36 
37 #include <iostream>
38 #include <fstream>
39 #include <iomanip>
40 
41 //#define EDM_ML_DEBUG
42 //#define plotDebug
43 
45  const SensitiveDetectorCatalog & clg,
46  edm::ParameterSet const & p, const SimTrackManager* manager) :
47  CaloSD(name, cpv, clg, p, manager,
48  (float)(p.getParameter<edm::ParameterSet>("HCalSD").getParameter<double>("TimeSliceUnit")),
49  p.getParameter<edm::ParameterSet>("HCalSD").getParameter<bool>("IgnoreTrackID")),
50  hcalConstants(nullptr), numberingFromDDD(nullptr), numberingScheme(nullptr),
51  showerLibrary(nullptr), hfshower(nullptr), showerParam(nullptr), showerPMT(nullptr),
52  showerBundle(nullptr), m_HBDarkening(nullptr), m_HEDarkening(nullptr),
53  m_HFDarkening(nullptr), hcalTestNS_(nullptr), depth_(1) {
54 
55  //static SimpleConfigurable<double> bk1(0.013, "HCalSD:BirkC1");
56  //static SimpleConfigurable<double> bk2(0.0568,"HCalSD:BirkC2");
57  //static SimpleConfigurable<double> bk3(1.75, "HCalSD:BirkC3");
58  // Values from NIM 80 (1970) 239-244: as implemented in Geant3
59 
61  useBirk = m_HC.getParameter<bool>("UseBirkLaw");
62  birk1 = m_HC.getParameter<double>("BirkC1")*(g/(MeV*cm2));
63  birk2 = m_HC.getParameter<double>("BirkC2");
64  birk3 = m_HC.getParameter<double>("BirkC3");
65  useShowerLibrary = m_HC.getParameter<bool>("UseShowerLibrary");
66  useParam = m_HC.getParameter<bool>("UseParametrize");
67  testNumber = m_HC.getParameter<bool>("TestNumberingScheme");
68  neutralDensity = m_HC.getParameter<bool>("doNeutralDensityFilter");
69  usePMTHit = m_HC.getParameter<bool>("UsePMTHits");
70  betaThr = m_HC.getParameter<double>("BetaThreshold");
71  eminHitHB = m_HC.getParameter<double>("EminHitHB")*MeV;
72  eminHitHE = m_HC.getParameter<double>("EminHitHE")*MeV;
73  eminHitHO = m_HC.getParameter<double>("EminHitHO")*MeV;
74  eminHitHF = m_HC.getParameter<double>("EminHitHF")*MeV;
75  useFibreBundle = m_HC.getParameter<bool>("UseFibreBundleHits");
76  deliveredLumi = m_HC.getParameter<double>("DelivLuminosity");
77  agingFlagHB = m_HC.getParameter<bool>("HBDarkening");
78  agingFlagHE = m_HC.getParameter<bool>("HEDarkening");
79  bool agingFlagHF = m_HC.getParameter<bool>("HFDarkening");
80  useHF = m_HC.getUntrackedParameter<bool>("UseHF",true);
81  bool forTBH2 = m_HC.getUntrackedParameter<bool>("ForTBH2",false);
82  useLayerWt = m_HC.getUntrackedParameter<bool>("UseLayerWt",false);
83  std::string file = m_HC.getUntrackedParameter<std::string>("WtFile","None");
84  testNS_ = m_HC.getUntrackedParameter<bool>("TestNS",false);
85  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
86  applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
87 
88 #ifdef EDM_ML_DEBUG
89  LogDebug("HcalSim") << "***************************************************"
90  << "\n"
91  << "* *"
92  << "\n"
93  << "* Constructing a HCalSD with name " << name << "\n"
94  << "* *"
95  << "\n"
96  << "***************************************************";
97 #endif
98  edm::LogInfo("HcalSim") << "HCalSD:: Use of HF code is set to " << useHF
99  << "\nUse of shower parametrization set to "
100  << useParam << "\nUse of shower library is set to "
101  << useShowerLibrary << "\nUse PMT Hit is set to "
102  << usePMTHit << " with beta Threshold "<< betaThr
103  << "\nUSe of FibreBundle Hit set to "<<useFibreBundle
104  << "\n Use of Birks law is set to "
105  << useBirk << " with three constants kB = "
106  << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3;
107  edm::LogInfo("HcalSim") << "HCalSD:: Suppression Flag " << suppressHeavy
108  << " protons below " << kmaxProton << " MeV,"
109  << " neutrons below " << kmaxNeutron << " MeV and"
110  << " ions below " << kmaxIon << " MeV\n"
111  << " Threshold for storing hits in HB: "
112  << eminHitHB << " HE: " << eminHitHE << " HO: "
113  << eminHitHO << " HF: " << eminHitHF << "\n"
114  << "Delivered luminosity for Darkening "
115  << deliveredLumi << " Flag (HE) " << agingFlagHE
116  << " Flag (HB) " << agingFlagHB
117  << " Flag (HF) " << agingFlagHF << "\n"
118  << "Application of Fiducial Cut " << applyFidCut
119  << "Flag for test number|neutral density filter "
120  << testNumber << " " << neutralDensity;
121 
122  HcalNumberingScheme* scheme;
123  if (testNumber || forTBH2)
124  scheme = dynamic_cast<HcalNumberingScheme*>(new HcalTestNumberingScheme(forTBH2));
125  else
126  scheme = new HcalNumberingScheme();
127  setNumberingScheme(scheme);
128 
129  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
130  std::vector<G4LogicalVolume *>::const_iterator lvcite;
131  const G4LogicalVolume* lv;
132  std::string attribute, value;
133  if (useHF) {
134  if (useParam) {
135  showerParam = new HFShowerParam(name, cpv, p);
136  } else {
137  if (useShowerLibrary) showerLibrary = new HFShowerLibrary(name, cpv, p);
138  hfshower = new HFShower(name, cpv, p, 0);
139  }
140 
141  // HF volume names
142  attribute = "Volume";
143  value = "HF";
144  DDSpecificsMatchesValueFilter filter0{DDValue(attribute,value,0)};
145  DDFilteredView fv0(cpv,filter0);
146  hfNames = getNames(fv0);
147  fv0.firstChild();
148  DDsvalues_type sv0(fv0.mergedSpecifics());
149  std::vector<double> temp = getDDDArray("Levels",sv0);
150  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
151  << " = " << value << " has " << hfNames.size()
152  << " elements";
153  for (unsigned int i=0; i < hfNames.size(); ++i) {
154  G4String namv = hfNames[i];
155  lv = nullptr;
156  for(lvcite=lvs->begin(); lvcite!=lvs->end(); lvcite++)
157  if((*lvcite)->GetName()==namv) {
158  lv = (*lvcite);
159  break;
160  }
161  hfLV.push_back(lv);
162  int level = static_cast<int>(temp[i]);
163  hfLevels.push_back(level);
164  edm::LogInfo("HcalSim") << "HCalSD: HF[" << i << "] = " << hfNames[i]
165  << " LV " << hfLV[i] << " at level "
166  << hfLevels[i];
167  }
168 
169  // HF Fibre volume names
170  value = "HFFibre";
171  DDSpecificsMatchesValueFilter filter1{DDValue(attribute,value,0)};
172  DDFilteredView fv1(cpv,filter1);
173  fibreNames = getNames(fv1);
174  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
175  << " = " << value << ":";
176  for (unsigned int i=0; i<fibreNames.size(); ++i) {
177  G4String namv = fibreNames[i];
178  lv = nullptr;
179  for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite) {
180  if ((*lvcite)->GetName() == namv) {
181  lv = (*lvcite);
182  break;
183  }
184  }
185  fibreLV.push_back(lv);
186  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
187  << " LV " << fibreLV[i];
188  }
189 
190  // HF PMT volume names
191  value = "HFPMT";
192  DDSpecificsMatchesValueFilter filter3{DDValue(attribute,value,0)};
193  DDFilteredView fv3(cpv,filter3);
194  std::vector<G4String> pmtNames = getNames(fv3);
195  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
196  << " = " << value << " have " << pmtNames.size()
197  << " entries";
198  for (unsigned int i=0; i<pmtNames.size(); ++i) {
199  G4String namv = pmtNames[i];
200  lv = nullptr;
201  for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite)
202  if ((*lvcite)->GetName() == namv) {
203  lv = (*lvcite);
204  break;
205  }
206  pmtLV.push_back(lv);
207  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << pmtNames[i]
208  << " LV " << pmtLV[i];
209  }
210  if (!pmtNames.empty()) showerPMT = new HFShowerPMT (name, cpv, p);
211 
212  // HF Fibre bundles
213  value = "HFFibreBundleStraight";
214  DDSpecificsMatchesValueFilter filter4{DDValue(attribute,value,0)};
215  DDFilteredView fv4(cpv,filter4);
216  std::vector<G4String> fibreNames = getNames(fv4);
217  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
218  << " = " << value << " have " << fibreNames.size()
219  << " entries";
220  for (unsigned int i=0; i<fibreNames.size(); ++i) {
221  G4String namv = fibreNames[i];
222  lv = nullptr;
223  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
224  if ((*lvcite)->GetName() == namv) {
225  lv = (*lvcite);
226  break;
227  }
228  fibre1LV.push_back(lv);
229  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
230  << " LV " << fibre1LV[i];
231  }
232 
233  // Geometry parameters for HF
234  value = "HFFibreBundleConical";
235  DDSpecificsMatchesValueFilter filter5{DDValue(attribute,value,0)};
236  DDFilteredView fv5(cpv,filter5);
237  fibreNames = getNames(fv5);
238  edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
239  << " = " << value << " have " << fibreNames.size()
240  << " entries";
241  for (unsigned int i=0; i<fibreNames.size(); ++i) {
242  G4String namv = fibreNames[i];
243  lv = nullptr;
244  for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite)
245  if ((*lvcite)->GetName() == namv) {
246  lv = (*lvcite);
247  break;
248  }
249  fibre2LV.push_back(lv);
250  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
251  << " LV " << fibre2LV[i];
252  }
253  if (!fibre1LV.empty() || !fibre2LV.empty())
254  showerBundle = new HFShowerFibreBundle (name, cpv, p);
255  }
256 
257  //Material list for HB/HE/HO sensitive detectors
258  const G4MaterialTable * matTab = G4Material::GetMaterialTable();
259  std::vector<G4Material*>::const_iterator matite;
260  attribute = "OnlyForHcalSimNumbering";
261  DDSpecificsHasNamedValueFilter filter2{attribute};
262  DDFilteredView fv2(cpv,filter2);
263  bool dodet = fv2.firstChild();
265 
266  while (dodet) {
267  const DDLogicalPart & log = fv2.logicalPart();
268  G4String namx = log.name().name();
269  if (!isItHF(namx) && !isItFibre(namx)) {
270  bool notIn = true;
271  for (unsigned int i=0; i<matNames.size(); ++i) {
272  if (!strcmp(matNames[i].c_str(),log.material().name().name().c_str())){
273  notIn = false;
274  break;
275  }
276  }
277  if (notIn) {
278  namx = log.material().name().name();
279  matNames.push_back(namx);
280  const G4Material* mat = nullptr;
281  for (matite = matTab->begin(); matite != matTab->end(); ++matite) {
282  if ((*matite)->GetName() == namx) {
283  mat = (*matite);
284  break;
285  }
286  }
287  materials.push_back(mat);
288  }
289  }
290  dodet = fv2.next();
291  }
292 
293  edm::LogInfo("HcalSim") << "HCalSD: Material names for " << attribute
294  << " = " << name << ":";
295  for (unsigned int i=0; i<matNames.size(); ++i)
296  edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << matNames[i]
297  << " pointer " << materials[i];
298 
299  mumPDG = mupPDG = 0;
300 
301  if (useLayerWt) readWeightFromFile(file);
302 
303  for (int i=0; i<9; ++i) hit_[i] = time_[i]= dist_[i] = nullptr;
304  hzvem = hzvhad = nullptr;
305 
306  if (agingFlagHF) m_HFDarkening.reset(new HFDarkening(m_HC.getParameter<edm::ParameterSet>("HFDarkeningParameterBlock")));
307 #ifdef plotDebug
309 
310  if ( tfile.isAvailable() ) {
311  static const char * const labels[] = {"HB", "HE", "HO", "HF Absorber", "HF PMT",
312  "HF Absorber Long", "HF Absorber Short",
313  "HF PMT Long", "HF PMT Short"};
314  TFileDirectory hcDir = tfile->mkdir("ProfileFromHCalSD");
315  char name[20], title[60];
316  for (int i=0; i<9; ++i) {
317  sprintf (title, "Hit energy in %s", labels[i]);
318  sprintf (name, "HCalSDHit%d", i);
319  hit_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
320  sprintf (title, "Energy (MeV)");
321  hit_[i]->GetXaxis()->SetTitle(title);
322  hit_[i]->GetYaxis()->SetTitle("Hits");
323  sprintf (title, "Time of the hit in %s", labels[i]);
324  sprintf (name, "HCalSDTime%d", i);
325  time_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
326  sprintf (title, "Time (ns)");
327  time_[i]->GetXaxis()->SetTitle(title);
328  time_[i]->GetYaxis()->SetTitle("Hits");
329  sprintf (title, "Longitudinal profile in %s", labels[i]);
330  sprintf (name, "HCalSDDist%d", i);
331  dist_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
332  sprintf (title, "Distance (mm)");
333  dist_[i]->GetXaxis()->SetTitle(title);
334  dist_[i]->GetYaxis()->SetTitle("Hits");
335  }
336  if (useHF && (!useParam)) {
337  hzvem = hcDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part)",330,0.0,1650.0);
338  hzvem->GetXaxis()->SetTitle("Longitudinal Profile (EM Part)");
339  hzvhad = hcDir.make<TH1F>("hzvhad","Longitudinal Profile (Had Part)",330,0.0,1650.0);
340  hzvhad->GetXaxis()->SetTitle("Longitudinal Profile (Hadronic Part)");
341  }
342  }
343 #endif
344 }
345 
347 
349  if (numberingScheme) delete numberingScheme;
350  if (showerLibrary) delete showerLibrary;
351  if (hfshower) delete hfshower;
352  if (showerParam) delete showerParam;
353  if (showerPMT) delete showerPMT;
354  if (showerBundle) delete showerBundle;
355  if (hcalTestNS_) delete hcalTestNS_;
356 }
357 
358 bool HCalSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
359 
360  NaNTrap( aStep ) ;
361 
362  if (aStep == nullptr) {
363  return true;
364  } else {
365  depth_ = (aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(0))%10;
366  const G4LogicalVolume* lv =
367  aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
368  G4String nameVolume = lv->GetName();
369  if (isItHF(aStep)) {
370  G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
371  double weight(1.0);
372  if (m_HFDarkening) {
373  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
374  double r = hitPoint.perp()/CLHEP::cm;
375  double z = std::abs(hitPoint.z())/CLHEP::cm;
376  double dose_acquired = 0.;
378  unsigned int hfZLayer = (int)((z - HFDarkening::lowZLimit)/5);
379  if (hfZLayer >= HFDarkening::upperZLimit) hfZLayer = (HFDarkening::upperZLimit-1);
380  float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
381  for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
382  dose_acquired = m_HFDarkening->dose(i,r);
383  weight *= m_HFDarkening->degradation(normalized_lumi*dose_acquired);
384  }
385  }
386 #ifdef EDM_ML_DEBUG
387  LogDebug("HcalSim") << "HCalSD: HFLumiDarkening at r = " << r
388  << ", z = " << z << " Dose " << dose_acquired
389  << " weight " << weight;
390 #endif
391  }
392  if (useParam) {
393 #ifdef EDM_ML_DEBUG
394  LogDebug("HcalSim") << "HCalSD: " << getNumberOfHits()
395  << " hits from parametrization in " << nameVolume
396  << " for Track " << aStep->GetTrack()->GetTrackID()
397  <<" (" << aStep->GetTrack()->GetDefinition()->GetParticleName()
398  <<")";
399 #endif
400  getFromParam(aStep, weight);
401 #ifdef EDM_ML_DEBUG
402  LogDebug("HcalSim") << "HCalSD: " << getNumberOfHits()
403  << " hits afterParamS*";
404 #endif
405  } else {
406  bool notaMuon = true;
407  if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
408  if (useShowerLibrary && notaMuon) {
409 #ifdef EDM_ML_DEBUG
410  LogDebug("HcalSim") << "HCalSD: Starts shower library from "
411  << nameVolume << " for Track "
412  << aStep->GetTrack()->GetTrackID() << " ("
413  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
414 #endif
415  getFromLibrary(aStep, weight);
416  } else if (isItFibre(lv)) {
417 #ifdef EDM_ML_DEBUG
418  LogDebug("HcalSim") << "HCalSD: Hit at Fibre in " << nameVolume
419  << " for Track "
420  << aStep->GetTrack()->GetTrackID() <<" ("
421  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
422 #endif
423  hitForFibre(aStep, weight);
424  }
425  }
426  } else if (isItPMT(lv)) {
427 #ifdef EDM_ML_DEBUG
428  LogDebug("HcalSim") << "HCalSD: Hit from PMT parametrization from "
429  << nameVolume << " for Track "
430  << aStep->GetTrack()->GetTrackID() << " ("
431  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
432 #endif
433  if (usePMTHit && showerPMT) getHitPMT(aStep);
434  } else if (isItStraightBundle(lv) || isItConicalBundle(lv)) {
435 #ifdef EDM_ML_DEBUG
436  LogDebug("HcalSim") << "HCalSD: Hit from FibreBundle from "
437  << nameVolume << " for Track "
438  << aStep->GetTrack()->GetTrackID() << " ("
439  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
440 #endif
443  } else {
444 #ifdef EDM_ML_DEBUG
445  LogDebug("HcalSim") << "HCalSD: Hit from standard path from "
446  << nameVolume << " for Track "
447  << aStep->GetTrack()->GetTrackID() << " ("
448  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
449 #endif
450  if (getStepInfo(aStep)) {
451 #ifdef plotDebug
452  if (edepositEM+edepositHAD > 0)
453  plotProfile(aStep, aStep->GetPreStepPoint()->GetPosition(),
454  edepositEM+edepositHAD,aStep->GetPostStepPoint()->GetGlobalTime(),0);
455 #endif
456  if (hitExists() == false && edepositEM+edepositHAD>0.) currentHit = createNewHit();
457  }
458  }
459  return true;
460  }
461 }
462 
463 double HCalSD::getEnergyDeposit(G4Step* aStep) {
464  double destep = aStep->GetTotalEnergyDeposit();
465  double weight = 1;
466  theTrack = aStep->GetTrack();
467 
468  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
469  uint32_t detid = setDetUnitId(aStep);
470  int det(0), ieta(0), phi(0), z(0), lay, depth(-1);
471  if (testNumber) {
472  HcalTestNumbering::unpackHcalIndex(detid,det,z,depth,ieta,phi,lay);
473  if (z==0) z = -1;
474  } else {
475  HcalDetId hcid(detid);
476  det = hcid.subdetId();
477  ieta = hcid.ietaAbs();
478  phi = hcid.iphi();
479  z = hcid.zside();
480  }
481  lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
482 #ifdef EDM_ML_DEBUG
483  edm::LogInfo("HcalSim") << "HCalSD: det: " << det << " ieta: "<< ieta
484  << " iphi: " << phi << " zside " << z << " lay: "
485  << lay-2;
486 #endif
487  if (depth_==0 && (det==1 || det==2) && ((!testNumber) || neutralDensity))
488  weight = hcalConstants->getLayer0Wt(det,phi,z);
489  if (useLayerWt) {
490  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
491  weight = layerWeight(det+2, hitPoint, depth_, lay);
492  }
493 
494  if (m_HBDarkening && det == 1) {
495  float dweight = m_HBDarkening->degradation(deliveredLumi,ieta,lay);
496  weight *= dweight;
497 #ifdef EDM_ML_DEBUG
498  edm::LogInfo("HcalSim") << "HCalSD: >>> HB Lumi: " << deliveredLumi
499  << " coefficient = " << dweight;
500 #endif
501  }
502 
503  if (m_HEDarkening && det == 2) {
504  float dweight = m_HEDarkening->degradation(deliveredLumi,ieta,lay);
505  weight *= dweight;
506 #ifdef EDM_ML_DEBUG
507  edm::LogInfo("HcalSim") << "HCalSD: >>> HE Lumi: " << deliveredLumi
508  << " coefficient = " << dweight;
509 #endif
510  }
511 
512  if (suppressHeavy) {
513  TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
514  if (trkInfo) {
515  int pdg = theTrack->GetDefinition()->GetPDGEncoding();
516  if (!(trkInfo->isPrimary())) { // Only secondary particles
517  double ke = theTrack->GetKineticEnergy()/MeV;
518  if ( pdg/1000000000 == 1 && (pdg/10000)%100 > 0 &&
519  (pdg/10)%100 > 0 && ke <kmaxIon ) weight = 0;
520  if ((pdg == 2212) && (ke < kmaxProton)) weight = 0;
521  if ((pdg == 2112) && (ke < kmaxNeutron)) weight = 0;
522 #ifdef EDM_ML_DEBUG
523  if (weight == 0)
524  edm::LogInfo("HcalSim") << "HCalSD:Ignore Track "
525  << theTrack->GetTrackID() << " Type "
526  << theTrack->GetDefinition()->GetParticleName()
527  << " Kinetic Energy " << ke << " MeV";
528 #endif
529  }
530  }
531  }
532 #ifdef EDM_ML_DEBUG
533  double weight0 = weight;
534 #endif
535  if (useBirk) {
536  const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
537  if (isItScintillator(mat)) weight *= getAttenuation(aStep, birk1, birk2, birk3);
538  }
539  double wt1 = getResponseWt(theTrack);
540  double wt2 = theTrack->GetWeight();
541 #ifdef EDM_ML_DEBUG
542  edm::LogInfo("HcalSim") << "HCalSD: Detector " << det+2 << " Depth " << depth_
543  << " weight " << weight0 << " " << weight << " " << wt1
544  << " " << wt2;
545 #endif
546  double edep = weight*wt1*destep;
547  if (wt2 > 0.0) { edep *= wt2; }
548  return edep;
549 }
550 
551 uint32_t HCalSD::setDetUnitId(const G4Step * aStep) {
552 
553  const G4StepPoint* prePoint = aStep->GetPreStepPoint();
554  const G4VTouchable* touch = prePoint->GetTouchable();
555  const G4ThreeVector& hitPoint = prePoint->GetPosition();
556 
557  int depth = (touch->GetReplicaNumber(0))%10 + 1;
558  int lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
559  int det = (touch->GetReplicaNumber(1))/1000;
560 
561  return setDetUnitId (det, hitPoint, depth, lay);
562 }
563 
565  if (scheme != nullptr) {
566  edm::LogInfo("HcalSim") << "HCalSD: updates numbering scheme for " << GetName();
567  if (numberingScheme) delete numberingScheme;
568  numberingScheme = scheme;
569  }
570 }
571 
572 void HCalSD::update(const BeginOfJob * job) {
573 
574  const edm::EventSetup* es = (*job)();
576  es->get<HcalSimNumberingRecord>().get(hdc);
577  if (hdc.isValid()) {
578  hcalConstants = (HcalDDDSimConstants*)(&(*hdc));
579  } else {
580  edm::LogError("HcalSim") << "HCalSD : Cannot find HcalDDDSimConstant";
581  throw cms::Exception("Unknown", "HCalSD") << "Cannot find HcalDDDSimConstant" << "\n";
582  }
583 
585 
586  edm::LogInfo("HcalSim") << "Maximum depth for HF " << hcalConstants->getMaxDepth(2);
587 
588  //Special Geometry parameters
590  edm::LogInfo("HcalSim") << "HCalSD: " << gpar.size()<< " gpar (cm)";
591  for (unsigned int ig=0; ig<gpar.size(); ig++)
592  edm::LogInfo("HcalSim") << "HCalSD: gpar[" << ig << "] = "
593  << gpar[ig]/cm << " cm";
594  //Test Hcal Numbering Scheme
595  if (testNS_) hcalTestNS_ = new HcalTestNS(es);
596 
597  if (agingFlagHB) {
599  es->get<HBHEDarkeningRecord>().get("HB",hdark);
600  m_HBDarkening = &*hdark;
601  }
602  if (agingFlagHE) {
604  es->get<HBHEDarkeningRecord>().get("HE",hdark);
605  m_HEDarkening = &*hdark;
606  }
607 }
608 
610  G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
611  G4String particleName;
612  mumPDG = theParticleTable->FindParticle(particleName="mu-")->GetPDGEncoding();
613  mupPDG = theParticleTable->FindParticle(particleName="mu+")->GetPDGEncoding();
614 #ifdef EDM_ML_DEBUG
615  LogDebug("HcalSim") << "HCalSD: Particle code for mu- = " << mumPDG
616  << " for mu+ = " << mupPDG;
617 #endif
618  if (showerLibrary) showerLibrary->initRun(theParticleTable,hcalConstants);
619  if (showerParam) showerParam->initRun(theParticleTable,hcalConstants);
620  if (hfshower) hfshower->initRun(theParticleTable,hcalConstants);
621  if (showerPMT) showerPMT->initRun(theParticleTable,hcalConstants);
622  if (showerBundle) showerBundle->initRun(theParticleTable,hcalConstants);
623 }
624 
625 bool HCalSD::filterHit(CaloG4Hit* aHit, double time) {
626  double threshold=0;
627  DetId theId(aHit->getUnitID());
628  switch (theId.subdetId()) {
629  case HcalBarrel:
630  threshold = eminHitHB; break;
631  case HcalEndcap:
632  threshold = eminHitHE; break;
633  case HcalOuter:
634  threshold = eminHitHO; break;
635  case HcalForward:
636  threshold = eminHitHF; break;
637  default:
638  break;
639  }
640  return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > threshold));
641 }
642 
643 uint32_t HCalSD::setDetUnitId (int det, const G4ThreeVector& pos, int depth, int lay=1) {
644  uint32_t id = 0;
645  if (numberingFromDDD) {
646  //get the ID's as eta, phi, depth, ... indices
647  HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, pos, depth, lay);
648  id = setDetUnitId(tmp);
649  }
650  return id;
651 }
652 
654  modifyDepth(tmp);
655  uint32_t id = (numberingScheme) ? numberingScheme->getUnitID(tmp) : 0;
656  if ((!testNumber) && hcalTestNS_) {
657  bool ok = hcalTestNS_->compare(tmp,id);
658  if (!ok)
659  edm::LogWarning("HcalSim") << "Det ID from HCalSD " << HcalDetId(id)
660  << " " << std::hex << id << std::dec
661  << " does not match one from relabller for "
662  << tmp.subdet << ":" << tmp.etaR << ":"
663  << tmp.phi << ":" << tmp.phis << ":"
664  << tmp.depth << ":" << tmp.lay << std::endl;
665  }
666  return id;
667 }
668 
669 std::vector<double> HCalSD::getDDDArray(const std::string & str,
670  const DDsvalues_type & sv) {
671 #ifdef EDM_ML_DEBUG
672  LogDebug("HcalSim") << "HCalSD:getDDDArray called for " << str;
673 #endif
674  DDValue value(str);
675  if (DDfetch(&sv,value)) {
676 #ifdef EDM_ML_DEBUG
677  LogDebug("HcalSim") << value;
678 #endif
679  const std::vector<double> & fvec = value.doubles();
680  int nval = fvec.size();
681  if (nval < 1) {
682  edm::LogError("HcalSim") << "HCalSD : # of " << str << " bins " << nval
683  << " < 2 ==> illegal";
684  throw cms::Exception("Unknown", "HCalSD") << "nval < 2 for array " << str << "\n";
685  }
686 
687  return fvec;
688  } else {
689  edm::LogError("HcalSim") << "HCalSD : cannot get array " << str;
690  throw cms::Exception("Unknown", "HCalSD") << "cannot get array " << str << "\n";
691  }
692 }
693 
694 std::vector<G4String> HCalSD::getNames(DDFilteredView& fv) {
695 
696  std::vector<G4String> tmp;
697  bool dodet = fv.firstChild();
698  while (dodet) {
699  const DDLogicalPart & log = fv.logicalPart();
700  bool ok = true;
701 
702  for (unsigned int i=0; i<tmp.size(); ++i) {
703  if (!strcmp(tmp[i].c_str(), log.name().name().c_str())) {
704  ok = false;
705  break;
706  }
707  }
708  if (ok) tmp.push_back(log.name().name());
709  dodet = fv.next();
710  }
711  return tmp;
712 }
713 
714 bool HCalSD::isItHF(const G4Step * aStep) {
715  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
716  int levels = (touch->GetHistoryDepth()) + 1;
717  for (unsigned int it=0; it < hfNames.size(); ++it) {
718  if (levels >= hfLevels[it]) {
719  const G4LogicalVolume* lv = touch->GetVolume(levels-hfLevels[it])->GetLogicalVolume();
720  if (lv == hfLV[it]) return true;
721  }
722  }
723  return false;
724 }
725 
726 bool HCalSD::isItHF (const G4String& name) {
727  std::vector<G4String>::const_iterator it = hfNames.begin();
728  for (; it != hfNames.end(); ++it) if (name == *it) return true;
729  return false;
730 }
731 
732 bool HCalSD::isItFibre (const G4LogicalVolume* lv) {
733  std::vector<const G4LogicalVolume*>::const_iterator ite = fibreLV.begin();
734  for (; ite != fibreLV.end(); ++ite) if (lv == *ite) return true;
735  return false;
736 }
737 
738 bool HCalSD::isItFibre (const G4String& name) {
739  std::vector<G4String>::const_iterator it = fibreNames.begin();
740  for (; it != fibreNames.end(); ++it) if (name == *it) return true;
741  return false;
742 }
743 
744 bool HCalSD::isItPMT (const G4LogicalVolume* lv) {
745  std::vector<const G4LogicalVolume*>::const_iterator ite = pmtLV.begin();
746  for (; ite != pmtLV.end(); ++ite) if (lv == *ite) return true;
747  return false;
748 }
749 
750 bool HCalSD::isItStraightBundle (const G4LogicalVolume* lv) {
751  std::vector<const G4LogicalVolume*>::const_iterator ite = fibre1LV.begin();
752  for (; ite != fibre1LV.end(); ++ite) if (lv == *ite) return true;
753  return false;
754 }
755 
756 bool HCalSD::isItConicalBundle (const G4LogicalVolume* lv) {
757  std::vector<const G4LogicalVolume*>::const_iterator ite = fibre2LV.begin();
758  for (; ite != fibre2LV.end(); ++ite) if (lv == *ite) return true;
759  return false;
760 }
761 
762 bool HCalSD::isItScintillator (const G4Material* mat) {
763  std::vector<const G4Material*>::const_iterator ite = materials.begin();
764  for (; ite != materials.end(); ++ite) if (mat == *ite) return true;
765  return false;
766 }
767 
768 bool HCalSD::isItinFidVolume (const G4ThreeVector& hitPoint) {
769  bool flag = true;
770  if (applyFidCut) {
771  int npmt = HFFibreFiducial:: PMTNumber(hitPoint);
772 #ifdef EDM_ML_DEBUG
773  edm::LogInfo("HcalSim") << "HCalSD::isItinFidVolume:#PMT= " << npmt
774  << " for hit point " << hitPoint;
775 #endif
776  if (npmt <= 0) flag = false;
777  }
778 #ifdef EDM_ML_DEBUG
779  edm::LogInfo("HcalSim") << "HCalSD::isItinFidVolume: point " << hitPoint
780  << " return flag " << flag;
781 #endif
782  return flag;
783 }
784 
785 void HCalSD::getFromLibrary (G4Step* aStep, double weight) {
786  preStepPoint = aStep->GetPreStepPoint();
787  theTrack = aStep->GetTrack();
788  int det = 5;
789  bool ok;
790 
791  std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, ok, weight, false);
792 
793  double etrack = preStepPoint->GetKineticEnergy();
794  int primaryID = setTrackID(aStep);
795 
796  // Reset entry point for new primary
797  posGlobal = preStepPoint->GetPosition();
798  resetForNewPrimary(posGlobal, etrack);
799 
800  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
801  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
802  edepositEM = 1.*GeV;
803  edepositHAD = 0.;
804  } else {
805  edepositEM = 0.;
806  edepositHAD = 1.*GeV;
807  }
808 #ifdef EDM_ML_DEBUG
809  edm::LogInfo("HcalSim") << "HCalSD::getFromLibrary " <<hits.size()
810  << " hits for " << GetName() << " of " << primaryID
811  << " with " << theTrack->GetDefinition()->GetParticleName()
812  << " of " << preStepPoint->GetKineticEnergy()/GeV << " GeV";
813 #endif
814  for (unsigned int i=0; i<hits.size(); ++i) {
815  G4ThreeVector hitPoint = hits[i].position;
816  if (isItinFidVolume (hitPoint)) {
817  int depth = hits[i].depth;
818  double time = hits[i].time;
819  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
820  currentID.setID(unitID, time, primaryID, 0);
821 #ifdef plotDebug
822  plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
823  bool emType = false;
824  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
825  emType = true;
826  plotHF(hitPoint,emType);
827 #endif
828 
829  // check if it is in the same unit and timeslice as the previous one
830  if (currentID == previousID) {
832  } else {
833  if (!checkHit()) currentHit = createNewHit();
834  }
835  }
836  }
837 
838  //Now kill the current track
839  if (ok) {
840  theTrack->SetTrackStatus(fStopAndKill);
841  G4TrackVector tv = *(aStep->GetSecondary());
842  for (unsigned int kk=0; kk<tv.size(); ++kk)
843  if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume())
844  tv[kk]->SetTrackStatus(fStopAndKill);
845  }
846 }
847 
848 void HCalSD::hitForFibre (const G4Step* aStep, double weight) { // if not ParamShower
849 
850  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
851  const G4Track* theTrack = aStep->GetTrack();
852  int primaryID = setTrackID(aStep);
853 
854  int det = 5;
855  std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight);
856 
857  G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
858  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
859  edepositEM = 1.*GeV;
860  edepositHAD = 0.;
861  } else {
862  edepositEM = 0.;
863  edepositHAD = 1.*GeV;
864  }
865 
866 #ifdef EDM_ML_DEBUG
867  edm::LogInfo("HcalSim") << "HCalSD::hitForFibre " << hits.size()
868  << " hits for " << GetName() << " of " << primaryID
869  << " with " << theTrack->GetDefinition()->GetParticleName()
870  << " of " << preStepPoint->GetKineticEnergy()/GeV
871  << " GeV in detector type " << det;
872 #endif
873  if (!hits.empty()) {
874  for (unsigned int i=0; i<hits.size(); ++i) {
875  G4ThreeVector hitPoint = hits[i].position;
876  if (isItinFidVolume (hitPoint)) {
877  int depth = hits[i].depth;
878  double time = hits[i].time;
879  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
880  currentID.setID(unitID, time, primaryID, 0);
881 #ifdef plotDebug
882  plotProfile(aStep, hitPoint, edepositEM, time, depth);
883  bool emType = false;
884  if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
885  emType = true;
886  plotHF(hitPoint,emType);
887 #endif
888  // check if it is in the same unit and timeslice as the previous one
889  if (currentID == previousID) {
891  } else {
892  posGlobal = preStepPoint->GetPosition();
893  if (!checkHit()) currentHit = createNewHit();
894  }
895  }
896  }
897  }
898 }
899 
900 void HCalSD::getFromParam (G4Step* aStep, double weight) {
901  std::vector<HFShowerParam::Hit> hits = showerParam->getHits(aStep, weight);
902  int nHit = static_cast<int>(hits.size());
903 
904  if (nHit > 0) {
905  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
906  int primaryID = setTrackID(aStep);
907 
908  int det = 5;
909 #ifdef EDM_ML_DEBUG
910  edm::LogInfo("HcalSim") << "HCalSD::getFromParam " << nHit << " hits for "
911  << GetName() << " of " << primaryID << " with "
912  << aStep->GetTrack()->GetDefinition()->GetParticleName()
913  << " of " << preStepPoint->GetKineticEnergy()/GeV
914  << " GeV in detector type " << det;
915 #endif
916  for (int i = 0; i<nHit; ++i) {
917  G4ThreeVector hitPoint = hits[i].position;
918  int depth = hits[i].depth;
919  double time = hits[i].time;
920  unsigned int unitID = setDetUnitId(det, hitPoint, depth);
921  currentID.setID(unitID, time, primaryID, 0);
922  edepositEM = hits[i].edep*GeV;
923  edepositHAD = 0.;
924 #ifdef plotDebug
925  plotProfile(aStep, hitPoint, edepositEM, time, depth);
926 #endif
927 
928  // check if it is in the same unit and timeslice as the previous one
929  if (currentID == previousID) {
931  } else {
932  posGlobal = preStepPoint->GetPosition();
933  if (!checkHit()) currentHit = createNewHit();
934  }
935  }
936  }
937 }
938 
939 void HCalSD::getHitPMT (const G4Step * aStep) {
940 
941  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
942  const G4Track* theTrack = aStep->GetTrack();
943  double edep = showerPMT->getHits(aStep);
944 
945  if (edep >= 0) {
946  double etrack = preStepPoint->GetKineticEnergy();
947  int primaryID = 0;
948  if (etrack >= energyCut) {
949  primaryID = theTrack->GetTrackID();
950  } else {
951  primaryID = theTrack->GetParentID();
952  if (primaryID == 0) primaryID = theTrack->GetTrackID();
953  }
954  // Reset entry point for new primary
955  posGlobal = preStepPoint->GetPosition();
956  resetForNewPrimary(posGlobal, etrack);
957 
958  //
959  int det = static_cast<int>(HcalForward);
960  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
961  double rr = (hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
962  double phi = (rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x()));
963  double etaR = showerPMT->getRadius();
964  int depth = 1;
965  if (etaR < 0) {
966  depth = 2;
967  etaR =-etaR;
968  }
969  if (hitPoint.z() < 0) etaR =-etaR;
970 #ifdef EDM_ML_DEBUG
971  edm::LogInfo("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR "
972  << etaR << " phi " << phi/deg << " depth " <<depth;
973 #endif
974  double time = (aStep->GetPostStepPoint()->GetGlobalTime());
975  uint32_t unitID = 0;
976  if (numberingFromDDD) {
978  depth,1);
979  unitID = setDetUnitId(tmp);
980  }
981  currentID.setID(unitID, time, primaryID, 1);
982 
983  edepositHAD = aStep->GetTotalEnergyDeposit();
984  edepositEM =-edepositHAD + (edep*GeV);
985 #ifdef plotDebug
986  plotProfile(aStep, hitPoint, edep*GeV, time, depth);
987 #endif
988 #ifdef EDM_ML_DEBUG
989  double beta = preStepPoint->GetBeta();
990  LogDebug("HcalSim") << "HCalSD::getHitPMT 1 hit for " << GetName()
991  << " of " << primaryID << " with "
992  << theTrack->GetDefinition()->GetParticleName()
993  << " of " << preStepPoint->GetKineticEnergy()/GeV
994  << " GeV with velocity " << beta << " UnitID "
995  << std::hex << unitID << std::dec;
996 #endif
997  // check if it is in the same unit and timeslice as the previous one
998  if (currentID == previousID) {
1000  } else {
1001  if (!checkHit()) currentHit = createNewHit();
1002  }
1003  }
1004 }
1005 
1006 void HCalSD::getHitFibreBundle (const G4Step* aStep, bool type) {
1007  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
1008  const G4Track* theTrack = aStep->GetTrack();
1009  double edep = showerBundle->getHits(aStep, type);
1010 
1011  if (edep >= 0) {
1012  double etrack = preStepPoint->GetKineticEnergy();
1013  int primaryID = 0;
1014  if (etrack >= energyCut) {
1015  primaryID = theTrack->GetTrackID();
1016  } else {
1017  primaryID = theTrack->GetParentID();
1018  if (primaryID == 0) primaryID = theTrack->GetTrackID();
1019  }
1020  // Reset entry point for new primary
1021  posGlobal = preStepPoint->GetPosition();
1022  resetForNewPrimary(posGlobal, etrack);
1023 
1024  //
1025  int det = static_cast<int>(HcalForward);
1026  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
1027  double rr = (hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
1028  double phi = (rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x()));
1029  double etaR = showerBundle->getRadius();
1030  int depth = 1;
1031  if (etaR < 0) {
1032  depth = 2;
1033  etaR =-etaR;
1034  }
1035  if (hitPoint.z() < 0) etaR =-etaR;
1036 #ifdef EDM_ML_DEBUG
1037  LogDebug("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR "
1038  << etaR << " phi " << phi/deg << " depth " <<depth;
1039 #endif
1040  double time = (aStep->GetPostStepPoint()->GetGlobalTime());
1041  uint32_t unitID = 0;
1042  if (numberingFromDDD) {
1043  HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det,etaR,phi,depth,1);
1044  unitID = setDetUnitId(tmp);
1045  }
1046  if (type) currentID.setID(unitID, time, primaryID, 3);
1047  else currentID.setID(unitID, time, primaryID, 2);
1048 
1049  edepositHAD = aStep->GetTotalEnergyDeposit();
1050  edepositEM =-edepositHAD + (edep*GeV);
1051 #ifdef plotDebug
1052  plotProfile(aStep, hitPoint, edep*GeV, time, depth);
1053 #endif
1054 #ifdef EDM_ML_DEBUG
1055  double beta = preStepPoint->GetBeta();
1056  LogDebug("HcalSim") << "HCalSD::getHitFibreBundle 1 hit for " << GetName()
1057  << " of " << primaryID << " with "
1058  << theTrack->GetDefinition()->GetParticleName()
1059  << " of " << preStepPoint->GetKineticEnergy()/GeV
1060  << " GeV with velocity " << beta << " UnitID "
1061  << std::hex << unitID << std::dec;
1062 #endif
1063  // check if it is in the same unit and timeslice as the previous one
1065  else if (!checkHit()) currentHit = createNewHit();
1066  } // non-zero energy deposit
1067 }
1068 
1069 int HCalSD::setTrackID (const G4Step* aStep) {
1070  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
1071  const G4Track* theTrack = aStep->GetTrack();
1072 
1073  double etrack = preStepPoint->GetKineticEnergy();
1074  TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
1075  int primaryID = trkInfo->getIDonCaloSurface();
1076  if (primaryID == 0) {
1077 #ifdef EDM_ML_DEBUG
1078  edm::LogInfo("HcalSim") << "HCalSD: Problem with primaryID **** set by "
1079  << "force to TkID **** " <<theTrack->GetTrackID();
1080 #endif
1081  primaryID = theTrack->GetTrackID();
1082  }
1083 
1084  if (primaryID != previousID.trackID())
1085  resetForNewPrimary(preStepPoint->GetPosition(), etrack);
1086 
1087  return primaryID;
1088 }
1089 
1091 
1092  std::ifstream infile;
1093  int entry=0;
1094  infile.open(fName.c_str(), std::ios::in);
1095  if (infile) {
1096  int det, zside, etaR, phi, lay;
1097  double wt;
1098  while (infile >> det >> zside >> etaR >> phi >> lay >> wt) {
1099  uint32_t id = HcalTestNumbering::packHcalIndex(det,zside,1,etaR,phi,lay);
1100  layerWeights.insert(std::pair<uint32_t,double>(id,wt));
1101  ++entry;
1102 #ifdef EDM_ML_DEBUG
1103  edm::LogInfo("HcalSim") << "HCalSD::readWeightFromFile:Entry " << entry
1104  << " ID " << std::hex << id << std::dec << " ("
1105  << det << "/" << zside << "/1/" << etaR << "/"
1106  << phi << "/" << lay << ") Weight " << wt;
1107 #endif
1108  }
1109  infile.close();
1110  }
1111  edm::LogInfo("HcalSim") << "HCalSD::readWeightFromFile: reads " << entry
1112  << " weights from " << fName;
1113  if (entry <= 0) useLayerWt = false;
1114 }
1115 
1116 double HCalSD::layerWeight(int det, const G4ThreeVector& pos, int depth, int lay) {
1117 
1118  double wt = 1.;
1119  if (numberingFromDDD) {
1120  //get the ID's as eta, phi, depth, ... indices
1122  depth, lay);
1123  modifyDepth(tmp);
1124  uint32_t id = HcalTestNumbering::packHcalIndex(tmp.subdet, tmp.zside, 1,
1125  tmp.etaR, tmp.phis,tmp.lay);
1126  std::map<uint32_t,double>::const_iterator ite = layerWeights.find(id);
1127  if (ite != layerWeights.end()) wt = ite->second;
1128 #ifdef EDM_ML_DEBUG
1129  edm::LogInfo("HcalSim") << "HCalSD::layerWeight: ID " << std::hex << id
1130  << std::dec << " (" << tmp.subdet << "/"
1131  << tmp.zside << "/1/" << tmp.etaR << "/"
1132  << tmp.phis << "/" << tmp.lay << ") Weight " <<wt;
1133 #endif
1134  }
1135  return wt;
1136 }
1137 
1138 void HCalSD::plotProfile(const G4Step* aStep,const G4ThreeVector& global, double edep,
1139  double time, int id) {
1140 
1141  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
1142  static const G4String modName[8] = {"HEModule", "HVQF" , "HBModule", "MBAT",
1143  "MBBT" , "MBBTC", "MBBT_R1P", "MBBT_R1M"};
1144  G4ThreeVector local;
1145  bool found=false;
1146  double depth=-2000;
1147  int idx = 4;
1148  for (int n=0; n<touch->GetHistoryDepth(); ++n) {
1149  G4String name = touch->GetVolume(n)->GetName();
1150 #ifdef EDM_ML_DEBUG
1151  LogDebug("HcalSim") << "plotProfile Depth " << n << " Name " << name;
1152 #endif
1153  for (unsigned int ii=0; ii<8; ++ii) {
1154  if (name == modName[ii]) {
1155  found = true;
1156  int dn = touch->GetHistoryDepth() - n;
1157  local = touch->GetHistory()->GetTransform(dn).TransformPoint(global);
1158  if (ii == 0) {depth = local.z() - 4006.5; idx = 1;}
1159  else if (ii == 1) {depth = local.z() + 825.0; idx = 3;}
1160  else if (ii == 2) {depth = local.x() - 1775.; idx = 0;}
1161  else {depth = local.y() + 15.; idx = 2;}
1162  break;
1163  }
1164  }
1165  if (found) break;
1166  }
1167  if (!found) depth = std::abs(global.z()) - 11500;
1168 #ifdef EDM_ML_DEBUG
1169  LogDebug("HcalSim") << "plotProfile Found " << found << " Global " << global
1170  << " Local " << local << " depth " << depth << " ID "
1171  << id << " EDEP " << edep << " Time " << time;
1172 #endif
1173  if (hit_[idx] != nullptr) hit_[idx]->Fill(edep);
1174  if (time_[idx] != nullptr) time_[idx]->Fill(time,edep);
1175  if (dist_[idx] != nullptr) dist_[idx]->Fill(depth,edep);
1176  int jd = 2*idx + id - 7;
1177  if (jd >= 0 && jd < 4) {
1178  jd += 5;
1179  if (hit_[jd] != nullptr) hit_[jd]->Fill(edep);
1180  if (time_[jd] != nullptr) time_[jd]->Fill(time,edep);
1181  if (dist_[jd] != nullptr) dist_[jd]->Fill(depth,edep);
1182  }
1183 }
1184 
1185 void HCalSD::plotHF(const G4ThreeVector& hitPoint, bool emType) {
1186  double zv = std::abs(hitPoint.z()) - gpar[4];
1187  if (emType) {
1188  if (hzvem != nullptr) hzvem->Fill(zv);
1189  } else {
1190  if (hzvhad != nullptr) hzvhad->Fill(zv);
1191  }
1192 }
1193 
1195  if (id.subdet == 4) {
1196  int ieta = (id.zside == 0) ? -id.etaR : id.etaR;
1197  if (hcalConstants->maxHFDepth(ieta,id.phis) > 2) {
1198  if (id.depth <= 2) {
1199  if (G4UniformRand() > 0.5) id.depth += 2;
1200  }
1201  }
1202  } else if ((id.subdet == 1 || id.subdet ==2) && testNumber) {
1203  if (depth_ == 0) id.depth = 1;
1204  else id.depth = 2;
1205  }
1206 }
#define LogDebug(id)
float edepositEM
Definition: CaloSD.h:122
const double beta
double energyCut
Definition: CaloSD.h:124
type
Definition: HCALResponse.h:21
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
double getAttenuation(const G4Step *aStep, double birk1, double birk2, double birk3)
Definition: CaloSD.cc:439
void readWeightFromFile(const std::string &)
Definition: HCalSD.cc:1090
std::vector< const G4LogicalVolume * > hfLV
Definition: HCalSD.h:107
bool useParam
Definition: HCalSD.h:99
double eminHitHE
Definition: HCalSD.h:100
HFShowerParam * showerParam
Definition: HCalSD.h:88
TH1F * time_[9]
Definition: HCalSD.h:109
G4int emPDG
Definition: CaloSD.h:137
double getRadius()
Definition: HFShowerPMT.cc:134
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:140
const double GeV
Definition: MathUtil.h:16
double kmaxNeutron
Definition: CaloSD.h:135
void updateHit(CaloG4Hit *)
Definition: CaloSD.cc:414
const N & name() const
Definition: DDBase.h:78
void getFromParam(G4Step *step, double weight)
Definition: HCalSD.cc:900
G4int depth_
Definition: HCalSD.h:102
int getIDonCaloSurface() const
bool useLayerWt
Definition: HCalSD.h:96
Definition: CaloSD.h:42
std::vector< double > gpar
Definition: HCalSD.h:103
bool useFibreBundle
Definition: HCalSD.h:96
int zside() const
get the z-side of the cell (1/-1)
Definition: HcalDetId.cc:93
double betaThr
Definition: HCalSD.h:98
void initRun(G4ParticleTable *, HcalDDDSimConstants *)
Definition: HFShower.cc:454
double deliveredLumi
Definition: HCalSD.h:101
std::vector< Hit > getHits(G4Step *aStep, bool &ok, double weight, bool onlyLong=false)
virtual uint32_t getUnitID(const HcalNumberingFromDDD::HcalID &id)
HFShowerFibreBundle * showerBundle
Definition: HCalSD.h:90
double eminHitHB
Definition: HCalSD.h:100
bool useShowerLibrary
Definition: HCalSD.h:99
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &)
Definition: HCalSD.cc:669
void plotProfile(const G4Step *step, const G4ThreeVector &pos, double edep, double time, int id)
Definition: HCalSD.cc:1138
Definition: weight.py:1
double birk2
Definition: HCalSD.h:98
void getFromLibrary(G4Step *step, double weight)
Definition: HCalSD.cc:785
bool usePMTHit
Definition: HCalSD.h:96
std::vector< const G4LogicalVolume * > fibre2LV
Definition: HCalSD.h:107
G4ThreeVector posGlobal
Definition: CaloSD.h:114
int zside(DetId const &)
static const unsigned int numberOfZLayers
Definition: HFDarkening.h:26
void modifyDepth(HcalNumberingFromDDD::HcalID &id)
Definition: HCalSD.cc:1194
void update(const BeginOfJob *) override
This routine will be called when the appropriate signal arrives.
Definition: HCalSD.cc:572
#define nullptr
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:564
type of data representation of DDCompactView
Definition: DDCompactView.h:90
bool agingFlagHE
Definition: HCalSD.h:91
double birk1
Definition: HCalSD.h:98
double kmaxProton
Definition: CaloSD.h:135
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:81
bool compare(HcalNumberingFromDDD::HcalID const &, uint32_t const &)
Definition: HcalTestNS.cc:25
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
static int PMTNumber(const G4ThreeVector &pe_effect)
G4bool checkHit()
Definition: CaloSD.cc:305
static TrackerG4SimHitNumberingScheme & numberingScheme(const DDCompactView &cpv, const GeometricDet &det)
TH1F * hzvhad
Definition: HCalSD.h:109
std::vector< Hit > getHits(const G4Step *aStep, double weight)
Definition: HFShower.cc:47
int setTrackID(const G4Step *step)
Definition: HCalSD.cc:1069
bool isItConicalBundle(const G4LogicalVolume *)
Definition: HCalSD.cc:756
void initRun(G4ParticleTable *, HcalDDDSimConstants *)
Definition: HFShowerPMT.cc:77
bool ProcessHits(G4Step *, G4TouchableHistory *) override
Definition: HCalSD.cc:358
double eminHitHF
Definition: HCalSD.h:100
const double MeV
static uint32_t packHcalIndex(int det, int z, int depth, int eta, int phi, int lay)
TH1F * hit_[9]
Definition: HCalSD.h:109
double getHits(const G4Step *aStep)
Definition: HFShowerPMT.cc:88
void resetForNewPrimary(const G4ThreeVector &, double)
Definition: CaloSD.cc:428
int getMaxDepth(const int &type) const
HCalSD(const std::string &, const DDCompactView &, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)
Definition: HCalSD.cc:44
bool agingFlagHB
Definition: HCalSD.h:91
bool testNumber
Definition: HCalSD.h:97
double kmaxIon
Definition: CaloSD.h:135
bool suppressHeavy
Definition: CaloSD.h:133
bool next()
set current node to the next node in the filtered tree
const std::vector< std::string > & getNames() const
HFShower * hfshower
Definition: HCalSD.h:87
bool useBirk
Definition: HCalSD.h:96
float edepositHAD
Definition: CaloSD.h:122
std::unique_ptr< HFDarkening > m_HFDarkening
Definition: HCalSD.h:94
double birk3
Definition: HCalSD.h:98
int trackID() const
Definition: CaloHitID.h:25
G4int epPDG
Definition: CaloSD.h:137
const HBHEDarkening * m_HEDarkening
Definition: HCalSD.h:93
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
std::maps an index to a DDValue. The index corresponds to the index assigned to the name of the std::...
Definition: DDsvalues.h:20
G4int gammaPDG
Definition: CaloSD.h:137
bool isAvailable() const
Definition: Service.h:46
bool applyFidCut
Definition: HCalSD.h:99
TH1F * dist_[9]
Definition: HCalSD.h:109
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void plotHF(const G4ThreeVector &pos, bool emType)
Definition: HCalSD.cc:1185
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
A DDLogicalPart aggregates information concerning material, solid and sensitveness ...
Definition: DDLogicalPart.h:92
HcalNumberingScheme * numberingScheme
Definition: HCalSD.h:85
T * make(const Args &...args) const
make new ROOT object
double getHits(const G4Step *aStep, bool type)
int getNumberOfHits()
Definition: CaloSD.cc:335
CaloHitID previousID
Definition: CaloSD.h:118
CaloG4Hit * currentHit
Definition: CaloSD.h:129
std::vector< int > hfLevels
Definition: HCalSD.h:104
void initRun(G4ParticleTable *, HcalDDDSimConstants *)
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:38
uint32_t setDetUnitId(const G4Step *step) override
Definition: HCalSD.cc:551
G4Track * theTrack
Definition: CaloSD.h:119
double getLayer0Wt(const int &det, const int &phi, const int &zside) const
virtual G4bool getStepInfo(G4Step *aStep)
Definition: CaloSD.cc:224
bool isItStraightBundle(const G4LogicalVolume *)
Definition: HCalSD.cc:750
bool testNS_
Definition: HCalSD.h:97
G4int mupPDG
Definition: HCalSD.h:102
ii
Definition: cuy.py:588
~HCalSD() override
Definition: HCalSD.cc:346
double eminHitHO
Definition: HCalSD.h:100
double tmaxHit
Definition: CaloSD.h:124
int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.cc:98
void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth=0)
Definition: CaloHitID.cc:44
TH1F * hzvem
Definition: HCalSD.h:109
G4StepPoint * preStepPoint
Definition: CaloSD.h:121
int iphi() const
get the cell iphi
Definition: HcalDetId.cc:103
const std::vector< double > & getGparHF() const
std::vector< const G4LogicalVolume * > pmtLV
Definition: HCalSD.h:107
Definition: DetId.h:18
static const unsigned int lowZLimit
Definition: HFDarkening.h:29
static const G4LogicalVolume * GetVolume(const std::string &name)
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
bool filterHit(CaloG4Hit *, double) override
Definition: HCalSD.cc:625
CaloHitID currentID
Definition: CaloSD.h:118
std::vector< Hit > getHits(G4Step *aStep, double weight)
static const unsigned int upperZLimit
Definition: HFDarkening.h:30
bool neutralDensity
Definition: HCalSD.h:97
void initRun(G4ParticleTable *, HcalDDDSimConstants *)
void getHitFibreBundle(const G4Step *step, bool type)
Definition: HCalSD.cc:1006
const T & get() const
Definition: EventSetup.h:59
void hitForFibre(const G4Step *step, double weight)
Definition: HCalSD.cc:848
bool isPrimary() const
bool isItPMT(const G4LogicalVolume *)
Definition: HCalSD.cc:744
DDsvalues_type mergedSpecifics() const
bool isItinFidVolume(const G4ThreeVector &)
Definition: HCalSD.cc:768
void getHitPMT(const G4Step *step)
Definition: HCalSD.cc:939
HcalNumberingFromDDD * numberingFromDDD
Definition: HCalSD.h:84
double layerWeight(int, const G4ThreeVector &, int, int)
Definition: HCalSD.cc:1116
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
HLT enums.
bool isItHF(const G4Step *)
Definition: HCalSD.cc:714
std::vector< G4String > fibreNames
Definition: HCalSD.h:105
G4bool hitExists()
Definition: CaloSD.cc:284
bool firstChild()
set the current node to the first child ...
bool isItScintillator(const G4Material *)
Definition: HCalSD.cc:762
#define EDM_ML_DEBUG
std::vector< G4String > hfNames
Definition: HCalSD.h:105
std::vector< const G4LogicalVolume * > fibreLV
Definition: HCalSD.h:107
G4int mumPDG
Definition: HCalSD.h:102
HcalTestNS * hcalTestNS_
Definition: HCalSD.h:95
std::map< uint32_t, double > layerWeights
Definition: HCalSD.h:108
HFShowerPMT * showerPMT
Definition: HCalSD.h:89
uint32_t getUnitID() const
Definition: CaloG4Hit.h:69
double getResponseWt(const G4Track *)
Definition: CaloSD.cc:581
std::vector< const G4LogicalVolume * > fibre1LV
Definition: HCalSD.h:107
HcalID unitID(int det, const CLHEP::Hep3Vector &pos, int depth, int lay=-1) const
HFShowerLibrary * showerLibrary
Definition: HCalSD.h:86
int maxHFDepth(const int &ieta, const int &iphi) const
std::vector< const G4Material * > materials
Definition: HCalSD.h:106
void initRun(G4ParticleTable *, HcalDDDSimConstants *)
void initRun() override
Definition: HCalSD.cc:609
const HBHEDarkening * m_HBDarkening
Definition: HCalSD.h:92
double getEnergyDeposit(G4Step *) override
Definition: HCalSD.cc:463
const std::string & name() const
Returns the name.
Definition: DDName.cc:90
const DDMaterial & material(void) const
Returns a reference object of the material this LogicalPart is made of.
CaloG4Hit * createNewHit()
Definition: CaloSD.cc:337
std::vector< G4String > matNames
Definition: HCalSD.h:105
float degradation(float intlumi, int ieta, int lay) const
double getEnergyDeposit() const
Definition: CaloG4Hit.h:81
bool useHF
Definition: HCalSD.h:99
bool isItFibre(const G4LogicalVolume *)
Definition: HCalSD.cc:732
void NaNTrap(const G4Step *step) const
HcalDDDSimConstants * hcalConstants
Definition: HCalSD.h:83