VHbbDataFormats/src/HbbCandidateFinderAlgo.cc

#include "VHbbAnalysis/VHbbDataFormats/interface/HbbCandidateFinderAlgo.h"
#include "VHbbAnalysis/VHbbDataFormats/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;
  }
};


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

  VHbbCandidateTools selector(verbose_);

  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;
  TVector3 higgsBoost;
  higgsBoost = (temp.H.fourMomentum).BoostVector();
  temp.H.helicities.push_back(selector.getHelicity(j1,higgsBoost));
  temp.H.helicities.push_back(selector.getHelicity(j2,higgsBoost));
  temp.H.deltaTheta = selector.getDeltaTheta(j1,j2);
  temp.additionalJets = addJets;
  temp.V.mets = met;
  temp.V.muons = mu;
  temp.V.electrons = ele;
  
  //
  // now see which kind of andidate this can be
  // 

  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;
 CompareJetPt ptComparator;
 
 if (verbose_){
   std::cout <<" CandidateFinder: Intermediate Jets = "<<tempJets.size()<<std::endl;
 }


 if (tempJets.size()<2) return false;
 
 std::sort(tempJets.begin(), tempJets.end(), bTagComparator);
 
if (tempJets[0].fourMomentum.Pt()>(tempJets[1].fourMomentum.Pt())){
   j1 = tempJets[0];
   j2 = tempJets[1];
 }else{
   j2 = tempJets[0];
   j1 = tempJets[1]; 
 }

//
 // additional jets
 //
 if (tempJets.size()>2){
   for (unsigned int i=2 ; i< tempJets.size(); ++i){
     addJets.push_back(tempJets[i]);
   }
 }

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