ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc
Revision: 1.13
Committed: Thu Apr 19 14:47:13 2012 UTC (13 years ago) by peiffer
Content type: text/plain
Branch: MAIN
Changes since 1.12: +3 -1 lines
Log Message:
isolation

File Contents

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