Physics/interface/TopDQMHelpers.h

0001 #ifndef TOPDQMHELPERS
0002 #define TOPDQMHELPERS
0003
0004 #include <string>
0005 #include <vector>
0006 #include <iostream>
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Common/interface/TriggerNames.h"
0009 #include "DataFormats/Common/interface/TriggerResults.h"
0010 #include "DataFormats/VertexReco/interface/Vertex.h"
0011 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0012 #include <DataFormats/METReco/interface/PFMET.h>
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "FWCore/Framework/interface/EDConsumerBase.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "DataFormats/PatCandidates/interface/Muon.h"
0018 #include "DataFormats/PatCandidates/interface/Electron.h"
0019 #include "DataFormats/PatCandidates/interface/Jet.h"
0020 #include "DataFormats/PatCandidates/interface/MET.h"
0021
0022 /**
0023    \fn      accept TopDQMHelpers.h "DQM/Physics/interface/TopDQMHelpers.h"
0024
0025    \brief   Helper function to determine trigger accepts.
0026
0027    Helper function to determine trigger accept for given TriggerResults and
0028    a given TriggerPath(s).
0029 */
0030
0031 inline bool accept(const edm::Event& event, const edm::TriggerResults& triggerTable, const std::string& triggerPath) {
0032   bool passed = false;
0033   const edm::TriggerNames& triggerNames = event.triggerNames(triggerTable);
0034   for (unsigned int i = 0; i < triggerNames.triggerNames().size(); ++i) {
0035     if (triggerNames.triggerNames()[i] == triggerPath) {
0036       if (triggerTable.accept(i)) {
0037         passed = true;
0038         break;
0039       }
0040     }
0041   }
0042   return passed;
0043 }
0044
0045 inline bool accept(const edm::Event& event,
0046                    const edm::TriggerResults& triggerTable,
0047                    const std::vector<std::string>& triggerPaths) {
0048   bool passed = false;
0049   for (unsigned int j = 0; j < triggerPaths.size(); ++j) {
0050     if (accept(event, triggerTable, triggerPaths[j])) {
0051       passed = true;
0052       break;
0053     }
0054   }
0055   return passed;
0056 }
0057
0058 #include "DataFormats/JetReco/interface/Jet.h"
0059 #include "FWCore/Framework/interface/EventSetup.h"
0060
0061 /**
0062    \class   Calculate TopDQMHelpers.h "DQM/Physics/interface/TopDQMHelpers.h"
0063
0064    \brief   Helper class for the calculation of a top and a W boson mass estime.
0065
0066    Helper class for the calculation of a top and a W boson mass estimate. The
0067    core implementation originates from the plugin TtSemiLepHypMaxSumPtWMass
0068    in TopQuarkAnalysis/TopJetCombination package. It may be extended to include
0069    b tag information.
0070 */
0071 class Calculate_miniAOD {
0072 public:
0073   /// default constructor
0074   Calculate_miniAOD(int maxNJets, double wMass);
0075   /// default destructor
0076   ~Calculate_miniAOD(){};
0077
0078   /// calculate W boson mass estimate
0079   double massWBoson(const std::vector<pat::Jet>& jets);
0080   /// calculate t-quark mass estimate
0081   double massTopQuark(const std::vector<pat::Jet>& jets);
0082   /// calculate b-tagged t-quark mass estimate
0083   // double massBTopQuark(const std::vector<reco::Jet>& jets, std::vector<bool>
0084   // bjet);
0085   double massBTopQuark(const std::vector<pat::Jet>& jets, std::vector<double> VbtagWP, double btagWP_);
0086
0087   /// calculate W boson transverse mass estimate
0088   double tmassWBoson(pat::Muon* lep, const pat::MET& met, const pat::Jet& b);
0089   /// calculate top quark transverse mass estimate
0090   double tmassTopQuark(pat::Muon* lep, const pat::MET& met, const pat::Jet& b);
0091
0092   /// calculate W boson transverse mass estimate
0093   double tmassWBoson(pat::Electron* lep, const pat::MET& met, const pat::Jet& b);
0094   /// calculate top quark transverse mass estimate
0095   double tmassTopQuark(pat::Electron* lep, const pat::MET& met, const pat::Jet& b);
0096
0097 private:
0098   /// do the calculation; this is called only once per event by the first
0099   /// function call to return a mass estimate. The once calculated values
0100   /// are cached afterwards
0101   void operator()(const std::vector<pat::Jet>& jets);
0102   /// do the calculation of the t-quark mass with one b-jet
0103   void operator2(const std::vector<pat::Jet>&, std::vector<double>, double);
0104   /// do the calculation of the transverse top and W masses
0105   void operator()(const pat::Jet& bJet, pat::Electron* lepton, const pat::MET& met);
0106
0107   void operator()(const pat::Jet& bJet, pat::Muon* lepton, const pat::MET& met);
0108
0109 private:
0110   /// indicate failed associations
0111   bool failed_;
0112   /// max. number of jets to be considered
0113   int maxNJets_;
0114   /// paramater of the w boson mass
0115   double wMass_;
0116   /// cache of w boson mass estimate
0117   double massWBoson_;
0118   /// cache of top quark mass estimate
0119   double massTopQuark_;
0120   /// cache of b-tagged top quark mass estimate
0121   double massBTopQuark_;
0122   /// cache of W boson transverse mass estimate
0123   double tmassWBoson_;
0124   /// cache of top quark transverse mass estimate
0125   double tmassTopQuark_;
0126 };
0127
0128 class Calculate {
0129 public:
0130   /// default constructor
0131   Calculate(int maxNJets, double wMass);
0132   /// default destructor
0133   ~Calculate(){};
0134
0135   /// calculate W boson mass estimate
0136   double massWBoson(const std::vector<reco::Jet>& jets);
0137   /// calculate t-quark mass estimate
0138   double massTopQuark(const std::vector<reco::Jet>& jets);
0139   /// calculate b-tagged t-quark mass estimate
0140   // double massBTopQuark(const std::vector<reco::Jet>& jets, std::vector<bool>
0141   // bjet);
0142   double massBTopQuark(const std::vector<reco::Jet>& jets, std::vector<double> VbtagWP, double btagWP_);
0143
0144   /// calculate W boson transverse mass estimate
0145   double tmassWBoson(reco::RecoCandidate* lep, const reco::MET& met, const reco::Jet& b);
0146   /// calculate top quark transverse mass estimate
0147   double tmassTopQuark(reco::RecoCandidate* lep, const reco::MET& met, const reco::Jet& b);
0148
0149 private:
0150   /// do the calculation; this is called only once per event by the first
0151   /// function call to return a mass estimate. The once calculated values
0152   /// are cached afterwards
0153   void operator()(const std::vector<reco::Jet>& jets);
0154   /// do the calculation of the t-quark mass with one b-jet
0155   void operator2(const std::vector<reco::Jet>&, std::vector<double>, double);
0156   /// do the calculation of the transverse top and W masses
0157   void operator()(const reco::Jet& bJet, reco::RecoCandidate* lepton, const reco::MET& met);
0158
0159 private:
0160   /// indicate failed associations
0161   bool failed_;
0162   /// max. number of jets to be considered
0163   int maxNJets_;
0164   /// paramater of the w boson mass
0165   double wMass_;
0166   /// cache of w boson mass estimate
0167   double massWBoson_;
0168   /// cache of top quark mass estimate
0169   double massTopQuark_;
0170   /// cache of b-tagged top quark mass estimate
0171   double massBTopQuark_;
0172   /// cache of W boson transverse mass estimate
0173   double tmassWBoson_;
0174   /// cache of top quark transverse mass estimate
0175   double tmassTopQuark_;
0176 };
0177
0178 #include "DataFormats/JetReco/interface/JetID.h"
0179 #include "DataFormats/JetReco/interface/PFJet.h"
0180 #include "DataFormats/JetReco/interface/CaloJet.h"
0181 #include "DataFormats/BTauReco/interface/JetTag.h"
0182 #include "DataFormats/Common/interface/Handle.h"
0183 #include "DataFormats/Common/interface/ValueMap.h"
0184 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0185 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0186 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0187 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0188 #include "FWCore/Utilities/interface/EDGetToken.h"
0189 #include "FWCore/Utilities/interface/InputTag.h"
0190
0191 /**
0192    \class   SelectionStep TopDQMHelpers.h
0193    "DQM/Physics/interface/TopDQMHelpers.h"
0194
0195    \brief   Templated helper class to allow a selection on a certain object
0196    collection.
0197
0198    Templated helper class to allow a selection on a certain object collection,
0199    which may be monitored by a separate class afterwards. The class wraps and
0200    slightly extends the features of the StringCutParser to allow also to apply
0201    event based selections, according to a minimal or maximal number of elements
0202    in the collection after the object selection has been applied. It takes an
0203    edm::ParameterSet in the constructor, which should contain the following
0204    elements:
0205     - src          : the input collection (mandatory).
0206     - select       : the selection string (mandatory).
0207     - min          : whether there is a min value on which to reject the whole
0208    event after
0209                      the selection (optional).
0210     - max          : whether there is a max value on which to reject the whole
0211    event after
0212                      the selection (optional).
0213     - electronId   : input tag of an electronId association map and cut value
0214                      (optional).
0215     - jetCorrector : label of jet corrector (optional).
0216     - jetBTagger   : parameters defining the btag algorithm and working point of
0217    choice
0218                      (optional).
0219     - jetID        : parameters defining the jetID value map and selection
0220    (optional).
0221
0222
0223    The parameters _src_ and _select_ are mandatory. The parameters _min_ and
0224    _max_ are optional. The parameters _electronId_ and _jetCorrector_ are
0225    optional. They are added to keep the possibility to apply selections on id'ed
0226    electrons or on corrected jets.  They may be omitted in the PSet for
0227    simplification reasons if not needed at any time.  They are not effiective
0228    for other object collections but electrons or jets.  If none of the two
0229    parameters _min_ or _max_ is found in the event the select function returns
0230    true if at least one object fullfilled the requirements.
0231
0232    The class has one template value, which is the object collection to apply the
0233    selection on. This has to be parsed to the StringCutParser class. The
0234    function select is overrided for jets to circumvent problems with the
0235    template specialisation. Note that for MET not type1 or muon corrections are
0236    supported on reco candidates.
0237 */
0238
0239 template <typename Object>
0240 class SelectionStep {
0241 public:
0242   /// default constructor
0243   SelectionStep(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC);
0244   /// default destructor
0245   ~SelectionStep(){};
0246
0247   /// apply selection
0248   bool select(const edm::Event& event);
0249   bool select(const edm::Event& event, const std::string& type);
0250   /// apply selection override for jets
0251   bool select(const edm::Event& event, const edm::EventSetup& setup);
0252   bool selectVertex(const edm::Event& event);
0253
0254 private:
0255   /// input collection
0256   edm::EDGetTokenT<edm::View<Object> > src_;
0257   /// min/max for object multiplicity
0258   int min_, max_;
0259   /// electronId label as extra selection type
0260   edm::EDGetTokenT<edm::ValueMap<float> > electronId_;
0261   /// electronId pattern we expect the following pattern:
0262   ///  0: fails
0263   ///  1: passes electron ID only
0264   ///  2: passes electron Isolation only
0265   ///  3: passes electron ID and Isolation only
0266   ///  4: passes conversion rejection
0267   ///  5: passes conversion rejection and ID
0268   ///  6: passes conversion rejection and Isolation
0269   ///  7: passes the whole selection
0270   /// As described on
0271   /// https://twiki.cern.ch/twiki/bin/view/CMS/SimpleCutBasedEleID
0272   double eidCutValue_;
0273   /// jet corrector as extra selection type
0274   edm::EDGetTokenT<reco::JetCorrector> jetCorrector_;
0275   /// choice for b-tag as extra selection type
0276   edm::EDGetTokenT<reco::JetTagCollection> btagLabel_;
0277   /// choice of b-tag working point as extra selection type
0278   double btagWorkingPoint_;
0279   /// jetID as an extra selection type
0280   edm::EDGetTokenT<reco::JetIDValueMap> jetIDLabel_;
0281
0282   edm::EDGetTokenT<edm::View<reco::Vertex> > pvs_;
0283   edm::EDGetTokenT<edm::View<reco::GsfElectron> > gsfEs_;
0284   /// string cut selector
0285   StringCutObjectSelector<Object> select_;
0286   /// selection string on the jetID
0287   StringCutObjectSelector<reco::JetID>* jetIDSelect_;
0288 };
0289
0290 /// default constructor
0291 template <typename Object>
0292 SelectionStep<Object>::SelectionStep(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
0293     : select_(cfg.getParameter<std::string>("select")), jetIDSelect_(nullptr) {
0294   src_ = iC.consumes<edm::View<Object> >(cfg.getParameter<edm::InputTag>("src"));
0295   // exist otherwise they are initialized with -1
0296   cfg.exists("min") ? min_ = cfg.getParameter<int>("min") : min_ = -1;
0297   cfg.exists("max") ? max_ = cfg.getParameter<int>("max") : max_ = -1;
0298   std::string mygSF = "gedGsfElectrons";
0299   gsfEs_ = iC.consumes<edm::View<reco::GsfElectron> >(cfg.getUntrackedParameter<edm::InputTag>("myGSF", mygSF));
0300   if (cfg.existsAs<edm::ParameterSet>("electronId")) {
0301     edm::ParameterSet elecId = cfg.getParameter<edm::ParameterSet>("electronId");
0302     electronId_ = iC.consumes<edm::ValueMap<float> >(elecId.getParameter<edm::InputTag>("src"));
0303     eidCutValue_ = elecId.getParameter<double>("cutValue");
0304   }
0305   if (cfg.exists("jetCorrector")) {
0306     jetCorrector_ = iC.consumes<reco::JetCorrector>(edm::InputTag(cfg.getParameter<std::string>("jetCorrector")));
0307   }
0308   if (cfg.existsAs<edm::ParameterSet>("jetBTagger")) {
0309     edm::ParameterSet jetBTagger = cfg.getParameter<edm::ParameterSet>("jetBTagger");
0310     btagLabel_ = iC.consumes<reco::JetTagCollection>(jetBTagger.getParameter<edm::InputTag>("label"));
0311     btagWorkingPoint_ = jetBTagger.getParameter<double>("workingPoint");
0312   }
0313   if (cfg.existsAs<edm::ParameterSet>("jetID")) {
0314     edm::ParameterSet jetID = cfg.getParameter<edm::ParameterSet>("jetID");
0315     jetIDLabel_ = iC.consumes<reco::JetIDValueMap>(jetID.getParameter<edm::InputTag>("label"));
0316     jetIDSelect_ = new StringCutObjectSelector<reco::JetID>(jetID.getParameter<std::string>("select"));
0317   }
0318 }
0319
0320 /// apply selection
0321 template <typename Object>
0322 bool SelectionStep<Object>::select(const edm::Event& event) {
0323   // fetch input collection
0324   edm::Handle<edm::View<Object> > src;
0325   if (!event.getByToken(src_, src))
0326     return false;
0327
0328   // load electronId value map if configured such
0329   edm::Handle<edm::ValueMap<float> > electronId;
0330   if (!electronId_.isUninitialized()) {
0331     if (!event.getByToken(electronId_, electronId))
0332       return false;
0333   }
0334
0335   // determine multiplicity of selected objects
0336   int n = 0;
0337   for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
0338     // special treatment for electrons
0339     if (dynamic_cast<const reco::GsfElectron*>(&*obj)) {
0340       unsigned int idx = obj - src->begin();
0341       if (electronId_.isUninitialized() ? true : ((double)(*electronId)[src->refAt(idx)] >= eidCutValue_)) {
0342         if (select_(*obj))
0343           ++n;
0344       }
0345     }
0346     // normal treatment
0347     else {
0348       if (select_(*obj))
0349         ++n;
0350     }
0351   }
0352   bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
0353   return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
0354 }
0355
0356 /// apply selection with special treatment for PFCandidates
0357 template <typename Object>
0358 bool SelectionStep<Object>::select(const edm::Event& event, const std::string& type) {
0359   // fetch input collection
0360   edm::Handle<edm::View<Object> > src;
0361   if (!event.getByToken(src_, src))
0362     return false;
0363
0364   // special for gsfElectron
0365   edm::Handle<edm::View<reco::GsfElectron> > elecs_gsf;
0366
0367   // load electronId value map if configured such
0368   edm::Handle<edm::ValueMap<float> > electronId;
0369   if (!electronId_.isUninitialized()) {
0370     if (!event.getByToken(electronId_, electronId))
0371       return false;
0372   }
0373
0374   // determine multiplicity of selected objects
0375   int n = 0;
0376   for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
0377     // special treatment for PF candidates
0378     if (dynamic_cast<const reco::PFCandidate*>(&*obj)) {
0379       reco::PFCandidate objtmp = dynamic_cast<const reco::PFCandidate&>(*obj);
0380
0381       if (objtmp.muonRef().isNonnull() && type == "muon") {
0382         if (select_(*obj)) {
0383           ++n;
0384         }
0385       } else if (objtmp.gsfElectronRef().isNonnull() && type == "electron") {
0386         if (select_(*obj)) {
0387           if (electronId_.isUninitialized()) {
0388             ++n;
0389           } else if (((double)(*electronId)[obj->gsfElectronRef()] >= eidCutValue_)) {
0390             ++n;
0391           }
0392         }
0393         //        idx_gsf++;
0394       }
0395     }
0396
0397     // special treatment for electrons
0398     else if (dynamic_cast<const reco::GsfElectron*>(&*obj)) {
0399       unsigned int idx = obj - src->begin();
0400       if (electronId_.isUninitialized() ? true : ((double)(*electronId)[src->refAt(idx)] >= eidCutValue_)) {
0401         if (select_(*obj))
0402           ++n;
0403       }
0404     }
0405
0406     // normal treatment
0407     else {
0408       if (select_(*obj))
0409         ++n;
0410     }
0411   }
0412   bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
0413   return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
0414 }
0415
0416 template <typename Object>
0417 bool SelectionStep<Object>::selectVertex(const edm::Event& event) {
0418   // fetch input collection
0419   edm::Handle<edm::View<Object> > src;
0420   if (!event.getByToken(src_, src))
0421     return false;
0422
0423   // load electronId value map if configured such
0424   edm::Handle<edm::ValueMap<float> > electronId;
0425   if (!electronId_.isUninitialized()) {
0426     if (!event.getByToken(electronId_, electronId))
0427       return false;
0428   }
0429
0430   // determine multiplicity of selected objects
0431   int n = 0;
0432   for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
0433     if (select_(*obj))
0434       ++n;
0435   }
0436   bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
0437   return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
0438 }
0439
0440 /// apply selection (w/o using the template class Object), override for jets
0441 template <typename Object>
0442 bool SelectionStep<Object>::select(const edm::Event& event, const edm::EventSetup& setup) {
0443   // fetch input collection
0444   edm::Handle<edm::View<Object> > src;
0445   if (!event.getByToken(src_, src))
0446     return false;
0447
0448   // load btag collection if configured such
0449   // NOTE that the JetTagCollection needs an
0450   // edm::View to reco::Jets; we have to add
0451   // another Handle bjets for this purpose
0452   edm::Handle<edm::View<reco::Jet> > bjets;
0453   edm::Handle<reco::JetTagCollection> btagger;
0454   edm::Handle<edm::View<reco::Vertex> > pvertex;
0455   if (!btagLabel_.isUninitialized()) {
0456     if (!event.getByToken(src_, bjets))
0457       return false;
0458     if (!event.getByToken(btagLabel_, btagger))
0459       return false;
0460     if (!event.getByToken(pvs_, pvertex))
0461       return false;
0462   }
0463
0464   // load jetID value map if configured such
0465   edm::Handle<reco::JetIDValueMap> jetID;
0466   if (jetIDSelect_) {
0467     if (!event.getByToken(jetIDLabel_, jetID))
0468       return false;
0469   }
0470
0471   // load jet corrector if configured such
0472   const reco::JetCorrector* corrector = nullptr;
0473   if (!jetCorrector_.isUninitialized()) {
0474     // check whether a jet corrector is in the event or not
0475     edm::Handle<reco::JetCorrector> correctorHandle = event.getHandle(jetCorrector_);
0476     if (correctorHandle.isValid()) {
0477       corrector = correctorHandle.product();
0478     } else {
0479       edm::LogVerbatim("TopDQMHelpers") << "\n"
0480                                         << "---------------------------------------------------------------\n"
0481                                         << " No reco::JetCorrrector available from Event:\n"
0482                                         << "  - Jets will not be corrected.\n"
0483                                         << "---------------------------------------------------------------"
0484                                            "\n";
0485     }
0486   }
0487   // determine multiplicity of selected objects
0488   int n = 0;
0489   for (typename edm::View<Object>::const_iterator obj = src->begin(); obj != src->end(); ++obj) {
0490     // check for chosen btag discriminator to be above the
0491     // corresponding working point if configured such
0492     unsigned int idx = obj - src->begin();
0493     if (btagLabel_.isUninitialized() ? true : (*btagger)[bjets->refAt(idx)] > btagWorkingPoint_) {
0494       bool passedJetID = true;
0495       // check jetID for calo jets
0496       if (jetIDSelect_ && dynamic_cast<const reco::CaloJet*>(src->refAt(idx).get())) {
0497         passedJetID = (*jetIDSelect_)((*jetID)[src->refAt(idx)]);
0498       }
0499       if (passedJetID) {
0500         // scale jet energy if configured such
0501         Object jet = *obj;
0502         jet.scaleEnergy(corrector ? corrector->correction(*obj) : 1.);
0503         if (select_(jet))
0504           ++n;
0505       }
0506     }
0507   }
0508   bool accept = (min_ >= 0 ? n >= min_ : true) && (max_ >= 0 ? n <= max_ : true);
0509   return (min_ < 0 && max_ < 0) ? (n > 0) : accept;
0510 }
0511
0512 #endif
0513
0514 /* Local Variables: */
0515 /* show-trailing-whitespace: t */
0516 /* truncate-lines: t */
0517 /* End: */