CMS 3D CMS Logo

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