GenFilters/plugins/PythiaFilterMotherSister.cc

0001
0002 #include "GeneratorInterface/GenFilters/plugins/PythiaFilterMotherSister.h"
0003 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0004
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006 #include <iostream>
0007
0008 using namespace edm;
0009 using namespace std;
0010
0011 PythiaFilterMotherSister::PythiaFilterMotherSister(const edm::ParameterSet& iConfig)
0012     : token_(consumes<edm::HepMCProduct>(
0013           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0014       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0015       minpcut(iConfig.getUntrackedParameter("MinP", 0.)),
0016       maxpcut(iConfig.getUntrackedParameter("MaxP", 10000.)),
0017       minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0018       maxptcut(iConfig.getUntrackedParameter("MaxPt", 10000.)),
0019       minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0020       maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
0021       minrapcut(iConfig.getUntrackedParameter("MinRapidity", -20.)),
0022       maxrapcut(iConfig.getUntrackedParameter("MaxRapidity", 20.)),
0023       minphicut(iConfig.getUntrackedParameter("MinPhi", -3.5)),
0024       maxphicut(iConfig.getUntrackedParameter("MaxPhi", 3.5)),
0025       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)),
0026       motherIDs(iConfig.getUntrackedParameter("MotherIDs", std::vector<int>{0})),
0027       sisterID(iConfig.getUntrackedParameter("SisterID", 0)),
0028       maxSisDisplacement(iConfig.getUntrackedParameter("MaxSisterDisplacement", -1.)),
0029       nephewIDs(iConfig.getUntrackedParameter("NephewIDs", std::vector<int>{0})),
0030       minNephewPts(iConfig.getUntrackedParameter("MinNephewPts", std::vector<double>{0.})) {
0031   if (nephewIDs.size() != minNephewPts.size()) {
0032     throw cms::Exception("BadConfig") << "PythiaFilterMotherSister: "
0033                                       << "'nephewIDs' and 'minNephewPts' need same length.";
0034   }
0035 }
0036
0037 PythiaFilterMotherSister::~PythiaFilterMotherSister() {
0038   // do anything here that needs to be done at desctruction time
0039   // (e.g. close files, deallocate resources etc.)
0040 }
0041
0042 //
0043 // member functions
0044 //
0045
0046 // ------------ method called to produce the data  ------------
0047 bool PythiaFilterMotherSister::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0048   using namespace edm;
0049   Handle<HepMCProduct> evt;
0050   iEvent.getByToken(token_, evt);
0051
0052   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0053
0054   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0055        ++p) {
0056     HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0057     double rapidity = 0.5 * log((mom.e() + mom.pz()) / (mom.e() - mom.pz()));
0058
0059     if (abs((*p)->pdg_id()) == particleID && mom.rho() > minpcut && mom.rho() < maxpcut &&
0060         (*p)->momentum().perp() > minptcut && (*p)->momentum().perp() < maxptcut && mom.eta() > minetacut &&
0061         mom.eta() < maxetacut && rapidity > minrapcut && rapidity < maxrapcut && (*p)->momentum().phi() > minphicut &&
0062         (*p)->momentum().phi() < maxphicut) {
0063       HepMC::GenParticle* mother = (*((*p)->production_vertex()->particles_in_const_begin()));
0064
0065       // check various possible mothers
0066       for (auto motherID : motherIDs) {
0067         if (abs(mother->pdg_id()) == abs(motherID)) {
0068           // loop over its daughters
0069           for (HepMC::GenVertex::particle_iterator dau = mother->end_vertex()->particles_begin(HepMC::children);
0070                dau != mother->end_vertex()->particles_end(HepMC::children);
0071                ++dau) {
0072             // find the daugther you're interested in
0073             if (abs((*dau)->pdg_id()) == abs(sisterID)) {
0074               int failNephewPt = 0;
0075               // check pt of the nephews
0076               for (HepMC::GenVertex::particle_iterator nephew = (*dau)->end_vertex()->particles_begin(HepMC::children);
0077                    nephew != (*dau)->end_vertex()->particles_end(HepMC::children);
0078                    ++nephew) {
0079                 int nephew_pdgId = abs((*nephew)->pdg_id());
0080                 for (unsigned int i = 0; i < nephewIDs.size(); i++) {
0081                   if (nephew_pdgId == abs(nephewIDs.at(i)))
0082                     failNephewPt += ((*nephew)->momentum().perp() < minNephewPts.at(i));
0083                 }
0084               }
0085               if (failNephewPt > 0)
0086                 return false;
0087               // calculate displacement of the sister particle, from production to decay
0088               HepMC::GenVertex* v1 = (*dau)->production_vertex();
0089               HepMC::GenVertex* v2 = (*dau)->end_vertex();
0090
0091               double lx12 = v1->position().x() - v2->position().x();
0092               double ly12 = v1->position().y() - v2->position().y();
0093               double lxy12 = sqrt(lx12 * lx12 + ly12 * ly12);
0094
0095               if (maxSisDisplacement != -1) {
0096                 if (lxy12 < maxSisDisplacement) {
0097                   return true;
0098                 }
0099               } else {
0100                 return true;
0101               }
0102             }
0103           }
0104         }
0105       }
0106     }
0107   }
0108
0109   return false;
0110 }