GenFilters/plugins/PythiaFilterMultiMother.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/global/EDFilter.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/EDGetToken.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0010 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0011
0012 #include <cmath>
0013 #include <cstdlib>
0014 #include <string>
0015 #include <vector>
0016
0017 class PythiaFilterMultiMother : public edm::global::EDFilter<> {
0018 public:
0019   explicit PythiaFilterMultiMother(const edm::ParameterSet&);
0020
0021   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0022
0023 private:
0024   // ----------member data ---------------------------
0025
0026   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0027   const int particleID;
0028   const double minpcut;
0029   const double maxpcut;
0030   const double minptcut;
0031   const double maxptcut;
0032   const double minetacut;
0033   const double maxetacut;
0034   const double minrapcut;
0035   const double maxrapcut;
0036   const double minphicut;
0037   const double maxphicut;
0038
0039   const int status;
0040   const std::vector<int> motherIDs;
0041   const int processID;
0042
0043   const double betaBoost;
0044 };
0045
0046 using namespace std;
0047
0048 PythiaFilterMultiMother::PythiaFilterMultiMother(const edm::ParameterSet& iConfig)
0049     : token_(consumes<edm::HepMCProduct>(
0050           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0051       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0052       minpcut(iConfig.getUntrackedParameter("MinP", 0.)),
0053       maxpcut(iConfig.getUntrackedParameter("MaxP", 10000.)),
0054       minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0055       maxptcut(iConfig.getUntrackedParameter("MaxPt", 10000.)),
0056       minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0057       maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
0058       minrapcut(iConfig.getUntrackedParameter("MinRapidity", -20.)),
0059       maxrapcut(iConfig.getUntrackedParameter("MaxRapidity", 20.)),
0060       minphicut(iConfig.getUntrackedParameter("MinPhi", -3.5)),
0061       maxphicut(iConfig.getUntrackedParameter("MaxPhi", 3.5)),
0062       status(iConfig.getUntrackedParameter("Status", 0)),
0063       motherIDs(iConfig.getUntrackedParameter("MotherIDs", std::vector<int>{0})),
0064       processID(iConfig.getUntrackedParameter("ProcessID", 0)),
0065       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {}
0066
0067 bool PythiaFilterMultiMother::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0068   bool accepted = false;
0069   edm::Handle<edm::HepMCProduct> evt;
0070   iEvent.getByToken(token_, evt);
0071
0072   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0073
0074   if (processID == 0 || processID == myGenEvent->signal_process_id()) {
0075     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0076          ++p) {
0077       HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0078       double rapidity = 0.5 * log((mom.e() + mom.pz()) / (mom.e() - mom.pz()));
0079
0080       if (abs((*p)->pdg_id()) == particleID && mom.rho() > minpcut && mom.rho() < maxpcut &&
0081           (*p)->momentum().perp() > minptcut && (*p)->momentum().perp() < maxptcut && mom.eta() > minetacut &&
0082           mom.eta() < maxetacut && rapidity > minrapcut && rapidity < maxrapcut && (*p)->momentum().phi() > minphicut &&
0083           (*p)->momentum().phi() < maxphicut) {
0084         for (std::vector<int>::const_iterator motherID = motherIDs.begin(); motherID != motherIDs.end(); ++motherID) {
0085           if (status == 0 && *motherID == 0) {
0086             accepted = true;
0087           }
0088           if (status != 0 && *motherID == 0) {
0089             if ((*p)->status() == status)
0090               accepted = true;
0091           }
0092
0093           HepMC::GenParticle* mother = (*((*p)->production_vertex()->particles_in_const_begin()));
0094
0095           if (status == 0 && *motherID != 0) {
0096             if (abs(mother->pdg_id()) == abs(*motherID)) {
0097               accepted = true;
0098             }
0099           }
0100           if (status != 0 && *motherID != 0) {
0101             if ((*p)->status() == status && abs(mother->pdg_id()) == abs(*motherID)) {
0102               accepted = true;
0103             }
0104           }
0105         }
0106       }
0107     }
0108   } else {
0109     accepted = true;
0110   }
0111   return accepted;
0112 }
0113
0114 DEFINE_FWK_MODULE(PythiaFilterMultiMother);