drawSystsShiftingNDdata.C

Go to the documentation of this file.

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "OscLib/OscCalcPMNSOpt.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"

#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Extrap/ModularExtrapComponent.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"

#include "Utilities/func/ROOTHelpers.h"

#include "TStyle.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TArrow.h"

using namespace ana;

const unsigned int kNumQuantiles = 4;
const std::vector<std::string> kOutputExts = {".png", ".eps", ".pdf", ".root"};


const std::map<const ISyst*, std::string> kSystRename
{
  { &kMAQEGenieReducedSyst2018,         "M_{A}^{QE}" },
  { &kRPACCQEEnhSyst2018,               "RPA CCQE (enh)" },
  { &kRPACCQESuppSyst2018,              "RPA CCQE (supp)" },
  { &kRPARESSyst2018,                   "RPA RES" },
  { &kMECShapeSyst2018Nu,               "MEC 4-mom. xfer shape" },
  // yeah -- use the same text as for nu.
  // these plots are illustrative only, and we won't make the WS versions.
  { &kMECShapeSyst2018AntiNu,           "MEC 4-mom. xfer shape" },
  { &kMECEnuShapeSyst2018Nu,            "MEC E_{#nu} shape" },
  { &kMECEnuShapeSyst2018AntiNu,        "MEC E_{#nu} shape" },  // etc.
  { &kMECInitStateNPFracSyst2018Nu,     "MEC it'l. state np fraction" },
  { &kMECInitStateNPFracSyst2018AntiNu, "MEC it'l. state np fraction" },

  { GetGenieKnobSyst(rwgt::fReweightMaCCRES),  "M_{A}^{RES}" },
  { GetGenieKnobSyst(rwgt::fReweightMvCCRES),  "M_{V}^{RES}" },
};

// Put a "NOvA Simulation" tag in the corner      .
// this was clashing with the .rootlogon...
namespace numu {
    void Simulation() {
        TLatex *prelim = new TLatex(.9, .95, "NOvA Simulation");
        prelim->SetTextColor(kGray + 1);
        prelim->SetNDC();
        prelim->SetTextSize(2 / 30.);
        prelim->SetTextAlign(32);
        prelim->Draw();
    }
}

void WriteBeam(Loaders::FluxType beam)
{
    auto beamLabel = new TLatex(.15, .92, (beam == Loaders::kFHC ? "Neutrino" : "Antineutrino") + TString(" beam"));
    beamLabel->SetNDC();
    beamLabel->SetTextSize(0.04);
    beamLabel->Draw();
}

void WriteQuartile(unsigned int q)
{
  if (q >= kNumQuantiles)
     return;

  TLatex* quart = new TLatex(.17, .75, Form("Quartile %i",q+1));
  quart->SetNDC();
  quart->SetTextSize(0.05);
  quart->Draw();
}

void WriteSystName(const ISyst * syst, bool top=false)
{
    std::string name = (kSystRename.find(syst) != kSystRename.end()) ? kSystRename.at(syst) : syst->LatexName();
    TLatex * text = new TLatex(0.17, top ? 0.83 : 0.2, name.c_str());
    text->SetNDC();
    text->SetTextSize(1.5/30.);
    text->Draw();
}

TH1* getResidualDiff(TH1* h1, TH1* h2)
{
  TH1* h = (TH1*) h1->Clone("h");
  h->Reset();

  for(int k = 1; k <= h1->GetNbinsX(); k++)
  {
    double val1 = (double) h1->GetBinContent(k);
    double val2 = (double) h2->GetBinContent(k);

    if(val1 * val2 > 0)
    {
      h->SetBinContent(k, (val1-val2) * 200./(val1 + val2));
    }
  }
  return h;
}

TDirectory * GetDir(TFile * f, const TString & dirName)
{
//    std::cout << "Loading directory: " << dirName << std::endl;
    return f->GetDirectory(dirName);
}

TH1 * GetQuantilePredictionHist(unsigned int quantNum, IPrediction * preds[], osc::IOscCalc * calc, double POT)
{
    if (quantNum < kNumQuantiles)
        return preds[quantNum]->Predict(calc).ToTH1(POT);

    // rather than taking the prediction that's extrapolated using all 4 quartiles to get the to get the sum,
    // sum up the predictions from each individually extrapolated quartile to get the total.
    TH1 * h = preds[0]->Predict(calc).ToTH1(POT);
    for (unsigned int i = 1; i < kNumQuantiles; i++)
        h->Add(preds[i]->Predict(calc).ToTH1(POT));
    return h;
}

void drawSystsShiftingNDdata(std::string hornCurrStr, std::string predFile, std::string outDir)
{
    Loaders::FluxType hornCurr;
    if (hornCurrStr == "FHC")
        hornCurr = Loaders::kFHC;
    else if (hornCurrStr == "RHC")
        hornCurr = Loaders::kRHC;
    else
    {
        std::cerr << "Unrecognized horn current string: '" << hornCurrStr << "'.  Options are FHC, RHC." << std::endl;
        return;
    }
    std::transform(hornCurrStr.begin(), hornCurrStr.end(), hornCurrStr.begin(), ::tolower);

  std::cout << "Loading predictions...";
  std::cout << std::flush;

  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
  gStyle->SetPaintTextFormat("6.3g");
  gStyle->SetFitFormat("6.3g");

  double POT = hornCurr == Loaders::kFHC ? kAna2018FHCPOT : kAna2018RHCPOT;

  osc::OscCalcPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc);

  calc.SetTh23(TMath::ASin(sqrt(0.58))); // 2018 joint fit result
  calc.SetDmsq32(2.51e-3); // Same

  std::vector<const ISyst*> systs = getAllDirectNumuSysts2018();

  const unsigned int kNumSysts = systs.size();

  IPrediction* predictionND[kNumSysts*2][kNumQuantiles];
  IPrediction* predictionFD[kNumSysts*2][kNumQuantiles];
  IPrediction* nominalPred[kNumQuantiles];

  TFile* fPred = new TFile( predFile.c_str() );

  std::string predNameND = "ccE_shiftNDdata";
  std::string predNameFD = "ccE_noExtrapFD";

  for(unsigned int i = 0; i < kNumQuantiles; i++)
  {
    nominalPred[i] = LoadFrom<PredictionNoExtrap>(GetDir(fPred, Form("nominalPred_%i",i))).release();

    for(unsigned int k = 0; k < kNumSysts; k++)
    {
      predictionND[2*k][i] = LoadFrom<PredictionInterp>(GetDir(fPred, Form( (predNameND+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();

      predictionND[2*k+1][i] = LoadFrom<PredictionInterp>(GetDir(fPred,  Form( (predNameND+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();

      // Now FD (no extrap)

      predictionFD[2*k][i] = LoadFrom<PredictionNoExtrap>(GetDir(fPred,  Form( (predNameFD+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();

      predictionFD[2*k+1][i] = LoadFrom<PredictionNoExtrap>(GetDir(fPred,  Form( (predNameFD+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();
    } // end loop ver systs
  } // end loop over quantiles

  std::cout << "  done." << std::endl;

  std::cout << "Loading went OK. Plotting!" << std::endl;

  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
  {
    TH1* hNominal = GetQuantilePredictionHist(i, nominalPred, &calc, POT);
    hNominal->Scale(0.1, "width");

    hNominal->SetLineColor(kBlack);
    double ymax = hNominal->GetBinContent(hNominal->GetMaximumBin()) * 2.0;
    hNominal->GetYaxis()->SetRangeUser(0, 15);
    hNominal->GetYaxis()->SetTitle("Events / (0.1 GeV)");
    hNominal->GetYaxis()->CenterTitle();
    hNominal->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
    hNominal->GetXaxis()->CenterTitle();

    for (unsigned int k = 0; k < kNumSysts; k++)
    {
      std::string fullName = hornCurrStr+"_ccE_"+systs[k]->ShortName();

      TCanvas* c = new TCanvas("c",(fullName).c_str());

      c->SetBottomMargin(0.15);
      c->SetLeftMargin(0.15);

      TH1* hShiftND_up = GetQuantilePredictionHist(i, predictionND[2*k], &calc, POT);
      hShiftND_up->Scale(0.1, "width");
      TH1* hShiftFD_up = GetQuantilePredictionHist(i, predictionFD[2*k], &calc, POT);
      hShiftFD_up->Scale(0.1, "width");

      TH1* hShiftND_dn = GetQuantilePredictionHist(i, predictionND[2*k+1], &calc, POT);
      hShiftND_dn->Scale(0.1, "width");
      TH1* hShiftFD_dn = GetQuantilePredictionHist(i, predictionFD[2*k+1], &calc, POT);
      hShiftFD_dn->Scale(0.1, "width");

      hShiftND_up->SetLineColor(kOrange+1);
      hShiftND_up->SetLineStyle(2);

      hShiftND_dn->SetLineColor(kViolet+1);
      hShiftND_dn->SetLineStyle(2);

      hShiftFD_up->SetMarkerStyle(20);
      hShiftFD_up->SetMarkerColor(kOrange+1);
      hShiftFD_up->SetLineWidth(0);
      hShiftFD_up->SetLineColor(0);
      hShiftFD_dn->SetMarkerStyle(21);
      hShiftFD_dn->SetMarkerColor(kViolet+1);
      hShiftFD_dn->SetLineWidth(0);
      hShiftFD_dn->SetLineColor(0);

      hNominal->Draw("hist");
      hShiftND_up->Draw("hist same");
      hShiftND_dn->Draw("hist same");
      hShiftFD_up->Draw("EX0 same");
      hShiftFD_dn->Draw("EX0 same");

      auto leg = std::make_unique<TLegend>(0.55,0.55,0.95,0.8);
      leg->AddEntry(hNominal, "Nominal with FD/ND","l");
      leg->AddEntry(hShiftND_up, "+ 1#sigma in ND","l");
      leg->AddEntry(hShiftFD_up, "+ 1#sigma in FD","p");
      leg->AddEntry(hShiftND_dn, "- 1#sigma in ND","l");
      leg->AddEntry(hShiftFD_dn, "- 1#sigma in FD","p");
      leg->SetFillStyle(0);
      leg->Draw();

      WriteSystName(systs[k], true);
      numu::Simulation();
      WriteQuartile(i);
      WriteBeam(hornCurr);

      auto qstr = (i < kNumQuantiles) ? std::to_string(i) : "all";
      util::SaveObj(c, fullName + "_Q" + qstr, outDir, kOutputExts, true);

      // Now the ratios

      c->Clear();

      hShiftND_up->Divide(hNominal);
      hShiftND_dn->Divide(hNominal);

      hShiftFD_up->Divide(hNominal);
      hShiftFD_dn->Divide(hNominal);

      hShiftND_up->GetYaxis()->SetRangeUser(0.8,1.2);
      //if(k == 0 && (kNumSysts <= 4)) hShiftND_up->GetYaxis()->SetRangeUser(0.5, 2.0); 
      hShiftND_up->GetYaxis()->SetTitle("Ratio to nominal MC");
      hShiftND_up->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
      hShiftND_up->GetYaxis()->CenterTitle();
      hShiftND_up->GetXaxis()->CenterTitle();

      hShiftND_up->SetMarkerColor(kOrange+1);
      hShiftND_up->SetLineWidth(2);
      hShiftND_up->SetLineStyle(1);
      hShiftND_up->SetLineColor(kOrange+1);
      hShiftND_dn->SetMarkerColor(kOrange+1);
      hShiftND_dn->SetLineWidth(2);
      hShiftND_dn->SetLineStyle(1);
      hShiftND_dn->SetLineColor(kOrange+1);

      hShiftFD_up->SetMarkerColor(kViolet+1);
      hShiftFD_up->SetLineWidth(2);
      hShiftFD_up->SetLineStyle(2);
      hShiftFD_up->SetLineColor(kViolet+1);
      hShiftFD_dn->SetMarkerColor(kViolet+1);
      hShiftFD_dn->SetLineWidth(2);
      hShiftFD_dn->SetLineStyle(2);
      hShiftFD_dn->SetLineColor(kViolet+1);

      hShiftND_up->Draw("hist");
      hShiftND_dn->Draw("hist same");
      hShiftFD_up->Draw("hist same");
      hShiftFD_dn->Draw("hist same");

      TGraph* grShade = ShadeBetweenHistograms(hShiftND_up, hShiftFD_up);
      grShade->SetFillColorAlpha(kRed-9, 0.35);
      grShade->SetLineColor(kRed-9);
      grShade->Draw("F");

      TGraph* grShade2 = ShadeBetweenHistograms(hShiftND_dn, hShiftFD_dn);
      grShade2->SetFillColorAlpha(kRed-9, 0.35);
      grShade2->SetLineColor(kRed-9);
      grShade2->Draw("F");

      TGraph* g1 = new TGraph;
      g1->SetPoint(0, -1, 1);
      g1->SetPoint(1, 10, 1);

      g1->SetLineWidth(1);
      g1->SetLineColor(kBlack);
      g1->SetLineStyle(2);

      g1->Draw("L");

      leg = std::make_unique<TLegend>(0.4,0.7,0.9,0.88);
      leg->SetFillStyle(0);
      leg->AddEntry(hShiftND_up, TString("#pm 1#sigma extrap. ND shift"),"l");
      leg->AddEntry(hShiftFD_up, TString("#pm 1#sigma FD shift"),"l");
      leg->SetFillColorAlpha(10,0.);
      leg->Draw();

      WriteSystName(systs[k]);
      numu::Simulation();
      WriteQuartile(i);
      WriteBeam(hornCurr);

      util::SaveObj(c, fullName + "_ratio_Q" + qstr, outDir, kOutputExts, true);

      // And finally the residual difference

      c->Clear();

      TH1* hResup = getResidualDiff(hShiftND_up, hShiftFD_up);

      TH1* hResdn = getResidualDiff(hShiftND_dn, hShiftFD_dn);

      hResup->GetYaxis()->SetTitle("Residual difference (%)");
      hResup->GetYaxis()->SetRangeUser( -10, 10 );
      //if(k == 0 && (kNumSysts <= 4)) hResup->GetYaxis()->SetRangeUser(-50, 50); 

      hResup->SetLineColor(kOrange+1);
      hResdn->SetLineColor(kViolet+1);

      hResup->Draw("hist");
      hResdn->Draw("hist same");

      leg->Clear();
      leg->AddEntry(hResup, TString("+1 #sigma shift in FD minus ND"),"l");
      leg->AddEntry(hResdn, TString("-1 #sigma shift in FD minus ND"),"l");
      leg->SetFillColorAlpha(10,0.);
      leg->Draw();

      TGraph* g0 = new TGraph;
      g0->SetPoint(0, -1, 0);
      g0->SetPoint(1, 10, 0);

      g0->SetLineWidth(1);
      g0->SetLineColor(kBlack);
      g0->SetLineStyle(2);

      g0->Draw("L");

      WriteSystName(systs[k]);
      numu::Simulation();
      WriteQuartile(i);
      WriteBeam(hornCurr);

      util::SaveObj(c, fullName + "_residual_Q" + qstr, outDir, kOutputExts, true);

      c->Close();
    } // end loop over systematics

  } // end loop over quantiles
  
  std::cout << "All plots finished!" << std::endl;

  fPred->Close();

} // End of function