PatBasicFWLiteJetAnalyzer_Selector.cc

Go to the documentation of this file.

/*   A macro for making a histogram of Jet Pt with cuts
This is a basic way to cut out jets of a certain Pt and Eta using an if statement
This example creates a histogram of Jet Pt, using Jets with Pt above 30 and ETA above -2.1 and below 2.1
*/

#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TSystem.h"
#include "TLorentzVector.h"

#include "PhysicsTools/SelectorUtils/interface/strbitset.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/MET.h"
#include "DataFormats/PatCandidates/interface/Photon.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Tau.h"
#include "PhysicsTools/SelectorUtils/interface/Selector.h"
#include "PhysicsTools/SelectorUtils/interface/JetIDSelectionFunctor.h"
#include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
#include "PhysicsTools/SelectorUtils/interface/RunLumiSelector.h"
#include "FWCore/FWLite/interface/FWLiteEnabler.h"
#include "PhysicsTools/FWLite/interface/TFileService.h"
#include "DataFormats/FWLite/interface/ChainEvent.h"
#include "FWCore/ParameterSetReader/interface/ProcessDescImpl.h"
#include "FWCore/ParameterSet/interface/ProcessDesc.h"

#include <iostream>
#include <memory>
#include <cmath>  //necessary for absolute function fabs()

using namespace std;


// -*- C++ -*-

// CMS includes
#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/Ptr.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "PhysicsTools/SelectorUtils/interface/EventSelector.h"

// Root includes
#include "TROOT.h"

using namespace std;

class JetIDStudiesSelector : public EventSelector {
public:
  JetIDStudiesSelector(edm::ParameterSet const& caloJetIdParams,
                       edm::ParameterSet const& pfJetIdParams,
                       edm::ParameterSet const& params)
      : jetSel_(new JetIDSelectionFunctor(caloJetIdParams)),
        pfJetSel_(new PFJetIDSelectionFunctor(pfJetIdParams)),
        jetSrc_(params.getParameter<edm::InputTag>("jetSrc")),
        pfJetSrc_(params.getParameter<edm::InputTag>("pfJetSrc")) {
    bool useCalo = params.getParameter<bool>("useCalo");

    push_back("Calo Cuts");
    push_back("Calo Kin Cuts");
    push_back("Calo Delta Phi");
    push_back("Calo Jet ID");
    push_back("PF Cuts");
    push_back("PF Kin Cuts");
    push_back("PF Delta Phi");
    push_back("PF Jet ID");

    set("Calo Cuts", useCalo);
    set("Calo Kin Cuts", useCalo);
    set("Calo Delta Phi", useCalo);
    set("Calo Jet ID", useCalo);
    set("PF Cuts", !useCalo);
    set("PF Kin Cuts", !useCalo);
    set("PF Delta Phi", !useCalo);
    set("PF Jet ID", !useCalo);

    // Indices for fast caching
    caloCuts_ = index_type(&bits_, std::string("Calo Cuts"));
    caloKin_ = index_type(&bits_, std::string("Calo Kin Cuts"));
    caloDeltaPhi_ = index_type(&bits_, std::string("Calo Delta Phi"));
    caloJetID_ = index_type(&bits_, std::string("Calo Jet ID"));
    pfCuts_ = index_type(&bits_, std::string("PF Cuts"));
    pfKin_ = index_type(&bits_, std::string("PF Kin Cuts"));
    pfDeltaPhi_ = index_type(&bits_, std::string("PF Delta Phi"));
    pfJetID_ = index_type(&bits_, std::string("PF Jet ID"));
  }

  ~JetIDStudiesSelector() override {}

  bool operator()(edm::EventBase const& event, pat::strbitset& ret) override {
    pat::strbitset retCaloJet = jetSel_->getBitTemplate();
    pat::strbitset retPFJet = pfJetSel_->getBitTemplate();

    if (considerCut(caloCuts_)) {
      passCut(ret, caloCuts_);
      event.getByLabel(jetSrc_, h_jets_);
      // Calo Cuts
      if (h_jets_->size() >= 2 || ignoreCut(caloKin_)) {
        passCut(ret, caloKin_);
        pat::Jet const& jet0 = h_jets_->at(0);
        pat::Jet const& jet1 = h_jets_->at(1);
        double dphi = fabs(deltaPhi<double>(jet0.phi(), jet1.phi()));

        if (fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(caloDeltaPhi_)) {
          passCut(ret, caloDeltaPhi_);

          retCaloJet.set(false);
          bool pass0 = (*jetSel_)(jet0, retCaloJet);
          retCaloJet.set(false);
          bool pass1 = (*jetSel_)(jet1, retCaloJet);
          if ((pass0 && pass1) || ignoreCut(caloJetID_)) {
            passCut(ret, caloJetID_);
            caloJet0_ = edm::Ptr<pat::Jet>(h_jets_, 0);
            caloJet1_ = edm::Ptr<pat::Jet>(h_jets_, 1);

            return (bool)ret;
          }  // end if found 2 "loose" jet ID jets
        }    // end if delta phi
      }      // end calo kin cuts
    }        // end if calo cuts

    if (considerCut(pfCuts_)) {
      passCut(ret, pfCuts_);
      event.getByLabel(pfJetSrc_, h_pfjets_);
      // PF Cuts
      if (h_pfjets_->size() >= 2 || ignoreCut(pfKin_)) {
        passCut(ret, pfKin_);
        pat::Jet const& jet0 = h_pfjets_->at(0);
        pat::Jet const& jet1 = h_pfjets_->at(1);
        double dphi = fabs(deltaPhi<double>(jet0.phi(), jet1.phi()));

        if (fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(pfDeltaPhi_)) {
          passCut(ret, pfDeltaPhi_);

          retPFJet.set(false);
          bool pass0 = (*pfJetSel_)(jet0, retPFJet);
          retPFJet.set(false);
          bool pass1 = (*pfJetSel_)(jet1, retPFJet);
          if ((pass0 && pass1) || ignoreCut(pfJetID_)) {
            passCut(ret, pfJetID_);
            pfJet0_ = edm::Ptr<pat::Jet>(h_pfjets_, 0);
            pfJet1_ = edm::Ptr<pat::Jet>(h_pfjets_, 1);

            return (bool)ret;
          }  // end if found 2 "loose" jet ID jets
        }    // end if delta phi
      }      // end pf kin cuts
    }        // end if pf cuts

    setIgnored(ret);

    return false;
  }  // end of method

  std::shared_ptr<JetIDSelectionFunctor> const& jetSel() const { return jetSel_; }
  std::shared_ptr<PFJetIDSelectionFunctor> const& pfJetSel() const { return pfJetSel_; }

  vector<pat::Jet> const& allCaloJets() const { return *h_jets_; }
  vector<pat::Jet> const& allPFJets() const { return *h_pfjets_; }

  pat::Jet const& caloJet0() const { return *caloJet0_; }
  pat::Jet const& caloJet1() const { return *caloJet1_; }

  pat::Jet const& pfJet0() const { return *pfJet0_; }
  pat::Jet const& pfJet1() const { return *pfJet1_; }

  // Fast caching indices
  index_type const& caloCuts() const { return caloCuts_; }
  index_type const& caloKin() const { return caloKin_; }
  index_type const& caloDeltaPhi() const { return caloDeltaPhi_; }
  index_type const& caloJetID() const { return caloJetID_; }
  index_type const& pfCuts() const { return pfCuts_; }
  index_type const& pfKin() const { return pfKin_; }
  index_type const& pfDeltaPhi() const { return pfDeltaPhi_; }
  index_type const& pfJetID() const { return pfJetID_; }

protected:
  std::shared_ptr<JetIDSelectionFunctor> jetSel_;
  std::shared_ptr<PFJetIDSelectionFunctor> pfJetSel_;
  edm::InputTag jetSrc_;
  edm::InputTag pfJetSrc_;

  edm::Handle<vector<pat::Jet> > h_jets_;
  edm::Handle<vector<pat::Jet> > h_pfjets_;

  edm::Ptr<pat::Jet> caloJet0_;
  edm::Ptr<pat::Jet> caloJet1_;

  edm::Ptr<pat::Jet> pfJet0_;
  edm::Ptr<pat::Jet> pfJet1_;

  // Fast caching indices
  index_type caloCuts_;
  index_type caloKin_;
  index_type caloDeltaPhi_;
  index_type caloJetID_;
  index_type pfCuts_;
  index_type pfKin_;
  index_type pfDeltaPhi_;
  index_type pfJetID_;
};

// ///////////////////// //
// // Main Subroutine // //
// ///////////////////// //

int main(int argc, char* argv[]) {
  if (argc < 2) {
    std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
    return 0;
  }

  cout << "Hello from " << argv[0] << "!" << endl;

  // load framework libraries
  gSystem->Load("libFWCoreFWLite");
  FWLiteEnabler::enable();

  cout << "Getting parameters" << endl;
  // Get the python configuration
  ProcessDescImpl builder(argv[1]);
  std::shared_ptr<edm::ProcessDesc> b = builder.processDesc();
  std::shared_ptr<edm::ParameterSet> parameters = b->getProcessPSet();
  parameters->registerIt();

  edm::ParameterSet const& jetStudiesParams = parameters->getParameter<edm::ParameterSet>("jetStudies");
  edm::ParameterSet const& pfJetStudiesParams = parameters->getParameter<edm::ParameterSet>("pfJetStudies");
  edm::ParameterSet const& caloJetIDParameters = parameters->getParameter<edm::ParameterSet>("jetIDSelector");
  edm::ParameterSet const& pfJetIDParameters = parameters->getParameter<edm::ParameterSet>("pfJetIDSelector");
  edm::ParameterSet const& plotParameters = parameters->getParameter<edm::ParameterSet>("plotParameters");
  edm::ParameterSet const& inputs = parameters->getParameter<edm::ParameterSet>("inputs");
  edm::ParameterSet const& outputs = parameters->getParameter<edm::ParameterSet>("outputs");

  cout << "Making RunLumiSelector" << endl;
  RunLumiSelector runLumiSel(inputs);

  cout << "setting up TFileService" << endl;
  // book a set of histograms
  fwlite::TFileService fs = fwlite::TFileService(outputs.getParameter<std::string>("outputName"));
  TFileDirectory theDir = fs.mkdir("histos");

  cout << "Setting up chain event" << endl;
  // This object 'event' is used both to get all information from the
  // event as well as to store histograms, etc.
  fwlite::ChainEvent ev(inputs.getParameter<std::vector<std::string> >("fileNames"));

  cout << "Booking histograms" << endl;
  // Book histograms

  std::map<std::string, TH1*> hists;

  hists["hist_nJet"] = theDir.make<TH1D>("hist_nJet", "Number of Calo Jets", 10, 0, 10);
  hists["hist_nPFJet"] = theDir.make<TH1D>("hist_nPFJet", "Number of PF Jets", 10, 0, 10);

  hists["hist_jetPt"] = theDir.make<TH1D>("hist_jetPt", "Jet p_{T}", 400, 0, 400);
  hists["hist_jetEtaVsPhi"] = theDir.make<TH2D>(
      "hist_jetEtaVsPhi", "Jet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi());
  hists["hist_jetNTracks"] = theDir.make<TH1D>("hist_jetNTracks", "Jet N_{TRACKS}", 20, 0, 20);
  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);
  hists["hist_jetEMF"] = theDir.make<TH1D>("hist_jetEMF", "Jet EMF", 200, 0, 1);
  hists["hist_jetPdgID"] = theDir.make<TH1D>("hist_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
  hists["hist_jetGenEmE"] = theDir.make<TH1D>("hist_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
  hists["hist_jetGenHadE"] = theDir.make<TH1D>("hist_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
  hists["hist_jetGenEMF"] = theDir.make<TH1D>("hist_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
  hists["hist_jetEoverGenE"] =
      theDir.make<TH1D>("hist_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
  hists["hist_jetCorr"] = theDir.make<TH1D>("hist_jetCorr", "Jet Correction Factor", 200, 0, 1.0);
  hists["hist_n90Hits"] = theDir.make<TH1D>("hist_n90Hits", "Jet n90Hits", 20, 0, 20);
  hists["hist_fHPD"] = theDir.make<TH1D>("hist_fHPD", "Jet fHPD", 200, 0, 1);
  hists["hist_nConstituents"] = theDir.make<TH1D>("hist_nConstituents", "Jet nConstituents", 20, 0, 20);
  hists["hist_jetCHF"] = theDir.make<TH1D>("hist_jetCHF", "Jet Charged Hadron Fraction", 200, 0, 1.0);

  hists["hist_good_jetPt"] = theDir.make<TH1D>("hist_good_jetPt", "Jet p_{T}", 400, 0, 400);
  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());
  hists["hist_good_jetNTracks"] = theDir.make<TH1D>("hist_good_jetNTracks", "Jet N_{TRACKS}", 20, 0, 20);
  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);
  hists["hist_good_jetEMF"] = theDir.make<TH1D>("hist_good_jetEMF", "Jet EMF", 200, 0, 1);
  hists["hist_good_jetPdgID"] = theDir.make<TH1D>("hist_good_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
  hists["hist_good_jetGenEmE"] = theDir.make<TH1D>("hist_good_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
  hists["hist_good_jetGenHadE"] = theDir.make<TH1D>("hist_good_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
  hists["hist_good_jetGenEMF"] = theDir.make<TH1D>("hist_good_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
  hists["hist_good_jetEoverGenE"] =
      theDir.make<TH1D>("hist_good_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
  hists["hist_good_jetCorr"] = theDir.make<TH1D>("hist_good_jetCorr", "Jet Correction Factor", 200, 0, 1.0);
  hists["hist_good_n90Hits"] = theDir.make<TH1D>("hist_good_n90Hits", "Jet n90Hits", 20, 0, 20);
  hists["hist_good_fHPD"] = theDir.make<TH1D>("hist_good_fHPD", "Jet fHPD", 200, 0, 1);
  hists["hist_good_nConstituents"] = theDir.make<TH1D>("hist_good_nConstituents", "Jet nConstituents", 20, 0, 20);
  hists["hist_good_jetCHF"] = theDir.make<TH1D>("hist_good_jetCHF", "Jet Charged Hadron Fraction", 200, 0, 1.0);

  hists["hist_pf_jetPt"] = theDir.make<TH1D>("hist_pf_jetPt", "PFJet p_{T}", 400, 0, 400);
  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());
  hists["hist_pf_jetNTracks"] = theDir.make<TH1D>("hist_pf_jetNTracks", "PFJet N_{TRACKS}", 20, 0, 20);
  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);
  hists["hist_pf_jetEMF"] = theDir.make<TH1D>("hist_pf_jetEMF", "PFJet EMF", 200, 0, 1);
  hists["hist_pf_jetCHF"] = theDir.make<TH1D>("hist_pf_jetCHF", "PFJet CHF", 200, 0, 1);
  hists["hist_pf_jetNHF"] = theDir.make<TH1D>("hist_pf_jetNHF", "PFJet NHF", 200, 0, 1);
  hists["hist_pf_jetCEF"] = theDir.make<TH1D>("hist_pf_jetCEF", "PFJet CEF", 200, 0, 1);
  hists["hist_pf_jetNEF"] = theDir.make<TH1D>("hist_pf_jetNEF", "PFJet NEF", 200, 0, 1);
  hists["hist_pf_jetPdgID"] = theDir.make<TH1D>("hist_pf_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
  hists["hist_pf_jetGenEmE"] = theDir.make<TH1D>("hist_pf_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
  hists["hist_pf_jetGenHadE"] = theDir.make<TH1D>("hist_pf_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
  hists["hist_pf_jetGenEMF"] = theDir.make<TH1D>("hist_pf_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
  hists["hist_pf_jetEoverGenE"] =
      theDir.make<TH1D>("hist_pf_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
  hists["hist_pf_jetCorr"] = theDir.make<TH1D>("hist_pf_jetCorr", "PFJet Correction Factor", 200, 0, 1.0);
  hists["hist_pf_nConstituents"] = theDir.make<TH1D>("hist_pf_nConstituents", "PFJet nConstituents", 20, 0, 20);

  hists["hist_pf_good_jetPt"] = theDir.make<TH1D>("hist_pf_good_jetPt", "PFJet p_{T}", 400, 0, 400);
  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());
  hists["hist_pf_good_jetNTracks"] = theDir.make<TH1D>("hist_pf_good_jetNTracks", "PFJet N_{TRACKS}", 20, 0, 20);
  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);
  hists["hist_pf_good_jetEMF"] = theDir.make<TH1D>("hist_pf_good_jetEMF", "PFJet EMF", 200, 0, 1);
  hists["hist_pf_good_jetCHF"] = theDir.make<TH1D>("hist_pf_good_jetCHF", "PFJet CHF", 200, 0, 1);
  hists["hist_pf_good_jetNHF"] = theDir.make<TH1D>("hist_pf_good_jetNHF", "PFJet NHF", 200, 0, 1);
  hists["hist_pf_good_jetCEF"] = theDir.make<TH1D>("hist_pf_good_jetCEF", "PFJet CEF", 200, 0, 1);
  hists["hist_pf_good_jetNEF"] = theDir.make<TH1D>("hist_pf_good_jetNEF", "PFJet NEF", 200, 0, 1);
  hists["hist_pf_good_jetPdgID"] =
      theDir.make<TH1D>("hist_pf_good_jetPdgID", "PDG Id of Jet Constituents", 10000, 0, 10000);
  hists["hist_pf_good_jetGenEmE"] = theDir.make<TH1D>("hist_pf_good_jetGenEmE", "Gen Jet EM Energy", 200, 0, 200);
  hists["hist_pf_good_jetGenHadE"] = theDir.make<TH1D>("hist_pf_good_jetGenHadE", "Gen Jet HAD Energy", 200, 0, 200);
  hists["hist_pf_good_jetGenEMF"] = theDir.make<TH1D>("hist_pf_good_jetGenEMF", "Gen Jet EMF", 200, 0, 1);
  hists["hist_pf_good_jetEoverGenE"] =
      theDir.make<TH1D>("hist_pf_good_jetEoverGenE", "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0);
  hists["hist_pf_good_jetCorr"] = theDir.make<TH1D>("hist_pf_good_jetCorr", "PFJet Correction Factor", 200, 0, 1.0);
  hists["hist_pf_good_nConstituents"] =
      theDir.make<TH1D>("hist_pf_good_nConstituents", "PFJet nConstituents", 20, 0, 20);

  hists["hist_mjj"] = theDir.make<TH1D>("hist_mjj", "Dijet mass", 300, 0, 300);
  hists["hist_dR_jj"] = theDir.make<TH1D>("hist_dR_jj", "#Delta R_{JJ}", 200, 0, TMath::TwoPi());
  hists["hist_imbalance_jj"] = theDir.make<TH1D>("hist_imbalance_jj", "Dijet imbalance", 200, -1.0, 1.0);

  hists["hist_pf_mjj"] = theDir.make<TH1D>("hist_pf_mjj", "Dijet mass", 300, 0, 300);
  hists["hist_pf_dR_jj"] = theDir.make<TH1D>("hist_pf_dR_jj", "#Delta R_{JJ}", 200, 0, TMath::TwoPi());
  hists["hist_pf_imbalance_jj"] = theDir.make<TH1D>("hist_pf_imbalance_jj", "Dijet imbalance", 200, -1.0, 1.0);

  cout << "Making functors" << endl;
  JetIDStudiesSelector caloSelector(caloJetIDParameters, pfJetIDParameters, jetStudiesParams);

  JetIDStudiesSelector pfSelector(caloJetIDParameters, pfJetIDParameters, pfJetStudiesParams);

  bool doTracks = plotParameters.getParameter<bool>("doTracks");
  bool useMC = plotParameters.getParameter<bool>("useMC");

  cout << "About to begin looping" << endl;

  int nev = 0;
  //loop through each event
  for (ev.toBegin(); !ev.atEnd(); ++ev, ++nev) {
    edm::EventBase const& event = ev;

    if (runLumiSel(ev) == false)
      continue;

    if (nev % 100 == 0)
      cout << "Processing run " << event.id().run() << ", lumi " << event.id().luminosityBlock() << ", event "
           << event.id().event() << endl;

    pat::strbitset retCalo = caloSelector.getBitTemplate();
    caloSelector(event, retCalo);

    pat::strbitset retPF = pfSelector.getBitTemplate();
    pfSelector(event, retPF);

    if (retCalo.test(caloSelector.caloKin())) {
      if (retCalo.test(caloSelector.caloDeltaPhi())) {
        vector<pat::Jet> const& allCaloJets = caloSelector.allCaloJets();

        for (std::vector<pat::Jet>::const_iterator jetBegin = allCaloJets.begin(),
                                                   jetEnd = jetBegin + 2,
                                                   ijet = jetBegin;
             ijet != jetEnd;
             ++ijet) {
          const pat::Jet& jet = *ijet;

          double pt = jet.pt();

          const reco::TrackRefVector& jetTracks = jet.associatedTracks();

          hists["hist_jetPt"]->Fill(pt);
          hists["hist_jetEtaVsPhi"]->Fill(jet.eta(), jet.phi());
          hists["hist_jetNTracks"]->Fill(jetTracks.size());
          hists["hist_jetNTracksVsPt"]->Fill(pt, jetTracks.size());
          hists["hist_jetEMF"]->Fill(jet.emEnergyFraction());
          hists["hist_jetCorr"]->Fill(jet.jecFactor("Uncorrected"));
          hists["hist_n90Hits"]->Fill(static_cast<int>(jet.jetID().n90Hits));
          hists["hist_fHPD"]->Fill(jet.jetID().fHPD);
          hists["hist_nConstituents"]->Fill(jet.nConstituents());

          if (useMC && jet.genJet() != nullptr) {
            hists["hist_jetGenEmE"]->Fill(jet.genJet()->emEnergy());
            hists["hist_jetGenHadE"]->Fill(jet.genJet()->hadEnergy());
            hists["hist_jetEoverGenE"]->Fill(jet.energy() / jet.genJet()->energy());
            hists["hist_jetGenEMF"]->Fill(jet.genJet()->emEnergy() / jet.genJet()->energy());
          }
          if (doTracks) {
            TLorentzVector p4_tracks(0, 0, 0, 0);
            for (reco::TrackRefVector::const_iterator itrk = jetTracks.begin(), itrkEnd = jetTracks.end();
                 itrk != itrkEnd;
                 ++itrk) {
              TLorentzVector p4_trk;
              double M_PION = 0.140;
              p4_trk.SetPtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), M_PION);
              p4_tracks += p4_trk;
            }
            hists["hist_jetCHF"]->Fill(p4_tracks.Energy() / jet.energy());
          }

        }  // end loop over jets

        if (retCalo.test(caloSelector.caloJetID())) {
          pat::Jet const& jet0 = caloSelector.caloJet0();
          pat::Jet const& jet1 = caloSelector.caloJet1();

          TLorentzVector p4_j0(jet0.px(), jet0.py(), jet0.pz(), jet0.energy());
          TLorentzVector p4_j1(jet1.px(), jet1.py(), jet1.pz(), jet1.energy());

          TLorentzVector p4_jj = p4_j0 + p4_j1;

          hists["hist_mjj"]->Fill(p4_jj.M());
          hists["hist_dR_jj"]->Fill(p4_j0.DeltaR(p4_j1));
          hists["hist_imbalance_jj"]->Fill((p4_j0.Perp() - p4_j1.Perp()) / (p4_j0.Perp() + p4_j1.Perp()));

          hists["hist_good_jetPt"]->Fill(jet0.pt());
          hists["hist_good_jetEtaVsPhi"]->Fill(jet0.eta(), jet0.phi());
          hists["hist_good_jetNTracks"]->Fill(jet0.associatedTracks().size());
          hists["hist_good_jetNTracksVsPt"]->Fill(jet0.pt(), jet0.associatedTracks().size());
          hists["hist_good_jetEMF"]->Fill(jet0.emEnergyFraction());
          hists["hist_good_jetCorr"]->Fill(jet0.jecFactor("Uncorrected"));
          hists["hist_good_n90Hits"]->Fill(static_cast<int>(jet0.jetID().n90Hits));
          hists["hist_good_fHPD"]->Fill(jet0.jetID().fHPD);
          hists["hist_good_nConstituents"]->Fill(jet0.nConstituents());

          hists["hist_good_jetPt"]->Fill(jet1.pt());
          hists["hist_good_jetEtaVsPhi"]->Fill(jet1.eta(), jet1.phi());
          hists["hist_good_jetNTracks"]->Fill(jet1.associatedTracks().size());
          hists["hist_good_jetNTracksVsPt"]->Fill(jet1.pt(), jet1.associatedTracks().size());
          hists["hist_good_jetEMF"]->Fill(jet1.emEnergyFraction());
          hists["hist_good_jetCorr"]->Fill(jet1.jecFactor("Uncorrected"));
          hists["hist_good_n90Hits"]->Fill(static_cast<int>(jet1.jetID().n90Hits));
          hists["hist_good_fHPD"]->Fill(jet1.jetID().fHPD);
          hists["hist_good_nConstituents"]->Fill(jet1.nConstituents());

        }  // end if passed calo jet id
      }    // end if passed dphi cuts
    }      // end if passed kin cuts

    if (retPF.test(pfSelector.pfDeltaPhi())) {
      vector<pat::Jet> const& allPFJets = pfSelector.allPFJets();

      for (std::vector<pat::Jet>::const_iterator jetBegin = allPFJets.begin(), jetEnd = jetBegin + 2, ijet = jetBegin;
           ijet != jetEnd;
           ++ijet) {
        const pat::Jet& jet = *ijet;

        double pt = jet.pt();

        hists["hist_pf_jetPt"]->Fill(pt);
        hists["hist_pf_jetEtaVsPhi"]->Fill(jet.eta(), jet.phi());
        hists["hist_pf_nConstituents"]->Fill(jet.nConstituents());
        hists["hist_pf_jetCEF"]->Fill(jet.chargedEmEnergyFraction());
        hists["hist_pf_jetNEF"]->Fill(jet.neutralEmEnergyFraction());
        hists["hist_pf_jetCHF"]->Fill(jet.chargedHadronEnergyFraction());
        hists["hist_pf_jetNHF"]->Fill(jet.neutralHadronEnergyFraction());

        if (useMC && jet.genJet() != nullptr) {
          hists["hist_pf_jetGenEmE"]->Fill(jet.genJet()->emEnergy());
          hists["hist_pf_jetGenHadE"]->Fill(jet.genJet()->hadEnergy());
          hists["hist_pf_jetEoverGenE"]->Fill(jet.energy() / jet.genJet()->energy());

          hists["hist_pf_jetGenEMF"]->Fill(jet.genJet()->emEnergy() / jet.genJet()->energy());
        }

      }  // end loop over jets

      if (retPF.test(pfSelector.pfJetID())) {
        pat::Jet const& jet0 = pfSelector.pfJet0();
        pat::Jet const& jet1 = pfSelector.pfJet1();

        TLorentzVector p4_j0(jet0.px(), jet0.py(), jet0.pz(), jet0.energy());
        TLorentzVector p4_j1(jet1.px(), jet1.py(), jet1.pz(), jet1.energy());

        TLorentzVector p4_jj = p4_j0 + p4_j1;

        hists["hist_pf_mjj"]->Fill(p4_jj.M());
        hists["hist_pf_dR_jj"]->Fill(p4_j0.DeltaR(p4_j1));
        hists["hist_pf_imbalance_jj"]->Fill((p4_j0.Perp() - p4_j1.Perp()) / (p4_j0.Perp() + p4_j1.Perp()));

        hists["hist_pf_good_jetPt"]->Fill(jet0.pt());
        hists["hist_pf_good_jetEtaVsPhi"]->Fill(jet0.eta(), jet0.phi());
        hists["hist_pf_good_nConstituents"]->Fill(jet0.nConstituents());
        hists["hist_pf_good_jetCEF"]->Fill(jet0.chargedEmEnergyFraction());
        hists["hist_pf_good_jetNEF"]->Fill(jet0.neutralEmEnergyFraction());
        hists["hist_pf_good_jetCHF"]->Fill(jet0.chargedHadronEnergyFraction());
        hists["hist_pf_good_jetNHF"]->Fill(jet0.neutralHadronEnergyFraction());

        hists["hist_pf_good_jetPt"]->Fill(jet1.pt());
        hists["hist_pf_good_jetEtaVsPhi"]->Fill(jet1.eta(), jet1.phi());
        hists["hist_pf_good_nConstituents"]->Fill(jet1.nConstituents());
        hists["hist_pf_good_jetCEF"]->Fill(jet1.chargedEmEnergyFraction());
        hists["hist_pf_good_jetNEF"]->Fill(jet1.neutralEmEnergyFraction());
        hists["hist_pf_good_jetCHF"]->Fill(jet1.chargedHadronEnergyFraction());
        hists["hist_pf_good_jetNHF"]->Fill(jet1.neutralHadronEnergyFraction());

      }  // end if 2 good PF jets

    }  // end if delta phi pf cuts

  }  // end loop over events

  cout << "Calo jet selection" << endl;
  caloSelector.print(std::cout);
  cout << "PF jet selection" << endl;
  pfSelector.print(std::cout);

  return 0;
}