FitSystEffectsAna.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 "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;

struct AngleValues{
  int nbins34;
  double min34;
  double max34;

  int nbins24;
  double min24;
  double max24;
};

double CalcSystEffect(double statonly, double withsyst);

double FindAngleFromChi2(double chi2, TH1* h);

void LoadMaps(TDirectory* dir,
              std::map<std::string, PredictionCombinePeriods*>* preds);

TH1* Plot1DSlice(IExperiment* expt, osc::OscCalcSterile* calc,
                 const IFitVar* var, std::vector<const IFitVar*> profVars,
                 TDirectory* rootOut, int nbins, double min, double max,
                 std::string name, Color_t color,
                 std::vector<const ISyst*> systs);

void Plot1DSlices(IExperiment* expt, osc::OscCalcSterile* calc,
                  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24, const IFitVar* vard24,
                  TDirectory* rootOut, AngleValues a,
                  std::string years, Color_t color,
                  std::string samplelabel,
                  std::map<std::string, TH1*>* h24s,
                  std::map<std::string, TH1*>* h34s,
                  std::vector<const ISyst*> systs = {});

void PrintAngles(FILE* text, std::string label,
                 double th2468, double th2490,
                 double th3468, double th3490);

void PrintSystEffect(FILE* text,
                     double th24stat68, double th2468,
                     double th24stat90, double th2490,
                     double th34stat68, double th3468,
                     double th34stat90, double th3490);

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

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

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

  // Set up path to file with filled CAFAna objects and for output of analysis
  std::string folder = "";
  std::string filenm = "pred_fitsysteffects_ana01_fix";

  std::string loadLocation = folder + filenm + ".root";
  std::string saveLocation = folder + filenm + "Ana.root";
  std::string textLocation = folder + 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);

  loadDir->cd((loadLocation + ":/sCosmic").c_str());
  Spectrum sCosmic = *Spectrum::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sData").c_str());
  Spectrum sData   = *Spectrum::LoadFrom(gDirectory);

  LoadMaps(rootL, &preds);

  rootL->Close(); // Don't forget to close the file!

  TH1* hCosmic = sCosmic.ToTH1(sCosmic.Livetime(), kLivetime);
  Spectrum sCos(hCosmic, sData.POT(), sData.Livetime());

  int nData = sData.Integral(sData.POT());
  std::cout << "Data spectrum: " << sData.POT()
            << " and " << nData << " events" << std::endl;
  double nCos = sCos.Integral(sCos.POT());
  std::cout << "Cos spectrum:  " << sCos.POT()
            << " and " << nCos  << " events" << std::endl;

  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc3f = DefaultSterileCalc(3);

  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  calc4f->SetNFlavors(4);
  calc4f->SetDm(4, 0.5); //Dm41 = 0.5 eV^2
  calc4f->SetAngle(2, 4, 0.);
  calc4f->SetAngle(3, 4, 0.);

  // Set up values for slice and surface plots
  AngleValues avals;
  avals.nbins34 = 500;
  avals.min34 = 0.;
  avals.max34 = 50.;
  avals.nbins24 = 450;
  avals.min24 = 0.;
  avals.max24 = 45.;

  const Color_t kFitColor = kRed;

  std::map<std::string, std::vector<const ISyst*> > systmap;
  systmap["OscParam"]   = std::vector<const ISyst*>();
  systmap["Beam"]   = std::vector<const ISyst*>();
  systmap["Birks"]  = std::vector<const ISyst*>();
  systmap["Calib"]  = std::vector<const ISyst*>();
  systmap["DataMC"] = std::vector<const ISyst*>();
  systmap["GENIE"]  = std::vector<const ISyst*>();
  systmap["NDCont"] = std::vector<const ISyst*>();
  systmap["NDRock"] = std::vector<const ISyst*>();
  systmap["Norm"]   = std::vector<const ISyst*>();
  systmap["Stats"]  = std::vector<const ISyst*>();

  systmap["OscParam"]  .push_back(&kNusOscParamSysts);
  systmap["Beam"]  .push_back(&kNusBeamSysts);
  systmap["Birks"] .push_back(&kNusBirksSyst);
  systmap["Calib"] .push_back(&kNusCalibFlatSyst);
  systmap["Calib"] .push_back(&kNusCalibRelSyst);
  systmap["Calib"] .push_back(&kNusCalibSlopeXSyst);
  systmap["Calib"] .push_back(&kNusCalibSlopeYSyst);
  systmap["DataMC"].push_back(&kNusNueCCSyst);
  systmap["DataMC"].push_back(&kNusNumuCCSyst);
  systmap["GENIE"] .push_back(&kNusGENIESmallSysts);
  systmap["NDCont"].push_back(&kNusNDContSyst);
  systmap["NDRock"].push_back(&kNusNDRockSyst);
  systmap["Norm"]  .push_back(&kNusNormSyst);
  systmap["Stats"] .push_back(&kNusMCStatsSyst);

  systmap["All"] = getAllNusSysts();

  std::map<std::string, TH1*> h24s;
  std::map<std::string, TH1*> h34s;

  // Open files for saving analysis output
  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
  FILE* textF;
  textF = fopen(textLocation.c_str(), "w");

  if(true) {
    std::cout << "NoSysts" << std::endl;
    ResetAngles(calc4f);

    GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.8, 3.2);
    CountingExperiment expt(&pred, sData, sCos);
    MultiExperiment multi({&th23Constraint, &expt});

    // Plot 1D fits for 1 data yr experiment
    Plot1DSlices(
      &multi, calc4f,
      &kFitTheta23InDegreesSterile,
      &kFitTheta34InDegreesSterile,
      &kFitTheta24InDegreesSterile,
      &kFitDelta24InPiUnitsSterile,
      rootF, avals, "1", kFitColor,
      "NoSysts", &h24s, &h34s
    );
  }

  double th24stat68 = FindAngleFromChi2(1.,    h24s["NoSysts"]);
  double th24stat90 = FindAngleFromChi2(2.706, h24s["NoSysts"]);
  double th34stat68 = FindAngleFromChi2(1.,    h34s["NoSysts"]);
  double th34stat90 = FindAngleFromChi2(2.706, h34s["NoSysts"]);
  PrintAngles(textF, "NoSysts", th24stat68, th24stat90, th34stat68, th34stat90);

  for(const auto& predpair : preds) {
    std::string label = predpair.first;
    std::cout << "Label: " << label << std::endl;
    ResetAngles(calc4f);

    GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.8, 3.2);
    CountingExperiment expt(predpair.second, sData, sCos);
    MultiExperiment multi({&th23Constraint, &expt});

    // Plot 1D fits for 1 data yr experiment
    Plot1DSlices(
      &multi, calc4f,
      &kFitTheta23InDegreesSterile,
      &kFitTheta34InDegreesSterile,
      &kFitTheta24InDegreesSterile,
      &kFitDelta24InPiUnitsSterile,
      rootF, avals, "1", kFitColor,
      label, &h24s, &h34s, systmap[label]
    );

    double th2468 = FindAngleFromChi2(1.,    h24s[label]);
    double th2490 = FindAngleFromChi2(2.706, h24s[label]);
    double th3468 = FindAngleFromChi2(1.,    h34s[label]);
    double th3490 = FindAngleFromChi2(2.706, h34s[label]);
    PrintAngles(textF, label, th2468, th2490, th3468, th3490);
    PrintSystEffect(textF,
                    th24stat68, th2468, th24stat90, th2490,
                    th34stat68, th3468, th34stat90, th3490);
  }
  
  fclose(textF); // Text file is done now, close it
  rootF->Close(); // Close the output root file
}

//------------------------------------------------------------------------------
double CalcSystEffect(double statonly, double withsyst)
{
  return 100.*(withsyst - statonly)/statonly;
}

//------------------------------------------------------------------------------
double FindAngleFromChi2(double chi2, TH1* h)
{
  for(int i = 1, n = h->GetNbinsX(); i < n; ++i) {
    if(h->GetBinContent(i, 1)    < chi2 &&
       h->GetBinContent(i+1, 1) >= chi2) {
      return h->GetXaxis()->GetBinCenter(i+1);
    }
  }

  return h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
}

//------------------------------------------------------------------------------
/// Load CAFAna objects from a TDirectory
/// The maps to contain these objects must already exist
void LoadMaps(TDirectory* dir,
              std::map<std::string, PredictionCombinePeriods*>* preds)
{
  TDirectory* tmp = gDirectory;
  TDirectory* keyDir = gDirectory;

  // Loop over all keys in this directory
  dir->cd();
  TIter nextkey(dir->GetListOfKeys());
  TKey* key;
  TKey* oldkey = 0;

  std::string sep = "__"; // This is the token separator
  std::string cafanaType = "";
  std::string sampleLabel = "";

  while((key = (TKey*)nextkey())) {
    // Keep only the highest cycle number for each key
    if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
    if(key->ReadObj()->IsFolder()) {
      keyDir = (TDirectory*)key->ReadObj();
    }
    else {
      continue;
    }

    std::string name(key->GetName());
    std::vector<std::string> tokens;
    std::size_t pos = 0, prev = 0;
    while(pos != std::string::npos) { // Loop over the TObject name
      pos = name.find(sep, prev); // Find the next separator

      if(pos != std::string::npos) {
        tokens.push_back(name.substr(prev, pos-prev));
        prev = pos+2; // Move past the current separator
      }
      else {
        tokens.push_back(name.substr(prev, name.length()-prev));
      }
    }

    if(tokens.size() > 1) {
      cafanaType = tokens[0];
      sampleLabel = tokens[1];
      if(cafanaType.compare("pred") == 0) {
        (*preds)[sampleLabel] = PredictionCombinePeriods::LoadFrom(keyDir).release();
      }
    }
    else { continue; }
  }

  tmp->cd();
  return;
}

//------------------------------------------------------------------------------
TH1* Plot1DSlice(IExperiment* expt, osc::OscCalcSterile* calc,
                 const IFitVar* var, std::vector<const IFitVar*> profVars,
                 TDirectory* rootOut, int nbins, double min, double max,
                 std::string name, Color_t color,
                 std::vector<const ISyst*> systs)
{
  // Create slice histogram with/without profiling based on input profVars vector size
  
  TH1* h = Profile(expt, calc, var, nbins, min, max, -1, profVars, systs);

  ResetAngles(calc); // Don't forget to reset calculator...

  // Formatting:
  h->SetLineColor(color); // Set the line color
  h->SetMaximum(5); // Want the plot to go from 0 to 5
  h->SetMinimum(0);
  h->SetName(name.c_str()); // Set the object name for saving
  rootOut->WriteTObject(h); // Save the raw histogram

  return h;
}

//------------------------------------------------------------------------------
void Plot1DSlices(IExperiment* expt, osc::OscCalcSterile* calc,
                  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24, const IFitVar* vard24,
                  TDirectory* rootOut, AngleValues a,
                  std::string years, Color_t color,
                  std::string samplelabel,
                  std::map<std::string, TH1*>* h24s,
                  std::map<std::string, TH1*>* h34s,
                  std::vector<const ISyst*> systs)
{
  // General form of the object name for saving
  std::string name = "h" + years + "yrTh";
  std::string fullname; // fullname adds 2 characters, the theta subscript indices
  
  // Plot profiled slices in same order as above
  // When slice is of one angle, the other two are profiled
  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
  fullname = name + "34Prof" + samplelabel;
  (*h34s)[samplelabel] = Plot1DSlice(
    expt, calc, var34, {var23, var24}, rootOut,
    a.nbins34, a.min34, a.max34, fullname, color, systs
  );
  fullname = name + "24Prof" + samplelabel;
  (*h24s)[samplelabel] = Plot1DSlice(
    expt, calc, var24, {var23, var34}, rootOut,
    a.nbins24, a.min24, a.max24, fullname, color, systs
  );
  
  return;
}

//------------------------------------------------------------------------------
void PrintAngles(FILE* text, std::string label,
                 double th2468, double th2490,
                 double th3468, double th3490)
{
  fprintf(text, "Angle limits for %s:\n", label.c_str());
  fprintf(text, "Th24: 68% (%): %4.2f, 90%%: %4.2f\n", th2468, th2490);
  fprintf(text, "Th34: 68% (%): %4.2f, 90%%: %4.2f\n", th3468, th3490);
  fprintf(text, "\n");
  return;
}

//------------------------------------------------------------------------------
void PrintSystEffect(FILE* text,
                     double th24stat68, double th2468,
                     double th24stat90, double th2490,
                     double th34stat68, double th3468,
                     double th34stat90, double th3490)
{
  fprintf(text, "Diff24: 68% (%) %4.2f, 90%%: %4.2f\n",
          CalcSystEffect(th24stat68, th2468),
          CalcSystEffect(th24stat90, th2490));
  fprintf(text, "Diff34: 68% (%) %4.2f, 90%%: %4.2f\n",
          CalcSystEffect(th34stat68, th3468),
          CalcSystEffect(th34stat90, th3490));
  fprintf(text, "\n");
  return;
}

//------------------------------------------------------------------------------
void ResetAngles(osc::OscCalcSterile* calc)
{
  ResetSterileCalcToDefault(calc);
  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
  calc->SetAngle(3, 4, 0.);
  calc->SetAngle(2, 4, 0.);

  return;
}