GenFilters/plugins/MCSingleParticleFilter.cc

0001
0002 #include "GeneratorInterface/GenFilters/plugins/MCSingleParticleFilter.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 MCSingleParticleFilter::MCSingleParticleFilter(const edm::ParameterSet& iConfig)
0012     : token_(consumes<edm::HepMCProduct>(
0013           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0014       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
0015   //here do whatever other initialization is needed
0016   vector<int> defpid;
0017   defpid.push_back(0);
0018   particleID = iConfig.getUntrackedParameter<vector<int> >("ParticleID", defpid);
0019   vector<double> defptmin;
0020   defptmin.push_back(0.);
0021   ptMin = iConfig.getUntrackedParameter<vector<double> >("MinPt", defptmin);
0022
0023   vector<double> defetamin;
0024   defetamin.push_back(-10.);
0025   etaMin = iConfig.getUntrackedParameter<vector<double> >("MinEta", defetamin);
0026   vector<double> defetamax;
0027   defetamax.push_back(10.);
0028   etaMax = iConfig.getUntrackedParameter<vector<double> >("MaxEta", defetamax);
0029   vector<int> defstat;
0030   defstat.push_back(0);
0031   status = iConfig.getUntrackedParameter<vector<int> >("Status", defstat);
0032
0033   // check for same size
0034   if ((ptMin.size() > 1 && particleID.size() != ptMin.size()) ||
0035       (etaMin.size() > 1 && particleID.size() != etaMin.size()) ||
0036       (etaMax.size() > 1 && particleID.size() != etaMax.size()) ||
0037       (status.size() > 1 && particleID.size() != status.size())) {
0038     edm::LogInfo("MCSingleParticleFilter")
0039         << "WARNING: size of MinPthat and/or MaxPthat not matching with ProcessID size!!" << endl;
0040   }
0041
0042   // if ptMin size smaller than particleID , fill up further with defaults
0043   if (particleID.size() > ptMin.size()) {
0044     vector<double> defptmin2;
0045     for (unsigned int i = 0; i < particleID.size(); i++) {
0046       defptmin2.push_back(0.);
0047     }
0048     ptMin = defptmin2;
0049   }
0050   // if etaMin size smaller than particleID , fill up further with defaults
0051   if (particleID.size() > etaMin.size()) {
0052     vector<double> defetamin2;
0053     for (unsigned int i = 0; i < particleID.size(); i++) {
0054       defetamin2.push_back(-10.);
0055     }
0056     etaMin = defetamin2;
0057   }
0058   // if etaMax size smaller than particleID , fill up further with defaults
0059   if (particleID.size() > etaMax.size()) {
0060     vector<double> defetamax2;
0061     for (unsigned int i = 0; i < particleID.size(); i++) {
0062       defetamax2.push_back(10.);
0063     }
0064     etaMax = defetamax2;
0065   }
0066   // if status size smaller than particleID , fill up further with defaults
0067   if (particleID.size() > status.size()) {
0068     vector<int> defstat2;
0069     for (unsigned int i = 0; i < particleID.size(); i++) {
0070       defstat2.push_back(0);
0071     }
0072     status = defstat2;
0073   }
0074
0075   // check if beta is smaller than 1
0076   if (std::abs(betaBoost) >= 1) {
0077     edm::LogError("MCSingleParticleFilter") << "Input beta boost is >= 1 !";
0078   }
0079 }
0080
0081 MCSingleParticleFilter::~MCSingleParticleFilter() {
0082   // do anything here that needs to be done at desctruction time
0083   // (e.g. close files, deallocate resources etc.)
0084 }
0085
0086 // ------------ method called to skim the data  ------------
0087 bool MCSingleParticleFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0088   using namespace edm;
0089   bool accepted = false;
0090   Handle<HepMCProduct> evt;
0091   iEvent.getByToken(token_, evt);
0092
0093   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0094
0095   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0096        ++p) {
0097     for (unsigned int i = 0; i < particleID.size(); i++) {
0098       if (particleID[i] == (*p)->pdg_id() || particleID[i] == 0) {
0099         if ((*p)->momentum().perp() > ptMin[i] && ((*p)->status() == status[i] || status[i] == 0)) {
0100           HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0101           if (mom.eta() > etaMin[i] && mom.eta() < etaMax[i]) {
0102             accepted = true;
0103           }
0104         }
0105       }
0106     }
0107   }
0108
0109   return accepted;
0110 }