ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc

Comparing UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc (file contents):
Revision 1.16 by peiffer, Mon Apr 30 09:01:16 2012 UTC vs.
Revision 1.29 by rkogler, Wed Jun 19 22:05:46 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 "UHHAnalysis/NtupleWriter/interface/JetProps.h"
22	>	#include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputer.h"
23	>	#include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputerWrapper.h"
24	>	#include "RecoBTau/JetTagComputer/interface/JetTagComputer.h"
25	>	#include "RecoBTau/JetTagComputer/interface/JetTagComputerRecord.h"
26	>	#include "RecoBTag/SecondaryVertex/interface/CombinedSVComputer.h"
27	>	#include "DataFormats/TrackReco/interface/Track.h"
28
29		//
30		// constants, enums and typedefs
#	Line 50 | Line 55 | NtupleWriter::NtupleWriter(const edm::Pa
55		doMuons = iConfig.getParameter<bool>("doMuons");
56		doTaus = iConfig.getParameter<bool>("doTaus");
57		doJets = iConfig.getParameter<bool>("doJets");
58	+	doGenJets = iConfig.getParameter<bool>("doGenJets");
59		doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");
60		doPhotons = iConfig.getParameter<bool>("doPhotons");
61		doMET = iConfig.getParameter<bool>("doMET");
#	Line 58 | Line 64 | NtupleWriter::NtupleWriter(const edm::Pa
64		doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
65		doPV = iConfig.getParameter<bool>("doPV");
66		doTopJets = iConfig.getParameter<bool>("doTopJets");
67	+	doTopJetsConstituents = iConfig.getParameter<bool>("doTopJetsConstituents");
68		doTrigger = iConfig.getParameter<bool>("doTrigger");
69	+	SVComputer_  = iConfig.getUntrackedParameter<edm::InputTag>("svComputer",edm::InputTag("combinedSecondaryVertex"));
70	+	doTagInfos = iConfig.getUntrackedParameter<bool>("doTagInfos",false);
71	+	storePFsAroundLeptons = iConfig.getUntrackedParameter<bool>("storePFsAroundLeptons",false);
72
73		// initialization of tree variables
74
#	Line 66 | Line 76 | NtupleWriter::NtupleWriter(const edm::Pa
76		tr->Branch("event",&event);
77		tr->Branch("luminosityBlock",&luminosityBlock);
78		tr->Branch("isRealData",&isRealData);
79	<	tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
79	>	tr->Branch("rho",&rho);
80	>	rho_source = iConfig.getParameter<edm::InputTag>("rho_source");
81	>
82	>	//tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
83		if(doLumiInfo){
84		tr->Branch("intgRecLumi",&intgRecLumi);
85		tr->Branch("intgDelLumi",&intgDelLumi);
#	Line 104 | Line 117 | NtupleWriter::NtupleWriter(const edm::Pa
117		tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
118		}
119		}
120	+	if(doGenJets){
121	+	genjet_sources = iConfig.getParameter<std::vector<std::string> >("genjet_sources");
122	+	genjet_ptmin = iConfig.getParameter<double> ("genjet_ptmin");
123	+	genjet_etamax = iConfig.getParameter<double> ("genjet_etamax");
124	+	for(size_t j=0; j< genjet_sources.size(); ++j){
125	+	tr->Branch( genjet_sources[j].c_str(), "std::vector<Particle>", &genjets[j]);
126	+	}
127	+	}
128		if(doTopJets){
129		topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
130		topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
#	Line 112 | Line 133 | NtupleWriter::NtupleWriter(const edm::Pa
133		tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
134		}
135		}
136	+	if(doTopJetsConstituents){
137	+	topjet_constituents_sources = iConfig.getParameter<std::vector<std::string> >("topjet_constituents_sources");
138	+	tr->Branch( "PFParticles", "std::vector<PFParticle>", &pfparticles);
139	+	}
140		if(doGenTopJets){
141		gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
142		gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
143		gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
144		for(size_t j=0; j< gentopjet_sources.size(); ++j){
145	<	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<TopJet>", &gentopjets[j]);
145	>	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<GenTopJet>", &gentopjets[j]);
146		}
147		}
148		if(doPhotons){
#	Line 139 | Line 164 | NtupleWriter::NtupleWriter(const edm::Pa
164		}
165		}
166		if(doGenInfo){
167	+	genparticle_source= iConfig.getParameter<edm::InputTag>("genparticle_source");
168		tr->Branch("genInfo","GenInfo",&genInfo);
169		tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
170		}
#	Line 150 | Line 176 | NtupleWriter::NtupleWriter(const edm::Pa
176		//tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
177		//tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
178		}
179	+	if (storePFsAroundLeptons){
180	+	pf_around_leptons_source = iConfig.getParameter<std::string>("pf_around_leptons_source");
181	+	}
182		newrun = true;
183		}
184
#	Line 181 | Line 210 | NtupleWriter::analyze(const edm::Event&
210		luminosityBlock = iEvent.luminosityBlock();
211		isRealData      = iEvent.isRealData();
212
213	<	if(isRealData){
214	<	edm::Handle<bool> bool_handle;
215	<	iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
216	<	HBHENoiseFilterResult = *(bool_handle.product());
217	<	}
218	<	else HBHENoiseFilterResult = false;
213	>	edm::Handle<double> m_rho;
214	>	iEvent.getByLabel(rho_source,m_rho);
215	>	rho=*m_rho;
216	>
217	>	//    if(isRealData){
218	>	//      edm::Handle<bool> bool_handle;
219	>	//      iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
220	>	//      HBHENoiseFilterResult = *(bool_handle.product());
221	>	//    }
222	>	//    else HBHENoiseFilterResult = false;
223
224		// ------------- primary vertices and beamspot  -------------
225
#	Line 202 | Line 235 | NtupleWriter::analyze(const edm::Event&
235		reco::Vertex reco_pv = reco_pvs[i];
236
237		PrimaryVertex pv;
238	<	pv.x =  reco_pv.x();
239	<	pv.y =  reco_pv.y();
240	<	pv.z =  reco_pv.z();
241	<	pv.nTracks =  reco_pv.nTracks();
242	<	//pv.isValid =  reco_pv.isValid();
243	<	pv.chi2 =  reco_pv.chi2();
244	<	pv.ndof =  reco_pv.ndof();
238	>	pv.set_x( reco_pv.x());
239	>	pv.set_y( reco_pv.y());
240	>	pv.set_z( reco_pv.z());
241	>	pv.set_nTracks( reco_pv.nTracks());
242	>	//pv.set_isValid( reco_pv.isValid());
243	>	pv.set_chi2( reco_pv.chi2());
244	>	pv.set_ndof( reco_pv.ndof());
245
246		pvs[j].push_back(pv);
247		}
#	Line 226 | Line 259 | NtupleWriter::analyze(const edm::Event&
259		// ------------- generator info -------------
260
261		if(doGenInfo){
262	<	genInfo.weights.clear();
263	<	genInfo.binningValues.clear();
262	>	genInfo.clear_weights();
263	>	genInfo.clear_binningValues();
264		genps.clear();
265
266		edm::Handle<GenEventInfoProduct> genEventInfoProduct;
267		iEvent.getByLabel("generator", genEventInfoProduct);
268		const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
269
270	<	genInfo.binningValues = genEventInfo.binningValues();
271	<	genInfo.weights = genEventInfo.weights();
272	<	genInfo.alphaQCD = genEventInfo.alphaQCD();
273	<	genInfo.alphaQED = genEventInfo.alphaQED();
274	<	genInfo.qScale = genEventInfo.qScale();
270	>	for(unsigned int k=0; k<genEventInfo.binningValues().size();++k){
271	>	genInfo.add_binningValue(genEventInfo.binningValues().at(k));
272	>	}
273	>	for(unsigned int k=0; k<genEventInfo.weights().size();++k){
274	>	genInfo.add_weight(genEventInfo.weights().at(k));
275	>	}
276	>	genInfo.set_alphaQCD(genEventInfo.alphaQCD());
277	>	genInfo.set_alphaQED(genEventInfo.alphaQED());
278	>	genInfo.set_qScale(genEventInfo.qScale());
279
280		const gen::PdfInfo* pdf = genEventInfo.pdf();
281		if(pdf){
282	<	genInfo.pdf_id1=pdf->id.first;
283	<	genInfo.pdf_id2=pdf->id.second;
284	<	genInfo.pdf_x1=pdf->x.first;
285	<	genInfo.pdf_x2=pdf->x.second;
286	<	genInfo.pdf_scalePDF=pdf->scalePDF;
287	<	genInfo.pdf_xPDF1=pdf->xPDF.first;
288	<	genInfo.pdf_xPDF2=pdf->xPDF.second;
282	>	genInfo.set_pdf_id1(pdf->id.first);
283	>	genInfo.set_pdf_id2(pdf->id.second);
284	>	genInfo.set_pdf_x1(pdf->x.first);
285	>	genInfo.set_pdf_x2(pdf->x.second);
286	>	genInfo.set_pdf_scalePDF(pdf->scalePDF);
287	>	genInfo.set_pdf_xPDF1(pdf->xPDF.first);
288	>	genInfo.set_pdf_xPDF2(pdf->xPDF.second);
289		}
290		else{
291	<	genInfo.pdf_id1=-999;
292	<	genInfo.pdf_id2=-999;
293	<	genInfo.pdf_x1=-999;
294	<	genInfo.pdf_x2=-999;
295	<	genInfo.pdf_scalePDF=-999;
296	<	genInfo.pdf_xPDF1=-999;
297	<	genInfo.pdf_xPDF2=-999;
291	>	genInfo.set_pdf_id1(-999);
292	>	genInfo.set_pdf_id2(-999);
293	>	genInfo.set_pdf_x1(-999);
294	>	genInfo.set_pdf_x2(-999);
295	>	genInfo.set_pdf_scalePDF(-999);
296	>	genInfo.set_pdf_xPDF1(-999);
297	>	genInfo.set_pdf_xPDF2(-999);
298		}
299
300		edm::Handle<std::vector<PileupSummaryInfo> > pus;
301		iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
302	<	genInfo.pileup_NumInteractions_intime=0;
303	<	genInfo.pileup_NumInteractions_ootbefore=0;
304	<	genInfo.pileup_NumInteractions_ootafter=0;
302	>	genInfo.set_pileup_NumInteractions_intime(0);
303	>	genInfo.set_pileup_NumInteractions_ootbefore(0);
304	>	genInfo.set_pileup_NumInteractions_ootafter(0);
305		if(pus.isValid()){
306	<	genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
306	>	genInfo.set_pileup_TrueNumInteractions ( (float) pus->at(0).getTrueNumInteractions());
307		for(size_t i=0; i<pus->size(); ++i){
308		if(pus->at(i).getBunchCrossing() == 0) // intime pileup
309	<	genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
309	>	genInfo.set_pileup_NumInteractions_intime( genInfo.pileup_NumInteractions_intime() + pus->at(i).getPU_NumInteractions());
310		else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
311	<	genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
311	>	genInfo.set_pileup_NumInteractions_ootbefore( genInfo.pileup_NumInteractions_ootbefore() + pus->at(i).getPU_NumInteractions());
312		}
313		else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
314	<	genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
314	>	genInfo.set_pileup_NumInteractions_ootafter( genInfo.pileup_NumInteractions_ootafter() + pus->at(i).getPU_NumInteractions());
315		}
316		}
317		}
318
319		edm::Handle<reco::GenParticleCollection> genPartColl;
320	<	iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
320	>	iEvent.getByLabel(genparticle_source, genPartColl);
321		int index=-1;
322		for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
323		index++;
324
325	<	//write out only top quarks and status 3 particles (works fine only for MadGraph)
326	<	if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
325	>	//write out only top quarks,final state leptons and status 3 particles (works fine only for MadGraph)
326	>	bool islepton = iter->status()==1 && abs(iter->pdgId())>=11 && abs(iter->pdgId())<=16 ;
327	>	if(abs(iter->pdgId())==6 || iter->status()==3 || islepton ||  doAllGenParticles){
328		GenParticle genp;
329	<	genp.charge = iter->charge();;
330	<	genp.pt = iter->p4().pt();
331	<	genp.eta = iter->p4().eta();
332	<	genp.phi = iter->p4().phi();
333	<	genp.energy = iter->p4().E();
334	<	genp.index =index;
335	<	genp.status = iter->status();
336	<	genp.pdgId = iter->pdgId();
337	<
338	<	genp.mother1=-1;
339	<	genp.mother2=-1;
340	<	genp.daughter1=-1;
341	<	genp.daughter2=-1;
329	>	genp.set_charge(iter->charge());
330	>	genp.set_pt(iter->p4().pt());
331	>	genp.set_eta(iter->p4().eta());
332	>	genp.set_phi(iter->p4().phi());
333	>	genp.set_energy(iter->p4().E());
334	>	genp.set_index(index);
335	>	genp.set_status( iter->status());
336	>	genp.set_pdgId( iter->pdgId());
337	>
338	>	genp.set_mother1(-1);
339	>	genp.set_mother2(-1);
340	>	genp.set_daughter1(-1);
341	>	genp.set_daughter2(-1);
342
343		int nm=iter->numberOfMothers();
344		int nd=iter->numberOfDaughters();
345
346
347	<	if (nm>0) genp.mother1 = iter->motherRef(0).key();
348	<	if (nm>1) genp.mother2 = iter->motherRef(1).key();
349	<	if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
350	<	if (nd>1) genp.daughter2 = iter->daughterRef(1).key();
351	<	//std::cout << genp.index <<"  pdgId = " << genp.pdgId << "  mo1 = " << genp.mother1 << "  mo2 = " << genp.mother2 <<"  da1 = " << genp.daughter1 << "  da2 = " << genp.daughter2 <<std::endl;
347	>	if (nm>0) genp.set_mother1( iter->motherRef(0).key());
348	>	if (nm>1) genp.set_mother2( iter->motherRef(1).key());
349	>	if (nd>0) genp.set_daughter1( iter->daughterRef(0).key());
350	>	if (nd>1) genp.set_daughter2( iter->daughterRef(1).key());
351	>
352		genps.push_back(genp);
353		}
354		}
355		}
356
357	+	// ------------- full PF collection -------------
358	+
359	+
360		// ------------- electrons -------------
361		if(doElectrons){
362
363	<	edm::Handle<reco::ConversionCollection> hConversions;
364	<	iEvent.getByLabel("allConversions", hConversions);
363	>	//      edm::Handle<reco::ConversionCollection> hConversions;
364	>	//      iEvent.getByLabel("allConversions", hConversions);
365
366	<	edm::Handle<reco::BeamSpot> beamSpot;
367	<	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
368	<	const reco::BeamSpot & bsp = *beamSpot;
366	>	//      edm::Handle<reco::BeamSpot> beamSpot;
367	>	//      iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
368	>	//      const reco::BeamSpot & bsp = *beamSpot;
369
370		for(size_t j=0; j< electron_sources.size(); ++j){
371		eles[j].clear();
#	Line 336 | Line 377 | NtupleWriter::analyze(const edm::Event&
377		pat::Electron pat_ele = pat_electrons[i];
378		Electron ele;
379
380	<	ele.charge =  pat_ele.charge();
381	<	ele.pt =  pat_ele.pt();
382	<	ele.eta =  pat_ele.eta();
383	<	ele.phi =  pat_ele.phi();
384	<	ele.energy =  pat_ele.energy();
385	<	ele.vertex_x = pat_ele.vertex().x();
386	<	ele.vertex_y = pat_ele.vertex().y();
387	<	ele.vertex_z = pat_ele.vertex().z();
388	<	ele.supercluster_eta = pat_ele.superCluster()->eta();
389	<	ele.supercluster_phi = pat_ele.superCluster()->phi();
390	<	ele.dB = pat_ele.dB();
391	<	//ele.particleIso = pat_ele.particleIso();
392	<	ele.neutralHadronIso = pat_ele.neutralHadronIso();
393	<	ele.chargedHadronIso = pat_ele.chargedHadronIso();
394	<	ele.trackIso = pat_ele.trackIso();
395	<	ele.photonIso = pat_ele.photonIso();
396	<	ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
397	<	ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
398	<	ele.gsfTrack_px= pat_ele.gsfTrack()->px();
399	<	ele.gsfTrack_py= pat_ele.gsfTrack()->py();
400	<	ele.gsfTrack_pz= pat_ele.gsfTrack()->pz();
401	<	ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
402	<	ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
403	<	ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
404	<	ele.passconversionveto = !ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position());
405	<	ele.dEtaIn=pat_ele.deltaEtaSuperClusterTrackAtVtx();
406	<	ele.dPhiIn=pat_ele.deltaPhiSuperClusterTrackAtVtx();
407	<	ele.sigmaIEtaIEta=pat_ele.sigmaIetaIeta();
408	<	ele.HoverE=pat_ele.hadronicOverEm();
409	<	ele.fbrem=pat_ele.fbrem();
410	<	ele.EoverPIn=pat_ele.eSuperClusterOverP();
411	<	ele.EcalEnergy=pat_ele.ecalEnergy();
412	<	//ele.mvaTrigV0=pat_ele.electronID("mvaTrigV0");
413	<	//ele.mvaNonTrigV0=pat_ele.electronID("mvaNonTrigV0");
380	>	ele.set_charge( pat_ele.charge());
381	>	ele.set_pt( pat_ele.pt());
382	>	ele.set_eta( pat_ele.eta());
383	>	ele.set_phi( pat_ele.phi());
384	>	ele.set_energy( pat_ele.energy());
385	>	ele.set_vertex_x(pat_ele.vertex().x());
386	>	ele.set_vertex_y(pat_ele.vertex().y());
387	>	ele.set_vertex_z(pat_ele.vertex().z());
388	>	ele.set_supercluster_eta(pat_ele.superCluster()->eta());
389	>	ele.set_supercluster_phi(pat_ele.superCluster()->phi());
390	>	ele.set_dB(pat_ele.dB());
391	>	//ele.set_particleIso(pat_ele.particleIso());
392	>	ele.set_neutralHadronIso(pat_ele.neutralHadronIso());
393	>	ele.set_chargedHadronIso(pat_ele.chargedHadronIso());
394	>	ele.set_trackIso(pat_ele.trackIso());
395	>	ele.set_photonIso(pat_ele.photonIso());
396	>	ele.set_puChargedHadronIso(pat_ele.puChargedHadronIso());
397	>	ele.set_gsfTrack_trackerExpectedHitsInner_numberOfLostHits(pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits());
398	>	ele.set_gsfTrack_px( pat_ele.gsfTrack()->px());
399	>	ele.set_gsfTrack_py( pat_ele.gsfTrack()->py());
400	>	ele.set_gsfTrack_pz( pat_ele.gsfTrack()->pz());
401	>	ele.set_gsfTrack_vx( pat_ele.gsfTrack()->vx());
402	>	ele.set_gsfTrack_vy( pat_ele.gsfTrack()->vy());
403	>	ele.set_gsfTrack_vz( pat_ele.gsfTrack()->vz());
404	>	//ele.set_passconversionveto(!ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position()));
405	>	ele.set_passconversionveto(pat_ele.passConversionVeto());
406	>	ele.set_dEtaIn(pat_ele.deltaEtaSuperClusterTrackAtVtx());
407	>	ele.set_dPhiIn(pat_ele.deltaPhiSuperClusterTrackAtVtx());
408	>	ele.set_sigmaIEtaIEta(pat_ele.sigmaIetaIeta());
409	>	ele.set_HoverE(pat_ele.hadronicOverEm());
410	>	ele.set_fbrem(pat_ele.fbrem());
411	>	ele.set_EoverPIn(pat_ele.eSuperClusterOverP());
412	>	ele.set_EcalEnergy(pat_ele.ecalEnergy());
413	>	ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
414	>	ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
415	>	float AEff03 = 0.00;
416	>	if(isRealData){
417	>	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2011);
418	>	}else{
419	>	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAFall11MC);
420	>	}
421	>	ele.set_AEff(AEff03);
422
423		eles[j].push_back(ele);
424	+
425	+	if (storePFsAroundLeptons){
426	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
427	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
428	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
429	+	StorePFCandsInCone(&ele, pf_cands, 0.4);
430	+	}
431	+
432		}
433		}
434		}
#	Line 389 | Line 446 | NtupleWriter::analyze(const edm::Event&
446		pat::Muon pat_mu = pat_muons[i];
447
448		Muon mu;
449	<	mu.charge =  pat_mu.charge();
450	<	mu.pt =  pat_mu.pt();
451	<	mu.eta =  pat_mu.eta();
452	<	mu.phi =  pat_mu.phi();
453	<	mu.energy =  pat_mu.energy();
454	<	mu.vertex_x = pat_mu.vertex().x();
455	<	mu.vertex_y = pat_mu.vertex().y();
456	<	mu.vertex_z = pat_mu.vertex().z();
457	<	mu.dB = pat_mu.dB();
458	<	//mu.particleIso = pat_mu.particleIso();
459	<	mu.neutralHadronIso = pat_mu.neutralHadronIso();
460	<	mu.chargedHadronIso = pat_mu.chargedHadronIso();
461	<	mu.trackIso = pat_mu.trackIso();
462	<	mu.photonIso = pat_mu.photonIso();
463	<	mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
464	<	mu.isGlobalMuon = pat_mu.isGlobalMuon();
465	<	mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
466	<	mu.isTrackerMuon = pat_mu.isTrackerMuon();
467	<	mu.numberOfMatchedStations = pat_mu.numberOfMatchedStations();
449	>	mu.set_charge( pat_mu.charge());
450	>	mu.set_pt( pat_mu.pt());
451	>	mu.set_eta( pat_mu.eta());
452	>	mu.set_phi( pat_mu.phi());
453	>	mu.set_energy( pat_mu.energy());
454	>	mu.set_vertex_x ( pat_mu.vertex().x());
455	>	mu.set_vertex_y ( pat_mu.vertex().y());
456	>	mu.set_vertex_z ( pat_mu.vertex().z());
457	>	mu.set_dB ( pat_mu.dB());
458	>	//mu.particleIso ( pat_mu.particleIso());
459	>	mu.set_neutralHadronIso ( pat_mu.neutralHadronIso());
460	>	mu.set_chargedHadronIso ( pat_mu.chargedHadronIso());
461	>	mu.set_trackIso ( pat_mu.trackIso());
462	>	mu.set_photonIso ( pat_mu.photonIso());
463	>	mu.set_puChargedHadronIso ( pat_mu.puChargedHadronIso());
464	>	mu.set_isGlobalMuon ( pat_mu.isGlobalMuon());
465	>	mu.set_isPFMuon ( pat_mu.isPFMuon());
466	>	mu.set_isStandAloneMuon ( pat_mu.isStandAloneMuon());
467	>	mu.set_isTrackerMuon ( pat_mu.isTrackerMuon());
468	>	mu.set_numberOfMatchedStations ( pat_mu.numberOfMatchedStations());
469		reco::TrackRef globalTrack = pat_mu.globalTrack();
470		if(!globalTrack.isNull()){
471	<	mu.globalTrack_chi2 = globalTrack->chi2();
472	<	mu.globalTrack_ndof = globalTrack->ndof();
473	<	mu.globalTrack_d0 = globalTrack->d0();
474	<	mu.globalTrack_d0Error = globalTrack->d0Error();
475	<	mu.globalTrack_numberOfValidHits = globalTrack->numberOfValidHits();
476	<	mu.globalTrack_numberOfLostHits = globalTrack->numberOfLostHits();
471	>	mu.set_globalTrack_chi2 ( globalTrack->chi2());
472	>	mu.set_globalTrack_ndof ( globalTrack->ndof());
473	>	mu.set_globalTrack_d0 ( globalTrack->d0());
474	>	mu.set_globalTrack_d0Error ( globalTrack->d0Error());
475	>	mu.set_globalTrack_numberOfValidHits ( globalTrack->numberOfValidHits());
476	>	mu.set_globalTrack_numberOfLostHits ( globalTrack->numberOfLostHits());
477	>	mu.set_globalTrack_numberOfValidMuonHits(globalTrack->hitPattern().numberOfValidMuonHits() );
478		}
479		else{
480	<	mu.globalTrack_chi2 = 0;
481	<	mu.globalTrack_ndof = 0;
482	<	mu.globalTrack_d0 = 0;
483	<	mu.globalTrack_d0Error = 0;
484	<	mu.globalTrack_numberOfValidHits = 0;
485	<	mu.globalTrack_numberOfLostHits = 0;
480	>	mu.set_globalTrack_chi2 ( 0);
481	>	mu.set_globalTrack_ndof ( 0);
482	>	mu.set_globalTrack_d0 ( 0);
483	>	mu.set_globalTrack_d0Error ( 0);
484	>	mu.set_globalTrack_numberOfValidHits ( 0);
485	>	mu.set_globalTrack_numberOfLostHits ( 0);
486		}
487		reco::TrackRef innerTrack = pat_mu.innerTrack();
488		if(!innerTrack.isNull()){
489	<	mu.innerTrack_chi2 = innerTrack->chi2();
490	<	mu.innerTrack_ndof = innerTrack->ndof();
491	<	mu.innerTrack_d0 = innerTrack->d0();
492	<	mu.innerTrack_d0Error = innerTrack->d0Error();
493	<	mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
494	<	mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
495	<	mu.innerTrack_trackerLayersWithMeasurement = innerTrack->hitPattern().trackerLayersWithMeasurement();
496	<	mu.innerTrack_numberOfValidPixelHits = innerTrack->hitPattern().numberOfValidPixelHits();
489	>	mu.set_innerTrack_chi2 ( innerTrack->chi2());
490	>	mu.set_innerTrack_ndof ( innerTrack->ndof());
491	>	mu.set_innerTrack_d0 ( innerTrack->d0());
492	>	mu.set_innerTrack_d0Error ( innerTrack->d0Error());
493	>	mu.set_innerTrack_numberOfValidHits ( innerTrack->numberOfValidHits());
494	>	mu.set_innerTrack_numberOfLostHits ( innerTrack->numberOfLostHits());
495	>	mu.set_innerTrack_trackerLayersWithMeasurement ( innerTrack->hitPattern().trackerLayersWithMeasurement());
496	>	mu.set_innerTrack_numberOfValidPixelHits ( innerTrack->hitPattern().numberOfValidPixelHits());
497		}
498		else{
499	<	mu.innerTrack_chi2 = 0;
500	<	mu.innerTrack_ndof = 0;
501	<	mu.innerTrack_d0 = 0;
502	<	mu.innerTrack_d0Error = 0;
503	<	mu.innerTrack_numberOfValidHits = 0;
504	<	mu.innerTrack_numberOfLostHits = 0;
505	<	mu.innerTrack_trackerLayersWithMeasurement = 0;
506	<	mu.innerTrack_numberOfValidPixelHits = 0;
499	>	mu.set_innerTrack_chi2 ( 0);
500	>	mu.set_innerTrack_ndof ( 0);
501	>	mu.set_innerTrack_d0 ( 0);
502	>	mu.set_innerTrack_d0Error ( 0);
503	>	mu.set_innerTrack_numberOfValidHits ( 0);
504	>	mu.set_innerTrack_numberOfLostHits ( 0);
505	>	mu.set_innerTrack_trackerLayersWithMeasurement ( 0);
506	>	mu.set_innerTrack_numberOfValidPixelHits ( 0);
507		}
508		reco::TrackRef outerTrack = pat_mu.outerTrack();
509		if(!outerTrack.isNull()){
510	<	mu.outerTrack_chi2 = outerTrack->chi2();
511	<	mu.outerTrack_ndof = outerTrack->ndof();
512	<	mu.outerTrack_d0 = outerTrack->d0();
513	<	mu.outerTrack_d0Error = outerTrack->d0Error();
514	<	mu.outerTrack_numberOfValidHits = outerTrack->numberOfValidHits();
515	<	mu.outerTrack_numberOfLostHits = outerTrack->numberOfLostHits();
510	>	mu.set_outerTrack_chi2 ( outerTrack->chi2());
511	>	mu.set_outerTrack_ndof ( outerTrack->ndof());
512	>	mu.set_outerTrack_d0 ( outerTrack->d0());
513	>	mu.set_outerTrack_d0Error ( outerTrack->d0Error());
514	>	mu.set_outerTrack_numberOfValidHits ( outerTrack->numberOfValidHits());
515	>	mu.set_outerTrack_numberOfLostHits ( outerTrack->numberOfLostHits());
516		}
517		else{
518	<	mu.outerTrack_chi2 = 0;
519	<	mu.outerTrack_ndof = 0;
520	<	mu.outerTrack_d0 = 0;
521	<	mu.outerTrack_d0Error = 0;
522	<	mu.outerTrack_numberOfValidHits = 0;
523	<	mu.outerTrack_numberOfLostHits = 0;
518	>	mu.set_outerTrack_chi2 ( 0);
519	>	mu.set_outerTrack_ndof ( 0);
520	>	mu.set_outerTrack_d0 ( 0);
521	>	mu.set_outerTrack_d0Error ( 0);
522	>	mu.set_outerTrack_numberOfValidHits ( 0);
523	>	mu.set_outerTrack_numberOfLostHits ( 0);
524		}
525
526		mus[j].push_back(mu);
527	+
528	+	if (storePFsAroundLeptons){
529	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
530	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
531	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
532	+	StorePFCandsInCone(&mu, pf_cands, 0.4);
533	+	}
534	+
535		}
536		}
537		}
#	Line 485 | Line 552 | NtupleWriter::analyze(const edm::Event&
552		if(fabs(pat_tau.eta()) > tau_etamax) continue;
553
554		Tau tau;
555	<	tau.charge =  pat_tau.charge();
556	<	tau.pt =  pat_tau.pt();
557	<	tau.eta =  pat_tau.eta();
558	<	tau.phi =  pat_tau.phi();
559	<	tau.energy =  pat_tau.energy();
560	<	tau.decayModeFinding = pat_tau.tauID("decayModeFinding")>0.5;
561	<	tau.byVLooseCombinedIsolationDeltaBetaCorr  = pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5;
562	<	tau.byLooseCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5;
563	<	tau.byMediumCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5;
564	<	tau.byTightCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5;
565	<	tau.againstElectronLoose  = pat_tau.tauID("againstElectronLoose")>0.5;
566	<	tau.againstElectronMedium = pat_tau.tauID("againstElectronMedium")>0.5;
567	<	tau.againstElectronTight = pat_tau.tauID("againstElectronTight")>0.5;
568	<	tau.againstElectronMVA  = pat_tau.tauID("againstElectronMVA")>0.5;
569	<	tau.againstMuonLoose = pat_tau.tauID("againstMuonLoose")>0.5;
570	<	tau.againstMuonMedium = pat_tau.tauID("againstMuonMedium")>0.5;
571	<	tau.againstMuonTight = pat_tau.tauID("againstMuonTight")>0.5;
572	<
573	<	reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
574	<	if(!leadPFCand.isNull()){
575	<	tau.leadPFCand_px = leadPFCand->px();
576	<	tau.leadPFCand_py = leadPFCand->py();
577	<	tau.leadPFCand_pz = leadPFCand->pz();
555	>	tau.set_charge( pat_tau.charge());
556	>	tau.set_pt( pat_tau.pt());
557	>	tau.set_eta( pat_tau.eta());
558	>	tau.set_phi( pat_tau.phi());
559	>	tau.set_energy( pat_tau.energy());
560	>	tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
561	>	//tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
562	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
563	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
564	>	tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
565	>	tau.set_byLooseIsolationMVA( pat_tau.tauID("byLooseIsolationMVA")>0.5);
566	>	tau.set_byMediumIsolationMVA( pat_tau.tauID("byMediumIsolationMVA")>0.5);
567	>	tau.set_byTightIsolationMVA( pat_tau.tauID("byTightIsolationMVA")>0.5);
568	>	tau.set_byLooseIsolationMVA2( pat_tau.tauID("byLooseIsolationMVA2")>0.5);
569	>	tau.set_byMediumIsolationMVA2( pat_tau.tauID("byMediumIsolationMVA2")>0.5);
570	>	tau.set_byTightIsolationMVA2( pat_tau.tauID("byTightIsolationMVA2")>0.5);
571	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr3Hits(  pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits")>0.5);
572	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits")>0.5);
573	>	tau.set_byTightCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr3Hits")>0.5);
574	>	tau.set_againstElectronLooseMVA3  ( pat_tau.tauID("againstElectronLooseMVA3")>0.5);
575	>	tau.set_againstElectronMediumMVA3 ( pat_tau.tauID("againstElectronMediumMVA3")>0.5);
576	>	tau.set_againstElectronTightMVA3 ( pat_tau.tauID("againstElectronTightMVA3")>0.5);
577	>	tau.set_againstElectronVTightMVA3 ( pat_tau.tauID("againstElectronVTightMVA3")>0.5);
578	>	tau.set_againstMuonLoose2 ( pat_tau.tauID("againstMuonLoose2")>0.5);
579	>	tau.set_againstMuonMedium2 ( pat_tau.tauID("againstMuonMedium2")>0.5);
580	>	tau.set_againstMuonTight2 ( pat_tau.tauID("againstMuonTight2")>0.5);
581	>	tau.set_byIsolationMVAraw(  pat_tau.tauID("byIsolationMVAraw"));
582	>	tau.set_byIsolationMVA2raw(  pat_tau.tauID("byIsolationMVA2raw"));
583	>	tau.set_decayMode( pat_tau.decayMode() );
584	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw"));
585	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw3Hits( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));
586	>
587	>	//       std::cout << pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") << std::endl;
588	>
589	>	//       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
590	>	//       if(!leadPFCand.isNull()){
591	>	//         tau.set_leadPFCand_px ( leadPFCand->px());
592	>	//         tau.set_leadPFCand_py ( leadPFCand->py());
593	>	//         tau.set_leadPFCand_pz ( leadPFCand->pz());
594	>	//       }
595	>	//       else{
596	>	//         tau.set_leadPFCand_px ( 0);
597	>	//         tau.set_leadPFCand_py ( 0);
598	>	//         tau.set_leadPFCand_pz ( 0);
599	>	//       }
600	>	taus[j].push_back(tau);
601	>
602	>	// store isolation info only for identified taus
603	>	if (pat_tau.tauID("decayModeFinding")>0.5){
604	>	bool storepfs = (pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5) || (pat_tau.tauID("byLooseIsolationMVA")>0.5);
605	>	if (storePFsAroundLeptons && storepfs){
606	>	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
607	>	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
608	>	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
609	>	StorePFCandsInCone(&tau, pf_cands, 0.4);
610	>	}
611		}
612	<	else{
613	<	tau.leadPFCand_px = 0;
614	<	tau.leadPFCand_py = 0;
615	<	tau.leadPFCand_pz = 0;
612	>	}
613	>	}
614	>	}
615	>
616	>	//-------------- gen jets -------------
617	>
618	>	if(doGenJets){
619	>	for(size_t j=0; j< genjet_sources.size(); ++j){
620	>
621	>	genjets[j].clear();
622	>
623	>	edm::Handle< std::vector<reco::GenJet> > genjet_handle;
624	>	iEvent.getByLabel(genjet_sources[j], genjet_handle);
625	>	const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
626	>
627	>	for (unsigned int i = 0; i < gen_jets.size(); ++i) {
628	>
629	>	const reco::GenJet* gen_jet = &gen_jets[i];
630	>	if(gen_jet->pt() < genjet_ptmin) continue;
631	>	if(fabs(gen_jet->eta()) > genjet_etamax) continue;
632	>
633	>	Particle jet;
634	>	jet.set_charge(gen_jet->charge());
635	>	jet.set_pt(gen_jet->pt());
636	>	jet.set_eta(gen_jet->eta());
637	>	jet.set_phi(gen_jet->phi());
638	>	jet.set_energy(gen_jet->energy());
639	>
640	>	// recalculate the jet charge
641	>	int jet_charge = 0;
642	>	std::vector<const reco::GenParticle * > jetgenps = gen_jet->getGenConstituents();
643	>	for(unsigned int l = 0; l<jetgenps.size(); ++l){
644	>	jet_charge +=  jetgenps[l]->charge();
645		}
646	<	taus[j].push_back(tau);
646	>	jet.set_charge(jet_charge);
647	>
648	>	genjets[j].push_back(jet);
649		}
650		}
651		}
#	Line 540 | Line 671 | NtupleWriter::analyze(const edm::Event&
671		//       }
672
673		Jet jet;
674	<	jet.charge = pat_jet.charge();
675	<	jet.pt = pat_jet.pt();
676	<	jet.eta = pat_jet.eta();
677	<	jet.phi = pat_jet.phi();
678	<	jet.energy = pat_jet.energy();
679	<	jet.numberOfDaughters =pat_jet.numberOfDaughters();
674	>	jet.set_charge(pat_jet.charge());
675	>	jet.set_pt(pat_jet.pt());
676	>	jet.set_eta(pat_jet.eta());
677	>	jet.set_phi(pat_jet.phi());
678	>	jet.set_energy(pat_jet.energy());
679	>	jet.set_flavor(pat_jet.partonFlavour());
680	>	jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
681		const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
682	<	jet.nTracks = jettracks.size();
683	<	jet.jetArea = pat_jet.jetArea();
552	<	jet.pileup = pat_jet.pileup();
682	>	jet.set_nTracks ( jettracks.size());
683	>	jet.set_jetArea(pat_jet.jetArea());
684		if(pat_jet.isPFJet()){
685	<	jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
686	<	jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
687	<	jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
688	<	jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
689	<	jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
690	<	jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
691	<	jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
692	<	jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
693	<	jet.muonMultiplicity =pat_jet.muonMultiplicity();
694	<	jet.electronMultiplicity =pat_jet.electronMultiplicity();
695	<	jet.photonMultiplicity =pat_jet.photonMultiplicity();
696	<	}
697	<
698	<	jecUnc->setJetEta(pat_jet.eta());
699	<	jecUnc->setJetPt(pat_jet.pt());
700	<	jet.JEC_uncertainty = jecUnc->getUncertainty(true);
701	<	jet.JEC_factor_raw = pat_jet.jecFactor("Uncorrected");
702	<
703	<	jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
704	<	jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
705	<	jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
575	<	jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
576	<	jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
577	<	jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
685	>	jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
686	>	jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
687	>	jet.set_chargedEmEnergyFraction (pat_jet.chargedEmEnergyFraction());
688	>	jet.set_chargedHadronEnergyFraction (pat_jet.chargedHadronEnergyFraction());
689	>	jet.set_muonEnergyFraction (pat_jet.muonEnergyFraction());
690	>	jet.set_photonEnergyFraction (pat_jet.photonEnergyFraction());
691	>	jet.set_chargedMultiplicity (pat_jet.chargedMultiplicity());
692	>	jet.set_neutralMultiplicity (pat_jet.neutralMultiplicity());
693	>	jet.set_muonMultiplicity (pat_jet.muonMultiplicity());
694	>	jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
695	>	jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
696	>	}
697	>
698	>	jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
699	>
700	>	jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
701	>	jet.set_btag_simpleSecondaryVertexHighPur(pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
702	>	jet.set_btag_combinedSecondaryVertex(pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags"));
703	>	jet.set_btag_combinedSecondaryVertexMVA(pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
704	>	jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
705	>	jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
706
707	+
708		const reco::GenJet *genj = pat_jet.genJet();
709		if(genj){
710	<	jet.genjet_pt = genj->pt();
711	<	jet.genjet_eta = genj->eta();
712	<	jet.genjet_phi = genj->phi();
713	<	jet.genjet_energy = genj->energy();
585	<	if(doAllGenParticles){
586	<	std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
587	<	for(unsigned int l = 0; l<jetgenps.size(); ++l){
588	<	for(unsigned int k=0; k< genps.size(); ++k){
589	<	if(jetgenps[l]->pt() == genps[k].pt && jetgenps[l]->pdgId() == genps[k].pdgId){
590	<	jet.genparticles_indices.push_back(genps[k].index);
591	<	}
592	<	}
710	>
711	>	for(unsigned int k=0; k<genjets->size(); ++k){
712	>	if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta()){
713	>	jet.set_genjet_index(k);
714		}
594	–	if(jet.genparticles_indices.size()!= jetgenps.size())
595	–	std::cout << "WARNING: Found only " << jet.genparticles_indices.size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
715		}
716	+	//         if( jet.genjet_index()<0){
717	+	//           std::cout<< "genjet not found for " << genj->pt() << "  " << genj->eta() << std::endl;
718	+	//         }
719	+
720		}
721	+
722		jets[j].push_back(jet);
723		}
724		}
725		}
726
727	+
728		// ------------- top jets -------------
729		if(doTopJets){
730	+
731		for(size_t j=0; j< topjet_sources.size(); ++j){
732
733		topjets[j].clear();
#	Line 617 | Line 743 | NtupleWriter::analyze(const edm::Event&
743		if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
744
745		TopJet topjet;
746	<	topjet.charge = pat_topjet.charge();
747	<	topjet.pt = pat_topjet.pt();
748	<	topjet.eta = pat_topjet.eta();
749	<	topjet.phi = pat_topjet.phi();
750	<	topjet.energy = pat_topjet.energy();
751	<	topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
746	>	topjet.set_charge(pat_topjet.charge());
747	>	topjet.set_pt(pat_topjet.pt());
748	>	topjet.set_eta(pat_topjet.eta());
749	>	topjet.set_phi(pat_topjet.phi());
750	>	topjet.set_energy(pat_topjet.energy());
751	>	topjet.set_flavor(pat_topjet.partonFlavour());
752	>	topjet.set_numberOfDaughters(pat_topjet.numberOfDaughters());
753		const reco::TrackRefVector&  topjettracks = pat_topjet.associatedTracks();
754	<	topjet.nTracks = topjettracks.size();
755	<	topjet.jetArea = pat_topjet.jetArea();
756	<	topjet.pileup = pat_topjet.pileup();
757	<	//       topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
758	<	//       topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
759	<	//       topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
760	<	//       topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
761	<	//       topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
762	<	//       topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
763	<	//       topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
764	<	//       topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
638	<	//       topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
639	<	//       topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
640	<	//       topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
641	<
642	<	jecUnc->setJetEta(pat_topjet.eta());
643	<	jecUnc->setJetPt(pat_topjet.pt());
644	<	topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
645	<	topjet.JEC_factor_raw = pat_topjet.jecFactor("Uncorrected");
646	<
647	<	topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
648	<	topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
649	<	topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
650	<	topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
651	<	topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
652	<	topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
754	>	topjet.set_nTracks( topjettracks.size());
755	>	topjet.set_jetArea( pat_topjet.jetArea());
756	>
757	>	topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
758	>
759	>	topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
760	>	topjet.set_btag_simpleSecondaryVertexHighPur(pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
761	>	topjet.set_btag_combinedSecondaryVertex(pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags"));
762	>	topjet.set_btag_combinedSecondaryVertexMVA(pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
763	>	topjet.set_btag_jetBProbability(pat_topjet.bDiscriminator("jetBProbabilityBJetTags"));
764	>	topjet.set_btag_jetProbability(pat_topjet.bDiscriminator("jetProbabilityBJetTags"));
765
766	+	/*
767		const reco::GenJet *genj = pat_topjet.genJet();
768		if(genj){
769	<	topjet.genjet_pt = genj->pt();
770	<	topjet.genjet_eta = genj->eta();
771	<	topjet.genjet_phi = genj->phi();
772	<	topjet.genjet_energy = genj->energy();
769	>	topjet.set_genjet_pt ( genj->pt());
770	>	topjet.set_genjet_eta ( genj->eta());
771	>	topjet.set_genjet_phi ( genj->phi());
772	>	topjet.set_genjet_energy ( genj->energy());
773		if(doAllGenParticles){
774		std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
775		for(unsigned int l = 0; l<jetgenps.size(); ++l){
776		for(unsigned int k=0; k< genps.size(); ++k){
777	<	if(jetgenps[l]->pt() == genps[k].pt && jetgenps[l]->pdgId() == genps[k].pdgId){
778	<	topjet.genparticles_indices.push_back(genps[k].index);
777	>	if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
778	>	topjet.add_genparticles_index(genps[k].index());
779	>	break;
780		}
781		}
782		}
783	<	if(topjet.genparticles_indices.size()!= jetgenps.size())
784	<	std::cout << "WARNING: Found only " << topjet.genparticles_indices.size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
783	>	if(topjet.genparticles_indices().size()!= jetgenps.size())
784	>	std::cout << "WARNING: Found only " << topjet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
785		}
786		}
787	+	*/
788	+
789	+	// add constituents to the jet, if requested
790	+	if (doTopJetsConstituents){
791	+
792	+	if (topjet_constituents_sources.size()>j){ //only add constituents if they are defined
793	+
794	+	edm::Handle<pat::JetCollection> pat_topjets_with_cands;
795	+	iEvent.getByLabel(topjet_constituents_sources[j], pat_topjets_with_cands);
796	+	pat::Jet* pat_topjet_wc = NULL;
797	+
798	+	for (unsigned int it = 0; it < pat_topjets_with_cands->size(); it++) {
799	+	const pat::Jet* cand =  dynamic_cast<pat::Jet const *>(&pat_topjets_with_cands->at(it));
800	+	double dphi = cand->phi() - pat_topjet.phi();
801	+	if (dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
802	+	if (dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
803	+	if (fabs(dphi)<0.5 && fabs(cand->eta()-pat_topjet.eta())<0.5){ // be generous: filtering, pruning... can change jet axis
804	+	pat_topjet_wc = const_cast<pat::Jet*>(cand);
805	+	break;
806	+	}
807	+	}
808	+
809	+	if (pat_topjet_wc){
810	+	StoreJetConstituents(pat_topjet_wc, &topjet);
811	+
812	+	// now run substructure information
813	+	JetProps substructure(&topjet);
814	+	substructure.set_pf_cands(&pfparticles);
815	+	double tau1 = substructure.GetNsubjettiness(1, Njettiness::onepass_kt_axes, 1., 2.0);
816	+	double tau2 = substructure.GetNsubjettiness(2, Njettiness::onepass_kt_axes, 1., 2.0);
817	+	double tau3 = substructure.GetNsubjettiness(3, Njettiness::onepass_kt_axes, 1., 2.0);
818	+	double qjets = substructure.GetQjetVolatility(iEvent.id().event(), 2.0);
819	+	topjet.set_tau1(tau1);
820	+	topjet.set_tau2(tau2);
821	+	topjet.set_tau3(tau3);
822	+	topjet.set_qjets_volatility(qjets);
823	+
824	+	}
825	+	}
826	+	}
827	+
828	+
829
830		for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
831		Particle subjet_v4;
832	<	subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
833	<	subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
834	<	subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
835	<	subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
836	<	topjet.subjets.push_back(subjet_v4);
832	>
833	>	reco::Candidate const * subjetd =  pat_topjet.daughter(k);
834	>	pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
835	>	if(patsubjetd)
836	>	{
837	>	subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
838	>	subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
839	>	subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
840	>	subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
841	>	topjet.add_subjet(subjet_v4);
842	>	//btag info
843	>	topjet.add_subFlavour(patsubjetd->partonFlavour());
844	>	topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
845	>	if (doTagInfos)
846	>	{
847	>	//ip tag info
848	>	reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
849	>	topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
850	>	topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
851	>	topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
852	>	topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
853	>	topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
854	>	topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
855	>	topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
856	>	topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
857	>	topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
858	>	topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
859	>	topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
860	>	topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));
861	>	topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
862	>	topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
863	>	topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
864	>	topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
865	>	topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
866	>	topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
867	>	topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
868	>	topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
869	>	topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
870	>	//sv tag info
871	>	reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
872	>	topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
873	>	topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
874	>	topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
875	>	topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
876	>	topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
877	>	topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
878	>	topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
879	>	std::vector<TLorentzVector> vp4; vp4.clear();
880	>	std::vector<float> vchi2; vchi2.clear();
881	>	std::vector<float> vndof; vndof.clear();
882	>	std::vector<float> vchi2ndof; vchi2ndof.clear();
883	>	std::vector<float> vtrsize; vtrsize.clear();
884	>	for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
885	>	{
886	>	ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
887	>	vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
888	>	vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());
889	>	vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());
890	>	vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());
891	>	vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());
892	>	}
893	>	topjet.add_subSecondaryVertex(vp4);
894	>	topjet.add_subVertexChi2(vchi2);
895	>	topjet.add_subVertexNdof(vndof);
896	>	topjet.add_subVertexNormalizedChi2(vchi2ndof);
897	>	topjet.add_subVertexTracksSize(vtrsize);
898	>	//try computer
899	>	edm::ESHandle<JetTagComputer> computerHandle;
900	>	iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
901	>	const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
902	>	if(computer)
903	>	{
904	>	computer->passEventSetup(iSetup);
905	>	std::vector<const reco::BaseTagInfo*>  baseTagInfos;
906	>	baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
907	>	baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );
908	>	JetTagComputer::TagInfoHelper helper(baseTagInfos);
909	>	reco::TaggingVariableList vars = computer->taggingVariables(helper);
910	>	topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
911	>	topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
912	>	topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
913	>	topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
914	>	}
915	>	}
916	>	}
917	>	else
918	>	{
919	>	//filling only standard information in case the subjet has not been pat-tified during the pattuples production
920	>	subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
921	>	subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
922	>	subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
923	>	subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
924	>	topjet.add_subjet(subjet_v4);
925	>	}
926	>
927	>
928		}
929		topjets[j].push_back(topjet);
930		}
931		}
932		}
933
934	<
934	>
935		// ------------- generator top jets -------------
936		if(doGenTopJets){
937		for(size_t j=0; j< gentopjet_sources.size(); ++j){
938
939		gentopjets[j].clear();
940	<
940	>
941		edm::Handle<reco::BasicJetCollection> reco_gentopjets;
695	–	//edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
942		iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
943
944		for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
#	Line 701 | Line 947 | NtupleWriter::analyze(const edm::Event&
947		if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
948		if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
949
950	<	TopJet gentopjet;
951	<	gentopjet.charge = reco_gentopjet.charge();
952	<	gentopjet.pt = reco_gentopjet.pt();
953	<	gentopjet.eta = reco_gentopjet.eta();
954	<	gentopjet.phi = reco_gentopjet.phi();
955	<	gentopjet.energy = reco_gentopjet.energy();
710	<	gentopjet.numberOfDaughters =reco_gentopjet.numberOfDaughters();
950	>	GenTopJet gentopjet;
951	>	gentopjet.set_charge(reco_gentopjet.charge());
952	>	gentopjet.set_pt(reco_gentopjet.pt());
953	>	gentopjet.set_eta(reco_gentopjet.eta());
954	>	gentopjet.set_phi(reco_gentopjet.phi());
955	>	gentopjet.set_energy(reco_gentopjet.energy());
956
957		for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
958		Particle subjet_v4;
959	<	subjet_v4.pt = reco_gentopjet.daughter(k)->p4().pt();
960	<	subjet_v4.eta = reco_gentopjet.daughter(k)->p4().eta();
961	<	subjet_v4.phi = reco_gentopjet.daughter(k)->p4().phi();
962	<	subjet_v4.energy = reco_gentopjet.daughter(k)->p4().E();
963	<	gentopjet.subjets.push_back(subjet_v4);
959	>	subjet_v4.set_pt(reco_gentopjet.daughter(k)->p4().pt());
960	>	subjet_v4.set_eta(reco_gentopjet.daughter(k)->p4().eta());
961	>	subjet_v4.set_phi(reco_gentopjet.daughter(k)->p4().phi());
962	>	subjet_v4.set_energy(reco_gentopjet.daughter(k)->p4().E());
963	>	gentopjet.add_subjet(subjet_v4);
964		}
965		gentopjets[j].push_back(gentopjet);
966		}
967		}
968		}
969
970	+
971		// ------------- photons -------------
972		if(doPhotons){
973		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 734 | Line 980 | NtupleWriter::analyze(const edm::Event&
980		for (unsigned int i = 0; i < pat_photons.size(); ++i) {
981		pat::Photon pat_photon = pat_photons[i];
982		Photon ph;
983	<	ph.charge = 0;
984	<	ph.pt =  pat_photon.pt();
985	<	ph.eta =  pat_photon.eta();
986	<	ph.phi =  pat_photon.phi();
987	<	ph.energy =  pat_photon.energy();
988	<	ph.vertex_x = pat_photon.vertex().x();
989	<	ph.vertex_y = pat_photon.vertex().y();
990	<	ph.vertex_z = pat_photon.vertex().z();
991	<	ph.supercluster_eta = pat_photon.superCluster()->eta();
992	<	ph.supercluster_phi = pat_photon.superCluster()->phi();
993	<	//       ph.neutralHadronIso = pat_photon.neutralHadronIso();
994	<	//       ph.chargedHadronIso = pat_photon.chargedHadronIso();
995	<	ph.trackIso = pat_photon.trackIso();
983	>	ph.set_charge(0);
984	>	ph.set_pt( pat_photon.pt());
985	>	ph.set_eta( pat_photon.eta());
986	>	ph.set_phi( pat_photon.phi());
987	>	ph.set_energy( pat_photon.energy());
988	>	ph.set_vertex_x(pat_photon.vertex().x());
989	>	ph.set_vertex_y(pat_photon.vertex().y());
990	>	ph.set_vertex_z(pat_photon.vertex().z());
991	>	ph.set_supercluster_eta(pat_photon.superCluster()->eta());
992	>	ph.set_supercluster_phi(pat_photon.superCluster()->phi());
993	>	//       ph.set_neutralHadronIso(pat_photon.neutralHadronIso());
994	>	//       ph.set_chargedHadronIso(pat_photon.chargedHadronIso());
995	>	ph.set_trackIso(pat_photon.trackIso());
996		phs[j].push_back(ph);
997		}
998		}
#	Line 766 | Line 1012 | NtupleWriter::analyze(const edm::Event&
1012		else{
1013		pat::MET pat_met = pat_mets[0];
1014
1015	<	met[j].pt=pat_met.pt();
1016	<	met[j].phi=pat_met.phi();
1017	<	met[j].mEtSig= pat_met.mEtSig();
1015	>	met[j].set_pt(pat_met.pt());
1016	>	met[j].set_phi(pat_met.phi());
1017	>	met[j].set_mEtSig(pat_met.mEtSig());
1018		}
1019
1020		}
#	Line 826 | Line 1072 | NtupleWriter::analyze(const edm::Event&
1072		tr->Fill();
1073		if(doLumiInfo)
1074		previouslumiblockwasfilled=true;
1075	+
1076	+	// clean up
1077	+	if(doTopJetsConstituents){
1078	+	pfparticles.clear();
1079	+	}
1080	+
1081	+	// need to do a check if necessary
1082	+	pfparticles.clear();
1083	+
1084		}
1085
1086
#	Line 854 | Line 1109 | void
1109		NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
1110		{
1111		if(doTrigger){
1112	<	bool setup_changed = false;
1113	<	hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1112	>	//bool setup_changed = false;
1113	>	//hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1114		newrun=true;
1115		}
1116
862	–	if(doJets || doTopJets){
863	–	edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
864	–	iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
865	–	JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
866	–	jecUnc = new JetCorrectionUncertainty(JetCorPar);
867	–	}
1117		}
1118
1119		// ------------ method called when ending the processing of a run  ------------
#	Line 917 | Line 1166 | NtupleWriter::fillDescriptions(edm::Conf
1166		descriptions.addDefault(desc);
1167		}
1168
1169	+	// ------------ method fills constituents of the pat_jet into the Ntuple and stores a reference
1170	+	// ------------ to those in the topjet
1171	+	// ------------ it is checked if the constituents have been stored already
1172	+	void
1173	+	NtupleWriter::StoreJetConstituents(pat::Jet* pat_jet, Jet* jet)
1174	+	{
1175	+	// checks if the pf cand has already been stored, only stores so far missing ones
1176	+	// PF candidates, then stores a reference to the pf candidate in the jet
1177	+	// also calculates the jet charge and sets it
1178	+
1179	+
1180	+	const std::vector<reco::PFCandidatePtr> jconstits = pat_jet->getPFConstituents();
1181	+
1182	+	unsigned int NstoredPFs = pfparticles.size();
1183	+
1184	+	// loop over all jet constituents and store them
1185	+	for (unsigned int i=0; i<jconstits.size(); ++i){
1186	+
1187	+	PFParticle part;
1188	+	const reco::PFCandidate* pf = jconstits[i].get();
1189	+
1190	+	// check if it has already been stored, omit if true
1191	+	int is_already_in_list = -1;
1192	+	for (unsigned int j=0; j<NstoredPFs; ++j){
1193	+	PFParticle spf = pfparticles[j];
1194	+	double r2 = pow(pf->eta()-spf.eta(),2) + pow(pf->phi()-spf.phi(),2);
1195	+	double dpt = fabs( pf->pt() - spf.pt() );
1196	+	if (r2<1e-10 && dpt<1e-10){
1197	+	is_already_in_list = j;
1198	+	break;
1199	+	}
1200	+	}
1201	+
1202	+	if (is_already_in_list>-1){
1203	+	continue;
1204	+	}
1205	+
1206	+
1207	+	part.set_pt(pf->pt());
1208	+	part.set_eta(pf->eta());
1209	+	part.set_phi(pf->phi());
1210	+	part.set_energy(pf->energy());
1211	+	part.set_charge(pf->charge());
1212	+
1213	+	part.set_ecal_en(pf->ecalEnergy());
1214	+	part.set_hcal_en(pf->hcalEnergy());
1215	+	reco::TrackRef trackref = pf->trackRef();
1216	+	if (!trackref.isNull()){
1217	+	part.set_track_mom(trackref->p());
1218	+	}
1219	+
1220	+	PFParticle::EParticleID id = PFParticle::eX;
1221	+	switch ( pf->translatePdgIdToType(pf->pdgId()) ){
1222	+	case reco::PFCandidate::X : id = PFParticle::eX; break;
1223	+	case reco::PFCandidate::h : id = PFParticle::eH; break;
1224	+	case reco::PFCandidate::e : id = PFParticle::eE; break;
1225	+	case reco::PFCandidate::mu : id = PFParticle::eMu; break;
1226	+	case reco::PFCandidate::gamma : id = PFParticle::eGamma; break;
1227	+	case reco::PFCandidate::h0 : id = PFParticle::eH0; break;
1228	+	case reco::PFCandidate::h_HF : id = PFParticle::eH_HF; break;
1229	+	case reco::PFCandidate::egamma_HF : id = PFParticle::eEgamma_HF; break;
1230	+	}
1231	+	part.set_particleID(id);
1232	+
1233	+	pfparticles.push_back(part);
1234	+
1235	+	// add a reference to the particle
1236	+	jet->add_pfconstituents_index(pfparticles.size()-1);
1237	+
1238	+	}
1239	+
1240	+	if(pat_jet->isPFJet()){
1241	+	jet->set_charge(pat_jet->charge());
1242	+	jet->set_neutralEmEnergyFraction(pat_jet->neutralEmEnergyFraction());
1243	+	jet->set_neutralHadronEnergyFraction (pat_jet->neutralHadronEnergyFraction());
1244	+	jet->set_chargedEmEnergyFraction (pat_jet->chargedEmEnergyFraction());
1245	+	jet->set_chargedHadronEnergyFraction (pat_jet->chargedHadronEnergyFraction());
1246	+	jet->set_muonEnergyFraction (pat_jet->muonEnergyFraction());
1247	+	jet->set_photonEnergyFraction (pat_jet->photonEnergyFraction());
1248	+	jet->set_chargedMultiplicity (pat_jet->chargedMultiplicity());
1249	+	jet->set_neutralMultiplicity (pat_jet->neutralMultiplicity());
1250	+	jet->set_muonMultiplicity (pat_jet->muonMultiplicity());
1251	+	jet->set_electronMultiplicity (pat_jet->electronMultiplicity());
1252	+	jet->set_photonMultiplicity (pat_jet->photonMultiplicity());
1253	+	}
1254	+
1255	+	return;
1256	+
1257	+	}
1258	+
1259	+	// ------------ method fills PF candidates in a cone of radius R around
1260	+	// ------------ a given particle (lepton, most likely)
1261	+	// ------------ it is checked if the constituents have been stored already
1262	+	void
1263	+	NtupleWriter::StorePFCandsInCone(Particle* inpart, const std::vector<reco::PFCandidate>& pf_cands, double R0)
1264	+	{
1265	+	// checks if the pf cand has already been stored, only stores so far missing ones
1266	+
1267	+	//cout << "\nStorePFCandsInCone: R0 = " << R0 << endl;
1268	+	//cout << "found " << pf_cands.size() << " PF candidates in the event" << endl;
1269	+	//cout << "Input inparticle: pt = " << inpart->pt() << " eta = " << inpart->eta() << " phi = " << inpart->phi() << endl;
1270	+
1271	+	unsigned int NstoredPFs = pfparticles.size();
1272	+	//cout << "got already " <<  NstoredPFs << " which need to be checked to avoid double counting" << endl;
1273	+
1274	+	// loop over all PF candidates and store the ones in a cone with radius R0
1275	+	for (unsigned int i=0; i<pf_cands.size(); ++i){
1276	+
1277	+	reco::PFCandidate pf = pf_cands[i];
1278	+
1279	+	// calculate distance in eta/phi
1280	+	double dphi = pf.phi() - inpart->phi();
1281	+	if (dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
1282	+	if (dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
1283	+	double dr2 = dphi*dphi + pow(pf.eta() - inpart->eta(),2);
1284	+
1285	+	// check if PF candidate is in cone around particle
1286	+	if (sqrt(dr2)>R0) continue;
1287	+
1288	+	//cout << "found particle with distance " << dr2 << " at position " << i << endl;
1289	+	//cout << "PF: pt = " << pf.pt() << " eta = " << pf.eta() << " phi = " << pf.phi() << endl;
1290	+
1291	+	// check if it's the same as the input particle
1292	+	double dpt = fabs( inpart->pt() - pf.pt() );
1293	+	if (dr2<1e-10 && dpt<1e-10){
1294	+	//cout << "same particle, don't store" << endl;
1295	+	continue;
1296	+	}
1297	+
1298	+	// check if it has already been stored, omit if true
1299	+	int is_already_in_list = -1;
1300	+	for (unsigned int j=0; j<NstoredPFs; ++j){
1301	+	PFParticle spf = pfparticles[j];
1302	+	double r2 = pow(pf.eta()-spf.eta(),2) + pow(pf.phi()-spf.phi(),2);
1303	+	double dpt = fabs( pf.pt() - spf.pt() );
1304	+	if (r2<1e-10 && dpt<1e-10){
1305	+	is_already_in_list = j;
1306	+	break;
1307	+	}
1308	+	}
1309	+
1310	+	if (is_already_in_list>-1){
1311	+	//cout << "is already in list, continue" << endl;
1312	+	continue;
1313	+	}
1314	+
1315	+
1316	+	PFParticle part;
1317	+	part.set_pt(pf.pt());
1318	+	part.set_eta(pf.eta());
1319	+	part.set_phi(pf.phi());
1320	+	part.set_energy(pf.energy());
1321	+	part.set_charge(pf.charge());
1322	+
1323	+	part.set_ecal_en(pf.ecalEnergy());
1324	+	part.set_hcal_en(pf.hcalEnergy());
1325	+	reco::TrackRef trackref = pf.trackRef();
1326	+	if (!trackref.isNull()){
1327	+	part.set_track_mom(trackref->p());
1328	+	}
1329	+
1330	+	PFParticle::EParticleID id = PFParticle::eX;
1331	+	switch ( pf.translatePdgIdToType(pf.pdgId()) ){
1332	+	case reco::PFCandidate::X : id = PFParticle::eX; break;
1333	+	case reco::PFCandidate::h : id = PFParticle::eH; break;
1334	+	case reco::PFCandidate::e : id = PFParticle::eE; break;
1335	+	case reco::PFCandidate::mu : id = PFParticle::eMu; break;
1336	+	case reco::PFCandidate::gamma : id = PFParticle::eGamma; break;
1337	+	case reco::PFCandidate::h0 : id = PFParticle::eH0; break;
1338	+	case reco::PFCandidate::h_HF : id = PFParticle::eH_HF; break;
1339	+	case reco::PFCandidate::egamma_HF : id = PFParticle::eEgamma_HF; break;
1340	+	}
1341	+	part.set_particleID(id);
1342	+
1343	+	pfparticles.push_back(part);
1344	+
1345	+	}
1346	+
1347	+	return;
1348	+
1349	+	}
1350	+
1351		//define this as a plug-in
1352		DEFINE_FWK_MODULE(NtupleWriter);

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines