GenFilters/plugins/MCSingleParticleYPt.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MCSingleParticleYPt
0004 // Class:      MCSingleParticleYPt
0005 //
0006 /*
0007
0008  Description: filter events based on the Pythia particleID, Pt, Y and status. It is based on MCSingleParticleFilter.
0009               It will used to filter a b-hadron with the given kinematics, only one b-hadron is required to match.
0010  Implementation: inherits from generic EDFilter
0011
0012 */
0013 //
0014 // Author: Alberto Sanchez-Hernandez
0015 // Adapted on: August 2016
0016 //
0017 //
0018 // Filter based on MCSingleParticleFilter.cc, but using rapidity instead of eta
0019
0020 #include "DataFormats/Common/interface/Handle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/global/EDFilter.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/Utilities/interface/EDGetToken.h"
0028 #include "FWCore/Utilities/interface/InputTag.h"
0029 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0030
0031 #include <cmath>
0032 #include <ostream>
0033 #include <string>
0034 #include <vector>
0035
0036 class MCSingleParticleYPt : public edm::global::EDFilter<> {
0037 public:
0038   explicit MCSingleParticleYPt(const edm::ParameterSet&);
0039
0040   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0041
0042 private:
0043   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0044   const int fVerbose;
0045   const bool fchekantiparticle;
0046   std::vector<int> particleID;
0047   std::vector<double> ptMin;
0048   std::vector<double> rapMin;
0049   std::vector<double> rapMax;
0050   std::vector<int> status;
0051 };
0052
0053 using namespace edm;
0054 using namespace std;
0055
0056 MCSingleParticleYPt::MCSingleParticleYPt(const edm::ParameterSet& iConfig)
0057     : token_(consumes<edm::HepMCProduct>(
0058           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0059       fVerbose(iConfig.getUntrackedParameter("verbose", 0)),
0060       fchekantiparticle(iConfig.getUntrackedParameter("CheckAntiparticle", true)) {
0061   vector<int> defpid;
0062   defpid.push_back(0);
0063   particleID = iConfig.getUntrackedParameter<vector<int> >("ParticleID", defpid);
0064   vector<double> defptmin;
0065   defptmin.push_back(0.);
0066   ptMin = iConfig.getUntrackedParameter<vector<double> >("MinPt", defptmin);
0067   vector<double> defrapmin;
0068   defrapmin.push_back(-10.);
0069   rapMin = iConfig.getUntrackedParameter<vector<double> >("MinY", defrapmin);
0070   vector<double> defrapmax;
0071   defrapmax.push_back(10.);
0072   rapMax = iConfig.getUntrackedParameter<vector<double> >("MaxY", defrapmax);
0073   vector<int> defstat;
0074   defstat.push_back(0);
0075   status = iConfig.getUntrackedParameter<vector<int> >("Status", defstat);
0076
0077   // check for same size
0078   if ((ptMin.size() > 1 && particleID.size() != ptMin.size()) ||
0079       (rapMin.size() > 1 && particleID.size() != rapMin.size()) ||
0080       (rapMax.size() > 1 && particleID.size() != rapMax.size()) ||
0081       (status.size() > 1 && particleID.size() != status.size())) {
0082     edm::LogWarning("MCSingleParticleYPt")
0083         << "WARNING: MCSingleParticleYPt : size of vector cuts do not match!!" << endl;
0084   }
0085
0086   // if ptMin size smaller than particleID , fill up further with defaults
0087   if (particleID.size() > ptMin.size()) {
0088     vector<double> defptmin2;
0089     for (unsigned int i = 0; i < particleID.size(); i++) {
0090       defptmin2.push_back(0.);
0091     }
0092     ptMin = defptmin2;
0093   }
0094   // if etaMin size smaller than particleID , fill up further with defaults
0095   if (particleID.size() > rapMin.size()) {
0096     vector<double> defrapmin2;
0097     for (unsigned int i = 0; i < particleID.size(); i++) {
0098       defrapmin2.push_back(-10.);
0099     }
0100     rapMin = defrapmin2;
0101   }
0102   // if etaMax size smaller than particleID , fill up further with defaults
0103   if (particleID.size() > rapMax.size()) {
0104     vector<double> defrapmax2;
0105     for (unsigned int i = 0; i < particleID.size(); i++) {
0106       defrapmax2.push_back(10.);
0107     }
0108     rapMax = defrapmax2;
0109   }
0110   // if status size smaller than particleID , fill up further with defaults
0111   if (particleID.size() > status.size()) {
0112     vector<int> defstat2;
0113     for (unsigned int i = 0; i < particleID.size(); i++) {
0114       defstat2.push_back(0);
0115     }
0116     status = defstat2;
0117   }
0118
0119   if (fVerbose > 0) {
0120     edm::LogInfo("MCSingleParticleYPt") << "----------------------------------------------------------------------"
0121                                         << std::endl;
0122     edm::LogInfo("MCSingleParticleYPt") << "----- MCSingleParticleYPt" << std::endl;
0123     for (unsigned int i = 0; i < particleID.size(); ++i) {
0124       edm::LogInfo("MCSingleParticleYPt") << " ID: " << particleID[i] << " pT > " << ptMin[i] << ",   " << rapMin[i]
0125                                           << " < y < " << rapMax[i] << ",   status = " << status[i] << std::endl;
0126     }
0127     if (fchekantiparticle)
0128       edm::LogInfo("MCSingleParticleYPt") << " anti-particles will be tested as well." << std::endl;
0129     edm::LogInfo("MCSingleParticleYPt") << "----------------------------------------------------------------------"
0130                                         << std::endl;
0131   }
0132 }
0133
0134 bool MCSingleParticleYPt::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0135   bool accepted = false;
0136   Handle<HepMCProduct> evt;
0137   iEvent.getByToken(token_, evt);
0138
0139   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0140   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0141        ++p) {
0142     if (fVerbose > 3)
0143       edm::LogInfo("MCSingleParticleYPt")
0144           << "Looking at particle : " << (*p)->pdg_id() << " status : " << (*p)->status() << std::endl;
0145
0146     for (unsigned int i = 0; i < particleID.size(); i++) {
0147       if (particleID[i] == (*p)->pdg_id() || (fchekantiparticle && (-particleID[i] == (*p)->pdg_id())) ||
0148           particleID[i] == 0) {
0149         // calculate rapidity just for the desired particle and make sure, this particles has enough energy
0150         double rapidity = ((*p)->momentum().e() - (*p)->momentum().pz()) > 0.
0151                               ? 0.5 * log(((*p)->momentum().e() + (*p)->momentum().pz()) /
0152                                           ((*p)->momentum().e() - (*p)->momentum().pz()))
0153                               : rapMax[i] + .1;
0154         if (fVerbose > 2)
0155           edm::LogInfo("MCSingleParticleYPt")
0156               << "Testing particle : " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp() << " y: " << rapidity
0157               << " status : " << (*p)->status() << endl;
0158         if ((*p)->momentum().perp() > ptMin[i] && rapidity > rapMin[i] && rapidity < rapMax[i] &&
0159             ((*p)->status() == status[i] || status[i] == 0)) {
0160           accepted = true;
0161           if (fVerbose > 1)
0162             edm::LogInfo("MCSingleParticleYPt")
0163                 << "Accepted particle : " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp() << " y: " << rapidity
0164                 << " status : " << (*p)->status() << endl;
0165           break;
0166         }
0167       }
0168     }
0169     if (accepted)
0170       break;
0171   }
0172   return accepted;
0173 }
0174
0175 DEFINE_FWK_MODULE(MCSingleParticleYPt);