LeptonSelection/src/HZZCiCElectronSelection.cc

#include <iostream>
#include <math.h>
#include "HZZCiCElectronSelection.h"
#include "ParseArgs.h"
#include "SelectionStatus.h"

using namespace std;

CICStruct eidVeryLoose; 
CICStruct eidLoose;
CICStruct eidMedium;
CICStruct eidTight; 
CICStruct eidSuperTight;
CICStruct eidHyperTight1;

CICStruct getCiCCuts(TString eleSeleAnality) {
  if(eleSeleAnality=="veryloose")        return eidVeryLoose;
  if(eleSeleAnality=="loose")        return eidLoose;
  else if(eleSeleAnality=="medium")  return eidMedium;
  else if(eleSeleAnality=="tight")   return eidTight;
  else {cout << "error! bad taughtness." << endl; assert(0); return eidLoose;}
}

void initCiCSelection() {
#include "CiCElectronSelectionDataTuning.i"
};

//
//https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideCategoryBasedElectronID
//
unsigned cicCategory(const mithep::TElectron *ele) {
  unsigned categories=0;
  float eta = fabs(ele->scEta);
  float eop = fabs(ele->EoverP); // this is what it should be ...
  float fbrem = ele->fBrem;

  if( ele->isEB ) { 
    if( eop>0.9 && eop<1.2 && fbrem>=0.12 ) {
      return ELECTRON_BREM_EB;
    }
    else if( (eta >  .445 && eta <  .45  ) ||
             (eta >  .79  && eta <  .81  ) ||
             (eta > 1.137 && eta < 1.157 ) ||
             (eta > 1.47285 && eta < 1.4744) ) { 
      return ELECTRON_CRACK_EB;
    }
    else if ( !ele->isEcalDriven ) { // ele->trackerDrivenSeed &&
      return ELECTRON_PURETRACKER;
    } 
    else if( fbrem<0.12 ) {
      return ELECTRON_LOWBREM_EB;
    }
    else return ELECTRON_BADTRACK_EB;

  } else { // isEE

    if( eop>0.82 && eop<1.22 && fbrem>=0.2 ) {
      return ELECTRON_BREM_EE;
    }
    else if ( eta>1.5 && eta<1.58 ) {
      return ELECTRON_CRACK_EE;
    }
    else if ( !ele->isEcalDriven ) { //ele->trackerDrivenSeed &&
      return ELECTRON_PURETRACKER;
    } 
    else if( fbrem<0.2 ) {
      return ELECTRON_LOWBREM_EE;
    }
    else return ELECTRON_BADTRACK_EE;
  }
  return 0xDEADDEAD;
};

SelectionStatus electronPreSelectionCic(ControlFlags &ctrl, const mithep::TElectron *electron) { 

  bool ret=true;

  //
  // this is baseline
  //
  /*
  ret &= ( fabs(electron->ip3dSig) < 4 );
  if( ctrl.debug ) cout << "elepresel : ret after IPsig: " << ret <<  endl;
  ret &= ( electron->pt > 7 );
  if( ctrl.debug ) cout << "elepresel : ret after scPt: " << ret <<  endl;
  ret &= ( fabs(electron->scEta) < 2.5 );
  if( ctrl.debug ) cout << "elepresel : ret after scEta: " << ret <<  endl;
  */

  //
  // this is for UF Z->4l
  //
  ret &= ( electron->pt > 7 );
  if( ctrl.debug ) cout << "elepresel : ret after scPt: " << ret <<  endl;
  ret &= ( fabs(electron->scEta) < 2.5 );
  if( ctrl.debug ) cout << "elepresel : ret after scEta: " << ret <<  endl;
  ret &= ( electron->corrNExpHitsInner < 2 );
  if( ctrl.debug ) cout << "elepresel : ret after nexp: " << ret <<  endl;
  ret &= ( !electron->isConv );
  if( ctrl.debug ) cout << "elepresel : ret after isconv: " << ret <<  endl;
  ret &= ( fabs(electron->ip3dSig) < 100 );
  if( ctrl.debug ) cout << "elepresel : ret after IPsig: " << ret <<  endl;

  SelectionStatus status;
  if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
  return status;
}


SelectionStatus failsCicSelection( ControlFlags &ctrl, const mithep::TElectron * ele, CICStruct &cuts ) {

  int level=0;
  unsigned failmask=0x0;
  unsigned category = cicCategory(ele);

  float scEt      = ele->scEt;
  float scEta     = ele->scEta;

  float deta      = fabs(ele->deltaEtaIn);
  float dphi      = fabs(ele->deltaPhiIn);
  float sigieie   = ele->sigiEtaiEta;
  float EsoPin    = ele->ESeedClusterOverPIn;
  float fBrem     = ele->fBrem;
  float HoE       = ele->HoverE;
  float EsoPinCor =  EsoPin + fBrem; 
  // weird, see http://cmslxr.fnal.gov/lxr/source/RecoEgamma/ElectronIdentification/src/CutBasedElectronID.cc
  if( fBrem < 0 ) { 
    EsoPinCor = EsoPin;
  } 

  float dist = (ele->partnerDist == -9999.? 9999:ele->partnerDist);
  float dcot = (ele->partnerDeltaCot == -9999.? 9999:ele->partnerDeltaCot);
  float dcotdistcomb = ((0.04 - std::max(dist, dcot)) > 0?(0.04 - std::max(dist, dcot)):0);
  //  cout << "debug is now set to : " << ctrl.debug << endl; 

  if( ctrl.debug ) cout << "--> scEt : " << scEt << endl;

  float max_deta, max_dphi, max_HoE, max_sigieie, max_missHits, max_conv;
  float max_IP;
  float min_EoP; // this is a min cut ... maxn_EoP  


  float etmin=10.;
  float etmax=40.;

  if( scEt < etmin ) {
    max_deta     = cuts.cutdetainl[category];
    max_dphi     = cuts.cutdphiinl[category];
    max_HoE      = cuts.cuthoel[category];
    max_sigieie  = cuts.cutseel[category];
    min_EoP      = cuts.cuteseedopcor[category];
    max_missHits = cuts.cutfmishits[category];
    max_conv     = cuts.cutdcotdist[category];
  }
  else if( scEt > etmin && scEt < etmax ) {
    //    float dET = scEt-10;
    max_deta     = cuts.cutdetainl[category]   + (1/etmin - 1/scEt)*(cuts.cutdetain[category]-cuts.cutdetainl[category])/(1/etmin - 1/etmax);
    max_dphi     = cuts.cutdphiinl[category]   + (1/etmin - 1/scEt)*(cuts.cutdphiin[category]-cuts.cutdphiinl[category])/(1/etmin - 1/etmax);
    max_HoE      = cuts.cuthoel[category]      + (1/etmin - 1/scEt)*(cuts.cuthoe[category]-cuts.cuthoel[category])/(1/etmin - 1/etmax);
    max_sigieie  = cuts.cutseel[category]      + (1/etmin - 1/scEt)*(cuts.cutsee[category]-cuts.cutseel[category])/(1/etmin - 1/etmax);
    min_EoP      = cuts.cuteseedopcor[category];
    max_missHits = cuts.cutfmishits[category];
    max_conv     = cuts.cutdcotdist[category];
  } else {
    max_deta     = cuts.cutdetain[category];   
    max_dphi     = cuts.cutdphiin[category];   
    max_HoE      = cuts.cuthoe[category];      
    max_sigieie  = cuts.cutsee[category];      
    min_EoP      = cuts.cuteseedopcor[category];
    max_missHits = cuts.cutfmishits[category];
    max_conv     = cuts.cutdcotdist[category];
  }

  if( ctrl.debug ) cout << "max_deta: " << max_deta     
       << "\tmax_dphi    : " << max_dphi    
       << "\tmax_HoE     : "<< max_HoE     
       << "\tmax_sigieie : " << max_sigieie  
       << "\tmin_EoP     : " << min_EoP     
       << endl;

  level++;
  if( deta > max_deta ) { 
    failmask |= (1<<level);
    if( ctrl.debug ) cout << "failing deta ... " << deta << ">" << max_deta << endl;
  }
  level++;
  if( dphi > max_dphi ) { 
    failmask |= (1<<level);
    if(ctrl.debug) cout << "failing dphi ... " << endl;
  }
  level++;
  if( EsoPinCor < min_EoP ) { 
    failmask |= (1<<level);
    if(ctrl.debug) cout << "failing eop ... " << endl;
  }
  level++;
  if( sigieie > max_sigieie ) {
    failmask |= (1<<level);
    if(ctrl.debug) cout << "failing sigieie: " << "sigieie: " << sigieie << " max_sigieie: " << max_sigieie << endl;
  }
  level++;
  if( HoE > max_HoE ) {
    failmask |= (1<<level);
    if( ctrl.debug ) cout << "failing hoe ... " << HoE << ">" << max_HoE << endl;
  }
  level++;
  if( ele->nExpHitsInner > max_missHits  ) {
    //    || dcotdistcomb > max_conv 
    failmask |= (1<<level);
    if(ctrl.debug) cout << "failing misshits ... " << endl;
  }

  // loose trkIso
  level++;  
  if( (ele->trkIso03/ele->pt) > 0.7 ) {
    failmask |= (1<<level);
    if(ctrl.debug) cout << "failing loose iso ... " << endl;
  }

  //   level++;   
  //   if( fabs(ele->ip) > max_IP ) { 
  //     failmask |= (1<<level); 
  //   } 


  SelectionStatus ret; 

  if( failmask  ) { ret.setStatus(SelectionStatus::FAIL); return ret; }

  if(ctrl.debug) { 
    cout << "going to return a positive cicID ... " << endl;
    cout << "\tscheme: " << cuts.scheme << endl;
  }
  if( cuts.scheme == CICLoose ) ret.setStatus(SelectionStatus::LOOSEID); 
  //  if( cuts.scheme == CiCMedium ) ret.setStatus(SelectionStatus::TIGHTID); 
  else if( cuts.scheme == CICTight ) ret.setStatus(SelectionStatus::TIGHTID); 
  else { 
    cout << "returning  undefined " << endl;
    ret.setStatus(SelectionStatus::UNDEFINED);
  }
                           
                           
  return ret;
};

Revision:	1.2
Committed:	Fri Feb 17 14:49:00 2012 UTC (13 years, 3 months ago) by khahn
Content type:	text/plain
Branch:	MAIN
CVS Tags:	compiled, synced_FSR_2, synced_FSR, synched2, synched, HEAD
Changes since 1.1:	+58 -20 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	khahn	1.1	#include <iostream>
2			#include <math.h>
3			#include "HZZCiCElectronSelection.h"
4	khahn	1.2	#include "ParseArgs.h"
5			#include "SelectionStatus.h"
6	khahn	1.1
7			using namespace std;
8
9			CICStruct eidVeryLoose;
10			CICStruct eidLoose;
11			CICStruct eidMedium;
12			CICStruct eidTight;
13			CICStruct eidSuperTight;
14			CICStruct eidHyperTight1;
15
16			CICStruct getCiCCuts(TString eleSeleAnality) {
17			if(eleSeleAnality=="veryloose") return eidVeryLoose;
18			if(eleSeleAnality=="loose") return eidLoose;
19			else if(eleSeleAnality=="medium") return eidMedium;
20			else if(eleSeleAnality=="tight") return eidTight;
21			else {cout << "error! bad taughtness." << endl; assert(0); return eidLoose;}
22			}
23
24			void initCiCSelection() {
25			#include "CiCElectronSelectionDataTuning.i"
26			};
27
28			//
29			//https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideCategoryBasedElectronID
30			//
31			unsigned cicCategory(const mithep::TElectron *ele) {
32			unsigned categories=0;
33			float eta = fabs(ele->scEta);
34			float eop = fabs(ele->EoverP); // this is what it should be ...
35			float fbrem = ele->fBrem;
36
37			if( ele->isEB ) {
38			if( eop>0.9 && eop<1.2 && fbrem>=0.12 ) {
39			return ELECTRON_BREM_EB;
40			}
41			else if( (eta > .445 && eta < .45 ) \|\|
42			(eta > .79 && eta < .81 ) \|\|
43			(eta > 1.137 && eta < 1.157 ) \|\|
44			(eta > 1.47285 && eta < 1.4744) ) {
45			return ELECTRON_CRACK_EB;
46			}
47			else if ( !ele->isEcalDriven ) { // ele->trackerDrivenSeed &&
48			return ELECTRON_PURETRACKER;
49			}
50			else if( fbrem<0.12 ) {
51			return ELECTRON_LOWBREM_EB;
52			}
53			else return ELECTRON_BADTRACK_EB;
54
55			} else { // isEE
56
57			if( eop>0.82 && eop<1.22 && fbrem>=0.2 ) {
58			return ELECTRON_BREM_EE;
59			}
60			else if ( eta>1.5 && eta<1.58 ) {
61			return ELECTRON_CRACK_EE;
62			}
63			else if ( !ele->isEcalDriven ) { //ele->trackerDrivenSeed &&
64			return ELECTRON_PURETRACKER;
65			}
66			else if( fbrem<0.2 ) {
67			return ELECTRON_LOWBREM_EE;
68			}
69			else return ELECTRON_BADTRACK_EE;
70			}
71			return 0xDEADDEAD;
72			};
73
74	khahn	1.2	SelectionStatus electronPreSelectionCic(ControlFlags &ctrl, const mithep::TElectron *electron) {
75
76			bool ret=true;
77
78			//
79			// this is baseline
80			//
81			/*
82			ret &= ( fabs(electron->ip3dSig) < 4 );
83			if( ctrl.debug ) cout << "elepresel : ret after IPsig: " << ret << endl;
84			ret &= ( electron->pt > 7 );
85			if( ctrl.debug ) cout << "elepresel : ret after scPt: " << ret << endl;
86			ret &= ( fabs(electron->scEta) < 2.5 );
87			if( ctrl.debug ) cout << "elepresel : ret after scEta: " << ret << endl;
88			*/
89
90			//
91			// this is for UF Z->4l
92			//
93			ret &= ( electron->pt > 7 );
94			if( ctrl.debug ) cout << "elepresel : ret after scPt: " << ret << endl;
95			ret &= ( fabs(electron->scEta) < 2.5 );
96			if( ctrl.debug ) cout << "elepresel : ret after scEta: " << ret << endl;
97			ret &= ( electron->corrNExpHitsInner < 2 );
98			if( ctrl.debug ) cout << "elepresel : ret after nexp: " << ret << endl;
99			ret &= ( !electron->isConv );
100			if( ctrl.debug ) cout << "elepresel : ret after isconv: " << ret << endl;
101			ret &= ( fabs(electron->ip3dSig) < 100 );
102			if( ctrl.debug ) cout << "elepresel : ret after IPsig: " << ret << endl;
103
104			SelectionStatus status;
105			if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
106			return status;
107			}
108
109
110			SelectionStatus failsCicSelection( ControlFlags &ctrl, const mithep::TElectron * ele, CICStruct &cuts ) {
111	khahn	1.1
112			int level=0;
113			unsigned failmask=0x0;
114			unsigned category = cicCategory(ele);
115
116			float scEt = ele->scEt;
117			float scEta = ele->scEta;
118
119			float deta = fabs(ele->deltaEtaIn);
120			float dphi = fabs(ele->deltaPhiIn);
121			float sigieie = ele->sigiEtaiEta;
122			float EsoPin = ele->ESeedClusterOverPIn;
123			float fBrem = ele->fBrem;
124			float HoE = ele->HoverE;
125			float EsoPinCor = EsoPin + fBrem;
126			// weird, see http://cmslxr.fnal.gov/lxr/source/RecoEgamma/ElectronIdentification/src/CutBasedElectronID.cc
127			if( fBrem < 0 ) {
128			EsoPinCor = EsoPin;
129			}
130
131			float dist = (ele->partnerDist == -9999.? 9999:ele->partnerDist);
132			float dcot = (ele->partnerDeltaCot == -9999.? 9999:ele->partnerDeltaCot);
133			float dcotdistcomb = ((0.04 - std::max(dist, dcot)) > 0?(0.04 - std::max(dist, dcot)):0);
134			// cout << "debug is now set to : " << ctrl.debug << endl;
135
136			if( ctrl.debug ) cout << "--> scEt : " << scEt << endl;
137
138			float max_deta, max_dphi, max_HoE, max_sigieie, max_missHits, max_conv;
139			float max_IP;
140			float min_EoP; // this is a min cut ... maxn_EoP
141
142
143			float etmin=10.;
144			float etmax=40.;
145
146			if( scEt < etmin ) {
147			max_deta = cuts.cutdetainl[category];
148			max_dphi = cuts.cutdphiinl[category];
149			max_HoE = cuts.cuthoel[category];
150			max_sigieie = cuts.cutseel[category];
151			min_EoP = cuts.cuteseedopcor[category];
152			max_missHits = cuts.cutfmishits[category];
153			max_conv = cuts.cutdcotdist[category];
154			}
155			else if( scEt > etmin && scEt < etmax ) {
156			// float dET = scEt-10;
157			max_deta = cuts.cutdetainl[category] + (1/etmin - 1/scEt)*(cuts.cutdetain[category]-cuts.cutdetainl[category])/(1/etmin - 1/etmax);
158			max_dphi = cuts.cutdphiinl[category] + (1/etmin - 1/scEt)*(cuts.cutdphiin[category]-cuts.cutdphiinl[category])/(1/etmin - 1/etmax);
159			max_HoE = cuts.cuthoel[category] + (1/etmin - 1/scEt)*(cuts.cuthoe[category]-cuts.cuthoel[category])/(1/etmin - 1/etmax);
160			max_sigieie = cuts.cutseel[category] + (1/etmin - 1/scEt)*(cuts.cutsee[category]-cuts.cutseel[category])/(1/etmin - 1/etmax);
161			min_EoP = cuts.cuteseedopcor[category];
162			max_missHits = cuts.cutfmishits[category];
163			max_conv = cuts.cutdcotdist[category];
164			} else {
165			max_deta = cuts.cutdetain[category];
166			max_dphi = cuts.cutdphiin[category];
167			max_HoE = cuts.cuthoe[category];
168			max_sigieie = cuts.cutsee[category];
169			min_EoP = cuts.cuteseedopcor[category];
170			max_missHits = cuts.cutfmishits[category];
171			max_conv = cuts.cutdcotdist[category];
172			}
173
174			if( ctrl.debug ) cout << "max_deta: " << max_deta
175			<< "\tmax_dphi : " << max_dphi
176			<< "\tmax_HoE : "<< max_HoE
177			<< "\tmax_sigieie : " << max_sigieie
178			<< "\tmin_EoP : " << min_EoP
179			<< endl;
180
181			level++;
182			if( deta > max_deta ) {
183			failmask \|= (1<<level);
184			if( ctrl.debug ) cout << "failing deta ... " << deta << ">" << max_deta << endl;
185			}
186			level++;
187			if( dphi > max_dphi ) {
188			failmask \|= (1<<level);
189			if(ctrl.debug) cout << "failing dphi ... " << endl;
190			}
191			level++;
192			if( EsoPinCor < min_EoP ) {
193			failmask \|= (1<<level);
194			if(ctrl.debug) cout << "failing eop ... " << endl;
195			}
196			level++;
197			if( sigieie > max_sigieie ) {
198			failmask \|= (1<<level);
199			if(ctrl.debug) cout << "failing sigieie: " << "sigieie: " << sigieie << " max_sigieie: " << max_sigieie << endl;
200			}
201			level++;
202			if( HoE > max_HoE ) {
203			failmask \|= (1<<level);
204			if( ctrl.debug ) cout << "failing hoe ... " << HoE << ">" << max_HoE << endl;
205			}
206			level++;
207			if( ele->nExpHitsInner > max_missHits ) {
208			// \|\| dcotdistcomb > max_conv
209			failmask \|= (1<<level);
210			if(ctrl.debug) cout << "failing misshits ... " << endl;
211			}
212
213			// loose trkIso
214			level++;
215			if( (ele->trkIso03/ele->pt) > 0.7 ) {
216			failmask \|= (1<<level);
217			if(ctrl.debug) cout << "failing loose iso ... " << endl;
218			}
219
220			// level++;
221			// if( fabs(ele->ip) > max_IP ) {
222			// failmask \|= (1<<level);
223			// }
224
225	khahn	1.2
226			SelectionStatus ret;
227
228			if( failmask ) { ret.setStatus(SelectionStatus::FAIL); return ret; }
229
230			if(ctrl.debug) {
231			cout << "going to return a positive cicID ... " << endl;
232			cout << "\tscheme: " << cuts.scheme << endl;
233			}
234			if( cuts.scheme == CICLoose ) ret.setStatus(SelectionStatus::LOOSEID);
235			// if( cuts.scheme == CiCMedium ) ret.setStatus(SelectionStatus::TIGHTID);
236			else if( cuts.scheme == CICTight ) ret.setStatus(SelectionStatus::TIGHTID);
237			else {
238			cout << "returning undefined " << endl;
239			ret.setStatus(SelectionStatus::UNDEFINED);
240			}
241
242
243			return ret;
244	khahn	1.1	};
245