|
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
|
// system include files
#include <memory>
#include <iostream>
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDFilter.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
using namespace edm;
using namespace std;
//
// class declaration
//
class LHEGenericMassFilter : public edm::global::EDFilter<> {
public:
explicit LHEGenericMassFilter(const edm::ParameterSet&);
~LHEGenericMassFilter() override = default;
static void fillDescriptions(edm::ConfigurationDescriptions&);
private:
bool filter(edm::StreamID, edm::Event&, edm::EventSetup const&) const override;
// ----------member data ---------------------------
const edm::EDGetTokenT<LHEEventProduct> src_;
const int numRequired_; // number of particles required to pass filter
const std::vector<int> particleID_; // vector of particle IDs to look for
const double minMass_;
const double maxMass_;
const bool requiredOutgoingStatus_; // Whether particles required to pass filter must have outgoing status
};
LHEGenericMassFilter::LHEGenericMassFilter(const edm::ParameterSet& iConfig)
: src_(consumes<LHEEventProduct>(iConfig.getParameter<edm::InputTag>("src"))),
numRequired_(iConfig.getParameter<int>("NumRequired")),
particleID_(iConfig.getParameter<std::vector<int>>("ParticleID")),
minMass_(iConfig.getParameter<double>("MinMass")),
maxMass_(iConfig.getParameter<double>("MaxMass")),
requiredOutgoingStatus_(iConfig.getParameter<bool>("RequiredOutgoingStatus")) {}
// ------------ method called to skim the data ------------
bool LHEGenericMassFilter::filter(edm::StreamID iID, edm::Event& iEvent, edm::EventSetup const& iSetup) const {
edm::Handle<LHEEventProduct> EvtHandle;
iEvent.getByToken(src_, EvtHandle);
int nFound = 0;
double Px = 0.;
double Py = 0.;
double Pz = 0.;
double E = 0.;
for (int i = 0; i < EvtHandle->hepeup().NUP; ++i) {
// if requiredOutgoingStatus_ keep only outgoing particles, otherwise keep them all
if (requiredOutgoingStatus_ && EvtHandle->hepeup().ISTUP[i] != 1) {
continue;
}
for (unsigned int j = 0; j < particleID_.size(); ++j) {
if (abs(particleID_[j]) == abs(EvtHandle->hepeup().IDUP[i])) {
nFound++;
Px = Px + EvtHandle->hepeup().PUP[i][0];
Py = Py + EvtHandle->hepeup().PUP[i][1];
Pz = Pz + EvtHandle->hepeup().PUP[i][2];
E = E + EvtHandle->hepeup().PUP[i][3];
break; // only match a given particle once!
}
} // loop over targets
} // loop over particles
// event accept/reject logic
if (nFound == numRequired_) {
double sqrdMass = E * E - (Px * Px + Py * Py + Pz * Pz);
if (sqrdMass > minMass_ * minMass_ && sqrdMass < maxMass_ * maxMass_) {
return true;
}
}
return false;
}
// ------------ method fills 'descriptions' with the allowed parameters for the module ------------
void LHEGenericMassFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("src", edm::InputTag("externalLHEProducer"));
desc.add<int>("NumRequired", 1);
desc.add<vector<int>>("ParticleID", std::vector<int>{1});
desc.add<double>("MinMass", 0.0);
desc.add<double>("MaxMass", 1.0);
desc.add<bool>("RequiredOutgoingStatus", true);
descriptions.addDefault(desc);
}
//define this as a plug-in
DEFINE_FWK_MODULE(LHEGenericMassFilter);
|