CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10/src/DQMOffline/Ecal/src/EERecoSummary.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    EERecoSummary
00004 // Class:      EERecoSummary
00005 // Original Author:  Martina Malberti
00006 // 
00007 // system include files
00008 #include <memory>
00009 
00010 // user include files
00011 #include "FWCore/Framework/interface/Frameworkfwd.h"
00012 #include "FWCore/Framework/interface/EDAnalyzer.h"
00013 
00014 #include "FWCore/Framework/interface/Event.h"
00015 #include "FWCore/Framework/interface/MakerMacros.h"
00016 #include "FWCore/Common/interface/EventBase.h"
00017 #include "DataFormats/Provenance/interface/EventAuxiliary.h"
00018 
00019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00020 #include "FWCore/ServiceRegistry/interface/Service.h"
00021 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00022 
00023 
00024 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00025 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00026 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00027 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00028 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
00029 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00030 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
00031 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
00032 
00033 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00034 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00035 #include "DataFormats/EcalDetId/interface/ESDetId.h"
00036 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00037 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00038 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00039 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00040 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00041 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
00042 #include "DataFormats/EgammaReco/interface/PreshowerClusterFwd.h"
00043 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00044 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
00045 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00046 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalCleaningAlgo.h"
00047 #include "RecoEcal/EgammaCoreTools/interface/EcalRecHitLess.h"
00048 
00049 #include "DataFormats/TrackReco/interface/Track.h"
00050 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00051 
00052 #include "DataFormats/JetReco/interface/CaloJet.h"
00053 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00054 
00055 #include "DQMOffline/Ecal/interface/EERecoSummary.h"
00056 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00057 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00058 #include "MagneticField/Engine/interface/MagneticField.h"
00059 
00060 #include "TVector3.h"
00061 
00062 #include <iostream>
00063 #include <cmath>
00064 #include <fstream>
00065 
00066 //
00067 // constructors and destructor
00068 //
00069 EERecoSummary::EERecoSummary(const edm::ParameterSet& ps)
00070 {
00071 
00072   prefixME_ = ps.getUntrackedParameter<std::string>("prefixME", "");
00073 
00074   //now do what ever initialization is needed
00075   recHitCollection_EE_       = ps.getParameter<edm::InputTag>("recHitCollection_EE");
00076   redRecHitCollection_EE_    = ps.getParameter<edm::InputTag>("redRecHitCollection_EE");
00077   basicClusterCollection_EE_ = ps.getParameter<edm::InputTag>("basicClusterCollection_EE");
00078   superClusterCollection_EE_ = ps.getParameter<edm::InputTag>("superClusterCollection_EE");
00079 
00080   ethrEE_                    = ps.getParameter<double>("ethrEE");
00081 
00082   scEtThrEE_                 = ps.getParameter<double>("scEtThrEE");
00083 
00084   // DQM Store -------------------
00085   dqmStore_ = edm::Service<DQMStore>().operator->();
00086 
00087   // Monitor Elements (ex THXD)
00088   dqmStore_->setCurrentFolder(prefixME_ + "/EERecoSummary"); // to organise the histos in folders
00089      
00090   // ReducedRecHits ----------------------------------------------
00091   // ... endcap 
00092   h_redRecHits_EE_recoFlag = dqmStore_->book1D("redRecHits_EE_recoFlag","redRecHits_EE_recoFlag",16,-0.5,15.5);  
00093 
00094   // RecHits ---------------------------------------------- 
00095   // ... endcap
00096   h_recHits_EE_recoFlag = dqmStore_->book1D("recHits_EE_recoFlag","recHits_EE_recoFlag",16,-0.5,15.5);  
00097 
00098   // ... endcap +
00099   h_recHits_EEP_energyMax     = dqmStore_->book1D("recHits_EEP_energyMax","recHitsEEP_energyMax",110,-10,100);
00100   h_recHits_EEP_Chi2          = dqmStore_->book1D("recHits_EEP_Chi2","recHits_EEP_Chi2",200,0,100);
00101   h_recHits_EEP_time          = dqmStore_->book1D("recHits_EEP_time","recHits_EEP_time",200,-50,50);
00102 
00103   // ... endcap -
00104   h_recHits_EEM_energyMax     = dqmStore_->book1D("recHits_EEM_energyMax","recHits_EEM_energyMax",110,-10,100);
00105   h_recHits_EEM_Chi2          = dqmStore_->book1D("recHits_EEM_Chi2","recHits_EEM_Chi2",200,0,100);
00106   h_recHits_EEM_time          = dqmStore_->book1D("recHits_EEM_time","recHits_EEM_time",200,-50,50);
00107 
00108   // Basic Clusters ----------------------------------------------    
00109   // ... associated endcap rec hits
00110   h_basicClusters_recHits_EE_recoFlag = dqmStore_->book1D("basicClusters_recHits_EE_recoFlag","basicClusters_recHits_EE_recoFlag",16,-0.5,15.5);  
00111 
00112   // Super Clusters ----------------------------------------------
00113   // ... endcap
00114   h_superClusters_EEP_nBC    = dqmStore_->book1D("superClusters_EEP_nBC","superClusters_EEP_nBC",100,0.,100.);
00115   h_superClusters_EEM_nBC    = dqmStore_->book1D("superClusters_EEM_nBC","superClusters_EEM_nBC",100,0.,100.);
00116 
00117   h_superClusters_eta        = dqmStore_->book1D("superClusters_eta","superClusters_eta",150,-3.,3.);
00118   h_superClusters_EE_phi     = dqmStore_->book1D("superClusters_EE_phi","superClusters_EE_phi",360,-3.1415927,3.1415927);
00119   
00120 }
00121 
00122 
00123 
00124 EERecoSummary::~EERecoSummary()
00125 {
00126         // do anything here that needs to be done at desctruction time
00127         // (e.g. close files, deallocate resources etc.)
00128 }
00129 
00130 
00131 //
00132 // member functions
00133 //
00134 
00135 // ------------ method called to for each event  ------------
00136 void EERecoSummary::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
00137 {
00138   
00139   //Get the magnetic field
00140   edm::ESHandle<MagneticField> theMagField;
00141   iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
00142 
00143   // --- REDUCED REC HITS ------------------------------------------------------------------------------------- 
00144   // ... endcap
00145   edm::Handle<EcalRecHitCollection> redRecHitsEE;
00146   ev.getByLabel( redRecHitCollection_EE_, redRecHitsEE );
00147   const EcalRecHitCollection* theEndcapEcalredRecHits = redRecHitsEE.product () ;
00148   if ( ! redRecHitsEE.isValid() ) {
00149     edm::LogWarning("EERecoSummary") << "redRecHitsEE not found"; 
00150   }
00151   
00152   for ( EcalRecHitCollection::const_iterator itr = theEndcapEcalredRecHits->begin () ;
00153         itr != theEndcapEcalredRecHits->end () ; ++itr)
00154   {
00155       
00156     EEDetId eeid( itr -> id() );
00157 
00158     h_redRecHits_EE_recoFlag->Fill( itr -> recoFlag() );
00159 
00160   }
00161 
00162   // --- REC HITS ------------------------------------------------------------------------------------- 
00163   
00164   // ... endcap
00165   edm::Handle<EcalRecHitCollection> recHitsEE;
00166   ev.getByLabel( recHitCollection_EE_, recHitsEE );
00167   const EcalRecHitCollection* theEndcapEcalRecHits = recHitsEE.product () ;
00168   if ( ! recHitsEE.isValid() ) {
00169     edm::LogWarning("EERecoSummary") << "recHitsEE not found"; 
00170   }
00171 
00172   float maxRecHitEnergyEEP = -999.;
00173   float maxRecHitEnergyEEM = -999.;
00174 
00175   EEDetId eeid_MrecHitEEM;
00176   EEDetId eeid_MrecHitEEP;
00177 
00178   for ( EcalRecHitCollection::const_iterator itr = theEndcapEcalRecHits->begin () ;
00179         itr != theEndcapEcalRecHits->end () ; ++itr)
00180     {
00181       
00182       EEDetId eeid( itr -> id() );
00183 
00184       // EE+
00185       if ( eeid.zside() > 0 ){
00186 
00187         h_recHits_EE_recoFlag       -> Fill( itr -> recoFlag() );
00188 
00189         // max E rec hit
00190         if (itr -> energy() > maxRecHitEnergyEEP && 
00191             !(eeid.ix()>=41 && eeid.ix()<=60 && eeid.iy()>=41 && eeid.iy()<=60) ) {
00192           maxRecHitEnergyEEP = itr -> energy() ;
00193         }
00194         
00195         // only channels above noise
00196         if (  itr -> energy() > ethrEE_ ){
00197           h_recHits_EEP_Chi2          -> Fill( itr -> chi2() );
00198           h_recHits_EEP_time          -> Fill( itr -> time() );
00199         }
00200       }
00201       
00202       // EE-
00203       if ( eeid.zside() < 0 ){
00204         
00205         h_recHits_EE_recoFlag       -> Fill( itr -> recoFlag() );
00206         
00207         // max E rec hit
00208         if (itr -> energy() > maxRecHitEnergyEEM && 
00209             !(eeid.ix()>=41 && eeid.ix()<=60 && eeid.iy()>=41 && eeid.iy()<=60) ) {
00210           maxRecHitEnergyEEM = itr -> energy() ;
00211         }
00212         
00213         // only channels above noise
00214         if (  itr -> energy() > ethrEE_ ) {
00215           h_recHits_EEM_Chi2          -> Fill( itr -> chi2() );
00216           h_recHits_EEM_time          -> Fill( itr -> time() );
00217         }
00218 
00219       }
00220     } // end loop over EE rec hits
00221 
00222   // energy
00223   h_recHits_EEP_energyMax -> Fill( maxRecHitEnergyEEP );
00224   h_recHits_EEM_energyMax -> Fill( maxRecHitEnergyEEM );
00225 
00226   //--- BASIC CLUSTERS --------------------------------------------------------------
00227 
00228   // ... endcap
00229   edm::Handle<reco::BasicClusterCollection> basicClusters_EE_h;
00230   ev.getByLabel( basicClusterCollection_EE_, basicClusters_EE_h );
00231   if ( ! basicClusters_EE_h.isValid() ) {
00232     edm::LogWarning("EERecoSummary") << "basicClusters_EE_h not found"; 
00233   }
00234 
00235   for (unsigned int icl = 0; icl < basicClusters_EE_h->size(); ++icl) {
00236     
00237     //Get the associated RecHits
00238     const std::vector<std::pair<DetId,float> > & hits= (*basicClusters_EE_h)[icl].hitsAndFractions();
00239     for (std::vector<std::pair<DetId,float> > ::const_iterator rh = hits.begin(); rh!=hits.end(); ++rh){
00240       
00241       EBRecHitCollection::const_iterator itrechit = theEndcapEcalRecHits->find((*rh).first);
00242       if (itrechit==theEndcapEcalRecHits->end()) continue;
00243       h_basicClusters_recHits_EE_recoFlag -> Fill ( itrechit -> recoFlag() );
00244     }
00245   
00246   }
00247 
00248   // Super Clusters
00249   // ... endcap
00250   edm::Handle<reco::SuperClusterCollection> superClusters_EE_h;
00251   ev.getByLabel( superClusterCollection_EE_, superClusters_EE_h );
00252   const reco::SuperClusterCollection* theEndcapSuperClusters = superClusters_EE_h.product () ;
00253   if ( ! superClusters_EE_h.isValid() ) {
00254     edm::LogWarning("EERecoSummary") << "superClusters_EE_h not found"; 
00255   }
00256 
00257   for (reco::SuperClusterCollection::const_iterator itSC = theEndcapSuperClusters->begin(); 
00258        itSC != theEndcapSuperClusters->end(); ++itSC ) {
00259 
00260     double scEt = itSC -> energy() * sin(2.*atan( exp(- itSC->position().eta() )));
00261 
00262     if (scEt < scEtThrEE_ ) continue;
00263 
00264     h_superClusters_eta       -> Fill( itSC -> eta() );
00265     h_superClusters_EE_phi    -> Fill( itSC -> phi() );
00266     
00267     if  ( itSC -> z() > 0 ){
00268       h_superClusters_EEP_nBC    -> Fill( itSC -> clustersSize() );      
00269     }
00270 
00271     if  ( itSC -> z() < 0 ){
00272       h_superClusters_EEM_nBC    -> Fill( itSC -> clustersSize() );      
00273     }
00274   }
00275 
00276   //--------------------------------------------------------
00277  
00278 }
00279 
00280 
00281 // ------------ method called once each job just before starting event loop  ------------
00282 void 
00283 EERecoSummary::beginJob()
00284 {
00285 }
00286 
00287 // ------------ method called once each job just after ending the event loop  ------------
00288 void 
00289 EERecoSummary::endJob() 
00290 {}
00291 
00292 
00293 //define this as a plug-in
00294 DEFINE_FWK_MODULE(EERecoSummary);