ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc
Revision: 1.1
Committed: Mon Apr 2 12:40:12 2012 UTC (13 years, 1 month ago) by peiffer
Content type: text/plain
Branch: MAIN
Branch point for: INITIAL
Log Message:
new

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