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;

  VHbbCandidate result;
  bool ok = false;
  //
  // first: hZmumu
  //
  
  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);
  
}


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


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

Revision:	1.1
Committed:	Tue Jun 28 12:03:54 2011 UTC (13 years, 10 months ago) by tboccali
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Jun28th2011
Log Message:	change candidate preselection
#	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			VHbbCandidateTools selector;
105
106			VHbbCandidate result;
107			bool ok = false;
108			//
109			// first: hZmumu
110			//
111
112			result = selector.getHZmumuCandidate(temp,ok);
113			if ( ok == true ){
114			result.setCandidateType(VHbbCandidate::Zmumu);
115			candidates.push_back(result);
116			}else{
117			// HZee
118			result = selector. getHZeeCandidate(temp,ok);
119			if ( ok == true ){
120			result.setCandidateType(VHbbCandidate::Zee);
121			candidates.push_back(result);
122			return;
123			}else{
124			//HWmunu
125			result = selector. getHWmunCandidate(temp,ok);
126			if ( ok == true ){
127			result.setCandidateType(VHbbCandidate::Wmun);
128			candidates.push_back(result);
129			return;
130			}else{
131			// HWenu
132			result = selector. getHWenCandidate(temp,ok);
133			if ( ok == true ){
134			result.setCandidateType(VHbbCandidate::Wen);
135			candidates.push_back(result);
136			return;
137			}else{
138			// HZnn
139			result = selector. getHZnnCandidate(temp,ok);
140			if ( ok == true ){
141			result.setCandidateType(VHbbCandidate::Znn);
142			candidates.push_back(result);
143			return;
144			}
145			}
146			}
147			}
148			}
149			return;
150			}
151
152			void HbbCandidateFinderAlgo::findMET(const VHbbEvent::METInfo & met, std::vector<VHbbEvent::METInfo>& out){
153			//
154			// just preselection: met significance > 2
155			if (met.metSig >2 ) out.push_back(met);
156
157			}
158
159
160			bool HbbCandidateFinderAlgo::findDiJets (const std::vector<VHbbEvent::SimpleJet>& jets, VHbbEvent::SimpleJet& j1, VHbbEvent::SimpleJet& j2,std::vector<VHbbEvent::SimpleJet>& addJets){
161
162			std::vector<VHbbEvent::SimpleJet> tempJets;
163
164			for (unsigned int i=0 ; i< jets.size(); ++i){
165			if (jets[i].fourMomentum.Pt()> jetPtThreshold)
166			tempJets.push_back(jets[i]);
167			}
168
169			CompareBTag bTagComparator;
170
171			if (tempJets.size()<2) return false;
172
173			std::sort(tempJets.begin(), tempJets.end(), bTagComparator);
174
175			j1 = tempJets[0];
176			j2 = tempJets[1];
177			//
178			// additional jets
179			//
180			if (tempJets.size()>2){
181			for (unsigned int i=2 ; i< tempJets.size(); ++i){
182			addJets.push_back(tempJets[i]);
183			}
184			}
185			CompareJetPt ptComparator;
186
187			std::sort(addJets.begin(), addJets.end(), ptComparator);
188			return true;
189			}
190
191			void HbbCandidateFinderAlgo::findMuons(const std::vector<VHbbEvent::MuonInfo>& muons, std::vector<VHbbEvent::MuonInfo>& out){
192			/* Use:
193			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:
194
195			We use RECO Muons that are both Global and Tracker
196			chi2/ndof < 10 for the global muon fit
197			The track associated to the muon must have
198			>= 1 pixel hits
199			>= 10 pixel + strip hits
200			>= 1 valid hit in the muon chambers
201			>= 2 muon stations
202			\|dxy\| < 0.2
203			\|eta\| < 2.4
204			Relative combined isolation (R) is required to be < 0.15
205			R = [Sum pT(trks) + Et(ECAL) + Et(HCAL)] / pT(mu) computed in a cone of radius 0.3 in eta-phi
206			pT(mu) > 20 GeV
207			*/
208			for (std::vector<VHbbEvent::MuonInfo>::const_iterator it = muons.begin(); it!= muons.end(); ++it){
209			if (
210			(*it). globChi2<10 &&
211			(*it).nPixelHits>= 1 &&
212			(*it).globNHits >= 1 &&
213			(*it).nHits >= 10 &&
214			(*it).cat ==1 &&
215			(*it).validMuStations >=2 &&
216			(*it).ipDb<.2 &&
217			((it).hIso+(it).eIso+(it).tIso)/(it).fourMomentum.Pt()<.15 &&
218			(*it).fourMomentum.Eta()<2.4 &&
219			(*it).fourMomentum.Pt()>15) {
220			out.push_back(*it);
221			}
222			}
223			}
224
225
226			void HbbCandidateFinderAlgo::findElectrons(const std::vector<VHbbEvent::ElectronInfo>& electrons, std::vector<VHbbEvent::ElectronInfo>& out){
227			/*
228			We adopt the standard cut-based selection from VBTF described in detail here.
229
230			Z -> ee
231			gsf electrons
232			VBTF WP95
233			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
234			pT(e) > 20
235
236			W -> e nu
237			gsf electrons
238			VBTF WP80
239			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
240			pT(e) > 30
241			*/
242
243			for (std::vector<VHbbEvent::ElectronInfo>::const_iterator it = electrons.begin(); it!= electrons.end(); ++it){
244			if (
245			// fake
246			// (*it).id95> &&
247			fabs((*it).fourMomentum.Eta()) < 2.5 &&
248			!( fabs((it).fourMomentum.Eta()) < 1.57 && fabs((it).fourMomentum.Eta()) > 1.44) &&
249			(*it).fourMomentum.Pt()>15 // I use the minimum ok for both Z and W
250			){
251			out.push_back(*it);
252			}
253			}
254			}
255