sensitivity2020.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Registry.h"
#include "CAFAna/Experiment/SolarConstraints.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/FC/FCSurface.h"
#include "CAFAna/Fit/IFitter.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "3FlavorAna/Systs/NueAcceptSysts.h"
#include "CAFAna/Systs/Systs.h"
#include "3FlavorAna/Systs/WrongSignNueSyst.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "3FlavorAna/Prediction/PredictionAddRock.h"
#include "CAFAna/Vars/FitVars.h"

#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"
#include "OscLib/IOscCalc.h"
#include "Utilities/rootlogon.C"

#include "TFile.h"
#include "TStyle.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TSystem.h"
#include "TH2.h"

using namespace ana;

// MakeFile: Do you need to make the intermediate files, or just load from them?
// MakeSurf: true false
// whatsurf: "th23dcp", "dmsqth23"
// MakeSens: true false
// whatplot: "mh", "dcp", "maxmix", "oct"
// isSyst  : use syst, not only nominals
// NuMuOnly: true false --> Run over only numu, or include nue too?
// horn    : "both", "fhc", "rhc"
// isNH    : true false
void sensitivity2020(bool MakeFile        = false,
                     bool MakeSurf        = false,
                     std::string whatsurf = "th23dcp",
                     bool MakeSens        = false, 
                     bool Varyth23        = false,
                     std::string whatplot = "mh",
                     bool isSyst          = false,
                     bool isPtExtrap      = false,
                     bool NuMuOnly        = false,
                     std::string horn     = "both",
                     bool isNH            = true,
                     std::string option   = "2019BestFit",
                     bool OnGrid          = false,
                     bool isFutureSens    = false ) {

  // --- Set some constants that I'll want to use.
  const int nnumu = 4;
  const std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
  const std::vector<double> th23_seeds  = {0.45, 0.55};
  ////////////////////////////////////////////////////////////////
  /////////////////////// Load nue and numu stuff ///////////////
  ///////////////////////////////////////////////////////////////

  std::string SystTag = isSyst    ?"with_syst":"just_stat";
  std::string sNuMu   = NuMuOnly  ?"numuonly" :"numu_nue";
  std::string sExtrap = isPtExtrap?"PtExtrap" :"nonPtExtrap";
  std::string sHier   = isNH      ?"NH"       :"IH";
  std::string TagName = sNuMu+"_"+horn+"_"+sHier+"_"+SystTag+"_"+sExtrap+"_"+option;
  std::string outfilename = "run_sensitivity_"+TagName;

  std::cout << "\n Greetings! You requested  to " << (MakeFile?" make files." : " only plot." ) 
            << "\nIt will use " << horn << " files and " << SystTag << " for " << sNuMu << ", plotting: "
            << (MakeSurf?"contours ":"no contours ")
            << (isNH    ?" NH "     :" IH "        ) << " ; " << whatsurf
            << (MakeSens?" slices " :" no slices"  )
            <<" in case of slices it will be "<< whatplot 
            << "\nAll using option "          << option 
            << "\n\t TagName is thus " << TagName
            << std::endl;
  
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  // Some options for different best fits.
  double BASE_dcp    = 0.82*M_PI;
  double BASE_th23   = 0.568;
  double BASE_dmsq32 = 2.41e-3;
  std::cout << "\n Default calc is set to 2020 BF." << std::endl;
  // --- 2019 best fit
  if ( option.find("2019BestFit") != std::string::npos ) {
    BASE_dcp    = 2*M_PI;
    BASE_th23   = 0.565;
    BASE_dmsq32 = 2.48e-3;
    std::cout << "\t Updated calc to 2019 BF." << std::endl;
  }
  // --- 2019 Best Fit - NH, LO
  if ( option.find("NHLO") != std::string::npos ) {
    BASE_dcp    = 1.9*M_PI;
    BASE_th23   = 0.480;
    BASE_dmsq32 = 2.49e-3;
    std::cout << "\t Updated calc to NH LO 2019 BF." << std::endl;
  } 
  // --- 2019 Best Fit - IH UO
  if ( option.find("IHUO") != std::string::npos ) {
    BASE_dcp    =  1.5*M_PI;
    BASE_th23   =  0.560;
    BASE_dmsq32 = -2.54e-3;
    std::cout << "\t Updated calc to IH UO 2019 BF." << std::endl;
  }
  calc  ->SetdCP   ( BASE_dcp    );
  calc  ->SetTh23  ( asin(sqrt(BASE_th23)) );
  calc  ->SetDmsq32( BASE_dmsq32 );
  // --- 2018 best fit
  //calc->SetdCP(1.21*M_PI);
  //calc->SetTh23(asin(sqrt(0.558)));
  //calc->SetDmsq32(2.44e-3);

  std::vector<double> pot;
  std::vector<double> livetime;

  double temp_potFD_rhc = kAna2020RHCPOT, temp_cosFD_rhc = kAna2020RHCLivetime;
  double temp_potFD_fhc = kAna2020FHCPOT, temp_cosFD_fhc = kAna2020FHCLivetime;
  if (isFutureSens) {
    temp_potFD_rhc = kFutureRHCPOT; temp_cosFD_rhc = kFutureRHCLivetime;
    temp_potFD_fhc = kFutureFHCPOT; temp_cosFD_fhc = kFutureFHCLivetime;
  }
  const double potFD_rhc = temp_potFD_rhc;
  const double potFD_fhc = temp_potFD_fhc;
  const double cosFD_rhc = temp_cosFD_rhc;
  const double cosFD_fhc = temp_cosFD_fhc;

  std::vector<const IPrediction*> predsvector;
  std::vector<std::pair <Spectrum*, double>> cosmicsvector;

  std::vector<const IPrediction*> predictionsNue;
  std::vector<const IPrediction*> predictionsNumu;
  std::vector<std::pair <Spectrum*, double>> cosmicsNue;
  std::vector<std::pair <Spectrum*, double>> cosmicsNumu;

  string decomp_fhc = "combo"; //should be combo
  string decomp_rhc = "prop";
  if (!isPtExtrap) {
    decomp_fhc  += "_noPt";
    decomp_rhc  += "_noPt";
    outfilename += "_noPtExtrap";
  }

  bool FakeNDData = false;

  // --- Load the Nu_e things.
  if (!NuMuOnly) {
    std::cout << "\n\t Not NuMuOnly, so pushing back the nue stuff...\n" << std::endl;
    if (horn=="fhc" || horn=="both") {
      predictionsNue  .push_back( GetNuePrediction2020(decomp_fhc, DefaultOscCalc(), isSyst, "fhc", FakeNDData) );
      cosmicsNue      .push_back( GetNueCosmics2020("fhc") );
    }
    if(horn=="rhc" || horn=="both") {
      predictionsNue.push_back( GetNuePrediction2020(decomp_rhc, DefaultOscCalc(), isSyst, "rhc", FakeNDData) );
      cosmicsNue    .push_back( GetNueCosmics2020("rhc") );
    }
    
    // --- Push the Nue stuff into the Cosmics and Preds vectors
    cosmicsvector.insert(cosmicsvector.end(), cosmicsNue    .begin(), cosmicsNue    .end());
    predsvector  .insert(predsvector  .end(), predictionsNue.begin(), predictionsNue.end());
  }

  // --- Load the Nu_mu things
  if (horn=="fhc" || horn=="both") {
    auto numu_preds_fhc   = GetNumuPredictions2020(nnumu, decomp_fhc, DefaultOscCalc(), isSyst, "fhc", FakeNDData);
    auto numu_cosmics_fhc = GetNumuCosmics2020    (nnumu, "fhc");
    predictionsNumu.insert(predictionsNumu.end(), numu_preds_fhc  .begin(), numu_preds_fhc  .end() );
    cosmicsNumu    .insert(cosmicsNumu    .end(), numu_cosmics_fhc.begin(), numu_cosmics_fhc.end() );
  }
  if (horn=="rhc" || horn=="both") {
    auto numu_preds_rhc   = GetNumuPredictions2020(nnumu, decomp_rhc, DefaultOscCalc(), isSyst, "rhc", FakeNDData);
    auto numu_cosmics_rhc = GetNumuCosmics2020    (nnumu, "rhc");
    predictionsNumu.insert( predictionsNumu.end(), numu_preds_rhc  .begin(), numu_preds_rhc  .end() );
    cosmicsNumu    .insert( cosmicsNumu    .end(), numu_cosmics_rhc.begin(), numu_cosmics_rhc.end() );
  }
  // --- Push the Numu stuff into the Cosmics and Preds vectors
  cosmicsvector.insert(cosmicsvector.end(), cosmicsNumu    .begin(), cosmicsNumu    .end());
  predsvector  .insert(predsvector  .end(), predictionsNumu.begin(), predictionsNumu.end()); 
 
  std::cout<<"\n\n Good job, you've loaded everything! \n\n"<< std::endl;

  //////////////////////////////////////////////////////////////
  ///////////////////// Check Predictions at first /////////////
  /////////////////////////////////////////////////////////////
  if (horn=="fhc" || horn=="both") {
    pot.push_back(potFD_fhc);
    livetime.push_back(cosFD_fhc);
  }
  if (horn=="rhc" || horn=="both") {
    pot.push_back(potFD_rhc);
    livetime.push_back(cosFD_rhc);
  }

  std::vector<TH1*> hNue; // nue signal
  std::vector<TH1*> hMC; //total MC
  std::vector<TH1*> hTotbkg; //total sum of all bkg
  std::vector<TH1*> hBkg; // all bkg except cosmics
  std::vector<TH1*> hNC;  //1
  std::vector<TH1*> hCC;  //2
  std::vector<TH1*> hBeam; //3
  std::vector<TH1*> hTau;  //4
  std::vector<TH1*> hCos;  //5
  std::vector<TH1*> hRock;  //6    rock events
  std::vector<TH1*> hNue_bg;  //7  wrong sign nue
  
  std::vector<TH1*> hNumu;
  std::vector<TH1*> hCosNumu;
  std::vector<TH1*> hMCNumu;

  std::cout<<">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Here are the predictions for 2020 nue ana with 2019 Best fit values: "<< std::endl;
  for(int i = 0; i < (int)predictionsNue.size(); i++){
    std::string sFHC;
    if (i==0){
      sFHC = "FHC";
      hNue   .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu    ).ToTH1(pot[i])); //signal FHC
      hNue_bg.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot[i])); //wrong sign bkg for FHC
    }
    if (i==1){
      sFHC = "RHC";
      hNue   .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot[i])); //signal RHC
      hNue_bg.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu    ).ToTH1(pot[i])); //wrong sign bkg for RHC
    }
    
    hTau .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(pot[i])); // tau CC bkg (2)
    hBeam.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(pot[i])); // beam bkg (3)
    hNC  .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAll     , Current::kNC, Sign::kBoth).ToTH1(pot[i])); // NC bkg (4)
    hCC  .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAllNuMu , Current::kCC, Sign::kBoth).ToTH1(pot[i])); // numu CC bkg (5)
    hMC  .push_back(predictionsNue[i]->Predict(calc).ToTH1(pot[i]));
    TH1D* temp = (TH1D*)hNue_bg[i]->Clone(UniqueName().c_str());
    temp->Add(hTau [i]);
    temp->Add(hBeam[i]);
    temp->Add(hNC  [i]);
    temp->Add(hCC  [i]);
    hBkg.push_back( (TH1D*)temp->Clone(UniqueName().c_str()));
    TH1D *tempCos = cosmicsNue[i].first->ToTH1(livetime[i], kBlack, kSolid, kLivetime);
    hCos.push_back( tempCos ); //Clone(UniqueName().c_str())); //cosmic bkg (6)
    std::cerr << "Done the cosmics..." << std::endl;
    temp  ->Add(hCos[i]);
    hMC[i]->Add(hCos[i]);
    hTotbkg.push_back(temp);
                   
    std::cout << sFHC << " mode, with " << pot[i] <<" POT, and " << livetime[i] << " s livetime." 
              <<"\n\t Total MC       " << hMC    [i]->Integral() 
              <<"\n\t Nu_e Signal    " << hNue   [i]->Integral() 
              <<"\n\t Total Bkg      " << hTotbkg[i]->Integral()
              <<"\n\t Total Beam bkg " << hBkg   [i]->Integral()
              <<"\n\t Wrong Sign bkg " << hNue_bg[i]->Integral()
              <<"\n\t Beam Nu_e      " << hBeam  [i]->Integral()
              <<"\n\t NC             " << hNC    [i]->Integral()
              <<"\n\t Nu_mu CC       " << hCC    [i]->Integral()
              <<"\n\t Nu_tau CC      " << hTau   [i]->Integral()
              <<"\n\t Cosmic         " << hCos   [i]->Integral()
              << std::endl;
  }

  std::cout<<">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Here are the predictions for 2020 numu ana with 2019 Best fit values: "<< std::endl;
  double AllQ_NuMu_Full, AllQ_Beam_Full, AllQ_Cosm_Full;
  double AllQ_NuMu_Dip , AllQ_Beam_Dip , AllQ_Cosm_Dip;
  AllQ_NuMu_Full = AllQ_Beam_Full = AllQ_Cosm_Full = AllQ_NuMu_Dip = AllQ_Beam_Dip = AllQ_Cosm_Dip = 0;
  for(int i = 0; i < (int)predictionsNumu.size(); i++){
    std::string sFHC = "FHC";
    double      pot  = potFD_fhc;
    double      live = cosFD_fhc;
    int quant = i+1;
    if (i>=4) {
      sFHC = "RHC";
      pot  = potFD_rhc;
      live = cosFD_rhc;
      quant= quant-4;
    }
    
    hNumu   .push_back(predictionsNumu[i]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));  
    hMCNumu .push_back(predictionsNumu[i]->Predict(calc).ToTH1(pot));
    TH1D *tempCos = cosmicsNumu[i].first->ToTH1(live, kBlack, kSolid, kLivetime);
    hCosNumu.push_back( tempCos ); //Clone(UniqueName().c_str()));

    double NuMu_Full = hNumu   [i]->Integral();
    double Beam_Full = hMCNumu [i]->Integral() - NuMu_Full;
    double Cosm_Full = hCosNumu[i]->Integral();
    double FOM_Full  = NuMu_Full / pow( NuMu_Full+Beam_Full+Cosm_Full, 0.5 );
    int bmin = hNumu[i]->GetXaxis()->FindBin(1.);
    int bmax = hNumu[i]->GetXaxis()->FindBin(2.);
    double NuMu_Dip  = hNumu   [i]->Integral( bmin, bmax );
    double Beam_Dip  = hMCNumu [i]->Integral( bmin, bmax ) - NuMu_Dip;
    double Cosm_Dip  = hCosNumu[i]->Integral( bmin, bmax );
    double FOM_Dip   = NuMu_Dip / pow( NuMu_Dip+Beam_Dip+Cosm_Dip, 0.5 );
    std::cout << sFHC << " mode, Quartile " << quant << " with " << pot <<" POT, and " << live << " s livetime."
              << "\n\t Nu_mu Signal Full: " << NuMu_Full << "\t Dip: " << NuMu_Dip 
              << "\n\t Beam Bkg     Full: " << Beam_Full << "\t Dip: " << Beam_Dip
              << "\n\t Cosmic       Full: " << Cosm_Full << "\t Dip: " << Cosm_Dip
              << "\n\t FOM_Full     Full: " << FOM_Full  << "\t Dip: " << FOM_Dip
              << std::endl;
    AllQ_NuMu_Full += NuMu_Full; AllQ_NuMu_Dip += NuMu_Dip;
    AllQ_Beam_Full += Beam_Full; AllQ_Beam_Dip += Beam_Dip;
    AllQ_Cosm_Full += Cosm_Full; AllQ_Cosm_Dip += Cosm_Dip;
    if (quant%4==0) {
      double AllQ_FOM_Full = AllQ_NuMu_Full / pow( AllQ_NuMu_Full+AllQ_Beam_Full+AllQ_Cosm_Full, 0.5 );
      double AllQ_FOM_Dip  = AllQ_NuMu_Dip  / pow( AllQ_NuMu_Dip +AllQ_Beam_Dip +AllQ_Cosm_Dip , 0.5 );
      std::cout << sFHC << " mode, All Quartile with " << pot <<" POT, and " << live << " s livetime."
                << "\n\t Nu_mu Signal Full: " << AllQ_NuMu_Full << "\t Dip: " << AllQ_NuMu_Dip 
                << "\n\t Beam Bkg     Full: " << AllQ_Beam_Full << "\t Dip: " << AllQ_Beam_Dip
                << "\n\t Cosmic       Full: " << AllQ_Cosm_Full << "\t Dip: " << AllQ_Cosm_Dip
                << "\n\t FOM_Full     Full: " << AllQ_FOM_Full  << "\t Dip: " << AllQ_FOM_Dip
                << "\n" << std::endl;
      AllQ_NuMu_Full = AllQ_Beam_Full = AllQ_Cosm_Full = AllQ_NuMu_Dip = AllQ_Beam_Dip = AllQ_Cosm_Dip = 0;
    }
  }
  

  std::cout<<"\n\n Good job, you've got your integrals, now to make sensitivites! \n\n"<< std::endl;
  /////////////////////////////////////////////////////////////
  /////////////////// Sensitivities ///////////////////////////
  ////////////////////////////////////////////////////////////

  // TODO Extrap, Extrap + Syst; also nue, numu, both; RHC, FHC, both
  // species from 2019 TN: th23 vs delta; th23 vs dmsq23; slices: hie, dcp, octant, max mixing

  //make fakedata at 2019 BF(make choosable)
  //ResetOscCalcToDefault(calc);
  //calc->SetdCP(1.21*M_PI);
  //calc->SetTh23(asin(sqrt(0.558)));
  //calc->SetDmsq32(2.44e-3);

  std::vector <const IExperiment * > expts;
  std::cout<<"Create fake data with systematic pred. and multiexp"<< std::endl;
  for(int i = 0; i < (int) predsvector.size(); ++i){
    double pot  = potFD_fhc;
    double live = cosFD_fhc;
    // Cases for when I need RHC pot/livetime.
    if(horn=="both" && ( (!NuMuOnly && (i==1 || i>5) ) || // Nue+NuMu  then vector is (nueFHC, nueRHC, numuFHC*4, numuRHC*4)
                         ( NuMuOnly && i>3 )              // NuMu only then vector is (numuFHC*4, numuRHC*4)
                         ) ) {
      pot  = potFD_rhc;
      live = cosFD_rhc;
    }
    std::cout << "\t My POT is " << pot << ", My Livetime is " << live << std::endl;       
    auto thisdata = GetFakeData ( predsvector[i], calc, pot, cosmicsvector[i].first );
    std::cout << "\t Fake data integral is: "<< thisdata->Integral(pot) << std::endl;
    expts.push_back(new SingleSampleExperiment( predsvector[i], *thisdata, *cosmicsvector[i].first, cosmicsvector[i].second, true) );
  }

  expts.push_back(WorldReactorConstraint2019());
  auto exptsAll = new MultiExperiment(expts);   
  std::cout<<"Experiments are ready, number of items in vector is "<<expts.size()<<" (the above plus WordReactorConstraint)" << std::endl;

  // --- Set systematic correlations if requested --- blindly copying out for now...
  std::cout<<"\n Now I want to set your systematics"<< std::endl;
  std::vector< const ISyst * > systs = {};
  SystShifts auxShifts = SystShifts::Nominal();
  std::vector <SystShifts> seedShifts = {};

  if(isSyst){
    std::cout << "\t Getting syst list and the correlations." << std::endl;
    if (NuMuOnly) {
      std::cout << "\t\t NuMuOnly, so get a reduced list of systs..." << std::endl;
      systs = GetJointFitSystematicList(isSyst, false, true, true, true, true, isPtExtrap);
      auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
      auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
      for (size_t i=0; i<expts.size(); ++i) {
        if (i<4) {
          exptsAll->SetSystCorrelations(i, notfornumuFHC);
        } else {
          if (i==8) continue;
          exptsAll->SetSystCorrelations(i, notfornumuRHC);
        }
      }
    } else {
      std::cout << "\t\t NuMu and Nue, so get all of the systs..." << std::endl;
      systs = GetJointFitSystematicList(isSyst, false, false, true, true, true, isPtExtrap);
      auto notfornueFHC  = GetCorrelations(true , true , isPtExtrap);
      auto notfornueRHC  = GetCorrelations(true , false, isPtExtrap);
      auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
      auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
      if(horn=="fhc"){
        exptsAll->SetSystCorrelations(0, notfornueFHC);
        for(int i=1; i<5;i++) {
          exptsAll->SetSystCorrelations(i, notfornumuFHC);
        }
      }
      if(horn=="both"){
        exptsAll->SetSystCorrelations(0, notfornueFHC);
        exptsAll->SetSystCorrelations(1, notfornueRHC);
        for(int i=0; i < 4; ++i) exptsAll->SetSystCorrelations(i+2, notfornumuFHC);
        for(int i=0; i < 4; ++i) exptsAll->SetSystCorrelations(i+6, notfornumuRHC);
      }
      
      auxShifts = SystShifts::Nominal();
      //Registry::Print();
    }
  }
  else { 
    std::cout<<"\t You're not setting any systematics"<< std::endl;
    systs={}; 
    //auxShifts = junkShifts;
    auxShifts = SystShifts();
  }

  // --- find BF here:
  std::cout << "\n Now to find the Best Fit." << std::endl;
  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13}; 
  MinuitFitter fit(exptsAll, fitvars, systs); 

  double minchi23 = 1E20;
  ResetOscCalcToDefault(calc);
  for(auto & th23seed:th23_seeds){
    for(int hie:{-1,1}){
      std::cout << "Finding best fit for " << (th23seed<0.5?"LO ":"UO ") << (hie>0? "NH " : "IH ") << std::endl;
      auto thisminchi = fit.Fit( calc, auxShifts ,
                                 SeedList({ 
                                     { &kFitDmSq32        , {hie*fabs(calc->GetDmsq32())} },
                                     { &kFitSinSqTheta23  , {th23seed} },
                                     { &kFitDeltaInPiUnits, delta_seeds  }
                                   }), 
                                 seedShifts, IFitter::kVerbose
                                 );
      minchi23= min(minchi23, thisminchi);
      ResetOscCalcToDefault(calc);
      if(isSyst){
        auxShifts = SystShifts::Nominal(); 
        std::cout<<"reset syst as well"<< std::endl;
      }
    }//end hie
  }//end th23
  std::cout << "Found overall min chi " << minchi23 << "\n" << std::endl;

  std::string FileDir = "";
  //std::string FileDir = "/nova/ana/3flavor/Ana2020/FutureSensitivities/";
  std::string PlotDir = "";
  if (!OnGrid) {
    FileDir = "files/";
    PlotDir = "plots/";
    gSystem->MakeDirectory( FileDir.c_str() );
    gSystem->MakeDirectory( PlotDir.c_str() );
  }
  // -----------------------------------
  // contours from 2019 TN --> Make Surf
  if(!MakeSurf){
    std::cout << "\n NOT making surfaces.\n" << std::endl;
  } else {

    int hie = 1;
    if(!isNH) { hie = -1; }

    std::string FileName = FileDir+outfilename+"_Surfonly_"+whatsurf+".root";
    std::string PlotName;
    if (whatsurf == "th23dcp" ) PlotName = TagName+"-surf_delta_th23";
    if (whatsurf == "dmsqth23") PlotName = TagName+"-surf_th23_dmsq23";

    std::cout << "\n You now want to make the surfaces called " << PlotName << ". They will be saved/loaded from " << FileName << std::endl;

    if(MakeFile){
      std::cout << "You are going to make the file containing the surfaces." << std::endl;
      TFile * outfile = new TFile( FileName.c_str(), "recreate" );
      outfile->cd();
      TVectorD v(1);
      v[0] = minchi23;
      v.Write("overall_minchi");

      std::cout<<"\t Clean the vectors"<< std::endl;
      predsvector.clear();

      ResetOscCalcToDefault(calc);
      calc->SetDmsq32(hie*fabs(BASE_dmsq32));
      if(isSyst) auxShifts = SystShifts::Nominal();
      
      // --- th23dcp
      if (whatsurf=="th23dcp") {
        int Ysteps = 30; int Xsteps = 30;
        std::cout<<"\t Doing th23dcp --> minimize with grid "<< Ysteps<<" * "<<Xsteps<<" points"<< std::endl;
        
        FrequentistSurface surf_th23_dcp (exptsAll, calc,
                                          &kFitDeltaInPiUnits, Xsteps, 0  , 2,
                                          &kFitSinSqTheta23  , Ysteps, 0.2, 0.7,
                                          {&kFitDmSq32Scaled, &kFitSinSq2Theta13},
                                          systs,
                                          SeedList({{&kFitSinSq2Theta13, {calc->GetTh13()}}})
                                          //,seedShifts
                                          );
        outfile->cd();
        surf_th23_dcp.SaveTo(outfile,  PlotName.c_str() ) ;
      }
      // --- dmsqth23
      else if (whatsurf=="dmsqth23") {
        int Ysteps = 26; int Xsteps = 26;
        std::cout<<"\t Doing dmsqth23 --> minimize with grid "<< Ysteps<<" * "<<Xsteps<<" points"<< std::endl;

        FrequentistSurface surf_th23_dmsq23 (exptsAll, calc,
                                             &kFitSinSqTheta23, Xsteps, 0.3, 0.7,
                                             &kFitDmSq32Scaled, Ysteps, (hie>0?2.1:-2.8), (hie>0?2.8:-2.1),
                                             {&kFitSinSq2Theta13, &kFitDeltaInPiUnits},
                                             systs,
                                             SeedList({{&kFitDeltaInPiUnits, delta_seeds}})
                                             //,seedShifts
                                             );
        outfile->cd();
        auto surf_test=surf_th23_dmsq23.ToTH2(minchi23);
        surf_test ->Write( TString(PlotName)+TString("_2DTest") );
        surf_th23_dmsq23.SaveTo( outfile,  PlotName.c_str() ) ;
      }
      outfile->Close();
      std::cout<<"Done with surfaces!"<< std::endl;
    }

    // ----> If not Make File, or once finished the above loop.
    std::cout << "You are going to load premade surfaces." << std::endl;
    TFile * infile = new TFile( FileName.c_str(),"read");
    
    auto mins =* (TVectorD*)infile->Get("overall_minchi");
    double minchi23 = mins[0];
    std::cout<<"Min chi: "<<minchi23<< std::endl;
    
    Int_t kFillColor1 = TColor::GetColorTransparent(k1SigmaConfidenceColorNH, 0.70);
    Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColorNH, 0.35);
    Int_t kFillColor3 = TColor::GetColorTransparent(k3SigmaConfidenceColorNH, 0.20);
    if (!isNH) {
      kFillColor1 = TColor::GetColorTransparent(k1SigmaConfidenceColorIH, 0.70);
      kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColorIH, 0.35);
      kFillColor3 = TColor::GetColorTransparent(k3SigmaConfidenceColorIH, 0.20);
    }    
    auto leg = ContourLegend(hie, false, false, kFillColor1, kFillColor2, kFillColor3, kBlack, false);
    
    // --- th23dcp
    if (whatsurf=="th23dcp"){
      auto surf_th23_dcp = *FrequentistSurface::LoadFrom(infile,  PlotName.c_str() ) ;
      
      TH2 * surf_th23_dcp_1Sigma = Gaussian68Percent2D(surf_th23_dcp);
      TH2 * surf_th23_dcp_2Sigma = Gaussian2Sigma2D   (surf_th23_dcp);
      TH2 * surf_th23_dcp_3Sigma = Gaussian3Sigma2D   (surf_th23_dcp);
      DrawContourCanvas("delta_th23", 0., 2., 0.225, 0.725);
      surf_th23_dcp.DrawContour(surf_th23_dcp_1Sigma, kSolid                , kFillColor1, minchi23);
      surf_th23_dcp.DrawContour(surf_th23_dcp_2Sigma, k2SigmaConfidenceStyle, kFillColor2, minchi23);
      surf_th23_dcp.DrawContour(surf_th23_dcp_3Sigma, kSolid                , kFillColor3, minchi23);
      
      leg->Draw();
      Simulation();
      
      gPad->Print((TString)PlotDir+PlotName+".pdf");
    }
    // --- dmsqth23
    else if (whatsurf=="dmsqth23") {
      auto surf_th23_dmsq23 = *FrequentistSurface::LoadFrom(infile,  PlotName.c_str() );
      
      TH2 * surf_th23_dmsq23_1Sigma = Gaussian68Percent2D(surf_th23_dmsq23);
      TH2 * surf_th23_dmsq23_2Sigma = Gaussian2Sigma2D   (surf_th23_dmsq23);
      TH2 * surf_th23_dmsq23_3Sigma = Gaussian3Sigma2D   (surf_th23_dmsq23);
      
      if (hie>0) DrawContourCanvas("th23_dm32", 0.3, 0.7,  2.05,  2.75);        
      else       DrawContourCanvas("th23_dm32", 0.3, 0.7, -2.95, -2.25);
      
      surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_1Sigma, kSolid                , kFillColor1, minchi23);
      surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_2Sigma, k2SigmaConfidenceStyle, kFillColor2, minchi23);
      surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_3Sigma, kSolid                , kFillColor3, minchi23);
      
      leg->Draw();
      Simulation();
      
      gPad->Print((TString)PlotDir+PlotName+".pdf");
    }
  } //MakeSurf end
  // contours from 2019 TN --> Make Surf
  // -----------------------------------

  // -----------------------------------
  // Now for the Make Sens loop
  if(!MakeSens){
    std::cout << "\n NOT making any sensitivities.\n" << std::endl;
  } else {    
    std::vector<std::pair<TGraph*, TString>> slice;
    std::vector<std::pair<TGraph*, TString>> slice2019;
    slice.clear();
    slice2019.clear();

    struct th23helper{
      std::vector<double> seeds;
      const IFitVar * var;
      TString label;
    };
    std::vector <th23helper> th23seeds;
    th23seeds.push_back( { {0.4,0.6}, &kFitSinSqTheta23           , "2019" } );
    th23seeds.push_back( { {0.45}   , &kFitSinSqTheta23LowerOctant, "LO SA"} );
    th23seeds.push_back( { {0.55}   , &kFitSinSqTheta23UpperOctant, "UO SA"} ); 

    int id, revid;
    std::string sVary = Varyth23 ? "Varyth23":"VarydCP";
    std::string FileName = FileDir+outfilename+"_"+sVary+"_Slicesonly_"+whatplot+".root";
    
    if(MakeFile){
      std::cout << "\nYou are going to make the file containing the sensitivites." << std::endl;
      TFile * outfile = new TFile( FileName.c_str(), "recreate" );
      outfile->cd();

      TVectorD v(1);
      v[0] = minchi23;
      v.Write("overall_minchi");

      th23helper th23seed = th23seeds[0];
      
      double th23  = BASE_th23;
      double delta = BASE_dcp;

      for(int hie : {-1,1}) {
        std::cout << "\nNow looking at hie = " << hie << std::endl;
        int steps = 25; double chi2r[steps], chi2w[steps], chi[steps]; double Fits[steps];
        double step = 2.0/(steps-1);
        for (int i =0; i<steps; i++){
          if ( Varyth23 ) { th23  = (i*step*0.1 ) + 0.4; }
          else            { delta = (i*step*M_PI);       } 

          ResetOscCalcToDefault(calc);
          calc  ->SetdCP   ( delta );
          calc  ->SetTh23  ( asin(sqrt(th23)) );
          calc  ->SetDmsq32( hie*fabs(BASE_dmsq32) );
 
            std::vector <const IExperiment * > expts_temp;
            for(int i = 0; i < (int) predsvector.size(); ++i){

              double pot = potFD_fhc;
              if(horn=="both" && ( (!NuMuOnly && (i==1 || i>5) ) || // Nue+NuMu  then vector is (nueFHC, nueRHC, numuFHC*4, numuRHC*4)
                                   ( NuMuOnly && i>3 )              // NuMu only then vector is (numuFHC*4, numuRHC*4)
                                   ) ) {
                pot  = potFD_rhc;
              }
              auto thisdata = GetFakeData (predsvector[i],calc, pot, cosmicsvector[i].first);
              expts_temp.push_back(new SingleSampleExperiment(predsvector[i], *thisdata, *cosmicsvector[i].first, cosmicsvector[i].second, true));
            }
            expts_temp.push_back(WorldReactorConstraint2019()); 
            auto exptsAll_temp = new MultiExperiment(expts_temp);
            auxShifts = SystShifts::Nominal();

            std::cout << "\tStep " << i << " of " << steps << ", th23 = " << th23 << ", delta = " << delta << " ("<<delta/M_PI<<"*M_PI)" << std::endl;

            if(isSyst){
              std::cout << "\t Getting syst list and the correlations." << std::endl;
              if (NuMuOnly) {
                std::cout << "\t\t NuMuOnly, so get a reduced list of systs..." << std::endl;
                systs = GetJointFitSystematicList(isSyst, false, true, true, true, true, isPtExtrap);
                auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
                auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
                for (size_t i=0; i<expts.size(); ++i) {
                  if (i<4) {
                    exptsAll_temp->SetSystCorrelations(i, notfornumuFHC);
                  } else {
                    if (i==8) continue;
                    exptsAll_temp->SetSystCorrelations(i, notfornumuRHC);
                  }
                }
              } else {
                systs = GetJointFitSystematicList(isSyst, false, false, true, true, true, isPtExtrap);
                auto notfornueFHC = GetCorrelations(true , true , isPtExtrap);
                auto notfornueRHC = GetCorrelations(true , false, isPtExtrap);
                auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
                auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
                if(horn=="fhc"){
                  exptsAll_temp->SetSystCorrelations(0, notfornueFHC);
                  for(int i=1; i<5;i++) {
                    exptsAll_temp->SetSystCorrelations(i, notfornumuFHC);
                  }
                }
                if(horn=="both"){
                  exptsAll_temp->SetSystCorrelations(0, notfornueFHC);
                  exptsAll_temp->SetSystCorrelations(1, notfornueRHC);
                  for(int i=0; i < 4; ++i) exptsAll_temp->SetSystCorrelations(i+2, notfornumuFHC);
                  for(int i=0; i < 4; ++i) exptsAll_temp->SetSystCorrelations(i+6, notfornumuRHC);
                }
                //Registry::Print();
              }
            }
          
            // ----- Now to do the MH sensitivities
            if( whatplot=="mh" ){ 
              std::cout << "\t\t Now to make the MH sensitivities...First the Normal Hierarchy." << std::endl; 

              std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13}; 
              MinuitFitter fit(exptsAll_temp, fitvars, systs); 
              
              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              chi2r[i] = fit.Fit(calc, auxShifts,
                                 SeedList({
                                     { &kFitDeltaInPiUnits, delta_seeds },
                                     { &kFitSinSqTheta23  , th23_seeds  }
                                 }), {}, IFitter::kQuiet
                                 );
              std::cout << "\t\t And now the Inverse Hierarchy." << std::endl;
              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(-hie*fabs(BASE_dmsq32));
              chi2w[i] = fit.Fit(calc, auxShifts,
                                 SeedList({ 
                                     { &kFitDeltaInPiUnits, delta_seeds },
                                     { &kFitSinSqTheta23  , th23_seeds  }
                                 }), {}, IFitter::kQuiet);
            }
            
            // --- Now to do the Delta CP sensitivities
            if( whatplot=="dcp" ){ 
              std::cout << "\t\t Now to make the dCP sensitivities...First do DeltaCP =  delta = " << delta << std::endl; 

              std::vector <const IFitVar *> fitvars = {&kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13}; 
              MinuitFitter fit(exptsAll_temp, fitvars, systs); 

              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              calc->SetdCP(delta);
              chi2r[i] = fit.Fit(calc, auxShifts,
                                 SeedList({
                                     { &kFitSinSqTheta23, th23_seeds},
                                     { &kFitDmSq32      , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
                                 }), {}, IFitter::kQuiet);
              std::cout << "\t\t And now do DeltaCP = 0" << std::endl;
              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              calc->SetdCP(0);
              double chi2w0 = fit.Fit(calc, auxShifts, 
                                      SeedList({
                                          { &kFitSinSqTheta23, th23_seeds},
                                          { &kFitDmSq32      , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
                                      }), {}, IFitter::kQuiet);
              std::cout << "\t\t And now do DeltaCP = PI" << std::endl;
              ResetOscCalcToDefault(calc);
              calc->SetdCP(M_PI);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              double chi2wpi = fit.Fit(calc, auxShifts, 
                                       SeedList({
                                           { &kFitSinSqTheta23, th23_seeds},
                                           { &kFitDmSq32      , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
                                        }), {}, IFitter::kQuiet);
              // And find the minimum value.
              chi2w[i] = std::min(chi2w0, chi2wpi);
              std::cout << "The smallest chi2 was " << chi2w[i] << std::endl;
            }
          
            // --- Now to do the Maximal Mixing sensitivities
            if( whatplot=="maxmix" ){
              std::cout << "\t\t Now to make the Maximal Mixing sensitivities...First do theta23 = " << th23 << std::endl; 
              
              std::vector <const IFitVar *> fitvars =   {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13}; 
              std::vector <const IFitVar *> fitvars23 = {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, &kFitSinSqTheta23 }; 
              MinuitFitter fit  (exptsAll_temp, fitvars  , systs);
              MinuitFitter fit23(exptsAll_temp, fitvars23, systs);

              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              calc->SetTh23  (asin(sqrt(th23)));
              chi2r[i] = fit23.Fit(calc, auxShifts, 
                                 SeedList({
                                     { &kFitDeltaInPiUnits, delta_seeds },
                                     { &kFitDmSq32        , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }, 
                                     { &kFitSinSqTheta23  , th23_seeds }
                                 }), {}, IFitter::kQuiet);
              std::cout << "\t\t And now do theta23 = 0.5" << std::endl; 
              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              calc->SetTh23(asin(sqrt(0.5)));
              chi2w[i] = fit.Fit(calc, auxShifts, 
                                 SeedList({
                                     { &kFitDeltaInPiUnits, delta_seeds},
                                     { &kFitDmSq32        ,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
                                 }), {}, IFitter::kQuiet);
            }
   
            // --- Now to do the Octant sensitivities
            if ( whatplot=="oct" ) {
              std::cout << "\t\t Now to make the octant sensitivities...First use the correct octant for th23 " << th23 << std::endl; 
              // Primitive flip between octants: minimize LO with kFitSinSqTheta23LowerOctant, 
              //  then use wrong octant and minimize there with kFitSinSqTheta23UpperOctant
              id=2; 
              revid=1;

              std::vector <const IFitVar *> fitvars1 = {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, th23seeds[id   ].var};
              std::vector <const IFitVar *> fitvars2 = {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, th23seeds[revid].var};
              
              MinuitFitter fit1(exptsAll_temp, fitvars1, systs); // Correct oct
              MinuitFitter fit2(exptsAll_temp, fitvars2, systs); // Wrong oct

              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              calc->SetTh23(asin(sqrt(th23)));
              chi2r[i] = fit1.Fit(calc, auxShifts, 
                                  SeedList({
                                      { &kFitDeltaInPiUnits, delta_seeds },
                                      { &kFitDmSq32        ,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }, 
                                      { th23seeds[id].var  , th23seeds[id].seeds }
                                  }), {}, IFitter::kQuiet);
              std::cout << "\t\t And now, use the wrong octant " << 1-th23 << std::endl;
              ResetOscCalcToDefault(calc);
              calc->SetDmsq32(hie*fabs(BASE_dmsq32));
              calc->SetTh23(asin(sqrt(1-th23)));
              chi2w[i] = fit2.Fit(calc, auxShifts, 
                                  SeedList({
                                      { &kFitDeltaInPiUnits, delta_seeds },
                                      { &kFitDmSq32        ,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }, 
                                      { th23seeds[revid].var  , th23seeds[revid].seeds }
                                  }), {}, IFitter::kQuiet);
            } 
          
            if (Varyth23) Fits[i] = th23;
            else          Fits[i] = delta/M_PI;
            chi[i] = sqrt(fabs(chi2w[i]-chi2r[i]));
   
            //std::cout<<delt<<" chi2 "<<chi[i]<<" right "<<chi2r[i]<<" wrong assumption :"<<chi2w[i]<< std::endl;
          }
          TString name = th23seed.label+" "+(hie>0?"NH":"IH");
          slice.push_back({new TGraph(steps, Fits, chi), name}); //0 is ih, 1 is nh

      } // end hie
      outfile->cd();

      TVectorD ss(1);
      ss[0] = slice.size();
      ss.Write("numsl");

      for(int i=0; i<(int) slice.size(); i++){
        slice[i].first->Write(slice[i].second);
      }
    
      outfile->Close();
      std::cout<<"Done with slices"<< std::endl;
    }//end MakeFile

    // ----> If not Make File, or when that loop finishes.
    std::cout << "You are going to read the file containing the sensitivities -- " << FileName << std::endl;
    TFile * infile = new TFile (FileName.c_str(),"read");
    auto mins =* (TVectorD*)infile->Get("overall_minchi");
    double minchi23 = mins[0];
    std::cout<<"Min chi: "<<minchi23<< std::endl;
    
    auto ss =* (TVectorD*)infile->Get("numsl");
    double numsl = ss[0];
    std::cout<<"Num slices: "<<numsl<< std::endl;
    
    // Clear the slice vectors.
    slice.clear();
    slice2019.clear();
    
    for (double th23 : {0.565}){ //2019 BF
      
      th23helper th23seed;
      if(th23 == 0.565) {th23seed = th23seeds[0];}
      if(th23 == 0.404) {th23seed = th23seeds[1];}
      if(th23 == 0.623) {th23seed = th23seeds[2];}
      
      for(int hie : {-1,1}){
        TString name = th23seed.label+" "+(hie>0?"NH":"IH");    
        TGraph* get_slice;
        infile->GetObject(name, get_slice);
        if(name.Contains("2019")) slice2019.push_back({get_slice,name});
        else slice.push_back({get_slice,name});
      }
    }

    new TCanvas;    
    double max = 5.1, xmin = 0, xmax = 2;
    if (Varyth23) {
      xmin = 0.4;
      xmax = 0.6;
      if(whatplot=="oct"                   ) max = 4.1;
      if(whatplot=="maxmix" && isFutureSens) max = 6.1;
    } else {
      if(whatplot=="dcp"   )              max = 3.1;
      if(whatplot=="maxmix")              max = 1.5;
      if(whatplot=="oct" && NuMuOnly    ) max = 0.5;
      if(whatplot=="mh"  && isFutureSens) max = 6.1;
    }
    DrawSliceCanvas(whatplot, max, xmin, xmax, false, false, isFutureSens );
    for(unsigned int i =0; i < slice2019.size(); i++){
      
      if(slice2019[i].second.Contains("IH")) slice2019[i].first->SetLineColor(kPrimColorIH);
      if(slice2019[i].second.Contains("NH")) slice2019[i].first->SetLineColor(k3SigmaConfidenceColorNH);
      
      slice2019[i].first->SetLineWidth(3);
      slice2019[i].first->Draw("same");
    }
    
    auto leg2 = SliceLegend(slice2019, (whatplot!="dcp"&&whatplot!="mh"), 0.53, whatplot=="mh", false, isSyst, isFutureSens);
    leg2->Draw();
    gPad->SetGrid();
    if (isFutureSens) {
      latex(0.138, 0.94, "NOvA 2020 Best Fit"                                     , 0, 12, 1.3/30.);
      latex(0.9  , 0.94, "#nu 31.5#times10^{20} + #bar{#nu} 31.5#times10^{20} POT", 0, 32, 1.3/30.);
      SimulationSide();
    } else {
      Simulation();
    }
    gPad->Print((TString)PlotDir+TagName+"-Sens-"+sVary+"-"+whatplot+".pdf");
    
  }//end MakeSens
}