GenFilters/plugins/LHEIntermediateParticlePtFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    LHEIntermediateParticlePtFilter
0004 // Class:      LHEIntermediateParticlePtFilter
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 <set>
0018
0019 // user include files
0020 #include "Math/Vector4D.h"
0021 #include "Math/Vector4Dfwd.h"
0022
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/global/EDFilter.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0031
0032 //
0033 // class declaration
0034 //
0035
0036 class LHEIntermediateParticlePtFilter : public edm::global::EDFilter<> {
0037 public:
0038   explicit LHEIntermediateParticlePtFilter(const edm::ParameterSet&);
0039   ~LHEIntermediateParticlePtFilter() override = default;
0040
0041   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0042
0043 private:
0044   // ----------member data ---------------------------
0045
0046   edm::EDGetTokenT<LHEEventProduct> src_;
0047   const std::vector<int> pdgIdVec_;
0048   std::set<int> pdgIds_;  // Set of PDG Ids to include
0049   const double ptMin_;    // number of particles required to pass filter
0050   const double ptMax_;    // number of particles required to pass filter
0051 };
0052
0053 using namespace edm;
0054 using namespace std;
0055
0056 LHEIntermediateParticlePtFilter::LHEIntermediateParticlePtFilter(const edm::ParameterSet& iConfig)
0057     : pdgIdVec_(iConfig.getParameter<std::vector<int>>("selectedPdgIds")),
0058       ptMin_(iConfig.getParameter<double>("ptMin")),
0059       ptMax_(iConfig.getParameter<double>("ptMax")) {
0060   //here do whatever other initialization is needed
0061   src_ = consumes<LHEEventProduct>(iConfig.getParameter<edm::InputTag>("src"));
0062   pdgIds_ = std::set<int>(pdgIdVec_.begin(), pdgIdVec_.end());
0063 }
0064
0065 // ------------ method called to skim the data  ------------
0066 bool LHEIntermediateParticlePtFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0067   edm::Handle<LHEEventProduct> EvtHandle;
0068   iEvent.getByToken(src_, EvtHandle);
0069
0070   std::vector<lhef::HEPEUP::FiveVector> lheParticles = EvtHandle->hepeup().PUP;
0071   std::vector<ROOT::Math::PxPyPzEVector> cands;
0072
0073   for (unsigned int i = 0; i < lheParticles.size(); ++i) {
0074     if (EvtHandle->hepeup().ISTUP[i] != 2) {  // look at intermediate resonances
0075       continue;
0076     }
0077     int pdgId = EvtHandle->hepeup().IDUP[i];
0078     if (pdgIds_.count(pdgId)) {
0079       cands.push_back(
0080           ROOT::Math::PxPyPzEVector(lheParticles[i][0], lheParticles[i][1], lheParticles[i][2], lheParticles[i][3]));
0081     }
0082   }
0083   double vpt_ = -1;
0084   if (!cands.empty()) {
0085     ROOT::Math::PxPyPzEVector tot = cands.at(0);
0086     for (unsigned icand = 1; icand < cands.size(); ++icand) {
0087       tot += cands.at(icand);
0088     }
0089     vpt_ = tot.pt();
0090   }
0091
0092   return (ptMax_ < 0. || vpt_ <= ptMax_) && (vpt_ > ptMin_);
0093 }
0094
0095 //define this as a plug-in
0096 DEFINE_FWK_MODULE(LHEIntermediateParticlePtFilter);