CMS 3D CMS Logo
DQMSourcePi0.cc
Go to the documentation of this file.
1 #include "FWCore/Framework/interface/Event.h"
2 #include "FWCore/Framework/interface/EventSetup.h"
3 #include "FWCore/Framework/interface/Frameworkfwd.h"
4 #include "FWCore/Framework/interface/MakerMacros.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 #include "FWCore/ParameterSet/interface/ParameterSet.h"
7 #include "FWCore/ServiceRegistry/interface/Service.h"
8 #include "FWCore/Utilities/interface/InputTag.h"
9 
10 // DQM include files
11 
12 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
13 #include "DQMServices/Core/interface/DQMStore.h"
14 
15 // work on collections
16 #include "DataFormats/DetId/interface/DetId.h"
17 #include "DataFormats/EcalDetId/interface/EBDetId.h"
18 #include "DataFormats/EcalDetId/interface/EEDetId.h"
19 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
20 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
21 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
22 #include "DataFormats/Math/interface/Point3D.h"
23 #include "DataFormats/Common/interface/Handle.h"
24 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
25 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
26 
27 // Geometry
28 #include "Geometry/Records/interface/CaloTopologyRecord.h"
29 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
30 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
31 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
32 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
33 #include "Geometry/CaloTopology/interface/CaloTopology.h"
34 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
35 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
36 #include "Geometry/Records/interface/IdealGeometryRecord.h"
37 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
38 
39 #include "TVector3.h"
40 
41 // Less than operator for sorting EcalRecHits according to energy.
42 inline bool ecalRecHitGreater(EcalRecHit x, EcalRecHit y) { return (x.energy() > y.energy()); }
43 
44 class DQMSourcePi0 : public DQMEDAnalyzer {
45 public:
46  DQMSourcePi0(const edm::ParameterSet &);
47  ~DQMSourcePi0() override;
48 
49 protected:
50  void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
51  void analyze(const edm::Event &e, const edm::EventSetup &c) override;
52 
53  void convxtalid(int &, int &);
54  int diff_neta_s(int, int);
55  int diff_nphi_s(int, int);
56 
57 private:
58  int eventCounter_;
59  PositionCalc posCalculator_;
60 
62  MonitorElement *hiPhiDistrEBpi0_;
63 
65  MonitorElement *hiXDistrEEpi0_;
66 
68  MonitorElement *hiPhiDistrEBeta_;
69 
71  MonitorElement *hiXDistrEEeta_;
72 
74  MonitorElement *hiEtaDistrEBpi0_;
75 
77  MonitorElement *hiYDistrEEpi0_;
78 
80  MonitorElement *hiEtaDistrEBeta_;
81 
83  MonitorElement *hiYDistrEEeta_;
84 
86  MonitorElement *hRechitEnergyEBpi0_;
87 
89  MonitorElement *hRechitEnergyEEpi0_;
90 
92  MonitorElement *hRechitEnergyEBeta_;
93 
95  MonitorElement *hRechitEnergyEEeta_;
96 
98  MonitorElement *hEventEnergyEBpi0_;
99 
101  MonitorElement *hEventEnergyEEpi0_;
102 
104  MonitorElement *hEventEnergyEBeta_;
105 
107  MonitorElement *hEventEnergyEEeta_;
108 
110  MonitorElement *hNRecHitsEBpi0_;
111 
113  MonitorElement *hNRecHitsEEpi0_;
114 
116  MonitorElement *hNRecHitsEBeta_;
117 
119  MonitorElement *hNRecHitsEEeta_;
120 
122  MonitorElement *hMeanRecHitEnergyEBpi0_;
123 
125  MonitorElement *hMeanRecHitEnergyEEpi0_;
126 
128  MonitorElement *hMeanRecHitEnergyEBeta_;
129 
131  MonitorElement *hMeanRecHitEnergyEEeta_;
132 
134  MonitorElement *hMinvPi0EB_;
135 
137  MonitorElement *hMinvPi0EE_;
138 
140  MonitorElement *hMinvEtaEB_;
141 
143  MonitorElement *hMinvEtaEE_;
144 
146  MonitorElement *hPt1Pi0EB_;
147 
149  MonitorElement *hPt1Pi0EE_;
150 
152  MonitorElement *hPt1EtaEB_;
153 
155  MonitorElement *hPt1EtaEE_;
156 
158  MonitorElement *hPt2Pi0EB_;
159 
161  MonitorElement *hPt2Pi0EE_;
162 
164  MonitorElement *hPt2EtaEB_;
165 
167  MonitorElement *hPt2EtaEE_;
168 
170  MonitorElement *hPtPi0EB_;
171 
173  MonitorElement *hPtPi0EE_;
174 
176  MonitorElement *hPtEtaEB_;
177 
179  MonitorElement *hPtEtaEE_;
180 
182  MonitorElement *hIsoPi0EB_;
183 
185  MonitorElement *hIsoPi0EE_;
186 
188  MonitorElement *hIsoEtaEB_;
189 
191  MonitorElement *hIsoEtaEE_;
192 
194  MonitorElement *hS4S91Pi0EB_;
195 
197  MonitorElement *hS4S91Pi0EE_;
198 
200  MonitorElement *hS4S91EtaEB_;
201 
203  MonitorElement *hS4S91EtaEE_;
204 
206  MonitorElement *hS4S92Pi0EB_;
207 
209  MonitorElement *hS4S92Pi0EE_;
210 
212  MonitorElement *hS4S92EtaEB_;
213 
215  MonitorElement *hS4S92EtaEE_;
216 
218  edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEBpi0_;
219  edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEBeta_;
220 
222  edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEEpi0_;
223  edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEEeta_;
224 
225  edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopoToken_;
226  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
227 
228  int gammaCandEtaSize_;
229  int gammaCandPhiSize_;
230 
231  double seleXtalMinEnergy_;
232  double seleXtalMinEnergyEndCap_;
233 
234  double clusSeedThr_;
235  int clusEtaSize_;
236  int clusPhiSize_;
237 
238  double clusSeedThrEndCap_;
239 
241  double selePtGamma_;
242  double selePtPi0_;
243  double seleMinvMaxPi0_;
244  double seleMinvMinPi0_;
245  double seleS4S9Gamma_;
246  double selePi0BeltDR_;
247  double selePi0BeltDeta_;
248  double selePi0Iso_;
249  double ptMinForIsolation_;
250 
252  double selePtGammaEndCap_;
253  double selePtPi0EndCap_;
254  double seleMinvMaxPi0EndCap_;
255  double seleMinvMinPi0EndCap_;
256  double seleS4S9GammaEndCap_;
257  double selePi0IsoEndCap_;
258  double selePi0BeltDREndCap_;
259  double selePi0BeltDetaEndCap_;
260  double ptMinForIsolationEndCap_;
261 
263  double selePtGammaEta_;
264  double selePtEta_;
265  double seleS4S9GammaEta_;
266  double seleS9S25GammaEta_;
267  double seleMinvMaxEta_;
268  double seleMinvMinEta_;
269  double ptMinForIsolationEta_;
270  double seleEtaIso_;
271  double seleEtaBeltDR_;
272  double seleEtaBeltDeta_;
273 
275  double selePtGammaEtaEndCap_;
276  double seleS4S9GammaEtaEndCap_;
277  double seleS9S25GammaEtaEndCap_;
278  double selePtEtaEndCap_;
279  double seleMinvMaxEtaEndCap_;
280  double seleMinvMinEtaEndCap_;
281  double ptMinForIsolationEtaEndCap_;
282  double seleEtaIsoEndCap_;
283  double seleEtaBeltDREndCap_;
284  double seleEtaBeltDetaEndCap_;
285 
286  bool ParameterLogWeighted_;
287  double ParameterX0_;
288  double ParameterT0_barl_;
289  double ParameterT0_endc_;
290  double ParameterT0_endcPresh_;
291  double ParameterW0_;
292 
293  std::vector<EBDetId> detIdEBRecHits;
294  std::vector<EcalRecHit> EBRecHits;
295 
296  std::vector<EEDetId> detIdEERecHits;
297  std::vector<EcalRecHit> EERecHits;
298 
300  unsigned int prescaleFactor_;
301 
303  std::string folderName_;
304 
306  bool saveToFile_;
307 
309  bool isMonEBpi0_;
310  bool isMonEBeta_;
311  bool isMonEEpi0_;
312  bool isMonEEeta_;
313 
315  std::string fileName_;
316 };
317 
318 #define TWOPI 6.283185308
319 
320 // ******************************************
321 // constructors
322 // *****************************************
323 
324 DQMSourcePi0::DQMSourcePi0(const edm::ParameterSet &ps) : eventCounter_(0) {
325  folderName_ = ps.getUntrackedParameter<std::string>("FolderName", "HLT/AlCaEcalPi0");
326  prescaleFactor_ = ps.getUntrackedParameter<int>("prescaleFactor", 1);
327  productMonitoredEBpi0_ =
328  consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEBpi0Tag"));
329  productMonitoredEBeta_ =
330  consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEBetaTag"));
331  productMonitoredEEpi0_ =
332  consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEEpi0Tag"));
333  productMonitoredEEeta_ =
334  consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEEetaTag"));
335  caloTopoToken_ = esConsumes();
336  caloGeomToken_ = esConsumes();
337 
338  isMonEBpi0_ = ps.getUntrackedParameter<bool>("isMonEBpi0", false);
339  isMonEBeta_ = ps.getUntrackedParameter<bool>("isMonEBeta", false);
340  isMonEEpi0_ = ps.getUntrackedParameter<bool>("isMonEEpi0", false);
341  isMonEEeta_ = ps.getUntrackedParameter<bool>("isMonEEeta", false);
342 
343  saveToFile_ = ps.getUntrackedParameter<bool>("SaveToFile", false);
344  fileName_ = ps.getUntrackedParameter<std::string>("FileName", "MonitorAlCaEcalPi0.root");
345 
346  clusSeedThr_ = ps.getParameter<double>("clusSeedThr");
347  clusSeedThrEndCap_ = ps.getParameter<double>("clusSeedThrEndCap");
348  clusEtaSize_ = ps.getParameter<int>("clusEtaSize");
349  clusPhiSize_ = ps.getParameter<int>("clusPhiSize");
350  if (clusPhiSize_ % 2 == 0 || clusEtaSize_ % 2 == 0)
351  edm::LogError("AlCaPi0RecHitsProducerError") << "Size of eta/phi for simple clustering should be odd numbers";
352 
353  seleXtalMinEnergy_ = ps.getParameter<double>("seleXtalMinEnergy");
354  seleXtalMinEnergyEndCap_ = ps.getParameter<double>("seleXtalMinEnergyEndCap");
355 
357  selePtGamma_ = ps.getParameter<double>("selePtGamma");
358  selePtPi0_ = ps.getParameter<double>("selePtPi0");
359  seleMinvMaxPi0_ = ps.getParameter<double>("seleMinvMaxPi0");
360  seleMinvMinPi0_ = ps.getParameter<double>("seleMinvMinPi0");
361  seleS4S9Gamma_ = ps.getParameter<double>("seleS4S9Gamma");
362  selePi0Iso_ = ps.getParameter<double>("selePi0Iso");
363  ptMinForIsolation_ = ps.getParameter<double>("ptMinForIsolation");
364  selePi0BeltDR_ = ps.getParameter<double>("selePi0BeltDR");
365  selePi0BeltDeta_ = ps.getParameter<double>("selePi0BeltDeta");
366 
368  selePtGammaEndCap_ = ps.getParameter<double>("selePtGammaEndCap");
369  selePtPi0EndCap_ = ps.getParameter<double>("selePtPi0EndCap");
370  seleS4S9GammaEndCap_ = ps.getParameter<double>("seleS4S9GammaEndCap");
371  seleMinvMaxPi0EndCap_ = ps.getParameter<double>("seleMinvMaxPi0EndCap");
372  seleMinvMinPi0EndCap_ = ps.getParameter<double>("seleMinvMinPi0EndCap");
373  ptMinForIsolationEndCap_ = ps.getParameter<double>("ptMinForIsolationEndCap");
374  selePi0BeltDREndCap_ = ps.getParameter<double>("selePi0BeltDREndCap");
375  selePi0BeltDetaEndCap_ = ps.getParameter<double>("selePi0BeltDetaEndCap");
376  selePi0IsoEndCap_ = ps.getParameter<double>("selePi0IsoEndCap");
377 
379  selePtGammaEta_ = ps.getParameter<double>("selePtGammaEta");
380  selePtEta_ = ps.getParameter<double>("selePtEta");
381  seleS4S9GammaEta_ = ps.getParameter<double>("seleS4S9GammaEta");
382  seleS9S25GammaEta_ = ps.getParameter<double>("seleS9S25GammaEta");
383  seleMinvMaxEta_ = ps.getParameter<double>("seleMinvMaxEta");
384  seleMinvMinEta_ = ps.getParameter<double>("seleMinvMinEta");
385  ptMinForIsolationEta_ = ps.getParameter<double>("ptMinForIsolationEta");
386  seleEtaIso_ = ps.getParameter<double>("seleEtaIso");
387  seleEtaBeltDR_ = ps.getParameter<double>("seleEtaBeltDR");
388  seleEtaBeltDeta_ = ps.getParameter<double>("seleEtaBeltDeta");
389 
391  selePtGammaEtaEndCap_ = ps.getParameter<double>("selePtGammaEtaEndCap");
392  selePtEtaEndCap_ = ps.getParameter<double>("selePtEtaEndCap");
393  seleS4S9GammaEtaEndCap_ = ps.getParameter<double>("seleS4S9GammaEtaEndCap");
394  seleS9S25GammaEtaEndCap_ = ps.getParameter<double>("seleS9S25GammaEtaEndCap");
395  seleMinvMaxEtaEndCap_ = ps.getParameter<double>("seleMinvMaxEtaEndCap");
396  seleMinvMinEtaEndCap_ = ps.getParameter<double>("seleMinvMinEtaEndCap");
397  ptMinForIsolationEtaEndCap_ = ps.getParameter<double>("ptMinForIsolationEtaEndCap");
398  seleEtaIsoEndCap_ = ps.getParameter<double>("seleEtaIsoEndCap");
399  seleEtaBeltDREndCap_ = ps.getParameter<double>("seleEtaBeltDREndCap");
400  seleEtaBeltDetaEndCap_ = ps.getParameter<double>("seleEtaBeltDetaEndCap");
401 
402  // Parameters for the position calculation:
403  edm::ParameterSet posCalcParameters = ps.getParameter<edm::ParameterSet>("posCalcParameters");
404  posCalculator_ = PositionCalc(posCalcParameters);
405 }
406 
407 DQMSourcePi0::~DQMSourcePi0() {}
408 
409 //--------------------------------------------------------
410 void DQMSourcePi0::bookHistograms(DQMStore::IBooker &ibooker, edm::Run const &irun, edm::EventSetup const &isetup) {
411  // create and cd into new folder
412  ibooker.setCurrentFolder(folderName_);
413 
414  // book some histograms 1D
415 
416  hiPhiDistrEBpi0_ = ibooker.book1D("iphiDistributionEBpi0", "RechitEB pi0 iphi", 361, 1, 361);
417  hiPhiDistrEBpi0_->setAxisTitle("i#phi ", 1);
418  hiPhiDistrEBpi0_->setAxisTitle("# rechits", 2);
419 
420  hiXDistrEEpi0_ = ibooker.book1D("iXDistributionEEpi0", "RechitEE pi0 ix", 100, 0, 100);
421  hiXDistrEEpi0_->setAxisTitle("ix ", 1);
422  hiXDistrEEpi0_->setAxisTitle("# rechits", 2);
423 
424  hiPhiDistrEBeta_ = ibooker.book1D("iphiDistributionEBeta", "RechitEB eta iphi", 361, 1, 361);
425  hiPhiDistrEBeta_->setAxisTitle("i#phi ", 1);
426  hiPhiDistrEBeta_->setAxisTitle("# rechits", 2);
427 
428  hiXDistrEEeta_ = ibooker.book1D("iXDistributionEEeta", "RechitEE eta ix", 100, 0, 100);
429  hiXDistrEEeta_->setAxisTitle("ix ", 1);
430  hiXDistrEEeta_->setAxisTitle("# rechits", 2);
431 
432  hiEtaDistrEBpi0_ = ibooker.book1D("iEtaDistributionEBpi0", "RechitEB pi0 ieta", 171, -85, 86);
433  hiEtaDistrEBpi0_->setAxisTitle("i#eta", 1);
434  hiEtaDistrEBpi0_->setAxisTitle("#rechits", 2);
435 
436  hiYDistrEEpi0_ = ibooker.book1D("iYDistributionEEpi0", "RechitEE pi0 iY", 100, 0, 100);
437  hiYDistrEEpi0_->setAxisTitle("iy", 1);
438  hiYDistrEEpi0_->setAxisTitle("#rechits", 2);
439 
440  hiEtaDistrEBeta_ = ibooker.book1D("iEtaDistributionEBeta", "RechitEB eta ieta", 171, -85, 86);
441  hiEtaDistrEBeta_->setAxisTitle("i#eta", 1);
442  hiEtaDistrEBeta_->setAxisTitle("#rechits", 2);
443 
444  hiYDistrEEeta_ = ibooker.book1D("iYDistributionEEeta", "RechitEE eta iY", 100, 0, 100);
445  hiYDistrEEeta_->setAxisTitle("iy", 1);
446  hiYDistrEEeta_->setAxisTitle("#rechits", 2);
447 
448  hRechitEnergyEBpi0_ = ibooker.book1D("rhEnergyEBpi0", "Pi0 rechits energy EB", 160, 0., 2.0);
449  hRechitEnergyEBpi0_->setAxisTitle("energy (GeV) ", 1);
450  hRechitEnergyEBpi0_->setAxisTitle("#rechits", 2);
451 
452  hRechitEnergyEEpi0_ = ibooker.book1D("rhEnergyEEpi0", "Pi0 rechits energy EE", 160, 0., 3.0);
453  hRechitEnergyEEpi0_->setAxisTitle("energy (GeV) ", 1);
454  hRechitEnergyEEpi0_->setAxisTitle("#rechits", 2);
455 
456  hRechitEnergyEBeta_ = ibooker.book1D("rhEnergyEBeta", "Eta rechits energy EB", 160, 0., 2.0);
457  hRechitEnergyEBeta_->setAxisTitle("energy (GeV) ", 1);
458  hRechitEnergyEBeta_->setAxisTitle("#rechits", 2);
459 
460  hRechitEnergyEEeta_ = ibooker.book1D("rhEnergyEEeta", "Eta rechits energy EE", 160, 0., 3.0);
461  hRechitEnergyEEeta_->setAxisTitle("energy (GeV) ", 1);
462  hRechitEnergyEEeta_->setAxisTitle("#rechits", 2);
463 
464  hEventEnergyEBpi0_ = ibooker.book1D("eventEnergyEBpi0", "Pi0 event energy EB", 100, 0., 20.0);
465  hEventEnergyEBpi0_->setAxisTitle("energy (GeV) ", 1);
466 
467  hEventEnergyEEpi0_ = ibooker.book1D("eventEnergyEEpi0", "Pi0 event energy EE", 100, 0., 50.0);
468  hEventEnergyEEpi0_->setAxisTitle("energy (GeV) ", 1);
469 
470  hEventEnergyEBeta_ = ibooker.book1D("eventEnergyEBeta", "Eta event energy EB", 100, 0., 20.0);
471  hEventEnergyEBeta_->setAxisTitle("energy (GeV) ", 1);
472 
473  hEventEnergyEEeta_ = ibooker.book1D("eventEnergyEEeta", "Eta event energy EE", 100, 0., 50.0);
474  hEventEnergyEEeta_->setAxisTitle("energy (GeV) ", 1);
475 
476  hNRecHitsEBpi0_ = ibooker.book1D("nRechitsEBpi0", "#rechits in pi0 collection EB", 100, 0., 250.);
477  hNRecHitsEBpi0_->setAxisTitle("rechits ", 1);
478 
479  hNRecHitsEEpi0_ = ibooker.book1D("nRechitsEEpi0", "#rechits in pi0 collection EE", 100, 0., 250.);
480  hNRecHitsEEpi0_->setAxisTitle("rechits ", 1);
481 
482  hNRecHitsEBeta_ = ibooker.book1D("nRechitsEBeta", "#rechits in eta collection EB", 100, 0., 250.);
483  hNRecHitsEBeta_->setAxisTitle("rechits ", 1);
484 
485  hNRecHitsEEeta_ = ibooker.book1D("nRechitsEEeta", "#rechits in eta collection EE", 100, 0., 250.);
486  hNRecHitsEEeta_->setAxisTitle("rechits ", 1);
487 
488  hMeanRecHitEnergyEBpi0_ = ibooker.book1D("meanEnergyEBpi0", "Mean rechit energy in pi0 collection EB", 50, 0., 2.);
489  hMeanRecHitEnergyEBpi0_->setAxisTitle("Mean Energy [GeV] ", 1);
490 
491  hMeanRecHitEnergyEEpi0_ = ibooker.book1D("meanEnergyEEpi0", "Mean rechit energy in pi0 collection EE", 100, 0., 5.);
492  hMeanRecHitEnergyEEpi0_->setAxisTitle("Mean Energy [GeV] ", 1);
493 
494  hMeanRecHitEnergyEBeta_ = ibooker.book1D("meanEnergyEBeta", "Mean rechit energy in eta collection EB", 50, 0., 2.);
495  hMeanRecHitEnergyEBeta_->setAxisTitle("Mean Energy [GeV] ", 1);
496 
497  hMeanRecHitEnergyEEeta_ = ibooker.book1D("meanEnergyEEeta", "Mean rechit energy in eta collection EE", 100, 0., 5.);
498  hMeanRecHitEnergyEEeta_->setAxisTitle("Mean Energy [GeV] ", 1);
499 
500  hMinvPi0EB_ = ibooker.book1D("Pi0InvmassEB", "Pi0 Invariant Mass in EB", 100, 0., 0.5);
501  hMinvPi0EB_->setAxisTitle("Inv Mass [GeV] ", 1);
502 
503  hMinvPi0EE_ = ibooker.book1D("Pi0InvmassEE", "Pi0 Invariant Mass in EE", 100, 0., 0.5);
504  hMinvPi0EE_->setAxisTitle("Inv Mass [GeV] ", 1);
505 
506  hMinvEtaEB_ = ibooker.book1D("EtaInvmassEB", "Eta Invariant Mass in EB", 100, 0., 0.85);
507  hMinvEtaEB_->setAxisTitle("Inv Mass [GeV] ", 1);
508 
509  hMinvEtaEE_ = ibooker.book1D("EtaInvmassEE", "Eta Invariant Mass in EE", 100, 0., 0.85);
510  hMinvEtaEE_->setAxisTitle("Inv Mass [GeV] ", 1);
511 
512  hPt1Pi0EB_ = ibooker.book1D("Pt1Pi0EB", "Pt 1st most energetic Pi0 photon in EB", 100, 0., 20.);
513  hPt1Pi0EB_->setAxisTitle("1st photon Pt [GeV] ", 1);
514 
515  hPt1Pi0EE_ = ibooker.book1D("Pt1Pi0EE", "Pt 1st most energetic Pi0 photon in EE", 100, 0., 20.);
516  hPt1Pi0EE_->setAxisTitle("1st photon Pt [GeV] ", 1);
517 
518  hPt1EtaEB_ = ibooker.book1D("Pt1EtaEB", "Pt 1st most energetic Eta photon in EB", 100, 0., 20.);
519  hPt1EtaEB_->setAxisTitle("1st photon Pt [GeV] ", 1);
520 
521  hPt1EtaEE_ = ibooker.book1D("Pt1EtaEE", "Pt 1st most energetic Eta photon in EE", 100, 0., 20.);
522  hPt1EtaEE_->setAxisTitle("1st photon Pt [GeV] ", 1);
523 
524  hPt2Pi0EB_ = ibooker.book1D("Pt2Pi0EB", "Pt 2nd most energetic Pi0 photon in EB", 100, 0., 20.);
525  hPt2Pi0EB_->setAxisTitle("2nd photon Pt [GeV] ", 1);
526 
527  hPt2Pi0EE_ = ibooker.book1D("Pt2Pi0EE", "Pt 2nd most energetic Pi0 photon in EE", 100, 0., 20.);
528  hPt2Pi0EE_->setAxisTitle("2nd photon Pt [GeV] ", 1);
529 
530  hPt2EtaEB_ = ibooker.book1D("Pt2EtaEB", "Pt 2nd most energetic Eta photon in EB", 100, 0., 20.);
531  hPt2EtaEB_->setAxisTitle("2nd photon Pt [GeV] ", 1);
532 
533  hPt2EtaEE_ = ibooker.book1D("Pt2EtaEE", "Pt 2nd most energetic Eta photon in EE", 100, 0., 20.);
534  hPt2EtaEE_->setAxisTitle("2nd photon Pt [GeV] ", 1);
535 
536  hPtPi0EB_ = ibooker.book1D("PtPi0EB", "Pi0 Pt in EB", 100, 0., 20.);
537  hPtPi0EB_->setAxisTitle("Pi0 Pt [GeV] ", 1);
538 
539  hPtPi0EE_ = ibooker.book1D("PtPi0EE", "Pi0 Pt in EE", 100, 0., 20.);
540  hPtPi0EE_->setAxisTitle("Pi0 Pt [GeV] ", 1);
541 
542  hPtEtaEB_ = ibooker.book1D("PtEtaEB", "Eta Pt in EB", 100, 0., 20.);
543  hPtEtaEB_->setAxisTitle("Eta Pt [GeV] ", 1);
544 
545  hPtEtaEE_ = ibooker.book1D("PtEtaEE", "Eta Pt in EE", 100, 0., 20.);
546  hPtEtaEE_->setAxisTitle("Eta Pt [GeV] ", 1);
547 
548  hIsoPi0EB_ = ibooker.book1D("IsoPi0EB", "Pi0 Iso in EB", 50, 0., 1.);
549  hIsoPi0EB_->setAxisTitle("Pi0 Iso", 1);
550 
551  hIsoPi0EE_ = ibooker.book1D("IsoPi0EE", "Pi0 Iso in EE", 50, 0., 1.);
552  hIsoPi0EE_->setAxisTitle("Pi0 Iso", 1);
553 
554  hIsoEtaEB_ = ibooker.book1D("IsoEtaEB", "Eta Iso in EB", 50, 0., 1.);
555  hIsoEtaEB_->setAxisTitle("Eta Iso", 1);
556 
557  hIsoEtaEE_ = ibooker.book1D("IsoEtaEE", "Eta Iso in EE", 50, 0., 1.);
558  hIsoEtaEE_->setAxisTitle("Eta Iso", 1);
559 
560  hS4S91Pi0EB_ = ibooker.book1D("S4S91Pi0EB", "S4S9 1st most energetic Pi0 photon in EB", 50, 0., 1.);
561  hS4S91Pi0EB_->setAxisTitle("S4S9 of the 1st Pi0 Photon ", 1);
562 
563  hS4S91Pi0EE_ = ibooker.book1D("S4S91Pi0EE", "S4S9 1st most energetic Pi0 photon in EE", 50, 0., 1.);
564  hS4S91Pi0EE_->setAxisTitle("S4S9 of the 1st Pi0 Photon ", 1);
565 
566  hS4S91EtaEB_ = ibooker.book1D("S4S91EtaEB", "S4S9 1st most energetic Eta photon in EB", 50, 0., 1.);
567  hS4S91EtaEB_->setAxisTitle("S4S9 of the 1st Eta Photon ", 1);
568 
569  hS4S91EtaEE_ = ibooker.book1D("S4S91EtaEE", "S4S9 1st most energetic Eta photon in EE", 50, 0., 1.);
570  hS4S91EtaEE_->setAxisTitle("S4S9 of the 1st Eta Photon ", 1);
571 
572  hS4S92Pi0EB_ = ibooker.book1D("S4S92Pi0EB", "S4S9 2nd most energetic Pi0 photon in EB", 50, 0., 1.);
573  hS4S92Pi0EB_->setAxisTitle("S4S9 of the 2nd Pi0 Photon", 1);
574 
575  hS4S92Pi0EE_ = ibooker.book1D("S4S92Pi0EE", "S4S9 2nd most energetic Pi0 photon in EE", 50, 0., 1.);
576  hS4S92Pi0EE_->setAxisTitle("S4S9 of the 2nd Pi0 Photon", 1);
577 
578  hS4S92EtaEB_ = ibooker.book1D("S4S92EtaEB", "S4S9 2nd most energetic Pi0 photon in EB", 50, 0., 1.);
579  hS4S92EtaEB_->setAxisTitle("S4S9 of the 2nd Eta Photon", 1);
580 
581  hS4S92EtaEE_ = ibooker.book1D("S4S92EtaEE", "S4S9 2nd most energetic Pi0 photon in EE", 50, 0., 1.);
582  hS4S92EtaEE_->setAxisTitle("S4S9 of the 2nd Eta Photon", 1);
583 }
584 
585 //-------------------------------------------------------------
586 void DQMSourcePi0::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
587  if (eventCounter_ % prescaleFactor_)
588  return;
589  eventCounter_++;
590 
591  auto const &theCaloTopology = iSetup.getHandle(caloTopoToken_);
592 
593  std::vector<EcalRecHit> seeds;
594  seeds.clear();
595 
596  std::vector<EBDetId> usedXtals;
597  usedXtals.clear();
598 
599  detIdEBRecHits.clear();
600  EBRecHits.clear();
601 
602  edm::Handle<EcalRecHitCollection> rhEBpi0;
603  edm::Handle<EcalRecHitCollection> rhEBeta;
604  edm::Handle<EcalRecHitCollection> rhEEpi0;
605  edm::Handle<EcalRecHitCollection> rhEEeta;
606 
607  if (isMonEBpi0_)
608  iEvent.getByToken(productMonitoredEBpi0_, rhEBpi0);
609  if (isMonEBeta_)
610  iEvent.getByToken(productMonitoredEBeta_, rhEBeta);
611  if (isMonEEpi0_)
612  iEvent.getByToken(productMonitoredEEpi0_, rhEEpi0);
613  if (isMonEEeta_)
614  iEvent.getByToken(productMonitoredEEeta_, rhEEeta);
615 
616  // Initialize the Position Calc
617 
618  //edm::ESHandle<CaloGeometry> geoHandle;
619  //iSetup.get<CaloGeometryRecord>().get(geoHandle);
620  const auto &geoHandle = iSetup.getHandle(caloGeomToken_);
621  const CaloSubdetectorGeometry *geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
622  const CaloSubdetectorGeometry *geometryEE_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
623  const CaloSubdetectorGeometry *geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
624 
625  const CaloSubdetectorTopology *topology_p = theCaloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel);
626  const CaloSubdetectorTopology *topology_ee = theCaloTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap);
627 
628  EcalRecHitCollection::const_iterator itb;
629 
630  // fill EB pi0 histos
631  if (isMonEBpi0_) {
632  if (rhEBpi0.isValid() && (!rhEBpi0->empty())) {
633  const EcalRecHitCollection *hitCollection_p = rhEBpi0.product();
634  float etot = 0;
635  for (itb = rhEBpi0->begin(); itb != rhEBpi0->end(); ++itb) {
636  EBDetId id(itb->id());
637  double energy = itb->energy();
638  if (energy < seleXtalMinEnergy_)
639  continue;
640 
641  EBDetId det = itb->id();
642 
643  detIdEBRecHits.push_back(det);
644  EBRecHits.push_back(*itb);
645 
646  if (energy > clusSeedThr_)
647  seeds.push_back(*itb);
648 
649  hiPhiDistrEBpi0_->Fill(id.iphi());
650  hiEtaDistrEBpi0_->Fill(id.ieta());
651  hRechitEnergyEBpi0_->Fill(itb->energy());
652 
653  etot += itb->energy();
654  } // Eb rechits
655 
656  hNRecHitsEBpi0_->Fill(rhEBpi0->size());
657  hMeanRecHitEnergyEBpi0_->Fill(etot / rhEBpi0->size());
658  hEventEnergyEBpi0_->Fill(etot);
659 
660  // cout << " EB RH Pi0 collection: #, mean rh_e, event E
661  // "<<rhEBpi0->size()<<" "<<etot/rhEBpi0->size()<<" "<<etot<<endl;
662 
663  // Pi0 maker
664 
665  // cout<< " RH coll size: "<<rhEBpi0->size()<<endl;
666  // cout<< " Pi0 seeds: "<<seeds.size()<<endl;
667 
668  int nClus;
669  std::vector<float> eClus;
670  std::vector<float> etClus;
671  std::vector<float> etaClus;
672  std::vector<float> thetaClus;
673  std::vector<float> phiClus;
674  std::vector<EBDetId> max_hit;
675 
676  std::vector<std::vector<EcalRecHit>> RecHitsCluster;
677  std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5;
678  std::vector<float> s4s9Clus;
679  std::vector<float> s9s25Clus;
680 
681  nClus = 0;
682 
683  // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
684  std::sort(seeds.begin(), seeds.end(), ecalRecHitGreater);
685 
686  for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
687  EBDetId seed_id = itseed->id();
688  std::vector<EBDetId>::const_iterator usedIds;
689 
690  bool seedAlreadyUsed = false;
691  for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
692  if (*usedIds == seed_id) {
693  seedAlreadyUsed = true;
694  // cout<< " Seed with energy "<<itseed->energy()<<" was used
695  // !"<<endl;
696  break;
697  }
698  }
699  if (seedAlreadyUsed)
700  continue;
701  std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
702  std::vector<std::pair<DetId, float>> clus_used;
703  // Reject the seed if not able to build the cluster around it correctly
704  // if(clus_v.size() < clusEtaSize_*clusPhiSize_){cout<<" Not enough
705  // RecHits "<<endl; continue;}
706  std::vector<EcalRecHit> RecHitsInWindow;
707  std::vector<EcalRecHit> RecHitsInWindow5x5;
708 
709  double simple_energy = 0;
710 
711  for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
712  EBDetId EBdet = *det;
713  // cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi
714  // "<<EBdet.iphi()<<endl;
715  bool HitAlreadyUsed = false;
716  for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
717  if (*usedIds == *det) {
718  HitAlreadyUsed = true;
719  break;
720  }
721  }
722  if (HitAlreadyUsed)
723  continue;
724 
725  std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), EBdet);
726  if (itdet == detIdEBRecHits.end())
727  continue;
728 
729  int nn = int(itdet - detIdEBRecHits.begin());
730  usedXtals.push_back(*det);
731  RecHitsInWindow.push_back(EBRecHits[nn]);
732  clus_used.push_back(std::make_pair(*det, 1));
733  simple_energy = simple_energy + EBRecHits[nn].energy();
734  }
735 
736  if (simple_energy <= 0)
737  continue;
738 
739  math::XYZPoint clus_pos =
740  posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
741  // cout<< " Simple Clustering: Total energy for this simple
742  // cluster : "<<simple_energy<<endl; cout<< " Simple Clustering:
743  // eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<endl; cout<< "
744  // Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<"
745  // "<<clus_pos.z()<<endl;
746 
747  float theta_s = 2. * atan(exp(-clus_pos.eta()));
748  // float p0x_s = simple_energy * sin(theta_s) *
749  // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
750  // sin(clus_pos.phi());
751  // float p0z_s = simple_energy * cos(theta_s);
752  // float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
753 
754  float et_s = simple_energy * sin(theta_s);
755  // cout << " Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<"
756  // "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
757 
758  // Compute S4/S9 variable
759  // We are not sure to have 9 RecHits so need to check eta and phi:
761  float s4s9_tmp[4];
762  for (int i = 0; i < 4; i++)
763  s4s9_tmp[i] = 0;
764 
765  int seed_ieta = seed_id.ieta();
766  int seed_iphi = seed_id.iphi();
767 
768  convxtalid(seed_iphi, seed_ieta);
769 
770  float e3x3 = 0;
771  float e5x5 = 0;
772 
773  for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
774  EBDetId det = (EBDetId)RecHitsInWindow[j].id();
775 
776  int ieta = det.ieta();
777  int iphi = det.iphi();
778 
779  convxtalid(iphi, ieta);
780 
781  float en = RecHitsInWindow[j].energy();
782 
783  int dx = diff_neta_s(seed_ieta, ieta);
784  int dy = diff_nphi_s(seed_iphi, iphi);
785 
786  if (dx <= 0 && dy <= 0)
787  s4s9_tmp[0] += en;
788  if (dx >= 0 && dy <= 0)
789  s4s9_tmp[1] += en;
790  if (dx <= 0 && dy >= 0)
791  s4s9_tmp[2] += en;
792  if (dx >= 0 && dy >= 0)
793  s4s9_tmp[3] += en;
794 
795  if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
796  e3x3 += en;
797  if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
798  e5x5 += en;
799  }
800 
801  if (e3x3 <= 0)
802  continue;
803 
804  float s4s9_max = *std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3;
805 
807  std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id, 5, 5);
808  for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
809  EBDetId det = *idItr;
810 
811  std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(), usedXtals.end(), det);
812 
814  if (itdet0 != usedXtals.end())
815  continue;
816 
817  // inside collections
818  std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), det);
819  if (itdet == detIdEBRecHits.end())
820  continue;
821 
822  int nn = int(itdet - detIdEBRecHits.begin());
823 
824  RecHitsInWindow5x5.push_back(EBRecHits[nn]);
825  e5x5 += EBRecHits[nn].energy();
826  }
827 
828  if (e5x5 <= 0)
829  continue;
830 
831  eClus.push_back(simple_energy);
832  etClus.push_back(et_s);
833  etaClus.push_back(clus_pos.eta());
834  thetaClus.push_back(theta_s);
835  phiClus.push_back(clus_pos.phi());
836  s4s9Clus.push_back(s4s9_max);
837  s9s25Clus.push_back(e3x3 / e5x5);
838  RecHitsCluster.push_back(RecHitsInWindow);
839  RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
840 
841  // std::cout<<" EB pi0 cluster (n,nxt,e,et eta,phi,s4s9)
842  //"<<nClus<<" "<<int(RecHitsInWindow.size())<<" "<<eClus[nClus]<<" "<<"
843  //"<<etClus[nClus]<<" "<<etaClus[nClus]<<" "<<phiClus[nClus]<<"
844  //"<<s4s9Clus[nClus]<<std::endl;
845 
846  nClus++;
847  }
848 
849  // cout<< " Pi0 clusters: "<<nClus<<endl;
850 
851  // Selection, based on Simple clustering
852  // pi0 candidates
853 
854  // if (nClus <= 1) return;
855  for (Int_t i = 0; i < nClus; i++) {
856  for (Int_t j = i + 1; j < nClus; j++) {
857  // cout<<" i "<<i<<" etClus[i] "<<etClus[i]<<" j "<<j<<"
858  // etClus[j] "<<etClus[j]<<endl;
859  if (etClus[i] > selePtGamma_ && etClus[j] > selePtGamma_ && s4s9Clus[i] > seleS4S9Gamma_ &&
860  s4s9Clus[j] > seleS4S9Gamma_) {
861  float p0x = etClus[i] * cos(phiClus[i]);
862  float p1x = etClus[j] * cos(phiClus[j]);
863  float p0y = etClus[i] * sin(phiClus[i]);
864  float p1y = etClus[j] * sin(phiClus[j]);
865  float p0z = eClus[i] * cos(thetaClus[i]);
866  float p1z = eClus[j] * cos(thetaClus[j]);
867 
868  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
869 
870  if (pt_pair < selePtPi0_)
871  continue;
872 
873  float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
874  (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
875  if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
876  // New Loop on cluster to measure isolation:
877  std::vector<int> IsoClus;
878  IsoClus.clear();
879  float Iso = 0;
880  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
881  for (Int_t k = 0; k < nClus; k++) {
882  if (etClus[k] < ptMinForIsolation_)
883  continue;
884 
885  if (k == i || k == j)
886  continue;
887  TVector3 ClusVect =
888  TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
889 
890  float dretacl = fabs(etaClus[k] - pairVect.Eta());
891  float drcl = ClusVect.DeltaR(pairVect);
892  // cout<< " Iso: k, E, drclpi0, detaclpi0, dphiclpi0
893  // "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
894  // "<<dretaclpi0<<endl;
895  if ((drcl < selePi0BeltDR_) && (dretacl < selePi0BeltDeta_)) {
896  // cout<< " ... good iso cluster #: "<<k<<"
897  // etClus[k] "<<etClus[k] <<endl;
898  Iso = Iso + etClus[k];
899  IsoClus.push_back(k);
900  }
901  }
902 
903  // cout<<" Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
904  if (Iso / pt_pair < selePi0Iso_) {
905  // for(unsigned int
906  // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
907  // for(unsigned int
908  // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
909 
910  hMinvPi0EB_->Fill(m_inv);
911  hPt1Pi0EB_->Fill(etClus[i]);
912  hPt2Pi0EB_->Fill(etClus[j]);
913  hPtPi0EB_->Fill(pt_pair);
914  hIsoPi0EB_->Fill(Iso / pt_pair);
915  hS4S91Pi0EB_->Fill(s4s9Clus[i]);
916  hS4S92Pi0EB_->Fill(s4s9Clus[j]);
917 
918  // cout <<" EB Simple Clustering: pi0 Candidate
919  // pt, eta, phi, Iso, m_inv, i, j : "<<pt_pair<<"
920  //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
921  //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
922  }
923  }
924  }
925  } // End of the "j" loop over Simple Clusters
926  } // End of the "i" loop over Simple Clusters
927 
928  } // rhEBpi0.valid() ends
929 
930  } // isMonEBpi0 ends
931 
932  //------------------ End of pi0 in EB --------------------------//
933 
934  // fill EB eta histos
935  if (isMonEBeta_) {
936  if (rhEBeta.isValid() && (!rhEBeta->empty())) {
937  const EcalRecHitCollection *hitCollection_p = rhEBeta.product();
938  float etot = 0;
939  for (itb = rhEBeta->begin(); itb != rhEBeta->end(); ++itb) {
940  EBDetId id(itb->id());
941  double energy = itb->energy();
942  if (energy < seleXtalMinEnergy_)
943  continue;
944 
945  EBDetId det = itb->id();
946 
947  detIdEBRecHits.push_back(det);
948  EBRecHits.push_back(*itb);
949 
950  if (energy > clusSeedThr_)
951  seeds.push_back(*itb);
952 
953  hiPhiDistrEBeta_->Fill(id.iphi());
954  hiEtaDistrEBeta_->Fill(id.ieta());
955  hRechitEnergyEBeta_->Fill(itb->energy());
956 
957  etot += itb->energy();
958  } // Eb rechits
959 
960  hNRecHitsEBeta_->Fill(rhEBeta->size());
961  hMeanRecHitEnergyEBeta_->Fill(etot / rhEBeta->size());
962  hEventEnergyEBeta_->Fill(etot);
963 
964  // cout << " EB RH Eta collection: #, mean rh_e, event E
965  // "<<rhEBeta->size()<<" "<<etot/rhEBeta->size()<<" "<<etot<<endl;
966 
967  // Eta maker
968 
969  // cout<< " RH coll size: "<<rhEBeta->size()<<endl;
970  // cout<< " Eta seeds: "<<seeds.size()<<endl;
971 
972  int nClus;
973  std::vector<float> eClus;
974  std::vector<float> etClus;
975  std::vector<float> etaClus;
976  std::vector<float> thetaClus;
977  std::vector<float> phiClus;
978  std::vector<EBDetId> max_hit;
979 
980  std::vector<std::vector<EcalRecHit>> RecHitsCluster;
981  std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5;
982  std::vector<float> s4s9Clus;
983  std::vector<float> s9s25Clus;
984 
985  nClus = 0;
986 
987  // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
988  std::sort(seeds.begin(), seeds.end(), ecalRecHitGreater);
989 
990  for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
991  EBDetId seed_id = itseed->id();
992  std::vector<EBDetId>::const_iterator usedIds;
993 
994  bool seedAlreadyUsed = false;
995  for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
996  if (*usedIds == seed_id) {
997  seedAlreadyUsed = true;
998  // cout<< " Seed with energy "<<itseed->energy()<<" was used
999  // !"<<endl;
1000  break;
1001  }
1002  }
1003  if (seedAlreadyUsed)
1004  continue;
1005  std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1006  std::vector<std::pair<DetId, float>> clus_used;
1007  // Reject the seed if not able to build the cluster around it correctly
1008  // if(clus_v.size() < clusEtaSize_*clusPhiSize_){cout<<" Not enough
1009  // RecHits "<<endl; continue;}
1010  std::vector<EcalRecHit> RecHitsInWindow;
1011  std::vector<EcalRecHit> RecHitsInWindow5x5;
1012 
1013  double simple_energy = 0;
1014 
1015  for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1016  EBDetId EBdet = *det;
1017  // cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi
1018  // "<<EBdet.iphi()<<endl;
1019  bool HitAlreadyUsed = false;
1020  for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
1021  if (*usedIds == *det) {
1022  HitAlreadyUsed = true;
1023  break;
1024  }
1025  }
1026  if (HitAlreadyUsed)
1027  continue;
1028 
1029  std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), EBdet);
1030  if (itdet == detIdEBRecHits.end())
1031  continue;
1032 
1033  int nn = int(itdet - detIdEBRecHits.begin());
1034  usedXtals.push_back(*det);
1035  RecHitsInWindow.push_back(EBRecHits[nn]);
1036  clus_used.push_back(std::make_pair(*det, 1));
1037  simple_energy = simple_energy + EBRecHits[nn].energy();
1038  }
1039 
1040  if (simple_energy <= 0)
1041  continue;
1042 
1043  math::XYZPoint clus_pos =
1044  posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
1045  // cout<< " Simple Clustering: Total energy for this simple
1046  // cluster : "<<simple_energy<<endl; cout<< " Simple Clustering:
1047  // eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<endl; cout<< "
1048  // Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<"
1049  // "<<clus_pos.z()<<endl;
1050 
1051  float theta_s = 2. * atan(exp(-clus_pos.eta()));
1052  // float p0x_s = simple_energy * sin(theta_s) *
1053  // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1054  // sin(clus_pos.phi());
1055  // float p0z_s = simple_energy * cos(theta_s);
1056  // float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1057 
1058  float et_s = simple_energy * sin(theta_s);
1059  // cout << " Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<"
1060  // "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
1061 
1062  // Compute S4/S9 variable
1063  // We are not sure to have 9 RecHits so need to check eta and phi:
1065  float s4s9_tmp[4];
1066  for (int i = 0; i < 4; i++)
1067  s4s9_tmp[i] = 0;
1068 
1069  int seed_ieta = seed_id.ieta();
1070  int seed_iphi = seed_id.iphi();
1071 
1072  convxtalid(seed_iphi, seed_ieta);
1073 
1074  float e3x3 = 0;
1075  float e5x5 = 0;
1076 
1077  for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1078  EBDetId det = (EBDetId)RecHitsInWindow[j].id();
1079 
1080  int ieta = det.ieta();
1081  int iphi = det.iphi();
1082 
1083  convxtalid(iphi, ieta);
1084 
1085  float en = RecHitsInWindow[j].energy();
1086 
1087  int dx = diff_neta_s(seed_ieta, ieta);
1088  int dy = diff_nphi_s(seed_iphi, iphi);
1089 
1090  if (dx <= 0 && dy <= 0)
1091  s4s9_tmp[0] += en;
1092  if (dx >= 0 && dy <= 0)
1093  s4s9_tmp[1] += en;
1094  if (dx <= 0 && dy >= 0)
1095  s4s9_tmp[2] += en;
1096  if (dx >= 0 && dy >= 0)
1097  s4s9_tmp[3] += en;
1098 
1099  if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1100  e3x3 += en;
1101  if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1102  e5x5 += en;
1103  }
1104 
1105  if (e3x3 <= 0)
1106  continue;
1107 
1108  float s4s9_max = *std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3;
1109 
1111  std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id, 5, 5);
1112  for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
1113  EBDetId det = *idItr;
1114 
1115  std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(), usedXtals.end(), det);
1116 
1118  if (itdet0 != usedXtals.end())
1119  continue;
1120 
1121  // inside collections
1122  std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), det);
1123  if (itdet == detIdEBRecHits.end())
1124  continue;
1125 
1126  int nn = int(itdet - detIdEBRecHits.begin());
1127 
1128  RecHitsInWindow5x5.push_back(EBRecHits[nn]);
1129  e5x5 += EBRecHits[nn].energy();
1130  }
1131 
1132  if (e5x5 <= 0)
1133  continue;
1134 
1135  eClus.push_back(simple_energy);
1136  etClus.push_back(et_s);
1137  etaClus.push_back(clus_pos.eta());
1138  thetaClus.push_back(theta_s);
1139  phiClus.push_back(clus_pos.phi());
1140  s4s9Clus.push_back(s4s9_max);
1141  s9s25Clus.push_back(e3x3 / e5x5);
1142  RecHitsCluster.push_back(RecHitsInWindow);
1143  RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
1144 
1145  // std::cout<<" EB Eta cluster (n,nxt,e,et eta,phi,s4s9)
1146  //"<<nClus<<" "<<int(RecHitsInWindow.size())<<" "<<eClus[nClus]<<" "<<"
1147  //"<<etClus[nClus]<<" "<<etaClus[nClus]<<" "<<phiClus[nClus]<<"
1148  //"<<s4s9Clus[nClus]<<std::endl;
1149 
1150  nClus++;
1151  }
1152 
1153  // cout<< " Eta clusters: "<<nClus<<endl;
1154 
1155  // Selection, based on Simple clustering
1156  // eta candidates
1157 
1158  // if (nClus <= 1) return;
1159  for (Int_t i = 0; i < nClus; i++) {
1160  for (Int_t j = i + 1; j < nClus; j++) {
1161  // cout<<" i "<<i<<" etClus[i] "<<etClus[i]<<" j "<<j<<"
1162  // etClus[j] "<<etClus[j]<<endl;
1163  if (etClus[i] > selePtGammaEta_ && etClus[j] > selePtGammaEta_ && s4s9Clus[i] > seleS4S9GammaEta_ &&
1164  s4s9Clus[j] > seleS4S9GammaEta_) {
1165  float p0x = etClus[i] * cos(phiClus[i]);
1166  float p1x = etClus[j] * cos(phiClus[j]);
1167  float p0y = etClus[i] * sin(phiClus[i]);
1168  float p1y = etClus[j] * sin(phiClus[j]);
1169  float p0z = eClus[i] * cos(thetaClus[i]);
1170  float p1z = eClus[j] * cos(thetaClus[j]);
1171 
1172  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1173 
1174  if (pt_pair < selePtEta_)
1175  continue;
1176 
1177  float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
1178  (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1179  if ((m_inv < seleMinvMaxEta_) && (m_inv > seleMinvMinEta_)) {
1180  // New Loop on cluster to measure isolation:
1181  std::vector<int> IsoClus;
1182  IsoClus.clear();
1183  float Iso = 0;
1184  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1185  for (Int_t k = 0; k < nClus; k++) {
1186  if (etClus[k] < ptMinForIsolationEta_)
1187  continue;
1188 
1189  if (k == i || k == j)
1190  continue;
1191  TVector3 ClusVect =
1192  TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
1193 
1194  float dretacl = fabs(etaClus[k] - pairVect.Eta());
1195  float drcl = ClusVect.DeltaR(pairVect);
1196  // cout<< " Iso: k, E, drclpi0, detaclpi0, dphiclpi0
1197  // "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
1198  // "<<dretaclpi0<<endl;
1199  if ((drcl < seleEtaBeltDR_) && (dretacl < seleEtaBeltDeta_)) {
1200  // cout<< " ... good iso cluster #: "<<k<<"
1201  // etClus[k] "<<etClus[k] <<endl;
1202  Iso = Iso + etClus[k];
1203  IsoClus.push_back(k);
1204  }
1205  }
1206 
1207  // cout<<" Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
1208  if (Iso / pt_pair < seleEtaIso_) {
1209  // for(unsigned int
1210  // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
1211  // for(unsigned int
1212  // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
1213 
1214  hMinvEtaEB_->Fill(m_inv);
1215  hPt1EtaEB_->Fill(etClus[i]);
1216  hPt2EtaEB_->Fill(etClus[j]);
1217  hPtEtaEB_->Fill(pt_pair);
1218  hIsoEtaEB_->Fill(Iso / pt_pair);
1219  hS4S91EtaEB_->Fill(s4s9Clus[i]);
1220  hS4S92EtaEB_->Fill(s4s9Clus[j]);
1221 
1222  // cout <<" EB Simple Clustering: Eta Candidate
1223  // pt, eta, phi, Iso, m_inv, i, j : "<<pt_pair<<"
1224  //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
1225  //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
1226  }
1227  }
1228  }
1229  } // End of the "j" loop over Simple Clusters
1230  } // End of the "i" loop over Simple Clusters
1231 
1232  } // rhEBeta.valid() ends
1233 
1234  } // isMonEBeta ends
1235 
1236  //------------------ End of Eta in EB --------------------------//
1237 
1238  //----------------- End of the EB --------------------------//
1239 
1240  //----------------- Start of the EE --------------------//
1241 
1242  // fill pi0 EE histos
1243  if (isMonEEpi0_) {
1244  if (rhEEpi0.isValid() && (!rhEEpi0->empty())) {
1245  const EcalRecHitCollection *hitCollection_ee = rhEEpi0.product();
1246  float etot = 0;
1247 
1248  detIdEERecHits.clear();
1249  EERecHits.clear();
1250 
1251  std::vector<EcalRecHit> seedsEndCap;
1252  seedsEndCap.clear();
1253 
1254  std::vector<EEDetId> usedXtalsEndCap;
1255  usedXtalsEndCap.clear();
1256 
1258  EERecHitCollection::const_iterator ite;
1259  for (ite = rhEEpi0->begin(); ite != rhEEpi0->end(); ite++) {
1260  double energy = ite->energy();
1261  if (energy < seleXtalMinEnergyEndCap_)
1262  continue;
1263 
1264  EEDetId det = ite->id();
1265  EEDetId id(ite->id());
1266 
1267  detIdEERecHits.push_back(det);
1268  EERecHits.push_back(*ite);
1269 
1270  hiXDistrEEpi0_->Fill(id.ix());
1271  hiYDistrEEpi0_->Fill(id.iy());
1272  hRechitEnergyEEpi0_->Fill(ite->energy());
1273 
1274  if (energy > clusSeedThrEndCap_)
1275  seedsEndCap.push_back(*ite);
1276 
1277  etot += ite->energy();
1278  } // EE rechits
1279 
1280  hNRecHitsEEpi0_->Fill(rhEEpi0->size());
1281  hMeanRecHitEnergyEEpi0_->Fill(etot / rhEEpi0->size());
1282  hEventEnergyEEpi0_->Fill(etot);
1283 
1284  // cout << " EE RH Pi0 collection: #, mean rh_e, event E
1285  // "<<rhEEpi0->size()<<" "<<etot/rhEEpi0->size()<<" "<<etot<<endl;
1286 
1287  int nClusEndCap;
1288  std::vector<float> eClusEndCap;
1289  std::vector<float> etClusEndCap;
1290  std::vector<float> etaClusEndCap;
1291  std::vector<float> thetaClusEndCap;
1292  std::vector<float> phiClusEndCap;
1293  std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1294  std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1295  std::vector<float> s4s9ClusEndCap;
1296  std::vector<float> s9s25ClusEndCap;
1297 
1298  nClusEndCap = 0;
1299 
1300  // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1301  std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1302 
1303  for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1304  EEDetId seed_id = itseed->id();
1305  std::vector<EEDetId>::const_iterator usedIds;
1306 
1307  bool seedAlreadyUsed = false;
1308  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1309  if (*usedIds == seed_id) {
1310  seedAlreadyUsed = true;
1311  break;
1312  }
1313  }
1314 
1315  if (seedAlreadyUsed)
1316  continue;
1317  std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1318  std::vector<std::pair<DetId, float>> clus_used;
1319 
1320  std::vector<EcalRecHit> RecHitsInWindow;
1321  std::vector<EcalRecHit> RecHitsInWindow5x5;
1322 
1323  float simple_energy = 0;
1324 
1325  for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1326  EEDetId EEdet = *det;
1327 
1328  bool HitAlreadyUsed = false;
1329  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1330  if (*usedIds == *det) {
1331  HitAlreadyUsed = true;
1332  break;
1333  }
1334  }
1335 
1336  if (HitAlreadyUsed)
1337  continue;
1338 
1339  std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1340  if (itdet == detIdEERecHits.end())
1341  continue;
1342 
1343  int nn = int(itdet - detIdEERecHits.begin());
1344  usedXtalsEndCap.push_back(*det);
1345  RecHitsInWindow.push_back(EERecHits[nn]);
1346  clus_used.push_back(std::make_pair(*det, 1));
1347  simple_energy = simple_energy + EERecHits[nn].energy();
1348  }
1349 
1350  if (simple_energy <= 0)
1351  continue;
1352 
1353  math::XYZPoint clus_pos =
1354  posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1355 
1356  float theta_s = 2. * atan(exp(-clus_pos.eta()));
1357  float et_s = simple_energy * sin(theta_s);
1358  // float p0x_s = simple_energy * sin(theta_s) *
1359  // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1360  // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1361 
1362  // Compute S4/S9 variable
1363  // We are not sure to have 9 RecHits so need to check eta and phi:
1364  float s4s9_tmp[4];
1365  for (int i = 0; i < 4; i++)
1366  s4s9_tmp[i] = 0;
1367 
1368  int ixSeed = seed_id.ix();
1369  int iySeed = seed_id.iy();
1370  float e3x3 = 0;
1371  float e5x5 = 0;
1372 
1373  for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1374  EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1375  int dx = ixSeed - det_this.ix();
1376  int dy = iySeed - det_this.iy();
1377 
1378  float en = RecHitsInWindow[j].energy();
1379 
1380  if (dx <= 0 && dy <= 0)
1381  s4s9_tmp[0] += en;
1382  if (dx >= 0 && dy <= 0)
1383  s4s9_tmp[1] += en;
1384  if (dx <= 0 && dy >= 0)
1385  s4s9_tmp[2] += en;
1386  if (dx >= 0 && dy >= 0)
1387  s4s9_tmp[3] += en;
1388 
1389  if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1390  e3x3 += en;
1391  if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1392  e5x5 += en;
1393  }
1394 
1395  if (e3x3 <= 0)
1396  continue;
1397 
1398  eClusEndCap.push_back(simple_energy);
1399  etClusEndCap.push_back(et_s);
1400  etaClusEndCap.push_back(clus_pos.eta());
1401  thetaClusEndCap.push_back(theta_s);
1402  phiClusEndCap.push_back(clus_pos.phi());
1403  s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1404  s9s25ClusEndCap.push_back(e3x3 / e5x5);
1405  RecHitsClusterEndCap.push_back(RecHitsInWindow);
1406  RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1407 
1408  // std::cout<<" EE pi0 cluster (n,nxt,e,et eta,phi,s4s9)
1409  //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1410  //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1411  //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1412  //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1413 
1414  nClusEndCap++;
1415  }
1416 
1417  // Selection, based on Simple clustering
1418  // pi0 candidates
1419 
1420  for (Int_t i = 0; i < nClusEndCap; i++) {
1421  for (Int_t j = i + 1; j < nClusEndCap; j++) {
1422  if (etClusEndCap[i] > selePtGammaEndCap_ && etClusEndCap[j] > selePtGammaEndCap_ &&
1423  s4s9ClusEndCap[i] > seleS4S9GammaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEndCap_) {
1424  float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1425  float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1426  float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1427  float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1428  float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1429  float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1430 
1431  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1432  if (pt_pair < selePtPi0EndCap_)
1433  continue;
1434  float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1435  (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1436 
1437  if ((m_inv < seleMinvMaxPi0EndCap_) && (m_inv > seleMinvMinPi0EndCap_)) {
1438  // New Loop on cluster to measure isolation:
1439  std::vector<int> IsoClus;
1440  IsoClus.clear();
1441  float Iso = 0;
1442  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1443  for (Int_t k = 0; k < nClusEndCap; k++) {
1444  if (etClusEndCap[k] < ptMinForIsolationEndCap_)
1445  continue;
1446 
1447  if (k == i || k == j)
1448  continue;
1449 
1450  TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1451  etClusEndCap[k] * sin(phiClusEndCap[k]),
1452  eClusEndCap[k] * cos(thetaClusEndCap[k]));
1453  float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1454  float drcl = clusVect.DeltaR(pairVect);
1455 
1456  if (drcl < selePi0BeltDREndCap_ && dretacl < selePi0BeltDetaEndCap_) {
1457  Iso = Iso + etClusEndCap[k];
1458  IsoClus.push_back(k);
1459  }
1460  }
1461 
1462  if (Iso / pt_pair < selePi0IsoEndCap_) {
1463  // cout <<" EE Simple Clustering: pi0 Candidate pt, eta, phi,
1464  // Iso, m_inv, i, j : "<<pt_pair<<" "<<pairVect.Eta()<<"
1465  // "<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1466  // "<<endl;
1467 
1468  hMinvPi0EE_->Fill(m_inv);
1469  hPt1Pi0EE_->Fill(etClusEndCap[i]);
1470  hPt2Pi0EE_->Fill(etClusEndCap[j]);
1471  hPtPi0EE_->Fill(pt_pair);
1472  hIsoPi0EE_->Fill(Iso / pt_pair);
1473  hS4S91Pi0EE_->Fill(s4s9ClusEndCap[i]);
1474  hS4S92Pi0EE_->Fill(s4s9ClusEndCap[j]);
1475  }
1476  }
1477  }
1478  } // End of the "j" loop over Simple Clusters
1479  } // End of the "i" loop over Simple Clusters
1480 
1481  } // rhEEpi0
1482  } // isMonEEpi0
1483 
1484  //================End of Pi0 endcap=======================//
1485 
1486  //================ Eta in EE===============================//
1487 
1488  // fill pi0 EE histos
1489  if (isMonEEeta_) {
1490  if (rhEEeta.isValid() && (!rhEEeta->empty())) {
1491  const EcalRecHitCollection *hitCollection_ee = rhEEeta.product();
1492  float etot = 0;
1493 
1494  detIdEERecHits.clear();
1495  EERecHits.clear();
1496 
1497  std::vector<EcalRecHit> seedsEndCap;
1498  seedsEndCap.clear();
1499 
1500  std::vector<EEDetId> usedXtalsEndCap;
1501  usedXtalsEndCap.clear();
1502 
1504  EERecHitCollection::const_iterator ite;
1505  for (ite = rhEEeta->begin(); ite != rhEEeta->end(); ite++) {
1506  double energy = ite->energy();
1507  if (energy < seleXtalMinEnergyEndCap_)
1508  continue;
1509 
1510  EEDetId det = ite->id();
1511  EEDetId id(ite->id());
1512 
1513  detIdEERecHits.push_back(det);
1514  EERecHits.push_back(*ite);
1515 
1516  hiXDistrEEeta_->Fill(id.ix());
1517  hiYDistrEEeta_->Fill(id.iy());
1518  hRechitEnergyEEeta_->Fill(ite->energy());
1519 
1520  if (energy > clusSeedThrEndCap_)
1521  seedsEndCap.push_back(*ite);
1522 
1523  etot += ite->energy();
1524  } // EE rechits
1525 
1526  hNRecHitsEEeta_->Fill(rhEEeta->size());
1527  hMeanRecHitEnergyEEeta_->Fill(etot / rhEEeta->size());
1528  hEventEnergyEEeta_->Fill(etot);
1529 
1530  // cout << " EE RH Eta collection: #, mean rh_e, event E
1531  // "<<rhEEeta->size()<<" "<<etot/rhEEeta->size()<<" "<<etot<<endl;
1532 
1533  int nClusEndCap;
1534  std::vector<float> eClusEndCap;
1535  std::vector<float> etClusEndCap;
1536  std::vector<float> etaClusEndCap;
1537  std::vector<float> thetaClusEndCap;
1538  std::vector<float> phiClusEndCap;
1539  std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1540  std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1541  std::vector<float> s4s9ClusEndCap;
1542  std::vector<float> s9s25ClusEndCap;
1543 
1544  nClusEndCap = 0;
1545 
1546  // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1547  std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1548 
1549  for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1550  EEDetId seed_id = itseed->id();
1551  std::vector<EEDetId>::const_iterator usedIds;
1552 
1553  bool seedAlreadyUsed = false;
1554  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1555  if (*usedIds == seed_id) {
1556  seedAlreadyUsed = true;
1557  break;
1558  }
1559  }
1560 
1561  if (seedAlreadyUsed)
1562  continue;
1563  std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1564  std::vector<std::pair<DetId, float>> clus_used;
1565 
1566  std::vector<EcalRecHit> RecHitsInWindow;
1567  std::vector<EcalRecHit> RecHitsInWindow5x5;
1568 
1569  float simple_energy = 0;
1570 
1571  for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1572  EEDetId EEdet = *det;
1573 
1574  bool HitAlreadyUsed = false;
1575  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1576  if (*usedIds == *det) {
1577  HitAlreadyUsed = true;
1578  break;
1579  }
1580  }
1581 
1582  if (HitAlreadyUsed)
1583  continue;
1584 
1585  std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1586  if (itdet == detIdEERecHits.end())
1587  continue;
1588 
1589  int nn = int(itdet - detIdEERecHits.begin());
1590  usedXtalsEndCap.push_back(*det);
1591  RecHitsInWindow.push_back(EERecHits[nn]);
1592  clus_used.push_back(std::make_pair(*det, 1));
1593  simple_energy = simple_energy + EERecHits[nn].energy();
1594  }
1595 
1596  if (simple_energy <= 0)
1597  continue;
1598 
1599  math::XYZPoint clus_pos =
1600  posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1601 
1602  float theta_s = 2. * atan(exp(-clus_pos.eta()));
1603  float et_s = simple_energy * sin(theta_s);
1604  // float p0x_s = simple_energy * sin(theta_s) *
1605  // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1606  // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1607 
1608  // Compute S4/S9 variable
1609  // We are not sure to have 9 RecHits so need to check eta and phi:
1610  float s4s9_tmp[4];
1611  for (int i = 0; i < 4; i++)
1612  s4s9_tmp[i] = 0;
1613 
1614  int ixSeed = seed_id.ix();
1615  int iySeed = seed_id.iy();
1616  float e3x3 = 0;
1617  float e5x5 = 0;
1618 
1619  for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1620  EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1621  int dx = ixSeed - det_this.ix();
1622  int dy = iySeed - det_this.iy();
1623 
1624  float en = RecHitsInWindow[j].energy();
1625 
1626  if (dx <= 0 && dy <= 0)
1627  s4s9_tmp[0] += en;
1628  if (dx >= 0 && dy <= 0)
1629  s4s9_tmp[1] += en;
1630  if (dx <= 0 && dy >= 0)
1631  s4s9_tmp[2] += en;
1632  if (dx >= 0 && dy >= 0)
1633  s4s9_tmp[3] += en;
1634 
1635  if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1636  e3x3 += en;
1637  if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1638  e5x5 += en;
1639  }
1640 
1641  if (e3x3 <= 0)
1642  continue;
1643 
1644  eClusEndCap.push_back(simple_energy);
1645  etClusEndCap.push_back(et_s);
1646  etaClusEndCap.push_back(clus_pos.eta());
1647  thetaClusEndCap.push_back(theta_s);
1648  phiClusEndCap.push_back(clus_pos.phi());
1649  s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1650  s9s25ClusEndCap.push_back(e3x3 / e5x5);
1651  RecHitsClusterEndCap.push_back(RecHitsInWindow);
1652  RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1653 
1654  // std::cout<<" EE Eta cluster (n,nxt,e,et eta,phi,s4s9)
1655  //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1656  //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1657  //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1658  //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1659 
1660  nClusEndCap++;
1661  }
1662 
1663  // Selection, based on Simple clustering
1664  // pi0 candidates
1665 
1666  for (Int_t i = 0; i < nClusEndCap; i++) {
1667  for (Int_t j = i + 1; j < nClusEndCap; j++) {
1668  if (etClusEndCap[i] > selePtGammaEtaEndCap_ && etClusEndCap[j] > selePtGammaEtaEndCap_ &&
1669  s4s9ClusEndCap[i] > seleS4S9GammaEtaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEtaEndCap_) {
1670  float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1671  float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1672  float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1673  float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1674  float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1675  float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1676 
1677  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1678  if (pt_pair < selePtEtaEndCap_)
1679  continue;
1680  float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1681  (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1682 
1683  if ((m_inv < seleMinvMaxEtaEndCap_) && (m_inv > seleMinvMinEtaEndCap_)) {
1684  // New Loop on cluster to measure isolation:
1685  std::vector<int> IsoClus;
1686  IsoClus.clear();
1687  float Iso = 0;
1688  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1689  for (Int_t k = 0; k < nClusEndCap; k++) {
1690  if (etClusEndCap[k] < ptMinForIsolationEtaEndCap_)
1691  continue;
1692 
1693  if (k == i || k == j)
1694  continue;
1695 
1696  TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1697  etClusEndCap[k] * sin(phiClusEndCap[k]),
1698  eClusEndCap[k] * cos(thetaClusEndCap[k]));
1699  float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1700  float drcl = clusVect.DeltaR(pairVect);
1701 
1702  if (drcl < seleEtaBeltDREndCap_ && dretacl < seleEtaBeltDetaEndCap_) {
1703  Iso = Iso + etClusEndCap[k];
1704  IsoClus.push_back(k);
1705  }
1706  }
1707 
1708  if (Iso / pt_pair < seleEtaIsoEndCap_) {
1709  // cout <<" EE Simple Clustering: Eta Candidate pt, eta,
1710  // phi, Iso, m_inv, i, j : "<<pt_pair<<" "<<pairVect.Eta()<<"
1711  //"<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1712  //"<<endl;
1713 
1714  hMinvEtaEE_->Fill(m_inv);
1715  hPt1EtaEE_->Fill(etClusEndCap[i]);
1716  hPt2EtaEE_->Fill(etClusEndCap[j]);
1717  hPtEtaEE_->Fill(pt_pair);
1718  hIsoEtaEE_->Fill(Iso / pt_pair);
1719  hS4S91EtaEE_->Fill(s4s9ClusEndCap[i]);
1720  hS4S92EtaEE_->Fill(s4s9ClusEndCap[j]);
1721  }
1722  }
1723  }
1724  } // End of the "j" loop over Simple Clusters
1725  } // End of the "i" loop over Simple Clusters
1726 
1727  } // rhEEeta
1728  } // isMonEEeta
1729 
1730  //================End of Pi0 endcap=======================//
1731 
1733 }
1734 
1735 void DQMSourcePi0::convxtalid(Int_t &nphi, Int_t &neta) {
1736  // Barrel only
1737  // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
1738  // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
1739 
1740  if (neta > 0)
1741  neta -= 1;
1742  if (nphi > 359)
1743  nphi = nphi - 360;
1744 
1745 } // end of convxtalid
1746 
1747 int DQMSourcePi0::diff_neta_s(Int_t neta1, Int_t neta2) {
1748  Int_t mdiff;
1749  mdiff = (neta1 - neta2);
1750  return mdiff;
1751 }
1752 
1753 // Calculate the distance in xtals taking into account the periodicity of the
1754 // Barrel
1755 int DQMSourcePi0::diff_nphi_s(Int_t nphi1, Int_t nphi2) {
1756  Int_t mdiff;
1757  if (std::abs(nphi1 - nphi2) < (360 - std::abs(nphi1 - nphi2))) {
1758  mdiff = nphi1 - nphi2;
1759  } else {
1760  mdiff = 360 - std::abs(nphi1 - nphi2);
1761  if (nphi1 > nphi2)
1762  mdiff = -mdiff;
1763  }
1764  return mdiff;
1765 }
1766 
1767 DEFINE_FWK_MODULE(DQMSourcePi0);
DQMSourcePi0::ParameterT0_barl_
double ParameterT0_barl_
Definition: DQMSourcePi0.cc:288
DQMSourcePi0::productMonitoredEEeta_
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEEeta_
Definition: DQMSourcePi0.cc:223
DQMSourcePi0::seleEtaBeltDREndCap_
double seleEtaBeltDREndCap_
Definition: DQMSourcePi0.cc:283
nphi
const int nphi
Definition: CMTRawAnalyzer.cc:418
DQMSourcePi0::analyze
void analyze(const edm::Event &e, const edm::EventSetup &c) override
Definition: DQMSourcePi0.cc:586
DQMSourcePi0::hiYDistrEEeta_
MonitorElement * hiYDistrEEeta_
Distribution of rechits in iy EE (eta)
Definition: DQMSourcePi0.cc:83
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:30
DQMSourcePi0::seleMinvMinPi0EndCap_
double seleMinvMinPi0EndCap_
Definition: DQMSourcePi0.cc:255
DQMSourcePi0::hMinvPi0EE_
MonitorElement * hMinvPi0EE_
Pi0 invariant mass in EE.
Definition: DQMSourcePi0.cc:137
deDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
DQMSourcePi0::hiXDistrEEeta_
MonitorElement * hiXDistrEEeta_
Distribution of rechits in ix EE (eta)
Definition: DQMSourcePi0.cc:71
DQMSourcePi0::ParameterT0_endcPresh_
double ParameterT0_endcPresh_
Definition: DQMSourcePi0.cc:290
DQMSourcePi0::hiPhiDistrEBpi0_
MonitorElement * hiPhiDistrEBpi0_
Distribution of rechits in iPhi (pi0)
Definition: DQMSourcePi0.cc:62
mps_fire.i
i
Definition: mps_fire.py:429
DQMSourcePi0::hPt1Pi0EB_
MonitorElement * hPt1Pi0EB_
Pt of the 1st most energetic Pi0 photon in EB.
Definition: DQMSourcePi0.cc:146
DQMSourcePi0::hPt2Pi0EE_
MonitorElement * hPt2Pi0EE_
Pt of the 2nd most energetic Pi0 photon in EE.
Definition: DQMSourcePi0.cc:161
DQMSourcePi0::hiYDistrEEpi0_
MonitorElement * hiYDistrEEpi0_
Distribution of rechits in iy EE (pi0)
Definition: DQMSourcePi0.cc:77
DQMSourcePi0::ptMinForIsolationEtaEndCap_
double ptMinForIsolationEtaEndCap_
Definition: DQMSourcePi0.cc:281
dqm::implementation::NavigatorBase::setCurrentFolder
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
DQMSourcePi0::hPtPi0EE_
MonitorElement * hPtPi0EE_
Pi0 Pt in EE.
Definition: DQMSourcePi0.cc:173
DQMSourcePi0::hMinvEtaEB_
MonitorElement * hMinvEtaEB_
Eta invariant mass in EB.
Definition: DQMSourcePi0.cc:140
ecalRecHitGreater
bool ecalRecHitGreater(EcalRecHit x, EcalRecHit y)
Definition: DQMSourcePi0.cc:42
DQMSourcePi0::seleEtaBeltDR_
double seleEtaBeltDR_
Definition: DQMSourcePi0.cc:271
DQMSourcePi0::detIdEERecHits
std::vector< EEDetId > detIdEERecHits
Definition: DQMSourcePi0.cc:296
DQMSourcePi0::bookHistograms
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: DQMSourcePi0.cc:410
edm::SortedCollection::size
size_type size() const
Definition: SortedCollection.h:215
DQMSourcePi0::posCalculator_
PositionCalc posCalculator_
Definition: DQMSourcePi0.cc:59
DQMSourcePi0::seleMinvMinEtaEndCap_
double seleMinvMinEtaEndCap_
Definition: DQMSourcePi0.cc:280
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
DQMSourcePi0::hS4S92EtaEB_
MonitorElement * hS4S92EtaEB_
S4S9 of the 2nd most energetic eta photon.
Definition: DQMSourcePi0.cc:212
DQMSourcePi0::detIdEBRecHits
std::vector< EBDetId > detIdEBRecHits
Definition: DQMSourcePi0.cc:293
DQMSourcePi0::productMonitoredEBeta_
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEBeta_
Definition: DQMSourcePi0.cc:219
DQMSourcePi0::selePtGammaEndCap_
double selePtGammaEndCap_
for pi0->gg endcap
Definition: DQMSourcePi0.cc:252
DQMSourcePi0::hMeanRecHitEnergyEBeta_
MonitorElement * hMeanRecHitEnergyEBeta_
Distribution of Mean energy per rechit EB (eta)
Definition: DQMSourcePi0.cc:128
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
edm::Handle::product
T const * product() const
Definition: Handle.h:70
DQMSourcePi0::ParameterT0_endc_
double ParameterT0_endc_
Definition: DQMSourcePi0.cc:289
DQMSourcePi0::seleMinvMaxPi0EndCap_
double seleMinvMaxPi0EndCap_
Definition: DQMSourcePi0.cc:254
EEDetId::ix
int ix() const
Definition: EEDetId.h:77
DQMSourcePi0::folderName_
std::string folderName_
DQM folder name.
Definition: DQMSourcePi0.cc:303
edm::SortedCollection< EcalRecHit >::const_iterator
std::vector< EcalRecHit >::const_iterator const_iterator
Definition: SortedCollection.h:80
DQMSourcePi0::hMeanRecHitEnergyEEpi0_
MonitorElement * hMeanRecHitEnergyEEpi0_
Distribution of Mean energy per rechit EE (pi0)
Definition: DQMSourcePi0.cc:125
DQMSourcePi0::hNRecHitsEBpi0_
MonitorElement * hNRecHitsEBpi0_
Distribution of number of RecHits EB (pi0)
Definition: DQMSourcePi0.cc:110
DQMSourcePi0::hPt2EtaEE_
MonitorElement * hPt2EtaEE_
Pt of the 2nd most energetic Eta photon in EE.
Definition: DQMSourcePi0.cc:167
DQMSourcePi0::hRechitEnergyEBpi0_
MonitorElement * hRechitEnergyEBpi0_
Energy Distribution of rechits EB (pi0)
Definition: DQMSourcePi0.cc:86
DQMSourcePi0::fileName_
std::string fileName_
Output file name if required.
Definition: DQMSourcePi0.cc:315
DQMSourcePi0::clusSeedThr_
double clusSeedThr_
Definition: DQMSourcePi0.cc:234
DQMSourcePi0::hRechitEnergyEEpi0_
MonitorElement * hRechitEnergyEEpi0_
Energy Distribution of rechits EE (pi0)
Definition: DQMSourcePi0.cc:89
DQMSourcePi0::seleMinvMaxEtaEndCap_
double seleMinvMaxEtaEndCap_
Definition: DQMSourcePi0.cc:279
DQMSourcePi0::EERecHits
std::vector< EcalRecHit > EERecHits
Definition: DQMSourcePi0.cc:297
DQMSourcePi0::hPtEtaEB_
MonitorElement * hPtEtaEB_
Eta Pt in EB.
Definition: DQMSourcePi0.cc:176
DQMSourcePi0::selePi0Iso_
double selePi0Iso_
Definition: DQMSourcePi0.cc:248
DQMSourcePi0::hIsoEtaEB_
MonitorElement * hIsoEtaEB_
Eta Iso EB.
Definition: DQMSourcePi0.cc:188
DQMSourcePi0::selePtEtaEndCap_
double selePtEtaEndCap_
Definition: DQMSourcePi0.cc:278
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
DQMSourcePi0::saveToFile_
bool saveToFile_
Write to file.
Definition: DQMSourcePi0.cc:306
l1ctLayer2EG_cff.id
id
Definition: l1ctLayer2EG_cff.py:85
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:129
DQMSourcePi0::seleS4S9GammaEtaEndCap_
double seleS4S9GammaEtaEndCap_
Definition: DQMSourcePi0.cc:276
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
DQMSourcePi0::hIsoPi0EB_
MonitorElement * hIsoPi0EB_
Pi0 Iso EB.
Definition: DQMSourcePi0.cc:182
DQMSourcePi0::hS4S91Pi0EE_
MonitorElement * hS4S91Pi0EE_
S4S9 of the 1st most energetic pi0 photon EE.
Definition: DQMSourcePi0.cc:197
DQMSourcePi0::hS4S91EtaEB_
MonitorElement * hS4S91EtaEB_
S4S9 of the 1st most energetic eta photon.
Definition: DQMSourcePi0.cc:200
hcalRecHitTable_cff.energy
energy
Definition: hcalRecHitTable_cff.py:13
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:35
DQMSourcePi0::selePi0BeltDeta_
double selePi0BeltDeta_
Definition: DQMSourcePi0.cc:247
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
DQMSourcePi0::seleS4S9Gamma_
double seleS4S9Gamma_
Definition: DQMSourcePi0.cc:245
DQMSourcePi0::ptMinForIsolationEta_
double ptMinForIsolationEta_
Definition: DQMSourcePi0.cc:269
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
dqm::impl::MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:292
DQMSourcePi0::seleMinvMinPi0_
double seleMinvMinPi0_
Definition: DQMSourcePi0.cc:244
DQMSourcePi0::hMinvPi0EB_
MonitorElement * hMinvPi0EB_
Pi0 invariant mass in EB.
Definition: DQMSourcePi0.cc:134
DQMSourcePi0::seleMinvMaxEta_
double seleMinvMaxEta_
Definition: DQMSourcePi0.cc:267
DQMSourcePi0::hS4S92Pi0EE_
MonitorElement * hS4S92Pi0EE_
S4S9 of the 2nd most energetic pi0 photon EE.
Definition: DQMSourcePi0.cc:209
DQMSourcePi0::convxtalid
void convxtalid(int &, int &)
Definition: DQMSourcePi0.cc:1735
DQMSourcePi0::ptMinForIsolationEndCap_
double ptMinForIsolationEndCap_
Definition: DQMSourcePi0.cc:260
hcalRecHitTable_cff.ieta
ieta
Definition: hcalRecHitTable_cff.py:15
DQMSourcePi0::ParameterLogWeighted_
bool ParameterLogWeighted_
Definition: DQMSourcePi0.cc:286
iEvent
int iEvent
Definition: GenABIO.cc:224
DQMSourcePi0::seleS4S9GammaEta_
double seleS4S9GammaEta_
Definition: DQMSourcePi0.cc:265
DQMSourcePi0::hiEtaDistrEBpi0_
MonitorElement * hiEtaDistrEBpi0_
Distribution of rechits in iEta (pi0)
Definition: DQMSourcePi0.cc:74
DQMSourcePi0::selePi0BeltDetaEndCap_
double selePi0BeltDetaEndCap_
Definition: DQMSourcePi0.cc:259
neta
const int neta
Definition: CMTRawAnalyzer.cc:417
DQMSourcePi0::hS4S91Pi0EB_
MonitorElement * hS4S91Pi0EB_
S4S9 of the 1st most energetic pi0 photon.
Definition: DQMSourcePi0.cc:194
DQMSourcePi0::hiXDistrEEpi0_
MonitorElement * hiXDistrEEpi0_
Distribution of rechits in ix EE (pi0)
Definition: DQMSourcePi0.cc:65
DQMSourcePi0::hiEtaDistrEBeta_
MonitorElement * hiEtaDistrEBeta_
Distribution of rechits in iEta (eta)
Definition: DQMSourcePi0.cc:80
DQMSourcePi0::productMonitoredEBpi0_
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEBpi0_
object to monitor
Definition: DQMSourcePi0.cc:218
DQMSourcePi0::isMonEBpi0_
bool isMonEBpi0_
which subdet will be monitored
Definition: DQMSourcePi0.cc:309
DQMSourcePi0::ParameterX0_
double ParameterX0_
Definition: DQMSourcePi0.cc:287
DQMSourcePi0::selePtPi0EndCap_
double selePtPi0EndCap_
Definition: DQMSourcePi0.cc:253
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
DQMSourcePi0::seleMinvMinEta_
double seleMinvMinEta_
Definition: DQMSourcePi0.cc:268
DQMSourcePi0::selePtGammaEtaEndCap_
double selePtGammaEtaEndCap_
for eta->gg endcap
Definition: DQMSourcePi0.cc:275
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
DQMSourcePi0::diff_nphi_s
int diff_nphi_s(int, int)
Definition: DQMSourcePi0.cc:1755
DQMSourcePi0::ptMinForIsolation_
double ptMinForIsolation_
Definition: DQMSourcePi0.cc:249
DQMSourcePi0::diff_neta_s
int diff_neta_s(int, int)
Definition: DQMSourcePi0.cc:1747
DQMSourcePi0::seleXtalMinEnergy_
double seleXtalMinEnergy_
Definition: DQMSourcePi0.cc:231
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
DQMSourcePi0::DQMSourcePi0
DQMSourcePi0(const edm::ParameterSet &)
Definition: DQMSourcePi0.cc:324
PositionCalc::Calculate_Location
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=nullptr)
Definition: PositionCalc.h:65
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:262
CaloSubdetectorTopology::getWindow
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
Definition: CaloSubdetectorTopology.cc:4
DQMSourcePi0::seleEtaBeltDeta_
double seleEtaBeltDeta_
Definition: DQMSourcePi0.cc:272
DQMSourcePi0::hPtEtaEE_
MonitorElement * hPtEtaEE_
Eta Pt in EE.
Definition: DQMSourcePi0.cc:179
DQMSourcePi0::seleXtalMinEnergyEndCap_
double seleXtalMinEnergyEndCap_
Definition: DQMSourcePi0.cc:232
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
hcalRecHitTable_cff.iphi
iphi
Definition: hcalRecHitTable_cff.py:16
DQMSourcePi0::seleS9S25GammaEtaEndCap_
double seleS9S25GammaEtaEndCap_
Definition: DQMSourcePi0.cc:277
DQMSourcePi0::selePi0BeltDR_
double selePi0BeltDR_
Definition: DQMSourcePi0.cc:246
DQMSourcePi0::hEventEnergyEEpi0_
MonitorElement * hEventEnergyEEpi0_
Distribution of total event energy EE (pi0)
Definition: DQMSourcePi0.cc:101
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:267
DQMSourcePi0::caloGeomToken_
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeomToken_
Definition: DQMSourcePi0.cc:226
DQMSourcePi0::hNRecHitsEEeta_
MonitorElement * hNRecHitsEEeta_
Distribution of number of RecHits EE (eta)
Definition: DQMSourcePi0.cc:119
ALPAKA_ACCELERATOR_NAMESPACE::ecal::reconstruction::internal::endcap::ix
ALPAKA_FN_ACC ALPAKA_FN_INLINE uint32_t ix(uint32_t id)
Definition: KernelHelpers.dev.cc:128
DQMSourcePi0::selePtGammaEta_
double selePtGammaEta_
for eta->gg barrel
Definition: DQMSourcePi0.cc:263
DQMSourcePi0::productMonitoredEEpi0_
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEEpi0_
object to monitor
Definition: DQMSourcePi0.cc:222
DQMSourcePi0::caloTopoToken_
edm::ESGetToken< CaloTopology, CaloTopologyRecord > caloTopoToken_
Definition: DQMSourcePi0.cc:225
DQMSourcePi0::seleEtaBeltDetaEndCap_
double seleEtaBeltDetaEndCap_
Definition: DQMSourcePi0.cc:284
ecalPiZeroTask_cfi.posCalcParameters
posCalcParameters
Definition: ecalPiZeroTask_cfi.py:20
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
DQMSourcePi0::hiPhiDistrEBeta_
MonitorElement * hiPhiDistrEBeta_
Distribution of rechits in iPhi (eta)
Definition: DQMSourcePi0.cc:68
DQMSourcePi0::selePtEta_
double selePtEta_
Definition: DQMSourcePi0.cc:264
DQMSourcePi0::hRechitEnergyEBeta_
MonitorElement * hRechitEnergyEBeta_
Energy Distribution of rechits EB (eta)
Definition: DQMSourcePi0.cc:92
HLT_FULL_cff.seeds
seeds
Definition: HLT_FULL_cff.py:77699
DQMSourcePi0::hPt1EtaEE_
MonitorElement * hPt1EtaEE_
Pt of the 1st most energetic Eta photon in EE.
Definition: DQMSourcePi0.cc:155
DQMSourcePi0::hEventEnergyEBpi0_
MonitorElement * hEventEnergyEBpi0_
Distribution of total event energy EB (pi0)
Definition: DQMSourcePi0.cc:98
DQMSourcePi0::seleEtaIso_
double seleEtaIso_
Definition: DQMSourcePi0.cc:270
DQMSourcePi0::hPt1EtaEB_
MonitorElement * hPt1EtaEB_
Pt of the 1st most energetic Eta photon in EB.
Definition: DQMSourcePi0.cc:152
DQMSourcePi0::hEventEnergyEBeta_
MonitorElement * hEventEnergyEBeta_
Distribution of total event energy EB (eta)
Definition: DQMSourcePi0.cc:104
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
DQMSourcePi0::hIsoEtaEE_
MonitorElement * hIsoEtaEE_
Eta Iso EE.
Definition: DQMSourcePi0.cc:191
DQMSourcePi0::hS4S92Pi0EB_
MonitorElement * hS4S92Pi0EB_
S4S9 of the 2nd most energetic pi0 photon.
Definition: DQMSourcePi0.cc:206
DQMSourcePi0::selePi0IsoEndCap_
double selePi0IsoEndCap_
Definition: DQMSourcePi0.cc:257
DQMSourcePi0::ParameterW0_
double ParameterW0_
Definition: DQMSourcePi0.cc:291
DQMSourcePi0::hPt2EtaEB_
MonitorElement * hPt2EtaEB_
Pt of the 2nd most energetic Eta photon in EB.
Definition: DQMSourcePi0.cc:164
DQMSourcePi0::seleMinvMaxPi0_
double seleMinvMaxPi0_
Definition: DQMSourcePi0.cc:243
DQMSourcePi0::EBRecHits
std::vector< EcalRecHit > EBRecHits
Definition: DQMSourcePi0.cc:294
DQMSourcePi0::hMinvEtaEE_
MonitorElement * hMinvEtaEE_
Eta invariant mass in EE.
Definition: DQMSourcePi0.cc:143
DQMSourcePi0::hPtPi0EB_
MonitorElement * hPtPi0EB_
Pi0 Pt in EB.
Definition: DQMSourcePi0.cc:170
DQMSourcePi0::seleEtaIsoEndCap_
double seleEtaIsoEndCap_
Definition: DQMSourcePi0.cc:282
DQMSourcePi0::hS4S92EtaEE_
MonitorElement * hS4S92EtaEE_
S4S9 of the 2nd most energetic eta photon EE.
Definition: DQMSourcePi0.cc:215
DQMSourcePi0::prescaleFactor_
unsigned int prescaleFactor_
Monitor every prescaleFactor_ events.
Definition: DQMSourcePi0.cc:300
DQMSourcePi0::hMeanRecHitEnergyEEeta_
MonitorElement * hMeanRecHitEnergyEEeta_
Distribution of Mean energy per rechit EE (eta)
Definition: DQMSourcePi0.cc:131
DQMSourcePi0::hIsoPi0EE_
MonitorElement * hIsoPi0EE_
Pi0 Iso EE.
Definition: DQMSourcePi0.cc:185
DQMSourcePi0::selePtPi0_
double selePtPi0_
Definition: DQMSourcePi0.cc:242
x
float x
Definition: beamSpotDipStandalone.cc:55
DQMSourcePi0::seleS4S9GammaEndCap_
double seleS4S9GammaEndCap_
Definition: DQMSourcePi0.cc:256
DQMSourcePi0::hNRecHitsEBeta_
MonitorElement * hNRecHitsEBeta_
Distribution of number of RecHits EB (eta)
Definition: DQMSourcePi0.cc:116
DQMSourcePi0::~DQMSourcePi0
~DQMSourcePi0() override
Definition: DQMSourcePi0.cc:407
ALPAKA_ACCELERATOR_NAMESPACE::ecal::reconstruction::internal::endcap::iy
ALPAKA_FN_ACC ALPAKA_FN_INLINE uint32_t iy(uint32_t id)
Definition: KernelHelpers.dev.cc:130
DQMSourcePi0::hPt1Pi0EE_
MonitorElement * hPt1Pi0EE_
Pt of the 1st most energetic Pi0 photon in EE.
Definition: DQMSourcePi0.cc:149
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
DQMSourcePi0::hS4S91EtaEE_
MonitorElement * hS4S91EtaEE_
S4S9 of the 1st most energetic eta photon EE.
Definition: DQMSourcePi0.cc:203
dqm::implementation::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
DQMSourcePi0::hRechitEnergyEEeta_
MonitorElement * hRechitEnergyEEeta_
Energy Distribution of rechits EE (eta)
Definition: DQMSourcePi0.cc:95
DQMSourcePi0::hMeanRecHitEnergyEBpi0_
MonitorElement * hMeanRecHitEnergyEBpi0_
Distribution of Mean energy per rechit EB (pi0)
Definition: DQMSourcePi0.cc:122
DQMSourcePi0::hEventEnergyEEeta_
MonitorElement * hEventEnergyEEeta_
Distribution of total event energy EE (eta)
Definition: DQMSourcePi0.cc:107
DQMSourcePi0::hNRecHitsEEpi0_
MonitorElement * hNRecHitsEEpi0_
Distribution of number of RecHits EE (pi0)
Definition: DQMSourcePi0.cc:113
DQMSourcePi0::selePi0BeltDREndCap_
double selePi0BeltDREndCap_
Definition: DQMSourcePi0.cc:258
edm::Run
Definition: Run.h:45
DQMSourcePi0::seleS9S25GammaEta_
double seleS9S25GammaEta_
Definition: DQMSourcePi0.cc:266
DQMSourcePi0::hPt2Pi0EB_
MonitorElement * hPt2Pi0EB_
Pt of the 2nd most energetic Pi0 photon in EB.
Definition: DQMSourcePi0.cc:158
DQMSourcePi0::selePtGamma_
double selePtGamma_
Definition: DQMSourcePi0.cc:241
EEDetId::iy
int iy() const
Definition: EEDetId.h:83
dqm::impl::MonitorElement::setAxisTitle
virtual void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
Definition: MonitorElement.cc:801
DQMSourcePi0::clusSeedThrEndCap_
double clusSeedThrEndCap_
Definition: DQMSourcePi0.cc:238
Generated for CMSSW Reference Manual by doxygen 1.8.14