accum_nue_numu_equivalent_fhc.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "TCanvas.h"
#include "TH2.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TFile.h"
#include "TProfile.h"
#include "TString.h"
#include "TPaveText.h"
#include "TGaxis.h"
#include "TROOT.h"
#include "TStyle.h"

#include "/nova/app/users/howard/Dev_Nue2018_Commit/CAFAna_includes/fhc_results.h"

#include "Utilities/rootlogon.C"

using namespace ana;

int countEvents(int iRun, const long eventList[], int totalEvents);
int countFakeEvents(int iRun, long eventList[], int totalEvents);
double ksTest(TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials);

void accum_nue_numu_equivalent_fhc()
{
  // Following Evan's guidelines more or less for time axis
  gSystem->Setenv("TZ","UTC");
  gStyle->SetTimeOffset(0);
  std::string taxis_labels = "%b '%y";

  // File with the exposure, livetime, POT plots
  TFile *fAccounting = new TFile("/nova/ana/users/howard/Ana2018-pot/nue_ana2018_pot-fhc.root","read");
  TH1D* hPOT = (TH1D*)fAccounting->Get("hPOTMassT"); // switch to using kT*POT
  hPOT->Scale(1./10.43); // scale by 1./10.43kT to get equivalent POT (see docdb-23417)
  int timeNBins = hPOT->GetXaxis()->GetNbins();
  double loTime = hPOT->GetXaxis()->GetBinLowEdge(1);
  double dxTime = hPOT->GetXaxis()->GetBinLowEdge(2) - hPOT->GetXaxis()->GetBinLowEdge(1);
  double hiTime = hPOT->GetXaxis()->GetBinLowEdge(timeNBins) + dxTime;

  // Set up time axis
  hPOT->SetTitle("");
  hPOT->GetXaxis()->SetTitle("Time");
  CenterTitles(hPOT);
  hPOT->GetXaxis()->SetTimeDisplay(1);
  hPOT->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
  hPOT->GetXaxis()->SetTitle("Date");
  hPOT->GetXaxis()->SetNdivisions(506,kFALSE);

  // Now make a plot of accumulated POT and accumulated event as function of time
  // ---------------------------------------------------------------------------------
  gStyle->SetPadTickY(0);
  std::cout << "Making new histogram with min,max,nbins = " << loTime << "," << hiTime << "," << timeNBins << std::endl;
  TH1D *hPOTAccumulated = new TH1D("hPOTAccumulated",";Time;Accumulated equivalent POT (x10^{20})",timeNBins,loTime,hiTime);
  TH1D *hEvtAccumulated = new TH1D("hEvtAccumulated",";Time;Accumulated Events",timeNBins,loTime,hiTime);
  TH1D *hEvtAccumulatedNumu = new TH1D("hEvtAccumulatedNumu",";Time;Accumulated Events",timeNBins,loTime,hiTime);

  for( int iRBin = 1; iRBin <= timeNBins; ++iRBin ){
    if( iRBin==1 )
      hPOTAccumulated->SetBinContent( iRBin, hPOT->GetBinContent(iRBin)/1.0e20 );
    else
      hPOTAccumulated->SetBinContent( iRBin, hPOTAccumulated->GetBinContent(iRBin-1)+(hPOT->GetBinContent(iRBin)/1.0e20) );
    hEvtAccumulated->SetBinContent( iRBin, countEvents(hEvtAccumulated->GetBinLowEdge(iRBin)+dxTime,nueTime,nEvts) );
    hEvtAccumulatedNumu->SetBinContent( iRBin, countEvents(hEvtAccumulatedNumu->GetBinLowEdge(iRBin)+dxTime,numuTime,nEvtsNumu) );
  }

  TCanvas *tcAccumulated = new TCanvas();

  // Blank plot to draw on
  TH1D *hPOTAccumBlank = new TH1D("hPOTAccumBlank",";Time;Accumulated equivalent POT (x10^{20})",timeNBins,loTime,hiTime);
  hPOTAccumBlank->GetXaxis()->SetRangeUser(loTime,1490000000);
  hPOTAccumBlank->GetYaxis()->SetRangeUser(0.,1.1*hPOTAccumulated->GetMaximum());
  CenterTitles(hPOTAccumBlank);

  hPOTAccumBlank->GetXaxis()->SetTimeDisplay(1);
  hPOTAccumBlank->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
  hPOTAccumBlank->GetXaxis()->SetTitle("Date");
  hPOTAccumBlank->GetXaxis()->SetNdivisions(506,kFALSE);

  hPOTAccumulated->GetXaxis()->SetTimeDisplay(1);
  hPOTAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
  hPOTAccumulated->GetXaxis()->SetTitle("Date");
  hPOTAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);

  hPOTAccumBlank->Draw("hist");
  hPOTAccumulated->Draw("hist same");

  // Set up the right axis and overlay the events
  tcAccumulated->Update();
  Float_t rightmax = 1.1*hEvtAccumulated->GetMaximum();
  Float_t rightmaxNumu = 1.1*hEvtAccumulatedNumu->GetMaximum();
  Float_t scale = gPad->GetUymax()/rightmax;
  Float_t scaleNumu = gPad->GetUymax()/rightmaxNumu;
  hEvtAccumulated->SetLineColor(kGreen+3);
  hEvtAccumulated->SetMarkerColor(kGreen+3);
  hEvtAccumulated->Scale(scale);

  hEvtAccumulated->GetXaxis()->SetTimeDisplay(1);
  hEvtAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
  hEvtAccumulated->GetXaxis()->SetTitle("Date");
  hEvtAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);

  hEvtAccumulated->Draw("hist same");
  hEvtAccumulatedNumu->SetLineColor(kRed);
  hEvtAccumulatedNumu->SetMarkerColor(kRed);
  hEvtAccumulatedNumu->Scale(scaleNumu);
  hEvtAccumulatedNumu->Draw("hist same");
  gStyle->SetLineWidth(3);
  TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(), gPad->GetUxmax(),gPad->GetUymax(),0,1.1,506,"+L");
  axis->SetLineColor(kBlack);
  axis->SetLabelColor(kBlack);
  // NOvA style axis options (copied from rootlogon file and set up to work)
  // --Axis titles
  axis->SetTitleSize(.055);
  axis->SetTitleOffset(.8);
  // --More space for y-axis to avoid clashing with big numbers
  axis->SetTitleOffset(.9);
  // --Axis labels (numbering)
  axis->SetLabelSize(.04);
  axis->SetLabelOffset(.005);
  // --Font
  axis->SetLabelFont(42);
  axis->SetTitleFont(42);
  // -----------------------------------------------------------------------
  axis->SetTitleColor(kBlack);
  axis->SetTitle("Fraction of Selected Events");
  axis->CenterTitle();
  axis->Draw();

  // Get ks test values
  double ksTestNue = ksTest(hPOTAccumulated,hEvtAccumulated,hPOT,nEvts,10000);
  double ksTestNumu = ksTest(hPOTAccumulated,hEvtAccumulatedNumu,hPOT,nEvtsNumu,10000);

  // Legend
  tcAccumulated->cd();
  TLegend *lAccum = new TLegend(0.12,0.55,0.55,0.87);
  lAccum->SetFillStyle(0);
  lAccum->AddEntry("hPOTAccumBlank","Accumulated POT","l");
  lAccum->AddEntry(hEvtAccumulated->GetName(),TString::Format("#nu_{e} candidates (KS prob = %.0f%%)",ksTestNue*100.),"l");
  lAccum->AddEntry(hEvtAccumulatedNumu->GetName(),TString::Format("#nu_{#mu} candidates (KS prob = %.0f%%)",ksTestNumu*100.),"l");
  lAccum->Draw();

  tcAccumulated->cd();
  TPaveText *beam = new TPaveText(0.55,0.13,0.75,0.2,"NB NDC");
  beam->AddText("Neutrino Beam");
  beam->Draw();

  Preliminary();
  tcAccumulated->Print("AccumulatedEvents_Ana2018_fhc_equivalent_time.pdf");
}

int countEvents( int iRun, const long eventList[], int totEvents ){
  int accumulatedEvts = 0;
  for( int iEvtCt=0; iEvtCt<totEvents; ++iEvtCt ){
    if( eventList[iEvtCt] < iRun ) accumulatedEvts++; // for time version we pass bin high edge, so want explicitly lower than.
  }
  return accumulatedEvts;
}

int countFakeEvents( int iRun, long eventList[], int totEvents ){
  int accumulatedEvts = 0;
  for( int iEvtCt=0; iEvtCt<totEvents; ++iEvtCt ){
    if( eventList[iEvtCt] < iRun ) accumulatedEvts++; // for time version we pass bin high edge, so want explicitly lower than.
  }
  return accumulatedEvts;
}

// Run pseudo-experiments with POT histograms (normalized here if needed) to test how likely the event distribution is
double ksTest( TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials ){
  if( nTrials==0 ) return -999.;

  // Set up histograms
  int nCDFBins = expCDF->GetXaxis()->GetNbins();
  double loBin = expCDF->GetXaxis()->GetBinLowEdge(1);
  double dxBin = expCDF->GetXaxis()->GetBinLowEdge(2) - expCDF->GetXaxis()->GetBinLowEdge(1);
  double hiBin = expCDF->GetXaxis()->GetBinLowEdge(nCDFBins);

  // Clone histograms to be acted upon
  TH1D *expCDFcpy = (TH1D*)expCDF->Clone();
  TH1D *expPDFcpy = (TH1D*)expPDF->Clone();
  TH1D *obsCDFcpy = (TH1D*)obsCDF->Clone();
  expCDFcpy->SetDirectory(0);
  expPDFcpy->SetDirectory(0);
  obsCDFcpy->SetDirectory(0);

  double expPDFcpyArea = expPDFcpy->Integral("width");
  expPDFcpy->Scale(1./expPDFcpyArea); // area normalized
  std::cout << "Check norm: " << fabs(expCDFcpy->GetMaximum()-1.0) << std::endl;
  if( fabs(expCDFcpy->GetMaximum()-1.0)>1.0e-9 ){
    // only normalize if not already normalized to have max at 1 (CDF)
    expCDFcpy->Scale(1.0e20/expPDF->Integral()); // factor 1e20 is because it's factored out of POT histogram
  }
  obsCDFcpy->Scale(1./obsCDFcpy->GetMaximum());

  std::cout << "Check PDF: " << expPDFcpy->Integral("width") << std::endl;

  // -----------------------------
  // Make some debug plots
  TCanvas *cExpPDF = new TCanvas();
  expPDFcpy->Draw("hist");
  if(nEvents==nEvts) cExpPDF->Print("debug_nue_expPDF_equivalent_time_fhc.pdf");
  else if(nEvents==nEvtsNumu) cExpPDF->Print("debug_numu_expPDF_equivalent_time_fhc.pdf");

  TCanvas *cExpCDF = new TCanvas();
  expCDFcpy->Draw("hist");
  if(nEvents==nEvts) cExpCDF->Print("debug_nue_expCDF_equivalent_time_fhc.pdf");
  else if(nEvents==nEvtsNumu) cExpCDF->Print("debug_numu_expCDF_equivalent_time_fhc.pdf");

  TCanvas *cObsCDF = new TCanvas();
  obsCDFcpy->Draw("hist");
  if(nEvents==nEvts) cObsCDF->Print("debug_nue_obsCDF_equivalent_time_fhc.pdf");
  else if(nEvents==nEvtsNumu) cObsCDF->Print("debug_numu_obsCDF_equivalent_time_fhc.pdf");
  // -----------------------------

  // Maximum difference between event CDF and POT CDF (value for KS test)
  double obsMaxDiff = 0.;
  double locMaxDiff = 0.;
  int    binMaxDiff = 0;
  for( int iBin=1; iBin<=nCDFBins; ++iBin ){
    if( fabs(expCDFcpy->GetBinContent(iBin)-obsCDFcpy->GetBinContent(iBin)) > obsMaxDiff){
      obsMaxDiff = fabs(expCDFcpy->GetBinContent(iBin)-obsCDFcpy->GetBinContent(iBin));
      locMaxDiff = expCDFcpy->GetBinCenter(iBin);
      binMaxDiff = iBin;
    }
  }

  // -----------------------------
  // Make some overlay debug plots
  TCanvas *cDBOverlay = new TCanvas();
  expCDFcpy->Draw("hist");
  obsCDFcpy->Draw("hist same");
  double minOfDiff = std::min(obsCDFcpy->GetBinContent(binMaxDiff),expCDFcpy->GetBinContent(binMaxDiff));
  double maxOfDiff = std::max(obsCDFcpy->GetBinContent(binMaxDiff),expCDFcpy->GetBinContent(binMaxDiff));
  TLine *overlayLine = new TLine(locMaxDiff,minOfDiff,locMaxDiff,maxOfDiff);
  overlayLine->SetLineColor(kRed);
  overlayLine->SetLineWidth(2);
  overlayLine->SetLineStyle(7);
  overlayLine->Draw();
  TPaveText *txtOverlay = new TPaveText(0.15,0.6,0.4,0.87,"NB NDC");
  txtOverlay->AddText( TString::Format("Max difference: %.3f",maxOfDiff-minOfDiff) );
  txtOverlay->Draw();
  if(nEvents==nEvts) cDBOverlay->Print("debug_nue_overlay_equivalent_time_fhc.pdf");
  else if(nEvents==nEvtsNumu) cDBOverlay->Print("debug_numu_overlay_equivalent_time_fhc.pdf");
  // -----------------------------

  std::cout << "Max Difference: " << obsMaxDiff << std::endl;

  unsigned int iPassTrial=0;

  TH1D *ksDistances;
  if( nEvents==nEvts ) ksDistances = new TH1D("ksDistancesNue",";max separation;fake trials",50,0,1);
  else if( nEvents==nEvtsNumu ) ksDistances = new TH1D("ksDistancesNumu",";max separation;fake trials",50,0,1);
  else return -999.;
  long fakeEvents[nEvtsNumu];
  TH1D *trialCDF;
  if( nEvents==nEvts ) trialCDF = new TH1D("trialCDFNue","",nCDFBins,loBin,hiBin+dxBin);
  else if( nEvents==nEvtsNumu ) trialCDF = new TH1D("trialCDFNumu","",nCDFBins,loBin,hiBin+dxBin);
  else return -999.;

  // Pseudo-experiments
  for( unsigned int iTrial=0; iTrial<nTrials; ++iTrial ){
    // Pick event times
    for( unsigned int iEvt=0; iEvt<nEvents; ++iEvt ){
      fakeEvents[iEvt] = expPDFcpy->GetRandom();
    }
    // Fill fake CDF
    for( int iRBin = 1; iRBin <= nCDFBins; ++iRBin ){
      trialCDF->SetBinContent( iRBin, countFakeEvents(trialCDF->GetBinCenter(iRBin),fakeEvents,nEvents) );
    }
    // Normalize CDF
    if( nEvents>0 ) trialCDF->Scale(1./double(nEvents));
    // Get the ks distance
    double trialMaxDiff = 0.;
    for( int iBin=1; iBin<=nCDFBins; ++iBin ){
      if( fabs(expCDFcpy->GetBinContent(iBin)-trialCDF->GetBinContent(iBin)) > trialMaxDiff)
        trialMaxDiff = fabs(expCDFcpy->GetBinContent(iBin)-trialCDF->GetBinContent(iBin));
    }
    ksDistances->Fill(trialMaxDiff);
    // If greater then observed max difference, count
    if(trialMaxDiff>=obsMaxDiff) iPassTrial+=1;
    // Reset histogram
    for( int iBin=1; iBin<=nCDFBins; ++iBin ){
      trialCDF->SetBinContent(iBin,0);
    }
  }

  // -----------------------------
  // Area normalize KS histogram and draw
  double ksArea = ksDistances->Integral();
  ksDistances->Scale(1./ksArea);

  TCanvas *cKSTest = new TCanvas();
  ksDistances->Draw();
  TLine *obsMaxDiffLine = new TLine(obsMaxDiff,0.,obsMaxDiff,0.95*ksDistances->GetMaximum());
  obsMaxDiffLine->SetLineColor(kRed);
  obsMaxDiffLine->Draw();
  if( nEvents==nEvts )
    cKSTest->Print("ksTestNue_equivalent_time_fhc.pdf");
  if( nEvents==nEvtsNumu )
    cKSTest->Print("ksTestNumu_equivalent_time_fhc.pdf");
  // -----------------------------

  // Result of the test is the fraction of pseudo-
  // experiments with KS distance >= observed KS distance
  return double(iPassTrial)/double(nTrials);
}