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