make_fc_surfaces_2020.C

Go to the documentation of this file.

#include <fstream>

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Utilities.h"

#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Exposures.h"

#include "CAFAna/Prediction/IPrediction.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"

#include "CAFAna/Vars/FitVars.h"

#include "CAFAna/FC/FCPoint.h"
#include "CAFAna/FC/FCCollection.h"
#include "CAFAna/FC/FCSurface.h"

#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"

#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"

#include "CAFAna/Analysis/Calcs.h"
#include "OscLib/IOscCalc.h"

#include "TRandom3.h"
#include "TH1.h"

namespace fc {
  enum kBeam_t {
    kRHC,
    kFHC
  };
  enum kFlavor_t {
    kNue,
    kNumu
  };
}

using namespace ana;

void SetSeeds(osc::IOscCalcAdjustable* calc, double delta, double th23, double dmsq)
{
  ResetOscCalcToDefault(calc);
  kFitDeltaInPiUnits.SetValue(calc, delta);
  kFitSinSqTheta23.SetValue(calc, th23);
  kFitDmSq32.SetValue(calc, dmsq);
}

void make_fc_surfaces_2020(int NPts, int startidx, int endidx, bool nh,
                           int N, int jid, int npts_flush = 0,
                           std::string plot = "deltassth23" ,
                           std::string options = "joint_fake2019" ,
                           std::string decomp = "",
                           TString beam = "both", bool corrSysts = false,
                           bool readFromFile = false, bool runOnGrid=false)
{
  assert(plot == "deltath13" || plot=="ssth23dmsq32" ||
         plot == "deltassth23" && "Plot must be ssth23dmsq32 or deltassth23");

  assert(beam == "FHCOnly" || beam == "RHCOnly" || beam == "both");

  if(readFromFile && !corrSysts){
    std::cout << "No bin list for stats only contours" << endl;
    //    exit(1);
  }

  //also treat the decomp option to get the right FCInput name

  bool fake2019 = options.find("fake2019") != std::string::npos;
  bool realData = options.find("realData") != std::string::npos;

  bool RHCOnly = beam.Contains("RHCOnly");
  bool FHCOnly = beam.Contains("FHCOnly");
  bool both = beam.Contains("both");

  bool nueOnly = options.find("nueOnly") != std::string::npos;
  bool numuOnly = options.find("numuOnly") != std::string::npos;
  bool joint = options.find("joint") != std::string::npos;

  bool isFHC = true;
  bool isRHC = true;

  // ssth23dmsq32 has 400 total bins while deltassth23 has 900
  int TotBins = 0;
  TotBins = (plot == "deltassth23") ? 900 : 400;
  int bins[TotBins];
  std::memset(bins, 0, TotBins);
  std::string hier = nh? "_NH_" : "_IH_";
  std::string stat = corrSysts? "_systs": "_stats";


  // If we don't want to run full corrections, we'll read from a text file which holds
  // a list of bins to be corrected based on deltachisq in each bin. This list is
  // generated by generate_fc_binlists.C

  ifstream binlist;
  // std::string path = getenv("FCHELPERANA2019_LIB_PATH");
  std::string path = "/nova/app/users/scalvez/fchelper/";
  path += "/binlists/";
  // TString binlist_fname = path+plot+hier+beam+stat+"_binlist.txt";
  TString binlist_fname = path+plot+hier+"binlist.txt";

  binlist.open(binlist_fname);
  int j = 0;
  // check if the binlist exists
  if(!binlist.good()){
    std::cout << "Could not find "<< binlist_fname << endl;
    exit(1);
  }
  if(readFromFile){
    do
      {
        binlist >> bins[j];
        j++;
      }while(!binlist.eof());
    binlist.close();
  }
  else for (int ibin = 0; ibin < TotBins; ibin++) bins[ibin]=ibin;

  int ibin = bins[startidx];
  int fbin = bins[endidx];

  std::string binname = "bins_"+std::to_string(ibin) + "_to_" + std::to_string(fbin);
  std::string number = std::to_string(N);
  std::string fcout = "FCCol_"+plot+"_"+hier+"_";
  if(both) fcout+="both";
  else fcout += RHCOnly? "RHCOnly":"FHCOnly";
  fcout += corrSysts? "_systs_":"_stats_";
  fcout += binname+"_v"+number+'_'+std::to_string(jid)+".root";

  std::string hierStr = nh ? "NH" : "IH";

  FCCollection fccol;

  // Relevant exposures, to be used later
  double thisPOT;
  double thisLivetime;

  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  if(fake2019) SetFakeCalc(calc, 0.564956, 2.48304e-3, 1.9997);
  // if(fake2019) SetFakeCalc(calc, 0.565, 2.48e-3, 0);
  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}

  // struct predictions {
  //   const string name;
  //   const IPrediction * pred;
  //   std::pair <Spectrum*, double> cos;
  //   double pot;
  //   double livetime;
  //   const fc::kBeam_t             beam;
  //   const fc::kFlavor_t           flavor;
  // };

  // std::vector <predictions> preds;
  // std::vector <predictions> preds_numu_only;
  // std::vector <Spectrum * > data;
  // std::vector <Spectrum * > data_numu_only;

  std::vector <ana::predictions> preds = LoadPredictions (corrSysts, runOnGrid, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  std::vector <ana::predictions> preds_numu_only = LoadPredictions (corrSysts, runOnGrid, decomp, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  // std::vector <Spectrum * > data = LoadRealData (runOnGrid, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  // std::vector <Spectrum * > data_numu_only = LoadRealData (runOnGrid, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);

  // We're setting up the requisite inputs to the nue experiment
  /*
  //////////////////////////////////////////////////////////////
  //  Nue Predictions
  /////////////////////////////////////////////////////////////

  if(!numuOnly) {
    if(FHCOnly || both ) {
      preds.push_back({"nue_fhc",
            GetNuePrediction2020(decomp+"combo", DefaultOscCalc(), corrSysts, "fhc",
                                 true, runOnGrid),
            GetNueCosmics2020("fhc", runOnGrid),
            GetPOT(), GetLT(),
            fc::kFHC, fc::kNue});
    }
    if(RHCOnly || both ) {
      preds.push_back({"nue_rhc",
            GetNuePrediction2020("prop", DefaultOscCalc(), corrSysts, "rhc",
                                 false, runOnGrid),
            GetNueCosmics2020("rhc", runOnGrid),
            GetPOT(false), GetLT(false),
            fc::kRHC, fc::kNue});
    }
  }
  std::cout << "\n------------------Done Getting nue preds ---------------------\n";

  ///////////////////////////////////////////////////////////////
  // Numu predictions
  //////////////////////////////////////////////////////////////

  // Construct numu inputs
  int nnumu = 4;
  if(!nueOnly) {
    for(fc::kBeam_t _beam : {fc::kFHC, fc::kRHC}) {
      if(((_beam == fc::kRHC) && FHCOnly) ||
         ((_beam == fc::kFHC) && RHCOnly))
        continue;
      double _pot      = GetPOT(_beam == fc::kFHC);
      double _livetime = GetLT (_beam == fc::kFHC);
      std::string beam_name = _beam == fc::kFHC ? "fhc" : "rhc";
      auto numu_preds   = GetNumuPredictions2020(nnumu, decomp, DefaultOscCalc(), corrSysts, beam_name, false, false);
      auto numu_cosmics = GetNumuCosmics2020    (nnumu,            beam_name, false);
      for (int i = 0; i< nnumu; i++ ){
        preds.push_back({"numu_" + beam_name + "_" +std::to_string(i+1),
              numu_preds[i],
              numu_cosmics[i],
              _pot,
              _livetime,
              _beam,
              fc::kNumu});
        preds_numu_only.push_back({"numu_" + beam_name + "_"+std::to_string(i+1),
              numu_preds[i],
              numu_cosmics[i],
              _pot,
              _livetime,
              _beam,
              fc::kNumu});
      }
    }
  }
  */
  // Grab the systs
  bool ptExtrap = false;
  // std::vector<const ISyst*> systs = {};
  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, ptExtrap);
  // if(corrSysts) systs = GetJointFitSystematicList(true, !nueOnly,
  //                                              !numuOnly, isFHC, isRHC);
  // Equivalent for 2020?
  // need to grab different systs for nue rhc and fhc
  std::vector<const ISyst*> snue_fhc = {};
  // if(corrSysts) snue_fhc = getAllAna2018Systs(kNueAna2018, true, kFHC);
  std::vector<const ISyst*> snue_rhc = {};
  // if(corrSysts) snue_rhc = getAllAna2018Systs(kNueAna2018, true, kRHC);

  // numu systs are independent of beam mode
  std::vector<const ISyst*> snumu = {};
  // snumu = GetJointFitSystematicList(true, !nueOnly,
  //                                   numuOnly,
  //                                   isFHC,
  //                                   isRHC);

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

  // Grab analysis contours to find bin center
  std::string p = ""; //delta && !corrSysts
  //  if(plot == "ssth23dmsq32" && corrSysts) p = "ssth23dmsq32_";
  //  if(plot == "deltassth23" && corrSysts) p = "deltassth23_";
  if(plot == "ssth23dmsq32") p = "dmsq_";
  if(plot == "deltath13") p = "th13_";

  // std::string helperupsname = std::string(getenv("FCHELPERANA2019_LIB_PATH"));
  std::string helperupsname = "/nova/app/users/scalvez/fchelper/";
  TString anasurfname = helperupsname
    + "/Surfaces/hist_contours_2019_joint_realData_"+beam+"_only";
  anasurfname += nh? "NH_" : "IH_";
  anasurfname += p;
  anasurfname += corrSysts? "systs.root": "stats.root";

  TFile *anasurffile = TFile::Open(anasurfname);
  anasurffile->cd();
  std::string plotShortName="delta_th23";
  if(plot == "ssth23dmsq32") plotShortName = "th23_dm32";
  if(nh) plotShortName += "_NH";
  else plotShortName += "_IH";
  TH2* anasurf = (TH2F*)anasurffile->Get(plotShortName.c_str());

  // Grab the FCHelper2018 file
  //std::string helperupsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
  TString helpername =
    helperupsname+"FCInputs/contours_FCInput_2019_joint_realData_"+beam+"_only";
  helpername += nh? "NH_" : "IH_";
  helpername += p;
  helpername += corrSysts? "systs.root":"stats.root";

  TFile *fchelp = TFile::Open(helpername);

  // Find the bin center for this job, probably should be arguments
  // dCP / pi

  // The surfaces are made with ana::ExpandedHistrogram() which sets the bin edges so that
  // the input min/max x and y fall in the center of the bins. Grab these boundaries.
  int NX = anasurf->GetXaxis()->GetNbins();
  int NY = anasurf->GetYaxis()->GetNbins();
  double minX = anasurf->GetXaxis()->GetBinCenter(1);       // root bins start at 1
  double maxX = anasurf->GetXaxis()->GetBinCenter(NX);
  double minY = anasurf->GetYaxis()->GetBinCenter(1);
  double maxY = anasurf->GetYaxis()->GetBinCenter(NY);

  double widthX = anasurf->GetXaxis()->GetBinWidth(1);
  double widthY = anasurf->GetYaxis()->GetBinWidth(1);
  std::cout << "widthX = " << widthX << endl;
  std::cout << "widthY = " << widthY << endl;
  assert(fbin < NX*NY && "Invalid bin number");

  int bin = 0;
  for (int ibin = startidx; ibin <= endidx; ibin++) {
    bin = bins[ibin];
    std::cout<<"- "<<ibin<<std::endl;
    std::cout << "bin " << bin << endl;
    // Once we run out of important bins, get out.
    if(bins[ibin] == 0 && ibin != 0) continue;

    //    bin = bins[i];

    int binY = bin/NX;
    int binX = bin - NX*binY;
    std::cout << "binX: " << binX << endl;
    std::cout << "binY: " << binY << endl;
    anasurffile->cd();
    double centerX = anasurf->GetXaxis()->GetBinCenter(binX+1);
    double centerY = anasurf->GetYaxis()->GetBinCenter(binY+1);
    std::cout <<  "centerX = " << centerX << endl;
    std::cout << "centerY = " << centerY << endl;
    fchelp->cd();

    // Pull out hists and get values at current bin
    double seedssth23 = 0;
    if(plot == "deltath13"){
      TH2* hssth23 = (TH2F*) fchelp->Get((hierStr+"_SinSqTheta23").c_str());
      seedssth23 = hssth23->GetBinContent(binX+1,binY+1);
      delete hssth23;
    }
    double seeddmsq32 = 0;
    if (plot == "deltassth23"){
      TH2* hdmsq = (TH2F*)fchelp->Get((hierStr+"_DmSq32").c_str());
      seeddmsq32 = hdmsq->GetBinContent(binX+1,binY+1); // root bins off-by-one
      delete hdmsq;
    }
    double seeddelta = 0;
    if (plot == "ssth23dmsq32"){
      TH2* hdelta = (TH2F*)fchelp->Get((hierStr+"_DeltaCPpi").c_str());
      seeddelta = hdelta->GetBinContent(binX+1,binY+1);
      delete hdelta;
    }
    TH2* hss2th13 = (TH2F*)fchelp->Get((hierStr+"_SinSq2Theta13").c_str());
    double seedss2th13 = hss2th13->GetBinContent(binX+1,binY+1);
    delete hss2th13;

    // Get current syst shifts
    std::map<std::string,double> seedsyst;
    for (const ISyst* syst : systs){
      TH2* h = (TH2F*)fchelp->Get((hierStr+"_"+syst->ShortName()).c_str());
      if (!h){
        std::cout << "Don't see a prof hist for " << syst->ShortName() << ".  Continuing, but check your ups version." << std::endl;
        continue;
      }
      double seedval = h->GetBinContent(binX+1,binY+1);
      seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
    }

    std::cout << "Parameters from anasurf " << endl;
    std::cout << "ssth23 = " <<  seedssth23 << endl;
    std::cout << "dmsq = " << seeddmsq32 << endl;
    std::cout << "delta = " << seeddelta << endl;
    std::cout << "th13 = " << seedss2th13 << endl;

    // Need random numbers to sample the bin
    TRandom3 rnd(0);

    for (int i = 0; i < NPts; i++){
      ResetOscCalcToDefault(calc);

      const IFitVar* dmsq_var;
      if(nh) dmsq_var = &kFitDmSq32ScaledNH;
      else dmsq_var = &kFitDmSq32ScaledIH;

      // backwards compatability
      double trueX = centerX;
      double trueY = centerY;
      std::cout << "trueX: " << trueX << endl;
      std::cout << "trueY: " << trueY << endl;
      // Plot = deltassth23: x=delta, y=ssth23
      // Plot = ssth23dmsq: x=ssth23, y=dmsq32
      // Plot = deltassth13: x=th13, y=delta
      if(plot != "deltath13"){
        kFitDeltaInPiUnits.SetValue(calc, plot=="deltassth23"?trueX:seeddelta);
        kFitSinSqTheta23.SetValue(calc, plot=="deltassth23"?trueY:trueX);
        dmsq_var->SetValue(calc, plot=="deltassth23"?seeddmsq32:trueY);
        kFitSinSq2Theta13.SetValue(calc, seedss2th13);
      }
      else{
        kFitDeltaInPiUnits.SetValue(calc, trueY);
        kFitSinSqTheta23.SetValue(calc, seedssth23);
        dmsq_var->SetValue(calc, seeddmsq32);
        kFitSinSq2Theta13.SetValue(calc, trueX);
      }

      // Need to seed psuedo-experiments at profiled syst shifts read from
      // FCInput
      SystShifts seednue_fhc;
      for (const ISyst *syst : snue_fhc){
        if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
        seednue_fhc.SetShift(syst,seedsyst[syst->ShortName()]);
      }
      SystShifts seednue_rhc;
      for (const ISyst *syst : snue_rhc) {
        if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
        seednue_rhc.SetShift(syst, seedsyst[syst->ShortName()]);

      }
      SystShifts seednumu;
      for (const ISyst *syst : snumu){
        if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
        seednumu.SetShift(syst,seedsyst[syst->ShortName()]);
      }


      // Make a mock data spectrum for nue
      // Use mock data + cosmic as predictions

      std::cout << "----------Making mock data at-----------------\n";
      std::cout << "ssth23 = " << kFitSinSqTheta23.GetValue(calc);
      std::cout << ", delta = " << kFitDeltaInPiUnits.GetValue(calc);
      std::cout << ", dmsq =  " << dmsq_var->GetValue(calc);
      std::cout << ", th13 = " << kFitSinSq2Theta13.GetValue(calc);
      std::cout << endl;

            // std::vector<Spectrum>                  mockdata;
      // std::vector<Spectrum>                  mockdata_numu;

      // Vestige?
      // SystShifts                             this_seednue;

      // for (auto ipred = 0; ipred < preds.size(); ipred++) {

      //   //   SystShifts tmp_systs;
      // //   if(preds[ipred].flavor == fc::kNue) {
      // //     if(preds[ipred].beam == fc::kFHC) tmp_systs = seednue_fhc;
      // //     else                              tmp_systs = seednue_rhc;
      // //   }
      // //   else                                tmp_systs = seednumu;
      //   std::cout << "___debug ipred " << ipred << std::endl;
      //   Spectrum      fakedata = preds[ipred].pred->Predict    (calc           ); // stats only
      //   // if(corrSysts) fakedata = preds[ipred].pred->PredictSyst(calc, tmp_systs); // with systs
      //   std::cout << "___debug cos first " << preds[ipred].cos.first << std::endl;

      //   fakedata += *preds[ipred].cos.first;
      //   std::cout << "___debug mock pot " << preds[ipred].pot << std::endl;
      //   mockdata.push_back(fakedata.MockData(preds[ipred].pot));
      // }

      std::vector<const IExperiment*> expts;
      std::vector<const IExperiment*> numu_expts;

      int nnumu = 4;
      // auto calc = DefaultOscCalc();
      //if(fake2019) SetFakeCalc(calc, 0.564956, 2.48304e-3, 1.9997);
      // if(fake2019) SetFakeCalc(calc, 0.565, 2.48e-3, 0);
      // else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}

      for(int ipred = 0; ipred < int(preds.size()); ++ipred) {
        double POT = preds[ipred].pot;
        // auto fakedata = GetFakeData (preds[ipred].pred, calc, POT, preds[ipred].cos.first, preds[ipred].livetime);
        auto fakedata = preds[ipred].pred->Predict(calc);
        auto hcos = preds[ipred].cos.first->ToTH1(preds[i].livetime, kBlack, kSolid, kLivetime);
        fakedata += Spectrum(hcos,POT,0);
        // cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc).Integral(POT)<<
        //   " cos "<<hcos->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
        auto mockdata = fakedata.MockData(POT);
        expts.push_back(new SingleSampleExperiment(preds[ipred].pred, mockdata, *preds[ipred].cos.first, preds[ipred].cos.second));
      }

      for(int ipred = 0; ipred < (int) preds_numu_only.size(); ++ipred){
        double POT = preds_numu_only[ipred].pot;
        // auto fakedata = GetFakeData (preds_numu_only[ipred].pred, calc, POT, preds_numu_only[ipred].cos.first);
        auto fakedata = preds_numu_only[ipred].pred->Predict(calc);
        auto hcos = preds[ipred].cos.first->ToTH1(preds[i].livetime, kBlack, kSolid, kLivetime);
        // cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc).Integral(POT)<<
        //   " cos "<<hcos->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
        fakedata += Spectrum(hcos,POT,0);
        auto mockdata = fakedata.MockData(POT);
        numu_expts.push_back(new SingleSampleExperiment(preds_numu_only[ipred].pred, mockdata, *preds_numu_only[ipred].cos.first, preds_numu_only[ipred].cos.second));
      }

      // Set up experiment, nova nue + nova numu + 2017 reactor
      // Get true chi-square before fit

      if (plot != "th13") {
        expts     .push_back( WorldReactorConstraint2019() );
        numu_expts.push_back( WorldReactorConstraint2019() );
      }
      std::cout << "Creating multiexperiment from " << expts.size()
                << " experiments" << std::endl;

      MultiExperiment expt(expts);
      MultiExperiment numu_expt(numu_expts);


      // Needed?
      std::cout<<"Clean verctors"<<std::endl;
      preds.clear();
      preds_numu_only.clear();

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


      // To fix when systs are available
      // MultiExperiment expt(expts);
      // MultiExperiment numu_expt(numu_expts);
      // // Correlate your systs
      // // TODO GetCorrelation function calls are different
      // for(int j = 0 ; j < (int)nue_preds.size(); j++){
      //   if(j==0 && !RHCOnly)
      //     expt.SetSystCorrelations(j, GetCorrelations(true, true));     // nue FHC
      //   else
      //     expt.SetSystCorrelations(j, GetCorrelations(true, false));    // nue RHC
      //   auto notfornumu = GetCorrelations(false, false);              // numu not dependent on beam
      //   for(int j = 0; j < (int)numu_preds.size(); ++j)
      //     expt.SetSystCorrelations(j+nue_preds.size(), notfornumu);
      // }



      /////////////////////////////////////////////////////////////////
      // Numu only fit to figure out dmsq seeds
      // (may be able to do without this for contours)
      /////////////////////////////////////////////////////////////////
      std::cout << "Fitting numu to find pre seeds" << endl;

      double maxmixing = .514;
      osc::IOscCalcAdjustable* calc = DefaultOscCalc();
      // osc::IOscCalcAdjustable* numu_calc = DefaultOscCalc();
      SystShifts auxShifts = SystShifts::Nominal();
      ResetOscCalcToDefault(calc);

      MinuitFitter fitnumu_only(&numu_expt, {&kFitSinSqTheta23, &kFitDmSq32});

      double minchi_numu =
        fitnumu_only.Fit(calc, auxShifts,
                         SeedList({{&kFitSinSqTheta23, {0.4}},
                               {&kFitDmSq32, {nh? 2.35e-3 : -2.35e-3}}}),
                         {},
                         IFitter::kQuiet);

      double pre_seed_dmsq = dmsq_var->GetValue(calc);
      double pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);

      double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
        (2*maxmixing - pre_seed_th23) : pre_seed_th23;
      double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
        pre_seed_th23 : (2*maxmixing - pre_seed_th23);
      if(pre_seed_th23_LO > .5) pre_seed_th23_LO = .45;

      std::cout << "Preseeds ---> ssth23_LO = " << pre_seed_th23_LO;
      std::cout << "\tdmsq_NH = " << fabs(pre_seed_dmsq) << endl;
      // done getting pre seeds
      // --------------------------------------------------------------
      ////////////////////////////////////////////////////////////////////
      // True fit
      ////////////////////////////////////////////////////////////////////

      std::vector<const IFitVar*> trueVars;
      if (plot=="deltassth23"){
        trueVars.push_back(dmsq_var);
        trueVars.push_back(&kFitSinSq2Theta13);
      }
      if (plot=="ssth23dmsq32"){
        trueVars.push_back(&kFitDeltaInPiUnits);
        trueVars.push_back(&kFitSinSq2Theta13);
      }


      MinuitFitter fittrue(&expt, trueVars, systs);
      double tb_th13 = 0;
      double tb_dmsq = 0;
      double tb_ssth23 = 0;
      double tb_delta = 0;
      SystShifts systshift = {};
      // if not correcting deltassth23, then put one element in this
      // vector so the following loop will run
      std::vector <double> tf_delta_seeds = {0, .5, 1., 1.5};
      std::vector <double> tf_ssth23_seeds = {pre_seed_th23_LO,
                                              maxmixing,
                                              pre_seed_th23_UO};

      // make a map of all of the seeds for the true fit
      std::map<const IFitVar*, std::vector<double>> tf_seedPts;
      if(plot == "deltassth23")
        tf_seedPts = {{dmsq_var, {nh? fabs(pre_seed_dmsq) : -1*fabs(pre_seed_dmsq)}},
                      {&kFitSinSq2Theta13, {0.082}}};
      else if(plot == "ssth23dmsq32")
        tf_seedPts = {{&kFitDeltaInPiUnits, tf_delta_seeds},
                      {&kFitSinSq2Theta13, {0.082}}};
      else if(plot == "deltath13")
        tf_seedPts = {{dmsq_var, {nh? fabs(pre_seed_dmsq) : -1*fabs(pre_seed_dmsq)}},
                      {&kFitSinSqTheta23, tf_ssth23_seeds}};

      std::cout << "True fit parameters (pre fit): ssth23 = " << kFitSinSqTheta23.GetValue(calc);
      std::cout << ", delta = " << kFitDeltaInPiUnits.GetValue(calc);
      std::cout << ", ss2th13 = " << kFitSinSq2Theta13.GetValue(calc);
      std::cout << ", dmsq = " << dmsq_var->GetValue(calc) << " x10^-3 eV" << endl;


      std::cout << "Starting true fit --->" << endl;
      auxShifts.ResetToNominal();
      double trueChi = fittrue.Fit(calc,
                                   auxShifts,
                                   tf_seedPts,
                                   {},
                                   IFitter::kQuiet);

      // save best parameters for posterity
      tb_delta = kFitDeltaInPiUnits.GetValue(calc);
      tb_ssth23 = kFitSinSqTheta23.GetValue(calc);
      tb_dmsq = dmsq_var->GetValue(calc);
      tb_th13 = kFitSinSq2Theta13.GetValue(calc);

      std::cout << "True fit parameters (post fit): ssth23 = " << tb_ssth23;
      std::cout << ", delta = " << tb_delta;
      std::cout << ", ss2th13 = " << tb_th13;
      std::cout << ", dmsq = " << tb_dmsq << " x10^-3 eV" << endl;

      // Set up best fitter
      MinuitFitter fit(&expt,
                 {&kFitDeltaInPiUnits, &kFitSinSqTheta23,
                     &kFitDmSq32Scaled, &kFitSinSq2Theta13},
                 systs);

      std::vector <double> bf_dmsq_seeds = {fabs(pre_seed_dmsq),
                                            -1*fabs(pre_seed_dmsq)};
      std::vector <double> bf_delta_seeds = {0., 0.5, 1., 1.5};
      std::vector <double> bf_th23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO,
                                            maxmixing};
      std::vector <double> bf_th13_seeds = {0.082};
      std::map <const IFitVar*, std::vector<double>> bf_seedPts =
        {{&kFitDeltaInPiUnits, bf_delta_seeds},
         {&kFitSinSqTheta23, bf_th23_seeds},
         {&kFitSinSq2Theta13, bf_th13_seeds},
         {&kFitDmSq32Scaled, bf_dmsq_seeds}};

      double bestdelta = 0;
      double bestssth23 = 0;
      double bestth13 = 0;
      double bestdmsq = 0;
      std::cout << "Starting best fit ----> " << endl;
      ResetOscCalcToDefault(calc);
      auxShifts.ResetToNominal();
      double bestChi = fit.Fit(calc,
                               auxShifts,
                               bf_seedPts,
                               {},
                               IFitter::kQuiet);

      bestdelta = kFitDeltaInPiUnits.GetValue(calc);
      bestssth23 = kFitSinSqTheta23.GetValue(calc);
      bestth13 = kFitSinSq2Theta13.GetValue(calc);
      bestdmsq = kFitDmSq32Scaled.GetValue(calc);
      std::cout << "Best fit parameters: ssth23 = " << bestssth23;
      std::cout << ", delta = " << bestdelta;
      std::cout << ", ss2th13 = " << bestth13;
      std::cout << ", dmsq = " << bestdmsq << " x10^-3 eV"<< endl;


      // Make pseudoexperiment dst
      // with FCPoint in its current state (r31771) the names of the branches are
      // not correct. We want to save best fit information
      FCPoint pt(seeddelta, seedssth23, seedss2th13, seeddmsq32,
                 pre_seed_th23, pre_seed_dmsq,
                 tb_delta, tb_ssth23, tb_th13, tb_dmsq,
                 trueChi,
                 bestdelta, bestssth23, bestth13, bestdmsq,
                 bestChi,
                 0, 0);

      fccol.AddPoint(pt);

      unsigned int nfc_pts = fccol.NPoints();
      if(npts_flush != 0) {
        if(nfc_pts % npts_flush == 0) {
          std::cout << "Saving " + std::to_string(nfc_pts) + "th point to " << fcout << endl;
          fccol.SaveToFile(fcout);
        }
      }
    }
  }// end loop over bins
  // Save to file
  fccol.SaveToFile(fcout);
}