GenFilters/plugins/LHEmttFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    LHEmttFilter
0004 // Class:      LHEmttFilter
0005 //
0006 /*
0007
0008  Description: Filter to select ttbar events with invariant mass over a certain threshold.
0009   (Based on LHEPtFilter)
0010  */
0011
0012 // system include files
0013 #include <memory>
0014 #include <iostream>
0015 #include <set>
0016
0017 // user include files
0018 #include "Math/Vector4D.h"
0019 #include "Math/Vector4Dfwd.h"
0020
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/global/EDFilter.h"
0023
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0029
0030 //
0031 // class declaration
0032 //
0033 class LHEmttFilter : public edm::global::EDFilter<> {
0034 public:
0035   explicit LHEmttFilter(const edm::ParameterSet&);
0036   ~LHEmttFilter() override;
0037
0038   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0039
0040 private:
0041   // ----------member data ---------------------------
0042
0043   edm::EDGetTokenT<LHEEventProduct> src_;
0044   double ptMin_ = 0;
0045   double MinInvMass_ = -1;
0046   double MaxInvMass_ = -1;
0047 };
0048
0049 using namespace edm;
0050 using namespace std;
0051
0052 LHEmttFilter::LHEmttFilter(const edm::ParameterSet& iConfig)
0053     : ptMin_(iConfig.getParameter<double>("ptMin")),
0054       MinInvMass_(iConfig.getParameter<double>("MinInvMass")),
0055       MaxInvMass_(iConfig.getParameter<double>("MaxInvMass")) {
0056   src_ = consumes<LHEEventProduct>(iConfig.getParameter<edm::InputTag>("src"));
0057 }
0058
0059 LHEmttFilter::~LHEmttFilter() {
0060   // do anything here that needs to be done at destruction time
0061   // (e.g. close files, deallocate resources etc.)
0062 }
0063
0064 // ------------ method called to skim the data  ------------
0065 bool LHEmttFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0066   edm::Handle<LHEEventProduct> EvtHandle;
0067   iEvent.getByToken(src_, EvtHandle);
0068
0069   std::vector<lhef::HEPEUP::FiveVector> lheParticles = EvtHandle->hepeup().PUP;
0070   std::vector<ROOT::Math::PxPyPzMVector> cands;
0071   std::vector<int> pdgId_cands;
0072
0073   for (unsigned int i = 0; i < lheParticles.size(); i++) {
0074     if (EvtHandle->hepeup().ISTUP[i] != 2) {  // keep only intermediate particles
0075       continue;
0076     }
0077     int pdgId = EvtHandle->hepeup().IDUP[i];
0078     if (std::abs(pdgId) != 6) {  // keep only top quarks
0079       continue;
0080     }
0081     pdgId_cands.push_back(pdgId);
0082     cands.push_back(
0083         ROOT::Math::PxPyPzMVector(lheParticles[i][0], lheParticles[i][1], lheParticles[i][2], lheParticles[i][4]));
0084   }
0085
0086   if (cands.size() != 2) {
0087     edm::LogWarning("LHEmttFilter Error") << "Number of top quarks found in the event != 2" << endl;
0088   }
0089
0090   double ttmass_ = -1;
0091   if (!cands.empty() && (cands.size() == 2)) {                           //two top quarks
0092     if (pdgId_cands.at(0) + pdgId_cands.at(1) == 0) {                    //exactly one t and one tbar
0093       if ((cands.at(0).pt() > ptMin_) && (cands.at(1).pt() > ptMin_)) {  //requiring minimum pT
0094         ROOT::Math::PxPyPzMVector tot = cands.at(0);
0095         for (unsigned icand = 1; icand < cands.size(); ++icand) {
0096           tot += cands.at(icand);
0097         }
0098         ttmass_ = tot.mass();
0099       }
0100     } else if (pdgId_cands.at(0) + pdgId_cands.at(1) != 0) {
0101       edm::LogWarning("LHEmttFilter Error") << "Found two t/tbar quarks instead of a ttbar pair" << endl;
0102     }
0103   }
0104
0105   if ((MinInvMass_ > -1 && ttmass_ < MinInvMass_) || (MaxInvMass_ > -1 && ttmass_ > MaxInvMass_)) {
0106     return false;
0107   } else {
0108     return true;
0109   }
0110 }
0111
0112 //define this as a plug-in
0113 DEFINE_FWK_MODULE(LHEmttFilter);