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.5 by peiffer, Thu Apr 5 09:48:20 2012 UTC vs.
Revision 1.21 by peiffer, Tue Jun 26 08:13:28 2012 UTC

# Line 50 | Line 50 | NtupleWriter::NtupleWriter(const edm::Pa
50    doMuons = iConfig.getParameter<bool>("doMuons");
51    doTaus = iConfig.getParameter<bool>("doTaus");
52    doJets = iConfig.getParameter<bool>("doJets");
53 +  doJECUncertainty = iConfig.getParameter<bool>("doJECUncertainty");
54 +  doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");  
55    doPhotons = iConfig.getParameter<bool>("doPhotons");
56    doMET = iConfig.getParameter<bool>("doMET");
57    doGenInfo = iConfig.getParameter<bool>("doGenInfo");
58 +  doAllGenParticles = iConfig.getParameter<bool>("doAllGenParticles");
59 +  doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
60    doPV = iConfig.getParameter<bool>("doPV");
61    doTopJets = iConfig.getParameter<bool>("doTopJets");
62 +  doTrigger = iConfig.getParameter<bool>("doTrigger");
63  
64    // initialization of tree variables
65  
# Line 62 | Line 67 | NtupleWriter::NtupleWriter(const edm::Pa
67    tr->Branch("event",&event);
68    tr->Branch("luminosityBlock",&luminosityBlock);
69    tr->Branch("isRealData",&isRealData);
70 <  tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
71 <  tr->Branch("beamspot_x0",&beamspot_x0);
67 <  tr->Branch("beamspot_y0",&beamspot_y0);
68 <  tr->Branch("beamspot_z0",&beamspot_z0);
70 >  tr->Branch("rho",&rho);
71 >  rho_source = iConfig.getParameter<edm::InputTag>("rho_source");
72  
73 +  //tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
74 +  if(doLumiInfo){
75 +    tr->Branch("intgRecLumi",&intgRecLumi);
76 +    tr->Branch("intgDelLumi",&intgDelLumi);
77 +  }
78 +  if(doPV){
79 +    tr->Branch("beamspot_x0",&beamspot_x0);
80 +    tr->Branch("beamspot_y0",&beamspot_y0);
81 +    tr->Branch("beamspot_z0",&beamspot_z0);
82 +  }
83    if(doElectrons){
84      electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
85      for(size_t j=0; j< electron_sources.size(); ++j){  
# Line 103 | Line 116 | NtupleWriter::NtupleWriter(const edm::Pa
116        tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
117      }
118    }
119 +  if(doGenTopJets){
120 +    gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
121 +    gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
122 +    gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
123 +    for(size_t j=0; j< gentopjet_sources.size(); ++j){  
124 +      tr->Branch( gentopjet_sources[j].c_str(), "std::vector<TopJet>", &gentopjets[j]);
125 +    }
126 +  }
127    if(doPhotons){
128      photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
129      for(size_t j=0; j< photon_sources.size(); ++j){  
# Line 125 | Line 146 | NtupleWriter::NtupleWriter(const edm::Pa
146      tr->Branch("genInfo","GenInfo",&genInfo);
147      tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
148    }
149 <
150 <  trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
151 <  //tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
152 <  tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);  
153 <  tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
154 <  tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
155 <  tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
156 <
149 >  if(doTrigger){
150 >    trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
151 >    //tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
152 >    tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);  
153 >    tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
154 >    //tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
155 >    //tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
156 >  }
157    newrun = true;
158   }
159  
# Line 164 | Line 185 | NtupleWriter::analyze(const edm::Event&
185     luminosityBlock = iEvent.luminosityBlock();
186     isRealData      = iEvent.isRealData();
187  
188 <   if(isRealData){
189 <     edm::Handle<bool> bool_handle;
190 <     iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
191 <     HBHENoiseFilterResult = *(bool_handle.product());
192 <   }
193 <   else HBHENoiseFilterResult = false;
188 >   edm::Handle<double> m_rho;
189 >   iEvent.getByLabel(rho_source,m_rho);
190 >   rho=*m_rho;
191 >
192 > //    if(isRealData){
193 > //      edm::Handle<bool> bool_handle;
194 > //      iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
195 > //      HBHENoiseFilterResult = *(bool_handle.product());
196 > //    }
197 > //    else HBHENoiseFilterResult = false;
198  
199     // ------------- primary vertices and beamspot  -------------
200  
# Line 185 | Line 210 | NtupleWriter::analyze(const edm::Event&
210           reco::Vertex reco_pv = reco_pvs[i];
211  
212           PrimaryVertex pv;
213 <         pv.x =  reco_pv.x();
214 <         pv.y =  reco_pv.y();
215 <         pv.z =  reco_pv.z();
216 <         pv.nTracks =  reco_pv.nTracks();
217 <         //pv.isValid =  reco_pv.isValid();
218 <         pv.chi2 =  reco_pv.chi2();
219 <         pv.ndof =  reco_pv.ndof();      
213 >         pv.set_x( reco_pv.x());
214 >         pv.set_y( reco_pv.y());
215 >         pv.set_z( reco_pv.z());
216 >         pv.set_nTracks( reco_pv.nTracks());
217 >         //pv.set_isValid( reco_pv.isValid());
218 >         pv.set_chi2( reco_pv.chi2());
219 >         pv.set_ndof( reco_pv.ndof());  
220  
221           pvs[j].push_back(pv);
222         }
223       }
224 +      
225 +     edm::Handle<reco::BeamSpot> beamSpot;
226 +     iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
227 +     const reco::BeamSpot & bsp = *beamSpot;
228 +    
229 +     beamspot_x0 = bsp.x0();
230 +     beamspot_y0 = bsp.y0();
231 +     beamspot_z0 = bsp.z0();
232     }
233 +
234 + // ------------- generator info -------------
235    
236 <   edm::Handle<reco::BeamSpot> beamSpot;
237 <   iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
238 <   const reco::BeamSpot & bsp = *beamSpot;
239 <  
240 <   beamspot_x0 = bsp.x0();
241 <   beamspot_y0 = bsp.y0();
242 <   beamspot_z0 = bsp.z0();
236 >   if(doGenInfo){
237 >     genInfo.clear_weights();
238 >     genInfo.clear_binningValues();
239 >     genps.clear();
240 >
241 >     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
242 >     iEvent.getByLabel("generator", genEventInfoProduct);
243 >     const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
244 >  
245 >     for(unsigned int k=0; k<genEventInfo.binningValues().size();++k){
246 >       genInfo.add_binningValue(genEventInfo.binningValues().at(k));
247 >     }
248 >     for(unsigned int k=0; k<genEventInfo.weights().size();++k){
249 >       genInfo.add_weight(genEventInfo.weights().at(k));
250 >     }
251 >     genInfo.set_alphaQCD(genEventInfo.alphaQCD());
252 >     genInfo.set_alphaQED(genEventInfo.alphaQED());
253 >     genInfo.set_qScale(genEventInfo.qScale());
254 >    
255 >     const gen::PdfInfo* pdf = genEventInfo.pdf();
256 >     if(pdf){
257 >       genInfo.set_pdf_id1(pdf->id.first);
258 >       genInfo.set_pdf_id2(pdf->id.second);
259 >       genInfo.set_pdf_x1(pdf->x.first);
260 >       genInfo.set_pdf_x2(pdf->x.second);
261 >       genInfo.set_pdf_scalePDF(pdf->scalePDF);
262 >       genInfo.set_pdf_xPDF1(pdf->xPDF.first);
263 >       genInfo.set_pdf_xPDF2(pdf->xPDF.second);
264 >     }
265 >     else{
266 >       genInfo.set_pdf_id1(-999);
267 >       genInfo.set_pdf_id2(-999);
268 >       genInfo.set_pdf_x1(-999);
269 >       genInfo.set_pdf_x2(-999);
270 >       genInfo.set_pdf_scalePDF(-999);
271 >       genInfo.set_pdf_xPDF1(-999);
272 >       genInfo.set_pdf_xPDF2(-999);
273 >     }
274 >
275 >     edm::Handle<std::vector<PileupSummaryInfo> > pus;
276 >     iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
277 >     genInfo.set_pileup_NumInteractions_intime(0);
278 >     genInfo.set_pileup_NumInteractions_ootbefore(0);
279 >     genInfo.set_pileup_NumInteractions_ootafter(0);
280 >     if(pus.isValid()){
281 >       genInfo.set_pileup_TrueNumInteractions ( (float) pus->at(0).getTrueNumInteractions());
282 >       for(size_t i=0; i<pus->size(); ++i){
283 >         if(pus->at(i).getBunchCrossing() == 0) // intime pileup
284 >           genInfo.set_pileup_NumInteractions_intime( genInfo.pileup_NumInteractions_intime() + pus->at(i).getPU_NumInteractions());
285 >         else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
286 >           genInfo.set_pileup_NumInteractions_ootbefore( genInfo.pileup_NumInteractions_ootbefore() + pus->at(i).getPU_NumInteractions());
287 >         }
288 >         else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
289 >           genInfo.set_pileup_NumInteractions_ootafter( genInfo.pileup_NumInteractions_ootafter() + pus->at(i).getPU_NumInteractions());
290 >         }
291 >       }
292 >     }
293 >
294 >     edm::Handle<reco::GenParticleCollection> genPartColl;
295 >     iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
296 >     int index=-1;
297 >     for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
298 >       index++;
299 >      
300 >       //write out only top quarks and status 3 particles (works fine only for MadGraph)
301 >       if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
302 >         GenParticle genp;
303 >         genp.set_charge(iter->charge());
304 >         genp.set_pt(iter->p4().pt());
305 >         genp.set_eta(iter->p4().eta());
306 >         genp.set_phi(iter->p4().phi());
307 >         genp.set_energy(iter->p4().E());
308 >         genp.set_index(index);
309 >         genp.set_status( iter->status());
310 >         genp.set_pdgId( iter->pdgId());
311 >
312 >         genp.set_mother1(-1);
313 >         genp.set_mother2(-1);
314 >         genp.set_daughter1(-1);
315 >         genp.set_daughter2(-1);
316 >        
317 >         int nm=iter->numberOfMothers();
318 >         int nd=iter->numberOfDaughters();
319 >
320 >        
321 >         if (nm>0) genp.set_mother1( iter->motherRef(0).key());
322 >         if (nm>1) genp.set_mother2( iter->motherRef(1).key());
323 >         if (nd>0) genp.set_daughter1( iter->daughterRef(0).key());
324 >         if (nd>1) genp.set_daughter2( iter->daughterRef(1).key());
325 >
326 >         genps.push_back(genp);
327 >       }
328 >     }
329 >   }
330  
331     // ------------- electrons -------------  
332     if(doElectrons){
333 +
334 + //      edm::Handle<reco::ConversionCollection> hConversions;
335 + //      iEvent.getByLabel("allConversions", hConversions);
336 +
337 + //      edm::Handle<reco::BeamSpot> beamSpot;
338 + //      iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
339 + //      const reco::BeamSpot & bsp = *beamSpot;
340 +
341       for(size_t j=0; j< electron_sources.size(); ++j){
342         eles[j].clear();
343         edm::Handle< std::vector<pat::Electron> > ele_handle;
# Line 218 | Line 348 | NtupleWriter::analyze(const edm::Event&
348           pat::Electron pat_ele = pat_electrons[i];
349           Electron ele;
350          
351 <         ele.charge =  pat_ele.charge();
352 <         ele.pt =  pat_ele.pt();
353 <         ele.eta =  pat_ele.eta();
354 <         ele.phi =  pat_ele.phi();
355 <         ele.energy =  pat_ele.energy();
356 <         ele.vertex_x = pat_ele.vertex().x();
357 <         ele.vertex_y = pat_ele.vertex().y();
358 <         ele.vertex_z = pat_ele.vertex().z();
359 <         ele.supercluster_eta = pat_ele.superCluster()->eta();
360 <         ele.supercluster_phi = pat_ele.superCluster()->phi();
361 <         ele.dB = pat_ele.dB();
362 <         //ele.particleIso = pat_ele.particleIso();
363 <         ele.neutralHadronIso = pat_ele.neutralHadronIso();
364 <         ele.chargedHadronIso = pat_ele.chargedHadronIso();
365 <         ele.trackIso = pat_ele.trackIso();
366 <         ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
367 <         ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
368 <         ele.gsfTrack_px= pat_ele.gsfTrack()->px();
369 <         ele.gsfTrack_py= pat_ele.gsfTrack()->py();
370 <         ele.gsfTrack_pz= pat_ele.gsfTrack()->pz();
371 <         ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
372 <         ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
373 <         ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
351 >         ele.set_charge( pat_ele.charge());
352 >         ele.set_pt( pat_ele.pt());
353 >         ele.set_eta( pat_ele.eta());
354 >         ele.set_phi( pat_ele.phi());
355 >         ele.set_energy( pat_ele.energy());
356 >         ele.set_vertex_x(pat_ele.vertex().x());
357 >         ele.set_vertex_y(pat_ele.vertex().y());
358 >         ele.set_vertex_z(pat_ele.vertex().z());
359 >         ele.set_supercluster_eta(pat_ele.superCluster()->eta());
360 >         ele.set_supercluster_phi(pat_ele.superCluster()->phi());
361 >         ele.set_dB(pat_ele.dB());
362 >         //ele.set_particleIso(pat_ele.particleIso());
363 >         ele.set_neutralHadronIso(pat_ele.neutralHadronIso());
364 >         ele.set_chargedHadronIso(pat_ele.chargedHadronIso());
365 >         ele.set_trackIso(pat_ele.trackIso());
366 >         ele.set_photonIso(pat_ele.photonIso());
367 >         ele.set_puChargedHadronIso(pat_ele.puChargedHadronIso());
368 >         ele.set_gsfTrack_trackerExpectedHitsInner_numberOfLostHits(pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits());
369 >         ele.set_gsfTrack_px( pat_ele.gsfTrack()->px());
370 >         ele.set_gsfTrack_py( pat_ele.gsfTrack()->py());
371 >         ele.set_gsfTrack_pz( pat_ele.gsfTrack()->pz());
372 >         ele.set_gsfTrack_vx( pat_ele.gsfTrack()->vx());
373 >         ele.set_gsfTrack_vy( pat_ele.gsfTrack()->vy());
374 >         ele.set_gsfTrack_vz( pat_ele.gsfTrack()->vz());
375 >         //ele.set_passconversionveto(!ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position()));
376 >         ele.set_passconversionveto(pat_ele.passConversionVeto());
377 >         ele.set_dEtaIn(pat_ele.deltaEtaSuperClusterTrackAtVtx());
378 >         ele.set_dPhiIn(pat_ele.deltaPhiSuperClusterTrackAtVtx());
379 >         ele.set_sigmaIEtaIEta(pat_ele.sigmaIetaIeta());
380 >         ele.set_HoverE(pat_ele.hadronicOverEm());
381 >         ele.set_fbrem(pat_ele.fbrem());
382 >         ele.set_EoverPIn(pat_ele.eSuperClusterOverP());
383 >         ele.set_EcalEnergy(pat_ele.ecalEnergy());
384 >         ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
385 >         ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
386 >
387           eles[j].push_back(ele);
388         }
389       }
# Line 259 | Line 402 | NtupleWriter::analyze(const edm::Event&
402           pat::Muon pat_mu = pat_muons[i];
403  
404           Muon mu;
405 <         mu.charge =  pat_mu.charge();
406 <         mu.pt =  pat_mu.pt();
407 <         mu.eta =  pat_mu.eta();
408 <         mu.phi =  pat_mu.phi();
409 <         mu.energy =  pat_mu.energy();
410 <         mu.vertex_x = pat_mu.vertex().x();
411 <         mu.vertex_y = pat_mu.vertex().y();
412 <         mu.vertex_z = pat_mu.vertex().z();
413 <         mu.dB = pat_mu.dB();
414 <         //mu.particleIso = pat_mu.particleIso();
415 <         mu.neutralHadronIso = pat_mu.neutralHadronIso();
416 <         mu.chargedHadronIso = pat_mu.chargedHadronIso();
417 <         mu.trackIso = pat_mu.trackIso();
418 <         mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
419 <         mu.isGlobalMuon = pat_mu.isGlobalMuon();
420 <         mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
421 <         mu.isTrackerMuon = pat_mu.isTrackerMuon();
422 <         mu.numberOfMatchedStations = pat_mu.numberOfMatchedStations();
405 >         mu.set_charge( pat_mu.charge());
406 >         mu.set_pt( pat_mu.pt());
407 >         mu.set_eta( pat_mu.eta());
408 >         mu.set_phi( pat_mu.phi());
409 >         mu.set_energy( pat_mu.energy());
410 >         mu.set_vertex_x ( pat_mu.vertex().x());
411 >         mu.set_vertex_y ( pat_mu.vertex().y());
412 >         mu.set_vertex_z ( pat_mu.vertex().z());
413 >         mu.set_dB ( pat_mu.dB());
414 >         //mu.particleIso ( pat_mu.particleIso());
415 >         mu.set_neutralHadronIso ( pat_mu.neutralHadronIso());
416 >         mu.set_chargedHadronIso ( pat_mu.chargedHadronIso());
417 >         mu.set_trackIso ( pat_mu.trackIso());
418 >         mu.set_photonIso ( pat_mu.photonIso());
419 >         mu.set_puChargedHadronIso ( pat_mu.puChargedHadronIso());
420 >         mu.set_isGlobalMuon ( pat_mu.isGlobalMuon());
421 >         mu.set_isPFMuon ( pat_mu.isPFMuon());
422 >         mu.set_isStandAloneMuon ( pat_mu.isStandAloneMuon());
423 >         mu.set_isTrackerMuon ( pat_mu.isTrackerMuon());
424 >         mu.set_numberOfMatchedStations ( pat_mu.numberOfMatchedStations());
425           reco::TrackRef globalTrack = pat_mu.globalTrack();
426           if(!globalTrack.isNull()){
427 <           mu.globalTrack_chi2 = globalTrack->chi2();
428 <           mu.globalTrack_ndof = globalTrack->ndof();
429 <           mu.globalTrack_d0 = globalTrack->d0();        
430 <           mu.globalTrack_d0Error = globalTrack->d0Error();
431 <           mu.globalTrack_numberOfValidHits = globalTrack->numberOfValidHits();
432 <           mu.globalTrack_numberOfLostHits = globalTrack->numberOfLostHits();
427 >           mu.set_globalTrack_chi2 ( globalTrack->chi2());
428 >           mu.set_globalTrack_ndof ( globalTrack->ndof());
429 >           mu.set_globalTrack_d0 ( globalTrack->d0());  
430 >           mu.set_globalTrack_d0Error ( globalTrack->d0Error());
431 >           mu.set_globalTrack_numberOfValidHits ( globalTrack->numberOfValidHits());
432 >           mu.set_globalTrack_numberOfLostHits ( globalTrack->numberOfLostHits());
433 >           mu.set_globalTrack_numberOfValidMuonHits(globalTrack->hitPattern().numberOfValidMuonHits() );
434           }
435           else{
436 <           mu.globalTrack_chi2 = 0;
437 <           mu.globalTrack_ndof = 0;
438 <           mu.globalTrack_d0 = 0;
439 <           mu.globalTrack_d0Error = 0;
440 <           mu.globalTrack_numberOfValidHits = 0;
441 <           mu.globalTrack_numberOfLostHits = 0;
436 >           mu.set_globalTrack_chi2 ( 0);
437 >           mu.set_globalTrack_ndof ( 0);
438 >           mu.set_globalTrack_d0 ( 0);
439 >           mu.set_globalTrack_d0Error ( 0);
440 >           mu.set_globalTrack_numberOfValidHits ( 0);
441 >           mu.set_globalTrack_numberOfLostHits ( 0);
442           }
443           reco::TrackRef innerTrack = pat_mu.innerTrack();
444           if(!innerTrack.isNull()){
445 <           mu.innerTrack_chi2 = innerTrack->chi2();
446 <           mu.innerTrack_ndof = innerTrack->ndof();
447 <           mu.innerTrack_d0 = innerTrack->d0();  
448 <           mu.innerTrack_d0Error = innerTrack->d0Error();
449 <           mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
450 <           mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
445 >           mu.set_innerTrack_chi2 ( innerTrack->chi2());
446 >           mu.set_innerTrack_ndof ( innerTrack->ndof());
447 >           mu.set_innerTrack_d0 ( innerTrack->d0());    
448 >           mu.set_innerTrack_d0Error ( innerTrack->d0Error());
449 >           mu.set_innerTrack_numberOfValidHits ( innerTrack->numberOfValidHits());
450 >           mu.set_innerTrack_numberOfLostHits ( innerTrack->numberOfLostHits());
451 >           mu.set_innerTrack_trackerLayersWithMeasurement ( innerTrack->hitPattern().trackerLayersWithMeasurement());
452 >           mu.set_innerTrack_numberOfValidPixelHits ( innerTrack->hitPattern().numberOfValidPixelHits());
453           }
454           else{
455 <           mu.innerTrack_chi2 = 0;
456 <           mu.innerTrack_ndof = 0;
457 <           mu.innerTrack_d0 = 0;
458 <           mu.innerTrack_d0Error = 0;
459 <           mu.innerTrack_numberOfValidHits = 0;
460 <           mu.innerTrack_numberOfLostHits = 0;
455 >           mu.set_innerTrack_chi2 ( 0);
456 >           mu.set_innerTrack_ndof ( 0);
457 >           mu.set_innerTrack_d0 ( 0);
458 >           mu.set_innerTrack_d0Error ( 0);
459 >           mu.set_innerTrack_numberOfValidHits ( 0);
460 >           mu.set_innerTrack_numberOfLostHits ( 0);
461 >           mu.set_innerTrack_trackerLayersWithMeasurement ( 0);
462 >           mu.set_innerTrack_numberOfValidPixelHits ( 0);
463           }
464           reco::TrackRef outerTrack = pat_mu.outerTrack();
465           if(!outerTrack.isNull()){
466 <           mu.outerTrack_chi2 = outerTrack->chi2();
467 <           mu.outerTrack_ndof = outerTrack->ndof();
468 <           mu.outerTrack_d0 = outerTrack->d0();  
469 <           mu.outerTrack_d0Error = outerTrack->d0Error();
470 <           mu.outerTrack_numberOfValidHits = outerTrack->numberOfValidHits();
471 <           mu.outerTrack_numberOfLostHits = outerTrack->numberOfLostHits();
466 >           mu.set_outerTrack_chi2 ( outerTrack->chi2());
467 >           mu.set_outerTrack_ndof ( outerTrack->ndof());
468 >           mu.set_outerTrack_d0 ( outerTrack->d0());    
469 >           mu.set_outerTrack_d0Error ( outerTrack->d0Error());
470 >           mu.set_outerTrack_numberOfValidHits ( outerTrack->numberOfValidHits());
471 >           mu.set_outerTrack_numberOfLostHits ( outerTrack->numberOfLostHits());
472           }
473           else{
474 <           mu.outerTrack_chi2 = 0;
475 <           mu.outerTrack_ndof = 0;
476 <           mu.outerTrack_d0 = 0;
477 <           mu.outerTrack_d0Error = 0;
478 <           mu.outerTrack_numberOfValidHits = 0;
479 <           mu.outerTrack_numberOfLostHits = 0;
474 >           mu.set_outerTrack_chi2 ( 0);
475 >           mu.set_outerTrack_ndof ( 0);
476 >           mu.set_outerTrack_d0 ( 0);
477 >           mu.set_outerTrack_d0Error ( 0);
478 >           mu.set_outerTrack_numberOfValidHits ( 0);
479 >           mu.set_outerTrack_numberOfLostHits ( 0);
480           }
481  
482           mus[j].push_back(mu);
# Line 350 | Line 500 | NtupleWriter::analyze(const edm::Event&
500           if(fabs(pat_tau.eta()) > tau_etamax) continue;
501  
502           Tau tau;
503 <         tau.charge =  pat_tau.charge();
504 <         tau.pt =  pat_tau.pt();
505 <         tau.eta =  pat_tau.eta();
506 <         tau.phi =  pat_tau.phi();
507 <         tau.energy =  pat_tau.energy();
508 <
503 >         tau.set_charge( pat_tau.charge());
504 >         tau.set_pt( pat_tau.pt());
505 >         tau.set_eta( pat_tau.eta());
506 >         tau.set_phi( pat_tau.phi());
507 >         tau.set_energy( pat_tau.energy());
508 >         tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
509 >         tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
510 >         tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
511 >         tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
512 >         tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
513 >         tau.set_againstElectronLoose  ( pat_tau.tauID("againstElectronLoose")>0.5);
514 >         tau.set_againstElectronMedium ( pat_tau.tauID("againstElectronMedium")>0.5);
515 >         tau.set_againstElectronTight ( pat_tau.tauID("againstElectronTight")>0.5);
516 >         tau.set_againstElectronMVA  ( pat_tau.tauID("againstElectronMVA")>0.5);
517 >         tau.set_againstMuonLoose ( pat_tau.tauID("againstMuonLoose")>0.5);
518 >         tau.set_againstMuonMedium ( pat_tau.tauID("againstMuonMedium")>0.5);
519 >         tau.set_againstMuonTight ( pat_tau.tauID("againstMuonTight")>0.5);
520 >
521 > //       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
522 > //       if(!leadPFCand.isNull()){
523 > //         tau.set_leadPFCand_px ( leadPFCand->px());
524 > //         tau.set_leadPFCand_py ( leadPFCand->py());
525 > //         tau.set_leadPFCand_pz ( leadPFCand->pz());
526 > //       }
527 > //       else{
528 > //         tau.set_leadPFCand_px ( 0);
529 > //         tau.set_leadPFCand_py ( 0);
530 > //         tau.set_leadPFCand_pz ( 0);
531 > //       }
532           taus[j].push_back(tau);
533         }
534       }
# Line 382 | Line 555 | NtupleWriter::analyze(const edm::Event&
555   //       }
556  
557           Jet jet;
558 <         jet.charge = pat_jet.charge();
559 <         jet.pt = pat_jet.pt();
560 <         jet.eta = pat_jet.eta();
561 <         jet.phi = pat_jet.phi();
562 <         jet.energy = pat_jet.energy();
563 <         jet.numberOfDaughters =pat_jet.numberOfDaughters();
558 >         jet.set_charge(pat_jet.charge());
559 >         jet.set_pt(pat_jet.pt());
560 >         jet.set_eta(pat_jet.eta());
561 >         jet.set_phi(pat_jet.phi());
562 >         jet.set_energy(pat_jet.energy());
563 >         jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
564           const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
565 <         jet.nTracks = jettracks.size();
566 <         jet.jetArea = pat_jet.jetArea();
567 <         jet.pileup = pat_jet.pileup();
565 >         jet.set_nTracks ( jettracks.size());
566 >         jet.set_jetArea(pat_jet.jetArea());
567 >         jet.set_pileup(pat_jet.pileup());
568           if(pat_jet.isPFJet()){
569 <           jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
570 <           jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
571 <           jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
572 <           jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
573 <           jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
574 <           jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
575 <           jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
576 <           jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
577 <           jet.muonMultiplicity =pat_jet.muonMultiplicity();
578 <           jet.electronMultiplicity =pat_jet.electronMultiplicity();
579 <           jet.photonMultiplicity =pat_jet.photonMultiplicity();
580 <         }
581 <
582 <         jecUnc->setJetEta(pat_jet.eta());
583 <         jecUnc->setJetPt(pat_jet.pt());
584 <         jet.JEC_uncertainty = jecUnc->getUncertainty(true);
585 <
586 <         jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
587 <         jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
588 <         jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
589 <         jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
590 <         jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
591 <         jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
592 <
569 >           jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
570 >           jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
571 >           jet.set_chargedEmEnergyFraction (pat_jet.chargedEmEnergyFraction());
572 >           jet.set_chargedHadronEnergyFraction (pat_jet.chargedHadronEnergyFraction());
573 >           jet.set_muonEnergyFraction (pat_jet.muonEnergyFraction());
574 >           jet.set_photonEnergyFraction (pat_jet.photonEnergyFraction());
575 >           jet.set_chargedMultiplicity (pat_jet.chargedMultiplicity());
576 >           jet.set_neutralMultiplicity (pat_jet.neutralMultiplicity());
577 >           jet.set_muonMultiplicity (pat_jet.muonMultiplicity());
578 >           jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
579 >           jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
580 >         }
581 >         if(doJECUncertainty){
582 >           jecUnc->setJetEta(pat_jet.eta());
583 >           jecUnc->setJetPt(pat_jet.pt());
584 >           jet.set_JEC_uncertainty(jecUnc->getUncertainty(true));
585 >         }
586 >         jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
587 >        
588 >         jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
589 >         jet.set_btag_simpleSecondaryVertexHighPur(pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
590 >         jet.set_btag_combinedSecondaryVertex(pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags"));
591 >         jet.set_btag_combinedSecondaryVertexMVA(pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
592 >         jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
593 >         jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
594 >
595 >         const reco::GenJet *genj = pat_jet.genJet();
596 >         if(genj){
597 >           jet.set_genjet_pt(genj->pt());
598 >           jet.set_genjet_eta(genj->eta());
599 >           jet.set_genjet_phi(genj->phi());
600 >           jet.set_genjet_energy(genj->energy());
601 >           if(doAllGenParticles){
602 >             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
603 >             for(unsigned int l = 0; l<jetgenps.size(); ++l){
604 >               for(unsigned int k=0; k< genps.size(); ++k){
605 >                 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
606 >                   jet.add_genparticles_index(genps[k].index());
607 >                   break;
608 >                 }
609 >               }
610 >             }
611 >             if(jet.genparticles_indices().size()!= jetgenps.size())
612 >               std::cout << "WARNING: Found only " << jet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
613 >           }
614 >         }
615           jets[j].push_back(jet);
616         }
617       }
# Line 433 | Line 628 | NtupleWriter::analyze(const edm::Event&
628         iEvent.getByLabel(topjet_sources[j], pat_topjets);
629  
630         for (unsigned int i = 0; i < pat_topjets->size(); i++) {
631 +
632           const pat::Jet  pat_topjet =  * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
633           if(pat_topjet.pt() < topjet_ptmin) continue;
634           if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
635  
636           TopJet topjet;
637 <         topjet.charge = pat_topjet.charge();
638 <         topjet.pt = pat_topjet.pt();
639 <         topjet.eta = pat_topjet.eta();
640 <         topjet.phi = pat_topjet.phi();
641 <         topjet.energy = pat_topjet.energy();
642 <         topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
637 >         topjet.set_charge(pat_topjet.charge());
638 >         topjet.set_pt(pat_topjet.pt());
639 >         topjet.set_eta(pat_topjet.eta());
640 >         topjet.set_phi(pat_topjet.phi());
641 >         topjet.set_energy(pat_topjet.energy());
642 >         topjet.set_numberOfDaughters(pat_topjet.numberOfDaughters());
643           const reco::TrackRefVector&  topjettracks = pat_topjet.associatedTracks();
644 <         topjet.nTracks = topjettracks.size();
645 <         topjet.jetArea = pat_topjet.jetArea();
646 <         topjet.pileup = pat_topjet.pileup();
647 < //       topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
648 < //       topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
649 < //       topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
650 < //       topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
651 < //       topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
652 < //       topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
653 < //       topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
654 < //       topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
655 < //       topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
656 < //       topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
657 < //       topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
658 <
659 <         jecUnc->setJetEta(pat_topjet.eta());
660 <         jecUnc->setJetPt(pat_topjet.pt());
661 <         topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
662 <
663 <         topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
664 <         topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
665 <         topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
666 <         topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
667 <         topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
668 <         topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
644 >         topjet.set_nTracks( topjettracks.size());
645 >         topjet.set_jetArea( pat_topjet.jetArea());
646 >         topjet.set_pileup( pat_topjet.pileup());
647 > //       topjet.set_neutralEmEnergyFraction(pat_topjet.neutralEmEnergyFraction());
648 > //       topjet.set_neutralHadronEnergyFraction(pat_topjet.neutralHadronEnergyFraction());
649 > //       topjet.set_chargedEmEnergyFraction(pat_topjet.chargedEmEnergyFraction());
650 > //       topjet.set_chargedHadronEnergyFraction(pat_topjet.chargedHadronEnergyFraction());
651 > //       topjet.set_muonEnergyFraction(pat_topjet.muonEnergyFraction());
652 > //       topjet.set_photonEnergyFraction(pat_topjet.photonEnergyFraction());
653 > //       topjet.set_chargedMultiplicity(pat_topjet.chargedMultiplicity());
654 > //       topjet.set_neutralMultiplicity(pat_topjet.neutralMultiplicity());
655 > //       topjet.set_muonMultiplicity(pat_topjet.muonMultiplicity());
656 > //       topjet.set_electronMultiplicity(pat_topjet.electronMultiplicity());
657 > //       topjet.set_photonMultiplicity(pat_topjet.photonMultiplicity());
658 >         if(doJECUncertainty){
659 >           jecUnc->setJetEta(pat_topjet.eta());
660 >           jecUnc->setJetPt(pat_topjet.pt());
661 >           topjet.set_JEC_uncertainty( jecUnc->getUncertainty(true));
662 >         }
663 >         topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
664 >
665 >         topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
666 >         topjet.set_btag_simpleSecondaryVertexHighPur(pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
667 >         topjet.set_btag_combinedSecondaryVertex(pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags"));
668 >         topjet.set_btag_combinedSecondaryVertexMVA(pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
669 >         topjet.set_btag_jetBProbability(pat_topjet.bDiscriminator("jetBProbabilityBJetTags"));
670 >         topjet.set_btag_jetProbability(pat_topjet.bDiscriminator("jetProbabilityBJetTags"));
671 >
672 >         /*
673 >         const reco::GenJet *genj = pat_topjet.genJet();
674 >         if(genj){
675 >           topjet.set_genjet_pt ( genj->pt());
676 >           topjet.set_genjet_eta ( genj->eta());
677 >           topjet.set_genjet_phi ( genj->phi());
678 >           topjet.set_genjet_energy ( genj->energy());
679 >           if(doAllGenParticles){
680 >             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
681 >             for(unsigned int l = 0; l<jetgenps.size(); ++l){
682 >               for(unsigned int k=0; k< genps.size(); ++k){
683 >                 if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
684 >                   topjet.add_genparticles_index(genps[k].index());
685 >                   break;
686 >                 }
687 >               }
688 >             }
689 >             if(topjet.genparticles_indices().size()!= jetgenps.size())
690 >               std::cout << "WARNING: Found only " << topjet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
691 >           }
692 >         }
693 >         */
694  
695           for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
696             Particle subjet_v4;
697 <           subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
698 <           subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
699 <           subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
700 <           subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
701 <           topjet.subjets.push_back(subjet_v4);
697 >           subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
698 >           subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
699 >           subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
700 >           subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
701 >           topjet.add_subjet(subjet_v4);
702           }
703           topjets[j].push_back(topjet);
704         }
705       }
706     }
707  
708 +
709 +   // ------------- generator top jets -------------
710 +   if(doGenTopJets){
711 +     for(size_t j=0; j< gentopjet_sources.size(); ++j){
712 +      
713 +       gentopjets[j].clear();
714 +
715 +       edm::Handle<reco::BasicJetCollection> reco_gentopjets;
716 +       //edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
717 +       iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
718 +
719 +       for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
720 +        
721 +         const reco::BasicJet  reco_gentopjet =  reco_gentopjets->at(i);
722 +         if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
723 +         if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
724 +
725 +         TopJet gentopjet;
726 +         gentopjet.set_charge(reco_gentopjet.charge());
727 +         gentopjet.set_pt(reco_gentopjet.pt());
728 +         gentopjet.set_eta(reco_gentopjet.eta());
729 +         gentopjet.set_phi(reco_gentopjet.phi());
730 +         gentopjet.set_energy(reco_gentopjet.energy());
731 +         gentopjet.set_numberOfDaughters(reco_gentopjet.numberOfDaughters());
732 +
733 +         for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
734 +           Particle subjet_v4;
735 +           subjet_v4.set_pt(reco_gentopjet.daughter(k)->p4().pt());
736 +           subjet_v4.set_eta(reco_gentopjet.daughter(k)->p4().eta());
737 +           subjet_v4.set_phi(reco_gentopjet.daughter(k)->p4().phi());
738 +           subjet_v4.set_energy(reco_gentopjet.daughter(k)->p4().E());
739 +           gentopjet.add_subjet(subjet_v4);
740 +         }
741 +         gentopjets[j].push_back(gentopjet);
742 +       }
743 +     }
744 +   }
745 +
746     // ------------- photons -------------
747     if(doPhotons){
748       for(size_t j=0; j< photon_sources.size(); ++j){
# Line 496 | Line 755 | NtupleWriter::analyze(const edm::Event&
755         for (unsigned int i = 0; i < pat_photons.size(); ++i) {
756           pat::Photon pat_photon = pat_photons[i];
757           Photon ph;
758 <         ph.charge = 0;
759 <         ph.pt =  pat_photon.pt();
760 <         ph.eta =  pat_photon.eta();
761 <         ph.phi =  pat_photon.phi();
762 <         ph.energy =  pat_photon.energy();
763 <         ph.vertex_x = pat_photon.vertex().x();
764 <         ph.vertex_y = pat_photon.vertex().y();
765 <         ph.vertex_z = pat_photon.vertex().z();
766 <         ph.supercluster_eta = pat_photon.superCluster()->eta();
767 <         ph.supercluster_phi = pat_photon.superCluster()->phi();
768 < //       ph.neutralHadronIso = pat_photon.neutralHadronIso();
769 < //       ph.chargedHadronIso = pat_photon.chargedHadronIso();
770 <         ph.trackIso = pat_photon.trackIso();
758 >         ph.set_charge(0);
759 >         ph.set_pt( pat_photon.pt());
760 >         ph.set_eta( pat_photon.eta());
761 >         ph.set_phi( pat_photon.phi());
762 >         ph.set_energy( pat_photon.energy());
763 >         ph.set_vertex_x(pat_photon.vertex().x());
764 >         ph.set_vertex_y(pat_photon.vertex().y());
765 >         ph.set_vertex_z(pat_photon.vertex().z());
766 >         ph.set_supercluster_eta(pat_photon.superCluster()->eta());
767 >         ph.set_supercluster_phi(pat_photon.superCluster()->phi());
768 > //       ph.set_neutralHadronIso(pat_photon.neutralHadronIso());
769 > //       ph.set_chargedHadronIso(pat_photon.chargedHadronIso());
770 >         ph.set_trackIso(pat_photon.trackIso());
771           phs[j].push_back(ph);
772         }
773       }
# Line 528 | Line 787 | NtupleWriter::analyze(const edm::Event&
787         else{
788           pat::MET pat_met = pat_mets[0];
789          
790 <         met[j].pt=pat_met.pt();
791 <         met[j].phi=pat_met.phi();
792 <         met[j].mEtSig= pat_met.mEtSig();
790 >         met[j].set_pt(pat_met.pt());
791 >         met[j].set_phi(pat_met.phi());
792 >         met[j].set_mEtSig(pat_met.mEtSig());
793         }
794        
795       }
796     }
797  
798     // ------------- trigger -------------
799 <  
800 <   edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
801 <   edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
802 <   iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
544 <  
545 <   const edm::Provenance *meta = dummy_TriggerEvent.provenance();
546 <   std::string nameProcess = meta->processName();
547 <   edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
548 <   triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
549 <  
550 <   edm::Handle<edm::TriggerResults> trigger;
551 <   iEvent.getByLabel(triggerResultTag, trigger);
552 <   const edm::TriggerResults& trig = *(trigger.product());
553 <  
554 <   triggerResults.clear();
555 <   triggerNames.clear();
556 <   L1_prescale.clear();
557 <   HLT_prescale.clear();
558 <
559 <   edm::Service<edm::service::TriggerNamesService> tns;
560 <   std::vector<std::string> triggerNames_all;
561 <   tns->getTrigPaths(trig,triggerNames_all);
562 <
563 <   if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
564 <   for(unsigned int i=0; i<trig.size(); ++i){
565 <     std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
566 <     for(; it!=trigger_prefixes.end(); ++it){
567 <       if(triggerNames_all[i].substr(0, it->size()) == *it)break;
568 <     }
569 <     if(it==trigger_prefixes.end()) continue;
570 <
571 <     //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
572 <     triggerResults.push_back(trig.accept(i));
573 <     if(newrun) triggerNames.push_back(triggerNames_all[i]);
574 <     if(isRealData){
575 <       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
576 <       L1_prescale.push_back(pre.first);
577 <       HLT_prescale.push_back(pre.second);
578 <     }
579 <   }
580 < //    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
581 < //      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
582 < //    }
583 <   newrun=false;
584 <
585 <   // ------------- generator info -------------
586 <  
587 <   if(doGenInfo){
588 <     genInfo.weights.clear();
589 <     genInfo.binningValues.clear();
590 <     genps.clear();
591 <
592 <     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
593 <     iEvent.getByLabel("generator", genEventInfoProduct);
594 <     const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
595 <  
596 <     genInfo.binningValues = genEventInfo.binningValues();
597 <     genInfo.weights = genEventInfo.weights();
598 <     genInfo.alphaQCD = genEventInfo.alphaQCD();
599 <     genInfo.alphaQED = genEventInfo.alphaQED();
600 <     genInfo.qScale = genEventInfo.qScale();
799 >   if(doTrigger){
800 >     edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
801 >     edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
802 >     iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
803      
804 <     const gen::PdfInfo* pdf = genEventInfo.pdf();
805 <     if(pdf){
806 <       genInfo.pdf_id1=pdf->id.first;
807 <       genInfo.pdf_id2=pdf->id.second;
808 <       genInfo.pdf_x1=pdf->x.first;
809 <       genInfo.pdf_x2=pdf->x.second;
810 <       genInfo.pdf_scalePDF=pdf->scalePDF;
811 <       genInfo.pdf_xPDF1=pdf->xPDF.first;
812 <       genInfo.pdf_xPDF2=pdf->xPDF.second;
813 <     }
814 <     else{
815 <       genInfo.pdf_id1=-999;
816 <       genInfo.pdf_id2=-999;
817 <       genInfo.pdf_x1=-999;
818 <       genInfo.pdf_x2=-999;
819 <       genInfo.pdf_scalePDF=-999;
820 <       genInfo.pdf_xPDF1=-999;
821 <       genInfo.pdf_xPDF2=-999;
822 <     }
823 <
824 <     edm::Handle<std::vector<PileupSummaryInfo> > pus;
825 <     iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
826 <     genInfo.pileup_NumInteractions_intime=0;
827 <     genInfo.pileup_NumInteractions_ootbefore=0;
828 <     genInfo.pileup_NumInteractions_ootafter=0;
829 <     if(pus.isValid()){
830 <       genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
831 <       for(size_t i=0; i<pus->size(); ++i){
832 <         if(pus->at(i).getBunchCrossing() == 0) // intime pileup
833 <            genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
834 <         else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
835 <            genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
836 <         }
837 <         else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
838 <           genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
839 <         }
840 <       }
841 <     }
842 <
843 <     edm::Handle<reco::GenParticleCollection> genPartColl;
642 <     iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
643 <     int index=-1;
644 <     for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
645 <       index++;
646 <      
647 <       //write out only top quarks and status 3 particles (works fine only for MadGraph)
648 <       if(abs(iter->pdgId())==6 || iter->status()==3){
649 <         GenParticle genp;
650 <         genp.charge = iter->charge();;
651 <         genp.pt = iter->p4().pt();
652 <         genp.eta = iter->p4().eta();
653 <         genp.phi = iter->p4().phi();
654 <         genp.energy = iter->p4().E();
655 <         genp.index =index;
656 <         genp.status = iter->status();
657 <         genp.pdgId = iter->pdgId();
658 <
659 <         genp.mother1=-1;
660 <         genp.mother2=-1;
661 <         genp.daughter1=-1;
662 <         genp.daughter2=-1;
663 <        
664 <         int nm=iter->numberOfMothers();
665 <         int nd=iter->numberOfDaughters();
666 <        
667 <         if (nm>0) genp.mother1 = iter->motherRef(0).key();
668 <         if (nm>1) genp.mother2 = iter->motherRef(nm-1).key();
669 <         if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
670 <         if (nd>1) genp.daughter2 = iter->daughterRef(nd-1).key();
671 <        
672 <         genps.push_back(genp);
673 <       }
674 <     }
675 <
804 >     const edm::Provenance *meta = dummy_TriggerEvent.provenance();
805 >     std::string nameProcess = meta->processName();
806 >     edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
807 >     triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
808 >    
809 >     edm::Handle<edm::TriggerResults> trigger;
810 >     iEvent.getByLabel(triggerResultTag, trigger);
811 >     const edm::TriggerResults& trig = *(trigger.product());
812 >    
813 >     triggerResults.clear();
814 >     triggerNames.clear();
815 > //      L1_prescale.clear();
816 > //      HLT_prescale.clear();
817 >    
818 >     edm::Service<edm::service::TriggerNamesService> tns;
819 >     std::vector<std::string> triggerNames_all;
820 >     tns->getTrigPaths(trig,triggerNames_all);
821 >    
822 >     if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
823 >     for(unsigned int i=0; i<trig.size(); ++i){
824 >       std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
825 >       for(; it!=trigger_prefixes.end(); ++it){
826 >         if(triggerNames_all[i].substr(0, it->size()) == *it)break;
827 >       }
828 >       if(it==trigger_prefixes.end()) continue;
829 >      
830 >       //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
831 >       triggerResults.push_back(trig.accept(i));
832 >       if(newrun) triggerNames.push_back(triggerNames_all[i]);
833 > //        if(isRealData){
834 > //       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
835 > //       L1_prescale.push_back(pre.first);
836 > //       HLT_prescale.push_back(pre.second);
837 > //       //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
838 > //        }
839 >     }
840 >     //    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
841 >     //      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
842 >     //    }
843 >     newrun=false;
844     }
845  
846  
847     tr->Fill();
848 <
848 >   if(doLumiInfo)
849 >     previouslumiblockwasfilled=true;
850   }
851  
852  
# Line 685 | Line 854 | NtupleWriter::analyze(const edm::Event&
854   void
855   NtupleWriter::beginJob()
856   {
857 +  if(doLumiInfo){
858 +    totalRecLumi=0;
859 +    totalDelLumi=0;
860 +    previouslumiblockwasfilled=false;
861 +  }
862   }
863  
864   // ------------ method called once each job just after ending the event loop  ------------
# Line 700 | Line 874 | NtupleWriter::endJob()
874   void
875   NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
876   {
877 <  bool setup_changed = false;
878 <  hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
879 <  newrun=true;
880 <
881 <  edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
708 <  iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
709 <  JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
710 <  jecUnc = new JetCorrectionUncertainty(JetCorPar);
877 >  if(doTrigger){
878 >    //bool setup_changed = false;
879 >    //hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
880 >    newrun=true;
881 >  }
882  
883 +  if(doJets || doTopJets){
884 +    if(doJECUncertainty){
885 +      edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
886 +      iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
887 +      JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
888 +      jecUnc = new JetCorrectionUncertainty(JetCorPar);
889 +    }
890 +  }
891   }
892  
893   // ------------ method called when ending the processing of a run  ------------
894   void
895   NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
896   {
897 +  if(doLumiInfo)
898 +    std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
899   }
900  
901   // ------------ method called when starting to processes a luminosity block  ------------
902   void
903 < NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
903 > NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
904   {
905 +  if(doLumiInfo){
906 +    edm::Handle<LumiSummary> l;
907 +    lumi.getByLabel("lumiProducer", l);
908 +    
909 +    //add lumi of lumi blocks without any event to next lumiblock
910 +    if(previouslumiblockwasfilled){
911 +      intgRecLumi=0;
912 +      intgDelLumi=0;
913 +    }
914 +    previouslumiblockwasfilled=false;
915 +    
916 +    if (l.isValid()){;
917 +      intgRecLumi+=l->intgRecLumi()*6.37;
918 +      intgDelLumi+=l->intgDelLumi()*6.37;
919 +      totalRecLumi+=l->intgRecLumi()*6.37;
920 +      totalDelLumi+=l->intgDelLumi()*6.37;
921 +    }
922 +    //std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<"   " << l->intgDelLumi()*6.37<<std::endl;
923 +    //std::cout << "summed:  "<< intgRecLumi << "   " << intgDelLumi << std::endl;
924 +  }
925   }
926  
927   // ------------ method called when ending the processing of a luminosity block  ------------

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines