GenFilters/plugins/MCProcessFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MCProcessFilter
0004 // Class:      MCProcessFilter
0005 //
0006 /*
0007
0008  Description: filter events based on the Pythia ProcessID and the Pt_hat
0009
0010  Implementation: inherits from generic EDFilter
0011
0012 */
0013 //
0014 // Original Author:  Filip Moortgat
0015 //         Created:  Mon Sept 11 10:57:54 CET 2006
0016 //
0017
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/global/EDFilter.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/EDGetToken.h"
0025 #include "FWCore/Utilities/interface/InputTag.h"
0026 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0027
0028 #include <iostream>
0029 #include <string>
0030 #include <vector>
0031
0032 class MCProcessFilter : public edm::global::EDFilter<> {
0033 public:
0034   explicit MCProcessFilter(const edm::ParameterSet&);
0035
0036   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0037
0038 private:
0039   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0040   std::vector<int> processID;
0041   std::vector<double> pthatMin;
0042   std::vector<double> pthatMax;
0043 };
0044
0045 using namespace std;
0046
0047 MCProcessFilter::MCProcessFilter(const edm::ParameterSet& iConfig)
0048     : token_(consumes<edm::HepMCProduct>(
0049           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))) {
0050   vector<int> defproc;
0051   defproc.push_back(0);
0052   processID = iConfig.getUntrackedParameter<vector<int> >("ProcessID", defproc);
0053   vector<double> defpthatmin;
0054   defpthatmin.push_back(0.);
0055   pthatMin = iConfig.getUntrackedParameter<vector<double> >("MinPthat", defpthatmin);
0056   vector<double> defpthatmax;
0057   defpthatmax.push_back(10000.);
0058   pthatMax = iConfig.getUntrackedParameter<vector<double> >("MaxPthat", defpthatmax);
0059
0060   // checkin size of phthat vectors -- default is allowed
0061   if ((pthatMin.size() > 1 && processID.size() != pthatMin.size()) ||
0062       (pthatMax.size() > 1 && processID.size() != pthatMax.size())) {
0063     cout << "WARNING: MCPROCESSFILTER : size of MinPthat and/or MaxPthat not matching with ProcessID size!!" << endl;
0064   }
0065
0066   // if pthatMin size smaller than processID , fill up further with defaults
0067   if (processID.size() > pthatMin.size()) {
0068     vector<double> defpthatmin2;
0069     for (unsigned int i = 0; i < processID.size(); i++) {
0070       defpthatmin2.push_back(0.);
0071     }
0072     pthatMin = defpthatmin2;
0073   }
0074   // if pthatMax size smaller than processID , fill up further with defaults
0075   if (processID.size() > pthatMax.size()) {
0076     vector<double> defpthatmax2;
0077     for (unsigned int i = 0; i < processID.size(); i++) {
0078       defpthatmax2.push_back(10000.);
0079     }
0080     pthatMax = defpthatmax2;
0081   }
0082 }
0083
0084 bool MCProcessFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0085   bool accepted = false;
0086   edm::Handle<edm::HepMCProduct> evt;
0087   iEvent.getByToken(token_, evt);
0088
0089   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0090
0091   // do the selection -- processID 0 is always accepted
0092   for (unsigned int i = 0; i < processID.size(); i++) {
0093     if (processID[i] == myGenEvent->signal_process_id() || processID[i] == 0) {
0094       if (myGenEvent->event_scale() > pthatMin[i] && myGenEvent->event_scale() < pthatMax[i]) {
0095         accepted = true;
0096       }
0097     }
0098   }
0099   return accepted;
0100 }
0101
0102 DEFINE_FWK_MODULE(MCProcessFilter);