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