demoFitSlices.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 "CAFAna/nue/Ana2018/joint_fit_2018_tools.h"
#include "Utilities/rootlogon.C"

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

void POTtag(bool both, bool FHCOnly, bool RHCOnly)
{
  TLatex* ltx;
  auto ltx1 = new TLatex(.16, 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("%1.2g#times10^{20} POT of #bar{#nu} data only", kAna2018RHCPOT/1E20).Data());}
  if(both) {ltx = new TLatex(.35, 0.93, TString::Format("%g#times10^{20} POT equiv #nu + %1.2g#times10^{20} POT #bar{#nu}", kAna2017FullDetEquivPOT/1E20, kAna2018RHCPOT/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 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 demoFitSlices(bool createFile = false,         /// Do you want to make fit and write it into the file or just plot existing slice
                   bool corrSysts = true,           /// syst or stat only?
                   TString options="joint_both",    /// what kind of analysis? numuOnly||nueOnly||joint ana for RHCOnly||RHCOnly||both modes
                   bool th23slice = false,          /// is it th23 slice?
                   bool octantSlice = true,         /// dou you want x2 curves split by octant?
                   bool dmsqSlice = false,          /// is it dmsq32 slice?
                   bool th13Slice = false)          /// is it th13 slice?
                                                    /// in case of th23slice = false && dmsqSlice = false && th13Slice = false it will be deltaCP slice

{
  ///First part follows demoFitContour logic, so I won't comment the loading and BF searches much
  
  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);

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

  TString outdir = "";
  TString outfilename (outdir +"hist_slices_2017_" + suffix);
  TString outFCfilename (outdir + "slices_FCInput_2018_" + suffix);   /// need for syst. debug profiles
  if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice) 
    outFCfilename = outdir + "slices_FCInput_2018_"+ suffix;


  if(createFile){
    
    //////////////////////////////////////////////////
    // Load Nue and Numu experiments
    //////////////////////////////////////////////////
    //need numu only for prestage seeds
    std::vector <const IPrediction * > preds;
    std::vector <const IPrediction * > preds_numu_only;
    std::vector <std::pair <TH1D*, double > > cosmics;
    std::vector <std::pair <TH1D*, double > > cosmics_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();
    SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17);

    if(!numuOnly) {
      if(FHCOnly || both ) {
         preds.push_back(GetNuePrediction2018("combo", DefaultOscCalc(), corrSysts, "fhc", false));
         cosmics.push_back(GetNueCosmics2018("fhc"));   
      } 
      if(RHCOnly || both ) {
         preds.push_back(GetNuePrediction2018("prop", DefaultOscCalc(), corrSysts, "rhc", false));
         cosmics.push_back(GetNueCosmics2018("rhc"));   
      } 
    }

    int nnumu = 4;
    if(!nueOnly) {
      if(FHCOnly || both ) {

        auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "fhc");
        preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
        preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());

        auto numu_cosmics = GetNumuCosmics2018(nnumu, "fhc");
        cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end()); 
        cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
      }
      if(RHCOnly || both ) {

        auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "rhc");
        preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
        preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());

        auto numu_cosmics = GetNumuCosmics2018(nnumu, "rhc");
        cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end()); 
        cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
      } 
    }

    cout<<"\n\n Predictions used for the fit:\n";

    for(int i = 0; i < (int) preds.size(); ++i){
      double POT;
      if(FHCOnly) POT = kAna2018FHCPOT;
      if(RHCOnly) POT = kAna2018RHCPOT;
      if(both) {
        if(nueOnly || both) {
           if(i==0) POT = kAna2018FHCPOT;
           if(i==1) POT = kAna2018RHCPOT;
        }
        if(numuOnly) {
           if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
           if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
         }
         if(both) {
           if(i >= 2 && i < 6) POT = kAna2018FHCPOT;
           if(i >= 6 && i < 10) POT = kAna2018RHCPOT;
        }       
      }
      auto thisdata = GetFakeData (preds[i],calc_fake, POT, cosmics[i].first);
      cout<<i<<" POT "<<POT<<" tot MC "<<preds[i]->Predict(calc_fake).Integral(POT)<<" cos "<<cosmics[i].first->Integral()<<" cos er "<<cosmics[i].second<<" analyze data "<<thisdata->Integral(POT)<<endl;
      expts.push_back(new SingleSampleExperiment(preds[i],*thisdata, cosmics[i].first,cosmics[i].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;
      if(FHCOnly) POT = kAna2018FHCPOT;
      if(RHCOnly) POT = kAna2018RHCPOT;
      if(both) {
           if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
           if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
      }
      auto thisdata = GetFakeData (preds_numu_only[i],calc_fake, POT, cosmics_numu_only[i].first);
      cout<<i<<" POT "<<POT<<" tot MC "<<preds_numu_only[i]->Predict(calc_fake).Integral(POT)<<" cos "<<cosmics_numu_only[i].first->Integral()<<" cos er "<<cosmics_numu_only[i].second<<" analyze data "<<thisdata->Integral(POT)<<endl;
      expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i],*thisdata, cosmics_numu_only[i].first,cosmics_numu_only[i].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";

    std::vector< const ISyst * > systs =  {};
    std::vector< const ISyst * > systs_numu_only =  {};
    if (corrSysts){      // syst case
       systs =  {&kAnaCalibrationSyst};                //here will be your own syst, which you made previously
       systs_numu_only =  {&kAnaCalibrationSyst};      //here will be your own syst, which you made previously  
    }

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

    osc::IOscCalcAdjustable* calc = DefaultOscCalc();

    SystShifts auxShifts = SystShifts::Nominal();

    // In case some systs need different seeds
    // If you run script with 2018 syst list you'll need seeds for some systematics
    // Comment for the workshop, but leave here for the education purposes
    
    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;

    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 - "dot" plot for systematics
          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";
/////////////////////////////////////////////////////////////////////////////
/////////////////////// 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"); 
    /// again, save BF into the file
    outfile->cd();
    TVectorD v(1);
    v[0] = minchi23;
    v.Write("overall_minchi");
    // how many points do you want?
    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";
      ///Structure defining the slice 
      struct ProfDef{
        const IFitVar * fitvar;   // which var to fix in each point?
        double minX;
        double maxX;
        std::vector <const IFitVar * > profvars; /// which vars to profile?
        std::vector <TString> profvarnames;      ///names of these vars
        std::map  <const IFitVar *, std::vector <double> >profseeds; /// seeds for osc. params
        TString shortname;
      };

       /// will loop once if you selected w/o octant split and twice if you want LO and UO curves
      for(auto const & thseed:th23seeds){  
        ProfDef profdef;
        const IFitVar* kCurTheta23 =  thseed.var; 
         
        /// deltaCP slice
        if(!th23slice && !dmsqSlice && !th13Slice){ 
          profdef.fitvar =  &kFitDeltaInPiUnits;
          profdef.minX = 0;
          profdef.maxX = 2;
          profdef.profvars =  {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23};               ///takes th23 vars from th23seeds (defined in seed controller) - it can be th23 varconstrained in LO or UO, or ordinary th23
          profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
          profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}} };  /// th23 seed from the seed controller
          profdef.shortname="delta";
        }
 
        ///th23 slice - no split into octants   
        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";
        }
        ///th13 slice
        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";
        }
        ///dmsq32 slice
        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";
        }
        // these three lines are arguments for SqrtProfile, they are used for FC
        std::map<const IFitVar*, TGraph*> profVarsMap;   
        std::map<const ISyst*, TGraph*> systMap;
        MakeMaps (profdef.profvars, profVarsMap, systs, systMap);

        
        std::cout << " Profile "<< thseed.label << "\n";
        ///make a profile       
        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;
        //save to the file
        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


///Done with fitting!
/// Now stile work

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

     //read minichi
    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

    /// loop over all slices in this ploe (IH, NH, LO, UO)
    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({new TGraph(h), hie + " " + thseed.label});
            //      graphs.push_back({HistoToCyclicGraph(h),hie + " " + thseed.label});
          }
        }
        DrawSliceCanvas(sliceName, 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, 3, 
                          (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
        }
        
        if(sliceName.Contains("dmsq")) {
          DrawSliceCanvas(sliceName, 3, 2.15, 2.75);
        }
        if(sliceName.Contains("th13")) {
          DrawSliceCanvas(sliceName, 3, 0.05, 0.2);
        }
      }

      //loop over all slices and set colors
      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("LO")) 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");
      }
      
      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 +
                    ext);
      }//end ext

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

      std::vector< const ISyst * > systs =  {};
      if (corrSysts){      // syst case
          systs =  {&kAnaCalibrationSyst}; // replace by your new syst
      }
      
      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"?k3SigmaConfidenceColorNH:kPrimColorIH);
            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