SpectrumParamEffectsAna.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/Prediction/PredictionCombinePeriods.h"
#include "NuXAna/macros/NuSPlotFunctions.h"

#include "OscLib/OscCalcSterile.h"

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

#include <string>

using namespace ana;

double CalcParamEffect(double nom, double shift);

void PlotEvents(TDirectory* out, std::string label,
                TH1* nom, TH1* hup, TH1* hdn);

void PrintEvents(FILE* text, std::string label,
                 double nom, double up, double dn);

void SetCalcParam(osc::OscCalcSterile* calc, std::string label, double shift);

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

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

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

  std::string loadLocation = "/nova/ana/steriles/Ana01/FitSystEffects.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 sCos = *Spectrum::LoadFrom(gDirectory);
  sCos.OverridePOT(kSecondAnaPOT);

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

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

  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  calc4f->SetNFlavors(4);
  ResetAngles(calc4f);

  std::map<std::string, std::pair<double, double> > parammap;
  parammap["t13"] = {-0.3675*M_PI/180., 0.3714*M_PI/180.};
  parammap["d13"] = {-0.27, 0.38};
  parammap["m32"] = {-0.00006, 0.00006};

  std::map<std::string, TH1*> hups;
  std::map<std::string, TH1*> hdns;
  std::map<std::string, double> nups;
  std::map<std::string, double> ndns;

  Spectrum snom = pred.Predict(calc4f) + sCos;
  TH1* hnom = snom.ToTH1(kSecondAnaPOT);
  double nnom = snom.Integral(kSecondAnaPOT);
  double scale = 83.509/nnom;
  nnom *= scale;

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

  for(const auto& parampair : parammap) {
    std::string label = parampair.first;
    double upshift = parampair.second.second;
    double dnshift = parampair.second.first;
    ResetAngles(calc4f);

    SetCalcParam(calc4f, label, upshift);
    Spectrum sup = pred.Predict(calc4f) + sCos;
    hups[label] = sup.ToTH1(kSecondAnaPOT);
    nups[label] = sup.Integral(kSecondAnaPOT)*scale;

    ResetAngles(calc4f);

    SetCalcParam(calc4f, label, dnshift);
    Spectrum sdn = pred.Predict(calc4f) + sCos;
    hdns[label] = sdn.ToTH1(kSecondAnaPOT);
    ndns[label] = sdn.Integral(kSecondAnaPOT)*scale;

    PrintEvents(textF, label, nnom, nups[label], ndns[label]);
    PlotEvents( rootF, label, hnom, hups[label], hdns[label]);
  }
  
  fclose(textF); // Text file is done now, close it
  rootF->Close(); // Close the output root file
}

//------------------------------------------------------------------------------
double CalcParamEffect(double nom, double shift)
{
  return 100.*(shift - nom)/nom;
}

//------------------------------------------------------------------------------
void PlotEvents(TDirectory* out, std::string label,
                TH1* nom, TH1* hup, TH1* hdn)
{
  double max = nom->GetBinContent(nom->GetMaximumBin());
  max = std::max(max, hup->GetBinContent(hup->GetMaximumBin()));
  max = std::max(max, hdn->GetBinContent(hdn->GetMaximumBin()));

  std::string title = ";Calorimetric Energy (GeV);Events / 0.25 GeV";

  SetHistOptions(nom, max, title, -404, kBlack, false);
  SetHistOptions(hup, max, title, -404, kRed,   false);
  SetHistOptions(hdn, max, title, -404, kBlue,  false);

  hup->SetLineStyle(2);
  hdn->SetLineStyle(3);

  std::string ptype = label.substr(0,1);

  std::string param = "";
  if(ptype.compare("t") == 0) {
    param = "#theta_{" + label.substr(1,2) + "} ";
  }
  if(ptype.compare("d") == 0) {
    param = "#delta_{" + label.substr(1,2) + "} ";
  }
  if(ptype.compare("m") == 0) {
    param = "#Deltam^{2}_{" + label.substr(1,2) + "} ";
  }

  double xL = 0.6, xR = 0.8, yB = 0.65, yT = 0.85;
  TLegend* leg = new TLegend(xL, yB, xR, yT);
  SetLegendOptions(leg);
  leg->AddEntry(nom, "Nominal", "l");
  leg->AddEntry(hup, (param + "+1#sigma").c_str(), "l");
  leg->AddEntry(hdn, (param + "-1#sigma").c_str(), "l");
  leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);

  std::string ctitle = "c" + label;
  TCanvas* c = new TCanvas(ctitle.c_str(), ctitle.c_str(), 800, 500);
  nom->Draw("hist");
  hup->Draw("hist same");
  hdn->Draw("hist same");
  leg->Draw();

  gPad->RedrawAxis();
  out->WriteTObject(c);
  return;
}

//------------------------------------------------------------------------------
void PrintEvents(FILE* text, std::string label,
                 double nom, double up, double dn)
{
  fprintf(text, "Number of events for %s:\n", label.c_str());
  fprintf(text, "Nom: %4.2f, Up: %4.2f, Down: %4.2f\n", nom, up, dn);
  fprintf(text, "Shifts: Up: %4.2f%%, Down: %4.2f%%\n\n",
          CalcParamEffect(nom, up),
          CalcParamEffect(nom, dn));
  return;
}

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

  return;
}

//------------------------------------------------------------------------------
void SetCalcParam(osc::OscCalcSterile* calc, std::string label, double shift)
{
  std::string ptype = label.substr(0,1);
  int val1 = std::stoi(label.substr(1,1));
  int val2 = std::stoi(label.substr(2,1));

  if(val1 > val2) { std::swap(val1, val2); }

  if(ptype.compare("t") == 0) {
    calc->SetAngle(val1, val2, calc->GetAngle(val1, val2) + shift);
  }
  else if(ptype.compare("d") == 0) {
    calc->SetDelta(val1, val2, calc->GetDelta(val1, val2) + shift);
  }
  else if(ptype.compare("m") == 0) {
    if(val2 != 1) {
      calc->SetDm(val2, calc->GetDm(val2) + shift);
    }
  }
  else {
    std::cout << "Unknown parameter type." << std::endl;
  }

  return;
}