CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
MaterialBudgetHcalAnalysis.cc
Go to the documentation of this file.
1 #include "FWCore/Framework/interface/Frameworkfwd.h"
2 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
3 #include "FWCore/Framework/interface/Event.h"
4 #include "FWCore/Framework/interface/MakerMacros.h"
5 #include "FWCore/Utilities/interface/Exception.h"
6 #include "FWCore/MessageLogger/interface/MessageLogger.h"
7 #include "FWCore/ParameterSet/interface/ParameterSet.h"
8 #include "FWCore/ServiceRegistry/interface/Service.h"
9 #include "CommonTools/UtilAlgos/interface/TFileService.h"
10 
11 #include "DataFormats/Math/interface/GeantUnits.h"
12 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingCalo.h"
13 
14 #include <TH1F.h>
15 #include <TH2F.h>
16 #include <TProfile.h>
17 #include <TProfile2D.h>
18 
19 #include <string>
20 #include <vector>
21 
22 using namespace geant_units::operators;
23 
24 class MaterialBudgetHcalAnalysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
25 public:
26  MaterialBudgetHcalAnalysis(const edm::ParameterSet &p);
27  ~MaterialBudgetHcalAnalysis() override = default;
28 
29  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
30 
31 private:
32  void analyze(edm::Event const &, edm::EventSetup const &) override;
33  void beginJob() override;
34 
35  static const uint32_t maxSet_ = 25, maxSet2_ = 9;
36  const int binEta_, binPhi_;
37  const double maxEta_, etaLow_, etaHigh_, etaLowMin_, etaLowMax_, etaMidMin_;
38  const double etaMidMax_, etaHighMin_, etaHighMax_, etaMinP_, etaMaxP_;
39  const edm::InputTag labelMBCalo_;
40  const edm::EDGetTokenT<MaterialAccountingCaloCollection> tokMBCalo_;
41  TH1F *me400_[maxSet_], *me800_[maxSet_], *me1300_[maxSet2_];
42  TH2F *me1200_[maxSet_], *me1400_[maxSet2_];
43  TProfile *me100_[maxSet_], *me200_[maxSet_], *me300_[maxSet_];
44  TProfile *me500_[maxSet_], *me600_[maxSet_], *me700_[maxSet_];
45  TProfile *me1500_[maxSet2_];
46  TProfile *me1600_[maxSet_], *me1700_[maxSet_], *me1800_[maxSet_];
47  TProfile *me1900_[maxSet_], *me2000_[maxSet_], *me2100_[maxSet_];
48  TProfile *me2200_[maxSet_], *me2300_[maxSet_], *me2400_[maxSet_];
49  TProfile2D *me900_[maxSet_], *me1000_[maxSet_], *me1100_[maxSet_];
50 };
51 
52 MaterialBudgetHcalAnalysis::MaterialBudgetHcalAnalysis(const edm::ParameterSet &p)
53  : binEta_(p.getParameter<int>("nBinEta")),
54  binPhi_(p.getParameter<int>("nBinPhi")),
55  maxEta_(p.getParameter<double>("maxEta")),
56  etaLow_(p.getParameter<double>("etaLow")),
57  etaHigh_(p.getParameter<double>("etaHigh")),
58  etaLowMin_(p.getParameter<double>("etaLowMin")),
59  etaLowMax_(p.getParameter<double>("etaLowMax")),
60  etaMidMin_(p.getParameter<double>("etaMidMin")),
61  etaMidMax_(p.getParameter<double>("etaMidMax")),
62  etaHighMin_(p.getParameter<double>("etaHighMin")),
63  etaHighMax_(p.getParameter<double>("etaHighMax")),
64  etaMinP_(p.getParameter<double>("etaMinP")),
65  etaMaxP_(p.getParameter<double>("etaMaxP")),
66  labelMBCalo_(p.getParameter<edm::InputTag>("labelMBCaloLabel")),
67  tokMBCalo_(consumes<MaterialAccountingCaloCollection>(labelMBCalo_)) {
68  usesResource(TFileService::kSharedResource);
69  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalAnalysis: == Eta plot: NX " << binEta_ << " Range "
70  << -maxEta_ << ":" << maxEta_ << " Phi plot: NX " << binPhi_ << " Range " << -1._pi
71  << ":" << 1._pi << " (Eta limit " << etaLow_ << ":" << etaHigh_ << ")"
72  << " Eta range (" << etaLowMin_ << ":" << etaLowMax_ << "), (" << etaMidMin_ << ":"
73  << etaMidMax_ << "), (" << etaHighMin_ << ":" << etaHighMax_
74  << ") Debug for eta range " << etaMinP_ << ":" << etaMaxP_;
75 }
76 
77 void MaterialBudgetHcalAnalysis::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
78  edm::ParameterSetDescription desc;
79  desc.add<int>("nBinEta", 260);
80  desc.add<int>("nBinPhi", 180);
81  desc.add<double>("maxEta", 5.2);
82  desc.add<double>("etaLow", -5.2);
83  desc.add<double>("etaHigh", 5.2);
84  desc.add<double>("etaMinP", 5.2);
85  desc.add<double>("etaMaxP", 0.0);
86  desc.add<double>("etaLowMin", 0.783);
87  desc.add<double>("etaLowMax", 0.870);
88  desc.add<double>("etaMidMin", 2.650);
89  desc.add<double>("etaMidMax", 2.868);
90  desc.add<double>("etaHighMin", 2.868);
91  desc.add<double>("etaHighMax", 3.000);
92  desc.add<edm::InputTag>("labelMBCaloLabel", edm::InputTag("g4SimHits", "HcalMatBCalo"));
93  descriptions.add("materialBudgetHcalAnalysis", desc);
94 }
95 
96 void MaterialBudgetHcalAnalysis::beginJob() {
97  // Book histograms
98  edm::Service<TFileService> tfile;
99 
100  if (!tfile.isAvailable())
101  throw cms::Exception("BadConfig") << "TFileService unavailable: "
102  << "please add it to config file";
103 
104  double maxPhi = 1._pi;
105  edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalAnalysis: Booking user histos === with " << binEta_
106  << " bins in eta from " << -maxEta_ << " to " << maxEta_ << " and " << binPhi_
107  << " bins in phi from " << -maxPhi << " to " << maxPhi;
108 
109  std::string iter;
110  std::string range0 = "(" + std::to_string(etaMidMin_) + ":" + std::to_string(etaMidMax_) + ") ";
111  std::string range1 = "(" + std::to_string(etaHighMin_) + ":" + std::to_string(etaHighMax_) + ") ";
112  std::string range2 = "(" + std::to_string(etaLowMin_) + ":" + std::to_string(etaLowMax_) + ") ";
113  // total X0
114  for (uint32_t i = 0; i < maxSet_; i++) {
115  iter = std::to_string(i);
116  me100_[i] = tfile->make<TProfile>(
117  std::to_string(i + 100).c_str(), ("MB(X0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
118  me200_[i] = tfile->make<TProfile>(
119  std::to_string(i + 200).c_str(), ("MB(L0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
120  me300_[i] = tfile->make<TProfile>(
121  std::to_string(i + 300).c_str(), ("MB(Step) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
122  me400_[i] = tfile->make<TH1F>(
123  std::to_string(i + 400).c_str(), ("Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
124  me500_[i] = tfile->make<TProfile>(
125  std::to_string(i + 500).c_str(), ("MB(X0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
126  me600_[i] = tfile->make<TProfile>(
127  std::to_string(i + 600).c_str(), ("MB(L0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
128  me700_[i] = tfile->make<TProfile>(
129  std::to_string(i + 700).c_str(), ("MB(Step) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
130  me800_[i] =
131  tfile->make<TH1F>(std::to_string(i + 800).c_str(), ("Phi in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
132  me900_[i] = tfile->make<TProfile2D>(std::to_string(i + 900).c_str(),
133  ("MB(X0) prof Eta Phi in region " + iter).c_str(),
134  binEta_ / 2,
135  -maxEta_,
136  maxEta_,
137  binPhi_ / 2,
138  -maxPhi,
139  maxPhi);
140  me1000_[i] = tfile->make<TProfile2D>(std::to_string(i + 1000).c_str(),
141  ("MB(L0) prof Eta Phi in region " + iter).c_str(),
142  binEta_ / 2,
143  -maxEta_,
144  maxEta_,
145  binPhi_ / 2,
146  -maxPhi,
147  maxPhi);
148  me1100_[i] = tfile->make<TProfile2D>(std::to_string(i + 1100).c_str(),
149  ("MB(Step) prof Eta Phi in region " + iter).c_str(),
150  binEta_ / 2,
151  -maxEta_,
152  maxEta_,
153  binPhi_ / 2,
154  -maxPhi,
155  maxPhi);
156  me1200_[i] = tfile->make<TH2F>(std::to_string(i + 1200).c_str(),
157  ("Eta vs Phi in region " + iter).c_str(),
158  binEta_ / 2,
159  -maxEta_,
160  maxEta_,
161  binPhi_ / 2,
162  -maxPhi,
163  maxPhi);
164  me1600_[i] = tfile->make<TProfile>(std::to_string(i + 1600).c_str(),
165  ("MB(X0) prof Ph in region " + range0 + iter).c_str(),
166  binPhi_,
167  -maxPhi,
168  maxPhi);
169  me1700_[i] = tfile->make<TProfile>(std::to_string(i + 1700).c_str(),
170  ("MB(L0) prof Ph in region " + range0 + iter).c_str(),
171  binPhi_,
172  -maxPhi,
173  maxPhi);
174  me1800_[i] = tfile->make<TProfile>(std::to_string(i + 1800).c_str(),
175  ("MB(Step) prof Ph in region " + range0 + iter).c_str(),
176  binPhi_,
177  -maxPhi,
178  maxPhi);
179  me1900_[i] = tfile->make<TProfile>(std::to_string(i + 1900).c_str(),
180  ("MB(X0) prof Ph in region " + range1 + iter).c_str(),
181  binPhi_,
182  -maxPhi,
183  maxPhi);
184  me2000_[i] = tfile->make<TProfile>(std::to_string(i + 2000).c_str(),
185  ("MB(L0) prof Ph in region " + range1 + iter).c_str(),
186  binPhi_,
187  -maxPhi,
188  maxPhi);
189  me2100_[i] = tfile->make<TProfile>(std::to_string(i + 2100).c_str(),
190  ("MB(Step) prof Ph in region " + range1 + iter).c_str(),
191  binPhi_,
192  -maxPhi,
193  maxPhi);
194  me2200_[i] = tfile->make<TProfile>(std::to_string(i + 2200).c_str(),
195  ("MB(X0) prof Ph in region " + range2 + iter).c_str(),
196  binPhi_,
197  -maxPhi,
198  maxPhi);
199  me2300_[i] = tfile->make<TProfile>(std::to_string(i + 2300).c_str(),
200  ("MB(L0) prof Ph in region " + range2 + iter).c_str(),
201  binPhi_,
202  -maxPhi,
203  maxPhi);
204  me2400_[i] = tfile->make<TProfile>(std::to_string(i + 2400).c_str(),
205  ("MB(Step) prof Ph in region " + range2 + iter).c_str(),
206  binPhi_,
207  -maxPhi,
208  maxPhi);
209  }
210  for (uint32_t i = 0; i < maxSet2_; i++) {
211  iter = std::to_string(i);
212  me1300_[i] = tfile->make<TH1F>(std::to_string(i + 1300).c_str(),
213  ("Events with layers Hit (0 all, 1 HB, ..) for " + iter).c_str(),
214  binEta_,
215  -maxEta_,
216  maxEta_);
217  me1400_[i] = tfile->make<TH2F>(std::to_string(i + 1400).c_str(),
218  ("Eta vs Phi for layers hit in " + iter).c_str(),
219  binEta_ / 2,
220  -maxEta_,
221  maxEta_,
222  binPhi_ / 2,
223  -maxPhi,
224  maxPhi);
225  me1500_[i] = tfile->make<TProfile>(std::to_string(i + 1500).c_str(),
226  ("Number of layers crossed (0 all, 1 HB, ..) for " + iter).c_str(),
227  binEta_,
228  -maxEta_,
229  maxEta_);
230  }
231 
232  edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalAnalysis: Booking user histos done ===";
233 }
234 
235 void MaterialBudgetHcalAnalysis::analyze(edm::Event const &iEvent, edm::EventSetup const &iSetup) {
236 #ifdef EDM_ML_DEBUG
237  edm::LogVerbatim("MaterialBudgetFull") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event()
238  << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
239  << iEvent.bunchCrossing();
240 #endif
241 
242  // Fill from the MB collection
243  auto const &hcalMBColl = iEvent.getHandle(tokMBCalo_);
244  if (hcalMBColl.isValid()) {
245  auto hcalMB = hcalMBColl.product();
246 #ifdef EDM_ML_DEBUG
247  edm::LogVerbatim("MaterialBudgetFull")
248  << "Finds HcalMaterialBudgetCollection with " << hcalMB->size() << " entries";
249 #endif
250 
251  for (auto itr = hcalMB->begin(); itr != hcalMB->end(); ++itr) {
252  for (uint32_t ii = 0; ii < itr->m_stepLen.size(); ++ii) {
253 #ifdef EDM_ML_DEBUG
254  if ((std::abs(itr->m_eta) >= etaMinP_) && (std::abs(itr->m_eta) <= etaMaxP_))
255  edm::LogVerbatim("MaterialBudget")
256  << "MaterialBudgetHcalAnalysis:FillHisto called with index " << ii << " integrated step "
257  << itr->m_stepLen[ii] << " X0 " << itr->m_radLen[ii] << " Lamda " << itr->m_intLen[ii];
258 #endif
259  if (ii < maxSet_) {
260  me100_[ii]->Fill(itr->m_eta, itr->m_radLen[ii]);
261  me200_[ii]->Fill(itr->m_eta, itr->m_intLen[ii]);
262  me300_[ii]->Fill(itr->m_eta, itr->m_stepLen[ii]);
263  me400_[ii]->Fill(itr->m_eta);
264 
265  if (itr->m_eta >= etaLow_ && itr->m_eta <= etaHigh_) {
266  me500_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
267  me600_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
268  me700_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
269  me800_[ii]->Fill(itr->m_phi);
270  }
271 
272  me900_[ii]->Fill(itr->m_eta, itr->m_phi, itr->m_radLen[ii]);
273  me1000_[ii]->Fill(itr->m_eta, itr->m_phi, itr->m_intLen[ii]);
274  me1100_[ii]->Fill(itr->m_eta, itr->m_phi, itr->m_stepLen[ii]);
275  me1200_[ii]->Fill(itr->m_eta, itr->m_phi);
276 
277  if ((std::abs(itr->m_eta) >= etaMidMin_) && (std::abs(itr->m_eta) <= etaMidMax_)) {
278  me1600_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
279  me1700_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
280  me1800_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
281  }
282 
283  if ((std::abs(itr->m_eta) >= etaHighMin_) && (std::abs(itr->m_eta) <= etaHighMax_)) {
284  me1900_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
285  me2000_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
286  me2100_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
287  }
288 
289  if ((std::abs(itr->m_eta) >= etaLowMin_) && (std::abs(itr->m_eta) <= etaLowMax_)) {
290  me2200_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
291  me2300_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
292  me2400_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
293  }
294  }
295  }
296 
297  me1300_[0]->Fill(itr->m_eta);
298  me1400_[0]->Fill(itr->m_eta, itr->m_phi);
299  if (itr->m_layers[0] > 0) {
300  me1300_[1]->Fill(itr->m_eta);
301  me1400_[1]->Fill(itr->m_eta, itr->m_phi);
302  }
303  if (itr->m_layers[0] >= 16) {
304  me1300_[2]->Fill(itr->m_eta);
305  me1400_[2]->Fill(itr->m_eta, itr->m_phi);
306  }
307  if (itr->m_layers[1] > 0) {
308  me1300_[3]->Fill(itr->m_eta);
309  me1400_[3]->Fill(itr->m_eta, itr->m_phi);
310  }
311  if (itr->m_layers[1] >= 16) {
312  me1300_[4]->Fill(itr->m_eta);
313  me1400_[4]->Fill(itr->m_eta, itr->m_phi);
314  }
315  if (itr->m_layers[2] > 0) {
316  me1300_[5]->Fill(itr->m_eta);
317  me1400_[5]->Fill(itr->m_eta, itr->m_phi);
318  }
319  if (itr->m_layers[2] >= 2) {
320  me1300_[6]->Fill(itr->m_eta);
321  me1400_[6]->Fill(itr->m_eta, itr->m_phi);
322  }
323  if (itr->m_layers[3] > 0) {
324  me1300_[7]->Fill(itr->m_eta);
325  me1400_[7]->Fill(itr->m_eta, itr->m_phi);
326  }
327  if (itr->m_layers[0] > 0 || itr->m_layers[1] > 0 || (itr->m_layers[3] > 0 && std::abs(itr->m_eta) > 3.0)) {
328  me1300_[8]->Fill(itr->m_eta);
329  me1400_[8]->Fill(itr->m_eta, itr->m_phi);
330  }
331  me1500_[0]->Fill(itr->m_eta, (double)(itr->m_layers[0] + itr->m_layers[1] + itr->m_layers[2] + itr->m_layers[3]));
332  me1500_[1]->Fill(itr->m_eta, (double)(itr->m_layers[0]));
333  me1500_[2]->Fill(itr->m_eta, (double)(itr->m_layers[1]));
334  me1500_[4]->Fill(itr->m_eta, (double)(itr->m_layers[3]));
335  }
336  }
337 }
338 
339 //define this as a plug-in
340 DEFINE_FWK_MODULE(MaterialBudgetHcalAnalysis);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:38
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:40
MaterialBudgetHcalAnalysis::tokMBCalo_
const edm::EDGetTokenT< MaterialAccountingCaloCollection > tokMBCalo_
Definition: MaterialBudgetHcalAnalysis.cc:40
mps_fire.i
i
Definition: mps_fire.py:428
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
MaterialAccountingCaloCollection
std::vector< MaterialAccountingCalo > MaterialAccountingCaloCollection
Definition: MaterialAccountingCalo.h:25
edm::EventBase::bunchCrossing
int bunchCrossing() const
Definition: EventBase.h:64
MaterialBudgetHcalAnalysis::analyze
void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: MaterialBudgetHcalAnalysis.cc:235
edm::EventBase::luminosityBlock
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:61
cond::impl::to_string
std::string to_string(const V &value)
Definition: OMSAccess.h:71
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
cuy.ii
int ii
Definition: cuy.py:589
analyze
example_stream void analyze(const edm::Event &, const edm::EventSetup &) override
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
bk::beginJob
void beginJob()
Definition: Breakpoints.cc:14
MaterialBudgetHcalAnalysis::MaterialBudgetHcalAnalysis
MaterialBudgetHcalAnalysis(const edm::ParameterSet &p)
Definition: MaterialBudgetHcalAnalysis.cc:52
edm::Event::getHandle
Handle< PROD > getHandle(EDGetTokenT< PROD > token) const
Definition: Event.h:563
iEvent
int iEvent
Definition: GenABIO.cc:224
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:40
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::ParameterSetDescription::add
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:95
fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
HLT_FULL_cff.maxPhi
tuple maxPhi
Definition: HLT_FULL_cff.py:39863
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HLT_FULL_cff.InputTag
tuple InputTag
Definition: HLT_FULL_cff.py:72505
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
MaterialBudgetHcalAnalysis::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: MaterialBudgetHcalAnalysis.cc:77
compare.tfile
tuple tfile
Definition: compare.py:324
Generated for CMSSW Reference Manual by doxygen 1.8.5