GenFilters/plugins/PythiaDauFilter.cc

0001
0002 #include "GeneratorInterface/GenFilters/plugins/PythiaDauFilter.h"
0003
0004 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0005 #include <iostream>
0006 #include <memory>
0007
0008 using namespace edm;
0009 using namespace std;
0010 using namespace Pythia8;
0011
0012 PythiaDauFilter::PythiaDauFilter(const edm::ParameterSet& iConfig)
0013     : token_(consumes<edm::HepMCProduct>(
0014           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0015       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0016       chargeconju(iConfig.getUntrackedParameter("ChargeConjugation", true)),
0017       ndaughters(iConfig.getUntrackedParameter("NumberDaughters", 0)),
0018       minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0019       maxptcut(iConfig.getUntrackedParameter("MaxPt", 14000.)),
0020       minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0021       maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)) {
0022   //now do what ever initialization is needed
0023   vector<int> defdauID;
0024   defdauID.push_back(0);
0025   dauIDs = iConfig.getUntrackedParameter<vector<int> >("DaughterIDs", defdauID);
0026   // create pythia8 instance to access particle data
0027   edm::LogInfo("PythiaDauFilter::PythiaDauFilter") << "Creating pythia8 instance for particle properties" << endl;
0028   if (!fLookupGen.get())
0029     fLookupGen = std::make_unique<Pythia>();
0030 }
0031
0032 PythiaDauFilter::~PythiaDauFilter() {
0033   // do anything here that needs to be done at desctruction time
0034   // (e.g. close files, deallocate resources etc.)
0035 }
0036
0037 //
0038 // member functions
0039 //
0040
0041 // ------------ method called to produce the data  ------------
0042 bool PythiaDauFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0043   using namespace edm;
0044   bool accepted = false;
0045   Handle<HepMCProduct> evt;
0046   iEvent.getByToken(token_, evt);
0047
0048   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0049
0050   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0051        ++p) {
0052     if ((*p)->pdg_id() != particleID)
0053       continue;
0054     int ndauac = 0;
0055     int ndau = 0;
0056     if ((*p)->end_vertex()) {
0057       for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0058            des != (*p)->end_vertex()->particles_end(HepMC::children);
0059            ++des) {
0060         ++ndau;
0061         for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0062           if ((*des)->pdg_id() != dauIDs[i])
0063             continue;
0064           if ((*des)->momentum().perp() > minptcut && (*des)->momentum().perp() < maxptcut &&
0065               (*des)->momentum().eta() > minetacut && (*des)->momentum().eta() < maxetacut) {
0066             ++ndauac;
0067             break;
0068           }
0069         }
0070       }
0071     }
0072     if (ndau == ndaughters && ndauac == ndaughters) {
0073       accepted = true;
0074       break;
0075     }
0076   }
0077
0078   if (!accepted && chargeconju) {
0079     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0080          ++p) {
0081       if ((*p)->pdg_id() != -particleID)
0082         continue;
0083       int ndauac = 0;
0084       int ndau = 0;
0085       if ((*p)->end_vertex()) {
0086         for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0087              des != (*p)->end_vertex()->particles_end(HepMC::children);
0088              ++des) {
0089           ++ndau;
0090           for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0091             int IDanti = -dauIDs[i];
0092             if (!(fLookupGen->particleData.isParticle(IDanti)))
0093               IDanti = dauIDs[i];
0094             if ((*des)->pdg_id() != IDanti)
0095               continue;
0096             if ((*des)->momentum().perp() > minptcut && (*des)->momentum().perp() < maxptcut &&
0097                 (*des)->momentum().eta() > minetacut && (*des)->momentum().eta() < maxetacut) {
0098               ++ndauac;
0099               break;
0100             }
0101           }
0102         }
0103       }
0104       if (ndau == ndaughters && ndauac == ndaughters) {
0105         accepted = true;
0106         break;
0107       }
0108     }
0109   }
0110
0111   if (accepted) {
0112     return true;
0113   } else {
0114     return false;
0115   }
0116 }