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.5 by peiffer, Thu Apr 5 09:48:20 2012 UTC vs.
Revision 1.27 by rkogler, Wed Jun 19 13:22:06 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");
62		doGenInfo = iConfig.getParameter<bool>("doGenInfo");
63	+	doAllGenParticles = iConfig.getParameter<bool>("doAllGenParticles");
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
72		// initialization of tree variables
73
#	Line 62 | Line 75 | NtupleWriter::NtupleWriter(const edm::Pa
75		tr->Branch("event",&event);
76		tr->Branch("luminosityBlock",&luminosityBlock);
77		tr->Branch("isRealData",&isRealData);
78	<	tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
79	<	tr->Branch("beamspot_x0",&beamspot_x0);
67	<	tr->Branch("beamspot_y0",&beamspot_y0);
68	<	tr->Branch("beamspot_z0",&beamspot_z0);
78	>	tr->Branch("rho",&rho);
79	>	rho_source = iConfig.getParameter<edm::InputTag>("rho_source");
80
81	+	//tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
82	+	if(doLumiInfo){
83	+	tr->Branch("intgRecLumi",&intgRecLumi);
84	+	tr->Branch("intgDelLumi",&intgDelLumi);
85	+	}
86	+	if(doPV){
87	+	tr->Branch("beamspot_x0",&beamspot_x0);
88	+	tr->Branch("beamspot_y0",&beamspot_y0);
89	+	tr->Branch("beamspot_z0",&beamspot_z0);
90	+	}
91		if(doElectrons){
92		electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
93		for(size_t j=0; j< electron_sources.size(); ++j){
#	Line 95 | Line 116 | NtupleWriter::NtupleWriter(const edm::Pa
116		tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
117		}
118		}
119	+	if(doGenJets){
120	+	genjet_sources = iConfig.getParameter<std::vector<std::string> >("genjet_sources");
121	+	genjet_ptmin = iConfig.getParameter<double> ("genjet_ptmin");
122	+	genjet_etamax = iConfig.getParameter<double> ("genjet_etamax");
123	+	for(size_t j=0; j< genjet_sources.size(); ++j){
124	+	tr->Branch( genjet_sources[j].c_str(), "std::vector<Particle>", &genjets[j]);
125	+	}
126	+	}
127		if(doTopJets){
128		topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
129		topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
#	Line 103 | Line 132 | NtupleWriter::NtupleWriter(const edm::Pa
132		tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
133		}
134		}
135	+	if(doTopJetsConstituents){
136	+	topjet_constituents_sources = iConfig.getParameter<std::vector<std::string> >("topjet_constituents_sources");
137	+	tr->Branch( "PFParticles", "std::vector<PFParticle>", &pfparticles);
138	+	}
139	+	if(doGenTopJets){
140	+	gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
141	+	gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
142	+	gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
143	+	for(size_t j=0; j< gentopjet_sources.size(); ++j){
144	+	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<GenTopJet>", &gentopjets[j]);
145	+	}
146	+	}
147		if(doPhotons){
148		photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
149		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 122 | Line 163 | NtupleWriter::NtupleWriter(const edm::Pa
163		}
164		}
165		if(doGenInfo){
166	+	genparticle_source= iConfig.getParameter<edm::InputTag>("genparticle_source");
167		tr->Branch("genInfo","GenInfo",&genInfo);
168		tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
169		}
170	<
171	<	trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
172	<	//tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
173	<	tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
174	<	tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
175	<	tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
176	<	tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
177	<
170	>	if(doTrigger){
171	>	trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
172	>	//tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
173	>	tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
174	>	tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
175	>	//tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
176	>	//tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
177	>	}
178		newrun = true;
179		}
180
#	Line 164 | Line 206 | NtupleWriter::analyze(const edm::Event&
206		luminosityBlock = iEvent.luminosityBlock();
207		isRealData      = iEvent.isRealData();
208
209	<	if(isRealData){
210	<	edm::Handle<bool> bool_handle;
211	<	iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
212	<	HBHENoiseFilterResult = *(bool_handle.product());
213	<	}
214	<	else HBHENoiseFilterResult = false;
209	>	edm::Handle<double> m_rho;
210	>	iEvent.getByLabel(rho_source,m_rho);
211	>	rho=*m_rho;
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;
219
220		// ------------- primary vertices and beamspot  -------------
221
#	Line 185 | Line 231 | NtupleWriter::analyze(const edm::Event&
231		reco::Vertex reco_pv = reco_pvs[i];
232
233		PrimaryVertex pv;
234	<	pv.x =  reco_pv.x();
235	<	pv.y =  reco_pv.y();
236	<	pv.z =  reco_pv.z();
237	<	pv.nTracks =  reco_pv.nTracks();
238	<	//pv.isValid =  reco_pv.isValid();
239	<	pv.chi2 =  reco_pv.chi2();
240	<	pv.ndof =  reco_pv.ndof();
234	>	pv.set_x( reco_pv.x());
235	>	pv.set_y( reco_pv.y());
236	>	pv.set_z( reco_pv.z());
237	>	pv.set_nTracks( reco_pv.nTracks());
238	>	//pv.set_isValid( reco_pv.isValid());
239	>	pv.set_chi2( reco_pv.chi2());
240	>	pv.set_ndof( reco_pv.ndof());
241
242		pvs[j].push_back(pv);
243		}
244		}
245	+
246	+	edm::Handle<reco::BeamSpot> beamSpot;
247	+	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
248	+	const reco::BeamSpot & bsp = *beamSpot;
249	+
250	+	beamspot_x0 = bsp.x0();
251	+	beamspot_y0 = bsp.y0();
252	+	beamspot_z0 = bsp.z0();
253		}
254	+
255	+	// ------------- generator info -------------
256
257	<	edm::Handle<reco::BeamSpot> beamSpot;
258	<	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
259	<	const reco::BeamSpot & bsp = *beamSpot;
260	<
261	<	beamspot_x0 = bsp.x0();
262	<	beamspot_y0 = bsp.y0();
263	<	beamspot_z0 = bsp.z0();
257	>	if(doGenInfo){
258	>	genInfo.clear_weights();
259	>	genInfo.clear_binningValues();
260	>	genps.clear();
261	>
262	>	edm::Handle<GenEventInfoProduct> genEventInfoProduct;
263	>	iEvent.getByLabel("generator", genEventInfoProduct);
264	>	const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
265	>
266	>	for(unsigned int k=0; k<genEventInfo.binningValues().size();++k){
267	>	genInfo.add_binningValue(genEventInfo.binningValues().at(k));
268	>	}
269	>	for(unsigned int k=0; k<genEventInfo.weights().size();++k){
270	>	genInfo.add_weight(genEventInfo.weights().at(k));
271	>	}
272	>	genInfo.set_alphaQCD(genEventInfo.alphaQCD());
273	>	genInfo.set_alphaQED(genEventInfo.alphaQED());
274	>	genInfo.set_qScale(genEventInfo.qScale());
275	>
276	>	const gen::PdfInfo* pdf = genEventInfo.pdf();
277	>	if(pdf){
278	>	genInfo.set_pdf_id1(pdf->id.first);
279	>	genInfo.set_pdf_id2(pdf->id.second);
280	>	genInfo.set_pdf_x1(pdf->x.first);
281	>	genInfo.set_pdf_x2(pdf->x.second);
282	>	genInfo.set_pdf_scalePDF(pdf->scalePDF);
283	>	genInfo.set_pdf_xPDF1(pdf->xPDF.first);
284	>	genInfo.set_pdf_xPDF2(pdf->xPDF.second);
285	>	}
286	>	else{
287	>	genInfo.set_pdf_id1(-999);
288	>	genInfo.set_pdf_id2(-999);
289	>	genInfo.set_pdf_x1(-999);
290	>	genInfo.set_pdf_x2(-999);
291	>	genInfo.set_pdf_scalePDF(-999);
292	>	genInfo.set_pdf_xPDF1(-999);
293	>	genInfo.set_pdf_xPDF2(-999);
294	>	}
295	>
296	>	edm::Handle<std::vector<PileupSummaryInfo> > pus;
297	>	iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
298	>	genInfo.set_pileup_NumInteractions_intime(0);
299	>	genInfo.set_pileup_NumInteractions_ootbefore(0);
300	>	genInfo.set_pileup_NumInteractions_ootafter(0);
301	>	if(pus.isValid()){
302	>	genInfo.set_pileup_TrueNumInteractions ( (float) pus->at(0).getTrueNumInteractions());
303	>	for(size_t i=0; i<pus->size(); ++i){
304	>	if(pus->at(i).getBunchCrossing() == 0) // intime pileup
305	>	genInfo.set_pileup_NumInteractions_intime( genInfo.pileup_NumInteractions_intime() + pus->at(i).getPU_NumInteractions());
306	>	else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
307	>	genInfo.set_pileup_NumInteractions_ootbefore( genInfo.pileup_NumInteractions_ootbefore() + pus->at(i).getPU_NumInteractions());
308	>	}
309	>	else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
310	>	genInfo.set_pileup_NumInteractions_ootafter( genInfo.pileup_NumInteractions_ootafter() + pus->at(i).getPU_NumInteractions());
311	>	}
312	>	}
313	>	}
314	>
315	>	edm::Handle<reco::GenParticleCollection> genPartColl;
316	>	iEvent.getByLabel(genparticle_source, genPartColl);
317	>	int index=-1;
318	>	for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
319	>	index++;
320	>
321	>	//write out only top quarks,final state leptons and status 3 particles (works fine only for MadGraph)
322	>	bool islepton = iter->status()==1 && abs(iter->pdgId())>=11 && abs(iter->pdgId())<=16 ;
323	>	if(abs(iter->pdgId())==6 || iter->status()==3 || islepton ||  doAllGenParticles){
324	>	GenParticle genp;
325	>	genp.set_charge(iter->charge());
326	>	genp.set_pt(iter->p4().pt());
327	>	genp.set_eta(iter->p4().eta());
328	>	genp.set_phi(iter->p4().phi());
329	>	genp.set_energy(iter->p4().E());
330	>	genp.set_index(index);
331	>	genp.set_status( iter->status());
332	>	genp.set_pdgId( iter->pdgId());
333	>
334	>	genp.set_mother1(-1);
335	>	genp.set_mother2(-1);
336	>	genp.set_daughter1(-1);
337	>	genp.set_daughter2(-1);
338	>
339	>	int nm=iter->numberOfMothers();
340	>	int nd=iter->numberOfDaughters();
341	>
342	>
343	>	if (nm>0) genp.set_mother1( iter->motherRef(0).key());
344	>	if (nm>1) genp.set_mother2( iter->motherRef(1).key());
345	>	if (nd>0) genp.set_daughter1( iter->daughterRef(0).key());
346	>	if (nd>1) genp.set_daughter2( iter->daughterRef(1).key());
347	>
348	>	genps.push_back(genp);
349	>	}
350	>	}
351	>	}
352
353		// ------------- electrons -------------
354		if(doElectrons){
355	+
356	+	//      edm::Handle<reco::ConversionCollection> hConversions;
357	+	//      iEvent.getByLabel("allConversions", hConversions);
358	+
359	+	//      edm::Handle<reco::BeamSpot> beamSpot;
360	+	//      iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
361	+	//      const reco::BeamSpot & bsp = *beamSpot;
362	+
363		for(size_t j=0; j< electron_sources.size(); ++j){
364		eles[j].clear();
365		edm::Handle< std::vector<pat::Electron> > ele_handle;
#	Line 218 | Line 370 | NtupleWriter::analyze(const edm::Event&
370		pat::Electron pat_ele = pat_electrons[i];
371		Electron ele;
372
373	<	ele.charge =  pat_ele.charge();
374	<	ele.pt =  pat_ele.pt();
375	<	ele.eta =  pat_ele.eta();
376	<	ele.phi =  pat_ele.phi();
377	<	ele.energy =  pat_ele.energy();
378	<	ele.vertex_x = pat_ele.vertex().x();
379	<	ele.vertex_y = pat_ele.vertex().y();
380	<	ele.vertex_z = pat_ele.vertex().z();
381	<	ele.supercluster_eta = pat_ele.superCluster()->eta();
382	<	ele.supercluster_phi = pat_ele.superCluster()->phi();
383	<	ele.dB = pat_ele.dB();
384	<	//ele.particleIso = pat_ele.particleIso();
385	<	ele.neutralHadronIso = pat_ele.neutralHadronIso();
386	<	ele.chargedHadronIso = pat_ele.chargedHadronIso();
387	<	ele.trackIso = pat_ele.trackIso();
388	<	ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
389	<	ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
390	<	ele.gsfTrack_px= pat_ele.gsfTrack()->px();
391	<	ele.gsfTrack_py= pat_ele.gsfTrack()->py();
392	<	ele.gsfTrack_pz= pat_ele.gsfTrack()->pz();
393	<	ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
394	<	ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
395	<	ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
373	>	ele.set_charge( pat_ele.charge());
374	>	ele.set_pt( pat_ele.pt());
375	>	ele.set_eta( pat_ele.eta());
376	>	ele.set_phi( pat_ele.phi());
377	>	ele.set_energy( pat_ele.energy());
378	>	ele.set_vertex_x(pat_ele.vertex().x());
379	>	ele.set_vertex_y(pat_ele.vertex().y());
380	>	ele.set_vertex_z(pat_ele.vertex().z());
381	>	ele.set_supercluster_eta(pat_ele.superCluster()->eta());
382	>	ele.set_supercluster_phi(pat_ele.superCluster()->phi());
383	>	ele.set_dB(pat_ele.dB());
384	>	//ele.set_particleIso(pat_ele.particleIso());
385	>	ele.set_neutralHadronIso(pat_ele.neutralHadronIso());
386	>	ele.set_chargedHadronIso(pat_ele.chargedHadronIso());
387	>	ele.set_trackIso(pat_ele.trackIso());
388	>	ele.set_photonIso(pat_ele.photonIso());
389	>	ele.set_puChargedHadronIso(pat_ele.puChargedHadronIso());
390	>	ele.set_gsfTrack_trackerExpectedHitsInner_numberOfLostHits(pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits());
391	>	ele.set_gsfTrack_px( pat_ele.gsfTrack()->px());
392	>	ele.set_gsfTrack_py( pat_ele.gsfTrack()->py());
393	>	ele.set_gsfTrack_pz( pat_ele.gsfTrack()->pz());
394	>	ele.set_gsfTrack_vx( pat_ele.gsfTrack()->vx());
395	>	ele.set_gsfTrack_vy( pat_ele.gsfTrack()->vy());
396	>	ele.set_gsfTrack_vz( pat_ele.gsfTrack()->vz());
397	>	//ele.set_passconversionveto(!ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position()));
398	>	ele.set_passconversionveto(pat_ele.passConversionVeto());
399	>	ele.set_dEtaIn(pat_ele.deltaEtaSuperClusterTrackAtVtx());
400	>	ele.set_dPhiIn(pat_ele.deltaPhiSuperClusterTrackAtVtx());
401	>	ele.set_sigmaIEtaIEta(pat_ele.sigmaIetaIeta());
402	>	ele.set_HoverE(pat_ele.hadronicOverEm());
403	>	ele.set_fbrem(pat_ele.fbrem());
404	>	ele.set_EoverPIn(pat_ele.eSuperClusterOverP());
405	>	ele.set_EcalEnergy(pat_ele.ecalEnergy());
406	>	ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
407	>	ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
408	>	float AEff03 = 0.00;
409	>	if(isRealData){
410	>	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2011);
411	>	}else{
412	>	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAFall11MC);
413	>	}
414	>	ele.set_AEff(AEff03);
415	>
416		eles[j].push_back(ele);
417		}
418		}
#	Line 259 | Line 431 | NtupleWriter::analyze(const edm::Event&
431		pat::Muon pat_mu = pat_muons[i];
432
433		Muon mu;
434	<	mu.charge =  pat_mu.charge();
435	<	mu.pt =  pat_mu.pt();
436	<	mu.eta =  pat_mu.eta();
437	<	mu.phi =  pat_mu.phi();
438	<	mu.energy =  pat_mu.energy();
439	<	mu.vertex_x = pat_mu.vertex().x();
440	<	mu.vertex_y = pat_mu.vertex().y();
441	<	mu.vertex_z = pat_mu.vertex().z();
442	<	mu.dB = pat_mu.dB();
443	<	//mu.particleIso = pat_mu.particleIso();
444	<	mu.neutralHadronIso = pat_mu.neutralHadronIso();
445	<	mu.chargedHadronIso = pat_mu.chargedHadronIso();
446	<	mu.trackIso = pat_mu.trackIso();
447	<	mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
448	<	mu.isGlobalMuon = pat_mu.isGlobalMuon();
449	<	mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
450	<	mu.isTrackerMuon = pat_mu.isTrackerMuon();
451	<	mu.numberOfMatchedStations = pat_mu.numberOfMatchedStations();
434	>	mu.set_charge( pat_mu.charge());
435	>	mu.set_pt( pat_mu.pt());
436	>	mu.set_eta( pat_mu.eta());
437	>	mu.set_phi( pat_mu.phi());
438	>	mu.set_energy( pat_mu.energy());
439	>	mu.set_vertex_x ( pat_mu.vertex().x());
440	>	mu.set_vertex_y ( pat_mu.vertex().y());
441	>	mu.set_vertex_z ( pat_mu.vertex().z());
442	>	mu.set_dB ( pat_mu.dB());
443	>	//mu.particleIso ( pat_mu.particleIso());
444	>	mu.set_neutralHadronIso ( pat_mu.neutralHadronIso());
445	>	mu.set_chargedHadronIso ( pat_mu.chargedHadronIso());
446	>	mu.set_trackIso ( pat_mu.trackIso());
447	>	mu.set_photonIso ( pat_mu.photonIso());
448	>	mu.set_puChargedHadronIso ( pat_mu.puChargedHadronIso());
449	>	mu.set_isGlobalMuon ( pat_mu.isGlobalMuon());
450	>	mu.set_isPFMuon ( pat_mu.isPFMuon());
451	>	mu.set_isStandAloneMuon ( pat_mu.isStandAloneMuon());
452	>	mu.set_isTrackerMuon ( pat_mu.isTrackerMuon());
453	>	mu.set_numberOfMatchedStations ( pat_mu.numberOfMatchedStations());
454		reco::TrackRef globalTrack = pat_mu.globalTrack();
455		if(!globalTrack.isNull()){
456	<	mu.globalTrack_chi2 = globalTrack->chi2();
457	<	mu.globalTrack_ndof = globalTrack->ndof();
458	<	mu.globalTrack_d0 = globalTrack->d0();
459	<	mu.globalTrack_d0Error = globalTrack->d0Error();
460	<	mu.globalTrack_numberOfValidHits = globalTrack->numberOfValidHits();
461	<	mu.globalTrack_numberOfLostHits = globalTrack->numberOfLostHits();
456	>	mu.set_globalTrack_chi2 ( globalTrack->chi2());
457	>	mu.set_globalTrack_ndof ( globalTrack->ndof());
458	>	mu.set_globalTrack_d0 ( globalTrack->d0());
459	>	mu.set_globalTrack_d0Error ( globalTrack->d0Error());
460	>	mu.set_globalTrack_numberOfValidHits ( globalTrack->numberOfValidHits());
461	>	mu.set_globalTrack_numberOfLostHits ( globalTrack->numberOfLostHits());
462	>	mu.set_globalTrack_numberOfValidMuonHits(globalTrack->hitPattern().numberOfValidMuonHits() );
463		}
464		else{
465	<	mu.globalTrack_chi2 = 0;
466	<	mu.globalTrack_ndof = 0;
467	<	mu.globalTrack_d0 = 0;
468	<	mu.globalTrack_d0Error = 0;
469	<	mu.globalTrack_numberOfValidHits = 0;
470	<	mu.globalTrack_numberOfLostHits = 0;
465	>	mu.set_globalTrack_chi2 ( 0);
466	>	mu.set_globalTrack_ndof ( 0);
467	>	mu.set_globalTrack_d0 ( 0);
468	>	mu.set_globalTrack_d0Error ( 0);
469	>	mu.set_globalTrack_numberOfValidHits ( 0);
470	>	mu.set_globalTrack_numberOfLostHits ( 0);
471		}
472		reco::TrackRef innerTrack = pat_mu.innerTrack();
473		if(!innerTrack.isNull()){
474	<	mu.innerTrack_chi2 = innerTrack->chi2();
475	<	mu.innerTrack_ndof = innerTrack->ndof();
476	<	mu.innerTrack_d0 = innerTrack->d0();
477	<	mu.innerTrack_d0Error = innerTrack->d0Error();
478	<	mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
479	<	mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
474	>	mu.set_innerTrack_chi2 ( innerTrack->chi2());
475	>	mu.set_innerTrack_ndof ( innerTrack->ndof());
476	>	mu.set_innerTrack_d0 ( innerTrack->d0());
477	>	mu.set_innerTrack_d0Error ( innerTrack->d0Error());
478	>	mu.set_innerTrack_numberOfValidHits ( innerTrack->numberOfValidHits());
479	>	mu.set_innerTrack_numberOfLostHits ( innerTrack->numberOfLostHits());
480	>	mu.set_innerTrack_trackerLayersWithMeasurement ( innerTrack->hitPattern().trackerLayersWithMeasurement());
481	>	mu.set_innerTrack_numberOfValidPixelHits ( innerTrack->hitPattern().numberOfValidPixelHits());
482		}
483		else{
484	<	mu.innerTrack_chi2 = 0;
485	<	mu.innerTrack_ndof = 0;
486	<	mu.innerTrack_d0 = 0;
487	<	mu.innerTrack_d0Error = 0;
488	<	mu.innerTrack_numberOfValidHits = 0;
489	<	mu.innerTrack_numberOfLostHits = 0;
484	>	mu.set_innerTrack_chi2 ( 0);
485	>	mu.set_innerTrack_ndof ( 0);
486	>	mu.set_innerTrack_d0 ( 0);
487	>	mu.set_innerTrack_d0Error ( 0);
488	>	mu.set_innerTrack_numberOfValidHits ( 0);
489	>	mu.set_innerTrack_numberOfLostHits ( 0);
490	>	mu.set_innerTrack_trackerLayersWithMeasurement ( 0);
491	>	mu.set_innerTrack_numberOfValidPixelHits ( 0);
492		}
493		reco::TrackRef outerTrack = pat_mu.outerTrack();
494		if(!outerTrack.isNull()){
495	<	mu.outerTrack_chi2 = outerTrack->chi2();
496	<	mu.outerTrack_ndof = outerTrack->ndof();
497	<	mu.outerTrack_d0 = outerTrack->d0();
498	<	mu.outerTrack_d0Error = outerTrack->d0Error();
499	<	mu.outerTrack_numberOfValidHits = outerTrack->numberOfValidHits();
500	<	mu.outerTrack_numberOfLostHits = outerTrack->numberOfLostHits();
495	>	mu.set_outerTrack_chi2 ( outerTrack->chi2());
496	>	mu.set_outerTrack_ndof ( outerTrack->ndof());
497	>	mu.set_outerTrack_d0 ( outerTrack->d0());
498	>	mu.set_outerTrack_d0Error ( outerTrack->d0Error());
499	>	mu.set_outerTrack_numberOfValidHits ( outerTrack->numberOfValidHits());
500	>	mu.set_outerTrack_numberOfLostHits ( outerTrack->numberOfLostHits());
501		}
502		else{
503	<	mu.outerTrack_chi2 = 0;
504	<	mu.outerTrack_ndof = 0;
505	<	mu.outerTrack_d0 = 0;
506	<	mu.outerTrack_d0Error = 0;
507	<	mu.outerTrack_numberOfValidHits = 0;
508	<	mu.outerTrack_numberOfLostHits = 0;
503	>	mu.set_outerTrack_chi2 ( 0);
504	>	mu.set_outerTrack_ndof ( 0);
505	>	mu.set_outerTrack_d0 ( 0);
506	>	mu.set_outerTrack_d0Error ( 0);
507	>	mu.set_outerTrack_numberOfValidHits ( 0);
508	>	mu.set_outerTrack_numberOfLostHits ( 0);
509		}
510
511		mus[j].push_back(mu);
#	Line 350 | Line 529 | NtupleWriter::analyze(const edm::Event&
529		if(fabs(pat_tau.eta()) > tau_etamax) continue;
530
531		Tau tau;
532	<	tau.charge =  pat_tau.charge();
533	<	tau.pt =  pat_tau.pt();
534	<	tau.eta =  pat_tau.eta();
535	<	tau.phi =  pat_tau.phi();
536	<	tau.energy =  pat_tau.energy();
532	>	tau.set_charge( pat_tau.charge());
533	>	tau.set_pt( pat_tau.pt());
534	>	tau.set_eta( pat_tau.eta());
535	>	tau.set_phi( pat_tau.phi());
536	>	tau.set_energy( pat_tau.energy());
537	>	tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
538	>	//tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
539	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
540	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
541	>	tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
542	>	tau.set_byLooseIsolationMVA( pat_tau.tauID("byLooseIsolationMVA")>0.5);
543	>	tau.set_byMediumIsolationMVA( pat_tau.tauID("byMediumIsolationMVA")>0.5);
544	>	tau.set_byTightIsolationMVA( pat_tau.tauID("byTightIsolationMVA")>0.5);
545	>	tau.set_byLooseIsolationMVA2( pat_tau.tauID("byLooseIsolationMVA2")>0.5);
546	>	tau.set_byMediumIsolationMVA2( pat_tau.tauID("byMediumIsolationMVA2")>0.5);
547	>	tau.set_byTightIsolationMVA2( pat_tau.tauID("byTightIsolationMVA2")>0.5);
548	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr3Hits(  pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits")>0.5);
549	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits")>0.5);
550	>	tau.set_byTightCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr3Hits")>0.5);
551	>	tau.set_againstElectronLooseMVA3  ( pat_tau.tauID("againstElectronLooseMVA3")>0.5);
552	>	tau.set_againstElectronMediumMVA3 ( pat_tau.tauID("againstElectronMediumMVA3")>0.5);
553	>	tau.set_againstElectronTightMVA3 ( pat_tau.tauID("againstElectronTightMVA3")>0.5);
554	>	tau.set_againstElectronVTightMVA3 ( pat_tau.tauID("againstElectronVTightMVA3")>0.5);
555	>	tau.set_againstMuonLoose2 ( pat_tau.tauID("againstMuonLoose2")>0.5);
556	>	tau.set_againstMuonMedium2 ( pat_tau.tauID("againstMuonMedium2")>0.5);
557	>	tau.set_againstMuonTight2 ( pat_tau.tauID("againstMuonTight2")>0.5);
558	>	tau.set_byIsolationMVAraw(  pat_tau.tauID("byIsolationMVAraw"));
559	>	tau.set_byIsolationMVA2raw(  pat_tau.tauID("byIsolationMVA2raw"));
560	>	tau.set_decayMode( pat_tau.decayMode() );
561	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw"));
562	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw3Hits( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));
563
564	+	//       std::cout << pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") << std::endl;
565	+
566	+	//       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
567	+	//       if(!leadPFCand.isNull()){
568	+	//         tau.set_leadPFCand_px ( leadPFCand->px());
569	+	//         tau.set_leadPFCand_py ( leadPFCand->py());
570	+	//         tau.set_leadPFCand_pz ( leadPFCand->pz());
571	+	//       }
572	+	//       else{
573	+	//         tau.set_leadPFCand_px ( 0);
574	+	//         tau.set_leadPFCand_py ( 0);
575	+	//         tau.set_leadPFCand_pz ( 0);
576	+	//       }
577		taus[j].push_back(tau);
578		}
579		}
580		}
581
582	+	//-------------- gen jets -------------
583	+
584	+	if(doGenJets){
585	+	for(size_t j=0; j< genjet_sources.size(); ++j){
586	+
587	+	genjets[j].clear();
588	+
589	+	edm::Handle< std::vector<reco::GenJet> > genjet_handle;
590	+	iEvent.getByLabel(genjet_sources[j], genjet_handle);
591	+	const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
592	+
593	+	for (unsigned int i = 0; i < gen_jets.size(); ++i) {
594	+
595	+	const reco::GenJet* gen_jet = &gen_jets[i];
596	+	if(gen_jet->pt() < genjet_ptmin) continue;
597	+	if(fabs(gen_jet->eta()) > genjet_etamax) continue;
598	+
599	+	Particle jet;
600	+	jet.set_charge(gen_jet->charge());
601	+	jet.set_pt(gen_jet->pt());
602	+	jet.set_eta(gen_jet->eta());
603	+	jet.set_phi(gen_jet->phi());
604	+	jet.set_energy(gen_jet->energy());
605	+
606	+	// recalculate the jet charge
607	+	int jet_charge = 0;
608	+	std::vector<const reco::GenParticle * > jetgenps = gen_jet->getGenConstituents();
609	+	for(unsigned int l = 0; l<jetgenps.size(); ++l){
610	+	jet_charge +=  jetgenps[l]->charge();
611	+	}
612	+	jet.set_charge(jet_charge);
613	+
614	+	genjets[j].push_back(jet);
615	+	}
616	+	}
617	+	}
618	+
619		// ------------- jets -------------
620		if(doJets){
621		for(size_t j=0; j< jet_sources.size(); ++j){
#	Line 382 | Line 637 | NtupleWriter::analyze(const edm::Event&
637		//       }
638
639		Jet jet;
640	<	jet.charge = pat_jet.charge();
641	<	jet.pt = pat_jet.pt();
642	<	jet.eta = pat_jet.eta();
643	<	jet.phi = pat_jet.phi();
644	<	jet.energy = pat_jet.energy();
645	<	jet.numberOfDaughters =pat_jet.numberOfDaughters();
640	>	jet.set_charge(pat_jet.charge());
641	>	jet.set_pt(pat_jet.pt());
642	>	jet.set_eta(pat_jet.eta());
643	>	jet.set_phi(pat_jet.phi());
644	>	jet.set_energy(pat_jet.energy());
645	>	jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
646		const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
647	<	jet.nTracks = jettracks.size();
648	<	jet.jetArea = pat_jet.jetArea();
649	<	jet.pileup = pat_jet.pileup();
647	>	jet.set_nTracks ( jettracks.size());
648	>	jet.set_jetArea(pat_jet.jetArea());
649	>	jet.set_pileup(pat_jet.pileup());
650		if(pat_jet.isPFJet()){
651	<	jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
652	<	jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
653	<	jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
654	<	jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
655	<	jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
656	<	jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
657	<	jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
658	<	jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
659	<	jet.muonMultiplicity =pat_jet.muonMultiplicity();
660	<	jet.electronMultiplicity =pat_jet.electronMultiplicity();
661	<	jet.photonMultiplicity =pat_jet.photonMultiplicity();
662	<	}
408	<
409	<	jecUnc->setJetEta(pat_jet.eta());
410	<	jecUnc->setJetPt(pat_jet.pt());
411	<	jet.JEC_uncertainty = jecUnc->getUncertainty(true);
412	<
413	<	jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
414	<	jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
415	<	jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
416	<	jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
417	<	jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
418	<	jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
651	>	jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
652	>	jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
653	>	jet.set_chargedEmEnergyFraction (pat_jet.chargedEmEnergyFraction());
654	>	jet.set_chargedHadronEnergyFraction (pat_jet.chargedHadronEnergyFraction());
655	>	jet.set_muonEnergyFraction (pat_jet.muonEnergyFraction());
656	>	jet.set_photonEnergyFraction (pat_jet.photonEnergyFraction());
657	>	jet.set_chargedMultiplicity (pat_jet.chargedMultiplicity());
658	>	jet.set_neutralMultiplicity (pat_jet.neutralMultiplicity());
659	>	jet.set_muonMultiplicity (pat_jet.muonMultiplicity());
660	>	jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
661	>	jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
662	>	}
663
664	+	jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
665	+
666	+	jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
667	+	jet.set_btag_simpleSecondaryVertexHighPur(pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
668	+	jet.set_btag_combinedSecondaryVertex(pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags"));
669	+	jet.set_btag_combinedSecondaryVertexMVA(pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
670	+	jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
671	+	jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
672	+
673	+
674	+	const reco::GenJet *genj = pat_jet.genJet();
675	+	if(genj){
676	+
677	+	for(unsigned int k=0; k<genjets->size(); ++k){
678	+	if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta()){
679	+	jet.set_genjet_index(k);
680	+	}
681	+	}
682	+	//         if( jet.genjet_index()<0){
683	+	//           std::cout<< "genjet not found for " << genj->pt() << "  " << genj->eta() << std::endl;
684	+	//         }
685	+
686	+	if(doAllGenParticles){
687	+	std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
688	+	for(unsigned int l = 0; l<jetgenps.size(); ++l){
689	+	for(unsigned int k=0; k< genps.size(); ++k){
690	+	if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
691	+	jet.add_genparticles_index(genps[k].index());
692	+	break;
693	+	}
694	+	}
695	+	}
696	+	if(jet.genparticles_indices().size()!= jetgenps.size())
697	+	std::cout << "WARNING: Found only " << jet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
698	+	}
699	+
700	+	}
701	+
702		jets[j].push_back(jet);
703		}
704		}
705		}
706
707	+
708		// ------------- top jets -------------
709		if(doTopJets){
710	+
711		for(size_t j=0; j< topjet_sources.size(); ++j){
712
713		topjets[j].clear();
#	Line 433 | Line 717 | NtupleWriter::analyze(const edm::Event&
717		iEvent.getByLabel(topjet_sources[j], pat_topjets);
718
719		for (unsigned int i = 0; i < pat_topjets->size(); i++) {
720	+
721		const pat::Jet  pat_topjet =  * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
722		if(pat_topjet.pt() < topjet_ptmin) continue;
723		if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
724
725		TopJet topjet;
726	<	topjet.charge = pat_topjet.charge();
727	<	topjet.pt = pat_topjet.pt();
728	<	topjet.eta = pat_topjet.eta();
729	<	topjet.phi = pat_topjet.phi();
730	<	topjet.energy = pat_topjet.energy();
731	<	topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
726	>	topjet.set_charge(pat_topjet.charge());
727	>	topjet.set_pt(pat_topjet.pt());
728	>	topjet.set_eta(pat_topjet.eta());
729	>	topjet.set_phi(pat_topjet.phi());
730	>	topjet.set_energy(pat_topjet.energy());
731	>	topjet.set_numberOfDaughters(pat_topjet.numberOfDaughters());
732		const reco::TrackRefVector&  topjettracks = pat_topjet.associatedTracks();
733	<	topjet.nTracks = topjettracks.size();
734	<	topjet.jetArea = pat_topjet.jetArea();
735	<	topjet.pileup = pat_topjet.pileup();
736	<	//       topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
737	<	//       topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
738	<	//       topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
739	<	//       topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
740	<	//       topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
741	<	//       topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
742	<	//       topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
743	<	//       topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
744	<	//       topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
745	<	//       topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
746	<	//       topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
747	<
748	<	jecUnc->setJetEta(pat_topjet.eta());
749	<	jecUnc->setJetPt(pat_topjet.pt());
750	<	topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
751	<
752	<	topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
753	<	topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
754	<	topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
755	<	topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
756	<	topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
757	<	topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
733	>	topjet.set_nTracks( topjettracks.size());
734	>	topjet.set_jetArea( pat_topjet.jetArea());
735	>	topjet.set_pileup( pat_topjet.pileup());
736	>
737	>	topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
738	>
739	>	topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
740	>	topjet.set_btag_simpleSecondaryVertexHighPur(pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
741	>	topjet.set_btag_combinedSecondaryVertex(pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags"));
742	>	topjet.set_btag_combinedSecondaryVertexMVA(pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
743	>	topjet.set_btag_jetBProbability(pat_topjet.bDiscriminator("jetBProbabilityBJetTags"));
744	>	topjet.set_btag_jetProbability(pat_topjet.bDiscriminator("jetProbabilityBJetTags"));
745	>
746	>	/*
747	>	const reco::GenJet *genj = pat_topjet.genJet();
748	>	if(genj){
749	>	topjet.set_genjet_pt ( genj->pt());
750	>	topjet.set_genjet_eta ( genj->eta());
751	>	topjet.set_genjet_phi ( genj->phi());
752	>	topjet.set_genjet_energy ( genj->energy());
753	>	if(doAllGenParticles){
754	>	std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
755	>	for(unsigned int l = 0; l<jetgenps.size(); ++l){
756	>	for(unsigned int k=0; k< genps.size(); ++k){
757	>	if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
758	>	topjet.add_genparticles_index(genps[k].index());
759	>	break;
760	>	}
761	>	}
762	>	}
763	>	if(topjet.genparticles_indices().size()!= jetgenps.size())
764	>	std::cout << "WARNING: Found only " << topjet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
765	>	}
766	>	}
767	>	*/
768	>
769	>	// add constituents to the jet, if requested
770	>	if (doTopJetsConstituents){
771	>
772	>	if (topjet_constituents_sources.size()>j){ //only add constituents if they are defined
773	>
774	>	edm::Handle<pat::JetCollection> pat_topjets_with_cands;
775	>	iEvent.getByLabel(topjet_constituents_sources[j], pat_topjets_with_cands);
776	>	pat::Jet* pat_topjet_wc = NULL;
777	>
778	>	for (unsigned int it = 0; it < pat_topjets_with_cands->size(); it++) {
779	>	const pat::Jet* cand =  dynamic_cast<pat::Jet const *>(&pat_topjets_with_cands->at(it));
780	>	double dphi = cand->phi() - pat_topjet.phi();
781	>	if (dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
782	>	if (dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
783	>	if (fabs(dphi)<0.5 && fabs(cand->eta()-pat_topjet.eta())<0.5){ // be generous: filtering, pruning... can change jet axis
784	>	pat_topjet_wc = const_cast<pat::Jet*>(cand);
785	>	break;
786	>	}
787	>	}
788	>
789	>	if (pat_topjet_wc){
790	>	StoreJetConstituents(pat_topjet_wc, &topjet);
791	>
792	>	// now run substructure information
793	>	JetProps substructure(&topjet);
794	>	substructure.set_pf_cands(&pfparticles);
795	>	double tau1 = substructure.GetNsubjettiness(1, Njettiness::onepass_kt_axes, 1., 2.0);
796	>	double tau2 = substructure.GetNsubjettiness(2, Njettiness::onepass_kt_axes, 1., 2.0);
797	>	double tau3 = substructure.GetNsubjettiness(3, Njettiness::onepass_kt_axes, 1., 2.0);
798	>	double qjets = substructure.GetQjetVolatility(iEvent.id().event(), 2.0);
799	>	topjet.set_tau1(tau1);
800	>	topjet.set_tau2(tau2);
801	>	topjet.set_tau3(tau3);
802	>	topjet.set_qjets_volatility(qjets);
803	>
804	>	}
805	>	}
806	>	}
807	>
808	>
809
810		for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
811		Particle subjet_v4;
812	<	subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
813	<	subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
814	<	subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
815	<	subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
816	<	topjet.subjets.push_back(subjet_v4);
812	>
813	>	reco::Candidate const * subjetd =  pat_topjet.daughter(k);
814	>	pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
815	>	if(patsubjetd)
816	>	{
817	>	subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
818	>	subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
819	>	subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
820	>	subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
821	>	topjet.add_subjet(subjet_v4);
822	>	//btag info
823	>	topjet.add_subFlavour(patsubjetd->partonFlavour());
824	>	topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
825	>	if (doTagInfos)
826	>	{
827	>	//ip tag info
828	>	reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
829	>	topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
830	>	topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
831	>	topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
832	>	topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
833	>	topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
834	>	topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
835	>	topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
836	>	topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
837	>	topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
838	>	topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
839	>	topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
840	>	topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));
841	>	topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
842	>	topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
843	>	topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
844	>	topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
845	>	topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
846	>	topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
847	>	topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
848	>	topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
849	>	topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
850	>	//sv tag info
851	>	reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
852	>	topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
853	>	topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
854	>	topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
855	>	topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
856	>	topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
857	>	topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
858	>	topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
859	>	std::vector<TLorentzVector> vp4; vp4.clear();
860	>	std::vector<float> vchi2; vchi2.clear();
861	>	std::vector<float> vndof; vndof.clear();
862	>	std::vector<float> vchi2ndof; vchi2ndof.clear();
863	>	std::vector<float> vtrsize; vtrsize.clear();
864	>	for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
865	>	{
866	>	ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
867	>	vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
868	>	vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());
869	>	vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());
870	>	vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());
871	>	vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());
872	>	}
873	>	topjet.add_subSecondaryVertex(vp4);
874	>	topjet.add_subVertexChi2(vchi2);
875	>	topjet.add_subVertexNdof(vndof);
876	>	topjet.add_subVertexNormalizedChi2(vchi2ndof);
877	>	topjet.add_subVertexTracksSize(vtrsize);
878	>	//try computer
879	>	edm::ESHandle<JetTagComputer> computerHandle;
880	>	iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
881	>	const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
882	>	if(computer)
883	>	{
884	>	computer->passEventSetup(iSetup);
885	>	std::vector<const reco::BaseTagInfo*>  baseTagInfos;
886	>	baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
887	>	baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );
888	>	JetTagComputer::TagInfoHelper helper(baseTagInfos);
889	>	reco::TaggingVariableList vars = computer->taggingVariables(helper);
890	>	topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
891	>	topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
892	>	topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
893	>	topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
894	>	}
895	>	}
896	>	}
897	>	else
898	>	{
899	>	//filling only standard information in case the subjet has not been pat-tified during the pattuples production
900	>	subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
901	>	subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
902	>	subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
903	>	subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
904	>	topjet.add_subjet(subjet_v4);
905	>	}
906	>
907	>
908		}
909		topjets[j].push_back(topjet);
910		}
911		}
912		}
913
914	+
915	+	// ------------- generator top jets -------------
916	+	if(doGenTopJets){
917	+	for(size_t j=0; j< gentopjet_sources.size(); ++j){
918	+
919	+	gentopjets[j].clear();
920	+
921	+	edm::Handle<reco::BasicJetCollection> reco_gentopjets;
922	+	iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
923	+
924	+	for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
925	+
926	+	const reco::BasicJet  reco_gentopjet =  reco_gentopjets->at(i);
927	+	if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
928	+	if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
929	+
930	+	GenTopJet gentopjet;
931	+	gentopjet.set_charge(reco_gentopjet.charge());
932	+	gentopjet.set_pt(reco_gentopjet.pt());
933	+	gentopjet.set_eta(reco_gentopjet.eta());
934	+	gentopjet.set_phi(reco_gentopjet.phi());
935	+	gentopjet.set_energy(reco_gentopjet.energy());
936	+
937	+	for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
938	+	Particle subjet_v4;
939	+	subjet_v4.set_pt(reco_gentopjet.daughter(k)->p4().pt());
940	+	subjet_v4.set_eta(reco_gentopjet.daughter(k)->p4().eta());
941	+	subjet_v4.set_phi(reco_gentopjet.daughter(k)->p4().phi());
942	+	subjet_v4.set_energy(reco_gentopjet.daughter(k)->p4().E());
943	+	gentopjet.add_subjet(subjet_v4);
944	+	}
945	+	gentopjets[j].push_back(gentopjet);
946	+	}
947	+	}
948	+	}
949	+
950	+
951		// ------------- photons -------------
952		if(doPhotons){
953		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 496 | Line 960 | NtupleWriter::analyze(const edm::Event&
960		for (unsigned int i = 0; i < pat_photons.size(); ++i) {
961		pat::Photon pat_photon = pat_photons[i];
962		Photon ph;
963	<	ph.charge = 0;
964	<	ph.pt =  pat_photon.pt();
965	<	ph.eta =  pat_photon.eta();
966	<	ph.phi =  pat_photon.phi();
967	<	ph.energy =  pat_photon.energy();
968	<	ph.vertex_x = pat_photon.vertex().x();
969	<	ph.vertex_y = pat_photon.vertex().y();
970	<	ph.vertex_z = pat_photon.vertex().z();
971	<	ph.supercluster_eta = pat_photon.superCluster()->eta();
972	<	ph.supercluster_phi = pat_photon.superCluster()->phi();
973	<	//       ph.neutralHadronIso = pat_photon.neutralHadronIso();
974	<	//       ph.chargedHadronIso = pat_photon.chargedHadronIso();
975	<	ph.trackIso = pat_photon.trackIso();
963	>	ph.set_charge(0);
964	>	ph.set_pt( pat_photon.pt());
965	>	ph.set_eta( pat_photon.eta());
966	>	ph.set_phi( pat_photon.phi());
967	>	ph.set_energy( pat_photon.energy());
968	>	ph.set_vertex_x(pat_photon.vertex().x());
969	>	ph.set_vertex_y(pat_photon.vertex().y());
970	>	ph.set_vertex_z(pat_photon.vertex().z());
971	>	ph.set_supercluster_eta(pat_photon.superCluster()->eta());
972	>	ph.set_supercluster_phi(pat_photon.superCluster()->phi());
973	>	//       ph.set_neutralHadronIso(pat_photon.neutralHadronIso());
974	>	//       ph.set_chargedHadronIso(pat_photon.chargedHadronIso());
975	>	ph.set_trackIso(pat_photon.trackIso());
976		phs[j].push_back(ph);
977		}
978		}
#	Line 528 | Line 992 | NtupleWriter::analyze(const edm::Event&
992		else{
993		pat::MET pat_met = pat_mets[0];
994
995	<	met[j].pt=pat_met.pt();
996	<	met[j].phi=pat_met.phi();
997	<	met[j].mEtSig= pat_met.mEtSig();
995	>	met[j].set_pt(pat_met.pt());
996	>	met[j].set_phi(pat_met.phi());
997	>	met[j].set_mEtSig(pat_met.mEtSig());
998		}
999
1000		}
1001		}
1002
1003		// ------------- trigger -------------
1004	<
1005	<	edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
1006	<	edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
1007	<	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();
1004	>	if(doTrigger){
1005	>	edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
1006	>	edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
1007	>	iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
1008
1009	<	const gen::PdfInfo* pdf = genEventInfo.pdf();
1010	<	if(pdf){
1011	<	genInfo.pdf_id1=pdf->id.first;
1012	<	genInfo.pdf_id2=pdf->id.second;
1013	<	genInfo.pdf_x1=pdf->x.first;
1014	<	genInfo.pdf_x2=pdf->x.second;
1015	<	genInfo.pdf_scalePDF=pdf->scalePDF;
1016	<	genInfo.pdf_xPDF1=pdf->xPDF.first;
1017	<	genInfo.pdf_xPDF2=pdf->xPDF.second;
1018	<	}
1019	<	else{
1020	<	genInfo.pdf_id1=-999;
1021	<	genInfo.pdf_id2=-999;
1022	<	genInfo.pdf_x1=-999;
1023	<	genInfo.pdf_x2=-999;
1024	<	genInfo.pdf_scalePDF=-999;
1025	<	genInfo.pdf_xPDF1=-999;
1026	<	genInfo.pdf_xPDF2=-999;
1027	<	}
1028	<
1029	<	edm::Handle<std::vector<PileupSummaryInfo> > pus;
1030	<	iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
1031	<	genInfo.pileup_NumInteractions_intime=0;
625	<	genInfo.pileup_NumInteractions_ootbefore=0;
626	<	genInfo.pileup_NumInteractions_ootafter=0;
627	<	if(pus.isValid()){
628	<	genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
629	<	for(size_t i=0; i<pus->size(); ++i){
630	<	if(pus->at(i).getBunchCrossing() == 0) // intime pileup
631	<	genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
632	<	else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
633	<	genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
634	<	}
635	<	else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
636	<	genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
637	<	}
1009	>	const edm::Provenance *meta = dummy_TriggerEvent.provenance();
1010	>	std::string nameProcess = meta->processName();
1011	>	edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
1012	>	triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
1013	>
1014	>	edm::Handle<edm::TriggerResults> trigger;
1015	>	iEvent.getByLabel(triggerResultTag, trigger);
1016	>	const edm::TriggerResults& trig = *(trigger.product());
1017	>
1018	>	triggerResults.clear();
1019	>	triggerNames.clear();
1020	>	//      L1_prescale.clear();
1021	>	//      HLT_prescale.clear();
1022	>
1023	>	edm::Service<edm::service::TriggerNamesService> tns;
1024	>	std::vector<std::string> triggerNames_all;
1025	>	tns->getTrigPaths(trig,triggerNames_all);
1026	>
1027	>	if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
1028	>	for(unsigned int i=0; i<trig.size(); ++i){
1029	>	std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
1030	>	for(; it!=trigger_prefixes.end(); ++it){
1031	>	if(triggerNames_all[i].substr(0, it->size()) == *it)break;
1032		}
1033	<	}
640	<
641	<	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++;
1033	>	if(it==trigger_prefixes.end()) continue;
1034
1035	<	//write out only top quarks and status 3 particles (works fine only for MadGraph)
1036	<	if(abs(iter->pdgId())==6 || iter->status()==3){
1037	<	GenParticle genp;
1038	<	genp.charge = iter->charge();;
1039	<	genp.pt = iter->p4().pt();
1040	<	genp.eta = iter->p4().eta();
1041	<	genp.phi = iter->p4().phi();
1042	<	genp.energy = iter->p4().E();
1043	<	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	<	}
1035	>	//triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
1036	>	triggerResults.push_back(trig.accept(i));
1037	>	if(newrun) triggerNames.push_back(triggerNames_all[i]);
1038	>	//        if(isRealData){
1039	>	//       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
1040	>	//       L1_prescale.push_back(pre.first);
1041	>	//       HLT_prescale.push_back(pre.second);
1042	>	//       //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
1043	>	//        }
1044		}
1045	<
1045	>	//    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
1046	>	//      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
1047	>	//    }
1048	>	newrun=false;
1049		}
1050
1051
1052		tr->Fill();
1053	+	if(doLumiInfo)
1054	+	previouslumiblockwasfilled=true;
1055	+
1056	+	// clean up
1057	+	m_stored_pfs.clear();
1058	+	if(doTopJetsConstituents){
1059	+	pfparticles.clear();
1060	+	}
1061
1062		}
1063
#	Line 685 | Line 1066 | NtupleWriter::analyze(const edm::Event&
1066		void
1067		NtupleWriter::beginJob()
1068		{
1069	+	if(doLumiInfo){
1070	+	totalRecLumi=0;
1071	+	totalDelLumi=0;
1072	+	previouslumiblockwasfilled=false;
1073	+	}
1074		}
1075
1076		// ------------ method called once each job just after ending the event loop  ------------
#	Line 700 | Line 1086 | NtupleWriter::endJob()
1086		void
1087		NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
1088		{
1089	<	bool setup_changed = false;
1090	<	hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1091	<	newrun=true;
1092	<
1093	<	edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
708	<	iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
709	<	JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
710	<	jecUnc = new JetCorrectionUncertainty(JetCorPar);
1089	>	if(doTrigger){
1090	>	//bool setup_changed = false;
1091	>	//hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1092	>	newrun=true;
1093	>	}
1094
1095		}
1096
#	Line 715 | Line 1098 | NtupleWriter::beginRun(edm::Run const& i
1098		void
1099		NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
1100		{
1101	+	if(doLumiInfo)
1102	+	std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
1103		}
1104
1105		// ------------ method called when starting to processes a luminosity block  ------------
1106		void
1107	<	NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
1107	>	NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
1108		{
1109	+	if(doLumiInfo){
1110	+	edm::Handle<LumiSummary> l;
1111	+	lumi.getByLabel("lumiProducer", l);
1112	+
1113	+	//add lumi of lumi blocks without any event to next lumiblock
1114	+	if(previouslumiblockwasfilled){
1115	+	intgRecLumi=0;
1116	+	intgDelLumi=0;
1117	+	}
1118	+	previouslumiblockwasfilled=false;
1119	+
1120	+	if (l.isValid()){;
1121	+	intgRecLumi+=l->intgRecLumi()*6.37;
1122	+	intgDelLumi+=l->intgDelLumi()*6.37;
1123	+	totalRecLumi+=l->intgRecLumi()*6.37;
1124	+	totalDelLumi+=l->intgDelLumi()*6.37;
1125	+	}
1126	+	//std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<"   " << l->intgDelLumi()*6.37<<std::endl;
1127	+	//std::cout << "summed:  "<< intgRecLumi << "   " << intgDelLumi << std::endl;
1128	+	}
1129		}
1130
1131		// ------------ method called when ending the processing of a luminosity block  ------------
#	Line 739 | Line 1144 | NtupleWriter::fillDescriptions(edm::Conf
1144		descriptions.addDefault(desc);
1145		}
1146
1147	+	// ------------ method fills constituents of the pat_jet into the Ntuple and stores a reference
1148	+	// ------------ to those in the topjet
1149	+	// ------------ it is checked if the constituents have been stored already
1150	+	void
1151	+	NtupleWriter::StoreJetConstituents(pat::Jet* pat_jet, Jet* jet)
1152	+	{
1153	+	// checks if the pf cand has already been stored, only stores so far missing
1154	+	// PF candidates, then stores a reference to the pf candidate in the jet
1155	+	// also calculates the jet charge and sets it
1156	+
1157	+
1158	+	const std::vector<reco::PFCandidatePtr> jconstits = pat_jet->getPFConstituents();
1159	+
1160	+	// loop over all jet constituents and store them
1161	+	for (unsigned int i=0; i<jconstits.size(); ++i){
1162	+
1163	+	PFParticle part;
1164	+	const reco::PFCandidate* pf = jconstits[i].get();
1165	+
1166	+	// check if it has already been stored, omit if true
1167	+	int is_already_in_list = -1;
1168	+	for (unsigned int j=0; j<m_stored_pfs.size(); ++j){
1169	+
1170	+	const reco::PFCandidate* spf = m_stored_pfs[j];
1171	+	double r2 = pow(pf->eta()-spf->eta(),2) + pow(pf->phi()-spf->phi(),2);
1172	+	double dpt = fabs( pf->pt() - spf->pt() );
1173	+	if (r2<1e-10 && dpt<1e-10){
1174	+	is_already_in_list = j;
1175	+	break;
1176	+	}
1177	+	}
1178	+
1179	+	if (is_already_in_list>-1){
1180	+	jet->add_pfconstituents_index(is_already_in_list);
1181	+	continue;
1182	+	}
1183	+
1184	+	part.set_pt(pf->pt());
1185	+	part.set_eta(pf->eta());
1186	+	part.set_phi(pf->phi());
1187	+	part.set_energy(pf->energy());
1188	+	part.set_charge(pf->charge());
1189	+
1190	+	part.set_ecal_en(pf->ecalEnergy());
1191	+	part.set_hcal_en(pf->hcalEnergy());
1192	+	reco::TrackRef trackref = pf->trackRef();
1193	+	if (!trackref.isNull()){
1194	+	part.set_track_mom(trackref->p());
1195	+	}
1196	+
1197	+	PFParticle::EParticleID id = PFParticle::eX;
1198	+	switch ( pf->translatePdgIdToType(pf->pdgId()) ){
1199	+	case reco::PFCandidate::X : id = PFParticle::eX; break;
1200	+	case reco::PFCandidate::h : id = PFParticle::eH; break;
1201	+	case reco::PFCandidate::e : id = PFParticle::eE; break;
1202	+	case reco::PFCandidate::mu : id = PFParticle::eMu; break;
1203	+	case reco::PFCandidate::gamma : id = PFParticle::eGamma; break;
1204	+	case reco::PFCandidate::h0 : id = PFParticle::eH0; break;
1205	+	case reco::PFCandidate::h_HF : id = PFParticle::eH_HF; break;
1206	+	case reco::PFCandidate::egamma_HF : id = PFParticle::eEgamma_HF; break;
1207	+	}
1208	+	part.set_particleID(id);
1209	+
1210	+	pfparticles.push_back(part);
1211	+	m_stored_pfs.push_back(pf);
1212	+
1213	+	// add a reference to the particle
1214	+	jet->add_pfconstituents_index(pfparticles.size()-1);
1215	+
1216	+	}
1217	+
1218	+	if(pat_jet->isPFJet()){
1219	+	jet->set_charge(pat_jet->charge());
1220	+	jet->set_neutralEmEnergyFraction (pat_jet->neutralEmEnergyFraction());
1221	+	jet->set_neutralHadronEnergyFraction (pat_jet->neutralHadronEnergyFraction());
1222	+	jet->set_chargedEmEnergyFraction (pat_jet->chargedEmEnergyFraction());
1223	+	jet->set_chargedHadronEnergyFraction (pat_jet->chargedHadronEnergyFraction());
1224	+	jet->set_muonEnergyFraction (pat_jet->muonEnergyFraction());
1225	+	jet->set_photonEnergyFraction (pat_jet->photonEnergyFraction());
1226	+	jet->set_chargedMultiplicity (pat_jet->chargedMultiplicity());
1227	+	jet->set_neutralMultiplicity (pat_jet->neutralMultiplicity());
1228	+	jet->set_muonMultiplicity (pat_jet->muonMultiplicity());
1229	+	jet->set_electronMultiplicity (pat_jet->electronMultiplicity());
1230	+	jet->set_photonMultiplicity (pat_jet->photonMultiplicity());
1231	+	}
1232	+
1233	+	return;
1234	+
1235	+	}
1236	+
1237		//define this as a plug-in
1238		DEFINE_FWK_MODULE(NtupleWriter);

Diff Legend

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