make_cosmic_eff_table.C
Go to the documentation of this file.
1 #include "TArrow.h"
2 #include "TCanvas.h"
3 #include "TFile.h"
4 #include "TH1.h"
5 #include "TH2.h"
6 #include "TLegend.h"
7 #include "TROOT.h"
8 
9 #include <cassert>
10 #include <iostream>
11 #include <fstream>
12 #include <set>
13 
14 const int kThreshRun = 11572; // Threshold change
15 const int kMysteryRun = 12450;
16 
17 const std::string kRelease = "S12.02.14";
18 
19 #ifndef __CINT__
20 #include <glob.h>
21 
22 std::vector<std::string> GetListOfFilesByWildcard(TString wildcardString)
23 {
24  std::vector<std::string> fileList;
25  glob_t g;
26  glob(wildcardString.Data(),
27  // GLOB_NOCHECK | // If no match return pattern
28  GLOB_TILDE, // Expand ~'s
29  0, &g);
30  for(unsigned int i = 0; i < g.gl_pathc; ++i)
31  fileList.push_back(g.gl_pathv[i]);
32  globfree(&g);
33 
34  // glob output is sorted by default
35  return fileList;
36 }
37 #endif
38 
40 {
41  OfflineChan(int p, int c) : plane(p), cell(c) {}
42 
43  int plane, cell;
44 
45  bool operator<(const OfflineChan& rhs) const
46  {
47  if(plane == rhs.plane) return cell < rhs.cell;
48  return plane < rhs.plane;
49  }
50 };
51 
52 struct State
53 {
54  State() {}
55  State(int r, int s, bool g) : run(r), subrun(s), good(g) {}
56 
57  void Write(std::ostream& os, int plane, int cell, int endrun, int endsubrun) const
58  {
59  os << run << ", " << subrun << ", " << endrun << ", " << endsubrun << ", ndos, " << plane << ", " << cell << ", " << good << ", " << kRelease << ", dummy, dummy" << std::endl;
60  }
61 
62  int run, subrun;
63  bool good;
64 };
65 
67 {
68  // This is to implement a 1000 event cut. Should really be in the module itself.
69  std::set<std::pair<int, int> > goodRuns;
70  std::ifstream is("good_stats.txt");
71  while(is){
72  int run, subrun;
73  is >> run >> subrun;
74  goodRuns.insert(std::make_pair(run, subrun));
75  }
76 
77  std::map<OfflineChan, State> states;
78 
79  std::ofstream table("cosmic_eff_table.csv");
80 
81  // Assuming file location and naming scheme only true in my particular setup
82  std::vector<std::string> fs = GetListOfFilesByWildcard("~/pbs/cosmic_eff/out/eff_r*_t02_cosmic_S12.02.14.root");
83 
84  TH1* hCutFrac = new TH1F("", "", fs.size(), 0, fs.size());
85  TH1* hCutFracCosmic = new TH1F("", "", fs.size(), 0, fs.size());
86  TH1* hCutFracXYTight = new TH1F("", "", fs.size(), 0, fs.size());
87  TH1* hCutFracXYLoose = new TH1F("", "", fs.size(), 0, fs.size());
88 
89  int prevsubrun = 0;
90  int prevrun = 0;
91 
92  long int rowswritten = 0;
93  long int totrows = 0;
94 
95  for(unsigned int n = 0; n < fs.size(); ++n){
96  TFile f(fs[n].c_str());
97  if(f.IsZombie()) return;
98 
99  const int p0 = fs[n].find("eff_r")+5;
100  const int p1 = fs[n].find("_s");
101  const int p2 = p1+2;
102  const int p3 = fs[n].find("_t02");
103  const int run = atoi(fs[n].substr(p0, p1-p0).c_str());
104  const int subrun = atoi(fs[n].substr(p2, p3-p2).c_str());
105 
106  if(!goodRuns.count(std::make_pair(run, subrun))){
107  std::cout << "Skipping for low stats..." << std::endl;
108  continue;
109  }
110 
111  // TODO: ensure monotonic
112 
113  std::cout << n << "/" << fs.size() << " : " << run << " " << subrun << std::endl;
114 
115  TH2* frac = (TH2*)f.Get("eff/frac");
116  TH2* htot = (TH2*)f.Get("eff/tot");
117  TH2* fracXYTight = (TH2*)f.Get("eff/fracXYTight");
118  TH2* fracXYLoose = (TH2*)f.Get("eff/fracXYLoose");
119  TH2* totXYTight = (TH2*)f.Get("eff/totXYTight");
120  TH2* totXYLoose = (TH2*)f.Get("eff/totXYLoose");
121 
122  TH2* known = (TH2*)f.Get("eff/known");
123 
124  if(!frac || !htot || !fracXYTight || !fracXYLoose || !totXYTight || !totXYLoose) continue;
125 
126  double totkept = 0;
127  double totcut = 0;
128  double totcutp = 0;
129  double totcutxyt = 0;
130  double totcutxyl = 0;
131 
132  int totcells = 0;
133 
134  const int X = frac->GetNbinsX()+1;
135  const int Y = frac->GetNbinsY()+1;
136  for(int x = 1; x < X; ++x){
137  for(int y = 1; y < Y; ++y){
138  const double t = htot->GetBinContent(x, y);
139  const double txyt = totXYTight->GetBinContent(x, y);
140  const double txyl = totXYLoose->GetBinContent(x, y);
141 
142  // Not enough stats. Just use whatever the old value was?
143  if(txyt < 20 && txyl < 20 && t < 20) continue;
144 
145  ++totcells;
146 
147  // Underflow bins
148  const int plane = x-1;
149  const int cell = y-1;
150 
151  bool good = true;
152 
153  if(txyt >= 20){
154  const double frXY = fracXYTight->GetBinContent(x, y);
155  if((run <= kThreshRun && frXY > .4) || (run > kThreshRun && frXY > .15)){
156  good = false;
157  ++totcutxyt;
158  }
159  }
160 
161  if(good && txyl >= 20){
162  const double frXY = fracXYLoose->GetBinContent(x, y);
163  if((run <= kThreshRun && frXY > .7) || (run > kThreshRun && frXY > .5)){
164  good = false;
165  ++totcutxyl;
166  }
167  }
168 
169  if(good && t >= 20){
170  const double fr = frac->GetBinContent(x, y);
171  if((run <= kMysteryRun && fr > .8) || (run > kMysteryRun && fr > .5)){
172  good = false;
173  ++totcutp;
174  }
175  }
176 
177  if(known && !known->GetBinContent(x, y)){
178  if(good)
179  ++totkept;
180  else
181  ++totcut;
182  }
183 
184  ++totrows;
185  OfflineChan chan(plane, cell);
186  std::map<OfflineChan, State>::iterator it = states.find(chan);
187  if(it == states.end()){
188  // Begin
189  states[chan] = State(run, subrun, good);
190  }
191  else{
192  const State oldstate = it->second;
193  if(oldstate.good != good){
194  // Write out the old one
195  oldstate.Write(table, plane, cell, prevrun, prevsubrun);
196  ++rowswritten;
197  // And start a new one
198  states[chan] = State(run, subrun, good);
199  }
200  }
201  } // end for y
202  } // end for x
203 
204  std::cout << (100*totcut)/(totkept+totcut) << "% cut" << std::endl;
205  const double tot = totkept+totcut;
206  if(tot && totcells > 1000){
207  hCutFrac->SetBinContent(n, 100*totcut/tot);
208  hCutFracCosmic->SetBinContent(n, 100*totcutp/tot);
209  hCutFracXYTight->SetBinContent(n, 100*totcutxyt/tot);
210  hCutFracXYLoose->SetBinContent(n, 100*totcutxyl/tot);
211  }
212 
213  if(run%50 == 0 && subrun == 0){
214  hCutFrac->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
215  }
216 
217  prevrun = run;
218  prevsubrun = subrun;
219  } // end for n
220 
221  // Finish off everything else
222  for(std::map<OfflineChan, State>::iterator it = states.begin(); it != states.end(); ++it){
223  it->second.Write(table, it->first.plane, it->first.cell, prevrun, prevsubrun);
224  ++rowswritten;
225  }
226  std::cout << "Wrote " << rowswritten << "/" << totrows << " = " << (10000*rowswritten/totrows)/100. << "%" << std::endl;
227 
228  hCutFrac->Draw();
229  hCutFrac->GetYaxis()->SetTitle("Percentage of channels cut");
230  hCutFrac->GetXaxis()->SetTickLength(0);
231  hCutFracCosmic->SetLineColor(kBlue);
232  hCutFracCosmic->Draw("same");
233  hCutFracXYTight->SetLineColor(kRed);
234  hCutFracXYTight->Draw("same");
235  hCutFracXYLoose->SetLineColor(kGreen+2);
236  hCutFracXYLoose->Draw("same");
237 
238  gPad->Print("plots/frac_cut.eps");
239  hCutFrac->GetYaxis()->SetRangeUser(0, 10);
240  gPad->Print("plots/frac_cut_zoom.eps");
241 }
bool operator<(const OfflineChan &rhs) const
const int kThreshRun
enum BeamMode kRed
set< int >::iterator it
const char * p
Definition: xmltok.h:285
void make_cosmic_eff_table()
std::vector< std::string > GetListOfFilesByWildcard(TString wildcardString)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
list fileList
Definition: cafExposure.py:30
Float_t Y
Definition: plot.C:38
const XML_Char * s
Definition: expat.h:262
void Write(std::ostream &os, int plane, int cell, int endrun, int endsubrun) const
const std::string kRelease
const int kMysteryRun
double frac(double x)
Fractional part.
OfflineChan(int p, int c)
State(int r, int s, bool g)
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
TRandom3 r(0)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Float_t X
Definition: plot.C:38
enum BeamMode kGreen
enum BeamMode kBlue
enum BeamMode string