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.16 by peiffer, Mon Apr 30 09:01:16 2012 UTC

#	Line 50 | Line 50 | NtupleWriter::NtupleWriter(const edm::Pa
50		doMuons = iConfig.getParameter<bool>("doMuons");
51		doTaus = iConfig.getParameter<bool>("doTaus");
52		doJets = iConfig.getParameter<bool>("doJets");
53	+	doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");
54		doPhotons = iConfig.getParameter<bool>("doPhotons");
55		doMET = iConfig.getParameter<bool>("doMET");
56		doGenInfo = iConfig.getParameter<bool>("doGenInfo");
57	+	doAllGenParticles = iConfig.getParameter<bool>("doAllGenParticles");
58	+	doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
59		doPV = iConfig.getParameter<bool>("doPV");
60		doTopJets = iConfig.getParameter<bool>("doTopJets");
61	+	doTrigger = iConfig.getParameter<bool>("doTrigger");
62
63		// initialization of tree variables
64
#	Line 63 | Line 67 | NtupleWriter::NtupleWriter(const edm::Pa
67		tr->Branch("luminosityBlock",&luminosityBlock);
68		tr->Branch("isRealData",&isRealData);
69		tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
70	<	tr->Branch("beamspot_x0",&beamspot_x0);
71	<	tr->Branch("beamspot_y0",&beamspot_y0);
72	<	tr->Branch("beamspot_z0",&beamspot_z0);
73	<
70	>	if(doLumiInfo){
71	>	tr->Branch("intgRecLumi",&intgRecLumi);
72	>	tr->Branch("intgDelLumi",&intgDelLumi);
73	>	}
74	>	if(doPV){
75	>	tr->Branch("beamspot_x0",&beamspot_x0);
76	>	tr->Branch("beamspot_y0",&beamspot_y0);
77	>	tr->Branch("beamspot_z0",&beamspot_z0);
78	>	}
79		if(doElectrons){
80		electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
81		for(size_t j=0; j< electron_sources.size(); ++j){
#	Line 103 | Line 112 | NtupleWriter::NtupleWriter(const edm::Pa
112		tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
113		}
114		}
115	+	if(doGenTopJets){
116	+	gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
117	+	gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
118	+	gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
119	+	for(size_t j=0; j< gentopjet_sources.size(); ++j){
120	+	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<TopJet>", &gentopjets[j]);
121	+	}
122	+	}
123		if(doPhotons){
124		photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
125		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 125 | Line 142 | NtupleWriter::NtupleWriter(const edm::Pa
142		tr->Branch("genInfo","GenInfo",&genInfo);
143		tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
144		}
145	<
146	<	trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
147	<	//tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
148	<	tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
149	<	tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
150	<	tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
151	<	tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
152	<
145	>	if(doTrigger){
146	>	trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
147	>	//tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
148	>	tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
149	>	tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
150	>	//tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
151	>	//tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
152	>	}
153		newrun = true;
154		}
155
#	Line 196 | Line 213 | NtupleWriter::analyze(const edm::Event&
213		pvs[j].push_back(pv);
214		}
215		}
216	+
217	+	edm::Handle<reco::BeamSpot> beamSpot;
218	+	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
219	+	const reco::BeamSpot & bsp = *beamSpot;
220	+
221	+	beamspot_x0 = bsp.x0();
222	+	beamspot_y0 = bsp.y0();
223	+	beamspot_z0 = bsp.z0();
224		}
225	+
226	+	// ------------- generator info -------------
227
228	<	edm::Handle<reco::BeamSpot> beamSpot;
229	<	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
230	<	const reco::BeamSpot & bsp = *beamSpot;
231	<
232	<	beamspot_x0 = bsp.x0();
233	<	beamspot_y0 = bsp.y0();
234	<	beamspot_z0 = bsp.z0();
228	>	if(doGenInfo){
229	>	genInfo.weights.clear();
230	>	genInfo.binningValues.clear();
231	>	genps.clear();
232	>
233	>	edm::Handle<GenEventInfoProduct> genEventInfoProduct;
234	>	iEvent.getByLabel("generator", genEventInfoProduct);
235	>	const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
236	>
237	>	genInfo.binningValues = genEventInfo.binningValues();
238	>	genInfo.weights = genEventInfo.weights();
239	>	genInfo.alphaQCD = genEventInfo.alphaQCD();
240	>	genInfo.alphaQED = genEventInfo.alphaQED();
241	>	genInfo.qScale = genEventInfo.qScale();
242	>
243	>	const gen::PdfInfo* pdf = genEventInfo.pdf();
244	>	if(pdf){
245	>	genInfo.pdf_id1=pdf->id.first;
246	>	genInfo.pdf_id2=pdf->id.second;
247	>	genInfo.pdf_x1=pdf->x.first;
248	>	genInfo.pdf_x2=pdf->x.second;
249	>	genInfo.pdf_scalePDF=pdf->scalePDF;
250	>	genInfo.pdf_xPDF1=pdf->xPDF.first;
251	>	genInfo.pdf_xPDF2=pdf->xPDF.second;
252	>	}
253	>	else{
254	>	genInfo.pdf_id1=-999;
255	>	genInfo.pdf_id2=-999;
256	>	genInfo.pdf_x1=-999;
257	>	genInfo.pdf_x2=-999;
258	>	genInfo.pdf_scalePDF=-999;
259	>	genInfo.pdf_xPDF1=-999;
260	>	genInfo.pdf_xPDF2=-999;
261	>	}
262	>
263	>	edm::Handle<std::vector<PileupSummaryInfo> > pus;
264	>	iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
265	>	genInfo.pileup_NumInteractions_intime=0;
266	>	genInfo.pileup_NumInteractions_ootbefore=0;
267	>	genInfo.pileup_NumInteractions_ootafter=0;
268	>	if(pus.isValid()){
269	>	genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
270	>	for(size_t i=0; i<pus->size(); ++i){
271	>	if(pus->at(i).getBunchCrossing() == 0) // intime pileup
272	>	genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
273	>	else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
274	>	genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
275	>	}
276	>	else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
277	>	genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
278	>	}
279	>	}
280	>	}
281	>
282	>	edm::Handle<reco::GenParticleCollection> genPartColl;
283	>	iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
284	>	int index=-1;
285	>	for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
286	>	index++;
287	>
288	>	//write out only top quarks and status 3 particles (works fine only for MadGraph)
289	>	if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
290	>	GenParticle genp;
291	>	genp.charge = iter->charge();;
292	>	genp.pt = iter->p4().pt();
293	>	genp.eta = iter->p4().eta();
294	>	genp.phi = iter->p4().phi();
295	>	genp.energy = iter->p4().E();
296	>	genp.index =index;
297	>	genp.status = iter->status();
298	>	genp.pdgId = iter->pdgId();
299	>
300	>	genp.mother1=-1;
301	>	genp.mother2=-1;
302	>	genp.daughter1=-1;
303	>	genp.daughter2=-1;
304	>
305	>	int nm=iter->numberOfMothers();
306	>	int nd=iter->numberOfDaughters();
307	>
308	>
309	>	if (nm>0) genp.mother1 = iter->motherRef(0).key();
310	>	if (nm>1) genp.mother2 = iter->motherRef(1).key();
311	>	if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
312	>	if (nd>1) genp.daughter2 = iter->daughterRef(1).key();
313	>	//std::cout << genp.index <<"  pdgId = " << genp.pdgId << "  mo1 = " << genp.mother1 << "  mo2 = " << genp.mother2 <<"  da1 = " << genp.daughter1 << "  da2 = " << genp.daughter2 <<std::endl;
314	>	genps.push_back(genp);
315	>	}
316	>	}
317	>	}
318
319		// ------------- electrons -------------
320		if(doElectrons){
321	+
322	+	edm::Handle<reco::ConversionCollection> hConversions;
323	+	iEvent.getByLabel("allConversions", hConversions);
324	+
325	+	edm::Handle<reco::BeamSpot> beamSpot;
326	+	iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
327	+	const reco::BeamSpot & bsp = *beamSpot;
328	+
329		for(size_t j=0; j< electron_sources.size(); ++j){
330		eles[j].clear();
331		edm::Handle< std::vector<pat::Electron> > ele_handle;
#	Line 233 | Line 351 | NtupleWriter::analyze(const edm::Event&
351		ele.neutralHadronIso = pat_ele.neutralHadronIso();
352		ele.chargedHadronIso = pat_ele.chargedHadronIso();
353		ele.trackIso = pat_ele.trackIso();
354	+	ele.photonIso = pat_ele.photonIso();
355		ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
356		ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
357		ele.gsfTrack_px= pat_ele.gsfTrack()->px();
#	Line 241 | Line 360 | NtupleWriter::analyze(const edm::Event&
360		ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
361		ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
362		ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
363	+	ele.passconversionveto = !ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position());
364	+	ele.dEtaIn=pat_ele.deltaEtaSuperClusterTrackAtVtx();
365	+	ele.dPhiIn=pat_ele.deltaPhiSuperClusterTrackAtVtx();
366	+	ele.sigmaIEtaIEta=pat_ele.sigmaIetaIeta();
367	+	ele.HoverE=pat_ele.hadronicOverEm();
368	+	ele.fbrem=pat_ele.fbrem();
369	+	ele.EoverPIn=pat_ele.eSuperClusterOverP();
370	+	ele.EcalEnergy=pat_ele.ecalEnergy();
371	+	//ele.mvaTrigV0=pat_ele.electronID("mvaTrigV0");
372	+	//ele.mvaNonTrigV0=pat_ele.electronID("mvaNonTrigV0");
373	+
374		eles[j].push_back(ele);
375		}
376		}
#	Line 272 | Line 402 | NtupleWriter::analyze(const edm::Event&
402		mu.neutralHadronIso = pat_mu.neutralHadronIso();
403		mu.chargedHadronIso = pat_mu.chargedHadronIso();
404		mu.trackIso = pat_mu.trackIso();
405	+	mu.photonIso = pat_mu.photonIso();
406		mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
407		mu.isGlobalMuon = pat_mu.isGlobalMuon();
408		mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
#	Line 302 | Line 433 | NtupleWriter::analyze(const edm::Event&
433		mu.innerTrack_d0Error = innerTrack->d0Error();
434		mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
435		mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
436	+	mu.innerTrack_trackerLayersWithMeasurement = innerTrack->hitPattern().trackerLayersWithMeasurement();
437	+	mu.innerTrack_numberOfValidPixelHits = innerTrack->hitPattern().numberOfValidPixelHits();
438		}
439		else{
440		mu.innerTrack_chi2 = 0;
#	Line 310 | Line 443 | NtupleWriter::analyze(const edm::Event&
443		mu.innerTrack_d0Error = 0;
444		mu.innerTrack_numberOfValidHits = 0;
445		mu.innerTrack_numberOfLostHits = 0;
446	+	mu.innerTrack_trackerLayersWithMeasurement = 0;
447	+	mu.innerTrack_numberOfValidPixelHits = 0;
448		}
449		reco::TrackRef outerTrack = pat_mu.outerTrack();
450		if(!outerTrack.isNull()){
#	Line 355 | Line 490 | NtupleWriter::analyze(const edm::Event&
490		tau.eta =  pat_tau.eta();
491		tau.phi =  pat_tau.phi();
492		tau.energy =  pat_tau.energy();
493	<
493	>	tau.decayModeFinding = pat_tau.tauID("decayModeFinding")>0.5;
494	>	tau.byVLooseCombinedIsolationDeltaBetaCorr  = pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5;
495	>	tau.byLooseCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5;
496	>	tau.byMediumCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5;
497	>	tau.byTightCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5;
498	>	tau.againstElectronLoose  = pat_tau.tauID("againstElectronLoose")>0.5;
499	>	tau.againstElectronMedium = pat_tau.tauID("againstElectronMedium")>0.5;
500	>	tau.againstElectronTight = pat_tau.tauID("againstElectronTight")>0.5;
501	>	tau.againstElectronMVA  = pat_tau.tauID("againstElectronMVA")>0.5;
502	>	tau.againstMuonLoose = pat_tau.tauID("againstMuonLoose")>0.5;
503	>	tau.againstMuonMedium = pat_tau.tauID("againstMuonMedium")>0.5;
504	>	tau.againstMuonTight = pat_tau.tauID("againstMuonTight")>0.5;
505	>
506	>	reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
507	>	if(!leadPFCand.isNull()){
508	>	tau.leadPFCand_px = leadPFCand->px();
509	>	tau.leadPFCand_py = leadPFCand->py();
510	>	tau.leadPFCand_pz = leadPFCand->pz();
511	>	}
512	>	else{
513	>	tau.leadPFCand_px = 0;
514	>	tau.leadPFCand_py = 0;
515	>	tau.leadPFCand_pz = 0;
516	>	}
517		taus[j].push_back(tau);
518		}
519		}
#	Line 409 | Line 567 | NtupleWriter::analyze(const edm::Event&
567		jecUnc->setJetEta(pat_jet.eta());
568		jecUnc->setJetPt(pat_jet.pt());
569		jet.JEC_uncertainty = jecUnc->getUncertainty(true);
570	+	jet.JEC_factor_raw = pat_jet.jecFactor("Uncorrected");
571
572		jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
573		jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
#	Line 417 | Line 576 | NtupleWriter::analyze(const edm::Event&
576		jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
577		jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
578
579	+	const reco::GenJet *genj = pat_jet.genJet();
580	+	if(genj){
581	+	jet.genjet_pt = genj->pt();
582	+	jet.genjet_eta = genj->eta();
583	+	jet.genjet_phi = genj->phi();
584	+	jet.genjet_energy = genj->energy();
585	+	if(doAllGenParticles){
586	+	std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
587	+	for(unsigned int l = 0; l<jetgenps.size(); ++l){
588	+	for(unsigned int k=0; k< genps.size(); ++k){
589	+	if(jetgenps[l]->pt() == genps[k].pt && jetgenps[l]->pdgId() == genps[k].pdgId){
590	+	jet.genparticles_indices.push_back(genps[k].index);
591	+	}
592	+	}
593	+	}
594	+	if(jet.genparticles_indices.size()!= jetgenps.size())
595	+	std::cout << "WARNING: Found only " << jet.genparticles_indices.size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
596	+	}
597	+	}
598		jets[j].push_back(jet);
599		}
600		}
#	Line 433 | Line 611 | NtupleWriter::analyze(const edm::Event&
611		iEvent.getByLabel(topjet_sources[j], pat_topjets);
612
613		for (unsigned int i = 0; i < pat_topjets->size(); i++) {
614	+
615		const pat::Jet  pat_topjet =  * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
616		if(pat_topjet.pt() < topjet_ptmin) continue;
617		if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
#	Line 463 | Line 642 | NtupleWriter::analyze(const edm::Event&
642		jecUnc->setJetEta(pat_topjet.eta());
643		jecUnc->setJetPt(pat_topjet.pt());
644		topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
645	+	topjet.JEC_factor_raw = pat_topjet.jecFactor("Uncorrected");
646
647		topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
648		topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
#	Line 471 | Line 651 | NtupleWriter::analyze(const edm::Event&
651		topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
652		topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
653
654	+	const reco::GenJet *genj = pat_topjet.genJet();
655	+	if(genj){
656	+	topjet.genjet_pt = genj->pt();
657	+	topjet.genjet_eta = genj->eta();
658	+	topjet.genjet_phi = genj->phi();
659	+	topjet.genjet_energy = genj->energy();
660	+	if(doAllGenParticles){
661	+	std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
662	+	for(unsigned int l = 0; l<jetgenps.size(); ++l){
663	+	for(unsigned int k=0; k< genps.size(); ++k){
664	+	if(jetgenps[l]->pt() == genps[k].pt && jetgenps[l]->pdgId() == genps[k].pdgId){
665	+	topjet.genparticles_indices.push_back(genps[k].index);
666	+	}
667	+	}
668	+	}
669	+	if(topjet.genparticles_indices.size()!= jetgenps.size())
670	+	std::cout << "WARNING: Found only " << topjet.genparticles_indices.size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
671	+	}
672	+	}
673	+
674		for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
675		Particle subjet_v4;
676		subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
#	Line 484 | Line 684 | NtupleWriter::analyze(const edm::Event&
684		}
685		}
686
687	+
688	+	// ------------- generator top jets -------------
689	+	if(doGenTopJets){
690	+	for(size_t j=0; j< gentopjet_sources.size(); ++j){
691	+
692	+	gentopjets[j].clear();
693	+
694	+	edm::Handle<reco::BasicJetCollection> reco_gentopjets;
695	+	//edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
696	+	iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
697	+
698	+	for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
699	+
700	+	const reco::BasicJet  reco_gentopjet =  reco_gentopjets->at(i);
701	+	if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
702	+	if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
703	+
704	+	TopJet gentopjet;
705	+	gentopjet.charge = reco_gentopjet.charge();
706	+	gentopjet.pt = reco_gentopjet.pt();
707	+	gentopjet.eta = reco_gentopjet.eta();
708	+	gentopjet.phi = reco_gentopjet.phi();
709	+	gentopjet.energy = reco_gentopjet.energy();
710	+	gentopjet.numberOfDaughters =reco_gentopjet.numberOfDaughters();
711	+
712	+	for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
713	+	Particle subjet_v4;
714	+	subjet_v4.pt = reco_gentopjet.daughter(k)->p4().pt();
715	+	subjet_v4.eta = reco_gentopjet.daughter(k)->p4().eta();
716	+	subjet_v4.phi = reco_gentopjet.daughter(k)->p4().phi();
717	+	subjet_v4.energy = reco_gentopjet.daughter(k)->p4().E();
718	+	gentopjet.subjets.push_back(subjet_v4);
719	+	}
720	+	gentopjets[j].push_back(gentopjet);
721	+	}
722	+	}
723	+	}
724	+
725		// ------------- photons -------------
726		if(doPhotons){
727		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 537 | Line 775 | NtupleWriter::analyze(const edm::Event&
775		}
776
777		// ------------- trigger -------------
778	<
779	<	edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
780	<	edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
781	<	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();
778	>	if(doTrigger){
779	>	edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
780	>	edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
781	>	iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
782
783	<	const gen::PdfInfo* pdf = genEventInfo.pdf();
784	<	if(pdf){
785	<	genInfo.pdf_id1=pdf->id.first;
786	<	genInfo.pdf_id2=pdf->id.second;
787	<	genInfo.pdf_x1=pdf->x.first;
788	<	genInfo.pdf_x2=pdf->x.second;
789	<	genInfo.pdf_scalePDF=pdf->scalePDF;
790	<	genInfo.pdf_xPDF1=pdf->xPDF.first;
791	<	genInfo.pdf_xPDF2=pdf->xPDF.second;
792	<	}
793	<	else{
794	<	genInfo.pdf_id1=-999;
795	<	genInfo.pdf_id2=-999;
796	<	genInfo.pdf_x1=-999;
797	<	genInfo.pdf_x2=-999;
798	<	genInfo.pdf_scalePDF=-999;
799	<	genInfo.pdf_xPDF1=-999;
800	<	genInfo.pdf_xPDF2=-999;
801	<	}
802	<
803	<	edm::Handle<std::vector<PileupSummaryInfo> > pus;
804	<	iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
805	<	genInfo.pileup_NumInteractions_intime=0;
806	<	genInfo.pileup_NumInteractions_ootbefore=0;
807	<	genInfo.pileup_NumInteractions_ootafter=0;
808	<	if(pus.isValid()){
809	<	genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
810	<	for(size_t i=0; i<pus->size(); ++i){
811	<	if(pus->at(i).getBunchCrossing() == 0) // intime pileup
812	<	genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
813	<	else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
814	<	genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
815	<	}
816	<	else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
817	<	genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
818	<	}
819	<	}
820	<	}
821	<
822	<	edm::Handle<reco::GenParticleCollection> genPartColl;
642	<	iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
643	<	int index=-1;
644	<	for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
645	<	index++;
646	<
647	<	//write out only top quarks and status 3 particles (works fine only for MadGraph)
648	<	if(abs(iter->pdgId())==6 || iter->status()==3){
649	<	GenParticle genp;
650	<	genp.charge = iter->charge();;
651	<	genp.pt = iter->p4().pt();
652	<	genp.eta = iter->p4().eta();
653	<	genp.phi = iter->p4().phi();
654	<	genp.energy = iter->p4().E();
655	<	genp.index =index;
656	<	genp.status = iter->status();
657	<	genp.pdgId = iter->pdgId();
658	<
659	<	genp.mother1=-1;
660	<	genp.mother2=-1;
661	<	genp.daughter1=-1;
662	<	genp.daughter2=-1;
663	<
664	<	int nm=iter->numberOfMothers();
665	<	int nd=iter->numberOfDaughters();
666	<
667	<	if (nm>0) genp.mother1 = iter->motherRef(0).key();
668	<	if (nm>1) genp.mother2 = iter->motherRef(nm-1).key();
669	<	if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
670	<	if (nd>1) genp.daughter2 = iter->daughterRef(nd-1).key();
671	<
672	<	genps.push_back(genp);
673	<	}
674	<	}
675	<
783	>	const edm::Provenance *meta = dummy_TriggerEvent.provenance();
784	>	std::string nameProcess = meta->processName();
785	>	edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
786	>	triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
787	>
788	>	edm::Handle<edm::TriggerResults> trigger;
789	>	iEvent.getByLabel(triggerResultTag, trigger);
790	>	const edm::TriggerResults& trig = *(trigger.product());
791	>
792	>	triggerResults.clear();
793	>	triggerNames.clear();
794	>	//      L1_prescale.clear();
795	>	//      HLT_prescale.clear();
796	>
797	>	edm::Service<edm::service::TriggerNamesService> tns;
798	>	std::vector<std::string> triggerNames_all;
799	>	tns->getTrigPaths(trig,triggerNames_all);
800	>
801	>	if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
802	>	for(unsigned int i=0; i<trig.size(); ++i){
803	>	std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
804	>	for(; it!=trigger_prefixes.end(); ++it){
805	>	if(triggerNames_all[i].substr(0, it->size()) == *it)break;
806	>	}
807	>	if(it==trigger_prefixes.end()) continue;
808	>
809	>	//triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
810	>	triggerResults.push_back(trig.accept(i));
811	>	if(newrun) triggerNames.push_back(triggerNames_all[i]);
812	>	//        if(isRealData){
813	>	//       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
814	>	//       L1_prescale.push_back(pre.first);
815	>	//       HLT_prescale.push_back(pre.second);
816	>	//       //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
817	>	//        }
818	>	}
819	>	//    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
820	>	//      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
821	>	//    }
822	>	newrun=false;
823		}
824
825
826		tr->Fill();
827	<
827	>	if(doLumiInfo)
828	>	previouslumiblockwasfilled=true;
829		}
830
831
#	Line 685 | Line 833 | NtupleWriter::analyze(const edm::Event&
833		void
834		NtupleWriter::beginJob()
835		{
836	+	if(doLumiInfo){
837	+	totalRecLumi=0;
838	+	totalDelLumi=0;
839	+	previouslumiblockwasfilled=false;
840	+	}
841		}
842
843		// ------------ method called once each job just after ending the event loop  ------------
#	Line 700 | Line 853 | NtupleWriter::endJob()
853		void
854		NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
855		{
856	<	bool setup_changed = false;
857	<	hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
858	<	newrun=true;
859	<
860	<	edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
708	<	iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
709	<	JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
710	<	jecUnc = new JetCorrectionUncertainty(JetCorPar);
856	>	if(doTrigger){
857	>	bool setup_changed = false;
858	>	hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
859	>	newrun=true;
860	>	}
861
862	+	if(doJets || doTopJets){
863	+	edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
864	+	iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
865	+	JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
866	+	jecUnc = new JetCorrectionUncertainty(JetCorPar);
867	+	}
868		}
869
870		// ------------ method called when ending the processing of a run  ------------
871		void
872		NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
873		{
874	+	if(doLumiInfo)
875	+	std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
876		}
877
878		// ------------ method called when starting to processes a luminosity block  ------------
879		void
880	<	NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
880	>	NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
881		{
882	+	if(doLumiInfo){
883	+	edm::Handle<LumiSummary> l;
884	+	lumi.getByLabel("lumiProducer", l);
885	+
886	+	//add lumi of lumi blocks without any event to next lumiblock
887	+	if(previouslumiblockwasfilled){
888	+	intgRecLumi=0;
889	+	intgDelLumi=0;
890	+	}
891	+	previouslumiblockwasfilled=false;
892	+
893	+	if (l.isValid()){;
894	+	intgRecLumi+=l->intgRecLumi()*6.37;
895	+	intgDelLumi+=l->intgDelLumi()*6.37;
896	+	totalRecLumi+=l->intgRecLumi()*6.37;
897	+	totalDelLumi+=l->intgDelLumi()*6.37;
898	+	}
899	+	//std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<"   " << l->intgDelLumi()*6.37<<std::endl;
900	+	//std::cout << "summed:  "<< intgRecLumi << "   " << intgDelLumi << std::endl;
901	+	}
902		}
903
904		// ------------ method called when ending the processing of a luminosity block  ------------

Diff Legend

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