CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
MaterialBudgetHcalHistos.cc
Go to the documentation of this file.
1 #include "Validation/Geometry/interface/MaterialBudgetHcalHistos.h"
2 
3 #include "DataFormats/Math/interface/GeantUnits.h"
4 #include "DetectorDescription/Core/interface/DDFilter.h"
5 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
6 #include "DetectorDescription/Core/interface/DDSplit.h"
7 #include "DetectorDescription/Core/interface/DDValue.h"
8 #include "FWCore/Utilities/interface/Exception.h"
9 #include "FWCore/MessageLogger/interface/MessageLogger.h"
10 #include "FWCore/ServiceRegistry/interface/Service.h"
11 #include "CommonTools/UtilAlgos/interface/TFileService.h"
12 
13 #include "DD4hep/Filter.h"
14 
15 #include <string>
16 #include <vector>
17 
18 using namespace geant_units::operators;
19 
20 MaterialBudgetHcalHistos::MaterialBudgetHcalHistos(const edm::ParameterSet& p) {
21  binEta_ = p.getUntrackedParameter<int>("NBinEta", 260);
22  binPhi_ = p.getUntrackedParameter<int>("NBinPhi", 180);
23  maxEta_ = p.getUntrackedParameter<double>("MaxEta", 5.2);
24  etaLow_ = p.getUntrackedParameter<double>("EtaLow", -5.2);
25  etaHigh_ = p.getUntrackedParameter<double>("EtaHigh", 5.2);
26  fillHistos_ = p.getUntrackedParameter<bool>("FillHisto", true);
27  printSum_ = p.getUntrackedParameter<bool>("PrintSummary", false);
28  fromdd4hep_ = p.getUntrackedParameter<bool>("Fromdd4hep", false);
29  etaMinP_ = p.getUntrackedParameter<double>("EtaMinP", 5.2);
30  etaMaxP_ = p.getUntrackedParameter<double>("EtaMaxP", 0.0);
31  etaLowMin_ = p.getUntrackedParameter<double>("EtaLowMin", 0.783);
32  etaLowMax_ = p.getUntrackedParameter<double>("EtaLowMax", 0.870);
33  etaMidMin_ = p.getUntrackedParameter<double>("EtaMidMin", 2.650);
34  etaMidMax_ = p.getUntrackedParameter<double>("EtaMidMax", 2.868);
35  etaHighMin_ = p.getUntrackedParameter<double>("EtaHighMin", 2.868);
36  etaHighMax_ = p.getUntrackedParameter<double>("EtaHighMax", 3.000);
37  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: FillHisto : " << fillHistos_ << " PrintSummary "
38  << printSum_ << " == Eta plot: NX " << binEta_ << " Range " << -maxEta_ << ":"
39  << maxEta_ << " Phi plot: NX " << binPhi_ << " Range " << -1._pi << ":" << 1._pi
40  << " (Eta limit " << etaLow_ << ":" << etaHigh_ << ")"
41  << " Eta range (" << etaLowMin_ << ":" << etaLowMax_ << "), (" << etaMidMin_ << ":"
42  << etaMidMax_ << "), (" << etaHighMin_ << ":" << etaHighMax_
43  << ") Debug for eta range " << etaMinP_ << ":" << etaMaxP_ << " FromDD4hep "
44  << fromdd4hep_;
45  if (fillHistos_)
46  book();
47 }
48 
49 void MaterialBudgetHcalHistos::fillBeginJob(const DDCompactView& cpv) {
50  constexpr int32_t addLevel = 0;
51  if (fillHistos_) {
52  std::string attribute = "ReadOutName";
53  std::string value = "HcalHits";
54  DDSpecificsMatchesValueFilter filter1{DDValue(attribute, value, 0)};
55  DDFilteredView fv1(cpv, filter1);
56  std::vector<std::string> names = getNames(fv1);
57  for (auto& name : names) {
58  std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
59  if (std::find(sensitives_.begin(), sensitives_.end(), namx) == sensitives_.end())
60  sensitives_.emplace_back(namx);
61  }
62  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
63  << value << " has " << sensitives_.size() << " elements";
64  for (unsigned int i = 0; i < sensitives_.size(); i++)
65  edm::LogVerbatim("MaterialBudgetFull")
66  << "MaterialBudgetHcalHistos: sensitives[" << i << "] = " << sensitives_[i];
67  attribute = "Volume";
68  value = "HF";
69  DDSpecificsMatchesValueFilter filter2{DDValue(attribute, value, 0)};
70  DDFilteredView fv2(cpv, filter2);
71  hfNames_ = getNames(fv2);
72  fv2.firstChild();
73  DDsvalues_type sv(fv2.mergedSpecifics());
74  std::vector<double> temp = getDDDArray("Levels", sv);
75  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
76  << value << " has " << hfNames_.size() << " elements";
77  for (unsigned int i = 0; i < hfNames_.size(); i++) {
78  int level = static_cast<int>(temp[i]);
79  hfLevels_.push_back(level + addLevel);
80  edm::LogVerbatim("MaterialBudgetFull")
81  << "MaterialBudgetHcalHistos: HF[" << i << "] = " << hfNames_[i] << " at level " << hfLevels_[i];
82  }
83 
84  const std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
85  attribute = "ReadOutName";
86  for (int k = 0; k < 2; k++) {
87  value = ecalRO[k];
88  DDSpecificsMatchesValueFilter filter3{DDValue(attribute, value, 0)};
89  DDFilteredView fv3(cpv, filter3);
90  std::vector<std::string> senstmp = getNames(fv3);
91  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute
92  << " = " << value << " has " << senstmp.size() << " elements";
93  for (unsigned int i = 0; i < senstmp.size(); i++) {
94  std::string name = senstmp[i].substr(0, 4);
95  if (std::find(sensitiveEC_.begin(), sensitiveEC_.end(), name) == sensitiveEC_.end())
96  sensitiveEC_.push_back(name);
97  }
98  }
99  for (unsigned int i = 0; i < sensitiveEC_.size(); i++)
100  edm::LogVerbatim("MaterialBudgetFull")
101  << "MaterialBudgetHcalHistos:sensitiveEC[" << i << "] = " << sensitiveEC_[i];
102  }
103 }
104 
105 void MaterialBudgetHcalHistos::fillBeginJob(const cms::DDCompactView& cpv) {
106  constexpr int32_t addLevel = 1;
107  if (fillHistos_) {
108  std::string attribute = "ReadOutName";
109  std::string value = "HcalHits";
110  const cms::DDFilter filter1(attribute, value);
111  cms::DDFilteredView fv1(cpv, filter1);
112  std::vector<std::string> names = getNames(fv1);
113  for (auto& name : names) {
114  std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
115  if (std::find(sensitives_.begin(), sensitives_.end(), namx) == sensitives_.end())
116  sensitives_.emplace_back(namx);
117  }
118  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
119  << value << " has " << sensitives_.size() << " elements";
120  for (unsigned int i = 0; i < sensitives_.size(); i++)
121  edm::LogVerbatim("MaterialBudgetFull")
122  << "MaterialBudgetHcalHistos: sensitives[" << i << "] = " << sensitives_[i];
123  attribute = "Volume";
124  value = "HF";
125  const cms::DDFilter filter2(attribute, value);
126  cms::DDFilteredView fv2(cpv, filter2);
127  std::vector<int> temp = fv2.get<std::vector<int> >("hf", "Levels");
128  hfNames_ = getNames(fv2);
129  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
130  << value << " has " << hfNames_.size() << " elements";
131  for (unsigned int i = 0; i < hfNames_.size(); i++) {
132  hfLevels_.push_back(temp[i] + addLevel);
133  edm::LogVerbatim("MaterialBudgetFull")
134  << "MaterialBudgetHcalHistos: HF[" << i << "] = " << hfNames_[i] << " at level " << hfLevels_[i];
135  }
136 
137  const std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
138  attribute = "ReadOutName";
139  for (int k = 0; k < 2; k++) {
140  value = ecalRO[k];
141  const cms::DDFilter filter(attribute, value);
142  cms::DDFilteredView fv(cpv, filter);
143  std::vector<std::string> senstmp = getNames(fv);
144  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute
145  << " = " << value << " has " << senstmp.size() << " elements";
146  for (unsigned int i = 0; i < senstmp.size(); i++) {
147  std::string name = senstmp[i].substr(0, 4);
148  if (std::find(sensitiveEC_.begin(), sensitiveEC_.end(), name) == sensitiveEC_.end())
149  sensitiveEC_.push_back(name);
150  }
151  }
152  for (unsigned int i = 0; i < sensitiveEC_.size(); i++)
153  edm::LogVerbatim("MaterialBudgetFull")
154  << "MaterialBudgetHcalHistos:sensitiveEC[" << i << "] = " << sensitiveEC_[i];
155  }
156 }
157 
158 void MaterialBudgetHcalHistos::fillStartTrack(const G4Track* aTrack) {
159  id_ = layer_ = steps_ = 0;
160  radLen_ = intLen_ = stepLen_ = 0;
161  nlayHB_ = nlayHE_ = nlayHF_ = nlayHO_ = 0;
162 
163  const G4ThreeVector& dir = aTrack->GetMomentum();
164  if (dir.theta() != 0) {
165  eta_ = dir.eta();
166  } else {
167  eta_ = -99;
168  }
169  phi_ = dir.phi();
170  double theEnergy = aTrack->GetTotalEnergy();
171  int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
172 
173  if (printSum_) {
174  matList_.clear();
175  stepLength_.clear();
176  radLength_.clear();
177  intLength_.clear();
178  }
179 
180  if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
181  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID() << " Code "
182  << theID << " Energy " << convertUnitsTo(1._GeV, theEnergy) << " GeV; Eta "
183  << eta_ << " Phi " << convertRadToDeg(phi_) << " PT "
184  << convertUnitsTo(1._GeV, dir.perp()) << " GeV *****";
185 }
186 
187 void MaterialBudgetHcalHistos::fillPerStep(const G4Step* aStep) {
188  G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
189  double step = aStep->GetStepLength();
190  double radl = material->GetRadlen();
191  double intl = material->GetNuclearInterLength();
192  double density = convertUnitsTo(1._g_per_cm3, material->GetDensity());
193 
194  int idOld = id_;
195  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
196  std::string name = (static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(0)->GetName())));
197  std::string matName = (static_cast<std::string>(dd4hep::dd::noNamespace(material->GetName())));
198  if (printSum_) {
199  bool found = false;
200  for (unsigned int ii = 0; ii < matList_.size(); ii++) {
201  if (matList_[ii] == matName) {
202  stepLength_[ii] += step;
203  radLength_[ii] += (step / radl);
204  intLength_[ii] += (step / intl);
205  found = true;
206  break;
207  }
208  }
209  if (!found) {
210  matList_.push_back(matName);
211  stepLength_.push_back(step);
212  radLength_.push_back(step / radl);
213  intLength_.push_back(step / intl);
214  }
215  if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
216  edm::LogVerbatim("MaterialBudget") << "Volume " << name << " id " << id_ << ":" << idOld << " Step " << step
217  << " Material " << matName << " Old Length " << stepLen_ << " X0 "
218  << step / radl << ":" << radLen_ << " Lambda " << step / intl << ":"
219  << intLen_;
220  } else {
221  if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
222  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Step at " << name << " id " << id_ << ":"
223  << idOld << " Length " << step << " in " << matName << " of density "
224  << density << " g/cc; Radiation Length " << radl << " mm; Interaction Length "
225  << intl << " mm\n Position "
226  << aStep->GetPreStepPoint()->GetPosition() << " Cylindrical R "
227  << aStep->GetPreStepPoint()->GetPosition().perp() << " Length (so far) "
228  << stepLen_ << " L/X0 " << step / radl << "/" << radLen_ << " L/Lambda "
229  << step / intl << "/" << intLen_;
230  }
231 
232  int det = 0, lay = 0;
233  double abseta = std::abs(eta_);
234  if (fillHistos_) {
235  edm::LogVerbatim("MaterialBudgetFull")
236  << "Volume " << name << ":" << matName << " EC:Sensitive:HF " << isItEC(name) << ":" << isSensitive(name) << ":"
237  << isItHF(touch) << " Eta " << abseta << " HC " << ((touch->GetReplicaNumber(1)) / 1000) << ":"
238  << ((touch->GetReplicaNumber(0) / 10) % 100 + 3) << " X0 " << (radLen_ + (step / radl)) << " Lambda "
239  << (intLen_ + (step / intl));
240  ;
241  if (isItEC(name)) {
242  det = 1;
243  lay = 1;
244  } else {
245  if (isSensitive(name)) {
246  if (isItHF(touch)) {
247  det = 5;
248  lay = 21;
249  if (lay != layer_)
250  ++nlayHF_;
251  } else {
252  det = (touch->GetReplicaNumber(1)) / 1000;
253  lay = (touch->GetReplicaNumber(0) / 10) % 100 + 3;
254  if (det == 4) {
255  if (abseta < 1.479)
256  lay = layer_ + 1;
257  else
258  lay--;
259  if (lay < 3)
260  lay = 3;
261  if (lay == layer_)
262  lay++;
263  if (lay > 20)
264  lay = 20;
265  if (lay != layer_)
266  ++nlayHE_;
267  } else if (lay != layer_) {
268  if (lay >= 20)
269  ++nlayHO_;
270  else
271  ++nlayHB_;
272  }
273  }
274  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Det " << det << " Layer " << lay << " Eta "
275  << eta_ << " Phi " << convertRadToDeg(phi_);
276  } else if (layer_ == 1) {
277  det = -1;
278  lay = 2;
279  }
280  }
281  if (det != 0) {
282  if (lay != layer_) {
283  id_ = lay;
284  layer_ = lay;
285  }
286  }
287 
288  if (id_ > idOld) {
289  if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
290  edm::LogVerbatim("MaterialBudget")
291  << "MaterialBudgetHcalHistos: Step at " << name << " calls filHisto with " << (id_ - 1);
292  fillHisto(id_ - 1);
293  }
294  }
295 
296  stepLen_ += step;
297  radLen_ += (step / radl);
298  intLen_ += (step / intl);
299  if (fillHistos_) {
300  if (id_ == 21) {
301  if (!isItHF(aStep->GetPostStepPoint()->GetTouchable())) {
302  if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
303  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: After HF in " << name << ":"
304  << static_cast<std::string>(dd4hep::dd::noNamespace(
305  aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)->GetName()))
306  << " calls fillHisto with " << id_;
307  fillHisto(idOld);
308  ++id_;
309  layer_ = 0;
310  }
311  }
312  }
313 }
314 
315 void MaterialBudgetHcalHistos::fillEndTrack() {
316  if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
317  edm::LogVerbatim("MaterialBudget") << "Number of layers hit in HB:" << nlayHB_ << " HE:" << nlayHE_
318  << " HO:" << nlayHO_ << " HF:" << nlayHF_;
319  if (fillHistos_) {
320  fillHisto(maxSet_ - 1);
321  fillLayer();
322  }
323  if (printSum_) {
324  for (unsigned int ii = 0; ii < matList_.size(); ii++) {
325  edm::LogVerbatim("MaterialBudget") << matList_[ii] << "\t" << stepLength_[ii] << "\t" << radLength_[ii] << "\t"
326  << intLength_[ii];
327  }
328  }
329 }
330 
331 void MaterialBudgetHcalHistos::book() {
332  // Book histograms
333  edm::Service<TFileService> tfile;
334 
335  if (!tfile.isAvailable())
336  throw cms::Exception("BadConfig") << "TFileService unavailable: "
337  << "please add it to config file";
338 
339  double maxPhi = 1._pi;
340  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Booking user histos === with " << binEta_
341  << " bins in eta from " << -maxEta_ << " to " << maxEta_ << " and " << binPhi_
342  << " bins in phi from " << -maxPhi << " to " << maxPhi;
343 
344  std::string iter;
345  std::string range0 = "(" + std::to_string(etaMidMin_) + ":" + std::to_string(etaMidMax_) + ") ";
346  std::string range1 = "(" + std::to_string(etaHighMin_) + ":" + std::to_string(etaHighMax_) + ") ";
347  std::string range2 = "(" + std::to_string(etaLowMin_) + ":" + std::to_string(etaLowMax_) + ") ";
348  // total X0
349  for (int i = 0; i < maxSet_; i++) {
350  iter = std::to_string(i);
351  me100[i] = tfile->make<TProfile>(
352  std::to_string(i + 100).c_str(), ("MB(X0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
353  me200[i] = tfile->make<TProfile>(
354  std::to_string(i + 200).c_str(), ("MB(L0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
355  me300[i] = tfile->make<TProfile>(
356  std::to_string(i + 300).c_str(), ("MB(Step) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
357  me400[i] = tfile->make<TH1F>(
358  std::to_string(i + 400).c_str(), ("Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
359  me500[i] = tfile->make<TProfile>(
360  std::to_string(i + 500).c_str(), ("MB(X0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
361  me600[i] = tfile->make<TProfile>(
362  std::to_string(i + 600).c_str(), ("MB(L0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
363  me700[i] = tfile->make<TProfile>(
364  std::to_string(i + 700).c_str(), ("MB(Step) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
365  me800[i] =
366  tfile->make<TH1F>(std::to_string(i + 800).c_str(), ("Phi in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
367  me900[i] = tfile->make<TProfile2D>(std::to_string(i + 900).c_str(),
368  ("MB(X0) prof Eta Phi in region " + iter).c_str(),
369  binEta_ / 2,
370  -maxEta_,
371  maxEta_,
372  binPhi_ / 2,
373  -maxPhi,
374  maxPhi);
375  me1000[i] = tfile->make<TProfile2D>(std::to_string(i + 1000).c_str(),
376  ("MB(L0) prof Eta Phi in region " + iter).c_str(),
377  binEta_ / 2,
378  -maxEta_,
379  maxEta_,
380  binPhi_ / 2,
381  -maxPhi,
382  maxPhi);
383  me1100[i] = tfile->make<TProfile2D>(std::to_string(i + 1100).c_str(),
384  ("MB(Step) prof Eta Phi in region " + iter).c_str(),
385  binEta_ / 2,
386  -maxEta_,
387  maxEta_,
388  binPhi_ / 2,
389  -maxPhi,
390  maxPhi);
391  me1200[i] = tfile->make<TH2F>(std::to_string(i + 1200).c_str(),
392  ("Eta vs Phi in region " + iter).c_str(),
393  binEta_ / 2,
394  -maxEta_,
395  maxEta_,
396  binPhi_ / 2,
397  -maxPhi,
398  maxPhi);
399  me1600[i] = tfile->make<TProfile>(std::to_string(i + 1600).c_str(),
400  ("MB(X0) prof Ph in region " + range0 + iter).c_str(),
401  binPhi_,
402  -maxPhi,
403  maxPhi);
404  me1700[i] = tfile->make<TProfile>(std::to_string(i + 1700).c_str(),
405  ("MB(L0) prof Ph in region " + range0 + iter).c_str(),
406  binPhi_,
407  -maxPhi,
408  maxPhi);
409  me1800[i] = tfile->make<TProfile>(std::to_string(i + 1800).c_str(),
410  ("MB(Step) prof Ph in region " + range0 + iter).c_str(),
411  binPhi_,
412  -maxPhi,
413  maxPhi);
414  me1900[i] = tfile->make<TProfile>(std::to_string(i + 1900).c_str(),
415  ("MB(X0) prof Ph in region " + range1 + iter).c_str(),
416  binPhi_,
417  -maxPhi,
418  maxPhi);
419  me2000[i] = tfile->make<TProfile>(std::to_string(i + 2000).c_str(),
420  ("MB(L0) prof Ph in region " + range1 + iter).c_str(),
421  binPhi_,
422  -maxPhi,
423  maxPhi);
424  me2100[i] = tfile->make<TProfile>(std::to_string(i + 2100).c_str(),
425  ("MB(Step) prof Ph in region " + range1 + iter).c_str(),
426  binPhi_,
427  -maxPhi,
428  maxPhi);
429  me2200[i] = tfile->make<TProfile>(std::to_string(i + 2200).c_str(),
430  ("MB(X0) prof Ph in region " + range2 + iter).c_str(),
431  binPhi_,
432  -maxPhi,
433  maxPhi);
434  me2300[i] = tfile->make<TProfile>(std::to_string(i + 2300).c_str(),
435  ("MB(L0) prof Ph in region " + range2 + iter).c_str(),
436  binPhi_,
437  -maxPhi,
438  maxPhi);
439  me2400[i] = tfile->make<TProfile>(std::to_string(i + 2400).c_str(),
440  ("MB(Step) prof Ph in region " + range2 + iter).c_str(),
441  binPhi_,
442  -maxPhi,
443  maxPhi);
444  }
445  for (int i = 0; i < maxSet2_; i++) {
446  iter = std::to_string(i);
447  me1300[i] = tfile->make<TH1F>(std::to_string(i + 1300).c_str(),
448  ("Events with layers Hit (0 all, 1 HB, ..) for " + iter).c_str(),
449  binEta_,
450  -maxEta_,
451  maxEta_);
452  me1400[i] = tfile->make<TH2F>(std::to_string(i + 1400).c_str(),
453  ("Eta vs Phi for layers hit in " + iter).c_str(),
454  binEta_ / 2,
455  -maxEta_,
456  maxEta_,
457  binPhi_ / 2,
458  -maxPhi,
459  maxPhi);
460  me1500[i] = tfile->make<TProfile>(std::to_string(i + 1500).c_str(),
461  ("Number of layers crossed (0 all, 1 HB, ..) for " + iter).c_str(),
462  binEta_,
463  -maxEta_,
464  maxEta_);
465  }
466 
467  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user histos done ===";
468 }
469 
470 void MaterialBudgetHcalHistos::fillHisto(int ii) {
471  if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
472  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos:FillHisto called with index " << ii
473  << " integrated step " << stepLen_ << " X0 " << radLen_ << " Lamda " << intLen_;
474 
475  if (ii >= 0 && ii < maxSet_) {
476  me100[ii]->Fill(eta_, radLen_);
477  me200[ii]->Fill(eta_, intLen_);
478  me300[ii]->Fill(eta_, stepLen_);
479  me400[ii]->Fill(eta_);
480 
481  if (eta_ >= etaLow_ && eta_ <= etaHigh_) {
482  me500[ii]->Fill(phi_, radLen_);
483  me600[ii]->Fill(phi_, intLen_);
484  me700[ii]->Fill(phi_, stepLen_);
485  me800[ii]->Fill(phi_);
486  }
487 
488  me900[ii]->Fill(eta_, phi_, radLen_);
489  me1000[ii]->Fill(eta_, phi_, intLen_);
490  me1100[ii]->Fill(eta_, phi_, stepLen_);
491  me1200[ii]->Fill(eta_, phi_);
492 
493  if ((std::abs(eta_) >= etaMidMin_) && (std::abs(eta_) <= etaMidMax_)) {
494  me1600[ii]->Fill(phi_, radLen_);
495  me1700[ii]->Fill(phi_, intLen_);
496  me1800[ii]->Fill(phi_, stepLen_);
497  }
498 
499  if ((std::abs(eta_) >= etaHighMin_) && (std::abs(eta_) <= etaHighMax_)) {
500  me1900[ii]->Fill(phi_, radLen_);
501  me2000[ii]->Fill(phi_, intLen_);
502  me2100[ii]->Fill(phi_, stepLen_);
503  }
504 
505  if ((std::abs(eta_) >= etaLowMin_) && (std::abs(eta_) <= etaLowMax_)) {
506  me2200[ii]->Fill(phi_, radLen_);
507  me2300[ii]->Fill(phi_, intLen_);
508  me2400[ii]->Fill(phi_, stepLen_);
509  }
510  }
511 }
512 
513 void MaterialBudgetHcalHistos::fillLayer() {
514  me1300[0]->Fill(eta_);
515  me1400[0]->Fill(eta_, phi_);
516  if (nlayHB_ > 0) {
517  me1300[1]->Fill(eta_);
518  me1400[1]->Fill(eta_, phi_);
519  }
520  if (nlayHB_ >= 16) {
521  me1300[2]->Fill(eta_);
522  me1400[2]->Fill(eta_, phi_);
523  }
524  if (nlayHE_ > 0) {
525  me1300[3]->Fill(eta_);
526  me1400[3]->Fill(eta_, phi_);
527  }
528  if (nlayHE_ >= 16) {
529  me1300[4]->Fill(eta_);
530  me1400[4]->Fill(eta_, phi_);
531  }
532  if (nlayHO_ > 0) {
533  me1300[5]->Fill(eta_);
534  me1400[5]->Fill(eta_, phi_);
535  }
536  if (nlayHO_ >= 2) {
537  me1300[6]->Fill(eta_);
538  me1400[6]->Fill(eta_, phi_);
539  }
540  if (nlayHF_ > 0) {
541  me1300[7]->Fill(eta_);
542  me1400[7]->Fill(eta_, phi_);
543  }
544  if (nlayHB_ > 0 || nlayHE_ > 0 || (nlayHF_ > 0 && std::abs(eta_) > 3.0)) {
545  me1300[8]->Fill(eta_);
546  me1400[8]->Fill(eta_, phi_);
547  }
548  me1500[0]->Fill(eta_, (double)(nlayHB_ + nlayHO_ + nlayHE_ + nlayHF_));
549  me1500[1]->Fill(eta_, (double)(nlayHB_));
550  me1500[2]->Fill(eta_, (double)(nlayHE_));
551  me1500[4]->Fill(eta_, (double)(nlayHF_));
552 }
553 
554 void MaterialBudgetHcalHistos::hend() {
555  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Save user histos ===";
556 }
557 
558 std::vector<std::string> MaterialBudgetHcalHistos::getNames(DDFilteredView& fv) {
559  std::vector<std::string> tmp;
560  bool dodet = fv.firstChild();
561  while (dodet) {
562  const DDLogicalPart& log = fv.logicalPart();
563  std::string namx = log.name().name();
564  if (std::find(tmp.begin(), tmp.end(), namx) == tmp.end())
565  tmp.push_back(namx);
566  dodet = fv.next();
567  }
568  return tmp;
569 }
570 
571 std::vector<std::string> MaterialBudgetHcalHistos::getNames(cms::DDFilteredView& fv) {
572  std::vector<std::string> tmp;
573  const std::vector<std::string> notIn = {
574  "CALO", "HCal", "MBBTL", "MBBTR", "MBBTC", "MBAT", "MBBT_R1M", "MBBT_R1P", "MBBT_R1MX", "MBBT_R1PX", "VCAL"};
575  while (fv.firstChild()) {
576  const std::string n{fv.name().data(), fv.name().size()};
577  if (std::find(notIn.begin(), notIn.end(), n) == notIn.end()) {
578  std::string::size_type pos = n.find(':');
579  const std::string namx = (pos == std::string::npos) ? n : std::string(n, pos + 1, n.size() - 1);
580  if (std::find(tmp.begin(), tmp.end(), namx) == tmp.end())
581  tmp.push_back(namx);
582  }
583  }
584  return tmp;
585 }
586 
587 std::vector<double> MaterialBudgetHcalHistos::getDDDArray(const std::string& str, const DDsvalues_type& sv) {
588  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos:getDDDArray called for " << str;
589  DDValue value(str);
590  if (DDfetch(&sv, value)) {
591  edm::LogVerbatim("MaterialBudgetFull") << value;
592  const std::vector<double>& fvec = value.doubles();
593  int nval = fvec.size();
594  if (nval < 1) {
595  throw cms::Exception("MaterialBudgetHcalHistos") << "nval = " << nval << " < 1 for array " << str << "\n";
596  }
597 
598  return fvec;
599  } else {
600  throw cms::Exception("MaterialBudgetHcalHistos") << "cannot get array " << str << "\n";
601  }
602 }
603 
604 bool MaterialBudgetHcalHistos::isSensitive(const std::string& name) {
605  std::vector<std::string>::const_iterator it = sensitives_.begin();
606  std::vector<std::string>::const_iterator itEnd = sensitives_.end();
607  std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
608  for (; it != itEnd; ++it)
609  if (namx == *it)
610  return true;
611  return false;
612 }
613 
614 bool MaterialBudgetHcalHistos::isItHF(const G4VTouchable* touch) {
615  int levels = ((touch->GetHistoryDepth()) + 1);
616  for (unsigned int it = 0; it < hfNames_.size(); it++) {
617  if (levels >= hfLevels_[it]) {
618  std::string name =
619  (static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(levels - hfLevels_[it])->GetName())))
620  .substr(0, 4);
621  if (name == hfNames_[it]) {
622  return true;
623  }
624  }
625  }
626  return false;
627 }
628 
629 bool MaterialBudgetHcalHistos::isItEC(const std::string& name) {
630  std::vector<std::string>::const_iterator it = sensitiveEC_.begin();
631  std::vector<std::string>::const_iterator itEnd = sensitiveEC_.end();
632  for (; it != itEnd; ++it)
633  if (name.substr(0, 4) == *it)
634  return true;
635  return false;
636 }
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
MaterialBudgetHcalHistos::isItEC
bool isItEC(const std::string &)
Definition: MaterialBudgetHcalHistos.cc:629
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
sistrip::layer_
static const char layer_[]
Definition: ConstantsForGranularity.h:30
DDFilteredView::logicalPart
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
Definition: DDFilteredView.cc:16
DDValue::doubles
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:111
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
mps_fire.i
i
Definition: mps_fire.py:428
DDBase::name
const N & name() const
Definition: DDBase.h:59
angle_units::operators::convertRadToDeg
constexpr NumType convertRadToDeg(NumType radians)
Definition: angle_units.h:21
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
alcazmumu_cfi.filter
tuple filter
Definition: alcazmumu_cfi.py:34
cuy.ii
int ii
Definition: cuy.py:589
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
DDCompactView
Compact representation of the geometrical detector hierarchy.
Definition: DDCompactView.h:81
DDfetch
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:79
trigger::size_type
uint16_t size_type
Definition: TriggerTypeDefs.h:18
names
const std::string names[nVars_]
Definition: PhotonIDValueMapProducer.cc:124
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
geant_units::operators::convertUnitsTo
constexpr NumType convertUnitsTo(double desiredUnits, NumType val)
Definition: GeantUnits.h:100
MaterialBudgetHcalHistos::getNames
std::vector< std::string > getNames(DDFilteredView &fv)
Definition: MaterialBudgetHcalHistos.cc:558
MaterialBudgetHcalHistos::fillBeginJob
void fillBeginJob(const DDCompactView &)
Definition: MaterialBudgetHcalHistos.cc:49
DDsvalues_type
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
Definition: DDsvalues.h:12
MaterialBudgetHcalHistos::fillStartTrack
void fillStartTrack(const G4Track *)
Definition: MaterialBudgetHcalHistos.cc:158
DDFilteredView::next
bool next()
set current node to the next node in the filtered tree
Definition: DDFilteredView.cc:67
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:40
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
MaterialBudgetHcalHistos::isSensitive
bool isSensitive(const std::string &)
Definition: MaterialBudgetHcalHistos.cc:604
DDLogicalPart
A DDLogicalPart aggregates information concerning material, solid and sensitveness ...
Definition: DDLogicalPart.h:93
cms::DDFilteredView::name
std::string_view name() const
Definition: DDFilteredView.cc:857
cms::DDFilteredView::firstChild
bool firstChild()
set the current node to the first child
Definition: DDFilteredView.cc:268
cms::DDFilteredView::get
T get(const std::string &)
extract attribute value
HLT_FULL_cff.maxPhi
tuple maxPhi
Definition: HLT_FULL_cff.py:53293
jetCorrFactors_cfi.levels
list levels
Definition: jetCorrFactors_cfi.py:49
DDFilteredView::mergedSpecifics
DDsvalues_type mergedSpecifics() const
Definition: DDFilteredView.cc:42
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
DDFilteredView::firstChild
bool firstChild()
set the current node to the first child ...
Definition: DDFilteredView.cc:86
MaterialBudgetHcalHistos::MaterialBudgetHcalHistos
MaterialBudgetHcalHistos(const edm::ParameterSet &p)
Definition: MaterialBudgetHcalHistos.cc:20
step
step
Definition: StallMonitor.cc:94
testEve_cfg.level
tuple level
Definition: testEve_cfg.py:47
compare.tfile
tuple tfile
Definition: compare.py:324
MsgTools::getNames
void getNames(uint8 *from, uint32 from_len, Strings &to)
Definition: MsgTools.h:85
str
#define str(s)
Definition: TestProcessor.cc:53
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
MaterialBudgetHcalHistos::getDDDArray
std::vector< double > getDDDArray(const std::string &str, const DDsvalues_type &sv)
Definition: MaterialBudgetHcalHistos.cc:587
DDName::name
const std::string & name() const
Returns the name.
Definition: DDName.cc:41
MaterialBudgetHcalHistos::isItHF
bool isItHF(const G4VTouchable *)
Definition: MaterialBudgetHcalHistos.cc:614
Generated for CMSSW Reference Manual by doxygen 1.8.5