GenFilters/plugins/BsJpsiPhiFilter.cc

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/global/EDFilter.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Utilities/interface/EDGetToken.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 #include "HepMC/GenVertex.h"
0009 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0010
0011 #include <cmath>
0012 #include <cstdlib>
0013 #include <string>
0014 #include <vector>
0015
0016 class BsJpsiPhiFilter : public edm::global::EDFilter<> {
0017 public:
0018   explicit BsJpsiPhiFilter(const edm::ParameterSet&);
0019
0020   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0021
0022 private:
0023   struct CutStruct {
0024     int type;
0025     double etaMin, etaMax, ptMin;
0026   };
0027
0028   HepMC::GenParticle* findParticle(HepMC::GenVertex*, const int requested_id) const;
0029
0030   HepMC::GenEvent::particle_const_iterator getNextBs(const HepMC::GenEvent::particle_const_iterator start,
0031                                                      const HepMC::GenEvent::particle_const_iterator end) const;
0032
0033   bool cuts(const HepMC::GenParticle* jpsi, const CutStruct& cut) const;
0034   bool etaInRange(float eta, float etamin, float etamax) const;
0035
0036   CutStruct leptonCuts, hadronCuts;
0037
0038   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0039 };
0040
0041 using namespace std;
0042
0043 BsJpsiPhiFilter::BsJpsiPhiFilter(const edm::ParameterSet& iConfig)
0044     : token_(consumes<edm::HepMCProduct>(
0045           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))) {
0046   hadronCuts.type = iConfig.getParameter<int>("hadronType");
0047   hadronCuts.etaMin = iConfig.getParameter<double>("hadronEtaMin");
0048   hadronCuts.etaMax = iConfig.getParameter<double>("hadronEtaMax");
0049   hadronCuts.ptMin = iConfig.getParameter<double>("hadronPtMin");
0050   leptonCuts.type = iConfig.getParameter<int>("leptonType");
0051   leptonCuts.etaMin = iConfig.getParameter<double>("leptonEtaMin");
0052   leptonCuts.etaMax = iConfig.getParameter<double>("leptonEtaMax");
0053   leptonCuts.ptMin = iConfig.getParameter<double>("leptonPtMin");
0054 }
0055
0056 HepMC::GenParticle* BsJpsiPhiFilter::findParticle(HepMC::GenVertex* vertex, const int requested_id) const {
0057   for (std::vector<HepMC::GenParticle*>::const_iterator p = vertex->particles_out_const_begin();
0058        p != vertex->particles_out_const_end();
0059        p++) {
0060     int event_particle_id = abs((*p)->pdg_id());
0061     if (requested_id == event_particle_id)
0062       return *p;
0063   }
0064   return nullptr;
0065 }
0066
0067 HepMC::GenEvent::particle_const_iterator BsJpsiPhiFilter::getNextBs(
0068     const HepMC::GenEvent::particle_const_iterator start, const HepMC::GenEvent::particle_const_iterator end) const {
0069   HepMC::GenEvent::particle_const_iterator p;
0070   for (p = start; p != end; p++) {
0071     int event_particle_id = abs((*p)->pdg_id());
0072     if (event_particle_id == 531)
0073       return p;
0074   }
0075   return p;
0076 }
0077
0078 bool BsJpsiPhiFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0079   edm::Handle<edm::HepMCProduct> evt;
0080   iEvent.getByToken(token_, evt);
0081
0082   const HepMC::GenEvent* generated_event = evt->GetEvent();
0083
0084   bool event_passed = false;
0085   HepMC::GenEvent::particle_const_iterator bs =
0086       getNextBs(generated_event->particles_begin(), generated_event->particles_end());
0087   while (bs != generated_event->particles_end()) {
0088     HepMC::GenVertex* outVertex = (*bs)->end_vertex();
0089
0090     HepMC::GenParticle* jpsi = nullptr;
0091     HepMC::GenParticle* phi = nullptr;
0092     int numChildren = outVertex->particles_out_size();
0093
0094     if ((numChildren == 2) && ((jpsi = findParticle(outVertex, 443)) != nullptr) &&
0095         ((phi = findParticle(outVertex, 333)) != nullptr)) {
0096       if (cuts(phi, hadronCuts) && cuts(jpsi, leptonCuts)) {
0097         event_passed = true;
0098         break;
0099       }
0100     }
0101     bs = getNextBs(++bs, generated_event->particles_end());
0102   }
0103
0104   delete generated_event;
0105
0106   return event_passed;
0107 }
0108
0109 bool BsJpsiPhiFilter::cuts(const HepMC::GenParticle* jpsi, const CutStruct& cut) const {
0110   HepMC::GenVertex* myVertex = jpsi->end_vertex();
0111   int numChildren = myVertex->particles_out_size();
0112   std::vector<HepMC::GenParticle*> psiChild;
0113   for (std::vector<HepMC::GenParticle*>::const_iterator p = myVertex->particles_out_const_begin();
0114        p != myVertex->particles_out_const_end();
0115        p++)
0116     psiChild.push_back((*p));
0117
0118   if (numChildren > 1) {
0119     if (psiChild.size() == 2 && (abs(psiChild[0]->pdg_id()) == cut.type) && (abs(psiChild[1]->pdg_id()) == cut.type)) {
0120       return ((etaInRange(psiChild[0]->momentum().eta(), cut.etaMin, cut.etaMax)) &&
0121               (etaInRange(psiChild[1]->momentum().eta(), cut.etaMin, cut.etaMax)) &&
0122               (psiChild[0]->momentum().perp() > cut.ptMin) && (psiChild[1]->momentum().perp() > cut.ptMin));
0123     }
0124     return false;
0125   }
0126   return false;
0127 }
0128
0129 bool BsJpsiPhiFilter::etaInRange(float eta, float etamin, float etamax) const {
0130   return ((etamin < eta) && (eta < etamax));
0131 }
0132
0133 DEFINE_FWK_MODULE(BsJpsiPhiFilter);