Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:42 

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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!