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 "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_;
  std::string 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_(0) {

  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_ = 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 = 0;
  if (!jetCorrector_.empty()) {
    // check whether a jet correcto is in the event setup or not
    if (setup.find(edm::eventsetup::EventSetupRecordKey::makeKey<
            JetCorrectionsRecord>())) {
      corrector = JetCorrector::getJetCorrector(jetCorrector_, setup);
    } 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: */