|
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
// -*- C++ -*-
//
// Package: MCSingleParticleYPt
// Class: MCSingleParticleYPt
//
/*
Description: filter events based on the Pythia particleID, Pt, Y and status. It is based on MCSingleParticleFilter.
It will used to filter a b-hadron with the given kinematics, only one b-hadron is required to match.
Implementation: inherits from generic EDFilter
*/
//
// Author: Alberto Sanchez-Hernandez
// Adapted on: August 2016
//
//
// Filter based on MCSingleParticleFilter.cc, but using rapidity instead of eta
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDFilter.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/EDGetToken.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include <cmath>
#include <ostream>
#include <string>
#include <vector>
class MCSingleParticleYPt : public edm::global::EDFilter<> {
public:
explicit MCSingleParticleYPt(const edm::ParameterSet&);
bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
private:
const edm::EDGetTokenT<edm::HepMCProduct> token_;
const int fVerbose;
const bool fchekantiparticle;
std::vector<int> particleID;
std::vector<double> ptMin;
std::vector<double> rapMin;
std::vector<double> rapMax;
std::vector<int> status;
};
using namespace edm;
using namespace std;
MCSingleParticleYPt::MCSingleParticleYPt(const edm::ParameterSet& iConfig)
: token_(consumes<edm::HepMCProduct>(
edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
fVerbose(iConfig.getUntrackedParameter("verbose", 0)),
fchekantiparticle(iConfig.getUntrackedParameter("CheckAntiparticle", true)) {
vector<int> defpid;
defpid.push_back(0);
particleID = iConfig.getUntrackedParameter<vector<int> >("ParticleID", defpid);
vector<double> defptmin;
defptmin.push_back(0.);
ptMin = iConfig.getUntrackedParameter<vector<double> >("MinPt", defptmin);
vector<double> defrapmin;
defrapmin.push_back(-10.);
rapMin = iConfig.getUntrackedParameter<vector<double> >("MinY", defrapmin);
vector<double> defrapmax;
defrapmax.push_back(10.);
rapMax = iConfig.getUntrackedParameter<vector<double> >("MaxY", defrapmax);
vector<int> defstat;
defstat.push_back(0);
status = iConfig.getUntrackedParameter<vector<int> >("Status", defstat);
// check for same size
if ((ptMin.size() > 1 && particleID.size() != ptMin.size()) ||
(rapMin.size() > 1 && particleID.size() != rapMin.size()) ||
(rapMax.size() > 1 && particleID.size() != rapMax.size()) ||
(status.size() > 1 && particleID.size() != status.size())) {
edm::LogWarning("MCSingleParticleYPt")
<< "WARNING: MCSingleParticleYPt : size of vector cuts do not match!!" << endl;
}
// if ptMin size smaller than particleID , fill up further with defaults
if (particleID.size() > ptMin.size()) {
vector<double> defptmin2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defptmin2.push_back(0.);
}
ptMin = defptmin2;
}
// if etaMin size smaller than particleID , fill up further with defaults
if (particleID.size() > rapMin.size()) {
vector<double> defrapmin2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defrapmin2.push_back(-10.);
}
rapMin = defrapmin2;
}
// if etaMax size smaller than particleID , fill up further with defaults
if (particleID.size() > rapMax.size()) {
vector<double> defrapmax2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defrapmax2.push_back(10.);
}
rapMax = defrapmax2;
}
// if status size smaller than particleID , fill up further with defaults
if (particleID.size() > status.size()) {
vector<int> defstat2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defstat2.push_back(0);
}
status = defstat2;
}
if (fVerbose > 0) {
edm::LogInfo("MCSingleParticleYPt") << "----------------------------------------------------------------------"
<< std::endl;
edm::LogInfo("MCSingleParticleYPt") << "----- MCSingleParticleYPt" << std::endl;
for (unsigned int i = 0; i < particleID.size(); ++i) {
edm::LogInfo("MCSingleParticleYPt") << " ID: " << particleID[i] << " pT > " << ptMin[i] << ", " << rapMin[i]
<< " < y < " << rapMax[i] << ", status = " << status[i] << std::endl;
}
if (fchekantiparticle)
edm::LogInfo("MCSingleParticleYPt") << " anti-particles will be tested as well." << std::endl;
edm::LogInfo("MCSingleParticleYPt") << "----------------------------------------------------------------------"
<< std::endl;
}
}
bool MCSingleParticleYPt::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
bool accepted = false;
Handle<HepMCProduct> evt;
iEvent.getByToken(token_, evt);
const HepMC::GenEvent* myGenEvent = evt->GetEvent();
for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
++p) {
if (fVerbose > 3)
edm::LogInfo("MCSingleParticleYPt")
<< "Looking at particle : " << (*p)->pdg_id() << " status : " << (*p)->status() << std::endl;
for (unsigned int i = 0; i < particleID.size(); i++) {
if (particleID[i] == (*p)->pdg_id() || (fchekantiparticle && (-particleID[i] == (*p)->pdg_id())) ||
particleID[i] == 0) {
// calculate rapidity just for the desired particle and make sure, this particles has enough energy
double rapidity = ((*p)->momentum().e() - (*p)->momentum().pz()) > 0.
? 0.5 * log(((*p)->momentum().e() + (*p)->momentum().pz()) /
((*p)->momentum().e() - (*p)->momentum().pz()))
: rapMax[i] + .1;
if (fVerbose > 2)
edm::LogInfo("MCSingleParticleYPt")
<< "Testing particle : " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp() << " y: " << rapidity
<< " status : " << (*p)->status() << endl;
if ((*p)->momentum().perp() > ptMin[i] && rapidity > rapMin[i] && rapidity < rapMax[i] &&
((*p)->status() == status[i] || status[i] == 0)) {
accepted = true;
if (fVerbose > 1)
edm::LogInfo("MCSingleParticleYPt")
<< "Accepted particle : " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp() << " y: " << rapidity
<< " status : " << (*p)->status() << endl;
break;
}
}
}
if (accepted)
break;
}
return accepted;
}
DEFINE_FWK_MODULE(MCSingleParticleYPt);
|