efficiency.C
Go to the documentation of this file.
1 // Macro for making plots of individual cell efficiencies vs w,
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. This will
5 // print a .png for every cell that an efficiency was calculated for in
6 // HitEfficiency_module.cc. This is best run in batch mode, unless you
7 // want a gizillion canvas's popping up in your face and then disappearing.
8
9 #include "TCanvas.h"
10 #include "TFile.h"
11 #include "TGraph.h"
12 #include "TH1.h"
13 #include "TH2.h"
14 #include "TMath.h"
15
16 #include <cmath>
17 #include <string>
18 #include <iostream>
19 #include <iomanip>
20 #include <cstdio>
21 #include <sstream>
22 #include <fstream>
23
24
25 string planestring(int plane)
26 {
27  std::ostringstream planestr;
28  if(plane<10){
29  planestr<<0;
30  planestr<<0;
31  planestr<<plane;
32  }
33  else if(plane<100){
34  planestr<<0;
35  planestr<<plane;
36  }
37  else{
38  planestr<<plane;
39  }
40  return planestr.str();
41 }
42
43 string cellstring(int cell)
44 {
45  std::ostringstream cellstr;
46  if(cell<10){
47  cellstr<<0;
48  cellstr<<0;
49  cellstr<<cell;
50  }
51  else if(cell<100){
52  cellstr<<0;
53  cellstr<<cell;
54  }
55  else{
56  cellstr<<cell;
57  }
58  return cellstr.str();
59 }
60
61
62 void efficiency()
63 {
64
66  // get the numerator and denominator for the efficiency
67  TH2F* hitsVwcell = tfile->Get("hiteff/hitsVwcelllencut");
68  TH2F* cellsVwcell = tfile->Get("hiteff/cellsVwcelllencut");
69
70  // Define some histogram binning
71  int nwbins = hitsVwcell->GetNbinsY();
72  double wmin = hitsVwcell->GetYaxis()->GetXmin();
73  double wmax = hitsVwcell->GetYaxis()->GetXmax();
74  int lowwbin = (int)((double)nwbins/3.0);
75  int highwbin = (int)(2.0*(double)nwbins/3.0);
76
77  // make a vector that holds all the cell detector states
78  const int npl = 32*28;
79  const int ncells = 32*12;
80
81  double warn = 0.80;
82  double good = 0.90;
83
84  unsigned short cellStatus[npl][ncells];
85  // default the status to nonexistent plot
86  for(unsigned int i = 0; i< npl; ++i){
87  for(unsigned int j = 0; j<ncells; ++j){
88  cellStatus[i][j] = 0;
89  }
90  }
91
92  // loop over every cell in the first diblock
93  for(int icell = 1; icell <= hitsVwcell->GetNbinsX(); ++icell){
94  // convert this to a plane and cell number
95  int cell = (icell-1)%384;
96  int plane = (icell-1-cell)/384;
97  int view = plane%2; // view = 0 y, view = 1 x
98
99  int celltot = cellsVwcell->ProjectionY("",icell,icell,"")->GetEntries();
100
101  // make the histogram and axis titles
102  string planename = planestring(plane);
103  string cellname = cellstring(cell);
104  TString axistitle("FD Cosmic Data - Plane "+planename+" Cell "+cellname+"; Distance From Center (cm);Efficiency");
105  TString histotitle("FD Cosmic Data - Plane "+planename+" Cell "+cellname);
106
107  // Make the efficiency histogram for this cell
108  TH1F* eff = new TH1F(histotitle.Data(),axistitle.Data(),nwbins,wmin,wmax);
109
110  // can't make a histogram for this cell, no data for the denominator
111  if(celltot == 0){ continue; }
112
113  // loop over every w bin
114  for(int iw = 1; iw <= nwbins; ++iw){
115  // get the efficiency at this cell/w, eff = k/n. k = number of hits at this cell/w, n = number of tricells at this cell/w
116  double k = hitsVwcell->GetBinContent(icell,iw);
117  double n = cellsVwcell->GetBinContent(icell,iw);
118
119  // calculate the variance of this efficiency
120  double var = (k+1)*(k+2)/((n+2)*(n+3)) - (k+1)*(k+1)/((n+2)*(n+2));
121
122  // Set the efficiency and error on this bin
123  if(n != 0){ eff->SetBinContent(iw,(k/n)); }
124  else{ eff->SetBinContent(iw,0); }
125  eff->SetBinError(iw,sqrt(var));
126  }
127
128  // Define some average efficiency metrics over different w ranges
129  double lowedgeeff = eff->Integral(1,lowwbin)/((double)(lowwbin - 1 + 1));
130  double mideff = eff->Integral(lowwbin+1,highwbin-1)/((double)(highwbin-1 - (lowwbin+1) + 1));
131  double highedgeeff = eff->Integral(highwbin,nwbins)/((double)(nwbins - (highwbin) + 1));
132
133  // calculate the cell state here. Assume bad until proven otherwise
134  cellStatus[plane][cell] = 1;
135  if(lowedgeeff > warn && mideff > warn && highedgeeff > warn){ cellStatus[plane][cell] = 2; }
136  if(lowedgeeff > good && mideff > good && highedgeeff > good){ cellStatus[plane][cell] = 3; }
137
138  // now make the efficiency histogram for this cell
139  // first get the plane and cell strings
140  TCanvas* c0 = new TCanvas("c0","c0");
141  eff->Draw();
142  eff->GetYaxis()->SetRangeUser(0.0,1.1);
143  TString viewstring("X");
144  if(!view){ viewstring = "Y"; }
145 // std::cout<<"Making printstring: planename: "<<planename<<" cellname: "<<cellname<<std::endl;
147  printstring.Append(viewstring);
148  printstring.Append("_"+planename+"_"+cellname+".png");
149  c0->Print(printstring.Data());
150 // std::cout<<"closing canvas"<<std::endl;
151  c0->Close();
152 // if(icell > 2){
153 // break;
154 // }
155
156  }
157
158
159  // loop over all the detector cells and write out the status
160  std::ofstream os("states.FD.js");
161  os << "states['FD'] = [\n";
162  for(int ipl = 0; ipl < npl; ++ipl){
163  os << "[";
164  for(int icl = 0; icl < ncells; ++icl){
165  int stat = cellStatus[ipl][icl];
166  os << stat;
167  if(icl != ncells-1){ os << "," ; }
168  }
169  os << "]";
170  if( ipl != npl-1){ os <<","; }
171  os << "\n";
172  }
173
174  os <<"];\n";
175
176
177 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double j
Definition: BetheBloch.cxx:29
string cellstring(int cell)
Definition: efficiency.C:43
def warn(msg)
print a warning message, from: https://cdcvs.fnal.gov/redmine/projects/novaart/repository/entry/trunk...
Definition: common_tools.py:16
string planestring(int plane)
Definition: efficiency.C:25
int ncells
Definition: geom.C:124
void efficiency()
Definition: efficiency.C:58