cana_time_peak.C

Go to the documentation of this file.

#include <iostream>
#include <cmath>
#include "TArrow.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TCut.h"
#include "TH1F.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TStyle.h"
#include "TText.h"
#include "TLegend.h"
#include "TLatex.h"
//......................................................................
// Usage:
// First load a chain ch of cana ntuples using something like:
//
// root [0] .x cana.C
//
// Then load this plot macro, project the histograms and make the plots
//
// root [1] .L cana_time_peak.C+
// root [2] proj_all(ch);
// root [3] plot_all(ch);
//
// For plot tweaking it should be possible to repeat these steps:
//
// root [4] .L cana_time_peak.C
// root [5] plot_all(ch);
//
//......................................................................
// Cuts optimized for seeing a time peak. Maximize S/sqrt(B)
//
const char* label = "NDOS Data, 09/12/13 - 09/27/13";
TCut run("run>=15071&&run<=15131");
TCut activity("nhit>=25&&nhitx>=3&&nhity>=3");  

TCut y1cut("y1<150");
TCut x1cut("abs(x1)<100");
TCut z1cut("280<z1 && z1<1200");
TCut fiducial(x1cut&&y1cut&&z1cut);

TCut cosnumi("cosnumi>0.95");
TCut canacuts(run&&fiducial&&activity&&cosnumi);

//......................................................................
// Where is the time peak?
//
double t1       = 217.0;
double t2       = 227.8;
TCut intime("tave>=217.0 && tave<=227.8");
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)");
double intimet  = (t2-t1);
double outtimet = 400.0;

void NDOSLabel()
{
  TLatex* tx = new TLatex(0.1, .93, label);
  tx->SetTextColor(kBlack);
  tx->SetNDC();
  tx->SetTextSize(1/30.);
  tx->SetTextAlign(12);
  tx->Draw();
}

// Put a "NOvA Preliminary" tag in the corner
void Preliminary()
{
  TLatex* prelim = new TLatex(.9, .93, "NOvA Preliminary");
  prelim->SetTextColor(kBlue);
  prelim->SetNDC();
  prelim->SetTextSize(1/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

void CenterTitles(TH1* histo)
{
  histo->GetXaxis()->CenterTitle();
  histo->GetYaxis()->CenterTitle();
  histo->GetZaxis()->CenterTitle();  
}

//......................................................................
// Make plots documenting the cuts used
//
TH1F* nhith;
void nhit_proj(TChain& ch) 
{
  nhith = new TH1F("nhith",";N_{hit};Events",101,-0.5,100.5);
  nhith->SetFillColor(kGreen-10);
  ch.Project("nhith","nhit", run);
}

void nhit_plot()
{
  gStyle->SetOptStat(0);

  TCanvas* nhitc = new TCanvas("nhitc","nhitc",800,600);
  nhitc->SetLogy();
  
  nhith->SetFillColor(kGreen-10);
  CenterTitles(nhith);
  nhith->Draw();
  Preliminary();

  TLine* nhitl = new TLine(25.0, 0, 25.0, nhith->GetMaximum());
  nhitl->SetLineStyle(3);
  nhitl->Draw();
  double ya1 = 0.75*log10(nhith->GetMaximum());
  double ya2 = 0.75*log10(nhith->GetMaximum())+0.9;
  double ya3 = pow(10,0.5*(ya1+ya2));
  ya1 = pow(10,ya1);
  ya2 = pow(10,ya2);
  TArrow* nhita1 = new TArrow(25.0, ya1, 30.0, ya1, 0.03);
  TArrow* nhita2 = new TArrow(25.0, ya2, 30.0, ya2, 0.03);
  nhita1->Draw();
  nhita2->Draw();

  TText* tx2 = new TText(30.0, ya3, "selected");
  tx2->SetTextFont(42);
  tx2->SetTextAlign(12);
  tx2->SetTextSize(0.035);
  tx2->Draw();

  NDOSLabel();

  nhitc->SaveAs("cana_nhit.pdf");
  nhitc->SaveAs("cana_nhit.png");
}

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

TH1F* y1h;
void y1_proj(TChain& ch) 
{
  y1h = new TH1F("y1h",";upper end y location (cm); Events / cm",
                       500,-250,250);
  ch.Project("y1h","y1", run&&activity);
}
void y1_plot() 
{
  gStyle->SetOptStat(0);

  TCanvas* y1c = new TCanvas("y1c","y1c",800,600);
  
  y1h->SetFillColor(kGreen-10);
  CenterTitles(y1h);
  y1h->Draw();
  Preliminary();

  TLine* y1l = new TLine(150.0, 0, 150.0, y1h->GetMaximum());
  y1l->SetLineStyle(3);
  y1l->Draw();
  double ya1 = 0.66*y1h->GetMaximum();
  double ya2 = 0.60*y1h->GetMaximum();
  double ya3 = 0.5*(ya1+ya2);
  
  TArrow* y1a1 = new TArrow(150.0, ya1, 120.0, ya1, 0.03);
  TArrow* y1a2 = new TArrow(150.0, ya2, 120.0, ya2, 0.03);
  y1a1->Draw();
  y1a2->Draw();

  NDOSLabel();

  TText* tx2 = new TText(110.0,ya3,"selected");
  tx2->SetTextAlign(32);
  tx2->SetTextFont(42);
  tx2->SetTextSize(0.035);
  tx2->Draw();

  y1c->SaveAs("cana_y1c.pdf");
  y1c->SaveAs("cana_y1c.png");
}

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

TH1F* x1h;
void x1_proj(TChain& ch) 
{
  x1h = new TH1F("x1h",";upper end x location (cm); Events / cm",
                 500,-250,250);
  ch.Project("x1h","x1", run&&activity&&y1cut);
}
void x1_plot() 
{
  gStyle->SetOptStat(0);
  
  TCanvas* x1c = new TCanvas("x1c","x1c",800,600);
  
  x1h->SetFillColor(kGreen-10);
  CenterTitles(x1h);
  x1h->Draw();
  Preliminary();
  
  TLine* x1l1 = new TLine(-100.0, 0, -100.0, x1h->GetMaximum());
  TLine* x1l2 = new TLine( 100.0, 0,  100.0, x1h->GetMaximum());
  x1l1->SetLineStyle(3);
  x1l2->SetLineStyle(3);
  x1l1->Draw();
  x1l2->Draw();
  double ya1 = 0.66*x1h->GetMaximum();
  double ya2 = 0.60*x1h->GetMaximum();
  double ya3 = 0.5*(ya1+ya2);
  
  TArrow* x1a1 = new TArrow(-100.0, ya1, -80.0, ya1, 0.02);
  TArrow* x1a2 = new TArrow(-100.0, ya2, -80.0, ya2, 0.02);
  x1a1->Draw();
  x1a2->Draw();

  TText* tx2 = new TText(0.0,ya3,"selected");
  tx2->SetTextFont(42);
  tx2->SetTextAlign(22);
  tx2->SetTextSize(0.035);
  tx2->Draw();

  TArrow* x1a3 = new TArrow(100.0, ya1, 80.0, ya1, 0.025);
  TArrow* x1a4 = new TArrow(100.0, ya2, 80.0, ya2, 0.025);
  x1a3->Draw();
  x1a4->Draw();

  NDOSLabel();

  x1c->SaveAs("cana_x1c.pdf");
  x1c->SaveAs("cana_x1c.png");
}

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

TH1F* z1h;
void z1_proj(TChain& ch) 
{
  z1h = new TH1F("z1h",";upper end z location (cm); Events / 10 cm",
                 140,0,1400);
  ch.Project("z1h","z1",run&&activity&&y1cut&&x1cut);
}
void z1_plot() 
{
  gStyle->SetOptStat(0);

  TCanvas* z1c = new TCanvas("z1c","z1c",800,600);

  z1h->SetFillColor(kGreen-10);
  CenterTitles(z1h);
  z1h->Draw();
  Preliminary();

  TLine* z1l1 = new TLine( 280.0, 0,  280.0, z1h->GetMaximum());
  TLine* z1l2 = new TLine(1200.0, 0, 1200.0, z1h->GetMaximum());
  z1l1->SetLineStyle(3);
  z1l2->SetLineStyle(3);
  z1l1->Draw();
  z1l2->Draw();
  double ya1 = 0.66*z1h->GetMaximum();
  double ya2 = 0.60*z1h->GetMaximum();
  double ya3 = 0.5*(ya1+ya2);
  
  TArrow* z1a1 = new TArrow(280.0, ya1, 330.0, ya1, 0.025);
  TArrow* z1a2 = new TArrow(280.0, ya2, 330.0, ya2, 0.025);
  z1a1->Draw();
  z1a2->Draw();

  TText* tx2 = new TText(700.0,ya3,"selected");
  tx2->SetTextFont(42);
  tx2->SetTextAlign(22);
  tx2->SetTextSize(0.035);
  tx2->Draw();

  TArrow* z1a3 = new TArrow(1200.0, ya1, 1150.0, ya1, 0.025);
  TArrow* z1a4 = new TArrow(1200.0, ya2, 1150.0, ya2, 0.025);
  z1a3->Draw();
  z1a4->Draw();

  NDOSLabel();

  z1c->SaveAs("cana_z1c.pdf");
  z1c->SaveAs("cana_z1c.png");
}

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

TH1F* cosnumiinh;
TH1F* cosnumiouth;
void cosnumi_proj(TChain& ch) 
{
  cosnumiinh = new TH1F("cosnumiinh",
                        ";cos #theta_{NuMI}; Events / 0.05",
                        20, 0.0, 1.0);
  cosnumiouth = new TH1F("cosnumiouth",
                         ";cos #theta_{NuMI}; Events / 0.05",
                         20, 0.0, 1.0);
  ch.Project("cosnumiinh", "cosnumi",run&&activity&&fiducial&&intime);
  ch.Project("cosnumiouth","cosnumi",run&&activity&&fiducial&&outtime);

  cosnumiouth->Scale(intimet/outtimet);
}
void cosnumi_plot() 
{
  gStyle->SetOptStat(0);
  
  TCanvas* cosnumic = new TCanvas("cosnumic","cosnumic",800,600);

  cosnumiinh->SetMarkerStyle(20);
  cosnumiouth->SetLineColor(kRed);
  cosnumiouth->SetLineWidth(2);

  cosnumiouth->SetMaximum(cosnumiinh->GetMaximum()+
                          1.5*sqrt(cosnumiinh->GetMaximum()));

  CenterTitles(cosnumiouth);
  cosnumiouth->Draw("hist");
  cosnumiinh-> Draw("same,e");

  Preliminary();

  double x1 = cosnumiouth->GetBinCenter(4);
  double y1 = cosnumiouth->GetBinContent(4);
  
  TLine* a1 = new TLine(x1+0.08, 1.12*y1, x1+0.03, 1.01*y1);
  a1->SetLineWidth(1);
  a1->SetLineColor(kRed);
  a1->Draw();

  TText* tx1 = new TText(x1+0.09, 1.14*y1,"out of time tracks");
  tx1->SetTextColor(kRed);
  tx1->Draw();

  x1 = cosnumiinh->GetBinCenter(18);
  y1 = cosnumiinh->GetBinContent(18);

  TLine* a2 = new TLine(x1-0.04, 1.01*y1, x1-0.08, 1.12*y1);
  a2->SetLineWidth(1);
  a2->SetLineColor(kBlack);
  a2->Draw();

  TText* tx2 = new TText(x1-0.06, 1.20*y1,"in time tracks");
  tx2->SetTextAlign(31);
  tx2->SetTextColor(kBlack);
  tx2->Draw();

  NDOSLabel();

  cosnumic->SaveAs("cana_cosnumic.pdf");
  cosnumic->SaveAs("cana_cosnumic.png");
}

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

TH1F* tavewideh;
void time_wide_proj(TChain& ch) 
{
  tavewideh = new TH1F("tavewideh",
                       ";NDOS event time [#musec];Events / 2 #musec",
                       250,0,500);
  ch.Project("tavewideh", "tave", canacuts);
}
void time_wide_plot()
{  
  gStyle->SetOptStat(0);
  
  TCanvas* cana_tave_wide = 
    new TCanvas("cana_tave_wide","cana_tave_wide",800,600);
  tavewideh->SetFillColor(kGreen-10);
  CenterTitles(tavewideh);
  tavewideh->SetMaximum(floor(tavewideh->GetMaximum())+1);
  tavewideh->Draw();
  Preliminary();
  
  unsigned int i1 = tavewideh->FindBin(t1);
  unsigned int i2 = tavewideh->FindBin(t2);
  double splusb = tavewideh->Integral(i1,i2);
  double b = 
    tavewideh->Integral(i1-25,i1-1) + 
    tavewideh->Integral(i2+1,i2+25);
  b *= (float)(i2-i1+1)/50.0;
  double fom = (splusb-b)/sqrt(b);
  std::cout << "FOM: " << fom << std::endl;

  NDOSLabel();

  /*
  TText* tx4 = new TText(240, 15, "Best fit for NuMI beam spill window");
  tx4->SetTextFont(32);
  tx4->SetTextSize(0.035);
  tx4->SetTextColor(kRed);
  tx4->Draw();
  */

  cana_tave_wide->SaveAs("cana_tave_wide.pdf");
  cana_tave_wide->SaveAs("cana_tave_wide.png");
}

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

TH1F*  tavezoomh;
int    ntbin_per_usec = 20;
double zoom_tbin_sz   = 1.0/(float)ntbin_per_usec;
double zoom_tlo  = floor(t1) - 5;
double zoom_thi  = floor(t2) + 6;
int    zoom_nbin = ntbin_per_usec*(zoom_thi-zoom_tlo);
void time_zoom_proj(TChain& ch)
{  
  gStyle->SetOptStat(0);
  
  tavezoomh = new TH1F("tavezoomh",
                       ";NDOS event time [#musec];Events / 50 nsec",
                       zoom_nbin,
                       zoom_tlo,
                       zoom_thi);
  ch.Project("tavezoomh", "tave", canacuts);
}

//......................................................................
// The NuMI spill time structure for time t in usec. Based on e-mail
// exchange with Phil Adamson. This applies to the period starting
// September 2013 -- ?
//
static const int nbatch = 6;
static const double dt_batch = 81*18.83E-3;
static const double dt_gap   =  5*18.83E-3;
static const double dt_cycle = dt_batch + dt_gap;
static const double spill_length = nbatch*dt_batch + (nbatch-1)*dt_gap;

double spill_profile(double t) 
{
  if (t<0.0)          return 0.0;
  if (t>spill_length) return 0.0;
  
  // Rewind cycle by cycle
  double ttmp = t;
  while (ttmp>dt_cycle) ttmp -= dt_cycle;
  
  // Check it we are inside the "on" part of the cycle
  if (ttmp<0.0)      return 0.0;
  if (ttmp>dt_batch) return 0.0;

  return 1.0;
}

//......................................................................
// Convolve the spill profile with a gaussian to get a template of what
// the spill profile should look like
//
TH1F* template_spillh = 0;
double template_tlo;
double template_thi;
int    template_nbin;
void make_template(double sigma) 
{
  if (template_spillh!=0) {
    delete template_spillh;
    template_spillh = 0;
  }
  if (template_spillh==0) {
    template_tlo =
      zoom_tbin_sz*floor(-5.0*sigma/zoom_tbin_sz);
    template_thi =
      zoom_tbin_sz*
      floor((spill_length+5.0*sigma)/zoom_tbin_sz+1);
    template_nbin = 10*(template_thi-template_tlo)/zoom_tbin_sz;
    template_spillh = new TH1F("tempalte_spillh",";usec",
                               template_nbin,
                               template_tlo,
                               template_thi);
  }
  
  int i, j;
  double ti, tj;
  double s, w;
  for (i=1; i<=template_nbin; ++i) {
    ti = template_spillh->GetBinCenter(i);
    s = spill_profile(ti);
    for (j=1; j<template_nbin; ++j) {
      tj = template_spillh->GetBinCenter(j);
      w = TMath::Gaus(ti,tj,sigma);
      
      template_spillh->Fill(tj, w*s);
    }
  }
  template_spillh->Scale(template_spillh->GetMaximum());
}

//......................................................................
// Make the fit function for the spill profile. b is the number of
// background evnets per bin, p is a scale factor on the signal
// strength, and t0 is where the spill should start.
TH1F* time_profileh;
void make_time_profile(double b,
                       double f,
                       double t0,
                       double sigma) 
{
  int i;
  
  static int ifirst = 1;
  if (ifirst) {
    time_profileh = new TH1F("time_profileh",
                             ";usec",
                             zoom_nbin,
                             zoom_tlo,
                             zoom_thi);
    ifirst = 0;
  }
  time_profileh->Reset();

  static double last_sigma = 0.0;
  if (last_sigma != sigma) {
    last_sigma = sigma;
    make_template(sigma);
  }
  
  double s, t;
  //
  // Put the background in
  //
  for (i=1; i<=time_profileh->GetNbinsX(); ++i) {
    t = time_profileh->GetBinCenter(i);
    time_profileh->Fill(t, b);
  }
  //
  // Put the signal in
  //
  for (i=1; i<=template_nbin; ++i) {
    t = template_spillh->GetBinCenter(i);
    s = template_spillh->GetBinContent(i);
    time_profileh->Fill(t+t0, f*s);
  }
}

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

double chi2()
{
  int i;
  double sum = 0;
  double ni, nui;
  for (i=1; i<=zoom_nbin; ++i) {
    nui = time_profileh->GetBinContent(i);
    ni  = tavezoomh->    GetBinContent(i);
    sum += nui;
    if (ni>0.0) {
      sum += ni*(log(ni/nui)-1);
    }
  }
  return 2.0*sum;
}

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

void fcn(int&, double*, double& f, double* p, int)
{
  make_time_profile(p[0], p[1], p[2], p[3]);
  f = chi2();
}

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

void minimize() 
{
  int ierr;
  TMinuit minuit;
  minuit.SetPrintLevel(-99);
  minuit.SetFCN(fcn);
  minuit.mnparm(0, "bkg",  1e-3, 1e-4, 0, 10.0, ierr);
  minuit.mnparm(1, "sf",   1e-3, 1e-4, 0, 10.0, ierr);
  minuit.mnparm(2, "t0",   1e-3, 1e-4, 0, 10.0, ierr);
  minuit.mnparm(3, "sigt", 1e-3, 1e-4, 0, 10.0, ierr);
  
  double chi2best    = 99E99;
  double p[4];
  double sigp[4];
  
  int i;
  double t0;
  double x[zoom_nbin];
  double y[zoom_nbin];
  for (i=0; i<zoom_nbin; ++i) {
    t0 = tavezoomh->GetBinCenter(i+1);

    minuit.mnparm(2, "t0", t0, 1e-3, 0.0, 500.0, ierr);
    minuit.FixParameter(2);

    minuit.mnparm(3, "sigt", 0.512/sqrt(12)/sqrt(25), 1e-3, 0.0, 500.0, ierr);
    minuit.FixParameter(3);
    
    minuit.Migrad();
    
    double fmin, fedm, errdef;
    int npari, nparx, istat;
    minuit.mnstat(fmin,fedm,errdef,npari,nparx,istat);
    x[i] = t0;
    y[i] = fmin;
    std::cerr << i << "\t" << t0 << "\t" << fmin << std::endl;
    if (fmin<chi2best) {
      chi2best = fmin;
      
      double val, err, xlolim, xuplim;
      int iuint;

      static TString bkg("bkg");
      minuit.mnpout(0, bkg, val, err, xlolim, xuplim, iuint);
      p[0]    = val;
      sigp[1] = err;

      static TString sf("sf");
      minuit.mnpout(1, sf, val, err, xlolim, xuplim, iuint);
      p[1]    = val;
      sigp[1] = err;
      
      static TString t0s("t0");
      minuit.mnpout(2, t0s, val, err, xlolim, xuplim, iuint);
      p[2]    = val;
      sigp[2] = err;

      static TString sigt("sigt");
      minuit.mnpout(3, sigt, val, err, xlolim, xuplim, iuint);
      p[3]    = val;
      sigp[3] = err;
    }
  }
  make_time_profile(p[0], p[1], p[2], p[3]);

  std::cerr << "chi2 = " << chi2best << std::endl;
  std::cerr << "Best b   = " << p[0] <<" +/- " << sigp[0] << std::endl;
  std::cerr << "Best sf  = " << p[1] <<" +/- " << sigp[1] << std::endl;
  std::cerr << "Best t0  = " << p[2] <<" +/- " << sigp[2] << std::endl;
  std::cerr << "Best sigt= " << p[3] <<" +/- " << sigp[3] << std::endl;
}

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

void time_zoom_plot() 
{
  TCanvas* time_zoom_drawc = 
    new TCanvas("time_zoom_drawc",
                "time_zoom_drawc",800,600);
  tavezoomh->SetFillColor(kGreen-10);
  CenterTitles(tavezoomh);
  tavezoomh->SetLineWidth(1);
  tavezoomh->Draw();
  Preliminary();

  if (time_profileh) {
    time_profileh->SetLineWidth(2);
    time_profileh->SetLineColor(kRed);
    time_profileh->Draw("same");
  }

  TLegend *tleg = new TLegend(0.15,0.75,0.43,0.85);
  tleg->SetTextFont(42);
  tleg->SetTextSize(0.04);
  tleg->AddEntry(time_profileh,"Fit to NuMI Beam Structure");
  tleg->Draw();
  
  NDOSLabel();

  time_zoom_drawc->SaveAs("cana_time_zoom.png");
  time_zoom_drawc->SaveAs("cana_time_zoom.pdf");
}

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

void proj_all(TChain& ch) 
{

  std::cout<<"starting proj"<<std::endl;
  nhit_proj(ch);
  std::cout<<"nhit done"<<std::endl;
  y1_proj(ch);
  std::cout<<"y1 done"<<std::endl;
  x1_proj(ch);
  std::cout<<"x1 done"<<std::endl;
  z1_proj(ch);
  std::cout<<"z1 done"<<std::endl;
  cosnumi_proj(ch);
  std::cout<<"cosnumi done"<<std::endl;
  time_wide_proj(ch);
  std::cout<<"time wide done"<<std::endl;
  time_zoom_proj(ch);
  std::cout<<"time zoom done"<<std::endl;
  minimize();
  std::cout<<"minimize done"<<std::endl;
}

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

void plot_all() 
{
  nhit_plot();
  y1_plot();
  x1_plot();
  z1_plot();
  cosnumi_plot();
  time_wide_plot();
  time_zoom_plot();
}



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