Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:01

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 }