CMS 3D CMS Logo

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