pass1_plots.C

Go to the documentation of this file.

const char* gPrefix = 0; // Plot files will carry this prefix
int         gMode   = 1; // We are comparing Trident MC to GENIE MC
// int      gMode   = 2; // ................ ND Data to GENIE MC (?)
int         gPrintPlots = 1; // Should we save plots to .pdf?
int         gPlotIt     = 1; // Turns plots on/off during development

//......................................................................

TFile* fpSIG = 0;
TFile* fpBKG = 0;

TTree* ntSIG = 0;
TTree* ntBKG = 0;

void open_files() 
{
  if (gMode==1) {
    std::cout << "Comparing trident MC to GENIE MC" << std::endl;
    gPrefix = "p1_T2G";
    fpSIG = TFile::Open("dimuon_signal_hist.root");
    fpBKG = TFile::Open("dimuon_background_hist.root");
    
    ntSIG = (TTree*)fpSIG->FindObjectAny("p1");
    ntBKG = (TTree*)fpBKG->FindObjectAny("p1");
  }
  else {
    std::cout << "Ugh oh - this mode isn't supported..." << std::endl;
  }
}

//......................................................................
//
// A lot of the plotting is repetitive. Make some functions to do the
// routine stuff and to enforce a good style.
//
TCanvas* canvas(const char* nm, const char* ti, const char* a)
{
  if (a == TString("square"))    return new TCanvas(nm, nm, 600,600);
  if (a == TString("portrait"))  return new TCanvas(nm, nm, 500,700);
  if (a == TString("landscape")) return new TCanvas(nm, nm, 700,500);
  return 0;
}
void area_normalize(TH1F* h1, TH1F* h2) 
{
  h2->Scale(h1->GetEntries()/h2->GetEntries());
}
float plot_max(TH1F* h1, TH1F* h2, float r, int islog=0)
{
  float m1 = h1->GetMaximum();
  float m2 = h2->GetMaximum();
  if (islog) {
    m1 = log10(m1);
    m2 = log10(m2);
  }

  float mx = 0;
  if (m1>m2) mx = (std::floor(m1/r)+1)*r;
  else       mx = (std::floor(m2/r)+1)*r;
  if (islog) mx = pow(10.0, mx);

  h1->SetMaximum(mx);
  h2->SetMaximum(mx);

  if (islog) {
    m1 = h1->GetMinimum();
    m2 = h2->GetMinimum();
    if (m1<=1 || m2<=1) h1->SetMinimum(0.1);
    else                h1->SetMinimum(1.0);
  }
  else {
    h1->SetMinimum(0);
  }

  return mx;
}
void show_x_cut(float x, float dx, float y1, float y2, float fy=0.82) 
{
  TLine* l = new TLine(x, y1, x, y2);
  l->SetLineColor(kBlack);
  l->SetLineStyle(3);
  l->Draw();

  double y = y1 + fy*(y2-y1);
  TArrow* a = new TArrow(x, y, x+dx, y, 0.02);
  
  a->SetLineColor(kBlack);
  a->Draw();
}
void sig_style(TH1F* h) 
{
  h->SetFillColor(kAzure-9);
  h->GetXaxis()->SetNdivisions(505);
  h->GetYaxis()->SetNdivisions(505);
}
void bkg_style(TH1F* h) 
{
  h->SetLineColor(kRed);
  h->SetLineWidth(1);
}

//......................................................................
//
// Again, the projecting and plotting of variables is repetitive. Code
// up the common features.
//
void proj_1D(TH1F**      h1,
             TH1F**      h2, 
             const char* base_name,
             const char* title,
             int         n,
             float       x1,
             float       x2,
             const char* var,
             TCut        cut)
{
  char sig_name[256];
  char bkg_name[256];
  sprintf(sig_name, "sig%s", base_name);
  sprintf(bkg_name, "bkg%s", base_name);

  TH1F* sig = new TH1F(sig_name, title, n, x1, x2);
  TH1F* bkg = new TH1F(bkg_name, title, n, x1, x2);
  (*h1) = sig;
  (*h2) = bkg;
  
  ntSIG->Project(sig_name, var, cut);
  ntBKG->Project(bkg_name, var, cut);

  std::cout << "|--* " 
            << base_name << " & "
            << sig->GetEntries() << " & "
            << bkg->GetEntries() << " & "
            << cut << " \\\\"
            << std::endl;

  area_normalize(sig, bkg);
}

//......................................................................

void plot_1D(TH1F*       h1,
             TH1F*       h2,
             const char* title,
             const char* aspect,
             int         islogy,
             float       r,
             float       cut_left=-99,
             float       cut_right=-99)
{
  char ti[256];
  sprintf(ti, "%s_%s", gPrefix, title);
  TCanvas* c = canvas(ti, 0, aspect);
  if (islogy==1) c->SetLogy();

  float mx = plot_max(h1, h2, r, islogy);
  
  sig_style(h1);
  bkg_style(h2);
  
  float xlo = h1->GetXaxis()->GetXmin();
  float xhi = h1->GetXaxis()->GetXmax();
  if (fabs((xlo+xhi)/(xhi-xlo))<0.05) {
    h1->GetXaxis()->CenterTitle();
  }
  h1->Draw();
  h2->Draw("hist,same");

  float dx = 0.05*(xhi-xlo);
  if (cut_left !=-99) { show_x_cut(cut_left,   dx, 0, mx); }
  if (cut_right!=-99) { show_x_cut(cut_right, -dx, 0, mx); }

  if (gPrintPlots) {
    char fname[256];
    sprintf(fname,"%s.pdf",c->GetName());
    c->SaveAs(fname);
  }
}

//......................................................................
TCut gCutNHIT;
TCut gCutVETO;
TCut gCutVTX;
TCut gCutSHAPE;
TCut gCutPASS1;

//......................................................................
//
// Here's where I walk through the distributions and make some cut
// decisions.
//
// The cuts proceed in roughly two big passes. First, are some basic
// quality cuts which ensure that we have events with enough
// information to have a good chance at a correct reconstruction and
// containment. We really don't have much choice about these. The
// second pass is to try to do some signal/background separation based
// on the event topology. Here I've chosen to try to be ~99%
// efficient (as a fraction of the selections above) while rejecting
// as much background as possible.
//
// Note: If you turn the plotting off, this will set all the global
// cuts.
//
void make_plots() 
{
  int liny   = 0; // Put the plot on linear y scale
  int logy   = 1; // Put the plot on log y scale

  /*====================================================================
    First level of cuts is to make sure we have enough information to
    work with in the event.

    We are searching for two muon tracks. Backgrounds are likely to be
    long muons (rock and cosmics) and events with long pion tracks.

    The interaction length in the detector is lam_I = ~40 cm so pions
    are very unlikely to travel more than 3 lam_I ~= 120 cm or about
    20 planes in the detecor. Rock muons and cosmics may traverse the
    entire detector length, ~220 planes.

    If I have two muon tracks traveling across 20 planes I should
    expect ~40 hits in the event. If the events are oriented along the
    beam direction, there should be roughly equal number of hits in
    each view.
  */
  TCut cutNONE = TCut("1");
  if (gPlotIt) {
  TH1F* sigNpl = 0;
  TH1F* bkgNpl = 0;
  proj_1D(&sigNpl, &bkgNpl,    "Nplane", ";number of hit planes;slice counts/5 planes", 44, 0, 220, "np", cutNONE);
  plot_1D( sigNpl,  bkgNpl, "01_Nplane", "portrait", logy, 0.5, 20, 200);
  }
  TCut cutNp("np>=20&&np<=200");

  if (gPlotIt) {
  TH1F* sigNhit = 0;
  TH1F* bkgNhit = 0;
  proj_1D(&sigNhit, &bkgNhit,    "Nhit", ";number of hit cells;slice counts/5 hits", 80, 0, 400, "(nhitx+nhity)", cutNONE);
  plot_1D( sigNhit,  bkgNhit, "01_Nhit", "portrait", logy, 0.5, 30, -99);
  }
  TCut cutNhit("(nhitx+nhity)>=30");

  if (gPlotIt) {
  TH1F* sigAnhit = 0;
  TH1F* bkgAnhit = 0;
  proj_1D(&sigAnhit, &bkgAnhit,    "Anhit", ";(N_{x}-N_{y})/(N_{x}+N_{y});slice counts/0.05", 20, -1, 1, "(nhitx-nhity)/(nhitx+nhity)", cutNp&&cutNhit);
  plot_1D( sigAnhit,  bkgAnhit, "01_Anhit", "portrait", logy, 0.5, -0.5, 0.5);
  }
  TCut cutAnhit("fabs(nhitx-nhity)/(nhitx+nhity)<=0.5");

  if (gPlotIt) {
  TH1F* sigN2p = 0;
  TH1F* bkgN2p = 0;
  proj_1D(&sigN2p, &bkgN2p,    "N2p", ";number of planes with 2 hits;slice counts", 100, 0, 100, "n2p", cutNp&&cutNhit&&cutAnhit);
  plot_1D( sigN2p,  bkgN2p, "01_N2p", "portrait", liny, 20, 8, -99);
  }
  TCut cutN2p("n2p>=8");

  if (gPlotIt) {
  TH1F* sigN3p = 0;
  TH1F* bkgN3p = 0;
  proj_1D(&sigN3p, &bkgN3p,    "N3p", ";number of planes with 3+ hits;slice counts", 100, 0, 100, "n3p", cutNp&&cutNhit&&cutAnhit);
  plot_1D( sigN3p,  bkgN3p, "01_N3p", "portrait", logy, 0.5, -99, 40);
  }
  TCut cutN3p("n3p<=40");

  //
  // Put this level of selection into a single cut to be used later.
  //
  gCutNHIT = TCut(cutNp&&cutNhit&&cutAnhit&&cutN2p&&cutN3p);
  //====================================================================

  /*====================================================================
    Next look at containment. This done by drawing a buffer zone of 5
    planes on the front face, and 5 cells on the top, bottom, left,
    and right and counting how many hits occur in those zones. Then we
    look where we guessed the event vertex to be and require it to be
    inside a fiducial volume.
  */
  if (gPlotIt) {
  TH1F* sigNp5 = 0;
  TH1F* bkgNp5 = 0;
  proj_1D(&sigNp5, &bkgNp5,    "Np5", ";number of hits on front face;slice counts/hit", 51, -0.5, 50.5, "np5", gCutNHIT);
  plot_1D( sigNp5,  bkgNp5, "02_Np5", "portrait", logy, 0.5, -99, 1.5);
  }
  TCut cutNp5("np5<=1");

  if (gPlotIt) {
  TH1F* sigNwall = 0;
  TH1F* bkgNwall = 0;
  proj_1D(&sigNwall, &bkgNwall,    "Nwall", ";number of side walls hit;slice counts", 5, -0.5, 4.5, "(nc5t>1)+(nc5b>1)+(nc5l>1)+(nc5r>1)", gCutNHIT);
  plot_1D( sigNwall,  bkgNwall, "02_Nwall", "portrait", logy, 1.0, -99, 2.5);
  }
  TCut cutNwall("((nc5t>1)+(nc5b>1)+(nc5l>1)+(nc5r>1))<=2");

  //
  // Put this level of selection into a single cut to be used later.
  //
  gCutVETO = TCut(cutNp5&&cutNwall);
  //====================================================================

  /*====================================================================
    Select verticies inside a fiducial region.
  */
  if (gPlotIt) {
  TH1F* sigVx = 0;
  TH1F* bkgVx = 0;
  proj_1D(&sigVx, &bkgVx,    "Vx", ";vertex x (cm);slice counts/ 10 cm", 40, -200, 200, "x", gCutNHIT&&gCutVETO);
  plot_1D( sigVx,  bkgVx, "03_Vx", "portrait", liny, 20, -160, 160);

  TH1F* sigVy = 0;
  TH1F* bkgVy = 0;
  proj_1D(&sigVy, &bkgVy,    "Vy", ";vertex y (cm);slice counts/ 10 cm", 40, -200, 200, "y", gCutNHIT&&gCutVETO);
  plot_1D( sigVy,  bkgVy, "03_Vy", "portrait", liny, 20, -160, 160);

  TH1F* sigVz = 0;
  TH1F* bkgVz = 0;
  proj_1D(&sigVz, &bkgVz,    "Vz", ";vertex z (cm);slice counts/ 50 cm", 28, 0, 1400, "z", gCutNHIT&&gCutVETO);
  plot_1D( sigVz,  bkgVz, "03_Vz", "landscape", liny, 20, 50, 1000);
  }
  //
  // Summary of vertex cut
  //
  gCutVTX = TCut("abs(x)<=160&&abs(y)<=160&&z>=50&&z<=1000");
  //====================================================================

  /*====================================================================
    Make some assessment of the overall shape of the event.  An
    event with two muon tracks and no hadronic activity should have ~2
    hits per plane consistent along its length. Other neutrino events
    will be "front loaded" with hadronic activity.
  */
  if (gPlotIt) {
  TH1F* sigHperP = 0;
  TH1F* bkgHperP = 0;
  proj_1D(&sigHperP, &bkgHperP,    "HperP", ";Nhit/Nplane;slice counts/0.2", 25, 0, 5, "(nhitx+nhity)/np", gCutNHIT&&gCutVETO&&gCutVTX);
  plot_1D( sigHperP,  bkgHperP, "04_HperP", "portrait", liny, 450, -99, 2.5);
  }
  TCut cutHperP("(nhitx+nhity)/np<=2.5");

  if (gPlotIt) {
  TH1F* sigNc = 0;
  TH1F* bkgNc = 0;
  proj_1D(&sigNc, &bkgNc,    "NcNp", ";(N^{2}_{cx}+N^{2}_{cy})^{1/2}/N_{planes};slice counts/0.1", 25, 0, 2.5, "sqrt(ncx**2+ncy**2)/np", gCutNHIT&&gCutVETO&&gCutVTX);
  plot_1D( sigNc,  bkgNc, "04_NcNp", "portrait", logy, 0.5, -99, 1.2);
  }
  TCut cutNcNp("sqrt(ncx**2+ncy**2)/np<=1.2");

  TH1F* sigN10 = 0;
  TH1F* bkgN10 = 0;
  if (gPlotIt) {
  proj_1D(&sigN10, &bkgN10,    "N10", ";number of hits in first ten planes;slice counts", 61, -0.5, 60.5, "n10", gCutNHIT&&gCutVETO&&gCutVTX);
  plot_1D( sigN10,  bkgN10, "04_N10", "portrait", logy, 0.5, -99, 25.5);
  }
  TCut cutN10("n10<=25");

  if (gPlotIt) {
  TH1F* sigPks = 0;
  TH1F* bkgPks = 0;
  proj_1D(&sigPks, &bkgPks,    "Pks", ";log_{10}(K-S probability);slice counts/0.2", 40, -8, 0, "log10(pks)", gCutNHIT&&gCutVETO&&gCutVTX);
  plot_1D( sigPks,  bkgPks, "04_Pks", "portrait", logy, 0.5, -4.5, -99);
  }
  TCut cutPks("log10(pks)>=-4.5");
  //
  // Put all the shape cuts together...
  //
  gCutSHAPE = TCut(cutHperP&&cutNcNp&&cutN3p&&cutN2p&&cutN10&&cutPks);
  //====================================================================

  //
  // Make a plot of the N plane distribution we started with now with
  // all cuts just to close out the counts
  //
  gCutPASS1 = TCut(gCutNHIT&&gCutVETO&&gCutVTX&&gCutSHAPE);
  if (gPlotIt) {
  TH1F* sigNpFinal = 0;
  TH1F* bkgNpFinal = 0;
  proj_1D(&sigNpFinal, &bkgNpFinal,    "NpFinal", ";number of hit planes;slice counts/5 planes", 44, 0, 220, "np", gCutPASS1);
  plot_1D( sigNpFinal,  bkgNpFinal, "05_NpFinal", "portrait", logy, 1.0, 20, 200);
  }
}

//......................................................................
//
// For each level of cuts, print out some events numbers and relevant
// information for signal and background events that pass and
// fail. Use these to check that things are working as we expect by
// looking at the events in the event display.
//
void scan(const char* sample,   // Signal or background?
          const char* cut_name, // What cut are we inspecting?
          TTree*      nt,       // Which ntuple?
          int         n0,       // Where to start the scan
          int         n1,       // How many entries to scan for "passed"
          int         n2,       // How many entries to scan for "failed"
          const char* var,      // Variables to print
          TCut        c,        // Base set of cuts
          TCut        s)         // The cut to toggle on and off
{
  //
  // Show all rows, no prompt to continue, during Scan()
  //
  nt->SetScanField(0);

  std::cout << "|== " << sample << " passing " << cut_name << std::endl;
  nt->Scan(var, c&&s, "", n1, n0);
  
  std::cout << "|== " << sample << " failing " << cut_name << std::endl;
  nt->Scan(var, c&&!s, "", n2, n0);  
}
void nt_scan()
{
  TCut noCut("1");

  scan("SIGNAL",     "NHIT", ntSIG, 0, 50, 50, "run:event:slice:np:nhitx:nhity:n2p:n3p", noCut, gCutNHIT);
  scan("BACKGROUND", "NHIT", ntBKG, 0, 50, 50, "run:event:slice:np:nhitx:nhity:n2p:n3p", noCut, gCutNHIT);

  scan("SIGNAL",     "VETO", ntSIG, 50, 30,  30, "run:event:slice:np5:nc5t:nc5b:nc5l:nc5r", gCutNHIT, gCutVETO);
  scan("BACKGROUND", "VETO", ntBKG, 50, 400, 30, "run:event:slice:np5:nc5t:nc5b:nc5l:nc5r", gCutNHIT, gCutVETO);

  scan("SIGNAL",     "VTX", ntSIG, 500, 40,  400,  "run:event:slice:x:y:z", gCutNHIT&&gCutVETO, gCutVTX);
  scan("BACKGROUND", "VTX", ntBKG, 500, 400, 1000, "run:event:slice:x:y:z", gCutNHIT&&gCutVETO, gCutVTX);

  scan("SIGNAL",     "SHAPE", ntSIG, 1000, 40,   5000, "run:event:slice:nhitx:nhity:np:ncx:ncy:log10(pks)", gCutNHIT&&gCutVETO&&gCutVTX, gCutSHAPE);
  scan("BACKGROUND", "SHAPE", ntBKG, 1000, 1000, 2000, "run:event:slice:nhitx:nhity:np:ncx:ncy:log10(pks)", gCutNHIT&&gCutVETO&&gCutVTX, gCutSHAPE);
}

//......................................................................

void pass1_plots()
{
  gStyle->SetOptStat(0);
  gErrorIgnoreLevel = kWarning;
  
  open_files();
  make_plots();
  nt_scan();
} 

////////////////////////////////////////////////////////////////////////