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.22 by peiffer, Wed Jul 25 09:56:57 2012 UTC vs.
Revision 1.27 by rkogler, Wed Jun 19 13:22:06 2013 UTC

#	Line 17 | Line 17
17		//
18		//
19
20	<	#include "UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
21	<
22	<
20	>	#include "UHHAnalysis/NtupleWriter/interface/NtupleWriter.h"
21	>	#include "UHHAnalysis/NtupleWriter/interface/JetProps.h"
22	>	#include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputer.h"
23	>	#include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputerWrapper.h"
24	>	#include "RecoBTau/JetTagComputer/interface/JetTagComputer.h"
25	>	#include "RecoBTau/JetTagComputer/interface/JetTagComputerRecord.h"
26	>	#include "RecoBTag/SecondaryVertex/interface/CombinedSVComputer.h"
27	>	#include "DataFormats/TrackReco/interface/Track.h"
28
29		//
30		// constants, enums and typedefs
#	Line 50 | Line 55 | NtupleWriter::NtupleWriter(const edm::Pa
55		doMuons = iConfig.getParameter<bool>("doMuons");
56		doTaus = iConfig.getParameter<bool>("doTaus");
57		doJets = iConfig.getParameter<bool>("doJets");
58	<	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
72		// initialization of tree variables
73
#	Line 108 | Line 116 | NtupleWriter::NtupleWriter(const edm::Pa
116		tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
117		}
118		}
119	+	if(doGenJets){
120	+	genjet_sources = iConfig.getParameter<std::vector<std::string> >("genjet_sources");
121	+	genjet_ptmin = iConfig.getParameter<double> ("genjet_ptmin");
122	+	genjet_etamax = iConfig.getParameter<double> ("genjet_etamax");
123	+	for(size_t j=0; j< genjet_sources.size(); ++j){
124	+	tr->Branch( genjet_sources[j].c_str(), "std::vector<Particle>", &genjets[j]);
125	+	}
126	+	}
127		if(doTopJets){
128		topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
129		topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
#	Line 116 | Line 132 | NtupleWriter::NtupleWriter(const edm::Pa
132		tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
133		}
134		}
135	+	if(doTopJetsConstituents){
136	+	topjet_constituents_sources = iConfig.getParameter<std::vector<std::string> >("topjet_constituents_sources");
137	+	tr->Branch( "PFParticles", "std::vector<PFParticle>", &pfparticles);
138	+	}
139		if(doGenTopJets){
140		gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
141		gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
142		gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
143		for(size_t j=0; j< gentopjet_sources.size(); ++j){
144	<	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<TopJet>", &gentopjets[j]);
144	>	tr->Branch( gentopjet_sources[j].c_str(), "std::vector<GenTopJet>", &gentopjets[j]);
145		}
146		}
147		if(doPhotons){
#	Line 298 | Line 318 | NtupleWriter::analyze(const edm::Event&
318		for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
319		index++;
320
321	<	//write out only top quarks and status 3 particles (works fine only for MadGraph)
322	<	if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
321	>	//write out only top quarks,final state leptons and status 3 particles (works fine only for MadGraph)
322	>	bool islepton = iter->status()==1 && abs(iter->pdgId())>=11 && abs(iter->pdgId())<=16 ;
323	>	if(abs(iter->pdgId())==6 || iter->status()==3 || islepton ||  doAllGenParticles){
324		GenParticle genp;
325		genp.set_charge(iter->charge());
326		genp.set_pt(iter->p4().pt());
#	Line 514 | Line 535 | NtupleWriter::analyze(const edm::Event&
535		tau.set_phi( pat_tau.phi());
536		tau.set_energy( pat_tau.energy());
537		tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
538	<	tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
538	>	//tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
539		tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
540		tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
541		tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
542	<	tau.set_againstElectronLoose  ( pat_tau.tauID("againstElectronLoose")>0.5);
543	<	tau.set_againstElectronMedium ( pat_tau.tauID("againstElectronMedium")>0.5);
544	<	tau.set_againstElectronTight ( pat_tau.tauID("againstElectronTight")>0.5);
545	<	tau.set_againstElectronMVA  ( pat_tau.tauID("againstElectronMVA")>0.5);
546	<	tau.set_againstMuonLoose ( pat_tau.tauID("againstMuonLoose")>0.5);
547	<	tau.set_againstMuonMedium ( pat_tau.tauID("againstMuonMedium")>0.5);
548	<	tau.set_againstMuonTight ( pat_tau.tauID("againstMuonTight")>0.5);
542	>	tau.set_byLooseIsolationMVA( pat_tau.tauID("byLooseIsolationMVA")>0.5);
543	>	tau.set_byMediumIsolationMVA( pat_tau.tauID("byMediumIsolationMVA")>0.5);
544	>	tau.set_byTightIsolationMVA( pat_tau.tauID("byTightIsolationMVA")>0.5);
545	>	tau.set_byLooseIsolationMVA2( pat_tau.tauID("byLooseIsolationMVA2")>0.5);
546	>	tau.set_byMediumIsolationMVA2( pat_tau.tauID("byMediumIsolationMVA2")>0.5);
547	>	tau.set_byTightIsolationMVA2( pat_tau.tauID("byTightIsolationMVA2")>0.5);
548	>	tau.set_byLooseCombinedIsolationDeltaBetaCorr3Hits(  pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits")>0.5);
549	>	tau.set_byMediumCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits")>0.5);
550	>	tau.set_byTightCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr3Hits")>0.5);
551	>	tau.set_againstElectronLooseMVA3  ( pat_tau.tauID("againstElectronLooseMVA3")>0.5);
552	>	tau.set_againstElectronMediumMVA3 ( pat_tau.tauID("againstElectronMediumMVA3")>0.5);
553	>	tau.set_againstElectronTightMVA3 ( pat_tau.tauID("againstElectronTightMVA3")>0.5);
554	>	tau.set_againstElectronVTightMVA3 ( pat_tau.tauID("againstElectronVTightMVA3")>0.5);
555	>	tau.set_againstMuonLoose2 ( pat_tau.tauID("againstMuonLoose2")>0.5);
556	>	tau.set_againstMuonMedium2 ( pat_tau.tauID("againstMuonMedium2")>0.5);
557	>	tau.set_againstMuonTight2 ( pat_tau.tauID("againstMuonTight2")>0.5);
558	>	tau.set_byIsolationMVAraw(  pat_tau.tauID("byIsolationMVAraw"));
559	>	tau.set_byIsolationMVA2raw(  pat_tau.tauID("byIsolationMVA2raw"));
560	>	tau.set_decayMode( pat_tau.decayMode() );
561	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw"));
562	>	tau.set_byCombinedIsolationDeltaBetaCorrRaw3Hits( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));
563
564	+	//       std::cout << pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") << std::endl;
565	+
566		//       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
567		//       if(!leadPFCand.isNull()){
568		//         tau.set_leadPFCand_px ( leadPFCand->px());
#	Line 542 | Line 579 | NtupleWriter::analyze(const edm::Event&
579		}
580		}
581
582	+	//-------------- gen jets -------------
583	+
584	+	if(doGenJets){
585	+	for(size_t j=0; j< genjet_sources.size(); ++j){
586	+
587	+	genjets[j].clear();
588	+
589	+	edm::Handle< std::vector<reco::GenJet> > genjet_handle;
590	+	iEvent.getByLabel(genjet_sources[j], genjet_handle);
591	+	const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
592	+
593	+	for (unsigned int i = 0; i < gen_jets.size(); ++i) {
594	+
595	+	const reco::GenJet* gen_jet = &gen_jets[i];
596	+	if(gen_jet->pt() < genjet_ptmin) continue;
597	+	if(fabs(gen_jet->eta()) > genjet_etamax) continue;
598	+
599	+	Particle jet;
600	+	jet.set_charge(gen_jet->charge());
601	+	jet.set_pt(gen_jet->pt());
602	+	jet.set_eta(gen_jet->eta());
603	+	jet.set_phi(gen_jet->phi());
604	+	jet.set_energy(gen_jet->energy());
605	+
606	+	// recalculate the jet charge
607	+	int jet_charge = 0;
608	+	std::vector<const reco::GenParticle * > jetgenps = gen_jet->getGenConstituents();
609	+	for(unsigned int l = 0; l<jetgenps.size(); ++l){
610	+	jet_charge +=  jetgenps[l]->charge();
611	+	}
612	+	jet.set_charge(jet_charge);
613	+
614	+	genjets[j].push_back(jet);
615	+	}
616	+	}
617	+	}
618	+
619		// ------------- jets -------------
620		if(doJets){
621		for(size_t j=0; j< jet_sources.size(); ++j){
#	Line 586 | Line 660 | NtupleWriter::analyze(const edm::Event&
660		jet.set_electronMultiplicity (pat_jet.electronMultiplicity());
661		jet.set_photonMultiplicity (pat_jet.photonMultiplicity());
662		}
663	<	if(doJECUncertainty){
590	<	jecUnc->setJetEta(pat_jet.eta());
591	<	jecUnc->setJetPt(pat_jet.pt());
592	<	jet.set_JEC_uncertainty(jecUnc->getUncertainty(true));
593	<	}
663	>
664		jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
665
666		jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
#	Line 603 | Line 673 | NtupleWriter::analyze(const edm::Event&
673
674		const reco::GenJet *genj = pat_jet.genJet();
675		if(genj){
676	<	jet.set_genjet_pt(genj->pt());
677	<	jet.set_genjet_eta(genj->eta());
678	<	jet.set_genjet_phi(genj->phi());
679	<	jet.set_genjet_energy(genj->energy());
676	>
677	>	for(unsigned int k=0; k<genjets->size(); ++k){
678	>	if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta()){
679	>	jet.set_genjet_index(k);
680	>	}
681	>	}
682	>	//         if( jet.genjet_index()<0){
683	>	//           std::cout<< "genjet not found for " << genj->pt() << "  " << genj->eta() << std::endl;
684	>	//         }
685	>
686		if(doAllGenParticles){
687		std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
688		for(unsigned int l = 0; l<jetgenps.size(); ++l){
#	Line 628 | Line 704 | NtupleWriter::analyze(const edm::Event&
704		}
705		}
706
707	+
708		// ------------- top jets -------------
709		if(doTopJets){
710	+
711		for(size_t j=0; j< topjet_sources.size(); ++j){
712
713		topjets[j].clear();
#	Line 655 | Line 733 | NtupleWriter::analyze(const edm::Event&
733		topjet.set_nTracks( topjettracks.size());
734		topjet.set_jetArea( pat_topjet.jetArea());
735		topjet.set_pileup( pat_topjet.pileup());
736	<	//       topjet.set_neutralEmEnergyFraction(pat_topjet.neutralEmEnergyFraction());
659	<	//       topjet.set_neutralHadronEnergyFraction(pat_topjet.neutralHadronEnergyFraction());
660	<	//       topjet.set_chargedEmEnergyFraction(pat_topjet.chargedEmEnergyFraction());
661	<	//       topjet.set_chargedHadronEnergyFraction(pat_topjet.chargedHadronEnergyFraction());
662	<	//       topjet.set_muonEnergyFraction(pat_topjet.muonEnergyFraction());
663	<	//       topjet.set_photonEnergyFraction(pat_topjet.photonEnergyFraction());
664	<	//       topjet.set_chargedMultiplicity(pat_topjet.chargedMultiplicity());
665	<	//       topjet.set_neutralMultiplicity(pat_topjet.neutralMultiplicity());
666	<	//       topjet.set_muonMultiplicity(pat_topjet.muonMultiplicity());
667	<	//       topjet.set_electronMultiplicity(pat_topjet.electronMultiplicity());
668	<	//       topjet.set_photonMultiplicity(pat_topjet.photonMultiplicity());
669	<	if(doJECUncertainty){
670	<	jecUnc->setJetEta(pat_topjet.eta());
671	<	jecUnc->setJetPt(pat_topjet.pt());
672	<	topjet.set_JEC_uncertainty( jecUnc->getUncertainty(true));
673	<	}
736	>
737		topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
738
739		topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
#	Line 703 | Line 766 | NtupleWriter::analyze(const edm::Event&
766		}
767		*/
768
769	+	// add constituents to the jet, if requested
770	+	if (doTopJetsConstituents){
771	+
772	+	if (topjet_constituents_sources.size()>j){ //only add constituents if they are defined
773	+
774	+	edm::Handle<pat::JetCollection> pat_topjets_with_cands;
775	+	iEvent.getByLabel(topjet_constituents_sources[j], pat_topjets_with_cands);
776	+	pat::Jet* pat_topjet_wc = NULL;
777	+
778	+	for (unsigned int it = 0; it < pat_topjets_with_cands->size(); it++) {
779	+	const pat::Jet* cand =  dynamic_cast<pat::Jet const *>(&pat_topjets_with_cands->at(it));
780	+	double dphi = cand->phi() - pat_topjet.phi();
781	+	if (dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
782	+	if (dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
783	+	if (fabs(dphi)<0.5 && fabs(cand->eta()-pat_topjet.eta())<0.5){ // be generous: filtering, pruning... can change jet axis
784	+	pat_topjet_wc = const_cast<pat::Jet*>(cand);
785	+	break;
786	+	}
787	+	}
788	+
789	+	if (pat_topjet_wc){
790	+	StoreJetConstituents(pat_topjet_wc, &topjet);
791	+
792	+	// now run substructure information
793	+	JetProps substructure(&topjet);
794	+	substructure.set_pf_cands(&pfparticles);
795	+	double tau1 = substructure.GetNsubjettiness(1, Njettiness::onepass_kt_axes, 1., 2.0);
796	+	double tau2 = substructure.GetNsubjettiness(2, Njettiness::onepass_kt_axes, 1., 2.0);
797	+	double tau3 = substructure.GetNsubjettiness(3, Njettiness::onepass_kt_axes, 1., 2.0);
798	+	double qjets = substructure.GetQjetVolatility(iEvent.id().event(), 2.0);
799	+	topjet.set_tau1(tau1);
800	+	topjet.set_tau2(tau2);
801	+	topjet.set_tau3(tau3);
802	+	topjet.set_qjets_volatility(qjets);
803	+
804	+	}
805	+	}
806	+	}
807	+
808	+
809	+
810		for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
811		Particle subjet_v4;
812	<	subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
813	<	subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
814	<	subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
815	<	subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
816	<	topjet.add_subjet(subjet_v4);
812	>
813	>	reco::Candidate const * subjetd =  pat_topjet.daughter(k);
814	>	pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
815	>	if(patsubjetd)
816	>	{
817	>	subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
818	>	subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
819	>	subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
820	>	subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
821	>	topjet.add_subjet(subjet_v4);
822	>	//btag info
823	>	topjet.add_subFlavour(patsubjetd->partonFlavour());
824	>	topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
825	>	if (doTagInfos)
826	>	{
827	>	//ip tag info
828	>	reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
829	>	topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
830	>	topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
831	>	topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
832	>	topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
833	>	topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
834	>	topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
835	>	topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
836	>	topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
837	>	topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
838	>	topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
839	>	topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
840	>	topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));
841	>	topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
842	>	topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
843	>	topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
844	>	topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
845	>	topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
846	>	topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
847	>	topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
848	>	topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
849	>	topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
850	>	//sv tag info
851	>	reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
852	>	topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
853	>	topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
854	>	topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
855	>	topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
856	>	topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
857	>	topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
858	>	topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
859	>	std::vector<TLorentzVector> vp4; vp4.clear();
860	>	std::vector<float> vchi2; vchi2.clear();
861	>	std::vector<float> vndof; vndof.clear();
862	>	std::vector<float> vchi2ndof; vchi2ndof.clear();
863	>	std::vector<float> vtrsize; vtrsize.clear();
864	>	for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
865	>	{
866	>	ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
867	>	vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
868	>	vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());
869	>	vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());
870	>	vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());
871	>	vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());
872	>	}
873	>	topjet.add_subSecondaryVertex(vp4);
874	>	topjet.add_subVertexChi2(vchi2);
875	>	topjet.add_subVertexNdof(vndof);
876	>	topjet.add_subVertexNormalizedChi2(vchi2ndof);
877	>	topjet.add_subVertexTracksSize(vtrsize);
878	>	//try computer
879	>	edm::ESHandle<JetTagComputer> computerHandle;
880	>	iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
881	>	const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
882	>	if(computer)
883	>	{
884	>	computer->passEventSetup(iSetup);
885	>	std::vector<const reco::BaseTagInfo*>  baseTagInfos;
886	>	baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
887	>	baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );
888	>	JetTagComputer::TagInfoHelper helper(baseTagInfos);
889	>	reco::TaggingVariableList vars = computer->taggingVariables(helper);
890	>	topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
891	>	topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
892	>	topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
893	>	topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
894	>	}
895	>	}
896	>	}
897	>	else
898	>	{
899	>	//filling only standard information in case the subjet has not been pat-tified during the pattuples production
900	>	subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
901	>	subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
902	>	subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
903	>	subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
904	>	topjet.add_subjet(subjet_v4);
905	>	}
906	>
907	>
908		}
909		topjets[j].push_back(topjet);
910		}
911		}
912		}
913
914	<
914	>
915		// ------------- generator top jets -------------
916		if(doGenTopJets){
917		for(size_t j=0; j< gentopjet_sources.size(); ++j){
918
919		gentopjets[j].clear();
920	<
920	>
921		edm::Handle<reco::BasicJetCollection> reco_gentopjets;
727	–	//edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
922		iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
923
924		for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
#	Line 733 | Line 927 | NtupleWriter::analyze(const edm::Event&
927		if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
928		if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
929
930	<	TopJet gentopjet;
930	>	GenTopJet gentopjet;
931		gentopjet.set_charge(reco_gentopjet.charge());
932		gentopjet.set_pt(reco_gentopjet.pt());
933		gentopjet.set_eta(reco_gentopjet.eta());
934		gentopjet.set_phi(reco_gentopjet.phi());
935		gentopjet.set_energy(reco_gentopjet.energy());
742	–	gentopjet.set_numberOfDaughters(reco_gentopjet.numberOfDaughters());
936
937		for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
938		Particle subjet_v4;
#	Line 754 | Line 947 | NtupleWriter::analyze(const edm::Event&
947		}
948		}
949
950	+
951		// ------------- photons -------------
952		if(doPhotons){
953		for(size_t j=0; j< photon_sources.size(); ++j){
#	Line 858 | Line 1052 | NtupleWriter::analyze(const edm::Event&
1052		tr->Fill();
1053		if(doLumiInfo)
1054		previouslumiblockwasfilled=true;
1055	+
1056	+	// clean up
1057	+	m_stored_pfs.clear();
1058	+	if(doTopJetsConstituents){
1059	+	pfparticles.clear();
1060	+	}
1061	+
1062		}
1063
1064
#	Line 891 | Line 1092 | NtupleWriter::beginRun(edm::Run const& i
1092		newrun=true;
1093		}
1094
894	–	if(doJets || doTopJets){
895	–	if(doJECUncertainty){
896	–	edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
897	–	iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
898	–	JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
899	–	jecUnc = new JetCorrectionUncertainty(JetCorPar);
900	–	}
901	–	}
1095		}
1096
1097		// ------------ method called when ending the processing of a run  ------------
#	Line 951 | Line 1144 | NtupleWriter::fillDescriptions(edm::Conf
1144		descriptions.addDefault(desc);
1145		}
1146
1147	+	// ------------ method fills constituents of the pat_jet into the Ntuple and stores a reference
1148	+	// ------------ to those in the topjet
1149	+	// ------------ it is checked if the constituents have been stored already
1150	+	void
1151	+	NtupleWriter::StoreJetConstituents(pat::Jet* pat_jet, Jet* jet)
1152	+	{
1153	+	// checks if the pf cand has already been stored, only stores so far missing
1154	+	// PF candidates, then stores a reference to the pf candidate in the jet
1155	+	// also calculates the jet charge and sets it
1156	+
1157	+
1158	+	const std::vector<reco::PFCandidatePtr> jconstits = pat_jet->getPFConstituents();
1159	+
1160	+	// loop over all jet constituents and store them
1161	+	for (unsigned int i=0; i<jconstits.size(); ++i){
1162	+
1163	+	PFParticle part;
1164	+	const reco::PFCandidate* pf = jconstits[i].get();
1165	+
1166	+	// check if it has already been stored, omit if true
1167	+	int is_already_in_list = -1;
1168	+	for (unsigned int j=0; j<m_stored_pfs.size(); ++j){
1169	+
1170	+	const reco::PFCandidate* spf = m_stored_pfs[j];
1171	+	double r2 = pow(pf->eta()-spf->eta(),2) + pow(pf->phi()-spf->phi(),2);
1172	+	double dpt = fabs( pf->pt() - spf->pt() );
1173	+	if (r2<1e-10 && dpt<1e-10){
1174	+	is_already_in_list = j;
1175	+	break;
1176	+	}
1177	+	}
1178	+
1179	+	if (is_already_in_list>-1){
1180	+	jet->add_pfconstituents_index(is_already_in_list);
1181	+	continue;
1182	+	}
1183	+
1184	+	part.set_pt(pf->pt());
1185	+	part.set_eta(pf->eta());
1186	+	part.set_phi(pf->phi());
1187	+	part.set_energy(pf->energy());
1188	+	part.set_charge(pf->charge());
1189	+
1190	+	part.set_ecal_en(pf->ecalEnergy());
1191	+	part.set_hcal_en(pf->hcalEnergy());
1192	+	reco::TrackRef trackref = pf->trackRef();
1193	+	if (!trackref.isNull()){
1194	+	part.set_track_mom(trackref->p());
1195	+	}
1196	+
1197	+	PFParticle::EParticleID id = PFParticle::eX;
1198	+	switch ( pf->translatePdgIdToType(pf->pdgId()) ){
1199	+	case reco::PFCandidate::X : id = PFParticle::eX; break;
1200	+	case reco::PFCandidate::h : id = PFParticle::eH; break;
1201	+	case reco::PFCandidate::e : id = PFParticle::eE; break;
1202	+	case reco::PFCandidate::mu : id = PFParticle::eMu; break;
1203	+	case reco::PFCandidate::gamma : id = PFParticle::eGamma; break;
1204	+	case reco::PFCandidate::h0 : id = PFParticle::eH0; break;
1205	+	case reco::PFCandidate::h_HF : id = PFParticle::eH_HF; break;
1206	+	case reco::PFCandidate::egamma_HF : id = PFParticle::eEgamma_HF; break;
1207	+	}
1208	+	part.set_particleID(id);
1209	+
1210	+	pfparticles.push_back(part);
1211	+	m_stored_pfs.push_back(pf);
1212	+
1213	+	// add a reference to the particle
1214	+	jet->add_pfconstituents_index(pfparticles.size()-1);
1215	+
1216	+	}
1217	+
1218	+	if(pat_jet->isPFJet()){
1219	+	jet->set_charge(pat_jet->charge());
1220	+	jet->set_neutralEmEnergyFraction (pat_jet->neutralEmEnergyFraction());
1221	+	jet->set_neutralHadronEnergyFraction (pat_jet->neutralHadronEnergyFraction());
1222	+	jet->set_chargedEmEnergyFraction (pat_jet->chargedEmEnergyFraction());
1223	+	jet->set_chargedHadronEnergyFraction (pat_jet->chargedHadronEnergyFraction());
1224	+	jet->set_muonEnergyFraction (pat_jet->muonEnergyFraction());
1225	+	jet->set_photonEnergyFraction (pat_jet->photonEnergyFraction());
1226	+	jet->set_chargedMultiplicity (pat_jet->chargedMultiplicity());
1227	+	jet->set_neutralMultiplicity (pat_jet->neutralMultiplicity());
1228	+	jet->set_muonMultiplicity (pat_jet->muonMultiplicity());
1229	+	jet->set_electronMultiplicity (pat_jet->electronMultiplicity());
1230	+	jet->set_photonMultiplicity (pat_jet->photonMultiplicity());
1231	+	}
1232	+
1233	+	return;
1234	+
1235	+	}
1236	+
1237		//define this as a plug-in
1238		DEFINE_FWK_MODULE(NtupleWriter);

Diff Legend

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