joint_fit_2018_tools.h

Go to the documentation of this file.

#pragma once

// TODO: find a better home for this header file; functions are duplicated

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
//#include "CAFAna/Prediction/PredictionSystJoint2018.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "3FlavorAna/Prediction/PredictionAddRock.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Analysis/CountingExperimentStatistics.h"
#include "3FlavorAna/Systs/JointAna2018Systs.h"
#include "CAFAna/Systs/XSecSysts.h"
#include "CAFAna/Vars/FitVars.h"

#include "CAFAna/Analysis/Style.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"

#include "Utilities/func/MathUtil.h"

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

using namespace ana;

#include "OscLib/IOscCalc.h"

#include "TFile.h"

TFile* GetGaussianSurface(TString par, TString hier)
{
  TString subdir = par.Contains("ssth23dmsq32")? "th23_dmsq" : "delta_th23";
  TString gaussDir = "/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/contours/"+subdir+"/syst/";
  TString gaussName = "hist_contours_2018_joint_realData_both_only";
  gaussName += hier;
  if(par.Contains("ssth23dmsq32"))
    gaussName += "combo_dmsq_systs.root";
  else
    gaussName += "combo_systs.root";
  

  std::cout << "Opening " << gaussDir + gaussName << std::endl;
  TFile* ret = TFile::Open(gaussDir + gaussName);
  return ret;
}


void CylindricalSmoothing(TH2* h)
{
  TH2* temp = (TH2*)h->Clone();

  int NX = h->GetXaxis()->GetNbins();
  int NY = h->GetYaxis()->GetNbins();
  double minX = h->GetXaxis()->GetBinLowEdge(1);
  double maxX = h->GetXaxis()->GetBinUpEdge(NX);
  double widthX = h->GetXaxis()->GetBinWidth(1);
  double minY = h->GetYaxis()->GetBinLowEdge(1);
  double maxY = h->GetYaxis()->GetBinUpEdge(NY);
  double widthY = h->GetYaxis()->GetBinWidth(1);

  //  assert(NX % 2 == 0 && 
  //     "Histogram must have an even number of X bins for this to work");

  // Save the location we expect the end column of the original
  // TH2 to be after we do the folding
  int foldedBin = 1 + NX/2;
 
  // Fold the initial histogram to join the ends of the X axis
  for(int i = 1; i <= NX/2; i++) {
    for(int j = 1; j <= NY/2; j++) {
      // x1, y1 denotes a position in the left-hand side half plane
      // x2, y2 denotes a position in the right-hand side half plane
      int x1 = i+NX/2;
      int x2 = i%(NX/2+1)+2;
      int y1 = j;
      int y2 = j;
      double z1 = h->GetBinContent(x1, y1);
      double z2 = h->GetBinContent(x2, y2);
      temp->SetBinContent(x1, y1, z2);
      temp->SetBinContent(x2, y2, z1);
    }
  }
  // smooth both histograms
  temp->Smooth();
  h->Smooth();
  
  // Replace end columns
  for(int y = 1; y <= NY; y++) {
    double sZ1 = temp->GetBinContent(foldedBin, y);
    double sZ2 = temp->GetBinContent(foldedBin, y);
    h->SetBinContent(1, y, sZ1);
    h->SetBinContent(NX, y, sZ2);
  }
}

std::string GetNuePredPath(std::string beam, bool isFakeData, bool NERSC) {
  std::string ret;
  bool isFHC = true;
  if(beam == "rhc") isFHC = false;

  // default path name for Nue Preds constructor
  // to point to official files
  if(!NERSC) return "";

  // NERSC paths
  else {
    std::string dir = std::string(getenv("FCHELPERANA2018_LIB_PATH"))
      + "/Predictions/provisional_singles";
    if(isFHC && !isFakeData)
      ret = dir + "/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_FHCAllSyst_nueconcat_realNDData_2018-05-22.root";
    else if(!isFHC && !isFakeData)
      ret = dir + "/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_RHCAllSyst_nueconcat_realNDData_2018-05-22.root";
    else if(isFHC && isFakeData)
      ret = dir + "/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_FHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
    else if(!isFHC & isFakeData)
      ret = dir + "/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_RHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
  }
  return ret;
}

const IPrediction * GetNuePrediction2018(std::string decomp,
                                         osc::IOscCalc *calc,
                                         bool corrSysts, 
                                         std::string beam, 
                                         bool isFakeData,
                                         bool GetFromUPS=false,
                                         bool minimizeMemory = false,
                                         bool NERSC = false){

  ENu2018ExtrapType extrap = kProportional;
  if(decomp == "noextrap") extrap = kNoExtrap;
  if(decomp == "prop") extrap = kProportional;
  if(decomp == "combo") extrap = kCombo;
  if(decomp == "fakexp") extrap = kFake;//Hack 

  PredictionSystJoint2018 * pred_fid;
  std::string nue_pred_path = GetNuePredPath(beam, isFakeData, NERSC);
  // Interactive running
  if (!GetFromUPS){
    if(corrSysts){
      std::cout << "Getting Nue Prediction from " + nue_pred_path << endl;
      if(beam=="rhc") pred_fid = new PredictionSystJoint2018 (extrap, calc, 
                                                              getAllAna2018Systs( kNueAna2018, true, kRHC), 
                                                              beam, isFakeData, true, 
                                                              nue_pred_path, minimizeMemory);
      else if(beam=="fhc") pred_fid = new PredictionSystJoint2018 (extrap, calc,
                                                              getAllAna2018Systs( kNueAna2018, true, kFHC),
                                                              beam, isFakeData, true,
                                                              nue_pred_path, minimizeMemory);
      else{ std::cout << "Beam mode needs to be either fhc or rhc" << endl; exit(1);}
    }
    else pred_fid = new PredictionSystJoint2018 (extrap, calc,{}, beam, isFakeData,
                                                 true, nue_pred_path, minimizeMemory);
  }
  // Grid running for FC
  else{
    std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));                       // need an update
    std::string name = "/pred_nue_reduced_v9.root";                                              //need an update
    if(corrSysts){
      if(beam=="rhc") pred_fid = new PredictionSystJoint2018(extrap, calc,  
                                                             getAllAna2018Systs( kNueAna2018, true, kRHC),
                                                             beam, isFakeData, true,
                                                             upsname+name, minimizeMemory);
      else if(beam=="fhc") pred_fid = new PredictionSystJoint2018(extrap, calc,
                                                                  getAllAna2018Systs( kNueAna2018, true, kFHC),
                                                                  beam, isFakeData, true,
                                                                  upsname+name, minimizeMemory);
      else{ std::cout << "Beam mode should be either fhc or rhc" << endl; exit(1);}
    }
    else
      pred_fid = new PredictionSystJoint2018(extrap, calc, {},beam, isFakeData,
                                             true, upsname+name, minimizeMemory);
  }
  


  double pot;
  if(beam=="rhc") pot = kAna2018RHCPOT;
  else pot = kAna2018FHCPOT;
  std::cout << "Scale to "<<pot<<std::endl;

  std::cout << "Nue prediction (no rock)"
            << pred_fid->Predict(calc).Integral(pot) <<std::endl;
  // Load rock

  double nonsScale = 1./10.45;
  double swapScale = 1./12.91;
  //fix to the normal path later
  std::string dirName = "/nova/ana/nu_e_ana/Ana2018/Predictions/FDRock/";
  if(NERSC) dirName = std::string(getenv("FCHELPERANA2018_LIB_PATH"))+"/Predictions/FDRock/";
  std::string fName2 = "fdrock_nue2018.root";
  std::string rName;

  if(beam=="rhc") rName = "pred_FDRock_merg_rhc";
  else rName = "pred_FDRock_merg_fhc";
  if( decomp == "noextrap"){
    if(beam=="rhc") rName = "pred_FDRock_rhc";
    else rName = "pred_FDRock_fhc";
  }
  // Point to the rock file in the ups product for FC'ing
  if (GetFromUPS){
    dirName = std::string(getenv("FCHELPERANA2017_LIB_PATH"));                                  //need an update
    fName2  = "/fdrock_nue2017.root";                                                           //need an update
    rName = "pred_FDRock";
  }

  auto rock = LoadFromFile<PredictionNoExtrap>(dirName + fName2, rName).release();

  std::cout << "Nue prediction (rock only - not scaled) "
            << rock->Predict(calc).Integral(pot) <<std::endl;
  std::cout << "Adding Rock" << endl;
  auto nue_pred = new PredictionAddRock (pred_fid, rock, nonsScale, swapScale);

  std::cout << "Nue prediction (rock added) "
            <<  nue_pred->Predict(calc).Integral(pot) << std::endl;
  std::cout << "Pure rock events:  "
            <<  nue_pred->Predict(calc).Integral(pot) - pred_fid->Predict(calc).Integral(pot) << std::endl;
  
  
  return nue_pred;
}

//////////////////////////////////////////////////////////////////////

std::pair <TH1D*, double > GetNueCosmics2018 (std::string beam, bool GetFromUPS=false, bool NERSC=false){
  Spectrum* outtime;
  if (GetFromUPS){
    std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));            //update
    std::string fname = "/nue2017_cospred.root";                                      //update
    std::string name  = "cos";       
    outtime = LoadFromFile<Spectrum>(upsname+fname,name).release();
  }
  else{
    std::string name;
    std::string anaDir = "/nova/ana/nu_e_ana/Ana2018/";
    if(NERSC) anaDir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
    if(beam=="rhc") name = "cosmic_spect_rhc";
    else name="cosmic_spect_fhc";
    outtime = LoadFromFile<Spectrum>
      (anaDir+"/Predictions/cosmic/cosmic_prediction_real_data.root", name ).release();
  }

  double livetime;
  if(beam=="rhc") livetime = kAna2018RHCLivetime ;
  else livetime = kAna2018FHCLivetime ;

  //double outN = outtime->Integral(outtime->Livetime(), 0, kLivetime);
  double outN = outtime->Integral(livetime, 0, kLivetime);

  std::cout << "\nAdding " << outtime->Integral(livetime, 0, kLivetime)
            << " cosmics from out-of-time." << std::endl;
  std::cout << "Which scale to " << outtime->Integral(livetime, 0,
                                                      kLivetime)
            << " in the spill window." << std::endl;

  std::cout << "uncertainty " << 1/sqrt(outN) << std::endl;

  TH1 *h = outtime->ToTH1(livetime, kBlack, kSolid, kLivetime);
//  assert( h->Integral() > 4 && "Nue cosmic prediction < 4. Something is wrong");    //RHC has ~0.6 cosmic events

  return {outtime->ToTH1(livetime, kBlack, kSolid, kLivetime), 1/sqrt(outN)};
}

//////////////////////////////////////////////////////////////////////

Spectrum* GetNueData2018 (std::string beam){
  Spectrum* intime;
  
  if(beam=="fhc"){
    intime = LoadFromFile<Spectrum>
      ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_fhc").release();
  }
  else if(beam=="rhc"){
    intime = LoadFromFile<Spectrum>
      ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_rhc").release();
  }
  else{ std::cout << "Beam mode != fhc or rhc" << std::endl; exit(1);}
  double POT, livetime;
  if(beam == "rhc") {livetime = kAna2018RHCLivetime; POT = kAna2018RHCPOT;}
  if(beam == "fhc") {livetime = kAna2018FHCLivetime; POT = kAna2018FHCPOT;}

  std::cout << "Loaded data spectrum:\n"
            << "POT " << intime->POT()
            << " (expected " << POT  << ")\n"
            << "Livetime " << intime->Livetime()
            << " (expected " << livetime << ")\n"
            << "Integral " << intime->Integral(intime->POT())
            << std::endl;

  return intime;
}

//////////////////////////////////////////////////////////////////////

std::vector <const IPrediction*> GetNumuPredictions2018(const int nq = 4,
                                                        bool useSysts = true,
                                                        std::string beam="fhc",
                                                        bool GetFromUPS=false,
                                                        ENu2018ExtrapType numuExtrap = kNuMu, 
                                                        bool minimizeMemory = false,
                                                        bool NERSC = false)
{
  auto calc = DefaultOscCalc();

  std::vector <const IPrediction*> numu_preds;
 
  std::string filename;
  std::string anaDir = "/nova/ana/nu_mu_ana/Ana2018/";
  if(NERSC) anaDir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
  if(beam=="rhc") filename=anaDir + "/Predictions/pred_numuconcat_rhc__numu2018.root";
  else filename= anaDir + "/Predictions/pred_numuconcat_fhc__numu2018.root";

  std::cout << "\nLoading Numu Predictions\n" ;
  for (int quant = 0; quant < nq; ++quant){
    if (GetFromUPS){
      std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));     //update
      std::string fname = "/pred_numu_reduced_v5.root";                          //update
      if(!useSysts) {      
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap,
                                                         calc, beam,
                                                         quant+1,{},
                                                         upsname+fname,
                                                         minimizeMemory));         //update
      }
      else{
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam,        
                                                         quant+1,
                                                         getAllAna2018Systs( kNumuAna2018, true),
                                                         upsname+fname,
                                                         minimizeMemory));          //update 
      }
    }
    else{
      if(!useSysts) {
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1, 
                                                         {}, filename, minimizeMemory));
      }
      else {
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1,
                                                         getAllAna2018Systs( kNumuAna2018, true), 
                                                         filename, minimizeMemory));
      }
    }
  }
  return numu_preds;
}
//////////////////////////////////////////////////////////////////////

std::vector <std::pair <TH1D*, double> > GetNumuCosmics2018(const int nq = 4, std::string beam="fhc",
                                                            bool GetFromUPS=false, bool NERSC = false)
{
   cout<<"Beam is "<<beam<<endl;
  //cosmics
  //the versionn in cafana has float not double and it's a pain
  TFile *fcosm;
  if (GetFromUPS){
    std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));          //update
    fcosm = TFile::Open((upsname+"/cos_bkg_hists_v2.root").c_str());                //update
  }
  else{
    std::string dir = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
    if(NERSC) dir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
    std::string filename;
    if (beam=="rhc") filename="cosmics_rhc__numu2018.root";
    else filename="cosmics_fhc__numu2018.root";
    fcosm = TFile::Open((dir+filename).c_str());
  }

  if(fcosm->IsZombie()) {std::cerr<< "bad cosmics\n"; exit(1);}
  std::vector <std::pair <TH1D*, double > > numu_cosmics;

  for (int i = 1; i <=nq; ++i) {
    auto h = (TH1D*) fcosm->Get((TString)"cosmics_q"+std::to_string(i))
      ->Clone(UniqueName().c_str());
    //auto hs = (TH1D*) fcosm->Get((TString)"q"+std::to_string(i)+"allsyst");        //update. no such stuff in numu files
    //numu_cosmics.push_back({h,hs->GetBinContent(1)});
    numu_cosmics.push_back({h,1/sqrt(h->Integral())}); 
    std::cout << "Cosmics all periods quantile " << i << " "
              << h->Integral()
              << " scale error "
              //<< hs->GetBinContent(1)
              << 1/sqrt(h->Integral())
              << "\n";
  }//end quantiles
  return numu_cosmics;
}
//////////////////////////////////////////////////////////////////////

std::vector <Spectrum * > GetNumuData2018(const int nq = 4, std::string beam = "fhc"){
 
  std::string file;
  if(beam == "fhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_fhc_full__numu2018.root";}
  if(beam == "rhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_rhc_full__numu2018.root";}

  std::vector <Spectrum *> numu_data;

  for (int i = 1; i <=nq; ++i) {

    auto f= new TFile(file.c_str());
    if(f->IsZombie()) {std::cerr<< "bad data\n"; exit(1);}
    numu_data.push_back(Spectrum::LoadFrom(f, ("fd_data_q"+std::to_string(i)).c_str())).release();

  }//end quantiles
  return numu_data;
}

//////////////////////////////////////////////////////////////////////

Spectrum * GetFakeData(const IPrediction* pred, osc::IOscCalc * calc,
                       const double pot, TH1D* cosmics = 0)
{
  auto temp = pred->Predict(calc).FakeData(pot);
  if(cosmics) temp += Spectrum(cosmics, pot,0);
  return new Spectrum(temp);
}

Spectrum * GetMockData(const IPrediction* pred, osc::IOscCalc * calc,
                       const double pot, TH1D* cosmics = 0)
{
  auto temp = pred->Predict(calc).MockData(pot);
  if(cosmics) temp += Spectrum(cosmics, pot,0);
  return new Spectrum(temp);
}

//////////////////////////////////////////////////////////////////////

void SetFakeCalc (osc::IOscCalcAdjustable * calc, double ssth23 = -5,
                  double dmsq32 = -5, double dCP_Pi = -5)
{
  if(dCP_Pi > -1) kFitDeltaInPiUnits.SetValue(calc, dCP_Pi);
  if(ssth23 > -1) kFitSinSqTheta23.SetValue(calc, ssth23);
  if(dmsq32 > -1) kFitDmSq32.SetValue(calc, dmsq32);
}

//////////////////////////////////////////////////////////////////////

std::vector<const ISyst*> GetJointFitSystematicList(bool corrSysts, bool nueExclusive = false, bool numuExclusive=false, bool isFHC = true, bool isRHC = true , bool intersection = true){
  std::vector<const ISyst*> ret ={} ;
  if(corrSysts){
    if(! intersection){
      if(nueExclusive){
        
         ret.push_back(&kRadCorrNue);
         ret.push_back(&kRadCorrNuebar);
         ret.push_back(&k2ndClassCurrs);
         AddJointAna2018NotCorrelSysts(ret, kNueAna2018, kBoth);
         AddNonLoadable2018Systs(ret, kNueAna2018);

      }
      if(numuExclusive){
        
         AddJointAna2018NotCorrelSysts(ret, kNumuAna2018, kBoth);
         if(isFHC) AddJointAna2018NotCorrelSysts(ret, kNueAna2018, kRHC);
         if(isRHC) AddJointAna2018NotCorrelSysts(ret, kNueAna2018, kFHC);
         if(isFHC) ret.push_back(&kAna2018NormRHC);
         if(isRHC) ret.push_back(&kAna2018NormFHC);
         AddNonLoadable2018Systs(ret, kNumuAna2018);

      }
    }
    else{
      if(nueExclusive) ret = getAllAna2018Systs(kNueAna2018);
      if(numuExclusive) ret = getAllAna2018Systs(kNumuAna2018);
      if(!nueExclusive && !numuExclusive) ret = getAllAna2018Systs(kJointAna2018, true, kBoth, true);
      //std::cout <<"I'm in else" <<std::endl;
      //if(!nueExclusive && !numuExclusive) ret = getAllAna2018Systs(kJointAna2018, true, true);
    }
  }
  std::cout << "Return list of  systematics: " ;
  for(auto & s:ret) std::cout << s->ShortName() << ", ";
  std::cout << std::endl;
  return ret;
}
std::vector<std::pair<const ISyst*,const ISyst*> > GetCorrelations (bool isNue, bool isFHC){
  std::vector<std::pair<const ISyst*,const ISyst*> > ret ={};
  std::cout<<(isNue?"not nue ":"not numu ")<<(isFHC&&!isNue?" ":(isFHC?" FHC ":" RHC "));
  //to first approximation if they are not the same they are not correlated.
  //might have to be more complicated later
  auto badsyst = GetJointFitSystematicList(true, !isNue, isNue, isFHC, !isFHC, false);
  for (auto & s:badsyst) ret.push_back ({s,NULL});

  return ret;
}

// For slices plots, plot by true octant, and make separate plots of each
class FitSinSqTheta23LowerOctant: public IConstrainedFitVar
{
public:
  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
    return util::sqr(sin(osc->GetTh23()));};
  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
    osc->SetTh23(asin(sqrt(Clamp(val))));};
  virtual std::string ShortName() const {return "ssth23LO";}
  virtual std::string LatexName() const {return "sin^{2}#theta_{23} LO";}
  virtual double LowLimit() const {return 0;}
  virtual double HighLimit() const {return 0.5;}
};

const IConstrainedFitVar * kFitSinSqTheta23LowerOctant = new FitSinSqTheta23LowerOctant();

class FitSinSqTheta23UpperOctant: public IConstrainedFitVar
{
public:
  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
    return util::sqr(sin(osc->GetTh23()));};
  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
    osc->SetTh23(asin(sqrt(Clamp(val))));};
  virtual std::string ShortName() const {return "ssth23UO";}
  virtual std::string LatexName() const {return "sin^{2}#theta_{23} UO";}
  virtual double LowLimit() const {return 0.5;} 
  virtual double HighLimit() const {return 1;}
};
const IConstrainedFitVar *  kFitSinSqTheta23UpperOctant = new FitSinSqTheta23UpperOctant();


////////////////////////////////////////////////////////////////////////////
//////////////////////// Slice, Contours Style /////////////////////////////
///////////////////////////////////////////////////////////////////////////
//move functions from plotting script here

void DrawSliceCanvas(TString slicename, const double ymax,
                     const double xmin = 0,  const double xmax = 2., bool overlay = false){
  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}})";
  }
  if(slicename.Contains("th13")) {
    title = ";sin^{2}2#theta_{13};Significance (#sqrt{#Delta#chi^{2}})";
  }
  auto axes = new TH2F("",title, 100, xmin, xmax, 100, 0, ymax);
  if(slicename.Contains("fccorr") || overlay) axes->GetYaxis()->SetTitle("Significance (#sigma)");
  // CenterTitles(axes);
  axes->GetXaxis()->CenterTitle();
  axes->GetYaxis()->CenterTitle();
  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, double low=0.7, bool isDmsq = false, bool overlay = false)
{
  TLegend * leg;
  if(isDelta) leg = new TLegend(0.6, 0.6, 0.9,0.89);
  else {
    leg = new TLegend(0.18, 0.2, 0.4, 0.49);
  }
  leg->SetTextSize(0.8*kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  for (auto &gr:graphs){
       TString str="";
    if(gr.second.Contains("O")){
      if(isDelta) leg->SetX1(0.57);
      if(isDmsq) leg->SetX1(0.18);
      if(!isDelta && !isDmsq) {leg->SetX1(0.6); leg->SetX2(0.89);}
      leg->SetMargin(0.18);
      str += gr.second.Contains("NH") ? "NH":"IH";
      str +=  gr.second.Contains("LO") ? " Lower": " Upper";
      str += " octant";
    }
    else{
      str="#splitline{";
      str += gr.second.Contains("NH") ? "Normal":"Inverted";
      str += "}{hierarchy}";    
      if(overlay) {
        str = "Uncorr. ";
        str += gr.second.Contains("NH") ? "NH" : "IH";
      }
    }
    leg->AddEntry(gr.first, str, "l");
  }
  return leg;
}


void  DrawContourCanvas (TString surfName, double minx = 0,
                         double maxx = 2, double miny = 0, double maxy = 1){
  auto c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetFillStyle(0);
  c1->SetLeftMargin(0.14);
  c1->SetBottomMargin(0.15);
  TH2* axes = new TH2F();
  TString title;
  if(surfName.Contains("delta_th23")) title=";#delta_{CP};sin^{2}#theta_{23}";
  if(surfName.Contains("th13_delta")) title=";sin^{2}2#theta_{13};#delta_{CP}/#pi";
  if(surfName.Contains("th23_dm32"))  title=";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
  axes = new TH2F("",title,100,minx,maxx,100,miny,maxy);
//  CenterTitles(axes);
  axes->GetXaxis()->CenterTitle();
  axes->GetYaxis()->CenterTitle(); 
  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.8);
  TGaxis::SetMaxDigits(3);
  if(surfName.Contains("th23_dm32")) {
    axes->GetYaxis()->SetDecimals() ;
    axes->GetYaxis()->SetTitleOffset(0.85);
    axes->GetYaxis()->SetLabelOffset(0.002);
    axes->GetYaxis()->SetLabelSize(0.058);
    axes->GetYaxis()->SetTitleSize(0.078);
    axes->GetXaxis()->SetTitleSize(.95*kBlessedTitleFontSize);
    axes->GetXaxis()->SetTitleOffset(0.86);
  }
  if(surfName.Contains("delta_th23"))  XAxisDeltaCPLabels(axes);
  c1->RedrawAxis();
}


TLegend * ContourLegend(int hie, bool fillContour, bool fccorr,
                        Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor,
                        bool bestFit){
  TLegend * leg = new TLegend(0.16,fccorr?0.19:0.17,0.65,0.28);
  leg->SetTextSize(kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  leg->SetMargin(1.3*leg->GetMargin());
  leg->SetNColumns(3);

  TGraph * dummy = new TGraph;
  dummy->SetLineWidth(3);

  if(fillContour){
    dummy->SetLineColor(kDarkColor);
    dummy->SetFillColor(kFillColor1);
    leg->AddEntry(dummy->Clone(),"1#sigma","f");
    dummy->SetFillColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","f");
    dummy->SetFillColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","f");
  }
  else{
    leg->SetMargin(1.4*leg->GetMargin());
    dummy->SetLineColor(kFillColor1);
    dummy->SetLineStyle(kSolid);
    leg->AddEntry(dummy->Clone(),"1#sigma","l");
    dummy->SetLineStyle(k2SigmaConfidenceStyle);
    dummy->SetLineColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","l");
    dummy->SetLineStyle(kSolid);
    dummy->SetLineColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","l");
  }
  if(bestFit){
    leg->SetNColumns(4);
    //leg->SetTextSize(0.8*kBlessedLabelFontSize);
//  leg->SetColumnSeparation(0.08);
    leg->SetX2(0.75);
    dummy->SetMarkerStyle(43);
    dummy->SetMarkerSize(2);
    dummy->SetMarkerColor(kDarkColor);
    dummy->SetLineColor(kDarkColor);
    dummy->SetLineStyle(7);
    if(hie>0) leg->AddEntry(dummy->Clone(),"Best Fit","p");
    if(hie<0) leg->AddEntry((TObject*)0, "#color[0]{Best Fit}", "");
    dummy->SetLineStyle(1);
  }
  if(!fccorr) leg->SetHeader("No Feldman-Cousins");
  else leg->SetHeader("Feldman-Cousins");
  return leg;
}


TLegend * FCOverlayContourLegend(int hie, bool fillContour,
                                 Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor)
{
  TLegend * leg = new TLegend(0.16,hie>0?.78:0.77,0.65,0.86);
  leg->SetTextSize(kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  leg->SetMargin(1.3*leg->GetMargin());
  leg->SetNColumns(3);

  TGraph * dummy = new TGraph;
  dummy->SetLineWidth(3);

  if(fillContour){
    dummy->SetLineColor(kDarkColor);
    dummy->SetFillColor(kFillColor1);
    leg->AddEntry(dummy->Clone(),"1#sigma","f");
    dummy->SetFillColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","f");
    dummy->SetFillColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","f");
  }
  else{
    leg->SetMargin(1.4*leg->GetMargin());
    dummy->SetLineColor(kFillColor1);
    dummy->SetLineStyle(kSolid);
    leg->AddEntry(dummy->Clone(),"1#sigma","l");
    dummy->SetLineStyle(k2SigmaConfidenceStyle);
    dummy->SetLineColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","l");
    dummy->SetLineStyle(kSolid);
    dummy->SetLineColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","l");
  }

  leg->SetHeader("No Feldman-Cousins");
  return leg;
}