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.
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/ParameterSet.h"
0031 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0032
0033 //
0034 // class declaration
0035 //
0036
0037 class LHEPtFilter : public edm::global::EDFilter<> {
0038 public:
0039   explicit LHEPtFilter(const edm::ParameterSet&);
0040   ~LHEPtFilter() override;
0041
0042   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0043
0044 private:
0045   // ----------member data ---------------------------
0046
0047   edm::EDGetTokenT<LHEEventProduct> src_;
0048   std::vector<int> pdgIdVec_;
0049   std::set<int> pdgIds_;  // Set of PDG Ids to include
0050   double ptMin_;          // number of particles required to pass filter
0051   double ptMax_;          // number of particles required to pass filter
0052 };
0053
0054 using namespace edm;
0055 using namespace std;
0056
0057 LHEPtFilter::LHEPtFilter(const edm::ParameterSet& iConfig)
0058     : pdgIdVec_(iConfig.getParameter<std::vector<int>>("selectedPdgIds")),
0059       ptMin_(iConfig.getParameter<double>("ptMin")),
0060       ptMax_(iConfig.getParameter<double>("ptMax")) {
0061   //here do whatever other initialization is needed
0062   src_ = consumes<LHEEventProduct>(iConfig.getParameter<edm::InputTag>("src"));
0063   pdgIds_ = std::set<int>(pdgIdVec_.begin(), pdgIdVec_.end());
0064 }
0065
0066 LHEPtFilter::~LHEPtFilter() {
0067   // do anything here that needs to be done at destruction time
0068   // (e.g. close files, deallocate resources etc.)
0069 }
0070
0071 // ------------ method called to skim the data  ------------
0072 bool LHEPtFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0073   edm::Handle<LHEEventProduct> EvtHandle;
0074   iEvent.getByToken(src_, EvtHandle);
0075
0076   std::vector<lhef::HEPEUP::FiveVector> lheParticles = EvtHandle->hepeup().PUP;
0077   std::vector<ROOT::Math::PxPyPzEVector> cands;
0078
0079   for (unsigned int i = 0; i < lheParticles.size(); ++i) {
0080     if (EvtHandle->hepeup().ISTUP[i] != 1) {  // keep only outgoing particles
0081       continue;
0082     }
0083     int pdgId = EvtHandle->hepeup().IDUP[i];
0084     if (pdgIds_.count(pdgId)) {
0085       cands.push_back(
0086           ROOT::Math::PxPyPzEVector(lheParticles[i][0], lheParticles[i][1], lheParticles[i][2], lheParticles[i][3]));
0087     }
0088   }
0089   double vpt_ = -1;
0090   if (!cands.empty()) {
0091     ROOT::Math::PxPyPzEVector tot = cands.at(0);
0092     for (unsigned icand = 1; icand < cands.size(); ++icand) {
0093       tot += cands.at(icand);
0094     }
0095     vpt_ = tot.pt();
0096   }
0097   if ((ptMax_ < 0. || vpt_ <= ptMax_) && vpt_ > ptMin_) {
0098     return true;
0099   } else {
0100     return false;
0101   }
0102 }
0103
0104 //define this as a plug-in
0105 DEFINE_FWK_MODULE(LHEPtFilter);