|
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
127
128
129
130
131
132
133
134
135
136
|
// -*- C++ -*-
//
// Package: MCDisplacementFilter
// Class: MCDisplacementFilter
//
// system include files
//#include <memory>
//#include <iostream>
#include <vector>
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDFilter.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
using namespace edm;
using namespace std;
//Filter particles based on their minimum and/or maximum displacement on the transverse plane and optionally on their pdgIds
//To run independently of pdgId, do not insert the particleIDs entry in filter declaration
// class decleration
//
namespace edm {
class HepMCProduct;
class ConfigurationDescriptions;
} // namespace edm
class MCDisplacementFilter : public edm::global::EDFilter<> {
public:
explicit MCDisplacementFilter(const edm::ParameterSet&);
~MCDisplacementFilter() override = default;
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
private:
// ----------member data ---------------------------
edm::InputTag hepMCProductTag_;
const edm::EDGetTokenT<edm::HepMCProduct> token_;
// To chose on which pdgIds the filter is applied - if ParticleIDs.at(0)==0 runs on all particles
std::vector<int> particleIDs_;
const float theUpperCut_; // Maximum displacement accepted
const float theLowerCut_; // Minimum displacement accepted
};
//methods implementation
//
//Class initialization
MCDisplacementFilter::MCDisplacementFilter(const edm::ParameterSet& iConfig)
: hepMCProductTag_(iConfig.getParameter<edm::InputTag>("hepMCProductTag")),
token_(consumes<edm::HepMCProduct>(hepMCProductTag_)),
particleIDs_(iConfig.getParameter<std::vector<int>>("ParticleIDs")),
theUpperCut_(iConfig.getParameter<double>("LengMax")),
theLowerCut_(iConfig.getParameter<double>("LengMin")) {}
//Filter description
void MCDisplacementFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("hepMCProductTag", edm::InputTag("generator", "unsmeared"));
desc.add<std::vector<int>>("ParticleIDs", std::vector<int>{0});
desc.add<double>("LengMax", -1.);
desc.add<double>("LengMin", -1.);
descriptions.addDefault(desc);
}
// ------------ method called to skim the data ------------
bool MCDisplacementFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
using namespace edm;
Handle<HepMCProduct> evt;
iEvent.getByToken(token_, evt);
bool pass = false;
bool matchedID = true;
const float theUpperCut2 = theUpperCut_ * theUpperCut_;
const float theLowerCut2 = theLowerCut_ * theLowerCut_;
const HepMC::GenEvent* generated_event = evt->GetEvent();
HepMC::GenEvent::particle_const_iterator p;
for (p = generated_event->particles_begin(); p != generated_event->particles_end(); p++) {
//matchedID might be moved to false to true for a particle in the event, it needs to be resetted everytime
if (particleIDs_.at(0) != 0)
matchedID = false;
//if a list of pdgId is provided, loop only on particles with those pdgId
for (unsigned int idx = 0; idx < particleIDs_.size(); idx++) {
if (abs((*p)->pdg_id()) == abs(particleIDs_.at(idx))) { //compares absolute values of pdgIds
matchedID = true;
break;
}
}
if (matchedID) {
if (theLowerCut_ <= 0. && theUpperCut_ <= 0.) {
pass = true;
break;
}
if (((*p)->production_vertex() != nullptr) && ((*p)->end_vertex() != nullptr)) {
float dist2 = (((*p)->production_vertex())->position().x() - ((*p)->end_vertex())->position().x()) *
(((*p)->production_vertex())->position().x() - ((*p)->end_vertex())->position().x()) +
(((*p)->production_vertex())->position().y() - ((*p)->end_vertex())->position().y()) *
(((*p)->production_vertex())->position().y() - ((*p)->end_vertex())->position().y());
// lower/upper cuts can be also <= 0 - prompt particle needs to be accepted in that case
if ((theLowerCut_ <= 0. || dist2 >= theLowerCut2) && (theUpperCut_ <= 0. || dist2 < theUpperCut2)) {
pass = true;
break;
}
}
if (((*p)->production_vertex() == nullptr) && (((*p)->end_vertex() != nullptr))) {
// lower/upper cuts can be also 0 - prompt particle needs to be accepted in that case
float distEndVert = (*p)->end_vertex()->position().perp();
if ((theLowerCut_ <= 0. || distEndVert >= theLowerCut_) && (theUpperCut_ <= 0. || distEndVert < theUpperCut_)) {
pass = true;
break;
}
}
}
}
return pass;
}
DEFINE_FWK_MODULE(MCDisplacementFilter);
|