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