reach_2018_dCPfractions.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Progress.h"
//#include "CAFAna/Core/Utilities.h"
//#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Ana2018/nue/joint_fit_2018_tools.h"

#include "OscLib/IOscCalc.h"
#include "Utilities/rootlogon.C"

using namespace ana;

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH2.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TMath.h"
#include "TString.h"

#include "TColor.h"

#include <iostream>
#include <iomanip>
#include <vector>


////////////////////////////////////////////////////////////////////////////////
void PrintOscParameters(const osc::IOscCalcAdjustable * calc)
{
  std::vector <const IFitVar*> vars = {&kFitDeltaInPiUnits,
                                       &kFitSinSqTheta23,
                                       &kFitDmSq32Scaled,
                                       &kFitSinSq2Theta13};
  std::cout << "Oscillation parameters ";
  for (const auto & v:vars){
    std::cout << v->LatexName() << " " << v->GetValue(calc) << ", ";
  }
  std::cout <<  std::endl;
}

void PrintSeparator(){
  std::cout << "\n" << std::string(100, '*')<< std::endl;
}

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

struct SuperPred{
  bool isNue;
  std::string flux;
  double pot;
  const IPrediction* prediction;
};

std::vector <SuperPred> GetAllPredictions(const double potFHC,
                                          const double potRHC,
                                          const bool stats = false,
                                          const bool nueOnly = false,
                                          const bool numuOnly = false,
                                          const std::string nueXP = "combo",
                                          const std::string numuXP = "prop")

{
  std::vector <SuperPred> allpreds;
  if(!numuOnly){
    if(potFHC > 0) {
      allpreds.push_back({true, "fhc", potFHC,
            GetNuePrediction2018(nueXP, DefaultOscCalc(), !stats, "fhc",
                                 true, false, true, false)});
    }

    std::cout << "\n" << std::string(100, '*')<< std::endl;

    if(potRHC > 0) {
      allpreds.push_back({true, "rhc", potRHC,
            GetNuePrediction2018("prop", DefaultOscCalc(), !stats, "rhc",
                                 true, false, true, false)});
    }
  } //end add nue predictions

  std::cout << "\n" << std::string(100, '*')<< std::endl;

  const int nnumu = 4;
  if(!nueOnly){
    if(potFHC > 0){
      auto numu_preds = GetNumuPredictions2018(nnumu, !stats, "fhc");
      for (int i = 0; i < nnumu; ++ i) {
        allpreds.push_back({false,"fhc", potFHC, numu_preds[i]});
      }
    }

    std::cout << "\n" << std::string(100, '*')<< std::endl;

    if(potRHC > 0){
      auto numu_preds = GetNumuPredictions2018(nnumu, !stats, "rhc");
      for (int i = 0; i < nnumu; ++ i) {
        allpreds.push_back({false,"rhc", potRHC, numu_preds[i]});
      }
    }
  }//end add numu predictions

  return allpreds;
} //end get all predictions

////////////////////////////////////////////////////////////////////////////////
struct potCombo{
  int year;
  double potFHC;
  double potRHC;
};

std::vector <potCombo> GetPOTCombination(const TString potOption)
{
  std::vector <potCombo> ret;
  if(potOption == "default_short"){
    ret =  {
      //{2017,  9.1,  0  },
//          {2018,  9.1,  8.1},
//          {2019, 11.6, 11.6},
            {2019,12.61,12.61},
//          {2020, 14.6, 14.6}
    };
  }

  if(potOption == "nominal_outdated"){
    ret = { //{2016, 6.05,0.00},
            {2017, 8.85,0.00},
            {2018, 8.85,8.10},
            {2019, 11.48,11.48},
            {2020, 14.48,14.48},
            {2021, 17.48,17.48},
            {2022, 20.48,20.48},
            {2023, 23.48,23.48},
            {2024, 26.48,26.48}};
  }
  if(potOption == "beampower_updated"){
    ret = { //{2016, 6.05,0.00},
            {2017, 9.47,0.00},
            {2018, 9.47,6.91},
            {2019,12.61,12.61},
            {2020,16.66,16.66},
            {2021,21.22,21.22},
            {2022,25.79,25.79},
            {2023,30.85,30.85},
            {2024,35.92,35.92} }; 
    }
    
  std::cout << "\n POT combination is: " << potOption << "\n";
  for(auto const & combo:ret){
    std::cout << combo.year << ": \t"
              << combo.potFHC<< " e20 POT FHC, \t"
              << combo.potRHC<< " e20 POT RHC"
              << std::endl;
  }
  return ret;
}

TString GetPOTComboTitle(const TString potOption, bool stats = true){
  TString ret; 

  if(potOption=="nominal_outdated") 
    ret = "#splitline{DUMMY 2018 ana, stats only }{old nom POT combo } ";
  if(potOption=="beampower_updated") 
    ret = "#splitline{2018 analysis techniques and}{projected beam exposure improvements} ";
  
  return ret;
}

TLatex * PrintParamLabel (TString label){
  TLatex * ltx = new TLatex();
  ltx->SetTextAlign(13);
  ltx->SetTextSize(0.04);
  ltx->DrawLatexNDC(0.12,0.97,label.Data());
  return ltx;
}

////////////////////////////////////////////////////////////////////////////////
TCanvas * PaintReachCanvas(TString param_label, 
                           TString pot_option,
                           bool pot_label = true,
                           bool significance = false,
                           double maxY = 100,
                           TString yaxistitle = "")
{

  auto ret = new TCanvas;
  gPad->SetTopMargin(0.10);

  auto combos = GetPOTCombination(pot_option);
  int minX = combos.begin()->year;
  int maxX = (combos.end()-1)->year;
  int nBins = (combos.size() + 1) * 2;

  auto axes = new TH2D ("ax",";Year;Fraction of #delta_{CP} values (%)",
                        nBins, minX - 0.5, maxX + 0.5, 100, 0, maxY);
  if(significance) axes->GetYaxis()->SetTitle("Significance #left(#sigma=#sqrt{#Delta#chi^{2}}#right)"); 
  if(yaxistitle != "") axes->GetYaxis()->SetTitle(yaxistitle);

  axes->GetXaxis()->CenterTitle();
  axes->GetYaxis()->CenterTitle();
  axes->GetYaxis()->SetTitleOffset(0.6);
  axes->Draw();

  PrintParamLabel(param_label);

  TLatex* ltx2 = new TLatex();
  ltx2->SetTextSize(1.2*axes->GetXaxis()->GetLabelSize());
  if(pot_label)  ltx2->DrawLatexNDC(0.35, 0.16, GetPOTComboTitle(pot_option));

//  gPad->SetGridx();
//  gPad->SetGridy();

  TLine * l1 = new TLine();
  l1->SetLineStyle(3);
  if(significance)  l1->DrawLine(minX-0.5,3,maxX+0.5,3.); //3 sigma line

  SimulationSide();
  return ret;
}

double GetFractiondCP(TGraph * gr, const double level){

  //integration above level makes more sense. no time to learn
  //transform the graph into a list of sigmas
  std::vector <double>sigmas;
  int nsteps = 1000;
  for(int i=0;i<nsteps;++i){
    sigmas.push_back(gr->Eval((double) i * 2./nsteps));
  }
  int numAbove = std::count_if(sigmas.begin(),sigmas.end(),
                               [level](double s){return s>=level;});

  return (double) 100* numAbove / sigmas.size();
}

// New for 2018: Fitter hates us. Shiqi wrote this one; is it worth committing?
class FitSinSqTheta23MaxMix: 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 "ssth23MaxMix";}
  virtual std::string LatexName() const {return "sin^{2}#theta_{23} MaxMix";}
  virtual double LowLimit() const {return 0.49;}
  virtual double HighLimit() const {return 0.51;}
};
const IConstrainedFitVar *  kFitSinSqTheta23MaxMix = new FitSinSqTheta23MaxMix();

////////////////////////////////////////////////////////////////////////////////
MultiExperiment * MakeFakeExperiment(const std::vector<SuperPred> allpreds,
                                     osc::IOscCalcAdjustable* calcFake,
                                     std::vector < std::pair <TH1D*, double> > cosmics = {},
                                     std::vector <const IExperiment * > constraints = {},
                                     bool verbose =false)
{

  bool hascosmics = cosmics.size()>0;
  if(hascosmics && cosmics.size()!=allpreds.size()){
    std::cerr << "\n Size of cosmics doesn't match predictions \n";
    exit(1);
  }
  std::vector <const IExperiment * > expts;

  if(verbose) {
    std::cout << "Creating Fake experiment \n";
    PrintOscParameters(calcFake);
  }
  
  for(int i = 0; i < (int) allpreds.size(); ++i){

    auto thisdata = GetFakeData (allpreds[i].prediction, calcFake, allpreds[i].pot, 
                                 hascosmics?cosmics[i].first:0);
    if(allpreds[i].pot > 0){
      if(hascosmics){
        expts.push_back(new SingleSampleExperiment(allpreds[i].prediction, *thisdata, 
                                                   cosmics[i].first, cosmics[i].second) );
      }
      else {
        expts.push_back(new SingleSampleExperiment(allpreds[i].prediction, *thisdata) ) ;
      }
      
      if(verbose){
        auto temp = thisdata->ToTH1(thisdata->POT());
        std::cout << "Added experiment " 
                  << allpreds[i].pot << " POT "
                  << allpreds[i].flux << " "
                  << (allpreds[i].isNue?"nue ":"numu ")
                  << temp->Integral() 
                  << " events\n";
        delete temp;
      }
    }
  }

  // Add additional constraints and create expt
  if(verbose) {
    std::cout << "Returning Multiexperiment from " << expts.size() 
              << " SingleSample experiments and "
              << constraints.size() << " external constraints" << std::endl;
  }
  expts.insert(expts.end(), constraints.begin(), constraints.end());
  
  MultiExperiment * expt = new MultiExperiment (expts);
  
  //Maybe set sys correlations here? or add a new function
  return expt;
}


////////////////////////////////////////////////////////////////////////////////
void  PlotDeltaFractionVsYear(TString meas, 
                              std::vector <TString> hierarchies,
                              std::vector <TString> thetas, 
                              std::vector <TString> combo_names, 
                              std::vector <Int_t >  colors,
                              bool react = true, 
                              std::map < TString, TString > ssth23_label = {},
                              std::map < TString, TString> combo_label={}) 
{
  double limit = 3.84; //dchisq for 95%o
  
  TLegend * leg = new TLegend (0.12, 0.6, 0.5, 0.89);
  leg->SetFillStyle(0);
  leg->SetMargin(1.3*leg->GetMargin());
  if( combo_names.size()>1 ) leg->SetNColumns(2);
  leg->SetTextSize(0.04);

  TString header = "";//"#splitline{NOvA joint #nu_{e}+#nu_{#mu}}{";
  if(meas == "hie") header += "Hierarchy at 95% CL";
  if(meas == "cpv") header += "CPV at 95% CL";
//  header+="}";
  leg->SetHeader(header);

  
  for(uint i = 0; i < combo_names.size(); ++i){
    
    auto combos = GetPOTCombination (combo_names[i]);
    
    int col = 0;
    for (auto const & hie: hierarchies){      
      for (auto const & theta:thetas){
        
        TString fname = ("dCPfraction/hist_" + meas+ "_" + hie + "_" + 
                         theta + "_" + (react ? "react_" : "") + 
                         combo_names[i] + ".root") ;
        
        auto file = new TFile (fname,"read");
        if(!file) {std::cerr << "No file " << fname << std::endl; exit(1);}
        auto gr_frac = new TGraph();
        
        for(auto & combo:combos){       
          
          auto gr = (TGraph*) file->Get(meas + "_" + hie + "_" + theta + "_" + 
                                        (react ? "react_" : "") + 
                                        combo_names[i] +
                                        "_" + std::to_string(combo.year) + "/gr_dchisq");
          
          if(!gr) {std::cerr << "Problem w graph \n"; continue;}
          auto thisfrac = GetFractiondCP(gr, limit);
          gr_frac->SetPoint(gr_frac->GetN(), combo.year, thisfrac);       
        }
        gr_frac->SetLineWidth(2);
        gr_frac->SetLineStyle(i+1);
        gr_frac->Draw("L");
        gr_frac->SetLineColor(colors[col]);
        ++col;
        leg->AddEntry(gr_frac, hie + " " + ssth23_label[theta] + " " + 
                      combo_label[combo_names[i]],"l");
      }//end ssth23
    }//end hie
  }//end combo names
  leg->Draw();
  gPad->RedrawAxis();
}


void  PlotDeltaFractionVsSignificance(TString meas, 
                                      std::vector <TString> hierarchies,
                                      std::vector <TString> thetas, 
                                      std::vector <TString> combo_names, 
                                      std::vector <Int_t >  colors,
                                      bool react = true, 
                                      std::map < TString, TString > ssth23_label = {},
                                      std::map < TString, TString> combo_label={}) 
{
  gPad->SetBottomMargin(0.12);

  TString head= ""; 
  double maxx = 6.;
  if(meas == "hie") head = "Hierarchy resolution";
  if(meas == "cpv") { head = "CP violation"; maxx =3.5;} 
  auto leg2 = new TLegend(0.6,0.6,0.9,0.9);
  leg2->SetFillStyle(0);
  leg2->SetHeader(head);

  auto ax2 = new TH2D("ax2",
                      ";Significance #left(#sigma=#sqrt{#Delta#chi^{2}}#right);Fraction of #delta_{CP} values (%)",
                      100,0,maxx,100,0,100);
  ax2->Draw();
  ax2->GetXaxis()->CenterTitle();
  ax2->GetYaxis()->CenterTitle();
  ax2->GetYaxis()->SetTitleOffset(0.8);

  int col =0;
  for (auto const & hie: hierarchies){  
    for (int th = 0; th < 2; ++th){
      
      TString fname = "dCPfraction/hist_" + meas + "_" + hie +"_"+thetas[th] + "_" + 
        (react? "react_":"") + combo_names[0] + ".root";
      auto file = new TFile (fname,"read");
      if(!fname) exit(1);
      
      for(int year:{2024}){
        auto gr = (TGraph*) file->Get(meas + "_"+ hie +"_" + thetas[th] + "_" + 
                                      (react?"react_":"") + combo_names[0] +
                                      "_" + std::to_string(year) + "/gr_dchisq");
        if(!gr) {std::cerr << "Problem w graph \n"; continue;}
        
        auto gr_frac = new TGraph();
        
        for(double lim = 0; lim < 6; lim += 0.05 ){
          auto thisfrac = GetFractiondCP(gr, lim*lim);
          gr_frac->SetPoint(gr_frac->GetN(), lim,thisfrac);       
        }
        gr_frac->SetLineWidth(2);
        gr_frac->SetLineStyle(th+1);
        gr_frac->Draw("L");
        gr_frac->SetLineColor(colors[col]); 
        leg2->AddEntry (gr_frac, hie + " " + ssth23_label[thetas[th]], "l");
//      if (th==0) leg2->AddEntry(gr_frac, std::to_string(year).c_str(), "l");
//      ++col;
      }//end years
    }//end th
    ++col;
  }//end hie
/*    auto dummy = new TGraph();
    dummy->SetLineColor(kGray+1);
    dummy->SetLineWidth(2);
    leg2->AddEntry(dummy->Clone(),"sin^{2}#theta_{23}=0.6","l");
    dummy->SetLineStyle(2);
    leg2->AddEntry(dummy->Clone(),"sin^{2}#theta_{23}=0.4","l");
*/
    leg2->SetTextSize(ax2->GetXaxis()->GetLabelSize());
    leg2->Draw();
//    gPad->SetGrid();
    gPad->RedrawAxis();
}
////////////////////////////////////////////////////////////////////////////////
void  PlotSignificanceVsYear_FixDelta(TString meas, 
                                      std::vector <TString> hierarchies,
                                      std::vector <double> deltas,
                                      std::vector <TString> thetas, 
                                      std::vector <TString> combo_names, 
                                      std::vector <Int_t >  colors,
                                      bool react = true, 
                                      std::map < double, TString> delta_label = {}
  )
{
  TString head= "";
  if(meas == "hie") head = "Hierarchy resolution";
  if(meas == "cpv") head = "CP violation";
  if(meas == "oct") head = "Octant determination";
  if(meas == "mix") head = "Max. mixing rejection";

  auto leg3 = new TLegend(0.1,0.6,0.5,0.89);
  leg3->SetHeader(head);
  leg3->SetTextSize(0.04);
  
  for(uint i = 0; i < combo_names.size(); ++i){
      
    auto combos = GetPOTCombination (combo_names[i]);
    int col = 0;
    for (auto const & hie: hierarchies){   
      for(double delta:deltas){
        
        auto gr_max = new TGraph();
        auto gr_min = new TGraph();
        
        for(auto & combo:combos){
          
          double maxchi = -999.; double minchi = 999.;
          
          for (auto const & theta:thetas){
            TString fname = "dCPfraction/hist_" + meas + "_" + hie + "_" + theta + "_" + 
              (react ? "react_":"") + combo_names[i] + ".root";
          
          TFile file (fname,"read");
          if(file.IsZombie()) {std::cerr << "No file " << fname << std::endl; exit(1);}
          
          auto gr = (TGraph*) file.Get(meas + "_" + hie + "_" + theta + "_" + 
                                       (react ? "react_":"") + combo_names[i] +
                                       "_" + std::to_string(combo.year) + "/gr_dchisq");
          
          if(!gr) {std::cerr << "Problem w graph \n"; continue;}
          
            maxchi = std::max(maxchi, gr->Eval(delta)); 
            minchi = std::min(minchi, gr->Eval(delta)); 
          } //end thetas
          gr_max->SetPoint(gr_max->GetN(),combo.year, sqrt(maxchi));
          gr_min->SetPoint(gr_min->GetN(),combo.year, sqrt(minchi));
        }//end years
        auto tempN = gr_min->GetN();
        for(int i = 0; i < tempN; ++i ){
          gr_max->SetPoint(gr_max->GetN(), 
                           gr_min->GetX()[tempN - i - 1], 
                           gr_min->GetY()[tempN - i - 1]);
        }
      
        gr_max->Draw("f");
        gr_max->SetFillColorAlpha(colors[col],.33);
        ++col;
        gr_max->SetLineColor(gr_max->GetFillColor());
        if(delta_label.size()==0) {
          leg3->AddEntry(gr_max, hie + TString::Format(" #delta_{CP}= %.2f#pi",delta),"f");}
        else{
          leg3->AddEntry(gr_max, hie + " #delta_{CP}="+delta_label[delta],"f");
        }
      }//end deltas
    }//end hie
  }//end combo names
  leg3->SetFillStyle(0);
  leg3->Draw();
}



////////////////////////////////////////////////////////////////////////////////
void  PlotSignificanceVsYear_FixTh23(TString meas, 
                                     std::vector <TString> hierarchies,
                                     std::vector <TString> thetas, 
                                     std::vector <TString> combo_names, 
                                     std::vector <Int_t >  colors,
                                     bool react = true, 
                                     std::map < TString, TString > ssth23_label = {})
{
  TString head= "";
  if(meas == "hie") head = "Hierarchy resolution";
  if(meas == "cpv") head = "CP violation";
  if(meas == "oct") head = "Octant determination";
  if(meas == "mix") head = "Max. mixing rejection";

  auto leg3 = new TLegend(0.1,0.6,0.5,0.88);
  leg3->SetHeader(head);
  leg3->SetTextSize(0.04);
  
  for(uint i = 0; i < combo_names.size(); ++i){
      
    auto combos = GetPOTCombination (combo_names[i]);
    int col = 0;
    for (auto const & hie: hierarchies){
      for (auto const & theta:thetas){
        
        TString fname = "dCPfraction/hist_" + meas + "_" + hie + "_" + theta + "_" + 
          (react ? "react_":"") + combo_names[i] + ".root";
        
        TFile file (fname,"read");
        if(file.IsZombie()) {std::cerr << "No file " << fname << std::endl; exit(1);}
        
        auto gr_max = new TGraph();
        auto gr_min = new TGraph();
        
        for(auto & combo:combos){
          
          double maxchi = -999.; double minchi = 999.;
          
          auto gr = (TGraph*) file.Get(meas + "_" + hie +"_" + theta + "_" + 
                                       (react ? "react_":"") + combo_names[i] +
                                       "_" + std::to_string(combo.year) + "/gr_dchisq");
          
          if(!gr) {std::cerr << "Problem w graph \n"; continue;}
          
          minchi = TMath::MinElement(gr->GetN(),gr->GetY());      
          maxchi = TMath::MaxElement(gr->GetN(),gr->GetY());      
          
          gr_max->SetPoint(gr_max->GetN(),combo.year, sqrt(maxchi));
          gr_min->SetPoint(gr_min->GetN(),combo.year, sqrt(minchi));
        } //end years
        auto tempN = gr_min->GetN();
        for(int i = 0; i < tempN; ++i ){
          gr_max->SetPoint(gr_max->GetN(), 
                           gr_min->GetX()[tempN - i - 1], 
                           gr_min->GetY()[tempN - i - 1]);
        }
        
        gr_max->Draw("f");
        gr_max->SetFillColorAlpha(colors[col],.33);
        ++col;
        gr_max->SetLineColor(gr_max->GetFillColor());
        if(ssth23_label.size()==0) {
          leg3->AddEntry(gr_max, hie + " " + theta, "f");}
        else{
          leg3->AddEntry(gr_max, hie + " "+ssth23_label[theta],"f");
        }
      }//end theta
    }//end hier
  }//end combo names
  leg3->SetFillStyle(0);
  leg3->Draw();
}


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

void reach_2018_dCPfractions(const bool createFile = false,
                             const TString options = "",
                             const TString pot_combo = "default_short")
{
  // Shared info
  TString filename = "dCPfraction/hist_"+ options + "_" + pot_combo + ".root";
  PrintSeparator();

  auto combos = GetPOTCombination(pot_combo);
  PrintSeparator();

  //////////////////////////////////////////////////////////////////////
  bool isNH = options.Contains("NH");
  bool isIH = options.Contains("IH");
  bool isUO = options.Contains("UO");  
  bool isLO = options.Contains("LO");
  bool isMAX = options.Contains("MAX");

  assert(isNH || isIH || !createFile);
  assert(isUO || isLO || isMAX || !createFile);

  double ssth23 = (isUO ? 0.6 : 0.4);
  if(isMAX) ssth23 = 0.5;

  double dmsq32 = 2.50e-3;
  if(!isNH) dmsq32 = -2.55e-3;

  double dcp_pi = 0.2;

  bool run_hie = options.Contains("hie");
  bool run_cpv = options.Contains("cpv");
  bool run_oct = options.Contains("oct");
  bool run_mix = options.Contains("mix");

  assert(run_hie || run_cpv || run_oct || run_mix);

  bool reactor = options.Contains("react");

  if( (run_mix || run_oct) && isMAX ) {
    std::cerr << "No wrong octant or rejection of maximal mixing if MAX is true"
             << std::endl;
    exit(1);
  }

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

  if(createFile){
    
    auto file = new TFile(filename, "recreate");
    auto preds = GetAllPredictions(1E20, 1E20, true);
    PrintSeparator();
    //////////////////////////////////////////////////////////////////////
    //Setup for hierarchy,cpv, dcp fractions
    double dstep = 2./32.;

    std::vector <const IFitVar *> fitvars_hie = {&kFitDeltaInPiUnits, 
                                                 &kFitSinSqTheta23,
                                                 &kFitDmSq32
    };
    std::vector <const IFitVar *> fitvars_cpv = {
                                                 &kFitSinSqTheta23,
                                                 &kFitDmSq32 
    };
    std::vector <const IFitVar *> fitvars_oct = {&kFitDeltaInPiUnits,
                                                 isUO ? 
                                                   kFitSinSqTheta23LowerOctant : 
                                                   kFitSinSqTheta23UpperOctant,
                                                 &kFitDmSq32 
    }; 
    std::vector <const IFitVar *> fitvars_mix = {&kFitDeltaInPiUnits,
                                                 kFitSinSqTheta23MaxMix,
                                                 &kFitDmSq32
    };

    std::map <const IFitVar*, std::vector <double> > seeds_hie =
      {
        {&kFitDeltaInPiUnits, {0.,0.5,1.,1.5}},
        {&kFitSinSqTheta23,   {0.3,0.7}}
      };
    std::map <const IFitVar*, std::vector <double> > seeds_cpv =
      {
        {&kFitSinSqTheta23,   {0.3,0.7}}
      };
    std::map <const IFitVar*, std::vector <double> > seeds_oct =
      {
        {&kFitDeltaInPiUnits, {0.,0.5,1.,1.5}}
      };
    std::map <const IFitVar*, std::vector <double> > seeds_mix =
      {
        {&kFitDeltaInPiUnits, {0.,0.5,1.,1.5}},
      };
    if(isUO) seeds_mix.insert({kFitSinSqTheta23LowerOctant,{0.49}});
    if(isLO) seeds_mix.insert({kFitSinSqTheta23UpperOctant,{0.51,0.6}});

    std::vector <const IExperiment* > constraints = {};
    if(reactor) {
      for (auto fv:{&fitvars_hie, &fitvars_cpv, &fitvars_oct, &fitvars_mix}){
        fv->push_back(&kFitSinSq2Theta13);
      }
      constraints.push_back(WorldReactorConstraint2017());
    }
    else {
      std::cerr << "\nWARNING: Not using reactor constraint\n";
    }
    //////////////////////////////////////////////////////////////////////
    for (auto const & combo:combos){
      PrintSeparator();
      std::cout << "\n\n Start POT combination "
                << combo.year << ": \t"
                << combo.potFHC<< " e20 POT FHC, \t"
                << combo.potRHC<< " e20 POT RHC \n"
                << std::endl;

      //Change the pot in the super predictions
      for(auto & pred:preds){
        if(pred.flux=="fhc") pred.pot = combo.potFHC * 1E20;
        if(pred.flux=="rhc") pred.pot = combo.potRHC * 1E20;
      }
      // Create an identifier and folder based on the pot combo
      // Store the POT info just in case
      TString dirYear = ( options + "_" + pot_combo + "_" + 
                          std::to_string(combo.year) );

      auto dir = file->mkdir(dirYear);
      dir->WriteObject(&combo,"potCombo");

      // Store the fake data calculator just in case - dCP will change
      auto calc_fake = DefaultOscCalc();
      SetFakeCalc(calc_fake, ssth23, dmsq32, dcp_pi);
      ana::SaveTo(* (osc::IOscCalc*) calc_fake, dir, "calc");

      //One graph = one dCP loop. Set steps earlier
      TGraph* gr = new TGraph();

      Progress prog ("Starting dCP loop");

      for(double del = 0; del < 2; del += dstep){
        prog.SetProgress(del/2);

        SetFakeCalc(calc_fake, ssth23, dmsq32, del);

        // Create the experiment
        std::vector < std::pair <TH1D*, double > > cosmics = {};

        //Some verbosity to check the fake experiments, but not all
        bool verb = (del == 0); 
        auto thisexpt = MakeFakeExperiment(preds, calc_fake, 
                                           cosmics, constraints, verb);

        //maybe add correlations here

        //placeholder
        std::vector <const ISyst* > mySysts = {};

        double minchi_true = 0;
        double minchi_wrong = 1e20;
        SystShifts systSeed = SystShifts::Nominal();
        auto calc = DefaultOscCalc();

        if(run_hie){
          // Get fits in wrong hierarchy
          MinuitFitter fit_hie(thisexpt, fitvars_hie, mySysts);
          calc->SetDmsq32(-1. * calc_fake->GetDmsq32()); //todo
          minchi_wrong = fit_hie.Fit(calc,systSeed, seeds_hie, {},
                                     verb ? IFitter::kVerbose : IFitter::kQuiet)->EvalMetricVal();
        }
        if(run_cpv){
          MinuitFitter fit_cpv(thisexpt, fitvars_cpv, mySysts);
          for(int hie:{-1,1}) { //test both hierarchies
            for(int dcp:{0,1}) { //test cp conserving values
              systSeed.ResetToNominal();
              calc = DefaultOscCalc();
              calc->SetDmsq32(hie * calc_fake->GetDmsq32()); 
              calc->SetdCP(dcp*M_PI);
              minchi_wrong = std::min(minchi_wrong,
                                      fit_cpv.Fit(calc,systSeed, seeds_cpv, {},
                                                  verb ? IFitter::kVerbose : 
                                                  IFitter::kQuiet)->EvalMetricVal());
            }//end dcp
          }//end hie
        } //end cpv

        if(run_oct){
          MinuitFitter fit_oct(thisexpt, fitvars_oct, mySysts);
          for(int hie:{-1,1}) { //test both hierarchies
            systSeed.ResetToNominal();
            calc = DefaultOscCalc();
            calc->SetDmsq32(hie * calc_fake->GetDmsq32()); 
            calc->SetTh23(isUO? asin(sqrt(0.4)): asin(sqrt(0.6))); 
            //wrong octant is set through seeds and fitvars
            minchi_wrong = std::min(minchi_wrong,
                                    fit_oct.Fit(calc,systSeed, seeds_oct, {},
                                                verb ? IFitter::kVerbose : 
                                                IFitter::kQuiet)->EvalMetricVal());
          }//end hie
        } //end oct
        
        if(run_mix){
          MinuitFitter fit_mix(thisexpt, fitvars_mix, mySysts);
          for(int hie:{-1,1}) { //test both hierarchies
            systSeed.ResetToNominal();
            calc = DefaultOscCalc();
            calc->SetDmsq32(hie * calc_fake->GetDmsq32()); 
            calc->SetTh23(asin(sqrt(0.5)));
            minchi_wrong = std::min(minchi_wrong,
                                    fit_mix.Fit(calc,systSeed, seeds_mix, {},
                                                verb ? IFitter::kVerbose : 
                                                IFitter::kQuiet)->EvalMetricVal());
          }//end hie
        } //end mix

        // Check best fit overall sometimes; todo: some cases actually want 
        // to store, not set to 0
        if(verb && 1){
          MinuitFitter fit_best(thisexpt, fitvars_hie, mySysts);          
          for(int hie:{-1,1}) { 
            systSeed.ResetToNominal();
            calc = DefaultOscCalc();
            calc->SetDmsq32(hie * calc_fake->GetDmsq32()); 
            minchi_true = fit_best.Fit(calc,systSeed, seeds_hie)->EvalMetricVal();
            std::cout << "\nFound true for " << (hie > 0 ? "NH" : "IH")
                      << " dCP/pi" << del << " minchi " << minchi_true
                      << std::endl;
            minchi_true = 0; //bad bad hack hack
          }
        }
        gr->SetPoint(gr->GetN(),del,minchi_wrong - minchi_true);
      }//end dCP loop
      prog.Done();
      dir->cd();
      gr->Write("gr_dchisq");
    }//end combos
    std::cout << "Wrote graphs to file " << filename << std::endl;  
  }//end create file

////////////////////////////////////////////////////////////////////////////////
  else{ //make plots

    std::vector <TString> thetas = {};
    if(isUO || 1) thetas.push_back("UO");
    if(isMAX || 1) thetas.push_back("MAX");
    if(isLO || 1) thetas.push_back("LO");

    std::map < TString, TString > ssth23_label = {{"UO","sin^{2}#theta_{23}=0.6"},
                                                  {"MAX","sin^{2}#theta_{23}=0.5"},
                                                  {"LO","sin^{2}#theta_{23}=0.4"}
    };

    std::vector <TString> combo_names = {}; //sue me
    if(0) combo_names.push_back("nominal_outdated");
    if(1) combo_names.push_back("beampower_updated");

    std::map < TString, TString > combo_label = {{"beampower_updated",""}};

    std::vector <Int_t > Colors = { 
      //violet blue cyan green yellow orange red magenta
      TColor::GetColor("#6c71c4"), TColor::GetColor("#268bd2"),
      TColor::GetColor("#2aa198"), TColor::GetColor("#859900"),
      TColor::GetColor("#b58900"), TColor::GetColor("#cb4b16"),
      TColor::GetColor("#dc322f"), TColor::GetColor("#d33682")
    };
    std::vector <Int_t > Colors_alt = {  
      //cyan  violet yellow magenta blue red orange green
      TColor::GetColor("#2aa198"), 
      TColor::GetColor("#6c71c4"), 
      TColor::GetColor("#b58900"),
      TColor::GetColor("#d33682"),
      TColor::GetColor("#268bd2"), 
      TColor::GetColor("#dc322f"),
      TColor::GetColor("#cb4b16"),
      TColor::GetColor("#859900"),
    };


    ////////////////////////////////////////////////////////////////////////////////
    // dCP fraction plots
    ////////////////////////////////////////////////////////////////////////////////
    TString labelShort = TString::Format("|#Deltam^{2}_{32}|=2.5#times10^{-3}eV^{2},  sin^{2}2#theta_{13}=0.082"); //todo: not generic
    //HIE
    //Having these here gives me more control. Should add a check that canvas exist or sth
    PaintReachCanvas(labelShort,combo_names[0],false,false,100);
    PlotDeltaFractionVsYear("hie", {"NH"}, thetas, combo_names, 
                            Colors_alt,
                            reactor, ssth23_label, combo_label);
    gPad->Print("reach2018_hie_dCPfrac_year_trueNH.pdf");
    
    PaintReachCanvas(labelShort,combo_names[0],false,false,100);
    PlotDeltaFractionVsYear("hie", {"IH"}, thetas, combo_names, 
                            Colors_alt,
                            reactor, ssth23_label, combo_label);
    gPad->Print("reach2018_hie_dCPfrac_year_trueIH.pdf");

    // CPV
    PaintReachCanvas(labelShort,combo_names[0],false,false,60);
    
    PlotDeltaFractionVsYear("cpv", {"NH"}, thetas, combo_names, 
                            Colors_alt, reactor, ssth23_label, combo_label);
    gPad->Print("reach2018_cpv_dCPfrac_year_trueNH.pdf");

    PaintReachCanvas(labelShort,combo_names[0],false,false,60);
    PlotDeltaFractionVsYear("cpv", {"IH"}, thetas, combo_names, 
                            Colors_alt, reactor, ssth23_label, combo_label);
    gPad->Print("reach2018_cpv_dCPfrac_year_trueIH.pdf");


    ////////////////////////////////////////////////////////////////////////////////
    // Delta fraction vs significance for year 2024
    TString label2024 = TString::Format("#Deltam^{2}_{32}=2.5#times10^{-3}eV^{2},  sin^{2}2#theta_{13}=0.082. #kern[1.5]{ } 36#times10^{20} POT(#nu)+ 36#times10^{20} POT(#bar{#nu})");

    TLine * l1 = new TLine();
    l1->SetLineStyle(3);

    //HIE
    new TCanvas();
    PlotDeltaFractionVsSignificance("hie", {"NH","IH"}, {"UO","LO"}, combo_names, 
                                    {Colors[1],Colors[6]},
                                    reactor, ssth23_label, combo_label);
    PrintParamLabel(label2024);
    SimulationSide();
    l1->DrawLine(3,0,3,100); //3 sigma line
    gPad->Print("reach2018_hie_frac_sigma.pdf");

    //CPV
    new TCanvas();
    PlotDeltaFractionVsSignificance("cpv", {"NH","IH"}, {"UO","LO"}, combo_names, 
                                    {Colors[1],Colors[6]},
                                    reactor, ssth23_label, combo_label);
    PrintParamLabel(label2024);
    SimulationSide();
    l1->DrawLine(3,0,3,60); //3 sigma line
    gPad->Print("reach2018_cpv_frac_sigma.pdf");
    
    //////////////////////////////////////////////////////////////////////////
    // XX vs year plots
    // Fixed delta: hie and cpv
    std::vector <double> deltas = {1.5,1.,0.,0.5};
    std::map <double, TString> delta_label = {{1.5,"3#pi/2"},{0.5," #pi/2"},
                                              {1.,"#pi"},{0.,"0"} };

    TString toplabel = TString::Format("sin^{2}#theta_{23} #in #left{0.4, 0.5, 0.6#right}, |#Deltam^{2}_{32}|=2.5#times10^{-3}eV^{2},  sin^{2}2#theta_{13}=0.082");

    //HIE
    PaintReachCanvas(toplabel,combo_names[0],true,true,5);
    PlotSignificanceVsYear_FixDelta("hie", {"NH"}, deltas, thetas, combo_names,
                                    Colors_alt, reactor, delta_label);
    gPad->Print("reach2018_hie_sigma_year_trueNH.pdf");

    PaintReachCanvas(toplabel,combo_names[0],true,true,5);
    PlotSignificanceVsYear_FixDelta("hie", {"IH"}, deltas, thetas, combo_names,
                                    Colors_alt, reactor, delta_label);
    gPad->Print("reach2018_hie_sigma_year_trueIH.pdf");

    //CPV
    PaintReachCanvas(toplabel,combo_names[0],true,true,5);
    PlotSignificanceVsYear_FixDelta("cpv", {"NH"}, {1.5,0.5}, thetas, combo_names,
                                    Colors_alt, reactor, delta_label);
    SimulationSide();
    gPad->Print("reach2018_cpv_sigma_year_trueNH.pdf");

    PaintReachCanvas(toplabel,combo_names[0],true,true,5);
    PlotSignificanceVsYear_FixDelta("cpv", {"IH"}, {1.5,0.5}, thetas, combo_names,
                                    Colors_alt, reactor, delta_label);
    SimulationSide();
    gPad->Print("reach2018_cpv_sigma_year_trueIH.pdf");
    
    // Fixed ssth23: oct and maxmix
    TString toplabel2 = TString::Format("#delta_{CP} #in #left[0,2#pi#right), |#Deltam^{2}_{32}=2.5#times10^{-3}eV^{2}|,  sin^{2}2#theta_{13}=0.082");
    
    //OCT
    PaintReachCanvas(toplabel2,combo_names[0],true,true,5);
    PlotSignificanceVsYear_FixTh23("oct", {"NH","IH"}, {"UO","LO"}, combo_names,
                                    Colors_alt, reactor, ssth23_label);
    
    gPad->Print("reach2018_oct_sigma_year.pdf");

    //MIX
    PaintReachCanvas(toplabel2,combo_names[0],true,true,7);
    PlotSignificanceVsYear_FixTh23("mix", {"NH","IH"}, {"UO","LO"}, combo_names,
                                   Colors_alt, reactor, ssth23_label);
    SimulationSide();
    gPad->Print("reach2018_mix_sigma_year.pdf");
  }//end make plots
} //end all