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