Reference Manual

Go to the documentation of this file.
00001 #include "Validation/EcalClusters/interface/EgammaSuperClusters.h"
00002 
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "FWCore/Utilities/interface/Exception.h"
00005 #include "FWCore/ServiceRegistry/interface/Service.h"
00006 
00007 #include "DataFormats/Common/interface/Handle.h"
00008 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00009 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00010 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
00011 
00012 #include "DataFormats/Math/interface/deltaPhi.h"
00013 
00014 #include "DataFormats/GeometryVector/interface/Pi.h"
00015 #include "DQMServices/Core/interface/DQMStore.h"
00016 
00017 EgammaSuperClusters::EgammaSuperClusters( const edm::ParameterSet& ps )
00018 {
00019         outputFile_ = ps.getUntrackedParameter<std::string>("outputFile", "");
00020         //CMSSW_Version_ = ps.getUntrackedParameter<std::string>("CMSSW_Version", "");
00021 
00022         verboseDBE_ = ps.getUntrackedParameter<bool>("verboseDBE", false);
00023 
00024         hist_min_Size_ = ps.getParameter<double>("hist_min_Size");
00025         hist_max_Size_ = ps.getParameter<double>("hist_max_Size");
00026         hist_bins_Size_ = ps.getParameter<int>   ("hist_bins_Size");
00027 
00028         hist_min_NumBC_ = ps.getParameter<double>("hist_min_NumBC");
00029         hist_max_NumBC_ = ps.getParameter<double>("hist_max_NumBC");
00030         hist_bins_NumBC_ = ps.getParameter<int>   ("hist_bins_NumBC");
00031 
00032         hist_min_ET_ = ps.getParameter<double>("hist_min_ET");
00033         hist_max_ET_ = ps.getParameter<double>("hist_max_ET");
00034         hist_bins_ET_ = ps.getParameter<int>   ("hist_bins_ET");
00035 
00036         hist_min_Eta_ = ps.getParameter<double>("hist_min_Eta");
00037         hist_max_Eta_ = ps.getParameter<double>("hist_max_Eta");
00038         hist_bins_Eta_ = ps.getParameter<int>   ("hist_bins_Eta");
00039 
00040         hist_min_Phi_ = ps.getParameter<double>("hist_min_Phi");
00041         hist_max_Phi_ = ps.getParameter<double>("hist_max_Phi");
00042         hist_bins_Phi_ = ps.getParameter<int>   ("hist_bins_Phi");
00043 
00044         hist_min_S1toS9_ = ps.getParameter<double>("hist_min_S1toS9");
00045         hist_max_S1toS9_ = ps.getParameter<double>("hist_max_S1toS9");
00046         hist_bins_S1toS9_ = ps.getParameter<int>   ("hist_bins_S1toS9");
00047 
00048         hist_min_S25toE_ = ps.getParameter<double>("hist_min_S25toE");
00049         hist_max_S25toE_ = ps.getParameter<double>("hist_max_S25toE");
00050         hist_bins_S25toE_ = ps.getParameter<int>   ("hist_bins_S25toE");
00051 
00052         hist_min_EoverTruth_ = ps.getParameter<double>("hist_min_EoverTruth");
00053         hist_max_EoverTruth_ = ps.getParameter<double>("hist_max_EoverTruth");
00054         hist_bins_EoverTruth_ = ps.getParameter<int>   ("hist_bins_EoverTruth");
00055         
00056         hist_min_deltaR_ = ps.getParameter<double>("hist_min_deltaR");
00057         hist_max_deltaR_ = ps.getParameter<double>("hist_max_deltaR");
00058         hist_bins_deltaR_ = ps.getParameter<int>   ("hist_bins_deltaR");
00059 
00060         hist_min_phiWidth_ = ps.getParameter<double>("hist_min_phiWidth");
00061         hist_max_phiWidth_ = ps.getParameter<double>("hist_max_phiWidth");
00062         hist_bins_phiWidth_ = ps.getParameter<int>("hist_bins_phiWidth");
00063 
00064         hist_min_etaWidth_ = ps.getParameter<double>("hist_min_etaWidth");
00065         hist_max_etaWidth_ = ps.getParameter<double>("hist_max_etaWidth");
00066         hist_bins_etaWidth_ = ps.getParameter<int>("hist_bins_etaWidth");
00067 
00068         hist_bins_preshowerE_ = ps.getParameter<int>("hist_bins_preshowerE");
00069         hist_min_preshowerE_ = ps.getParameter<double>("hist_min_preshowerE");
00070         hist_max_preshowerE_ = ps.getParameter<double>("hist_max_preshowerE");
00071 
00072         hist_min_R_ = ps.getParameter<double>("hist_min_R");
00073         hist_max_R_ = ps.getParameter<double>("hist_max_R");
00074         hist_bins_R_ = ps.getParameter<int>   ("hist_bins_R");
00075 
00076         MCTruthCollection_ = ps.getParameter<edm::InputTag>("MCTruthCollection");
00077 
00078         barrelRawSuperClusterCollection_ = ps.getParameter<edm::InputTag>("barrelRawSuperClusterCollection");
00079         barrelCorSuperClusterCollection_ = ps.getParameter<edm::InputTag>("barrelCorSuperClusterCollection");
00080 
00081         endcapRawSuperClusterCollection_ = ps.getParameter<edm::InputTag>("endcapRawSuperClusterCollection");
00082         endcapPreSuperClusterCollection_ = ps.getParameter<edm::InputTag>("endcapPreSuperClusterCollection");
00083         endcapCorSuperClusterCollection_ = ps.getParameter<edm::InputTag>("endcapCorSuperClusterCollection");
00084 
00085         barrelRecHitCollection_ = ps.getParameter<edm::InputTag>("barrelRecHitCollection");
00086         endcapRecHitCollection_ = ps.getParameter<edm::InputTag>("endcapRecHitCollection");
00087 }
00088 
00089 EgammaSuperClusters::~EgammaSuperClusters() {}
00090 
00091 void EgammaSuperClusters::beginJob() 
00092 {
00093         dbe_ = edm::Service<DQMStore>().operator->();                   
00094 
00095         if ( verboseDBE_ )
00096         {
00097         dbe_->setVerbose(1);
00098                 dbe_->showDirStructure();
00099         }
00100         else 
00101                 dbe_->setVerbose(0);
00102 
00103         //dbe_->setCurrentFolder("Ecal/CMSSW_"+CMSSW_Version_+"/EcalClusters/SuperClusters/");
00104         dbe_->setCurrentFolder("EcalClusterV/EcalSuperClusters/");
00105 
00106         // Number of SuperClusters
00107         //
00108         hist_EB_RawSC_Size_ 
00109                 = dbe_->book1D("hist_EB_RawSC_Size_","# Raw SuperClusters in Barrel",
00110                         hist_bins_Size_,hist_min_Size_,hist_max_Size_);
00111         hist_EE_RawSC_Size_ 
00112                 = dbe_->book1D("hist_EE_RawSC_Size_","# Raw SuperClusters in Endcap",
00113                         hist_bins_Size_,hist_min_Size_,hist_max_Size_);
00114         hist_EB_CorSC_Size_
00115                 = dbe_->book1D("hist_EB_CorSC_Size_","# Corrected SuperClusters in Barrel",
00116                         hist_bins_Size_,hist_min_Size_,hist_max_Size_);
00117         hist_EE_CorSC_Size_
00118                 = dbe_->book1D("hist_EE_CorSC_Size_","# Corrected SuperClusters in Endcap",
00119                         hist_bins_Size_,hist_min_Size_,hist_max_Size_);
00120         hist_EE_PreSC_Size_
00121                 = dbe_->book1D("hist_EE_PreSC_Size_","# SuperClusters with Preshower in Endcap",
00122                         hist_bins_Size_,hist_min_Size_,hist_max_Size_);
00123 
00124         // Number of BasicClusters in SuperCluster
00125         //
00126         hist_EB_RawSC_NumBC_ 
00127                 = dbe_->book1D("hist_EB_RawSC_NumBC_","# of Basic Clusters in Raw Super Clusters in Barrel",
00128                         hist_bins_NumBC_,hist_min_NumBC_,hist_max_NumBC_);
00129         hist_EE_RawSC_NumBC_ 
00130                 = dbe_->book1D("hist_EE_RawSC_NumBC_","# of Basic Clusters in Raw Super Clusters in Endcap",
00131                 hist_bins_NumBC_,hist_min_NumBC_,hist_max_NumBC_);
00132         hist_EB_CorSC_NumBC_
00133                 = dbe_->book1D("hist_EB_CorSC_NumBC_","# of Basic Clusters in Corrected SuperClusters in Barrel",
00134                         hist_bins_NumBC_,hist_min_NumBC_,hist_max_NumBC_);
00135         hist_EE_CorSC_NumBC_
00136                 = dbe_->book1D("hist_EE_CorSC_NumBC_","# of Basic Clusters in Corrected SuperClusters in Endcap",
00137                         hist_bins_NumBC_,hist_min_NumBC_,hist_max_NumBC_);
00138         hist_EE_PreSC_NumBC_
00139                 = dbe_->book1D("hist_EE_PreSC_NumBC_","# of Basic Clusters in SuperClusters with Preshower in Endcap",
00140                         hist_bins_NumBC_,hist_min_NumBC_,hist_max_NumBC_);
00141 
00142         // ET distribution of SuperClusters
00143         //
00144         hist_EB_RawSC_ET_ 
00145                 = dbe_->book1D("hist_EB_RawSC_ET_","ET of Raw SuperClusters in Barrel",
00146                         hist_bins_ET_,hist_min_ET_,hist_max_ET_);
00147         hist_EE_RawSC_ET_ 
00148                 = dbe_->book1D("hist_EE_RawSC_ET_","ET of Raw SuperClusters in Endcap",
00149                         hist_bins_ET_,hist_min_ET_,hist_max_ET_);
00150         hist_EB_CorSC_ET_
00151                 = dbe_->book1D("hist_EB_CorSC_ET_","ET of Corrected SuperClusters in Barrel",
00152                         hist_bins_ET_,hist_min_ET_,hist_max_ET_);
00153         hist_EE_CorSC_ET_
00154                 = dbe_->book1D("hist_EE_CorSC_ET_","ET of Corrected SuperClusters in Endcap",
00155                         hist_bins_ET_,hist_min_ET_,hist_max_ET_);
00156         hist_EE_PreSC_ET_
00157                 = dbe_->book1D("hist_EE_PreSC_ET_","ET of SuperClusters with Preshower in Endcap",
00158                         hist_bins_ET_,hist_min_ET_,hist_max_ET_);
00159 
00160         // Eta distribution of SuperClusters
00161         //
00162         hist_EB_RawSC_Eta_ 
00163                 = dbe_->book1D("hist_EB_RawSC_Eta_","Eta of Raw SuperClusters in Barrel",
00164                         hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_);
00165         hist_EE_RawSC_Eta_ 
00166                 = dbe_->book1D("hist_EE_RawSC_Eta_","Eta of Raw SuperClusters in Endcap",
00167                         hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_);
00168         hist_EB_CorSC_Eta_
00169                 = dbe_->book1D("hist_EB_CorSC_Eta_","Eta of Corrected SuperClusters in Barrel",
00170                         hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_);
00171         hist_EE_CorSC_Eta_
00172                 = dbe_->book1D("hist_EE_CorSC_Eta_","Eta of Corrected SuperClusters in Endcap",
00173                         hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_);
00174         hist_EE_PreSC_Eta_
00175                 = dbe_->book1D("hist_EE_PreSC_Eta_","Eta of SuperClusters with Preshower in Endcap",
00176                         hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_);
00177 
00178         // Phi distribution of SuperClusters
00179         //
00180         hist_EB_RawSC_Phi_
00181                 = dbe_->book1D("hist_EB_RawSC_Phi_","Phi of Raw SuperClusters in Barrel",
00182                         hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_);
00183         hist_EE_RawSC_Phi_
00184                 = dbe_->book1D("hist_EE_RawSC_Phi_","Phi of Raw SuperClusters in Endcap",
00185                         hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_);
00186         hist_EB_CorSC_Phi_
00187                 = dbe_->book1D("hist_EB_CorSC_Phi_","Phi of Corrected SuperClusters in Barrel",
00188                         hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_);
00189         hist_EE_CorSC_Phi_
00190                 = dbe_->book1D("hist_EE_CorSC_Phi_","Phi of Corrected SuperClusters in Endcap",
00191                         hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_);
00192         hist_EE_PreSC_Phi_
00193                 = dbe_->book1D("hist_EE_PreSC_Phi_","Phi of SuperClusters with Preshower in Endcap",
00194                         hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_);
00195 
00196         // S1/S9 distribution of SuperClusters
00197         //
00198         hist_EB_RawSC_S1toS9_ 
00199                 = dbe_->book1D("hist_EB_RawSC_S1toS9_","S1/S9 of Raw Super Clusters in Barrel",
00200                         hist_bins_S1toS9_,hist_min_S1toS9_,hist_max_S1toS9_);
00201         hist_EE_RawSC_S1toS9_ 
00202                 = dbe_->book1D("hist_EE_RawSC_S1toS9_","S1/S9 of Raw Super Clusters in Endcap",
00203                         hist_bins_S1toS9_,hist_min_S1toS9_,hist_max_S1toS9_);
00204         hist_EB_CorSC_S1toS9_
00205                 = dbe_->book1D("hist_EB_CorSC_S1toS9_","S1/S9 of Corrected SuperClusters in Barrel",
00206                         hist_bins_S1toS9_,hist_min_S1toS9_,hist_max_S1toS9_);
00207         hist_EE_CorSC_S1toS9_
00208                 = dbe_->book1D("hist_EE_CorSC_S1toS9_","S1/S9 of Corrected SuperClusters in Endcap",
00209                         hist_bins_S1toS9_,hist_min_S1toS9_,hist_max_S1toS9_);
00210         hist_EE_PreSC_S1toS9_
00211                 = dbe_->book1D("hist_EE_PreSC_S1toS9_","S1/S9 of SuperClusters with Preshower in Endcap",
00212                         hist_bins_S1toS9_,hist_min_S1toS9_,hist_max_S1toS9_);
00213 
00214         // S25/E distribution of SuperClusters
00215         //
00216         hist_EB_RawSC_S25toE_ 
00217                 = dbe_->book1D("hist_EB_RawSC_S25toE_","S25/E of Raw Super Clusters in Barrel",
00218                         hist_bins_S25toE_,hist_min_S25toE_,hist_max_S25toE_);
00219         hist_EE_RawSC_S25toE_ 
00220                 = dbe_->book1D("hist_EE_RawSC_S25toE_","S25/E of Raw Super Clusters in Endcap",
00221                         hist_bins_S25toE_,hist_min_S25toE_,hist_max_S25toE_);
00222         hist_EB_CorSC_S25toE_
00223                 = dbe_->book1D("hist_EB_CorSC_S25toE_","S25/E of Corrected SuperClusters in Barrel",
00224                         hist_bins_S25toE_,hist_min_S25toE_,hist_max_S25toE_);
00225         hist_EE_CorSC_S25toE_
00226                 = dbe_->book1D("hist_EE_CorSC_S25toE_","S25/E of Corrected SuperClusters in Endcap",
00227                         hist_bins_S25toE_,hist_min_S25toE_,hist_max_S25toE_);
00228         hist_EE_PreSC_S25toE_
00229                 = dbe_->book1D("hist_EE_PreSC_S25toE_","S25/E of SuperClusters with Preshower in Endcap",
00230                         hist_bins_S25toE_,hist_min_S25toE_,hist_max_S25toE_);
00231 
00232         // E/E(true) distribution of SuperClusters
00233         //
00234         hist_EB_RawSC_EoverTruth_ 
00235                 = dbe_->book1D("hist_EB_RawSC_EoverTruth_","E/True E of Raw SuperClusters in Barrel",   
00236                         hist_bins_EoverTruth_,hist_min_EoverTruth_,hist_max_EoverTruth_);
00237         hist_EE_RawSC_EoverTruth_ 
00238                 = dbe_->book1D("hist_EE_RawSC_EoverTruth_","E/True E of Raw SuperClusters in Endcap",
00239                         hist_bins_EoverTruth_,hist_min_EoverTruth_,hist_max_EoverTruth_);
00240         hist_EB_CorSC_EoverTruth_
00241                 = dbe_->book1D("hist_EB_CorSC_EoverTruth_","E/True E of Corrected SuperClusters in Barrel",
00242                         hist_bins_EoverTruth_,hist_min_EoverTruth_,hist_max_EoverTruth_);
00243         hist_EE_CorSC_EoverTruth_
00244                 = dbe_->book1D("hist_EE_CorSC_EoverTruth_","E/True E of Corrected SuperClusters in Endcap",
00245                         hist_bins_EoverTruth_,hist_min_EoverTruth_,hist_max_EoverTruth_);
00246         hist_EE_PreSC_EoverTruth_
00247                 = dbe_->book1D("hist_EE_PreSC_EoverTruth_","E/True E of SuperClusters with Preshower in Endcap",
00248                         hist_bins_EoverTruth_,hist_min_EoverTruth_,hist_max_EoverTruth_);
00249 
00250         // dR distribution of SuperClusters from truth
00251         //
00252         hist_EB_RawSC_deltaR_ 
00253                 = dbe_->book1D("hist_EB_RawSC_deltaR_","dR to MC truth of Raw Super Clusters in Barrel",
00254                         hist_bins_deltaR_,hist_min_deltaR_,hist_max_deltaR_);
00255         hist_EE_RawSC_deltaR_ 
00256                 = dbe_->book1D("hist_EE_RawSC_deltaR_","dR to MC truth of Raw Super Clusters in Endcap",
00257                         hist_bins_deltaR_,hist_min_deltaR_,hist_max_deltaR_);
00258         hist_EB_CorSC_deltaR_
00259                 = dbe_->book1D("hist_EB_CorSC_deltaR_","dR to MC truth of Corrected SuperClusters in Barrel",
00260                         hist_bins_deltaR_,hist_min_deltaR_,hist_max_deltaR_);
00261         hist_EE_CorSC_deltaR_
00262                 = dbe_->book1D("hist_EE_CorSC_deltaR_","dR to MC truth of Corrected SuperClusters in Endcap",
00263                         hist_bins_deltaR_,hist_min_deltaR_,hist_max_deltaR_);
00264         hist_EE_PreSC_deltaR_
00265                 = dbe_->book1D("hist_EE_PreSC_deltaR_","dR to MC truth of SuperClusters with Preshower in Endcap",
00266                         hist_bins_deltaR_,hist_min_deltaR_,hist_max_deltaR_);
00267 
00268         // phi width stored in corrected SuperClusters
00269         hist_EB_CorSC_phiWidth_
00270                 = dbe_->book1D("hist_EB_CorSC_phiWidth_","phiWidth of Corrected Super Clusters in Barrel",
00271                         hist_bins_phiWidth_,hist_min_phiWidth_,hist_max_phiWidth_);
00272         hist_EE_CorSC_phiWidth_
00273                 = dbe_->book1D("hist_EE_CorSC_phiWidth_","phiWidth of Corrected Super Clusters in Endcap",
00274                         hist_bins_phiWidth_,hist_min_phiWidth_,hist_max_phiWidth_);
00275 
00276         // eta width stored in corrected SuperClusters
00277         hist_EB_CorSC_etaWidth_
00278                 = dbe_->book1D("hist_EB_CorSC_etaWidth_","etaWidth of Corrected Super Clusters in Barrel",
00279                         hist_bins_etaWidth_,hist_min_etaWidth_,hist_max_etaWidth_);
00280         hist_EE_CorSC_etaWidth_
00281                 = dbe_->book1D("hist_EE_CorSC_etaWidth_","etaWidth of Corrected Super Clusters in Endcap",
00282                         hist_bins_etaWidth_,hist_min_etaWidth_,hist_max_etaWidth_);
00283 
00284 
00285         // preshower energy
00286         hist_EE_PreSC_preshowerE_
00287                 = dbe_->book1D("hist_EE_PreSC_preshowerE_","preshower energy in Super Clusters with Preshower in Endcap",
00288                         hist_bins_preshowerE_,hist_min_preshowerE_,hist_max_preshowerE_);
00289         hist_EE_CorSC_preshowerE_
00290                 = dbe_->book1D("hist_EE_CorSC_preshowerE_","preshower energy in Corrected Super Clusters with Preshower in Endcap",
00291                         hist_bins_preshowerE_,hist_min_preshowerE_,hist_max_preshowerE_);
00292 
00293 
00294         //
00295         hist_EB_CorSC_ET_vs_Eta_ = dbe_->book2D( "hist_EB_CorSC_ET_vs_Eta_", "Corr Super Cluster ET versus Eta in Barrel", 
00296                                               hist_bins_ET_, hist_min_ET_, hist_max_ET_,
00297                                               hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_ );
00298 
00299         hist_EB_CorSC_ET_vs_Phi_ = dbe_->book2D( "hist_EB_CorSC_ET_vs_Phi_", "Corr Super Cluster ET versus Phi in Barrel", 
00300                                               hist_bins_ET_, hist_min_ET_, hist_max_ET_,
00301                                               hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_ );
00302 
00303         hist_EE_CorSC_ET_vs_Eta_ = dbe_->book2D( "hist_EE_CorSC_ET_vs_Eta_", "Corr Super Cluster ET versus Eta in Endcap", 
00304                                               hist_bins_ET_, hist_min_ET_, hist_max_ET_,
00305                                               hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_ );
00306 
00307         hist_EE_CorSC_ET_vs_Phi_ = dbe_->book2D( "hist_EE_CorSC_ET_vs_Phi_", "Corr Super Cluster ET versus Phi in Endcap", 
00308                                               hist_bins_ET_, hist_min_ET_, hist_max_ET_,
00309                                               hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_ );
00310 
00311         hist_EE_CorSC_ET_vs_R_ = dbe_->book2D( "hist_EE_CorSC_ET_vs_R_", "Corr Super Cluster ET versus Radius in Endcap", 
00312                                               hist_bins_ET_, hist_min_ET_, hist_max_ET_,
00313                                               hist_bins_R_,hist_min_R_,hist_max_R_ );
00314 
00315 
00316 }
00317 
00318 void EgammaSuperClusters::analyze( const edm::Event& evt, const edm::EventSetup& es )
00319 {
00320 
00321         bool skipMC = false;
00322         bool skipBarrel = false;
00323         bool skipEndcap = false;
00324 
00325         //
00326         // Get MCTRUTH
00327         //
00328         edm::Handle<edm::HepMCProduct> pMCTruth ;
00329         evt.getByLabel(MCTruthCollection_, pMCTruth);
00330         if (!pMCTruth.isValid()) {
00331           edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label "
00332                                                << MCTruthCollection_.label();
00333           skipMC = true;
00334         }
00335         const HepMC::GenEvent* genEvent = pMCTruth->GetEvent();
00336 
00337         if( skipMC ) return;
00338 
00339         //
00340         // Get the BARREL products 
00341         //
00342         edm::Handle<reco::SuperClusterCollection> pBarrelRawSuperClusters;
00343         evt.getByLabel(barrelRawSuperClusterCollection_, pBarrelRawSuperClusters);
00344         if (!pBarrelRawSuperClusters.isValid()) {
00345           edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label " 
00346                                                << barrelRawSuperClusterCollection_.label();
00347           skipBarrel = true;
00348         }
00349 
00350         edm::Handle<reco::SuperClusterCollection> pBarrelCorSuperClusters;
00351         evt.getByLabel(barrelCorSuperClusterCollection_, pBarrelCorSuperClusters);
00352         if (!pBarrelCorSuperClusters.isValid()) {
00353           edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label "
00354                                                << barrelCorSuperClusterCollection_.label();
00355           skipBarrel = true;
00356         }
00357 
00358         edm::Handle< EBRecHitCollection > pBarrelRecHitCollection;
00359         evt.getByLabel( barrelRecHitCollection_, pBarrelRecHitCollection );
00360         if ( ! pBarrelRecHitCollection.isValid() ) {
00361           skipBarrel = true;
00362         }
00363         edm::Handle< EERecHitCollection > pEndcapRecHitCollection;
00364         evt.getByLabel( endcapRecHitCollection_, pEndcapRecHitCollection );
00365         if ( ! pEndcapRecHitCollection.isValid() ) {
00366           skipEndcap = true;
00367         }
00368 
00369         if( skipBarrel || skipEndcap ) return;
00370 
00371         EcalClusterLazyTools lazyTool( evt, es, barrelRecHitCollection_, endcapRecHitCollection_ );
00372 
00373         // Get the BARREL collections        
00374         const reco::SuperClusterCollection* barrelRawSuperClusters = pBarrelRawSuperClusters.product();
00375         const reco::SuperClusterCollection* barrelCorSuperClusters = pBarrelCorSuperClusters.product();
00376 
00377         // Number of entries in collections
00378         hist_EB_RawSC_Size_->Fill(barrelRawSuperClusters->size());
00379         hist_EB_CorSC_Size_->Fill(barrelCorSuperClusters->size());
00380 
00381         // Do RAW BARREL SuperClusters
00382         for(reco::SuperClusterCollection::const_iterator aClus = barrelRawSuperClusters->begin(); 
00383                 aClus != barrelRawSuperClusters->end(); aClus++)
00384         {
00385                 // kinematics
00386                 hist_EB_RawSC_NumBC_->Fill(aClus->clustersSize());
00387                 hist_EB_RawSC_ET_->Fill(aClus->energy()/std::cosh(aClus->position().eta()));
00388                 hist_EB_RawSC_Eta_->Fill(aClus->position().eta());
00389                 hist_EB_RawSC_Phi_->Fill(aClus->position().phi());
00390 
00391                 // cluster shape
00392                 const reco::CaloClusterPtr seed = aClus->seed();
00393                 hist_EB_RawSC_S1toS9_->Fill( lazyTool.eMax( *seed ) / lazyTool.e3x3( *seed ) );
00394                 hist_EB_RawSC_S25toE_->Fill( lazyTool.e5x5( *seed ) / aClus->energy() );
00395 
00396                 // truth
00397                 double dRClosest = 999.9;
00398                 double energyClosest = 0;
00399                 closestMCParticle(genEvent, *aClus, dRClosest, energyClosest);
00400 
00401                 if (dRClosest < 0.1)
00402                 {
00403                         hist_EB_RawSC_EoverTruth_->Fill(aClus->energy()/energyClosest);         
00404                         hist_EB_RawSC_deltaR_->Fill(dRClosest);
00405 
00406                 }
00407 
00408         }
00409 
00410         // Do CORRECTED BARREL SuperClusters
00411         for(reco::SuperClusterCollection::const_iterator aClus = barrelCorSuperClusters->begin();
00412                 aClus != barrelCorSuperClusters->end(); aClus++)
00413         {
00414                 // kinematics
00415                 hist_EB_CorSC_NumBC_->Fill(aClus->clustersSize());
00416                 hist_EB_CorSC_ET_->Fill(aClus->energy()/std::cosh(aClus->position().eta()));
00417                 hist_EB_CorSC_Eta_->Fill(aClus->position().eta());
00418                 hist_EB_CorSC_Phi_->Fill(aClus->position().phi());
00419 
00420                 hist_EB_CorSC_ET_vs_Eta_->Fill( aClus->energy()/std::cosh(aClus->position().eta()), aClus->eta() );
00421                 hist_EB_CorSC_ET_vs_Phi_->Fill( aClus->energy()/std::cosh(aClus->position().eta()), aClus->phi() );
00422 
00423 
00424                 // cluster shape
00425                 const reco::CaloClusterPtr seed = aClus->seed();
00426                 hist_EB_CorSC_S1toS9_->Fill( lazyTool.eMax( *seed ) / lazyTool.e3x3( *seed ) );
00427                 hist_EB_CorSC_S25toE_->Fill( lazyTool.e5x5( *seed ) / aClus->energy() );
00428 
00429                 // correction variables
00430                 hist_EB_CorSC_phiWidth_->Fill(aClus->phiWidth());
00431                 hist_EB_CorSC_etaWidth_->Fill(aClus->etaWidth());
00432 
00433                 // truth
00434                 double dRClosest = 999.9;
00435                 double energyClosest = 0;
00436                 closestMCParticle(genEvent, *aClus, dRClosest, energyClosest);
00437 
00438                 if (dRClosest < 0.1)
00439                 {
00440                         hist_EB_CorSC_EoverTruth_->Fill(aClus->energy()/energyClosest);
00441                         hist_EB_CorSC_deltaR_->Fill(dRClosest);
00442 
00443                 }
00444 
00445         }
00446 
00447         //
00448         // Get the ENDCAP products
00449         //
00450         edm::Handle<reco::SuperClusterCollection> pEndcapRawSuperClusters;
00451         evt.getByLabel(endcapRawSuperClusterCollection_, pEndcapRawSuperClusters);
00452         if (!pEndcapRawSuperClusters.isValid()) {
00453           edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label " 
00454                                                << endcapRawSuperClusterCollection_.label();
00455         }
00456 
00457         edm::Handle<reco::SuperClusterCollection> pEndcapPreSuperClusters;
00458         evt.getByLabel(endcapPreSuperClusterCollection_, pEndcapPreSuperClusters);
00459         if (!pEndcapPreSuperClusters.isValid()) {
00460           edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label "
00461                                                << endcapPreSuperClusterCollection_.label();
00462         }
00463 
00464         edm::Handle<reco::SuperClusterCollection> pEndcapCorSuperClusters;
00465         evt.getByLabel(endcapCorSuperClusterCollection_, pEndcapCorSuperClusters);
00466         if (!pEndcapCorSuperClusters.isValid()) {
00467           edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label "
00468                                                << endcapCorSuperClusterCollection_.label();
00469         }
00470 
00471         // Get the ENDCAP collections
00472         const reco::SuperClusterCollection* endcapRawSuperClusters = pEndcapRawSuperClusters.product();
00473         const reco::SuperClusterCollection* endcapPreSuperClusters = pEndcapPreSuperClusters.product();
00474         const reco::SuperClusterCollection* endcapCorSuperClusters = pEndcapCorSuperClusters.product();
00475 
00476         // Number of entries in collections
00477         hist_EE_RawSC_Size_->Fill(endcapRawSuperClusters->size());
00478         hist_EE_PreSC_Size_->Fill(endcapPreSuperClusters->size());
00479         hist_EE_CorSC_Size_->Fill(endcapCorSuperClusters->size());
00480 
00481         // Do RAW ENDCAP SuperClusters
00482         for(reco::SuperClusterCollection::const_iterator aClus = endcapRawSuperClusters->begin(); 
00483                 aClus != endcapRawSuperClusters->end(); aClus++)
00484         {
00485                 hist_EE_RawSC_NumBC_->Fill(aClus->clustersSize());
00486                 hist_EE_RawSC_ET_->Fill(aClus->energy()/std::cosh(aClus->position().eta()));
00487                 hist_EE_RawSC_Eta_->Fill(aClus->position().eta());
00488                 hist_EE_RawSC_Phi_->Fill(aClus->position().phi());
00489 
00490                 const reco::CaloClusterPtr seed = aClus->seed();
00491                 hist_EE_RawSC_S1toS9_->Fill( lazyTool.eMax( *seed ) / lazyTool.e3x3( *seed ) );
00492                 hist_EE_RawSC_S25toE_->Fill( lazyTool.e5x5( *seed ) / aClus->energy() );
00493 
00494                 // truth
00495                 double dRClosest = 999.9;
00496                 double energyClosest = 0;
00497                 closestMCParticle(genEvent, *aClus, dRClosest, energyClosest);
00498 
00499                 if (dRClosest < 0.1)
00500                 {
00501                         hist_EE_RawSC_EoverTruth_->Fill(aClus->energy()/energyClosest);
00502                         hist_EE_RawSC_deltaR_->Fill(dRClosest);
00503                 }
00504 
00505         }
00506 
00507         // Do ENDCAP SuperClusters with PRESHOWER
00508         for(reco::SuperClusterCollection::const_iterator aClus = endcapPreSuperClusters->begin();
00509                 aClus != endcapPreSuperClusters->end(); aClus++)
00510         {
00511                 hist_EE_PreSC_NumBC_->Fill(aClus->clustersSize());
00512                 hist_EE_PreSC_ET_->Fill(aClus->energy()/std::cosh(aClus->position().eta()));
00513                 hist_EE_PreSC_Eta_->Fill(aClus->position().eta());
00514                 hist_EE_PreSC_Phi_->Fill(aClus->position().phi());
00515                 hist_EE_PreSC_preshowerE_->Fill(aClus->preshowerEnergy());
00516 
00517                 const reco::CaloClusterPtr seed = aClus->seed();
00518                 hist_EE_PreSC_S1toS9_->Fill( lazyTool.eMax( *seed ) / lazyTool.e3x3( *seed ) );
00519                 hist_EE_PreSC_S25toE_->Fill( lazyTool.e5x5( *seed ) / aClus->energy() );
00520 
00521                 // truth
00522                 double dRClosest = 999.9;
00523                 double energyClosest = 0;
00524                 closestMCParticle(genEvent, *aClus, dRClosest, energyClosest);
00525 
00526                 if (dRClosest < 0.1)
00527                 {
00528                         hist_EE_PreSC_EoverTruth_->Fill(aClus->energy()/energyClosest);
00529                         hist_EE_PreSC_deltaR_->Fill(dRClosest);
00530                 }
00531 
00532         }
00533 
00534         // Do CORRECTED ENDCAP SuperClusters
00535         for(reco::SuperClusterCollection::const_iterator aClus = endcapCorSuperClusters->begin();
00536                 aClus != endcapCorSuperClusters->end(); aClus++)
00537         {
00538                 hist_EE_CorSC_NumBC_->Fill(aClus->clustersSize());
00539                 hist_EE_CorSC_ET_->Fill(aClus->energy()/std::cosh(aClus->position().eta()));
00540                 hist_EE_CorSC_Eta_->Fill(aClus->position().eta());
00541                 hist_EE_CorSC_Phi_->Fill(aClus->position().phi());
00542                 hist_EE_CorSC_preshowerE_->Fill(aClus->preshowerEnergy());
00543 
00544                 hist_EE_CorSC_ET_vs_Eta_->Fill( aClus->energy()/std::cosh(aClus->position().eta()), aClus->eta() );
00545                 hist_EE_CorSC_ET_vs_Phi_->Fill( aClus->energy()/std::cosh(aClus->position().eta()), aClus->phi() );
00546                 hist_EE_CorSC_ET_vs_R_->Fill( aClus->energy()/std::cosh(aClus->position().eta()), 
00547                                               std::sqrt( std::pow(aClus->x(),2) + std::pow(aClus->y(),2) ) );
00548 
00549 
00550                 // correction variables
00551                 hist_EE_CorSC_phiWidth_->Fill(aClus->phiWidth());
00552                 hist_EE_CorSC_etaWidth_->Fill(aClus->etaWidth());
00553 
00554                 const reco::CaloClusterPtr seed = aClus->seed();
00555                 hist_EE_CorSC_S1toS9_->Fill( lazyTool.eMax( *seed ) / lazyTool.e3x3( *seed ) );
00556                 hist_EE_CorSC_S25toE_->Fill( lazyTool.e5x5( *seed ) / aClus->energy() );
00557 
00558                 // truth
00559                 double dRClosest = 999.9;
00560                 double energyClosest = 0;
00561                 closestMCParticle(genEvent, *aClus, dRClosest, energyClosest);
00562 
00563                 if (dRClosest < 0.1)
00564                 {
00565                         hist_EE_CorSC_EoverTruth_->Fill(aClus->energy()/energyClosest);
00566                         hist_EE_CorSC_deltaR_->Fill(dRClosest);
00567                 }
00568 
00569         }
00570 
00571 }
00572 
00573 void EgammaSuperClusters::endJob()
00574 {
00575         if (outputFile_.size() != 0 && dbe_) dbe_->save(outputFile_);
00576 }
00577 
00578 //
00579 // Closest MC Particle
00580 //
00581 void EgammaSuperClusters::closestMCParticle(const HepMC::GenEvent *genEvent, const reco::SuperCluster &sc, 
00582                                                 double &dRClosest, double &energyClosest)
00583 {
00584 
00585         // SuperCluster eta, phi
00586         double scEta = sc.eta();
00587         double scPhi = sc.phi();
00588 
00589         // initialize dRClosest to a large number
00590         dRClosest = 999.9;
00591 
00592         // loop over the MC truth particles to find the
00593         // closest to the superCluster in dR space
00594         for(HepMC::GenEvent::particle_const_iterator currentParticle = genEvent->particles_begin();
00595                 currentParticle != genEvent->particles_end(); currentParticle++ )
00596         {
00597                 if((*currentParticle)->status() == 1)
00598                 {
00599                         // need GenParticle in ECAL co-ordinates
00600                         HepMC::FourVector vtx = (*currentParticle)->production_vertex()->position();
00601                         double phiTrue = (*currentParticle)->momentum().phi();
00602                         double etaTrue = ecalEta((*currentParticle)->momentum().eta(), vtx.z()/10., vtx.perp()/10.);
00603 
00604                         double dPhi = reco::deltaPhi(phiTrue, scPhi);
00605                         double dEta = scEta - etaTrue;
00606                         double deltaR = std::sqrt(dPhi*dPhi + dEta*dEta);
00607 
00608                         if(deltaR < dRClosest)
00609                         {
00610                                 dRClosest = deltaR;
00611                                 energyClosest = (*currentParticle)->momentum().e();
00612                         }
00613 
00614                 } // end if stable particle     
00615 
00616         } // end loop on get particles
00617 
00618 }
00619 
00620 
00621 //
00622 // Compute Eta in the ECAL co-ordinate system
00623 //
00624 float EgammaSuperClusters::ecalEta(float EtaParticle , float Zvertex, float plane_Radius)
00625 {  
00626         const float R_ECAL           = 136.5;
00627         const float Z_Endcap         = 328.0;
00628         const float etaBarrelEndcap  = 1.479;
00629 
00630         if(EtaParticle != 0.)
00631         {
00632                 float Theta = 0.0  ;
00633                 float ZEcal = (R_ECAL-plane_Radius)*sinh(EtaParticle)+Zvertex;
00634 
00635                 if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
00636                 if(Theta<0.0) Theta = Theta+Geom::pi() ;
00637 
00638                 float ETA = - log(tan(0.5*Theta));
00639 
00640                 if( fabs(ETA) > etaBarrelEndcap )
00641                 {
00642                         float Zend = Z_Endcap ;
00643                         if(EtaParticle<0.0 )  Zend = -Zend ;
00644                         float Zlen = Zend - Zvertex ;
00645                         float RR = Zlen/sinh(EtaParticle);
00646                         Theta = atan((RR+plane_Radius)/Zend);
00647                         if(Theta<0.0) Theta = Theta+Geom::pi() ;
00648                         ETA = - log(tan(0.5*Theta));
00649                 }
00650                 
00651                 return ETA;
00652         }
00653         else
00654         {
00655                 edm::LogWarning("")  << "[EgammaSuperClusters::ecalEta] Warning: Eta equals to zero, not correcting" ;
00656                 return EtaParticle;
00657         }
00658 }
00659 
00660
CMSSW_4_4_3_patch1/src/Validation/EcalClusters/src/EgammaSuperClusters.cc
