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