plot_DataMCComp_numu.C

Go to the documentation of this file.

//------------------------------------------------------------------------------
// This plots the noextrap predictions made by make_DataMCComp_numu.C.
// The map of variables are defined at the bottom of the header.
// FHC pot-equiv is put in by hand at line 145
//
// The fakeData bool should be false if you are using real data. I just had that
//  there for testing while the data files were still blinded
//
// Andrew Sutton - ats2mk@virginia.edu
//------------------------------------------------------------------------------

#pragma once

#include "DataMCComp_header.h"

using namespace ana;

TH1D* Load1DHistFromPred(TFile* inFile, std::string name, PredSel predsel, double pot, osc::IOscCalc *calc);

TH1D* LoadFakeDataFromPred(TFile* inFile, std::string name, double pot);

TH1D* Load1DHistFromSpec(TFile* inFile, std::string name, bool cosmic, double  pot, double lt);

void plot_DataMCComp_numu(std::string inName, bool fakeData = false)
{
  TFile *inFile = new TFile(inName.c_str(), "READ");
  //TFile *dataFile = new TFile("dataonly_fhc_period1235.root", "READ");
  //  TFile *dataFile = new TFile("dataonly_fhc_period910.root", "READ");

  Sign::Sign_t WS;
  std::string horn, HORN;
  double pot, lt;
  
  if (inName.find("fhc") != std::string::npos) {
    WS   = Sign::kAntiNu;
    horn = "Neutrino Beam";
    HORN = "FHC";
    pot  = kAna2020FHCPOT;
    lt   = kAna2020FHCLivetime;
  }
  else if (inName.find("rhc") != std::string::npos) {
    WS   = Sign::kNu;
    HORN = "RHC";
    horn = "Antineutrino Beam";
    pot  = kAna2020RHCPOT; 
    lt   = kAna2020RHCLivetime; 
  }
  else std::abort();

  // get the BF
  std::string sHieOct  = "NH_UO";
  std::string sFitDir  = "/nova/ana/3flavor/Ana2020/Fits/BestFits/"; 
  std::string sFitFile = sFitDir  + "bestfits_joint_realData_both_systs.root";
  TFile *fBestFit = new TFile(sFitFile.c_str(), "READ");
  auto calc = (osc::IOscCalcAdjustable*)(LoadFrom<osc::IOscCalc>(fBestFit, (sHieOct + "_calc").c_str())).release();

  // Get POT and LT from the datafile
  double dpot = pot;
  double dlt  = lt;
  if (inFile->GetDirectory("dtSpec_muonE_AllQuants")) {
    std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom(inFile, "dtSpec_muonE_AllQuants");
    dpot = TempSpec->POT();
    dlt  = TempSpec->Livetime();
  }
  
  // Quantiles and desired
  const std::vector<std::string> quantNames =
    { "AllQuants", "Quant1", "Quant2", "Quant3", "Quant4"};

  std::string WSH = "WS" + HORN;
  
  std::vector< TH1D* > hPred[5];
  std::vector< TH1D* > hBkg[5];
  std::vector< TH1D* > hNumu[5];
  std::vector< TH1D* > hWS[5];
  std::vector< TH1D* > hCosm[5];
  std::vector< TH1D* > hData[5];
  std::vector< TH1D* > hDataUS[5];
  
  for (uint q = 0; q < quantNames.size(); ++q) {
    for (auto var : vars) {

      std::string varName = var.first.first;
      
      hPred[q].push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at("All") , pot, calc) );
      hBkg[q] .push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at("All") , pot, calc) );
      hNumu[q].push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at("Numu"), pot, calc) );
      hWS[q]  .push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at(WSH)   , pot, calc) );

      hCosm[q].push_back( Load1DHistFromSpec(inFile, "csSpec_" + varName + "_" + quantNames[q], true , pot, lt) );

      if (fakeData) {
        hData  [q].push_back( LoadFakeDataFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], pot) );
        hDataUS[q].push_back( LoadFakeDataFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], pot) );
      }
      else {
        hData  [q].push_back( Load1DHistFromSpec(inFile, "dtSpec_" + varName + "_" + quantNames[q], false, dpot, dlt) );
        hDataUS[q].push_back( Load1DHistFromSpec(inFile, "dtSpec_" + varName + "_" + quantNames[q], false, dpot, dlt) );
      }
      
      // adding together for custom hist stack
      hBkg [q].back()->Add(hNumu[q].back(), -1); // remove the signal
      hBkg [q].back()->Add(hCosm[q].back(), +1); // add on the cosmics for stack
      hWS  [q].back()->Add(hBkg [q].back(), +1); // add bkg to WS for stack
      hPred[q].back()->Add(hCosm[q].back(), +1); // add cosmics to full prediction
      
      // Beautification
      PimpHist(hPred [q].back(), kSolid, kFullPredColor        , 3, kFullCircle, kFullPredColor        , 1);
      PimpHist(hWS   [q].back(), kSolid, kNumuWSColor          , 1, kFullCircle, kNumuWSColor          , 1);
      PimpHist(hBkg  [q].back(), kSolid, kTotalBkgColor        , 1, kFullCircle, kTotalBkgColor        , 1);
      PimpHist(hCosm [q].back(), kSolid, kCosmicBackgroundColor, 1, kFullCircle, kCosmicBackgroundColor, 1);

      FillWithDimColor(hWS[q].back()  , 0);
      FillWithDimColor(hCosm[q].back(), 0);
      FillWithDimColor(hBkg[q].back() , 0);

      if (varName == "recoE"){
        hData[q].back()->Scale(0.1, "width");
        hPred[q].back()->Scale(0.1, "width");
        hCosm[q].back()->Scale(0.1, "width");
        hBkg [q].back()->Scale(0.1, "width");
        hWS  [q].back()->Scale(0.1, "width");
        hNumu[q].back()->Scale(0.1, "width");
      }
    }// end var loop
  }// end quant loop

  // Plotting
  int v = 0;
  for (auto var: vars) {
    std::ofstream f(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.txt").c_str());
    std::ofstream f_ByQ(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.txt").c_str());
    TCanvas *c = new TCanvas("c", "c", 960, 800);
    TCanvas *c_ByQ;
    TPad *p1, *p2, *p3, *p4;
    SplitCanvasQuant(c_ByQ, p1, p2, p3, p4);
    c_ByQ->SetTitle("c_ByQ");

    TLatex* xTitle = new TLatex(0.5, 0.03, var.first.second.c_str());
    xTitle->SetNDC();
    xTitle->SetTextAlign(22);
    TString eventAppend = Form("%.2f #times 10^{20} POT", dpot/1e20);
    if (HORN == "FHC") {
      // put in by hand for equiv
      dpot = 13.6e20;
      eventAppend = Form("%.2f #times 10^{20} POT-equiv.", dpot/1e20);
    }
    if (var.first.first == "recoE")
      eventAppend = "0.1 GeV";
    TLatex* yTitle = new TLatex(0.023, 0.5, "Events / " + eventAppend);
    yTitle->SetNDC();
    yTitle->SetTextAlign(22);
    yTitle->SetTextAngle(90.);
 
    // setup the legends
    TLegend *leg   = new TLegend(0.60,0.60,0.88,0.88);
    TLegend *leg_l = new TLegend(0.10,0.75,0.25,0.88);
    TLegend *leg_r = new TLegend(0.52,0.75,0.65,0.88);
    leg->AddEntry(hData[0][v], "FD data"       ,"LEP"); leg_l->AddEntry(hData[0][v], "FD data"       ,"LEP");
    leg->AddEntry(hPred[0][v], "2020 Best-fit" ,"L"  ); leg_l->AddEntry(hPred[0][v], "2020 Best-fit" ,"L"  );
    leg->AddEntry(hBkg [0][v], "Background"    ,"BF" ); leg_l->AddEntry(hBkg [0][v], "Background"    ,"BF" );
    leg  ->SetTextSize  (0.05); leg  ->SetBorderSize(0); leg  ->SetFillStyle (0);
    leg_l->SetTextSize  (0.03); leg_l->SetBorderSize(0); leg_l->SetFillStyle (0);
    leg_r->SetTextSize  (0.03); leg_r->SetBorderSize(0); leg_r->SetFillStyle (0);

    // Get the maximum quant bin for plotting
    float maxDataByQ = 0;
    float maxPredByQ = 0;
    for (uint q = 1; q < quantNames.size(); ++q) {
      float maxData = hData[q][v]->GetMaximum();
      float maxPred = hPred[q][v]->GetMaximum();

      if (maxData > maxDataByQ) maxDataByQ = maxData;
      if (maxPred > maxPredByQ) maxPredByQ = maxPred;
    }

    // quantile loop
    for (uint q = 0; q < quantNames.size(); ++q) {
      // Making plots look better and changing pads
      if (q == 0)
        MakeHistCanvasReady_Quant(q, hPred[q][v], hData[q][v]->GetMaximum()*2.2);
      else
        MakeHistCanvasReady_Quant(q, hPred[q][v], maxDataByQ*1.8);
      
      if (q == 0) c ->cd();
      if (q == 1) p1->cd();
      if (q == 2) p2->cd();
      if (q == 3) p3->cd();
      if (q == 4) p4->cd();
      
      // Set the data error bars and draw everything
      TGraph * grData = GraphWithPoissonErrors2(hData[q][v], hPred[q][v], false, true);
      if (var.first.first == "recoE") grData = graphAsymmErrorScaled(hData[q][v], hDataUS[q][v]);
      grData->SetMarkerStyle(kFullCircle);
      
      hPred[q][v]->Draw("hist same");
      hBkg[q][v] ->Draw("hist same");
      if ( !(var.first.first == "trueE" || var.first.first == "resE") )
        grData->Draw("P");
      
      DrawQuantLabel(q);
    }// end quant loop

    c->cd();
    leg->Draw();
    xTitle->Draw();
    yTitle->Draw();
    Preliminary();
    DrawBeamLabel(HORN=="FHC");
    c->Update();
    c->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.pdf").c_str());
    c->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.png").c_str());
    c->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.eps").c_str());
    f << HORN << " " << var.first.second
      << " with an un-extrapolated prediction set to the 2020 joint analysis best fit point"
      << std::endl;

    c_ByQ->cd();
    Preliminary();
    DrawBeamLabel(HORN=="FHC");
    xTitle->SetTextSize(0.04);
    yTitle->SetTextSize(0.04);
    xTitle->Draw();
    yTitle->Draw();
    p1->cd();
    leg_l->Draw();
    p2->cd();
    leg_r->Draw();
    c_ByQ->Update();
    c_ByQ->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.pdf").c_str());
    c_ByQ->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.png").c_str());
    c_ByQ->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.eps").c_str());
    f_ByQ << HORN << " " << var.first.second
          << " split by hadronic energy fraction quartile with an"
          << " un-extrapolated prediction set to the 2020 joint analysis best fit point"
          << std::endl;

    f.close();
    f_ByQ.close();
      
    delete c;
    delete c_ByQ;

    ++v;
  }// end var loop
}


//*******************************************************************************
TH1D* Load1DHistFromPred(TFile* inFile, std::string name, PredSel predsel, double pot, osc::IOscCalc* calc)
{
  // -- Here you can change the best fit point
  // osc::OscCalcDumb* calc = new osc::OscCalcDumb();
  // osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  // calc->SetdCP (2*M_PI);
  // calc->SetTh23(asin(sqrt(0.565)));
  // calc->SetDmsq32(2.48e-3);

  std::unique_ptr<PredictionNoExtrap> TempPred = 
    PredictionNoExtrap::LoadFrom(inFile, (name).c_str());

  return TempPred->PredictComponent(calc, predsel.flav, predsel.curr, predsel.sign).ToTH1(pot);
}

//*******************************************************************************
TH1D* LoadFakeDataFromPred(TFile* inFile, std::string name, double pot)
{
  osc::OscCalcDumb* calc = new osc::OscCalcDumb();

  std::unique_ptr<PredictionNoExtrap> TempPred = 
    PredictionNoExtrap::LoadFrom(inFile, (name).c_str());

  return TempPred->Predict(calc).MockData(pot).ToTH1(pot);

}

//*******************************************************************************
TH1D* Load1DHistFromSpec(TFile* inFile, std::string name, bool cosmic, double pot, double lt)
{
  std::unique_ptr<Spectrum> TempSpec = 
    Spectrum::LoadFrom(inFile, (name).c_str());
 
  if(cosmic)
    return TempSpec->ToTH1(lt, kLivetime);
  else
    return TempSpec->ToTH1(pot);
}