HTMonitor.cc

Go to the documentation of this file.

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DQMOffline/Trigger/plugins/HTMonitor.h"

#include "CommonTools/TriggerUtils/interface/GenericTriggerEventFlag.h"

#include "DataFormats/Math/interface/deltaPhi.h"

// -----------------------------
//  constructors and destructor
// -----------------------------

HTMonitor::HTMonitor(const edm::ParameterSet& iConfig)
    : folderName_(iConfig.getParameter<std::string>("FolderName")),
      metInputTag_(iConfig.getParameter<edm::InputTag>("met")),
      jetInputTag_(iConfig.getParameter<edm::InputTag>("jets")),
      eleInputTag_(iConfig.getParameter<edm::InputTag>("electrons")),
      muoInputTag_(iConfig.getParameter<edm::InputTag>("muons")),
      vtxInputTag_(iConfig.getParameter<edm::InputTag>("vertices")),
      metToken_(consumes<reco::PFMETCollection>(metInputTag_)),
      jetToken_(mayConsume<reco::JetView>(jetInputTag_)),
      eleToken_(mayConsume<reco::GsfElectronCollection>(eleInputTag_)),
      muoToken_(mayConsume<reco::MuonCollection>(muoInputTag_)),
      vtxToken_(mayConsume<reco::VertexCollection>(vtxInputTag_)),
      ht_variable_binning_(
          iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<std::vector<double> >("htBinning")),
      ht_binning_(
          getHistoPSet(iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet>("htPSet"))),
      ls_binning_(
          getHistoPSet(iConfig.getParameter<edm::ParameterSet>("histoPSet").getParameter<edm::ParameterSet>("lsPSet"))),
      num_genTriggerEventFlag_(new GenericTriggerEventFlag(
          iConfig.getParameter<edm::ParameterSet>("numGenericTriggerEventPSet"), consumesCollector(), *this)),
      den_genTriggerEventFlag_(new GenericTriggerEventFlag(
          iConfig.getParameter<edm::ParameterSet>("denGenericTriggerEventPSet"), consumesCollector(), *this)),
      metSelection_(iConfig.getParameter<std::string>("metSelection")),
      jetSelection_(iConfig.getParameter<std::string>("jetSelection")),
      eleSelection_(iConfig.getParameter<std::string>("eleSelection")),
      muoSelection_(iConfig.getParameter<std::string>("muoSelection")),
      jetSelection_HT_(iConfig.getParameter<std::string>("jetSelection_HT")),
      njets_(iConfig.getParameter<unsigned>("njets")),
      nelectrons_(iConfig.getParameter<unsigned>("nelectrons")),
      nmuons_(iConfig.getParameter<unsigned>("nmuons")),
      dEtaCut_(iConfig.getParameter<double>("dEtaCut")) {
  /* mia: THIS CODE SHOULD BE DELETED !!!! */
  string quantity = iConfig.getParameter<std::string>("quantity");
  if (quantity == "HT") {
    quantity_ = HT;
  } else if (quantity == "Mjj") {
    quantity_ = MJJ;
  } else if (quantity == "softdrop") {
    quantity_ = SOFTDROP;
  } else {
    throw cms::Exception("quantity not defined")
        << "the quantity '" << quantity << "' is undefined. Please check your config!" << std::endl;
  }

  // this vector has to be alligned to the the number of Tokens accessed by this module
  warningPrinted4token_.push_back(false);  // PFMETCollection
  warningPrinted4token_.push_back(false);  // JetCollection
  warningPrinted4token_.push_back(false);  // GsfElectronCollection
  warningPrinted4token_.push_back(false);  // MuonCollection
  warningPrinted4token_.push_back(false);  // VertexCollection
}

HTMonitor::~HTMonitor() = default;

HTMonitor::MEHTbinning HTMonitor::getHistoPSet(const edm::ParameterSet& pset) {
  return HTMonitor::MEHTbinning{
      pset.getParameter<unsigned>("nbins"),
      pset.getParameter<double>("xmin"),
      pset.getParameter<double>("xmax"),
  };
}

HTMonitor::MEHTbinning HTMonitor::getHistoLSPSet(const edm::ParameterSet& pset) {
  return HTMonitor::MEHTbinning{pset.getParameter<unsigned>("nbins"), 0., double(pset.getParameter<unsigned>("nbins"))};
}

void HTMonitor::setHTitle(HTME& me, const std::string& titleX, const std::string& titleY) {
  me.numerator->setAxisTitle(titleX, 1);
  me.numerator->setAxisTitle(titleY, 2);
  me.denominator->setAxisTitle(titleX, 1);
  me.denominator->setAxisTitle(titleY, 2);
}

void HTMonitor::bookME(DQMStore::IBooker& ibooker,
                       HTME& me,
                       const std::string& histname,
                       const std::string& histtitle,
                       int nbins,
                       double min,
                       double max) {
  me.numerator = ibooker.book1D(histname + "_numerator", histtitle + " (numerator)", nbins, min, max);
  me.denominator = ibooker.book1D(histname + "_denominator", histtitle + " (denominator)", nbins, min, max);
}
void HTMonitor::bookME(DQMStore::IBooker& ibooker,
                       HTME& me,
                       const std::string& histname,
                       const std::string& histtitle,
                       const std::vector<double>& binning) {
  int nbins = binning.size() - 1;
  std::vector<float> fbinning(binning.begin(), binning.end());
  float* arr = &fbinning[0];
  me.numerator = ibooker.book1D(histname + "_numerator", histtitle + " (numerator)", nbins, arr);
  me.denominator = ibooker.book1D(histname + "_denominator", histtitle + " (denominator)", nbins, arr);
}
void HTMonitor::bookME(DQMStore::IBooker& ibooker,
                       HTME& me,
                       const std::string& histname,
                       const std::string& histtitle,
                       int nbinsX,
                       double xmin,
                       double xmax,
                       double ymin,
                       double ymax) {
  me.numerator =
      ibooker.bookProfile(histname + "_numerator", histtitle + " (numerator)", nbinsX, xmin, xmax, ymin, ymax);
  me.denominator =
      ibooker.bookProfile(histname + "_denominator", histtitle + " (denominator)", nbinsX, xmin, xmax, ymin, ymax);
}
void HTMonitor::bookME(DQMStore::IBooker& ibooker,
                       HTME& me,
                       const std::string& histname,
                       const std::string& histtitle,
                       int nbinsX,
                       double xmin,
                       double xmax,
                       int nbinsY,
                       double ymin,
                       double ymax) {
  me.numerator =
      ibooker.book2D(histname + "_numerator", histtitle + " (numerator)", nbinsX, xmin, xmax, nbinsY, ymin, ymax);
  me.denominator =
      ibooker.book2D(histname + "_denominator", histtitle + " (denominator)", nbinsX, xmin, xmax, nbinsY, ymin, ymax);
}
void HTMonitor::bookME(DQMStore::IBooker& ibooker,
                       HTME& me,
                       const std::string& histname,
                       const std::string& histtitle,
                       const std::vector<double>& binningX,
                       const std::vector<double>& binningY) {
  int nbinsX = binningX.size() - 1;
  std::vector<float> fbinningX(binningX.begin(), binningX.end());
  float* arrX = &fbinningX[0];
  int nbinsY = binningY.size() - 1;
  std::vector<float> fbinningY(binningY.begin(), binningY.end());
  float* arrY = &fbinningY[0];

  me.numerator = ibooker.book2D(histname + "_numerator", histtitle + " (numerator)", nbinsX, arrX, nbinsY, arrY);
  me.denominator = ibooker.book2D(histname + "_denominator", histtitle + " (denominator)", nbinsX, arrX, nbinsY, arrY);
}

void HTMonitor::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) {
  std::string histname, histtitle;

  std::string currentFolder = folderName_;
  ibooker.setCurrentFolder(currentFolder);

  switch (quantity_) {
    case HT: {
      histname = "ht_variable";
      histtitle = "HT";
      bookME(ibooker, qME_variableBinning_, histname, histtitle, ht_variable_binning_);
      setHTitle(qME_variableBinning_, "HT [GeV]", "events / [GeV]");

      histname = "htVsLS";
      histtitle = "HT vs LS";
      bookME(ibooker,
             htVsLS_,
             histname,
             histtitle,
             ls_binning_.nbins,
             ls_binning_.xmin,
             ls_binning_.xmax,
             ht_binning_.xmin,
             ht_binning_.xmax);
      setHTitle(htVsLS_, "LS", "HT [GeV]");

      histname = "deltaphi_metjet1";
      histtitle = "DPHI_METJ1";
      bookME(ibooker, deltaphimetj1ME_, histname, histtitle, phi_binning_.nbins, phi_binning_.xmin, phi_binning_.xmax);
      setHTitle(deltaphimetj1ME_, "delta phi (met, j1)", "events / 0.1 rad");

      histname = "deltaphi_jet1jet2";
      histtitle = "DPHI_J1J2";
      bookME(ibooker, deltaphij1j2ME_, histname, histtitle, phi_binning_.nbins, phi_binning_.xmin, phi_binning_.xmax);
      setHTitle(deltaphij1j2ME_, "delta phi (j1, j2)", "events / 0.1 rad");
      break;
    }

    case MJJ: {
      histname = "mjj_variable";
      histtitle = "Mjj";
      bookME(ibooker, qME_variableBinning_, histname, histtitle, ht_variable_binning_);
      setHTitle(qME_variableBinning_, "Mjj [GeV]", "events / [GeV]");
      break;
    }

    case SOFTDROP: {
      histname = "softdrop_variable";
      histtitle = "softdropmass";
      bookME(ibooker, qME_variableBinning_, histname, histtitle, ht_variable_binning_);
      setHTitle(qME_variableBinning_, "leading jet softdropmass [GeV]", "events / [GeV]");
      break;
    }
  }

  // Initialize the GenericTriggerEventFlag
  if (num_genTriggerEventFlag_ && num_genTriggerEventFlag_->on())
    num_genTriggerEventFlag_->initRun(iRun, iSetup);
  if (den_genTriggerEventFlag_ && den_genTriggerEventFlag_->on())
    den_genTriggerEventFlag_->initRun(iRun, iSetup);
}

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
void HTMonitor::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
  // Filter out events if Trigger Filtering is requested
  if (den_genTriggerEventFlag_->on() && !den_genTriggerEventFlag_->accept(iEvent, iSetup))
    return;

  edm::Handle<reco::PFMETCollection> metHandle;
  iEvent.getByToken(metToken_, metHandle);
  if (!metHandle.isValid()) {
    if (!warningPrinted4token_[0]) {
      edm::LogWarning("HTMonitor") << "skipping events because the collection " << metInputTag_.label().c_str()
                                   << " is not available";
      warningPrinted4token_[0] = true;
    }
    return;
  }
  reco::PFMET pfmet = metHandle->front();
  if (!metSelection_(pfmet))
    return;

  edm::Handle<reco::JetView> jetHandle;  //add a configurable jet collection & jet pt selection
  iEvent.getByToken(jetToken_, jetHandle);
  if (!jetHandle.isValid()) {
    if (!warningPrinted4token_[1]) {
      edm::LogWarning("HTMonitor") << "skipping events because the collection " << jetInputTag_.label().c_str()
                                   << " is not available";
      warningPrinted4token_[1] = true;
    }
    return;
  }
  std::vector<reco::Jet> jets;
  if (jetHandle->size() < njets_)
    return;
  for (auto const& j : *jetHandle) {
    if (jetSelection_(j)) {
      jets.push_back(j);
    }
  }

  if (jets.size() < njets_)
    return;

  float deltaPhi_met_j1 = 10.0;
  float deltaPhi_j1_j2 = 10.0;

  if (!jets.empty())
    deltaPhi_met_j1 = fabs(deltaPhi(pfmet.phi(), jets[0].phi()));
  if (jets.size() >= 2)
    deltaPhi_j1_j2 = fabs(deltaPhi(jets[0].phi(), jets[1].phi()));

  std::vector<reco::GsfElectron> electrons;
  edm::Handle<reco::GsfElectronCollection> eleHandle;
  iEvent.getByToken(eleToken_, eleHandle);
  if (eleHandle.isValid()) {
    if (eleHandle->size() < nelectrons_)
      return;
    for (auto const& e : *eleHandle) {
      if (eleSelection_(e))
        electrons.push_back(e);
    }
    if (electrons.size() < nelectrons_)
      return;
  } else {
    if (!warningPrinted4token_[2]) {
      warningPrinted4token_[2] = true;
      if (eleInputTag_.label().empty())
        edm::LogWarning("HTMonitor") << "GsfElectronCollection not set";
      else
        edm::LogWarning("HTMonitor") << "skipping events because the collection " << eleInputTag_.label().c_str()
                                     << " is not available";
    }
    if (!eleInputTag_.label().empty())
      return;
  }

  reco::Vertex vtx;
  edm::Handle<reco::VertexCollection> vtxHandle;
  iEvent.getByToken(vtxToken_, vtxHandle);
  if (vtxHandle.isValid()) {
    for (auto const& v : *vtxHandle) {
      bool isFake = v.isFake();

      if (!isFake) {
        vtx = v;
        break;
      }
    }
  } else {
    if (!warningPrinted4token_[3]) {
      warningPrinted4token_[3] = true;
      if (vtxInputTag_.label().empty())
        edm::LogWarning("HTMonitor") << "VertexCollection not set";
      else
        edm::LogWarning("HTMonitor") << "skipping events because the collection " << vtxInputTag_.label().c_str()
                                     << " is not available";
    }
    if (!vtxInputTag_.label().empty())
      return;
  }

  std::vector<reco::Muon> muons;
  edm::Handle<reco::MuonCollection> muoHandle;
  iEvent.getByToken(muoToken_, muoHandle);
  if (muoHandle.isValid()) {
    if (muoHandle->size() < nmuons_)
      return;
    for (auto const& m : *muoHandle) {
      if (muoSelection_(m) && m.isGlobalMuon() && m.isPFMuon() && m.globalTrack()->normalizedChi2() < 10. &&
          m.globalTrack()->hitPattern().numberOfValidMuonHits() > 0 && m.numberOfMatchedStations() > 1 &&
          fabs(m.muonBestTrack()->dxy(vtx.position())) < 0.2 && fabs(m.muonBestTrack()->dz(vtx.position())) < 0.5 &&
          m.innerTrack()->hitPattern().numberOfValidPixelHits() > 0 &&
          m.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5)
        muons.push_back(m);
    }
    if (muons.size() < nmuons_)
      return;
  } else {
    if (!warningPrinted4token_[4]) {
      warningPrinted4token_[4] = true;
      if (muoInputTag_.label().empty())
        edm::LogWarning("HTMonitor") << "MuonCollection not set";
      else
        edm::LogWarning("HTMonitor") << "skipping events because the collection " << muoInputTag_.label().c_str()
                                     << " is not available";
    }
    if (!muoInputTag_.label().empty())
      return;
  }

  // fill histograms
  switch (quantity_) {
    case HT: {
      float ht = 0.0;
      for (auto const& j : *jetHandle) {
        if (jetSelection_HT_(j))
          ht += j.pt();
      }

      // filling histograms (denominator)
      qME_variableBinning_.denominator->Fill(ht);

      deltaphimetj1ME_.denominator->Fill(deltaPhi_met_j1);
      deltaphij1j2ME_.denominator->Fill(deltaPhi_j1_j2);

      int ls = iEvent.id().luminosityBlock();
      htVsLS_.denominator->Fill(ls, ht);

      // applying selection for numerator
      if (num_genTriggerEventFlag_->on() && !num_genTriggerEventFlag_->accept(iEvent, iSetup))
        return;

      // filling histograms (num_genTriggerEventFlag_)
      qME_variableBinning_.numerator->Fill(ht);

      htVsLS_.numerator->Fill(ls, ht);
      deltaphimetj1ME_.numerator->Fill(deltaPhi_met_j1);
      deltaphij1j2ME_.numerator->Fill(deltaPhi_j1_j2);
      break;
    }

    case MJJ: {
      if (jets.size() < 2)
        return;

      // deltaEta cut
      if (fabs(jets[0].p4().Eta() - jets[1].p4().Eta()) >= dEtaCut_)
        return;
      float mjj = (jets[0].p4() + jets[1].p4()).M();

      qME_variableBinning_.denominator->Fill(mjj);

      // applying selection for numerator
      if (num_genTriggerEventFlag_->on() && !num_genTriggerEventFlag_->accept(iEvent, iSetup))
        return;

      qME_variableBinning_.numerator->Fill(mjj);
      break;
    }

    case SOFTDROP: {
      if (jets.size() < 2)
        return;

      // deltaEta cut
      if (fabs(jets[0].p4().Eta() - jets[1].p4().Eta()) >= dEtaCut_)
        return;

      float softdrop = jets[0].p4().M();

      qME_variableBinning_.denominator->Fill(softdrop);

      // applying selection for numerator
      if (num_genTriggerEventFlag_->on() && !num_genTriggerEventFlag_->accept(iEvent, iSetup))
        return;

      qME_variableBinning_.numerator->Fill(softdrop);
      break;
    }
  }
}

void HTMonitor::fillHistoPSetDescription(edm::ParameterSetDescription& pset) {
  pset.add<unsigned int>("nbins");
  pset.add<double>("xmin");
  pset.add<double>("xmax");
}

void HTMonitor::fillHistoLSPSetDescription(edm::ParameterSetDescription& pset) {
  pset.add<unsigned int>("nbins", 2500);
  pset.add<double>("xmin", 0.);
  pset.add<double>("xmax", 2500.);
}

void HTMonitor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<std::string>("FolderName", "HLT/HT");
  desc.add<std::string>("quantity", "HT");

  desc.add<edm::InputTag>("met", edm::InputTag("pfMet"));
  desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
  desc.add<edm::InputTag>("electrons", edm::InputTag("gedGsfElectrons"));
  desc.add<edm::InputTag>("muons", edm::InputTag("muons"));
  desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
  desc.add<std::string>("metSelection", "pt > 0");
  desc.add<std::string>("jetSelection", "pt > 0");
  desc.add<std::string>("eleSelection", "pt > 0");
  desc.add<std::string>("muoSelection", "pt > 0");
  desc.add<std::string>("jetSelection_HT", "pt > 30 && eta < 2.5");
  desc.add<unsigned>("njets", 0);
  desc.add<unsigned>("nelectrons", 0);
  desc.add<unsigned>("nmuons", 0);
  desc.add<double>("dEtaCut", 1.3);

  edm::ParameterSetDescription genericTriggerEventPSet;
  genericTriggerEventPSet.add<bool>("andOr");
  genericTriggerEventPSet.add<edm::InputTag>("dcsInputTag", edm::InputTag("scalersRawToDigi"));
  genericTriggerEventPSet.add<std::vector<int> >("dcsPartitions", {});
  genericTriggerEventPSet.add<bool>("andOrDcs", false);
  genericTriggerEventPSet.add<bool>("errorReplyDcs", true);
  genericTriggerEventPSet.add<std::string>("dbLabel", "");
  genericTriggerEventPSet.add<bool>("andOrHlt", true);
  genericTriggerEventPSet.add<edm::InputTag>("hltInputTag", edm::InputTag("TriggerResults::HLT"));
  genericTriggerEventPSet.add<std::vector<std::string> >("hltPaths", {});
  genericTriggerEventPSet.add<std::string>("hltDBKey", "");
  genericTriggerEventPSet.add<bool>("errorReplyHlt", false);
  genericTriggerEventPSet.add<unsigned int>("verbosityLevel", 1);

  desc.add<edm::ParameterSetDescription>("numGenericTriggerEventPSet", genericTriggerEventPSet);
  desc.add<edm::ParameterSetDescription>("denGenericTriggerEventPSet", genericTriggerEventPSet);

  edm::ParameterSetDescription histoPSet;
  edm::ParameterSetDescription htPSet;
  fillHistoPSetDescription(htPSet);
  histoPSet.add<edm::ParameterSetDescription>("htPSet", htPSet);
  std::vector<double> bins = {0.,   20.,  40.,   60.,   80.,   90.,   100.,  110.,  120.,  130.,  140.,
                              150., 160., 170.,  180.,  190.,  200.,  220.,  240.,  260.,  280.,  300.,
                              350., 400., 450.,  500.,  550.,  600.,  650.,  700.,  750.,  800.,  850.,
                              900., 950., 1000., 1050., 1100., 1200., 1300., 1400., 1500., 2000., 2500.};
  histoPSet.add<std::vector<double> >("htBinning", bins);

  edm::ParameterSetDescription lsPSet;
  fillHistoLSPSetDescription(lsPSet);
  histoPSet.add<edm::ParameterSetDescription>("lsPSet", lsPSet);

  desc.add<edm::ParameterSetDescription>("histoPSet", histoPSet);

  descriptions.add("htMonitoring", desc);
}

// Define this as a plug-in
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(HTMonitor);