plot_cosmic_eff.C

Go to the documentation of this file.

#include "TArrow.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TLegend.h"
#include "TROOT.h"

#include <cassert>
#include <fstream>
#include <iostream>
#include <set>

const int kThreshRun = 11572; // Threshold change
const int kMysteryRun = 12450;

#ifndef __CINT__
#include <glob.h>

std::vector<std::string> GetListOfFilesByWildcard(TString wildcardString)
{
  std::vector<std::string> fileList;
  glob_t g;
  glob(wildcardString.Data(),
       // GLOB_NOCHECK | // If no match return pattern
       GLOB_TILDE, // Expand ~'s
       0, &g);
  for(unsigned int i = 0; i < g.gl_pathc; ++i)
    fileList.push_back(g.gl_pathv[i]);
  globfree(&g);

  // glob output is sorted by default
  return fileList;
}
#endif

void make_plots(TFile* f, TH2* hVsRun, TGraph* g, int run, std::string suffix, double xpos,
                TH2*& cut)
{
  TH2F* hits = (TH2F*)f->Get(("eff/hits"+suffix).c_str());
  TH2F* miss = (TH2F*)f->Get(("eff/miss"+suffix).c_str());
  TH2F* known = (TH2F*)f->Get("eff/known");

  if(!hits || !miss || !known) return;

  if(!cut){
    cut = new TH2F(*hits);
    cut->Reset();
  }

  TH2F tot = *hits;
  tot.Add(miss);

  TH2F frac = *miss;
  frac.Divide(&tot);

  TH1F proj("", "", 100, 0, 1);

  int totCells = 0;

  const int X = frac.GetNbinsX()+1;
  const int Y = frac.GetNbinsY()+1;
  for(int x = 1; x < X; ++x){
    for(int y = 1; y < Y; ++y){
      const bool isBad = known->GetBinContent(x, y);
      if(isBad) continue;

      // The occupancy cuts got it already
      if(cut->GetBinContent(x, y)) continue;

      // Not enough stats to trust our efficiency fraction
      if(tot.GetBinContent(x, y) < 20) continue;

      ++totCells;

      const double fr = frac.GetBinContent(x, y);

      if(suffix == "XYTight"){
        if(run <= kThreshRun){
          if(fr > .4) cut->SetBinContent(x, y, 1);
        }
        else{
          if(fr > .15) cut->SetBinContent(x, y, 1);
        }
      }
      if(suffix == "XYTight"){
        if(run <= kThreshRun){
          if(fr > .7) cut->SetBinContent(x, y, 1);
        }
        else{
          if(fr > .5) cut->SetBinContent(x, y, 1);
        }
      }
      if(suffix == ""){
        if(run <= kMysteryRun){
          if(fr > .8) cut->SetBinContent(x, y, 1);
        }
        else{
          if(fr > .5) cut->SetBinContent(x, y, 1);
        }
      }

      if(tot.GetBinContent(x, y)){
        proj.Fill(fr);
        hVsRun->Fill(xpos, fr);
      }
    } // end for y
  } // end for x

  // Don't include runs with few channels: noisy
  if(g && totCells > 1000){
    int iy;
    // If it's worse than .95 we're cutting it regardless
    int topbin = proj.GetXaxis()->FindBin(.95);
    // Come down from there to find what would cut 1% of what's left
    for(iy = topbin; iy >= 0; --iy){
      if(proj.Integral(iy, topbin) > .01*totCells) break;
    }
    g->SetPoint(g->GetN(), xpos, proj.GetXaxis()->GetBinCenter(iy));
  }
}

void plot_cosmic_eff()
{
  assert(gROOT->IsBatch());

  // This is to implement a 1000 event cut. Should really be in the module itself.
  std::set<std::pair<int, int> > goodRuns;
  std::ifstream is("good_stats.txt");
  while(is){
    int run, subrun;
    is >> run >> subrun;
    goodRuns.insert(std::make_pair(run, subrun));
  }

  // Assuming file location and naming scheme only true in my particular setup
  std::vector<std::string> fs = GetListOfFilesByWildcard("~/pbs/cosmic_eff/out/eff_r*_t02_cosmic_S12.02.14.root");

  TH2* hVsRun = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
  TH2* hVsRunXYTight = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
  TH2* hVsRunZTight = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
  TH2* hVsRunXYLoose = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);
  TH2* hVsRunZLoose = new TH2F("", ";Run;Fraction of times missed from track", fs.size(), 0, fs.size(), 101, 0, 1.01);

  TGraph* gXYTight = new TGraph;
  TGraph* gXYLoose = new TGraph;
  TGraph* gPoint = new TGraph;

  TGraph* gPointHand = new TGraph;
  TGraph* gXYTightHand = new TGraph;
  TGraph* gXYLooseHand = new TGraph;

  gPointHand->SetPoint(0, 0, .8);
  gXYTightHand->SetPoint(0, 0, .4);
  gXYLooseHand->SetPoint(0, 0, .7);

  for(unsigned int n = 0; n < fs.size(); ++n){
    TFile f(fs[n].c_str());
    if(f.IsZombie()) continue;

    const int p0 = fs[n].find("eff_r")+5;
    const int p1 = fs[n].find("_s");
    const int p2 = p1+2;
    const int p3 = fs[n].find("_t02");
    const int run = atoi(fs[n].substr(p0, p1-p0).c_str());
    const int subrun = atoi(fs[n].substr(p2, p3-p2).c_str());

    if(!goodRuns.count(std::make_pair(run, subrun))){
      std::cout << "Skipping for low stats..." << std::endl;
      continue;
    }

    if(run > kThreshRun && gXYTightHand->GetN() == 1){
      gXYTightHand->SetPoint(1, n+.5, .4);
      gXYTightHand->SetPoint(2, n+.5, .15);
      gXYLooseHand->SetPoint(1, n+.5, .7);
      gXYLooseHand->SetPoint(2, n+.5, .5);
    }
    if(run > kMysteryRun && gPointHand->GetN() == 1){
      gPointHand->SetPoint(1, n+.5, .8);
      gPointHand->SetPoint(2, n+.5, .5);
    }

    std::cout << n << "/" << fs.size() << " : " << run << " " << subrun << std::endl;

    // Keep track of which ones have been cut already
    TH2* cut = 0;

    make_plots(&f, hVsRunXYTight, gXYTight, run, "XYTight", n+.5, cut);
    make_plots(&f, hVsRunXYLoose, gXYLoose, run, "XYLoose", n+.5, cut);
    make_plots(&f, hVsRun, gPoint, run, "", n+.5, cut);

    make_plots(&f, hVsRunZTight, 0, run, "ZTight", n+.5, cut);
    make_plots(&f, hVsRunZLoose, 0, run, "ZLoose", n+.5, cut);

    if(run%50 == 0 && subrun == 0){
      hVsRun->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
      hVsRunXYTight->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
      hVsRunZTight->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
      hVsRunXYLoose->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
      hVsRunZLoose->GetXaxis()->SetBinLabel(n+1, TString::Format("%d", run));
    }

    delete cut;
  } // end for n

  gXYTightHand->SetPoint(3, fs.size(), .15);
  gXYLooseHand->SetPoint(3, fs.size(), .5);
  gPointHand->SetPoint(3, fs.size(), .5);

  hVsRun->GetXaxis()->SetTickLength(0);
  hVsRun->Draw("colz");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_log.eps");
  gPoint->SetLineColor(kRed);
  gPoint->Draw("l");
  gPointHand->SetLineColor(kGreen);
  gPointHand->Draw("l");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run_trace.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_trace_log.eps");

  hVsRunXYTight->GetXaxis()->SetTickLength(0);
  hVsRunXYTight->Draw("colz");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run_xytight.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_xytight_log.eps");
  gXYTight->SetLineColor(kRed);
  gXYTight->Draw("l");
  gXYTightHand->SetLineColor(kGreen);
  gXYTightHand->Draw("l");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run_xytight_trace.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_xytight_trace_log.eps");

  hVsRunZTight->GetXaxis()->SetTickLength(0);
  hVsRunZTight->Draw("colz");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run_ztight.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_ztight_log.eps");

  hVsRunXYLoose->GetXaxis()->SetTickLength(0);
  hVsRunXYLoose->Draw("colz");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run_xyloose.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_xyloose_log.eps");
  gXYLoose->SetLineColor(kRed);
  gXYLoose->Draw("l");
  gXYLooseHand->SetLineColor(kGreen);
  gXYLooseHand->Draw("l");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run_xyloose_trace.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_xyloose_trace_log.eps");

  hVsRunZLoose->GetXaxis()->SetTickLength(0);
  hVsRunZLoose->Draw("colz");
  gPad->SetLogz(false);
  gPad->Print("plots/proj_vs_run_zloose.eps");
  gPad->SetLogz();
  gPad->Print("plots/proj_vs_run_zloose_log.eps");

  TFile fout("cosmic_eff.root", "RECREATE");
  hVsRun->Write("hVsRun");
  hVsRunXYLoose->Write("hVsRunXYLoose");
  hVsRunZLoose->Write("hVsRunZLoose");
  hVsRunXYLoose->Write("hVsRunXYTight");
  hVsRunZLoose->Write("hVsRunZTight");
  fout.Close();
}