avgefficiency.C
Go to the documentation of this file.
1 // Macro for making plots of average cell efficiencies by view
2 // Runs on the output histogram file from HitEfficiency_module.cc,
3 // though you will want to either TChain or hadd a bunch of these
4 // files together to get any sort of reasonable statistics
5 
6 #include "TCanvas.h"
7 #include "TFile.h"
8 #include "TGraph.h"
9 #include "TH1.h"
10 #include "TH2.h"
11 #include "TMath.h"
12 #include "TLegend.h"
13 #include "TProfile.h"
14 
15 #include <cmath>
16 #include <string>
17 #include <iostream>
18 
19 
21 {
22 
23  double nmincell = 500; // mininum number of entries in denominator of cell efficiency for a cell to be "high" stat
24 
25  // Get the ana histo file and grab the relevant histograms
26  TFile* tfile = new TFile("nova/ana/users/raddatz/raddatz/LowThreshold/haddlowthres.root","READ");
27 
28  // numerator of efficiency vs w (lencut histos are for path lengths 6-7 cm)
29  TH2F* hitsVwcell = (TH2F*)tfile->Get("hiteff/hitsVwcelllencut");
30  // denominator of efficiency vs w (lencut histos are for path lengths 6-7 cm)
31  TH2F* cellsVwcell = (TH2F*)tfile->Get("hiteff/cellsVwcelllencut");
32  // numerator of efficiency vs path length vs w for x cells
33  TH2F* hitsVwVlenx = (TH2F*)tfile->Get("hiteff/hitsVwVlenx");
34  // denominator of efficiency vs path length vs w for x cells
35  TH2F* cellsVwVlenx = (TH2F*)tfile->Get("hiteff/cellsVwVlenx");
36  // numerator of efficiency vs path length vs w for y cells
37  TH2F* hitsVwVleny = (TH2F*)tfile->Get("hiteff/hitsVwVleny");
38  // denominator of efficiency vs path length vs w for y cells
39  TH2F* cellsVwVleny = (TH2F*)tfile->Get("hiteff/cellsVwVleny");
40 
41  // Define a bunch of histogram binning
42  // First grab predefined binning from histograms
43  int ncellbins = hitsVwcell->GetNbinsX();
44  double cellmax = hitsVwcell->GetXaxis()->GetXmax();
45  int nwbins = hitsVwcell->GetNbinsY();
46  double wmin = hitsVwcell->GetYaxis()->GetXmin();
47  double wmax = hitsVwcell->GetYaxis()->GetXmax();
48  int npathbins = hitsVwVlenx->GetNbinsY();
49  double pathmin = hitsVwVlenx->GetYaxis()->GetXmin();
50  double pathmax = hitsVwVlenx->GetYaxis()->GetXmax();
51  int neffbins = 101;
52  double effmin = -0.05;
53  double effmax = 1.05;
54  double plmin = -0.5;
55  int cellmaxplanecell = ((int)cellmax-1)%384;
56  double plmax = ((int)cellmax-1-cellmaxplanecell)/384 + 0.5;
57  int nplbins = (int)(plmax-plmin);
58  int nclbins = 385;
59  double clmin = -0.5;
60  double clmax = 384.5;
61 
62 
63  // Define a bunch of histograms
64  // efficiency as a function of path length and w histograms as well as ratio between views
65  TH2F* effVwVlenx = new TH2F("effVwVlenx",";W (cm); Path Length (cm)",nwbins,wmin,wmax,npathbins,pathmin,pathmax);
66  TH2F* effVwVleny = new TH2F("effVwVleny",";W (cm); Path Length (cm)",nwbins,wmin,wmax,npathbins,pathmin,pathmax);
67  TH2F* effVwVlenratio = new TH2F("effVwVlenratio",";W (cm);Path Length (cm)",
68  effVwVlenx->GetNbinsX(),effVwVlenx->GetXaxis()->GetXmin(),effVwVlenx->GetXaxis()->GetXmax(),
69  effVwVlenx->GetNbinsY(),effVwVlenx->GetYaxis()->GetXmin(),effVwVlenx->GetYaxis()->GetXmax());
70 
71  // efficiency vs. w averaged over cells for each view
72  TH1F* xeffVw = new TH1F("xeffVw",";W (cm);Efficiency",nwbins,wmin,wmax);
73  TH1F* yeffVw = new TH1F("yeffVw",";W (cm);Efficiency",nwbins,wmin,wmax);
74 
75  // efficiency vs w averaged over cells for each view using only cells with high stats
76  TH1F* xeffVwcut = new TH1F("xeffVwcut",";W (cm);Efficiency",nwbins,wmin,wmax);
77  TH1F* yeffVwcut = new TH1F("yeffVwcut",";W (cm);Efficiency",nwbins,wmin,wmax);
78 
79  // average efficiency Vs Cell Vs Plane for the lowest third, mid third, and highest third of the cell for x cells
80  TH2F* effslowCellVPlanex = new TH2F("effslowCellVPlanex",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
81  TH2F* effsmidCellVPlanex = new TH2F("effsmidCellVPlanex",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
82  TH2F* effshighCellVPlanex = new TH2F("effshighCellVPlanex",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
83 
84  // averaged efficiency of the lowest third, mid third, and highest third of the cell for x cells
85  TH1F* effslowx = new TH1F("effslowx",";Efficiency;NCells",neffbins,effmin,effmax);
86  TH1F* effsmidx = new TH1F("effsmidx",";Efficiency;NCells",neffbins,effmin,effmax);
87  TH1F* effshighx = new TH1F("effshighx",";Efficiency;NCells",neffbins,effmin,effmax);
88 
89  // average efficiency Vs Cell Vs Plane for the lowest third, mid third, and highest third of the cell for y cells
90  TH2F* effslowCellVPlaney = new TH2F("effslowCellVPlaney",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
91  TH2F* effsmidCellVPlaney = new TH2F("effsmidCellVPlaney",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
92  TH2F* effshighCellVPlaney = new TH2F("effshighCellVPlaney",";Plane;Cell",nplbins,plmin,plmax,nclbins,clmin,clmax);
93 
94  // averaged efficiency of the lowest third, mid third, and highest third of the cell for y cells
95  TH1F* effslowy = new TH1F("effslowy",";Efficiency;NCells",neffbins,effmin,effmax);
96  TH1F* effsmidy = new TH1F("effsmidy",";Efficiency;NCells",neffbins,effmin,effmax);
97  TH1F* effshighy = new TH1F("effshighy",";Efficiency;NCells",neffbins,effmin,effmax);
98 
99  // helper histogram to get the average efficiency of a cell over varius w ranges
100  TH1F* eff = new TH1F("eff","",nwbins,wmin,wmax);
101 
102  int lowwbin = (int)((double)nwbins/3.0);
103  int highwbin = (int)(2.0*(double)nwbins/3.0);
104  // helper histograms for getting the average and variation in average cell efficiencies by view
105  TH2F* effVwx = new TH2F("effVwx",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);
106  TH2F* effVwy = new TH2F("effVwy",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);
107  TH2F* effVwxcut = new TH2F("effVwxcut",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);
108  TH2F* effVwycut = new TH2F("effVwycut",";W (cm);Efficiency",nwbins,wmin,wmax,neffbins,effmin,effmax);
109 
110 
111  TCanvas* c1 = new TCanvas("c1","c1");
112  effVwVlenx->Divide(hitsVwVlenx,cellsVwVlenx);
113  effVwVlenx->Draw("colz");
114  effVwVlenx->GetZaxis()->SetRangeUser(0.5,1.0);
115 
116  TCanvas* c2 = new TCanvas("c2","c2");
117  effVwVleny->Divide(hitsVwVleny,cellsVwVleny);
118  effVwVleny->Draw("colz");
119  effVwVleny->GetZaxis()->SetRangeUser(0.5,1.0);
120 
121  TCanvas* c3 = new TCanvas("c3","c3");
122  effVwVlenratio->Divide(effVwVleny,effVwVlenx);
123  effVwVlenratio->Draw("colz");
124  effVwVlenratio->GetZaxis()->SetRangeUser(0.8,1.2);
125 
126 
127  double nxcell = 0.0;
128  double nycell = 0.0;
129 
130  double nxcellcut = 0.0;
131  double nycellcut = 0.0;
132 
133  // loop over every cell in the first diblock
134  for(int icell = 1; icell <= ncellbins; ++icell){
135  // convert this to a plane and cell number
136  int cell = (icell-1)%384;
137  int plane = (icell-1-cell)/384;
138  int modcell = cell%32;
139 // int module = (cell-modcell)/32;
140 
141  int view = plane%2; // view = 0 y, view = 1 x
142  // need to switch the module in cell for x view to keep consistent cell numbering within module between views
143  if(view){ modcell = 31-modcell; }
144 
145  int celltot = cellsVwcell->ProjectionY("",icell,icell,"")->GetEntries();
146  // If never get any tri-cells on this cell skip it
147  if(celltot == 0){ continue; }
148  if(plane > (int)plmax){ continue; }
149 
150  if(view){ ++nxcell; }
151  else{ ++nycell; }
152  if(celltot > nmincell){
153  if(view){ ++nxcellcut; }
154  else{ ++nycellcut; }
155  }
156  // loop over every w bin
157  for(int iw = 1; iw <= nwbins; ++iw){
158  // get the efficiency at this cell/w, eff = k/n
159  double k = hitsVwcell->GetBinContent(icell,iw);
160  double n = cellsVwcell->GetBinContent(icell,iw);
161  double ef = k/n;
162  if(n == 0){ ef = 0; }
163  // calculate the variance of this
164  double var = (k+1)*(k+2)/((n+2)*(n+3)) - (k+1)*(k+1)/((n+2)*(n+2));
165 
166  double w = hitsVwcell->GetYaxis()->GetBinCenter(iw);
167 
168  eff->SetBinContent(iw,ef);
169  if(view){
170  xeffVw->SetBinContent(iw,(ef+xeffVw->GetBinContent(iw)));
171  effVwx->Fill(w,ef);
172  if(celltot > nmincell){
173  xeffVwcut->SetBinContent(iw,(ef+xeffVwcut->GetBinContent(iw)));
174  effVwxcut->Fill(w,ef);
175  }
176  }
177  else{
178  yeffVw->SetBinContent(iw,(ef+yeffVw->GetBinContent(iw)));
179  effVwy->Fill(w,ef);
180  if(celltot > nmincell){
181  yeffVwcut->SetBinContent(iw,(ef+yeffVwcut->GetBinContent(iw)));
182  effVwycut->Fill(w,ef);
183  }
184  }
185 
186  eff->SetBinError(iw,sqrt(var));
187  } // end loop over w bins
188 
189  // now get the average efficiency of this cell over different w ranges
190  double lowedgeeff = eff->Integral(1,lowwbin)/((double)(lowwbin - 1 + 1));
191  double mideff = eff->Integral(lowwbin+1,highwbin-1)/((double)(highwbin-1 - (lowwbin+1) + 1));
192  double highedgeeff = eff->Integral(highwbin,nwbins)/((double)(nwbins - (highwbin) + 1));
193 
194  int planebin = effslowCellVPlanex->GetXaxis()->FindFixBin(plane);
195  int cellbin = effslowCellVPlanex->GetYaxis()->FindFixBin(cell);
196 
197  if(view){
198  effslowCellVPlanex->SetBinContent(planebin,cellbin,lowedgeeff);
199  effslowx->Fill(lowedgeeff);
200  effsmidCellVPlanex->SetBinContent(planebin,cellbin,mideff);
201  effsmidx->Fill(mideff);
202  effshighCellVPlanex->SetBinContent(planebin,cellbin,highedgeeff);
203  effshighx->Fill(highedgeeff);
204  }
205  else{
206  effslowCellVPlaney->SetBinContent(planebin,cellbin,lowedgeeff);
207  effslowy->Fill(lowedgeeff);
208  effsmidCellVPlaney->SetBinContent(planebin,cellbin,mideff);
209  effsmidy->Fill(mideff);
210  effshighCellVPlaney->SetBinContent(planebin,cellbin,highedgeeff);
211  effshighy->Fill(highedgeeff);
212  }
213 
214  }
215 
216 
217  // scale all the histograms based on the number of entries
218  xeffVw->Scale(1.0/nxcell);
219  yeffVw->Scale(1.0/nycell);
220  xeffVwcut->Scale(1.0/nxcellcut);
221  yeffVwcut->Scale(1.0/nycellcut);
222 
223 
224 
225  // set errors
226  for(int iw = 1; iw <= xeffVw->GetNbinsX(); ++iw){
227  double err = xeffVw->GetBinError(iw);
228  xeffVw->SetBinError(iw,(sqrt(err)/nxcell));
229  }
230  for(int iw = 1; iw <= yeffVw->GetNbinsX(); ++iw){
231  double err = yeffVw->GetBinError(iw);
232  yeffVw->SetBinError(iw,(sqrt(err)/nycell));
233  }
234  TCanvas* c4 = new TCanvas("c4","c4");
235  xeffVw->Draw();
236  yeffVw->Draw("same");
237  xeffVw->GetYaxis()->SetRangeUser(0.0,1.1);
238  xeffVw->SetLineColor(kBlue);
239  yeffVw->SetLineColor(kRed);
240  xeffVw->SetMarkerColor(kBlue);
241  yeffVw->SetMarkerColor(kRed);
242  TLegend* leg4 = new TLegend(0.55,0.15,0.85,0.35);
243  leg4->AddEntry(xeffVw,"Vertical Cells","LE");
244  leg4->AddEntry(yeffVw,"Horizontal Cells","LE");
245  leg4->Draw("same");
246 
247  TCanvas* c5 = new TCanvas("c5","c5");
248  TProfile* effVwprofx = effVwx->ProfileX("effVwprofx",1,-1,"s");
249  effVwprofx->GetYaxis()->SetRangeUser(0.0,1.1);
250  effVwprofx->GetYaxis()->SetTitle("Efficiency");
251  effVwprofx->Draw();
252  TProfile* effVwprofy = effVwy->ProfileX("effVwprofy",1,-1,"s");
253  effVwprofy->Draw("same");
254  effVwprofx->SetLineColor(kBlue);
255  effVwprofx->SetMarkerColor(kBlue);
256  effVwprofy->SetLineColor(kRed);
257  effVwprofy->SetMarkerColor(kRed);
258  TLegend* leg5 = new TLegend(0.55,0.15,0.85,0.35);
259  leg5->AddEntry(effVwprofx,"Vertical Cells","LE");
260  leg5->AddEntry(effVwprofy,"Horizontal Cells","LE");
261  leg5->Draw("same");
262 
263 
264  for(int iw = 1; iw <= xeffVwcut->GetNbinsX(); ++iw){
265  double err = xeffVwcut->GetBinError(iw);
266  xeffVwcut->SetBinError(iw,(sqrt(err)/nxcellcut));
267  }
268  for(int iw = 1; iw <= yeffVwcut->GetNbinsX(); ++iw){
269  double err = yeffVwcut->GetBinError(iw);
270  yeffVwcut->SetBinError(iw,(sqrt(err)/nycellcut));
271  }
272  TCanvas* c6 = new TCanvas("c6","c6");
273  xeffVwcut->Draw();
274  yeffVwcut->Draw("same");
275  xeffVwcut->GetYaxis()->SetRangeUser(0.0,1.1);
276  xeffVwcut->SetLineColor(kBlue);
277  yeffVwcut->SetLineColor(kRed);
278  xeffVwcut->SetMarkerColor(kBlue);
279  yeffVwcut->SetMarkerColor(kRed);
280  TLegend* leg6 = new TLegend(0.55,0.15,0.85,0.35);
281  leg6->AddEntry(xeffVwcut,"Vertical Cells","LE");
282  leg6->AddEntry(yeffVwcut,"Horizontal Cells","LE");
283  leg6->Draw("same");
284 
285  TCanvas* c7 = new TCanvas("c7","c7");
286  TProfile* effVwprofxcut = effVwxcut->ProfileX("effVwprofxcut",1,-1,"s");
287  effVwprofxcut->GetYaxis()->SetRangeUser(0.0,1.1);
288  effVwprofxcut->GetYaxis()->SetTitle("Efficiency");
289  effVwprofxcut->Draw();
290  TProfile* effVwprofycut = effVwycut->ProfileX("effVwprofycut",1,-1,"s");
291  effVwprofycut->Draw("same");
292  effVwprofxcut->SetLineColor(kBlue);
293  effVwprofxcut->SetMarkerColor(kBlue);
294  effVwprofycut->SetLineColor(kRed);
295  effVwprofycut->SetMarkerColor(kRed);
296  TLegend* leg7 = new TLegend(0.55,0.15,0.85,0.35);
297  leg7->AddEntry(effVwprofxcut,"Vertical Cells","LE");
298  leg7->AddEntry(effVwprofycut,"Horizontal Cells","LE");
299  leg7->Draw("same");
300 
301  TCanvas* c8 = new TCanvas("c8","c8");
302  effslowCellVPlanex->Draw("colz");
303  effslowCellVPlanex->GetZaxis()->SetRangeUser(0.0,1.0);
304 
305  TCanvas* c9 = new TCanvas("c9","c9");
306  effslowCellVPlaney->Draw("colz");
307  effslowCellVPlaney->GetZaxis()->SetRangeUser(0.0,1.0);
308 
309  TCanvas* c10 = new TCanvas("c10","c10");
310  effslowx->Draw();
311  effslowy->Draw("same");
312  effslowx->SetLineColor(kBlue);
313  effslowy->SetLineColor(kRed);
314  TLegend* leg10 = new TLegend(0.15,0.65,0.35,0.85);
315  leg10->AddEntry(effslowx,"X Cells","LE");
316  leg10->AddEntry(effslowy,"Y Cells","LE");
317  leg10->Draw("same");
318 
319  TCanvas* c11 = new TCanvas("c11","c11");
320  effsmidCellVPlanex->Draw("colz");
321  effsmidCellVPlanex->GetZaxis()->SetRangeUser(0.0,1.0);
322 
323  TCanvas* c12 = new TCanvas("c12","c12");
324  effsmidCellVPlaney->Draw("colz");
325  effsmidCellVPlaney->GetZaxis()->SetRangeUser(0.0,1.0);
326 
327  TCanvas* c13 = new TCanvas("c13","c13");
328  effsmidx->Draw();
329  effsmidy->Draw("same");
330  effsmidx->SetLineColor(kBlue);
331  effsmidy->SetLineColor(kRed);
332  TLegend* leg13 = new TLegend(0.15,0.65,0.35,0.85);
333  leg13->AddEntry(effsmidx,"X Cells","LE");
334  leg13->AddEntry(effsmidy,"Y Cells","LE");
335  leg13->Draw("same");
336 
337  TCanvas* c14 = new TCanvas("c14","c14");
338  effshighCellVPlanex->Draw("colz");
339  effshighCellVPlanex->GetZaxis()->SetRangeUser(0.0,1.0);
340 
341  TCanvas* c15 = new TCanvas("c15","c15");
342  effshighCellVPlaney->Draw("colz");
343  effshighCellVPlaney->GetZaxis()->SetRangeUser(0.0,1.0);
344 
345  TCanvas* c16 = new TCanvas("c16","c16");
346  effshighx->Draw();
347  effshighy->Draw("same");
348  effshighx->SetLineColor(kBlue);
349  effshighy->SetLineColor(kRed);
350  TLegend* leg16 = new TLegend(0.15,0.65,0.35,0.85);
351  leg16->AddEntry(effshighx,"X Cells","LE");
352  leg16->AddEntry(effshighy,"Y Cells","LE");
353  leg16->Draw("same");
354 
355 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
c2
Definition: demo5.py:33
void avgefficiency()
Definition: avgefficiency.C:20
c1
Definition: demo5.py:24
Float_t w
Definition: plot.C:20