GenFilters/plugins/PythiaFilterMultiAncestor.cc

0001 /*
0002 This filter allows to select events where a given particle originates *either directly
0003 (daughter) or indirectly ((grand)^n-daughter)*  from a list of possible ancestors.
0004
0005 It also allows to filter on the immediate daughters of the said particle.
0006
0007 Kinematic selections can also be applied on the particle and its daughters.
0008
0009 The example below shows how to select a jpsi->mumu coming from *any* b hadron (notice
0010 that MotherIDs is just 5, i.e. b-quark) however long the decay chain is.
0011 This includes, for example, B->Jpsi + X as well as B->Psi(2S)(->JPsi) X
0012
0013 process.jpsi_from_bhadron_filter = cms.EDFilter("PythiaFilterMultiAncestor",
0014     DaughterIDs = cms.untracked.vint32(-13, 13),
0015     DaughterMaxEtas = cms.untracked.vdouble(3., 3.),
0016     DaughterMaxPts = cms.untracked.vdouble(100000.0, 100000.0),
0017     DaughterMinEtas = cms.untracked.vdouble(-2.6, -2.6),
0018     DaughterMinPts = cms.untracked.vdouble(2.5, 2.5),
0019     MaxEta = cms.untracked.double(3.0),
0020     MinEta = cms.untracked.double(-3.0),
0021     MinPt = cms.untracked.double(6.0),
0022     MotherIDs = cms.untracked.vint32(5),
0023     ParticleID = cms.untracked.int32(443)
0024 )
0025 */
0026
0027 // system include files
0028 #include <memory>
0029 #include <iostream>
0030
0031 // user include files
0032 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0033
0034 #include "FWCore/Framework/interface/global/EDFilter.h"
0035 #include "FWCore/Framework/interface/Event.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039
0040 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0041
0042 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0043
0044 using namespace edm;
0045 using namespace std;
0046
0047 namespace edm {
0048   class HepMCProduct;
0049 }
0050
0051 class PythiaFilterMultiAncestor : public edm::global::EDFilter<> {
0052 public:
0053   explicit PythiaFilterMultiAncestor(const edm::ParameterSet&);
0054
0055   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0056
0057 private:
0058   bool isAncestor(HepMC::GenParticle* particle, int IDtoMatch) const;
0059
0060   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0061   const int particleID;
0062   const double minpcut;
0063   const double maxpcut;
0064   const double minptcut;
0065   const double maxptcut;
0066   const double minetacut;
0067   const double maxetacut;
0068   const double minrapcut;
0069   const double maxrapcut;
0070   const double minphicut;
0071   const double maxphicut;
0072
0073   const int status;
0074   const std::vector<int> motherIDs;
0075   const std::vector<int> daughterIDs;
0076   const std::vector<double> daughterMinPts;
0077   const std::vector<double> daughterMaxPts;
0078   const std::vector<double> daughterMinEtas;
0079   const std::vector<double> daughterMaxEtas;
0080
0081   const int processID;
0082
0083   const double betaBoost;
0084 };
0085
0086 PythiaFilterMultiAncestor::PythiaFilterMultiAncestor(const edm::ParameterSet& iConfig)
0087     : token_(consumes<edm::HepMCProduct>(
0088           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0089       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0090       minpcut(iConfig.getUntrackedParameter("MinP", 0.)),
0091       maxpcut(iConfig.getUntrackedParameter("MaxP", 10000.)),
0092       minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0093       maxptcut(iConfig.getUntrackedParameter("MaxPt", 10000.)),
0094       minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0095       maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
0096       minrapcut(iConfig.getUntrackedParameter("MinRapidity", -20.)),
0097       maxrapcut(iConfig.getUntrackedParameter("MaxRapidity", 20.)),
0098       minphicut(iConfig.getUntrackedParameter("MinPhi", -3.5)),
0099       maxphicut(iConfig.getUntrackedParameter("MaxPhi", 3.5)),
0100       status(iConfig.getUntrackedParameter("Status", 0)),
0101       motherIDs(iConfig.getUntrackedParameter("MotherIDs", std::vector<int>{0})),
0102       daughterIDs(iConfig.getUntrackedParameter("DaughterIDs", std::vector<int>{0})),
0103       daughterMinPts(iConfig.getUntrackedParameter("DaughterMinPts", std::vector<double>{0.})),
0104       daughterMaxPts(iConfig.getUntrackedParameter("DaughterMaxPts", std::vector<double>{10000.})),
0105       daughterMinEtas(iConfig.getUntrackedParameter("DaughterMinEtas", std::vector<double>{-10.})),
0106       daughterMaxEtas(iConfig.getUntrackedParameter("DaughterMaxEtas", std::vector<double>{10.})),
0107       processID(iConfig.getUntrackedParameter("ProcessID", 0)),
0108       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
0109   //now do what ever initialization is needed
0110 }
0111
0112 // ------------ access the full genealogy ---------
0113 bool PythiaFilterMultiAncestor::isAncestor(HepMC::GenParticle* particle, int IDtoMatch) const {
0114   for (HepMC::GenVertex::particle_iterator ancestor = particle->production_vertex()->particles_begin(HepMC::ancestors);
0115        ancestor != particle->production_vertex()->particles_end(HepMC::ancestors);
0116        ++ancestor) {
0117     // std::cout << __LINE__ << "]\t particle's PDG ID " << particle->pdg_id()
0118     //                       << " \t particle's ancestor's PDG ID " << (*ancestor)->pdg_id()
0119     //                       << " \t ID to match " << IDtoMatch << std::endl;
0120
0121     if (abs((*ancestor)->pdg_id()) == abs(IDtoMatch)) {
0122       //  std::cout << __LINE__ << "]\t found!" << std::endl;
0123       return true;
0124     }
0125   }
0126
0127   // std::cout << __LINE__ << "]\t nope, no luck" << std::endl;
0128   return false;
0129 }
0130
0131 // ------------ method called to produce the data  ------------
0132 bool PythiaFilterMultiAncestor::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0133   using namespace edm;
0134   bool accepted = false;
0135   Handle<HepMCProduct> evt;
0136   iEvent.getByToken(token_, evt);
0137
0138   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0139
0140   if (processID == 0 || processID == myGenEvent->signal_process_id()) {
0141     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0142          ++p) {
0143       HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0144       double rapidity = 0.5 * log((mom.e() + mom.pz()) / (mom.e() - mom.pz()));
0145
0146       if (abs((*p)->pdg_id()) == particleID && mom.rho() > minpcut && mom.rho() < maxpcut &&
0147           (*p)->momentum().perp() > minptcut && (*p)->momentum().perp() < maxptcut && mom.eta() > minetacut &&
0148           mom.eta() < maxetacut && rapidity > minrapcut && rapidity < maxrapcut && (*p)->momentum().phi() > minphicut &&
0149           (*p)->momentum().phi() < maxphicut) {
0150         // find the mother
0151         for (std::vector<int>::const_iterator motherID = motherIDs.begin(); motherID != motherIDs.end(); ++motherID) {
0152           // check status if no mother's pdgID is specified
0153           if (status == 0 && *motherID == 0) {
0154             accepted = true;
0155           }
0156           if (status != 0 && *motherID == 0) {
0157             if ((*p)->status() == status)
0158               accepted = true;
0159           }
0160
0161           // check the mother's pdgID
0162           if (status == 0 && *motherID != 0) {
0163             // if (abs(mother->pdg_id()) == abs(*motherID)) {
0164             if (isAncestor(*p, *motherID)) {
0165               accepted = true;
0166             }
0167           }
0168           if (status != 0 && *motherID != 0) {
0169             // if ((*p)->status() == status && abs(mother->pdg_id()) == abs(*motherID)){
0170             if ((*p)->status() == status && isAncestor(*p, *motherID)) {
0171               accepted = true;
0172             }
0173           }
0174         }
0175
0176         // find the daughters
0177         if (accepted & (!daughterIDs.empty())) {
0178           // if you got this far it means that the mother was found
0179           // now let's check the daughters
0180           // use a counter, if there's enough daughters that match the pdg and kinematic
0181           // criteria accept the event
0182           uint good_dau = 0;
0183           uint good_dau_cc = 0;
0184           for (HepMC::GenVertex::particle_iterator dau = (*p)->end_vertex()->particles_begin(HepMC::children);
0185                dau != (*p)->end_vertex()->particles_end(HepMC::children);
0186                ++dau) {
0187             for (unsigned int i = 0; i < daughterIDs.size(); ++i) {
0188               // if a daughter has its pdgID among the desired ones, apply kin cuts on it
0189               // if it survives, add a notch to the counter
0190               if ((*dau)->pdg_id() == daughterIDs[i]) {
0191                 if ((*dau)->momentum().perp() < daughterMinPts[i])
0192                   continue;
0193                 if ((*dau)->momentum().perp() > daughterMaxPts[i])
0194                   continue;
0195                 if ((*dau)->momentum().eta() < daughterMinEtas[i])
0196                   continue;
0197                 if ((*dau)->momentum().eta() > daughterMaxEtas[i])
0198                   continue;
0199                 ++good_dau;
0200               }
0201               // check charge conjugation
0202               if (-(*dau)->pdg_id() == daughterIDs[i]) {  // notice minus sign
0203                 if ((*dau)->momentum().perp() < daughterMinPts[i])
0204                   continue;
0205                 if ((*dau)->momentum().perp() > daughterMaxPts[i])
0206                   continue;
0207                 if ((*dau)->momentum().eta() < daughterMinEtas[i])
0208                   continue;
0209                 if ((*dau)->momentum().eta() > daughterMaxEtas[i])
0210                   continue;
0211                 ++good_dau_cc;
0212               }
0213             }
0214           }
0215           if (good_dau < daughterIDs.size() && good_dau_cc < daughterIDs.size())
0216             accepted = false;
0217         }
0218       }
0219       // only need to satisfy the conditions _once_
0220       if (accepted)
0221         break;
0222     }
0223
0224   } else {
0225     accepted = true;
0226   }
0227
0228   if (accepted) {
0229     return true;
0230   } else {
0231     return false;
0232   }
0233 }
0234
0235 DEFINE_FWK_MODULE(PythiaFilterMultiAncestor);