GenFilters/plugins/LHEPtFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    LHEPtFilter
0004 // Class:      LHEPtFilter
0005 //
0006 /*
0007
0008  Description: Filter to select events with pT in a given range; includes a switch for sum type (vector or scalar).
0009  (Based on LHEGenericFilter)
0010
0011
0012 */
0013 //
0014
0015 // system include files
0016 #include <memory>
0017 #include <iostream>
0018 #include <set>
0019
0020 // user include files
0021 #include "Math/Vector4D.h"
0022 #include "Math/Vector4Dfwd.h"
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/global/EDFilter.h"
0026
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029
0030 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0033
0034 //
0035 // class declaration
0036 //
0037
0038 class LHEPtFilter : public edm::global::EDFilter<> {
0039 public:
0040   explicit LHEPtFilter(const edm::ParameterSet&);
0041   ~LHEPtFilter() override;
0042
0043   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0044   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0045
0046 private:
0047   // ----------member data ---------------------------
0048
0049   edm::EDGetTokenT<LHEEventProduct> src_;
0050   std::vector<int> pdgIdVec_;
0051   std::set<int> pdgIds_;  // Set of PDG Ids to include
0052   double ptMin_;          // number of particles required to pass filter
0053   double ptMax_;          // number of particles required to pass filter
0054   bool isScalar_;         // true for scalar sum or false for vector sum
0055 };
0056
0057 using namespace edm;
0058 using namespace std;
0059
0060 LHEPtFilter::LHEPtFilter(const edm::ParameterSet& iConfig)
0061     : pdgIdVec_(iConfig.getParameter<std::vector<int>>("selectedPdgIds")),
0062       ptMin_(iConfig.getParameter<double>("ptMin")),
0063       ptMax_(iConfig.getParameter<double>("ptMax")),
0064       isScalar_(iConfig.getParameter<bool>("isScalar")) {
0065   //here do whatever other initialization is needed
0066   src_ = consumes<LHEEventProduct>(iConfig.getParameter<edm::InputTag>("src"));
0067   pdgIds_ = std::set<int>(pdgIdVec_.begin(), pdgIdVec_.end());
0068 }
0069
0070 LHEPtFilter::~LHEPtFilter() {
0071   // do anything here that needs to be done at destruction time
0072   // (e.g. close files, deallocate resources etc.)
0073 }
0074
0075 // ------------ method called to skim the data  ------------
0076 bool LHEPtFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0077   edm::Handle<LHEEventProduct> EvtHandle;
0078   iEvent.getByToken(src_, EvtHandle);
0079
0080   std::vector<lhef::HEPEUP::FiveVector> lheParticles = EvtHandle->hepeup().PUP;
0081   std::vector<ROOT::Math::PxPyPzEVector> cands;
0082
0083   for (unsigned int i = 0; i < lheParticles.size(); ++i) {
0084     if (EvtHandle->hepeup().ISTUP[i] != 1) {  // keep only outgoing particles
0085       continue;
0086     }
0087     int pdgId = EvtHandle->hepeup().IDUP[i];
0088     if (pdgIds_.count(pdgId)) {
0089       cands.push_back(
0090           ROOT::Math::PxPyPzEVector(lheParticles[i][0], lheParticles[i][1], lheParticles[i][2], lheParticles[i][3]));
0091     }
0092   }
0093   double vpt_ = -1;
0094   if (isScalar_) {  // do the scalar sum
0095     if (!cands.empty()) {
0096       double tot = cands.at(0).pt();
0097       for (unsigned icand = 1; icand < cands.size(); ++icand) {
0098         tot += cands.at(icand).pt();
0099       }
0100       vpt_ = tot;
0101     }
0102     if ((ptMax_ < 0. || vpt_ <= ptMax_) && vpt_ > ptMin_) {
0103       return true;
0104     } else {
0105       return false;
0106     }
0107   } else {  // else do the vector sum
0108     if (!cands.empty()) {
0109       ROOT::Math::PxPyPzEVector tot = cands.at(0);
0110       for (unsigned icand = 1; icand < cands.size(); ++icand) {
0111         tot += cands.at(icand);
0112       }
0113       vpt_ = tot.pt();
0114     }
0115     if ((ptMax_ < 0. || vpt_ <= ptMax_) && vpt_ > ptMin_) {
0116       return true;
0117     } else {
0118       return false;
0119     }
0120   }
0121 }
0122
0123 void LHEPtFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0124   edm::ParameterSetDescription desc;
0125   desc.add<bool>("isScalar", false);  // default is false
0126 }
0127
0128 //define this as a plug-in
0129 DEFINE_FWK_MODULE(LHEPtFilter);