joint_fit_2019_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 "3FlavorAna/Ana2019/nue/joint_fit_2019_loader_tools.h"
//#include "./joint_fit_2019_loader_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 <algorithm>
#include <vector>
#include <string>

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 MakeMaps ( std::vector <const IFitVar*> profvars,
                std::map<const IFitVar*, TGraph*> &profMap,
                std::vector <const ISyst* > systs,
                std::map<const ISyst*, TGraph*> &systMap);
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);
void POTtag(bool both, bool FHCOnly, bool RHCOnly);
void AverageCyclicPoints(TGraph* g);
void AddCyclicPoints(TGraph* g);
TGraph* HistoToCyclicGraph(TH1* h, 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);



void joint_fit_2019_slices(bool createFile = false,
                           bool corrSysts = true,
                           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 fake2018 = options.Contains("fake2018");
  bool realData = options.Contains("realData");

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

  TString tag;
  if(FHCOnly) tag="FHCOnly/";
  if(RHCOnly) tag="RHCOnly/";
  if(both) tag = "RHC_and_FHC/";

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

  // official path commented so things aren't being overwritten
  TString outdir = "/nova/ana/nu_e_ana/Ana2019/Results/"+tag+"slices/FCInputs/";
  if (runOnGrid) outdir = "./";

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

  // outdir = "/nova/ana/nu_e_ana/Ana2019/Results/"+tag+"slices/";
  outdir="./results/blessed_plots/slices/";
  outdir += corrSysts?"syst/":"stat/";
  if (runOnGrid) outdir = "./";

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

  if(createFile){

   //////////////////////////////////////////////////
   // Load Nue and Numu experiments
   //////////////////////////////////////////////////
   //need numu only for prestage seeds

   struct predictions {
           const string name;
           const IPrediction * pred;
           std::pair <TH1D*, double> cos;
           double pot;
   };

   std::vector <predictions> preds;
   std::vector <predictions> preds_numu_only;
   std::vector <Spectrum * > data;
   std::vector <Spectrum * > data_numu_only;
   std::vector <const IExperiment * > expts;
   std::vector <const IExperiment * > expts_numu_only;

    auto calc_fake = DefaultOscCalc();
    if(fake2018) SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17);
    else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}

    if(!numuOnly) {
      if(FHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false, runOnGrid), GetNueCosmics2019("fhc", runOnGrid), GetPOT()});
      }
      if(RHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false, runOnGrid), GetNueCosmics2019("rhc", runOnGrid), GetPOT(false)});
      }
    }

   int nnumu = 4;
   if(!nueOnly) {
      if(FHCOnly || both ) {
        auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc", runOnGrid);
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc", runOnGrid);
        for (int i = 0; i< nnumu; i++ ){
            preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
            preds_numu_only.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
        }
      }
      if(RHCOnly || both ) {
        auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc", runOnGrid);
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc", runOnGrid);
        for (int i = 0; i< nnumu; i++ ){
            preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
            preds_numu_only.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
        }
       }
    }


   if(realData){
      if(!numuOnly){
        if(FHCOnly || both ) data.push_back(GetNueData2019("fhc", runOnGrid));
        if(RHCOnly || both ) data.push_back(GetNueData2019("rhc", runOnGrid));
      }
      if(!nueOnly){
        if(FHCOnly || both ){
          auto numu_data = GetNumuData2019(nnumu, "fhc", runOnGrid);
          data.insert(data.end(),numu_data.begin(), numu_data.end());
          data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
        if(RHCOnly || both ){
          auto numu_data = GetNumuData2019(nnumu, "rhc", runOnGrid);
          data.insert(data.end(),numu_data.begin(), numu_data.end());
          data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
      }
    }

    for(int i = 0; i < int(preds.size()); ++i){
        double POT = preds[i].pot;
        auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first);
        if(realData) thisdata = data[i];
        cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds[i].cos.first->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
        expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
    }

    cout<<"Make numu only experiment to get the seeds"<<endl;
    //make numu only experiment for seeds:

    for(int i = 0; i < (int) preds_numu_only.size(); ++i){
        double POT = preds_numu_only[i].pot;
        auto thisdata = GetFakeData (preds_numu_only[i].pred, calc_fake, POT, preds_numu_only[i].cos.first);
        if(realData) thisdata = data_numu_only[i];
        cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds_numu_only[i].cos.first->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
        expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i].pred, *thisdata, preds_numu_only[i].cos.first, preds_numu_only[i].cos.second));
    }
    ////////////////////////////////////////////////////////////
    // Add constraints, make experiments
    ////////////////////////////////////////////////////////////
    std::cout << "\nCreating multiexperiment\n" << std::endl;

    if (!th13Slice) {expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";}
    else {std::cout << "No reactor constraint, that's th13 slice\n";}
    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);

   std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
              << expts_numu_only.size() << " experiments\n\n" << std::endl;
    auto exptThis_numu_only = new MultiExperiment(expts_numu_only);

    ////////////////////////////////////////////////////////////
    // Systematics
    ////////////////////////////////////////////////////////////

    std::cout << "Systematics for the fit:\n";
    auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true);

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

    std::cout << "Systematics for the numu only fit:\n";
    auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, true, true, true);

    ////////////////////////////////////////////////////////////
    // 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 && !th23slice && !dmsqSlice && !th13Slice && octantSlice && (both || FHCOnly)){   //Cherenkov seed for deltaCP and neutron as well :(
      cout<<"Add Cherenkov seeds \n";
      for (double systshift:{-0.5, +0.5}){                                     //for the FHC Only or FHC+RHC ana
        SystShifts tempShifts (&kAnaCherenkovSyst, systshift);
        seedShifts.push_back(tempShifts);
      }
      if(both){
        cout<<"Add neutron systematic seeds \n";                               //For the joint RHC+FHC ana only
        for (double systshift:{-0.5, +0.5}){
        SystShifts tempShifts (&kNeutronVisEScalePrimariesSyst2018, systshift);
        seedShifts.push_back(tempShifts);
        }
      }
    }


  //////////////////////////////////////////////////////////////////////
  ///////////////////////// Seed controller ////////////////////////////
  //////////////////////////////////////////////////////////////////////

   cout<<"------------------- Start prestage seeds --------------------------"<<endl;

    double minchi_numu = 1E20;
    double pre_seed_th23;
    double pre_seed_dmsq;
    ResetOscCalcToDefault(calc);
    auxShifts.ResetToNominal();

   double maxmix = 0.514;  // from the numu results
   double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;

   for(int hie:{1}){
      std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
      Fitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32}, {});

      auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
                                        {{&kFitDmSq32,{hie*2.5e-3}},
                                        {&kFitSinSqTheta23, {0.4} }}, {}, //seedShifts,
                                        Fitter::kVerbose);

     if(thisminchi < minchi_numu) minchi_numu = thisminchi;

     pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
     pre_seed_dmsq = kFitDmSq32.GetValue(calc);
   }

  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));

  ResetOscCalcToDefault(calc);
  auxShifts.ResetToNominal();

  cout<<"------------------- End prestage seeds --------------------------"<<endl;


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

   //for NH:
    if(octantSlice) {
    th23seeds.push_back( { {0.499, numu_pre_seedLONH}, kFitSinSqTheta23LowerOctant, "LO"});
    th23seeds.push_back( { {0.501, numu_pre_seedUONH}, kFitSinSqTheta23UpperOctant, "UO"});
   }
   else th23seeds.push_back({ {numu_pre_seedLONH, 0.5, numu_pre_seedUONH}, &kFitSinSqTheta23, ""});

   std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
   std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};

   /////////////////////////////////////////////////////////////////////////
   ////////////////////////// Different best fit searches //////////////////
   /////////////////////////////////////////////////////////////////////////

    // Find the best fit points
    double minchi23 = 1E20;
    for(int hie:{-1, 1}){
      for (auto & thseed:th23seeds){

        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,
                                                &kFitDmSq32,
                                                &kFitSinSq2Theta13};

        Fitter fit23(exptThis, fitvars, systs);
        cout<<" pre dmsq seed is "<<pre_seed_dmsq<<endl;
        auto thisminchi = fit23.Fit(calc, auxShifts,
                                    {
                                      { &kFitDmSq32,   {hie*pre_seed_dmsq} },
                                      { thseed.var,          thseed.seeds },
                                      { &kFitDeltaInPiUnits, delta_seeds }
                                    }, seedShifts,
                                    Fitter::kVerbose);
        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=%.6f", thisminchi);
          shifts->SetTitle(str);
          gPad->Update();
          TLine *l=new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
          l->Draw("same");
          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 << "\n\n";

    cout<<"Check with oldstyle seeds "<<endl;
    double minchi23test = 1E20;
      for(int hie:{-1,1}){
        std::cout << "\n\nFinding test best fit "
                  << (hie>0? "NH " : "IH ")
                  << "\n\n";
        std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                &kFitSinSqTheta23,
                                                &kFitDmSq32,
                                                &kFitSinSq2Theta13};

        Fitter fit23(exptThis, fitvars, systs);
        auto thisminchi = fit23.Fit(calc,auxShifts ,{
            {&kFitDmSq32,{hie*fabs(kFitDmSq32.GetValue(calc))}},
            {&kFitSinSqTheta23, th23_old_seeds},
            {&kFitDeltaInPiUnits, delta_seeds}
          },  seedShifts
          );
        minchi23test= min(minchi23test, thisminchi);
        if(corrSysts){
        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=%.6f", thisminchi);
          shifts->SetTitle(str);
          gPad->Print(outdir  + "debug_slice_shifts_bestfit_" + suffix +
                      (hie>0? "NH": "IH") + "test_chisq_with_no_octant_differentiation.pdf");
        }

        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }//end hie
    std::cout << "\nFound overall test minchi " << minchi23test << std::endl;


    cout<<"<<<<<<<<<<<<<<<<<<<<<<<<<< Test with only systematics free"<<endl;
    minchi23test = 1E20;
    ResetOscCalcToDefault(calc);
    auxShifts.ResetToNominal();

    calc->SetdCP(0.16576*M_PI);
    calc->SetTh23(asin(sqrt(0.584509)));
    calc->SetDmsq32(2.50515e-3);
    calc->SetTh13(asin(sqrt(0.0824663))/2);

    Fitter fit23(exptThis, {}, systs);
    auto thisminchi = fit23.Fit(calc,auxShifts ,{}, seedShifts,
          Fitter::kVerbose);
   ofstream txtfile(outdir  + "debug_slice_shifts_bestfit_" + suffix +
                      + "test_chisq_with_no_octant_test_with_osc_par_fixed.txt");
   for (auto &v:systs){
       txtfile<<v->ShortName().c_str()<<"="<<auxShifts.GetShift(v)<<endl;
    }
    txtfile.close();
    minchi23test= min(minchi23test, thisminchi);
    auto shifts = PlotSystShifts(auxShifts);
          TString str = "Best fit " ;
          str+= TString::Format(" LL=%.6f", thisminchi);
          shifts->SetTitle(str);
          gPad->Print(outdir +  "debug_slice_shifts_bestfit_" + suffix +
                      + "test_chisq_with_no_octant_test_with_osc_par_fixed.pdf");
   ResetOscCalcToDefault(calc);
   auxShifts.ResetToNominal();
   cout<<"<<<<<<<<<<<< This is not true BF : "<< minchi23test<<endl;

  /////////////////////////////////////////////////////////////////////////////
  /////////////////////// Now, slices! ////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////

    //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");

    int steps = 60;// 80;

    //Now create all the slices
    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 && !th13Slice){
          profdef.fitvar =  &kFitDeltaInPiUnits;
          profdef.minX = 0;
          profdef.maxX = 2;
          profdef.profvars =  {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23};
          profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
          profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}}, {&kFitSinSq2Theta13, {0.082}} };
          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,delta_seeds}};
          profdef.shortname="th23";
        }

        if(th13Slice){
          profdef.fitvar = &kFitSinSq2Theta13;
          profdef.minX = 0.01;
          profdef.maxX = 0.2;
          profdef.profvars =  {&kFitDmSq32, kCurTheta23, &kFitDeltaInPiUnits};
          profdef.profvarnames =  {"DmSq32","SinSqTheta23","DeltaCPpi"};
          profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}} };
          profdef.shortname="th13";
        }

        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,delta_seeds}};
          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";

        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{
    TString infilename = "/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/slices/";
    infilename += corrSysts?"syst/":"stat/";
    infilename += "hist_slices_2019_";
    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");

    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.45}, kFitSinSqTheta23LowerOctant, "LO"});
       th23seeds.push_back( { {0.55}, kFitSinSqTheta23UpperOctant, "UO"});
     }
    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:{"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);
            }
            if(fccorr)
              graphs.push_back({new TGraph(h), hie + " " + thseed.label});
            else
              graphs.push_back({HistoToCyclicGraph(h, true),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(sliceName.Contains("th13")) {
          DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 0.05, 0.2);
        }
      }
     
      TString fcFolder="/nova/ana/nu_e_ana/Ana2019/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
          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]);
          }
        }
        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(outdir + "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(outdir + "points_slice_" +  suffix +
          // gPad->Print(TString("./test/") + [>(corrSysts?"syst/":"stat/")+<] "points_slice_" +  suffix +
                      (fccorr?"_fccorr":"") +
                      (fourierFit?"_smooth":"") +
                      ".pdf");
      }
    }//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, numuOnly, nueOnly, true, 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(outdir + "debug_slice_prof_" + suffix + "_"+ str + ".pdf");
      }//end profile vars
    }//end debug

  }//end plots
}//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 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 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", kAna2017FullDetEquivPOT/1E20).Data());}
  if (RHCOnly){ ltx = new TLatex(.52, 0.93, TString::Format("%2.4g#times10^{20} POT of #bar{#nu} data only", kAna2019RHCPOT/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}", kAna2017FullDetEquivPOT/1E20, kAna2019RHCPOT/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(TH1* h, 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*/)
    return HistoToCyclicGraph(hist, useRoot);
  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);
  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(fourierFit) {
    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