JetAnalyzer_HeavyIons_matching.cc

Go to the documentation of this file.

//
// Jet Analyzer class for heavy ion jets. for DQM jet analysis monitoring
// For CMSSW_7_5_X, especially reading background subtracted jets - jet matching 1 - 1
// author: Raghav Kunnawalkam Elayavalli,
//         April 6th 2015
//         Rutgers University, email: raghav.k.e at CERN dot CH
//

#include "DQMOffline/JetMET/interface/JetAnalyzer_HeavyIons_matching.h"

using namespace edm;
using namespace reco;
using namespace std;

// declare the constructors:

JetAnalyzer_HeavyIons_matching::JetAnalyzer_HeavyIons_matching(const edm::ParameterSet& iConfig)
    : mInputJet1Collection(iConfig.getParameter<edm::InputTag>("src_Jet1")),
      mInputJet2Collection(iConfig.getParameter<edm::InputTag>("src_Jet2")),
      JetType1(iConfig.getUntrackedParameter<std::string>("Jet1")),
      JetType2(iConfig.getUntrackedParameter<std::string>("Jet2")),
      mRecoJetPtThreshold(iConfig.getParameter<double>("recoJetPtThreshold")),
      mRecoDelRMatch(iConfig.getParameter<double>("recoDelRMatch")),
      mRecoJetEtaCut(iConfig.getParameter<double>("recoJetEtaCut")) {
  std::string inputCollectionLabelJet1(mInputJet1Collection.label());
  std::string inputCollectionLabelJet2(mInputJet2Collection.label());

  //consumes

  if (std::string("VsCalo") == JetType1)
    caloJet1Token_ = consumes<reco::CaloJetCollection>(mInputJet1Collection);
  if (std::string("VsPF") == JetType1)
    pfJetsToken_ = consumes<reco::PFJetCollection>(mInputJet1Collection);
  if (std::string("PuCalo") == JetType1)
    caloJet2Token_ = consumes<reco::CaloJetCollection>(mInputJet1Collection);
  if (std::string("PuPF") == JetType1)
    basicJetsToken_ = consumes<reco::BasicJetCollection>(mInputJet1Collection);

  if (std::string("VsCalo") == JetType2)
    caloJet1Token_ = consumes<reco::CaloJetCollection>(mInputJet2Collection);
  if (std::string("VsPF") == JetType2)
    pfJetsToken_ = consumes<reco::PFJetCollection>(mInputJet2Collection);
  if (std::string("PuCalo") == JetType2)
    caloJet2Token_ = consumes<reco::CaloJetCollection>(mInputJet2Collection);
  if (std::string("PuPF") == JetType2)
    basicJetsToken_ = consumes<reco::BasicJetCollection>(mInputJet2Collection);

  // initialize the Jet matching histograms

  mpT_ratio_Jet1Jet2 = nullptr;
  mpT_Jet1_matched = nullptr;
  mpT_Jet2_matched = nullptr;
  mpT_Jet1_unmatched = nullptr;
  mpT_Jet2_unmatched = nullptr;

  // we need to add histograms which will hold the hadronic and electromagnetic energy content for the unmatched Jets.
  if (std::string("VsCalo") == JetType1 || std::string("PuCalo") == JetType1) {
    mHadEnergy_Jet1_unmatched = nullptr;
    mEmEnergy_Jet1_unmatched = nullptr;
  }
  if (std::string("VsCalo") == JetType2 || std::string("PuCalo") == JetType2) {
    mHadEnergy_Jet2_unmatched = nullptr;
    mEmEnergy_Jet2_unmatched = nullptr;
  }

  if (std::string("VsPF") == JetType1) {
    mChargedHadronEnergy_Jet1_unmatched = nullptr;
    mNeutralHadronEnergy_Jet1_unmatched = nullptr;
    mChargedEmEnergy_Jet1_unmatched = nullptr;
    mNeutralEmEnergy_Jet1_unmatched = nullptr;
    mChargedMuEnergy_Jet1_unmatched = nullptr;
    mChargedHadEnergyFraction_Jet1_unmatched = nullptr;
    mNeutralHadEnergyFraction_Jet1_unmatched = nullptr;
    mPhotonEnergyFraction_Jet1_unmatched = nullptr;
    mElectronEnergyFraction_Jet1_unmatched = nullptr;
    mMuonEnergyFraction_Jet1_unmatched = nullptr;
  }

  if (std::string("VsPF") == JetType2) {
    mChargedHadronEnergy_Jet2_unmatched = nullptr;
    mNeutralHadronEnergy_Jet2_unmatched = nullptr;
    mChargedEmEnergy_Jet2_unmatched = nullptr;
    mNeutralEmEnergy_Jet2_unmatched = nullptr;
    mChargedMuEnergy_Jet2_unmatched = nullptr;
    mChargedHadEnergyFraction_Jet2_unmatched = nullptr;
    mNeutralHadEnergyFraction_Jet2_unmatched = nullptr;
    mPhotonEnergyFraction_Jet2_unmatched = nullptr;
    mElectronEnergyFraction_Jet2_unmatched = nullptr;
    mMuonEnergyFraction_Jet2_unmatched = nullptr;
  }
}

void JetAnalyzer_HeavyIons_matching::bookHistograms(DQMStore::IBooker& ibooker,
                                                    edm::Run const& iRun,
                                                    edm::EventSetup const&) {
  ibooker.setCurrentFolder("JetMET/HIJetValidation/" + mInputJet1Collection.label() + "_DeltaRMatched_" +
                           mInputJet2Collection.label());

  mpT_ratio_Jet1Jet2 = ibooker.book1D(
      "Ratio_Jet1pT_vs_Jet2pT", Form(";Matched %s Jet pT / %s Jet pT;Counts", "JetType1", "JetType2"), 100, 0, 10);
  mpT_Jet1_matched =
      ibooker.book1D("Jet1_matched_jet_Spectra", Form(";Matched %s Spectra; counts", "JetType1"), 100, 0, 1000);
  mpT_Jet2_matched =
      ibooker.book1D("Jet2_matched_jet_Spectra", Form(";Matched %s Spectra; counts", "JetType2"), 100, 0, 1000);
  mpT_Jet1_unmatched =
      ibooker.book1D("Jet1_unmatched_jet_Spectra", Form(";Unmatched %s spectra;counts", "JetType1"), 100, 0, 1000);
  mpT_Jet2_unmatched =
      ibooker.book1D("Jet2_unmatched_jet_Spectra", Form(";Unmatched %s spectra;counts", "JetType2"), 100, 0, 1000);

  if (std::string("VsCalo") == JetType1 || std::string("PuCalo") == JetType1) {
    mHadEnergy_Jet1_unmatched =
        ibooker.book1D("HadEnergy_Jet1_unmatched",
                       Form("HadEnergy_Jet1_unmatched;HadEnergy unmatched %s;counts", "JetType1"),
                       50,
                       0,
                       200);
    mEmEnergy_Jet1_unmatched = ibooker.book1D(
        "EmEnergy_Jet1_unmatched", Form("EmEnergy_Jet1_unmatched;EMEnergy unmatched %s;counts", "JetType1"), 50, 0, 200);
  }

  if (std::string("VsCalo") == JetType2 || std::string("PuCalo") == JetType2) {
    mHadEnergy_Jet2_unmatched =
        ibooker.book1D("HadEnergy_Jet2_unmatched",
                       Form("HadEnergy_Jet2_unmatched;HadEnergy unmatched %s;counts", "JetType2"),
                       50,
                       0,
                       200);
    mEmEnergy_Jet2_unmatched = ibooker.book1D(
        "EmEnergy_Jet2_unmatched", Form("EmEnergy_Jet2_unmatched;EMEnergy unmatched %s;counts", "JetType2"), 50, 0, 200);
  }

  if (std::string("VsPF") == JetType1) {
    mChargedHadronEnergy_Jet1_unmatched = ibooker.book1D(
        "ChargedHadronEnergy_Jet1_unmatched", Form(";charged HAD energy unmatched %s;counts", "JetType1"), 100, 0, 300);
    mNeutralHadronEnergy_Jet1_unmatched = ibooker.book1D(
        "neutralHadronEnergy_Jet1_unmatched", Form(";neutral HAD energy unmatched %s;counts", "JetType1"), 100, 0, 300);
    mChargedEmEnergy_Jet1_unmatched = ibooker.book1D(
        "ChargedEmEnergy_Jet1_unmatched", Form(";charged EM energy unmatched %s;counts", "JetType1"), 100, 0, 300);
    mNeutralEmEnergy_Jet1_unmatched = ibooker.book1D(
        "neutralEmEnergy_Jet1_unmatched", Form(";neutral EM energy unmatched %s;counts", "JetType1"), 100, 0, 300);
    mChargedMuEnergy_Jet1_unmatched = ibooker.book1D(
        "ChargedMuEnergy_Jet1_unmatched", Form(";charged Mu energy unmatched %s;counts", "JetType1"), 100, 0, 300);

    mChargedHadEnergyFraction_Jet1_unmatched = ibooker.book1D(
        "ChargedHadEnergyFraction_Jet1_unmatched", Form(";h^{+/-} Energy Fraction %s;counts", "JetType1"), 50, 0, 1);
    mNeutralHadEnergyFraction_Jet1_unmatched = ibooker.book1D(
        "NeutralHadEnergyFraction_Jet1_unmatched", Form(";h^{0} Energy Fraction %s;counts", "JetType1"), 50, 0, 1);
    mPhotonEnergyFraction_Jet1_unmatched = ibooker.book1D(
        "PhotonEnergyFraction_Jet1_unmatched", Form(";#gamma Energy Fraction %s;counts", "JetType1"), 50, 0, 1);
    mElectronEnergyFraction_Jet1_unmatched = ibooker.book1D(
        "ElectronEnergyFraction_Jet1_unmatched", Form(";e Energy Fraction %s;counts", "JetType1"), 50, 0, 1);
    mMuonEnergyFraction_Jet1_unmatched = ibooker.book1D(
        "MuonoEnergyFraction_Jet1_unmatched", Form(";#mu Energy Fraction %s;counts", "JetType1"), 50, 0, 1);
  }

  if (std::string("VsPF") == JetType2) {
    mChargedHadronEnergy_Jet2_unmatched = ibooker.book1D(
        "ChargedHadronEnergy_Jet2_unmatched", Form(";charged HAD energy unmatched %s;counts", "JetType2"), 100, 0, 300);
    mNeutralHadronEnergy_Jet2_unmatched = ibooker.book1D(
        "neutralHadronEnergy_Jet2_unmatched", Form(";neutral HAD energy unmatched %s;counts", "JetType2"), 100, 0, 300);
    mChargedEmEnergy_Jet2_unmatched = ibooker.book1D(
        "ChargedEmEnergy_Jet2_unmatched", Form(";charged EM energy unmatched %s;counts", "JetType2"), 100, 0, 300);
    mNeutralEmEnergy_Jet2_unmatched = ibooker.book1D(
        "neutralEmEnergy_Jet2_unmatched", Form(";neutral EM energy unmatched %s;counts", "JetType2"), 100, 0, 300);
    mChargedMuEnergy_Jet2_unmatched = ibooker.book1D(
        "ChargedMuEnergy_Jet2_unmatched", Form(";charged Mu energy unmatched %s;counts", "JetType2"), 100, 0, 300);

    mChargedHadEnergyFraction_Jet2_unmatched = ibooker.book1D(
        "ChargedHadEnergyFraction_Jet2_unmatched", Form(";h^{+/-} Energy Fraction %s;counts", "JetType2"), 50, 0, 1);
    mNeutralHadEnergyFraction_Jet2_unmatched = ibooker.book1D(
        "NeutralHadEnergyFraction_Jet2_unmatched", Form(";h^{0} Energy Fraction %s;counts", "JetType2"), 50, 0, 1);
    mPhotonEnergyFraction_Jet2_unmatched = ibooker.book1D(
        "PhotonEnergyFraction_Jet2_unmatched", Form(";#gamma Energy Fraction %s;counts", "JetType2"), 50, 0, 1);
    mElectronEnergyFraction_Jet2_unmatched = ibooker.book1D(
        "ElectronEnergyFraction_Jet2_unmatched", Form(";e Energy Fraction %s;counts", "JetType2"), 50, 0, 1);
    mMuonEnergyFraction_Jet2_unmatched = ibooker.book1D(
        "MuonoEnergyFraction_Jet2_unmatched", Form(";#mu Energy Fraction %s;counts", "JetType2"), 50, 0, 1);
  }

  if (mOutputFile.empty())
    LogInfo("OutputInfo") << " Histograms will NOT be saved";
  else
    LogInfo("OutputInfo") << " Histograms will be saved to file:" << mOutputFile;
}

//------------------------------------------------------------------------------
// ~JetAnalyzer_HeavyIons
//------------------------------------------------------------------------------
JetAnalyzer_HeavyIons_matching::~JetAnalyzer_HeavyIons_matching() {}

//------------------------------------------------------------------------------
// beginJob
//------------------------------------------------------------------------------
//void JetAnalyzer_HeavyIons_matching::beginJob() {
//  std::cout<<"inside the begin job function"<<endl;
//}

//------------------------------------------------------------------------------
// endJob
//------------------------------------------------------------------------------
//void JetAnalyzer_HeavyIons_matching::endJob()
//{
//  if (!mOutputFile.empty() && &*edm::Service<DQMStore>())
//    {
//      edm::Service<DQMStore>()->save(mOutputFile);
//    }
//}

//------------------------------------------------------------------------------
// analyze
//------------------------------------------------------------------------------
void JetAnalyzer_HeavyIons_matching::analyze(const edm::Event& mEvent, const edm::EventSetup& mSetup) {
  // Get the Jet collection
  //----------------------------------------------------------------------------

  std::vector<const Jet*> recoJet1;
  recoJet1.clear();
  std::vector<const Jet*> recoJet2;
  recoJet2.clear();

  edm::Handle<CaloJetCollection> caloJet1;
  edm::Handle<CaloJetCollection> caloJet2;
  edm::Handle<JPTJetCollection> jptJets;
  edm::Handle<PFJetCollection> pfJets;
  edm::Handle<BasicJetCollection> basicJets;

  if (std::string("VsCalo") == JetType1) {
    mEvent.getByToken(caloJet1Token_, caloJet1);
    for (unsigned ijet = 0; ijet < caloJet1->size(); ++ijet)
      recoJet1.push_back(&(*caloJet1)[ijet]);
  }
  if (std::string("PuCalo") == JetType1) {
    mEvent.getByToken(caloJet2Token_, caloJet1);
    for (unsigned ijet = 0; ijet < caloJet1->size(); ++ijet)
      recoJet1.push_back(&(*caloJet1)[ijet]);
  }
  if (std::string("VsPF") == JetType1) {
    mEvent.getByToken(pfJetsToken_, pfJets);
    for (unsigned ijet = 0; ijet < pfJets->size(); ++ijet)
      recoJet1.push_back(&(*pfJets)[ijet]);
  }
  if (std::string("PuPF") == JetType1) {
    mEvent.getByToken(basicJetsToken_, basicJets);
    for (unsigned ijet = 0; ijet < basicJets->size(); ++ijet)
      recoJet1.push_back(&(*basicJets)[ijet]);
  }

  if (std::string("VsCalo") == JetType2) {
    mEvent.getByToken(caloJet1Token_, caloJet2);
    for (unsigned ijet = 0; ijet < caloJet2->size(); ++ijet)
      recoJet2.push_back(&(*caloJet2)[ijet]);
  }
  if (std::string("PuCalo") == JetType2) {
    mEvent.getByToken(caloJet2Token_, caloJet2);
    for (unsigned ijet = 0; ijet < caloJet2->size(); ++ijet)
      recoJet2.push_back(&(*caloJet2)[ijet]);
  }
  if (std::string("VsPF") == JetType2) {
    mEvent.getByToken(pfJetsToken_, pfJets);
    for (unsigned ijet = 0; ijet < pfJets->size(); ++ijet)
      recoJet2.push_back(&(*pfJets)[ijet]);
  }
  if (std::string("PuPF") == JetType2) {
    mEvent.getByToken(basicJetsToken_, basicJets);
    for (unsigned ijet = 0; ijet < basicJets->size(); ++ijet)
      recoJet2.push_back(&(*basicJets)[ijet]);
  }

  // start to perform the matching - between recoJet1 and recoJet2.

  Int_t Jet1_nref = recoJet1.size();
  Int_t Jet2_nref = recoJet2.size();

  int jet1 = 0;
  int jet2 = 0;

  std::vector<MyJet> vJet1, vJet2;
  std::vector<int> Jet1_ID(Jet1_nref), Jet2_ID(Jet2_nref);

  if (Jet1_nref == 0 || Jet2_nref == 0)
    return;

  for (unsigned ijet1 = 0; ijet1 < recoJet1.size(); ++ijet1) {
    if (recoJet1[ijet1]->pt() < mRecoJetPtThreshold)
      continue;
    if (fabs(recoJet1[ijet1]->eta()) < mRecoJetEtaCut)
      continue;

    MyJet JET1;
    JET1.eta = recoJet1[ijet1]->eta();
    JET1.phi = recoJet1[ijet1]->phi();
    JET1.pt = recoJet1[ijet1]->pt();
    JET1.id = ijet1;

    vJet1.push_back(JET1);
    jet1++;

  }  // first jet loop

  for (unsigned ijet2 = 0; ijet2 < recoJet2.size(); ++ijet2) {
    if (recoJet2[ijet2]->pt() < mRecoJetPtThreshold)
      continue;
    if (fabs(recoJet2[ijet2]->eta()) < mRecoJetEtaCut)
      continue;

    MyJet JET2;
    JET2.eta = recoJet2[ijet2]->eta();
    JET2.phi = recoJet2[ijet2]->phi();
    JET2.pt = recoJet2[ijet2]->pt();
    JET2.id = ijet2;

    vJet2.push_back(JET2);
    jet2++;

  }  // second jet loop

  bool onlyJet1 = (jet1 > 0 && jet2 == 0) ? true : false;
  bool onlyJet2 = (jet1 == 0 && jet2 > 0) ? true : false;
  bool bothJet1Jet2 = (jet1 > 0 && jet2 > 0) ? true : false;

  int matchedJets = 0;
  int unmatchedJet1 = 0;
  int unmatchedJet2 = 0;

  std::vector<MyJet>::const_iterator iJet, jJet;

  if (onlyJet1) {
    for (iJet = vJet1.begin(); iJet != vJet1.end(); ++iJet) {
      int pj = (*iJet).id;

      mpT_Jet1_unmatched->Fill(recoJet1[pj]->pt());
    }

  } else if (onlyJet2) {
    for (iJet = vJet2.begin(); iJet != vJet2.end(); ++iJet) {
      int cj = (*iJet).id;

      mpT_Jet2_unmatched->Fill(recoJet2[cj]->pt());
    }

  } else if (bothJet1Jet2) {
    ABMatchedJets mABMatchedJets;

    for (iJet = vJet1.begin(); iJet != vJet1.end(); ++iJet) {
      for (jJet = vJet2.begin(); jJet != vJet2.end(); ++jJet) {
        mABMatchedJets.insert(std::make_pair(*iJet, *jJet));
      }
    }

    ABItr itr;
    // matched Jets matching Jet 1 to Jet 2
    for (itr = mABMatchedJets.begin(); itr != mABMatchedJets.end(); ++itr) {
      ABJetPair jetpair = (*itr);
      MyJet Aj = jetpair.first;
      MyJet Bj = jetpair.second;

      float delr = JetAnalyzer_HeavyIons_matching::deltaRR(Bj.eta, Bj.phi, Aj.eta, Aj.phi);

      if (delr < mRecoDelRMatch && Jet1_ID[Aj.id] == 0) {
        mpT_ratio_Jet1Jet2->Fill((Float_t)recoJet2[Bj.id]->pt() / recoJet1[Aj.id]->pt());

        mpT_Jet1_matched->Fill(recoJet1[Aj.id]->pt());
        mpT_Jet2_matched->Fill(recoJet2[Bj.id]->pt());

        Jet1_ID[Aj.id] = 1;
        Jet2_ID[Bj.id] = 1;

        matchedJets++;
      }
    }

    // for unmatched Jets
    for (itr = mABMatchedJets.begin(); itr != mABMatchedJets.end(); ++itr) {
      ABJetPair jetpair = (*itr);

      MyJet Aj = jetpair.first;
      MyJet Bj = jetpair.second;

      if (Jet1_ID[Aj.id] == 0) {
        mpT_Jet1_unmatched->Fill(recoJet1[Aj.id]->pt());
        unmatchedJet1++;
        Jet1_ID[Aj.id] = 1;

        if (std::string("VsCalo") == JetType1 || std::string("PuCalo") == JetType1) {
          mHadEnergy_Jet1_unmatched->Fill((*caloJet1)[Aj.id].hadEnergyInHO() + (*caloJet1)[Aj.id].hadEnergyInHB() +
                                          (*caloJet1)[Aj.id].hadEnergyInHF() + (*caloJet1)[Aj.id].hadEnergyInHE());
          mEmEnergy_Jet1_unmatched->Fill((*caloJet1)[Aj.id].emEnergyInEB() + (*caloJet1)[Aj.id].emEnergyInEE() +
                                         (*caloJet1)[Aj.id].emEnergyInHF());
        }

        if (std::string("VsPF") == JetType1) {
          mChargedHadronEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].chargedHadronEnergy());
          mNeutralHadronEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].neutralHadronEnergy());
          mChargedEmEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].chargedEmEnergy());
          mNeutralEmEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].neutralEmEnergy());
          mChargedMuEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].chargedMuEnergy());

          mChargedHadEnergyFraction_Jet1_unmatched->Fill((*pfJets)[Aj.id].chargedHadronEnergyFraction());
          mNeutralHadEnergyFraction_Jet1_unmatched->Fill((*pfJets)[Aj.id].neutralHadronEnergyFraction());
          mPhotonEnergyFraction_Jet1_unmatched->Fill((*pfJets)[Aj.id].photonEnergyFraction());
          mElectronEnergyFraction_Jet1_unmatched->Fill((*pfJets)[Aj.id].electronEnergyFraction());
          mMuonEnergyFraction_Jet1_unmatched->Fill((*pfJets)[Aj.id].muonEnergyFraction());
        }
      }

      if (Jet2_ID[Bj.id] == 0) {
        mpT_Jet2_unmatched->Fill(recoJet2[Bj.id]->pt());
        unmatchedJet2++;
        Jet2_ID[Bj.id] = 2;
        if (std::string("VsCalo") == JetType2 || std::string("PuCalo") == JetType2) {
          mHadEnergy_Jet2_unmatched->Fill((*caloJet2)[Bj.id].hadEnergyInHO() + (*caloJet2)[Bj.id].hadEnergyInHB() +
                                          (*caloJet2)[Bj.id].hadEnergyInHF() + (*caloJet2)[Bj.id].hadEnergyInHE());
          mEmEnergy_Jet2_unmatched->Fill((*caloJet2)[Bj.id].emEnergyInEB() + (*caloJet2)[Bj.id].emEnergyInEE() +
                                         (*caloJet2)[Bj.id].emEnergyInHF());
        }

        if (std::string("VsPF") == JetType2) {
          mChargedHadronEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].chargedHadronEnergy());
          mNeutralHadronEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].neutralHadronEnergy());
          mChargedEmEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].chargedEmEnergy());
          mNeutralEmEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].neutralEmEnergy());
          mChargedMuEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].chargedMuEnergy());

          mChargedHadEnergyFraction_Jet2_unmatched->Fill((*pfJets)[Bj.id].chargedHadronEnergyFraction());
          mNeutralHadEnergyFraction_Jet2_unmatched->Fill((*pfJets)[Bj.id].neutralHadronEnergyFraction());
          mPhotonEnergyFraction_Jet2_unmatched->Fill((*pfJets)[Bj.id].photonEnergyFraction());
          mElectronEnergyFraction_Jet2_unmatched->Fill((*pfJets)[Bj.id].electronEnergyFraction());
          mMuonEnergyFraction_Jet2_unmatched->Fill((*pfJets)[Bj.id].muonEnergyFraction());
        }
      }
    }

  }  // both Jet1 and Jet2 in the event.
}