Reference Manual

Go to the documentation of this file.
00001 #include "FWCore/Framework/interface/Event.h"
00002 #include "FWCore/Framework/interface/EventSetup.h"
00003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00004 #include "FWCore/ServiceRegistry/interface/Service.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 // DQM include files
00008 
00009 #include "DQMServices/Core/interface/MonitorElement.h"
00010 
00011 // work on collections
00012 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00013 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00014 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00015 
00016 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00017 
00018 #include "DQMServices/Core/interface/DQMStore.h"
00019 
00020 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00021 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00022 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00023 #include "DataFormats/DetId/interface/DetId.h"
00024 #include "DataFormats/Math/interface/Point3D.h"
00025 
00026 #include "FWCore/Framework/interface/Event.h"
00027 #include "FWCore/Framework/interface/EventSetup.h"
00028 #include "DataFormats/Common/interface/Handle.h"
00029 #include "FWCore/Framework/interface/ESHandle.h"
00030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00031 
00032 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
00033 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
00035 #include "RecoEcal/EgammaCoreTools/interface/EcalEtaPhiRegion.h"
00036 
00037 #include "TVector3.h"
00038 
00039 #include "DQM/HLTEvF/interface/HLTAlCaMonPi0.h"
00040 
00041 
00042 #define TWOPI 6.283185308
00043 
00044 
00045 
00046 using namespace std;
00047 using namespace edm;
00048 
00049 
00050 // ******************************************
00051 // constructors
00052 // *****************************************
00053 
00054 HLTAlCaMonPi0::HLTAlCaMonPi0( const edm::ParameterSet& ps ) :
00055 eventCounter_(0)
00056 {
00057   dbe_ = Service<DQMStore>().operator->();
00058   folderName_ = ps.getUntrackedParameter<std::string>("FolderName","HLT/AlCaEcalPi0");
00059   prescaleFactor_ = ps.getUntrackedParameter<int>("prescaleFactor",1);
00060   productMonitoredEBpi0_= ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEBpi0Tag");
00061   productMonitoredEBeta_= ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEBetaTag");
00062   productMonitoredEEpi0_= ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEEpi0Tag");
00063   productMonitoredEEeta_= ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEEetaTag");
00064 
00065   isMonEBpi0_ = ps.getUntrackedParameter<bool>("isMonEBpi0",false);
00066   isMonEBeta_ = ps.getUntrackedParameter<bool>("isMonEBeta",false);
00067   isMonEEpi0_ = ps.getUntrackedParameter<bool>("isMonEEpi0",false);
00068   isMonEEeta_ = ps.getUntrackedParameter<bool>("isMonEEeta",false);
00069 
00070   saveToFile_=ps.getUntrackedParameter<bool>("SaveToFile",false);
00071   fileName_=  ps.getUntrackedParameter<std::string>("FileName","MonitorAlCaEcalPi0.root");
00072 
00073   clusSeedThr_ = ps.getParameter<double> ("clusSeedThr");
00074   clusSeedThrEndCap_ = ps.getParameter<double> ("clusSeedThrEndCap");
00075   clusEtaSize_ = ps.getParameter<int> ("clusEtaSize");
00076   clusPhiSize_ = ps.getParameter<int> ("clusPhiSize");
00077   if ( clusPhiSize_ % 2 == 0 ||  clusEtaSize_ % 2 == 0)
00078     edm::LogError("AlCaPi0RecHitsProducerError") << "Size of eta/phi for simple clustering should be odd numbers";
00079 
00080 
00081   seleXtalMinEnergy_ = ps.getParameter<double>("seleXtalMinEnergy");
00082   seleXtalMinEnergyEndCap_ = ps.getParameter<double>("seleXtalMinEnergyEndCap");
00083 
00085     selePtGamma_ = ps.getParameter<double> ("selePtGamma");  
00086     selePtPi0_ = ps.getParameter<double> ("selePtPi0");  
00087     seleMinvMaxPi0_ = ps.getParameter<double> ("seleMinvMaxPi0");  
00088     seleMinvMinPi0_ = ps.getParameter<double> ("seleMinvMinPi0");  
00089     seleS4S9Gamma_ = ps.getParameter<double> ("seleS4S9Gamma");  
00090     selePi0Iso_ = ps.getParameter<double> ("selePi0Iso");  
00091     ptMinForIsolation_ = ps.getParameter<double> ("ptMinForIsolation");
00092     selePi0BeltDR_ = ps.getParameter<double> ("selePi0BeltDR");  
00093     selePi0BeltDeta_ = ps.getParameter<double> ("selePi0BeltDeta");  
00094 
00095   
00097     selePtGammaEndCap_ = ps.getParameter<double> ("selePtGammaEndCap");  
00098     selePtPi0EndCap_ = ps.getParameter<double> ("selePtPi0EndCap");   
00099     seleS4S9GammaEndCap_ = ps.getParameter<double> ("seleS4S9GammaEndCap");  
00100     seleMinvMaxPi0EndCap_ = ps.getParameter<double> ("seleMinvMaxPi0EndCap");  
00101     seleMinvMinPi0EndCap_ = ps.getParameter<double> ("seleMinvMinPi0EndCap");  
00102     ptMinForIsolationEndCap_ = ps.getParameter<double> ("ptMinForIsolationEndCap");
00103     selePi0BeltDREndCap_ = ps.getParameter<double> ("selePi0BeltDREndCap");  
00104     selePi0BeltDetaEndCap_ = ps.getParameter<double> ("selePi0BeltDetaEndCap");  
00105     selePi0IsoEndCap_ = ps.getParameter<double> ("selePi0IsoEndCap");  
00106 
00107 
00109     selePtGammaEta_ = ps.getParameter<double> ("selePtGammaEta");  
00110     selePtEta_ = ps.getParameter<double> ("selePtEta");   
00111     seleS4S9GammaEta_ = ps.getParameter<double> ("seleS4S9GammaEta");  
00112     seleS9S25GammaEta_ = ps.getParameter<double> ("seleS9S25GammaEta");  
00113     seleMinvMaxEta_ = ps.getParameter<double> ("seleMinvMaxEta");  
00114     seleMinvMinEta_ = ps.getParameter<double> ("seleMinvMinEta");  
00115     ptMinForIsolationEta_ = ps.getParameter<double> ("ptMinForIsolationEta");
00116     seleEtaIso_ = ps.getParameter<double> ("seleEtaIso");  
00117     seleEtaBeltDR_ = ps.getParameter<double> ("seleEtaBeltDR");  
00118     seleEtaBeltDeta_ = ps.getParameter<double> ("seleEtaBeltDeta");  
00119 
00120   
00122     selePtGammaEtaEndCap_ = ps.getParameter<double> ("selePtGammaEtaEndCap");  
00123     selePtEtaEndCap_ = ps.getParameter<double> ("selePtEtaEndCap");   
00124     seleS4S9GammaEtaEndCap_ = ps.getParameter<double> ("seleS4S9GammaEtaEndCap");  
00125     seleS9S25GammaEtaEndCap_ = ps.getParameter<double> ("seleS9S25GammaEtaEndCap");  
00126     seleMinvMaxEtaEndCap_ = ps.getParameter<double> ("seleMinvMaxEtaEndCap");  
00127     seleMinvMinEtaEndCap_ = ps.getParameter<double> ("seleMinvMinEtaEndCap");  
00128     ptMinForIsolationEtaEndCap_ = ps.getParameter<double> ("ptMinForIsolationEtaEndCap");
00129     seleEtaIsoEndCap_ = ps.getParameter<double> ("seleEtaIsoEndCap");  
00130     seleEtaBeltDREndCap_ = ps.getParameter<double> ("seleEtaBeltDREndCap");  
00131     seleEtaBeltDetaEndCap_ = ps.getParameter<double> ("seleEtaBeltDetaEndCap");  
00132 
00133 
00134   // Parameters for the position calculation:
00135   edm::ParameterSet posCalcParameters = 
00136     ps.getParameter<edm::ParameterSet>("posCalcParameters");
00137   posCalculator_ = PositionCalc(posCalcParameters);
00138 
00139 }
00140 
00141 
00142 HLTAlCaMonPi0::~HLTAlCaMonPi0()
00143 {}
00144 
00145 
00146 //--------------------------------------------------------
00147 void HLTAlCaMonPi0::beginJob(){
00148 
00149 
00150   // create and cd into new folder
00151   dbe_->setCurrentFolder(folderName_);
00152 
00153   // book some histograms 1D
00154 
00155   hiPhiDistrEBpi0_ = dbe_->book1D("iphiDistributionEBpi0", "RechitEB pi0 iphi", 361, 1,361);
00156   hiPhiDistrEBpi0_->setAxisTitle("i#phi ", 1);
00157   hiPhiDistrEBpi0_->setAxisTitle("# rechits", 2);
00158 
00159   hiXDistrEEpi0_ = dbe_->book1D("iXDistributionEEpi0", "RechitEE pi0 ix", 100, 0,100);
00160   hiXDistrEEpi0_->setAxisTitle("ix ", 1);
00161   hiXDistrEEpi0_->setAxisTitle("# rechits", 2);
00162 
00163   hiPhiDistrEBeta_ = dbe_->book1D("iphiDistributionEBeta", "RechitEB eta iphi", 361, 1,361);
00164   hiPhiDistrEBeta_->setAxisTitle("i#phi ", 1);
00165   hiPhiDistrEBeta_->setAxisTitle("# rechits", 2);
00166 
00167   hiXDistrEEeta_ = dbe_->book1D("iXDistributionEEeta", "RechitEE eta ix", 100, 0,100);
00168   hiXDistrEEeta_->setAxisTitle("ix ", 1);
00169   hiXDistrEEeta_->setAxisTitle("# rechits", 2);
00170 
00171 
00172   hiEtaDistrEBpi0_ = dbe_->book1D("iEtaDistributionEBpi0", "RechitEB pi0 ieta", 171, -85, 86);
00173   hiEtaDistrEBpi0_->setAxisTitle("i#eta", 1);
00174   hiEtaDistrEBpi0_->setAxisTitle("#rechits", 2);
00175 
00176   hiYDistrEEpi0_ = dbe_->book1D("iYDistributionEBpi0", "RechitEE pi0 iY", 100, 0, 100);
00177   hiYDistrEEpi0_->setAxisTitle("iy", 1);
00178   hiYDistrEEpi0_->setAxisTitle("#rechits", 2);
00179 
00180   hiEtaDistrEBeta_ = dbe_->book1D("iEtaDistributionEBeta", "RechitEB eta ieta", 171, -85, 86);
00181   hiEtaDistrEBeta_->setAxisTitle("i#eta", 1);
00182   hiEtaDistrEBeta_->setAxisTitle("#rechits", 2);
00183 
00184   hiYDistrEEeta_ = dbe_->book1D("iYDistributionEBeta", "RechitEE eta iY", 100, 0, 100);
00185   hiYDistrEEeta_->setAxisTitle("iy", 1);
00186   hiYDistrEEeta_->setAxisTitle("#rechits", 2);
00187 
00188 
00189   hRechitEnergyEBpi0_ = dbe_->book1D("rhEnergyEBpi0","Pi0 rechits energy EB",160,0.,2.0);
00190   hRechitEnergyEBpi0_->setAxisTitle("energy (GeV) ",1);
00191   hRechitEnergyEBpi0_->setAxisTitle("#rechits",2);
00192 
00193   hRechitEnergyEEpi0_ = dbe_->book1D("rhEnergyEEpi0","Pi0 rechits energy EE",160,0.,3.0);
00194   hRechitEnergyEEpi0_->setAxisTitle("energy (GeV) ",1);
00195   hRechitEnergyEEpi0_->setAxisTitle("#rechits",2);
00196 
00197   hRechitEnergyEBeta_ = dbe_->book1D("rhEnergyEBeta","Eta rechits energy EB",160,0.,2.0);
00198   hRechitEnergyEBeta_->setAxisTitle("energy (GeV) ",1);
00199   hRechitEnergyEBeta_->setAxisTitle("#rechits",2);
00200 
00201   hRechitEnergyEEeta_ = dbe_->book1D("rhEnergyEEeta","Eta rechits energy EE",160,0.,3.0);
00202   hRechitEnergyEEeta_->setAxisTitle("energy (GeV) ",1);
00203   hRechitEnergyEEeta_->setAxisTitle("#rechits",2);
00204 
00205   hEventEnergyEBpi0_ = dbe_->book1D("eventEnergyEBpi0","Pi0 event energy EB",100,0.,20.0);
00206   hEventEnergyEBpi0_->setAxisTitle("energy (GeV) ",1);
00207 
00208   hEventEnergyEEpi0_ = dbe_->book1D("eventEnergyEEpi0","Pi0 event energy EE",100,0.,50.0);
00209   hEventEnergyEEpi0_->setAxisTitle("energy (GeV) ",1);
00210 
00211   hEventEnergyEBeta_ = dbe_->book1D("eventEnergyEBeta","Eta event energy EB",100,0.,20.0);
00212   hEventEnergyEBeta_->setAxisTitle("energy (GeV) ",1);
00213 
00214   hEventEnergyEEeta_ = dbe_->book1D("eventEnergyEEeta","Eta event energy EE",100,0.,50.0);
00215   hEventEnergyEEeta_->setAxisTitle("energy (GeV) ",1);
00216 
00217   hNRecHitsEBpi0_ = dbe_->book1D("nRechitsEBpi0","#rechits in pi0 collection EB",100,0.,250.);
00218   hNRecHitsEBpi0_->setAxisTitle("rechits ",1);
00219   
00220   hNRecHitsEEpi0_ = dbe_->book1D("nRechitsEEpi0","#rechits in pi0 collection EE",100,0.,250.);
00221   hNRecHitsEEpi0_->setAxisTitle("rechits ",1);
00222   
00223   hNRecHitsEBeta_ = dbe_->book1D("nRechitsEBeta","#rechits in eta collection EB",100,0.,250.);
00224   hNRecHitsEBeta_->setAxisTitle("rechits ",1);
00225   
00226   hNRecHitsEEeta_ = dbe_->book1D("nRechitsEEeta","#rechits in eta collection EE",100,0.,250.);
00227   hNRecHitsEEeta_->setAxisTitle("rechits ",1);
00228   
00229   hMeanRecHitEnergyEBpi0_ = dbe_->book1D("meanEnergyEBpi0","Mean rechit energy in pi0 collection EB",50,0.,2.);
00230   hMeanRecHitEnergyEBpi0_->setAxisTitle("Mean Energy [GeV] ",1);
00231   
00232   hMeanRecHitEnergyEEpi0_ = dbe_->book1D("meanEnergyEEpi0","Mean rechit energy in pi0 collection EE",100,0.,5.);
00233   hMeanRecHitEnergyEEpi0_->setAxisTitle("Mean Energy [GeV] ",1);
00234   
00235   hMeanRecHitEnergyEBeta_ = dbe_->book1D("meanEnergyEBeta","Mean rechit energy in eta collection EB",50,0.,2.);
00236   hMeanRecHitEnergyEBeta_->setAxisTitle("Mean Energy [GeV] ",1);
00237   
00238   hMeanRecHitEnergyEEeta_ = dbe_->book1D("meanEnergyEEeta","Mean rechit energy in eta collection EE",100,0.,5.);
00239   hMeanRecHitEnergyEEeta_->setAxisTitle("Mean Energy [GeV] ",1);
00240   
00241   hMinvPi0EB_ = dbe_->book1D("Pi0InvmassEB","Pi0 Invariant Mass in EB",100,0.,0.5);
00242   hMinvPi0EB_->setAxisTitle("Inv Mass [GeV] ",1);
00243 
00244   hMinvPi0EE_ = dbe_->book1D("Pi0InvmassEE","Pi0 Invariant Mass in EE",100,0.,0.5);
00245   hMinvPi0EE_->setAxisTitle("Inv Mass [GeV] ",1);
00246   
00247   hMinvEtaEB_ = dbe_->book1D("EtaInvmassEB","Eta Invariant Mass in EB",100,0.,0.85);
00248   hMinvEtaEB_->setAxisTitle("Inv Mass [GeV] ",1);
00249 
00250   hMinvEtaEE_ = dbe_->book1D("EtaInvmassEE","Eta Invariant Mass in EE",100,0.,0.85);
00251   hMinvEtaEE_->setAxisTitle("Inv Mass [GeV] ",1);
00252 
00253   
00254   hPt1Pi0EB_ = dbe_->book1D("Pt1Pi0EB","Pt 1st most energetic Pi0 photon in EB",100,0.,20.);
00255   hPt1Pi0EB_->setAxisTitle("1st photon Pt [GeV] ",1);
00256   
00257   hPt1Pi0EE_ = dbe_->book1D("Pt1Pi0EE","Pt 1st most energetic Pi0 photon in EE",100,0.,20.);
00258   hPt1Pi0EE_->setAxisTitle("1st photon Pt [GeV] ",1);
00259 
00260   hPt1EtaEB_ = dbe_->book1D("Pt1EtaEB","Pt 1st most energetic Eta photon in EB",100,0.,20.);
00261   hPt1EtaEB_->setAxisTitle("1st photon Pt [GeV] ",1);
00262   
00263   hPt1EtaEE_ = dbe_->book1D("Pt1EtaEE","Pt 1st most energetic Eta photon in EE",100,0.,20.);
00264   hPt1EtaEE_->setAxisTitle("1st photon Pt [GeV] ",1);
00265   
00266   hPt2Pi0EB_ = dbe_->book1D("Pt2Pi0EB","Pt 2nd most energetic Pi0 photon in EB",100,0.,20.);
00267   hPt2Pi0EB_->setAxisTitle("2nd photon Pt [GeV] ",1);
00268 
00269   hPt2Pi0EE_ = dbe_->book1D("Pt2Pi0EE","Pt 2nd most energetic Pi0 photon in EE",100,0.,20.);
00270   hPt2Pi0EE_->setAxisTitle("2nd photon Pt [GeV] ",1);
00271 
00272   hPt2EtaEB_ = dbe_->book1D("Pt2EtaEB","Pt 2nd most energetic Eta photon in EB",100,0.,20.);
00273   hPt2EtaEB_->setAxisTitle("2nd photon Pt [GeV] ",1);
00274 
00275   hPt2EtaEE_ = dbe_->book1D("Pt2EtaEE","Pt 2nd most energetic Eta photon in EE",100,0.,20.);
00276   hPt2EtaEE_->setAxisTitle("2nd photon Pt [GeV] ",1);
00277 
00278   
00279   hPtPi0EB_ = dbe_->book1D("PtPi0EB","Pi0 Pt in EB",100,0.,20.);
00280   hPtPi0EB_->setAxisTitle("Pi0 Pt [GeV] ",1);
00281 
00282   hPtPi0EE_ = dbe_->book1D("PtPi0EE","Pi0 Pt in EE",100,0.,20.);
00283   hPtPi0EE_->setAxisTitle("Pi0 Pt [GeV] ",1);
00284 
00285   hPtEtaEB_ = dbe_->book1D("PtEtaEB","Eta Pt in EB",100,0.,20.);
00286   hPtEtaEB_->setAxisTitle("Eta Pt [GeV] ",1);
00287 
00288   hPtEtaEE_ = dbe_->book1D("PtEtaEE","Eta Pt in EE",100,0.,20.);
00289   hPtEtaEE_->setAxisTitle("Eta Pt [GeV] ",1);
00290 
00291   hIsoPi0EB_ = dbe_->book1D("IsoPi0EB","Pi0 Iso in EB",50,0.,1.);
00292   hIsoPi0EB_->setAxisTitle("Pi0 Iso",1);
00293 
00294   hIsoPi0EE_ = dbe_->book1D("IsoPi0EE","Pi0 Iso in EE",50,0.,1.);
00295   hIsoPi0EE_->setAxisTitle("Pi0 Iso",1);
00296 
00297   hIsoEtaEB_ = dbe_->book1D("IsoEtaEB","Eta Iso in EB",50,0.,1.);
00298   hIsoEtaEB_->setAxisTitle("Eta Iso",1);
00299 
00300   hIsoEtaEE_ = dbe_->book1D("IsoEtaEE","Eta Iso in EE",50,0.,1.);
00301   hIsoEtaEE_->setAxisTitle("Eta Iso",1);
00302 
00303   hS4S91Pi0EB_ = dbe_->book1D("S4S91Pi0EB","S4S9 1st most energetic Pi0 photon in EB",50,0.,1.);
00304   hS4S91Pi0EB_->setAxisTitle("S4S9 of the 1st Pi0 Photon ",1);
00305 
00306   hS4S91Pi0EE_ = dbe_->book1D("S4S91Pi0EE","S4S9 1st most energetic Pi0 photon in EE",50,0.,1.);
00307   hS4S91Pi0EE_->setAxisTitle("S4S9 of the 1st Pi0 Photon ",1);
00308 
00309   hS4S91EtaEB_ = dbe_->book1D("S4S91EtaEB","S4S9 1st most energetic Eta photon in EB",50,0.,1.);
00310   hS4S91EtaEB_->setAxisTitle("S4S9 of the 1st Eta Photon ",1);
00311 
00312   hS4S91EtaEE_ = dbe_->book1D("S4S91EtaEE","S4S9 1st most energetic Eta photon in EE",50,0.,1.);
00313   hS4S91EtaEE_->setAxisTitle("S4S9 of the 1st Eta Photon ",1);
00314 
00315   hS4S92Pi0EB_ = dbe_->book1D("S4S92Pi0EB","S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
00316   hS4S92Pi0EB_->setAxisTitle("S4S9 of the 2nd Pi0 Photon",1);
00317 
00318   hS4S92Pi0EE_ = dbe_->book1D("S4S92Pi0EE","S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
00319   hS4S92Pi0EE_->setAxisTitle("S4S9 of the 2nd Pi0 Photon",1);
00320 
00321   hS4S92EtaEB_ = dbe_->book1D("S4S92EtaEB","S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
00322   hS4S92EtaEB_->setAxisTitle("S4S9 of the 2nd Eta Photon",1);
00323 
00324   hS4S92EtaEE_ = dbe_->book1D("S4S92EtaEE","S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
00325   hS4S92EtaEE_->setAxisTitle("S4S9 of the 2nd Eta Photon",1);
00326 
00327   
00328 
00329 
00330 }
00331 
00332 //--------------------------------------------------------
00333 void HLTAlCaMonPi0::beginRun(const edm::Run& r, const EventSetup& context) {
00334 
00335 }
00336 
00337 //--------------------------------------------------------
00338 void HLTAlCaMonPi0::beginLuminosityBlock(const LuminosityBlock& lumiSeg, 
00339      const EventSetup& context) {
00340   
00341 }
00342 
00343 //-------------------------------------------------------------
00344 
00345 void HLTAlCaMonPi0::analyze(const Event& iEvent, 
00346                                const EventSetup& iSetup ){  
00347  
00348   if (eventCounter_% prescaleFactor_ ) return; 
00349   eventCounter_++;
00350 
00351   edm::ESHandle<CaloTopology> theCaloTopology;
00352   iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
00353 
00354 
00355   std::vector<EcalRecHit> seeds;
00356   seeds.clear();
00357 
00358   std::vector<EBDetId> usedXtals;
00359   usedXtals.clear();
00360 
00361   detIdEBRecHits.clear(); 
00362   EBRecHits.clear();  
00363 
00364 
00365 
00366 
00367 
00368     
00369   edm::Handle<EcalRecHitCollection> rhEBpi0;
00370   edm::Handle<EcalRecHitCollection> rhEBeta;
00371   edm::Handle<EcalRecHitCollection> rhEEpi0;
00372   edm::Handle<EcalRecHitCollection> rhEEeta;
00373 
00374  bool GetRecHitsCollectionEBpi0 = true;
00375  if(isMonEBpi0_) {
00376    try { 
00377      iEvent.getByLabel(productMonitoredEBpi0_, rhEBpi0); 
00378    }catch( cms::Exception& exception ) {
00379      LogDebug("HLTAlCaPi0DQMSource") << "no EcalRecHits EB_pi0, can not run all stuffs" << iEvent.eventAuxiliary ().event() <<" run: "<<iEvent.eventAuxiliary ().run()  << std::endl;
00380      GetRecHitsCollectionEBpi0 = false;
00381    }
00382  }
00383 
00384  bool GetRecHitsCollectionEBeta = true;
00385  if(isMonEBeta_) {
00386    try { 
00387      if(isMonEBeta_) iEvent.getByLabel(productMonitoredEBeta_, rhEBeta); 
00388      if(isMonEEpi0_) iEvent.getByLabel(productMonitoredEEpi0_, rhEEpi0);
00389      if(isMonEEeta_) iEvent.getByLabel(productMonitoredEEeta_, rhEEeta);
00390    }catch( cms::Exception& exception ) {
00391      LogDebug("HLTAlCaPi0DQMSource") << "no EcalRecHits EB_Eta, can not run all stuffs" << iEvent.eventAuxiliary ().event() <<" run: "<<iEvent.eventAuxiliary ().run()  << std::endl;
00392      GetRecHitsCollectionEBeta = false;
00393    }
00394  }
00395 
00396  bool GetRecHitsCollectionEEpi0 = true;
00397  if(isMonEEpi0_) {
00398    try { 
00399      if(isMonEEpi0_) iEvent.getByLabel(productMonitoredEEpi0_, rhEEpi0);
00400    }catch( cms::Exception& exception ) {
00401      LogDebug("HLTAlCaPi0DQMSource") << "no EcalRecHits EE_Pi0, can not run all stuffs" << iEvent.eventAuxiliary ().event() <<" run: "<<iEvent.eventAuxiliary ().run()  << std::endl;
00402      GetRecHitsCollectionEEpi0 = false;
00403    }
00404  }
00405 
00406  bool GetRecHitsCollectionEEeta = true;
00407  if(isMonEEeta_) {
00408    try { 
00409      if(isMonEEeta_) iEvent.getByLabel(productMonitoredEEeta_, rhEEeta);
00410    }catch( cms::Exception& exception ) {
00411      LogDebug("HLTAlCaPi0DQMSource") << "no EcalRecHits EE_Eta, can not run all stuffs" << iEvent.eventAuxiliary ().event() <<" run: "<<iEvent.eventAuxiliary ().run()  << std::endl;
00412      GetRecHitsCollectionEEeta = false;
00413    }
00414  }
00415 
00416 
00417 
00418 
00419   // Initialize the Position Calc
00420   const CaloSubdetectorGeometry *geometry_p; 
00421   const CaloSubdetectorGeometry *geometryEE_p;    
00422   const CaloSubdetectorGeometry *geometryES_p;
00423   
00424   const CaloSubdetectorTopology *topology_p;
00425   const CaloSubdetectorTopology *topology_ee;
00426   
00427 
00428 
00429   edm::ESHandle<CaloGeometry> geoHandle;
00430   iSetup.get<CaloGeometryRecord>().get(geoHandle);     
00431   geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00432   geometryEE_p = geoHandle->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00433   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00434   topology_p = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalBarrel);
00435   topology_ee = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalEndcap);
00436   
00437   
00438   
00439   EcalRecHitCollection::const_iterator itb;
00440   
00441   // fill EB pi0 histos 
00442   if(isMonEBpi0_ ){
00443     if (rhEBpi0.isValid() && (rhEBpi0->size() > 0) && GetRecHitsCollectionEBpi0){
00444 
00445 
00446       const EcalRecHitCollection *hitCollection_p = rhEBpi0.product();
00447       float etot =0;
00448       for(itb=rhEBpi0->begin(); itb!=rhEBpi0->end(); ++itb){
00449         
00450         EBDetId id(itb->id());
00451         double energy = itb->energy();
00452         if( energy < seleXtalMinEnergy_) continue; 
00453 
00454         EBDetId det = itb->id();
00455 
00456 
00457         detIdEBRecHits.push_back(det);
00458         EBRecHits.push_back(*itb);
00459 
00460         if (energy > clusSeedThr_) seeds.push_back(*itb);
00461 
00462         hiPhiDistrEBpi0_->Fill(id.iphi());
00463         hiEtaDistrEBpi0_->Fill(id.ieta());
00464         hRechitEnergyEBpi0_->Fill(itb->energy());
00465         
00466         etot+= itb->energy();    
00467       } // Eb rechits
00468       
00469       hNRecHitsEBpi0_->Fill(rhEBpi0->size());
00470       hMeanRecHitEnergyEBpi0_->Fill(etot/rhEBpi0->size());
00471       hEventEnergyEBpi0_->Fill(etot);
00472 
00473       //      std::cout << " EB RH Pi0 collection: #, mean rh_e, event E "<<rhEBpi0->size()<<" "<<etot/rhEBpi0->size()<<" "<<etot<<std::endl;   
00474 
00475 
00476       // Pi0 maker
00477 
00478       //std::cout<< " RH coll size: "<<rhEBpi0->size()<<std::endl;
00479       //std::cout<< " Pi0 seeds: "<<seeds.size()<<std::endl;
00480 
00481       int nClus;
00482       std::vector<float> eClus;
00483       std::vector<float> etClus;
00484       std::vector<float> etaClus; 
00485       std::vector<float> thetaClus;
00486       std::vector<float> phiClus;
00487       std::vector<EBDetId> max_hit;
00488 
00489       std::vector< std::vector<EcalRecHit> > RecHitsCluster;
00490       std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5;
00491       std::vector<float> s4s9Clus;
00492       std::vector<float> s9s25Clus;
00493 
00494 
00495       nClus=0;
00496 
00497 
00498       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
00499       sort(seeds.begin(), seeds.end(), ecalRecHitLess());
00500 
00501       for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
00502         EBDetId seed_id = itseed->id();
00503         std::vector<EBDetId>::const_iterator usedIds;
00504 
00505         bool seedAlreadyUsed=false;
00506         for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
00507           if(*usedIds==seed_id){
00508             seedAlreadyUsed=true;
00509             //std::cout<< " Seed with energy "<<itseed->energy()<<" was used !"<<std::endl;
00510             break; 
00511           }
00512         }
00513         if(seedAlreadyUsed)continue;
00514         std::vector<DetId> clus_v = topology_p->getWindow(seed_id,clusEtaSize_,clusPhiSize_);       
00515         std::vector< std::pair<DetId, float> > clus_used;
00516         //Reject the seed if not able to build the cluster around it correctly
00517         //if(clus_v.size() < clusEtaSize_*clusPhiSize_){std::cout<<" Not enough RecHits "<<std::endl; continue;}
00518         vector<EcalRecHit> RecHitsInWindow;
00519         vector<EcalRecHit> RecHitsInWindow5x5;
00520 
00521         double simple_energy = 0; 
00522 
00523         for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
00524           EBDetId EBdet = *det;
00525           //      std::cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi "<<EBdet.iphi()<<std::endl;
00526           bool  HitAlreadyUsed=false;
00527           for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
00528             if(*usedIds==*det){
00529               HitAlreadyUsed=true;
00530               break;
00531             }
00532           }
00533           if(HitAlreadyUsed)continue;
00534 
00535 
00536           std::vector<EBDetId>::iterator itdet = find( detIdEBRecHits.begin(),detIdEBRecHits.end(),EBdet);
00537           if(itdet == detIdEBRecHits.end()) continue; 
00538       
00539           int nn = int(itdet - detIdEBRecHits.begin());
00540           usedXtals.push_back(*det);
00541           RecHitsInWindow.push_back(EBRecHits[nn]);
00542           clus_used.push_back(std::pair<DetId, float>(*det, 1) );
00543           simple_energy = simple_energy + EBRecHits[nn].energy();
00544 
00545 
00546         }
00547 
00548         if(simple_energy <= 0) continue;
00549    
00550         math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hitCollection_p,geometry_p,geometryES_p);
00551         //std::cout<< "       Simple Clustering: Total energy for this simple cluster : "<<simple_energy<<std::endl; 
00552         //std::cout<< "       Simple Clustering: eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<std::endl; 
00553         //std::cout<< "       Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<" "<<clus_pos.z()<<std::endl; 
00554 
00555         float theta_s = 2. * atan(exp(-clus_pos.eta()));
00556         //      float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
00557         //float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
00558         //    float p0z_s = simple_energy * cos(theta_s);
00559         //float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
00560         
00561         float et_s = simple_energy * sin(theta_s);
00562         //std::cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<" "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<std::endl;
00563 
00564         //Compute S4/S9 variable
00565         //We are not sure to have 9 RecHits so need to check eta and phi:
00567         float s4s9_tmp[4];
00568         for(int i=0;i<4;i++)s4s9_tmp[i]= 0;
00569 
00570         int seed_ieta = seed_id.ieta();
00571         int seed_iphi = seed_id.iphi();
00572     
00573         convxtalid( seed_iphi,seed_ieta);
00574     
00575         float e3x3 = 0; 
00576         float e5x5 = 0; 
00577 
00578         for(unsigned int j=0; j<RecHitsInWindow.size();j++){
00579           EBDetId det = (EBDetId)RecHitsInWindow[j].id(); 
00580       
00581           int ieta = det.ieta();
00582           int iphi = det.iphi();
00583       
00584           convxtalid(iphi,ieta);
00585       
00586           float en = RecHitsInWindow[j].energy(); 
00587       
00588           int dx = diff_neta_s(seed_ieta,ieta);
00589           int dy = diff_nphi_s(seed_iphi,iphi);
00590       
00591           if(dx <= 0 && dy <=0) s4s9_tmp[0] += en; 
00592           if(dx >= 0 && dy <=0) s4s9_tmp[1] += en; 
00593           if(dx <= 0 && dy >=0) s4s9_tmp[2] += en; 
00594           if(dx >= 0 && dy >=0) s4s9_tmp[3] += en; 
00595 
00596           if(abs(dx)<=1 && abs(dy)<=1) e3x3 += en; 
00597           if(abs(dx)<=2 && abs(dy)<=2) e5x5 += en; 
00598 
00599       
00600         }
00601     
00602         if(e3x3 <= 0) continue;
00603 
00604         float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3; 
00605 
00606 
00608     std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id,5,5); 
00609     for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
00610       EBDetId det = *idItr;
00611       
00612 
00613       std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(),usedXtals.end(),det);
00614 
00616       if(itdet0 != usedXtals.end()) continue; 
00617       
00618       //inside collections
00619       std::vector<EBDetId>::iterator itdet = find( detIdEBRecHits.begin(),detIdEBRecHits.end(),det);
00620       if(itdet == detIdEBRecHits.end()) continue; 
00621 
00622       int nn = int(itdet - detIdEBRecHits.begin());
00623       
00624       RecHitsInWindow5x5.push_back(EBRecHits[nn]);
00625       e5x5 += EBRecHits[nn].energy();
00626       
00627     }
00628 
00629 
00630 
00631         if(e5x5 <= 0) continue;
00632 
00633             eClus.push_back(simple_energy);
00634             etClus.push_back(et_s);
00635             etaClus.push_back(clus_pos.eta());
00636             thetaClus.push_back(theta_s);
00637             phiClus.push_back(clus_pos.phi());
00638             s4s9Clus.push_back(s4s9_max);
00639             s9s25Clus.push_back(e3x3/e5x5);
00640             RecHitsCluster.push_back(RecHitsInWindow);
00641             RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
00642             
00643             //      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;
00644 
00645             nClus++;
00646 
00647 
00648   }
00649     
00650       // std::cout<< " Pi0 clusters: "<<nClus<<std::endl;
00651 
00652       // Selection, based on Simple clustering
00653       //pi0 candidates
00654       int npi0_s=0;
00655 
00656 
00657       //      if (nClus <= 1) return;
00658       for(Int_t i=0 ; i<nClus ; i++){
00659         for(Int_t j=i+1 ; j<nClus ; j++){
00660           //      std::cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"  etClus[j] "<<etClus[j]<<std::endl;
00661           if( etClus[i]>selePtGamma_ && etClus[j]>selePtGamma_ && s4s9Clus[i]>seleS4S9Gamma_ && s4s9Clus[j]>seleS4S9Gamma_){
00662 
00663 
00664           float p0x = etClus[i] * cos(phiClus[i]);
00665           float p1x = etClus[j] * cos(phiClus[j]);
00666           float p0y = etClus[i] * sin(phiClus[i]);
00667           float p1y = etClus[j] * sin(phiClus[j]);
00668           float p0z = eClus[i] * cos(thetaClus[i]);
00669           float p1z = eClus[j] * cos(thetaClus[j]);
00670         
00671         
00672           float pt_pair = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
00673 
00674           if (pt_pair < selePtPi0_)continue;
00675 
00676           float m_inv = sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );  
00677             if ( (m_inv<seleMinvMaxPi0_) && (m_inv>seleMinvMinPi0_) ){
00678 
00679               //New Loop on cluster to measure isolation:
00680               std::vector<int> IsoClus;
00681               IsoClus.clear();
00682               float Iso = 0;
00683               TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
00684               for(Int_t k=0 ; k<nClus ; k++){
00685 
00686 
00687                 if(etClus[k] < ptMinForIsolation_) continue;
00688 
00689                 if(k==i || k==j)continue;
00690               TVector3 ClusVect = TVector3(etClus[k] *cos(phiClus[k]), etClus[k] * sin(phiClus[k]) , eClus[k] * cos(thetaClus[k]));
00691 
00692                 float dretacl = fabs(etaClus[k] - pairVect.Eta());
00693                 float drcl = ClusVect.DeltaR(pairVect);
00694                 //      std::cout<< "   Iso: k, E, drclpi0, detaclpi0, dphiclpi0 "<<k<<" "<<eClus[k]<<" "<<drclpi0<<" "<<dretaclpi0<<std::endl;
00695                 if((drcl<selePi0BeltDR_) && (dretacl<selePi0BeltDeta_) ){
00696                   //              std::cout<< "   ... good iso cluster #: "<<k<<" etClus[k] "<<etClus[k] <<std::endl;
00697                   Iso = Iso + etClus[k];
00698                   IsoClus.push_back(k);
00699                 }
00700               }
00701 
00702               //      std::cout<<"  Iso/pt_pi0 "<<Iso/pt_pi0<<std::endl;
00703               if(Iso/pt_pair<selePi0Iso_){
00704                 //for(unsigned int Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
00705                 //for(unsigned int Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
00706 
00707 
00708                 hMinvPi0EB_->Fill(m_inv);
00709                 hPt1Pi0EB_->Fill(etClus[i]);
00710                 hPt2Pi0EB_->Fill(etClus[j]);
00711                 hPtPi0EB_->Fill(pt_pair);
00712                 hIsoPi0EB_->Fill(Iso/pt_pair);
00713                 hS4S91Pi0EB_->Fill(s4s9Clus[i]);
00714                 hS4S92Pi0EB_->Fill(s4s9Clus[j]);
00715                 
00716                 //              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;  
00717 
00718                 npi0_s++;
00719                 
00720               }
00721             
00722 
00723 
00724           
00725             }
00726           }
00727         } // End of the "j" loop over Simple Clusters
00728       } // End of the "i" loop over Simple Clusters
00729 
00730       //      std::cout<<"  (Simple Clustering) EB Pi0 candidates #: "<<npi0_s<<std::endl;
00731 
00732     } // rhEBpi0.valid() ends
00733 
00734   } //  isMonEBpi0 ends
00735 
00736   //------------------ End of pi0 in EB --------------------------//
00737 
00738   // fill EB eta histos 
00739   if(isMonEBeta_ ){
00740     if (rhEBeta.isValid() && (rhEBeta->size() > 0) && GetRecHitsCollectionEBeta){
00741 
00742 
00743       const EcalRecHitCollection *hitCollection_p = rhEBeta.product();
00744       float etot =0;
00745       for(itb=rhEBeta->begin(); itb!=rhEBeta->end(); ++itb){
00746         
00747         EBDetId id(itb->id());
00748         double energy = itb->energy();
00749         if( energy < seleXtalMinEnergy_) continue; 
00750 
00751         EBDetId det = itb->id();
00752 
00753 
00754         detIdEBRecHits.push_back(det);
00755         EBRecHits.push_back(*itb);
00756 
00757         if (energy > clusSeedThr_) seeds.push_back(*itb);
00758 
00759         hiPhiDistrEBeta_->Fill(id.iphi());
00760         hiEtaDistrEBeta_->Fill(id.ieta());
00761         hRechitEnergyEBeta_->Fill(itb->energy());
00762         
00763         etot+= itb->energy();    
00764       } // Eb rechits
00765       
00766       hNRecHitsEBeta_->Fill(rhEBeta->size());
00767       hMeanRecHitEnergyEBeta_->Fill(etot/rhEBeta->size());
00768       hEventEnergyEBeta_->Fill(etot);
00769 
00770       //      std::cout << " EB RH Eta collection: #, mean rh_e, event E "<<rhEBeta->size()<<" "<<etot/rhEBeta->size()<<" "<<etot<<std::endl;   
00771 
00772 
00773       // Eta maker
00774 
00775       //std::cout<< " RH coll size: "<<rhEBeta->size()<<std::endl;
00776       //std::cout<< " Eta seeds: "<<seeds.size()<<std::endl;
00777 
00778       int nClus;
00779       std::vector<float> eClus;
00780       std::vector<float> etClus;
00781       std::vector<float> etaClus; 
00782       std::vector<float> thetaClus;
00783       std::vector<float> phiClus;
00784       std::vector<EBDetId> max_hit;
00785 
00786       std::vector< std::vector<EcalRecHit> > RecHitsCluster;
00787       std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5;
00788       std::vector<float> s4s9Clus;
00789       std::vector<float> s9s25Clus;
00790 
00791 
00792       nClus=0;
00793 
00794       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
00795       sort(seeds.begin(), seeds.end(), ecalRecHitLess());
00796 
00797       for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
00798         EBDetId seed_id = itseed->id();
00799         std::vector<EBDetId>::const_iterator usedIds;
00800 
00801         bool seedAlreadyUsed=false;
00802         for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
00803           if(*usedIds==seed_id){
00804             seedAlreadyUsed=true;
00805             //std::cout<< " Seed with energy "<<itseed->energy()<<" was used !"<<std::endl;
00806             break; 
00807           }
00808         }
00809         if(seedAlreadyUsed)continue;
00810         std::vector<DetId> clus_v = topology_p->getWindow(seed_id,clusEtaSize_,clusPhiSize_);       
00811         std::vector< std::pair<DetId, float> > clus_used;
00812         //Reject the seed if not able to build the cluster around it correctly
00813         //if(clus_v.size() < clusEtaSize_*clusPhiSize_){std::cout<<" Not enough RecHits "<<std::endl; continue;}
00814         vector<EcalRecHit> RecHitsInWindow;
00815         vector<EcalRecHit> RecHitsInWindow5x5;
00816 
00817         double simple_energy = 0; 
00818 
00819         for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
00820           EBDetId EBdet = *det;
00821           //      std::cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi "<<EBdet.iphi()<<std::endl;
00822           bool  HitAlreadyUsed=false;
00823           for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
00824             if(*usedIds==*det){
00825               HitAlreadyUsed=true;
00826               break;
00827             }
00828           }
00829           if(HitAlreadyUsed)continue;
00830 
00831 
00832           std::vector<EBDetId>::iterator itdet = find( detIdEBRecHits.begin(),detIdEBRecHits.end(),EBdet);
00833           if(itdet == detIdEBRecHits.end()) continue; 
00834       
00835           int nn = int(itdet - detIdEBRecHits.begin());
00836           usedXtals.push_back(*det);
00837           RecHitsInWindow.push_back(EBRecHits[nn]);
00838           clus_used.push_back(std::pair<DetId, float>(*det, 1));
00839           simple_energy = simple_energy + EBRecHits[nn].energy();
00840 
00841 
00842         }
00843 
00844         if(simple_energy <= 0) continue;
00845    
00846         math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hitCollection_p,geometry_p,geometryES_p);
00847         //std::cout<< "       Simple Clustering: Total energy for this simple cluster : "<<simple_energy<<std::endl; 
00848         //std::cout<< "       Simple Clustering: eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<std::endl; 
00849         //std::cout<< "       Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<" "<<clus_pos.z()<<std::endl; 
00850 
00851         float theta_s = 2. * atan(exp(-clus_pos.eta()));
00852         //      float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
00853         //float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
00854         //    float p0z_s = simple_energy * cos(theta_s);
00855         //float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
00856         
00857         float et_s = simple_energy * sin(theta_s);
00858         //std::cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<" "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<std::endl;
00859 
00860         //Compute S4/S9 variable
00861         //We are not sure to have 9 RecHits so need to check eta and phi:
00863         float s4s9_tmp[4];
00864         for(int i=0;i<4;i++)s4s9_tmp[i]= 0;
00865 
00866         int seed_ieta = seed_id.ieta();
00867         int seed_iphi = seed_id.iphi();
00868     
00869         convxtalid( seed_iphi,seed_ieta);
00870     
00871         float e3x3 = 0; 
00872         float e5x5 = 0; 
00873 
00874         for(unsigned int j=0; j<RecHitsInWindow.size();j++){
00875           EBDetId det = (EBDetId)RecHitsInWindow[j].id(); 
00876       
00877           int ieta = det.ieta();
00878           int iphi = det.iphi();
00879       
00880           convxtalid(iphi,ieta);
00881       
00882           float en = RecHitsInWindow[j].energy(); 
00883       
00884           int dx = diff_neta_s(seed_ieta,ieta);
00885           int dy = diff_nphi_s(seed_iphi,iphi);
00886       
00887           if(dx <= 0 && dy <=0) s4s9_tmp[0] += en; 
00888           if(dx >= 0 && dy <=0) s4s9_tmp[1] += en; 
00889           if(dx <= 0 && dy >=0) s4s9_tmp[2] += en; 
00890           if(dx >= 0 && dy >=0) s4s9_tmp[3] += en; 
00891 
00892           if(abs(dx)<=1 && abs(dy)<=1) e3x3 += en; 
00893           if(abs(dx)<=2 && abs(dy)<=2) e5x5 += en; 
00894 
00895       
00896         }
00897     
00898         if(e3x3 <= 0) continue;
00899 
00900         float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3; 
00901 
00902 
00904     std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id,5,5); 
00905     for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
00906       EBDetId det = *idItr;
00907       
00908 
00909       std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(),usedXtals.end(),det);
00910 
00912       if(itdet0 != usedXtals.end()) continue; 
00913       
00914       //inside collections
00915       std::vector<EBDetId>::iterator itdet = find( detIdEBRecHits.begin(),detIdEBRecHits.end(),det);
00916       if(itdet == detIdEBRecHits.end()) continue; 
00917 
00918       int nn = int(itdet - detIdEBRecHits.begin());
00919       
00920       RecHitsInWindow5x5.push_back(EBRecHits[nn]);
00921       e5x5 += EBRecHits[nn].energy();
00922       
00923     }
00924 
00925 
00926 
00927         if(e5x5 <= 0) continue;
00928 
00929             eClus.push_back(simple_energy);
00930             etClus.push_back(et_s);
00931             etaClus.push_back(clus_pos.eta());
00932             thetaClus.push_back(theta_s);
00933             phiClus.push_back(clus_pos.phi());
00934             s4s9Clus.push_back(s4s9_max);
00935             s9s25Clus.push_back(e3x3/e5x5);
00936             RecHitsCluster.push_back(RecHitsInWindow);
00937             RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
00938             
00939             //      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;
00940 
00941             nClus++;
00942             
00943 
00944   }
00945     
00946       // std::cout<< " Eta clusters: "<<nClus<<std::endl;
00947 
00948       // Selection, based on Simple clustering
00949       //eta candidates
00950       int npi0_s=0;
00951 
00952 
00953       //      if (nClus <= 1) return;
00954       for(Int_t i=0 ; i<nClus ; i++){
00955         for(Int_t j=i+1 ; j<nClus ; j++){
00956           //      std::cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"  etClus[j] "<<etClus[j]<<std::endl;
00957           if( etClus[i]>selePtGammaEta_ && etClus[j]>selePtGammaEta_ && s4s9Clus[i]>seleS4S9GammaEta_ && s4s9Clus[j]>seleS4S9GammaEta_){
00958 
00959 
00960           float p0x = etClus[i] * cos(phiClus[i]);
00961           float p1x = etClus[j] * cos(phiClus[j]);
00962           float p0y = etClus[i] * sin(phiClus[i]);
00963           float p1y = etClus[j] * sin(phiClus[j]);
00964           float p0z = eClus[i] * cos(thetaClus[i]);
00965           float p1z = eClus[j] * cos(thetaClus[j]);
00966         
00967         
00968           float pt_pair = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
00969 
00970           if (pt_pair < selePtEta_)continue;
00971 
00972           float m_inv = sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );  
00973             if ( (m_inv<seleMinvMaxEta_) && (m_inv>seleMinvMinEta_) ){
00974 
00975               //New Loop on cluster to measure isolation:
00976               std::vector<int> IsoClus;
00977               IsoClus.clear();
00978               float Iso = 0;
00979               TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
00980               for(Int_t k=0 ; k<nClus ; k++){
00981 
00982 
00983                 if(etClus[k] < ptMinForIsolationEta_) continue;
00984 
00985                 if(k==i || k==j)continue;
00986               TVector3 ClusVect = TVector3(etClus[k] *cos(phiClus[k]), etClus[k] * sin(phiClus[k]) , eClus[k] * cos(thetaClus[k]));
00987 
00988                 float dretacl = fabs(etaClus[k] - pairVect.Eta());
00989                 float drcl = ClusVect.DeltaR(pairVect);
00990                 //      std::cout<< "   Iso: k, E, drclpi0, detaclpi0, dphiclpi0 "<<k<<" "<<eClus[k]<<" "<<drclpi0<<" "<<dretaclpi0<<std::endl;
00991                 if((drcl<seleEtaBeltDR_) && (dretacl<seleEtaBeltDeta_) ){
00992                   //              std::cout<< "   ... good iso cluster #: "<<k<<" etClus[k] "<<etClus[k] <<std::endl;
00993                   Iso = Iso + etClus[k];
00994                   IsoClus.push_back(k);
00995                 }
00996               }
00997 
00998               //      std::cout<<"  Iso/pt_pi0 "<<Iso/pt_pi0<<std::endl;
00999               if(Iso/pt_pair<seleEtaIso_){
01000                 //for(unsigned int Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
01001                 //for(unsigned int Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
01002 
01003 
01004                 hMinvEtaEB_->Fill(m_inv);
01005                 hPt1EtaEB_->Fill(etClus[i]);
01006                 hPt2EtaEB_->Fill(etClus[j]);
01007                 hPtEtaEB_->Fill(pt_pair);
01008                 hIsoEtaEB_->Fill(Iso/pt_pair);
01009                 hS4S91EtaEB_->Fill(s4s9Clus[i]);
01010                 hS4S92EtaEB_->Fill(s4s9Clus[j]);
01011                 
01012                 //              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;  
01013 
01014                 npi0_s++;
01015                 
01016               }
01017             
01018 
01019 
01020           
01021             }
01022           }
01023         } // End of the "j" loop over Simple Clusters
01024       } // End of the "i" loop over Simple Clusters
01025 
01026       //      std::cout<<"  (Simple Clustering) EB Eta candidates #: "<<npi0_s<<std::endl;
01027 
01028     } // rhEBeta.valid() ends
01029 
01030   } //  isMonEBeta ends
01031 
01032   //------------------ End of Eta in EB --------------------------//
01033 
01034 
01035 
01036       //----------------- End of the EB --------------------------//
01037 
01038 
01039 
01040 
01041       //----------------- Start of the EE --------------------//
01042 
01043       // fill pi0 EE histos
01044       if(isMonEEpi0_){
01045         if (rhEEpi0.isValid() && (rhEEpi0->size() > 0) && GetRecHitsCollectionEEpi0){
01046 
01047           const EcalRecHitCollection *hitCollection_ee = rhEEpi0.product();
01048           float etot =0;
01049 
01050 
01051       detIdEERecHits.clear(); 
01052       EERecHits.clear();  
01053 
01054 
01055       std::vector<EcalRecHit> seedsEndCap;
01056       seedsEndCap.clear();
01057 
01058       std::vector<EEDetId> usedXtalsEndCap;
01059       usedXtalsEndCap.clear();
01060 
01061 
01063       EERecHitCollection::const_iterator ite;
01064       for (ite=rhEEpi0->begin(); ite!=rhEEpi0->end(); ite++) {
01065         double energy = ite->energy();
01066         if( energy < seleXtalMinEnergyEndCap_) continue; 
01067     
01068         EEDetId det = ite->id();
01069         EEDetId id(ite->id());
01070 
01071     
01072         detIdEERecHits.push_back(det);
01073         EERecHits.push_back(*ite);
01074     
01075         hiXDistrEEpi0_->Fill(id.ix());
01076         hiYDistrEEpi0_->Fill(id.iy());
01077         hRechitEnergyEEpi0_->Fill(ite->energy());
01078 
01079         if (energy > clusSeedThrEndCap_) seedsEndCap.push_back(*ite);
01080         
01081         etot+= ite->energy();    
01082       } // EE rechits
01083       
01084       hNRecHitsEEpi0_->Fill(rhEEpi0->size());
01085       hMeanRecHitEnergyEEpi0_->Fill(etot/rhEEpi0->size());
01086       hEventEnergyEEpi0_->Fill(etot);
01087 
01088       //      std::cout << " EE RH Pi0 collection: #, mean rh_e, event E "<<rhEEpi0->size()<<" "<<etot/rhEEpi0->size()<<" "<<etot<<std::endl;   
01089   
01090       int nClusEndCap;
01091       std::vector<float> eClusEndCap;
01092       std::vector<float> etClusEndCap;
01093       std::vector<float> etaClusEndCap;
01094       std::vector<float> thetaClusEndCap;
01095       std::vector<float> phiClusEndCap;
01096       std::vector< std::vector<EcalRecHit> > RecHitsClusterEndCap;
01097       std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5EndCap;
01098       std::vector<float> s4s9ClusEndCap;
01099       std::vector<float> s9s25ClusEndCap;
01100 
01101 
01102       nClusEndCap=0;
01103   
01104   
01105       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
01106       sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitLess());
01107   
01108       for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
01109         EEDetId seed_id = itseed->id();
01110         std::vector<EEDetId>::const_iterator usedIds;
01111     
01112         bool seedAlreadyUsed=false;
01113         for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
01114           if(*usedIds==seed_id){
01115             seedAlreadyUsed=true;
01116             break; 
01117           }
01118         }
01119 
01120         if(seedAlreadyUsed)continue;
01121         std::vector<DetId> clus_v = topology_ee->getWindow(seed_id,clusEtaSize_,clusPhiSize_);      
01122         std::vector< std::pair<DetId, float> > clus_used;   
01123 
01124         vector<EcalRecHit> RecHitsInWindow;
01125         vector<EcalRecHit> RecHitsInWindow5x5;
01126     
01127         float simple_energy = 0; 
01128 
01129         for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
01130           EEDetId EEdet = *det;
01131       
01132           bool  HitAlreadyUsed=false;
01133           for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
01134             if(*usedIds==*det){
01135               HitAlreadyUsed=true;
01136               break;
01137             }
01138           }
01139      
01140           if(HitAlreadyUsed)continue;
01141       
01142           std::vector<EEDetId>::iterator itdet = find( detIdEERecHits.begin(),detIdEERecHits.end(),EEdet);
01143           if(itdet == detIdEERecHits.end()) continue; 
01144       
01145       
01146           int nn = int(itdet - detIdEERecHits.begin());
01147           usedXtalsEndCap.push_back(*det);
01148           RecHitsInWindow.push_back(EERecHits[nn]);
01149           clus_used.push_back(std::pair<DetId, float>(*det, 1));
01150           simple_energy = simple_energy + EERecHits[nn].energy();
01151       
01152       
01153         }
01154     
01155         if( simple_energy <= 0) continue;
01156         
01157         math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hitCollection_ee,geometryEE_p,geometryES_p);
01158     
01159 
01160         float theta_s = 2. * atan(exp(-clus_pos.eta()));
01161         float et_s = simple_energy * sin(theta_s);
01162         //      float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
01163         //float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
01164         //float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
01165     
01166     
01167         //Compute S4/S9 variable
01168         //We are not sure to have 9 RecHits so need to check eta and phi:
01169         float s4s9_tmp[4];
01170         for(int i=0;i<4;i++) s4s9_tmp[i]= 0; 
01171     
01172         int ixSeed = seed_id.ix();
01173         int iySeed = seed_id.iy();
01174         float e3x3 = 0; 
01175         float e5x5 = 0;
01176 
01177         for(unsigned int j=0; j<RecHitsInWindow.size();j++){
01178           EEDetId det_this = (EEDetId)RecHitsInWindow[j].id(); 
01179           int dx = ixSeed - det_this.ix();
01180           int dy = iySeed - det_this.iy();
01181       
01182           float en = RecHitsInWindow[j].energy(); 
01183       
01184           if(dx <= 0 && dy <=0) s4s9_tmp[0] += en; 
01185           if(dx >= 0 && dy <=0) s4s9_tmp[1] += en; 
01186           if(dx <= 0 && dy >=0) s4s9_tmp[2] += en; 
01187           if(dx >= 0 && dy >=0) s4s9_tmp[3] += en; 
01188 
01189           if( abs(dx)<=1 && abs(dy)<=1) e3x3 += en; 
01190           if( abs(dx)<=2 && abs(dy)<=2) e5x5 += en; 
01191 
01192         }
01193 
01194       
01195         if(e3x3 <= 0) continue;
01196 
01197         eClusEndCap.push_back(simple_energy);
01198         etClusEndCap.push_back(et_s);
01199         etaClusEndCap.push_back(clus_pos.eta());
01200         thetaClusEndCap.push_back(theta_s);
01201         phiClusEndCap.push_back(clus_pos.phi());
01202         s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
01203         s9s25ClusEndCap.push_back(e3x3/e5x5);
01204         RecHitsClusterEndCap.push_back(RecHitsInWindow);
01205         RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
01206         
01207 
01208         //      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;
01209 
01210 
01211         nClusEndCap++;
01212 
01213 
01214       }
01215 
01216 
01217 
01218       // Selection, based on Simple clustering
01219       //pi0 candidates
01220       int npi0_se=0;
01221 
01222 
01223       for(Int_t i=0 ; i<nClusEndCap ; i++){
01224         for(Int_t j=i+1 ; j<nClusEndCap ; j++){
01225       
01226           if( etClusEndCap[i]>selePtGammaEndCap_ && etClusEndCap[j]>selePtGammaEndCap_ && s4s9ClusEndCap[i]>seleS4S9GammaEndCap_ && s4s9ClusEndCap[j]>seleS4S9GammaEndCap_){
01227 
01228           float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
01229           float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
01230           float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
01231           float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
01232           float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
01233           float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
01234 
01235         
01236             float pt_pair = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
01237             if (pt_pair < selePtPi0EndCap_)continue;
01238             float m_inv = sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );  
01239  
01240 
01241             if ( (m_inv<seleMinvMaxPi0EndCap_) && (m_inv>seleMinvMinPi0EndCap_) ){
01242           
01243 
01244               //New Loop on cluster to measure isolation:
01245               std::vector<int> IsoClus;
01246               IsoClus.clear();
01247               float Iso = 0;
01248               TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
01249               for(Int_t k=0 ; k<nClusEndCap ; k++){
01250 
01251                 if(etClusEndCap[k] < ptMinForIsolationEndCap_) continue; 
01252             
01253             
01254                 if(k==i || k==j)continue;
01255             
01256             
01257               TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]), etClusEndCap[k] * sin(phiClusEndCap[k]) , eClusEndCap[k] * cos(thetaClusEndCap[k]) ) ;
01258               float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
01259               float drcl = clusVect.DeltaR(pairVect);
01260               
01261               if(drcl < selePi0BeltDREndCap_ && dretacl < selePi0BeltDetaEndCap_ ){
01262                 Iso = Iso + etClusEndCap[k];
01263                 IsoClus.push_back(k);
01264               }
01265             }
01266 
01267           
01268               if(Iso/pt_pair<selePi0IsoEndCap_){
01269                 //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;  
01270 
01271                 hMinvPi0EE_->Fill(m_inv);
01272                 hPt1Pi0EE_->Fill(etClusEndCap[i]);
01273                 hPt2Pi0EE_->Fill(etClusEndCap[j]);
01274                 hPtPi0EE_->Fill(pt_pair);
01275                 hIsoPi0EE_->Fill(Iso/pt_pair);
01276                 hS4S91Pi0EE_->Fill(s4s9ClusEndCap[i]);
01277                 hS4S92Pi0EE_->Fill(s4s9ClusEndCap[j]);
01278 
01279                 npi0_se++;
01280               }
01281           
01282             }
01283           }
01284         } // End of the "j" loop over Simple Clusters
01285       } // End of the "i" loop over Simple Clusters
01286 
01287       //      std::cout<<"  (Simple Clustering) EE Pi0 candidates #: "<<npi0_se<<std::endl;
01288 
01289         } //rhEEpi0
01290       } //isMonEEpi0
01291 
01292       //================End of Pi0 endcap=======================//
01293 
01294 
01295       //================ Eta in EE===============================//
01296 
01297       // fill pi0 EE histos
01298       if(isMonEEeta_){
01299         if (rhEEeta.isValid() && (rhEEeta->size() > 0) && GetRecHitsCollectionEEeta){
01300 
01301           const EcalRecHitCollection *hitCollection_ee = rhEEeta.product();
01302           float etot =0;
01303 
01304       detIdEERecHits.clear(); 
01305       EERecHits.clear();  
01306 
01307 
01308       std::vector<EcalRecHit> seedsEndCap;
01309       seedsEndCap.clear();
01310 
01311       std::vector<EEDetId> usedXtalsEndCap;
01312       usedXtalsEndCap.clear();
01313 
01314 
01316       EERecHitCollection::const_iterator ite;
01317       for (ite=rhEEeta->begin(); ite!=rhEEeta->end(); ite++) {
01318         double energy = ite->energy();
01319         if( energy < seleXtalMinEnergyEndCap_) continue; 
01320     
01321         EEDetId det = ite->id();
01322         EEDetId id(ite->id());
01323 
01324     
01325         detIdEERecHits.push_back(det);
01326         EERecHits.push_back(*ite);
01327     
01328         hiXDistrEEeta_->Fill(id.ix());
01329         hiYDistrEEeta_->Fill(id.iy());
01330         hRechitEnergyEEeta_->Fill(ite->energy());
01331 
01332         if (energy > clusSeedThrEndCap_) seedsEndCap.push_back(*ite);
01333         
01334         etot+= ite->energy();    
01335       } // EE rechits
01336       
01337       hNRecHitsEEeta_->Fill(rhEEeta->size());
01338       hMeanRecHitEnergyEEeta_->Fill(etot/rhEEeta->size());
01339       hEventEnergyEEeta_->Fill(etot);
01340 
01341       //      std::cout << " EE RH Eta collection: #, mean rh_e, event E "<<rhEEeta->size()<<" "<<etot/rhEEeta->size()<<" "<<etot<<std::endl;   
01342   
01343       int nClusEndCap;
01344       std::vector<float> eClusEndCap;
01345       std::vector<float> etClusEndCap;
01346       std::vector<float> etaClusEndCap;
01347       std::vector<float> thetaClusEndCap;
01348       std::vector<float> phiClusEndCap;
01349       std::vector< std::vector<EcalRecHit> > RecHitsClusterEndCap;
01350       std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5EndCap;
01351       std::vector<float> s4s9ClusEndCap;
01352       std::vector<float> s9s25ClusEndCap;
01353 
01354 
01355       nClusEndCap=0;
01356   
01357   
01358       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
01359       sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitLess());
01360   
01361       for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
01362         EEDetId seed_id = itseed->id();
01363         std::vector<EEDetId>::const_iterator usedIds;
01364     
01365         bool seedAlreadyUsed=false;
01366         for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
01367           if(*usedIds==seed_id){
01368             seedAlreadyUsed=true;
01369             break; 
01370           }
01371         }
01372 
01373         if(seedAlreadyUsed)continue;
01374         std::vector<DetId> clus_v = topology_ee->getWindow(seed_id,clusEtaSize_,clusPhiSize_);      
01375         std::vector< std::pair<DetId, float> > clus_used;
01376 
01377         vector<EcalRecHit> RecHitsInWindow;
01378         vector<EcalRecHit> RecHitsInWindow5x5;
01379     
01380         float simple_energy = 0; 
01381 
01382         for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
01383           EEDetId EEdet = *det;
01384       
01385           bool  HitAlreadyUsed=false;
01386           for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
01387             if(*usedIds==*det){
01388               HitAlreadyUsed=true;
01389               break;
01390             }
01391           }
01392      
01393           if(HitAlreadyUsed)continue;
01394       
01395           std::vector<EEDetId>::iterator itdet = find( detIdEERecHits.begin(),detIdEERecHits.end(),EEdet);
01396           if(itdet == detIdEERecHits.end()) continue; 
01397       
01398       
01399           int nn = int(itdet - detIdEERecHits.begin());
01400           usedXtalsEndCap.push_back(*det);
01401           RecHitsInWindow.push_back(EERecHits[nn]);
01402           clus_used.push_back(std::pair<DetId, float>(*det, 1));
01403           simple_energy = simple_energy + EERecHits[nn].energy();
01404       
01405       
01406         }
01407     
01408         if( simple_energy <= 0) continue;
01409         
01410         math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hitCollection_ee,geometryEE_p,geometryES_p);
01411     
01412 
01413         float theta_s = 2. * atan(exp(-clus_pos.eta()));
01414         float et_s = simple_energy * sin(theta_s);
01415         //      float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
01416         //float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
01417         //float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
01418     
01419     
01420         //Compute S4/S9 variable
01421         //We are not sure to have 9 RecHits so need to check eta and phi:
01422         float s4s9_tmp[4];
01423         for(int i=0;i<4;i++) s4s9_tmp[i]= 0; 
01424     
01425         int ixSeed = seed_id.ix();
01426         int iySeed = seed_id.iy();
01427         float e3x3 = 0; 
01428         float e5x5 = 0;
01429 
01430         for(unsigned int j=0; j<RecHitsInWindow.size();j++){
01431           EEDetId det_this = (EEDetId)RecHitsInWindow[j].id(); 
01432           int dx = ixSeed - det_this.ix();
01433           int dy = iySeed - det_this.iy();
01434       
01435           float en = RecHitsInWindow[j].energy(); 
01436       
01437           if(dx <= 0 && dy <=0) s4s9_tmp[0] += en; 
01438           if(dx >= 0 && dy <=0) s4s9_tmp[1] += en; 
01439           if(dx <= 0 && dy >=0) s4s9_tmp[2] += en; 
01440           if(dx >= 0 && dy >=0) s4s9_tmp[3] += en; 
01441 
01442           if( abs(dx)<=1 && abs(dy)<=1) e3x3 += en; 
01443           if( abs(dx)<=2 && abs(dy)<=2) e5x5 += en; 
01444 
01445         }
01446 
01447       
01448         if(e3x3 <= 0) continue;
01449 
01450         eClusEndCap.push_back(simple_energy);
01451         etClusEndCap.push_back(et_s);
01452         etaClusEndCap.push_back(clus_pos.eta());
01453         thetaClusEndCap.push_back(theta_s);
01454         phiClusEndCap.push_back(clus_pos.phi());
01455         s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
01456         s9s25ClusEndCap.push_back(e3x3/e5x5);
01457         RecHitsClusterEndCap.push_back(RecHitsInWindow);
01458         RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
01459         
01460         //      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;
01461 
01462         nClusEndCap++;
01463 
01464       }
01465 
01466 
01467 
01468       // Selection, based on Simple clustering
01469       //pi0 candidates
01470       int npi0_se=0;
01471 
01472 
01473       for(Int_t i=0 ; i<nClusEndCap ; i++){
01474         for(Int_t j=i+1 ; j<nClusEndCap ; j++){
01475       
01476           if( etClusEndCap[i]>selePtGammaEtaEndCap_ && etClusEndCap[j]>selePtGammaEtaEndCap_ && s4s9ClusEndCap[i]>seleS4S9GammaEtaEndCap_ && s4s9ClusEndCap[j]>seleS4S9GammaEtaEndCap_){
01477 
01478           float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
01479           float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
01480           float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
01481           float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
01482           float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
01483           float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
01484 
01485         
01486             float pt_pair = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
01487             if (pt_pair < selePtEtaEndCap_)continue;
01488             float m_inv = sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );  
01489  
01490 
01491             if ( (m_inv<seleMinvMaxEtaEndCap_) && (m_inv>seleMinvMinEtaEndCap_) ){
01492           
01493 
01494               //New Loop on cluster to measure isolation:
01495               std::vector<int> IsoClus;
01496               IsoClus.clear();
01497               float Iso = 0;
01498               TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
01499               for(Int_t k=0 ; k<nClusEndCap ; k++){
01500 
01501                 if(etClusEndCap[k] < ptMinForIsolationEtaEndCap_) continue; 
01502             
01503             
01504                 if(k==i || k==j)continue;
01505             
01506             
01507               TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]), etClusEndCap[k] * sin(phiClusEndCap[k]) , eClusEndCap[k] * cos(thetaClusEndCap[k]) ) ;
01508               float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
01509               float drcl = clusVect.DeltaR(pairVect);
01510               
01511               if(drcl < seleEtaBeltDREndCap_ && dretacl < seleEtaBeltDetaEndCap_ ){
01512                 Iso = Iso + etClusEndCap[k];
01513                 IsoClus.push_back(k);
01514               }
01515             }
01516 
01517           
01518               if(Iso/pt_pair<seleEtaIsoEndCap_){
01519                 //      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;  
01520 
01521                 hMinvEtaEE_->Fill(m_inv);
01522                 hPt1EtaEE_->Fill(etClusEndCap[i]);
01523                 hPt2EtaEE_->Fill(etClusEndCap[j]);
01524                 hPtEtaEE_->Fill(pt_pair);
01525                 hIsoEtaEE_->Fill(Iso/pt_pair);
01526                 hS4S91EtaEE_->Fill(s4s9ClusEndCap[i]);
01527                 hS4S92EtaEE_->Fill(s4s9ClusEndCap[j]);
01528 
01529                 npi0_se++;
01530               }
01531           
01532             }
01533           }
01534         } // End of the "j" loop over Simple Clusters
01535       } // End of the "i" loop over Simple Clusters
01536 
01537       //      std::cout<<"  (Simple Clustering) EE Eta candidates #: "<<npi0_se<<std::endl;
01538 
01539         } //rhEEeta
01540       } //isMonEEeta
01541 
01542       //================End of Pi0 endcap=======================//
01543 
01544 
01545 
01546 
01548 
01549 
01550 } 
01551 
01552 
01553 
01554 
01555 //--------------------------------------------------------
01556 void HLTAlCaMonPi0::endLuminosityBlock(const LuminosityBlock& lumiSeg, 
01557                                           const EventSetup& context) {
01558 }
01559 //--------------------------------------------------------
01560 void HLTAlCaMonPi0::endRun(const Run& r, const EventSetup& context){
01561 
01562 }
01563 //--------------------------------------------------------
01564 void HLTAlCaMonPi0::endJob(){
01565 
01566   if(dbe_) {  
01567     if (saveToFile_) {
01568       dbe_->save(fileName_);
01569     }
01570   }
01571 }
01572 
01573 
01574 void HLTAlCaMonPi0::convxtalid(Int_t &nphi,Int_t &neta)
01575 {
01576   // Barrel only
01577   // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
01578   // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
01579   
01580   if(neta > 0) neta -= 1;
01581   if(nphi > 359) nphi=nphi-360;
01582   
01583   // final check
01584   if(nphi >359 || nphi <0 || neta< -85 || neta > 84)
01585     {
01586       std::cout <<" unexpected fatal error in HLTPi0RecHitsFilter::convxtalid "<<  nphi <<  " " << neta <<  " " 
01587                 <<std::endl;
01588       //exit(1);
01589     }
01590 } //end of convxtalid
01591 
01592 int HLTAlCaMonPi0::diff_neta_s(Int_t neta1, Int_t neta2){
01593   Int_t mdiff;
01594   mdiff=(neta1-neta2);
01595   return mdiff;
01596 }
01597 
01598 // Calculate the distance in xtals taking into account the periodicity of the Barrel
01599 int HLTAlCaMonPi0::diff_nphi_s(Int_t nphi1,Int_t nphi2) {
01600   Int_t mdiff;
01601   if(abs(nphi1-nphi2) < (360-abs(nphi1-nphi2))) {
01602     mdiff=nphi1-nphi2;
01603   }
01604   else {
01605     mdiff=360-abs(nphi1-nphi2);
01606     if(nphi1>nphi2) mdiff=-mdiff;
01607   }
01608   return mdiff;
01609 }
01610 
01611
/data/refman/pasoursint/CMSSW_5_3_0/src/DQM/HLTEvF/plugins/HLTAlCaMonPi0.cc
