File indexing completed on 2024-04-06 12:13:44
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
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
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
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
0058
0059
0060
0061
0062
0063
0064
0065
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
0074 nAccepted = 0;
0075 }
0076
0077
0078
0079
0080
0081
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
0118 void HighMultiplicityGenFilter::beginJob() {}
0119
0120
0121 void HighMultiplicityGenFilter::endJob() {
0122 std::cout << "There are " << nAccepted << " events with multiplicity greater than " << nMin << std::endl;
0123 }
0124
0125
0126 DEFINE_FWK_MODULE(HighMultiplicityGenFilter);