EgammaCoreTools/test/testEcalClusterToolsFWLite.C

0001 #if !defined(__CINT__) && !defined(__MAKECINT__)
0002 #include "DataFormats/FWLite/interface/Handle.h"
0003 #include "DataFormats/FWLite/interface/Event.h"
0004
0005 //Headers for ROOT
0006 #include "TFile.h"
0007 #include "TH1F.h"
0008 #include "TMath.h"
0009
0010 //Headers for the topology and cluster tools
0011 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0012 #include "Geometry/CaloTopology/interface/EcalBarrelHardcodedTopology.h"
0013 #include "Geometry/CaloTopology/interface/EcalEndcapHardcodedTopology.h"
0014 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0015
0016 //Headers for the data items
0017 #include "DataFormats/DetId/interface/DetId.h"
0018 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0019 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0020 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0021 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0022 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0023 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0024 #endif
0025
0026
0027 void testEcalClusterToolsFWLite() {
0028   std::cout <<"opening file"<<std::endl;
0029   TFile *file=TFile::Open("rfio:///castor/cern.ch/cms/store/user/meridian/meridian/SingleGammaPt35_DSZS_V1/SingleGammaPt35_DSZS_V1/00b02d884670d693cb397a1e0af88088/SingleGammaPt35_cfi_py_GEN_SIM_DIGI_L1_DIGI2RAW_RAW2DIGI_RECO_1.root");
0030
0031   TH1F iEtaiEtaEB("iEtaiEtaEB","iEtaiEtaEB",100,0.,0.03);
0032   TH1F iEtaiPhiEB("iEtaiPhiEB","iEtaiPhiEB",100,0.,0.03);
0033   TH1F iPhiiPhiEB("iPhiiPhiEB","iPhiiPhiEB",100,0.,0.03);
0034
0035   TH1F iEtaiEtaEE("iEtaiEtaEE","iEtaiEtaEE",100,0.,0.03);
0036   TH1F iEtaiPhiEE("iEtaiPhiEE","iEtaiPhiEE",100,0.,0.03);
0037   TH1F iPhiiPhiEE("iPhiiPhiEE","iPhiiPhiEE",100,0.,0.03);
0038
0039   fwlite::Event ev(file);
0040
0041   CaloTopology *topology=new CaloTopology();
0042   EcalBarrelHardcodedTopology* ebTopology=new EcalBarrelHardcodedTopology();
0043   EcalEndcapHardcodedTopology* eeTopology=new EcalEndcapHardcodedTopology();
0044   topology->setSubdetTopology(DetId::Ecal,EcalBarrel,ebTopology);
0045   topology->setSubdetTopology(DetId::Ecal,EcalEndcap,eeTopology);
0046
0047   for( ev.toBegin();
0048        ! ev.atEnd();
0049        ++ev) {
0050
0051     fwlite::Handle<reco::BasicClusterCollection > pEBClusters;
0052     pEBClusters.getByLabel(ev,"hybridSuperClusters","hybridBarrelBasicClusters");
0053     const reco::BasicClusterCollection *ebClusters = pEBClusters.ptr();
0054
0055     fwlite::Handle<reco::BasicClusterCollection > pEEClusters;
0056     pEEClusters.getByLabel(ev,"multi5x5BasicClusters","multi5x5EndcapBasicClusters");
0057     const reco::BasicClusterCollection *eeClusters = pEEClusters.ptr();
0058
0059     fwlite::Handle< EcalRecHitCollection > pEBRecHits;
0060     pEBRecHits.getByLabel( ev, "reducedEcalRecHitsEB","");
0061     const EcalRecHitCollection *ebRecHits = pEBRecHits.ptr();
0062
0063     fwlite::Handle< EcalRecHitCollection > pEERecHits;
0064     pEERecHits.getByLabel( ev, "reducedEcalRecHitsEE","");
0065     const EcalRecHitCollection *eeRecHits = pEERecHits.ptr();
0066
0067     std::cout << "========== BARREL ==========" << std::endl;
0068     for (reco::BasicClusterCollection::const_iterator it = ebClusters->begin(); it != ebClusters->end(); ++it ) {
0069       std::cout << "----- new cluster -----" << std::endl;
0070       std::cout << "----------------- size: " << (*it).size() << " energy: " << (*it).energy() << std::endl;
0071
0072       std::cout << "e1x3..................... " << EcalClusterTools::e1x3( *it, ebRecHits, topology ) << std::endl;
0073       std::cout << "e3x1..................... " << EcalClusterTools::e3x1( *it, ebRecHits, topology ) << std::endl;
0074       std::cout << "e1x5..................... " << EcalClusterTools::e1x5( *it, ebRecHits, topology ) << std::endl;
0075       //std::cout << "e5x1..................... " << EcalClusterTools::e5x1( *it, ebRecHits, topology ) << std::endl;
0076       std::cout << "e2x2..................... " << EcalClusterTools::e2x2( *it, ebRecHits, topology ) << std::endl;
0077       std::cout << "e3x3..................... " << EcalClusterTools::e3x3( *it, ebRecHits, topology ) << std::endl;
0078       std::cout << "e4x4..................... " << EcalClusterTools::e4x4( *it, ebRecHits, topology ) << std::endl;
0079       std::cout << "e5x5..................... " << EcalClusterTools::e5x5( *it, ebRecHits, topology ) << std::endl;
0080       std::cout << "e2x5Right................ " << EcalClusterTools::e2x5Right( *it, ebRecHits, topology ) << std::endl;
0081       std::cout << "e2x5Left................. " << EcalClusterTools::e2x5Left( *it, ebRecHits, topology ) << std::endl;
0082       std::cout << "e2x5Top.................. " << EcalClusterTools::e2x5Top( *it, ebRecHits, topology ) << std::endl;
0083       std::cout << "e2x5Bottom............... " << EcalClusterTools::e2x5Bottom( *it, ebRecHits, topology ) << std::endl;
0084       std::cout << "e2x5Max.................. " << EcalClusterTools::e2x5Max( *it, ebRecHits, topology ) << std::endl;
0085       std::cout << "eMax..................... " << EcalClusterTools::eMax( *it, ebRecHits ) << std::endl;
0086       std::cout << "e2nd..................... " << EcalClusterTools::e2nd( *it, ebRecHits ) << std::endl;
0087       std::vector<float> vEta = EcalClusterTools::energyBasketFractionEta( *it, ebRecHits );
0088       std::cout << "energyBasketFractionEta..";
0089       for (size_t i = 0; i < vEta.size(); ++i ) {
0090     std::cout << " " << vEta[i];
0091       }
0092       std::cout << std::endl;
0093       std::vector<float> vPhi = EcalClusterTools::energyBasketFractionPhi( *it, ebRecHits );
0094       std::cout << "energyBasketFractionPhi..";
0095       for (size_t i = 0; i < vPhi.size(); ++i ) {
0096     std::cout << " " << vPhi[i];
0097       }
0098       std::cout << std::endl;
0099       std::vector<float> vLocCov = EcalClusterTools::localCovariances( *it, ebRecHits, topology );
0100       std::cout << "local covariances........ " << vLocCov[0] << " " << vLocCov[1] << " " << vLocCov[2] << std::endl;
0101       if ((*it).energy() < 10)
0102     continue;
0103       iEtaiEtaEB.Fill(TMath::Sqrt(vLocCov[0]));
0104       iEtaiPhiEB.Fill(TMath::Sqrt(vLocCov[1]));
0105       iPhiiPhiEB.Fill(TMath::Sqrt(vLocCov[2]));
0106     }
0107
0108     std::cout << "========== ENDCAPS ==========" << std::endl;
0109     for (reco::BasicClusterCollection::const_iterator it = eeClusters->begin(); it != eeClusters->end(); ++it ) {
0110       std::cout << "----- new cluster -----" << std::endl;
0111       std::cout << "----------------- size: " << (*it).size() << " energy: " << (*it).energy() << std::endl;
0112
0113       std::cout << "e1x3..................... " << EcalClusterTools::e1x3( *it, eeRecHits, topology ) << std::endl;
0114       std::cout << "e3x1..................... " << EcalClusterTools::e3x1( *it, eeRecHits, topology ) << std::endl;
0115       std::cout << "e1x5..................... " << EcalClusterTools::e1x5( *it, eeRecHits, topology ) << std::endl;
0116       //std::cout << "e5x1..................... " << EcalClusterTools::e5x1( *it, eeRecHits, topology ) << std::endl;
0117       std::cout << "e2x2..................... " << EcalClusterTools::e2x2( *it, eeRecHits, topology ) << std::endl;
0118       std::cout << "e3x3..................... " << EcalClusterTools::e3x3( *it, eeRecHits, topology ) << std::endl;
0119       std::cout << "e4x4..................... " << EcalClusterTools::e4x4( *it, eeRecHits, topology ) << std::endl;
0120       std::cout << "e5x5..................... " << EcalClusterTools::e5x5( *it, eeRecHits, topology ) << std::endl;
0121       std::cout << "e2x5Right................ " << EcalClusterTools::e2x5Right( *it, eeRecHits, topology ) << std::endl;
0122       std::cout << "e2x5Left................. " << EcalClusterTools::e2x5Left( *it, eeRecHits, topology ) << std::endl;
0123       std::cout << "e2x5Top.................. " << EcalClusterTools::e2x5Top( *it, eeRecHits, topology ) << std::endl;
0124       std::cout << "e2x5Bottom............... " << EcalClusterTools::e2x5Bottom( *it, eeRecHits, topology ) << std::endl;
0125       std::cout << "eMax..................... " << EcalClusterTools::eMax( *it, eeRecHits ) << std::endl;
0126       std::cout << "e2nd..................... " << EcalClusterTools::e2nd( *it, eeRecHits ) << std::endl;
0127       std::vector<float> vLocCov = EcalClusterTools::localCovariances( *it, eeRecHits, topology );
0128       std::cout << "local covariances........ " << vLocCov[0] << " " << vLocCov[1] << " " << vLocCov[2] << std::endl;
0129       if ((*it).energy() < 10)
0130     continue;
0131       iEtaiEtaEE.Fill(TMath::Sqrt(vLocCov[0]));
0132       iEtaiPhiEE.Fill(TMath::Sqrt(vLocCov[1]));
0133       iPhiiPhiEE.Fill(TMath::Sqrt(vLocCov[2]));
0134
0135     }
0136
0137
0138   }
0139
0140   //Writing OutFile
0141   TFile outFile("locCov.root","RECREATE");
0142   outFile.cd();
0143   iEtaiEtaEB.Write();
0144   iEtaiPhiEB.Write();
0145   iPhiiPhiEB.Write();
0146   iEtaiEtaEE.Write();
0147   iEtaiPhiEE.Write();
0148   iPhiiPhiEE.Write();
0149   outFile.Write();
0150   outFile.Close();
0151
0152   delete topology;
0153 }