HbbAnalyzer/plugins/HbbCandidateFinderAlgo.cc

#include "VHbbAnalysis/HbbAnalyzer/interface/HbbCandidateFinderAlgo.h"
#include "VHbbAnalysis/HbbAnalyzer/interface/VHbbCandidateTools.h"

#include <iostream>
#include<cstdlib>

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;
  }
};

float HbbCandidateFinderAlgo::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 HbbCandidateFinderAlgo::run (const VHbbEvent* event, 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;

  //
  // now see which kind of andidate this can be
  // 
  VHbbCandidateTools selector(verbose_);

  VHbbCandidate result;
  bool ok = false;
  //
  // first: hZmumu
  //
  
  if (verbose_){
    std::cout <<" START SELECTION "<<std::endl;
  }
  
  result = selector.getHZmumuCandidate(temp,ok);
  if ( ok == true ){
    result.setCandidateType(VHbbCandidate::Zmumu);
    candidates.push_back(result);
  }else{
    //      HZee
    result = selector. getHZeeCandidate(temp,ok);
    if ( ok == true ){
      result.setCandidateType(VHbbCandidate::Zee);
      candidates.push_back(result);
      return;
    }else{
      //HWmunu
      result = selector. getHWmunCandidate(temp,ok);
      if ( ok == true ){
        result.setCandidateType(VHbbCandidate::Wmun);
        candidates.push_back(result);
        return;
      }else{
        // HWenu
        result = selector. getHWenCandidate(temp,ok);
        if ( ok == true ){
          result.setCandidateType(VHbbCandidate::Wen);
          candidates.push_back(result);
          return;
        }else{
          // HZnn
          result = selector. getHZnnCandidate(temp,ok);
          if ( ok == true ){
            result.setCandidateType(VHbbCandidate::Znn);
            candidates.push_back(result);
            return;
          }
        }
      }
    }
  }
  return;
}

void HbbCandidateFinderAlgo::findMET(const VHbbEvent::METInfo & met, std::vector<VHbbEvent::METInfo>& out){
  //

  //  just preselection: met significance > 2 

    if (met.metSig >2 ) out.push_back(met);
    if (verbose_){
      std::cout <<" CandidateFinder: Input MET = "<<met.metSig<<" Output MET = "<<out.size()<<std::endl;
    }
    
}


bool HbbCandidateFinderAlgo::findDiJets (const std::vector<VHbbEvent::SimpleJet>& jets, VHbbEvent::SimpleJet& j1, VHbbEvent::SimpleJet& j2,std::vector<VHbbEvent::SimpleJet>& addJets){
  
 std::vector<VHbbEvent::SimpleJet> tempJets;

 if (verbose_){
   std::cout <<" CandidateFinder: Input Jets = "<<jets.size()<<std::endl;
 }

 for (unsigned int i=0 ; i< jets.size(); ++i){
   if (jets[i].fourMomentum.Pt()> jetPtThreshold)
     tempJets.push_back(jets[i]);
 }
 
 CompareBTag  bTagComparator;

 
 if (verbose_){
   std::cout <<" CandidateFinder: Intermediate Jets = "<<tempJets.size()<<std::endl;
 }


 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;

 if (verbose_){
   std::cout <<" CandidateFinder: Output Jets = "<<2<<" Additional = "<<addJets.size()<<std::endl;
 }


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


}

void HbbCandidateFinderAlgo::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()>15) {
      out.push_back(*it);
    }
  }

    if (verbose_){
      std::cout <<" CandidateFinder: Input Muons = "<<muons.size()<<" Output Muons = "<<out.size()<<std::endl;
    }


}


void HbbCandidateFinderAlgo::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>  &&
        fabs((*it).fourMomentum.Eta()) < 2.5 &&
        !(      fabs((*it).fourMomentum.Eta()) < 1.57 &&        fabs((*it).fourMomentum.Eta()) > 1.44) &&
        (*it).fourMomentum.Pt()>15 //  I use the minimum ok for both Z and W
        ){
      out.push_back(*it);
    }  
  }
    if (verbose_){
      std::cout <<" CandidateFinder: Input Electrons = "<<electrons.size()<<" Output Electrons = "<<out.size()<<std::endl;
    }


}

Revision:	1.3
Committed:	Thu Jun 30 09:05:39 2011 UTC (13 years, 10 months ago) by tboccali
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Jun30th2011
Changes since 1.2:	+1 -1 lines
Log Message:	fix met
#	User	Rev	Content
1	tboccali	1.1	#include "VHbbAnalysis/HbbAnalyzer/interface/HbbCandidateFinderAlgo.h"
2			#include "VHbbAnalysis/HbbAnalyzer/interface/VHbbCandidateTools.h"
3
4			#include <iostream>
5			#include<cstdlib>
6
7			struct CompareJetPt {
8			bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
9			return j1.fourMomentum.Pt() > j2.fourMomentum.Pt();
10			}
11			};
12
13			struct CompareBTag {
14			bool operator()(const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
15			return j1.csv > j2.csv;
16			}
17			};
18
19			float HbbCandidateFinderAlgo::getDeltaTheta( const VHbbEvent::SimpleJet & j1, const VHbbEvent::SimpleJet & j2 ) const{
20
21			double deltaTheta = 1e10;
22			TLorentzVector pi(0,0,0,0);
23			TLorentzVector v_j1 = j1.chargedTracksFourMomentum;
24			TLorentzVector v_j2 = j2.chargedTracksFourMomentum;
25
26			if( v_j2.Mag() == 0
27			\|\| v_j1.Mag() == 0 )
28			return deltaTheta = 1e10;
29
30			//use j1 to calculate the pull vector
31			TVector2 t = j1.tVector;
32
33			if( t.Mod() == 0 )
34			return deltaTheta = 1e10;
35
36			Double_t dphi = v_j2.Phi()- v_j1.Phi();
37			if ( dphi > M_PI ) {
38			dphi -= 2.0*M_PI;
39			} else if ( dphi <= -M_PI ) {
40			dphi += 2.0*M_PI;
41			}
42			Double_t deltaeta = v_j2.Rapidity() - v_j1.Rapidity();
43			TVector2 BBdir( deltaeta, dphi );
44
45			deltaTheta = t.DeltaPhi(BBdir);
46
47			return deltaTheta;
48
49			}
50
51
52
53
54
55			void HbbCandidateFinderAlgo::run (const VHbbEvent* event, std::vector<VHbbCandidate> & candidates){
56			//
57			// first find the jets
58			//
59
60			VHbbEvent::SimpleJet j1,j2;
61			std::vector<VHbbEvent::SimpleJet> addJets;
62			bool foundJets = findDiJets(event->simpleJets2,j1,j2,addJets) ;
63
64			if (verbose_){
65			std::cout <<" Found Dijets: "<<foundJets<< " Additional: "<<addJets.size()<< std::endl;
66			}
67
68			if (foundJets == false) return;
69			//
70			// search for leptons
71			//
72			std::vector<VHbbEvent::MuonInfo> mu;
73			findMuons(event->muInfo,mu);
74			std::vector<VHbbEvent::ElectronInfo> ele;
75			findElectrons(event->eleInfo,ele);
76
77			std::vector<VHbbEvent::METInfo> met;
78			findMET(event->pfmet, met);
79
80			if (verbose_){
81			std::cout <<" Electrons: "<< ele.size()<<std::endl;
82			std::cout <<" Muons : "<< mu.size()<<std::endl;
83			std::cout <<" MET : "<< met.size()<<std::endl;
84			}
85			if (ele.size()<1 && mu.size() < 1 && met.size()<1) return;
86
87			//
88			// fill!
89			//
90			VHbbCandidate temp;
91			temp.H.jets.push_back(j1);
92			temp.H.jets.push_back(j2);
93			temp.H.fourMomentum = (j1).fourMomentum+(j2).fourMomentum;
94			temp.H.deltaTheta = getDeltaTheta(j1,j2);
95			// temp.H.deltaTheta = getDeltaTheta()
96			temp.additionalJets = addJets;
97			temp.V.mets = met;
98			temp.V.muons = mu;
99			temp.V.electrons = ele;
100
101			//
102			// now see which kind of andidate this can be
103			//
104	tboccali	1.3	VHbbCandidateTools selector(verbose_);
105	tboccali	1.1
106			VHbbCandidate result;
107			bool ok = false;
108			//
109			// first: hZmumu
110			//
111
112	tboccali	1.2	if (verbose_){
113			std::cout <<" START SELECTION "<<std::endl;
114			}
115
116	tboccali	1.1	result = selector.getHZmumuCandidate(temp,ok);
117	tboccali	1.2	if ( ok == true ){
118			result.setCandidateType(VHbbCandidate::Zmumu);
119			candidates.push_back(result);
120			}else{
121			// HZee
122			result = selector. getHZeeCandidate(temp,ok);
123	tboccali	1.1	if ( ok == true ){
124	tboccali	1.2	result.setCandidateType(VHbbCandidate::Zee);
125	tboccali	1.1	candidates.push_back(result);
126	tboccali	1.2	return;
127	tboccali	1.1	}else{
128	tboccali	1.2	//HWmunu
129			result = selector. getHWmunCandidate(temp,ok);
130			if ( ok == true ){
131			result.setCandidateType(VHbbCandidate::Wmun);
132			candidates.push_back(result);
133			return;
134			}else{
135			// HWenu
136			result = selector. getHWenCandidate(temp,ok);
137	tboccali	1.1	if ( ok == true ){
138	tboccali	1.2	result.setCandidateType(VHbbCandidate::Wen);
139	tboccali	1.1	candidates.push_back(result);
140			return;
141			}else{
142	tboccali	1.2	// HZnn
143			result = selector. getHZnnCandidate(temp,ok);
144	tboccali	1.1	if ( ok == true ){
145	tboccali	1.2	result.setCandidateType(VHbbCandidate::Znn);
146	tboccali	1.1	candidates.push_back(result);
147			return;
148			}
149			}
150	tboccali	1.2	}
151	tboccali	1.1	}
152	tboccali	1.2	}
153	tboccali	1.1	return;
154			}
155
156			void HbbCandidateFinderAlgo::findMET(const VHbbEvent::METInfo & met, std::vector<VHbbEvent::METInfo>& out){
157			//
158	tboccali	1.2
159	tboccali	1.1	// just preselection: met significance > 2
160	tboccali	1.2
161	tboccali	1.1	if (met.metSig >2 ) out.push_back(met);
162	tboccali	1.2	if (verbose_){
163			std::cout <<" CandidateFinder: Input MET = "<<met.metSig<<" Output MET = "<<out.size()<<std::endl;
164			}
165
166	tboccali	1.1	}
167
168
169			bool HbbCandidateFinderAlgo::findDiJets (const std::vector<VHbbEvent::SimpleJet>& jets, VHbbEvent::SimpleJet& j1, VHbbEvent::SimpleJet& j2,std::vector<VHbbEvent::SimpleJet>& addJets){
170
171			std::vector<VHbbEvent::SimpleJet> tempJets;
172	tboccali	1.2
173			if (verbose_){
174			std::cout <<" CandidateFinder: Input Jets = "<<jets.size()<<std::endl;
175			}
176
177	tboccali	1.1	for (unsigned int i=0 ; i< jets.size(); ++i){
178			if (jets[i].fourMomentum.Pt()> jetPtThreshold)
179			tempJets.push_back(jets[i]);
180			}
181
182			CompareBTag bTagComparator;
183
184	tboccali	1.2
185			if (verbose_){
186			std::cout <<" CandidateFinder: Intermediate Jets = "<<tempJets.size()<<std::endl;
187			}
188
189
190	tboccali	1.1	if (tempJets.size()<2) return false;
191
192			std::sort(tempJets.begin(), tempJets.end(), bTagComparator);
193
194			j1 = tempJets[0];
195			j2 = tempJets[1];
196			//
197			// additional jets
198			//
199			if (tempJets.size()>2){
200			for (unsigned int i=2 ; i< tempJets.size(); ++i){
201			addJets.push_back(tempJets[i]);
202			}
203			}
204			CompareJetPt ptComparator;
205
206	tboccali	1.2	if (verbose_){
207			std::cout <<" CandidateFinder: Output Jets = "<<2<<" Additional = "<<addJets.size()<<std::endl;
208			}
209
210
211	tboccali	1.1	std::sort(addJets.begin(), addJets.end(), ptComparator);
212			return true;
213	tboccali	1.2
214
215	tboccali	1.1	}
216
217			void HbbCandidateFinderAlgo::findMuons(const std::vector<VHbbEvent::MuonInfo>& muons, std::vector<VHbbEvent::MuonInfo>& out){
218			/* Use:
219			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:
220
221			We use RECO Muons that are both Global and Tracker
222			chi2/ndof < 10 for the global muon fit
223			The track associated to the muon must have
224			>= 1 pixel hits
225			>= 10 pixel + strip hits
226			>= 1 valid hit in the muon chambers
227			>= 2 muon stations
228			\|dxy\| < 0.2
229			\|eta\| < 2.4
230			Relative combined isolation (R) is required to be < 0.15
231			R = [Sum pT(trks) + Et(ECAL) + Et(HCAL)] / pT(mu) computed in a cone of radius 0.3 in eta-phi
232			pT(mu) > 20 GeV
233			*/
234			for (std::vector<VHbbEvent::MuonInfo>::const_iterator it = muons.begin(); it!= muons.end(); ++it){
235			if (
236			(*it). globChi2<10 &&
237			(*it).nPixelHits>= 1 &&
238			(*it).globNHits >= 1 &&
239			(*it).nHits >= 10 &&
240			(*it).cat ==1 &&
241			(*it).validMuStations >=2 &&
242			(*it).ipDb<.2 &&
243			((it).hIso+(it).eIso+(it).tIso)/(it).fourMomentum.Pt()<.15 &&
244			(*it).fourMomentum.Eta()<2.4 &&
245			(*it).fourMomentum.Pt()>15) {
246			out.push_back(*it);
247			}
248			}
249	tboccali	1.2
250			if (verbose_){
251			std::cout <<" CandidateFinder: Input Muons = "<<muons.size()<<" Output Muons = "<<out.size()<<std::endl;
252			}
253
254
255
256	tboccali	1.1	}
257
258
259			void HbbCandidateFinderAlgo::findElectrons(const std::vector<VHbbEvent::ElectronInfo>& electrons, std::vector<VHbbEvent::ElectronInfo>& out){
260			/*
261			We adopt the standard cut-based selection from VBTF described in detail here.
262
263			Z -> ee
264			gsf electrons
265			VBTF WP95
266			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
267			pT(e) > 20
268
269			W -> e nu
270			gsf electrons
271			VBTF WP80
272			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
273			pT(e) > 30
274			*/
275
276			for (std::vector<VHbbEvent::ElectronInfo>::const_iterator it = electrons.begin(); it!= electrons.end(); ++it){
277			if (
278			// fake
279			// (*it).id95> &&
280			fabs((*it).fourMomentum.Eta()) < 2.5 &&
281			!( fabs((it).fourMomentum.Eta()) < 1.57 && fabs((it).fourMomentum.Eta()) > 1.44) &&
282			(*it).fourMomentum.Pt()>15 // I use the minimum ok for both Z and W
283			){
284			out.push_back(*it);
285			}
286			}
287	tboccali	1.2	if (verbose_){
288			std::cout <<" CandidateFinder: Input Electrons = "<<electrons.size()<<" Output Electrons = "<<out.size()<<std::endl;
289			}
290
291
292	tboccali	1.1	}
293