|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
// -*- C++ -*-
//
// Package: HighMultiplicityGenFilter
// Class: HighMultiplicityGenFilter
//
/**\class HighMultiplicityGenFilter HighMultiplicityGenFilter.cc davidlw/HighMultiplicityGenFilter/src/HighMultiplicityGenFilter.cc
Description: <one line class summary>
Implementation:
<Notes on implementation>
*/
//
// Original Author: Wei Li
// Created: Tue Dec 8 23:51:37 EST 2009
//
//
// system include files
#include <memory>
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDFilter.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Utilities/interface/ESGetToken.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
//
// class declaration
//
class HighMultiplicityGenFilter : public edm::one::EDFilter<> {
public:
explicit HighMultiplicityGenFilter(const edm::ParameterSet&);
~HighMultiplicityGenFilter() override = default;
private:
void beginJob() override;
bool filter(edm::Event&, const edm::EventSetup&) override;
void endJob() override;
// ----------member data ---------------------------
const edm::ESGetToken<ParticleDataTable, edm::DefaultRecord> pdtToken_;
const edm::EDGetTokenT<edm::HepMCProduct> hepmcSrc;
const double etaMax;
const double ptMin;
const int nMin;
int nAccepted;
};
//
// constants, enums and typedefs
//
//
// static data member definitions
//
//
// constructors and destructor
//
HighMultiplicityGenFilter::HighMultiplicityGenFilter(const edm::ParameterSet& iConfig)
: pdtToken_(esConsumes<ParticleDataTable, edm::DefaultRecord>()),
hepmcSrc(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("generatorSmeared"))),
etaMax(iConfig.getUntrackedParameter<double>("etaMax")),
ptMin(iConfig.getUntrackedParameter<double>("ptMin")),
nMin(iConfig.getUntrackedParameter<int>("nMin")) {
//now do what ever initialization is needed
nAccepted = 0;
}
//
// member functions
//
// ------------ method called on each new Event ------------
bool HighMultiplicityGenFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
bool accepted = false;
const edm::Handle<edm::HepMCProduct>& evt = iEvent.getHandle(hepmcSrc);
const edm::ESHandle<ParticleDataTable>& pdt = iSetup.getHandle(pdtToken_);
const HepMC::GenEvent* myGenEvent = evt->GetEvent();
int nMult = 0;
for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
++p) {
if ((*p)->status() != 1)
continue;
double charge = 0;
int pid = (*p)->pdg_id();
if (abs(pid) > 100000) {
std::cout << "pid=" << pid << " status=" << (*p)->status() << std::endl;
continue;
}
const ParticleData* part = pdt->particle(pid);
if (part)
charge = part->charge();
if (charge == 0)
continue;
if ((*p)->momentum().perp() > ptMin && fabs((*p)->momentum().eta()) < etaMax)
nMult++;
}
if (nMult >= nMin) {
nAccepted++;
accepted = true;
}
return accepted;
}
// ------------ method called once each job just before starting event loop ------------
void HighMultiplicityGenFilter::beginJob() {}
// ------------ method called once each job just after ending the event loop ------------
void HighMultiplicityGenFilter::endJob() {
std::cout << "There are " << nAccepted << " events with multiplicity greater than " << nMin << std::endl;
}
//define this as a plug-in
DEFINE_FWK_MODULE(HighMultiplicityGenFilter);
|