HepMCCandAlgos/plugins/FlavorHistoryFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FlavorHistoryFilter
0004 // Class:      FlavorHistoryFilter
0005 //
0006 /**\class FlavorHistoryFilter FlavorHistoryFilter.cc PhysicsTools/FlavorHistoryFilter/src/FlavorHistoryFilter.cc
0007
0008  Description:
0009
0010  This now filters events hierarchically. Previously this was done at the python configuration
0011  level, which was cumbersome for users to use.
0012
0013  Now, the hierarchy is:
0014
0015  Create prioritized paths to separate HF composition samples.
0016
0017  These are exclusive priorities, so sample "i" will not overlap with "i+1".
0018  Note that the "dr" values below correspond to the dr between the
0019  matched genjet, and the sister genjet.
0020
0021  1) W+bb with >= 2 jets from the ME (dr > 0.5)
0022  2) W+b or W+bb with 1 jet from the ME
0023  3) W+cc from the ME (dr > 0.5)
0024  4) W+c or W+cc with 1 jet from the ME
0025  5) W+bb with 1 jet from the parton shower (dr == 0.0)
0026  6) W+cc with 1 jet from the parton shower (dr == 0.0)
0027
0028  These are the "trash bin" samples that we're throwing away:
0029
0030  7) W+bb with >= 2 partons but 1 jet from the ME (dr == 0.0)
0031  8) W+cc with >= 2 partons but 1 jet from the ME (dr == 0.0)
0032  9) W+bb with >= 2 partons but 2 jets from the PS (dr > 0.5)
0033  10)W+cc with >= 2 partons but 2 jets from the PS (dr > 0.5)
0034
0035  And here is the true "light flavor" sample:
0036
0037  11) Veto of all the previous (W+ light jets)
0038
0039  Implementation:
0040      <Notes on implementation>
0041 */
0042 //
0043 // Original Author:  "Salvatore Rappoccio"
0044 //         Created:  Sat Jun 28 00:41:21 CDT 2008
0045 //
0046 //
0047
0048 // system include files
0049 #include <memory>
0050
0051 // user include files
0052 #include "FWCore/Framework/interface/Frameworkfwd.h"
0053 #include "FWCore/Framework/interface/global/EDFilter.h"
0054
0055 #include "FWCore/Framework/interface/Event.h"
0056 #include "FWCore/Framework/interface/MakerMacros.h"
0057
0058 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0059
0060 #include "DataFormats/HepMCCandidate/interface/FlavorHistoryEvent.h"
0061 #include "DataFormats/HepMCCandidate/interface/FlavorHistory.h"
0062 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0063
0064 #include "PhysicsTools/HepMCCandAlgos/interface/FlavorHistorySelectorUtil.h"
0065
0066 #include "FWCore/Utilities/interface/Exception.h"
0067 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0068 #include "DataFormats/Math/interface/deltaR.h"
0069
0070 #include "PhysicsTools/HepMCCandAlgos/interface/FlavorHistorySelectorUtil.h"
0071
0072 #include <vector>
0073
0074 //
0075 // class declaration
0076 //
0077
0078 class FlavorHistoryFilter : public edm::global::EDFilter<> {
0079 public:
0080   typedef reco::FlavorHistory::FLAVOR_T flavor_type;
0081   typedef std::vector<int> flavor_vector;
0082
0083   explicit FlavorHistoryFilter(const edm::ParameterSet&);
0084
0085   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0086
0087 private:
0088   bool filter(edm::StreamID, ::edm::Event&, const edm::EventSetup&) const override;
0089
0090   // ----------member data ---------------------------
0091   edm::EDGetTokenT<reco::FlavorHistoryEvent> bsrcToken_;  // Input b flavor history collection name
0092   edm::EDGetTokenT<reco::FlavorHistoryEvent> csrcToken_;  // Input c flavor history collection name
0093   int pathToSelect_;                                      // Select any of the following paths:
0094   double dr_;                                             // dr with which to cut off the events
0095   // Note! The "b" and "c" here refer to the number of matched b and c genjets, respectively
0096   std::unique_ptr<reco::FlavorHistorySelectorUtil const> bb_me_;  // To select bb->2 events from matrix element... Path 1
0097   std::unique_ptr<reco::FlavorHistorySelectorUtil const> b_me_;  // To select  b->1 events from matrix element... Path 2
0098   std::unique_ptr<reco::FlavorHistorySelectorUtil const> cc_me_;  // To select cc->2 events from matrix element... Path 3
0099   std::unique_ptr<reco::FlavorHistorySelectorUtil const> c_me_;  // To select  c->1 events from matrix element... Path 4
0100   std::unique_ptr<reco::FlavorHistorySelectorUtil const> b_ps_;  // To select bb->2 events from parton shower ... Path 5
0101   std::unique_ptr<reco::FlavorHistorySelectorUtil const> c_ps_;  // To select cc->2 events from parton shower ... Path 6
0102   std::unique_ptr<reco::FlavorHistorySelectorUtil const>
0103       bb_me_comp_;  // To select bb->1 events from matrix element... Path 7
0104   std::unique_ptr<reco::FlavorHistorySelectorUtil const>
0105       cc_me_comp_;  // To select cc->1 events from matrix element... Path 8
0106   std::unique_ptr<reco::FlavorHistorySelectorUtil const>
0107       b_ps_comp_;  // To select bb->2 events from parton shower ... Path 9
0108   std::unique_ptr<reco::FlavorHistorySelectorUtil const>
0109       c_ps_comp_;  // To select cc->1 events from parton shower ... Path 10
0110                    // The veto of all of these is               ... Path 11
0111 };
0112
0113 using namespace edm;
0114 using namespace reco;
0115 using namespace std;
0116
0117 //
0118 // constants, enums and typedefs
0119 //
0120
0121 //
0122 // static data member definitions
0123 //
0124
0125 //
0126 // constructors and destructor
0127 //
0128 FlavorHistoryFilter::FlavorHistoryFilter(const edm::ParameterSet& iConfig)
0129     : bsrcToken_(consumes<FlavorHistoryEvent>(iConfig.getParameter<edm::InputTag>("bsrc"))),
0130       csrcToken_(consumes<FlavorHistoryEvent>(iConfig.getParameter<edm::InputTag>("csrc"))) {
0131   pathToSelect_ = iConfig.getParameter<int>("pathToSelect");
0132
0133   // This is the "interface" delta R with which to decide
0134   // where to take the event from
0135   dr_ = iConfig.getParameter<double>("dr");
0136   bool verbose(iConfig.getParameter<bool>("verbose"));
0137
0138   // Set up the boundaries.
0139   // dr0 = 0.0
0140   // dr1 = set by user
0141   // dr2 = infinity
0142   double dr0 = 0.0;
0143   double dr1 = dr_;
0144   double dr2 = 99999.0;
0145
0146   // These are the processes that can come from the matrix element calculation
0147   std::vector<int> me_ids;
0148   me_ids.push_back(2);  // flavor excitation
0149   me_ids.push_back(3);  // flavor creation
0150
0151   // These are the processes that can come from the parton shower calculation
0152   std::vector<int> ps_ids;
0153   ps_ids.push_back(1);  // gluon splitting
0154
0155   // To select bb->2 events from matrix element... Path 1
0156   bb_me_ = std::make_unique<FlavorHistorySelectorUtil>(5, 2, me_ids, dr1, dr2, verbose);
0157
0158   // To select  b->1 events from matrix element... Path 2
0159   b_me_ = std::make_unique<FlavorHistorySelectorUtil>(5, 1, me_ids, dr0, dr0, verbose);
0160
0161   // To select cc->2 events from matrix element... Path 3
0162   cc_me_ = std::make_unique<FlavorHistorySelectorUtil>(4, 2, me_ids, dr1, dr2, verbose);
0163
0164   // To select  c->1 events from matrix element... Path 4
0165   c_me_ = std::make_unique<FlavorHistorySelectorUtil>(4, 1, me_ids, dr0, dr0, verbose);
0166
0167   // To select bb->2 events from parton shower ... Path 5
0168   b_ps_ = std::make_unique<FlavorHistorySelectorUtil>(5, 1, ps_ids, dr0, dr1, verbose);
0169
0170   // To select cc->2 events from parton shower ... Path 6
0171   c_ps_ = std::make_unique<FlavorHistorySelectorUtil>(4, 1, ps_ids, dr0, dr1, verbose);
0172
0173   // To select bb->1 events from matrix element... Path 7
0174   bb_me_comp_ = std::make_unique<FlavorHistorySelectorUtil>(5, 2, me_ids, dr0, dr1, verbose);
0175
0176   // To select cc->1 events from matrix element... Path 8
0177   cc_me_comp_ = std::make_unique<FlavorHistorySelectorUtil>(4, 2, me_ids, dr0, dr1, verbose);
0178
0179   // To select bb->2 events from parton shower ... Path 9
0180   b_ps_comp_ = std::make_unique<FlavorHistorySelectorUtil>(5, 2, ps_ids, dr1, dr2, verbose);
0181
0182   // To select cc->1 events from parton shower ... Path 10
0183   c_ps_comp_ = std::make_unique<FlavorHistorySelectorUtil>(4, 2, ps_ids, dr1, dr2, verbose);
0184
0185   // The veto of all of these is               ... Path 11
0186
0187   // This will write 1-11 (the path number), or 0 if error.
0188   produces<unsigned int>();
0189 }
0190
0191 void FlavorHistoryFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0192   edm::ParameterSetDescription desc;
0193
0194   desc.add<edm::InputTag>("bsrc");
0195   desc.add<edm::InputTag>("csrc");
0196   desc.add<int>("pathToSelect", -1);
0197   desc.add<double>("dr");
0198   desc.add<bool>("verbose");
0199
0200   descriptions.addDefault(desc);
0201 }
0202
0203 bool FlavorHistoryFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0204   // Get the flavor history
0205   Handle<FlavorHistoryEvent> bFlavorHistoryEvent;
0206   iEvent.getByToken(bsrcToken_, bFlavorHistoryEvent);
0207
0208   Handle<FlavorHistoryEvent> cFlavorHistoryEvent;
0209   iEvent.getByToken(csrcToken_, cFlavorHistoryEvent);
0210
0211   auto selection = std::make_unique<unsigned int>(0);
0212
0213   // Get the number of matched b-jets in the event
0214   unsigned int nb = bFlavorHistoryEvent->nb();
0215   // Get the number of matched c-jets in the event
0216   unsigned int nc = cFlavorHistoryEvent->nc();
0217   // Get the two flavor sources. The highest takes precedence
0218   // over the rest.
0219   FlavorHistory::FLAVOR_T bFlavorSource = bFlavorHistoryEvent->flavorSource();
0220   FlavorHistory::FLAVOR_T cFlavorSource = cFlavorHistoryEvent->flavorSource();
0221   FlavorHistory::FLAVOR_T flavorSource = FlavorHistory::FLAVOR_NULL;
0222   // Get the highest flavor in the event
0223   unsigned int highestFlavor = 0;
0224   // Get the delta r between the two heavy flavor matched jets.
0225   double dr = -1;
0226
0227   // Preference is in increasing priority:
0228   //  1: gluon splitting
0229   //  2: flavor excitation
0230   //  3: flavor creation (matrix element)
0231   //  4: flavor decay
0232   if (bFlavorSource >= cFlavorSource) {
0233     flavorSource = bFlavorHistoryEvent->flavorSource();
0234     highestFlavor = bFlavorHistoryEvent->highestFlavor();
0235     dr = bFlavorHistoryEvent->deltaR();
0236   } else {
0237     flavorSource = cFlavorHistoryEvent->flavorSource();
0238     highestFlavor = cFlavorHistoryEvent->highestFlavor();
0239     dr = cFlavorHistoryEvent->deltaR();
0240   }
0241
0242   *selection = 0;
0243   // Now make hierarchical determination
0244   if (bb_me_->select(nb, nc, highestFlavor, flavorSource, dr))
0245     *selection = 1;
0246   else if (b_me_->select(nb, nc, highestFlavor, flavorSource, dr))
0247     *selection = 2;
0248   else if (cc_me_->select(nb, nc, highestFlavor, flavorSource, dr))
0249     *selection = 3;
0250   else if (c_me_->select(nb, nc, highestFlavor, flavorSource, dr))
0251     *selection = 4;
0252   else if (b_ps_->select(nb, nc, highestFlavor, flavorSource, dr))
0253     *selection = 5;
0254   else if (c_ps_->select(nb, nc, highestFlavor, flavorSource, dr))
0255     *selection = 6;
0256   else if (bb_me_comp_->select(nb, nc, highestFlavor, flavorSource, dr))
0257     *selection = 7;
0258   else if (cc_me_comp_->select(nb, nc, highestFlavor, flavorSource, dr))
0259     *selection = 8;
0260   else if (b_ps_comp_->select(nb, nc, highestFlavor, flavorSource, dr))
0261     *selection = 9;
0262   else if (c_ps_comp_->select(nb, nc, highestFlavor, flavorSource, dr))
0263     *selection = 10;
0264   else
0265     *selection = 11;
0266
0267   bool pass = false;
0268   if (pathToSelect_ > 0) {
0269     pass = (*selection > 0 && *selection == static_cast<unsigned int>(pathToSelect_));
0270   } else {
0271     pass = true;
0272   }
0273
0274   iEvent.put(std::move(selection));
0275
0276   return pass;
0277 }
0278
0279 //define this as a plug-in
0280 DEFINE_FWK_MODULE(FlavorHistoryFilter);