ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc
(Generate patch)

Comparing UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc (file contents):
Revision 1.6 by peiffer, Wed Apr 11 15:15:44 2012 UTC vs.
Revision 1.24 by peiffer, Wed Jun 5 14:13:58 2013 UTC

# Line 17 | Line 17
17   //
18   //
19  
20 < #include "UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
21 <
22 <
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
# Line 50 | Line 53 | NtupleWriter::NtupleWriter(const edm::Pa
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 <
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("HBHENoiseFilterResult",&HBHENoiseFilterResult);
76 <  tr->Branch("beamspot_x0",&beamspot_x0);
68 <  tr->Branch("beamspot_y0",&beamspot_y0);
69 <  tr->Branch("beamspot_z0",&beamspot_z0);
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){  
# Line 96 | Line 113 | NtupleWriter::NtupleWriter(const edm::Pa
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");
# Line 104 | Line 129 | NtupleWriter::NtupleWriter(const edm::Pa
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){  
# Line 123 | Line 156 | NtupleWriter::NtupleWriter(const edm::Pa
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    }
# Line 131 | Line 165 | NtupleWriter::NtupleWriter(const edm::Pa
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);
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   }
# Line 165 | Line 199 | NtupleWriter::analyze(const edm::Event&
199     luminosityBlock = iEvent.luminosityBlock();
200     isRealData      = iEvent.isRealData();
201  
202 <   if(isRealData){
203 <     edm::Handle<bool> bool_handle;
204 <     iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
205 <     HBHENoiseFilterResult = *(bool_handle.product());
206 <   }
207 <   else HBHENoiseFilterResult = false;
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  
# Line 186 | Line 224 | NtupleWriter::analyze(const edm::Event&
224           reco::Vertex reco_pv = reco_pvs[i];
225  
226           PrimaryVertex pv;
227 <         pv.x =  reco_pv.x();
228 <         pv.y =  reco_pv.y();
229 <         pv.z =  reco_pv.z();
230 <         pv.nTracks =  reco_pv.nTracks();
231 <         //pv.isValid =  reco_pv.isValid();
232 <         pv.chi2 =  reco_pv.chi2();
233 <         pv.ndof =  reco_pv.ndof();      
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 <   edm::Handle<reco::BeamSpot> beamSpot;
251 <   iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
252 <   const reco::BeamSpot & bsp = *beamSpot;
253 <  
254 <   beamspot_x0 = bsp.x0();
255 <   beamspot_y0 = bsp.y0();
256 <   beamspot_z0 = bsp.z0();
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 and status 3 particles (works fine only for MadGraph)
315 >       if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
316 >         GenParticle genp;
317 >         genp.set_charge(iter->charge());
318 >         genp.set_pt(iter->p4().pt());
319 >         genp.set_eta(iter->p4().eta());
320 >         genp.set_phi(iter->p4().phi());
321 >         genp.set_energy(iter->p4().E());
322 >         genp.set_index(index);
323 >         genp.set_status( iter->status());
324 >         genp.set_pdgId( iter->pdgId());
325 >
326 >         genp.set_mother1(-1);
327 >         genp.set_mother2(-1);
328 >         genp.set_daughter1(-1);
329 >         genp.set_daughter2(-1);
330 >        
331 >         int nm=iter->numberOfMothers();
332 >         int nd=iter->numberOfDaughters();
333 >
334 >        
335 >         if (nm>0) genp.set_mother1( iter->motherRef(0).key());
336 >         if (nm>1) genp.set_mother2( iter->motherRef(1).key());
337 >         if (nd>0) genp.set_daughter1( iter->daughterRef(0).key());
338 >         if (nd>1) genp.set_daughter2( iter->daughterRef(1).key());
339 >
340 >         genps.push_back(genp);
341 >       }
342 >     }
343 >   }
344  
345     // ------------- electrons -------------  
346     if(doElectrons){
347 +
348 + //      edm::Handle<reco::ConversionCollection> hConversions;
349 + //      iEvent.getByLabel("allConversions", hConversions);
350 +
351 + //      edm::Handle<reco::BeamSpot> beamSpot;
352 + //      iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
353 + //      const reco::BeamSpot & bsp = *beamSpot;
354 +
355       for(size_t j=0; j< electron_sources.size(); ++j){
356         eles[j].clear();
357         edm::Handle< std::vector<pat::Electron> > ele_handle;
# Line 219 | Line 362 | NtupleWriter::analyze(const edm::Event&
362           pat::Electron pat_ele = pat_electrons[i];
363           Electron ele;
364          
365 <         ele.charge =  pat_ele.charge();
366 <         ele.pt =  pat_ele.pt();
367 <         ele.eta =  pat_ele.eta();
368 <         ele.phi =  pat_ele.phi();
369 <         ele.energy =  pat_ele.energy();
370 <         ele.vertex_x = pat_ele.vertex().x();
371 <         ele.vertex_y = pat_ele.vertex().y();
372 <         ele.vertex_z = pat_ele.vertex().z();
373 <         ele.supercluster_eta = pat_ele.superCluster()->eta();
374 <         ele.supercluster_phi = pat_ele.superCluster()->phi();
375 <         ele.dB = pat_ele.dB();
376 <         //ele.particleIso = pat_ele.particleIso();
377 <         ele.neutralHadronIso = pat_ele.neutralHadronIso();
378 <         ele.chargedHadronIso = pat_ele.chargedHadronIso();
379 <         ele.trackIso = pat_ele.trackIso();
380 <         ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
381 <         ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
382 <         ele.gsfTrack_px= pat_ele.gsfTrack()->px();
383 <         ele.gsfTrack_py= pat_ele.gsfTrack()->py();
384 <         ele.gsfTrack_pz= pat_ele.gsfTrack()->pz();
385 <         ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
386 <         ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
387 <         ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
365 >         ele.set_charge( pat_ele.charge());
366 >         ele.set_pt( pat_ele.pt());
367 >         ele.set_eta( pat_ele.eta());
368 >         ele.set_phi( pat_ele.phi());
369 >         ele.set_energy( pat_ele.energy());
370 >         ele.set_vertex_x(pat_ele.vertex().x());
371 >         ele.set_vertex_y(pat_ele.vertex().y());
372 >         ele.set_vertex_z(pat_ele.vertex().z());
373 >         ele.set_supercluster_eta(pat_ele.superCluster()->eta());
374 >         ele.set_supercluster_phi(pat_ele.superCluster()->phi());
375 >         ele.set_dB(pat_ele.dB());
376 >         //ele.set_particleIso(pat_ele.particleIso());
377 >         ele.set_neutralHadronIso(pat_ele.neutralHadronIso());
378 >         ele.set_chargedHadronIso(pat_ele.chargedHadronIso());
379 >         ele.set_trackIso(pat_ele.trackIso());
380 >         ele.set_photonIso(pat_ele.photonIso());
381 >         ele.set_puChargedHadronIso(pat_ele.puChargedHadronIso());
382 >         ele.set_gsfTrack_trackerExpectedHitsInner_numberOfLostHits(pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits());
383 >         ele.set_gsfTrack_px( pat_ele.gsfTrack()->px());
384 >         ele.set_gsfTrack_py( pat_ele.gsfTrack()->py());
385 >         ele.set_gsfTrack_pz( pat_ele.gsfTrack()->pz());
386 >         ele.set_gsfTrack_vx( pat_ele.gsfTrack()->vx());
387 >         ele.set_gsfTrack_vy( pat_ele.gsfTrack()->vy());
388 >         ele.set_gsfTrack_vz( pat_ele.gsfTrack()->vz());
389 >         //ele.set_passconversionveto(!ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position()));
390 >         ele.set_passconversionveto(pat_ele.passConversionVeto());
391 >         ele.set_dEtaIn(pat_ele.deltaEtaSuperClusterTrackAtVtx());
392 >         ele.set_dPhiIn(pat_ele.deltaPhiSuperClusterTrackAtVtx());
393 >         ele.set_sigmaIEtaIEta(pat_ele.sigmaIetaIeta());
394 >         ele.set_HoverE(pat_ele.hadronicOverEm());
395 >         ele.set_fbrem(pat_ele.fbrem());
396 >         ele.set_EoverPIn(pat_ele.eSuperClusterOverP());
397 >         ele.set_EcalEnergy(pat_ele.ecalEnergy());
398 >         ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
399 >         ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
400 >         float AEff03 = 0.00;
401 >         if(isRealData){
402 >           AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2011);
403 >         }else{
404 >           AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAFall11MC);
405 >         }
406 >         ele.set_AEff(AEff03);
407 >
408           eles[j].push_back(ele);
409         }
410       }
# Line 260 | Line 423 | NtupleWriter::analyze(const edm::Event&
423           pat::Muon pat_mu = pat_muons[i];
424  
425           Muon mu;
426 <         mu.charge =  pat_mu.charge();
427 <         mu.pt =  pat_mu.pt();
428 <         mu.eta =  pat_mu.eta();
429 <         mu.phi =  pat_mu.phi();
430 <         mu.energy =  pat_mu.energy();
431 <         mu.vertex_x = pat_mu.vertex().x();
432 <         mu.vertex_y = pat_mu.vertex().y();
433 <         mu.vertex_z = pat_mu.vertex().z();
434 <         mu.dB = pat_mu.dB();
435 <         //mu.particleIso = pat_mu.particleIso();
436 <         mu.neutralHadronIso = pat_mu.neutralHadronIso();
437 <         mu.chargedHadronIso = pat_mu.chargedHadronIso();
438 <         mu.trackIso = pat_mu.trackIso();
439 <         mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
440 <         mu.isGlobalMuon = pat_mu.isGlobalMuon();
441 <         mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
442 <         mu.isTrackerMuon = pat_mu.isTrackerMuon();
443 <         mu.numberOfMatchedStations = pat_mu.numberOfMatchedStations();
426 >         mu.set_charge( pat_mu.charge());
427 >         mu.set_pt( pat_mu.pt());
428 >         mu.set_eta( pat_mu.eta());
429 >         mu.set_phi( pat_mu.phi());
430 >         mu.set_energy( pat_mu.energy());
431 >         mu.set_vertex_x ( pat_mu.vertex().x());
432 >         mu.set_vertex_y ( pat_mu.vertex().y());
433 >         mu.set_vertex_z ( pat_mu.vertex().z());
434 >         mu.set_dB ( pat_mu.dB());
435 >         //mu.particleIso ( pat_mu.particleIso());
436 >         mu.set_neutralHadronIso ( pat_mu.neutralHadronIso());
437 >         mu.set_chargedHadronIso ( pat_mu.chargedHadronIso());
438 >         mu.set_trackIso ( pat_mu.trackIso());
439 >         mu.set_photonIso ( pat_mu.photonIso());
440 >         mu.set_puChargedHadronIso ( pat_mu.puChargedHadronIso());
441 >         mu.set_isGlobalMuon ( pat_mu.isGlobalMuon());
442 >         mu.set_isPFMuon ( pat_mu.isPFMuon());
443 >         mu.set_isStandAloneMuon ( pat_mu.isStandAloneMuon());
444 >         mu.set_isTrackerMuon ( pat_mu.isTrackerMuon());
445 >         mu.set_numberOfMatchedStations ( pat_mu.numberOfMatchedStations());
446           reco::TrackRef globalTrack = pat_mu.globalTrack();
447           if(!globalTrack.isNull()){
448 <           mu.globalTrack_chi2 = globalTrack->chi2();
449 <           mu.globalTrack_ndof = globalTrack->ndof();
450 <           mu.globalTrack_d0 = globalTrack->d0();        
451 <           mu.globalTrack_d0Error = globalTrack->d0Error();
452 <           mu.globalTrack_numberOfValidHits = globalTrack->numberOfValidHits();
453 <           mu.globalTrack_numberOfLostHits = globalTrack->numberOfLostHits();
448 >           mu.set_globalTrack_chi2 ( globalTrack->chi2());
449 >           mu.set_globalTrack_ndof ( globalTrack->ndof());
450 >           mu.set_globalTrack_d0 ( globalTrack->d0());  
451 >           mu.set_globalTrack_d0Error ( globalTrack->d0Error());
452 >           mu.set_globalTrack_numberOfValidHits ( globalTrack->numberOfValidHits());
453 >           mu.set_globalTrack_numberOfLostHits ( globalTrack->numberOfLostHits());
454 >           mu.set_globalTrack_numberOfValidMuonHits(globalTrack->hitPattern().numberOfValidMuonHits() );
455           }
456           else{
457 <           mu.globalTrack_chi2 = 0;
458 <           mu.globalTrack_ndof = 0;
459 <           mu.globalTrack_d0 = 0;
460 <           mu.globalTrack_d0Error = 0;
461 <           mu.globalTrack_numberOfValidHits = 0;
462 <           mu.globalTrack_numberOfLostHits = 0;
457 >           mu.set_globalTrack_chi2 ( 0);
458 >           mu.set_globalTrack_ndof ( 0);
459 >           mu.set_globalTrack_d0 ( 0);
460 >           mu.set_globalTrack_d0Error ( 0);
461 >           mu.set_globalTrack_numberOfValidHits ( 0);
462 >           mu.set_globalTrack_numberOfLostHits ( 0);
463           }
464           reco::TrackRef innerTrack = pat_mu.innerTrack();
465           if(!innerTrack.isNull()){
466 <           mu.innerTrack_chi2 = innerTrack->chi2();
467 <           mu.innerTrack_ndof = innerTrack->ndof();
468 <           mu.innerTrack_d0 = innerTrack->d0();  
469 <           mu.innerTrack_d0Error = innerTrack->d0Error();
470 <           mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
471 <           mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
466 >           mu.set_innerTrack_chi2 ( innerTrack->chi2());
467 >           mu.set_innerTrack_ndof ( innerTrack->ndof());
468 >           mu.set_innerTrack_d0 ( innerTrack->d0());    
469 >           mu.set_innerTrack_d0Error ( innerTrack->d0Error());
470 >           mu.set_innerTrack_numberOfValidHits ( innerTrack->numberOfValidHits());
471 >           mu.set_innerTrack_numberOfLostHits ( innerTrack->numberOfLostHits());
472 >           mu.set_innerTrack_trackerLayersWithMeasurement ( innerTrack->hitPattern().trackerLayersWithMeasurement());
473 >           mu.set_innerTrack_numberOfValidPixelHits ( innerTrack->hitPattern().numberOfValidPixelHits());
474           }
475           else{
476 <           mu.innerTrack_chi2 = 0;
477 <           mu.innerTrack_ndof = 0;
478 <           mu.innerTrack_d0 = 0;
479 <           mu.innerTrack_d0Error = 0;
480 <           mu.innerTrack_numberOfValidHits = 0;
481 <           mu.innerTrack_numberOfLostHits = 0;
476 >           mu.set_innerTrack_chi2 ( 0);
477 >           mu.set_innerTrack_ndof ( 0);
478 >           mu.set_innerTrack_d0 ( 0);
479 >           mu.set_innerTrack_d0Error ( 0);
480 >           mu.set_innerTrack_numberOfValidHits ( 0);
481 >           mu.set_innerTrack_numberOfLostHits ( 0);
482 >           mu.set_innerTrack_trackerLayersWithMeasurement ( 0);
483 >           mu.set_innerTrack_numberOfValidPixelHits ( 0);
484           }
485           reco::TrackRef outerTrack = pat_mu.outerTrack();
486           if(!outerTrack.isNull()){
487 <           mu.outerTrack_chi2 = outerTrack->chi2();
488 <           mu.outerTrack_ndof = outerTrack->ndof();
489 <           mu.outerTrack_d0 = outerTrack->d0();  
490 <           mu.outerTrack_d0Error = outerTrack->d0Error();
491 <           mu.outerTrack_numberOfValidHits = outerTrack->numberOfValidHits();
492 <           mu.outerTrack_numberOfLostHits = outerTrack->numberOfLostHits();
487 >           mu.set_outerTrack_chi2 ( outerTrack->chi2());
488 >           mu.set_outerTrack_ndof ( outerTrack->ndof());
489 >           mu.set_outerTrack_d0 ( outerTrack->d0());    
490 >           mu.set_outerTrack_d0Error ( outerTrack->d0Error());
491 >           mu.set_outerTrack_numberOfValidHits ( outerTrack->numberOfValidHits());
492 >           mu.set_outerTrack_numberOfLostHits ( outerTrack->numberOfLostHits());
493           }
494           else{
495 <           mu.outerTrack_chi2 = 0;
496 <           mu.outerTrack_ndof = 0;
497 <           mu.outerTrack_d0 = 0;
498 <           mu.outerTrack_d0Error = 0;
499 <           mu.outerTrack_numberOfValidHits = 0;
500 <           mu.outerTrack_numberOfLostHits = 0;
495 >           mu.set_outerTrack_chi2 ( 0);
496 >           mu.set_outerTrack_ndof ( 0);
497 >           mu.set_outerTrack_d0 ( 0);
498 >           mu.set_outerTrack_d0Error ( 0);
499 >           mu.set_outerTrack_numberOfValidHits ( 0);
500 >           mu.set_outerTrack_numberOfLostHits ( 0);
501           }
502  
503           mus[j].push_back(mu);
# Line 351 | Line 521 | NtupleWriter::analyze(const edm::Event&
521           if(fabs(pat_tau.eta()) > tau_etamax) continue;
522  
523           Tau tau;
524 <         tau.charge =  pat_tau.charge();
525 <         tau.pt =  pat_tau.pt();
526 <         tau.eta =  pat_tau.eta();
527 <         tau.phi =  pat_tau.phi();
528 <         tau.energy =  pat_tau.energy();
529 <
524 >         tau.set_charge( pat_tau.charge());
525 >         tau.set_pt( pat_tau.pt());
526 >         tau.set_eta( pat_tau.eta());
527 >         tau.set_phi( pat_tau.phi());
528 >         tau.set_energy( pat_tau.energy());
529 >         tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
530 >         tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
531 >         tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
532 >         tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
533 >         tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
534 >         tau.set_againstElectronLoose  ( pat_tau.tauID("againstElectronLoose")>0.5);
535 >         tau.set_againstElectronMedium ( pat_tau.tauID("againstElectronMedium")>0.5);
536 >         tau.set_againstElectronTight ( pat_tau.tauID("againstElectronTight")>0.5);
537 >         tau.set_againstElectronMVA  ( pat_tau.tauID("againstElectronMVA")>0.5);
538 >         tau.set_againstMuonLoose ( pat_tau.tauID("againstMuonLoose")>0.5);
539 >         tau.set_againstMuonMedium ( pat_tau.tauID("againstMuonMedium")>0.5);
540 >         tau.set_againstMuonTight ( pat_tau.tauID("againstMuonTight")>0.5);
541 >
542 > //       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
543 > //       if(!leadPFCand.isNull()){
544 > //         tau.set_leadPFCand_px ( leadPFCand->px());
545 > //         tau.set_leadPFCand_py ( leadPFCand->py());
546 > //         tau.set_leadPFCand_pz ( leadPFCand->pz());
547 > //       }
548 > //       else{
549 > //         tau.set_leadPFCand_px ( 0);
550 > //         tau.set_leadPFCand_py ( 0);
551 > //         tau.set_leadPFCand_pz ( 0);
552 > //       }
553           taus[j].push_back(tau);
554         }
555       }
# Line 383 | Line 576 | NtupleWriter::analyze(const edm::Event&
576   //       }
577  
578           Jet jet;
579 <         jet.charge = pat_jet.charge();
580 <         jet.pt = pat_jet.pt();
581 <         jet.eta = pat_jet.eta();
582 <         jet.phi = pat_jet.phi();
583 <         jet.energy = pat_jet.energy();
584 <         jet.numberOfDaughters =pat_jet.numberOfDaughters();
579 >         jet.set_charge(pat_jet.charge());
580 >         jet.set_pt(pat_jet.pt());
581 >         jet.set_eta(pat_jet.eta());
582 >         jet.set_phi(pat_jet.phi());
583 >         jet.set_energy(pat_jet.energy());
584 >         jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
585           const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
586 <         jet.nTracks = jettracks.size();
587 <         jet.jetArea = pat_jet.jetArea();
588 <         jet.pileup = pat_jet.pileup();
586 >         jet.set_nTracks ( jettracks.size());
587 >         jet.set_jetArea(pat_jet.jetArea());
588 >         jet.set_pileup(pat_jet.pileup());
589           if(pat_jet.isPFJet()){
590 <           jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
591 <           jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
592 <           jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
593 <           jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
594 <           jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
595 <           jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
596 <           jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
597 <           jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
598 <           jet.muonMultiplicity =pat_jet.muonMultiplicity();
599 <           jet.electronMultiplicity =pat_jet.electronMultiplicity();
600 <           jet.photonMultiplicity =pat_jet.photonMultiplicity();
601 <         }
602 <
603 <         jecUnc->setJetEta(pat_jet.eta());
604 <         jecUnc->setJetPt(pat_jet.pt());
605 <         jet.JEC_uncertainty = jecUnc->getUncertainty(true);
606 <
607 <         jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
608 <         jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
609 <         jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
610 <         jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
611 <         jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
612 <         jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
590 >           jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
591 >           jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
592 >           jet.set_chargedEmEnergyFraction (pat_jet.chargedEmEnergyFraction());
593 >           jet.set_chargedHadronEnergyFraction (pat_jet.chargedHadronEnergyFraction());
594 >           jet.set_muonEnergyFraction (pat_jet.muonEnergyFraction());
595 >           jet.set_photonEnergyFraction (pat_jet.photonEnergyFraction());
596 >           jet.set_chargedMultiplicity (pat_jet.chargedMultiplicity());
597 >           jet.set_neutralMultiplicity (pat_jet.neutralMultiplicity());
598 >           jet.set_muonMultiplicity (pat_jet.muonMultiplicity());
599 >           jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
600 >           jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
601 >         }
602 >         if(doJECUncertainty){
603 >           jecUnc->setJetEta(pat_jet.eta());
604 >           jecUnc->setJetPt(pat_jet.pt());
605 >           jet.set_JEC_uncertainty(jecUnc->getUncertainty(true));
606 >         }
607 >         jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
608 >        
609 >         jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
610 >         jet.set_btag_simpleSecondaryVertexHighPur(pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
611 >         jet.set_btag_combinedSecondaryVertex(pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags"));
612 >         jet.set_btag_combinedSecondaryVertexMVA(pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
613 >         jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
614 >         jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
615  
616 +        
617 +         const reco::GenJet *genj = pat_jet.genJet();
618 +         if(genj){
619 +           jet.set_genjet_pt(genj->pt());
620 +           jet.set_genjet_eta(genj->eta());
621 +           jet.set_genjet_phi(genj->phi());
622 +           jet.set_genjet_energy(genj->energy());
623 +           if(doAllGenParticles){
624 +             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
625 +             for(unsigned int l = 0; l<jetgenps.size(); ++l){
626 +               for(unsigned int k=0; k< genps.size(); ++k){
627 +                 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
628 +                   jet.add_genparticles_index(genps[k].index());
629 +                   break;
630 +                 }
631 +               }
632 +             }
633 +             if(jet.genparticles_indices().size()!= jetgenps.size())
634 +               std::cout << "WARNING: Found only " << jet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
635 +           }
636 +          
637 +         }
638 +        
639           jets[j].push_back(jet);
640         }
641       }
642     }
643  
644 +   //-------------- gen jets -------------
645 +
646 +   if(doGenJets){
647 +     for(size_t j=0; j< genjet_sources.size(); ++j){
648 +      
649 +       genjets[j].clear();
650 +
651 +       edm::Handle< std::vector<reco::GenJet> > genjet_handle;
652 +       iEvent.getByLabel(genjet_sources[j], genjet_handle);
653 +       const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
654 +  
655 +       for (unsigned int i = 0; i < gen_jets.size(); ++i) {
656 +         pat::Jet gen_jet = gen_jets[i];
657 +         if(gen_jet.pt() < genjet_ptmin) continue;
658 +         if(fabs(gen_jet.eta()) > genjet_etamax) continue;
659 +
660 +         Particle jet;
661 +         jet.set_charge(gen_jet.charge());
662 +         jet.set_pt(gen_jet.pt());
663 +         jet.set_eta(gen_jet.eta());
664 +         jet.set_phi(gen_jet.phi());
665 +         jet.set_energy(gen_jet.energy());
666 +
667 +         genjets[j].push_back(jet);
668 +
669 +       }
670 +     }
671 +   }
672 +
673     // ------------- top jets -------------
674     if(doTopJets){
675       for(size_t j=0; j< topjet_sources.size(); ++j){
# Line 434 | Line 681 | NtupleWriter::analyze(const edm::Event&
681         iEvent.getByLabel(topjet_sources[j], pat_topjets);
682  
683         for (unsigned int i = 0; i < pat_topjets->size(); i++) {
684 +
685           const pat::Jet  pat_topjet =  * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
686           if(pat_topjet.pt() < topjet_ptmin) continue;
687           if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
688  
689           TopJet topjet;
690 <         topjet.charge = pat_topjet.charge();
691 <         topjet.pt = pat_topjet.pt();
692 <         topjet.eta = pat_topjet.eta();
693 <         topjet.phi = pat_topjet.phi();
694 <         topjet.energy = pat_topjet.energy();
695 <         topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
690 >         topjet.set_charge(pat_topjet.charge());
691 >         topjet.set_pt(pat_topjet.pt());
692 >         topjet.set_eta(pat_topjet.eta());
693 >         topjet.set_phi(pat_topjet.phi());
694 >         topjet.set_energy(pat_topjet.energy());
695 >         topjet.set_numberOfDaughters(pat_topjet.numberOfDaughters());
696           const reco::TrackRefVector&  topjettracks = pat_topjet.associatedTracks();
697 <         topjet.nTracks = topjettracks.size();
698 <         topjet.jetArea = pat_topjet.jetArea();
699 <         topjet.pileup = pat_topjet.pileup();
700 < //       topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
701 < //       topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
702 < //       topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
703 < //       topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
704 < //       topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
705 < //       topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
706 < //       topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
707 < //       topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
708 < //       topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
709 < //       topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
710 < //       topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
711 <
712 <         jecUnc->setJetEta(pat_topjet.eta());
713 <         jecUnc->setJetPt(pat_topjet.pt());
714 <         topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
715 <
716 <         topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
717 <         topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
718 <         topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
719 <         topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
720 <         topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
721 <         topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
697 >         topjet.set_nTracks( topjettracks.size());
698 >         topjet.set_jetArea( pat_topjet.jetArea());
699 >         topjet.set_pileup( pat_topjet.pileup());
700 > //       topjet.set_neutralEmEnergyFraction(pat_topjet.neutralEmEnergyFraction());
701 > //       topjet.set_neutralHadronEnergyFraction(pat_topjet.neutralHadronEnergyFraction());
702 > //       topjet.set_chargedEmEnergyFraction(pat_topjet.chargedEmEnergyFraction());
703 > //       topjet.set_chargedHadronEnergyFraction(pat_topjet.chargedHadronEnergyFraction());
704 > //       topjet.set_muonEnergyFraction(pat_topjet.muonEnergyFraction());
705 > //       topjet.set_photonEnergyFraction(pat_topjet.photonEnergyFraction());
706 > //       topjet.set_chargedMultiplicity(pat_topjet.chargedMultiplicity());
707 > //       topjet.set_neutralMultiplicity(pat_topjet.neutralMultiplicity());
708 > //       topjet.set_muonMultiplicity(pat_topjet.muonMultiplicity());
709 > //       topjet.set_electronMultiplicity(pat_topjet.electronMultiplicity());
710 > //       topjet.set_photonMultiplicity(pat_topjet.photonMultiplicity());
711 >         if(doJECUncertainty){
712 >           jecUnc->setJetEta(pat_topjet.eta());
713 >           jecUnc->setJetPt(pat_topjet.pt());
714 >           topjet.set_JEC_uncertainty( jecUnc->getUncertainty(true));
715 >         }
716 >         topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
717 >
718 >         topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
719 >         topjet.set_btag_simpleSecondaryVertexHighPur(pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
720 >         topjet.set_btag_combinedSecondaryVertex(pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags"));
721 >         topjet.set_btag_combinedSecondaryVertexMVA(pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
722 >         topjet.set_btag_jetBProbability(pat_topjet.bDiscriminator("jetBProbabilityBJetTags"));
723 >         topjet.set_btag_jetProbability(pat_topjet.bDiscriminator("jetProbabilityBJetTags"));
724 >
725 >         /*
726 >         const reco::GenJet *genj = pat_topjet.genJet();
727 >         if(genj){
728 >           topjet.set_genjet_pt ( genj->pt());
729 >           topjet.set_genjet_eta ( genj->eta());
730 >           topjet.set_genjet_phi ( genj->phi());
731 >           topjet.set_genjet_energy ( genj->energy());
732 >           if(doAllGenParticles){
733 >             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
734 >             for(unsigned int l = 0; l<jetgenps.size(); ++l){
735 >               for(unsigned int k=0; k< genps.size(); ++k){
736 >                 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
737 >                   topjet.add_genparticles_index(genps[k].index());
738 >                   break;
739 >                 }
740 >               }
741 >             }
742 >             if(topjet.genparticles_indices().size()!= jetgenps.size())
743 >               std::cout << "WARNING: Found only " << topjet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
744 >           }
745 >         }
746 >         */
747  
748           for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
749             Particle subjet_v4;
750 <           subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
751 <           subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
752 <           subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
753 <           subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
754 <           topjet.subjets.push_back(subjet_v4);
750 >
751 >           reco::Candidate const * subjetd =  pat_topjet.daughter(k);
752 >           pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
753 >           if(patsubjetd)
754 >           {
755 >              subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
756 >              subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
757 >              subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
758 >              subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
759 >              topjet.add_subjet(subjet_v4);
760 >              //btag info
761 >              topjet.add_subFlavour(patsubjetd->partonFlavour());
762 >              topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
763 >              if (doTagInfos)
764 >                {
765 >                  //ip tag info
766 >                  reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
767 >                  topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
768 >                  topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
769 >                  topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
770 >                  topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
771 >                  topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
772 >                  topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
773 >                  topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
774 >                  topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
775 >                  topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
776 >                  topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
777 >                  topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
778 >                  topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));    
779 >                  topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
780 >                  topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
781 >                  topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
782 >                  topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
783 >                  topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
784 >                  topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
785 >                  topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
786 >                  topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
787 >                  topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
788 >                  //sv tag info
789 >                  reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
790 >                  topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
791 >                  topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
792 >                  topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
793 >                  topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
794 >                  topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
795 >                  topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
796 >                  topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
797 >                  std::vector<TLorentzVector> vp4; vp4.clear();
798 >                  std::vector<float> vchi2; vchi2.clear();
799 >                  std::vector<float> vndof; vndof.clear();
800 >                  std::vector<float> vchi2ndof; vchi2ndof.clear();
801 >                  std::vector<float> vtrsize; vtrsize.clear();
802 >                  for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
803 >                    {
804 >                      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
805 >                      vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
806 >                      vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());  
807 >                      vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());  
808 >                      vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());  
809 >                      vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());  
810 >                    }
811 >                  topjet.add_subSecondaryVertex(vp4);
812 >                  topjet.add_subVertexChi2(vchi2);
813 >                  topjet.add_subVertexNdof(vndof);
814 >                  topjet.add_subVertexNormalizedChi2(vchi2ndof);
815 >                  topjet.add_subVertexTracksSize(vtrsize);
816 >                  //try computer
817 >                  edm::ESHandle<JetTagComputer> computerHandle;
818 >                  iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
819 >                  const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
820 >                  if(computer)
821 >                    {
822 >                      computer->passEventSetup(iSetup);
823 >                      std::vector<const reco::BaseTagInfo*>  baseTagInfos;
824 >                      baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
825 >                      baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );      
826 >                      JetTagComputer::TagInfoHelper helper(baseTagInfos);
827 >                      reco::TaggingVariableList vars = computer->taggingVariables(helper);
828 >                      topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
829 >                      topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
830 >                      topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
831 >                      topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
832 >                    }
833 >                }
834 >           }
835 >           else
836 >             {
837 >               //filling only standard information in case the subjet has not been pat-tified during the pattuples production
838 >               subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
839 >               subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
840 >               subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
841 >               subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
842 >               topjet.add_subjet(subjet_v4);
843 >             }
844 >          
845 >          
846           }
847           topjets[j].push_back(topjet);
848         }
849       }
850     }
851 +  
852 +  
853 +   // ------------- generator top jets -------------
854 +   if(doGenTopJets){
855 +     for(size_t j=0; j< gentopjet_sources.size(); ++j){
856 +      
857 +       gentopjets[j].clear();
858 +      
859 +       edm::Handle<reco::BasicJetCollection> reco_gentopjets;
860 +       //edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
861 +       iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
862 +
863 +       for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
864 +        
865 +         const reco::BasicJet  reco_gentopjet =  reco_gentopjets->at(i);
866 +         if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
867 +         if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
868 +
869 +         TopJet gentopjet;
870 +         gentopjet.set_charge(reco_gentopjet.charge());
871 +         gentopjet.set_pt(reco_gentopjet.pt());
872 +         gentopjet.set_eta(reco_gentopjet.eta());
873 +         gentopjet.set_phi(reco_gentopjet.phi());
874 +         gentopjet.set_energy(reco_gentopjet.energy());
875 +         gentopjet.set_numberOfDaughters(reco_gentopjet.numberOfDaughters());
876 +
877 +         for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
878 +           Particle subjet_v4;
879 +           subjet_v4.set_pt(reco_gentopjet.daughter(k)->p4().pt());
880 +           subjet_v4.set_eta(reco_gentopjet.daughter(k)->p4().eta());
881 +           subjet_v4.set_phi(reco_gentopjet.daughter(k)->p4().phi());
882 +           subjet_v4.set_energy(reco_gentopjet.daughter(k)->p4().E());
883 +           gentopjet.add_subjet(subjet_v4);
884 +         }
885 +         gentopjets[j].push_back(gentopjet);
886 +       }
887 +     }
888 +   }
889  
890     // ------------- photons -------------
891     if(doPhotons){
# Line 497 | Line 899 | NtupleWriter::analyze(const edm::Event&
899         for (unsigned int i = 0; i < pat_photons.size(); ++i) {
900           pat::Photon pat_photon = pat_photons[i];
901           Photon ph;
902 <         ph.charge = 0;
903 <         ph.pt =  pat_photon.pt();
904 <         ph.eta =  pat_photon.eta();
905 <         ph.phi =  pat_photon.phi();
906 <         ph.energy =  pat_photon.energy();
907 <         ph.vertex_x = pat_photon.vertex().x();
908 <         ph.vertex_y = pat_photon.vertex().y();
909 <         ph.vertex_z = pat_photon.vertex().z();
910 <         ph.supercluster_eta = pat_photon.superCluster()->eta();
911 <         ph.supercluster_phi = pat_photon.superCluster()->phi();
912 < //       ph.neutralHadronIso = pat_photon.neutralHadronIso();
913 < //       ph.chargedHadronIso = pat_photon.chargedHadronIso();
914 <         ph.trackIso = pat_photon.trackIso();
902 >         ph.set_charge(0);
903 >         ph.set_pt( pat_photon.pt());
904 >         ph.set_eta( pat_photon.eta());
905 >         ph.set_phi( pat_photon.phi());
906 >         ph.set_energy( pat_photon.energy());
907 >         ph.set_vertex_x(pat_photon.vertex().x());
908 >         ph.set_vertex_y(pat_photon.vertex().y());
909 >         ph.set_vertex_z(pat_photon.vertex().z());
910 >         ph.set_supercluster_eta(pat_photon.superCluster()->eta());
911 >         ph.set_supercluster_phi(pat_photon.superCluster()->phi());
912 > //       ph.set_neutralHadronIso(pat_photon.neutralHadronIso());
913 > //       ph.set_chargedHadronIso(pat_photon.chargedHadronIso());
914 >         ph.set_trackIso(pat_photon.trackIso());
915           phs[j].push_back(ph);
916         }
917       }
# Line 529 | Line 931 | NtupleWriter::analyze(const edm::Event&
931         else{
932           pat::MET pat_met = pat_mets[0];
933          
934 <         met[j].pt=pat_met.pt();
935 <         met[j].phi=pat_met.phi();
936 <         met[j].mEtSig= pat_met.mEtSig();
934 >         met[j].set_pt(pat_met.pt());
935 >         met[j].set_phi(pat_met.phi());
936 >         met[j].set_mEtSig(pat_met.mEtSig());
937         }
938        
939       }
# Line 554 | Line 956 | NtupleWriter::analyze(const edm::Event&
956      
957       triggerResults.clear();
958       triggerNames.clear();
959 <     L1_prescale.clear();
960 <     HLT_prescale.clear();
959 > //      L1_prescale.clear();
960 > //      HLT_prescale.clear();
961      
962       edm::Service<edm::service::TriggerNamesService> tns;
963       std::vector<std::string> triggerNames_all;
# Line 572 | Line 974 | NtupleWriter::analyze(const edm::Event&
974         //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
975         triggerResults.push_back(trig.accept(i));
976         if(newrun) triggerNames.push_back(triggerNames_all[i]);
977 <       if(isRealData){
978 <         std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
979 <         L1_prescale.push_back(pre.first);
980 <         HLT_prescale.push_back(pre.second);
981 <       }
977 > //        if(isRealData){
978 > //       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
979 > //       L1_prescale.push_back(pre.first);
980 > //       HLT_prescale.push_back(pre.second);
981 > //       //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
982 > //        }
983       }
984       //    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
985       //      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
# Line 584 | Line 987 | NtupleWriter::analyze(const edm::Event&
987       newrun=false;
988     }
989  
587   // ------------- generator info -------------
588  
589   if(doGenInfo){
590     genInfo.weights.clear();
591     genInfo.binningValues.clear();
592     genps.clear();
593
594     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
595     iEvent.getByLabel("generator", genEventInfoProduct);
596     const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
597  
598     genInfo.binningValues = genEventInfo.binningValues();
599     genInfo.weights = genEventInfo.weights();
600     genInfo.alphaQCD = genEventInfo.alphaQCD();
601     genInfo.alphaQED = genEventInfo.alphaQED();
602     genInfo.qScale = genEventInfo.qScale();
603    
604     const gen::PdfInfo* pdf = genEventInfo.pdf();
605     if(pdf){
606       genInfo.pdf_id1=pdf->id.first;
607       genInfo.pdf_id2=pdf->id.second;
608       genInfo.pdf_x1=pdf->x.first;
609       genInfo.pdf_x2=pdf->x.second;
610       genInfo.pdf_scalePDF=pdf->scalePDF;
611       genInfo.pdf_xPDF1=pdf->xPDF.first;
612       genInfo.pdf_xPDF2=pdf->xPDF.second;
613     }
614     else{
615       genInfo.pdf_id1=-999;
616       genInfo.pdf_id2=-999;
617       genInfo.pdf_x1=-999;
618       genInfo.pdf_x2=-999;
619       genInfo.pdf_scalePDF=-999;
620       genInfo.pdf_xPDF1=-999;
621       genInfo.pdf_xPDF2=-999;
622     }
623
624     edm::Handle<std::vector<PileupSummaryInfo> > pus;
625     iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
626     genInfo.pileup_NumInteractions_intime=0;
627     genInfo.pileup_NumInteractions_ootbefore=0;
628     genInfo.pileup_NumInteractions_ootafter=0;
629     if(pus.isValid()){
630       genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
631       for(size_t i=0; i<pus->size(); ++i){
632         if(pus->at(i).getBunchCrossing() == 0) // intime pileup
633            genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
634         else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
635            genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
636         }
637         else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
638           genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
639         }
640       }
641     }
642
643     edm::Handle<reco::GenParticleCollection> genPartColl;
644     iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
645     int index=-1;
646     for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
647       index++;
648      
649       //write out only top quarks and status 3 particles (works fine only for MadGraph)
650       if(abs(iter->pdgId())==6 || iter->status()==3){
651         GenParticle genp;
652         genp.charge = iter->charge();;
653         genp.pt = iter->p4().pt();
654         genp.eta = iter->p4().eta();
655         genp.phi = iter->p4().phi();
656         genp.energy = iter->p4().E();
657         genp.index =index;
658         genp.status = iter->status();
659         genp.pdgId = iter->pdgId();
660
661         genp.mother1=-1;
662         genp.mother2=-1;
663         genp.daughter1=-1;
664         genp.daughter2=-1;
665        
666         int nm=iter->numberOfMothers();
667         int nd=iter->numberOfDaughters();
668        
669         if (nm>0) genp.mother1 = iter->motherRef(0).key();
670         if (nm>1) genp.mother2 = iter->motherRef(nm-1).key();
671         if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
672         if (nd>1) genp.daughter2 = iter->daughterRef(nd-1).key();
673        
674         genps.push_back(genp);
675       }
676     }
677
678   }
679
990  
991     tr->Fill();
992 <
992 >   if(doLumiInfo)
993 >     previouslumiblockwasfilled=true;
994   }
995  
996  
# Line 687 | Line 998 | NtupleWriter::analyze(const edm::Event&
998   void
999   NtupleWriter::beginJob()
1000   {
1001 +  if(doLumiInfo){
1002 +    totalRecLumi=0;
1003 +    totalDelLumi=0;
1004 +    previouslumiblockwasfilled=false;
1005 +  }
1006   }
1007  
1008   // ------------ method called once each job just after ending the event loop  ------------
# Line 702 | Line 1018 | NtupleWriter::endJob()
1018   void
1019   NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
1020   {
1021 <  bool setup_changed = false;
1022 <  hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1023 <  newrun=true;
1024 <
1025 <  edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
710 <  iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
711 <  JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
712 <  jecUnc = new JetCorrectionUncertainty(JetCorPar);
1021 >  if(doTrigger){
1022 >    //bool setup_changed = false;
1023 >    //hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1024 >    newrun=true;
1025 >  }
1026  
1027 +  if(doJets || doTopJets){
1028 +    if(doJECUncertainty){
1029 +      edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
1030 +      iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
1031 +      JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
1032 +      jecUnc = new JetCorrectionUncertainty(JetCorPar);
1033 +    }
1034 +  }
1035   }
1036  
1037   // ------------ method called when ending the processing of a run  ------------
1038   void
1039   NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
1040   {
1041 +  if(doLumiInfo)
1042 +    std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
1043   }
1044  
1045   // ------------ method called when starting to processes a luminosity block  ------------
1046   void
1047 < NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
1047 > NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
1048   {
1049 +  if(doLumiInfo){
1050 +    edm::Handle<LumiSummary> l;
1051 +    lumi.getByLabel("lumiProducer", l);
1052 +    
1053 +    //add lumi of lumi blocks without any event to next lumiblock
1054 +    if(previouslumiblockwasfilled){
1055 +      intgRecLumi=0;
1056 +      intgDelLumi=0;
1057 +    }
1058 +    previouslumiblockwasfilled=false;
1059 +    
1060 +    if (l.isValid()){;
1061 +      intgRecLumi+=l->intgRecLumi()*6.37;
1062 +      intgDelLumi+=l->intgDelLumi()*6.37;
1063 +      totalRecLumi+=l->intgRecLumi()*6.37;
1064 +      totalDelLumi+=l->intgDelLumi()*6.37;
1065 +    }
1066 +    //std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<"   " << l->intgDelLumi()*6.37<<std::endl;
1067 +    //std::cout << "summed:  "<< intgRecLumi << "   " << intgDelLumi << std::endl;
1068 +  }
1069   }
1070  
1071   // ------------ method called when ending the processing of a luminosity block  ------------

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines