TopDQMHelpers.h

Go to the documentation of this file.

#ifndef TOPDQMHELPERS
#define TOPDQMHELPERS

#include <string>
#include <vector>
#include <iostream>
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Common/interface/TriggerNames.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include <DataFormats/METReco/interface/PFMET.h>
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/EDConsumerBase.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/MET.h"

inline bool accept(const edm::Event& event, const edm::TriggerResults& triggerTable, const std::string& triggerPath) {
  bool passed = false;
  const edm::TriggerNames& triggerNames = event.triggerNames(triggerTable);
  for (unsigned int i = 0; i < triggerNames.triggerNames().size(); ++i) {
    if (triggerNames.triggerNames()[i] == triggerPath) {
      if (triggerTable.accept(i)) {
        passed = true;
        break;
      }
    }
  }
  return passed;
}

inline bool accept(const edm::Event& event,
                   const edm::TriggerResults& triggerTable,
                   const std::vector<std::string>& triggerPaths) {
  bool passed = false;
  for (unsigned int j = 0; j < triggerPaths.size(); ++j) {
    if (accept(event, triggerTable, triggerPaths[j])) {
      passed = true;
      break;
    }
  }
  return passed;
}

#include "DataFormats/JetReco/interface/Jet.h"
#include "FWCore/Framework/interface/EventSetup.h"

class Calculate_miniAOD {
public:
  Calculate_miniAOD(int maxNJets, double wMass);
  ~Calculate_miniAOD(){};

  double massWBoson(const std::vector<pat::Jet>& jets);
  double massTopQuark(const std::vector<pat::Jet>& jets);
  // double massBTopQuark(const std::vector<reco::Jet>& jets, std::vector<bool>
  // bjet);
  double massBTopQuark(const std::vector<pat::Jet>& jets, std::vector<double> VbtagWP, double btagWP_);

  double tmassWBoson(pat::Muon* lep, const pat::MET& met, const pat::Jet& b);
  double tmassTopQuark(pat::Muon* lep, const pat::MET& met, const pat::Jet& b);

  double tmassWBoson(pat::Electron* lep, const pat::MET& met, const pat::Jet& b);
  double tmassTopQuark(pat::Electron* lep, const pat::MET& met, const pat::Jet& b);

private:
  void operator()(const std::vector<pat::Jet>& jets);
  void operator2(const std::vector<pat::Jet>&, std::vector<double>, double);
  void operator()(const pat::Jet& bJet, pat::Electron* lepton, const pat::MET& met);

  void operator()(const pat::Jet& bJet, pat::Muon* lepton, const pat::MET& met);

private:
  bool failed_;
  int maxNJets_;
  double wMass_;
  double massWBoson_;
  double massTopQuark_;
  double massBTopQuark_;
  double tmassWBoson_;
  double tmassTopQuark_;
};

class Calculate {
public:
  Calculate(int maxNJets, double wMass);
  ~Calculate(){};

  double massWBoson(const std::vector<reco::Jet>& jets);
  double massTopQuark(const std::vector<reco::Jet>& jets);
  // double massBTopQuark(const std::vector<reco::Jet>& jets, std::vector<bool>
  // bjet);
  double massBTopQuark(const std::vector<reco::Jet>& jets, std::vector<double> VbtagWP, double btagWP_);

  double tmassWBoson(reco::RecoCandidate* lep, const reco::MET& met, const reco::Jet& b);
  double tmassTopQuark(reco::RecoCandidate* lep, const reco::MET& met, const reco::Jet& b);

private:
  void operator()(const std::vector<reco::Jet>& jets);
  void operator2(const std::vector<reco::Jet>&, std::vector<double>, double);
  void operator()(const reco::Jet& bJet, reco::RecoCandidate* lepton, const reco::MET& met);

private:
  bool failed_;
  int maxNJets_;
  double wMass_;
  double massWBoson_;
  double massTopQuark_;
  double massBTopQuark_;
  double tmassWBoson_;
  double tmassTopQuark_;
};

#include "DataFormats/JetReco/interface/JetID.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/BTauReco/interface/JetTag.h"
#include "DataFormats/Common/interface/ValueMap.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "JetMETCorrections/Objects/interface/JetCorrector.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"

template <typename Object>
class SelectionStep {
public:
  SelectionStep(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC);
  ~SelectionStep(){};

  bool select(const edm::Event& event);
  bool select(const edm::Event& event, const std::string& type);
  bool select(const edm::Event& event, const edm::EventSetup& setup);
  bool selectVertex(const edm::Event& event);

private:
  edm::EDGetTokenT<edm::View<Object> > src_;
  int min_, max_;
  edm::EDGetTokenT<edm::ValueMap<float> > electronId_;
  double eidCutValue_;
  edm::ESGetToken<JetCorrector, JetCorrectionsRecord> jetCorrector_;
  edm::EDGetTokenT<reco::JetTagCollection> btagLabel_;
  double btagWorkingPoint_;
  edm::EDGetTokenT<reco::JetIDValueMap> jetIDLabel_;

  edm::EDGetTokenT<edm::View<reco::Vertex> > pvs_;
  edm::EDGetTokenT<edm::View<reco::GsfElectron> > gsfEs_;
  StringCutObjectSelector<Object> select_;
  StringCutObjectSelector<reco::JetID>* jetIDSelect_;
};

template <typename Object>
SelectionStep<Object>::SelectionStep(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
    : select_(cfg.getParameter<std::string>("select")), jetIDSelect_(nullptr) {
  src_ = iC.consumes<edm::View<Object> >(cfg.getParameter<edm::InputTag>("src"));
  // exist otherwise they are initialized with -1
  cfg.exists("min") ? min_ = cfg.getParameter<int>("min") : min_ = -1;
  cfg.exists("max") ? max_ = cfg.getParameter<int>("max") : max_ = -1;
  std::string mygSF = "gedGsfElectrons";
  gsfEs_ = iC.consumes<edm::View<reco::GsfElectron> >(cfg.getUntrackedParameter<edm::InputTag>("myGSF", mygSF));
  if (cfg.existsAs<edm::ParameterSet>("electronId")) {
    edm::ParameterSet elecId = cfg.getParameter<edm::ParameterSet>("electronId");
    electronId_ = iC.consumes<edm::ValueMap<float> >(elecId.getParameter<edm::InputTag>("src"));
    eidCutValue_ = elecId.getParameter<double>("cutValue");
  }
  if (cfg.exists("jetCorrector")) {
    jetCorrector_ = iC.esConsumes(edm::ESInputTag("", cfg.getParameter<std::string>("jetCorrector")));
  }
  if (cfg.existsAs<edm::ParameterSet>("jetBTagger")) {
    edm::ParameterSet jetBTagger = cfg.getParameter<edm::ParameterSet>("jetBTagger");
    btagLabel_ = iC.consumes<reco::JetTagCollection>(jetBTagger.getParameter<edm::InputTag>("label"));
    btagWorkingPoint_ = jetBTagger.getParameter<double>("workingPoint");
  }
  if (cfg.existsAs<edm::ParameterSet>("jetID")) {
    edm::ParameterSet jetID = cfg.getParameter<edm::ParameterSet>("jetID");
    jetIDLabel_ = iC.consumes<reco::JetIDValueMap>(jetID.getParameter<edm::InputTag>("label"));
    jetIDSelect_ = new StringCutObjectSelector<reco::JetID>(jetID.getParameter<std::string>("select"));
  }
}

template <typename Object>
bool SelectionStep<Object>::select(const edm::Event& event) {
  // fetch input collection
  edm::Handle<edm::View<Object> > src;
  if (!event.getByToken(src_, src))
    return false;

  // load electronId value map if configured such
  edm::Handle<edm::ValueMap<float> > electronId;
  if (!electronId_.isUninitialized()) {
    if (!event.getByToken(electronId_, electronId))
      return false;
  }

  // determine multiplicity of selected objects
  int n = 0;
  for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
    // special treatment for electrons
    if (dynamic_cast<const reco::GsfElectron*>(&*obj)) {
      unsigned int idx = obj - src->begin();
      if (electronId_.isUninitialized() ? true : ((double)(*electronId)[src->refAt(idx)] >= eidCutValue_)) {
        if (select_(*obj))
          ++n;
      }
    }
    // normal treatment
    else {
      if (select_(*obj))
        ++n;
    }
  }
  bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
  return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
}

template <typename Object>
bool SelectionStep<Object>::select(const edm::Event& event, const std::string& type) {
  // fetch input collection
  edm::Handle<edm::View<Object> > src;
  if (!event.getByToken(src_, src))
    return false;

  // special for gsfElectron
  edm::Handle<edm::View<reco::GsfElectron> > elecs_gsf;

  // load electronId value map if configured such
  edm::Handle<edm::ValueMap<float> > electronId;
  if (!electronId_.isUninitialized()) {
    if (!event.getByToken(electronId_, electronId))
      return false;
  }

  // determine multiplicity of selected objects
  int n = 0;
  for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
    // special treatment for PF candidates
    if (dynamic_cast<const reco::PFCandidate*>(&*obj)) {
      reco::PFCandidate objtmp = dynamic_cast<const reco::PFCandidate&>(*obj);

      if (objtmp.muonRef().isNonnull() && type == "muon") {
        if (select_(*obj)) {
          ++n;
        }
      } else if (objtmp.gsfElectronRef().isNonnull() && type == "electron") {
        if (select_(*obj)) {
          if (electronId_.isUninitialized()) {
            ++n;
          } else if (((double)(*electronId)[obj->gsfElectronRef()] >= eidCutValue_)) {
            ++n;
          }
        }
        //        idx_gsf++;
      }
    }

    // special treatment for electrons
    else if (dynamic_cast<const reco::GsfElectron*>(&*obj)) {
      unsigned int idx = obj - src->begin();
      if (electronId_.isUninitialized() ? true : ((double)(*electronId)[src->refAt(idx)] >= eidCutValue_)) {
        if (select_(*obj))
          ++n;
      }
    }

    // normal treatment
    else {
      if (select_(*obj))
        ++n;
    }
  }
  bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
  return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
}

template <typename Object>
bool SelectionStep<Object>::selectVertex(const edm::Event& event) {
  // fetch input collection
  edm::Handle<edm::View<Object> > src;
  if (!event.getByToken(src_, src))
    return false;

  // load electronId value map if configured such
  edm::Handle<edm::ValueMap<float> > electronId;
  if (!electronId_.isUninitialized()) {
    if (!event.getByToken(electronId_, electronId))
      return false;
  }

  // determine multiplicity of selected objects
  int n = 0;
  for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
    if (select_(*obj))
      ++n;
  }
  bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
  return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
}

template <typename Object>
bool SelectionStep<Object>::select(const edm::Event& event, const edm::EventSetup& setup) {
  // fetch input collection
  edm::Handle<edm::View<Object> > src;
  if (!event.getByToken(src_, src))
    return false;

  // load btag collection if configured such
  // NOTE that the JetTagCollection needs an
  // edm::View to reco::Jets; we have to add
  // another Handle bjets for this purpose
  edm::Handle<edm::View<reco::Jet> > bjets;
  edm::Handle<reco::JetTagCollection> btagger;
  edm::Handle<edm::View<reco::Vertex> > pvertex;
  if (!btagLabel_.isUninitialized()) {
    if (!event.getByToken(src_, bjets))
      return false;
    if (!event.getByToken(btagLabel_, btagger))
      return false;
    if (!event.getByToken(pvs_, pvertex))
      return false;
  }

  // load jetID value map if configured such
  edm::Handle<reco::JetIDValueMap> jetID;
  if (jetIDSelect_) {
    if (!event.getByToken(jetIDLabel_, jetID))
      return false;
  }

  // load jet corrector if configured such
  const JetCorrector* corrector = nullptr;
  if (!jetCorrector_.isInitialized() && jetCorrector_.hasValidIndex()) {
    // check whether a jet correcto is in the event setup or not
    if (setup.find(edm::eventsetup::EventSetupRecordKey::makeKey<JetCorrectionsRecord>())) {
      corrector = &setup.getData(jetCorrector_);
    } else {
      edm::LogVerbatim("TopDQMHelpers") << "\n"
                                        << "---------------------------------------------------------------\n"
                                        << " No JetCorrectionsRecord available from EventSetup:\n"
                                        << "  - Jets will not be corrected.\n"
                                        << "  - If you want to change this add the following lines to your "
                                           "cfg file:\n"
                                        << "\n"
                                        << "  ## load jet corrections\n"
                                        << "  process.load(\""
                                           "JetMETCorrections.Configuration.JetCorrectionServicesAllAlgos_cff"
                                           "\")\n"
                                        << "  process.prefer(\"ak5CaloL2L3\")\n"
                                        << "---------------------------------------------------------------"
                                           "\n";
    }
  }
  // determine multiplicity of selected objects
  int n = 0;
  for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
    // check for chosen btag discriminator to be above the
    // corresponding working point if configured such
    unsigned int idx = obj - src->begin();
    if (btagLabel_.isUninitialized() ? true : (*btagger)[bjets->refAt(idx)] > btagWorkingPoint_) {
      bool passedJetID = true;
      // check jetID for calo jets
      if (jetIDSelect_ && dynamic_cast<const reco::CaloJet*>(src->refAt(idx).get())) {
        passedJetID = (*jetIDSelect_)((*jetID)[src->refAt(idx)]);
      }
      if (passedJetID) {
        // scale jet energy if configured such
        Object jet = *obj;
        jet.scaleEnergy(corrector ? corrector->correction(*obj) : 1.);
        if (select_(jet))
          ++n;
      }
    }
  }
  bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
  return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
}

#endif

/* Local Variables: */
/* show-trailing-whitespace: t */
/* truncate-lines: t */
/* End: */