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.14 by peiffer, Mon Apr 23 14:26:23 2012 UTC vs.
Revision 1.25 by peiffer, Thu Jun 6 07:38:31 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");
# Line 59 | Line 64 | NtupleWriter::NtupleWriter(const edm::Pa
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);
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);
# Line 104 | 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 139 | 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 147 | 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 181 | 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 202 | 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         }
# Line 223 | Line 245 | NtupleWriter::analyze(const edm::Event&
245       beamspot_z0 = bsp.z0();
246     }
247  
248 + // ------------- generator info -------------
249 +  
250 +   if(doGenInfo){
251 +     genInfo.clear_weights();
252 +     genInfo.clear_binningValues();
253 +     genps.clear();
254 +
255 +     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
256 +     iEvent.getByLabel("generator", genEventInfoProduct);
257 +     const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
258 +  
259 +     for(unsigned int k=0; k<genEventInfo.binningValues().size();++k){
260 +       genInfo.add_binningValue(genEventInfo.binningValues().at(k));
261 +     }
262 +     for(unsigned int k=0; k<genEventInfo.weights().size();++k){
263 +       genInfo.add_weight(genEventInfo.weights().at(k));
264 +     }
265 +     genInfo.set_alphaQCD(genEventInfo.alphaQCD());
266 +     genInfo.set_alphaQED(genEventInfo.alphaQED());
267 +     genInfo.set_qScale(genEventInfo.qScale());
268 +    
269 +     const gen::PdfInfo* pdf = genEventInfo.pdf();
270 +     if(pdf){
271 +       genInfo.set_pdf_id1(pdf->id.first);
272 +       genInfo.set_pdf_id2(pdf->id.second);
273 +       genInfo.set_pdf_x1(pdf->x.first);
274 +       genInfo.set_pdf_x2(pdf->x.second);
275 +       genInfo.set_pdf_scalePDF(pdf->scalePDF);
276 +       genInfo.set_pdf_xPDF1(pdf->xPDF.first);
277 +       genInfo.set_pdf_xPDF2(pdf->xPDF.second);
278 +     }
279 +     else{
280 +       genInfo.set_pdf_id1(-999);
281 +       genInfo.set_pdf_id2(-999);
282 +       genInfo.set_pdf_x1(-999);
283 +       genInfo.set_pdf_x2(-999);
284 +       genInfo.set_pdf_scalePDF(-999);
285 +       genInfo.set_pdf_xPDF1(-999);
286 +       genInfo.set_pdf_xPDF2(-999);
287 +     }
288 +
289 +     edm::Handle<std::vector<PileupSummaryInfo> > pus;
290 +     iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
291 +     genInfo.set_pileup_NumInteractions_intime(0);
292 +     genInfo.set_pileup_NumInteractions_ootbefore(0);
293 +     genInfo.set_pileup_NumInteractions_ootafter(0);
294 +     if(pus.isValid()){
295 +       genInfo.set_pileup_TrueNumInteractions ( (float) pus->at(0).getTrueNumInteractions());
296 +       for(size_t i=0; i<pus->size(); ++i){
297 +         if(pus->at(i).getBunchCrossing() == 0) // intime pileup
298 +           genInfo.set_pileup_NumInteractions_intime( genInfo.pileup_NumInteractions_intime() + pus->at(i).getPU_NumInteractions());
299 +         else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
300 +           genInfo.set_pileup_NumInteractions_ootbefore( genInfo.pileup_NumInteractions_ootbefore() + pus->at(i).getPU_NumInteractions());
301 +         }
302 +         else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
303 +           genInfo.set_pileup_NumInteractions_ootafter( genInfo.pileup_NumInteractions_ootafter() + pus->at(i).getPU_NumInteractions());
304 +         }
305 +       }
306 +     }
307 +
308 +     edm::Handle<reco::GenParticleCollection> genPartColl;
309 +     iEvent.getByLabel(genparticle_source, genPartColl);
310 +     int index=-1;
311 +     for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
312 +       index++;
313 +      
314 +       //write out only top quarks 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);
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;
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();
# Line 243 | 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.photonIso = pat_ele.photonIso();
381 <         ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
382 <         ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
383 <         ele.gsfTrack_px= pat_ele.gsfTrack()->px();
384 <         ele.gsfTrack_py= pat_ele.gsfTrack()->py();
385 <         ele.gsfTrack_pz= pat_ele.gsfTrack()->pz();
386 <         ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
387 <         ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
388 <         ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
389 <         ele.passconversionveto = !ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position());
390 <         ele.dEtaIn=pat_ele.deltaEtaSuperClusterTrackAtVtx();
391 <         ele.dPhiIn=pat_ele.deltaPhiSuperClusterTrackAtVtx();
392 <         ele.sigmaIEtaIEta=pat_ele.sigmaIetaIeta();
393 <         ele.HoverE=pat_ele.hadronicOverEm();
394 <         ele.fbrem=pat_ele.fbrem();
395 <         ele.EoverPIn=pat_ele.eSuperClusterOverP();
396 <         ele.EcalEnergy=pat_ele.ecalEnergy();
397 <         //ele.mvaTrigV0=pat_ele.electronID("mvaTrigV0");
398 <         //ele.mvaNonTrigV0=pat_ele.electronID("mvaNonTrigV0");
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         }
# Line 296 | 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.photonIso = pat_mu.photonIso();
440 <         mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
441 <         mu.isGlobalMuon = pat_mu.isGlobalMuon();
442 <         mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
443 <         mu.isTrackerMuon = pat_mu.isTrackerMuon();
444 <         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();
472 <           mu.innerTrack_trackerLayersWithMeasurement = innerTrack->hitPattern().trackerLayersWithMeasurement();
473 <           mu.innerTrack_numberOfValidPixelHits = innerTrack->hitPattern().numberOfValidPixelHits();
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;
482 <           mu.innerTrack_trackerLayersWithMeasurement = 0;
483 <           mu.innerTrack_numberOfValidPixelHits = 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 392 | 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 <         tau.decayModeFinding = pat_tau.tauID("decayModeFinding")>0.5;
530 <         tau.byVLooseCombinedIsolationDeltaBetaCorr  = pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5;
531 <         tau.byLooseCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5;
532 <         tau.byMediumCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5;
533 <         tau.byTightCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5;
534 <         tau.againstElectronLoose  = pat_tau.tauID("againstElectronLoose")>0.5;
535 <         tau.againstElectronMedium = pat_tau.tauID("againstElectronMedium")>0.5;
536 <         tau.againstElectronTight = pat_tau.tauID("againstElectronTight")>0.5;
537 <         tau.againstElectronMVA  = pat_tau.tauID("againstElectronMVA")>0.5;
538 <         tau.againstMuonLoose = pat_tau.tauID("againstMuonLoose")>0.5;
539 <         tau.againstMuonMedium = pat_tau.tauID("againstMuonMedium")>0.5;
540 <         tau.againstMuonTight = pat_tau.tauID("againstMuonTight")>0.5;
541 <
542 <         reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
543 <         if(!leadPFCand.isNull()){
544 <           tau.leadPFCand_px = leadPFCand->px();
545 <           tau.leadPFCand_py = leadPFCand->py();
546 <           tau.leadPFCand_pz = leadPFCand->pz();
547 <         }
548 <         else{
549 <           tau.leadPFCand_px = 0;
550 <           tau.leadPFCand_py = 0;
551 <           tau.leadPFCand_pz = 0;
552 <         }
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       }
556     }
557  
558 +   //-------------- gen jets -------------
559 +
560 +   if(doGenJets){
561 +     for(size_t j=0; j< genjet_sources.size(); ++j){
562 +      
563 +       genjets[j].clear();
564 +
565 +       edm::Handle< std::vector<reco::GenJet> > genjet_handle;
566 +       iEvent.getByLabel(genjet_sources[j], genjet_handle);
567 +       const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
568 +  
569 +       for (unsigned int i = 0; i < gen_jets.size(); ++i) {
570 +         const reco::GenJet* gen_jet = &gen_jets[i];
571 +         if(gen_jet->pt() < genjet_ptmin) continue;
572 +         if(fabs(gen_jet->eta()) > genjet_etamax) continue;
573 +
574 +         Particle jet;
575 +         jet.set_charge(gen_jet->charge());
576 +         jet.set_pt(gen_jet->pt());
577 +         jet.set_eta(gen_jet->eta());
578 +         jet.set_phi(gen_jet->phi());
579 +         jet.set_energy(gen_jet->energy());
580 +
581 +         genjets[j].push_back(jet);
582 +
583 +       }
584 +     }
585 +   }
586 +
587     // ------------- jets -------------
588     if(doJets){
589       for(size_t j=0; j< jet_sources.size(); ++j){
# Line 447 | Line 605 | NtupleWriter::analyze(const edm::Event&
605   //       }
606  
607           Jet jet;
608 <         jet.charge = pat_jet.charge();
609 <         jet.pt = pat_jet.pt();
610 <         jet.eta = pat_jet.eta();
611 <         jet.phi = pat_jet.phi();
612 <         jet.energy = pat_jet.energy();
613 <         jet.numberOfDaughters =pat_jet.numberOfDaughters();
608 >         jet.set_charge(pat_jet.charge());
609 >         jet.set_pt(pat_jet.pt());
610 >         jet.set_eta(pat_jet.eta());
611 >         jet.set_phi(pat_jet.phi());
612 >         jet.set_energy(pat_jet.energy());
613 >         jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
614           const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
615 <         jet.nTracks = jettracks.size();
616 <         jet.jetArea = pat_jet.jetArea();
617 <         jet.pileup = pat_jet.pileup();
615 >         jet.set_nTracks ( jettracks.size());
616 >         jet.set_jetArea(pat_jet.jetArea());
617 >         jet.set_pileup(pat_jet.pileup());
618           if(pat_jet.isPFJet()){
619 <           jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
620 <           jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
621 <           jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
622 <           jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
623 <           jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
624 <           jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
625 <           jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
626 <           jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
627 <           jet.muonMultiplicity =pat_jet.muonMultiplicity();
628 <           jet.electronMultiplicity =pat_jet.electronMultiplicity();
629 <           jet.photonMultiplicity =pat_jet.photonMultiplicity();
630 <         }
631 <
632 <         jecUnc->setJetEta(pat_jet.eta());
633 <         jecUnc->setJetPt(pat_jet.pt());
634 <         jet.JEC_uncertainty = jecUnc->getUncertainty(true);
635 <         jet.JEC_factor_raw = pat_jet.jecFactor("Uncorrected");
636 <
637 <         jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
638 <         jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
639 <         jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
640 <         jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
641 <         jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
642 <         jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
619 >           jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
620 >           jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
621 >           jet.set_chargedEmEnergyFraction (pat_jet.chargedEmEnergyFraction());
622 >           jet.set_chargedHadronEnergyFraction (pat_jet.chargedHadronEnergyFraction());
623 >           jet.set_muonEnergyFraction (pat_jet.muonEnergyFraction());
624 >           jet.set_photonEnergyFraction (pat_jet.photonEnergyFraction());
625 >           jet.set_chargedMultiplicity (pat_jet.chargedMultiplicity());
626 >           jet.set_neutralMultiplicity (pat_jet.neutralMultiplicity());
627 >           jet.set_muonMultiplicity (pat_jet.muonMultiplicity());
628 >           jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
629 >           jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
630 >         }
631 >         if(doJECUncertainty){
632 >           jecUnc->setJetEta(pat_jet.eta());
633 >           jecUnc->setJetPt(pat_jet.pt());
634 >           jet.set_JEC_uncertainty(jecUnc->getUncertainty(true));
635 >         }
636 >         jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
637 >        
638 >         jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
639 >         jet.set_btag_simpleSecondaryVertexHighPur(pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
640 >         jet.set_btag_combinedSecondaryVertex(pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags"));
641 >         jet.set_btag_combinedSecondaryVertexMVA(pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
642 >         jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
643 >         jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
644  
645 +        
646           const reco::GenJet *genj = pat_jet.genJet();
647           if(genj){
488           jet.genjet_pt = genj->pt();
489           jet.genjet_eta = genj->eta();
490           jet.genjet_phi = genj->phi();
491           jet.genjet_energy = genj->energy();
492         }
648  
649 +           for(unsigned int k=0; k<genjets->size(); ++k){
650 +             if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta())
651 +               jet.set_genjet_index(k);
652 +           }
653 + //         if( jet.genjet_index()<0){
654 + //           std::cout<< "genjet not found for " << genj->pt() << "  " << genj->eta() << std::endl;
655 + //         }
656 +
657 +           if(doAllGenParticles){
658 +             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
659 +             for(unsigned int l = 0; l<jetgenps.size(); ++l){
660 +               for(unsigned int k=0; k< genps.size(); ++k){
661 +                 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
662 +                   jet.add_genparticles_index(genps[k].index());
663 +                   break;
664 +                 }
665 +               }
666 +             }
667 +             if(jet.genparticles_indices().size()!= jetgenps.size())
668 +               std::cout << "WARNING: Found only " << jet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
669 +           }
670 +          
671 +         }
672 +        
673           jets[j].push_back(jet);
674         }
675       }
676     }
677  
678 +
679 +
680     // ------------- top jets -------------
681     if(doTopJets){
682       for(size_t j=0; j< topjet_sources.size(); ++j){
# Line 513 | Line 694 | NtupleWriter::analyze(const edm::Event&
694           if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
695  
696           TopJet topjet;
697 <         topjet.charge = pat_topjet.charge();
698 <         topjet.pt = pat_topjet.pt();
699 <         topjet.eta = pat_topjet.eta();
700 <         topjet.phi = pat_topjet.phi();
701 <         topjet.energy = pat_topjet.energy();
702 <         topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
697 >         topjet.set_charge(pat_topjet.charge());
698 >         topjet.set_pt(pat_topjet.pt());
699 >         topjet.set_eta(pat_topjet.eta());
700 >         topjet.set_phi(pat_topjet.phi());
701 >         topjet.set_energy(pat_topjet.energy());
702 >         topjet.set_numberOfDaughters(pat_topjet.numberOfDaughters());
703           const reco::TrackRefVector&  topjettracks = pat_topjet.associatedTracks();
704 <         topjet.nTracks = topjettracks.size();
705 <         topjet.jetArea = pat_topjet.jetArea();
706 <         topjet.pileup = pat_topjet.pileup();
707 < //       topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
708 < //       topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
709 < //       topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
710 < //       topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
711 < //       topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
712 < //       topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
713 < //       topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
714 < //       topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
715 < //       topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
716 < //       topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
717 < //       topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
718 <
719 <         jecUnc->setJetEta(pat_topjet.eta());
720 <         jecUnc->setJetPt(pat_topjet.pt());
721 <         topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
722 <         topjet.JEC_factor_raw = pat_topjet.jecFactor("Uncorrected");
723 <
543 <         topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
544 <         topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
545 <         topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
546 <         topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
547 <         topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
548 <         topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
704 >         topjet.set_nTracks( topjettracks.size());
705 >         topjet.set_jetArea( pat_topjet.jetArea());
706 >         topjet.set_pileup( pat_topjet.pileup());
707 > //       topjet.set_neutralEmEnergyFraction(pat_topjet.neutralEmEnergyFraction());
708 > //       topjet.set_neutralHadronEnergyFraction(pat_topjet.neutralHadronEnergyFraction());
709 > //       topjet.set_chargedEmEnergyFraction(pat_topjet.chargedEmEnergyFraction());
710 > //       topjet.set_chargedHadronEnergyFraction(pat_topjet.chargedHadronEnergyFraction());
711 > //       topjet.set_muonEnergyFraction(pat_topjet.muonEnergyFraction());
712 > //       topjet.set_photonEnergyFraction(pat_topjet.photonEnergyFraction());
713 > //       topjet.set_chargedMultiplicity(pat_topjet.chargedMultiplicity());
714 > //       topjet.set_neutralMultiplicity(pat_topjet.neutralMultiplicity());
715 > //       topjet.set_muonMultiplicity(pat_topjet.muonMultiplicity());
716 > //       topjet.set_electronMultiplicity(pat_topjet.electronMultiplicity());
717 > //       topjet.set_photonMultiplicity(pat_topjet.photonMultiplicity());
718 >         if(doJECUncertainty){
719 >           jecUnc->setJetEta(pat_topjet.eta());
720 >           jecUnc->setJetPt(pat_topjet.pt());
721 >           topjet.set_JEC_uncertainty( jecUnc->getUncertainty(true));
722 >         }
723 >         topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
724  
725 +         topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
726 +         topjet.set_btag_simpleSecondaryVertexHighPur(pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
727 +         topjet.set_btag_combinedSecondaryVertex(pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags"));
728 +         topjet.set_btag_combinedSecondaryVertexMVA(pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
729 +         topjet.set_btag_jetBProbability(pat_topjet.bDiscriminator("jetBProbabilityBJetTags"));
730 +         topjet.set_btag_jetProbability(pat_topjet.bDiscriminator("jetProbabilityBJetTags"));
731 +
732 +         /*
733           const reco::GenJet *genj = pat_topjet.genJet();
734           if(genj){
735 <           topjet.genjet_pt = genj->pt();
736 <           topjet.genjet_eta = genj->eta();
737 <           topjet.genjet_phi = genj->phi();
738 <           topjet.genjet_energy = genj->energy();
735 >           topjet.set_genjet_pt ( genj->pt());
736 >           topjet.set_genjet_eta ( genj->eta());
737 >           topjet.set_genjet_phi ( genj->phi());
738 >           topjet.set_genjet_energy ( genj->energy());
739 >           if(doAllGenParticles){
740 >             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
741 >             for(unsigned int l = 0; l<jetgenps.size(); ++l){
742 >               for(unsigned int k=0; k< genps.size(); ++k){
743 >                 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
744 >                   topjet.add_genparticles_index(genps[k].index());
745 >                   break;
746 >                 }
747 >               }
748 >             }
749 >             if(topjet.genparticles_indices().size()!= jetgenps.size())
750 >               std::cout << "WARNING: Found only " << topjet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
751 >           }
752           }
753 +         */
754  
755           for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
756             Particle subjet_v4;
757 <           subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
758 <           subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
759 <           subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
760 <           subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
761 <           topjet.subjets.push_back(subjet_v4);
757 >
758 >           reco::Candidate const * subjetd =  pat_topjet.daughter(k);
759 >           pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
760 >           if(patsubjetd)
761 >           {
762 >              subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
763 >              subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
764 >              subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
765 >              subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
766 >              topjet.add_subjet(subjet_v4);
767 >              //btag info
768 >              topjet.add_subFlavour(patsubjetd->partonFlavour());
769 >              topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
770 >              if (doTagInfos)
771 >                {
772 >                  //ip tag info
773 >                  reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
774 >                  topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
775 >                  topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
776 >                  topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
777 >                  topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
778 >                  topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
779 >                  topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
780 >                  topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
781 >                  topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
782 >                  topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
783 >                  topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
784 >                  topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
785 >                  topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));    
786 >                  topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
787 >                  topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
788 >                  topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
789 >                  topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
790 >                  topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
791 >                  topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
792 >                  topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
793 >                  topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
794 >                  topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
795 >                  //sv tag info
796 >                  reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
797 >                  topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
798 >                  topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
799 >                  topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
800 >                  topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
801 >                  topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
802 >                  topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
803 >                  topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
804 >                  std::vector<TLorentzVector> vp4; vp4.clear();
805 >                  std::vector<float> vchi2; vchi2.clear();
806 >                  std::vector<float> vndof; vndof.clear();
807 >                  std::vector<float> vchi2ndof; vchi2ndof.clear();
808 >                  std::vector<float> vtrsize; vtrsize.clear();
809 >                  for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
810 >                    {
811 >                      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
812 >                      vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
813 >                      vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());  
814 >                      vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());  
815 >                      vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());  
816 >                      vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());  
817 >                    }
818 >                  topjet.add_subSecondaryVertex(vp4);
819 >                  topjet.add_subVertexChi2(vchi2);
820 >                  topjet.add_subVertexNdof(vndof);
821 >                  topjet.add_subVertexNormalizedChi2(vchi2ndof);
822 >                  topjet.add_subVertexTracksSize(vtrsize);
823 >                  //try computer
824 >                  edm::ESHandle<JetTagComputer> computerHandle;
825 >                  iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
826 >                  const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
827 >                  if(computer)
828 >                    {
829 >                      computer->passEventSetup(iSetup);
830 >                      std::vector<const reco::BaseTagInfo*>  baseTagInfos;
831 >                      baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
832 >                      baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );      
833 >                      JetTagComputer::TagInfoHelper helper(baseTagInfos);
834 >                      reco::TaggingVariableList vars = computer->taggingVariables(helper);
835 >                      topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
836 >                      topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
837 >                      topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
838 >                      topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
839 >                    }
840 >                }
841 >           }
842 >           else
843 >             {
844 >               //filling only standard information in case the subjet has not been pat-tified during the pattuples production
845 >               subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
846 >               subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
847 >               subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
848 >               subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
849 >               topjet.add_subjet(subjet_v4);
850 >             }
851 >          
852 >          
853           }
854           topjets[j].push_back(topjet);
855         }
856       }
857     }
858 <
859 <
858 >  
859 >  
860     // ------------- generator top jets -------------
861     if(doGenTopJets){
862       for(size_t j=0; j< gentopjet_sources.size(); ++j){
863        
864         gentopjets[j].clear();
865 <
865 >      
866         edm::Handle<reco::BasicJetCollection> reco_gentopjets;
867         //edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
868         iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
# Line 586 | Line 874 | NtupleWriter::analyze(const edm::Event&
874           if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
875  
876           TopJet gentopjet;
877 <         gentopjet.charge = reco_gentopjet.charge();
878 <         gentopjet.pt = reco_gentopjet.pt();
879 <         gentopjet.eta = reco_gentopjet.eta();
880 <         gentopjet.phi = reco_gentopjet.phi();
881 <         gentopjet.energy = reco_gentopjet.energy();
882 <         gentopjet.numberOfDaughters =reco_gentopjet.numberOfDaughters();
877 >         gentopjet.set_charge(reco_gentopjet.charge());
878 >         gentopjet.set_pt(reco_gentopjet.pt());
879 >         gentopjet.set_eta(reco_gentopjet.eta());
880 >         gentopjet.set_phi(reco_gentopjet.phi());
881 >         gentopjet.set_energy(reco_gentopjet.energy());
882 >         gentopjet.set_numberOfDaughters(reco_gentopjet.numberOfDaughters());
883  
884           for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
885             Particle subjet_v4;
886 <           subjet_v4.pt = reco_gentopjet.daughter(k)->p4().pt();
887 <           subjet_v4.eta = reco_gentopjet.daughter(k)->p4().eta();
888 <           subjet_v4.phi = reco_gentopjet.daughter(k)->p4().phi();
889 <           subjet_v4.energy = reco_gentopjet.daughter(k)->p4().E();
890 <           gentopjet.subjets.push_back(subjet_v4);
886 >           subjet_v4.set_pt(reco_gentopjet.daughter(k)->p4().pt());
887 >           subjet_v4.set_eta(reco_gentopjet.daughter(k)->p4().eta());
888 >           subjet_v4.set_phi(reco_gentopjet.daughter(k)->p4().phi());
889 >           subjet_v4.set_energy(reco_gentopjet.daughter(k)->p4().E());
890 >           gentopjet.add_subjet(subjet_v4);
891           }
892           gentopjets[j].push_back(gentopjet);
893         }
# Line 618 | Line 906 | NtupleWriter::analyze(const edm::Event&
906         for (unsigned int i = 0; i < pat_photons.size(); ++i) {
907           pat::Photon pat_photon = pat_photons[i];
908           Photon ph;
909 <         ph.charge = 0;
910 <         ph.pt =  pat_photon.pt();
911 <         ph.eta =  pat_photon.eta();
912 <         ph.phi =  pat_photon.phi();
913 <         ph.energy =  pat_photon.energy();
914 <         ph.vertex_x = pat_photon.vertex().x();
915 <         ph.vertex_y = pat_photon.vertex().y();
916 <         ph.vertex_z = pat_photon.vertex().z();
917 <         ph.supercluster_eta = pat_photon.superCluster()->eta();
918 <         ph.supercluster_phi = pat_photon.superCluster()->phi();
919 < //       ph.neutralHadronIso = pat_photon.neutralHadronIso();
920 < //       ph.chargedHadronIso = pat_photon.chargedHadronIso();
921 <         ph.trackIso = pat_photon.trackIso();
909 >         ph.set_charge(0);
910 >         ph.set_pt( pat_photon.pt());
911 >         ph.set_eta( pat_photon.eta());
912 >         ph.set_phi( pat_photon.phi());
913 >         ph.set_energy( pat_photon.energy());
914 >         ph.set_vertex_x(pat_photon.vertex().x());
915 >         ph.set_vertex_y(pat_photon.vertex().y());
916 >         ph.set_vertex_z(pat_photon.vertex().z());
917 >         ph.set_supercluster_eta(pat_photon.superCluster()->eta());
918 >         ph.set_supercluster_phi(pat_photon.superCluster()->phi());
919 > //       ph.set_neutralHadronIso(pat_photon.neutralHadronIso());
920 > //       ph.set_chargedHadronIso(pat_photon.chargedHadronIso());
921 >         ph.set_trackIso(pat_photon.trackIso());
922           phs[j].push_back(ph);
923         }
924       }
# Line 650 | Line 938 | NtupleWriter::analyze(const edm::Event&
938         else{
939           pat::MET pat_met = pat_mets[0];
940          
941 <         met[j].pt=pat_met.pt();
942 <         met[j].phi=pat_met.phi();
943 <         met[j].mEtSig= pat_met.mEtSig();
941 >         met[j].set_pt(pat_met.pt());
942 >         met[j].set_phi(pat_met.phi());
943 >         met[j].set_mEtSig(pat_met.mEtSig());
944         }
945        
946       }
# Line 675 | Line 963 | NtupleWriter::analyze(const edm::Event&
963      
964       triggerResults.clear();
965       triggerNames.clear();
966 <     L1_prescale.clear();
967 <     HLT_prescale.clear();
966 > //      L1_prescale.clear();
967 > //      HLT_prescale.clear();
968      
969       edm::Service<edm::service::TriggerNamesService> tns;
970       std::vector<std::string> triggerNames_all;
# Line 693 | Line 981 | NtupleWriter::analyze(const edm::Event&
981         //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
982         triggerResults.push_back(trig.accept(i));
983         if(newrun) triggerNames.push_back(triggerNames_all[i]);
984 <       if(isRealData){
985 <         std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
986 <         L1_prescale.push_back(pre.first);
987 <         HLT_prescale.push_back(pre.second);
988 <         //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
989 <       }
984 > //        if(isRealData){
985 > //       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
986 > //       L1_prescale.push_back(pre.first);
987 > //       HLT_prescale.push_back(pre.second);
988 > //       //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
989 > //        }
990       }
991       //    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
992       //      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
# Line 706 | Line 994 | NtupleWriter::analyze(const edm::Event&
994       newrun=false;
995     }
996  
709   // ------------- generator info -------------
710  
711   if(doGenInfo){
712     genInfo.weights.clear();
713     genInfo.binningValues.clear();
714     genps.clear();
715
716     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
717     iEvent.getByLabel("generator", genEventInfoProduct);
718     const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
719  
720     genInfo.binningValues = genEventInfo.binningValues();
721     genInfo.weights = genEventInfo.weights();
722     genInfo.alphaQCD = genEventInfo.alphaQCD();
723     genInfo.alphaQED = genEventInfo.alphaQED();
724     genInfo.qScale = genEventInfo.qScale();
725    
726     const gen::PdfInfo* pdf = genEventInfo.pdf();
727     if(pdf){
728       genInfo.pdf_id1=pdf->id.first;
729       genInfo.pdf_id2=pdf->id.second;
730       genInfo.pdf_x1=pdf->x.first;
731       genInfo.pdf_x2=pdf->x.second;
732       genInfo.pdf_scalePDF=pdf->scalePDF;
733       genInfo.pdf_xPDF1=pdf->xPDF.first;
734       genInfo.pdf_xPDF2=pdf->xPDF.second;
735     }
736     else{
737       genInfo.pdf_id1=-999;
738       genInfo.pdf_id2=-999;
739       genInfo.pdf_x1=-999;
740       genInfo.pdf_x2=-999;
741       genInfo.pdf_scalePDF=-999;
742       genInfo.pdf_xPDF1=-999;
743       genInfo.pdf_xPDF2=-999;
744     }
745
746     edm::Handle<std::vector<PileupSummaryInfo> > pus;
747     iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
748     genInfo.pileup_NumInteractions_intime=0;
749     genInfo.pileup_NumInteractions_ootbefore=0;
750     genInfo.pileup_NumInteractions_ootafter=0;
751     if(pus.isValid()){
752       genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
753       for(size_t i=0; i<pus->size(); ++i){
754         if(pus->at(i).getBunchCrossing() == 0) // intime pileup
755            genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
756         else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
757            genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
758         }
759         else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
760           genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
761         }
762       }
763     }
764
765     edm::Handle<reco::GenParticleCollection> genPartColl;
766     iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
767     int index=-1;
768     for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
769       index++;
770      
771       //write out only top quarks and status 3 particles (works fine only for MadGraph)
772       if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
773         GenParticle genp;
774         genp.charge = iter->charge();;
775         genp.pt = iter->p4().pt();
776         genp.eta = iter->p4().eta();
777         genp.phi = iter->p4().phi();
778         genp.energy = iter->p4().E();
779         genp.index =index;
780         genp.status = iter->status();
781         genp.pdgId = iter->pdgId();
782
783         genp.mother1=-1;
784         genp.mother2=-1;
785         genp.daughter1=-1;
786         genp.daughter2=-1;
787        
788         int nm=iter->numberOfMothers();
789         int nd=iter->numberOfDaughters();
790
791        
792         if (nm>0) genp.mother1 = iter->motherRef(0).key();
793         if (nm>1) genp.mother2 = iter->motherRef(1).key();
794         if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
795         if (nd>1) genp.daughter2 = iter->daughterRef(1).key();
796         //std::cout << genp.index <<"  pdgId = " << genp.pdgId << "  mo1 = " << genp.mother1 << "  mo2 = " << genp.mother2 <<"  da1 = " << genp.daughter1 << "  da2 = " << genp.daughter2 <<std::endl;
797         genps.push_back(genp);
798       }
799     }
800
801   }
997  
998     tr->Fill();
999     if(doLumiInfo)
# Line 831 | Line 1026 | void
1026   NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
1027   {
1028    if(doTrigger){
1029 <    bool setup_changed = false;
1030 <    hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1029 >    //bool setup_changed = false;
1030 >    //hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1031      newrun=true;
1032    }
1033  
1034    if(doJets || doTopJets){
1035 <    edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
1036 <    iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
1037 <    JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
1038 <    jecUnc = new JetCorrectionUncertainty(JetCorPar);
1035 >    if(doJECUncertainty){
1036 >      edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
1037 >      iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
1038 >      JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
1039 >      jecUnc = new JetCorrectionUncertainty(JetCorPar);
1040 >    }
1041    }
1042   }
1043  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines