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.26 by peiffer, Mon Jun 10 13:33:05 2013 UTC vs.
Revision 1.27 by rkogler, Wed Jun 19 13:22:06 2013 UTC

#	Line 18 | Line 18
18		//
19
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 54 | Line 56 | NtupleWriter::NtupleWriter(const edm::Pa
56		doTaus = iConfig.getParameter<bool>("doTaus");
57		doJets = iConfig.getParameter<bool>("doJets");
58		doGenJets = iConfig.getParameter<bool>("doGenJets");
57	–	doJECUncertainty = iConfig.getParameter<bool>("doJECUncertainty");
59		doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");
60		doPhotons = iConfig.getParameter<bool>("doPhotons");
61		doMET = iConfig.getParameter<bool>("doMET");
#	Line 63 | 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
74		tr->Branch("run",&run);
#	Line 129 | 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 584 | Line 591 | NtupleWriter::analyze(const edm::Event&
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;
#	Line 595 | Line 603 | NtupleWriter::analyze(const edm::Event&
603		jet.set_phi(gen_jet->phi());
604		jet.set_energy(gen_jet->energy());
605
606	<	genjets[j].push_back(jet);
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		}
#	Line 645 | 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){
649	<	jecUnc->setJetEta(pat_jet.eta());
650	<	jecUnc->setJetPt(pat_jet.pt());
651	<	jet.set_JEC_uncertainty(jecUnc->getUncertainty(true));
652	<	}
663	>
664		jet.set_JEC_factor_raw(pat_jet.jecFactor("Uncorrected"));
665
666		jet.set_btag_simpleSecondaryVertexHighEff(pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
#	Line 664 | Line 675 | NtupleWriter::analyze(const edm::Event&
675		if(genj){
676
677		for(unsigned int k=0; k<genjets->size(); ++k){
678	<	if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta())
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;
#	Line 693 | Line 705 | NtupleWriter::analyze(const edm::Event&
705		}
706
707
696	–
708		// ------------- top jets -------------
709		if(doTopJets){
710	+
711		for(size_t j=0; j< topjet_sources.size(); ++j){
712
713		topjets[j].clear();
#	Line 721 | 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());
725	<	//       topjet.set_neutralHadronEnergyFraction(pat_topjet.neutralHadronEnergyFraction());
726	<	//       topjet.set_chargedEmEnergyFraction(pat_topjet.chargedEmEnergyFraction());
727	<	//       topjet.set_chargedHadronEnergyFraction(pat_topjet.chargedHadronEnergyFraction());
728	<	//       topjet.set_muonEnergyFraction(pat_topjet.muonEnergyFraction());
729	<	//       topjet.set_photonEnergyFraction(pat_topjet.photonEnergyFraction());
730	<	//       topjet.set_chargedMultiplicity(pat_topjet.chargedMultiplicity());
731	<	//       topjet.set_neutralMultiplicity(pat_topjet.neutralMultiplicity());
732	<	//       topjet.set_muonMultiplicity(pat_topjet.muonMultiplicity());
733	<	//       topjet.set_electronMultiplicity(pat_topjet.electronMultiplicity());
734	<	//       topjet.set_photonMultiplicity(pat_topjet.photonMultiplicity());
735	<	if(doJECUncertainty){
736	<	jecUnc->setJetEta(pat_topjet.eta());
737	<	jecUnc->setJetPt(pat_topjet.pt());
738	<	topjet.set_JEC_uncertainty( jecUnc->getUncertainty(true));
739	<	}
736	>
737		topjet.set_JEC_factor_raw( pat_topjet.jecFactor("Uncorrected"));
738
739		topjet.set_btag_simpleSecondaryVertexHighEff(pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
#	Line 769 | 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
#	Line 872 | Line 910 | NtupleWriter::analyze(const edm::Event&
910		}
911		}
912		}
913	<
913	>
914
915		// ------------- generator top jets -------------
916		if(doGenTopJets){
#	Line 881 | Line 919 | NtupleWriter::analyze(const edm::Event&
919		gentopjets[j].clear();
920
921		edm::Handle<reco::BasicJetCollection> reco_gentopjets;
884	–	//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 890 | 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());
899	–	gentopjet.set_numberOfDaughters(reco_gentopjet.numberOfDaughters());
936
937		for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
938		Particle subjet_v4;
#	Line 911 | 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 1015 | 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 1048 | Line 1092 | NtupleWriter::beginRun(edm::Run const& i
1092		newrun=true;
1093		}
1094
1051	–	if(doJets || doTopJets){
1052	–	if(doJECUncertainty){
1053	–	edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
1054	–	iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
1055	–	JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
1056	–	jecUnc = new JetCorrectionUncertainty(JetCorPar);
1057	–	}
1058	–	}
1095		}
1096
1097		// ------------ method called when ending the processing of a run  ------------
#	Line 1108 | 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