DataMCNDAna_nus17.C

Go to the documentation of this file.

// A macro for Near Detector N-1 analysis
//    cafe -bq DataMCNDAna_nus17.C <folder>
// where folder = the directory the input file lives
// and the output file will be saved too.
//
// Gavin S. Davies

#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "NuXAna/macros/NuSPlotFunctions.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TKey.h"
#include "TLegend.h"

#include <string>

using namespace ana;

void make_dataMC(TFile* SavedFile,
                 TFile* OutFile);

void make_plots(const ana::Spectrum& data,
                const ana::CheatDecomp& MC,
                std::string plotlabel,
                TFile* plotsOut);

PlotOptions GetPlotOptions(std::string name);

void DataMCNDAna_nus17(std::string folder)
{
  TH1::AddDirectory(0);
  // load in loaded specs
  std::string filein = folder + "/dataMCNDLoadOut.root";
  std::string fileout = folder + "/dataMCNDAnaOut.root";

  TFile* spSaved = new TFile(filein.c_str());
  TFile* plotsOut = new TFile(fileout.c_str(), "RECREATE");

  // make plots
  make_dataMC(spSaved, plotsOut);

  spSaved->Close();
  plotsOut->Close();
  std::cout << "plotted all" << std::endl;
}

//..............................................................................
void make_dataMC(TFile* SavedFile,
                 TFile* OutFile)
{
  TIter next(SavedFile->GetListOfKeys());
  TKey *key;
  while ((key=(TKey*)next()))
  {
    const char* fDatKey = key->GetName();
    std::string strDatKey(fDatKey);

    if(strDatKey.find("MC") != std::string::npos) { continue; }

    std::string strMCKey(fDatKey);
    strMCKey.replace(strMCKey.size()-4, 4, "MC");

    std::string fPlotLabel(fDatKey);
    fPlotLabel.erase(fPlotLabel.end()-4, fPlotLabel.end());

    std::cout << strDatKey << ", " << strMCKey << std::endl;
    std::cout << fPlotLabel << std::endl;

    CheatDecomp sMC = *CheatDecomp::LoadFrom(SavedFile, strMCKey.c_str());
    Spectrum sData = *Spectrum::LoadFrom(SavedFile, strDatKey.c_str());

    make_plots(sData, sMC, fPlotLabel.c_str(), OutFile);
  }
};

//..............................................................................
void make_plots(const ana::Spectrum& data,
                const ana::CheatDecomp& MC,
                std::string plotlabel,
                TFile* plotsOut)
{
  Spectrum empty = MC.NCTotalComponent();
  empty.Clear();

  Spectrum mcspecs[MAXSPEC] = {
    MC.NCTotalComponent()   +
    MC.NueComponent()  + MC.AntiNueComponent() +
    MC.NumuComponent() + MC.AntiNumuComponent(),
    MC.NCTotalComponent(),
    MC.NumuComponent() + MC.AntiNumuComponent(),
    MC.NueComponent()  + MC.AntiNueComponent(),
    empty,
    empty
  };

  double pot = data.POT();
  std::cout << "data POT: " << pot << std::endl;
  TH1* hDat    = data   .ToTH1(pot);
  TH1* hMCTot  = mcspecs[0].ToTH1(pot);
  TH1* hMCNC   = mcspecs[1].ToTH1(pot);
  TH1* hMCNuMu = mcspecs[2].ToTH1(pot);
  TH1* hMCNuE  = mcspecs[3].ToTH1(pot);

  std::string xLabel(hMCNC->GetXaxis()->GetTitle());
  if(xLabel.find("PT/P") != std::string::npos) {
    xLabel.replace(xLabel.find("PT/P"), 4, "p_{T}/p");
  }
  if(xLabel.find("RecoE") != std::string::npos) {
    xLabel = "Visible Energy (GeV)";
  }
  if(xLabel.find("RecoELong") != std::string::npos) {
    xLabel = "Visible Energy (GeV)";
  }

  std::string yLabel = "";
  if(plotlabel.find("RecoE") != std::string::npos) { yLabel = "10^{3} Events / 0.1 GeV / 8.04 #times10^{20} POT"; }
  if(plotlabel.find("RecoELong") != std::string::npos) { yLabel = "10^{3} Events / 0.1 GeV / 8.04 #times10^{20} POT"; }
  else yLabel = "10^{3} Events / 8.04 #times10^{20} POT";
  std::string histTitle = ";" + xLabel + ";" + yLabel;

  // Numu Stack = Numu
  TH1* hNumuStack = (TH1*)hMCNuMu->Clone();

  // Nue Stack = Numu Stack + Nue = Numu + Nue
  TH1* hNueStack  = (TH1*)hNumuStack->Clone();
  hNueStack->Add(hMCNuE);

  // NC Stack = Nue Stack + NC = Numu + Nue + NC
  TH1* hNCStack   = (TH1*)hNueStack->Clone();
  hNCStack->Add(hMCNC);

  double maxi = 0;
  double maxDat    = hDat   ->GetBinContent(hDat   ->GetMaximumBin());
  double maxMCTot  = hMCTot ->GetBinContent(hMCTot ->GetMaximumBin());
  double maxMCNC   = hMCNC  ->GetBinContent(hMCNC  ->GetMaximumBin());
  double maxMCNuE  = hMCNuE ->GetBinContent(hMCNuE ->GetMaximumBin());
  double maxMCNuMu = hMCNuMu->GetBinContent(hMCNuMu->GetMaximumBin());

  maxi = std::max(maxDat, maxMCTot);
  maxi = std::max(maxi, maxMCNC);
  maxi = std::max(maxi, maxMCNuE);
  maxi = std::max(maxi, maxMCNuMu);

  int nbins = hDat->GetNbinsX();
  for(int i = 1; i <= nbins; ++i) {
    hDat->SetBinError(i, sqrt(hDat->GetBinContent(i)));
  }

  // Scale down the ND by 10e3 (don't forget the max values!)
  double ndScale = 1./1000.;
  hDat   ->Scale(ndScale);
  hMCTot ->Scale(ndScale);
  hMCNC  ->Scale(ndScale);
  hMCNuE ->Scale(ndScale);
  hMCNuMu->Scale(ndScale);
  hNCStack  ->Scale(ndScale);
  hNueStack ->Scale(ndScale);
  hNumuStack->Scale(ndScale);
  maxi *= ndScale;

  SetHistOptions(hDat,    maxi, histTitle, 0, kBlack,                  false);
  SetHistOptions(hMCTot,  maxi, histTitle, 0, kTotalMCColor,           false);
  SetHistOptions(hMCNC,   maxi, histTitle, 0, kNCBackgroundColor,      false);
  SetHistOptions(hMCNuE,  maxi, histTitle, 0, kBeamNueBackgroundColor, false);
  SetHistOptions(hMCNuMu, maxi, histTitle, 0, kNumuBackgroundColor,    false);

  hDat->SetMarkerStyle(kFullCircle);

  SetHistOptions(hNueStack,  maxi, histTitle, 0, kBeamNueBackgroundColor, true);
  SetHistOptions(hNumuStack, maxi, histTitle, 0, kNumuBackgroundColor,    true);
  SetHistOptions(hNCStack,   maxi, histTitle, 0, kNCBackgroundColor,      true);

  // Create a legend for the histograms
  double xL = 0.5, xR = 0.7, yB = 0.64, yT = 0.84;
  TLegend* leg = new TLegend(xL, yB, xR, yT);
  SetLegendOptions(leg);
  leg->AddEntry(hDat,    "ND Data", "lep");
  leg->AddEntry(hMCTot,  "Total Prediction", "l");
  leg->AddEntry(hMCNC,   "NC Prediction", "l");
  leg->AddEntry(hMCNuE,  "#nu_{e} CC Background", "l");
  leg->AddEntry(hMCNuMu, "#nu_{#mu} CC Background", "l");
  leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);

  TH1* hRat  = (TH1*)hDat->Clone("hRat");
  hRat->Divide(hMCTot);
  TH1* hLine = (TH1*)hMCTot->Clone("");

  double maxvalRat = 1., minvalRat = 1.;
  for(int i = 1; i <= nbins; ++i) {
    if(hDat->GetBinContent(i) > 0) {
      double errpct = hDat->GetBinError(i)/hDat->GetBinContent(i);
      hRat->SetBinError(i, errpct);
    }

    double binval = hRat->GetBinContent(i);
    double errval = hRat->GetBinError(i);
    if(binval > 0) {
      maxvalRat = std::max(maxvalRat, binval + errval);
      minvalRat = std::min(minvalRat, binval - errval);
    }

    hLine->SetBinContent(i, 1.);
  }

  hRat->GetYaxis()->SetTitle("Data / Total MC");
  hRat->SetLineColor(kBlack);
  hLine->GetYaxis()->SetTitle("Data / Total MC");
  hLine->SetLineColor(kTotalMCColor);

  hRat ->SetMaximum(1.05*maxvalRat);
  hLine->SetMaximum(1.05*maxvalRat);
  hRat ->SetMinimum(0.95*minvalRat);
  hLine->SetMinimum(0.95*minvalRat);

  CenterTitles(hRat);

  if(plotlabel.find("NM1") == std::string::npos &&
     plotlabel.find("ECalLongAllE") == std::string::npos) {
    TCanvas* cnorat = new TCanvas(plotlabel.c_str(), plotlabel.c_str(), 800, 500);
    gPad->SetFillStyle(0);
    hDat->SetMinimum(0.001);
    hMCTot->SetMinimum(0.001);
    hMCNC->SetMinimum(0.001);
    hMCNuE->SetMinimum(0.001);
    hMCNuMu->SetMinimum(0.001);
    //gPad->SetLogy();
    hDat->Draw("E1");
    hMCTot->Draw("hist same");
    hMCNC->Draw("hist same");
    hMCNuE->Draw("hist same");
    hMCNuMu->Draw("hist same");
    leg->Draw();
    gPad->RedrawAxis();

    Preliminary();
    plotsOut->WriteTObject(cnorat);

    // Create a canvas for and draw the stacked histograms
    TCanvas* cStack = new TCanvas((plotlabel + "Stack").c_str(),
                                  (plotlabel + "Stack").c_str(), 800, 500);
    gPad->SetFillStyle(0);
    //gPad->SetLogy();
    TLegend* legStack = new TLegend(xL, yB, xR, yT);
    SetLegendOptions(legStack);
    legStack->AddEntry(hDat, "ND Data", "lep");
    legStack->AddEntry(hNCStack,   "NC 3 Flavor Prediction", "f");
    legStack->AddEntry(hNueStack,  "#nu_{e} CC Background", "f");
    legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
    legStack->SetY1(legStack->GetY2() - legStack->GetNRows()*0.05);

    hNCStack  ->Draw("hist");
    hNueStack ->Draw("hist same");
    hNumuStack->Draw("hist same");
    hDat->Draw("same");
    legStack->Draw();
    gPad->RedrawAxis();

    if(plotlabel.find("RecoELong") != std::string::npos) {
      FILE* textFO;
      textFO = fopen("out2.txt", "a+");
      //    Spectrum sData = hDat.FakeData(pot);
      fprintf(textFO, "Event numbers at the %2.2s for %s %.2fe20 POT:\n", "ND", plotlabel.c_str(), 8.04);
      if(hDat) { fprintf(textFO, "%12.12s, ", "Data"); }
      else     { fprintf(textFO, "%14.14s", ""); }
      fprintf(textFO, "%12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
              "All", "NC (3 Flav)", "Numu", "Nue", "FOM");
      
      double nAll  = hMCTot ->Integral();
      double nNC3F = hMCNC  ->Integral();
      double nNumu = hMCNuMu->Integral();
      double nNue  = hMCNuE ->Integral();
      double fom   = nNC3F/sqrt(nAll);
      if(hDat) {
        double nData = hDat->Integral(1, hDat->GetNbinsX());
           fprintf(textFO, "%10E, ", nData);
      }
      else {
      fprintf(textFO, "%14.14s", "");
      }
      fprintf(textFO, "%10E, %10E, %10E, %10E, %10E\n\n",
              nAll, nNC3F, nNumu, nNue, fom);
      fclose(textFO);
    }
    Preliminary();
    plotsOut->WriteTObject(cStack);
  }
  else {
    FILE* textF;
    textF = fopen("out.txt", "a+");
    Spectrum sData = data.FakeData(pot);
    std::cout << "POT: " << pot << std::endl;
    /// TODO: The pot values need to not be hard-coded here
    PlotStack(mcspecs, plotsOut, textF, plotlabel, "ND",
              8.04e20, 8.04, GetPlotOptions(plotlabel), &sData);
    fclose(textF);
  }

  TCanvas* canvas = new TCanvas((plotlabel + "Rat").c_str(),
                                (plotlabel + "Rat").c_str(), 800, 800);
  gPad->SetFillStyle(0);
  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
  TPad* pRatio = new TPad("pRatio", "", 0., 0.,    1., 0.375);
  pSpecs->Draw();
  pRatio->Draw();

  pSpecs->cd();
  hDat->Draw("E1");
  hDat->GetXaxis()->SetRangeUser(0,5);
  hDat->GetXaxis()->SetNdivisions(-405);
  hMCTot->Draw("hist same");
  hMCNC->Draw("hist same");
  hMCNuE->Draw("hist same");
  hMCNuMu->Draw("hist same");
  leg->Draw();
  gPad->RedrawAxis();

  Preliminary();
  pRatio->cd();
  hLine->Draw("hist ][");
  hRat->Draw("E1 same");

  plotsOut->WriteTObject(canvas);
  std::cout << "done\n";
};

//------------------------------------------------------------------------------
PlotOptions GetPlotOptions(std::string name)
{
  PlotOptions ret;

  int div = -404;
  double xL = 0.5, xR = 0.7, yB = 0.62, yT = 0.87;
  double texX = xL + .01;
  double texY = yB - 0.045;
  double cutL = 1.;
  double cutR = 0.;

  if(name.find("RecoELong")    != std::string::npos) {
    div = -406;
    xL = 0.35;
    xR = 0.55;
    texX = xL + .01;
    cutL = 0.5;
    cutR = 4.0;
  }
  if(name.find("VtxX")    != std::string::npos) {
    div = -404;
    xL = 0.4;
    xR = 0.6;
    texY = yT - 0.045;
    cutL = -100.;
    cutR =  100.;
  }
  if(name.find("VtxY")    != std::string::npos) {
    div = -404;
    xL = 0.4;
    xR = 0.6;
    texX = 0.2;
    texY = yT - 0.045;
    cutL = -100.;
    cutR =  100.;
  }
  if(name.find("VtxZ")    != std::string::npos) {
    div = -503;
    xL = 0.28;
    xR = 0.48;
    texX = xL + 0.01;
    cutL = 200.;
    cutR = 1000.;
  }
  if(name.find("CVN")     != std::string::npos) {
    div = -505;
    xL = 0.35;
    xR = 0.55;
    texX = xL + 0.01;
    cutL = 0.2;
  }
  if(name.find("NHit")    != std::string::npos) {
    div = -505;
    cutL = 20.;
    cutR = 200.;
  }
  if(name.find("PTP")     != std::string::npos) {
    div = -505;
    xL = 0.125;
    xR = 0.325;
    texX = 0.5;
    texY = yT - 0.045;
    cutR =  0.8;
  }
  if(name.find("EpH") != std::string::npos) {
    div = -505;
    cutL = 0.009;
  }
  if(name.find("DistTop") != std::string::npos) {
    div = -404;
    xL = 0.4;
    xR = 0.6;
    texX = 0.675;
    texY = yT - 0.045;
  }

  ret.ndiv = div;
  ret.xL = xL;
  ret.xR = xR;
  ret.yT = yT;
  ret.yB = yB;
  ret.texX = texX;
  ret.texY = texY;
  ret.cutL = cutL;
  ret.cutR = cutR;

  return ret;
}