CMS 3D CMS Logo

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