Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:00 

0001 /*   A macro for making a histogram of Jet Pt with cuts
0002 This is a basic way to cut out jets of a certain Pt and Eta using an if statement
0003 This example creates a histogram of Jet Pt, using Jets with Pt above 30 and ETA above -2.1 and below 2.1
0004 */
0005 
0006 #include "TFile.h"
0007 #include "TH1.h"
0008 #include "TH2.h"
0009 #include "TCanvas.h"
0010 #include "TLegend.h"
0011 #include "TSystem.h"
0012 #include "TLorentzVector.h"
0013 
0014 #include "PhysicsTools/SelectorUtils/interface/strbitset.h"
0015 #include "DataFormats/PatCandidates/interface/Jet.h"
0016 #include "DataFormats/PatCandidates/interface/MET.h"
0017 #include "DataFormats/PatCandidates/interface/Photon.h"
0018 #include "DataFormats/PatCandidates/interface/Muon.h"
0019 #include "DataFormats/PatCandidates/interface/Electron.h"
0020 #include "DataFormats/PatCandidates/interface/Tau.h"
0021 #include "PhysicsTools/SelectorUtils/interface/Selector.h"
0022 #include "PhysicsTools/SelectorUtils/interface/JetIDSelectionFunctor.h"
0023 #include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
0024 #include "PhysicsTools/SelectorUtils/interface/RunLumiSelector.h"
0025 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0026 #include "PhysicsTools/FWLite/interface/TFileService.h"
0027 #include "DataFormats/FWLite/interface/ChainEvent.h"
0028 #include "FWCore/ParameterSetReader/interface/ProcessDescImpl.h"
0029 #include "FWCore/ParameterSet/interface/ProcessDesc.h"
0030 
0031 #include <iostream>
0032 #include <memory>
0033 #include <cmath>  //necessary for absolute function fabs()
0034 
0035 using namespace std;
0036 
0037 ///-------------------------
0038 /// DRIVER FUNCTION
0039 ///-------------------------
0040 
0041 // -*- C++ -*-
0042 
0043 // CMS includes
0044 #include "FWCore/Utilities/interface/InputTag.h"
0045 #include "DataFormats/Common/interface/Handle.h"
0046 #include "DataFormats/Common/interface/Ptr.h"
0047 #include "DataFormats/PatCandidates/interface/Jet.h"
0048 #include "PhysicsTools/SelectorUtils/interface/EventSelector.h"
0049 
0050 // Root includes
0051 #include "TROOT.h"
0052 
0053 using namespace std;
0054 
0055 class JetIDStudiesSelector : public EventSelector {
0056 public:
0057   JetIDStudiesSelector(edm::ParameterSet const& caloJetIdParams,
0058                        edm::ParameterSet const& pfJetIdParams,
0059                        edm::ParameterSet const& params)
0060       : jetSel_(new JetIDSelectionFunctor(caloJetIdParams)),
0061         pfJetSel_(new PFJetIDSelectionFunctor(pfJetIdParams)),
0062         jetSrc_(params.getParameter<edm::InputTag>("jetSrc")),
0063         pfJetSrc_(params.getParameter<edm::InputTag>("pfJetSrc")) {
0064     bool useCalo = params.getParameter<bool>("useCalo");
0065 
0066     push_back("Calo Cuts");
0067     push_back("Calo Kin Cuts");
0068     push_back("Calo Delta Phi");
0069     push_back("Calo Jet ID");
0070     push_back("PF Cuts");
0071     push_back("PF Kin Cuts");
0072     push_back("PF Delta Phi");
0073     push_back("PF Jet ID");
0074 
0075     set("Calo Cuts", useCalo);
0076     set("Calo Kin Cuts", useCalo);
0077     set("Calo Delta Phi", useCalo);
0078     set("Calo Jet ID", useCalo);
0079     set("PF Cuts", !useCalo);
0080     set("PF Kin Cuts", !useCalo);
0081     set("PF Delta Phi", !useCalo);
0082     set("PF Jet ID", !useCalo);
0083 
0084     // Indices for fast caching
0085     caloCuts_ = index_type(&bits_, std::string("Calo Cuts"));
0086     caloKin_ = index_type(&bits_, std::string("Calo Kin Cuts"));
0087     caloDeltaPhi_ = index_type(&bits_, std::string("Calo Delta Phi"));
0088     caloJetID_ = index_type(&bits_, std::string("Calo Jet ID"));
0089     pfCuts_ = index_type(&bits_, std::string("PF Cuts"));
0090     pfKin_ = index_type(&bits_, std::string("PF Kin Cuts"));
0091     pfDeltaPhi_ = index_type(&bits_, std::string("PF Delta Phi"));
0092     pfJetID_ = index_type(&bits_, std::string("PF Jet ID"));
0093   }
0094 
0095   ~JetIDStudiesSelector() override {}
0096 
0097   bool operator()(edm::EventBase const& event, pat::strbitset& ret) override {
0098     pat::strbitset retCaloJet = jetSel_->getBitTemplate();
0099     pat::strbitset retPFJet = pfJetSel_->getBitTemplate();
0100 
0101     if (considerCut(caloCuts_)) {
0102       passCut(ret, caloCuts_);
0103       event.getByLabel(jetSrc_, h_jets_);
0104       // Calo Cuts
0105       if (h_jets_->size() >= 2 || ignoreCut(caloKin_)) {
0106         passCut(ret, caloKin_);
0107         pat::Jet const& jet0 = h_jets_->at(0);
0108         pat::Jet const& jet1 = h_jets_->at(1);
0109         double dphi = fabs(deltaPhi<double>(jet0.phi(), jet1.phi()));
0110 
0111         if (fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(caloDeltaPhi_)) {
0112           passCut(ret, caloDeltaPhi_);
0113 
0114           retCaloJet.set(false);
0115           bool pass0 = (*jetSel_)(jet0, retCaloJet);
0116           retCaloJet.set(false);
0117           bool pass1 = (*jetSel_)(jet1, retCaloJet);
0118           if ((pass0 && pass1) || ignoreCut(caloJetID_)) {
0119             passCut(ret, caloJetID_);
0120             caloJet0_ = edm::Ptr<pat::Jet>(h_jets_, 0);
0121             caloJet1_ = edm::Ptr<pat::Jet>(h_jets_, 1);
0122 
0123             return (bool)ret;
0124           }  // end if found 2 "loose" jet ID jets
0125         }    // end if delta phi
0126       }      // end calo kin cuts
0127     }        // end if calo cuts
0128 
0129     if (considerCut(pfCuts_)) {
0130       passCut(ret, pfCuts_);
0131       event.getByLabel(pfJetSrc_, h_pfjets_);
0132       // PF Cuts
0133       if (h_pfjets_->size() >= 2 || ignoreCut(pfKin_)) {
0134         passCut(ret, pfKin_);
0135         pat::Jet const& jet0 = h_pfjets_->at(0);
0136         pat::Jet const& jet1 = h_pfjets_->at(1);
0137         double dphi = fabs(deltaPhi<double>(jet0.phi(), jet1.phi()));
0138 
0139         if (fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(pfDeltaPhi_)) {
0140           passCut(ret, pfDeltaPhi_);
0141 
0142           retPFJet.set(false);
0143           bool pass0 = (*pfJetSel_)(jet0, retPFJet);
0144           retPFJet.set(false);
0145           bool pass1 = (*pfJetSel_)(jet1, retPFJet);
0146           if ((pass0 && pass1) || ignoreCut(pfJetID_)) {
0147             passCut(ret, pfJetID_);
0148             pfJet0_ = edm::Ptr<pat::Jet>(h_pfjets_, 0);
0149             pfJet1_ = edm::Ptr<pat::Jet>(h_pfjets_, 1);
0150 
0151             return (bool)ret;
0152           }  // end if found 2 "loose" jet ID jets
0153         }    // end if delta phi
0154       }      // end pf kin cuts
0155     }        // end if pf cuts
0156 
0157     setIgnored(ret);
0158 
0159     return false;
0160   }  // end of method
0161 
0162   std::shared_ptr<JetIDSelectionFunctor> const& jetSel() const { return jetSel_; }
0163   std::shared_ptr<PFJetIDSelectionFunctor> const& pfJetSel() const { return pfJetSel_; }
0164 
0165   vector<pat::Jet> const& allCaloJets() const { return *h_jets_; }
0166   vector<pat::Jet> const& allPFJets() const { return *h_pfjets_; }
0167 
0168   pat::Jet const& caloJet0() const { return *caloJet0_; }
0169   pat::Jet const& caloJet1() const { return *caloJet1_; }
0170 
0171   pat::Jet const& pfJet0() const { return *pfJet0_; }
0172   pat::Jet const& pfJet1() const { return *pfJet1_; }
0173 
0174   // Fast caching indices
0175   index_type const& caloCuts() const { return caloCuts_; }
0176   index_type const& caloKin() const { return caloKin_; }
0177   index_type const& caloDeltaPhi() const { return caloDeltaPhi_; }
0178   index_type const& caloJetID() const { return caloJetID_; }
0179   index_type const& pfCuts() const { return pfCuts_; }
0180   index_type const& pfKin() const { return pfKin_; }
0181   index_type const& pfDeltaPhi() const { return pfDeltaPhi_; }
0182   index_type const& pfJetID() const { return pfJetID_; }
0183 
0184 protected:
0185   std::shared_ptr<JetIDSelectionFunctor> jetSel_;
0186   std::shared_ptr<PFJetIDSelectionFunctor> pfJetSel_;
0187   edm::InputTag jetSrc_;
0188   edm::InputTag pfJetSrc_;
0189 
0190   edm::Handle<vector<pat::Jet> > h_jets_;
0191   edm::Handle<vector<pat::Jet> > h_pfjets_;
0192 
0193   edm::Ptr<pat::Jet> caloJet0_;
0194   edm::Ptr<pat::Jet> caloJet1_;
0195 
0196   edm::Ptr<pat::Jet> pfJet0_;
0197   edm::Ptr<pat::Jet> pfJet1_;
0198 
0199   // Fast caching indices
0200   index_type caloCuts_;
0201   index_type caloKin_;
0202   index_type caloDeltaPhi_;
0203   index_type caloJetID_;
0204   index_type pfCuts_;
0205   index_type pfKin_;
0206   index_type pfDeltaPhi_;
0207   index_type pfJetID_;
0208 };
0209 
0210 ///////////////////////////
0211 // ///////////////////// //
0212 // // Main Subroutine // //
0213 // ///////////////////// //
0214 ///////////////////////////
0215 
0216 int main(int argc, char* argv[]) {
0217   if (argc < 2) {
0218     std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
0219     return 0;
0220   }
0221 
0222   cout << "Hello from " << argv[0] << "!" << endl;
0223 
0224   // load framework libraries
0225   gSystem->Load("libFWCoreFWLite");
0226   FWLiteEnabler::enable();
0227 
0228   cout << "Getting parameters" << endl;
0229   // Get the python configuration
0230   ProcessDescImpl builder(argv[1], true);
0231   std::shared_ptr<edm::ProcessDesc> b = builder.processDesc();
0232   std::shared_ptr<edm::ParameterSet> parameters = b->getProcessPSet();
0233   parameters->registerIt();
0234 
0235   edm::ParameterSet const& jetStudiesParams = parameters->getParameter<edm::ParameterSet>("jetStudies");
0236   edm::ParameterSet const& pfJetStudiesParams = parameters->getParameter<edm::ParameterSet>("pfJetStudies");
0237   edm::ParameterSet const& caloJetIDParameters = parameters->getParameter<edm::ParameterSet>("jetIDSelector");
0238   edm::ParameterSet const& pfJetIDParameters = parameters->getParameter<edm::ParameterSet>("pfJetIDSelector");
0239   edm::ParameterSet const& plotParameters = parameters->getParameter<edm::ParameterSet>("plotParameters");
0240   edm::ParameterSet const& inputs = parameters->getParameter<edm::ParameterSet>("inputs");
0241   edm::ParameterSet const& outputs = parameters->getParameter<edm::ParameterSet>("outputs");
0242 
0243   cout << "Making RunLumiSelector" << endl;
0244   RunLumiSelector runLumiSel(inputs);
0245 
0246   cout << "setting up TFileService" << endl;
0247   // book a set of histograms
0248   fwlite::TFileService fs = fwlite::TFileService(outputs.getParameter<std::string>("outputName"));
0249   TFileDirectory theDir = fs.mkdir("histos");
0250 
0251   cout << "Setting up chain event" << endl;
0252   // This object 'event' is used both to get all information from the
0253   // event as well as to store histograms, etc.
0254   fwlite::ChainEvent ev(inputs.getParameter<std::vector<std::string> >("fileNames"));
0255 
0256   cout << "Booking histograms" << endl;
0257   // Book histograms
0258 
0259   std::map<std::string, TH1*> hists;
0260 
0261   hists["hist_nJet"] = theDir.make<TH1D>("hist_nJet", "Number of Calo Jets", 10, 0, 10);
0262   hists["hist_nPFJet"] = theDir.make<TH1D>("hist_nPFJet", "Number of PF Jets", 10, 0, 10);
0263 
0264   hists["hist_jetPt"] = theDir.make<TH1D>("hist_jetPt", "Jet p_{T}", 400, 0, 400);
0265   hists["hist_jetEtaVsPhi"] = theDir.make<TH2D>(
0266       "hist_jetEtaVsPhi", "Jet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi());
0267   hists["hist_jetNTracks"] = theDir.make<TH1D>("hist_jetNTracks", "Jet N_{TRACKS}", 20, 0, 20);
0268   hists["hist_jetNTracksVsPt"] = theDir.make<TH2D>(
0269       "hist_jetNTracksVsPt", "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}", 20, 0, 200, 20, 0, 20);
0270   hists["hist_jetEMF"] = theDir.make<TH1D>("hist_jetEMF", "Jet EMF", 200, 0, 1);
0271   hists["hist_jetPdgID"] = theDir.make<TH1D>("hist_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
0272   hists["hist_jetGenEmE"] = theDir.make<TH1D>("hist_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
0273   hists["hist_jetGenHadE"] = theDir.make<TH1D>("hist_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
0274   hists["hist_jetGenEMF"] = theDir.make<TH1D>("hist_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
0275   hists["hist_jetEoverGenE"] =
0276       theDir.make<TH1D>("hist_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
0277   hists["hist_jetCorr"] = theDir.make<TH1D>("hist_jetCorr", "Jet Correction Factor", 200, 0, 1.0);
0278   hists["hist_n90Hits"] = theDir.make<TH1D>("hist_n90Hits", "Jet n90Hits", 20, 0, 20);
0279   hists["hist_fHPD"] = theDir.make<TH1D>("hist_fHPD", "Jet fHPD", 200, 0, 1);
0280   hists["hist_nConstituents"] = theDir.make<TH1D>("hist_nConstituents", "Jet nConstituents", 20, 0, 20);
0281   hists["hist_jetCHF"] = theDir.make<TH1D>("hist_jetCHF", "Jet Charged Hadron Fraction", 200, 0, 1.0);
0282 
0283   hists["hist_good_jetPt"] = theDir.make<TH1D>("hist_good_jetPt", "Jet p_{T}", 400, 0, 400);
0284   hists["hist_good_jetEtaVsPhi"] = theDir.make<TH2D>(
0285       "hist_good_jetEtaVsPhi", "Jet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi());
0286   hists["hist_good_jetNTracks"] = theDir.make<TH1D>("hist_good_jetNTracks", "Jet N_{TRACKS}", 20, 0, 20);
0287   hists["hist_good_jetNTracksVsPt"] =
0288       theDir.make<TH2D>("hist_good_jetNTracksVsPt",
0289                         "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",
0290                         20,
0291                         0,
0292                         200,
0293                         20,
0294                         0,
0295                         20);
0296   hists["hist_good_jetEMF"] = theDir.make<TH1D>("hist_good_jetEMF", "Jet EMF", 200, 0, 1);
0297   hists["hist_good_jetPdgID"] = theDir.make<TH1D>("hist_good_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
0298   hists["hist_good_jetGenEmE"] = theDir.make<TH1D>("hist_good_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
0299   hists["hist_good_jetGenHadE"] = theDir.make<TH1D>("hist_good_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
0300   hists["hist_good_jetGenEMF"] = theDir.make<TH1D>("hist_good_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
0301   hists["hist_good_jetEoverGenE"] =
0302       theDir.make<TH1D>("hist_good_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
0303   hists["hist_good_jetCorr"] = theDir.make<TH1D>("hist_good_jetCorr", "Jet Correction Factor", 200, 0, 1.0);
0304   hists["hist_good_n90Hits"] = theDir.make<TH1D>("hist_good_n90Hits", "Jet n90Hits", 20, 0, 20);
0305   hists["hist_good_fHPD"] = theDir.make<TH1D>("hist_good_fHPD", "Jet fHPD", 200, 0, 1);
0306   hists["hist_good_nConstituents"] = theDir.make<TH1D>("hist_good_nConstituents", "Jet nConstituents", 20, 0, 20);
0307   hists["hist_good_jetCHF"] = theDir.make<TH1D>("hist_good_jetCHF", "Jet Charged Hadron Fraction", 200, 0, 1.0);
0308 
0309   hists["hist_pf_jetPt"] = theDir.make<TH1D>("hist_pf_jetPt", "PFJet p_{T}", 400, 0, 400);
0310   hists["hist_pf_jetEtaVsPhi"] = theDir.make<TH2D>(
0311       "hist_pf_jetEtaVsPhi", "PFJet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi());
0312   hists["hist_pf_jetNTracks"] = theDir.make<TH1D>("hist_pf_jetNTracks", "PFJet N_{TRACKS}", 20, 0, 20);
0313   hists["hist_pf_jetNTracksVsPt"] = theDir.make<TH2D>(
0314       "hist_pf_jetNTracksVsPt", "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}", 20, 0, 200, 20, 0, 20);
0315   hists["hist_pf_jetEMF"] = theDir.make<TH1D>("hist_pf_jetEMF", "PFJet EMF", 200, 0, 1);
0316   hists["hist_pf_jetCHF"] = theDir.make<TH1D>("hist_pf_jetCHF", "PFJet CHF", 200, 0, 1);
0317   hists["hist_pf_jetNHF"] = theDir.make<TH1D>("hist_pf_jetNHF", "PFJet NHF", 200, 0, 1);
0318   hists["hist_pf_jetCEF"] = theDir.make<TH1D>("hist_pf_jetCEF", "PFJet CEF", 200, 0, 1);
0319   hists["hist_pf_jetNEF"] = theDir.make<TH1D>("hist_pf_jetNEF", "PFJet NEF", 200, 0, 1);
0320   hists["hist_pf_jetPdgID"] = theDir.make<TH1D>("hist_pf_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
0321   hists["hist_pf_jetGenEmE"] = theDir.make<TH1D>("hist_pf_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
0322   hists["hist_pf_jetGenHadE"] = theDir.make<TH1D>("hist_pf_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
0323   hists["hist_pf_jetGenEMF"] = theDir.make<TH1D>("hist_pf_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
0324   hists["hist_pf_jetEoverGenE"] =
0325       theDir.make<TH1D>("hist_pf_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
0326   hists["hist_pf_jetCorr"] = theDir.make<TH1D>("hist_pf_jetCorr", "PFJet Correction Factor", 200, 0, 1.0);
0327   hists["hist_pf_nConstituents"] = theDir.make<TH1D>("hist_pf_nConstituents", "PFJet nConstituents", 20, 0, 20);
0328 
0329   hists["hist_pf_good_jetPt"] = theDir.make<TH1D>("hist_pf_good_jetPt", "PFJet p_{T}", 400, 0, 400);
0330   hists["hist_pf_good_jetEtaVsPhi"] = theDir.make<TH2D>(
0331       "hist_pf_good_jetEtaVsPhi", "PFJet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi());
0332   hists["hist_pf_good_jetNTracks"] = theDir.make<TH1D>("hist_pf_good_jetNTracks", "PFJet N_{TRACKS}", 20, 0, 20);
0333   hists["hist_pf_good_jetNTracksVsPt"] =
0334       theDir.make<TH2D>("hist_pf_good_jetNTracksVsPt",
0335                         "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",
0336                         20,
0337                         0,
0338                         200,
0339                         20,
0340                         0,
0341                         20);
0342   hists["hist_pf_good_jetEMF"] = theDir.make<TH1D>("hist_pf_good_jetEMF", "PFJet EMF", 200, 0, 1);
0343   hists["hist_pf_good_jetCHF"] = theDir.make<TH1D>("hist_pf_good_jetCHF", "PFJet CHF", 200, 0, 1);
0344   hists["hist_pf_good_jetNHF"] = theDir.make<TH1D>("hist_pf_good_jetNHF", "PFJet NHF", 200, 0, 1);
0345   hists["hist_pf_good_jetCEF"] = theDir.make<TH1D>("hist_pf_good_jetCEF", "PFJet CEF", 200, 0, 1);
0346   hists["hist_pf_good_jetNEF"] = theDir.make<TH1D>("hist_pf_good_jetNEF", "PFJet NEF", 200, 0, 1);
0347   hists["hist_pf_good_jetPdgID"] =
0348       theDir.make<TH1D>("hist_pf_good_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
0349   hists["hist_pf_good_jetGenEmE"] = theDir.make<TH1D>("hist_pf_good_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
0350   hists["hist_pf_good_jetGenHadE"] = theDir.make<TH1D>("hist_pf_good_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
0351   hists["hist_pf_good_jetGenEMF"] = theDir.make<TH1D>("hist_pf_good_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
0352   hists["hist_pf_good_jetEoverGenE"] =
0353       theDir.make<TH1D>("hist_pf_good_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
0354   hists["hist_pf_good_jetCorr"] = theDir.make<TH1D>("hist_pf_good_jetCorr", "PFJet Correction Factor", 200, 0, 1.0);
0355   hists["hist_pf_good_nConstituents"] =
0356       theDir.make<TH1D>("hist_pf_good_nConstituents", "PFJet nConstituents", 20, 0, 20);
0357 
0358   hists["hist_mjj"] = theDir.make<TH1D>("hist_mjj", "Dijet mass", 300, 0, 300);
0359   hists["hist_dR_jj"] = theDir.make<TH1D>("hist_dR_jj", "#Delta R_{JJ}", 200, 0, TMath::TwoPi());
0360   hists["hist_imbalance_jj"] = theDir.make<TH1D>("hist_imbalance_jj", "Dijet imbalance", 200, -1.0, 1.0);
0361 
0362   hists["hist_pf_mjj"] = theDir.make<TH1D>("hist_pf_mjj", "Dijet mass", 300, 0, 300);
0363   hists["hist_pf_dR_jj"] = theDir.make<TH1D>("hist_pf_dR_jj", "#Delta R_{JJ}", 200, 0, TMath::TwoPi());
0364   hists["hist_pf_imbalance_jj"] = theDir.make<TH1D>("hist_pf_imbalance_jj", "Dijet imbalance", 200, -1.0, 1.0);
0365 
0366   cout << "Making functors" << endl;
0367   JetIDStudiesSelector caloSelector(caloJetIDParameters, pfJetIDParameters, jetStudiesParams);
0368 
0369   JetIDStudiesSelector pfSelector(caloJetIDParameters, pfJetIDParameters, pfJetStudiesParams);
0370 
0371   bool doTracks = plotParameters.getParameter<bool>("doTracks");
0372   bool useMC = plotParameters.getParameter<bool>("useMC");
0373 
0374   cout << "About to begin looping" << endl;
0375 
0376   int nev = 0;
0377   //loop through each event
0378   for (ev.toBegin(); !ev.atEnd(); ++ev, ++nev) {
0379     edm::EventBase const& event = ev;
0380 
0381     if (runLumiSel(ev) == false)
0382       continue;
0383 
0384     if (nev % 100 == 0)
0385       cout << "Processing run " << event.id().run() << ", lumi " << event.id().luminosityBlock() << ", event "
0386            << event.id().event() << endl;
0387 
0388     pat::strbitset retCalo = caloSelector.getBitTemplate();
0389     caloSelector(event, retCalo);
0390 
0391     pat::strbitset retPF = pfSelector.getBitTemplate();
0392     pfSelector(event, retPF);
0393 
0394     ///------------------
0395     /// CALO JETS
0396     ///------------------
0397     if (retCalo.test(caloSelector.caloKin())) {
0398       if (retCalo.test(caloSelector.caloDeltaPhi())) {
0399         vector<pat::Jet> const& allCaloJets = caloSelector.allCaloJets();
0400 
0401         for (std::vector<pat::Jet>::const_iterator jetBegin = allCaloJets.begin(),
0402                                                    jetEnd = jetBegin + 2,
0403                                                    ijet = jetBegin;
0404              ijet != jetEnd;
0405              ++ijet) {
0406           const pat::Jet& jet = *ijet;
0407 
0408           double pt = jet.pt();
0409 
0410           const reco::TrackRefVector& jetTracks = jet.associatedTracks();
0411 
0412           hists["hist_jetPt"]->Fill(pt);
0413           hists["hist_jetEtaVsPhi"]->Fill(jet.eta(), jet.phi());
0414           hists["hist_jetNTracks"]->Fill(jetTracks.size());
0415           hists["hist_jetNTracksVsPt"]->Fill(pt, jetTracks.size());
0416           hists["hist_jetEMF"]->Fill(jet.emEnergyFraction());
0417           hists["hist_jetCorr"]->Fill(jet.jecFactor("Uncorrected"));
0418           hists["hist_n90Hits"]->Fill(static_cast<int>(jet.jetID().n90Hits));
0419           hists["hist_fHPD"]->Fill(jet.jetID().fHPD);
0420           hists["hist_nConstituents"]->Fill(jet.nConstituents());
0421 
0422           if (useMC && jet.genJet() != nullptr) {
0423             hists["hist_jetGenEmE"]->Fill(jet.genJet()->emEnergy());
0424             hists["hist_jetGenHadE"]->Fill(jet.genJet()->hadEnergy());
0425             hists["hist_jetEoverGenE"]->Fill(jet.energy() / jet.genJet()->energy());
0426             hists["hist_jetGenEMF"]->Fill(jet.genJet()->emEnergy() / jet.genJet()->energy());
0427           }
0428           if (doTracks) {
0429             TLorentzVector p4_tracks(0, 0, 0, 0);
0430             for (reco::TrackRefVector::const_iterator itrk = jetTracks.begin(), itrkEnd = jetTracks.end();
0431                  itrk != itrkEnd;
0432                  ++itrk) {
0433               TLorentzVector p4_trk;
0434               double M_PION = 0.140;
0435               p4_trk.SetPtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), M_PION);
0436               p4_tracks += p4_trk;
0437             }
0438             hists["hist_jetCHF"]->Fill(p4_tracks.Energy() / jet.energy());
0439           }
0440 
0441         }  // end loop over jets
0442 
0443         if (retCalo.test(caloSelector.caloJetID())) {
0444           pat::Jet const& jet0 = caloSelector.caloJet0();
0445           pat::Jet const& jet1 = caloSelector.caloJet1();
0446 
0447           TLorentzVector p4_j0(jet0.px(), jet0.py(), jet0.pz(), jet0.energy());
0448           TLorentzVector p4_j1(jet1.px(), jet1.py(), jet1.pz(), jet1.energy());
0449 
0450           TLorentzVector p4_jj = p4_j0 + p4_j1;
0451 
0452           hists["hist_mjj"]->Fill(p4_jj.M());
0453           hists["hist_dR_jj"]->Fill(p4_j0.DeltaR(p4_j1));
0454           hists["hist_imbalance_jj"]->Fill((p4_j0.Perp() - p4_j1.Perp()) / (p4_j0.Perp() + p4_j1.Perp()));
0455 
0456           hists["hist_good_jetPt"]->Fill(jet0.pt());
0457           hists["hist_good_jetEtaVsPhi"]->Fill(jet0.eta(), jet0.phi());
0458           hists["hist_good_jetNTracks"]->Fill(jet0.associatedTracks().size());
0459           hists["hist_good_jetNTracksVsPt"]->Fill(jet0.pt(), jet0.associatedTracks().size());
0460           hists["hist_good_jetEMF"]->Fill(jet0.emEnergyFraction());
0461           hists["hist_good_jetCorr"]->Fill(jet0.jecFactor("Uncorrected"));
0462           hists["hist_good_n90Hits"]->Fill(static_cast<int>(jet0.jetID().n90Hits));
0463           hists["hist_good_fHPD"]->Fill(jet0.jetID().fHPD);
0464           hists["hist_good_nConstituents"]->Fill(jet0.nConstituents());
0465 
0466           hists["hist_good_jetPt"]->Fill(jet1.pt());
0467           hists["hist_good_jetEtaVsPhi"]->Fill(jet1.eta(), jet1.phi());
0468           hists["hist_good_jetNTracks"]->Fill(jet1.associatedTracks().size());
0469           hists["hist_good_jetNTracksVsPt"]->Fill(jet1.pt(), jet1.associatedTracks().size());
0470           hists["hist_good_jetEMF"]->Fill(jet1.emEnergyFraction());
0471           hists["hist_good_jetCorr"]->Fill(jet1.jecFactor("Uncorrected"));
0472           hists["hist_good_n90Hits"]->Fill(static_cast<int>(jet1.jetID().n90Hits));
0473           hists["hist_good_fHPD"]->Fill(jet1.jetID().fHPD);
0474           hists["hist_good_nConstituents"]->Fill(jet1.nConstituents());
0475 
0476         }  // end if passed calo jet id
0477       }    // end if passed dphi cuts
0478     }      // end if passed kin cuts
0479 
0480     ///------------------
0481     /// PF JETS
0482     ///------------------
0483     if (retPF.test(pfSelector.pfDeltaPhi())) {
0484       vector<pat::Jet> const& allPFJets = pfSelector.allPFJets();
0485 
0486       for (std::vector<pat::Jet>::const_iterator jetBegin = allPFJets.begin(), jetEnd = jetBegin + 2, ijet = jetBegin;
0487            ijet != jetEnd;
0488            ++ijet) {
0489         const pat::Jet& jet = *ijet;
0490 
0491         double pt = jet.pt();
0492 
0493         hists["hist_pf_jetPt"]->Fill(pt);
0494         hists["hist_pf_jetEtaVsPhi"]->Fill(jet.eta(), jet.phi());
0495         hists["hist_pf_nConstituents"]->Fill(jet.nConstituents());
0496         hists["hist_pf_jetCEF"]->Fill(jet.chargedEmEnergyFraction());
0497         hists["hist_pf_jetNEF"]->Fill(jet.neutralEmEnergyFraction());
0498         hists["hist_pf_jetCHF"]->Fill(jet.chargedHadronEnergyFraction());
0499         hists["hist_pf_jetNHF"]->Fill(jet.neutralHadronEnergyFraction());
0500 
0501         if (useMC && jet.genJet() != nullptr) {
0502           hists["hist_pf_jetGenEmE"]->Fill(jet.genJet()->emEnergy());
0503           hists["hist_pf_jetGenHadE"]->Fill(jet.genJet()->hadEnergy());
0504           hists["hist_pf_jetEoverGenE"]->Fill(jet.energy() / jet.genJet()->energy());
0505 
0506           hists["hist_pf_jetGenEMF"]->Fill(jet.genJet()->emEnergy() / jet.genJet()->energy());
0507         }
0508 
0509       }  // end loop over jets
0510 
0511       if (retPF.test(pfSelector.pfJetID())) {
0512         pat::Jet const& jet0 = pfSelector.pfJet0();
0513         pat::Jet const& jet1 = pfSelector.pfJet1();
0514 
0515         TLorentzVector p4_j0(jet0.px(), jet0.py(), jet0.pz(), jet0.energy());
0516         TLorentzVector p4_j1(jet1.px(), jet1.py(), jet1.pz(), jet1.energy());
0517 
0518         TLorentzVector p4_jj = p4_j0 + p4_j1;
0519 
0520         hists["hist_pf_mjj"]->Fill(p4_jj.M());
0521         hists["hist_pf_dR_jj"]->Fill(p4_j0.DeltaR(p4_j1));
0522         hists["hist_pf_imbalance_jj"]->Fill((p4_j0.Perp() - p4_j1.Perp()) / (p4_j0.Perp() + p4_j1.Perp()));
0523 
0524         hists["hist_pf_good_jetPt"]->Fill(jet0.pt());
0525         hists["hist_pf_good_jetEtaVsPhi"]->Fill(jet0.eta(), jet0.phi());
0526         hists["hist_pf_good_nConstituents"]->Fill(jet0.nConstituents());
0527         hists["hist_pf_good_jetCEF"]->Fill(jet0.chargedEmEnergyFraction());
0528         hists["hist_pf_good_jetNEF"]->Fill(jet0.neutralEmEnergyFraction());
0529         hists["hist_pf_good_jetCHF"]->Fill(jet0.chargedHadronEnergyFraction());
0530         hists["hist_pf_good_jetNHF"]->Fill(jet0.neutralHadronEnergyFraction());
0531 
0532         hists["hist_pf_good_jetPt"]->Fill(jet1.pt());
0533         hists["hist_pf_good_jetEtaVsPhi"]->Fill(jet1.eta(), jet1.phi());
0534         hists["hist_pf_good_nConstituents"]->Fill(jet1.nConstituents());
0535         hists["hist_pf_good_jetCEF"]->Fill(jet1.chargedEmEnergyFraction());
0536         hists["hist_pf_good_jetNEF"]->Fill(jet1.neutralEmEnergyFraction());
0537         hists["hist_pf_good_jetCHF"]->Fill(jet1.chargedHadronEnergyFraction());
0538         hists["hist_pf_good_jetNHF"]->Fill(jet1.neutralHadronEnergyFraction());
0539 
0540       }  // end if 2 good PF jets
0541 
0542     }  // end if delta phi pf cuts
0543 
0544   }  // end loop over events
0545 
0546   cout << "Calo jet selection" << endl;
0547   caloSelector.print(std::cout);
0548   cout << "PF jet selection" << endl;
0549   pfSelector.print(std::cout);
0550 
0551   return 0;
0552 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!