ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/VHbbAnalysis/HbbAnalyzer/plugins/HbbAnalyzerNew.cc

Comparing UserCode/VHbbAnalysis/HbbAnalyzer/plugins/HbbAnalyzerNew.cc (file contents):
Revision 1.1 by tboccali, Wed Jun 8 17:25:32 2011 UTC vs.
Revision 1.9 by tboccali, Thu Jun 30 08:56:27 2011 UTC

#	Line 17 | Line 17 | Implementation:
17		//
18		//
19
20	<	#include "UserCode/HbbAnalyzer/interface/HbbAnalyzerNew.h"
21	<	#include "UserCode/HbbAnalyzer/interface/VHbbEvent.h"
20	>	#include "VHbbAnalysis/HbbAnalyzer/interface/HbbAnalyzerNew.h"
21	>	#include "VHbbAnalysis/HbbAnalyzer/interface/VHbbEvent.h"
22	>
23	>	#include "DataFormats/GeometryVector/interface/VectorUtil.h"
24	>
25
26		#define GENPTOLOR(a) TLorentzVector((a).px(), (a).py(), (a).pz(), (a).energy())
27		#define GENPTOLORP(a) TLorentzVector((a)->px(), (a)->py(), (a)->pz(), (a)->energy())
28
29		HbbAnalyzerNew::HbbAnalyzerNew(const edm::ParameterSet& iConfig):
30	<	eleLabel_(iConfig.getUntrackedParameter<edm::InputTag>("electronTag")),
31	<	muoLabel_(iConfig.getUntrackedParameter<edm::InputTag>("muonTag")),
32	<	jetLabel_(iConfig.getUntrackedParameter<edm::InputTag>("jetTag")),
33	<	subjetLabel_(iConfig.getUntrackedParameter<edm::InputTag>("subjetTag")),
34	<	simplejet1Label_(iConfig.getUntrackedParameter<edm::InputTag>("simplejet1Tag")),
35	<	simplejet2Label_(iConfig.getUntrackedParameter<edm::InputTag>("simplejet2Tag")),
36	<	simplejet3Label_(iConfig.getUntrackedParameter<edm::InputTag>("simplejet3Tag")),
37	<	simplejet4Label_(iConfig.getUntrackedParameter<edm::InputTag>("simplejet4Tag")),
38	<	tauLabel_(iConfig.getUntrackedParameter<edm::InputTag>("tauTag")),
39	<	metLabel_(iConfig.getUntrackedParameter<edm::InputTag>("metTag")),
40	<	phoLabel_(iConfig.getUntrackedParameter<edm::InputTag>("photonTag")),
41	<	dimuLabel_(iConfig.getUntrackedParameter<edm::InputTag>("dimuTag")),
42	<	dielecLabel_(iConfig.getUntrackedParameter<edm::InputTag>("dielecTag")),
43	<	hltResults_(iConfig.getUntrackedParameter<edm::InputTag>("hltResultsTag")){
30	>	eleLabel_(iConfig.getParameter<edm::InputTag>("electronTag")),
31	>	muoLabel_(iConfig.getParameter<edm::InputTag>("muonTag")),
32	>	jetLabel_(iConfig.getParameter<edm::InputTag>("jetTag")),
33	>	subjetLabel_(iConfig.getParameter<edm::InputTag>("subjetTag")),
34	>	simplejet1Label_(iConfig.getParameter<edm::InputTag>("simplejet1Tag")),
35	>	simplejet2Label_(iConfig.getParameter<edm::InputTag>("simplejet2Tag")),
36	>	simplejet3Label_(iConfig.getParameter<edm::InputTag>("simplejet3Tag")),
37	>	simplejet4Label_(iConfig.getParameter<edm::InputTag>("simplejet4Tag")),
38	>	tauLabel_(iConfig.getParameter<edm::InputTag>("tauTag")),
39	>	metLabel_(iConfig.getParameter<edm::InputTag>("metTag")),
40	>	phoLabel_(iConfig.getParameter<edm::InputTag>("photonTag")),
41	>	dimuLabel_(iConfig.getParameter<edm::InputTag>("dimuTag")),
42	>	dielecLabel_(iConfig.getParameter<edm::InputTag>("dielecTag")),
43	>	hltResults_(iConfig.getParameter<edm::InputTag>("hltResultsTag")),
44	>	runOnMC_(iConfig.getParameter<bool>("runOnMC")), verbose_(iConfig.getUntrackedParameter<bool>("verbose")) {
45	>
46		//
47		// put the setwhatproduced etc etc
48
#	Line 311 | Line 316 | HbbAnalyzerNew::produce(edm::Event& iEve
316
317		Handle<GenParticleCollection> genParticles;
318
314	–	bool isMC=1;
319		bool printJet=0;
320
321
322	<	if(isMC){
322	>	if(runOnMC_){
323
324		int Hin=-99,Win=-99,Zin=-99,bin=-99,bbarin=-99;
325
#	Line 513 | Line 517 | HbbAnalyzerNew::produce(edm::Event& iEve
517		VHbbEvent::SimpleJet sj;
518		sj.flavour = jet_iter->partonFlavour();
519
516	–
520		sj.tche=jet_iter->bDiscriminator("trackCountingHighEffBJetTags");
521		sj.tchp=jet_iter->bDiscriminator("trackCountingHighPurBJetTags");
522		sj.jp=jet_iter->bDiscriminator("jetProbabilityBJetTags");
#	Line 524 | Line 527 | HbbAnalyzerNew::produce(edm::Event& iEve
527		sj.charge=jet_iter->jetCharge();
528		sj.ntracks=jet_iter->associatedTracks().size();
529		sj.fourMomentum=GENPTOLORP(jet_iter);
530	+	sj.chargedTracksFourMomentum=(getChargedTracksMomentum(&*(jet_iter)));
531	+
532	+	//
533	+	// add tVector
534	+	//
535	+	sj.tVector = getTvect(&(*jet_iter));
536
537		Particle::LorentzVector p4Jet = jet_iter->p4();
538
539	<	if(isMC){
539	>	if(runOnMC_){
540		double minb1DR=9999.;
541		for(size_t i = 0; i < genParticles->size(); ++ i) {
542		const GenParticle & p = (*genParticles)[i];
#	Line 535 | Line 544 | HbbAnalyzerNew::produce(edm::Event& iEve
544		if(abs(id)<=6 || id==21 || id==23 || abs(id)==24){
545		double bb1DR=TMath::Sqrt((p.eta()-p4Jet.eta())*(p.eta()-p4Jet.eta())+
546		(p.phi()-p4Jet.phi())*(p.phi()-p4Jet.phi()));
547	<	if(bb1DR<minb1DR) {minb1DR=bb1DR; sj.b1BestMCid=id; if(p.mother()!=0) sj.b1BestMCmomid=p.mother()->pdgId();}
547	>	if(bb1DR<minb1DR) {minb1DR=bb1DR; sj.bestMCid=id; if(p.mother()!=0) sj.bestMCmomid=p.mother()->pdgId();}
548		}
549		}
550		} //isMC
#	Line 560 | Line 569 | HbbAnalyzerNew::produce(edm::Event& iEve
569		sj.charge=jet_iter->jetCharge();
570		sj.ntracks=jet_iter->associatedTracks().size();
571		sj.fourMomentum=GENPTOLORP(jet_iter);
572	+	sj.chargedTracksFourMomentum=(getChargedTracksMomentum(&*(jet_iter)));
573	+	sj.tVector = getTvect(&(*jet_iter));
574
575		Particle::LorentzVector p4Jet = jet_iter->p4();
576
577	<	if(isMC){
577	>	if(runOnMC_){
578		double minb2DR=9999.;
579		for(size_t i = 0; i < genParticles->size(); ++ i) {
580		const GenParticle & p = (*genParticles)[i];
#	Line 571 | Line 582 | HbbAnalyzerNew::produce(edm::Event& iEve
582		if(abs(id)<=6 || id==21 || id==23 || abs(id)==24){
583		double bb2DR=TMath::Sqrt((p.eta()-p4Jet.eta())*(p.eta()-p4Jet.eta())+
584		(p.phi()-p4Jet.phi())*(p.phi()-p4Jet.phi()));
585	<	if(bb2DR<minb2DR) {minb2DR=bb2DR; sj.b1BestMCid=id; if(p.mother()!=0) sj.b1BestMCmomid=p.mother()->pdgId();}
585	>	if(bb2DR<minb2DR) {minb2DR=bb2DR; sj.bestMCid=id; if(p.mother()!=0) sj.bestMCmomid=p.mother()->pdgId();}
586		}
587		}
588		}   //isMC
#	Line 623 | Line 634 | HbbAnalyzerNew::produce(edm::Event& iEve
634		hj.subFourMomentum.push_back(GENPTOLORP(icandJet));
635		hj.etaSub.push_back(icandJet->eta());
636		hj.phiSub.push_back(icandJet->phi());
626	–
637
638		}
639		hbbInfo->hardJets.push_back(hj);
#	Line 648 | Line 658 | HbbAnalyzerNew::produce(edm::Event& iEve
658		VHbbEvent::SimpleJet sj;
659
660		sj.flavour = subjet_iter->partonFlavour();
661	<
661	>	sj.tVector = getTvect(&(*subjet_iter));
662		sj.tche=subjet_iter->bDiscriminator("trackCountingHighEffBJetTags");
663		sj.tchp=subjet_iter->bDiscriminator("trackCountingHighPurBJetTags");
664		sj.jp=subjet_iter->bDiscriminator("jetProbabilityBJetTags");
#	Line 659 | Line 669 | HbbAnalyzerNew::produce(edm::Event& iEve
669		sj.charge=subjet_iter->jetCharge();
670		sj.ntracks=subjet_iter->associatedTracks().size();
671		sj.fourMomentum=GENPTOLORP(subjet_iter);
672	<
672	>	sj.fourMomentum=(getChargedTracksMomentum(&*(subjet_iter)));
673		hbbInfo->subJets.push_back(sj);
674
675		}
#	Line 669 | Line 679 | HbbAnalyzerNew::produce(edm::Event& iEve
679		// met is calomet
680		//
681
672	–
673	–	for(edm::View<pat::MET>::const_iterator met = mets.begin(); met!=mets.end(); ++met){
674	–	hbbInfo->calomet.sumEt=met->sumEt();
675	–	hbbInfo->calomet.metSig=met->mEtSig();
676	–	hbbInfo->calomet.eLong=met->e_longitudinal();
677	–	hbbInfo->calomet.fourMomentum=GENPTOLORP(met);
678	–	}
679	–
682		edm::Handle<edm::View<pat::MET> > metTCHandle;
683		iEvent.getByLabel("patMETsTC",metTCHandle);
684		edm::View<pat::MET> metsTC = *metTCHandle;
685	<	for(edm::View<pat::MET>::const_iterator metTC = metsTC.begin(); metTC!=metsTC.end(); ++metTC){
686	<	hbbInfo->calomet.sumEt=metTC->sumEt();
687	<	hbbInfo->calomet.metSig=metTC->mEtSig();
688	<	hbbInfo->calomet.eLong=metTC->e_longitudinal();
689	<	hbbInfo->calomet.fourMomentum=GENPTOLORP(metTC);
685	>	if(metsTC.size()){
686	>	hbbInfo->tcmet.sumEt=(metsTC[0]).sumEt();
687	>	hbbInfo->tcmet.metSig=(metsTC[0]).mEtSig();
688	>	hbbInfo->tcmet.eLong=(metsTC[0]).e_longitudinal();
689	>	hbbInfo->tcmet.fourMomentum=GENPTOLOR((metsTC[0]));
690	>	if (verbose_)     std::cout <<" METTC "<<     hbbInfo->tcmet.metSig <<" " <<     hbbInfo->tcmet.sumEt<<std::endl;
691		}
692	+
693	+	if(mets.size()){
694	+	hbbInfo->calomet.sumEt=(mets[0]).sumEt();
695	+	hbbInfo->calomet.metSig=(mets[0]).mEtSig();
696	+	hbbInfo->calomet.eLong=(mets[0]).e_longitudinal();
697	+	hbbInfo->calomet.fourMomentum=GENPTOLOR((mets[0]));
698	+	if (verbose_)     std::cout <<" METTC "<<     hbbInfo->calomet.metSig <<" " <<     hbbInfo->calomet.sumEt<<std::endl;
699	+	}
700	+
701	+
702		edm::Handle<edm::View<pat::MET> > metPFHandle;
703		iEvent.getByLabel("patMETsPF",metPFHandle);
704		edm::View<pat::MET> metsPF = *metPFHandle;
705	<	for(edm::View<pat::MET>::const_iterator metPF = metsPF.begin(); metPF!=metsPF.end(); ++metPF){
706	<	hbbInfo->calomet.sumEt=metPF->sumEt();
707	<	hbbInfo->calomet.metSig=metPF->mEtSig();
708	<	hbbInfo->calomet.eLong=metPF->e_longitudinal();
709	<	hbbInfo->calomet.fourMomentum=GENPTOLORP(metPF);
705	>
706	>	if(metsPF.size()){
707	>	hbbInfo->pfmet.sumEt=(metsPF[0]).sumEt();
708	>	hbbInfo->pfmet.metSig=(metsPF[0]).mEtSig();
709	>	hbbInfo->pfmet.eLong=(metsPF[0]).e_longitudinal();
710	>	hbbInfo->pfmet.fourMomentum=GENPTOLOR((metsPF[0]));
711	>	if (verbose_)     std::cout <<" METTC "<<     hbbInfo->pfmet.metSig <<" " <<     hbbInfo->pfmet.sumEt<<std::endl;
712		}
713
714
715	+	if(verbose_){
716	+	std::cout << "METs: calomet "<<mets.size()<<" tcmet "<<metsTC.size()<<" pfmet "<<metsPF.size()<<std::endl;
717	+	}
718	+
719	+	std::cout << " INPUT MUONS "<<muons.size()<<std::endl;
720
721		for(edm::View<pat::Muon>::const_iterator mu = muons.begin(); mu!=muons.end(); ++mu){
722		VHbbEvent::MuonInfo mf;
#	Line 730 | Line 750 | HbbAnalyzerNew::produce(edm::Event& iEve
750		const reco::HitPattern& q = gTrack->hitPattern();
751		mf.globChi2=gTrack.get()->normalizedChi2();
752		mf.globNHits=q.numberOfValidMuonHits();
753	+	mf. validMuStations = q. muonStationsWithValidHits();
754		}else{
755		mf.globChi2=-99;
756		mf.globNHits=-99;
757		}
758	+	//
759	+	// add stamuon
760	+
761	+	//    if (mu->isStandAloneMuon()) {
762	+	//      reco::TrackRef sta = mu->standAloneMuon();
763	+	//
764	+	//    }
765	+
766
767		//     int muInfo[12];
768		//     fillMuBlock(mu,  muInfo);
769	<	if(isMC){
769	>	if(runOnMC_){
770		const GenParticle* muMc = mu->genLepton();
771		if(muMc!=0){
772		mf.mcId=muMc->pdgId();
#	Line 748 | Line 777 | HbbAnalyzerNew::produce(edm::Event& iEve
777		hbbInfo->muInfo.push_back(mf);
778		}
779
780	<
780	>	std::cout << " INPUT electrons "<<electrons.size()<<std::endl;
781		for(edm::View<pat::Electron>::const_iterator elec = electrons.begin(); elec!=electrons.end(); ++elec){
782		VHbbEvent::ElectronInfo ef;
783		ef.fourMomentum=GENPTOLORP(elec);
#	Line 761 | Line 790 | HbbAnalyzerNew::produce(edm::Event& iEve
790		ef.hIso=elec->hcalIso();
791		Geom::Phi<double> deltaphi(elec->superCluster()->phi()-atan2(hbbInfo->calomet.fourMomentum.Py(),hbbInfo->calomet.fourMomentum.Px()));
792		ef.acop = deltaphi.value();
793	<
794	<	if(isMC){
793	>	//
794	>	// fill eleids
795	>	//
796	>	ef.id95 = elec->electronID("simpleEleId95cIso");
797	>	ef.id85 = elec->electronID("simpleEleId85cIso");
798	>	ef.id70 = elec->electronID("simpleEleId70cIso");
799	>	ef.id95r = elec->electronID("simpleEleId95relIso");
800	>	ef.id70r = elec->electronID("simpleEleId70relIso");
801	>	ef.id85r = elec->electronID("simpleEleId85relIso");
802	>
803	>	if(runOnMC_){
804		const GenParticle* elecMc = elec->genLepton();
805		if(elecMc!=0){
806		ef.mcId=elecMc->pdgId();
#	Line 773 | Line 811 | HbbAnalyzerNew::produce(edm::Event& iEve
811		hbbInfo->eleInfo.push_back(ef);
812		}
813
814	+
815	+	std::cout << " INPUT taus "<<taus.size()<<std::endl;
816		for(edm::View<pat::Tau>::const_iterator tau = taus.begin(); tau!=taus.end(); ++tau){
817		VHbbEvent::TauInfo tf;
818		tf.fourMomentum=GENPTOLORP(tau);
#	Line 862 | Line 902 | HbbAnalyzerNew::produce(edm::Event& iEve
902		const reco::HitPattern& q = gTrack->hitPattern();
903		df.daughter1.globNHits=q.numberOfValidMuonHits();
904		df.daughter1.globChi2=gTrack.get()->normalizedChi2();
905	+	df.daughter1.validMuStations = q. muonStationsWithValidHits();
906		}
907		if(muonDau1.isGlobalMuon()){
908		TrackRef gTrack = muonDau1.globalTrack();
909		const reco::HitPattern& q = gTrack->hitPattern();
910		df.daughter2.globNHits=q.numberOfValidMuonHits();
911		df.daughter2.globChi2=gTrack.get()->normalizedChi2();
912	+	df.daughter2.validMuStations = q. muonStationsWithValidHits();
913		}
914
915		}
#	Line 900 | Line 942 | HbbAnalyzerNew::produce(edm::Event& iEve
942		df.daughter1.hIso = elecDau0.hcalIso();
943		df.daughter2.hIso = elecDau1.hcalIso();
944
945	<
945	>	// ids
946	>	df.daughter1.id95 = elecDau0.electronID("simpleEleId95cIso");
947	>	df.daughter1.id85 = elecDau0.electronID  ("simpleEleId85cIso");
948	>	df.daughter1.id70 = elecDau0.electronID  ("simpleEleId70cIso");
949	>	df.daughter1.id95r = elecDau0.electronID ("simpleEleId95relIso");
950	>	df.daughter1.id85r = elecDau0.electronID ("simpleEleId85relIso");
951	>	df.daughter1.id70r = elecDau0.electronID ("simpleEleId70relIso");
952	>
953	>
954	>	df.daughter2.id95 = elecDau1.electronID("simpleEleId95cIso");
955	>	df.daughter2.id85 = elecDau1.electronID  ("simpleEleId85cIso");
956	>	df.daughter2.id70 = elecDau1.electronID  ("simpleEleId70cIso");
957	>	df.daughter2.id95r = elecDau1.electronID ("simpleEleId95relIso");
958	>	df.daughter2.id85r = elecDau1.electronID ("simpleEleId85relIso");
959	>	df.daughter2.id70r = elecDau1.electronID ("simpleEleId70relIso");
960	>
961		hbbInfo->diElectronInfo.push_back(df);
962
963		}
964	+	if (verbose_){
965	+	std::cout <<" Pushing hbbInfo "<<std::endl;
966	+	std::cout <<" SimpleJets1 = "<<hbbInfo->simpleJets.size()<<std::endl<<
967	+	" SimpleJets2 = "<<hbbInfo->simpleJets2.size()<<std::endl<<
968	+	" SubJets = "<<hbbInfo->subJets.size()<<std::endl<<
969	+	" HardJets = "<<hbbInfo->hardJets.size()<<std::endl<<
970	+	" Muons = "<<hbbInfo->muInfo.size()<<std::endl<<
971	+	" Electrons = "<<hbbInfo->eleInfo.size()<<std::endl<<
972	+	" Taus = "<<hbbInfo->tauInfo.size()<<std::endl<<
973	+	" Electrons = "<<hbbInfo->eleInfo.size()<<std::endl<<
974	+	"--------------------- "<<std::endl;
975	+	}
976	+
977	+
978		iEvent.put(hbbInfo);
979
980	+
981		}
982
983		void
#	Line 939 | Line 1011 | void
1011		HbbAnalyzerNew::endJob() {
1012		}
1013
1014	+	TVector2 HbbAnalyzerNew::getTvect( const pat::Jet* patJet ){
1015	+
1016	+	TVector2 t_Vect(0,0);
1017	+	TVector2 null(0,0);
1018	+	TVector2 ci(0,0);
1019	+	TLorentzVector pi(0,0,0,0);
1020	+	TLorentzVector J(0,0,0,0);
1021	+	TVector2 r(0,0);
1022	+	double patJetpfcPt = 1e10;
1023	+	double r_mag = 1e10;
1024	+	unsigned int nOfconst = 0;
1025	+
1026	+
1027	+	if (patJet->isPFJet() == false) {
1028	+	return t_Vect;
1029	+	}
1030	+
1031	+
1032	+	//re-reconstruct the jet direction with the charged tracks
1033	+	std::vector<reco::PFCandidatePtr>
1034	+	patJetpfc = patJet->getPFConstituents();
1035	+	for(size_t idx = 0; idx < patJetpfc.size(); idx++){
1036	+	if( patJetpfc.at(idx)->charge() != 0 ){
1037	+	pi.SetPtEtaPhiE( patJetpfc.at(idx)->pt(), patJetpfc.at(idx)->eta(), patJetpfc.at(idx)->phi(), patJetpfc.at(idx)->energy() );
1038	+	J += pi;
1039	+	nOfconst++;
1040	+	}
1041	+	}
1042	+	// if there are less than two charged tracks do not calculate the pull (there is not enough info). It returns a null vector
1043	+
1044	+	if( nOfconst < 2 )
1045	+	return null;
1046	+
1047	+
1048	+
1049	+	TVector2 v_J( J.Rapidity(), J.Phi() );
1050	+	//calculate TVector using only charged tracks
1051	+	for(size_t idx = 0; idx < patJetpfc.size(); idx++){
1052	+	if( patJetpfc.at(idx)->charge() != 0  ){
1053	+	patJetpfcPt = patJetpfc.at(idx)->pt();
1054	+	pi.SetPtEtaPhiE( patJetpfc.at(idx)->pt(), patJetpfc.at(idx)->eta(), patJetpfc.at(idx)->phi(), patJetpfc.at(idx)->energy() );
1055	+	r.Set( pi.Rapidity() - J.Rapidity(), Geom::deltaPhi( patJetpfc.at(idx)->phi(), J.Phi() ) );
1056	+	r_mag = r.Mod();
1057	+	t_Vect += ( patJetpfcPt / J.Pt() ) * r_mag * r;
1058	+	}
1059	+	}
1060	+
1061	+
1062	+	return t_Vect;
1063	+
1064	+	}
1065	+
1066	+	TLorentzVector HbbAnalyzerNew::getChargedTracksMomentum(const pat::Jet* patJet ){
1067	+	//  return TLorentzVector();
1068	+	TLorentzVector pi(0,0,0,0);
1069	+	TLorentzVector v_j1(0,0,0,0);
1070	+
1071	+
1072	+	//  std::cout <<"fff ECCCCCCOOOOO "<<patJet->isPFJet()<<std::endl;
1073	+
1074	+	if (patJet->isPFJet() == false ){
1075	+	v_j1 = GENPTOLORP(patJet);
1076	+	return v_j1;
1077	+	}
1078	+	std::vector<reco::PFCandidatePtr>
1079	+	j1pfc = patJet->getPFConstituents();
1080	+	for(size_t idx = 0; idx < j1pfc.size(); idx++){
1081	+	if( j1pfc.at(idx)->charge() != 0 ){
1082	+	pi.SetPtEtaPhiE( j1pfc.at(idx)->pt(), j1pfc.at(idx)->eta(), j1pfc.at(idx)->phi(), j1pfc.at(idx)->energy() );
1083	+	v_j1 += pi;
1084	+	}
1085	+	}
1086	+	return v_j1;
1087	+	//re-
1088	+	}
1089	+
1090		//define this as a plug-in
1091		DEFINE_FWK_MODULE(HbbAnalyzerNew);

Diff Legend

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