PlotNus17PredSystsData.C

Go to the documentation of this file.

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

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Core/IFitVar.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "NuXAna/Systs/Nus17Systs.h"
#include "OscLib/OscCalcSterile.h"

#include "NuXAna/macros/NuSPlotFunctions.h"

using namespace ana;

// Definitions for organization

// One struct to contain all values for angle fits
struct AngleValues{
  int nbins23;
  double min23;
  double max23;

  int nbins34;
  double min34;
  double max34;

  int nbins24;
  double min24;
  double max24;
};

void TexInfo();

TLegend* MakeLegend(double x1, double y1, double x2, double y2, double textsize = 0.03)
{
  TLegend* leg = new TLegend(x1, y1, x2, y2);
  leg->SetBorderSize(0);
  leg->SetFillColor(kWhite);
  leg->SetFillStyle(0);
  leg->SetFillStyle(0);
  leg->SetTextFont(42);
  leg->SetTextSize(textsize);
  return leg;
}

void SigmaLines(double xmin, double xmax);

void PlotText(const double xpos, const double start, const double step);

// Reset calculator values after each fit to preset defaults
void ResetAngles(osc::OscCalcSterile* calc);

// Plot a single 1D slice canvas -- this function does the actual plotting
// The macro does not call this function directly in the main function
// \param expt The CAFAna experiment object, contains FD prediction and "data"
// \param calc Sterile calculator for oscillating the spectra
// \param var Parameter to fit for!
// \param profVars Vector of parameters to profile over
// \param rootOut Already open file for root output
// \param name Base of the string for creating the name for saving the root objects

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


void CompareSlices(TH1* stats, TH1* systs, std::string angle, TDirectory* rootOut, double xmin, double xmax);

// Plot 2D slice surfaces
// The macro does not call this function directly
// \param rootOut Already open file for root output
// \param prof String indicating whether surfaces are profiled or not (optional)
void Plot2DSlice(IExperiment* expt1,
                 IExperiment* expt2,
                 osc::OscCalcSterile* calc,
                 const IFitVar* varA,
                 int nbinsA, double minA, double maxA,
                 const IFitVar* varB,
                 int nbinsB, double minB, double maxB,
                 TDirectory* rootOut, std::string angles,
                 std::vector<const IFitVar*> vars,
                 std::vector<const ISyst*> systs = {});




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

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

  std::string loadLocation = folder + filenm + ".root";
  std::string saveLocation = "plot_" + filenm + "_shape_2D_lower.root";
  std::string textLocation = "count_" + filenm + "_shape_2D_lower.txt";

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

  // Load in the NuMI data spectrum
  Spectrum* intime = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Nus17/fout_nus17_box_opening_restricted.root","spec_nus_vise_numi").release();

  TH1 *hSpectrum = intime->ToTH1(intime->POT());
  Spectrum data_spec(hSpectrum, intime->POT(), intime->Livetime());
  std::cout << "Data spectrum: " << data_spec.POT()
            <<" POT and "
            << data_spec.ToTH1(intime->POT())->GetSumOfWeights()
            << " events" << std::endl;
  delete hSpectrum;
  TFile *cos = new TFile("$NOVA_ANA/steriles/Nus17/nus17_cosmic_background.root","read");
  TH1   *hcosmic = (TH1*)cos->Get("cosall");
  Spectrum cosmic(hcosmic, kAna2017POT, kAna2017Livetime);
  std::cout << "Cosmic Data spectrum: " << cosmic.POT()
            <<" POT and "
            << cosmic.ToTH1(cosmic.POT())->GetSumOfWeights()
            << " events"
            << ", Livetime: "
            << cosmic.Livetime() << std::endl;
  delete hcosmic;
  cos->Close();
  
  // Load filled objects from file
  TFile* rootL = new TFile(loadLocation.c_str(), "UPDATE");
  TDirectory* tmpL = gDirectory;
  TDirectory* loadDir = gDirectory;

  loadDir->cd((loadLocation + ":/nus17_predI").c_str());
  std::unique_ptr<PredictionInterp> pred = PredictionInterp::LoadFrom(gDirectory);
  
  
  tmpL->cd();
  rootL->Close(); // Don't forget to close the file!

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

  Spectrum prediction       = pred.get()->Predict(calc3f);
    

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


  std::string labelEReco = "Visible Energy (GeV)";

  // // Open files for saving analysis output
  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
  
  // Create Experiment objects, one each for the 1 data yr and 3 data yr predictions
  GaussianConstraint s2th23Constraint(&kFitSinSqTheta23Sterile, 0.404, 0.030);
  GaussianConstraint dmsq32Constraint(&kFitDmSq32Sterile, 2.67e-3, 0.11e-3);
  // 0.404 +- 0.030 --> dCP = 1.48pi
  // 0.624 +- 0.030 --> dCP = 0.74pi
    
  SingleSampleExperiment expt9s(pred.get(), data_spec, cosmic);
  MultiExperiment multi9s({&s2th23Constraint, &dmsq32Constraint, &expt9s});
  //MultiExperiment multi9s({&s2th23Constraint, &expt9s});
    
  // Set up values for slice and surface plots
  AngleValues avals;
  avals.nbins23 = 200;
  avals.min23 = 20.;
  avals.max23 = 70.;
  avals.nbins34 = 200;
  avals.min34 = 0.;
  avals.max34 = 50.;
  avals.nbins24 = 180;
  avals.min24 = 0.;
  avals.max24 = 45.;
    
  std::string fullname, special;
  
  const Color_t kStatsFitCol = kBlack;
  const Color_t kRateCol = kAzure+2;
  const Color_t kShapeCol = kBlack;
  
  const Color_t kFitColor = kRed;
  
    // fullname = "34";
    // special  = "statC_24e20";
    // //Plot 1D fits for 9e20 experiment: Counting
    // TH1* h34C_1 = Plot1DSlice(&multi9s, calc4f,
    //                        &kFitTheta34InDegreesSterile,
    //                        {&kFitSinSqTheta23Sterile,
    //                            &kFitTheta24InDegreesSterile,
    //                            &kFitDelta24InPiUnitsSterile},
    //                        avals.nbins34, avals.min34, avals.max34,
    //                        rootF, fullname, special, kRateCol,
    //                        -1
    //                        );
                              
    
    // special = "systS_9e20";
    // //Plot 1D fits for 9e20 experiment: Sample
    // TH1* h34S_1 = Plot1DSlice(&multi9s, calc4f,
    //                        &kFitTheta34InDegreesSterile,
    //                        {&kFitSinSqTheta23Sterile,
    //                            &kFitTheta24InDegreesSterile,
    //                            &kFitDelta24InPiUnitsSterile},
    //                        avals.nbins34, avals.min34, avals.max34,
    //                        rootF, fullname, special, kShapeCol,
    //                        -1, systs
    //                        );
    
    // CompareSlices(h34C_1, h34S_1, "34_stats_9e20", rootF,
    //                    avals.min34, avals.max34);
    
    // fullname = "24";
    // special  = "statC_9e20";
    // //Plot 1D fits for 9e20 experiment: Counting
    // TH1* h24C_1 = Plot1DSlice(&multi9s, calc4f,
    //                        &kFitTheta24InDegreesSterile,
    //                        {&kFitSinSqTheta23Sterile,
    //                            &kFitTheta34InDegreesSterile,
    //                            &kFitDelta24InPiUnitsSterile},
    //                        avals.nbins24, avals.min24, avals.max24,
    //                        rootF, fullname, special, kRateCol,
    //                        -1
    //                        );
                              
    // //fullname = "24";
    // special = "systS_9e20";
    // //Plot 1D fits for 9e20 experiment: Sample
    // TH1* h24S_1 = Plot1DSlice(&multi9s, calc4f,
    //                        &kFitTheta24InDegreesSterile,
    //                        {&kFitSinSqTheta23Sterile,
    //                            &kFitTheta34InDegreesSterile,
    //                            &kFitDelta24InPiUnitsSterile},
    //                        avals.nbins24, avals.min24, avals.max24,
    //                        rootF, fullname, special, kShapeCol,
    //                        -1, systs
    //                        );
    
    // CompareSlices(h24C_1, h24S_1, "24_shape_9e20", rootF,
    //                    avals.min24, avals.max24);
    
    
    // //Lower the number of bins for surfaces
    
    avals.nbins34 = 50;
    avals.min34 = 0.;
    avals.max34 = 50.;
    avals.nbins24 = 45;
    avals.min24 = 0.;
    avals.max24 = 45.;
    
    
    special = "34vs24_9e20_shape";
    //Plot the surfaces for 9e20 rate vs. shape experiment
    Plot2DSlice(&multi9s, &multi9s, calc4f,
                &kFitTheta34InDegreesSterile,
                avals.nbins34, avals.min34, avals.max34,
                &kFitTheta24InDegreesSterile,
                avals.nbins24, avals.min24, avals.max24,
                rootF, special,
                {&kFitSinSqTheta23Sterile,
                    &kFitDelta24InPiUnitsSterile},
                systs
                );
    
    rootF->Close(); // Close the output root file
    
}


//........................................................................
TH1* Plot1DSlice(IExperiment* expt, osc::OscCalcSterile* calc,
                 const IFitVar* var, std::vector<const IFitVar*> profVars,
                 int nbins, double min, double max,
                 TDirectory* rootOut, std::string indices,
                 std::string special, Color_t color,
                 double minchi, 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);
  //TH1* h = systs.size() > 0 ?  Profile(expt, calc, var, nbins, min, max, -1, profVars, systs, { {&kFitDelta24InPiUnitsSterile,{1.5} } }) :
  //Profile(expt, calc, var, nbins, min, max, -1, profVars, {}, { {&kFitDelta24InPiUnitsSterile, {1.5} } });
  ResetAngles(calc); // Don't forget to reset calculator...
  TH1* h = systs.size() > 0 ?  Profile(expt, calc, var, nbins, min, max, -1, profVars, systs) :
  Profile(expt, calc, var, nbins, min, max, -1, profVars);


  std::string fullname = indices + special;
  // Formatting
  h->SetLineColor(color);
  h->SetMaximum(5);
  h->SetMinimum(0);
  h->SetName(fullname.c_str()); // Object name for saving
  rootOut->WriteTObject(h);     // Save raw histogram


  std::cout << fullname << std::endl;
  std::string cname  = "c1D" + fullname;
  std::string ctitle = "Theta" + indices + "_ChiSq1DSlice_" + special;
  std::string htitle = "#theta_{"+indices+"} (deg.)";
  h->GetXaxis()->SetTitle(htitle.c_str());

  // Actually draw things together
  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
  h->GetXaxis()->SetRangeUser(min+0.25,max);
  h->Draw("l");
  TexInfo();
  SigmaLines(min, max);
  PlotText(0.135, 0.82, 0.055);
  TLegend* leg = 0;
  // Special case legend position for each angle slice
  if(indices.compare("34") == 0){
    std::cout << "Making 34 legend" << std::endl;
    leg = MakeLegend(0.135, 0.35, 0.435, 0.50, 0.037);
  }
  if(indices.compare("24") == 0)
    leg = MakeLegend(0.600, 0.35, 0.900, 0.50, 0.037);
  if(indices.compare("23") == 0)
    leg = MakeLegend(0.350, 0.35, 0.650, 0.50, 0.037);

  if(systs.size() > 0)
    leg->AddEntry(h,"syst.","L");
  else
    leg->AddEntry(h,"stat.","L");
  
  leg->Draw();

  //Simulation();
  Preliminary();
  CenterTitles(h);
  rootOut->WriteTObject(c);

  return h;
}


void CompareSlices(TH1* stats, TH1* systs, std::string indices, TDirectory* rootOut, double min, double max)
{

  std::string cname  = "c1DComp" + indices;
  std::string ctitle = "Theta" + indices + " ChiSq1Dcomp";
  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
  stats->GetXaxis()->SetRangeUser(min+0.25,max);
  systs->GetXaxis()->SetRangeUser(min+0.25,max);
  stats->Draw("l");
  systs->Draw("l same");
  TexInfo();
  SigmaLines(min, max);
  //Simulation();
  Preliminary();
  PlotText(0.135, 0.82, 0.055);
  TLegend* leg = 0;
  // Special case legend for angle slice
  if(indices.compare("34") == 0)
    leg = MakeLegend(0.135, 0.35, 0.435, 0.50, 0.037);
  if(indices.compare("24") == 0)
    leg = MakeLegend(0.600, 0.35, 0.900, 0.50, 0.037);
  if(indices.compare("23") == 0)
    leg = MakeLegend(0.350, 0.35, 0.650, 0.50, 0.037);
  else leg = MakeLegend(0.135, 0.35, 0.435, 0.50, 0.037);
  leg->AddEntry(stats,"stat.","L");
  leg->AddEntry(systs,"syst.","L");
  leg->Draw();
  rootOut->WriteTObject(c);

  return;
}


void Plot2DSlice(IExperiment* expt1,
                 IExperiment* expt2,
                 osc::OscCalcSterile* calc,
                 const IFitVar* varA,
                 int nbinsA, double minA, double maxA,
                 const IFitVar* varB,
                 int nbinsB, double minB, double maxB,
                 TDirectory* rootOut, std::string indices,
                 std::vector<const IFitVar*> vars,
                 std::vector<const ISyst*> systs)
{

  ResetAngles(calc);
  std::string cname      = "c2D" + indices;
  std::string cnamestat  = "c2D" + indices + "statC";
  std::string cnamesyst  = "c2D" + indices + "statS";
  std::string hnamestat  = "h" + indices + "statC";
  std::string hnamesyst  = "h" + indices + "statS";
  std::string ctitle     = "Surface for Theta" + indices;
  std::string ctitlestat = "Surface for Theta" + indices + " statC";
  std::string ctitlesyst = "Surface for Theta" + indices + " systS";

  FrequentistSurface surfStat(expt1, calc,
                   varA, nbinsA, minA, maxA,
                   varB, nbinsB, minB, maxB,
                   vars//, {}, { {&kFitDelta24InPiUnitsSterile, {1.5} } }
                   );
  ResetAngles(calc);
  
  FrequentistSurface surfSyst(expt2, calc,
                   varA, nbinsA, minA, maxA,
                   varB, nbinsB, minB, maxB,
                   vars, systs//, { {&kFitDelta24InPiUnitsSterile, {1.5} } }
                   );
  ResetAngles(calc);

  TH1* hSurfStat = surfStat.ToTH2();     // Create a histogram from the surfaces
  hSurfStat->SetName(hnamestat.c_str()); // Set the name for saving
  rootOut->WriteTObject(hSurfStat);      // Save the raw histogram

  TH1* hSurfSyst = surfSyst.ToTH2();     // Create a histogram from the surfaces
  hSurfSyst->SetName(hnamesyst.c_str()); // Set the name for saving
  rootOut->WriteTObject(hSurfSyst);      // Save the raw histogram

  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str());

  surfStat.DrawContour(Gaussian68Percent1D(surfStat), kDashed, kAzure+2);
  surfStat.DrawContour(Gaussian90Percent1D(surfStat), kSolid,  kAzure+2);
  surfSyst.DrawContour(Gaussian68Percent1D(surfSyst), kDashed, kBlack);
  surfSyst.DrawContour(Gaussian90Percent1D(surfSyst), kSolid,  kBlack);

  PlotText(0.135, 0.82, 0.055);
  TLegend* legcomp = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);
  TH1* dummy   = new TH1F("","",1,0,1);

  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kDashed);
  dummy->SetLineWidth(2);
  legcomp->AddEntry(dummy->Clone(),"68% C.L. (stat.)","L");
  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kSolid);
  legcomp->AddEntry(dummy->Clone(),"90% C.L. (stat.)","L");
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kDashed);
  legcomp->AddEntry(dummy->Clone(),"68% C.L. (syst.)","L");
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kSolid);
  legcomp->AddEntry(dummy->Clone(),"90% C.L. (syst.)","L");
  legcomp->Draw();
  //Simulation();
  Preliminary();
  rootOut->WriteTObject(c);

  // Plot Experiment 1
  TCanvas* cStat = new TCanvas(cnamestat.c_str(), ctitlestat.c_str());

  surfStat.DrawContour(Gaussian68Percent1D(surfStat), kDashed, kAzure+2);
  surfStat.DrawContour(Gaussian90Percent1D(surfStat), kSolid,  kAzure+2);

  PlotText(0.135, 0.82, 0.055);
  TLegend* legStat = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);

  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kDashed);
  dummy->SetLineWidth(2);

  legStat->AddEntry(dummy->Clone(),"68% C.L. (stat.)","L");
  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kSolid);
  legStat->AddEntry(dummy->Clone(),"90% C.L. (stat.)","L");

  legStat->Draw();
  //Simulation();
  Preliminary();
  rootOut->WriteTObject(cStat);

  //Plot Experiment 2
  TCanvas* cSyst = new TCanvas(cnamesyst.c_str(), ctitlesyst.c_str());

  surfSyst.DrawContour(Gaussian68Percent1D(surfSyst), kDashed, kBlack);
  surfSyst.DrawContour(Gaussian90Percent1D(surfSyst), kSolid,  kBlack);

  PlotText(0.135, 0.82, 0.055);
  TLegend* legSyst = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kDashed);
  dummy->SetLineWidth(2);
  legSyst->AddEntry(dummy->Clone(),"68% C.L. (syst.)","L");
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kSolid);
  legSyst->AddEntry(dummy->Clone(),"90% C.L. (syst.)","L");
  legSyst->Draw();
  //Simulation();
  Preliminary();
  rootOut->WriteTObject(cSyst);

  return;
}

void PlotText(const double xpos, const double start, const double step)
{
  TLatex* tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(xpos, start, "NOvA 8.85 #times 10^{20} POT-equiv.");
  tex->DrawLatex(xpos, start - 1*step, "#Deltam^{2}_{32} = 2.67#times10^{-3} eV^{2}");
  tex->DrawLatex(xpos, start - 2*step, "#theta_{13} = 8.5#circ, ^{}sin^{2}#theta_{23} = 0.404");
  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{41} = 0.5 eV^{2}, ^{}#delta_{13} = 1.48#pi");
}

void TexInfo()
{
  const double xtex = 0.76;
  const double y68  = 0.28;
  const double y90  = 0.55;

  TLatex *tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(xtex, y68, "68% C.L."); // 1-sigma
  tex->DrawLatex(xtex, y90, "90% C.L."); // 2-sigma
}

void SigmaLines(double xmin, double xmax)
{
  TLine* line68 = new TLine();
  TLine* line90 = new TLine();
  // Fixes horizontal effects at max canvas value
  TLine* line05 = new TLine();
  // Dashed line for 1 and 2 sigma lines
  line68->SetLineStyle(kDashed);
  line90->SetLineStyle(kDashed);

  line68->SetLineWidth(2);
  line90->SetLineWidth(2);
  line05->SetLineWidth(2);
  line68->DrawLine(xmin, 1,    xmax, 1);
  line90->DrawLine(xmin, 2.71, xmax, 2.71);
  line05->DrawLine(xmin, 5,    xmax, 5);
}


//........................................................................
void ResetAngles(osc::OscCalcSterile* calc)
{
  ResetSterileCalcToDefault(calc);
  calc->SetDm(4, 0.5);          // #deltam41 = 0.5 eV^2
  //calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
  calc->SetAngle(3, 4, 0); // 0, icecube
  calc->SetAngle(2, 4, 0); // 0, sk
  calc->SetDelta(1, 3, 1.48*M_PI);
  //calc->SetDelta(2, 4, 3*M_PI/2);
  calc->SetDm(3, 7.53e-5 + 2.67e-3);
  calc->SetAngle(2, 3, (asin(sqrt(0.404))));

  return;
}