Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-22 06:25:02

0001 // -*- C++ -*-
0002 //
0003 // Package:    HighMultiplicityGenFilter
0004 // Class:      HighMultiplicityGenFilter
0005 //
0006 /**\class HighMultiplicityGenFilter HighMultiplicityGenFilter.cc davidlw/HighMultiplicityGenFilter/src/HighMultiplicityGenFilter.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Wei Li
0015 //         Created:  Tue Dec  8 23:51:37 EST 2009
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDFilter.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/ESHandle.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/Utilities/interface/ESGetToken.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0032 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0033 //
0034 // class declaration
0035 //
0036 
0037 class HighMultiplicityGenFilter : public edm::one::EDFilter<> {
0038 public:
0039   explicit HighMultiplicityGenFilter(const edm::ParameterSet&);
0040   ~HighMultiplicityGenFilter() override = default;
0041 
0042 private:
0043   void beginJob() override;
0044   bool filter(edm::Event&, const edm::EventSetup&) override;
0045   void endJob() override;
0046 
0047   // ----------member data ---------------------------
0048   const edm::ESGetToken<ParticleDataTable, edm::DefaultRecord> pdtToken_;
0049   const edm::EDGetTokenT<edm::HepMCProduct> hepmcSrc;
0050   const double etaMax;
0051   const double ptMin;
0052   const int nMin;
0053   int nAccepted;
0054 };
0055 
0056 //
0057 // constants, enums and typedefs
0058 //
0059 
0060 //
0061 // static data member definitions
0062 //
0063 
0064 //
0065 // constructors and destructor
0066 //
0067 HighMultiplicityGenFilter::HighMultiplicityGenFilter(const edm::ParameterSet& iConfig)
0068     : pdtToken_(esConsumes<ParticleDataTable, edm::DefaultRecord>()),
0069       hepmcSrc(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("generatorSmeared"))),
0070       etaMax(iConfig.getUntrackedParameter<double>("etaMax")),
0071       ptMin(iConfig.getUntrackedParameter<double>("ptMin")),
0072       nMin(iConfig.getUntrackedParameter<int>("nMin")) {
0073   //now do what ever initialization is needed
0074   nAccepted = 0;
0075 }
0076 
0077 //
0078 // member functions
0079 //
0080 
0081 // ------------ method called on each new Event  ------------
0082 bool HighMultiplicityGenFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0083   bool accepted = false;
0084   const edm::Handle<edm::HepMCProduct>& evt = iEvent.getHandle(hepmcSrc);
0085   const edm::ESHandle<ParticleDataTable>& pdt = iSetup.getHandle(pdtToken_);
0086 
0087   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0088 
0089   int nMult = 0;
0090   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0091        ++p) {
0092     if ((*p)->status() != 1)
0093       continue;
0094 
0095     double charge = 0;
0096     int pid = (*p)->pdg_id();
0097     if (abs(pid) > 100000) {
0098       std::cout << "pid=" << pid << " status=" << (*p)->status() << std::endl;
0099       continue;
0100     }
0101     const ParticleData* part = pdt->particle(pid);
0102     if (part)
0103       charge = part->charge();
0104     if (charge == 0)
0105       continue;
0106 
0107     if ((*p)->momentum().perp() > ptMin && fabs((*p)->momentum().eta()) < etaMax)
0108       nMult++;
0109   }
0110   if (nMult >= nMin) {
0111     nAccepted++;
0112     accepted = true;
0113   }
0114   return accepted;
0115 }
0116 
0117 // ------------ method called once each job just before starting event loop  ------------
0118 void HighMultiplicityGenFilter::beginJob() {}
0119 
0120 // ------------ method called once each job just after ending the event loop  ------------
0121 void HighMultiplicityGenFilter::endJob() {
0122   std::cout << "There are " << nAccepted << " events with multiplicity greater than " << nMin << std::endl;
0123 }
0124 
0125 //define this as a plug-in
0126 DEFINE_FWK_MODULE(HighMultiplicityGenFilter);