plot_cosmic_eff.C
Go to the documentation of this file.
1 #include "TArrow.h"
2 #include "TCanvas.h"
3 #include "TFile.h"
4 #include "TGraph.h"
5 #include "TH1.h"
6 #include "TH2.h"
7 #include "TLegend.h"
8 #include "TROOT.h"
9 
10 #include <cassert>
11 #include <fstream>
12 #include <iostream>
13 #include <set>
14 
15 const int kThreshRun = 11572; // Threshold change
16 const int kMysteryRun = 12450;
17 
18 #ifndef __CINT__
19 #include <glob.h>
20 
21 std::vector<std::string> GetListOfFilesByWildcard(TString wildcardString)
22 {
23  std::vector<std::string> fileList;
24  glob_t g;
25  glob(wildcardString.Data(),
26  // GLOB_NOCHECK | // If no match return pattern
27  GLOB_TILDE, // Expand ~'s
28  0, &g);
29  for(unsigned int i = 0; i < g.gl_pathc; ++i)
30  fileList.push_back(g.gl_pathv[i]);
31  globfree(&g);
32 
33  // glob output is sorted by default
34  return fileList;
35 }
36 #endif
37 
38 void make_plots(TFile* f, TH2* hVsRun, TGraph* g, int run, std::string suffix, double xpos,
39  TH2*& cut)
40 {
41  TH2F* hits = (TH2F*)f->Get(("eff/hits"+suffix).c_str());
42  TH2F* miss = (TH2F*)f->Get(("eff/miss"+suffix).c_str());
43  TH2F* known = (TH2F*)f->Get("eff/known");
44 
45  if(!hits || !miss || !known) return;
46 
47  if(!cut){
48  cut = new TH2F(*hits);
49  cut->Reset();
50  }
51 
52  TH2F tot = *hits;
53  tot.Add(miss);
54 
55  TH2F frac = *miss;
56  frac.Divide(&tot);
57 
58  TH1F proj("", "", 100, 0, 1);
59 
60  int totCells = 0;
61 
62  const int X = frac.GetNbinsX()+1;
63  const int Y = frac.GetNbinsY()+1;
64  for(int x = 1; x < X; ++x){
65  for(int y = 1; y < Y; ++y){
66  const bool isBad = known->GetBinContent(x, y);
67  if(isBad) continue;
68 
69  // The occupancy cuts got it already
70  if(cut->GetBinContent(x, y)) continue;
71 
72  // Not enough stats to trust our efficiency fraction
73  if(tot.GetBinContent(x, y) < 20) continue;
74 
75  ++totCells;
76 
77  const double fr = frac.GetBinContent(x, y);
78 
79  if(suffix == "XYTight"){
80  if(run <= kThreshRun){
81  if(fr > .4) cut->SetBinContent(x, y, 1);
82  }
83  else{
84  if(fr > .15) cut->SetBinContent(x, y, 1);
85  }
86  }
87  if(suffix == "XYTight"){
88  if(run <= kThreshRun){
89  if(fr > .7) cut->SetBinContent(x, y, 1);
90  }
91  else{
92  if(fr > .5) cut->SetBinContent(x, y, 1);
93  }
94  }
95  if(suffix == ""){
96  if(run <= kMysteryRun){
97  if(fr > .8) cut->SetBinContent(x, y, 1);
98  }
99  else{
100  if(fr > .5) cut->SetBinContent(x, y, 1);
101  }
102  }
103 
104  if(tot.GetBinContent(x, y)){
105  proj.Fill(fr);
106  hVsRun->Fill(xpos, fr);
107  }
108  } // end for y
109  } // end for x
110 
111  // Don't include runs with few channels: noisy
112  if(g && totCells > 1000){
113  int iy;
114  // If it's worse than .95 we're cutting it regardless
115  int topbin = proj.GetXaxis()->FindBin(.95);
116  // Come down from there to find what would cut 1% of what's left
117  for(iy = topbin; iy >= 0; --iy){
118  if(proj.Integral(iy, topbin) > .01*totCells) break;
119  }
120  g->SetPoint(g->GetN(), xpos, proj.GetXaxis()->GetBinCenter(iy));
121  }
122 }
123 
125 {
126  assert(gROOT->IsBatch());
127 
128  // This is to implement a 1000 event cut. Should really be in the module itself.
129  std::set<std::pair<int, int> > goodRuns;
130  std::ifstream is("good_stats.txt");
131  while(is){
132  int run, subrun;
133  is >> run >> subrun;
134  goodRuns.insert(std::make_pair(run, subrun));
135  }
136 
137  // Assuming file location and naming scheme only true in my particular setup
138  std::vector<std::string> fs = GetListOfFilesByWildcard("~/pbs/cosmic_eff/out/eff_r*_t02_cosmic_S12.02.14.root");
139 
140  TH2* hVsRun = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
141  TH2* hVsRunXYTight = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
142  TH2* hVsRunZTight = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
143  TH2* hVsRunXYLoose = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
144  TH2* hVsRunZLoose = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
145 
146  TGraph* gXYTight = new TGraph;
147  TGraph* gXYLoose = new TGraph;
148  TGraph* gPoint = new TGraph;
149 
150  TGraph* gPointHand = new TGraph;
151  TGraph* gXYTightHand = new TGraph;
152  TGraph* gXYLooseHand = new TGraph;
153 
154  gPointHand->SetPoint(0, 0, .8);
155  gXYTightHand->SetPoint(0, 0, .4);
156  gXYLooseHand->SetPoint(0, 0, .7);
157 
158  for(unsigned int n = 0; n < fs.size(); ++n){
159  TFile f(fs[n].c_str());
160  if(f.IsZombie()) continue;
161 
162  const int p0 = fs[n].find("eff_r")+5;
163  const int p1 = fs[n].find("_s");
164  const int p2 = p1+2;
165  const int p3 = fs[n].find("_t02");
166  const int run = atoi(fs[n].substr(p0, p1-p0).c_str());
167  const int subrun = atoi(fs[n].substr(p2, p3-p2).c_str());
168 
169  if(!goodRuns.count(std::make_pair(run, subrun))){
170  std::cout << "Skipping for low stats..." << std::endl;
171  continue;
172  }
173 
174  if(run > kThreshRun && gXYTightHand->GetN() == 1){
175  gXYTightHand->SetPoint(1, n+.5, .4);
176  gXYTightHand->SetPoint(2, n+.5, .15);
177  gXYLooseHand->SetPoint(1, n+.5, .7);
178  gXYLooseHand->SetPoint(2, n+.5, .5);
179  }
180  if(run > kMysteryRun && gPointHand->GetN() == 1){
181  gPointHand->SetPoint(1, n+.5, .8);
182  gPointHand->SetPoint(2, n+.5, .5);
183  }
184 
185  std::cout << n << "/" << fs.size() << " : " << run << " " << subrun << std::endl;
186 
187  // Keep track of which ones have been cut already
188  TH2* cut = 0;
189 
190  make_plots(&f, hVsRunXYTight, gXYTight, run, "XYTight", n+.5, cut);
191  make_plots(&f, hVsRunXYLoose, gXYLoose, run, "XYLoose", n+.5, cut);
192  make_plots(&f, hVsRun, gPoint, run, "", n+.5, cut);
193 
194  make_plots(&f, hVsRunZTight, 0, run, "ZTight", n+.5, cut);
195  make_plots(&f, hVsRunZLoose, 0, run, "ZLoose", n+.5, cut);
196 
197  if(run%50 == 0 && subrun == 0){
198  hVsRun->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
199  hVsRunXYTight->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
200  hVsRunZTight->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
201  hVsRunXYLoose->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
202  hVsRunZLoose->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
203  }
204 
205  delete cut;
206  } // end for n
207 
208  gXYTightHand->SetPoint(3, fs.size(), .15);
209  gXYLooseHand->SetPoint(3, fs.size(), .5);
210  gPointHand->SetPoint(3, fs.size(), .5);
211 
212  hVsRun->GetXaxis()->SetTickLength(0);
213  hVsRun->Draw("colz");
214  gPad->SetLogz(false);
215  gPad->Print("plots/proj_vs_run.eps");
216  gPad->SetLogz();
217  gPad->Print("plots/proj_vs_run_log.eps");
218  gPoint->SetLineColor(kRed);
219  gPoint->Draw("l");
220  gPointHand->SetLineColor(kGreen);
221  gPointHand->Draw("l");
222  gPad->SetLogz(false);
223  gPad->Print("plots/proj_vs_run_trace.eps");
224  gPad->SetLogz();
225  gPad->Print("plots/proj_vs_run_trace_log.eps");
226 
227  hVsRunXYTight->GetXaxis()->SetTickLength(0);
228  hVsRunXYTight->Draw("colz");
229  gPad->SetLogz(false);
230  gPad->Print("plots/proj_vs_run_xytight.eps");
231  gPad->SetLogz();
232  gPad->Print("plots/proj_vs_run_xytight_log.eps");
233  gXYTight->SetLineColor(kRed);
234  gXYTight->Draw("l");
235  gXYTightHand->SetLineColor(kGreen);
236  gXYTightHand->Draw("l");
237  gPad->SetLogz(false);
238  gPad->Print("plots/proj_vs_run_xytight_trace.eps");
239  gPad->SetLogz();
240  gPad->Print("plots/proj_vs_run_xytight_trace_log.eps");
241 
242  hVsRunZTight->GetXaxis()->SetTickLength(0);
243  hVsRunZTight->Draw("colz");
244  gPad->SetLogz(false);
245  gPad->Print("plots/proj_vs_run_ztight.eps");
246  gPad->SetLogz();
247  gPad->Print("plots/proj_vs_run_ztight_log.eps");
248 
249  hVsRunXYLoose->GetXaxis()->SetTickLength(0);
250  hVsRunXYLoose->Draw("colz");
251  gPad->SetLogz(false);
252  gPad->Print("plots/proj_vs_run_xyloose.eps");
253  gPad->SetLogz();
254  gPad->Print("plots/proj_vs_run_xyloose_log.eps");
255  gXYLoose->SetLineColor(kRed);
256  gXYLoose->Draw("l");
257  gXYLooseHand->SetLineColor(kGreen);
258  gXYLooseHand->Draw("l");
259  gPad->SetLogz(false);
260  gPad->Print("plots/proj_vs_run_xyloose_trace.eps");
261  gPad->SetLogz();
262  gPad->Print("plots/proj_vs_run_xyloose_trace_log.eps");
263 
264  hVsRunZLoose->GetXaxis()->SetTickLength(0);
265  hVsRunZLoose->Draw("colz");
266  gPad->SetLogz(false);
267  gPad->Print("plots/proj_vs_run_zloose.eps");
268  gPad->SetLogz();
269  gPad->Print("plots/proj_vs_run_zloose_log.eps");
270 
271  TFile fout("cosmic_eff.root", "RECREATE");
272  hVsRun->Write("hVsRun");
273  hVsRunXYLoose->Write("hVsRunXYLoose");
274  hVsRunZLoose->Write("hVsRunZLoose");
275  hVsRunXYLoose->Write("hVsRunXYTight");
276  hVsRunZLoose->Write("hVsRunZTight");
277  fout.Close();
278 }
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
InputPoint * gPoint
void hits()
Definition: readHits.C:15
double frac(double x)
Fractional part.
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Float_t proj
Definition: plot.C:35
assert(nhit_max >=nhit_nbins)
const int kThreshRun
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const int kMysteryRun
Float_t X
Definition: plot.C:38
void plot_cosmic_eff()
void make_plots(TFile *f, TH2 *hVsRun, TGraph *g, int run, std::string suffix, double xpos, TH2 *&cut)