joint_fit_2017_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/Fit/FrequentistSurface.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 "3FlavorAna/Prediction/PredictionSystNue2017.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "OscLib/IOscCalc.h"
#include "3FlavorAna/Ana2017/nue/joint_fit_2017_tools.h"

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

#include <algorithm>
#include <vector>
#include <string>

using namespace ana;


TGraph * signFuncGr;
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.55,0.7},
                 {2,0.3,0.5},{2,0.5,0.7}};
  }
  if(grname.Contains("delta")){
    sig_range = {{0.,1., 2.},
                 {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 MakeMaps ( std::vector <const IFitVar*> profvars,
                std::map<const IFitVar*, TGraph*> &profMap, 
                std::vector <const ISyst* > systs,
                std::map<const ISyst*, TGraph*> &systMap)
{
  
  for (const IFitVar * var:profvars) 
              profMap.insert(std::pair<const IFitVar*, TGraph*> (var, new TGraph()));
  for (const ISyst* syst : systs)
    systMap.insert(std::pair<const ISyst*,   TGraph*> (syst, new TGraph()));
}

void SaveMaps ( TDirectory * dir, 
                TString prefix, 
                std::vector <const IFitVar*> profvars,
                std::vector <TString > profnames,
                std::map<const IFitVar*, TGraph*> &profMap, 
                std::vector <const ISyst* > systs,
                std::map<const ISyst*, TGraph*> &systMap)
{
  TDirectory *tmp = gDirectory;
  dir->cd();
  for (int i = 0; i < (int) profvars.size(); ++i){
    profMap[profvars[i]]->Write((prefix + "_" + profnames[i]));
  }
  for (const ISyst* s : systs)
    systMap[s]->Write((prefix + "_" + s->ShortName()));
  tmp->cd();
}


void DrawSliceCanvas(TString slicename, const double ymax,
                     const double xmin = 0,  const double xmax = 2.){
  auto c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetFillStyle(0);
  c1->SetLeftMargin(0.14);
  c1->SetBottomMargin(0.15);
  TString title;
  if(slicename.Contains("delta")) {
    title = ";delta;Significance (#sqrt{#Delta#chi^{2}})";
  }
  if(slicename.Contains("th23")) {
    title = ";sin^{2}#theta_{23};Significance (#sqrt{#Delta#chi^{2}})";
  }
  if(slicename.Contains("dmsq")) {
    title = ";|#Deltam^{2}_{32}| (10^{-3} eV^{2});Significance (#sqrt{#Delta#chi^{2}})";
  }
  auto axes = new TH2F("",title, 100, xmin, xmax, 100, 0, ymax);
  if(slicename.Contains("fccorr")) axes->GetYaxis()->SetTitle("Significance (#sigma)");
  CenterTitles(axes);
  axes->Draw();
  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetXaxis()->SetTitleOffset(0.8);
  axes->GetYaxis()->SetTitleOffset(0.75);
  if(slicename.Contains("delta")) XAxisDeltaCPLabels(axes);
  c1->RedrawAxis();
}

TLegend * SliceLegend(std::vector <std::pair <TGraph*, TString > > & graphs, bool isDelta)
{
  TLegend * leg;
  if(isDelta) leg = new TLegend(0.7, 0.6, 0.89,0.89);
  else leg = new TLegend(0.7, 0.2, 0.89, 0.49);
  leg->SetTextSize(kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  for (auto &gr:graphs){
    leg->AddEntry(gr.first,gr.second, "l");
  }
  return leg;
}


void POTtag()
{
  TLatex* ltx = new TLatex(.16, .81, TString::Format("#splitline{NOvA FD}{%g#times10^{20} POT equiv.}", kAna2017FullDetEquivPOT/1E20).Data());
  ltx->SetNDC();
  ltx->SetTextAlign(11);
  ltx->SetTextSize(kBlessedLabelFontSize);
  ltx->Draw();
}

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() <<  x+2 <<  y << "\n\n";
}


TGraph* HistoToCyclicGraph(TH1* h)
{
  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);

  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 = 80;
    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);
  //  hist->Smooth(10);   

  //now return sqrt
  for (int i = 1; i <= nbins; ++i){
    hist->SetBinContent(i, sqrt(hist->GetBinContent(i)));
  }
  if(isDelta && !useRoot) return HistoToCyclicGraph(hist);
  else return new TGraph(hist);
}

TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName, 
                       bool fourierFit, bool useRoot = false, bool fccol=false)
{
  bool th23plot = fcFileName.Contains("th23");
  bool dmsqplot = fcFileName.Contains("dmsq");
  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(); 
    fcXH = new FCSurface(n,minX,maxX,3,0,1);
    fcXH->Add(*fccol); 
  }
  TGraph* sigXH = new TGraph;
  TGraph* pcXH = new TGraph;

  for(int i = 0; i < n; ++i){
    const double delta = sliceXH->GetX()[i];
    const double upXH = sliceXH->GetY()[i];

    //HACK HACK
    if((!fcFileName.Contains("delta") && upXH > 11)){
      sigXH->SetPoint(sigXH->GetN(), delta, sqrt(upXH));
    }
    else {
      double pXH;
      if(upXH > 0) pXH = fcXH->GetBin(i+1, 2)->SignificanceForUpValue(upXH);
      else pXH = 0; //HACK HACK
      sigXH->SetPoint(sigXH->GetN(), delta, PValueToSigma(1-pXH));
   }
  }
  if(fourierFit) {
    sigXH = SmoothWithFourierFit(sigXH, fcFileName.Contains("delta"), useRoot);
  }

  return sigXH; 
}

void joint_fit_2017_slices(std::string decomp = "prop",
                           bool createFile=false,
                           bool corrSysts = false,
                           TString options="joint",
                           bool th23slice = false,
                           bool octantSlice = true,
                           bool dmsqSlice = false,
                           bool fccorr=false,
                           bool fourierFit = false)
{
  bool nueOnly = options.Contains("nueOnly");
  bool numuOnly = options.Contains("numuOnly");
  bool joint = options.Contains("joint");
  assert (nueOnly || numuOnly || joint);

  bool fakeNHLO = options.Contains("fakeNHLO");
  bool fakeNHUO = options.Contains("fakeNHUO");
  bool fake2017 = options.Contains("fake2017");
  bool realData = options.Contains("realData");

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

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

  TString outdir = "/nova/ana/nu_e_ana/Ana2017/Results/slices/";
  //  TString outdir = "slices/";
  TString outfilename (outdir + "hist_slices_2017_" + suffix);
  TString outFCfilename (outdir + "slices_FCInput_2017_" + suffix);
  if(joint && corrSysts && !th23slice && ! dmsqSlice) 
    outFCfilename = outdir + "slices_FCInput_2017_" + decomp;

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

  if(createFile){
    
    //////////////////////////////////////////////////
    // Load Nue and Numu experiments
    //////////////////////////////////////////////////

    std::vector <const IPrediction * > preds;
    std::vector <std::pair <TH1D*, double > > cosmics;
    std::vector <Spectrum * > data;
    std::vector <const IExperiment * > expts;

    auto calc_fake = DefaultOscCalc();
    if(fakeNHLO) SetFakeCalc(calc_fake, 0.404, 2.7e-3, 1.48);
    else if(fakeNHUO) SetFakeCalc(calc_fake, 0.623, 2.7e-3, 0.74);
    else if(fake2017) SetFakeCalc(calc_fake, 0.56, 2.43e-3, 1.21);
    else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}

    if(!numuOnly) {
      preds.push_back(GetNuePrediction2017(decomp, DefaultOscCalc(), corrSysts));
      cosmics.push_back(GetNueCosmics2017());
    }

    int nnumu = 4;

    if(!nueOnly) {
      auto numu_preds = GetNumuPredictions2017(nnumu, corrSysts);
      preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
      auto numu_cosmics = GetNumuCosmics2017(nnumu);
      cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
    }

    if(realData){
      if(!numuOnly){
        data.push_back(GetNueData2017());
      }
      if(!nueOnly){
        auto numu_data = GetNumuData2017(nnumu);
        data.insert(data.end(),numu_data.begin(), numu_data.end());
      }
    }
    for(int i = 0; i < (int) preds.size(); ++i){
      auto thisdata = GetFakeData (preds[i],calc_fake, kAna2017POT, cosmics[i].first);
      if(realData) thisdata = data[i];
      expts.push_back(new SingleSampleExperiment(preds[i],*thisdata,
                                                 cosmics[i].first,cosmics[i].second));
    }

    ////////////////////////////////////////////////////////////
    // Add constraints, make experiments
    ////////////////////////////////////////////////////////////
    std::cout << "\nCreating multiexperiment\n" << std::endl;

    std::cout << "Adding WorldReactorConstraint2017\n";
    expts.push_back(WorldReactorConstraint2017());
    if(nueOnly) {
      std::cout << "Adding Dmsq32ConstraintPDG2017\n";
      expts.push_back(&kDmsq32ConstraintPDG2017);
    }
    std::cout << "Creating Multiexperiment with a total of "
              << expts.size() << " experiments\n\n" << std::endl;
    auto exptThis = new MultiExperiment(expts);

    ////////////////////////////////////////////////////////////
    // Systematics
    ////////////////////////////////////////////////////////////
    std::cout << "Systematics for the fit:\n";
    auto systs = GetJointFitSystematicList(corrSysts, !numuOnly, !nueOnly, true);

    std::cout << "\n\nSystematic correlations...\n";
    if(!nueOnly && ! numuOnly && corrSysts){
      exptThis->SetSystCorrelations(0, GetCorrelations(true));
      auto notfornumu = GetCorrelations(false);
      for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
    }

    ////////////////////////////////////////////////////////////
    // Fit
    ////////////////////////////////////////////////////////////
    std::cout << "Starting the fit" << std::endl;

    osc::IOscCalcAdjustable* calc = DefaultOscCalc();

    SystShifts auxShifts = SystShifts::Nominal();

    // In case some systs need different seeds
    std::vector <SystShifts> seedShifts = {};
    if(corrSysts){
      for (double systshift:{+0.5,-0.5}){
        SystShifts tempShifts (kMECq0ShapeSyst2017,systshift);
        seedShifts.push_back(tempShifts);
      }
    }

    // Find the best fit points
    double minchi23 = 1E20;

    for (auto & thseed: th23seeds){
      for(int hie:{-1, 1}){

        std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
                  << thseed.label << ", "
                  << "ssth23 seeds ";
        for (auto const & s:thseed.seeds) std::cout << s << ", ";
        std::cout << std::endl;

        std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                thseed.var, 
                                                &kFitDmSq32Scaled, 
                                                &kFitSinSq2Theta13};
        
        MinuitFitter fit23(exptThis, fitvars, systs);

        auto thisminchi = fit23.Fit(calc, auxShifts,
                                    {
                                      { &kFitDmSq32,        {hie * fabs(calc->GetDmsq32())} },
                                      { thseed.var,          thseed.seeds },
                                      { &kFitDeltaInPiUnits,{0.,0.5, 1., 1.5} }
                                    }, seedShifts);
        minchi23= min(minchi23, thisminchi);

        if(corrSysts){
          // Plot the systematic pulls and label with BFV
          auto shifts = PlotSystShifts(auxShifts);
          TString str = "Best fit " ;
          for (auto &v:fitvars){
            str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
          }
          str+= TString::Format(" LL=%.3f", thisminchi);
          shifts->SetTitle(str);
          gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
                      (hie>0? "NH": "IH") + thseed.label + ".pdf");
        }

        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }
    }
    std::cout << "\nFound overall minchi " << minchi23 << std::endl;

    //We'll save to this file, make pretty plots later
    TFile * outfile = new TFile(outfilename+".root","recreate");
    TFile * outFCfile = new TFile(outFCfilename+".root","recreate");

    outfile->cd();
    TVectorD v(1);
    v[0] = minchi23;
    v.Write("overall_minchi");

    //HACK HACK
    int steps = 80;// 5;

    //Now create all the surfaces
    for(int hie: {-1, +1}){

      ResetOscCalcToDefault(calc);
      calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

      std::cout << "Starting profile " << (hie>0? "NH ": "IH") << "\n\n";

      struct ProfDef{
        const IFitVar * fitvar;
        double minX;
        double maxX;
        std::vector <const IFitVar * > profvars;
        std::vector <TString> profvarnames;
        std::map  <const IFitVar *, std::vector <double> >profseeds;
        TString shortname;
      };

      for(auto const & thseed:th23seeds){  
        ProfDef profdef;
        const IFitVar* kCurTheta23 =  thseed.var; 

        if(!th23slice && !dmsqSlice){ 
          profdef.fitvar =  &kFitDeltaInPiUnits;
          profdef.minX = 0;
          profdef.maxX = 2;
          profdef.profvars =  {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23};
          profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
          profdef.profseeds = {{kCurTheta23, thseed.seeds}};
          profdef.shortname="delta";
        }

        if(th23slice){
          profdef.fitvar = &kFitSinSqTheta23;
          profdef.minX = (nueOnly ? 0:0.3);
          profdef.maxX = (nueOnly ? 1:0.7);
          profdef.profvars =  {&kFitDmSq32, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
          profdef.profvarnames =  {"DmSq32","SinSq2Theta13","DeltaCPpi"};
          profdef.profseeds = {{&kFitDeltaInPiUnits,{0.5,1.,1.5}}};
          profdef.shortname="th23";
        }

        if(dmsqSlice){
          profdef.fitvar = &kFitDmSq32;
          profdef.minX = hie > 0 ? 2e-3 : -3e-3;
          profdef.maxX = hie > 0 ? 3e-3 : -2e-3;
          profdef.profvars =   {kCurTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
          profdef.profvarnames = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
          profdef.profseeds = { {kCurTheta23, thseed.seeds},
                                {&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}, };
          profdef.shortname="dmsq";
        }

        std::map<const IFitVar*, TGraph*> profVarsMap;
        std::map<const ISyst*, TGraph*> systMap;
        MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
        
        std::cout << " Profile "<< thseed.label << "\n\n";
        
        auto slice = SqrtProfile(exptThis, calc, 
                                 profdef.fitvar, steps, profdef.minX, profdef.maxX,
                                 minchi23,
                                 profdef.profvars,
                                 systs, 
                                 profdef.profseeds,
                                 seedShifts,
                                 profVarsMap, systMap);
        //Save!
        TString hieroctstr = (hie > 0 ? "NH" : "IH") + thseed.label;
        
        outfile->cd();
        slice->Write((TString ) "slice_" + profdef.shortname + "_" + hieroctstr);
        SaveMaps (outFCfile, hieroctstr, 
                  profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
                
      }//end seeds
    }//end hie
  
    delete outfile;
    std::cout << "\nProfiles saved to " << outfilename << std::endl;
    return;
  }//end createfile




  //////////////////////////////////////////////////////////////////////
  // Make plots
  //////////////////////////////////////////////////////////////////////
  else{ 
    std::cout << "Making plots from " << outfilename+".root" << std::endl;
    TFile * infile = new TFile (outfilename+".root","read");

    auto mins =* (TVectorD*)infile->Get("overall_minchi");
    double minchi23 = mins[0];
    
    std::vector <TString> sliceNames = {};
    if (!th23slice && ! dmsqSlice) sliceNames = {"slice_delta_"};
    if (th23slice) sliceNames = {"slice_th23_"};
    if (dmsqSlice) sliceNames = {"slice_dmsq_"};
    
    for(TString sliceName:sliceNames){
      std::vector < std::pair <TGraph*, TString> > graphs;
      
      if(sliceName.Contains("delta")) {
        for (TString hie:{"NH","IH"}){
          for(auto const & thseed:th23seeds){
            auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
            if(!h) {
              std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; 
              exit(1);
            }
            graphs.push_back({HistoToCyclicGraph(h),hie + " " + thseed.label});
          }
        }
        DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 5, 0, 2);
      }
      else {
        for (TString hie:{"NH","IH"}) {
          for(auto const & thseed:th23seeds){
            auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
            if(!h) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
            graphs.push_back({new TGraph(h),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(fccorr) {
        if(!joint || !corrSysts || decomp != "combo")
        { std::cerr << "\nWe don't have FC\n\n"; exit(1);}
        
        TString fcFolder="/nova/ana/nu_e_ana/Ana2017/FC/";
        for (auto &gr:graphs) {
          TString fcSliceN =  "FCCol_";
          fcSliceN += gr.second ;
          if(th23slice) fcSliceN += "_slice_ssth23";
          else if(dmsqSlice) fcSliceN += "_slice_dmsq32";
          else fcSliceN += "_slice_delta";
          fcSliceN.ReplaceAll (" LO","_LO"); //gah
          fcSliceN.ReplaceAll (" UO","_UO"); 
          fcSliceN.ReplaceAll (" ",""); 
          bool useFourier = (fcSliceN.Contains("delta") && 
                             (fcSliceN.Contains("IH") /*|| fcSliceN.Contains("LO")*/)) ;
          gr.first = FCCorrectSlice(gr.first, fcFolder + fcSliceN + ".root",
                                    fourierFit, !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("UO")) gr.first->SetLineStyle(7);
        gr.first->SetLineWidth(3);
        if(dmsqSlice){ //HACK
          for (int i =0; i < gr.first->GetN(); ++i){
            gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
          }
        }
        gr.first->Draw("l same");
      }
      
      gPad->RedrawAxis();
      
      auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice);
      leg->Draw();
      
      POTtag();
      //gPad->SetGrid();
      for(TString ext: {".pdf",".root"}) {
        gPad->Print(outdir + "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(outdir + "points_slice_" +  suffix +
                      (fccorr?"_fccorr":"") +
                      (fourierFit?"_smooth":"") + 
                      ".pdf");
      }
    }//end slicename

    
    bool debug = true;

    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(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};

      auto systs = GetJointFitSystematicList(corrSysts, !numuOnly, !nueOnly, 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.8;}
        if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
        if(str=="SinSqTheta23")  {miny=0.4; maxy=0.6;}
        
        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);
        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"?k3SigmaConfidenceColorNH:kPrimColorIH);
            if(thseed.label.Contains("UO")) gr->SetLineStyle(7);
            gr->SetLineWidth(3);
            gr->SetMarkerColor(gr->GetLineColor());
            gr->Draw("l same");
          }//end th23
        }//end hie
        gPad->SetGrid();
        gPad->Print(outdir + "debug_slice_prof_" + suffix + "_"+ str + ".pdf");
      }//end profile vars
    }//end debug
    
  }//end plots
}//end all