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.1.1.1 by peiffer, Mon Apr 2 15:18:11 2012 UTC vs.
Revision 1.30 by jott, Thu Jun 20 15:03:51 2013 UTC

#	Line 17 | Line 17
17		//
18		//
19
20	<	#include "UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
21	<
22	<
23	<
24	<	//
25	<	// constants, enums and typedefs
26	<	//
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	>	namespace{
30	>
31	>	/// Use info from pat_jet to fill info in jet
32	>	void fill_jet_info(const pat::Jet & pat_jet, Jet & jet)
33	>	{
34	>	jet.set_charge(pat_jet.charge());
35	>	jet.set_pt(pat_jet.pt());
36	>	jet.set_eta(pat_jet.eta());
37	>	jet.set_phi(pat_jet.phi());
38	>	jet.set_energy(pat_jet.energy());
39	>	jet.set_flavor(pat_jet.partonFlavour());
40	>	jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
41	>	const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
42	>	jet.set_nTracks ( jettracks.size());
43	>	jet.set_jetArea(pat_jet.jetArea());
44	>	if(pat_jet.isPFJet()){
45	>	jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
46	>	jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
47	>	jet.set_chargedEmEnergyFraction (pat_jet.chargedEmEnergyFraction());
48	>	jet.set_chargedHadronEnergyFraction (pat_jet.chargedHadronEnergyFraction());
49	>	jet.set_muonEnergyFraction (pat_jet.muonEnergyFraction());
50	>	jet.set_photonEnergyFraction (pat_jet.photonEnergyFraction());
51	>	jet.set_chargedMultiplicity (pat_jet.chargedMultiplicity());
52	>	jet.set_neutralMultiplicity (pat_jet.neutralMultiplicity());
53	>	jet.set_muonMultiplicity (pat_jet.muonMultiplicity());
54	>	jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
55	>	jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
56	>	}
57	>	jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
58	>	jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
59	>	jet.set_btag_simpleSecondaryVertexHighPur(pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
60	>	jet.set_btag_combinedSecondaryVertex(pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags"));
61	>	jet.set_btag_combinedSecondaryVertexMVA(pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags"));
62	>	jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
63	>	jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
64	>	}
65
66	<	//
29	<	// static data member definitions
30	<	//
66	>	}
67
68		//
69		// constructors and destructor
70		//
71		NtupleWriter::NtupleWriter(const edm::ParameterSet& iConfig)
36	–
72		{
73	<	//now do what ever initialization is needed
39	<	//edm::Service<TFileService> fs;
40	<	//tr = fs->make<TTree>("AnalysisTree","AnalysisTree");
41	<
42	<	fileName =  iConfig.getParameter<std::string>("fileName");
73	>	fileName = iConfig.getParameter<std::string>("fileName");
74		outfile = new TFile(fileName,"RECREATE");
75		outfile->cd();
76		tr = new TTree("AnalysisTree","AnalysisTree");
77
47	–	std::string name;
48	–
78		doElectrons = iConfig.getParameter<bool>("doElectrons");
79		doMuons = iConfig.getParameter<bool>("doMuons");
80		doTaus = iConfig.getParameter<bool>("doTaus");
81		doJets = iConfig.getParameter<bool>("doJets");
82	+	doGenJets = iConfig.getParameter<bool>("doGenJets");
83	+	doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");
84		doPhotons = iConfig.getParameter<bool>("doPhotons");
85		doMET = iConfig.getParameter<bool>("doMET");
86		doGenInfo = iConfig.getParameter<bool>("doGenInfo");
87	+	doAllGenParticles = iConfig.getParameter<bool>("doAllGenParticles");
88	+	doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
89		doPV = iConfig.getParameter<bool>("doPV");
90		doTopJets = iConfig.getParameter<bool>("doTopJets");
91	+	doTopJetsConstituents = iConfig.getParameter<bool>("doTopJetsConstituents");
92	+	doTrigger = iConfig.getParameter<bool>("doTrigger");
93	+	SVComputer_  = iConfig.getUntrackedParameter<edm::InputTag>("svComputer",edm::InputTag("combinedSecondaryVertex"));
94	+	doTagInfos = iConfig.getUntrackedParameter<bool>("doTagInfos",false);
95	+	storePFsAroundLeptons = iConfig.getUntrackedParameter<bool>("storePFsAroundLeptons",false);
96
97		// initialization of tree variables
98
#	Line 62 | Line 100 | NtupleWriter::NtupleWriter(const edm::Pa
100		tr->Branch("event",&event);
101		tr->Branch("luminosityBlock",&luminosityBlock);
102		tr->Branch("isRealData",&isRealData);
103	<	tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
104	<	tr->Branch("beamspot_x0",&beamspot_x0);
67	<	tr->Branch("beamspot_y0",&beamspot_y0);
68	<	tr->Branch("beamspot_z0",&beamspot_z0);
103	>	tr->Branch("rho",&rho);
104	>	rho_source = iConfig.getParameter<edm::InputTag>("rho_source");
105
106	+	if(doLumiInfo){
107	+	tr->Branch("intgRecLumi",&intgRecLumi);
108	+	tr->Branch("intgDelLumi",&intgDelLumi);
109	+	}
110	+	if(doPV){
111	+	tr->Branch("beamspot_x0",&beamspot_x0);
112	+	tr->Branch("beamspot_y0",&beamspot_y0);
113	+	tr->Branch("beamspot_z0",&beamspot_z0);
114	+	}
115		if(doElectrons){
116		electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
117		for(size_t j=0; j< electron_sources.size(); ++j){
#	Line 95 | Line 140 | NtupleWriter::NtupleWriter(const edm::Pa
140		tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
141		}
142		}
143	+	if(doGenJets){
144	+	genjet_sources = iConfig.getParameter<std::vector<std::string> >("genjet_sources");
145	+	genjet_ptmin = iConfig.getParameter<double> ("genjet_ptmin");
146	+	genjet_etamax = iConfig.getParameter<double> ("genjet_etamax");
147	+	for(size_t j=0; j< genjet_sources.size(); ++j){
148	+	tr->Branch( genjet_sources[j].c_str(), "std::vector<Particle>", &genjets[j]);
149	+	}
150	+	}
151		if(doTopJets){
152		topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
153		topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
#	Line 103 | Line 156 | NtupleWriter::NtupleWriter(const edm::Pa
156		tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
157		}
158		}
159	+	if(doTopJetsConstituents){
160	+	topjet_constituents_sources = iConfig.getParameter<std::vector<std::string> >("topjet_constituents_sources");
161	+	tr->Branch( "PFParticles", "std::vector<PFParticle>", &pfparticles);
162	+	}
163	+	if(doGenTopJets){
164	+	gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
165	+	gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
166	+	gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
167	+	for(size_t j=0; j< gentopjet_sources.size(); ++j){
168	+	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<GenTopJet>", &gentopjets[j]);
169	+	}
170	+	}
171		if(doPhotons){
172		photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
173		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 122 | Line 187 | NtupleWriter::NtupleWriter(const edm::Pa
187		}
188		}
189		if(doGenInfo){
190	+	genparticle_source= iConfig.getParameter<edm::InputTag>("genparticle_source");
191		tr->Branch("genInfo","GenInfo",&genInfo);
192		tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
193		}
194	<
195	<	trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
196	<	//tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
197	<	tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
198	<	tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
199	<	tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
200	<	tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
201	<
194	>	if(doTrigger){
195	>	trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
196	>	//tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
197	>	tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
198	>	tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
199	>	//tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
200	>	//tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
201	>	}
202	>	if (storePFsAroundLeptons){
203	>	pf_around_leptons_source = iConfig.getParameter<std::string>("pf_around_leptons_source");
204	>	}
205		newrun = true;
206		}
207
#	Line 151 | Line 220 | NtupleWriter::~NtupleWriter()
220		//
221
222		// ------------ method called for each event  ------------
223	<	void
155	<	NtupleWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
223	>	void NtupleWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
224		{
157	–	//   using namespace edm;
158	–
159	–
225		// ------------- common variables -------------
226
227		run = iEvent.id().run();
#	Line 164 | Line 229 | NtupleWriter::analyze(const edm::Event&
229		luminosityBlock = iEvent.luminosityBlock();
230		isRealData      = iEvent.isRealData();
231
232	<	if(isRealData){
233	<	edm::Handle<bool> bool_handle;
234	<	iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
235	<	HBHENoiseFilterResult = *(bool_handle.product());
232	>	edm::Handle<double> m_rho;
233	>	iEvent.getByLabel(rho_source,m_rho);
234	>	rho=*m_rho;
235	>
236	>	// ------------- primary vertices and beamspot  -------------
237	>
238	>	if(doPV){
239	>	for(size_t j=0; j< pv_sources.size(); ++j){
240	>	pvs[j].clear();
241	>
242	>	edm::Handle< std::vector<reco::Vertex> > pv_handle;
243	>	iEvent.getByLabel(pv_sources[j], pv_handle);
244	>	const std::vector<reco::Vertex>& reco_pvs = *(pv_handle.product());
245	>
246	>	for (unsigned int i = 0; i <  reco_pvs.size(); ++i) {
247	>	reco::Vertex reco_pv = reco_pvs[i];
248	>
249	>	PrimaryVertex pv;
250	>	pv.set_x( reco_pv.x());
251	>	pv.set_y( reco_pv.y());
252	>	pv.set_z( reco_pv.z());
253	>	pv.set_nTracks( reco_pv.nTracks());
254	>	//pv.set_isValid( reco_pv.isValid());
255	>	pv.set_chi2( reco_pv.chi2());
256	>	pv.set_ndof( reco_pv.ndof());
257	>
258	>	pvs[j].push_back(pv);
259	>	}
260	>	}
261	>
262	>	edm::Handle<reco::BeamSpot> beamSpot;
263	>	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
264	>	const reco::BeamSpot & bsp = *beamSpot;
265	>
266	>	beamspot_x0 = bsp.x0();
267	>	beamspot_y0 = bsp.y0();
268	>	beamspot_z0 = bsp.z0();
269		}
270	<	else HBHENoiseFilterResult = false;
270	>
271	>	// ------------- generator info -------------
272	>
273	>	if(doGenInfo){
274	>	genInfo.clear_weights();
275	>	genInfo.clear_binningValues();
276	>	genps.clear();
277	>
278	>	edm::Handle<GenEventInfoProduct> genEventInfoProduct;
279	>	iEvent.getByLabel("generator", genEventInfoProduct);
280	>	const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
281	>
282	>	for(unsigned int k=0; k<genEventInfo.binningValues().size();++k){
283	>	genInfo.add_binningValue(genEventInfo.binningValues().at(k));
284	>	}
285	>	for(unsigned int k=0; k<genEventInfo.weights().size();++k){
286	>	genInfo.add_weight(genEventInfo.weights().at(k));
287	>	}
288	>	genInfo.set_alphaQCD(genEventInfo.alphaQCD());
289	>	genInfo.set_alphaQED(genEventInfo.alphaQED());
290	>	genInfo.set_qScale(genEventInfo.qScale());
291	>
292	>	const gen::PdfInfo* pdf = genEventInfo.pdf();
293	>	if(pdf){
294	>	genInfo.set_pdf_id1(pdf->id.first);
295	>	genInfo.set_pdf_id2(pdf->id.second);
296	>	genInfo.set_pdf_x1(pdf->x.first);
297	>	genInfo.set_pdf_x2(pdf->x.second);
298	>	genInfo.set_pdf_scalePDF(pdf->scalePDF);
299	>	genInfo.set_pdf_xPDF1(pdf->xPDF.first);
300	>	genInfo.set_pdf_xPDF2(pdf->xPDF.second);
301	>	}
302	>	else{
303	>	genInfo.set_pdf_id1(-999);
304	>	genInfo.set_pdf_id2(-999);
305	>	genInfo.set_pdf_x1(-999);
306	>	genInfo.set_pdf_x2(-999);
307	>	genInfo.set_pdf_scalePDF(-999);
308	>	genInfo.set_pdf_xPDF1(-999);
309	>	genInfo.set_pdf_xPDF2(-999);
310	>	}
311	>
312	>	edm::Handle<std::vector<PileupSummaryInfo> > pus;
313	>	iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
314	>	genInfo.set_pileup_NumInteractions_intime(0);
315	>	genInfo.set_pileup_NumInteractions_ootbefore(0);
316	>	genInfo.set_pileup_NumInteractions_ootafter(0);
317	>	if(pus.isValid()){
318	>	genInfo.set_pileup_TrueNumInteractions ( (float) pus->at(0).getTrueNumInteractions());
319	>	for(size_t i=0; i<pus->size(); ++i){
320	>	if(pus->at(i).getBunchCrossing() == 0) // intime pileup
321	>	genInfo.set_pileup_NumInteractions_intime( genInfo.pileup_NumInteractions_intime() + pus->at(i).getPU_NumInteractions());
322	>	else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
323	>	genInfo.set_pileup_NumInteractions_ootbefore( genInfo.pileup_NumInteractions_ootbefore() + pus->at(i).getPU_NumInteractions());
324	>	}
325	>	else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
326	>	genInfo.set_pileup_NumInteractions_ootafter( genInfo.pileup_NumInteractions_ootafter() + pus->at(i).getPU_NumInteractions());
327	>	}
328	>	}
329	>	}
330	>
331	>	edm::Handle<reco::GenParticleCollection> genPartColl;
332	>	iEvent.getByLabel(genparticle_source, genPartColl);
333	>	int index=-1;
334	>	for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
335	>	index++;
336	>
337	>	//write out only top quarks,final state leptons and status 3 particles (works fine only for MadGraph)
338	>	bool islepton = iter->status()==1 && abs(iter->pdgId())>=11 && abs(iter->pdgId())<=16 ;
339	>	if(abs(iter->pdgId())==6 || iter->status()==3 || islepton ||  doAllGenParticles){
340	>	GenParticle genp;
341	>	genp.set_charge(iter->charge());
342	>	genp.set_pt(iter->p4().pt());
343	>	genp.set_eta(iter->p4().eta());
344	>	genp.set_phi(iter->p4().phi());
345	>	genp.set_energy(iter->p4().E());
346	>	genp.set_index(index);
347	>	genp.set_status( iter->status());
348	>	genp.set_pdgId( iter->pdgId());
349	>
350	>	genp.set_mother1(-1);
351	>	genp.set_mother2(-1);
352	>	genp.set_daughter1(-1);
353	>	genp.set_daughter2(-1);
354	>
355	>	int nm=iter->numberOfMothers();
356	>	int nd=iter->numberOfDaughters();
357	>
358	>
359	>	if (nm>0) genp.set_mother1( iter->motherRef(0).key());
360	>	if (nm>1) genp.set_mother2( iter->motherRef(1).key());
361	>	if (nd>0) genp.set_daughter1( iter->daughterRef(0).key());
362	>	if (nd>1) genp.set_daughter2( iter->daughterRef(1).key());
363	>
364	>	genps.push_back(genp);
365	>	}
366	>	}
367	>	}
368	>
369	>	// ------------- full PF collection -------------
370	>
371
372		// ------------- electrons -------------
373		if(doElectrons){
374	+
375	+	//      edm::Handle<reco::ConversionCollection> hConversions;
376	+	//      iEvent.getByLabel("allConversions", hConversions);
377	+
378	+	//      edm::Handle<reco::BeamSpot> beamSpot;
379	+	//      iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
380	+	//      const reco::BeamSpot & bsp = *beamSpot;
381	+
382		for(size_t j=0; j< electron_sources.size(); ++j){
383		eles[j].clear();
384		edm::Handle< std::vector<pat::Electron> > ele_handle;
#	Line 183 | Line 389 | NtupleWriter::analyze(const edm::Event&
389		pat::Electron pat_ele = pat_electrons[i];
390		Electron ele;
391
392	<	ele.charge =  pat_ele.charge();
393	<	ele.pt =  pat_ele.pt();
394	<	ele.eta =  pat_ele.eta();
395	<	ele.phi =  pat_ele.phi();
396	<	ele.energy =  pat_ele.energy();
397	<	ele.vertex_x = pat_ele.vertex().x();
398	<	ele.vertex_y = pat_ele.vertex().y();
399	<	ele.vertex_z = pat_ele.vertex().z();
400	<	ele.supercluster_eta = pat_ele.superCluster()->eta();
401	<	ele.supercluster_phi = pat_ele.superCluster()->phi();
402	<	ele.dB = pat_ele.dB();
403	<	//ele.particleIso = pat_ele.particleIso();
404	<	ele.neutralHadronIso = pat_ele.neutralHadronIso();
405	<	ele.chargedHadronIso = pat_ele.chargedHadronIso();
406	<	ele.trackIso = pat_ele.trackIso();
407	<	ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
392	>	ele.set_charge( pat_ele.charge());
393	>	ele.set_pt( pat_ele.pt());
394	>	ele.set_eta( pat_ele.eta());
395	>	ele.set_phi( pat_ele.phi());
396	>	ele.set_energy( pat_ele.energy());
397	>	ele.set_vertex_x(pat_ele.vertex().x());
398	>	ele.set_vertex_y(pat_ele.vertex().y());
399	>	ele.set_vertex_z(pat_ele.vertex().z());
400	>	ele.set_supercluster_eta(pat_ele.superCluster()->eta());
401	>	ele.set_supercluster_phi(pat_ele.superCluster()->phi());
402	>	ele.set_dB(pat_ele.dB());
403	>	//ele.set_particleIso(pat_ele.particleIso());
404	>	ele.set_neutralHadronIso(pat_ele.neutralHadronIso());
405	>	ele.set_chargedHadronIso(pat_ele.chargedHadronIso());
406	>	ele.set_trackIso(pat_ele.trackIso());
407	>	ele.set_photonIso(pat_ele.photonIso());
408	>	ele.set_puChargedHadronIso(pat_ele.puChargedHadronIso());
409	>	ele.set_gsfTrack_trackerExpectedHitsInner_numberOfLostHits(pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits());
410	>	ele.set_gsfTrack_px( pat_ele.gsfTrack()->px());
411	>	ele.set_gsfTrack_py( pat_ele.gsfTrack()->py());
412	>	ele.set_gsfTrack_pz( pat_ele.gsfTrack()->pz());
413	>	ele.set_gsfTrack_vx( pat_ele.gsfTrack()->vx());
414	>	ele.set_gsfTrack_vy( pat_ele.gsfTrack()->vy());
415	>	ele.set_gsfTrack_vz( pat_ele.gsfTrack()->vz());
416	>	//ele.set_passconversionveto(!ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position()));
417	>	ele.set_passconversionveto(pat_ele.passConversionVeto());
418	>	ele.set_dEtaIn(pat_ele.deltaEtaSuperClusterTrackAtVtx());
419	>	ele.set_dPhiIn(pat_ele.deltaPhiSuperClusterTrackAtVtx());
420	>	ele.set_sigmaIEtaIEta(pat_ele.sigmaIetaIeta());
421	>	ele.set_HoverE(pat_ele.hadronicOverEm());
422	>	ele.set_fbrem(pat_ele.fbrem());
423	>	ele.set_EoverPIn(pat_ele.eSuperClusterOverP());
424	>	ele.set_EcalEnergy(pat_ele.ecalEnergy());
425	>	ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
426	>	ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
427	>	float AEff03 = 0.00;
428	>	if(isRealData){
429	>	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2011);
430	>	}else{
431	>	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAFall11MC);
432	>	}
433	>	ele.set_AEff(AEff03);
434
435		eles[j].push_back(ele);
436	+
437	+	if (storePFsAroundLeptons){
438	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
439	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
440	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
441	+	StorePFCandsInCone(&ele, pf_cands, 0.4);
442	+	}
443	+
444		}
445		}
446		}
#	Line 218 | Line 458 | NtupleWriter::analyze(const edm::Event&
458		pat::Muon pat_mu = pat_muons[i];
459
460		Muon mu;
461	<	mu.charge =  pat_mu.charge();
462	<	mu.pt =  pat_mu.pt();
463	<	mu.eta =  pat_mu.eta();
464	<	mu.phi =  pat_mu.phi();
465	<	mu.energy =  pat_mu.energy();
466	<	mu.vertex_x = pat_mu.vertex().x();
467	<	mu.vertex_y = pat_mu.vertex().y();
468	<	mu.vertex_z = pat_mu.vertex().z();
469	<	mu.dB = pat_mu.dB();
470	<	//mu.particleIso = pat_mu.particleIso();
471	<	mu.neutralHadronIso = pat_mu.neutralHadronIso();
472	<	mu.chargedHadronIso = pat_mu.chargedHadronIso();
473	<	mu.trackIso = pat_mu.trackIso();
474	<	mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
475	<	mu.isGlobalMuon = pat_mu.isGlobalMuon();
476	<	mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
477	<	mu.isTrackerMuon = pat_mu.isTrackerMuon();
478	<	mu.numberOfMatchedStations = pat_mu.numberOfMatchedStations();
461	>	mu.set_charge( pat_mu.charge());
462	>	mu.set_pt( pat_mu.pt());
463	>	mu.set_eta( pat_mu.eta());
464	>	mu.set_phi( pat_mu.phi());
465	>	mu.set_energy( pat_mu.energy());
466	>	mu.set_vertex_x ( pat_mu.vertex().x());
467	>	mu.set_vertex_y ( pat_mu.vertex().y());
468	>	mu.set_vertex_z ( pat_mu.vertex().z());
469	>	mu.set_dB ( pat_mu.dB());
470	>	//mu.particleIso ( pat_mu.particleIso());
471	>	mu.set_neutralHadronIso ( pat_mu.neutralHadronIso());
472	>	mu.set_chargedHadronIso ( pat_mu.chargedHadronIso());
473	>	mu.set_trackIso ( pat_mu.trackIso());
474	>	mu.set_photonIso ( pat_mu.photonIso());
475	>	mu.set_puChargedHadronIso ( pat_mu.puChargedHadronIso());
476	>	mu.set_isGlobalMuon ( pat_mu.isGlobalMuon());
477	>	mu.set_isPFMuon ( pat_mu.isPFMuon());
478	>	mu.set_isStandAloneMuon ( pat_mu.isStandAloneMuon());
479	>	mu.set_isTrackerMuon ( pat_mu.isTrackerMuon());
480	>	mu.set_numberOfMatchedStations ( pat_mu.numberOfMatchedStations());
481		reco::TrackRef globalTrack = pat_mu.globalTrack();
482		if(!globalTrack.isNull()){
483	<	mu.globalTrack_chi2 = globalTrack->chi2();
484	<	mu.globalTrack_ndof = globalTrack->ndof();
485	<	mu.globalTrack_d0 = globalTrack->d0();
486	<	mu.globalTrack_d0Error = globalTrack->d0Error();
487	<	mu.globalTrack_numberOfValidHits = globalTrack->numberOfValidHits();
488	<	mu.globalTrack_numberOfLostHits = globalTrack->numberOfLostHits();
483	>	mu.set_globalTrack_chi2 ( globalTrack->chi2());
484	>	mu.set_globalTrack_ndof ( globalTrack->ndof());
485	>	mu.set_globalTrack_d0 ( globalTrack->d0());
486	>	mu.set_globalTrack_d0Error ( globalTrack->d0Error());
487	>	mu.set_globalTrack_numberOfValidHits ( globalTrack->numberOfValidHits());
488	>	mu.set_globalTrack_numberOfLostHits ( globalTrack->numberOfLostHits());
489	>	mu.set_globalTrack_numberOfValidMuonHits(globalTrack->hitPattern().numberOfValidMuonHits() );
490		}
491		else{
492	<	mu.globalTrack_chi2 = 0;
493	<	mu.globalTrack_ndof = 0;
494	<	mu.globalTrack_d0 = 0;
495	<	mu.globalTrack_d0Error = 0;
496	<	mu.globalTrack_numberOfValidHits = 0;
497	<	mu.globalTrack_numberOfLostHits = 0;
492	>	mu.set_globalTrack_chi2 ( 0);
493	>	mu.set_globalTrack_ndof ( 0);
494	>	mu.set_globalTrack_d0 ( 0);
495	>	mu.set_globalTrack_d0Error ( 0);
496	>	mu.set_globalTrack_numberOfValidHits ( 0);
497	>	mu.set_globalTrack_numberOfLostHits ( 0);
498		}
499		reco::TrackRef innerTrack = pat_mu.innerTrack();
500		if(!innerTrack.isNull()){
501	<	mu.innerTrack_chi2 = innerTrack->chi2();
502	<	mu.innerTrack_ndof = innerTrack->ndof();
503	<	mu.innerTrack_d0 = innerTrack->d0();
504	<	mu.innerTrack_d0Error = innerTrack->d0Error();
505	<	mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
506	<	mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
501	>	mu.set_innerTrack_chi2 ( innerTrack->chi2());
502	>	mu.set_innerTrack_ndof ( innerTrack->ndof());
503	>	mu.set_innerTrack_d0 ( innerTrack->d0());
504	>	mu.set_innerTrack_d0Error ( innerTrack->d0Error());
505	>	mu.set_innerTrack_numberOfValidHits ( innerTrack->numberOfValidHits());
506	>	mu.set_innerTrack_numberOfLostHits ( innerTrack->numberOfLostHits());
507	>	mu.set_innerTrack_trackerLayersWithMeasurement ( innerTrack->hitPattern().trackerLayersWithMeasurement());
508	>	mu.set_innerTrack_numberOfValidPixelHits ( innerTrack->hitPattern().numberOfValidPixelHits());
509		}
510		else{
511	<	mu.innerTrack_chi2 = 0;
512	<	mu.innerTrack_ndof = 0;
513	<	mu.innerTrack_d0 = 0;
514	<	mu.innerTrack_d0Error = 0;
515	<	mu.innerTrack_numberOfValidHits = 0;
516	<	mu.innerTrack_numberOfLostHits = 0;
511	>	mu.set_innerTrack_chi2 ( 0);
512	>	mu.set_innerTrack_ndof ( 0);
513	>	mu.set_innerTrack_d0 ( 0);
514	>	mu.set_innerTrack_d0Error ( 0);
515	>	mu.set_innerTrack_numberOfValidHits ( 0);
516	>	mu.set_innerTrack_numberOfLostHits ( 0);
517	>	mu.set_innerTrack_trackerLayersWithMeasurement ( 0);
518	>	mu.set_innerTrack_numberOfValidPixelHits ( 0);
519		}
520		reco::TrackRef outerTrack = pat_mu.outerTrack();
521		if(!outerTrack.isNull()){
522	<	mu.outerTrack_chi2 = outerTrack->chi2();
523	<	mu.outerTrack_ndof = outerTrack->ndof();
524	<	mu.outerTrack_d0 = outerTrack->d0();
525	<	mu.outerTrack_d0Error = outerTrack->d0Error();
526	<	mu.outerTrack_numberOfValidHits = outerTrack->numberOfValidHits();
527	<	mu.outerTrack_numberOfLostHits = outerTrack->numberOfLostHits();
522	>	mu.set_outerTrack_chi2 ( outerTrack->chi2());
523	>	mu.set_outerTrack_ndof ( outerTrack->ndof());
524	>	mu.set_outerTrack_d0 ( outerTrack->d0());
525	>	mu.set_outerTrack_d0Error ( outerTrack->d0Error());
526	>	mu.set_outerTrack_numberOfValidHits ( outerTrack->numberOfValidHits());
527	>	mu.set_outerTrack_numberOfLostHits ( outerTrack->numberOfLostHits());
528		}
529		else{
530	<	mu.outerTrack_chi2 = 0;
531	<	mu.outerTrack_ndof = 0;
532	<	mu.outerTrack_d0 = 0;
533	<	mu.outerTrack_d0Error = 0;
534	<	mu.outerTrack_numberOfValidHits = 0;
535	<	mu.outerTrack_numberOfLostHits = 0;
530	>	mu.set_outerTrack_chi2 ( 0);
531	>	mu.set_outerTrack_ndof ( 0);
532	>	mu.set_outerTrack_d0 ( 0);
533	>	mu.set_outerTrack_d0Error ( 0);
534	>	mu.set_outerTrack_numberOfValidHits ( 0);
535	>	mu.set_outerTrack_numberOfLostHits ( 0);
536		}
537
538		mus[j].push_back(mu);
539	+
540	+	if (storePFsAroundLeptons){
541	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
542	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
543	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
544	+	StorePFCandsInCone(&mu, pf_cands, 0.4);
545	+	}
546	+
547		}
548		}
549		}
#	Line 309 | Line 564 | NtupleWriter::analyze(const edm::Event&
564		if(fabs(pat_tau.eta()) > tau_etamax) continue;
565
566		Tau tau;
567	<	tau.charge =  pat_tau.charge();
568	<	tau.pt =  pat_tau.pt();
569	<	tau.eta =  pat_tau.eta();
570	<	tau.phi =  pat_tau.phi();
571	<	tau.energy =  pat_tau.energy();
567	>	tau.set_charge( pat_tau.charge());
568	>	tau.set_pt( pat_tau.pt());
569	>	tau.set_eta( pat_tau.eta());
570	>	tau.set_phi( pat_tau.phi());
571	>	tau.set_energy( pat_tau.energy());
572	>	tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
573	>	//tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
574	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
575	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
576	>	tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
577	>	tau.set_byLooseIsolationMVA( pat_tau.tauID("byLooseIsolationMVA")>0.5);
578	>	tau.set_byMediumIsolationMVA( pat_tau.tauID("byMediumIsolationMVA")>0.5);
579	>	tau.set_byTightIsolationMVA( pat_tau.tauID("byTightIsolationMVA")>0.5);
580	>	tau.set_byLooseIsolationMVA2( pat_tau.tauID("byLooseIsolationMVA2")>0.5);
581	>	tau.set_byMediumIsolationMVA2( pat_tau.tauID("byMediumIsolationMVA2")>0.5);
582	>	tau.set_byTightIsolationMVA2( pat_tau.tauID("byTightIsolationMVA2")>0.5);
583	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr3Hits(  pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits")>0.5);
584	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits")>0.5);
585	>	tau.set_byTightCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr3Hits")>0.5);
586	>	tau.set_againstElectronLooseMVA3  ( pat_tau.tauID("againstElectronLooseMVA3")>0.5);
587	>	tau.set_againstElectronMediumMVA3 ( pat_tau.tauID("againstElectronMediumMVA3")>0.5);
588	>	tau.set_againstElectronTightMVA3 ( pat_tau.tauID("againstElectronTightMVA3")>0.5);
589	>	tau.set_againstElectronVTightMVA3 ( pat_tau.tauID("againstElectronVTightMVA3")>0.5);
590	>	tau.set_againstMuonLoose2 ( pat_tau.tauID("againstMuonLoose2")>0.5);
591	>	tau.set_againstMuonMedium2 ( pat_tau.tauID("againstMuonMedium2")>0.5);
592	>	tau.set_againstMuonTight2 ( pat_tau.tauID("againstMuonTight2")>0.5);
593	>	tau.set_byIsolationMVAraw(  pat_tau.tauID("byIsolationMVAraw"));
594	>	tau.set_byIsolationMVA2raw(  pat_tau.tauID("byIsolationMVA2raw"));
595	>	tau.set_decayMode( pat_tau.decayMode() );
596	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw"));
597	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw3Hits( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));
598
599	+	//       std::cout << pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") << std::endl;
600	+
601	+	//       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
602	+	//       if(!leadPFCand.isNull()){
603	+	//         tau.set_leadPFCand_px ( leadPFCand->px());
604	+	//         tau.set_leadPFCand_py ( leadPFCand->py());
605	+	//         tau.set_leadPFCand_pz ( leadPFCand->pz());
606	+	//       }
607	+	//       else{
608	+	//         tau.set_leadPFCand_px ( 0);
609	+	//         tau.set_leadPFCand_py ( 0);
610	+	//         tau.set_leadPFCand_pz ( 0);
611	+	//       }
612		taus[j].push_back(tau);
613	+
614	+	// store isolation info only for identified taus
615	+	if (pat_tau.tauID("decayModeFinding")>0.5){
616	+	bool storepfs = (pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5) || (pat_tau.tauID("byLooseIsolationMVA")>0.5);
617	+	if (storePFsAroundLeptons && storepfs){
618	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
619	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
620	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
621	+	StorePFCandsInCone(&tau, pf_cands, 0.4);
622	+	}
623	+	}
624	+	}
625	+	}
626	+	}
627	+
628	+	//-------------- gen jets -------------
629	+
630	+	if(doGenJets){
631	+	for(size_t j=0; j< genjet_sources.size(); ++j){
632	+
633	+	genjets[j].clear();
634	+
635	+	edm::Handle< std::vector<reco::GenJet> > genjet_handle;
636	+	iEvent.getByLabel(genjet_sources[j], genjet_handle);
637	+	const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
638	+
639	+	for (unsigned int i = 0; i < gen_jets.size(); ++i) {
640	+
641	+	const reco::GenJet* gen_jet = &gen_jets[i];
642	+	if(gen_jet->pt() < genjet_ptmin) continue;
643	+	if(fabs(gen_jet->eta()) > genjet_etamax) continue;
644	+
645	+	Particle jet;
646	+	jet.set_charge(gen_jet->charge());
647	+	jet.set_pt(gen_jet->pt());
648	+	jet.set_eta(gen_jet->eta());
649	+	jet.set_phi(gen_jet->phi());
650	+	jet.set_energy(gen_jet->energy());
651	+
652	+	// recalculate the jet charge
653	+	int jet_charge = 0;
654	+	std::vector<const reco::GenParticle * > jetgenps = gen_jet->getGenConstituents();
655	+	for(unsigned int l = 0; l<jetgenps.size(); ++l){
656	+	jet_charge +=  jetgenps[l]->charge();
657	+	}
658	+	jet.set_charge(jet_charge);
659	+
660	+	genjets[j].push_back(jet);
661		}
662		}
663		}
#	Line 339 | Line 681 | NtupleWriter::analyze(const edm::Event&
681		//       for(size_t k=0; k<dis.size(); ++k){
682		//         std::cout << dis[k].first << std::endl;
683		//       }
342	–
684		Jet jet;
685	<	jet.charge = pat_jet.charge();
686	<	jet.pt = pat_jet.pt();
687	<	jet.eta = pat_jet.eta();
688	<	jet.phi = pat_jet.phi();
689	<	jet.energy = pat_jet.energy();
690	<	jet.numberOfDaughters =pat_jet.numberOfDaughters();
691	<	const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
692	<	jet.nTracks = jettracks.size();
693	<	jet.jetArea = pat_jet.jetArea();
694	<	jet.pileup = pat_jet.pileup();
695	<	if(pat_jet.isPFJet()){
696	<	jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
356	<	jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
357	<	jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
358	<	jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
359	<	jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
360	<	jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
361	<	jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
362	<	jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
363	<	jet.muonMultiplicity =pat_jet.muonMultiplicity();
364	<	jet.electronMultiplicity =pat_jet.electronMultiplicity();
365	<	jet.photonMultiplicity =pat_jet.photonMultiplicity();
366	<	}
367	<	jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
368	<	jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
369	<	jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
370	<	jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
371	<	jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
372	<	jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
373	<
685	>	fill_jet_info(pat_jet, jet);
686	>	const reco::GenJet *genj = pat_jet.genJet();
687	>	if(genj){
688	>	for(unsigned int k=0; k<genjets->size(); ++k){
689	>	if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta()){
690	>	jet.set_genjet_index(k);
691	>	}
692	>	}
693	>	//         if( jet.genjet_index()<0){
694	>	//           std::cout<< "genjet not found for " << genj->pt() << "  " << genj->eta() << std::endl;
695	>	//         }
696	>	}
697		jets[j].push_back(jet);
698		}
699		}
700		}
701
702	+
703		// ------------- top jets -------------
704		if(doTopJets){
705	+
706		for(size_t j=0; j< topjet_sources.size(); ++j){
707
708		topjets[j].clear();
#	Line 387 | Line 712 | NtupleWriter::analyze(const edm::Event&
712		iEvent.getByLabel(topjet_sources[j], pat_topjets);
713
714		for (unsigned int i = 0; i < pat_topjets->size(); i++) {
715	<	const pat::Jet  pat_topjet =  * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
715	>	const pat::Jet & pat_topjet =  dynamic_cast<pat::Jet const&>(pat_topjets->at(i));
716		if(pat_topjet.pt() < topjet_ptmin) continue;
717		if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
718
719		TopJet topjet;
720	<	topjet.charge = pat_topjet.charge();
721	<	topjet.pt = pat_topjet.pt();
722	<	topjet.eta = pat_topjet.eta();
723	<	topjet.phi = pat_topjet.phi();
724	<	topjet.energy = pat_topjet.energy();
725	<	topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
726	<	const reco::TrackRefVector&  topjettracks = pat_topjet.associatedTracks();
727	<	topjet.nTracks = topjettracks.size();
728	<	topjet.jetArea = pat_topjet.jetArea();
729	<	topjet.pileup = pat_topjet.pileup();
730	<	//       topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
731	<	//       topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
732	<	//       topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
733	<	//       topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
734	<	//       topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
735	<	//       topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
736	<	//       topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
737	<	//       topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
738	<	//       topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
739	<	//       topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
740	<	//       topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
741	<
742	<	topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
743	<	topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
744	<	topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
745	<	topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
746	<	topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
747	<	topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
720	>	fill_jet_info(pat_topjet, topjet);
721	>
722	>	/*
723	>	const reco::GenJet *genj = pat_topjet.genJet();
724	>	if(genj){
725	>	topjet.set_genjet_pt ( genj->pt());
726	>	topjet.set_genjet_eta ( genj->eta());
727	>	topjet.set_genjet_phi ( genj->phi());
728	>	topjet.set_genjet_energy ( genj->energy());
729	>	if(doAllGenParticles){
730	>	std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
731	>	for(unsigned int l = 0; l<jetgenps.size(); ++l){
732	>	for(unsigned int k=0; k< genps.size(); ++k){
733	>	if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
734	>	topjet.add_genparticles_index(genps[k].index());
735	>	break;
736	>	}
737	>	}
738	>	}
739	>	if(topjet.genparticles_indices().size()!= jetgenps.size())
740	>	std::cout << "WARNING: Found only " << topjet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
741	>	}
742	>	}
743	>	*/
744	>
745	>	// add constituents to the jet, if requested
746	>	if (doTopJetsConstituents){
747	>	if (topjet_constituents_sources.size()>j){ //only add constituents if they are defined
748	>	edm::Handle<pat::JetCollection> pat_topjets_with_cands;
749	>	iEvent.getByLabel(topjet_constituents_sources[j], pat_topjets_with_cands);
750	>	const pat::Jet* pat_topjet_wc = NULL;
751	>
752	>	for (unsigned int it = 0; it < pat_topjets_with_cands->size(); it++) {
753	>	const pat::Jet* cand =  dynamic_cast<pat::Jet const *>(&pat_topjets_with_cands->at(it));
754	>	assert(cand);
755	>	double dphi = deltaPhi(cand->phi(), pat_topjet.phi());
756	>	if (fabs(dphi)<0.5 && fabs(cand->eta()-pat_topjet.eta())<0.5){ // be generous: filtering, pruning... can change jet axis
757	>	pat_topjet_wc = cand;
758	>	break;
759	>	}
760	>	}
761	>
762	>	if (pat_topjet_wc){
763	>	StoreJetConstituents(*pat_topjet_wc, topjet);
764	>	// now run substructure information
765	>	JetProps substructure(&topjet);
766	>	substructure.set_pf_cands(&pfparticles);
767	>	double tau1 = substructure.GetNsubjettiness(1, Njettiness::onepass_kt_axes, 1., 2.0);
768	>	double tau2 = substructure.GetNsubjettiness(2, Njettiness::onepass_kt_axes, 1., 2.0);
769	>	double tau3 = substructure.GetNsubjettiness(3, Njettiness::onepass_kt_axes, 1., 2.0);
770	>	double qjets = substructure.GetQjetVolatility(iEvent.id().event(), 2.0);
771	>	topjet.set_tau1(tau1);
772	>	topjet.set_tau2(tau2);
773	>	topjet.set_tau3(tau3);
774	>	topjet.set_qjets_volatility(qjets);
775	>
776	>	}
777	>	}
778	>	}
779	>
780	>
781
782		for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
783		Particle subjet_v4;
784	<	subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
785	<	subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
786	<	subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
787	<	subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
788	<	topjet.subjets.push_back(subjet_v4);
784	>
785	>	reco::Candidate const * subjetd =  pat_topjet.daughter(k);
786	>	pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
787	>	if(patsubjetd)
788	>	{
789	>	subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
790	>	subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
791	>	subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
792	>	subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
793	>	topjet.add_subjet(subjet_v4);
794	>	//btag info
795	>	topjet.add_subFlavour(patsubjetd->partonFlavour());
796	>	topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
797	>	if (doTagInfos)
798	>	{
799	>	//ip tag info
800	>	reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
801	>	topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
802	>	topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
803	>	topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
804	>	topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
805	>	topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
806	>	topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
807	>	topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
808	>	topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
809	>	topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
810	>	topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
811	>	topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
812	>	topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));
813	>	topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
814	>	topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
815	>	topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
816	>	topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
817	>	topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
818	>	topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
819	>	topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
820	>	topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
821	>	topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
822	>	//sv tag info
823	>	reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
824	>	topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
825	>	topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
826	>	topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
827	>	topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
828	>	topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
829	>	topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
830	>	topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
831	>	std::vector<TLorentzVector> vp4; vp4.clear();
832	>	std::vector<float> vchi2; vchi2.clear();
833	>	std::vector<float> vndof; vndof.clear();
834	>	std::vector<float> vchi2ndof; vchi2ndof.clear();
835	>	std::vector<float> vtrsize; vtrsize.clear();
836	>	for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
837	>	{
838	>	ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
839	>	vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
840	>	vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());
841	>	vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());
842	>	vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());
843	>	vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());
844	>	}
845	>	topjet.add_subSecondaryVertex(vp4);
846	>	topjet.add_subVertexChi2(vchi2);
847	>	topjet.add_subVertexNdof(vndof);
848	>	topjet.add_subVertexNormalizedChi2(vchi2ndof);
849	>	topjet.add_subVertexTracksSize(vtrsize);
850	>	//try computer
851	>	edm::ESHandle<JetTagComputer> computerHandle;
852	>	iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
853	>	const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
854	>	if(computer)
855	>	{
856	>	computer->passEventSetup(iSetup);
857	>	std::vector<const reco::BaseTagInfo*>  baseTagInfos;
858	>	baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
859	>	baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );
860	>	JetTagComputer::TagInfoHelper helper(baseTagInfos);
861	>	reco::TaggingVariableList vars = computer->taggingVariables(helper);
862	>	topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
863	>	topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
864	>	topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
865	>	topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
866	>	}
867	>	}
868	>	}
869	>	else
870	>	{
871	>	//filling only standard information in case the subjet has not been pat-tified during the pattuples production
872	>	subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
873	>	subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
874	>	subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
875	>	subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
876	>	topjet.add_subjet(subjet_v4);
877	>	}
878	>
879	>
880		}
881		topjets[j].push_back(topjet);
882		}
883		}
884		}
885
886	+
887	+	// ------------- generator top jets -------------
888	+	if(doGenTopJets){
889	+	for(size_t j=0; j< gentopjet_sources.size(); ++j){
890	+
891	+	gentopjets[j].clear();
892	+
893	+	edm::Handle<reco::BasicJetCollection> reco_gentopjets;
894	+	iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
895	+
896	+	for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
897	+
898	+	const reco::BasicJet  reco_gentopjet =  reco_gentopjets->at(i);
899	+	if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
900	+	if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
901	+
902	+	GenTopJet gentopjet;
903	+	gentopjet.set_charge(reco_gentopjet.charge());
904	+	gentopjet.set_pt(reco_gentopjet.pt());
905	+	gentopjet.set_eta(reco_gentopjet.eta());
906	+	gentopjet.set_phi(reco_gentopjet.phi());
907	+	gentopjet.set_energy(reco_gentopjet.energy());
908	+
909	+	for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
910	+	Particle subjet_v4;
911	+	subjet_v4.set_pt(reco_gentopjet.daughter(k)->p4().pt());
912	+	subjet_v4.set_eta(reco_gentopjet.daughter(k)->p4().eta());
913	+	subjet_v4.set_phi(reco_gentopjet.daughter(k)->p4().phi());
914	+	subjet_v4.set_energy(reco_gentopjet.daughter(k)->p4().E());
915	+	gentopjet.add_subjet(subjet_v4);
916	+	}
917	+	gentopjets[j].push_back(gentopjet);
918	+	}
919	+	}
920	+	}
921	+
922	+
923		// ------------- photons -------------
924		if(doPhotons){
925		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 446 | Line 932 | NtupleWriter::analyze(const edm::Event&
932		for (unsigned int i = 0; i < pat_photons.size(); ++i) {
933		pat::Photon pat_photon = pat_photons[i];
934		Photon ph;
935	<	ph.charge = 0;
936	<	ph.pt =  pat_photon.pt();
937	<	ph.eta =  pat_photon.eta();
938	<	ph.phi =  pat_photon.phi();
939	<	ph.energy =  pat_photon.energy();
940	<	ph.vertex_x = pat_photon.vertex().x();
941	<	ph.vertex_y = pat_photon.vertex().y();
942	<	ph.vertex_z = pat_photon.vertex().z();
943	<	ph.supercluster_eta = pat_photon.superCluster()->eta();
944	<	ph.supercluster_phi = pat_photon.superCluster()->phi();
945	<	//       ph.neutralHadronIso = pat_photon.neutralHadronIso();
946	<	//       ph.chargedHadronIso = pat_photon.chargedHadronIso();
947	<	ph.trackIso = pat_photon.trackIso();
935	>	ph.set_charge(0);
936	>	ph.set_pt( pat_photon.pt());
937	>	ph.set_eta( pat_photon.eta());
938	>	ph.set_phi( pat_photon.phi());
939	>	ph.set_energy( pat_photon.energy());
940	>	ph.set_vertex_x(pat_photon.vertex().x());
941	>	ph.set_vertex_y(pat_photon.vertex().y());
942	>	ph.set_vertex_z(pat_photon.vertex().z());
943	>	ph.set_supercluster_eta(pat_photon.superCluster()->eta());
944	>	ph.set_supercluster_phi(pat_photon.superCluster()->phi());
945	>	//       ph.set_neutralHadronIso(pat_photon.neutralHadronIso());
946	>	//       ph.set_chargedHadronIso(pat_photon.chargedHadronIso());
947	>	ph.set_trackIso(pat_photon.trackIso());
948		phs[j].push_back(ph);
949		}
950		}
#	Line 478 | Line 964 | NtupleWriter::analyze(const edm::Event&
964		else{
965		pat::MET pat_met = pat_mets[0];
966
967	<	met[j].pt=pat_met.pt();
968	<	met[j].phi=pat_met.phi();
969	<	met[j].mEtSig= pat_met.mEtSig();
967	>	met[j].set_pt(pat_met.pt());
968	>	met[j].set_phi(pat_met.phi());
969	>	met[j].set_mEtSig(pat_met.mEtSig());
970		}
971
972		}
973		}
974
975		// ------------- trigger -------------
976	<
977	<	edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
978	<	edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
979	<	iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
494	<
495	<	const edm::Provenance *meta = dummy_TriggerEvent.provenance();
496	<	std::string nameProcess = meta->processName();
497	<	edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
498	<	triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
499	<
500	<	edm::Handle<edm::TriggerResults> trigger;
501	<	iEvent.getByLabel(triggerResultTag, trigger);
502	<	const edm::TriggerResults& trig = *(trigger.product());
503	<
504	<	triggerResults.clear();
505	<	triggerNames.clear();
506	<	L1_prescale.clear();
507	<	HLT_prescale.clear();
508	<
509	<	edm::Service<edm::service::TriggerNamesService> tns;
510	<	std::vector<std::string> triggerNames_all;
511	<	tns->getTrigPaths(trig,triggerNames_all);
512	<
513	<	if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
514	<	for(unsigned int i=0; i<trig.size(); ++i){
515	<	std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
516	<	for(; it!=trigger_prefixes.end(); ++it){
517	<	if(triggerNames_all[i].substr(0, it->size()) == *it)break;
518	<	}
519	<	if(it==trigger_prefixes.end()) continue;
520	<
521	<	//triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
522	<	triggerResults.push_back(trig.accept(i));
523	<	if(newrun) triggerNames.push_back(triggerNames_all[i]);
524	<	if(isRealData){
525	<	std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
526	<	L1_prescale.push_back(pre.first);
527	<	HLT_prescale.push_back(pre.second);
528	<	}
529	<	}
530	<	//    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
531	<	//      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
532	<	//    }
533	<	newrun=false;
534	<
535	<	// ------------- generator info -------------
536	<
537	<	if(doGenInfo){
538	<	genInfo.weights.clear();
539	<	genInfo.binningValues.clear();
540	<	genps.clear();
541	<
542	<	edm::Handle<GenEventInfoProduct> genEventInfoProduct;
543	<	iEvent.getByLabel("generator", genEventInfoProduct);
544	<	const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
545	<
546	<	genInfo.binningValues = genEventInfo.binningValues();
547	<	genInfo.weights = genEventInfo.weights();
548	<	genInfo.alphaQCD = genEventInfo.alphaQCD();
549	<	genInfo.alphaQED = genEventInfo.alphaQED();
550	<	genInfo.qScale = genEventInfo.qScale();
976	>	if(doTrigger){
977	>	edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
978	>	edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
979	>	iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
980
981	<	const gen::PdfInfo* pdf = genEventInfo.pdf();
982	<	if(pdf){
983	<	genInfo.pdf_id1=pdf->id.first;
984	<	genInfo.pdf_id2=pdf->id.second;
985	<	genInfo.pdf_x1=pdf->x.first;
986	<	genInfo.pdf_x2=pdf->x.second;
987	<	genInfo.pdf_scalePDF=pdf->scalePDF;
988	<	genInfo.pdf_xPDF1=pdf->xPDF.first;
989	<	genInfo.pdf_xPDF2=pdf->xPDF.second;
990	<	}
991	<	else{
992	<	genInfo.pdf_id1=-999;
993	<	genInfo.pdf_id2=-999;
994	<	genInfo.pdf_x1=-999;
995	<	genInfo.pdf_x2=-999;
996	<	genInfo.pdf_scalePDF=-999;
997	<	genInfo.pdf_xPDF1=-999;
998	<	genInfo.pdf_xPDF2=-999;
999	<	}
1000	<
1001	<	edm::Handle<std::vector<PileupSummaryInfo> > pus;
1002	<	iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
1003	<	genInfo.pileup_NumInteractions_intime=0;
575	<	genInfo.pileup_NumInteractions_ootbefore=0;
576	<	genInfo.pileup_NumInteractions_ootafter=0;
577	<	if(pus.isValid()){
578	<	genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
579	<	for(size_t i=0; i<pus->size(); ++i){
580	<	if(pus->at(i).getBunchCrossing() == 0) // intime pileup
581	<	genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
582	<	else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
583	<	genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
584	<	}
585	<	else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
586	<	genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
587	<	}
981	>	const edm::Provenance *meta = dummy_TriggerEvent.provenance();
982	>	std::string nameProcess = meta->processName();
983	>	edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
984	>	triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
985	>
986	>	edm::Handle<edm::TriggerResults> trigger;
987	>	iEvent.getByLabel(triggerResultTag, trigger);
988	>	const edm::TriggerResults& trig = *(trigger.product());
989	>
990	>	triggerResults.clear();
991	>	triggerNames.clear();
992	>	//      L1_prescale.clear();
993	>	//      HLT_prescale.clear();
994	>
995	>	edm::Service<edm::service::TriggerNamesService> tns;
996	>	std::vector<std::string> triggerNames_all;
997	>	tns->getTrigPaths(trig,triggerNames_all);
998	>
999	>	if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
1000	>	for(unsigned int i=0; i<trig.size(); ++i){
1001	>	std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
1002	>	for(; it!=trigger_prefixes.end(); ++it){
1003	>	if(triggerNames_all[i].substr(0, it->size()) == *it)break;
1004		}
1005	<	}
590	<
591	<	edm::Handle<reco::GenParticleCollection> genPartColl;
592	<	iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
593	<	int index=-1;
594	<	for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
595	<	index++;
1005	>	if(it==trigger_prefixes.end()) continue;
1006
1007	<	//write out only top quarks and status 3 particles (works fine only for MadGraph)
1008	<	if(abs(iter->pdgId())==6 || iter->status()==3){
1009	<	GenParticle genp;
1010	<	genp.charge = iter->charge();;
1011	<	genp.pt = iter->p4().pt();
1012	<	genp.eta = iter->p4().eta();
1013	<	genp.phi = iter->p4().phi();
1014	<	genp.energy = iter->p4().E();
1015	<	genp.index =index;
606	<	genp.status = iter->status();
607	<	genp.pdgId = iter->pdgId();
608	<
609	<	genp.mother1=-1;
610	<	genp.mother2=-1;
611	<	genp.daughter1=-1;
612	<	genp.daughter2=-1;
613	<
614	<	int nm=iter->numberOfMothers();
615	<	int nd=iter->numberOfDaughters();
616	<
617	<	if (nm>0) genp.mother1 = iter->motherRef(0).key();
618	<	if (nm>1) genp.mother2 = iter->motherRef(nm-1).key();
619	<	if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
620	<	if (nd>1) genp.daughter2 = iter->daughterRef(nd-1).key();
621	<
622	<	genps.push_back(genp);
623	<	}
1007	>	//triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
1008	>	triggerResults.push_back(trig.accept(i));
1009	>	if(newrun) triggerNames.push_back(triggerNames_all[i]);
1010	>	//        if(isRealData){
1011	>	//       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
1012	>	//       L1_prescale.push_back(pre.first);
1013	>	//       HLT_prescale.push_back(pre.second);
1014	>	//       //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
1015	>	//        }
1016		}
1017	<
1017	>	//    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
1018	>	//      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
1019	>	//    }
1020	>	newrun=false;
1021		}
1022
628	–	// ------------- primary vertices and beamspot  -------------
1023
1024	<	if(doPV){
1025	<	for(size_t j=0; j< pv_sources.size(); ++j){
1026	<	pvs[j].clear();
633	<
634	<	edm::Handle< std::vector<reco::Vertex> > pv_handle;
635	<	iEvent.getByLabel(pv_sources[j], pv_handle);
636	<	const std::vector<reco::Vertex>& reco_pvs = *(pv_handle.product());
637	<
638	<	for (unsigned int i = 0; i <  reco_pvs.size(); ++i) {
639	<	reco::Vertex reco_pv = reco_pvs[i];
640	<
641	<	PrimaryVertex pv;
642	<	pv.x =  reco_pv.x();
643	<	pv.y =  reco_pv.y();
644	<	pv.z =  reco_pv.z();
645	<	pv.nTracks =  reco_pv.nTracks();
646	<	//pv.isValid =  reco_pv.isValid();
647	<	pv.chi2 =  reco_pv.chi2();
648	<	pv.ndof =  reco_pv.ndof();
1024	>	tr->Fill();
1025	>	if(doLumiInfo)
1026	>	previouslumiblockwasfilled=true;
1027
1028	<	pvs[j].push_back(pv);
1029	<	}
1030	<	}
1028	>	// clean up
1029	>	if(doTopJetsConstituents){
1030	>	pfparticles.clear();
1031		}
1032
1033	<	edm::Handle<reco::BeamSpot> beamSpot;
1034	<	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
657	<	const reco::BeamSpot & bsp = *beamSpot;
658	<
659	<	beamspot_x0 = bsp.x0();
660	<	beamspot_y0 = bsp.y0();
661	<	beamspot_z0 = bsp.z0();
662	<
663	<	tr->Fill();
1033	>	// need to do a check if necessary
1034	>	pfparticles.clear();
1035
1036		}
1037
#	Line 669 | Line 1040 | NtupleWriter::analyze(const edm::Event&
1040		void
1041		NtupleWriter::beginJob()
1042		{
1043	+	if(doLumiInfo){
1044	+	totalRecLumi=0;
1045	+	totalDelLumi=0;
1046	+	previouslumiblockwasfilled=false;
1047	+	}
1048		}
1049
1050		// ------------ method called once each job just after ending the event loop  ------------
#	Line 684 | Line 1060 | NtupleWriter::endJob()
1060		void
1061		NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
1062		{
1063	<	bool setup_changed = false;
1064	<	hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1065	<	newrun=true;
1063	>	if(doTrigger){
1064	>	//bool setup_changed = false;
1065	>	//hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
1066	>	newrun=true;
1067	>	}
1068	>
1069		}
1070
1071		// ------------ method called when ending the processing of a run  ------------
1072	<	void
694	<	NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
1072	>	void NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
1073		{
1074	+	if(doLumiInfo)
1075	+	std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
1076		}
1077
1078		// ------------ method called when starting to processes a luminosity block  ------------
1079	<	void
700	<	NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
701	<	{
702	<	}
703	<
704	<	// ------------ method called when ending the processing of a luminosity block  ------------
705	<	void
706	<	NtupleWriter::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
1079	>	void NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
1080		{
1081	+	if(doLumiInfo){
1082	+	edm::Handle<LumiSummary> l;
1083	+	lumi.getByLabel("lumiProducer", l);
1084	+
1085	+	//add lumi of lumi blocks without any event to next lumiblock
1086	+	if(previouslumiblockwasfilled){
1087	+	intgRecLumi=0;
1088	+	intgDelLumi=0;
1089	+	}
1090	+	previouslumiblockwasfilled=false;
1091	+
1092	+	if (l.isValid()){;
1093	+	intgRecLumi+=l->intgRecLumi()*6.37;
1094	+	intgDelLumi+=l->intgDelLumi()*6.37;
1095	+	totalRecLumi+=l->intgRecLumi()*6.37;
1096	+	totalDelLumi+=l->intgDelLumi()*6.37;
1097	+	}
1098	+	//std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<"   " << l->intgDelLumi()*6.37<<std::endl;
1099	+	//std::cout << "summed:  "<< intgRecLumi << "   " << intgDelLumi << std::endl;
1100	+	}
1101		}
1102
1103		// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
#	Line 717 | Line 1110 | NtupleWriter::fillDescriptions(edm::Conf
1110		descriptions.addDefault(desc);
1111		}
1112
1113	+
1114	+	namespace{
1115	+
1116	+	PFParticle PFCandidate2PFParticle(const reco::PFCandidate & pf){
1117	+	PFParticle part;
1118	+	part.set_pt(pf.pt());
1119	+	part.set_eta(pf.eta());
1120	+	part.set_phi(pf.phi());
1121	+	part.set_energy(pf.energy());
1122	+	part.set_charge(pf.charge());
1123	+	part.set_ecal_en(pf.ecalEnergy());
1124	+	part.set_hcal_en(pf.hcalEnergy());
1125	+	reco::TrackRef trackref = pf.trackRef();
1126	+	if (!trackref.isNull()){
1127	+	part.set_track_mom(trackref->p());
1128	+	}
1129	+	PFParticle::EParticleID id = PFParticle::eX;
1130	+	switch ( pf.translatePdgIdToType(pf.pdgId()) ){
1131	+	case reco::PFCandidate::X : id = PFParticle::eX; break;
1132	+	case reco::PFCandidate::h : id = PFParticle::eH; break;
1133	+	case reco::PFCandidate::e : id = PFParticle::eE; break;
1134	+	case reco::PFCandidate::mu : id = PFParticle::eMu; break;
1135	+	case reco::PFCandidate::gamma : id = PFParticle::eGamma; break;
1136	+	case reco::PFCandidate::h0 : id = PFParticle::eH0; break;
1137	+	case reco::PFCandidate::h_HF : id = PFParticle::eH_HF; break;
1138	+	case reco::PFCandidate::egamma_HF : id = PFParticle::eEgamma_HF; break;
1139	+	}
1140	+	part.set_particleID(id);
1141	+	return part;
1142	+	}
1143	+
1144	+	// add pf to pfs, ensuring there is no duplication. Retuns the index
1145	+	// of pf in pfs.
1146	+	size_t add_pfpart(const reco::PFCandidate & pf, vector<PFParticle> & pfs){
1147	+	for(size_t j=0; j<pfs.size(); ++j){
1148	+	PFParticle spf = pfs[j];
1149	+	// note: static_cast to float is to ensure the comparison is done with the same precision as these quantities
1150	+	// have been stored in spf. Otherwise, there could be a non-zero difference just because of conversion loss
1151	+	// from double to float.
1152	+	double r = fabs(static_cast<float>(pf.eta()) - spf.eta()) + fabs(static_cast<float>(pf.phi()) - spf.phi());
1153	+	double dpt = fabs(static_cast<float>(pf.pt()) - spf.pt());
1154	+	if (r == 0.0 && dpt == 0.0){
1155	+	return j;
1156	+	}
1157	+	}
1158	+	pfs.push_back(PFCandidate2PFParticle(pf));
1159	+	return pfs.size()-1;
1160	+	}
1161	+
1162	+
1163	+	}
1164	+
1165	+
1166	+	void NtupleWriter::StoreJetConstituents(const pat::Jet& pat_jet, Jet& jet)
1167	+	{
1168	+	const std::vector<reco::PFCandidatePtr> jconstits = pat_jet.getPFConstituents();
1169	+	for (unsigned int i=0; i<jconstits.size(); ++i){
1170	+	const reco::PFCandidate* pf = jconstits[i].get();
1171	+	size_t pfparticles_index = add_pfpart(*pf, pfparticles);
1172	+	jet.add_pfconstituents_index(pfparticles_index);
1173	+	}
1174	+	}
1175	+
1176	+	void NtupleWriter::StorePFCandsInCone(Particle* inpart, const std::vector<reco::PFCandidate>& pf_cands, double R0)
1177	+	{
1178	+	for (unsigned int i=0; i<pf_cands.size(); ++i){
1179	+	const reco::PFCandidate & pf = pf_cands[i];
1180	+	double dr = deltaR(*inpart, pf);
1181	+	if (dr>R0) continue;
1182	+
1183	+	// don't store particle itself:
1184	+	double dpt = fabs( inpart->pt() - pf.pt() );
1185	+	if (dr<1e-10 && dpt<1e-10){
1186	+	continue;
1187	+	}
1188	+	add_pfpart(pf, pfparticles);
1189	+	}
1190	+	}
1191	+
1192		//define this as a plug-in
1193		DEFINE_FWK_MODULE(NtupleWriter);

Diff Legend

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