GenFilters/plugins/PythiaAllDauVFilter.cc

0001 #include "PythiaAllDauVFilter.h"
0002 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0003 #include <iostream>
0004 #include <memory>
0005 #include <vector>
0006
0007 using namespace edm;
0008 using namespace std;
0009 using namespace Pythia8;
0010
0011 PythiaAllDauVFilter::PythiaAllDauVFilter(const edm::ParameterSet& iConfig)
0012     : fVerbose(iConfig.getUntrackedParameter("verbose", 0)),
0013       token_(consumes<edm::HepMCProduct>(iConfig.getUntrackedParameter<edm::InputTag>("moduleLabel"))),
0014       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0015       motherID(iConfig.getUntrackedParameter("MotherID", 0)),
0016       chargeconju(iConfig.getUntrackedParameter("ChargeConjugation", true)),
0017       ndaughters(iConfig.getUntrackedParameter("NumberDaughters", 0)),
0018       maxptcut(iConfig.getUntrackedParameter("MaxPt", 14000.)) {
0019   //now do what ever initialization is needed
0020   vector<int> defdauID;
0021   defdauID.push_back(0);
0022   dauIDs = iConfig.getUntrackedParameter<vector<int> >("DaughterIDs", defdauID);
0023   vector<double> defminptcut;
0024   defminptcut.push_back(0.);
0025   minptcut = iConfig.getUntrackedParameter<vector<double> >("MinPt", defminptcut);
0026   vector<double> defminetacut;
0027   defminetacut.push_back(-10.);
0028   minetacut = iConfig.getUntrackedParameter<vector<double> >("MinEta", defminetacut);
0029   vector<double> defmaxetacut;
0030   defmaxetacut.push_back(10.);
0031   maxetacut = iConfig.getUntrackedParameter<vector<double> >("MaxEta", defmaxetacut);
0032
0033   // create pythia8 instance to access particle data
0034   edm::LogInfo("PythiaAllDauVFilter") << "Creating pythia8 instance for particle properties" << endl;
0035   if (!fLookupGen.get())
0036     fLookupGen = std::make_unique<Pythia>();
0037
0038   if (chargeconju) {
0039     antiParticleID = -particleID;
0040     if (!(fLookupGen->particleData.isParticle(antiParticleID)))
0041       antiParticleID = particleID;
0042
0043     int antiId;
0044     for (size_t i = 0; i < dauIDs.size(); i++) {
0045       antiId = -dauIDs[i];
0046       if (!(fLookupGen->particleData.isParticle(antiId)))
0047         antiId = dauIDs[i];
0048
0049       antiDauIDs.push_back(antiId);
0050     }
0051   }
0052
0053   edm::LogInfo("PythiaAllDauVFilter") << "----------------------------------------------------------------------"
0054                                       << endl;
0055   edm::LogInfo("PythiaAllDauVFilter") << "--- PythiaAllDauVFilter" << endl;
0056   for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0057     edm::LogInfo("PythiaAllDauVFilter") << "ID: " << dauIDs[i] << " pT > " << minptcut[i] << " " << minetacut[i]
0058                                         << " eta < " << maxetacut[i] << endl;
0059   }
0060   if (chargeconju)
0061     for (unsigned int i = 0; i < antiDauIDs.size(); ++i) {
0062       edm::LogInfo("PythiaAllDauVFilter") << "ID: " << antiDauIDs[i] << " pT > " << minptcut[i] << " " << minetacut[i]
0063                                           << " eta < " << maxetacut[i] << endl;
0064     }
0065   edm::LogInfo("PythiaAllDauVFilter") << "maxptcut   = " << maxptcut << endl;
0066   edm::LogInfo("PythiaAllDauVFilter") << "particleID = " << particleID << endl;
0067   if (chargeconju)
0068     edm::LogInfo("PythiaAllDauVFilter") << "antiParticleID = " << antiParticleID << endl;
0069
0070   edm::LogInfo("PythiaAllDauVFilter") << "motherID   = " << motherID << endl;
0071 }
0072
0073 PythiaAllDauVFilter::~PythiaAllDauVFilter() {
0074   // do anything here that needs to be done at desctruction time
0075   // (e.g. close files, deallocate resources etc.)
0076 }
0077
0078 //
0079 // member functions
0080 //
0081
0082 // ------------ method called to produce the data  ------------
0083 bool PythiaAllDauVFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0084   using namespace edm;
0085   bool accepted = false;
0086   Handle<HepMCProduct> evt;
0087   iEvent.getByToken(token_, evt);
0088
0089   int OK(1);
0090   vector<int> vparticles;
0091   vector<bool> foundDaughter(dauIDs.size(), false);
0092   auto dauCollection = &dauIDs;
0093
0094   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0095
0096   if (fVerbose > 5) {
0097     edm::LogInfo("PythiaAllDauVFilter") << "looking for " << particleID << endl;
0098   }
0099
0100   for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0101     if ((*p)->pdg_id() == particleID) {
0102       dauCollection = &(dauIDs);
0103     } else if (chargeconju and ((*p)->pdg_id() == antiParticleID)) {
0104       dauCollection = &(antiDauIDs);
0105     } else {
0106       continue;
0107     }
0108
0109     // -- Check for mother of this particle
0110     if (0 != motherID) {
0111       OK = 0;
0112       for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0113            des != (*p)->production_vertex()->particles_in_const_end();
0114            ++des) {
0115         if (fVerbose > 10) {
0116           edm::LogInfo("PythiaAllDauVFilter") << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0117                                               << " eta: " << (*des)->momentum().eta() << endl;
0118         }
0119         if (abs(motherID) == abs((*des)->pdg_id())) {
0120           OK = 1;
0121           break;
0122         }
0123       }
0124     }
0125     if (0 == OK)
0126       continue;
0127
0128     // -- check for daugthers
0129     int ndau = 0;
0130     for (unsigned int i = 0; i < foundDaughter.size(); ++i) {
0131       foundDaughter[i] = false;
0132     }
0133     if (fVerbose > 5) {
0134       edm::LogInfo("PythiaAllDauVFilter") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0135                                           << " eta: " << (*p)->momentum().eta() << endl;
0136     }
0137     if ((*p)->end_vertex()) {
0138       for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0139            des != (*p)->end_vertex()->particles_end(HepMC::children);
0140            ++des) {
0141         ++ndau;
0142         if (fVerbose > 5) {
0143           edm::LogInfo("PythiaAllDauVFilter")
0144               << "\t daughter : ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0145               << " eta: " << (*des)->momentum().eta() << endl;
0146         }
0147         for (unsigned int i = 0; i < dauCollection->size(); ++i) {
0148           if ((*des)->pdg_id() != dauCollection->at(i))
0149             continue;
0150
0151           // possible to have more than one daughter of same pdgID and same/different kinematic constraints
0152           if (foundDaughter[i])
0153             continue;
0154
0155           if (fVerbose > 5) {
0156             edm::LogInfo("PythiaAllDauVFilter")
0157                 << "\t\t checking cuts of , daughter i = " << i << " pT = " << (*des)->momentum().perp()
0158                 << " eta = " << (*des)->momentum().eta() << endl;
0159           }
0160           if ((*des)->momentum().perp() > minptcut[i] && (*des)->momentum().perp() < maxptcut &&
0161               (*des)->momentum().eta() > minetacut[i] && (*des)->momentum().eta() < maxetacut[i]) {
0162             foundDaughter[i] = true;
0163             vparticles.push_back((*des)->pdg_id());
0164             if (fVerbose > 2) {
0165               edm::LogInfo("PythiaAllDauVFilter")
0166                   << "\t  accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0167                   << " eta = " << (*des)->momentum().eta() << endl;
0168             }
0169             break;
0170           }
0171         }
0172       }
0173     }
0174
0175     // -- ( number of daughtrs == daughters passing cut ) and ( all daughters specified are found)
0176     if (ndau == ndaughters) {
0177       accepted = true;
0178       for (unsigned int i = 0; i < foundDaughter.size(); ++i) {
0179         if (!foundDaughter[i]) {
0180           accepted = false;
0181         }
0182       }
0183       if (accepted and (fVerbose > 0)) {
0184         edm::LogInfo("PythiaAllDauVFilter") << "  accepted this decay from " << (*p)->pdg_id();
0185         for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0186           edm::LogInfo("PythiaAllDauVFilter") << vparticles[iv] << " ";
0187         edm::LogInfo("PythiaAllDauVFilter") << " from mother = " << motherID << endl;
0188       }
0189     }
0190
0191     if (accepted)
0192       break;
0193   }
0194
0195   delete myGenEvent;
0196   return accepted;
0197 }
0198
0199 //define this as a plug-in
0200 DEFINE_FWK_MODULE(PythiaAllDauVFilter);