SystsOscParamAna.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "NuXAna/Systs/NusSysts.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "NuXAna/macros/NuSPlotFunctions.h"
#include "NuXAna/macros/NuSSystFunctions.h"

#include "OscLib/OscCalcSterile.h"

#include "TArrow.h"
#include "TCanvas.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TH1.h"
#include "TKey.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TPad.h"

using namespace ana;

TH1* GetNC(IPrediction* specs, double POT, osc::OscCalcSterile* calc);
TH1* GetBG(IPrediction* specs, double POT, osc::OscCalcSterile* calc);

std::string MakeName(std::string chan, std::string det,
                     std::string syst, std::string sigma);

std::string MakeTitle(std::string chan, std::string det,
                      std::string syst, std::string sigma);

void PlotSyst(TH1* nom, TH1* hp1, TH1* hm1,
              TDirectory* out, FILE* text, strings strs, bool use_max);
void PlotSyst(TH1* h, TH1* hp1, TH1* hm1,
              TDirectory* out, FILE* text, strings strs);

// After running a fit, reset the values in a calculator to preset defaults
void ResetAngles(osc::OscCalcSterile* calc);

void SystsOscParamAna()
{
  TH1::AddDirectory(0);

  // A map of systematic labels
  std::map<std::string, std::string> shiftlabels;
  shiftlabels["theta13"]   = "#theta_{13}";
  shiftlabels["delta13"]   = "#delta_{13}";
  shiftlabels["Deltam232"] = "#Deltam^2_{32}";

  std::map<std::string, PredictionCombinePeriods*> preds;

  // Set up path to file with filled CAFAna objects and for output of analysis
  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
  std::string filenm = "SystsOscParam";

  std::string loadLocation = "/nova/ana/steriles/Ana01/FitSystEffects.root";
  std::string saveLocation = filenm + "Ana.root";
  std::string textLocation = filenm + "Ana.txt";

  // Load filled objects from file
  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
  TDirectory* tmpL = gDirectory;
  TDirectory* loadDir = gDirectory;
  loadDir->cd((loadLocation + ":/pred").c_str());
  PredictionCombinePeriods pred = *PredictionCombinePeriods::LoadFrom(gDirectory);
  rootL->Close();

  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc4fnom = DefaultSterileCalc(4);
  calc4fnom->SetNFlavors(4);
  ResetAngles(calc4fnom);
  
  osc::OscCalcSterile* calc4fplus = DefaultSterileCalc(4);
  calc4fplus->SetNFlavors(4);
  ResetAngles(calc4fplus);
  
  osc::OscCalcSterile* calc4fminus = DefaultSterileCalc(4);
  calc4fminus->SetNFlavors(4);
  ResetAngles(calc4fminus);


  
  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
  FILE* textF;
  textF = fopen(textLocation.c_str(), "w");

  strings strs;
  strs.fComponent = "NC";
  strs.fDet = "EX";
  strs.fPOT = StringFromDouble(NCSCALE);
  strs.fSample = "Ana01";
  strs.fSystType = "OscParam";
  strs.fXLabel = "Calorimetric Energy (GeV)";
  
  InitializeSystText(textF, strs);

  std::map<std::string, std::pair<double, double> > parammap;
  parammap["theta13"] = {8.475, 0.0353};
  parammap["delta13"] = {1.39, 0.325};
  parammap["Deltam232"] = {0.00244, 0.00006};

  TH1* hNomNC_predSt = GetNC(&pred, NCSCALE, calc4fnom);
  TH1* hNomBG_predSt = GetBG(&pred, NCSCALE, calc4fnom);

  for(const auto& shifts : shiftlabels) {
    ResetAngles(calc4fplus);
  
    strs.fSystS = shifts.first;
    strs.fSystL = shifts.second;

    if(shifts.first.compare("theta13") == 0) {
      calc4fplus->SetAngle(1, 3, calc4fplus->GetAngle(1, 3) + parammap[shifts.first].second);
    }
    else if(shifts.first.compare("delta13") == 0) {
      calc4fplus->SetDelta(1, 3, calc4fplus->GetDelta(1,3) + parammap[shifts.first].second);
    }
    else {
      calc4fplus->SetDm(3, calc4fplus->GetDm(3) + parammap[shifts.first].second);
    }
    
    TH1* hSp1NC_predSt = GetNC(&pred, NCSCALE, calc4fplus);
    TH1* hSp1BG_predSt = GetBG(&pred, NCSCALE, calc4fplus);
    
    ResetAngles(calc4fminus);

    if(shifts.first.compare("theta13") == 0) {
      calc4fminus->SetAngle(1, 3, calc4fminus->GetAngle(1, 3) - parammap[shifts.first].second);
    }
    else if(shifts.first.compare("delta13") == 0) {
      calc4fminus->SetDelta(1, 3, calc4fminus->GetDelta(1,3) - parammap[shifts.first].second);
    }
    else {
      calc4fminus->SetDm(3, calc4fminus->GetDm(3) - parammap[shifts.first].second);
    }
    
    TH1* hSm1NC_predSt = GetNC(&pred, NCSCALE, calc4fminus);
    TH1* hSm1BG_predSt = GetBG(&pred, NCSCALE, calc4fminus);

    strs.fComponent = "NC";
    strs.fDet = "EX";
    std::cout << "Plotting signal" << std::endl;
    PlotSyst(hNomNC_predSt, hSp1NC_predSt, hSm1NC_predSt,
             rootF, textF, strs, false);

    strs.fComponent = "BG";
    PlotSyst(hNomBG_predSt, hSp1BG_predSt, hSm1BG_predSt,
             rootF, textF, strs, false);
  }

  strs.fComponent = "NC";
  PlotSystBand(hNomNC_predSt, rootF, textF, strs);

  strs.fComponent = "BG";
  PlotSystBand(hNomBG_predSt, rootF, textF, strs);

  fclose(textF); // Text file is done now, close it
  rootF->Close(); // Close the output root file
}

//------------------------------------------------------------------------------
TH1* GetNC(IPrediction* specs, double POT, osc::OscCalcSterile* calc)
{
  Spectrum s = specs->PredictComponent(
    calc, Flavors::kAll, Current::kNC, Sign::kBoth
  );
  TH1* ret = s.ToTH1(POT);
  std::cout << "NC: " << ret->Integral(0,-1) << std::endl;
  ZeroError(ret);
  return ret;
}

//------------------------------------------------------------------------------
TH1* GetBG(IPrediction* specs, double POT, osc::OscCalcSterile* calc)
{
  Spectrum s = specs->PredictComponent(
    calc, Flavors::kAll, Current::kCC, Sign::kBoth
  );
  TH1* ret = s.ToTH1(POT);
  std::cout << " Bkg: " << ret->Integral(0,-1) << std::endl;
  ZeroError(ret);
  return ret;
}

//------------------------------------------------------------------------------
std::string MakeName(std::string chan, std::string det,
                     std::string syst, std::string sigma)
{
  // This pattern must match that in NusSystFromHist::LoadHists
  std::string ret = "h" + chan + "_" + det + "_" + syst + "_" + sigma;
  return ret;
}

//------------------------------------------------------------------------------
std::string MakeTitle(std::string chan, std::string det,
                      std::string syst, std::string sigma)
{
  std::string ret = chan + " " + det + " " + syst + " Systematic, "
                    + sigma + "#sigma";
  return ret;
}

//----------------------------------------------------------------------------
/// Plot FD syst shifted and nominal spectra
/// Should be able to change the input to IPrediction and make this redundant
/// This is for plotting an on/off style shift
void PlotSyst(TH1* nom, TH1* hp1, TH1* hm1,
              TDirectory* out, FILE* text, strings strs, bool use_max)
{
  AddErrorInQuadrature(nom, hp1, hm1, use_max);
  PlotSyst(nom, hp1, hm1, out, text, strs);
  return;
}

//----------------------------------------------------------------------------
/// Plot a nominal spectrum and a shifted spectra
/// This is for a systematic with multiple sigma levels
/// \param h The nominal spectrum
/// \param hp1 The shifted spectrum, +1 sigma
/// \param hm1 The shifted spectrum, -1 sigma
/// \param out The directory/file for root output
/// \param text The file for text output
/// \param strs Contains strings for pretty plots
void PlotSyst(TH1* h, TH1* hp1, TH1* hm1,
              TDirectory* out, FILE* text, strings strs)
{
  std::string xlabel = strs.fXLabel;
  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
  std::string title = strs.fDet + " " + strs.fComponent +
                      " Error for " + strs.fSystL + " Systematic";
  std::string fullTitle = ";" + xlabel + ";" + ylabel;
//  std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
  std::string name = "c" + strs.fComponent + strs.fDet +
                      strs.fSystS + strs.fSample;
  std::string ylabelRat = "Shifted / Nominal";
  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;

  TH1* hRat = (TH1*)h->Clone();
  TH1* hp1Rat = (TH1*)hp1->Clone();
  TH1* hm1Rat = (TH1*)hm1->Clone();

  int nBins = h->GetNbinsX();
  double nom_evts = h->Integral(1, nBins);
  double evts_diff = 0.;
  double percent_diff   = 100.;
  double percent_diff_p = 100.;
  double percent_diff_m = 100.;
  for(int i = 0; i <= nBins+1; ++i) {
    hRat->SetBinContent(i, 1.);

    if(i != 0 && i != nBins+1) {
      evts_diff += std::max(std::abs(hp1->GetBinContent(i) - h->GetBinContent(i)),
                            std::abs(hm1->GetBinContent(i) - h->GetBinContent(i)));
    }
  }
  percent_diff   *= (evts_diff/nom_evts);
  percent_diff_p *= ((hp1->Integral(1, nBins) - nom_evts)/nom_evts);
  percent_diff_m *= ((nom_evts - hm1->Integral(1, nBins))/nom_evts);

  fprintf(text, "%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
          strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
          percent_diff, percent_diff_p, percent_diff_m);

  hp1Rat->Divide(h);
  hm1Rat->Divide(h);

  double maxval = h->GetBinContent(h->GetMaximumBin());
  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
  maxval = std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));

  double minval = h->GetBinContent(h->GetMinimumBin());
  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
  minval = std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));

  if(maxval < minval) { std::swap(maxval, minval); }

  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
  maxvalRat = std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));

  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
  for(int i = 1; i <= nBins; ++i) {
    if(hp1Rat->GetBinContent(i) > 0.) {
      minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
    }
    if(hm1Rat->GetBinContent(i) > 0.) {
      minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
    }
  }

  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
    minvalRat = 1 - std::abs(maxvalRat - 1.);
    maxvalRat = 1 + std::abs(maxvalRat - 1.);
  }
  else {
    minvalRat = 1 - std::abs(1. - minvalRat);
    maxvalRat = 1 + std::abs(1. - minvalRat);
  }

  CenterTitles(h);
  h->SetLineColor(kBlack);
  h->SetMaximum(1.1*maxval);
  h->SetMinimum(0);
  h->SetTitle(fullTitle.c_str());

  CenterTitles(hp1);
  hp1->SetLineColor(kRed);
  hp1->SetLineStyle(2);
  hp1->SetMaximum(1.1*maxval);
  hp1->SetMinimum(0);
  hp1->SetTitle(fullTitle.c_str());

  CenterTitles(hm1);
  hm1->SetLineColor(kRed);
  hm1->SetLineStyle(2);
  hm1->SetMaximum(1.1*maxval);
  hm1->SetMinimum(0);
  hm1->SetTitle(fullTitle.c_str());

  CenterTitles(hRat);
  hRat->SetLineColor(kBlack);
  hRat->SetMaximum(1.1*maxvalRat);
  hRat->SetMinimum(0.9*minvalRat);
  hRat->SetTitle(fullTitleRat.c_str());

  CenterTitles(hp1Rat);
  hp1Rat->SetLineColor(kRed);
  hp1Rat->SetLineStyle(2);
  hp1Rat->SetMaximum(1.1*maxvalRat);
  hp1Rat->SetMinimum(0.9*minvalRat);
  hp1Rat->SetTitle(fullTitleRat.c_str());

  CenterTitles(hm1Rat);
  hm1Rat->SetLineColor(kRed);
  hm1Rat->SetLineStyle(2);
  hm1Rat->SetMaximum(1.1*maxvalRat);
  hm1Rat->SetMinimum(0.9*minvalRat);
  hm1Rat->SetTitle(fullTitleRat.c_str());

  double xL = 0.6, xR = 0.85;
  double yB = 0.6, yT = 0.85;
  TLegend* leg = new TLegend(xL, yB, xR, yT);
  leg->AddEntry(h,   "Nominal",      "l");
  leg->AddEntry(hp1, "+/- 1 #sigma", "l");

  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
  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();
  h->Draw("hist");
  hp1->Draw("hist same");
  hm1->Draw("hist same");
  leg->Draw();
  Simulation();

  pRatio->cd();
  hRat->Draw("hist");
  hp1Rat->Draw("hist same");
  hm1Rat->Draw("hist same");

  c->Update();
  out->WriteTObject(c);
  
  hRat->SetName( MakeName(strs.fComponent, strs.fDet, strs.fSystType, "0").c_str() );
  hRat->SetTitle(MakeTitle(strs.fComponent, strs.fDet, strs.fSystType, "0").c_str() );
  hp1Rat->SetName( MakeName (strs.fComponent, strs.fDet, strs.fSystType, "+1").c_str() );
  hp1Rat->SetTitle(MakeTitle(strs.fComponent, strs.fDet, strs.fSystType, "+1").c_str() );
  hm1Rat->SetName( MakeName (strs.fComponent, strs.fDet, strs.fSystType, "-1").c_str() );
  hm1Rat->SetTitle(MakeTitle(strs.fComponent, strs.fDet, strs.fSystType, "-1").c_str() );
  ZeroError(hRat);
  ZeroError(hp1Rat);
  ZeroError(hm1Rat);
  out->WriteTObject(hRat);
  out->WriteTObject(hp1Rat);
  out->WriteTObject(hm1Rat);
  return;
}

//------------------------------------------------------------------------------
void ResetAngles(osc::OscCalcSterile* calc)
{
  ResetSterileCalcToDefault(calc);
  calc->SetDm(4, 0.5);
  calc->SetAngle(2, 4, 0.17453293); // 10 degrees
  calc->SetAngle(3, 4, 0.61086524); // 35 degrees

  return;
}