Reference Manual

Go to the documentation of this file.
00001 /*   A macro for making a histogram of Jet Pt with cuts
00002 This is a basic way to cut out jets of a certain Pt and Eta using an if statement
00003 This example creates a histogram of Jet Pt, using Jets with Pt above 30 and ETA above -2.1 and below 2.1
00004 */
00005 
00006 #include "TFile.h"
00007 #include "TH1.h"
00008 #include "TH2.h"
00009 #include "TCanvas.h"
00010 #include "TLegend.h"
00011 #include "TSystem.h"
00012 #include "TLorentzVector.h"
00013 
00014 #include "PhysicsTools/SelectorUtils/interface/strbitset.h"
00015 #include "DataFormats/PatCandidates/interface/Jet.h"
00016 #include "DataFormats/PatCandidates/interface/MET.h"
00017 #include "DataFormats/PatCandidates/interface/Photon.h"
00018 #include "DataFormats/PatCandidates/interface/Muon.h"
00019 #include "DataFormats/PatCandidates/interface/Electron.h"
00020 #include "DataFormats/PatCandidates/interface/Tau.h"
00021 #include "PhysicsTools/SelectorUtils/interface/Selector.h"
00022 #include "PhysicsTools/SelectorUtils/interface/JetIDSelectionFunctor.h"
00023 #include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
00024 #include "PhysicsTools/SelectorUtils/interface/RunLumiSelector.h"
00025 #include "FWCore/FWLite/interface/AutoLibraryLoader.h"
00026 #include "PhysicsTools/FWLite/interface/TFileService.h"
00027 #include "DataFormats/FWLite/interface/ChainEvent.h"
00028 #include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
00029 #include "FWCore/ParameterSet/interface/ProcessDesc.h"
00030 
00031 
00032 #include <iostream>
00033 #include <cmath>      //necessary for absolute function fabs()
00034 
00035 using namespace std;
00036 
00040 
00041 // -*- C++ -*-
00042 
00043 // CMS includes
00044 #include "FWCore/Utilities/interface/InputTag.h"
00045 #include "DataFormats/Common/interface/Handle.h"
00046 #include "DataFormats/Common/interface/Ptr.h"
00047 #include "DataFormats/PatCandidates/interface/Jet.h"
00048 #include "PhysicsTools/SelectorUtils/interface/EventSelector.h"
00049 
00050 
00051 // Root includes
00052 #include "TROOT.h"
00053 
00054 using namespace std;
00055 
00056 
00057 class JetIDStudiesSelector : public EventSelector {
00058 public:
00059   JetIDStudiesSelector( edm::ParameterSet const & caloJetIdParams,
00060                         edm::ParameterSet const & pfJetIdParams,
00061                         edm::ParameterSet const & params ) :
00062     jetSel_    (new JetIDSelectionFunctor(caloJetIdParams)), 
00063     pfJetSel_  (new PFJetIDSelectionFunctor(pfJetIdParams)),
00064     jetSrc_    (params.getParameter<edm::InputTag>("jetSrc")), 
00065     pfJetSrc_  (params.getParameter<edm::InputTag>("pfJetSrc")) 
00066   {
00067     bool useCalo = params.getParameter<bool>("useCalo");
00068                
00069     push_back("Calo Cuts");
00070     push_back("Calo Kin Cuts");
00071     push_back("Calo Delta Phi");
00072     push_back("Calo Jet ID");
00073     push_back("PF Cuts");
00074     push_back("PF Kin Cuts");
00075     push_back("PF Delta Phi");
00076     push_back("PF Jet ID");
00077 
00078     set("Calo Cuts", useCalo);
00079     set("Calo Kin Cuts", useCalo);
00080     set("Calo Delta Phi", useCalo);
00081     set("Calo Jet ID", useCalo);
00082     set("PF Cuts", !useCalo);
00083     set("PF Kin Cuts", !useCalo);
00084     set("PF Delta Phi", !useCalo);    
00085     set("PF Jet ID", !useCalo);
00086 
00087     // Indices for fast caching
00088     caloCuts_      = index_type(&bits_, std::string("Calo Cuts") );
00089     caloKin_       = index_type(&bits_, std::string("Calo Kin Cuts"));
00090     caloDeltaPhi_  = index_type(&bits_, std::string("Calo Delta Phi"));
00091     caloJetID_     = index_type(&bits_, std::string("Calo Jet ID"));
00092     pfCuts_        = index_type(&bits_, std::string("PF Cuts"));
00093     pfKin_         = index_type(&bits_, std::string("PF Kin Cuts"));
00094     pfDeltaPhi_    = index_type(&bits_, std::string("PF Delta Phi"));
00095     pfJetID_       = index_type(&bits_, std::string("PF Jet ID"));
00096     
00097   }
00098 
00099   virtual ~JetIDStudiesSelector() {}
00100 
00101   virtual bool operator()( edm::EventBase const & event, pat::strbitset & ret){
00102 
00103     pat::strbitset retCaloJet = jetSel_->getBitTemplate();
00104     pat::strbitset retPFJet = pfJetSel_->getBitTemplate();
00105 
00106 
00107     if ( considerCut(caloCuts_) ) {
00108       passCut(ret, caloCuts_);
00109       event.getByLabel( jetSrc_, h_jets_ );
00110       // Calo Cuts
00111       if ( h_jets_->size() >= 2 || ignoreCut(caloKin_) ) {
00112         passCut(ret, caloKin_);
00113         pat::Jet const & jet0 = h_jets_->at(0);
00114         pat::Jet const & jet1 = h_jets_->at(1);
00115         double dphi = fabs(deltaPhi<double>( jet0.phi(),
00116                                              jet1.phi() ) );
00117           
00118         if ( fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(caloDeltaPhi_) ) {
00119           passCut(ret, caloDeltaPhi_);
00120 
00121             
00122           retCaloJet.set(false);
00123           bool pass0 = (*jetSel_)( jet0, retCaloJet );
00124           retCaloJet.set(false);
00125           bool pass1 = (*jetSel_)( jet1, retCaloJet );
00126           if ( (pass0 && pass1) || ignoreCut(caloJetID_) ) {
00127             passCut(ret, caloJetID_);
00128             caloJet0_ = edm::Ptr<pat::Jet>( h_jets_, 0);
00129             caloJet1_ = edm::Ptr<pat::Jet>( h_jets_, 1);
00130 
00131             return (bool)ret;
00132           }// end if found 2 "loose" jet ID jets
00133         }// end if delta phi
00134       }// end calo kin cuts
00135     }// end if calo cuts
00136 
00137 
00138     if ( considerCut(pfCuts_) ) {
00139 
00140       passCut(ret, pfCuts_);
00141       event.getByLabel( pfJetSrc_, h_pfjets_ );
00142       // PF Cuts
00143       if ( h_pfjets_->size() >= 2 || ignoreCut(pfKin_) ) {
00144         passCut( ret, pfKin_);
00145         pat::Jet const & jet0 = h_pfjets_->at(0);
00146         pat::Jet const & jet1 = h_pfjets_->at(1);
00147         double dphi = fabs(deltaPhi<double>( jet0.phi(),
00148                                              jet1.phi() ) );
00149           
00150         if ( fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(pfDeltaPhi_) ) {
00151           passCut(ret, pfDeltaPhi_);
00152 
00153             
00154           retPFJet.set(false);
00155           bool pass0 = (*pfJetSel_)( jet0, retPFJet );
00156           retPFJet.set(false);
00157           bool pass1 = (*pfJetSel_)( jet1, retPFJet );
00158           if ( (pass0 && pass1) || ignoreCut(pfJetID_) ) {
00159             passCut(ret, pfJetID_);
00160             pfJet0_ = edm::Ptr<pat::Jet>( h_pfjets_, 0);
00161             pfJet1_ = edm::Ptr<pat::Jet>( h_pfjets_, 1);
00162 
00163             return (bool)ret;
00164           }// end if found 2 "loose" jet ID jets
00165         }// end if delta phi
00166       }// end pf kin cuts
00167     }// end if pf cuts
00168 
00169 
00170     setIgnored(ret);
00171 
00172     return false;
00173   }// end of method
00174 
00175   boost::shared_ptr<JetIDSelectionFunctor> const &   jetSel()     const { return jetSel_;}
00176   boost::shared_ptr<PFJetIDSelectionFunctor> const & pfJetSel()   const { return pfJetSel_;}
00177 
00178   vector<pat::Jet>            const &   allCaloJets () const { return *h_jets_; }
00179   vector<pat::Jet>            const &   allPFJets   () const { return *h_pfjets_; }
00180 
00181   pat::Jet                    const &   caloJet0() const { return *caloJet0_; }
00182   pat::Jet                    const &   caloJet1() const { return *caloJet1_; }
00183 
00184   pat::Jet                    const &   pfJet0() const { return *pfJet0_; }
00185   pat::Jet                    const &   pfJet1() const { return *pfJet1_; }
00186 
00187 
00188   // Fast caching indices
00189   index_type const &   caloCuts()        const {return caloCuts_;}      
00190   index_type const &   caloKin()         const {return caloKin_;}       
00191   index_type const &   caloDeltaPhi()    const {return caloDeltaPhi_;}  
00192   index_type const &   caloJetID()       const {return caloJetID_;}     
00193   index_type const &   pfCuts()          const {return pfCuts_;}        
00194   index_type const &   pfKin()           const {return pfKin_;}         
00195   index_type const &   pfDeltaPhi()      const {return pfDeltaPhi_;}    
00196   index_type const &   pfJetID()         const {return pfJetID_;}       
00197 
00198 
00199 protected:
00200   boost::shared_ptr<JetIDSelectionFunctor>   jetSel_;
00201   boost::shared_ptr<PFJetIDSelectionFunctor> pfJetSel_;
00202   edm::InputTag                              jetSrc_;
00203   edm::InputTag                              pfJetSrc_;
00204   
00205   edm::Handle<vector<pat::Jet> >             h_jets_;
00206   edm::Handle<vector<pat::Jet> >             h_pfjets_;
00207 
00208   edm::Ptr<pat::Jet>                         caloJet0_;
00209   edm::Ptr<pat::Jet>                         caloJet1_;
00210 
00211   edm::Ptr<pat::Jet>                         pfJet0_;
00212   edm::Ptr<pat::Jet>                         pfJet1_;
00213 
00214   // Fast caching indices
00215   index_type   caloCuts_;      
00216   index_type   caloKin_;       
00217   index_type   caloDeltaPhi_;  
00218   index_type   caloJetID_;     
00219   index_type   pfCuts_;        
00220   index_type   pfKin_;         
00221   index_type   pfDeltaPhi_;    
00222   index_type   pfJetID_;       
00223 
00224 };
00225 
00227 // ///////////////////// //
00228 // // Main Subroutine // //
00229 // ///////////////////// //
00231 
00232 int main (int argc, char* argv[]) 
00233 {
00234 
00235 
00236   if ( argc < 2 ) {
00237     std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
00238     return 0;
00239   }
00240 
00241   cout << "Hello from " << argv[0] << "!" << endl;
00242 
00243 
00244   // load framework libraries
00245   gSystem->Load( "libFWCoreFWLite" );
00246   AutoLibraryLoader::enable();
00247 
00248 
00249   cout << "Getting parameters" << endl;
00250   // Get the python configuration
00251   PythonProcessDesc builder(argv[1]);
00252   boost::shared_ptr<edm::ProcessDesc> b = builder.processDesc();
00253   boost::shared_ptr<edm::ParameterSet> parameters = b->getProcessPSet();
00254   parameters->registerIt(); 
00255 
00256   edm::ParameterSet const& jetStudiesParams    = parameters->getParameter<edm::ParameterSet>("jetStudies");
00257   edm::ParameterSet const& pfJetStudiesParams  = parameters->getParameter<edm::ParameterSet>("pfJetStudies");
00258   edm::ParameterSet const& caloJetIDParameters = parameters->getParameter<edm::ParameterSet>("jetIDSelector");
00259   edm::ParameterSet const& pfJetIDParameters   = parameters->getParameter<edm::ParameterSet>("pfJetIDSelector");
00260   edm::ParameterSet const& plotParameters      = parameters->getParameter<edm::ParameterSet>("plotParameters");
00261   edm::ParameterSet const& inputs              = parameters->getParameter<edm::ParameterSet>("inputs");
00262   edm::ParameterSet const& outputs             = parameters->getParameter<edm::ParameterSet>("outputs");
00263 
00264   cout << "Making RunLumiSelector" << endl;
00265   RunLumiSelector runLumiSel( inputs );
00266   
00267   cout << "setting up TFileService" << endl;
00268   // book a set of histograms
00269   fwlite::TFileService fs = fwlite::TFileService( outputs.getParameter<std::string>("outputName") );
00270   TFileDirectory theDir = fs.mkdir( "histos" ); 
00271     
00272   cout << "Setting up chain event" << endl;
00273   // This object 'event' is used both to get all information from the
00274   // event as well as to store histograms, etc.
00275   fwlite::ChainEvent ev ( inputs.getParameter<std::vector<std::string> > ("fileNames") );
00276 
00277 
00278   cout << "Booking histograms" << endl;
00279    // Book histograms
00280 
00281   std::map<std::string, TH1*> hists;
00282 
00283    hists["hist_nJet"                     ] = theDir.make<TH1D>( "hist_nJet"                    ,"Number of Calo Jets", 10, 0, 10 ) ;
00284    hists["hist_nPFJet"                   ] = theDir.make<TH1D>( "hist_nPFJet"                  ,"Number of PF Jets", 10, 0, 10 ) ;
00285 
00286    hists["hist_jetPt"                    ] = theDir.make<TH1D>( "hist_jetPt"                   , "Jet p_{T}", 400, 0, 400 ) ;
00287    hists["hist_jetEtaVsPhi"              ] = theDir.make<TH2D>( "hist_jetEtaVsPhi"             , "Jet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
00288    hists["hist_jetNTracks"               ] = theDir.make<TH1D>( "hist_jetNTracks"              , "Jet N_{TRACKS}", 20, 0, 20 ) ;
00289    hists["hist_jetNTracksVsPt"           ] = theDir.make<TH2D>( "hist_jetNTracksVsPt"          , "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
00290    hists["hist_jetEMF"                   ] = theDir.make<TH1D>( "hist_jetEMF"                  , "Jet EMF", 200, 0, 1) ;
00291    hists["hist_jetPdgID"                 ] = theDir.make<TH1D>( "hist_jetPdgID"                , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
00292    hists["hist_jetGenEmE"                ] = theDir.make<TH1D>( "hist_jetGenEmE"               , "Gen Jet EM Energy", 200, 0, 200 ) ;
00293    hists["hist_jetGenHadE"               ] = theDir.make<TH1D>( "hist_jetGenHadE"              , "Gen Jet HAD Energy", 200, 0, 200 ) ;
00294    hists["hist_jetGenEMF"                ] = theDir.make<TH1D>( "hist_jetGenEMF"               , "Gen Jet EMF", 200, 0, 1) ;
00295    hists["hist_jetEoverGenE"             ] = theDir.make<TH1D>( "hist_jetEoverGenE"            , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
00296    hists["hist_jetCorr"                  ] = theDir.make<TH1D>( "hist_jetCorr"                 , "Jet Correction Factor", 200, 0, 1.0 ) ;
00297    hists["hist_n90Hits"                  ] = theDir.make<TH1D>( "hist_n90Hits"                 , "Jet n90Hits", 20, 0, 20) ;
00298    hists["hist_fHPD"                     ] = theDir.make<TH1D>( "hist_fHPD"                    , "Jet fHPD", 200, 0, 1) ;
00299    hists["hist_nConstituents"            ] = theDir.make<TH1D>( "hist_nConstituents"           , "Jet nConstituents", 20, 0, 20 ) ;
00300    hists["hist_jetCHF"                   ] = theDir.make<TH1D>( "hist_jetCHF"                  , "Jet Charged Hadron Fraction", 200, 0, 1.0) ;
00301                                               
00302    hists["hist_good_jetPt"             ] = theDir.make<TH1D>( "hist_good_jetPt"            , "Jet p_{T}", 400, 0, 400 ) ;
00303    hists["hist_good_jetEtaVsPhi"       ] = theDir.make<TH2D>( "hist_good_jetEtaVsPhi"      , "Jet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
00304    hists["hist_good_jetNTracks"        ] = theDir.make<TH1D>( "hist_good_jetNTracks"       , "Jet N_{TRACKS}", 20, 0, 20 ) ;
00305    hists["hist_good_jetNTracksVsPt"    ] = theDir.make<TH2D>( "hist_good_jetNTracksVsPt"   , "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
00306    hists["hist_good_jetEMF"            ] = theDir.make<TH1D>( "hist_good_jetEMF"           , "Jet EMF", 200, 0, 1) ;
00307    hists["hist_good_jetPdgID"          ] = theDir.make<TH1D>( "hist_good_jetPdgID"         , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
00308    hists["hist_good_jetGenEmE"         ] = theDir.make<TH1D>( "hist_good_jetGenEmE"        , "Gen Jet EM Energy", 200, 0, 200 ) ;
00309    hists["hist_good_jetGenHadE"        ] = theDir.make<TH1D>( "hist_good_jetGenHadE"       , "Gen Jet HAD Energy", 200, 0, 200 ) ;
00310    hists["hist_good_jetGenEMF"         ] = theDir.make<TH1D>( "hist_good_jetGenEMF"        , "Gen Jet EMF", 200, 0, 1) ;
00311    hists["hist_good_jetEoverGenE"      ] = theDir.make<TH1D>( "hist_good_jetEoverGenE"     , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
00312    hists["hist_good_jetCorr"           ] = theDir.make<TH1D>( "hist_good_jetCorr"          , "Jet Correction Factor", 200, 0, 1.0 ) ;
00313    hists["hist_good_n90Hits"           ] = theDir.make<TH1D>( "hist_good_n90Hits"          , "Jet n90Hits", 20, 0, 20) ;
00314    hists["hist_good_fHPD"              ] = theDir.make<TH1D>( "hist_good_fHPD"             , "Jet fHPD", 200, 0, 1) ;
00315    hists["hist_good_nConstituents"     ] = theDir.make<TH1D>( "hist_good_nConstituents"    , "Jet nConstituents", 20, 0, 20 ) ;
00316    hists["hist_good_jetCHF"            ] = theDir.make<TH1D>( "hist_good_jetCHF"           , "Jet Charged Hadron Fraction", 200, 0, 1.0) ;
00317                                                           
00318                                               
00319    hists["hist_pf_jetPt"                 ] = theDir.make<TH1D>( "hist_pf_jetPt"                , "PFJet p_{T}", 400, 0, 400 ) ;
00320    hists["hist_pf_jetEtaVsPhi"           ] = theDir.make<TH2D>( "hist_pf_jetEtaVsPhi"          , "PFJet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
00321    hists["hist_pf_jetNTracks"            ] = theDir.make<TH1D>( "hist_pf_jetNTracks"           , "PFJet N_{TRACKS}", 20, 0, 20 ) ;
00322    hists["hist_pf_jetNTracksVsPt"        ] = theDir.make<TH2D>( "hist_pf_jetNTracksVsPt"       , "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
00323    hists["hist_pf_jetEMF"                ] = theDir.make<TH1D>( "hist_pf_jetEMF"               , "PFJet EMF", 200, 0, 1) ;
00324    hists["hist_pf_jetCHF"                ] = theDir.make<TH1D>( "hist_pf_jetCHF"               , "PFJet CHF", 200, 0, 1) ;
00325    hists["hist_pf_jetNHF"                ] = theDir.make<TH1D>( "hist_pf_jetNHF"               , "PFJet NHF", 200, 0, 1) ;
00326    hists["hist_pf_jetCEF"                ] = theDir.make<TH1D>( "hist_pf_jetCEF"               , "PFJet CEF", 200, 0, 1) ;
00327    hists["hist_pf_jetNEF"                ] = theDir.make<TH1D>( "hist_pf_jetNEF"               , "PFJet NEF", 200, 0, 1) ;
00328    hists["hist_pf_jetPdgID"              ] = theDir.make<TH1D>( "hist_pf_jetPdgID"             , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
00329    hists["hist_pf_jetGenEmE"             ] = theDir.make<TH1D>( "hist_pf_jetGenEmE"            , "Gen Jet EM Energy", 200, 0, 200 ) ;
00330    hists["hist_pf_jetGenHadE"            ] = theDir.make<TH1D>( "hist_pf_jetGenHadE"           , "Gen Jet HAD Energy", 200, 0, 200 ) ;
00331    hists["hist_pf_jetGenEMF"             ] = theDir.make<TH1D>( "hist_pf_jetGenEMF"            , "Gen Jet EMF", 200, 0, 1) ;
00332    hists["hist_pf_jetEoverGenE"          ] = theDir.make<TH1D>( "hist_pf_jetEoverGenE"         , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
00333    hists["hist_pf_jetCorr"               ] = theDir.make<TH1D>( "hist_pf_jetCorr"              , "PFJet Correction Factor", 200, 0, 1.0 ) ;
00334    hists["hist_pf_nConstituents"         ] = theDir.make<TH1D>( "hist_pf_nConstituents"        , "PFJet nConstituents", 20, 0, 20 ) ;
00335                                               
00336    hists["hist_pf_good_jetPt"          ] = theDir.make<TH1D>( "hist_pf_good_jetPt"         , "PFJet p_{T}", 400, 0, 400 ) ;
00337    hists["hist_pf_good_jetEtaVsPhi"    ] = theDir.make<TH2D>( "hist_pf_good_jetEtaVsPhi"   , "PFJet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
00338    hists["hist_pf_good_jetNTracks"     ] = theDir.make<TH1D>( "hist_pf_good_jetNTracks"    , "PFJet N_{TRACKS}", 20, 0, 20 ) ;
00339    hists["hist_pf_good_jetNTracksVsPt" ] = theDir.make<TH2D>( "hist_pf_good_jetNTracksVsPt", "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
00340    hists["hist_pf_good_jetEMF"         ] = theDir.make<TH1D>( "hist_pf_good_jetEMF"        , "PFJet EMF", 200, 0, 1) ;
00341    hists["hist_pf_good_jetCHF"         ] = theDir.make<TH1D>( "hist_pf_good_jetCHF"        , "PFJet CHF", 200, 0, 1) ;
00342    hists["hist_pf_good_jetNHF"         ] = theDir.make<TH1D>( "hist_pf_good_jetNHF"        , "PFJet NHF", 200, 0, 1) ;
00343    hists["hist_pf_good_jetCEF"         ] = theDir.make<TH1D>( "hist_pf_good_jetCEF"        , "PFJet CEF", 200, 0, 1) ;
00344    hists["hist_pf_good_jetNEF"         ] = theDir.make<TH1D>( "hist_pf_good_jetNEF"        , "PFJet NEF", 200, 0, 1) ;
00345    hists["hist_pf_good_jetPdgID"       ] = theDir.make<TH1D>( "hist_pf_good_jetPdgID"      , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
00346    hists["hist_pf_good_jetGenEmE"      ] = theDir.make<TH1D>( "hist_pf_good_jetGenEmE"     , "Gen Jet EM Energy", 200, 0, 200 ) ;
00347    hists["hist_pf_good_jetGenHadE"     ] = theDir.make<TH1D>( "hist_pf_good_jetGenHadE"    , "Gen Jet HAD Energy", 200, 0, 200 ) ;
00348    hists["hist_pf_good_jetGenEMF"      ] = theDir.make<TH1D>( "hist_pf_good_jetGenEMF"     , "Gen Jet EMF", 200, 0, 1) ;
00349    hists["hist_pf_good_jetEoverGenE"   ] = theDir.make<TH1D>( "hist_pf_good_jetEoverGenE"  , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
00350    hists["hist_pf_good_jetCorr"        ] = theDir.make<TH1D>( "hist_pf_good_jetCorr"       , "PFJet Correction Factor", 200, 0, 1.0 ) ;
00351    hists["hist_pf_good_nConstituents"  ] = theDir.make<TH1D>( "hist_pf_good_nConstituents" , "PFJet nConstituents", 20, 0, 20 ) ;
00352 
00353    hists["hist_mjj"                           ] = theDir.make<TH1D>( "hist_mjj"                          , "Dijet mass", 300, 0, 300 ) ;
00354    hists["hist_dR_jj"                         ] = theDir.make<TH1D>( "hist_dR_jj"                        , "#Delta R_{JJ}", 200, 0, TMath::TwoPi() ) ;
00355    hists["hist_imbalance_jj"                  ] = theDir.make<TH1D>( "hist_imbalance_jj"                 , "Dijet imbalance", 200, -1.0, 1.0 )  ;
00356                                               
00357    hists["hist_pf_mjj"                        ] = theDir.make<TH1D>( "hist_pf_mjj"                       , "Dijet mass", 300, 0, 300 ) ;
00358    hists["hist_pf_dR_jj"                      ] = theDir.make<TH1D>( "hist_pf_dR_jj"                     , "#Delta R_{JJ}", 200, 0, TMath::TwoPi() ) ;
00359    hists["hist_pf_imbalance_jj"               ] = theDir.make<TH1D>( "hist_pf_imbalance_jj"              , "Dijet imbalance", 200, -1.0, 1.0 )  ;
00360 
00361 
00362    
00363    cout << "Making functors" << endl;
00364    JetIDStudiesSelector caloSelector( caloJetIDParameters,
00365                                       pfJetIDParameters,
00366                                       jetStudiesParams );
00367 
00368    JetIDStudiesSelector pfSelector( caloJetIDParameters,
00369                                     pfJetIDParameters,
00370                                     pfJetStudiesParams );
00371 
00372    bool doTracks = plotParameters.getParameter<bool>("doTracks");
00373    bool useMC    = plotParameters.getParameter<bool>("useMC");
00374 
00375    cout << "About to begin looping" << endl;
00376 
00377   int nev = 0;
00378   //loop through each event
00379   for (ev.toBegin(); ! ev.atEnd(); ++ev, ++nev) {
00380 
00381 
00382     edm::EventBase const & event = ev;
00383 
00384     if ( runLumiSel(ev) == false ) continue;
00385 
00386     if ( nev % 100 == 0 ) cout << "Processing run " << event.id().run() << ", lumi " << event.id().luminosityBlock() << ", event " << event.id().event() << endl;
00387 
00388 
00389 
00390 
00391     pat::strbitset retCalo = caloSelector.getBitTemplate();
00392     caloSelector( event, retCalo );
00393 
00394 
00395     pat::strbitset retPF = pfSelector.getBitTemplate();
00396     pfSelector( event, retPF );
00397 
00401     if ( retCalo.test( caloSelector.caloKin() ) ) {
00402       if ( retCalo.test( caloSelector.caloDeltaPhi() ) ) {
00403         vector<pat::Jet>  const & allCaloJets = caloSelector.allCaloJets();
00404 
00405         for ( std::vector<pat::Jet>::const_iterator jetBegin = allCaloJets.begin(),
00406                 jetEnd = jetBegin + 2, ijet = jetBegin;
00407               ijet != jetEnd; ++ijet ) {
00408         
00409           const pat::Jet & jet = *ijet;
00410 
00411           double pt = jet.pt();
00412 
00413           const reco::TrackRefVector & jetTracks = jet.associatedTracks();
00414 
00415           hists["hist_jetPt"]->Fill( pt );
00416           hists["hist_jetEtaVsPhi"]->Fill( jet.eta(), jet.phi() );
00417           hists["hist_jetNTracks"]->Fill( jetTracks.size() );
00418           hists["hist_jetNTracksVsPt"]->Fill( pt, jetTracks.size() );
00419           hists["hist_jetEMF"]->Fill( jet.emEnergyFraction() ); 
00420           hists["hist_jetCorr"]->Fill( jet.jecFactor("Uncorrected") );
00421           hists["hist_n90Hits"]->Fill( static_cast<int>(jet.jetID().n90Hits) );
00422           hists["hist_fHPD"]->Fill( jet.jetID().fHPD );
00423           hists["hist_nConstituents"]->Fill( jet.nConstituents() );
00424 
00425           if ( useMC && jet.genJet() != 0 ) {
00426             hists["hist_jetGenEmE"]->Fill( jet.genJet()->emEnergy() );
00427             hists["hist_jetGenHadE"]->Fill( jet.genJet()->hadEnergy() );
00428             hists["hist_jetEoverGenE"]->Fill( jet.energy() / jet.genJet()->energy() );
00429             hists["hist_jetGenEMF"]->Fill( jet.genJet()->emEnergy() / jet.genJet()->energy() );
00430           }
00431           if ( doTracks ) {
00432             TLorentzVector p4_tracks(0,0,0,0);
00433             for ( reco::TrackRefVector::const_iterator itrk = jetTracks.begin(),
00434                     itrkEnd = jetTracks.end();
00435                   itrk != itrkEnd; ++itrk ) {
00436               TLorentzVector p4_trk;
00437               double M_PION = 0.140;
00438               p4_trk.SetPtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), M_PION );
00439               p4_tracks += p4_trk;
00440             }
00441             hists["hist_jetCHF"]->Fill( p4_tracks.Energy() / jet.energy() );
00442           }
00443 
00444       
00445         } // end loop over jets
00446 
00447     
00448 
00449         if ( retCalo.test( caloSelector.caloJetID() ) ) {
00450           pat::Jet const & jet0 = caloSelector.caloJet0();
00451           pat::Jet const & jet1 = caloSelector.caloJet1();
00452 
00453           TLorentzVector p4_j0( jet0.px(), jet0.py(), jet0.pz(), jet0.energy() );
00454           TLorentzVector p4_j1( jet1.px(), jet1.py(), jet1.pz(), jet1.energy() );
00455 
00456           TLorentzVector p4_jj = p4_j0 + p4_j1;
00457 
00458           hists["hist_mjj"]->Fill( p4_jj.M() );
00459           hists["hist_dR_jj"]->Fill( p4_j0.DeltaR( p4_j1 ) );
00460           hists["hist_imbalance_jj"]->Fill( (p4_j0.Perp() - p4_j1.Perp() ) /
00461                                                      (p4_j0.Perp() + p4_j1.Perp() ) );
00462 
00463           hists["hist_good_jetPt"]->Fill( jet0.pt() );
00464           hists["hist_good_jetEtaVsPhi"]->Fill( jet0.eta(), jet0.phi() );
00465           hists["hist_good_jetNTracks"]->Fill( jet0.associatedTracks().size() );
00466           hists["hist_good_jetNTracksVsPt"]->Fill( jet0.pt(), jet0.associatedTracks().size() );
00467           hists["hist_good_jetEMF"]->Fill( jet0.emEnergyFraction() );   
00468           hists["hist_good_jetCorr"]->Fill( jet0.jecFactor("Uncorrected") );
00469           hists["hist_good_n90Hits"]->Fill( static_cast<int>(jet0.jetID().n90Hits) );
00470           hists["hist_good_fHPD"]->Fill( jet0.jetID().fHPD );
00471           hists["hist_good_nConstituents"]->Fill( jet0.nConstituents() );
00472 
00473 
00474           hists["hist_good_jetPt"]->Fill( jet1.pt() );
00475           hists["hist_good_jetEtaVsPhi"]->Fill( jet1.eta(), jet1.phi() );
00476           hists["hist_good_jetNTracks"]->Fill( jet1.associatedTracks().size() );
00477           hists["hist_good_jetNTracksVsPt"]->Fill( jet1.pt(), jet1.associatedTracks().size() );
00478           hists["hist_good_jetEMF"]->Fill( jet1.emEnergyFraction() );   
00479           hists["hist_good_jetCorr"]->Fill( jet1.jecFactor("Uncorrected") );
00480           hists["hist_good_n90Hits"]->Fill( static_cast<int>(jet1.jetID().n90Hits) );
00481           hists["hist_good_fHPD"]->Fill( jet1.jetID().fHPD );
00482           hists["hist_good_nConstituents"]->Fill( jet1.nConstituents() );
00483 
00484         }// end if passed calo jet id
00485       }// end if passed dphi cuts
00486     }// end if passed kin cuts
00487 
00488 
00492     if ( retPF.test( pfSelector.pfDeltaPhi() ) ) {
00493 
00494       vector<pat::Jet> const & allPFJets = pfSelector.allPFJets();
00495 
00496       for ( std::vector<pat::Jet>::const_iterator jetBegin = allPFJets.begin(),
00497               jetEnd = jetBegin + 2, ijet = jetBegin;
00498             ijet != jetEnd; ++ijet ) {
00499         
00500         const pat::Jet & jet = *ijet;
00501         
00502         double pt = jet.pt();
00503       
00504         hists["hist_pf_jetPt"]->Fill( pt );
00505         hists["hist_pf_jetEtaVsPhi"]->Fill( jet.eta(), jet.phi() );
00506         hists["hist_pf_nConstituents"]->Fill( jet.nConstituents() );
00507         hists["hist_pf_jetCEF"]->Fill( jet.chargedEmEnergyFraction()  );
00508         hists["hist_pf_jetNEF"]->Fill( jet.neutralEmEnergyFraction()  );
00509         hists["hist_pf_jetCHF"]->Fill( jet.chargedHadronEnergyFraction()  );
00510         hists["hist_pf_jetNHF"]->Fill( jet.neutralHadronEnergyFraction()  );
00511 
00512 
00513         if ( useMC && jet.genJet() != 0 ) {
00514           hists["hist_pf_jetGenEmE"]->Fill( jet.genJet()->emEnergy() );
00515           hists["hist_pf_jetGenHadE"]->Fill( jet.genJet()->hadEnergy() );
00516           hists["hist_pf_jetEoverGenE"]->Fill( jet.energy() / jet.genJet()->energy() );
00517 
00518           hists["hist_pf_jetGenEMF"]->Fill( jet.genJet()->emEnergy() / jet.genJet()->energy() );
00519         }
00520  
00521       } // end loop over jets
00522 
00523     
00524 
00525       if ( retPF.test( pfSelector.pfJetID() ) ) {
00526         pat::Jet const & jet0 = pfSelector.pfJet0();
00527         pat::Jet const & jet1 = pfSelector.pfJet1();
00528 
00529         TLorentzVector p4_j0( jet0.px(), jet0.py(), jet0.pz(), jet0.energy() );
00530         TLorentzVector p4_j1( jet1.px(), jet1.py(), jet1.pz(), jet1.energy() );
00531 
00532         TLorentzVector p4_jj = p4_j0 + p4_j1;
00533 
00534         hists["hist_pf_mjj"]->Fill( p4_jj.M() );
00535         hists["hist_pf_dR_jj"]->Fill( p4_j0.DeltaR( p4_j1 ) );
00536         hists["hist_pf_imbalance_jj"]->Fill( (p4_j0.Perp() - p4_j1.Perp() ) /
00537                                                       (p4_j0.Perp() + p4_j1.Perp() ) );
00538 
00539         hists["hist_pf_good_jetPt"]->Fill( jet0.pt() );
00540         hists["hist_pf_good_jetEtaVsPhi"]->Fill( jet0.eta(), jet0.phi() );
00541         hists["hist_pf_good_nConstituents"]->Fill( jet0.nConstituents() );
00542         hists["hist_pf_good_jetCEF"]->Fill( jet0.chargedEmEnergyFraction()  );
00543         hists["hist_pf_good_jetNEF"]->Fill( jet0.neutralEmEnergyFraction()  );
00544         hists["hist_pf_good_jetCHF"]->Fill( jet0.chargedHadronEnergyFraction()  );
00545         hists["hist_pf_good_jetNHF"]->Fill( jet0.neutralHadronEnergyFraction()  );
00546 
00547 
00548         hists["hist_pf_good_jetPt"]->Fill( jet1.pt() );
00549         hists["hist_pf_good_jetEtaVsPhi"]->Fill( jet1.eta(), jet1.phi() );
00550         hists["hist_pf_good_nConstituents"]->Fill( jet1.nConstituents() );
00551         hists["hist_pf_good_jetCEF"]->Fill( jet1.chargedEmEnergyFraction()  );
00552         hists["hist_pf_good_jetNEF"]->Fill( jet1.neutralEmEnergyFraction()  );
00553         hists["hist_pf_good_jetCHF"]->Fill( jet1.chargedHadronEnergyFraction()  );
00554         hists["hist_pf_good_jetNHF"]->Fill( jet1.neutralHadronEnergyFraction()  );
00555 
00556       } // end if 2 good PF jets
00557     
00558     }// end if delta phi pf cuts
00559     
00560   } // end loop over events
00561     
00562   cout << "Calo jet selection" << endl;
00563   caloSelector.print(std::cout);
00564   cout << "PF jet selection" << endl;
00565   pfSelector.print(std::cout);
00566 
00567 
00568 
00569 
00570   return 0;
00571 }
00572
/data/refman/pasoursint/CMSSW_5_2_9/src/PhysicsTools/PatExamples/bin/PatBasicFWLiteJetAnalyzer_Selector.cc
