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.21 by peiffer, Tue Jun 26 08:13:28 2012 UTC vs.
Revision 1.29 by rkogler, Wed Jun 19 22:05:46 2013 UTC

#	Line 17 | Line 17
17		//
18		//
19
20	<	#include "UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
21	<
22	<
20	>	#include "UHHAnalysis/NtupleWriter/interface/NtupleWriter.h"
21	>	#include "UHHAnalysis/NtupleWriter/interface/JetProps.h"
22	>	#include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputer.h"
23	>	#include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputerWrapper.h"
24	>	#include "RecoBTau/JetTagComputer/interface/JetTagComputer.h"
25	>	#include "RecoBTau/JetTagComputer/interface/JetTagComputerRecord.h"
26	>	#include "RecoBTag/SecondaryVertex/interface/CombinedSVComputer.h"
27	>	#include "DataFormats/TrackReco/interface/Track.h"
28
29		//
30		// constants, enums and typedefs
#	Line 50 | Line 55 | NtupleWriter::NtupleWriter(const edm::Pa
55		doMuons = iConfig.getParameter<bool>("doMuons");
56		doTaus = iConfig.getParameter<bool>("doTaus");
57		doJets = iConfig.getParameter<bool>("doJets");
58	<	doJECUncertainty = iConfig.getParameter<bool>("doJECUncertainty");
58	>	doGenJets = iConfig.getParameter<bool>("doGenJets");
59		doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");
60		doPhotons = iConfig.getParameter<bool>("doPhotons");
61		doMET = iConfig.getParameter<bool>("doMET");
#	Line 59 | Line 64 | NtupleWriter::NtupleWriter(const edm::Pa
64		doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
65		doPV = iConfig.getParameter<bool>("doPV");
66		doTopJets = iConfig.getParameter<bool>("doTopJets");
67	+	doTopJetsConstituents = iConfig.getParameter<bool>("doTopJetsConstituents");
68		doTrigger = iConfig.getParameter<bool>("doTrigger");
69	+	SVComputer_  = iConfig.getUntrackedParameter<edm::InputTag>("svComputer",edm::InputTag("combinedSecondaryVertex"));
70	+	doTagInfos = iConfig.getUntrackedParameter<bool>("doTagInfos",false);
71	+	storePFsAroundLeptons = iConfig.getUntrackedParameter<bool>("storePFsAroundLeptons",false);
72
73		// initialization of tree variables
74
#	Line 108 | Line 117 | NtupleWriter::NtupleWriter(const edm::Pa
117		tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
118		}
119		}
120	+	if(doGenJets){
121	+	genjet_sources = iConfig.getParameter<std::vector<std::string> >("genjet_sources");
122	+	genjet_ptmin = iConfig.getParameter<double> ("genjet_ptmin");
123	+	genjet_etamax = iConfig.getParameter<double> ("genjet_etamax");
124	+	for(size_t j=0; j< genjet_sources.size(); ++j){
125	+	tr->Branch( genjet_sources[j].c_str(), "std::vector<Particle>", &genjets[j]);
126	+	}
127	+	}
128		if(doTopJets){
129		topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
130		topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
#	Line 116 | Line 133 | NtupleWriter::NtupleWriter(const edm::Pa
133		tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
134		}
135		}
136	+	if(doTopJetsConstituents){
137	+	topjet_constituents_sources = iConfig.getParameter<std::vector<std::string> >("topjet_constituents_sources");
138	+	tr->Branch( "PFParticles", "std::vector<PFParticle>", &pfparticles);
139	+	}
140		if(doGenTopJets){
141		gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
142		gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
143		gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
144		for(size_t j=0; j< gentopjet_sources.size(); ++j){
145	<	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<TopJet>", &gentopjets[j]);
145	>	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<GenTopJet>", &gentopjets[j]);
146		}
147		}
148		if(doPhotons){
#	Line 143 | Line 164 | NtupleWriter::NtupleWriter(const edm::Pa
164		}
165		}
166		if(doGenInfo){
167	+	genparticle_source= iConfig.getParameter<edm::InputTag>("genparticle_source");
168		tr->Branch("genInfo","GenInfo",&genInfo);
169		tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
170		}
#	Line 154 | Line 176 | NtupleWriter::NtupleWriter(const edm::Pa
176		//tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
177		//tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
178		}
179	+	if (storePFsAroundLeptons){
180	+	pf_around_leptons_source = iConfig.getParameter<std::string>("pf_around_leptons_source");
181	+	}
182		newrun = true;
183		}
184
#	Line 292 | Line 317 | NtupleWriter::analyze(const edm::Event&
317		}
318
319		edm::Handle<reco::GenParticleCollection> genPartColl;
320	<	iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
320	>	iEvent.getByLabel(genparticle_source, genPartColl);
321		int index=-1;
322		for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
323		index++;
324
325	<	//write out only top quarks and status 3 particles (works fine only for MadGraph)
326	<	if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
325	>	//write out only top quarks,final state leptons and status 3 particles (works fine only for MadGraph)
326	>	bool islepton = iter->status()==1 && abs(iter->pdgId())>=11 && abs(iter->pdgId())<=16 ;
327	>	if(abs(iter->pdgId())==6 || iter->status()==3 || islepton ||  doAllGenParticles){
328		GenParticle genp;
329		genp.set_charge(iter->charge());
330		genp.set_pt(iter->p4().pt());
#	Line 328 | Line 354 | NtupleWriter::analyze(const edm::Event&
354		}
355		}
356
357	+	// ------------- full PF collection -------------
358	+
359	+
360		// ------------- electrons -------------
361		if(doElectrons){
362
#	Line 383 | Line 412 | NtupleWriter::analyze(const edm::Event&
412		ele.set_EcalEnergy(pat_ele.ecalEnergy());
413		ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
414		ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
415	+	float AEff03 = 0.00;
416	+	if(isRealData){
417	+	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2011);
418	+	}else{
419	+	AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAFall11MC);
420	+	}
421	+	ele.set_AEff(AEff03);
422
423		eles[j].push_back(ele);
424	+
425	+	if (storePFsAroundLeptons){
426	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
427	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
428	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
429	+	StorePFCandsInCone(&ele, pf_cands, 0.4);
430	+	}
431	+
432		}
433		}
434		}
#	Line 480 | Line 524 | NtupleWriter::analyze(const edm::Event&
524		}
525
526		mus[j].push_back(mu);
527	+
528	+	if (storePFsAroundLeptons){
529	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
530	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
531	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
532	+	StorePFCandsInCone(&mu, pf_cands, 0.4);
533	+	}
534	+
535		}
536		}
537		}
#	Line 506 | Line 558 | NtupleWriter::analyze(const edm::Event&
558		tau.set_phi( pat_tau.phi());
559		tau.set_energy( pat_tau.energy());
560		tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
561	<	tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
561	>	//tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
562		tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
563		tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
564		tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
565	<	tau.set_againstElectronLoose  ( pat_tau.tauID("againstElectronLoose")>0.5);
566	<	tau.set_againstElectronMedium ( pat_tau.tauID("againstElectronMedium")>0.5);
567	<	tau.set_againstElectronTight ( pat_tau.tauID("againstElectronTight")>0.5);
568	<	tau.set_againstElectronMVA  ( pat_tau.tauID("againstElectronMVA")>0.5);
569	<	tau.set_againstMuonLoose ( pat_tau.tauID("againstMuonLoose")>0.5);
570	<	tau.set_againstMuonMedium ( pat_tau.tauID("againstMuonMedium")>0.5);
571	<	tau.set_againstMuonTight ( pat_tau.tauID("againstMuonTight")>0.5);
565	>	tau.set_byLooseIsolationMVA( pat_tau.tauID("byLooseIsolationMVA")>0.5);
566	>	tau.set_byMediumIsolationMVA( pat_tau.tauID("byMediumIsolationMVA")>0.5);
567	>	tau.set_byTightIsolationMVA( pat_tau.tauID("byTightIsolationMVA")>0.5);
568	>	tau.set_byLooseIsolationMVA2( pat_tau.tauID("byLooseIsolationMVA2")>0.5);
569	>	tau.set_byMediumIsolationMVA2( pat_tau.tauID("byMediumIsolationMVA2")>0.5);
570	>	tau.set_byTightIsolationMVA2( pat_tau.tauID("byTightIsolationMVA2")>0.5);
571	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr3Hits(  pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits")>0.5);
572	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits")>0.5);
573	>	tau.set_byTightCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr3Hits")>0.5);
574	>	tau.set_againstElectronLooseMVA3  ( pat_tau.tauID("againstElectronLooseMVA3")>0.5);
575	>	tau.set_againstElectronMediumMVA3 ( pat_tau.tauID("againstElectronMediumMVA3")>0.5);
576	>	tau.set_againstElectronTightMVA3 ( pat_tau.tauID("againstElectronTightMVA3")>0.5);
577	>	tau.set_againstElectronVTightMVA3 ( pat_tau.tauID("againstElectronVTightMVA3")>0.5);
578	>	tau.set_againstMuonLoose2 ( pat_tau.tauID("againstMuonLoose2")>0.5);
579	>	tau.set_againstMuonMedium2 ( pat_tau.tauID("againstMuonMedium2")>0.5);
580	>	tau.set_againstMuonTight2 ( pat_tau.tauID("againstMuonTight2")>0.5);
581	>	tau.set_byIsolationMVAraw(  pat_tau.tauID("byIsolationMVAraw"));
582	>	tau.set_byIsolationMVA2raw(  pat_tau.tauID("byIsolationMVA2raw"));
583	>	tau.set_decayMode( pat_tau.decayMode() );
584	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw"));
585	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw3Hits( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));
586
587	+	//       std::cout << pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") << std::endl;
588	+
589		//       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
590		//       if(!leadPFCand.isNull()){
591		//         tau.set_leadPFCand_px ( leadPFCand->px());
#	Line 530 | Line 598 | NtupleWriter::analyze(const edm::Event&
598		//         tau.set_leadPFCand_pz ( 0);
599		//       }
600		taus[j].push_back(tau);
601	+
602	+	// store isolation info only for identified taus
603	+	if (pat_tau.tauID("decayModeFinding")>0.5){
604	+	bool storepfs = (pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5) || (pat_tau.tauID("byLooseIsolationMVA")>0.5);
605	+	if (storePFsAroundLeptons && storepfs){
606	+	edm::Handle< std::vector<reco::PFCandidate> > pfcand_handle;
607	+	iEvent.getByLabel(pf_around_leptons_source, pfcand_handle);  // default: "pfNoPileUpPFlow"
608	+	const std::vector<reco::PFCandidate>& pf_cands = *(pfcand_handle.product());
609	+	StorePFCandsInCone(&tau, pf_cands, 0.4);
610	+	}
611	+	}
612	+	}
613	+	}
614	+	}
615	+
616	+	//-------------- gen jets -------------
617	+
618	+	if(doGenJets){
619	+	for(size_t j=0; j< genjet_sources.size(); ++j){
620	+
621	+	genjets[j].clear();
622	+
623	+	edm::Handle< std::vector<reco::GenJet> > genjet_handle;
624	+	iEvent.getByLabel(genjet_sources[j], genjet_handle);
625	+	const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
626	+
627	+	for (unsigned int i = 0; i < gen_jets.size(); ++i) {
628	+
629	+	const reco::GenJet* gen_jet = &gen_jets[i];
630	+	if(gen_jet->pt() < genjet_ptmin) continue;
631	+	if(fabs(gen_jet->eta()) > genjet_etamax) continue;
632	+
633	+	Particle jet;
634	+	jet.set_charge(gen_jet->charge());
635	+	jet.set_pt(gen_jet->pt());
636	+	jet.set_eta(gen_jet->eta());
637	+	jet.set_phi(gen_jet->phi());
638	+	jet.set_energy(gen_jet->energy());
639	+
640	+	// recalculate the jet charge
641	+	int jet_charge = 0;
642	+	std::vector<const reco::GenParticle * > jetgenps = gen_jet->getGenConstituents();
643	+	for(unsigned int l = 0; l<jetgenps.size(); ++l){
644	+	jet_charge +=  jetgenps[l]->charge();
645	+	}
646	+	jet.set_charge(jet_charge);
647	+
648	+	genjets[j].push_back(jet);
649		}
650		}
651		}
#	Line 560 | Line 676 | NtupleWriter::analyze(const edm::Event&
676		jet.set_eta(pat_jet.eta());
677		jet.set_phi(pat_jet.phi());
678		jet.set_energy(pat_jet.energy());
679	+	jet.set_flavor(pat_jet.partonFlavour());
680		jet.set_numberOfDaughters (pat_jet.numberOfDaughters());
681		const reco::TrackRefVector&  jettracks = pat_jet.associatedTracks();
682		jet.set_nTracks ( jettracks.size());
683		jet.set_jetArea(pat_jet.jetArea());
567	–	jet.set_pileup(pat_jet.pileup());
684		if(pat_jet.isPFJet()){
685		jet.set_neutralEmEnergyFraction (pat_jet.neutralEmEnergyFraction());
686		jet.set_neutralHadronEnergyFraction (pat_jet.neutralHadronEnergyFraction());
#	Line 578 | Line 694 | NtupleWriter::analyze(const edm::Event&
694		jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
695		jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
696		}
697	<	if(doJECUncertainty){
582	<	jecUnc->setJetEta(pat_jet.eta());
583	<	jecUnc->setJetPt(pat_jet.pt());
584	<	jet.set_JEC_uncertainty(jecUnc->getUncertainty(true));
585	<	}
697	>
698		jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
699
700		jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
#	Line 592 | Line 704 | NtupleWriter::analyze(const edm::Event&
704		jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
705		jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
706
707	+
708		const reco::GenJet *genj = pat_jet.genJet();
709		if(genj){
710	<	jet.set_genjet_pt(genj->pt());
711	<	jet.set_genjet_eta(genj->eta());
712	<	jet.set_genjet_phi(genj->phi());
713	<	jet.set_genjet_energy(genj->energy());
601	<	if(doAllGenParticles){
602	<	std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
603	<	for(unsigned int l = 0; l<jetgenps.size(); ++l){
604	<	for(unsigned int k=0; k< genps.size(); ++k){
605	<	if(jetgenps[l]->pt() == genps[k].pt() && jetgenps[l]->pdgId() == genps[k].pdgId()){
606	<	jet.add_genparticles_index(genps[k].index());
607	<	break;
608	<	}
609	<	}
710	>
711	>	for(unsigned int k=0; k<genjets->size(); ++k){
712	>	if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta()){
713	>	jet.set_genjet_index(k);
714		}
611	–	if(jet.genparticles_indices().size()!= jetgenps.size())
612	–	std::cout << "WARNING: Found only " << jet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
715		}
716	+	//         if( jet.genjet_index()<0){
717	+	//           std::cout<< "genjet not found for " << genj->pt() << "  " << genj->eta() << std::endl;
718	+	//         }
719	+
720		}
721	+
722		jets[j].push_back(jet);
723		}
724		}
725		}
726
727	+
728		// ------------- top jets -------------
729		if(doTopJets){
730	+
731		for(size_t j=0; j< topjet_sources.size(); ++j){
732
733		topjets[j].clear();
#	Line 639 | Line 748 | NtupleWriter::analyze(const edm::Event&
748		topjet.set_eta(pat_topjet.eta());
749		topjet.set_phi(pat_topjet.phi());
750		topjet.set_energy(pat_topjet.energy());
751	+	topjet.set_flavor(pat_topjet.partonFlavour());
752		topjet.set_numberOfDaughters(pat_topjet.numberOfDaughters());
753		const reco::TrackRefVector&  topjettracks = pat_topjet.associatedTracks();
754		topjet.set_nTracks( topjettracks.size());
755		topjet.set_jetArea( pat_topjet.jetArea());
756	<	topjet.set_pileup( pat_topjet.pileup());
647	<	//       topjet.set_neutralEmEnergyFraction(pat_topjet.neutralEmEnergyFraction());
648	<	//       topjet.set_neutralHadronEnergyFraction(pat_topjet.neutralHadronEnergyFraction());
649	<	//       topjet.set_chargedEmEnergyFraction(pat_topjet.chargedEmEnergyFraction());
650	<	//       topjet.set_chargedHadronEnergyFraction(pat_topjet.chargedHadronEnergyFraction());
651	<	//       topjet.set_muonEnergyFraction(pat_topjet.muonEnergyFraction());
652	<	//       topjet.set_photonEnergyFraction(pat_topjet.photonEnergyFraction());
653	<	//       topjet.set_chargedMultiplicity(pat_topjet.chargedMultiplicity());
654	<	//       topjet.set_neutralMultiplicity(pat_topjet.neutralMultiplicity());
655	<	//       topjet.set_muonMultiplicity(pat_topjet.muonMultiplicity());
656	<	//       topjet.set_electronMultiplicity(pat_topjet.electronMultiplicity());
657	<	//       topjet.set_photonMultiplicity(pat_topjet.photonMultiplicity());
658	<	if(doJECUncertainty){
659	<	jecUnc->setJetEta(pat_topjet.eta());
660	<	jecUnc->setJetPt(pat_topjet.pt());
661	<	topjet.set_JEC_uncertainty( jecUnc->getUncertainty(true));
662	<	}
756	>
757		topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
758
759		topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
#	Line 692 | Line 786 | NtupleWriter::analyze(const edm::Event&
786		}
787		*/
788
789	+	// add constituents to the jet, if requested
790	+	if (doTopJetsConstituents){
791	+
792	+	if (topjet_constituents_sources.size()>j){ //only add constituents if they are defined
793	+
794	+	edm::Handle<pat::JetCollection> pat_topjets_with_cands;
795	+	iEvent.getByLabel(topjet_constituents_sources[j], pat_topjets_with_cands);
796	+	pat::Jet* pat_topjet_wc = NULL;
797	+
798	+	for (unsigned int it = 0; it < pat_topjets_with_cands->size(); it++) {
799	+	const pat::Jet* cand =  dynamic_cast<pat::Jet const *>(&pat_topjets_with_cands->at(it));
800	+	double dphi = cand->phi() - pat_topjet.phi();
801	+	if (dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
802	+	if (dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
803	+	if (fabs(dphi)<0.5 && fabs(cand->eta()-pat_topjet.eta())<0.5){ // be generous: filtering, pruning... can change jet axis
804	+	pat_topjet_wc = const_cast<pat::Jet*>(cand);
805	+	break;
806	+	}
807	+	}
808	+
809	+	if (pat_topjet_wc){
810	+	StoreJetConstituents(pat_topjet_wc, &topjet);
811	+
812	+	// now run substructure information
813	+	JetProps substructure(&topjet);
814	+	substructure.set_pf_cands(&pfparticles);
815	+	double tau1 = substructure.GetNsubjettiness(1, Njettiness::onepass_kt_axes, 1., 2.0);
816	+	double tau2 = substructure.GetNsubjettiness(2, Njettiness::onepass_kt_axes, 1., 2.0);
817	+	double tau3 = substructure.GetNsubjettiness(3, Njettiness::onepass_kt_axes, 1., 2.0);
818	+	double qjets = substructure.GetQjetVolatility(iEvent.id().event(), 2.0);
819	+	topjet.set_tau1(tau1);
820	+	topjet.set_tau2(tau2);
821	+	topjet.set_tau3(tau3);
822	+	topjet.set_qjets_volatility(qjets);
823	+
824	+	}
825	+	}
826	+	}
827	+
828	+
829	+
830		for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
831		Particle subjet_v4;
832	<	subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
833	<	subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
834	<	subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
835	<	subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
836	<	topjet.add_subjet(subjet_v4);
832	>
833	>	reco::Candidate const * subjetd =  pat_topjet.daughter(k);
834	>	pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
835	>	if(patsubjetd)
836	>	{
837	>	subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
838	>	subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
839	>	subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
840	>	subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
841	>	topjet.add_subjet(subjet_v4);
842	>	//btag info
843	>	topjet.add_subFlavour(patsubjetd->partonFlavour());
844	>	topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
845	>	if (doTagInfos)
846	>	{
847	>	//ip tag info
848	>	reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
849	>	topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
850	>	topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
851	>	topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
852	>	topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
853	>	topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
854	>	topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
855	>	topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
856	>	topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
857	>	topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
858	>	topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
859	>	topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
860	>	topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));
861	>	topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
862	>	topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
863	>	topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
864	>	topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
865	>	topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
866	>	topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
867	>	topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
868	>	topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
869	>	topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
870	>	//sv tag info
871	>	reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
872	>	topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
873	>	topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
874	>	topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
875	>	topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
876	>	topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
877	>	topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
878	>	topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
879	>	std::vector<TLorentzVector> vp4; vp4.clear();
880	>	std::vector<float> vchi2; vchi2.clear();
881	>	std::vector<float> vndof; vndof.clear();
882	>	std::vector<float> vchi2ndof; vchi2ndof.clear();
883	>	std::vector<float> vtrsize; vtrsize.clear();
884	>	for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
885	>	{
886	>	ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
887	>	vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
888	>	vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());
889	>	vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());
890	>	vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());
891	>	vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());
892	>	}
893	>	topjet.add_subSecondaryVertex(vp4);
894	>	topjet.add_subVertexChi2(vchi2);
895	>	topjet.add_subVertexNdof(vndof);
896	>	topjet.add_subVertexNormalizedChi2(vchi2ndof);
897	>	topjet.add_subVertexTracksSize(vtrsize);
898	>	//try computer
899	>	edm::ESHandle<JetTagComputer> computerHandle;
900	>	iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
901	>	const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
902	>	if(computer)
903	>	{
904	>	computer->passEventSetup(iSetup);
905	>	std::vector<const reco::BaseTagInfo*>  baseTagInfos;
906	>	baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
907	>	baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );
908	>	JetTagComputer::TagInfoHelper helper(baseTagInfos);
909	>	reco::TaggingVariableList vars = computer->taggingVariables(helper);
910	>	topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
911	>	topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
912	>	topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
913	>	topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
914	>	}
915	>	}
916	>	}
917	>	else
918	>	{
919	>	//filling only standard information in case the subjet has not been pat-tified during the pattuples production
920	>	subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
921	>	subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
922	>	subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
923	>	subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
924	>	topjet.add_subjet(subjet_v4);
925	>	}
926	>
927	>
928		}
929		topjets[j].push_back(topjet);
930		}
931		}
932		}
933
934	<
934	>
935		// ------------- generator top jets -------------
936		if(doGenTopJets){
937		for(size_t j=0; j< gentopjet_sources.size(); ++j){
938
939		gentopjets[j].clear();
940	<
940	>
941		edm::Handle<reco::BasicJetCollection> reco_gentopjets;
716	–	//edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
942		iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
943
944		for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
#	Line 722 | Line 947 | NtupleWriter::analyze(const edm::Event&
947		if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
948		if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
949
950	<	TopJet gentopjet;
950	>	GenTopJet gentopjet;
951		gentopjet.set_charge(reco_gentopjet.charge());
952		gentopjet.set_pt(reco_gentopjet.pt());
953		gentopjet.set_eta(reco_gentopjet.eta());
954		gentopjet.set_phi(reco_gentopjet.phi());
955		gentopjet.set_energy(reco_gentopjet.energy());
731	–	gentopjet.set_numberOfDaughters(reco_gentopjet.numberOfDaughters());
956
957		for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
958		Particle subjet_v4;
#	Line 743 | Line 967 | NtupleWriter::analyze(const edm::Event&
967		}
968		}
969
970	+
971		// ------------- photons -------------
972		if(doPhotons){
973		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 847 | Line 1072 | NtupleWriter::analyze(const edm::Event&
1072		tr->Fill();
1073		if(doLumiInfo)
1074		previouslumiblockwasfilled=true;
1075	+
1076	+	// clean up
1077	+	if(doTopJetsConstituents){
1078	+	pfparticles.clear();
1079	+	}
1080	+
1081	+	// need to do a check if necessary
1082	+	pfparticles.clear();
1083	+
1084		}
1085
1086
#	Line 880 | Line 1114 | NtupleWriter::beginRun(edm::Run const& i
1114		newrun=true;
1115		}
1116
883	–	if(doJets || doTopJets){
884	–	if(doJECUncertainty){
885	–	edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
886	–	iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
887	–	JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
888	–	jecUnc = new JetCorrectionUncertainty(JetCorPar);
889	–	}
890	–	}
1117		}
1118
1119		// ------------ method called when ending the processing of a run  ------------
#	Line 940 | Line 1166 | NtupleWriter::fillDescriptions(edm::Conf
1166		descriptions.addDefault(desc);
1167		}
1168
1169	+	// ------------ method fills constituents of the pat_jet into the Ntuple and stores a reference
1170	+	// ------------ to those in the topjet
1171	+	// ------------ it is checked if the constituents have been stored already
1172	+	void
1173	+	NtupleWriter::StoreJetConstituents(pat::Jet* pat_jet, Jet* jet)
1174	+	{
1175	+	// checks if the pf cand has already been stored, only stores so far missing ones
1176	+	// PF candidates, then stores a reference to the pf candidate in the jet
1177	+	// also calculates the jet charge and sets it
1178	+
1179	+
1180	+	const std::vector<reco::PFCandidatePtr> jconstits = pat_jet->getPFConstituents();
1181	+
1182	+	unsigned int NstoredPFs = pfparticles.size();
1183	+
1184	+	// loop over all jet constituents and store them
1185	+	for (unsigned int i=0; i<jconstits.size(); ++i){
1186	+
1187	+	PFParticle part;
1188	+	const reco::PFCandidate* pf = jconstits[i].get();
1189	+
1190	+	// check if it has already been stored, omit if true
1191	+	int is_already_in_list = -1;
1192	+	for (unsigned int j=0; j<NstoredPFs; ++j){
1193	+	PFParticle spf = pfparticles[j];
1194	+	double r2 = pow(pf->eta()-spf.eta(),2) + pow(pf->phi()-spf.phi(),2);
1195	+	double dpt = fabs( pf->pt() - spf.pt() );
1196	+	if (r2<1e-10 && dpt<1e-10){
1197	+	is_already_in_list = j;
1198	+	break;
1199	+	}
1200	+	}
1201	+
1202	+	if (is_already_in_list>-1){
1203	+	continue;
1204	+	}
1205	+
1206	+
1207	+	part.set_pt(pf->pt());
1208	+	part.set_eta(pf->eta());
1209	+	part.set_phi(pf->phi());
1210	+	part.set_energy(pf->energy());
1211	+	part.set_charge(pf->charge());
1212	+
1213	+	part.set_ecal_en(pf->ecalEnergy());
1214	+	part.set_hcal_en(pf->hcalEnergy());
1215	+	reco::TrackRef trackref = pf->trackRef();
1216	+	if (!trackref.isNull()){
1217	+	part.set_track_mom(trackref->p());
1218	+	}
1219	+
1220	+	PFParticle::EParticleID id = PFParticle::eX;
1221	+	switch ( pf->translatePdgIdToType(pf->pdgId()) ){
1222	+	case reco::PFCandidate::X : id = PFParticle::eX; break;
1223	+	case reco::PFCandidate::h : id = PFParticle::eH; break;
1224	+	case reco::PFCandidate::e : id = PFParticle::eE; break;
1225	+	case reco::PFCandidate::mu : id = PFParticle::eMu; break;
1226	+	case reco::PFCandidate::gamma : id = PFParticle::eGamma; break;
1227	+	case reco::PFCandidate::h0 : id = PFParticle::eH0; break;
1228	+	case reco::PFCandidate::h_HF : id = PFParticle::eH_HF; break;
1229	+	case reco::PFCandidate::egamma_HF : id = PFParticle::eEgamma_HF; break;
1230	+	}
1231	+	part.set_particleID(id);
1232	+
1233	+	pfparticles.push_back(part);
1234	+
1235	+	// add a reference to the particle
1236	+	jet->add_pfconstituents_index(pfparticles.size()-1);
1237	+
1238	+	}
1239	+
1240	+	if(pat_jet->isPFJet()){
1241	+	jet->set_charge(pat_jet->charge());
1242	+	jet->set_neutralEmEnergyFraction(pat_jet->neutralEmEnergyFraction());
1243	+	jet->set_neutralHadronEnergyFraction (pat_jet->neutralHadronEnergyFraction());
1244	+	jet->set_chargedEmEnergyFraction (pat_jet->chargedEmEnergyFraction());
1245	+	jet->set_chargedHadronEnergyFraction (pat_jet->chargedHadronEnergyFraction());
1246	+	jet->set_muonEnergyFraction (pat_jet->muonEnergyFraction());
1247	+	jet->set_photonEnergyFraction (pat_jet->photonEnergyFraction());
1248	+	jet->set_chargedMultiplicity (pat_jet->chargedMultiplicity());
1249	+	jet->set_neutralMultiplicity (pat_jet->neutralMultiplicity());
1250	+	jet->set_muonMultiplicity (pat_jet->muonMultiplicity());
1251	+	jet->set_electronMultiplicity (pat_jet->electronMultiplicity());
1252	+	jet->set_photonMultiplicity (pat_jet->photonMultiplicity());
1253	+	}
1254	+
1255	+	return;
1256	+
1257	+	}
1258	+
1259	+	// ------------ method fills PF candidates in a cone of radius R around
1260	+	// ------------ a given particle (lepton, most likely)
1261	+	// ------------ it is checked if the constituents have been stored already
1262	+	void
1263	+	NtupleWriter::StorePFCandsInCone(Particle* inpart, const std::vector<reco::PFCandidate>& pf_cands, double R0)
1264	+	{
1265	+	// checks if the pf cand has already been stored, only stores so far missing ones
1266	+
1267	+	//cout << "\nStorePFCandsInCone: R0 = " << R0 << endl;
1268	+	//cout << "found " << pf_cands.size() << " PF candidates in the event" << endl;
1269	+	//cout << "Input inparticle: pt = " << inpart->pt() << " eta = " << inpart->eta() << " phi = " << inpart->phi() << endl;
1270	+
1271	+	unsigned int NstoredPFs = pfparticles.size();
1272	+	//cout << "got already " <<  NstoredPFs << " which need to be checked to avoid double counting" << endl;
1273	+
1274	+	// loop over all PF candidates and store the ones in a cone with radius R0
1275	+	for (unsigned int i=0; i<pf_cands.size(); ++i){
1276	+
1277	+	reco::PFCandidate pf = pf_cands[i];
1278	+
1279	+	// calculate distance in eta/phi
1280	+	double dphi = pf.phi() - inpart->phi();
1281	+	if (dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
1282	+	if (dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
1283	+	double dr2 = dphi*dphi + pow(pf.eta() - inpart->eta(),2);
1284	+
1285	+	// check if PF candidate is in cone around particle
1286	+	if (sqrt(dr2)>R0) continue;
1287	+
1288	+	//cout << "found particle with distance " << dr2 << " at position " << i << endl;
1289	+	//cout << "PF: pt = " << pf.pt() << " eta = " << pf.eta() << " phi = " << pf.phi() << endl;
1290	+
1291	+	// check if it's the same as the input particle
1292	+	double dpt = fabs( inpart->pt() - pf.pt() );
1293	+	if (dr2<1e-10 && dpt<1e-10){
1294	+	//cout << "same particle, don't store" << endl;
1295	+	continue;
1296	+	}
1297	+
1298	+	// check if it has already been stored, omit if true
1299	+	int is_already_in_list = -1;
1300	+	for (unsigned int j=0; j<NstoredPFs; ++j){
1301	+	PFParticle spf = pfparticles[j];
1302	+	double r2 = pow(pf.eta()-spf.eta(),2) + pow(pf.phi()-spf.phi(),2);
1303	+	double dpt = fabs( pf.pt() - spf.pt() );
1304	+	if (r2<1e-10 && dpt<1e-10){
1305	+	is_already_in_list = j;
1306	+	break;
1307	+	}
1308	+	}
1309	+
1310	+	if (is_already_in_list>-1){
1311	+	//cout << "is already in list, continue" << endl;
1312	+	continue;
1313	+	}
1314	+
1315	+
1316	+	PFParticle part;
1317	+	part.set_pt(pf.pt());
1318	+	part.set_eta(pf.eta());
1319	+	part.set_phi(pf.phi());
1320	+	part.set_energy(pf.energy());
1321	+	part.set_charge(pf.charge());
1322	+
1323	+	part.set_ecal_en(pf.ecalEnergy());
1324	+	part.set_hcal_en(pf.hcalEnergy());
1325	+	reco::TrackRef trackref = pf.trackRef();
1326	+	if (!trackref.isNull()){
1327	+	part.set_track_mom(trackref->p());
1328	+	}
1329	+
1330	+	PFParticle::EParticleID id = PFParticle::eX;
1331	+	switch ( pf.translatePdgIdToType(pf.pdgId()) ){
1332	+	case reco::PFCandidate::X : id = PFParticle::eX; break;
1333	+	case reco::PFCandidate::h : id = PFParticle::eH; break;
1334	+	case reco::PFCandidate::e : id = PFParticle::eE; break;
1335	+	case reco::PFCandidate::mu : id = PFParticle::eMu; break;
1336	+	case reco::PFCandidate::gamma : id = PFParticle::eGamma; break;
1337	+	case reco::PFCandidate::h0 : id = PFParticle::eH0; break;
1338	+	case reco::PFCandidate::h_HF : id = PFParticle::eH_HF; break;
1339	+	case reco::PFCandidate::egamma_HF : id = PFParticle::eEgamma_HF; break;
1340	+	}
1341	+	part.set_particleID(id);
1342	+
1343	+	pfparticles.push_back(part);
1344	+
1345	+	}
1346	+
1347	+	return;
1348	+
1349	+	}
1350	+
1351		//define this as a plug-in
1352		DEFINE_FWK_MODULE(NtupleWriter);

Diff Legend

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