ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HbbAnalysis/src/HistosElecsBase.cc
Revision: 1.6
Committed: Wed Jun 9 14:34:17 2010 UTC (14 years, 11 months ago) by amagnan
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.5: +0 -0 lines
State: FILE REMOVED
Error occurred while calculating annotation data.
Log Message:
clean up histogram code, moved to analysis in UserCode/amagnan

File Contents

# Content
1 #include <iostream>
2 #include <fstream>
3 #include <sstream>
4 #include <cmath>
5 #include <algorithm>
6
7 #include "UserCode/HbbAnalysis/interface/HistosElecsBase.hh"
8
9 namespace HbbAnalysis {//namespace
10
11 void HistosElecsBase::Initialise(TFileDirectory & aDir, std::string aName, bool aDoGenMatched){
12
13 doGenMatched_ = aDoGenMatched;
14 CreateHistos(aName,aDir);
15
16 p_pTVeryLoose = aDir.make<TH1F>("p_pTVeryLoose",";p_{T}^{VeryLoose};N_{entries}",200,0,200);
17 p_etaVeryLoose = aDir.make<TH1F>("p_etaVeryLoose",";|#eta^{VeryLoose}|;N_{entries}",80,0,4);
18 p_pTLoose = aDir.make<TH1F>("p_pTLoose",";p_{T}^{Loose};N_{entries}",200,0,200);
19 p_etaLoose = aDir.make<TH1F>("p_etaLoose",";|#eta^{Loose}|;N_{entries}",80,0,4);
20 p_pTMedium = aDir.make<TH1F>("p_pTMedium",";p_{T}^{Medium};N_{entries}",200,0,200);
21 p_etaMedium = aDir.make<TH1F>("p_etaMedium",";|#eta^{Medium}|;N_{entries}",80,0,4);
22 p_pTTight = aDir.make<TH1F>("p_pTTight",";p_{T}^{Tight};N_{entries}",200,0,200);
23 p_etaTight = aDir.make<TH1F>("p_etaTight",";|#eta^{Tight}|;N_{entries}",80,0,4);
24 p_pTSuperTight = aDir.make<TH1F>("p_pTSuperTight",";p_{T}^{SuperTight};N_{entries}",200,0,200);
25 p_etaSuperTight = aDir.make<TH1F>("p_etaSuperTight",";|#eta^{SuperTight}|;N_{entries}",80,0,4);
26
27 p_abseta = aDir.make<TH1F>("p_abseta",";|#eta|;N_{entries}",80,0,4);
28
29 p_deltaRJet[0] = aDir.make<TH1F>("p_deltaRJetFlav0",";#Delta R;N_{entries}",100,0,5);
30 p_deltaRJet[1] = aDir.make<TH1F>("p_deltaRJetFlav1",";#Delta R;N_{entries}",100,0,5);
31 p_deltaRJet[2] = aDir.make<TH1F>("p_deltaRJetFlav2",";#Delta R;N_{entries}",100,0,5);
32 p_deltaRJet[3] = aDir.make<TH1F>("p_deltaRJetFlav3",";#Delta R;N_{entries}",100,0,5);
33
34 p_recoOverGen = aDir.make<TH1F>("p_recoOverGen",";p_{T}^{reco}/p_{T}^{gen};N_{entries}",300,0,3);
35
36 for (unsigned int bin(0); bin<20; bin++){//loop on eta bins
37 std::ostringstream retatitle;
38 retatitle << "p_recoOverGen_vsEta_" << bin ;
39 p_recoOverGen_vsEta[bin] = aDir.make<TH1F>(retatitle.str().c_str(),";p_{T}^{reco}/p_{T}^{gen};N_{entries}",300,0,3);
40 }//loop on eta bins
41
42 for (unsigned int bin(0); bin<10; bin++){//loop on pt bins
43 std::ostringstream rpttitle;
44 rpttitle << "p_recoOverGen_vsPt_" << bin;
45 p_recoOverGen_vsPt[bin] = aDir.make<TH1F>(rpttitle.str().c_str(),";p_{T}^{reco}/p_{T}^{gen};N_{entries}",300,0,3);
46 }//loop on pt bins
47
48
49 p_nElectrons = aDir.make<TH1F>("p_nElectrons",";N_{electrons};N_{entries}",10,0,10);
50
51 p_electronID = aDir.make<TH1F>("p_electronID",";electronID;N_{entries}",20,0,20);
52
53 p_scSigmaEtaEta = aDir.make<TH1F>("p_scSigmaEtaEta",";#sigma_{#eta,#eta};N_{entries}",500,0,0.1);
54 p_scSigmaIEtaIEta = aDir.make<TH1F>("p_scSigmaIEtaIEta",";#sigma_{i#eta,i#eta};N_{entries}",500,0,0.1);
55 p_scE1x5 = aDir.make<TH1F>("p_scE1x5",";E_{sc}^{1#times 5} (GeV);N_{entries}",500,0,500);
56 p_scE2x5Max = aDir.make<TH1F>("p_scE2x5Max",";E_{sc}^{2#times 5} (GeV);N_{entries}",500,0,500);
57 p_scE5x5 = aDir.make<TH1F>("p_scE5x5",";E_{sc}^{5#times 5} (GeV);N_{entries}",500,0,500);
58 p_scE1x5OverscE5x5 = aDir.make<TH1F>("p_scE1x5OverscE5x5",";E_{sc}^{1#times 5}/E_{sc}^{5#times 5};N_{entries}",500,0,2);
59 p_scE2x5MaxOverscE5x5 = aDir.make<TH1F>("p_scE2x5MaxOverscE5x5",";E_{sc}^{2#times 5}/E_{sc}^{5#times 5};N_{entries}",500,0,2);
60
61 p_eSuperClusterOverP = aDir.make<TH1F>("p_eSuperClusterOverP",";E_{sc}/p_{track};N_{entries}",100,0.5,2);
62
63 p_HoverE = aDir.make<TH1F>("p_HoverE",";E_{had}/E_{em};N_{entries}",150,0,0.15);
64 p_deltaPhiIn = aDir.make<TH1F>("p_deltaPhiIn",";#Delta#phi_{in};N_{entries}",150,0,0.15);
65 p_deltaEtaIn = aDir.make<TH1F>("p_deltaEtaIn",";#Delta#eta_{in};N_{entries}",150,0,0.03);
66
67 p_caloIso = aDir.make<TH1F>("p_caloIso",";caloIso;N_{entries}",400,0,20);
68 p_hcalIso = aDir.make<TH1F>("p_hcalIso",";hcalIso;N_{entries}",400,0,20);
69
70 p_trackIso = aDir.make<TH1F>("p_trackIso",";trackIso (GeV);N_{entries}",400,0,20);
71 p_trackIsoOverEt = aDir.make<TH1F>("p_trackIsoOverEt",";trackIso/E_{T};N_{entries}",400,0,1);
72 p_ecalIso = aDir.make<TH1F>("p_ecalIso",";ecalIso (GeV);N_{entries}",400,0,20);
73 p_ecalIsoOverEt = aDir.make<TH1F>("p_ecalIsoOverEt",";ecalIso/E_{T};N_{entries}",400,0,1);
74 p_combIso = aDir.make<TH1F>("p_combIso",";combIso (GeV);N_{entries}",400,0,40);
75 p_combIsoOverEt = aDir.make<TH1F>("p_combIsoOverEt",";combIso/E_{T};N_{entries}",400,0,1);
76
77 p_gsfTrk_pT = aDir.make<TH1F>("p_gsfTrk_pT",";p_{T}^{trk} (GeV);N_{entries}",200,0,200);
78 p_gsfTrk_IPxy = aDir.make<TH1F>("p_gsfTrk_IPxy",";IP_{xy}^{trk};N_{entries}",100,-0.1,0.1);
79 p_gsfTrk_IPz = aDir.make<TH1F>("p_gsfTrk_IPz",";IP_{z}^{trk};N_{entries}",100,-1,1);
80
81 //efficiency curves vs cut value
82 peff_eIso[0] = aDir.make<TH1F>("peff_eIso_trkIso",";trkIso (GeV);N_{sel}/N_{tot}",50,0,5);
83 peff_eIso[1] = aDir.make<TH1F>("peff_eIso_trkIsoOverEt",";trkIso/E_{T};N_{sel}/N_{tot}",50,0,1);
84 peff_eIso[2] = aDir.make<TH1F>("peff_eIso_ecalIso",";ecalIso (GeV);N_{sel}/N_{tot}",50,0,5);
85 peff_eIso[3] = aDir.make<TH1F>("peff_eIso_ecalIsoOverEt",";ecalIso/E_{T};N_{sel}/N_{tot}",50,0,1);
86 peff_eIso[4] = aDir.make<TH1F>("peff_eIso_combIso",";combIso (GeV);N_{sel}/N_{tot}",50,0,10);
87 peff_eIso[5] = aDir.make<TH1F>("peff_eIso_combIsoOverEt",";combIso/E_{T};N_{sel}/N_{tot}",50,0,1);
88
89 //efficiency curves vs pT muon
90 peff_eEt[0] = aDir.make<TH1F>("peff_eEt_trkIso",";E_{T} (GeV);N_{sel}/N_{tot}",18,10,100);
91 const unsigned int lNBinsPt = peff_eEt[0]->GetNbinsX();
92 const double lMinPt = peff_eEt[0]->GetXaxis()->GetXmin();
93 const double lMaxPt = peff_eEt[0]->GetXaxis()->GetXmax();
94 peff_eEt[1] = aDir.make<TH1F>("peff_eEt_trkIsoOverEt",";E_{T} (GeV);N_{sel}/N_{tot}",lNBinsPt,lMinPt,lMaxPt);
95 peff_eEt[2] = aDir.make<TH1F>("peff_eEt_ecalIso",";E_{T} (GeV);N_{sel}/N_{tot}",lNBinsPt,lMinPt,lMaxPt);
96 peff_eEt[3] = aDir.make<TH1F>("peff_eEt_ecalIsoOverEt",";E_{T} (GeV);N_{sel}/N_{tot}",lNBinsPt,lMinPt,lMaxPt);
97 peff_eEt[4] = aDir.make<TH1F>("peff_eEt_combIso",";E_{T} (GeV);N_{sel}/N_{tot}",lNBinsPt,lMinPt,lMaxPt);
98 peff_eEt[5] = aDir.make<TH1F>("peff_eEt_combIsoOverEt",";E_{T} (GeV);N_{sel}/N_{tot}",lNBinsPt,lMinPt,lMaxPt);
99
100
101 //efficiency curves vs eta muon
102 peff_eEta[0] = aDir.make<TH1F>("peff_eEta_trkIso",";#eta;N_{sel}/N_{tot}",42,-2.1,2.1);
103 const unsigned int lNBinsEta = peff_eEta[0]->GetNbinsX();
104 const double lMinEta = peff_eEta[0]->GetXaxis()->GetXmin();
105 const double lMaxEta = peff_eEta[0]->GetXaxis()->GetXmax();
106 peff_eEta[1] = aDir.make<TH1F>("peff_eEta_trkIsoOverEt",";#eta;N_{sel}/N_{tot}",lNBinsEta,lMinEta,lMaxEta);
107 peff_eEta[2] = aDir.make<TH1F>("peff_eEta_ecalIso",";#eta;N_{sel}/N_{tot}",lNBinsEta,lMinEta,lMaxEta);
108 peff_eEta[3] = aDir.make<TH1F>("peff_eEta_ecalIsoOverEt",";#eta;N_{sel}/N_{tot}",lNBinsEta,lMinEta,lMaxEta);
109 peff_eEta[4] = aDir.make<TH1F>("peff_eEta_combIso",";#eta;N_{sel}/N_{tot}",lNBinsEta,lMinEta,lMaxEta);
110 peff_eEta[5] = aDir.make<TH1F>("peff_eEta_combIsoOverEt",";#eta;N_{sel}/N_{tot}",lNBinsEta,lMinEta,lMaxEta);
111
112 //efficiency calculation tool
113 for (unsigned int i(0); i<6; i++){
114 const unsigned int lNBins = peff_eIso[i]->GetNbinsX();
115 const double lMin = peff_eIso[i]->GetXaxis()->GetXmin();
116 const double lMax = peff_eIso[i]->GetXaxis()->GetXmax();
117 isoEff_[i].initialise(lNBins,lMin,lMax);
118 isoEffEt_[i].initialise(lNBinsPt,lMinPt,lMaxPt);
119 isoEffEta_[i].initialise(lNBinsEta,lMinEta,lMaxEta);
120 }
121
122
123 }
124
125 void HistosElecsBase::FillEventHistograms(const std::vector<HbbAnalysis::Electron> & aCol){
126
127 if (aCol.size() > 0) {
128 if (doGenMatched_) {
129 unsigned int nMatched = 0;
130 for (std::vector<HbbAnalysis::Electron>::const_iterator iElec = aCol.begin();
131 iElec != aCol.end();
132 iElec++)
133 {
134 if (MatchesGenElectron(*iElec)) nMatched++;
135 }
136 p_nElectrons->Fill(nMatched);
137 }
138 else p_nElectrons->Fill(aCol.size());
139
140 }
141 }
142
143 void HistosElecsBase::FillHistograms(const HbbAnalysis::Electron & aElec, const std::vector<float> & dR)
144 {
145 bool lIsGenMatched = !doGenMatched_ || (doGenMatched_ && MatchesGenElectron(aElec));
146 if (lIsGenMatched) {//genMatched
147 for (unsigned int lFlav(0); lFlav<dR.size(); lFlav++){
148 p_deltaRJet[lFlav]->Fill(dR[lFlav]);
149 }
150 }
151 }
152
153 void HistosElecsBase::FillHistograms(const HbbAnalysis::Electron & aElec, bool isLead){//FillHistograms
154
155 bool lIsGenMatched = !doGenMatched_ || (doGenMatched_ && MatchesGenElectron(aElec));
156 if (lIsGenMatched) {//genMatched
157 if (isLead) {//isLead
158 FillBaseHistograms(aElec.recoVars().pT,aElec.recoVars().eta,aElec.recoVars().phi,aElec.recoVars().charge);
159
160 // if ( aElec.gsfTrack().isAvailable() && !aElec.gsfTrack().isNull() ) {
161 // p_gsfTrk_pT->Fill(aElec.gsfTrack()->pt());
162 // if ( aRecoVertices->size() >= 1 ) {
163 // const reco::Vertex & thePrimaryEventVertex = (*aRecoVertices->begin());
164 // p_gsfTrk_IPxy->Fill(aElec.gsfTrack()->dxy(thePrimaryEventVertex.position()));
165 // p_gsfTrk_IPz->Fill(aElec.gsfTrack()->dz(thePrimaryEventVertex.position()));
166 // }
167 // }
168
169
170 if (aElec.genVars().valid && fabs(aElec.genVars().eta)<2.1 && aElec.genVars().pT>10) {
171 int binEta = static_cast<int>((aElec.genVars().eta+2.1)/0.21);
172 int binpt;
173 if (aElec.genVars().pT < 100) binpt = static_cast<int>((aElec.genVars().pT-10)/11.25);
174 else if (aElec.genVars().pT < 120) binpt = 8;
175 else binpt = 9;
176
177 p_recoOverGen_vsEta[binEta]->Fill(aElec.recoVars().pT/aElec.genVars().pT);
178 p_recoOverGen_vsPt[binpt]->Fill(aElec.recoVars().pT/aElec.genVars().pT);
179 p_recoOverGen->Fill(aElec.recoVars().pT/aElec.genVars().pT);
180 }
181
182
183 std::vector<std::pair<std::string,float> > lIds = aElec.idVars().electronIDs;
184 if (lIds.size() > 5) std::cout << "--- WARNING: histo will be out-of-range, number of electron IDs = " << lIds.size() << std::endl;
185 for (unsigned int iId(0); iId<lIds.size(); iId++){
186 float lVal = lIds.at(iId).second;
187 //const std::string & lName = lIds.at(iId).first;
188 //std::cout << "--- ElectronID : Bin " << iId << ", id name = " << lName << ", id val = " << lVal << std::endl;
189 //if (lIds.at(iId).second > 1) std::cout << "--- ElectronID : Bin " << iId << ", id name = " << lName << ", id val = " << lVal << std::endl;
190 //val = 0,1(id),2(iso),3(id+iso) in new ids
191 p_electronID->Fill(4*iId+lVal);
192
193 if (iId==0 && lVal>0.5) {
194 p_pTLoose->Fill(aElec.recoVars().pT);
195 p_etaLoose->Fill(fabs(aElec.recoVars().eta));
196 }
197 if (iId==1 && lVal>0.5) {
198 p_pTSuperTight->Fill(aElec.recoVars().pT);
199 p_etaSuperTight->Fill(fabs(aElec.recoVars().eta));
200 }
201 if (iId==2 && lVal>0.5) {
202 p_pTVeryLoose->Fill(aElec.recoVars().pT);
203 p_etaVeryLoose->Fill(fabs(aElec.recoVars().eta));
204 }
205 if (iId==3 && lVal>0.5) {
206 p_pTMedium->Fill(aElec.recoVars().pT);
207 p_etaMedium->Fill(fabs(aElec.recoVars().eta));
208 }
209 if (iId==4 && lVal>0.5) {
210 p_pTTight->Fill(aElec.recoVars().pT);
211 p_etaTight->Fill(fabs(aElec.recoVars().eta));
212 }
213
214
215 }
216
217 p_abseta->Fill(fabs(aElec.recoVars().eta));
218
219 p_scSigmaEtaEta->Fill(aElec.scVars().sigmaEtaEta);
220 p_scSigmaIEtaIEta->Fill(aElec.scVars().sigmaIEtaIEta);
221 p_scE1x5->Fill(aElec.scVars().e1x5);
222 p_scE2x5Max->Fill(aElec.scVars().e2x5Max);
223 p_scE5x5->Fill(aElec.scVars().e5x5);
224
225 p_eSuperClusterOverP->Fill(aElec.scVars().eOverP);
226
227 p_HoverE->Fill(aElec.idVars().hOverE);
228 p_deltaPhiIn->Fill(aElec.idVars().deltaPhiIn);
229 p_deltaEtaIn->Fill(aElec.idVars().deltaEtaIn);
230
231 if (aElec.scVars().e5x5 > 0) {
232 p_scE1x5OverscE5x5->Fill(aElec.scVars().e1x5/aElec.scVars().e5x5);
233 p_scE2x5MaxOverscE5x5->Fill(aElec.scVars().e2x5Max/aElec.scVars().e5x5);
234 }
235
236 p_caloIso->Fill(aElec.isoVars().calo);
237 p_hcalIso->Fill(aElec.isoVars().hcal);
238
239 double lIsoVar[6] = {
240 aElec.isoVars().track,
241 aElec.isoVars().track/aElec.recoVars().pT,
242 aElec.isoVars().ecal,
243 aElec.isoVars().ecal/aElec.recoVars().pT,
244 aElec.isoVars().track+aElec.isoVars().ecal+aElec.isoVars().hcal,
245 (aElec.isoVars().track+aElec.isoVars().ecal+aElec.isoVars().hcal)/aElec.recoVars().pT
246 };
247
248 p_trackIso->Fill(lIsoVar[0]);
249 p_trackIsoOverEt->Fill(lIsoVar[1]);
250 p_ecalIso->Fill(lIsoVar[2]);
251 p_ecalIsoOverEt->Fill(lIsoVar[3]);
252 p_combIso->Fill(lIsoVar[4]);
253 p_combIsoOverEt->Fill(lIsoVar[5]);
254
255 double lIsoCut[6] = {2.4,0.05,2.4,0.09,3.,0.1};
256
257 for (unsigned int i(0); i<6; i++){//loop on iso variables
258
259 for (unsigned int lBin(0); lBin<isoEff_[i].numberOfBins(); lBin++){
260 double lCut = isoEff_[i].xMin()+isoEff_[i].stepSize()*lBin;
261
262 isoEff_[i].incrementTotal(lBin);
263 if (lIsoVar[i] <= lCut){
264 isoEff_[i].incrementPass(lBin);
265 }
266 }
267
268 //eff vs pT
269 for (unsigned int lBin(0); lBin<isoEffEt_[i].numberOfBins(); lBin++){
270 double lCutMin = isoEffEt_[i].xMin()+isoEffEt_[i].stepSize()*lBin;
271 double lCutMax = isoEffEt_[i].xMin()+isoEffEt_[i].stepSize()*(lBin+1);
272
273 if (aElec.recoVars().pT >= lCutMin && aElec.recoVars().pT < lCutMax){
274 isoEffEt_[i].incrementTotal(lBin);
275
276 if (lIsoVar[i] <= lIsoCut[i]){
277 isoEffEt_[i].incrementPass(lBin);
278 }
279 }
280 }
281
282 //eff vs eta
283 for (unsigned int lBin(0); lBin<isoEffEta_[i].numberOfBins(); lBin++){
284 double lCutMin = isoEffEta_[i].xMin()+isoEffEta_[i].stepSize()*lBin;
285 double lCutMax = isoEffEta_[i].xMin()+isoEffEta_[i].stepSize()*(lBin+1);
286
287 if (aElec.recoVars().eta >= lCutMin && aElec.recoVars().eta < lCutMax){
288 isoEffEta_[i].incrementTotal(lBin);
289
290 if (lIsoVar[i] <= lIsoCut[i]){
291 isoEffEta_[i].incrementPass(lBin);
292 }
293 }
294 }
295
296
297 }//loop on iso variables
298
299
300
301 }//isLead
302 }//genMatched
303
304 }//FillHistograms
305
306
307
308 void HistosElecsBase::FillEffHistograms(){//FillEffHistograms
309
310 for (unsigned int i(0); i<6; i++){//loop on iso variables
311
312 for (unsigned int lBin(0); lBin<isoEff_[i].numberOfBins(); lBin++){
313 peff_eIso[i]->SetBinContent(lBin+1,isoEff_[i].getRatio(lBin));
314 peff_eIso[i]->SetBinError(lBin+1,isoEff_[i].getRatioError(lBin));
315 }
316 //eff vs pT
317 for (unsigned int lBin(0); lBin<isoEffEt_[i].numberOfBins(); lBin++){
318 peff_eEt[i]->SetBinContent(lBin+1,isoEffEt_[i].getRatio(lBin));
319 peff_eEt[i]->SetBinError(lBin+1,isoEffEt_[i].getRatioError(lBin));
320 }
321
322 //eff vs eta
323 for (unsigned int lBin(0); lBin<isoEffEta_[i].numberOfBins(); lBin++){
324 peff_eEta[i]->SetBinContent(lBin+1,isoEffEta_[i].getRatio(lBin));
325 peff_eEta[i]->SetBinError(lBin+1,isoEffEta_[i].getRatioError(lBin));
326 }
327
328 }//loop on iso variables
329
330 }//FillEffHistograms
331
332 bool HistosElecsBase::MatchesGenElectron(const HbbAnalysis::Electron& aElec)
333 {
334 return (aElec.genVars().valid && abs(aElec.genVars().pdgId) == 11);
335 }
336
337 }//namespace
338
339