plot_joint_fit_2020_slices.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/Plots.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Experiment/Dmsq32Constraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/FC/FCSurface.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "OscLib/IOscCalc.h"
#include "Utilities/rootlogon.C"

#include "../joint_fit_2020_loader_tools.h"
#include "../joint_fit_2020_datarelease_tools.h"

#include "TCanvas.h"
#include "TBox.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.h"
#include "TText.h"
#include "TLatex.h"
#include "TPad.h"
#include "TLine.h"
#include "TMarker.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TGaxis.h"
#include "TGraphSmooth.h"

#include <algorithm>
#include <vector>
#include <string>
#include <fstream>
#include <sstream>

using namespace ana;


TGraph * signFuncGr;
Double_t signFunc(Double_t *x, Double_t *t);
void DrawSignMarker(double x, double y, int color);
void FindSignPoint(TGraph* gr, double sigma, double lowx, double highx, int color);
void HighlightSignPoints(TGraph* gr, TString grname);
void POTtag(bool both, bool FHCOnly, bool RHCOnly);
void AverageCyclicPoints(TGraph* g);
void AddCyclicPoints(TGraph* g);
TGraph* HistoToCyclicGraph(TGraph* g, bool useRoot);
void SmoothWithFourierFit(TH1* hist,
                          int nbins, double minX, double maxX,
                          int fourierFitN);
//why is this function so ugly
TGraph * SmoothWithFourierFit(TGraph* gr0, bool isDelta,
                              bool useRoot, int fourierFitN);
TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
                       bool fourierFit, std::string plot_name,
                       bool useRoot, bool fccol);

bool is_file_exist(const char *fileName)
{
    std::ifstream infile(fileName);
    return infile.good();
}


void plot_joint_fit_2020_slices(bool corrSysts = false,
                                bool runOnGrid = false,
                                TString options="joint_realData_both",
                                bool th23slice = false,
                                bool octantSlice = true,
                                bool dmsqSlice = false,
                                bool th13Slice = false,
                                bool fccorr = false,
                                bool fourierFit = true)

{
  bool nueOnly = options.Contains("nueOnly");
  bool numuOnly = options.Contains("numuOnly");
  bool joint = options.Contains("joint");
  assert (nueOnly || numuOnly || joint);

  bool FHCOnly = options.Contains("FHCOnly");
  bool RHCOnly = options.Contains("RHCOnly");
  bool both = options.Contains("both");
  assert (FHCOnly || RHCOnly || both);

  bool fake2019 = options.Contains("fake2019");
  bool realData = options.Contains("realData");

  auto suffix = options;
  if(corrSysts) suffix+="_systs";
  else suffix+="_stats";

  if(th23slice) suffix+="_th23";
  if(dmsqSlice) suffix+="_dmsq";
  if(th13Slice) suffix+="_th13";
  if(!octantSlice)suffix+="_noOct";

  TString outdir = "/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/";
  if (runOnGrid) outdir = "./";

  TString outFCfilename (outdir + "slices_FCInput_2020_" + suffix);
  if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice) outFCfilename = outdir + "slices_FCInput_2020_" +suffix;

  TString outfilename (outdir + "hist_slices_2020_" + suffix);

  TString infilename = outdir;

  infilename += "hist_slices_2020_";
  infilename += options;
  infilename += corrSysts?"_systs":"_stats";
  if(th23slice) infilename+="_th23";
  if(dmsqSlice) infilename+="_dmsq";
  if(th13Slice) infilename+="_th13";
  if(!octantSlice)infilename+="_noOct";
  infilename += ".root";

  std::cout << "Making plots from " << infilename << std::endl;
  TFile * infile = new TFile (infilename,"read");
  TFile * official =0;
  TString officialpreffix = "./";
  bool makedatarelease = true;
  if(makedatarelease &&
     !(corrSysts && options.Contains("joint") && fccorr )) {
     std::cerr << "\nWARNING: "
               << "makedatarelease is set to true"
               << "for a weird combination of options\n";
  }
  if(makedatarelease) {
     official = MakeDataReleaseFileSlice(officialpreffix,
                                         th23slice, dmsqSlice);
  }
  auto mins =* (TVectorD*)infile->Get("overall_minchi");
  double minchi23 = mins[0];

  std::vector <TString> sliceNames = {};
  if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {"slice_delta_"};
  if (th23slice) sliceNames = {"slice_th23_"};
  if (dmsqSlice) sliceNames = {"slice_dmsq_"};
  if (th13Slice) sliceNames = {"slice_th13_"};

  // it shouldn't play any role, but copy old style seeding here for sake of naming, everywhere is just .lable

  struct th23helper{
      std::vector<double> seeds;
      const IFitVar * var;
      TString label;
  };
  std::vector <th23helper> th23seeds;
  if(octantSlice) {
     th23seeds.push_back( { {0.55}, &kFitSinSqTheta23UpperOctant, "UO"});
     th23seeds.push_back( { {0.45}, &kFitSinSqTheta23LowerOctant, "LO"});
  }
  else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});

  //end of old-style-seeding

  for(TString sliceName:sliceNames){
      std::vector < std::pair <TGraph*, TString> > graphs;

      if(sliceName.Contains("delta")) {
         for (TString hie:{"IH", "NH"}){
              for(auto const & thseed:th23seeds){
                  //auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label); // previously was TH1 saved into the file
                  auto gr = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
                  if(!gr) {
                    std::cerr << "Problem " << sliceName + hie + thseed.label << "\n";
                    exit(1);
                  }
                  if(fccorr)
                     graphs.push_back({gr, hie + " " + thseed.label});
                  else
                     graphs.push_back({HistoToCyclicGraph(gr, true), hie + " " + thseed.label});
              }
         }
         DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 5, 0, 2);
      }
      else {
           for (TString hie:{"NH","IH"}) {
                for(auto const & thseed:th23seeds){
                    auto gr = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
                    if(!gr) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
                    graphs.push_back({gr, hie + " " + thseed.label});
                }
           }

           if(sliceName.Contains("th23")) {
              DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3,
                             (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
           }

           if(sliceName.Contains("dmsq")) {
              DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 2.15, 2.75);
           }
           if(sliceName.Contains("th13")) {
              DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 0.05, 0.2);
           }
      }
     
      TString fcFolder="/nova/ana/nu_e_ana/Ana2020/FC/slices/";
      std::string plot_name = "";
      if(th23slice) {
         plot_name = "ssth23";
         fcFolder += "th23/";
      }
      else if(th13Slice) {
              plot_name = "th13";
              fcFolder += "th13/";
      }
      else if(dmsqSlice) {
              plot_name = "dmsq";
              fcFolder += "dmsq/";
      }
      else {
         plot_name = "delta";
         fcFolder += "delta/";
      }

      if(fccorr) {
         if(!joint || !corrSysts)
         { std::cerr << "\nWe don't have FC\n\n"; exit(1);}

         // TString fcFolder="./FC/";
         for (auto &gr:graphs) {
              TString fcSliceN = plot_name;
              fcSliceN += gr.second.Contains("NH")? "_nh" : "_ih";
              if(!th23slice)
                 fcSliceN += gr.second.Contains("UO")? "uo" : "lo";
              fcSliceN += ".root";

               // only use fourier fit for period boundary conditions
               // Fourier fit doesn't like delta_UO_nh. Would need
               // lots of harmonics to fit it correctly anyway
               if (is_file_exist(fcFolder+fcSliceN)){
                   std::cout<<"Found FC file for this curve"<<std::endl;
                   bool useFourier = (fcSliceN.Contains("delta") &&
                                      !fcSliceN.Contains("delta_nhuo"));

                   gr.first = FCCorrectSlice(gr.first, fcFolder + fcSliceN,
                                             fourierFit, plot_name, !useFourier, true);
               }
         }

      }//end fccorr

      for (auto &gr:graphs) {
           //if(gr.second.Contains("IH")) gr.first->SetLineColor(kPrimColorIH);
           //if(gr.second.Contains("NH")) gr.first->SetLineColor(k3SigmaConfidenceColorNH);
           if(gr.second.Contains("IH")) gr.first->SetLineColor(cih_dark);
           if(gr.second.Contains("NH")) gr.first->SetLineColor(cnh);
           if(gr.second.Contains("LO")) gr.first->SetLineStyle(7);
           gr.first->SetLineWidth(3);
           if(dmsqSlice){
              for (int i =0; i < gr.first->GetN(); ++i){
                   gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
              }
           }
           TGraphSmooth *test = new TGraphSmooth("normal");
           TGraph * test_sm = test->SmoothKern(gr.first,"normal", 0.01);
           test_sm->SetLineColor(gr.first->GetLineColor());
           test_sm->SetLineStyle(gr.first->GetLineStyle());
           test_sm->SetLineWidth(3);
           test_sm->Draw("l same");
           //gr.first->Draw("l same");
      }

      POTtag(both, FHCOnly, RHCOnly);
      PreliminarySide();
      gPad->RedrawAxis();

      auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
      leg->Draw();

      //      outdir = "test_delta_slice/";
      //gPad->SetGrid();
      for(TString ext: {".pdf",".root"}) {
          gPad->Print("./slice_" +  suffix +
          // gPad->Print(TString("./test/") + [>(corrSysts?"syst/":"stat/")+<] "slice_" +  suffix +
                     (fccorr?"_fccorr":"") +
                     (fourierFit?"_smooth":"") +
                     ext);
      }//end ext

      bool annotate = 1;
      if(annotate) {

         for (auto &gr:graphs) {
               HighlightSignPoints(gr.first, sliceName);
         }
         gPad->SetGrid();
         gPad->Print("./points_slice_" +  suffix +
          // gPad->Print(TString("./test/") + [>(corrSysts?"syst/":"stat/")+<] "points_slice_" +  suffix +
                    (fccorr?"_fccorr":"") +
                    (fourierFit?"_smooth":"") +
                    ".pdf");
      }
      if(official)  GraphsToDataRelease(official, graphs);
  }//end slicename


  bool debug = false;

  if(debug){

     TFile * inFCfile = new TFile(outFCfilename+".root","read");
     std::cout << "\n\nOpen profiled " << outFCfilename+".root\n" << std::endl;

     std::vector <TString> strprof = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
     if(th23slice) strprof = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
     if(th13Slice) strprof = {"DmSq32","DeltaCPpi", "SinSqTheta23"};
     if(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
     
     auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, true);;

     for (const ISyst* s : systs){
          strprof.push_back(s->ShortName());
     }

     for(auto const & str:strprof){
         new TCanvas;
         double miny=-1.5, maxy=1.5;
         if(str=="DeltaCPpi")     {miny=0; maxy=2;}
         if(str=="DmSq32")        {miny=2.2; maxy=2.9;}
         if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
         if(str=="SinSqTheta23")  {miny=0.3; maxy=0.8;}

         auto axes = new TH2F("",";#delta_{CP}/#pi;" + str, 100, 0, 2, 100, miny, maxy);
         if(dmsqSlice) axes = new TH2F("",";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" + str, 100, 2.2, 2.8, 100, miny, maxy);
         if(th23slice) axes = new TH2F("",";sin^{2}#theta_{23};" + str, 100, 0.3, 0.7, 100, miny, maxy);
         if(th13Slice) axes = new TH2F("",";sin^{2}#theta_{13};" + str, 100, 0.01, 0.2, 100, miny, maxy);
         axes->Draw("hist");

         for (TString hie:{"NH","IH"}){
              for(auto const & thseed:th23seeds){
                  auto gr = (TGraph*)  inFCfile->Get(hie + thseed.label + "_" + str);
                  if(!gr) {
                     std::cerr << "Problem " << hie + thseed.label + "_" + str << "\n";
                     exit(1);
                  }
                  if(str=="DmSq32"){
                     for (int i = 0; i < (int) gr->GetN(); ++i){
                          gr->SetPoint(i, gr->GetX()[i], 1000* fabs(gr->GetY()[i]));
                     }
                  }
                  if(dmsqSlice){
                     for (int i = 0; i < (int) gr->GetN(); ++i){
                          gr->SetPoint(i, 1000*fabs(gr->GetX()[i]), gr->GetY()[i]);
                     }
                  }
                  //gr->SetLineColor(hie=="NH"?k2SigmaConfidenceColorNH:kPrimColorIH);
                  gr->SetLineColor(hie=="NH"?cnh:cih_dark);
                  if(thseed.label.Contains("LO")) gr->SetLineStyle(7);
                  gr->SetLineWidth(3);
                  gr->SetMarkerColor(gr->GetLineColor());
                  gr->Draw("l same");
              }//end th23
         }//end hie
         gPad->SetGrid();
         gPad->Print("./debug_slice_prof_" + suffix + "_"+ str + ".pdf");
     }//end profile vars
  }//end debug

}//end all


/// And supporting functions:

Double_t signFunc(Double_t *x, Double_t *t){

  double y = signFuncGr->Eval(x[0]);
  if(t[0] > 0 && abs(y - t[0]) > 0.05) return 999.;
  return abs(y - t[0]);
}
void DrawSignMarker(double x, double y, int color)
{
  TMarker * m = new TMarker (x, y,29);
  m->SetMarkerColor(color);
  m->SetMarkerSize(3);
  TString pstr = TString::Format("#bf{(%.2f,%.2f)}",x,y);
  TLatex * ltx = new TLatex(x,y,pstr.Data());
  ltx->SetTextAngle(45);
  ltx->SetTextColor(color);
  ltx->Draw();
  m->Draw();
}
void FindSignPoint(TGraph* gr, double sigma, double lowx=1, double highx=2, int color = kMagenta)
{
  double minX = gr->GetX()[0];
  double maxX = gr->GetX()[gr->GetN()];
  signFuncGr = gr;
  TF1 * f1 = new TF1("sign",signFunc,minX,maxX,1);
  f1->SetParameter(0,sigma);
  double xmin = f1->GetMinimumX(lowx,highx);

  if(f1->Eval(xmin) < 5) DrawSignMarker(xmin, gr->Eval(xmin), color);
}



void HighlightSignPoints(TGraph* gr, TString grname)
{
  struct SPoint{
    double sigma;
    double xlow;
    double xhigh;
  };
  std::vector <SPoint > sig_range;
  //fish for local minima points / 123 sgma points
  if(grname.Contains("dmsq")){
    sig_range = {{0.,2.3,2.6},
                 {1,2.3,2.4},{1,2.4,2.6},
                 {2,2.2,2.5},{2,2.5,2.7}};
  }
  if(grname.Contains("th23")){
    sig_range = {{0.,0,0.5},{0,0.5,0.7},{0,0.5,0.5001},
                 {1,0.3,0.45},{1,0.45,0.5},{1,0.5,0.515},{1,0.515,0.55},{1,0.54,0.7},
                 {2,0.3,0.5},{2,0.5,0.7}};
  }
  if(grname.Contains("th13")){
    sig_range = {{0.,0,0.1},{0,0.1,0.15},{0,0.15,0.2},
                 {1,0.0,0.1},{1,0.1,0.15},{1,0.15,0.2},
                 {2,0.0,0.1},{2,0.1,0.2}};
  }
  if(grname.Contains("delta")){
    sig_range = {{0.,1., 2.}, {0.,0., 1.},
                 {1., 0, 0.5},{1.,0.5, 1.},{1.,1, 1.7},{1.,1.7, 2.},
                 {2., 0, 0.5},{2.,0.5, 1.},{2.,1, 1.3},{2.,1.5, 2.},
                 {3., 0, 0.5}, {3,0.5,1.5}, {3,1.7,2.}};

  }
  int color = kMagenta;
  for (auto &sr:sig_range){
    if(sr.sigma == 1 ) color = kGreen + 1;
    if(sr.sigma == 2 ) color = kAzure ;
    if(sr.sigma == 3 ) color = kOrange ;
    FindSignPoint(gr,sr.sigma,sr.xlow, sr.xhigh, color);
  }
}

void POTtag(bool both, bool FHCOnly, bool RHCOnly)
{
  TLatex* ltx = new TLatex();
  auto ltx1 = new TLatex(.15, 0.93, TString::Format("NOvA FD"));
  if (FHCOnly){ ltx = new TLatex(.42, 0.93, TString::Format("%g#times10^{20} POT equiv. of #nu data only", kAna2020FHCPOT/1E20).Data());}
  if (RHCOnly){ ltx = new TLatex(.52, 0.93, TString::Format("%2.4g#times10^{20} POT of #bar{#nu} data only", kAna2020RHCPOT/1E20).Data());}
  if(both) {ltx = new TLatex(.32, 0.93, TString::Format("%g#times10^{20} POT equiv #nu + %2.4g#times10^{20} POT #bar{#nu}", kAna2020FHCFullDetEquivPOT/1E20, kAna2020RHCPOT/double(1E20)).Data());}
  ltx1->SetNDC();
  ltx1->SetTextAlign(11);
  ltx1->SetTextSize(0.8*kBlessedLabelFontSize);
  ltx1->Draw();
  ltx->SetNDC();
  ltx->SetTextAlign(11);
  ltx->SetTextSize(0.8*kBlessedLabelFontSize);
  ltx->Draw("same");
}

void AverageCyclicPoints(TGraph* g)
{
  double x0, y0;
  double x2pi, y2pi;
  g->GetPoint(g->GetN()-1, x2pi, y2pi);
  g->GetPoint(1, x0, y0);
  std::cout << "(x0, y0) = (" << x0 << ", " << y0 << ")\n";
  std::cout << "(x2pi, y2pi) = (" << x2pi << ", " << y2pi << ")\n";

  double avgY = ( y0 + y2pi ) / 2;
  std::cout << "Average = " << avgY << std::endl;
  g->SetPoint(0, 0, avgY);
  g->SetPoint(g->GetN(), 2, avgY);

  std::cout << "\nFound average cyclic point: (0, " << avgY << ")\n";
  std::cout << "\nAdded average cyclic points " << 0 << ": " << g->GetX()[0] << " " << g->GetY()[0];
  std::cout << ", " << g->GetN() << ": " << g->GetX()[g->GetN()-1] << " " << g->GetY()[g->GetN()-1] << "\n\n";
}



void AddCyclicPoints(TGraph* g)
{
  double x, y;
  g->GetPoint(g->GetN()-1, x, y);
  g->SetPoint(0, x-2, y);
  g->GetPoint(1, x, y);
  g->SetPoint(g->GetN(), x+2, y);

  std::cout << "\n Added cyclic points " << 0 << " " << g->GetX()[0]<< " " <<  g->GetY()[0]
            << ", " << g->GetN()-1 <<  x+2 <<  y << "\n\n";
}


TGraph* HistoToCyclicGraph(TGraph* g, bool useRoot = false)
{
 // TGraph* g = new TGraph;
 // g->SetPoint(0, 0, 0); // temporary
 // for(int i = 1; i <= h->GetNbinsX(); ++i)
 //   g->SetPoint(i, h->GetXaxis()->GetBinCenter(i), h->GetBinContent(i));

  //g->SetLineColor(h->GetLineColor());
  g->SetLineWidth(2);

  if(useRoot)
    AverageCyclicPoints(g);
  else
    AddCyclicPoints(g);

  return g;
}


void SmoothWithFourierFit(TH1* hist,
                          int nbins = 80, double minX = 0, double maxX = 2,
                          int fourierFitN=3)
{
  FitToFourier fitF (hist, minX, maxX, fourierFitN);
  auto f1= fitF.Fit();
  hist->Reset();
  hist->Add(f1);
}
//why is this function so ugly
TGraph * SmoothWithFourierFit(TGraph* gr0,
                              bool isDelta = true,
                              bool useRoot = false,
                              int fourierFitN = 2)
{
  int nbins = gr0->GetN();
  double minX = gr0->GetX()[0];
  double maxX = gr0->GetX()[nbins-1];
  //TGraph to center bin
  minX = minX  - (gr0->GetX()[1] - gr0->GetX()[0])/2;
  maxX = maxX  + (gr0->GetX()[1] - gr0->GetX()[0])/2;

  //Fourier fit  relies on periodicity, ignore extra points
  if(isDelta && !useRoot){
    nbins = 60;
    minX = 0;
    maxX = 2;
  }
  std::cout << "Smooth " << nbins << " nbins " << minX << ", " << maxX << "\n";

  TH1D * hist = new TH1D(UniqueName().c_str(),"",nbins,minX,maxX);

  for (int i = 0; i < nbins; ++i){
    double val = gr0->Eval(hist->GetBinCenter(i+1));
    hist->SetBinContent(i + 1, val * val); //use dchisq not chi2 to smooth
  }
  if(useRoot) hist->Smooth(2);
  else  SmoothWithFourierFit(hist, nbins, minX, maxX, fourierFitN);
  //now return sqrt
  for (int i = 1; i <= nbins; ++i){
    hist->SetBinContent(i, sqrt(hist->GetBinContent(i)));
  }
  /*if(isDelta && !useRoot) // NB: commented, previously HistoToCyclicGraph used TH1D, now there is a graph
    return HistoToCyclicGraph(hist, useRoot); // some corrections are needed
  else */return new TGraph(hist);
}

TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
                       bool fourierFit, std::string plot_name,
                       bool useRoot = false, bool fccol=false)
{
  bool th23plot = fcFileName.Contains("th23");
  bool dmsqplot = fcFileName.Contains("dmsq");
  bool dmsqnhuo = fcFileName.Contains("dmsq_nhuo");
  TString hiername = fcFileName.Contains("nh")? "NH":"IH";

  double minX = 0;
  double maxX = 2;
  if (th23plot) {minX = 0.3; maxX = 0.7;}
  if (dmsqplot) {
    minX = hiername=="NH" ? 2e-3 : -3e-3;
    maxX = hiername=="NH" ? 3e-3 : -2e-3;
  }

  //We want slice, not sqrt
  auto n= sqrtslice->GetN();

  auto sliceXH = new TGraph(n);
  for (int i = 0; i < n;i++ ){
    sliceXH->SetPoint(i,sqrtslice->GetX()[i],util::sqr(sqrtslice->GetY()[i]));
  }
  FCSurface *fcXH = 0;
  if(!fccol)  fcXH = FCSurface::LoadFromFile(fcFileName.Data()).release();
  else {
    auto fccol = FCCollection::LoadFromFile(fcFileName.Data()).release();
    //Only one Y bin required for slices.
    // All slices this year had 60 bins. Need to find a better way
    // next year
    fcXH = new FCSurface(60,minX,maxX,1,0,1);
    fcXH->Add(*fccol, plot_name);
  }


  TGraph* sigXH = new TGraph;
  TGraph* pcXH = new TGraph;

  TFile* bestfitfile = new TFile("/nova/ana/nu_e_ana/Ana2019/Results/BestFits/bestfits_joint_realData_both_systs.root", "read");
  auto calcbf = (osc::IOscCalcAdjustable*)(LoadFrom<osc::IOscCalc>(bestfitfile, "NH_UO_calc")).release();
  double dmsq32bf = kFitDmSq32.GetValue(calcbf);
  double th23bf = kFitSinSqTheta23.GetValue(calcbf);
  double deltabf = kFitDeltaInPiUnits.GetValue(calcbf);
  for(int i = 0; i < n; ++i){
    const double delta = sliceXH->GetX()[i];
    const double upXH = sliceXH->GetY()[i];

    // Fall back to gaussian up if dchisq is too large or FCBin is empty (mostly for deltaCP IH slices)
    if(upXH > 13 || fcXH->GetBin(i+1,1)->Empty()){
      std::cout << " Falling back on gaussian up value: " << upXH
                << " for bin " << i << " , delta = " << delta << std::endl;
      sigXH->SetPoint(sigXH->GetN(), delta, sqrt(upXH));
    }
    else {
      double pXH;
      if(upXH > 0) pXH = fcXH->GetBin(i+1, 1)->SignificanceForUpValue(upXH);
      else pXH = 0; //HACK HACK
      sigXH->SetPoint(sigXH->GetN(), delta, PValueToSigma(1-pXH));
   }
   if(dmsqnhuo && (i+1 < n) && (delta < dmsq32bf) && (dmsq32bf < sliceXH->GetX()[i+1])){
      sigXH->SetPoint(sigXH->GetN(), dmsq32bf, 0.);
   } 
   if ( special_smoothing && (i+1 < n) && (delta < deltabf) && (deltabf < sliceXH->GetX()[i+1])){
      std::cout << "Setting BFP at deltaCP = " << deltabf << std::endl;
      sigXH->SetPoint(sigXH->GetN(), deltabf, 0.);
   }
   if ( fcFileName.Contains("th23_nh") && (i+1 < n) && (delta < th23bf) && (th23bf < sliceXH->GetX()[i+1])){
      std::cout << "Setting BFP at th23 = " << th23bf << std::endl;
      sigXH->SetPoint(sigXH->GetN(), th23bf, 0.);
   } 
  }
  if(fourierFit && !th23plot) {
    sigXH = SmoothWithFourierFit(sigXH, fcFileName.Contains("delta"), useRoot);
  }

  return sigXH;
}

// th23
// cafe joint_fit_2019_slices.C false true false joint_realData_both true false false false true true

// dmsq
// cafe joint_fit_2019_slices.C false true false joint_realData_both false true true false true true

// delta
// cafe joint_fit_2019_slices.C false true false joint_realData_both false true false false true true