plot_nue_data_mc.C

Go to the documentation of this file.

#ifdef __CINT__
void plot_nd_data_mc(std::string fname){}
#else

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Analysis/Style.h"
using namespace ana;

#include "OscLib/func/IOscCalculator.h"

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

void Legend(double x0, double y0, double x1, double y1)
{
  // Doesn't handle log-y well
  //  TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
  TLegend* leg = new TLegend(x0, y0, x1, y1);
  leg->SetFillStyle(0);

  TH1* dummy = new TH1F("", "", 1, 0, 1);
  dummy->SetMarkerStyle(kFullCircle);

  leg->AddEntry(dummy->Clone(), "ND data", "lep");
  dummy->SetLineColor(kGray+2);
  dummy->SetMarkerColor(kGray+2);
  //  leg->AddEntry(dummy->Clone(), "ND data stagger", "lep");
  dummy->SetLineColor(kTotalMCColor);
  leg->AddEntry(dummy->Clone(), "Total MC", "l");
  dummy->SetLineColor(kNCBackgroundColor);
  leg->AddEntry(dummy->Clone(), "NC", "l");
  dummy->SetLineColor(kBeamNueBackgroundColor);
  leg->AddEntry(dummy->Clone(), "Beam #nu_{e} CC", "l");
  dummy->SetLineColor(kNumuBackgroundColor);
  leg->AddEntry(dummy->Clone(), "#nu_{#mu} CC", "l");

  leg->Draw();
}

void plot_nd_data_mc(std::string fname)
{
  struct HistDef
  {
    std::string name;
  };

  // TODO definitely share. Merge macros?
  const int kNumVars = 52;
  const HistDef defs[kNumVars] = {
    "lem",
    "lem_zoom",
    "nhits",
    "lem_pidexp",
    "lem_meany",
    "lem_qfrac",
    "lem_ediff",
    "lem_enrich",
    "calE",
    "nvtx",
    "vtxx",
    "vtxy",
    "vtxz",

    "nshw",
    "shwstartx",
    "shwstarty",
    "shwstartz",
    "shwstopx",
    "shwstopy",
    "shwstopz",
    "longprong",
    "e_per_hit",
    "shw_nhit",
    "shw_len",
    "shw_width",
    "shw_gap",
    "shw_efrac",

    "shw_e",
    "had_e",

    "ptp",
    "costheta",
    "hits_per_plane",
    "pi0mass",

    "lid",
    "lid_zoom",

    "mulll",
    "elll",
    "emulll",
    "eglll",
    "eplll",
    "enlll",
    "epilll",
    "epi0lll",

    "mullt",
    "ellt",
    "emullt",
    "egllt",
    "epllt",
    "enllt",
    "epillt",
    "epi0llt",

    "confuse"
  };

  const int kNumSels = 4;
  const std::string selNames[kNumSels] = {"nocut", "presel", "lem", "lid"};

  Spectrum* spects[kNumSels][kNumVars];
  Spectrum* spectsStagger[kNumSels][kNumVars];
  IPrediction* preds[kNumSels][kNumVars];

  osc::NoOscillations noosc;

  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const char* selName = selNames[selIdx].c_str();
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const char* varName = defs[varIdx].name.c_str();

      spects[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/spect_%s", selName, varName).Data()).release();
      spectsStagger[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/spect_%s_stagger", selName, varName).Data()).release();
      preds[selIdx][varIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/pred_%s", selName, varName).Data()).release();

      if(selIdx == 0 && varIdx == 0)
        std::cout << "data, MC POT: "
                  << spects[selIdx][varIdx]->POT() << " "
                  << preds[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;


      if(defs[varIdx].name == "confuse" && selIdx == 1){
        const double pot = spects[selIdx][varIdx]->POT();
        TH1* hd = spects[selIdx][varIdx]->ToTH1(pot);
        TH1* hm = preds[selIdx][varIdx]->Predict(&noosc).ToTH1(pot);
        std::cout << "Neither, LID, LEM, both" << std::endl;
        std::cout << "DATA: " << hd->GetBinContent(1) << " " << hd->GetBinContent(2) << " " << hd->GetBinContent(3) << " " << hd->GetBinContent(4) << " " << std::endl;
        std::cout << "MC: " << hm->GetBinContent(1) << " " << hm->GetBinContent(2) << " " << hm->GetBinContent(3) << " " << hm->GetBinContent(4) << " " << std::endl;
      }
    } // end for varIdx
  } // end for selIdx

  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const std::string selName = selNames[selIdx];
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const std::string varName = defs[varIdx].name;

      new TCanvas;

      Spectrum* spect = spects[selIdx][varIdx];
      Spectrum* spectStagger = spectsStagger[selIdx][varIdx];
      IPrediction* pred = preds[selIdx][varIdx];

      TH1* h = DataMCComparisonComponents(*spect, pred, &noosc);
      h->SetTitle("NO#nuA Near Detector");
      h->GetXaxis()->CenterTitle();
      h->GetYaxis()->CenterTitle();

      TH1* hstag = spectStagger->ToTH1(spect->POT());
      hstag->SetMarkerStyle(kFullCircle);
      hstag->SetMarkerColor(kGray+2);
      hstag->SetLineColor(kGray+2);
      //      hstag->Draw("ep same");

      TH1* horig = spect->ToTH1(spect->POT());
      horig->SetMarkerStyle(kFullCircle);
      horig->Draw("ep same");


      Legend(.6, .6, .85, .85);

      gPad->Print(("plots/"+selName+"_"+varName+"_lin.eps").c_str());


      h = DataMCComparisonComponents(*spect, pred, &noosc);
      h->SetTitle("NO#nuA Near Detector");
      h->GetXaxis()->CenterTitle();
      h->GetYaxis()->CenterTitle();

      hstag->Draw("ep same");
      horig->Draw("ep same");

      gPad->SetLogy();

      Legend(.6, .6, .85, .85);

      gPad->Print(("plots/"+selName+"_"+varName+"_log.eps").c_str());


      new TCanvas;

      DataMCRatio(*spect, pred, &noosc);

      gPad->Print(("plots/"+selName+"_"+varName+"_ratio.eps").c_str());


      TCanvas* c = new TCanvas;

      TPad* p1 = new TPad("", "", 0, 0, 1, 1);
      TPad* p2 = new TPad("", "", 0, 0, 1, 1);
      p1->SetBottomMargin(.3);
      p2->SetTopMargin(.7);
      p1->SetFillStyle(0);
      p2->SetFillStyle(0);
      p1->Draw();
      p2->Draw();

      p1->cd();
      h = DataMCComparisonComponents(*spect, pred, &noosc);
      h->SetTitle("NO#nuA Near Detector");
      h->GetXaxis()->CenterTitle();
      h->GetYaxis()->CenterTitle();
      hstag->Draw("ep same");
      horig->Draw("ep same");
      Legend(.6, .6, .85, .85);

      p2->cd();
      DataMCRatio(*spect, pred, &noosc, 0, 2.99);
      Ratio ratioStagger(*spectStagger, pred->PredictUnoscillated());
      TH1* hratio = ratioStagger.ToTH1();
      hratio->SetLineColor(kGray+2);
      //      hratio->Draw("e same");

      c->cd();

      gPad->Print(("plots/"+selName+"_"+varName+"_lin_ratio.eps").c_str());

      p1->cd();
      h = DataMCComparisonComponents(*spect, pred, &noosc);
      h->SetTitle("NO#nuA Near Detector");
      h->GetXaxis()->CenterTitle();
      h->GetYaxis()->CenterTitle();
      //      hstag->Draw("ep same");
      horig->Draw("ep same");
      gPad->SetLogy();
      Legend(.6, .6, .85, .85);
      c->cd();
      gPad->Print(("plots/"+selName+"_"+varName+"_log_ratio.eps").c_str());
    } // end for varIdx
  } // end for selIdx
}

#endif