make_cosmic_eff_table.C

Go to the documentation of this file.

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

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

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

const std::string kRelease = "S12.02.14";

#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

struct OfflineChan
{
  OfflineChan(int p, int c) : plane(p), cell(c) {}

  int plane, cell;

  bool operator<(const OfflineChan& rhs) const
  {
    if(plane == rhs.plane) return cell < rhs.cell;
    return plane < rhs.plane;
  }
};

struct State
{
  State() {}
  State(int r, int s, bool g) : run(r), subrun(s), good(g) {}

  void Write(std::ostream& os, int plane, int cell, int endrun, int endsubrun) const
  {
    os << run << ", " << subrun << ", " << endrun << ", " << endsubrun << ", ndos, " << plane << ", " << cell << ", " << good << ", " << kRelease << ", dummy, dummy" << std::endl;
  }

  int run, subrun;
  bool good;
};

void make_cosmic_eff_table()
{
  // 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));
  }

  std::map<OfflineChan, State> states;

  std::ofstream table("cosmic_eff_table.csv");

  // 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");

  TH1* hCutFrac = new TH1F("", "", fs.size(), 0, fs.size());
  TH1* hCutFracCosmic = new TH1F("", "", fs.size(), 0, fs.size());
  TH1* hCutFracXYTight = new TH1F("", "", fs.size(), 0, fs.size());
  TH1* hCutFracXYLoose = new TH1F("", "", fs.size(), 0, fs.size());

  int prevsubrun = 0;
  int prevrun = 0;

  long int rowswritten = 0;
  long int totrows = 0;

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

    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;
    }

    // TODO: ensure monotonic

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

    TH2* frac = (TH2*)f.Get("eff/frac");
    TH2* htot = (TH2*)f.Get("eff/tot");
    TH2* fracXYTight = (TH2*)f.Get("eff/fracXYTight");
    TH2* fracXYLoose = (TH2*)f.Get("eff/fracXYLoose");
    TH2* totXYTight = (TH2*)f.Get("eff/totXYTight");
    TH2* totXYLoose = (TH2*)f.Get("eff/totXYLoose");

    TH2* known = (TH2*)f.Get("eff/known");

    if(!frac || !htot || !fracXYTight || !fracXYLoose || !totXYTight || !totXYLoose) continue;

    double totkept = 0;
    double totcut = 0;
    double totcutp = 0;
    double totcutxyt = 0;
    double totcutxyl = 0;

    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 double t = htot->GetBinContent(x, y);
        const double txyt = totXYTight->GetBinContent(x, y);
        const double txyl = totXYLoose->GetBinContent(x, y);

        // Not enough stats. Just use whatever the old value was?
        if(txyt < 20 && txyl < 20 && t < 20) continue;

        ++totcells;

        // Underflow bins
        const int plane = x-1;
        const int cell = y-1;

        bool good = true;

        if(txyt >= 20){
          const double frXY = fracXYTight->GetBinContent(x, y);
          if((run <= kThreshRun && frXY > .4) || (run > kThreshRun && frXY > .15)){
            good = false;
            ++totcutxyt;
          }
        }

        if(good && txyl >= 20){
          const double frXY = fracXYLoose->GetBinContent(x, y);
          if((run <= kThreshRun && frXY > .7) || (run > kThreshRun && frXY > .5)){
            good = false;
            ++totcutxyl;
          }
        }

        if(good && t >= 20){
          const double fr = frac->GetBinContent(x, y);
          if((run <= kMysteryRun && fr > .8) || (run > kMysteryRun && fr > .5)){
            good = false;
            ++totcutp;
          }
        }

        if(known && !known->GetBinContent(x, y)){
          if(good)
            ++totkept;
          else
            ++totcut;
        }

        ++totrows;
        OfflineChan chan(plane, cell);
        std::map<OfflineChan, State>::iterator it = states.find(chan);
        if(it == states.end()){
          // Begin
          states[chan] = State(run, subrun, good);
        }
        else{
          const State oldstate = it->second;
          if(oldstate.good != good){
            // Write out the old one
            oldstate.Write(table, plane, cell, prevrun, prevsubrun);
            ++rowswritten;
            // And start a new one
            states[chan] = State(run, subrun, good);
          }
        }
      } // end for y
    } // end for x

    std::cout << (100*totcut)/(totkept+totcut) << "% cut" << std::endl;
    const double tot = totkept+totcut;
    if(tot && totcells > 1000){
      hCutFrac->SetBinContent(n, 100*totcut/tot);
      hCutFracCosmic->SetBinContent(n, 100*totcutp/tot);
      hCutFracXYTight->SetBinContent(n, 100*totcutxyt/tot);
      hCutFracXYLoose->SetBinContent(n, 100*totcutxyl/tot);
    }

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

    prevrun = run;
    prevsubrun = subrun;
  } // end for n

  // Finish off everything else
  for(std::map<OfflineChan, State>::iterator it = states.begin(); it != states.end(); ++it){
    it->second.Write(table, it->first.plane, it->first.cell, prevrun, prevsubrun);
    ++rowswritten;
  }
  std::cout << "Wrote " << rowswritten << "/" << totrows << " = " << (10000*rowswritten/totrows)/100. << "%" << std::endl;

  hCutFrac->Draw();
  hCutFrac->GetYaxis()->SetTitle("Percentage of channels cut");
  hCutFrac->GetXaxis()->SetTickLength(0);
  hCutFracCosmic->SetLineColor(kBlue);
  hCutFracCosmic->Draw("same");
  hCutFracXYTight->SetLineColor(kRed);
  hCutFracXYTight->Draw("same");
  hCutFracXYLoose->SetLineColor(kGreen+2);
  hCutFracXYLoose->Draw("same");

  gPad->Print("plots/frac_cut.eps");
  hCutFrac->GetYaxis()->SetRangeUser(0, 10);
  gPad->Print("plots/frac_cut_zoom.eps");
}