1 |
peiffer |
1.1 |
// -*- C++ -*-
|
2 |
|
|
//
|
3 |
|
|
// Package: NtupleWriter
|
4 |
|
|
// Class: NtupleWriter
|
5 |
|
|
//
|
6 |
|
|
/**\class NtupleWriter NtupleWriter.cc NtupleWriter/src/NtupleWriter.cc
|
7 |
|
|
|
8 |
|
|
Description: [one line class summary]
|
9 |
|
|
|
10 |
|
|
Implementation:
|
11 |
|
|
[Notes on implementation]
|
12 |
|
|
*/
|
13 |
|
|
//
|
14 |
|
|
// Original Author: Thomas Peiffer,,,Uni Hamburg
|
15 |
|
|
// Created: Tue Mar 13 08:43:34 CET 2012
|
16 |
|
|
// $Id$
|
17 |
|
|
//
|
18 |
|
|
//
|
19 |
|
|
|
20 |
|
|
#include "UserCode/UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
|
21 |
|
|
|
22 |
|
|
|
23 |
|
|
|
24 |
|
|
//
|
25 |
|
|
// constants, enums and typedefs
|
26 |
|
|
//
|
27 |
|
|
|
28 |
|
|
//
|
29 |
|
|
// static data member definitions
|
30 |
|
|
//
|
31 |
|
|
|
32 |
|
|
//
|
33 |
|
|
// constructors and destructor
|
34 |
|
|
//
|
35 |
|
|
NtupleWriter::NtupleWriter(const edm::ParameterSet& iConfig)
|
36 |
|
|
|
37 |
|
|
{
|
38 |
|
|
//now do what ever initialization is needed
|
39 |
|
|
//edm::Service<TFileService> fs;
|
40 |
|
|
//tr = fs->make<TTree>("AnalysisTree","AnalysisTree");
|
41 |
|
|
|
42 |
|
|
fileName = iConfig.getParameter<std::string>("fileName");
|
43 |
|
|
outfile = new TFile(fileName,"RECREATE");
|
44 |
|
|
outfile->cd();
|
45 |
|
|
tr = new TTree("AnalysisTree","AnalysisTree");
|
46 |
|
|
|
47 |
|
|
std::string name;
|
48 |
|
|
|
49 |
|
|
doElectrons = iConfig.getParameter<bool>("doElectrons");
|
50 |
|
|
doMuons = iConfig.getParameter<bool>("doMuons");
|
51 |
|
|
doTaus = iConfig.getParameter<bool>("doTaus");
|
52 |
|
|
doJets = iConfig.getParameter<bool>("doJets");
|
53 |
|
|
doPhotons = iConfig.getParameter<bool>("doPhotons");
|
54 |
|
|
doMET = iConfig.getParameter<bool>("doMET");
|
55 |
|
|
doGenInfo = iConfig.getParameter<bool>("doGenInfo");
|
56 |
|
|
doPV = iConfig.getParameter<bool>("doPV");
|
57 |
|
|
doTopJets = iConfig.getParameter<bool>("doTopJets");
|
58 |
|
|
|
59 |
|
|
// initialization of tree variables
|
60 |
|
|
|
61 |
|
|
tr->Branch("run",&run);
|
62 |
|
|
tr->Branch("event",&event);
|
63 |
|
|
tr->Branch("luminosityBlock",&luminosityBlock);
|
64 |
|
|
tr->Branch("isRealData",&isRealData);
|
65 |
|
|
tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
|
66 |
|
|
tr->Branch("beamspot_x0",&beamspot_x0);
|
67 |
|
|
tr->Branch("beamspot_y0",&beamspot_y0);
|
68 |
|
|
tr->Branch("beamspot_z0",&beamspot_z0);
|
69 |
|
|
|
70 |
|
|
if(doElectrons){
|
71 |
|
|
electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
|
72 |
|
|
for(size_t j=0; j< electron_sources.size(); ++j){
|
73 |
|
|
tr->Branch( electron_sources[j].c_str(), "std::vector<Electron>", &eles[j]);
|
74 |
|
|
}
|
75 |
|
|
}
|
76 |
|
|
if(doMuons){
|
77 |
|
|
muon_sources = iConfig.getParameter<std::vector<std::string> >("muon_sources");
|
78 |
|
|
for(size_t j=0; j< muon_sources.size(); ++j){
|
79 |
|
|
tr->Branch( muon_sources[j].c_str(), "std::vector<Muon>", &mus[j]);
|
80 |
|
|
}
|
81 |
|
|
}
|
82 |
|
|
if(doTaus){
|
83 |
|
|
tau_sources = iConfig.getParameter<std::vector<std::string> >("tau_sources");
|
84 |
|
|
tau_ptmin = iConfig.getParameter<double> ("tau_ptmin");
|
85 |
|
|
tau_etamax = iConfig.getParameter<double> ("tau_etamax");
|
86 |
|
|
for(size_t j=0; j< tau_sources.size(); ++j){
|
87 |
|
|
tr->Branch( tau_sources[j].c_str(), "std::vector<Tau>", &taus[j]);
|
88 |
|
|
}
|
89 |
|
|
}
|
90 |
|
|
if(doJets){
|
91 |
|
|
jet_sources = iConfig.getParameter<std::vector<std::string> >("jet_sources");
|
92 |
|
|
jet_ptmin = iConfig.getParameter<double> ("jet_ptmin");
|
93 |
|
|
jet_etamax = iConfig.getParameter<double> ("jet_etamax");
|
94 |
|
|
for(size_t j=0; j< jet_sources.size(); ++j){
|
95 |
|
|
tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
|
96 |
|
|
}
|
97 |
|
|
}
|
98 |
|
|
if(doTopJets){
|
99 |
|
|
topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
|
100 |
|
|
topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
|
101 |
|
|
topjet_etamax = iConfig.getParameter<double> ("topjet_etamax");
|
102 |
|
|
for(size_t j=0; j< topjet_sources.size(); ++j){
|
103 |
|
|
tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
|
104 |
|
|
}
|
105 |
|
|
}
|
106 |
|
|
if(doPhotons){
|
107 |
|
|
photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
|
108 |
|
|
for(size_t j=0; j< photon_sources.size(); ++j){
|
109 |
|
|
tr->Branch( photon_sources[j].c_str(), "std::vector<Photon>", &phs[j]);
|
110 |
|
|
}
|
111 |
|
|
}
|
112 |
|
|
if(doMET){
|
113 |
|
|
met_sources = iConfig.getParameter<std::vector<std::string> >("met_sources");
|
114 |
|
|
for(size_t j=0; j< met_sources.size(); ++j){
|
115 |
|
|
tr->Branch(met_sources[j].c_str(), "MET", &met[j]);
|
116 |
|
|
}
|
117 |
|
|
}
|
118 |
|
|
if(doPV){
|
119 |
|
|
pv_sources = iConfig.getParameter<std::vector<std::string> >("pv_sources");
|
120 |
|
|
for(size_t j=0; j< pv_sources.size(); ++j){
|
121 |
|
|
tr->Branch( pv_sources[j].c_str(), "std::vector<PrimaryVertex>", &pvs[j]);
|
122 |
|
|
}
|
123 |
|
|
}
|
124 |
|
|
if(doGenInfo){
|
125 |
|
|
tr->Branch("genInfo","GenInfo",&genInfo);
|
126 |
|
|
tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
|
127 |
|
|
}
|
128 |
|
|
|
129 |
|
|
trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
|
130 |
|
|
//tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
|
131 |
|
|
tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
|
132 |
|
|
tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
|
133 |
|
|
tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
|
134 |
|
|
tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
|
135 |
|
|
|
136 |
|
|
newrun = true;
|
137 |
|
|
}
|
138 |
|
|
|
139 |
|
|
|
140 |
|
|
NtupleWriter::~NtupleWriter()
|
141 |
|
|
{
|
142 |
|
|
|
143 |
|
|
// do anything here that needs to be done at desctruction time
|
144 |
|
|
// (e.g. close files, deallocate resources etc.)
|
145 |
|
|
|
146 |
|
|
}
|
147 |
|
|
|
148 |
|
|
|
149 |
|
|
//
|
150 |
|
|
// member functions
|
151 |
|
|
//
|
152 |
|
|
|
153 |
|
|
// ------------ method called for each event ------------
|
154 |
|
|
void
|
155 |
|
|
NtupleWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
156 |
|
|
{
|
157 |
|
|
// using namespace edm;
|
158 |
|
|
|
159 |
|
|
|
160 |
|
|
// ------------- common variables -------------
|
161 |
|
|
|
162 |
|
|
run = iEvent.id().run();
|
163 |
|
|
event = iEvent.id().event();
|
164 |
|
|
luminosityBlock = iEvent.luminosityBlock();
|
165 |
|
|
isRealData = iEvent.isRealData();
|
166 |
|
|
|
167 |
|
|
if(isRealData){
|
168 |
|
|
edm::Handle<bool> bool_handle;
|
169 |
|
|
iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
|
170 |
|
|
HBHENoiseFilterResult = *(bool_handle.product());
|
171 |
|
|
}
|
172 |
|
|
else HBHENoiseFilterResult = false;
|
173 |
|
|
|
174 |
|
|
// ------------- electrons -------------
|
175 |
|
|
if(doElectrons){
|
176 |
|
|
for(size_t j=0; j< electron_sources.size(); ++j){
|
177 |
|
|
eles[j].clear();
|
178 |
|
|
edm::Handle< std::vector<pat::Electron> > ele_handle;
|
179 |
|
|
iEvent.getByLabel(electron_sources[j], ele_handle);
|
180 |
|
|
const std::vector<pat::Electron>& pat_electrons = *(ele_handle.product());
|
181 |
|
|
|
182 |
|
|
for (unsigned int i = 0; i < pat_electrons.size(); ++i) {
|
183 |
|
|
pat::Electron pat_ele = pat_electrons[i];
|
184 |
|
|
Electron ele;
|
185 |
|
|
|
186 |
|
|
ele.charge = pat_ele.charge();
|
187 |
|
|
ele.pt = pat_ele.pt();
|
188 |
|
|
ele.eta = pat_ele.eta();
|
189 |
|
|
ele.phi = pat_ele.phi();
|
190 |
|
|
ele.energy = pat_ele.energy();
|
191 |
|
|
ele.vertex_x = pat_ele.vertex().x();
|
192 |
|
|
ele.vertex_y = pat_ele.vertex().y();
|
193 |
|
|
ele.vertex_z = pat_ele.vertex().z();
|
194 |
|
|
ele.supercluster_eta = pat_ele.superCluster()->eta();
|
195 |
|
|
ele.supercluster_phi = pat_ele.superCluster()->phi();
|
196 |
|
|
ele.dB = pat_ele.dB();
|
197 |
|
|
//ele.particleIso = pat_ele.particleIso();
|
198 |
|
|
ele.neutralHadronIso = pat_ele.neutralHadronIso();
|
199 |
|
|
ele.chargedHadronIso = pat_ele.chargedHadronIso();
|
200 |
|
|
ele.trackIso = pat_ele.trackIso();
|
201 |
|
|
ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
|
202 |
|
|
|
203 |
|
|
eles[j].push_back(ele);
|
204 |
|
|
}
|
205 |
|
|
}
|
206 |
|
|
}
|
207 |
|
|
|
208 |
|
|
// ------------- muons -------------
|
209 |
|
|
if(doMuons){
|
210 |
|
|
for(size_t j=0; j< muon_sources.size(); ++j){
|
211 |
|
|
mus[j].clear();
|
212 |
|
|
|
213 |
|
|
edm::Handle< std::vector<pat::Muon> > mu_handle;
|
214 |
|
|
iEvent.getByLabel(muon_sources[j], mu_handle);
|
215 |
|
|
const std::vector<pat::Muon>& pat_muons = *(mu_handle.product());
|
216 |
|
|
|
217 |
|
|
for (unsigned int i = 0; i <pat_muons.size() ; ++i) {
|
218 |
|
|
pat::Muon pat_mu = pat_muons[i];
|
219 |
|
|
|
220 |
|
|
Muon mu;
|
221 |
|
|
mu.charge = pat_mu.charge();
|
222 |
|
|
mu.pt = pat_mu.pt();
|
223 |
|
|
mu.eta = pat_mu.eta();
|
224 |
|
|
mu.phi = pat_mu.phi();
|
225 |
|
|
mu.energy = pat_mu.energy();
|
226 |
|
|
mu.vertex_x = pat_mu.vertex().x();
|
227 |
|
|
mu.vertex_y = pat_mu.vertex().y();
|
228 |
|
|
mu.vertex_z = pat_mu.vertex().z();
|
229 |
|
|
mu.dB = pat_mu.dB();
|
230 |
|
|
//mu.particleIso = pat_mu.particleIso();
|
231 |
|
|
mu.neutralHadronIso = pat_mu.neutralHadronIso();
|
232 |
|
|
mu.chargedHadronIso = pat_mu.chargedHadronIso();
|
233 |
|
|
mu.trackIso = pat_mu.trackIso();
|
234 |
|
|
mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
|
235 |
|
|
mu.isGlobalMuon = pat_mu.isGlobalMuon();
|
236 |
|
|
mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
|
237 |
|
|
mu.isTrackerMuon = pat_mu.isTrackerMuon();
|
238 |
|
|
mu.numberOfMatchedStations = pat_mu.numberOfMatchedStations();
|
239 |
|
|
reco::TrackRef globalTrack = pat_mu.globalTrack();
|
240 |
|
|
if(!globalTrack.isNull()){
|
241 |
|
|
mu.globalTrack_chi2 = globalTrack->chi2();
|
242 |
|
|
mu.globalTrack_ndof = globalTrack->ndof();
|
243 |
|
|
mu.globalTrack_d0 = globalTrack->d0();
|
244 |
|
|
mu.globalTrack_d0Error = globalTrack->d0Error();
|
245 |
|
|
mu.globalTrack_numberOfValidHits = globalTrack->numberOfValidHits();
|
246 |
|
|
mu.globalTrack_numberOfLostHits = globalTrack->numberOfLostHits();
|
247 |
|
|
}
|
248 |
|
|
else{
|
249 |
|
|
mu.globalTrack_chi2 = 0;
|
250 |
|
|
mu.globalTrack_ndof = 0;
|
251 |
|
|
mu.globalTrack_d0 = 0;
|
252 |
|
|
mu.globalTrack_d0Error = 0;
|
253 |
|
|
mu.globalTrack_numberOfValidHits = 0;
|
254 |
|
|
mu.globalTrack_numberOfLostHits = 0;
|
255 |
|
|
}
|
256 |
|
|
reco::TrackRef innerTrack = pat_mu.innerTrack();
|
257 |
|
|
if(!innerTrack.isNull()){
|
258 |
|
|
mu.innerTrack_chi2 = innerTrack->chi2();
|
259 |
|
|
mu.innerTrack_ndof = innerTrack->ndof();
|
260 |
|
|
mu.innerTrack_d0 = innerTrack->d0();
|
261 |
|
|
mu.innerTrack_d0Error = innerTrack->d0Error();
|
262 |
|
|
mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
|
263 |
|
|
mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
|
264 |
|
|
}
|
265 |
|
|
else{
|
266 |
|
|
mu.innerTrack_chi2 = 0;
|
267 |
|
|
mu.innerTrack_ndof = 0;
|
268 |
|
|
mu.innerTrack_d0 = 0;
|
269 |
|
|
mu.innerTrack_d0Error = 0;
|
270 |
|
|
mu.innerTrack_numberOfValidHits = 0;
|
271 |
|
|
mu.innerTrack_numberOfLostHits = 0;
|
272 |
|
|
}
|
273 |
|
|
reco::TrackRef outerTrack = pat_mu.outerTrack();
|
274 |
|
|
if(!outerTrack.isNull()){
|
275 |
|
|
mu.outerTrack_chi2 = outerTrack->chi2();
|
276 |
|
|
mu.outerTrack_ndof = outerTrack->ndof();
|
277 |
|
|
mu.outerTrack_d0 = outerTrack->d0();
|
278 |
|
|
mu.outerTrack_d0Error = outerTrack->d0Error();
|
279 |
|
|
mu.outerTrack_numberOfValidHits = outerTrack->numberOfValidHits();
|
280 |
|
|
mu.outerTrack_numberOfLostHits = outerTrack->numberOfLostHits();
|
281 |
|
|
}
|
282 |
|
|
else{
|
283 |
|
|
mu.outerTrack_chi2 = 0;
|
284 |
|
|
mu.outerTrack_ndof = 0;
|
285 |
|
|
mu.outerTrack_d0 = 0;
|
286 |
|
|
mu.outerTrack_d0Error = 0;
|
287 |
|
|
mu.outerTrack_numberOfValidHits = 0;
|
288 |
|
|
mu.outerTrack_numberOfLostHits = 0;
|
289 |
|
|
}
|
290 |
|
|
|
291 |
|
|
mus[j].push_back(mu);
|
292 |
|
|
}
|
293 |
|
|
}
|
294 |
|
|
}
|
295 |
|
|
// ------------- taus -------------
|
296 |
|
|
|
297 |
|
|
if(doTaus){
|
298 |
|
|
for(size_t j=0; j< tau_sources.size(); ++j){
|
299 |
|
|
taus[j].clear();
|
300 |
|
|
|
301 |
|
|
|
302 |
|
|
edm::Handle< std::vector<pat::Tau> > tau_handle;
|
303 |
|
|
iEvent.getByLabel(tau_sources[j], tau_handle);
|
304 |
|
|
const std::vector<pat::Tau>& pat_taus = *(tau_handle.product());
|
305 |
|
|
|
306 |
|
|
for (unsigned int i = 0; i < pat_taus.size(); ++i) {
|
307 |
|
|
pat::Tau pat_tau = pat_taus[i];
|
308 |
|
|
if(pat_tau.pt() < tau_ptmin) continue;
|
309 |
|
|
if(fabs(pat_tau.eta()) > tau_etamax) continue;
|
310 |
|
|
|
311 |
|
|
Tau tau;
|
312 |
|
|
tau.charge = pat_tau.charge();
|
313 |
|
|
tau.pt = pat_tau.pt();
|
314 |
|
|
tau.eta = pat_tau.eta();
|
315 |
|
|
tau.phi = pat_tau.phi();
|
316 |
|
|
tau.energy = pat_tau.energy();
|
317 |
|
|
|
318 |
|
|
taus[j].push_back(tau);
|
319 |
|
|
}
|
320 |
|
|
}
|
321 |
|
|
}
|
322 |
|
|
|
323 |
|
|
// ------------- jets -------------
|
324 |
|
|
if(doJets){
|
325 |
|
|
for(size_t j=0; j< jet_sources.size(); ++j){
|
326 |
|
|
|
327 |
|
|
jets[j].clear();
|
328 |
|
|
|
329 |
|
|
edm::Handle< std::vector<pat::Jet> > jet_handle;
|
330 |
|
|
iEvent.getByLabel(jet_sources[j], jet_handle);
|
331 |
|
|
const std::vector<pat::Jet>& pat_jets = *(jet_handle.product());
|
332 |
|
|
|
333 |
|
|
for (unsigned int i = 0; i < pat_jets.size(); ++i) {
|
334 |
|
|
pat::Jet pat_jet = pat_jets[i];
|
335 |
|
|
if(pat_jet.pt() < jet_ptmin) continue;
|
336 |
|
|
if(fabs(pat_jet.eta()) > jet_etamax) continue;
|
337 |
|
|
// std::cout << "available btag discriminators: " << std::endl;
|
338 |
|
|
// const std::vector<std::pair<std::string, float> > & dis = pat_jets[i].getPairDiscri();
|
339 |
|
|
// for(size_t k=0; k<dis.size(); ++k){
|
340 |
|
|
// std::cout << dis[k].first << std::endl;
|
341 |
|
|
// }
|
342 |
|
|
|
343 |
|
|
Jet jet;
|
344 |
|
|
jet.charge = pat_jet.charge();
|
345 |
|
|
jet.pt = pat_jet.pt();
|
346 |
|
|
jet.eta = pat_jet.eta();
|
347 |
|
|
jet.phi = pat_jet.phi();
|
348 |
|
|
jet.energy = pat_jet.energy();
|
349 |
|
|
jet.numberOfDaughters =pat_jet.numberOfDaughters();
|
350 |
|
|
const reco::TrackRefVector& jettracks = pat_jet.associatedTracks();
|
351 |
|
|
jet.nTracks = jettracks.size();
|
352 |
|
|
jet.jetArea = pat_jet.jetArea();
|
353 |
|
|
jet.pileup = pat_jet.pileup();
|
354 |
|
|
if(pat_jet.isPFJet()){
|
355 |
|
|
jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
|
356 |
|
|
jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
|
357 |
|
|
jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
|
358 |
|
|
jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
|
359 |
|
|
jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
|
360 |
|
|
jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
|
361 |
|
|
jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
|
362 |
|
|
jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
|
363 |
|
|
jet.muonMultiplicity =pat_jet.muonMultiplicity();
|
364 |
|
|
jet.electronMultiplicity =pat_jet.electronMultiplicity();
|
365 |
|
|
jet.photonMultiplicity =pat_jet.photonMultiplicity();
|
366 |
|
|
}
|
367 |
|
|
jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
|
368 |
|
|
jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
|
369 |
|
|
jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
|
370 |
|
|
jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
|
371 |
|
|
jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
|
372 |
|
|
jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
|
373 |
|
|
|
374 |
|
|
jets[j].push_back(jet);
|
375 |
|
|
}
|
376 |
|
|
}
|
377 |
|
|
}
|
378 |
|
|
|
379 |
|
|
// ------------- top jets -------------
|
380 |
|
|
if(doTopJets){
|
381 |
|
|
for(size_t j=0; j< topjet_sources.size(); ++j){
|
382 |
|
|
|
383 |
|
|
topjets[j].clear();
|
384 |
|
|
|
385 |
|
|
edm::Handle<pat::JetCollection> pat_topjets;
|
386 |
|
|
//edm::Handle<std::vector<reco::Jet> > pat_topjets;
|
387 |
|
|
iEvent.getByLabel(topjet_sources[j], pat_topjets);
|
388 |
|
|
|
389 |
|
|
for (unsigned int i = 0; i < pat_topjets->size(); i++) {
|
390 |
|
|
const pat::Jet pat_topjet = * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
|
391 |
|
|
if(pat_topjet.pt() < topjet_ptmin) continue;
|
392 |
|
|
if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
|
393 |
|
|
|
394 |
|
|
TopJet topjet;
|
395 |
|
|
topjet.charge = pat_topjet.charge();
|
396 |
|
|
topjet.pt = pat_topjet.pt();
|
397 |
|
|
topjet.eta = pat_topjet.eta();
|
398 |
|
|
topjet.phi = pat_topjet.phi();
|
399 |
|
|
topjet.energy = pat_topjet.energy();
|
400 |
|
|
topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
|
401 |
|
|
const reco::TrackRefVector& topjettracks = pat_topjet.associatedTracks();
|
402 |
|
|
topjet.nTracks = topjettracks.size();
|
403 |
|
|
topjet.jetArea = pat_topjet.jetArea();
|
404 |
|
|
topjet.pileup = pat_topjet.pileup();
|
405 |
|
|
// topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
|
406 |
|
|
// topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
|
407 |
|
|
// topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
|
408 |
|
|
// topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
|
409 |
|
|
// topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
|
410 |
|
|
// topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
|
411 |
|
|
// topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
|
412 |
|
|
// topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
|
413 |
|
|
// topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
|
414 |
|
|
// topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
|
415 |
|
|
// topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
|
416 |
|
|
|
417 |
|
|
topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
|
418 |
|
|
topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
|
419 |
|
|
topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
|
420 |
|
|
topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
|
421 |
|
|
topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
|
422 |
|
|
topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
|
423 |
|
|
|
424 |
|
|
for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
|
425 |
|
|
Particle subjet_v4;
|
426 |
|
|
subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
|
427 |
|
|
subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
|
428 |
|
|
subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
|
429 |
|
|
subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
|
430 |
|
|
topjet.subjets.push_back(subjet_v4);
|
431 |
|
|
}
|
432 |
|
|
topjets[j].push_back(topjet);
|
433 |
|
|
}
|
434 |
|
|
}
|
435 |
|
|
}
|
436 |
|
|
|
437 |
|
|
// ------------- photons -------------
|
438 |
|
|
if(doPhotons){
|
439 |
|
|
for(size_t j=0; j< photon_sources.size(); ++j){
|
440 |
|
|
phs[j].clear();
|
441 |
|
|
|
442 |
|
|
edm::Handle< std::vector<pat::Photon> > photon_handle;
|
443 |
|
|
iEvent.getByLabel(photon_sources[j], photon_handle);
|
444 |
|
|
const std::vector<pat::Photon>& pat_photons = *(photon_handle.product());
|
445 |
|
|
|
446 |
|
|
for (unsigned int i = 0; i < pat_photons.size(); ++i) {
|
447 |
|
|
pat::Photon pat_photon = pat_photons[i];
|
448 |
|
|
Photon ph;
|
449 |
|
|
ph.charge = 0;
|
450 |
|
|
ph.pt = pat_photon.pt();
|
451 |
|
|
ph.eta = pat_photon.eta();
|
452 |
|
|
ph.phi = pat_photon.phi();
|
453 |
|
|
ph.energy = pat_photon.energy();
|
454 |
|
|
ph.vertex_x = pat_photon.vertex().x();
|
455 |
|
|
ph.vertex_y = pat_photon.vertex().y();
|
456 |
|
|
ph.vertex_z = pat_photon.vertex().z();
|
457 |
|
|
ph.supercluster_eta = pat_photon.superCluster()->eta();
|
458 |
|
|
ph.supercluster_phi = pat_photon.superCluster()->phi();
|
459 |
|
|
// ph.neutralHadronIso = pat_photon.neutralHadronIso();
|
460 |
|
|
// ph.chargedHadronIso = pat_photon.chargedHadronIso();
|
461 |
|
|
ph.trackIso = pat_photon.trackIso();
|
462 |
|
|
phs[j].push_back(ph);
|
463 |
|
|
}
|
464 |
|
|
}
|
465 |
|
|
}
|
466 |
|
|
|
467 |
|
|
// ------------- MET -------------
|
468 |
|
|
if(doMET){
|
469 |
|
|
for(size_t j=0; j< met_sources.size(); ++j){
|
470 |
|
|
|
471 |
|
|
edm::Handle< std::vector<pat::MET> > met_handle;
|
472 |
|
|
iEvent.getByLabel(met_sources[j], met_handle);
|
473 |
|
|
const std::vector<pat::MET>& pat_mets = *(met_handle.product());
|
474 |
|
|
|
475 |
|
|
if(pat_mets.size()!=1){
|
476 |
|
|
std::cout<< "WARNING: number of METs = " << pat_mets.size() <<", should be 1" << std::endl;
|
477 |
|
|
}
|
478 |
|
|
else{
|
479 |
|
|
pat::MET pat_met = pat_mets[0];
|
480 |
|
|
|
481 |
|
|
met[j].pt=pat_met.pt();
|
482 |
|
|
met[j].phi=pat_met.phi();
|
483 |
|
|
met[j].mEtSig= pat_met.mEtSig();
|
484 |
|
|
}
|
485 |
|
|
|
486 |
|
|
}
|
487 |
|
|
}
|
488 |
|
|
|
489 |
|
|
// ------------- trigger -------------
|
490 |
|
|
|
491 |
|
|
edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
|
492 |
|
|
edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
|
493 |
|
|
iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
|
494 |
|
|
|
495 |
|
|
const edm::Provenance *meta = dummy_TriggerEvent.provenance();
|
496 |
|
|
std::string nameProcess = meta->processName();
|
497 |
|
|
edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
|
498 |
|
|
triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
|
499 |
|
|
|
500 |
|
|
edm::Handle<edm::TriggerResults> trigger;
|
501 |
|
|
iEvent.getByLabel(triggerResultTag, trigger);
|
502 |
|
|
const edm::TriggerResults& trig = *(trigger.product());
|
503 |
|
|
|
504 |
|
|
triggerResults.clear();
|
505 |
|
|
triggerNames.clear();
|
506 |
|
|
L1_prescale.clear();
|
507 |
|
|
HLT_prescale.clear();
|
508 |
|
|
|
509 |
|
|
edm::Service<edm::service::TriggerNamesService> tns;
|
510 |
|
|
std::vector<std::string> triggerNames_all;
|
511 |
|
|
tns->getTrigPaths(trig,triggerNames_all);
|
512 |
|
|
|
513 |
|
|
if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
|
514 |
|
|
for(unsigned int i=0; i<trig.size(); ++i){
|
515 |
|
|
std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
|
516 |
|
|
for(; it!=trigger_prefixes.end(); ++it){
|
517 |
|
|
if(triggerNames_all[i].substr(0, it->size()) == *it)break;
|
518 |
|
|
}
|
519 |
|
|
if(it==trigger_prefixes.end()) continue;
|
520 |
|
|
|
521 |
|
|
//triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
|
522 |
|
|
triggerResults.push_back(trig.accept(i));
|
523 |
|
|
if(newrun) triggerNames.push_back(triggerNames_all[i]);
|
524 |
|
|
if(isRealData){
|
525 |
|
|
std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
|
526 |
|
|
L1_prescale.push_back(pre.first);
|
527 |
|
|
HLT_prescale.push_back(pre.second);
|
528 |
|
|
}
|
529 |
|
|
}
|
530 |
|
|
// for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
|
531 |
|
|
// std::cout << (*iter).first << " " << (*iter).second << std::endl;
|
532 |
|
|
// }
|
533 |
|
|
newrun=false;
|
534 |
|
|
|
535 |
|
|
// ------------- generator info -------------
|
536 |
|
|
|
537 |
|
|
if(doGenInfo){
|
538 |
|
|
genInfo.weights.clear();
|
539 |
|
|
genInfo.binningValues.clear();
|
540 |
|
|
genps.clear();
|
541 |
|
|
|
542 |
|
|
edm::Handle<GenEventInfoProduct> genEventInfoProduct;
|
543 |
|
|
iEvent.getByLabel("generator", genEventInfoProduct);
|
544 |
|
|
const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
|
545 |
|
|
|
546 |
|
|
genInfo.binningValues = genEventInfo.binningValues();
|
547 |
|
|
genInfo.weights = genEventInfo.weights();
|
548 |
|
|
genInfo.alphaQCD = genEventInfo.alphaQCD();
|
549 |
|
|
genInfo.alphaQED = genEventInfo.alphaQED();
|
550 |
|
|
genInfo.qScale = genEventInfo.qScale();
|
551 |
|
|
|
552 |
|
|
const gen::PdfInfo* pdf = genEventInfo.pdf();
|
553 |
|
|
if(pdf){
|
554 |
|
|
genInfo.pdf_id1=pdf->id.first;
|
555 |
|
|
genInfo.pdf_id2=pdf->id.second;
|
556 |
|
|
genInfo.pdf_x1=pdf->x.first;
|
557 |
|
|
genInfo.pdf_x2=pdf->x.second;
|
558 |
|
|
genInfo.pdf_scalePDF=pdf->scalePDF;
|
559 |
|
|
genInfo.pdf_xPDF1=pdf->xPDF.first;
|
560 |
|
|
genInfo.pdf_xPDF2=pdf->xPDF.second;
|
561 |
|
|
}
|
562 |
|
|
else{
|
563 |
|
|
genInfo.pdf_id1=-999;
|
564 |
|
|
genInfo.pdf_id2=-999;
|
565 |
|
|
genInfo.pdf_x1=-999;
|
566 |
|
|
genInfo.pdf_x2=-999;
|
567 |
|
|
genInfo.pdf_scalePDF=-999;
|
568 |
|
|
genInfo.pdf_xPDF1=-999;
|
569 |
|
|
genInfo.pdf_xPDF2=-999;
|
570 |
|
|
}
|
571 |
|
|
|
572 |
|
|
edm::Handle<std::vector<PileupSummaryInfo> > pus;
|
573 |
|
|
iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
|
574 |
|
|
genInfo.pileup_NumInteractions_intime=0;
|
575 |
|
|
genInfo.pileup_NumInteractions_ootbefore=0;
|
576 |
|
|
genInfo.pileup_NumInteractions_ootafter=0;
|
577 |
|
|
if(pus.isValid()){
|
578 |
|
|
genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
|
579 |
|
|
for(size_t i=0; i<pus->size(); ++i){
|
580 |
|
|
if(pus->at(i).getBunchCrossing() == 0) // intime pileup
|
581 |
|
|
genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
|
582 |
|
|
else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
|
583 |
|
|
genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
|
584 |
|
|
}
|
585 |
|
|
else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
|
586 |
|
|
genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
|
587 |
|
|
}
|
588 |
|
|
}
|
589 |
|
|
}
|
590 |
|
|
|
591 |
|
|
edm::Handle<reco::GenParticleCollection> genPartColl;
|
592 |
|
|
iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
|
593 |
|
|
int index=-1;
|
594 |
|
|
for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
|
595 |
|
|
index++;
|
596 |
|
|
|
597 |
|
|
//write out only top quarks and status 3 particles (works fine only for MadGraph)
|
598 |
|
|
if(abs(iter->pdgId())==6 || iter->status()==3){
|
599 |
|
|
GenParticle genp;
|
600 |
|
|
genp.charge = iter->charge();;
|
601 |
|
|
genp.pt = iter->p4().pt();
|
602 |
|
|
genp.eta = iter->p4().eta();
|
603 |
|
|
genp.phi = iter->p4().phi();
|
604 |
|
|
genp.energy = iter->p4().E();
|
605 |
|
|
genp.index =index;
|
606 |
|
|
genp.status = iter->status();
|
607 |
|
|
genp.pdgId = iter->pdgId();
|
608 |
|
|
|
609 |
|
|
genp.mother1=-1;
|
610 |
|
|
genp.mother2=-1;
|
611 |
|
|
genp.daughter1=-1;
|
612 |
|
|
genp.daughter2=-1;
|
613 |
|
|
|
614 |
|
|
int nm=iter->numberOfMothers();
|
615 |
|
|
int nd=iter->numberOfDaughters();
|
616 |
|
|
|
617 |
|
|
if (nm>0) genp.mother1 = iter->motherRef(0).key();
|
618 |
|
|
if (nm>1) genp.mother2 = iter->motherRef(nm-1).key();
|
619 |
|
|
if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
|
620 |
|
|
if (nd>1) genp.daughter2 = iter->daughterRef(nd-1).key();
|
621 |
|
|
|
622 |
|
|
genps.push_back(genp);
|
623 |
|
|
}
|
624 |
|
|
}
|
625 |
|
|
|
626 |
|
|
}
|
627 |
|
|
|
628 |
|
|
// ------------- primary vertices and beamspot -------------
|
629 |
|
|
|
630 |
|
|
if(doPV){
|
631 |
|
|
for(size_t j=0; j< pv_sources.size(); ++j){
|
632 |
|
|
pvs[j].clear();
|
633 |
|
|
|
634 |
|
|
edm::Handle< std::vector<reco::Vertex> > pv_handle;
|
635 |
|
|
iEvent.getByLabel(pv_sources[j], pv_handle);
|
636 |
|
|
const std::vector<reco::Vertex>& reco_pvs = *(pv_handle.product());
|
637 |
|
|
|
638 |
|
|
for (unsigned int i = 0; i < reco_pvs.size(); ++i) {
|
639 |
|
|
reco::Vertex reco_pv = reco_pvs[i];
|
640 |
|
|
|
641 |
|
|
PrimaryVertex pv;
|
642 |
|
|
pv.x = reco_pv.x();
|
643 |
|
|
pv.y = reco_pv.y();
|
644 |
|
|
pv.z = reco_pv.z();
|
645 |
|
|
pv.nTracks = reco_pv.nTracks();
|
646 |
|
|
//pv.isValid = reco_pv.isValid();
|
647 |
|
|
pv.chi2 = reco_pv.chi2();
|
648 |
|
|
pv.ndof = reco_pv.ndof();
|
649 |
|
|
|
650 |
|
|
pvs[j].push_back(pv);
|
651 |
|
|
}
|
652 |
|
|
}
|
653 |
|
|
}
|
654 |
|
|
|
655 |
|
|
edm::Handle<reco::BeamSpot> beamSpot;
|
656 |
|
|
iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
|
657 |
|
|
const reco::BeamSpot & bsp = *beamSpot;
|
658 |
|
|
|
659 |
|
|
beamspot_x0 = bsp.x0();
|
660 |
|
|
beamspot_y0 = bsp.y0();
|
661 |
|
|
beamspot_z0 = bsp.z0();
|
662 |
|
|
|
663 |
|
|
tr->Fill();
|
664 |
|
|
|
665 |
|
|
}
|
666 |
|
|
|
667 |
|
|
|
668 |
|
|
// ------------ method called once each job just before starting event loop ------------
|
669 |
|
|
void
|
670 |
|
|
NtupleWriter::beginJob()
|
671 |
|
|
{
|
672 |
|
|
}
|
673 |
|
|
|
674 |
|
|
// ------------ method called once each job just after ending the event loop ------------
|
675 |
|
|
void
|
676 |
|
|
NtupleWriter::endJob()
|
677 |
|
|
{
|
678 |
|
|
outfile->cd();
|
679 |
|
|
tr->Write();
|
680 |
|
|
outfile->Close();
|
681 |
|
|
}
|
682 |
|
|
|
683 |
|
|
// ------------ method called when starting to processes a run ------------
|
684 |
|
|
void
|
685 |
|
|
NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup)
|
686 |
|
|
{
|
687 |
|
|
bool setup_changed = false;
|
688 |
|
|
hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
|
689 |
|
|
newrun=true;
|
690 |
|
|
}
|
691 |
|
|
|
692 |
|
|
// ------------ method called when ending the processing of a run ------------
|
693 |
|
|
void
|
694 |
|
|
NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
|
695 |
|
|
{
|
696 |
|
|
}
|
697 |
|
|
|
698 |
|
|
// ------------ method called when starting to processes a luminosity block ------------
|
699 |
|
|
void
|
700 |
|
|
NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
|
701 |
|
|
{
|
702 |
|
|
}
|
703 |
|
|
|
704 |
|
|
// ------------ method called when ending the processing of a luminosity block ------------
|
705 |
|
|
void
|
706 |
|
|
NtupleWriter::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
|
707 |
|
|
{
|
708 |
|
|
}
|
709 |
|
|
|
710 |
|
|
// ------------ method fills 'descriptions' with the allowed parameters for the module ------------
|
711 |
|
|
void
|
712 |
|
|
NtupleWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
|
713 |
|
|
//The following says we do not know what parameters are allowed so do no validation
|
714 |
|
|
// Please change this to state exactly what you do use, even if it is no parameters
|
715 |
|
|
edm::ParameterSetDescription desc;
|
716 |
|
|
desc.setUnknown();
|
717 |
|
|
descriptions.addDefault(desc);
|
718 |
|
|
}
|
719 |
|
|
|
720 |
|
|
//define this as a plug-in
|
721 |
|
|
DEFINE_FWK_MODULE(NtupleWriter);
|