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  int npi0_s = 0;
854 
855  // if (nClus <= 1) return;
856  for (Int_t i = 0; i < nClus; i++) {
857  for (Int_t j = i + 1; j < nClus; j++) {
858  // cout<<" i "<<i<<" etClus[i] "<<etClus[i]<<" j "<<j<<"
859  // etClus[j] "<<etClus[j]<<endl;
860  if (etClus[i] > selePtGamma_ && etClus[j] > selePtGamma_ && s4s9Clus[i] > seleS4S9Gamma_ &&
861  s4s9Clus[j] > seleS4S9Gamma_) {
862  float p0x = etClus[i] * cos(phiClus[i]);
863  float p1x = etClus[j] * cos(phiClus[j]);
864  float p0y = etClus[i] * sin(phiClus[i]);
865  float p1y = etClus[j] * sin(phiClus[j]);
866  float p0z = eClus[i] * cos(thetaClus[i]);
867  float p1z = eClus[j] * cos(thetaClus[j]);
868 
869  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
870 
871  if (pt_pair < selePtPi0_)
872  continue;
873 
874  float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
875  (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
876  if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
877  // New Loop on cluster to measure isolation:
878  std::vector<int> IsoClus;
879  IsoClus.clear();
880  float Iso = 0;
881  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
882  for (Int_t k = 0; k < nClus; k++) {
883  if (etClus[k] < ptMinForIsolation_)
884  continue;
885 
886  if (k == i || k == j)
887  continue;
888  TVector3 ClusVect =
889  TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
890 
891  float dretacl = fabs(etaClus[k] - pairVect.Eta());
892  float drcl = ClusVect.DeltaR(pairVect);
893  // cout<< " Iso: k, E, drclpi0, detaclpi0, dphiclpi0
894  // "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
895  // "<<dretaclpi0<<endl;
896  if ((drcl < selePi0BeltDR_) && (dretacl < selePi0BeltDeta_)) {
897  // cout<< " ... good iso cluster #: "<<k<<"
898  // etClus[k] "<<etClus[k] <<endl;
899  Iso = Iso + etClus[k];
900  IsoClus.push_back(k);
901  }
902  }
903 
904  // cout<<" Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
905  if (Iso / pt_pair < selePi0Iso_) {
906  // for(unsigned int
907  // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
908  // for(unsigned int
909  // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
910 
911  hMinvPi0EB_->Fill(m_inv);
912  hPt1Pi0EB_->Fill(etClus[i]);
913  hPt2Pi0EB_->Fill(etClus[j]);
914  hPtPi0EB_->Fill(pt_pair);
915  hIsoPi0EB_->Fill(Iso / pt_pair);
916  hS4S91Pi0EB_->Fill(s4s9Clus[i]);
917  hS4S92Pi0EB_->Fill(s4s9Clus[j]);
918 
919  // cout <<" EB Simple Clustering: pi0 Candidate
920  // pt, eta, phi, Iso, m_inv, i, j : "<<pt_pair<<"
921  //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
922  //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
923 
924  npi0_s++;
925  }
926  }
927  }
928  } // End of the "j" loop over Simple Clusters
929  } // End of the "i" loop over Simple Clusters
930 
931  // cout<<" (Simple Clustering) EB Pi0 candidates #: "<<npi0_s<<endl;
932 
933  } // rhEBpi0.valid() ends
934 
935  } // isMonEBpi0 ends
936 
937  //------------------ End of pi0 in EB --------------------------//
938 
939  // fill EB eta histos
940  if (isMonEBeta_) {
941  if (rhEBeta.isValid() && (!rhEBeta->empty())) {
942  const EcalRecHitCollection *hitCollection_p = rhEBeta.product();
943  float etot = 0;
944  for (itb = rhEBeta->begin(); itb != rhEBeta->end(); ++itb) {
945  EBDetId id(itb->id());
946  double energy = itb->energy();
947  if (energy < seleXtalMinEnergy_)
948  continue;
949 
950  EBDetId det = itb->id();
951 
952  detIdEBRecHits.push_back(det);
953  EBRecHits.push_back(*itb);
954 
955  if (energy > clusSeedThr_)
956  seeds.push_back(*itb);
957 
958  hiPhiDistrEBeta_->Fill(id.iphi());
959  hiEtaDistrEBeta_->Fill(id.ieta());
960  hRechitEnergyEBeta_->Fill(itb->energy());
961 
962  etot += itb->energy();
963  } // Eb rechits
964 
965  hNRecHitsEBeta_->Fill(rhEBeta->size());
966  hMeanRecHitEnergyEBeta_->Fill(etot / rhEBeta->size());
967  hEventEnergyEBeta_->Fill(etot);
968 
969  // cout << " EB RH Eta collection: #, mean rh_e, event E
970  // "<<rhEBeta->size()<<" "<<etot/rhEBeta->size()<<" "<<etot<<endl;
971 
972  // Eta maker
973 
974  // cout<< " RH coll size: "<<rhEBeta->size()<<endl;
975  // cout<< " Eta seeds: "<<seeds.size()<<endl;
976 
977  int nClus;
978  std::vector<float> eClus;
979  std::vector<float> etClus;
980  std::vector<float> etaClus;
981  std::vector<float> thetaClus;
982  std::vector<float> phiClus;
983  std::vector<EBDetId> max_hit;
984 
985  std::vector<std::vector<EcalRecHit>> RecHitsCluster;
986  std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5;
987  std::vector<float> s4s9Clus;
988  std::vector<float> s9s25Clus;
989 
990  nClus = 0;
991 
992  // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
993  std::sort(seeds.begin(), seeds.end(), ecalRecHitGreater);
994 
995  for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
996  EBDetId seed_id = itseed->id();
997  std::vector<EBDetId>::const_iterator usedIds;
998 
999  bool seedAlreadyUsed = false;
1000  for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
1001  if (*usedIds == seed_id) {
1002  seedAlreadyUsed = true;
1003  // cout<< " Seed with energy "<<itseed->energy()<<" was used
1004  // !"<<endl;
1005  break;
1006  }
1007  }
1008  if (seedAlreadyUsed)
1009  continue;
1010  std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1011  std::vector<std::pair<DetId, float>> clus_used;
1012  // Reject the seed if not able to build the cluster around it correctly
1013  // if(clus_v.size() < clusEtaSize_*clusPhiSize_){cout<<" Not enough
1014  // RecHits "<<endl; continue;}
1015  std::vector<EcalRecHit> RecHitsInWindow;
1016  std::vector<EcalRecHit> RecHitsInWindow5x5;
1017 
1018  double simple_energy = 0;
1019 
1020  for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1021  EBDetId EBdet = *det;
1022  // cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi
1023  // "<<EBdet.iphi()<<endl;
1024  bool HitAlreadyUsed = false;
1025  for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
1026  if (*usedIds == *det) {
1027  HitAlreadyUsed = true;
1028  break;
1029  }
1030  }
1031  if (HitAlreadyUsed)
1032  continue;
1033 
1034  std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), EBdet);
1035  if (itdet == detIdEBRecHits.end())
1036  continue;
1037 
1038  int nn = int(itdet - detIdEBRecHits.begin());
1039  usedXtals.push_back(*det);
1040  RecHitsInWindow.push_back(EBRecHits[nn]);
1041  clus_used.push_back(std::make_pair(*det, 1));
1042  simple_energy = simple_energy + EBRecHits[nn].energy();
1043  }
1044 
1045  if (simple_energy <= 0)
1046  continue;
1047 
1048  math::XYZPoint clus_pos =
1049  posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
1050  // cout<< " Simple Clustering: Total energy for this simple
1051  // cluster : "<<simple_energy<<endl; cout<< " Simple Clustering:
1052  // eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<endl; cout<< "
1053  // Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<"
1054  // "<<clus_pos.z()<<endl;
1055 
1056  float theta_s = 2. * atan(exp(-clus_pos.eta()));
1057  // float p0x_s = simple_energy * sin(theta_s) *
1058  // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1059  // sin(clus_pos.phi());
1060  // float p0z_s = simple_energy * cos(theta_s);
1061  // float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1062 
1063  float et_s = simple_energy * sin(theta_s);
1064  // cout << " Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<"
1065  // "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
1066 
1067  // Compute S4/S9 variable
1068  // We are not sure to have 9 RecHits so need to check eta and phi:
1070  float s4s9_tmp[4];
1071  for (int i = 0; i < 4; i++)
1072  s4s9_tmp[i] = 0;
1073 
1074  int seed_ieta = seed_id.ieta();
1075  int seed_iphi = seed_id.iphi();
1076 
1077  convxtalid(seed_iphi, seed_ieta);
1078 
1079  float e3x3 = 0;
1080  float e5x5 = 0;
1081 
1082  for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1083  EBDetId det = (EBDetId)RecHitsInWindow[j].id();
1084 
1085  int ieta = det.ieta();
1086  int iphi = det.iphi();
1087 
1088  convxtalid(iphi, ieta);
1089 
1090  float en = RecHitsInWindow[j].energy();
1091 
1092  int dx = diff_neta_s(seed_ieta, ieta);
1093  int dy = diff_nphi_s(seed_iphi, iphi);
1094 
1095  if (dx <= 0 && dy <= 0)
1096  s4s9_tmp[0] += en;
1097  if (dx >= 0 && dy <= 0)
1098  s4s9_tmp[1] += en;
1099  if (dx <= 0 && dy >= 0)
1100  s4s9_tmp[2] += en;
1101  if (dx >= 0 && dy >= 0)
1102  s4s9_tmp[3] += en;
1103 
1104  if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1105  e3x3 += en;
1106  if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1107  e5x5 += en;
1108  }
1109 
1110  if (e3x3 <= 0)
1111  continue;
1112 
1113  float s4s9_max = *std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3;
1114 
1116  std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id, 5, 5);
1117  for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
1118  EBDetId det = *idItr;
1119 
1120  std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(), usedXtals.end(), det);
1121 
1123  if (itdet0 != usedXtals.end())
1124  continue;
1125 
1126  // inside collections
1127  std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), det);
1128  if (itdet == detIdEBRecHits.end())
1129  continue;
1130 
1131  int nn = int(itdet - detIdEBRecHits.begin());
1132 
1133  RecHitsInWindow5x5.push_back(EBRecHits[nn]);
1134  e5x5 += EBRecHits[nn].energy();
1135  }
1136 
1137  if (e5x5 <= 0)
1138  continue;
1139 
1140  eClus.push_back(simple_energy);
1141  etClus.push_back(et_s);
1142  etaClus.push_back(clus_pos.eta());
1143  thetaClus.push_back(theta_s);
1144  phiClus.push_back(clus_pos.phi());
1145  s4s9Clus.push_back(s4s9_max);
1146  s9s25Clus.push_back(e3x3 / e5x5);
1147  RecHitsCluster.push_back(RecHitsInWindow);
1148  RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
1149 
1150  // std::cout<<" EB Eta cluster (n,nxt,e,et eta,phi,s4s9)
1151  //"<<nClus<<" "<<int(RecHitsInWindow.size())<<" "<<eClus[nClus]<<" "<<"
1152  //"<<etClus[nClus]<<" "<<etaClus[nClus]<<" "<<phiClus[nClus]<<"
1153  //"<<s4s9Clus[nClus]<<std::endl;
1154 
1155  nClus++;
1156  }
1157 
1158  // cout<< " Eta clusters: "<<nClus<<endl;
1159 
1160  // Selection, based on Simple clustering
1161  // eta candidates
1162  int npi0_s = 0;
1163 
1164  // if (nClus <= 1) return;
1165  for (Int_t i = 0; i < nClus; i++) {
1166  for (Int_t j = i + 1; j < nClus; j++) {
1167  // cout<<" i "<<i<<" etClus[i] "<<etClus[i]<<" j "<<j<<"
1168  // etClus[j] "<<etClus[j]<<endl;
1169  if (etClus[i] > selePtGammaEta_ && etClus[j] > selePtGammaEta_ && s4s9Clus[i] > seleS4S9GammaEta_ &&
1170  s4s9Clus[j] > seleS4S9GammaEta_) {
1171  float p0x = etClus[i] * cos(phiClus[i]);
1172  float p1x = etClus[j] * cos(phiClus[j]);
1173  float p0y = etClus[i] * sin(phiClus[i]);
1174  float p1y = etClus[j] * sin(phiClus[j]);
1175  float p0z = eClus[i] * cos(thetaClus[i]);
1176  float p1z = eClus[j] * cos(thetaClus[j]);
1177 
1178  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1179 
1180  if (pt_pair < selePtEta_)
1181  continue;
1182 
1183  float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
1184  (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1185  if ((m_inv < seleMinvMaxEta_) && (m_inv > seleMinvMinEta_)) {
1186  // New Loop on cluster to measure isolation:
1187  std::vector<int> IsoClus;
1188  IsoClus.clear();
1189  float Iso = 0;
1190  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1191  for (Int_t k = 0; k < nClus; k++) {
1192  if (etClus[k] < ptMinForIsolationEta_)
1193  continue;
1194 
1195  if (k == i || k == j)
1196  continue;
1197  TVector3 ClusVect =
1198  TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
1199 
1200  float dretacl = fabs(etaClus[k] - pairVect.Eta());
1201  float drcl = ClusVect.DeltaR(pairVect);
1202  // cout<< " Iso: k, E, drclpi0, detaclpi0, dphiclpi0
1203  // "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
1204  // "<<dretaclpi0<<endl;
1205  if ((drcl < seleEtaBeltDR_) && (dretacl < seleEtaBeltDeta_)) {
1206  // cout<< " ... good iso cluster #: "<<k<<"
1207  // etClus[k] "<<etClus[k] <<endl;
1208  Iso = Iso + etClus[k];
1209  IsoClus.push_back(k);
1210  }
1211  }
1212 
1213  // cout<<" Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
1214  if (Iso / pt_pair < seleEtaIso_) {
1215  // for(unsigned int
1216  // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
1217  // for(unsigned int
1218  // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
1219 
1220  hMinvEtaEB_->Fill(m_inv);
1221  hPt1EtaEB_->Fill(etClus[i]);
1222  hPt2EtaEB_->Fill(etClus[j]);
1223  hPtEtaEB_->Fill(pt_pair);
1224  hIsoEtaEB_->Fill(Iso / pt_pair);
1225  hS4S91EtaEB_->Fill(s4s9Clus[i]);
1226  hS4S92EtaEB_->Fill(s4s9Clus[j]);
1227 
1228  // cout <<" EB Simple Clustering: Eta Candidate
1229  // pt, eta, phi, Iso, m_inv, i, j : "<<pt_pair<<"
1230  //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
1231  //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
1232 
1233  npi0_s++;
1234  }
1235  }
1236  }
1237  } // End of the "j" loop over Simple Clusters
1238  } // End of the "i" loop over Simple Clusters
1239 
1240  // cout<<" (Simple Clustering) EB Eta candidates #: "<<npi0_s<<endl;
1241 
1242  } // rhEBeta.valid() ends
1243 
1244  } // isMonEBeta ends
1245 
1246  //------------------ End of Eta in EB --------------------------//
1247 
1248  //----------------- End of the EB --------------------------//
1249 
1250  //----------------- Start of the EE --------------------//
1251 
1252  // fill pi0 EE histos
1253  if (isMonEEpi0_) {
1254  if (rhEEpi0.isValid() && (!rhEEpi0->empty())) {
1255  const EcalRecHitCollection *hitCollection_ee = rhEEpi0.product();
1256  float etot = 0;
1257 
1258  detIdEERecHits.clear();
1259  EERecHits.clear();
1260 
1261  std::vector<EcalRecHit> seedsEndCap;
1262  seedsEndCap.clear();
1263 
1264  std::vector<EEDetId> usedXtalsEndCap;
1265  usedXtalsEndCap.clear();
1266 
1268  EERecHitCollection::const_iterator ite;
1269  for (ite = rhEEpi0->begin(); ite != rhEEpi0->end(); ite++) {
1270  double energy = ite->energy();
1271  if (energy < seleXtalMinEnergyEndCap_)
1272  continue;
1273 
1274  EEDetId det = ite->id();
1275  EEDetId id(ite->id());
1276 
1277  detIdEERecHits.push_back(det);
1278  EERecHits.push_back(*ite);
1279 
1280  hiXDistrEEpi0_->Fill(id.ix());
1281  hiYDistrEEpi0_->Fill(id.iy());
1282  hRechitEnergyEEpi0_->Fill(ite->energy());
1283 
1284  if (energy > clusSeedThrEndCap_)
1285  seedsEndCap.push_back(*ite);
1286 
1287  etot += ite->energy();
1288  } // EE rechits
1289 
1290  hNRecHitsEEpi0_->Fill(rhEEpi0->size());
1291  hMeanRecHitEnergyEEpi0_->Fill(etot / rhEEpi0->size());
1292  hEventEnergyEEpi0_->Fill(etot);
1293 
1294  // cout << " EE RH Pi0 collection: #, mean rh_e, event E
1295  // "<<rhEEpi0->size()<<" "<<etot/rhEEpi0->size()<<" "<<etot<<endl;
1296 
1297  int nClusEndCap;
1298  std::vector<float> eClusEndCap;
1299  std::vector<float> etClusEndCap;
1300  std::vector<float> etaClusEndCap;
1301  std::vector<float> thetaClusEndCap;
1302  std::vector<float> phiClusEndCap;
1303  std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1304  std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1305  std::vector<float> s4s9ClusEndCap;
1306  std::vector<float> s9s25ClusEndCap;
1307 
1308  nClusEndCap = 0;
1309 
1310  // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1311  std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1312 
1313  for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1314  EEDetId seed_id = itseed->id();
1315  std::vector<EEDetId>::const_iterator usedIds;
1316 
1317  bool seedAlreadyUsed = false;
1318  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1319  if (*usedIds == seed_id) {
1320  seedAlreadyUsed = true;
1321  break;
1322  }
1323  }
1324 
1325  if (seedAlreadyUsed)
1326  continue;
1327  std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1328  std::vector<std::pair<DetId, float>> clus_used;
1329 
1330  std::vector<EcalRecHit> RecHitsInWindow;
1331  std::vector<EcalRecHit> RecHitsInWindow5x5;
1332 
1333  float simple_energy = 0;
1334 
1335  for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1336  EEDetId EEdet = *det;
1337 
1338  bool HitAlreadyUsed = false;
1339  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1340  if (*usedIds == *det) {
1341  HitAlreadyUsed = true;
1342  break;
1343  }
1344  }
1345 
1346  if (HitAlreadyUsed)
1347  continue;
1348 
1349  std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1350  if (itdet == detIdEERecHits.end())
1351  continue;
1352 
1353  int nn = int(itdet - detIdEERecHits.begin());
1354  usedXtalsEndCap.push_back(*det);
1355  RecHitsInWindow.push_back(EERecHits[nn]);
1356  clus_used.push_back(std::make_pair(*det, 1));
1357  simple_energy = simple_energy + EERecHits[nn].energy();
1358  }
1359 
1360  if (simple_energy <= 0)
1361  continue;
1362 
1363  math::XYZPoint clus_pos =
1364  posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1365 
1366  float theta_s = 2. * atan(exp(-clus_pos.eta()));
1367  float et_s = simple_energy * sin(theta_s);
1368  // float p0x_s = simple_energy * sin(theta_s) *
1369  // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1370  // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1371 
1372  // Compute S4/S9 variable
1373  // We are not sure to have 9 RecHits so need to check eta and phi:
1374  float s4s9_tmp[4];
1375  for (int i = 0; i < 4; i++)
1376  s4s9_tmp[i] = 0;
1377 
1378  int ixSeed = seed_id.ix();
1379  int iySeed = seed_id.iy();
1380  float e3x3 = 0;
1381  float e5x5 = 0;
1382 
1383  for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1384  EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1385  int dx = ixSeed - det_this.ix();
1386  int dy = iySeed - det_this.iy();
1387 
1388  float en = RecHitsInWindow[j].energy();
1389 
1390  if (dx <= 0 && dy <= 0)
1391  s4s9_tmp[0] += en;
1392  if (dx >= 0 && dy <= 0)
1393  s4s9_tmp[1] += en;
1394  if (dx <= 0 && dy >= 0)
1395  s4s9_tmp[2] += en;
1396  if (dx >= 0 && dy >= 0)
1397  s4s9_tmp[3] += en;
1398 
1399  if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1400  e3x3 += en;
1401  if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1402  e5x5 += en;
1403  }
1404 
1405  if (e3x3 <= 0)
1406  continue;
1407 
1408  eClusEndCap.push_back(simple_energy);
1409  etClusEndCap.push_back(et_s);
1410  etaClusEndCap.push_back(clus_pos.eta());
1411  thetaClusEndCap.push_back(theta_s);
1412  phiClusEndCap.push_back(clus_pos.phi());
1413  s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1414  s9s25ClusEndCap.push_back(e3x3 / e5x5);
1415  RecHitsClusterEndCap.push_back(RecHitsInWindow);
1416  RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1417 
1418  // std::cout<<" EE pi0 cluster (n,nxt,e,et eta,phi,s4s9)
1419  //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1420  //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1421  //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1422  //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1423 
1424  nClusEndCap++;
1425  }
1426 
1427  // Selection, based on Simple clustering
1428  // pi0 candidates
1429  int npi0_se = 0;
1430 
1431  for (Int_t i = 0; i < nClusEndCap; i++) {
1432  for (Int_t j = i + 1; j < nClusEndCap; j++) {
1433  if (etClusEndCap[i] > selePtGammaEndCap_ && etClusEndCap[j] > selePtGammaEndCap_ &&
1434  s4s9ClusEndCap[i] > seleS4S9GammaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEndCap_) {
1435  float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1436  float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1437  float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1438  float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1439  float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1440  float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1441 
1442  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1443  if (pt_pair < selePtPi0EndCap_)
1444  continue;
1445  float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1446  (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1447 
1448  if ((m_inv < seleMinvMaxPi0EndCap_) && (m_inv > seleMinvMinPi0EndCap_)) {
1449  // New Loop on cluster to measure isolation:
1450  std::vector<int> IsoClus;
1451  IsoClus.clear();
1452  float Iso = 0;
1453  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1454  for (Int_t k = 0; k < nClusEndCap; k++) {
1455  if (etClusEndCap[k] < ptMinForIsolationEndCap_)
1456  continue;
1457 
1458  if (k == i || k == j)
1459  continue;
1460 
1461  TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1462  etClusEndCap[k] * sin(phiClusEndCap[k]),
1463  eClusEndCap[k] * cos(thetaClusEndCap[k]));
1464  float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1465  float drcl = clusVect.DeltaR(pairVect);
1466 
1467  if (drcl < selePi0BeltDREndCap_ && dretacl < selePi0BeltDetaEndCap_) {
1468  Iso = Iso + etClusEndCap[k];
1469  IsoClus.push_back(k);
1470  }
1471  }
1472 
1473  if (Iso / pt_pair < selePi0IsoEndCap_) {
1474  // cout <<" EE Simple Clustering: pi0 Candidate pt, eta, phi,
1475  // Iso, m_inv, i, j : "<<pt_pair<<" "<<pairVect.Eta()<<"
1476  // "<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1477  // "<<endl;
1478 
1479  hMinvPi0EE_->Fill(m_inv);
1480  hPt1Pi0EE_->Fill(etClusEndCap[i]);
1481  hPt2Pi0EE_->Fill(etClusEndCap[j]);
1482  hPtPi0EE_->Fill(pt_pair);
1483  hIsoPi0EE_->Fill(Iso / pt_pair);
1484  hS4S91Pi0EE_->Fill(s4s9ClusEndCap[i]);
1485  hS4S92Pi0EE_->Fill(s4s9ClusEndCap[j]);
1486 
1487  npi0_se++;
1488  }
1489  }
1490  }
1491  } // End of the "j" loop over Simple Clusters
1492  } // End of the "i" loop over Simple Clusters
1493 
1494  // cout<<" (Simple Clustering) EE Pi0 candidates #:
1495  // "<<npi0_se<<endl;
1496 
1497  } // rhEEpi0
1498  } // isMonEEpi0
1499 
1500  //================End of Pi0 endcap=======================//
1501 
1502  //================ Eta in EE===============================//
1503 
1504  // fill pi0 EE histos
1505  if (isMonEEeta_) {
1506  if (rhEEeta.isValid() && (!rhEEeta->empty())) {
1507  const EcalRecHitCollection *hitCollection_ee = rhEEeta.product();
1508  float etot = 0;
1509 
1510  detIdEERecHits.clear();
1511  EERecHits.clear();
1512 
1513  std::vector<EcalRecHit> seedsEndCap;
1514  seedsEndCap.clear();
1515 
1516  std::vector<EEDetId> usedXtalsEndCap;
1517  usedXtalsEndCap.clear();
1518 
1520  EERecHitCollection::const_iterator ite;
1521  for (ite = rhEEeta->begin(); ite != rhEEeta->end(); ite++) {
1522  double energy = ite->energy();
1523  if (energy < seleXtalMinEnergyEndCap_)
1524  continue;
1525 
1526  EEDetId det = ite->id();
1527  EEDetId id(ite->id());
1528 
1529  detIdEERecHits.push_back(det);
1530  EERecHits.push_back(*ite);
1531 
1532  hiXDistrEEeta_->Fill(id.ix());
1533  hiYDistrEEeta_->Fill(id.iy());
1534  hRechitEnergyEEeta_->Fill(ite->energy());
1535 
1536  if (energy > clusSeedThrEndCap_)
1537  seedsEndCap.push_back(*ite);
1538 
1539  etot += ite->energy();
1540  } // EE rechits
1541 
1542  hNRecHitsEEeta_->Fill(rhEEeta->size());
1543  hMeanRecHitEnergyEEeta_->Fill(etot / rhEEeta->size());
1544  hEventEnergyEEeta_->Fill(etot);
1545 
1546  // cout << " EE RH Eta collection: #, mean rh_e, event E
1547  // "<<rhEEeta->size()<<" "<<etot/rhEEeta->size()<<" "<<etot<<endl;
1548 
1549  int nClusEndCap;
1550  std::vector<float> eClusEndCap;
1551  std::vector<float> etClusEndCap;
1552  std::vector<float> etaClusEndCap;
1553  std::vector<float> thetaClusEndCap;
1554  std::vector<float> phiClusEndCap;
1555  std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1556  std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1557  std::vector<float> s4s9ClusEndCap;
1558  std::vector<float> s9s25ClusEndCap;
1559 
1560  nClusEndCap = 0;
1561 
1562  // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1563  std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1564 
1565  for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1566  EEDetId seed_id = itseed->id();
1567  std::vector<EEDetId>::const_iterator usedIds;
1568 
1569  bool seedAlreadyUsed = false;
1570  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1571  if (*usedIds == seed_id) {
1572  seedAlreadyUsed = true;
1573  break;
1574  }
1575  }
1576 
1577  if (seedAlreadyUsed)
1578  continue;
1579  std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1580  std::vector<std::pair<DetId, float>> clus_used;
1581 
1582  std::vector<EcalRecHit> RecHitsInWindow;
1583  std::vector<EcalRecHit> RecHitsInWindow5x5;
1584 
1585  float simple_energy = 0;
1586 
1587  for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1588  EEDetId EEdet = *det;
1589 
1590  bool HitAlreadyUsed = false;
1591  for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1592  if (*usedIds == *det) {
1593  HitAlreadyUsed = true;
1594  break;
1595  }
1596  }
1597 
1598  if (HitAlreadyUsed)
1599  continue;
1600 
1601  std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1602  if (itdet == detIdEERecHits.end())
1603  continue;
1604 
1605  int nn = int(itdet - detIdEERecHits.begin());
1606  usedXtalsEndCap.push_back(*det);
1607  RecHitsInWindow.push_back(EERecHits[nn]);
1608  clus_used.push_back(std::make_pair(*det, 1));
1609  simple_energy = simple_energy + EERecHits[nn].energy();
1610  }
1611 
1612  if (simple_energy <= 0)
1613  continue;
1614 
1615  math::XYZPoint clus_pos =
1616  posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1617 
1618  float theta_s = 2. * atan(exp(-clus_pos.eta()));
1619  float et_s = simple_energy * sin(theta_s);
1620  // float p0x_s = simple_energy * sin(theta_s) *
1621  // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1622  // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1623 
1624  // Compute S4/S9 variable
1625  // We are not sure to have 9 RecHits so need to check eta and phi:
1626  float s4s9_tmp[4];
1627  for (int i = 0; i < 4; i++)
1628  s4s9_tmp[i] = 0;
1629 
1630  int ixSeed = seed_id.ix();
1631  int iySeed = seed_id.iy();
1632  float e3x3 = 0;
1633  float e5x5 = 0;
1634 
1635  for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1636  EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1637  int dx = ixSeed - det_this.ix();
1638  int dy = iySeed - det_this.iy();
1639 
1640  float en = RecHitsInWindow[j].energy();
1641 
1642  if (dx <= 0 && dy <= 0)
1643  s4s9_tmp[0] += en;
1644  if (dx >= 0 && dy <= 0)
1645  s4s9_tmp[1] += en;
1646  if (dx <= 0 && dy >= 0)
1647  s4s9_tmp[2] += en;
1648  if (dx >= 0 && dy >= 0)
1649  s4s9_tmp[3] += en;
1650 
1651  if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1652  e3x3 += en;
1653  if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1654  e5x5 += en;
1655  }
1656 
1657  if (e3x3 <= 0)
1658  continue;
1659 
1660  eClusEndCap.push_back(simple_energy);
1661  etClusEndCap.push_back(et_s);
1662  etaClusEndCap.push_back(clus_pos.eta());
1663  thetaClusEndCap.push_back(theta_s);
1664  phiClusEndCap.push_back(clus_pos.phi());
1665  s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1666  s9s25ClusEndCap.push_back(e3x3 / e5x5);
1667  RecHitsClusterEndCap.push_back(RecHitsInWindow);
1668  RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1669 
1670  // std::cout<<" EE Eta cluster (n,nxt,e,et eta,phi,s4s9)
1671  //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1672  //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1673  //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1674  //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1675 
1676  nClusEndCap++;
1677  }
1678 
1679  // Selection, based on Simple clustering
1680  // pi0 candidates
1681  int npi0_se = 0;
1682 
1683  for (Int_t i = 0; i < nClusEndCap; i++) {
1684  for (Int_t j = i + 1; j < nClusEndCap; j++) {
1685  if (etClusEndCap[i] > selePtGammaEtaEndCap_ && etClusEndCap[j] > selePtGammaEtaEndCap_ &&
1686  s4s9ClusEndCap[i] > seleS4S9GammaEtaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEtaEndCap_) {
1687  float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1688  float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1689  float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1690  float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1691  float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1692  float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1693 
1694  float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1695  if (pt_pair < selePtEtaEndCap_)
1696  continue;
1697  float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1698  (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1699 
1700  if ((m_inv < seleMinvMaxEtaEndCap_) && (m_inv > seleMinvMinEtaEndCap_)) {
1701  // New Loop on cluster to measure isolation:
1702  std::vector<int> IsoClus;
1703  IsoClus.clear();
1704  float Iso = 0;
1705  TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1706  for (Int_t k = 0; k < nClusEndCap; k++) {
1707  if (etClusEndCap[k] < ptMinForIsolationEtaEndCap_)
1708  continue;
1709 
1710  if (k == i || k == j)
1711  continue;
1712 
1713  TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1714  etClusEndCap[k] * sin(phiClusEndCap[k]),
1715  eClusEndCap[k] * cos(thetaClusEndCap[k]));
1716  float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1717  float drcl = clusVect.DeltaR(pairVect);
1718 
1719  if (drcl < seleEtaBeltDREndCap_ && dretacl < seleEtaBeltDetaEndCap_) {
1720  Iso = Iso + etClusEndCap[k];
1721  IsoClus.push_back(k);
1722  }
1723  }
1724 
1725  if (Iso / pt_pair < seleEtaIsoEndCap_) {
1726  // cout <<" EE Simple Clustering: Eta Candidate pt, eta,
1727  // phi, Iso, m_inv, i, j : "<<pt_pair<<" "<<pairVect.Eta()<<"
1728  //"<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1729  //"<<endl;
1730 
1731  hMinvEtaEE_->Fill(m_inv);
1732  hPt1EtaEE_->Fill(etClusEndCap[i]);
1733  hPt2EtaEE_->Fill(etClusEndCap[j]);
1734  hPtEtaEE_->Fill(pt_pair);
1735  hIsoEtaEE_->Fill(Iso / pt_pair);
1736  hS4S91EtaEE_->Fill(s4s9ClusEndCap[i]);
1737  hS4S92EtaEE_->Fill(s4s9ClusEndCap[j]);
1738 
1739  npi0_se++;
1740  }
1741  }
1742  }
1743  } // End of the "j" loop over Simple Clusters
1744  } // End of the "i" loop over Simple Clusters
1745 
1746  // cout<<" (Simple Clustering) EE Eta candidates #:
1747  // "<<npi0_se<<endl;
1748 
1749  } // rhEEeta
1750  } // isMonEEeta
1751 
1752  //================End of Pi0 endcap=======================//
1753 
1755 }
1756 
1757 void DQMSourcePi0::convxtalid(Int_t &nphi, Int_t &neta) {
1758  // Barrel only
1759  // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
1760  // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
1761 
1762  if (neta > 0)
1763  neta -= 1;
1764  if (nphi > 359)
1765  nphi = nphi - 360;
1766 
1767 } // end of convxtalid
1768 
1769 int DQMSourcePi0::diff_neta_s(Int_t neta1, Int_t neta2) {
1770  Int_t mdiff;
1771  mdiff = (neta1 - neta2);
1772  return mdiff;
1773 }
1774 
1775 // Calculate the distance in xtals taking into account the periodicity of the
1776 // Barrel
1777 int DQMSourcePi0::diff_nphi_s(Int_t nphi1, Int_t nphi2) {
1778  Int_t mdiff;
1779  if (std::abs(nphi1 - nphi2) < (360 - std::abs(nphi1 - nphi2))) {
1780  mdiff = nphi1 - nphi2;
1781  } else {
1782  mdiff = 360 - std::abs(nphi1 - nphi2);
1783  if (nphi1 > nphi2)
1784  mdiff = -mdiff;
1785  }
1786  return mdiff;
1787 }
1788 
1789 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:65
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
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
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
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
DetachedQuadStep_cff.seeds
seeds
Definition: DetachedQuadStep_cff.py:199
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:1757
DQMSourcePi0::ptMinForIsolationEndCap_
double ptMinForIsolationEndCap_
Definition: DQMSourcePi0.cc:260
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:1777
DQMSourcePi0::ptMinForIsolation_
double ptMinForIsolation_
Definition: DQMSourcePi0.cc:249
DQMSourcePi0::diff_neta_s
int diff_neta_s(int, int)
Definition: DQMSourcePi0.cc:1769
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
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
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
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
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
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
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