HbbAnalyzer/plugins/HbbCandidateFinder.cc

#include "VHbbAnalysis/HbbAnalyzer/interface/HbbCandidateFinder.h"

struct CompareJetPt {
  bool operator()( const VHbbEvent::SimpleJet& j1, const  VHbbEvent::SimpleJet& j2 ) const {
    return j1.fourMomentum.Pt() > j2.fourMomentum.Pt();
  }
};

  struct CompareBTag {
    bool operator()(const  VHbbEvent::SimpleJet& j1, const  VHbbEvent::SimpleJet& j2 ) const {
      return j1.csv > j2.csv;
    }
  };

HbbCandidateFinder::HbbCandidateFinder(const edm::ParameterSet& iConfig):   vhbbevent_(iConfig.getParameter<edm::InputTag>("VHbbEventLabel")), verbose_ (iConfig.getParameter<bool>("verbose")), jetPtThreshold(iConfig.getParameter<double>("jetPtThreshold")){
  produces<std::vector<VHbbCandidate > >();
}

HbbCandidateFinder::~HbbCandidateFinder(){}

void HbbCandidateFinder::beginJob(){}
void HbbCandidateFinder::endJob(){}


float HbbCandidateFinder::getDeltaTheta( const VHbbEvent::SimpleJet & j1, const VHbbEvent::SimpleJet & j2 ) const{

 double deltaTheta = 1e10;
 TLorentzVector pi(0,0,0,0);
 TLorentzVector v_j1 = j1.chargedTracksFourMomentum;
 TLorentzVector v_j2 = j2.chargedTracksFourMomentum;
 
 if( v_j2.Mag() == 0 
     || v_j1.Mag() == 0 )
   return deltaTheta = 1e10;
 
 //use j1 to calculate the pull vector
 TVector2 t = j1.tVector;
 
 if( t.Mod() == 0 )
   return deltaTheta = 1e10;
 
 Double_t dphi =  v_j2.Phi()- v_j1.Phi();
 if ( dphi > M_PI ) {
   dphi -= 2.0*M_PI;
 } else if ( dphi <= -M_PI ) {
   dphi += 2.0*M_PI;
 }
 Double_t deltaeta = v_j2.Rapidity() - v_j1.Rapidity();
 TVector2 BBdir( deltaeta, dphi );
 
 deltaTheta = t.DeltaPhi(BBdir);
 
 return deltaTheta;
 
}


void HbbCandidateFinder::produce( edm::Event& iEvent, const edm::EventSetup& iEventSetup){
  
  std::auto_ptr<std::vector<VHbbCandidate> >  vHbbCandidates( new std::vector<VHbbCandidate>  );

  edm::Handle<VHbbEvent>  vHbbEvent;
  //  iEvent.getByLabel(vhbbevent_, vHbbEvent);
  iEvent.getByType(vHbbEvent);
  

  //
  // start searching for candidates
  //

  //  hbbCandidateFinderAlgo(vHbbCandidates, vHbbEvent-> result());
  // do nothing for a test
  
  run(vHbbEvent.product(),vHbbCandidates);


  if (verbose_)
    std::cout <<" Pushing VHbb candidates: "<<vHbbCandidates->size()<<std::endl;

  iEvent.put(vHbbCandidates);  

}

void HbbCandidateFinder::run (const VHbbEvent* event, std::auto_ptr<std::vector<VHbbCandidate> > & candidates){
  //
  // first find the jets
  //

  VHbbEvent::SimpleJet j1,j2;
  std::vector<VHbbEvent::SimpleJet> addJets;
  bool foundJets = findDiJets(event->simpleJets2,j1,j2,addJets) ;

  if (verbose_){
    std::cout <<" Found Dijets: "<<foundJets<< " Additional: "<<addJets.size()<< std::endl;
  }
  
  if (foundJets == false) return;
  //
  // search for leptons
  //
  std::vector<VHbbEvent::MuonInfo> mu;
  findMuons(event->muInfo,mu);
  std::vector<VHbbEvent::ElectronInfo> ele;
 findElectrons(event->eleInfo,ele);
  
  std::vector<VHbbEvent::METInfo> met;
  findMET(event->pfmet, met);

  if (verbose_){
    std::cout <<" Electrons: "<< ele.size()<<std::endl;
    std::cout <<" Muons    : "<< mu.size()<<std::endl;
    std::cout <<" MET      : "<< met.size()<<std::endl;
  }
  if (ele.size()<1 && mu.size() < 1 && met.size()<1) return;

  //
  // fill!
  //
  VHbbCandidate temp;
  temp.H.jets.push_back(j1);
  temp.H.jets.push_back(j2);
  temp.H.fourMomentum = (j1).fourMomentum+(j2).fourMomentum;
  temp.H.deltaTheta = getDeltaTheta(j1,j2);
  //  temp.H.deltaTheta = getDeltaTheta()
  temp.additionalJets = addJets;
  temp.V.mets = met;
  temp.V.muons = mu;
  temp.V.electrons = ele;
  
  candidates->push_back(temp);
}

void HbbCandidateFinder::findMET(const VHbbEvent::METInfo & met, std::vector<VHbbEvent::METInfo>& out){
  //
  //  just preselection: met significance > 2 
    if (met.metSig >2 ) out.push_back(met);
  
}


bool HbbCandidateFinder::findDiJets (const std::vector<VHbbEvent::SimpleJet>& jets, VHbbEvent::SimpleJet& j1, VHbbEvent::SimpleJet& j2,std::vector<VHbbEvent::SimpleJet>& addJets){
  
 std::vector<VHbbEvent::SimpleJet> tempJets;
 
 for (unsigned int i=0 ; i< jets.size(); ++i){
   if (jets[i].fourMomentum.Pt()> jetPtThreshold)
     tempJets.push_back(jets[i]);
 }
 
 CompareBTag  bTagComparator;

 if (tempJets.size()<2) return false;
 
 std::sort(tempJets.begin(), tempJets.end(), bTagComparator);
 
 j1 = tempJets[0];
 j2 = tempJets[1];
 //
 // additional jets
 //
 if (tempJets.size()>2){
   for (unsigned int i=2 ; i< tempJets.size(); ++i){
     addJets.push_back(tempJets[i]);
   }
 }
  CompareJetPt ptComparator;

  std::sort(addJets.begin(), addJets.end(), ptComparator);
 return true;
}

void HbbCandidateFinder::findMuons(const std::vector<VHbbEvent::MuonInfo>& muons, std::vector<VHbbEvent::MuonInfo>& out){
  /* Use:
For both W -> mu nu and Z -> mu mu, we adopt the standard VBTF muon selection described in VbtfWmunuBaselineSelection. The explicit cuts are reproduced here:

    We use RECO Muons that are both Global and Tracker
    chi2/ndof < 10 for the global muon fit
    The track associated to the muon must have
        >= 1 pixel hits
        >= 10 pixel + strip hits
        >= 1 valid hit in the muon chambers
        >= 2 muon stations
        |dxy| < 0.2
        |eta| < 2.4 
    Relative combined isolation (R) is required to be < 0.15
        R = [Sum pT(trks) + Et(ECAL) + Et(HCAL)] / pT(mu) computed in a cone of radius 0.3 in eta-phi 
    pT(mu) > 20 GeV 
  */
  for (std::vector<VHbbEvent::MuonInfo>::const_iterator it = muons.begin(); it!= muons.end(); ++it){
    if (
        (*it). globChi2<10 &&
        (*it).nPixelHits>= 1 &&
        (*it).globNHits >= 1 &&
        (*it).nHits >= 10 &&
        (*it).cat ==1 && 
        (*it).validMuStations >=2 &&
        (*it).ipDb<.2 &&
        ((*it).hIso+(*it).eIso+(*it).tIso)/(*it).fourMomentum.Pt()<.15 &&
        (*it).fourMomentum.Eta()<2.4 &&
        (*it).fourMomentum.Pt()>20) {
      out.push_back(*it);
    }
  }
}


void HbbCandidateFinder::findElectrons(const std::vector<VHbbEvent::ElectronInfo>& electrons, std::vector<VHbbEvent::ElectronInfo>& out){
  /*
We adopt the standard cut-based selection from VBTF described in detail here.

    Z -> ee
        gsf electrons
        VBTF WP95
        |eta|<2.5, excluding the gap 1.44 < |eta| < 1.57
        pT(e) > 20 

    W -> e nu
        gsf electrons
        VBTF WP80
        |eta|<2.5, excluding the gap 1.44 < |eta| < 1.57
        pT(e) > 30 
  */

  for (std::vector<VHbbEvent::ElectronInfo>::const_iterator it = electrons.begin(); it!= electrons.end(); ++it){
    if (
        // fake
        //      (*it).id95>  &&
        std::abs((*it).fourMomentum.Eta()) < 2.5 &&
        !(      std::abs((*it).fourMomentum.Eta()) < 1.57 &&    std::abs((*it).fourMomentum.Eta()) > 1.44) &&
        (*it).fourMomentum.Pt()>20 //  I use the minimum ok for both Z and W
        ){
      out.push_back(*it);
    }  
  }
}


//define this as a plug-in
DEFINE_FWK_MODULE(HbbCandidateFinder);
Revision:	1.4
Committed:	Tue Jun 21 15:49:01 2011 UTC (13 years, 10 months ago) by tboccali
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+64 -33 lines
Log Message:	maybe working preselection?
#	User	Rev	Content
1	tboccali	1.1	#include "VHbbAnalysis/HbbAnalyzer/interface/HbbCandidateFinder.h"
2
3	tboccali	1.3	struct CompareJetPt {
4			bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
5			return j1.fourMomentum.Pt() > j2.fourMomentum.Pt();
6			}
7			};
8	tboccali	1.1
9	tboccali	1.3	struct CompareBTag {
10			bool operator()(const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
11			return j1.csv > j2.csv;
12			}
13			};
14	tboccali	1.1
15	tboccali	1.3	HbbCandidateFinder::HbbCandidateFinder(const edm::ParameterSet& iConfig): vhbbevent_(iConfig.getParameter<edm::InputTag>("VHbbEventLabel")), verbose_ (iConfig.getParameter<bool>("verbose")), jetPtThreshold(iConfig.getParameter<double>("jetPtThreshold")){
16	tboccali	1.1	produces<std::vector<VHbbCandidate > >();
17			}
18
19			HbbCandidateFinder::~HbbCandidateFinder(){}
20
21			void HbbCandidateFinder::beginJob(){}
22			void HbbCandidateFinder::endJob(){}
23
24
25	tboccali	1.4	float HbbCandidateFinder::getDeltaTheta( const VHbbEvent::SimpleJet & j1, const VHbbEvent::SimpleJet & j2 ) const{
26
27			double deltaTheta = 1e10;
28			TLorentzVector pi(0,0,0,0);
29			TLorentzVector v_j1 = j1.chargedTracksFourMomentum;
30			TLorentzVector v_j2 = j2.chargedTracksFourMomentum;
31
32			if( v_j2.Mag() == 0
33			\|\| v_j1.Mag() == 0 )
34			return deltaTheta = 1e10;
35
36			//use j1 to calculate the pull vector
37			TVector2 t = j1.tVector;
38
39			if( t.Mod() == 0 )
40			return deltaTheta = 1e10;
41
42			Double_t dphi = v_j2.Phi()- v_j1.Phi();
43			if ( dphi > M_PI ) {
44			dphi -= 2.0*M_PI;
45			} else if ( dphi <= -M_PI ) {
46			dphi += 2.0*M_PI;
47			}
48			Double_t deltaeta = v_j2.Rapidity() - v_j1.Rapidity();
49			TVector2 BBdir( deltaeta, dphi );
50
51			deltaTheta = t.DeltaPhi(BBdir);
52
53			return deltaTheta;
54
55			}
56
57	tboccali	1.1
58
59
60			void HbbCandidateFinder::produce( edm::Event& iEvent, const edm::EventSetup& iEventSetup){
61
62			std::auto_ptr<std::vector<VHbbCandidate> > vHbbCandidates( new std::vector<VHbbCandidate> );
63
64			edm::Handle<VHbbEvent> vHbbEvent;
65	tboccali	1.3	// iEvent.getByLabel(vhbbevent_, vHbbEvent);
66			iEvent.getByType(vHbbEvent);
67	tboccali	1.1
68
69			//
70			// start searching for candidates
71			//
72
73			// hbbCandidateFinderAlgo(vHbbCandidates, vHbbEvent-> result());
74			// do nothing for a test
75
76	tboccali	1.3	run(vHbbEvent.product(),vHbbCandidates);
77
78	tboccali	1.2
79	tboccali	1.3	if (verbose_)
80			std::cout <<" Pushing VHbb candidates: "<<vHbbCandidates->size()<<std::endl;
81	tboccali	1.2
82	tboccali	1.1	iEvent.put(vHbbCandidates);
83
84			}
85
86	tboccali	1.2	void HbbCandidateFinder::run (const VHbbEvent* event, std::auto_ptr<std::vector<VHbbCandidate> > & candidates){
87			//
88			// first find the jets
89			//
90	tboccali	1.3
91			VHbbEvent::SimpleJet j1,j2;
92			std::vector<VHbbEvent::SimpleJet> addJets;
93			bool foundJets = findDiJets(event->simpleJets2,j1,j2,addJets) ;
94
95			if (verbose_){
96			std::cout <<" Found Dijets: "<<foundJets<< " Additional: "<<addJets.size()<< std::endl;
97			}
98	tboccali	1.2
99	tboccali	1.3	if (foundJets == false) return;
100	tboccali	1.2	//
101	tboccali	1.4	// search for leptons
102	tboccali	1.2	//
103	tboccali	1.3	std::vector<VHbbEvent::MuonInfo> mu;
104			findMuons(event->muInfo,mu);
105			std::vector<VHbbEvent::ElectronInfo> ele;
106	tboccali	1.4	findElectrons(event->eleInfo,ele);
107	tboccali	1.3
108	tboccali	1.4	std::vector<VHbbEvent::METInfo> met;
109			findMET(event->pfmet, met);
110
111	tboccali	1.3	if (verbose_){
112			std::cout <<" Electrons: "<< ele.size()<<std::endl;
113			std::cout <<" Muons : "<< mu.size()<<std::endl;
114	tboccali	1.4	std::cout <<" MET : "<< met.size()<<std::endl;
115	tboccali	1.3	}
116	tboccali	1.4	if (ele.size()<1 && mu.size() < 1 && met.size()<1) return;
117	tboccali	1.2
118			//
119			// fill!
120			//
121	tboccali	1.4	VHbbCandidate temp;
122			temp.H.jets.push_back(j1);
123			temp.H.jets.push_back(j2);
124			temp.H.fourMomentum = (j1).fourMomentum+(j2).fourMomentum;
125			temp.H.deltaTheta = getDeltaTheta(j1,j2);
126			// temp.H.deltaTheta = getDeltaTheta()
127			temp.additionalJets = addJets;
128			temp.V.mets = met;
129			temp.V.muons = mu;
130			temp.V.electrons = ele;
131
132			candidates->push_back(temp);
133			}
134	tboccali	1.2
135	tboccali	1.4	void HbbCandidateFinder::findMET(const VHbbEvent::METInfo & met, std::vector<VHbbEvent::METInfo>& out){
136			//
137			// just preselection: met significance > 2
138			if (met.metSig >2 ) out.push_back(met);
139
140	tboccali	1.2	}
141
142	tboccali	1.4
143	tboccali	1.3	bool HbbCandidateFinder::findDiJets (const std::vector<VHbbEvent::SimpleJet>& jets, VHbbEvent::SimpleJet& j1, VHbbEvent::SimpleJet& j2,std::vector<VHbbEvent::SimpleJet>& addJets){
144	tboccali	1.2
145	tboccali	1.3	std::vector<VHbbEvent::SimpleJet> tempJets;
146
147			for (unsigned int i=0 ; i< jets.size(); ++i){
148			if (jets[i].fourMomentum.Pt()> jetPtThreshold)
149			tempJets.push_back(jets[i]);
150			}
151
152			CompareBTag bTagComparator;
153
154			if (tempJets.size()<2) return false;
155
156			std::sort(tempJets.begin(), tempJets.end(), bTagComparator);
157
158			j1 = tempJets[0];
159			j2 = tempJets[1];
160			//
161			// additional jets
162			//
163			if (tempJets.size()>2){
164			for (unsigned int i=2 ; i< tempJets.size(); ++i){
165			addJets.push_back(tempJets[i]);
166			}
167			}
168			CompareJetPt ptComparator;
169
170			std::sort(addJets.begin(), addJets.end(), ptComparator);
171			return true;
172			}
173
174			void HbbCandidateFinder::findMuons(const std::vector<VHbbEvent::MuonInfo>& muons, std::vector<VHbbEvent::MuonInfo>& out){
175			/* Use:
176			For both W -> mu nu and Z -> mu mu, we adopt the standard VBTF muon selection described in VbtfWmunuBaselineSelection. The explicit cuts are reproduced here:
177
178			We use RECO Muons that are both Global and Tracker
179			chi2/ndof < 10 for the global muon fit
180			The track associated to the muon must have
181			>= 1 pixel hits
182			>= 10 pixel + strip hits
183			>= 1 valid hit in the muon chambers
184			>= 2 muon stations
185			\|dxy\| < 0.2
186			\|eta\| < 2.4
187			Relative combined isolation (R) is required to be < 0.15
188			R = [Sum pT(trks) + Et(ECAL) + Et(HCAL)] / pT(mu) computed in a cone of radius 0.3 in eta-phi
189			pT(mu) > 20 GeV
190			*/
191			for (std::vector<VHbbEvent::MuonInfo>::const_iterator it = muons.begin(); it!= muons.end(); ++it){
192			if (
193	tboccali	1.4	(*it). globChi2<10 &&
194	tboccali	1.3	(*it).nPixelHits>= 1 &&
195	tboccali	1.4	(*it).globNHits >= 1 &&
196			(*it).nHits >= 10 &&
197			(*it).cat ==1 &&
198			(*it).validMuStations >=2 &&
199	tboccali	1.3	(*it).ipDb<.2 &&
200			((it).hIso+(it).eIso+(it).tIso)/(it).fourMomentum.Pt()<.15 &&
201			(*it).fourMomentum.Eta()<2.4 &&
202			(*it).fourMomentum.Pt()>20) {
203			out.push_back(*it);
204	tboccali	1.2	}
205			}
206			}
207
208
209	tboccali	1.3	void HbbCandidateFinder::findElectrons(const std::vector<VHbbEvent::ElectronInfo>& electrons, std::vector<VHbbEvent::ElectronInfo>& out){
210			/*
211			We adopt the standard cut-based selection from VBTF described in detail here.
212
213			Z -> ee
214			gsf electrons
215			VBTF WP95
216			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
217			pT(e) > 20
218
219			W -> e nu
220			gsf electrons
221			VBTF WP80
222			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
223			pT(e) > 30
224			*/
225
226			for (std::vector<VHbbEvent::ElectronInfo>::const_iterator it = electrons.begin(); it!= electrons.end(); ++it){
227			if (
228			// fake
229			// (*it).id95> &&
230			std::abs((*it).fourMomentum.Eta()) < 2.5 &&
231			!( std::abs((it).fourMomentum.Eta()) < 1.57 && std::abs((it).fourMomentum.Eta()) > 1.44) &&
232	tboccali	1.4	(*it).fourMomentum.Pt()>20 // I use the minimum ok for both Z and W
233	tboccali	1.3	){
234			out.push_back(*it);
235			}
236			}
237	tboccali	1.2	}
238
239	tboccali	1.1
240			//define this as a plug-in
241			DEFINE_FWK_MODULE(HbbCandidateFinder);