ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc
Revision: 1.26
Committed: Mon Jun 10 13:33:05 2013 UTC (11 years, 11 months ago) by peiffer
Content type: text/plain
Branch: MAIN
Changes since 1.25: +28 -11 lines
Log Message:
new tau variables

File Contents

# Content
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.25 2013/06/06 07:38:31 peiffer Exp $
17 //
18 //
19
20 #include "UHHAnalysis/NtupleWriter/interface/NtupleWriter.h"
21 #include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputer.h"
22 #include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputerWrapper.h"
23 #include "RecoBTau/JetTagComputer/interface/JetTagComputer.h"
24 #include "RecoBTau/JetTagComputer/interface/JetTagComputerRecord.h"
25 #include "RecoBTag/SecondaryVertex/interface/CombinedSVComputer.h"
26
27 //
28 // constants, enums and typedefs
29 //
30
31 //
32 // static data member definitions
33 //
34
35 //
36 // constructors and destructor
37 //
38 NtupleWriter::NtupleWriter(const edm::ParameterSet& iConfig)
39
40 {
41 //now do what ever initialization is needed
42 //edm::Service<TFileService> fs;
43 //tr = fs->make<TTree>("AnalysisTree","AnalysisTree");
44
45 fileName = iConfig.getParameter<std::string>("fileName");
46 outfile = new TFile(fileName,"RECREATE");
47 outfile->cd();
48 tr = new TTree("AnalysisTree","AnalysisTree");
49
50 std::string name;
51
52 doElectrons = iConfig.getParameter<bool>("doElectrons");
53 doMuons = iConfig.getParameter<bool>("doMuons");
54 doTaus = iConfig.getParameter<bool>("doTaus");
55 doJets = iConfig.getParameter<bool>("doJets");
56 doGenJets = iConfig.getParameter<bool>("doGenJets");
57 doJECUncertainty = iConfig.getParameter<bool>("doJECUncertainty");
58 doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");
59 doPhotons = iConfig.getParameter<bool>("doPhotons");
60 doMET = iConfig.getParameter<bool>("doMET");
61 doGenInfo = iConfig.getParameter<bool>("doGenInfo");
62 doAllGenParticles = iConfig.getParameter<bool>("doAllGenParticles");
63 doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
64 doPV = iConfig.getParameter<bool>("doPV");
65 doTopJets = iConfig.getParameter<bool>("doTopJets");
66 doTrigger = iConfig.getParameter<bool>("doTrigger");
67 SVComputer_ = iConfig.getUntrackedParameter<edm::InputTag>("svComputer",edm::InputTag("combinedSecondaryVertex"));
68 doTagInfos = iConfig.getUntrackedParameter<bool>("doTagInfos",false);
69 // initialization of tree variables
70
71 tr->Branch("run",&run);
72 tr->Branch("event",&event);
73 tr->Branch("luminosityBlock",&luminosityBlock);
74 tr->Branch("isRealData",&isRealData);
75 tr->Branch("rho",&rho);
76 rho_source = iConfig.getParameter<edm::InputTag>("rho_source");
77
78 //tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
79 if(doLumiInfo){
80 tr->Branch("intgRecLumi",&intgRecLumi);
81 tr->Branch("intgDelLumi",&intgDelLumi);
82 }
83 if(doPV){
84 tr->Branch("beamspot_x0",&beamspot_x0);
85 tr->Branch("beamspot_y0",&beamspot_y0);
86 tr->Branch("beamspot_z0",&beamspot_z0);
87 }
88 if(doElectrons){
89 electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
90 for(size_t j=0; j< electron_sources.size(); ++j){
91 tr->Branch( electron_sources[j].c_str(), "std::vector<Electron>", &eles[j]);
92 }
93 }
94 if(doMuons){
95 muon_sources = iConfig.getParameter<std::vector<std::string> >("muon_sources");
96 for(size_t j=0; j< muon_sources.size(); ++j){
97 tr->Branch( muon_sources[j].c_str(), "std::vector<Muon>", &mus[j]);
98 }
99 }
100 if(doTaus){
101 tau_sources = iConfig.getParameter<std::vector<std::string> >("tau_sources");
102 tau_ptmin = iConfig.getParameter<double> ("tau_ptmin");
103 tau_etamax = iConfig.getParameter<double> ("tau_etamax");
104 for(size_t j=0; j< tau_sources.size(); ++j){
105 tr->Branch( tau_sources[j].c_str(), "std::vector<Tau>", &taus[j]);
106 }
107 }
108 if(doJets){
109 jet_sources = iConfig.getParameter<std::vector<std::string> >("jet_sources");
110 jet_ptmin = iConfig.getParameter<double> ("jet_ptmin");
111 jet_etamax = iConfig.getParameter<double> ("jet_etamax");
112 for(size_t j=0; j< jet_sources.size(); ++j){
113 tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
114 }
115 }
116 if(doGenJets){
117 genjet_sources = iConfig.getParameter<std::vector<std::string> >("genjet_sources");
118 genjet_ptmin = iConfig.getParameter<double> ("genjet_ptmin");
119 genjet_etamax = iConfig.getParameter<double> ("genjet_etamax");
120 for(size_t j=0; j< genjet_sources.size(); ++j){
121 tr->Branch( genjet_sources[j].c_str(), "std::vector<Particle>", &genjets[j]);
122 }
123 }
124 if(doTopJets){
125 topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
126 topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
127 topjet_etamax = iConfig.getParameter<double> ("topjet_etamax");
128 for(size_t j=0; j< topjet_sources.size(); ++j){
129 tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
130 }
131 }
132 if(doGenTopJets){
133 gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
134 gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
135 gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
136 for(size_t j=0; j< gentopjet_sources.size(); ++j){
137 tr->Branch( gentopjet_sources[j].c_str(), "std::vector<TopJet>", &gentopjets[j]);
138 }
139 }
140 if(doPhotons){
141 photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
142 for(size_t j=0; j< photon_sources.size(); ++j){
143 tr->Branch( photon_sources[j].c_str(), "std::vector<Photon>", &phs[j]);
144 }
145 }
146 if(doMET){
147 met_sources = iConfig.getParameter<std::vector<std::string> >("met_sources");
148 for(size_t j=0; j< met_sources.size(); ++j){
149 tr->Branch(met_sources[j].c_str(), "MET", &met[j]);
150 }
151 }
152 if(doPV){
153 pv_sources = iConfig.getParameter<std::vector<std::string> >("pv_sources");
154 for(size_t j=0; j< pv_sources.size(); ++j){
155 tr->Branch( pv_sources[j].c_str(), "std::vector<PrimaryVertex>", &pvs[j]);
156 }
157 }
158 if(doGenInfo){
159 genparticle_source= iConfig.getParameter<edm::InputTag>("genparticle_source");
160 tr->Branch("genInfo","GenInfo",&genInfo);
161 tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
162 }
163 if(doTrigger){
164 trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
165 //tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
166 tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
167 tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
168 //tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
169 //tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
170 }
171 newrun = true;
172 }
173
174
175 NtupleWriter::~NtupleWriter()
176 {
177
178 // do anything here that needs to be done at desctruction time
179 // (e.g. close files, deallocate resources etc.)
180
181 }
182
183
184 //
185 // member functions
186 //
187
188 // ------------ method called for each event ------------
189 void
190 NtupleWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
191 {
192 // using namespace edm;
193
194
195 // ------------- common variables -------------
196
197 run = iEvent.id().run();
198 event = iEvent.id().event();
199 luminosityBlock = iEvent.luminosityBlock();
200 isRealData = iEvent.isRealData();
201
202 edm::Handle<double> m_rho;
203 iEvent.getByLabel(rho_source,m_rho);
204 rho=*m_rho;
205
206 // if(isRealData){
207 // edm::Handle<bool> bool_handle;
208 // iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
209 // HBHENoiseFilterResult = *(bool_handle.product());
210 // }
211 // else HBHENoiseFilterResult = false;
212
213 // ------------- primary vertices and beamspot -------------
214
215 if(doPV){
216 for(size_t j=0; j< pv_sources.size(); ++j){
217 pvs[j].clear();
218
219 edm::Handle< std::vector<reco::Vertex> > pv_handle;
220 iEvent.getByLabel(pv_sources[j], pv_handle);
221 const std::vector<reco::Vertex>& reco_pvs = *(pv_handle.product());
222
223 for (unsigned int i = 0; i < reco_pvs.size(); ++i) {
224 reco::Vertex reco_pv = reco_pvs[i];
225
226 PrimaryVertex pv;
227 pv.set_x( reco_pv.x());
228 pv.set_y( reco_pv.y());
229 pv.set_z( reco_pv.z());
230 pv.set_nTracks( reco_pv.nTracks());
231 //pv.set_isValid( reco_pv.isValid());
232 pv.set_chi2( reco_pv.chi2());
233 pv.set_ndof( reco_pv.ndof());
234
235 pvs[j].push_back(pv);
236 }
237 }
238
239 edm::Handle<reco::BeamSpot> beamSpot;
240 iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
241 const reco::BeamSpot & bsp = *beamSpot;
242
243 beamspot_x0 = bsp.x0();
244 beamspot_y0 = bsp.y0();
245 beamspot_z0 = bsp.z0();
246 }
247
248 // ------------- generator info -------------
249
250 if(doGenInfo){
251 genInfo.clear_weights();
252 genInfo.clear_binningValues();
253 genps.clear();
254
255 edm::Handle<GenEventInfoProduct> genEventInfoProduct;
256 iEvent.getByLabel("generator", genEventInfoProduct);
257 const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
258
259 for(unsigned int k=0; k<genEventInfo.binningValues().size();++k){
260 genInfo.add_binningValue(genEventInfo.binningValues().at(k));
261 }
262 for(unsigned int k=0; k<genEventInfo.weights().size();++k){
263 genInfo.add_weight(genEventInfo.weights().at(k));
264 }
265 genInfo.set_alphaQCD(genEventInfo.alphaQCD());
266 genInfo.set_alphaQED(genEventInfo.alphaQED());
267 genInfo.set_qScale(genEventInfo.qScale());
268
269 const gen::PdfInfo* pdf = genEventInfo.pdf();
270 if(pdf){
271 genInfo.set_pdf_id1(pdf->id.first);
272 genInfo.set_pdf_id2(pdf->id.second);
273 genInfo.set_pdf_x1(pdf->x.first);
274 genInfo.set_pdf_x2(pdf->x.second);
275 genInfo.set_pdf_scalePDF(pdf->scalePDF);
276 genInfo.set_pdf_xPDF1(pdf->xPDF.first);
277 genInfo.set_pdf_xPDF2(pdf->xPDF.second);
278 }
279 else{
280 genInfo.set_pdf_id1(-999);
281 genInfo.set_pdf_id2(-999);
282 genInfo.set_pdf_x1(-999);
283 genInfo.set_pdf_x2(-999);
284 genInfo.set_pdf_scalePDF(-999);
285 genInfo.set_pdf_xPDF1(-999);
286 genInfo.set_pdf_xPDF2(-999);
287 }
288
289 edm::Handle<std::vector<PileupSummaryInfo> > pus;
290 iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
291 genInfo.set_pileup_NumInteractions_intime(0);
292 genInfo.set_pileup_NumInteractions_ootbefore(0);
293 genInfo.set_pileup_NumInteractions_ootafter(0);
294 if(pus.isValid()){
295 genInfo.set_pileup_TrueNumInteractions ( (float) pus->at(0).getTrueNumInteractions());
296 for(size_t i=0; i<pus->size(); ++i){
297 if(pus->at(i).getBunchCrossing() == 0) // intime pileup
298 genInfo.set_pileup_NumInteractions_intime( genInfo.pileup_NumInteractions_intime() + pus->at(i).getPU_NumInteractions());
299 else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
300 genInfo.set_pileup_NumInteractions_ootbefore( genInfo.pileup_NumInteractions_ootbefore() + pus->at(i).getPU_NumInteractions());
301 }
302 else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
303 genInfo.set_pileup_NumInteractions_ootafter( genInfo.pileup_NumInteractions_ootafter() + pus->at(i).getPU_NumInteractions());
304 }
305 }
306 }
307
308 edm::Handle<reco::GenParticleCollection> genPartColl;
309 iEvent.getByLabel(genparticle_source, genPartColl);
310 int index=-1;
311 for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
312 index++;
313
314 //write out only top quarks,final state leptons and status 3 particles (works fine only for MadGraph)
315 bool islepton = iter->status()==1 && abs(iter->pdgId())>=11 && abs(iter->pdgId())<=16 ;
316 if(abs(iter->pdgId())==6 || iter->status()==3 || islepton || doAllGenParticles){
317 GenParticle genp;
318 genp.set_charge(iter->charge());
319 genp.set_pt(iter->p4().pt());
320 genp.set_eta(iter->p4().eta());
321 genp.set_phi(iter->p4().phi());
322 genp.set_energy(iter->p4().E());
323 genp.set_index(index);
324 genp.set_status( iter->status());
325 genp.set_pdgId( iter->pdgId());
326
327 genp.set_mother1(-1);
328 genp.set_mother2(-1);
329 genp.set_daughter1(-1);
330 genp.set_daughter2(-1);
331
332 int nm=iter->numberOfMothers();
333 int nd=iter->numberOfDaughters();
334
335
336 if (nm>0) genp.set_mother1( iter->motherRef(0).key());
337 if (nm>1) genp.set_mother2( iter->motherRef(1).key());
338 if (nd>0) genp.set_daughter1( iter->daughterRef(0).key());
339 if (nd>1) genp.set_daughter2( iter->daughterRef(1).key());
340
341 genps.push_back(genp);
342 }
343 }
344 }
345
346 // ------------- electrons -------------
347 if(doElectrons){
348
349 // edm::Handle<reco::ConversionCollection> hConversions;
350 // iEvent.getByLabel("allConversions", hConversions);
351
352 // edm::Handle<reco::BeamSpot> beamSpot;
353 // iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
354 // const reco::BeamSpot & bsp = *beamSpot;
355
356 for(size_t j=0; j< electron_sources.size(); ++j){
357 eles[j].clear();
358 edm::Handle< std::vector<pat::Electron> > ele_handle;
359 iEvent.getByLabel(electron_sources[j], ele_handle);
360 const std::vector<pat::Electron>& pat_electrons = *(ele_handle.product());
361
362 for (unsigned int i = 0; i < pat_electrons.size(); ++i) {
363 pat::Electron pat_ele = pat_electrons[i];
364 Electron ele;
365
366 ele.set_charge( pat_ele.charge());
367 ele.set_pt( pat_ele.pt());
368 ele.set_eta( pat_ele.eta());
369 ele.set_phi( pat_ele.phi());
370 ele.set_energy( pat_ele.energy());
371 ele.set_vertex_x(pat_ele.vertex().x());
372 ele.set_vertex_y(pat_ele.vertex().y());
373 ele.set_vertex_z(pat_ele.vertex().z());
374 ele.set_supercluster_eta(pat_ele.superCluster()->eta());
375 ele.set_supercluster_phi(pat_ele.superCluster()->phi());
376 ele.set_dB(pat_ele.dB());
377 //ele.set_particleIso(pat_ele.particleIso());
378 ele.set_neutralHadronIso(pat_ele.neutralHadronIso());
379 ele.set_chargedHadronIso(pat_ele.chargedHadronIso());
380 ele.set_trackIso(pat_ele.trackIso());
381 ele.set_photonIso(pat_ele.photonIso());
382 ele.set_puChargedHadronIso(pat_ele.puChargedHadronIso());
383 ele.set_gsfTrack_trackerExpectedHitsInner_numberOfLostHits(pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits());
384 ele.set_gsfTrack_px( pat_ele.gsfTrack()->px());
385 ele.set_gsfTrack_py( pat_ele.gsfTrack()->py());
386 ele.set_gsfTrack_pz( pat_ele.gsfTrack()->pz());
387 ele.set_gsfTrack_vx( pat_ele.gsfTrack()->vx());
388 ele.set_gsfTrack_vy( pat_ele.gsfTrack()->vy());
389 ele.set_gsfTrack_vz( pat_ele.gsfTrack()->vz());
390 //ele.set_passconversionveto(!ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position()));
391 ele.set_passconversionveto(pat_ele.passConversionVeto());
392 ele.set_dEtaIn(pat_ele.deltaEtaSuperClusterTrackAtVtx());
393 ele.set_dPhiIn(pat_ele.deltaPhiSuperClusterTrackAtVtx());
394 ele.set_sigmaIEtaIEta(pat_ele.sigmaIetaIeta());
395 ele.set_HoverE(pat_ele.hadronicOverEm());
396 ele.set_fbrem(pat_ele.fbrem());
397 ele.set_EoverPIn(pat_ele.eSuperClusterOverP());
398 ele.set_EcalEnergy(pat_ele.ecalEnergy());
399 ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
400 ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
401 float AEff03 = 0.00;
402 if(isRealData){
403 AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2011);
404 }else{
405 AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAFall11MC);
406 }
407 ele.set_AEff(AEff03);
408
409 eles[j].push_back(ele);
410 }
411 }
412 }
413
414 // ------------- muons -------------
415 if(doMuons){
416 for(size_t j=0; j< muon_sources.size(); ++j){
417 mus[j].clear();
418
419 edm::Handle< std::vector<pat::Muon> > mu_handle;
420 iEvent.getByLabel(muon_sources[j], mu_handle);
421 const std::vector<pat::Muon>& pat_muons = *(mu_handle.product());
422
423 for (unsigned int i = 0; i <pat_muons.size() ; ++i) {
424 pat::Muon pat_mu = pat_muons[i];
425
426 Muon mu;
427 mu.set_charge( pat_mu.charge());
428 mu.set_pt( pat_mu.pt());
429 mu.set_eta( pat_mu.eta());
430 mu.set_phi( pat_mu.phi());
431 mu.set_energy( pat_mu.energy());
432 mu.set_vertex_x ( pat_mu.vertex().x());
433 mu.set_vertex_y ( pat_mu.vertex().y());
434 mu.set_vertex_z ( pat_mu.vertex().z());
435 mu.set_dB ( pat_mu.dB());
436 //mu.particleIso ( pat_mu.particleIso());
437 mu.set_neutralHadronIso ( pat_mu.neutralHadronIso());
438 mu.set_chargedHadronIso ( pat_mu.chargedHadronIso());
439 mu.set_trackIso ( pat_mu.trackIso());
440 mu.set_photonIso ( pat_mu.photonIso());
441 mu.set_puChargedHadronIso ( pat_mu.puChargedHadronIso());
442 mu.set_isGlobalMuon ( pat_mu.isGlobalMuon());
443 mu.set_isPFMuon ( pat_mu.isPFMuon());
444 mu.set_isStandAloneMuon ( pat_mu.isStandAloneMuon());
445 mu.set_isTrackerMuon ( pat_mu.isTrackerMuon());
446 mu.set_numberOfMatchedStations ( pat_mu.numberOfMatchedStations());
447 reco::TrackRef globalTrack = pat_mu.globalTrack();
448 if(!globalTrack.isNull()){
449 mu.set_globalTrack_chi2 ( globalTrack->chi2());
450 mu.set_globalTrack_ndof ( globalTrack->ndof());
451 mu.set_globalTrack_d0 ( globalTrack->d0());
452 mu.set_globalTrack_d0Error ( globalTrack->d0Error());
453 mu.set_globalTrack_numberOfValidHits ( globalTrack->numberOfValidHits());
454 mu.set_globalTrack_numberOfLostHits ( globalTrack->numberOfLostHits());
455 mu.set_globalTrack_numberOfValidMuonHits(globalTrack->hitPattern().numberOfValidMuonHits() );
456 }
457 else{
458 mu.set_globalTrack_chi2 ( 0);
459 mu.set_globalTrack_ndof ( 0);
460 mu.set_globalTrack_d0 ( 0);
461 mu.set_globalTrack_d0Error ( 0);
462 mu.set_globalTrack_numberOfValidHits ( 0);
463 mu.set_globalTrack_numberOfLostHits ( 0);
464 }
465 reco::TrackRef innerTrack = pat_mu.innerTrack();
466 if(!innerTrack.isNull()){
467 mu.set_innerTrack_chi2 ( innerTrack->chi2());
468 mu.set_innerTrack_ndof ( innerTrack->ndof());
469 mu.set_innerTrack_d0 ( innerTrack->d0());
470 mu.set_innerTrack_d0Error ( innerTrack->d0Error());
471 mu.set_innerTrack_numberOfValidHits ( innerTrack->numberOfValidHits());
472 mu.set_innerTrack_numberOfLostHits ( innerTrack->numberOfLostHits());
473 mu.set_innerTrack_trackerLayersWithMeasurement ( innerTrack->hitPattern().trackerLayersWithMeasurement());
474 mu.set_innerTrack_numberOfValidPixelHits ( innerTrack->hitPattern().numberOfValidPixelHits());
475 }
476 else{
477 mu.set_innerTrack_chi2 ( 0);
478 mu.set_innerTrack_ndof ( 0);
479 mu.set_innerTrack_d0 ( 0);
480 mu.set_innerTrack_d0Error ( 0);
481 mu.set_innerTrack_numberOfValidHits ( 0);
482 mu.set_innerTrack_numberOfLostHits ( 0);
483 mu.set_innerTrack_trackerLayersWithMeasurement ( 0);
484 mu.set_innerTrack_numberOfValidPixelHits ( 0);
485 }
486 reco::TrackRef outerTrack = pat_mu.outerTrack();
487 if(!outerTrack.isNull()){
488 mu.set_outerTrack_chi2 ( outerTrack->chi2());
489 mu.set_outerTrack_ndof ( outerTrack->ndof());
490 mu.set_outerTrack_d0 ( outerTrack->d0());
491 mu.set_outerTrack_d0Error ( outerTrack->d0Error());
492 mu.set_outerTrack_numberOfValidHits ( outerTrack->numberOfValidHits());
493 mu.set_outerTrack_numberOfLostHits ( outerTrack->numberOfLostHits());
494 }
495 else{
496 mu.set_outerTrack_chi2 ( 0);
497 mu.set_outerTrack_ndof ( 0);
498 mu.set_outerTrack_d0 ( 0);
499 mu.set_outerTrack_d0Error ( 0);
500 mu.set_outerTrack_numberOfValidHits ( 0);
501 mu.set_outerTrack_numberOfLostHits ( 0);
502 }
503
504 mus[j].push_back(mu);
505 }
506 }
507 }
508 // ------------- taus -------------
509
510 if(doTaus){
511 for(size_t j=0; j< tau_sources.size(); ++j){
512 taus[j].clear();
513
514
515 edm::Handle< std::vector<pat::Tau> > tau_handle;
516 iEvent.getByLabel(tau_sources[j], tau_handle);
517 const std::vector<pat::Tau>& pat_taus = *(tau_handle.product());
518
519 for (unsigned int i = 0; i < pat_taus.size(); ++i) {
520 pat::Tau pat_tau = pat_taus[i];
521 if(pat_tau.pt() < tau_ptmin) continue;
522 if(fabs(pat_tau.eta()) > tau_etamax) continue;
523
524 Tau tau;
525 tau.set_charge( pat_tau.charge());
526 tau.set_pt( pat_tau.pt());
527 tau.set_eta( pat_tau.eta());
528 tau.set_phi( pat_tau.phi());
529 tau.set_energy( pat_tau.energy());
530 tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
531 //tau.set_byVLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
532 tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
533 tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
534 tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
535 tau.set_byLooseIsolationMVA( pat_tau.tauID("byLooseIsolationMVA")>0.5);
536 tau.set_byMediumIsolationMVA( pat_tau.tauID("byMediumIsolationMVA")>0.5);
537 tau.set_byTightIsolationMVA( pat_tau.tauID("byTightIsolationMVA")>0.5);
538 tau.set_byLooseIsolationMVA2( pat_tau.tauID("byLooseIsolationMVA2")>0.5);
539 tau.set_byMediumIsolationMVA2( pat_tau.tauID("byMediumIsolationMVA2")>0.5);
540 tau.set_byTightIsolationMVA2( pat_tau.tauID("byTightIsolationMVA2")>0.5);
541 tau.set_byLooseCombinedIsolationDeltaBetaCorr3Hits( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits")>0.5);
542 tau.set_byMediumCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits")>0.5);
543 tau.set_byTightCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr3Hits")>0.5);
544 tau.set_againstElectronLooseMVA3 ( pat_tau.tauID("againstElectronLooseMVA3")>0.5);
545 tau.set_againstElectronMediumMVA3 ( pat_tau.tauID("againstElectronMediumMVA3")>0.5);
546 tau.set_againstElectronTightMVA3 ( pat_tau.tauID("againstElectronTightMVA3")>0.5);
547 tau.set_againstElectronVTightMVA3 ( pat_tau.tauID("againstElectronVTightMVA3")>0.5);
548 tau.set_againstMuonLoose2 ( pat_tau.tauID("againstMuonLoose2")>0.5);
549 tau.set_againstMuonMedium2 ( pat_tau.tauID("againstMuonMedium2")>0.5);
550 tau.set_againstMuonTight2 ( pat_tau.tauID("againstMuonTight2")>0.5);
551 tau.set_byIsolationMVAraw( pat_tau.tauID("byIsolationMVAraw"));
552 tau.set_byIsolationMVA2raw( pat_tau.tauID("byIsolationMVA2raw"));
553 tau.set_decayMode( pat_tau.decayMode() );
554 tau.set_byCombinedIsolationDeltaBetaCorrRaw( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw"));
555 tau.set_byCombinedIsolationDeltaBetaCorrRaw3Hits( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));
556
557 // std::cout << pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") << std::endl;
558
559 // reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
560 // if(!leadPFCand.isNull()){
561 // tau.set_leadPFCand_px ( leadPFCand->px());
562 // tau.set_leadPFCand_py ( leadPFCand->py());
563 // tau.set_leadPFCand_pz ( leadPFCand->pz());
564 // }
565 // else{
566 // tau.set_leadPFCand_px ( 0);
567 // tau.set_leadPFCand_py ( 0);
568 // tau.set_leadPFCand_pz ( 0);
569 // }
570 taus[j].push_back(tau);
571 }
572 }
573 }
574
575 //-------------- gen jets -------------
576
577 if(doGenJets){
578 for(size_t j=0; j< genjet_sources.size(); ++j){
579
580 genjets[j].clear();
581
582 edm::Handle< std::vector<reco::GenJet> > genjet_handle;
583 iEvent.getByLabel(genjet_sources[j], genjet_handle);
584 const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
585
586 for (unsigned int i = 0; i < gen_jets.size(); ++i) {
587 const reco::GenJet* gen_jet = &gen_jets[i];
588 if(gen_jet->pt() < genjet_ptmin) continue;
589 if(fabs(gen_jet->eta()) > genjet_etamax) continue;
590
591 Particle jet;
592 jet.set_charge(gen_jet->charge());
593 jet.set_pt(gen_jet->pt());
594 jet.set_eta(gen_jet->eta());
595 jet.set_phi(gen_jet->phi());
596 jet.set_energy(gen_jet->energy());
597
598 genjets[j].push_back(jet);
599
600 }
601 }
602 }
603
604 // ------------- jets -------------
605 if(doJets){
606 for(size_t j=0; j< jet_sources.size(); ++j){
607
608 jets[j].clear();
609
610 edm::Handle< std::vector<pat::Jet> > jet_handle;
611 iEvent.getByLabel(jet_sources[j], jet_handle);
612 const std::vector<pat::Jet>& pat_jets = *(jet_handle.product());
613
614 for (unsigned int i = 0; i < pat_jets.size(); ++i) {
615 pat::Jet pat_jet = pat_jets[i];
616 if(pat_jet.pt() < jet_ptmin) continue;
617 if(fabs(pat_jet.eta()) > jet_etamax) continue;
618 // std::cout << "available btag discriminators: " << std::endl;
619 // const std::vector<std::pair<std::string, float> > & dis = pat_jets[i].getPairDiscri();
620 // for(size_t k=0; k<dis.size(); ++k){
621 // std::cout << dis[k].first << std::endl;
622 // }
623
624 Jet jet;
625 jet.set_charge(pat_jet.charge());
626 jet.set_pt(pat_jet.pt());
627 jet.set_eta(pat_jet.eta());
628 jet.set_phi(pat_jet.phi());
629 jet.set_energy(pat_jet.energy());
630 jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
631 const reco::TrackRefVector& jettracks = pat_jet.associatedTracks();
632 jet.set_nTracks ( jettracks.size());
633 jet.set_jetArea(pat_jet.jetArea());
634 jet.set_pileup(pat_jet.pileup());
635 if(pat_jet.isPFJet()){
636 jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
637 jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
638 jet.set_chargedEmEnergyFraction (pat_jet.chargedEmEnergyFraction());
639 jet.set_chargedHadronEnergyFraction (pat_jet.chargedHadronEnergyFraction());
640 jet.set_muonEnergyFraction (pat_jet.muonEnergyFraction());
641 jet.set_photonEnergyFraction (pat_jet.photonEnergyFraction());
642 jet.set_chargedMultiplicity (pat_jet.chargedMultiplicity());
643 jet.set_neutralMultiplicity (pat_jet.neutralMultiplicity());
644 jet.set_muonMultiplicity (pat_jet.muonMultiplicity());
645 jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
646 jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
647 }
648 if(doJECUncertainty){
649 jecUnc->setJetEta(pat_jet.eta());
650 jecUnc->setJetPt(pat_jet.pt());
651 jet.set_JEC_uncertainty(jecUnc->getUncertainty(true));
652 }
653 jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
654
655 jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
656 jet.set_btag_simpleSecondaryVertexHighPur(pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
657 jet.set_btag_combinedSecondaryVertex(pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags"));
658 jet.set_btag_combinedSecondaryVertexMVA(pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
659 jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
660 jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
661
662
663 const reco::GenJet *genj = pat_jet.genJet();
664 if(genj){
665
666 for(unsigned int k=0; k<genjets->size(); ++k){
667 if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta())
668 jet.set_genjet_index(k);
669 }
670 // if( jet.genjet_index()<0){
671 // std::cout<< "genjet not found for " << genj->pt() << " " << genj->eta() << std::endl;
672 // }
673
674 if(doAllGenParticles){
675 std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
676 for(unsigned int l = 0; l<jetgenps.size(); ++l){
677 for(unsigned int k=0; k< genps.size(); ++k){
678 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
679 jet.add_genparticles_index(genps[k].index());
680 break;
681 }
682 }
683 }
684 if(jet.genparticles_indices().size()!= jetgenps.size())
685 std::cout << "WARNING: Found only " << jet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
686 }
687
688 }
689
690 jets[j].push_back(jet);
691 }
692 }
693 }
694
695
696
697 // ------------- top jets -------------
698 if(doTopJets){
699 for(size_t j=0; j< topjet_sources.size(); ++j){
700
701 topjets[j].clear();
702
703 edm::Handle<pat::JetCollection> pat_topjets;
704 //edm::Handle<std::vector<reco::Jet> > pat_topjets;
705 iEvent.getByLabel(topjet_sources[j], pat_topjets);
706
707 for (unsigned int i = 0; i < pat_topjets->size(); i++) {
708
709 const pat::Jet pat_topjet = * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
710 if(pat_topjet.pt() < topjet_ptmin) continue;
711 if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
712
713 TopJet topjet;
714 topjet.set_charge(pat_topjet.charge());
715 topjet.set_pt(pat_topjet.pt());
716 topjet.set_eta(pat_topjet.eta());
717 topjet.set_phi(pat_topjet.phi());
718 topjet.set_energy(pat_topjet.energy());
719 topjet.set_numberOfDaughters(pat_topjet.numberOfDaughters());
720 const reco::TrackRefVector& topjettracks = pat_topjet.associatedTracks();
721 topjet.set_nTracks( topjettracks.size());
722 topjet.set_jetArea( pat_topjet.jetArea());
723 topjet.set_pileup( pat_topjet.pileup());
724 // topjet.set_neutralEmEnergyFraction(pat_topjet.neutralEmEnergyFraction());
725 // topjet.set_neutralHadronEnergyFraction(pat_topjet.neutralHadronEnergyFraction());
726 // topjet.set_chargedEmEnergyFraction(pat_topjet.chargedEmEnergyFraction());
727 // topjet.set_chargedHadronEnergyFraction(pat_topjet.chargedHadronEnergyFraction());
728 // topjet.set_muonEnergyFraction(pat_topjet.muonEnergyFraction());
729 // topjet.set_photonEnergyFraction(pat_topjet.photonEnergyFraction());
730 // topjet.set_chargedMultiplicity(pat_topjet.chargedMultiplicity());
731 // topjet.set_neutralMultiplicity(pat_topjet.neutralMultiplicity());
732 // topjet.set_muonMultiplicity(pat_topjet.muonMultiplicity());
733 // topjet.set_electronMultiplicity(pat_topjet.electronMultiplicity());
734 // topjet.set_photonMultiplicity(pat_topjet.photonMultiplicity());
735 if(doJECUncertainty){
736 jecUnc->setJetEta(pat_topjet.eta());
737 jecUnc->setJetPt(pat_topjet.pt());
738 topjet.set_JEC_uncertainty( jecUnc->getUncertainty(true));
739 }
740 topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
741
742 topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
743 topjet.set_btag_simpleSecondaryVertexHighPur(pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
744 topjet.set_btag_combinedSecondaryVertex(pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags"));
745 topjet.set_btag_combinedSecondaryVertexMVA(pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
746 topjet.set_btag_jetBProbability(pat_topjet.bDiscriminator("jetBProbabilityBJetTags"));
747 topjet.set_btag_jetProbability(pat_topjet.bDiscriminator("jetProbabilityBJetTags"));
748
749 /*
750 const reco::GenJet *genj = pat_topjet.genJet();
751 if(genj){
752 topjet.set_genjet_pt ( genj->pt());
753 topjet.set_genjet_eta ( genj->eta());
754 topjet.set_genjet_phi ( genj->phi());
755 topjet.set_genjet_energy ( genj->energy());
756 if(doAllGenParticles){
757 std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
758 for(unsigned int l = 0; l<jetgenps.size(); ++l){
759 for(unsigned int k=0; k< genps.size(); ++k){
760 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
761 topjet.add_genparticles_index(genps[k].index());
762 break;
763 }
764 }
765 }
766 if(topjet.genparticles_indices().size()!= jetgenps.size())
767 std::cout << "WARNING: Found only " << topjet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
768 }
769 }
770 */
771
772 for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
773 Particle subjet_v4;
774
775 reco::Candidate const * subjetd = pat_topjet.daughter(k);
776 pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
777 if(patsubjetd)
778 {
779 subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
780 subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
781 subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
782 subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
783 topjet.add_subjet(subjet_v4);
784 //btag info
785 topjet.add_subFlavour(patsubjetd->partonFlavour());
786 topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
787 if (doTagInfos)
788 {
789 //ip tag info
790 reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
791 topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
792 topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
793 topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
794 topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
795 topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
796 topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
797 topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
798 topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
799 topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
800 topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
801 topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
802 topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));
803 topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
804 topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
805 topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
806 topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
807 topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
808 topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
809 topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
810 topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
811 topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
812 //sv tag info
813 reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
814 topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
815 topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
816 topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
817 topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
818 topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
819 topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
820 topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
821 std::vector<TLorentzVector> vp4; vp4.clear();
822 std::vector<float> vchi2; vchi2.clear();
823 std::vector<float> vndof; vndof.clear();
824 std::vector<float> vchi2ndof; vchi2ndof.clear();
825 std::vector<float> vtrsize; vtrsize.clear();
826 for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
827 {
828 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
829 vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
830 vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());
831 vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());
832 vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());
833 vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());
834 }
835 topjet.add_subSecondaryVertex(vp4);
836 topjet.add_subVertexChi2(vchi2);
837 topjet.add_subVertexNdof(vndof);
838 topjet.add_subVertexNormalizedChi2(vchi2ndof);
839 topjet.add_subVertexTracksSize(vtrsize);
840 //try computer
841 edm::ESHandle<JetTagComputer> computerHandle;
842 iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
843 const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
844 if(computer)
845 {
846 computer->passEventSetup(iSetup);
847 std::vector<const reco::BaseTagInfo*> baseTagInfos;
848 baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
849 baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );
850 JetTagComputer::TagInfoHelper helper(baseTagInfos);
851 reco::TaggingVariableList vars = computer->taggingVariables(helper);
852 topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
853 topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
854 topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
855 topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
856 }
857 }
858 }
859 else
860 {
861 //filling only standard information in case the subjet has not been pat-tified during the pattuples production
862 subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
863 subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
864 subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
865 subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
866 topjet.add_subjet(subjet_v4);
867 }
868
869
870 }
871 topjets[j].push_back(topjet);
872 }
873 }
874 }
875
876
877 // ------------- generator top jets -------------
878 if(doGenTopJets){
879 for(size_t j=0; j< gentopjet_sources.size(); ++j){
880
881 gentopjets[j].clear();
882
883 edm::Handle<reco::BasicJetCollection> reco_gentopjets;
884 //edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
885 iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
886
887 for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
888
889 const reco::BasicJet reco_gentopjet = reco_gentopjets->at(i);
890 if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
891 if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
892
893 TopJet gentopjet;
894 gentopjet.set_charge(reco_gentopjet.charge());
895 gentopjet.set_pt(reco_gentopjet.pt());
896 gentopjet.set_eta(reco_gentopjet.eta());
897 gentopjet.set_phi(reco_gentopjet.phi());
898 gentopjet.set_energy(reco_gentopjet.energy());
899 gentopjet.set_numberOfDaughters(reco_gentopjet.numberOfDaughters());
900
901 for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
902 Particle subjet_v4;
903 subjet_v4.set_pt(reco_gentopjet.daughter(k)->p4().pt());
904 subjet_v4.set_eta(reco_gentopjet.daughter(k)->p4().eta());
905 subjet_v4.set_phi(reco_gentopjet.daughter(k)->p4().phi());
906 subjet_v4.set_energy(reco_gentopjet.daughter(k)->p4().E());
907 gentopjet.add_subjet(subjet_v4);
908 }
909 gentopjets[j].push_back(gentopjet);
910 }
911 }
912 }
913
914 // ------------- photons -------------
915 if(doPhotons){
916 for(size_t j=0; j< photon_sources.size(); ++j){
917 phs[j].clear();
918
919 edm::Handle< std::vector<pat::Photon> > photon_handle;
920 iEvent.getByLabel(photon_sources[j], photon_handle);
921 const std::vector<pat::Photon>& pat_photons = *(photon_handle.product());
922
923 for (unsigned int i = 0; i < pat_photons.size(); ++i) {
924 pat::Photon pat_photon = pat_photons[i];
925 Photon ph;
926 ph.set_charge(0);
927 ph.set_pt( pat_photon.pt());
928 ph.set_eta( pat_photon.eta());
929 ph.set_phi( pat_photon.phi());
930 ph.set_energy( pat_photon.energy());
931 ph.set_vertex_x(pat_photon.vertex().x());
932 ph.set_vertex_y(pat_photon.vertex().y());
933 ph.set_vertex_z(pat_photon.vertex().z());
934 ph.set_supercluster_eta(pat_photon.superCluster()->eta());
935 ph.set_supercluster_phi(pat_photon.superCluster()->phi());
936 // ph.set_neutralHadronIso(pat_photon.neutralHadronIso());
937 // ph.set_chargedHadronIso(pat_photon.chargedHadronIso());
938 ph.set_trackIso(pat_photon.trackIso());
939 phs[j].push_back(ph);
940 }
941 }
942 }
943
944 // ------------- MET -------------
945 if(doMET){
946 for(size_t j=0; j< met_sources.size(); ++j){
947
948 edm::Handle< std::vector<pat::MET> > met_handle;
949 iEvent.getByLabel(met_sources[j], met_handle);
950 const std::vector<pat::MET>& pat_mets = *(met_handle.product());
951
952 if(pat_mets.size()!=1){
953 std::cout<< "WARNING: number of METs = " << pat_mets.size() <<", should be 1" << std::endl;
954 }
955 else{
956 pat::MET pat_met = pat_mets[0];
957
958 met[j].set_pt(pat_met.pt());
959 met[j].set_phi(pat_met.phi());
960 met[j].set_mEtSig(pat_met.mEtSig());
961 }
962
963 }
964 }
965
966 // ------------- trigger -------------
967 if(doTrigger){
968 edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
969 edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
970 iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
971
972 const edm::Provenance *meta = dummy_TriggerEvent.provenance();
973 std::string nameProcess = meta->processName();
974 edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
975 triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
976
977 edm::Handle<edm::TriggerResults> trigger;
978 iEvent.getByLabel(triggerResultTag, trigger);
979 const edm::TriggerResults& trig = *(trigger.product());
980
981 triggerResults.clear();
982 triggerNames.clear();
983 // L1_prescale.clear();
984 // HLT_prescale.clear();
985
986 edm::Service<edm::service::TriggerNamesService> tns;
987 std::vector<std::string> triggerNames_all;
988 tns->getTrigPaths(trig,triggerNames_all);
989
990 if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
991 for(unsigned int i=0; i<trig.size(); ++i){
992 std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
993 for(; it!=trigger_prefixes.end(); ++it){
994 if(triggerNames_all[i].substr(0, it->size()) == *it)break;
995 }
996 if(it==trigger_prefixes.end()) continue;
997
998 //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
999 triggerResults.push_back(trig.accept(i));
1000 if(newrun) triggerNames.push_back(triggerNames_all[i]);
1001 // if(isRealData){
1002 // std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
1003 // L1_prescale.push_back(pre.first);
1004 // HLT_prescale.push_back(pre.second);
1005 // //std::cout << triggerNames_all[i] << " " << pre.first << " " <<pre.second << " " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
1006 // }
1007 }
1008 // for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
1009 // std::cout << (*iter).first << " " << (*iter).second << std::endl;
1010 // }
1011 newrun=false;
1012 }
1013
1014
1015 tr->Fill();
1016 if(doLumiInfo)
1017 previouslumiblockwasfilled=true;
1018 }
1019
1020
1021 // ------------ method called once each job just before starting event loop ------------
1022 void
1023 NtupleWriter::beginJob()
1024 {
1025 if(doLumiInfo){
1026 totalRecLumi=0;
1027 totalDelLumi=0;
1028 previouslumiblockwasfilled=false;
1029 }
1030 }
1031
1032 // ------------ method called once each job just after ending the event loop ------------
1033 void
1034 NtupleWriter::endJob()
1035 {
1036 outfile->cd();
1037 tr->Write();
1038 outfile->Close();
1039 }
1040
1041 // ------------ method called when starting to processes a run ------------
1042 void
1043 NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup)
1044 {
1045 if(doTrigger){
1046 //bool setup_changed = false;
1047 //hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1048 newrun=true;
1049 }
1050
1051 if(doJets || doTopJets){
1052 if(doJECUncertainty){
1053 edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
1054 iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
1055 JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
1056 jecUnc = new JetCorrectionUncertainty(JetCorPar);
1057 }
1058 }
1059 }
1060
1061 // ------------ method called when ending the processing of a run ------------
1062 void
1063 NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
1064 {
1065 if(doLumiInfo)
1066 std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
1067 }
1068
1069 // ------------ method called when starting to processes a luminosity block ------------
1070 void
1071 NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
1072 {
1073 if(doLumiInfo){
1074 edm::Handle<LumiSummary> l;
1075 lumi.getByLabel("lumiProducer", l);
1076
1077 //add lumi of lumi blocks without any event to next lumiblock
1078 if(previouslumiblockwasfilled){
1079 intgRecLumi=0;
1080 intgDelLumi=0;
1081 }
1082 previouslumiblockwasfilled=false;
1083
1084 if (l.isValid()){;
1085 intgRecLumi+=l->intgRecLumi()*6.37;
1086 intgDelLumi+=l->intgDelLumi()*6.37;
1087 totalRecLumi+=l->intgRecLumi()*6.37;
1088 totalDelLumi+=l->intgDelLumi()*6.37;
1089 }
1090 //std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<" " << l->intgDelLumi()*6.37<<std::endl;
1091 //std::cout << "summed: "<< intgRecLumi << " " << intgDelLumi << std::endl;
1092 }
1093 }
1094
1095 // ------------ method called when ending the processing of a luminosity block ------------
1096 void
1097 NtupleWriter::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
1098 {
1099 }
1100
1101 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
1102 void
1103 NtupleWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1104 //The following says we do not know what parameters are allowed so do no validation
1105 // Please change this to state exactly what you do use, even if it is no parameters
1106 edm::ParameterSetDescription desc;
1107 desc.setUnknown();
1108 descriptions.addDefault(desc);
1109 }
1110
1111 //define this as a plug-in
1112 DEFINE_FWK_MODULE(NtupleWriter);