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;
  if (useHighestPtHiggs_ == false){
    foundJets = findDiJets(event->simpleJets2,j1,j2,addJets) ;
  }else{
    foundJets= findDiJetsHighestPt(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 
  // removed on request by A. Rizzi - no filter at all
  //    if (met.metSig >2 ) out.push_back(met);
  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;


}


bool HbbCandidateFinderAlgo::findDiJetsHighestPt (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]);
 }

 if (tempJets.size()<2) return false;

 //loop over the dijets and save the one with highest Pt

 CompareJetPt ptComparator;
 std::sort(tempJets.begin(), tempJets.end(), ptComparator);
 //
 // so if i<j, pt(i)>pt(j)
 //
 
 float highestPt = -1000;
 unsigned  int highesti,highestj;
 for (unsigned int i =0; i< tempJets.size()-1; ++i){
   for (unsigned int j =i+1; j< tempJets.size(); ++j){
     float pt = (tempJets[i].fourMomentum+tempJets[j].fourMomentum).Pt();
     if (pt> highestPt){
       highestPt = pt;
       highesti=i;
       highestj=j;
     }
   }
 }
 j1 = tempJets[highesti]; 
 j2 = tempJets[highestj]; 

 for (unsigned int i=0; i<tempJets.size(); ++i){
   if (i!= highesti && i!= highestj)
     addJets.push_back(tempJets[i]);
 }

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

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