make_fc_nus_surfs_nersc_2019.C

Go to the documentation of this file.

#include "TFile.h"
#include "TH2.h"
#include "TRandom3.h"
#include "TDirectory.h"
#include "TH1.h"

#include <ctime>
#include <fstream>
#include <sstream>

#include "StandardRecord/Proxy/SRProxy.h"



#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/Sample.h"

#include "CAFAna/Prediction/CovarianceMatrix.h"

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

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

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

#include "NuXAna/macros/Nus19FHC/Utilities.h"

using namespace ana;

std::vector<double> LogBins(int nbins,double xlo, double xhi)
{
  double log_spacing = exp( (log(xhi) - log(xlo))/(nbins) );
  std::vector<double> binning;
  for (int i = 0; i <= nbins; ++i) {
    if (i == 0) binning.push_back(xlo);
    else        binning.push_back(log_spacing*binning[i-1]);
  }
  return binning;
}

std::vector<double> LinBins(int nbins,double xlo, double xhi)
{
  double lin_spacing = (xhi - xlo)/(nbins);
  std::vector<double> binning;
  for (int i = 0; i <= nbins; ++i) {
    if (i == 0) binning.push_back(xlo);
    else        binning.push_back(lin_spacing+binning[i-1]);
  }
  return binning;
}

void make_fc_nus_surfs_nersc_2019(int NPts = 1, int startidx = 0, int endidx = 0, bool nh = true, 
         int N = 1, int jid = 1, int npts_flush = 1, 
         std::string plot = "th24vsth34_UO" ,
         TString beam = "neardetfardetFHCcovmx", bool corrSysts = true,
         bool readFromFile = true, double max_time_in_sec = 3600)
{
  //clock starts
  std::clock_t tstart = std::clock();
  
  DontAddDirectory guard;
  
  FCCollection fccol;

  std::vector<std::string> polarity = { "fhc", "rhc" };
  std::vector<bool> use_polarity = GetOptionConfig(beam, polarity);
  std::vector<std::string> detector = { "neardet", "fardet" };
  std::vector<bool> use_detector = GetOptionConfig(beam, detector);

  bool covmx = beam.Contains("covmx", TString::kIgnoreCase);
  bool extrap = beam.Contains("extrap", TString::kIgnoreCase);
  assert(covmx xor extrap);

  if (covmx && use_polarity[1])
    throw std::runtime_error("RHC covariance matrix not supported yet.");

  if (extrap && use_polarity[0])
    throw std::runtime_error("Extrapolated fit not supported for FHC.");

  std::cout << "beam:: " << beam << std::endl;
  std::cout << "systs:: " << corrSysts << std::endl;
  if (covmx && corrSysts) {
    corrSysts = false;
    std::cout << "Using covariance matrix, so disabling traditional systematics." << std::endl;
  }

  bool upperOctant = false;
  if(plot == "th24vsth34_UO" || 
     plot == "th24_UO" || 
     plot == "th34_UO" ||
     plot == "dm41vsth24_UO" ||
     plot == "dm41vsth34_UO")
    {
       upperOctant = true;
    }

  bool th24 = false;
  bool th34 = false;
  bool th24vsth34 = false;
  bool dm41vsth24 = false;
  bool dm41vsth34 = false;
  if(plot == "th24_UO" || plot == "th24_LO"){
    th24 = true;
  }
  else if(plot == "th34_UO" || plot == "th34_LO"){
    th34 = true;
  }
  else if(plot == "th24vsth34_UO" || plot == "th24vsth34_LO") th24vsth34 = true;
  else if(plot == "dm41vsth24_UO" || plot == "dm41vsth24_LO") dm41vsth24 = true;
  else if(plot == "dm41vsth34_UO" || plot == "dm41vsth34_LO") dm41vsth34 = true;

  beam += "extrap";
  if(corrSysts) beam += "syst";

  // Get POT & livetimes
  //double nd_pot = use_polarity[1] ? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
  double fd_pot = use_polarity[1] ? kAna2018RHCPOT : kAna2018FHCPOT;
  double fd_livetime = use_polarity[1] ? kAna2018RHCLivetime : kAna2018FHCLivetime;

  IPrediction* pred = nullptr;
  Spectrum* cosmic_ext = nullptr;

  std::vector<covmx::Sample> samples;
  std::vector<IPrediction*> preds;
  std::vector<Spectrum*> cosmic;
  std::vector<Spectrum*> data;
  std::vector<Spectrum*>::iterator iter;
  
  CovarianceMatrix* mx = nullptr;

  std::string pull_syst = "none";
  std::vector<const ISyst*> systs = {};
  SystConcat* sc = nullptr;
  
  // OscCalc for Sterile 3+1 Model
  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
  SetAngles(calc);
  
  if (covmx) {

    for (size_t d = 0; d < 2; ++d) {   // for detector
      for (size_t p = 0; p < 2; ++p) { // for beam polarity
        if (!use_polarity[p] || !use_detector[d]) continue;
        samples.push_back(covmx::Sample(covmx::Analysis::kNC,
          (covmx::Polarity)p, (covmx::Detector)d));
        preds.push_back(LoadPrediction(samples.back(), pull_syst == "none" ? "nosysts" : "systs"));
      } // for beam polarity
    } // for detector

    for (auto sample : samples)
      cosmic.push_back(LoadCosmic(sample).release());
 
    std::vector<const ISyst*> isysts = GetSystsFromFile(samples[0]);
 
    if (pull_syst == "all") {
      // All pull systs, single detector
      if (samples.size() == 1)
        systs = isysts;
      // All pull systs, both detectors
      else {
        sc = new SystConcat(GetSystConcat(samples, isysts, preds));
        systs = sc->GetActiveSysts();
      }
    }
    else if (pull_syst != "none") {
      // Single pull syst, single detector
      if (samples.size() == 1) {
        for (const ISyst* syst : isysts) {
          if (GetSystBaseName(syst->ShortName()) == pull_syst) {
            systs = { syst };
            break;
          }
        }
        if (systs.empty()) {
          std::ostringstream err;
          err << "Couldn't find requested pull systematic " << pull_syst;
          assert (false and err.str().c_str());
        }
      }
      // Single pull syst, both detectors
      else {
        std::vector<SystConcat> all_scs = GetSystConcats(samples, preds, isysts);
        for (auto it_sc : all_scs) {
          if (GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pull_syst) {
            sc = new SystConcat(it_sc);
            systs = sc->GetActiveSysts();
            break;
          }
        }
      }
    }
  
    for (size_t i = 0; i < preds.size(); ++i) {
      data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(12e20)));
      if(samples[i].detector == covmx::Detector::kFarDet)
        data.back()->OverrideLivetime(kNus19SensLivetime);
      if (cosmic[i]) *data[i] += *cosmic[i];
    }

    std::vector<Binning> bins;
    for (auto d : data) bins.push_back(d->GetBinnings()[0]);
    mx = new CovarianceMatrix(bins);
    std::cout << std::endl;
    for (CovarianceMatrix* it_mx : GetMatrices(samples)) {
      if (GetSystBaseName(it_mx->GetName()) != pull_syst) {
        std::cout << "Adding " << it_mx->GetName()
          << " to covariance matrix." << std::endl;
        mx->AddMatrix(it_mx);
      }
      else
        std::cout << "Omitting pull term " << pull_syst
          << " from matrix." << std::endl;
    }

  } // if covmx
  else {

    // load prediction for Extrap fit
    pred = LoadExtrapPrediction(beam);
      
    // load cosmics for Extrap fit
    //cosmic = new Spectrum(LoadCosmicSpectrum(rhc));
    cosmic_ext = new Spectrum(LoadCosmicSpectrum(polarity[1]));
    cosmic_ext->Scale(fd_livetime/cosmic_ext->Livetime());
    cosmic_ext->OverridePOT(fd_pot);

    if (corrSysts)
      //systs = rhc ? GetNus18RHCExtrapSysts() : GetNus18FHCFDSysts("all");
      systs = LoadISysts();

  }

  std::cout << "Number of Systs:: " << systs.size() << std::endl; 
 
  // Binnings
  std::vector<double> th24Bins;
  std::vector<double> th34Bins;
  std::vector<double> dm2Bins;

  if(th24){
    th24Bins = LinBins(180, 0., 45.);
  }
  else if(th34) {
    th34Bins = LinBins(360, 0., 90.);
  }
  else if (th24vsth34) {
    th24Bins = LinBins(45, 0., 45.);
    th34Bins = LinBins(90, 0., 90.);
  }
  else if (dm41vsth24) {
    th24Bins = LogBins(50, 1.0e-6, 1);
    dm2Bins  = LogBins(50, 1.0e-2, 100);
  }
  else if (dm41vsth34) {
    th34Bins = LogBins(50, 1.0e-6, 1);
    dm2Bins  = LogBins(50, 1.0e-2, 100);
  }

  double binxlowedge = 0.;
  double binxuppedge = 0.; 
  double binylowedge = 0.;
  double binyuppedge = 0.;
  double binxcenter = 0.;
  double binycenter = 0.;

  int TotBins = 4050;
  int bins[TotBins];

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

  std::string binname = "_bins_"+std::to_string(ibin) + "_to_" + std::to_string(fbin);
  std::string hierStr = nh ? "NH" : "IH";
  std::string number = std::to_string(N);
  std::string fcout = "FCCol_nus18_"+plot+"_"+hierStr+"_";
  fcout += use_polarity[0] ? "FHC":"RHC";
  fcout += corrSysts ? "_systs":"_stats";
  fcout += binname+"_v"+number+'_'+std::to_string(jid)+".root";

  // finished loading
  std::clock_t tload = std::clock();
  double loadtime = ( tload - tstart ) / (double) CLOCKS_PER_SEC;  

  TRandom3 rnd(0); // zero seeds randomly  

  for (int i = startidx; i <= endidx; i++) {
    
    int iter_bin = bins[i];

    if(th24){
      binxlowedge = th34Bins.front();
      binxuppedge = th34Bins.back(); 
      binylowedge = th24Bins[iter_bin]; 
      binyuppedge = th24Bins[iter_bin+1];
      binxcenter = 0.5*(binxlowedge+binxuppedge);
      binycenter = 0.5*(binylowedge+binyuppedge);
    }
    else if(th34){
      binxlowedge = th34Bins[iter_bin];
      binxuppedge = th34Bins[iter_bin+1]; 
      binylowedge = th24Bins.front(); 
      binyuppedge = th24Bins.back();
      binxcenter = 0.5*(binxlowedge+binxuppedge);
      binycenter = 0.5*(binylowedge+binyuppedge);
    }
    else if (th24vsth34) {
      int x_bins = th34Bins.size() - 1;
      int x_idx = iter_bin % x_bins;
      int y_idx = iter_bin / x_bins;
      binxlowedge = th34Bins[x_idx];
      binxuppedge = th34Bins[x_idx+1]; 
      binylowedge = th24Bins[y_idx]; 
      binyuppedge = th24Bins[y_idx+1];
      binxcenter = 0.5*(binxlowedge+binxuppedge);
      binycenter = 0.5*(binylowedge+binyuppedge);
    }
    else if (dm41vsth24) {
      int x_bins = th24Bins.size() - 1;
      int x_idx = iter_bin % x_bins;
      int y_idx = iter_bin / x_bins;
      binxlowedge = th24Bins[x_idx];
      binxuppedge = th24Bins[x_idx+1]; 
      binylowedge = dm2Bins[y_idx]; 
      binyuppedge = dm2Bins[y_idx+1];
      binxcenter = 0.5*(binxlowedge+binxuppedge);
      binycenter = 0.5*(binylowedge+binyuppedge);
    }
    else if (dm41vsth34) {
      int x_bins = th34Bins.size() - 1;
      int x_idx = iter_bin % x_bins;
      int y_idx = iter_bin / x_bins;
      binxlowedge = th34Bins[x_idx];
      binxuppedge = th34Bins[x_idx+1]; 
      binylowedge = dm2Bins[y_idx]; 
      binyuppedge = dm2Bins[y_idx+1];
      binxcenter = 0.5*(binxlowedge+binxuppedge);
      binycenter = 0.5*(binylowedge+binyuppedge);
    }

    std::cout << "Bin X (low-center-up): " << binxlowedge << "-" << binxcenter << "-" << binxuppedge << std::endl;
    std::cout << "Bin Y (low-center-up): " << binylowedge << "-" << binycenter << "-" << binyuppedge << std::endl;
    
    double profDelta13_PiUnits = 1.46;
    double profssTh13 = pow( sin( kNus19Th13 ) , 2);
    double profssTh23 = 0.470;
    double profDmSq = kNus19Dm21 + kNus19Dm32;
    double profDelta24_PiUnits = 1.;
    if(upperOctant){
      profDelta13_PiUnits = kNus19Delta13/M_PI;
      profssTh23 = pow( sin( kNus19Th23 ) , 2);
    }
    double profTh24 = 0.;
    double profTh34 = 0.;
    double profDm41 = 0.;  
    if(th24 || th34 || th24vsth34) {
      profTh24 = binycenter;
      profTh34 = binxcenter;
      profDm41 = 0.5;
    }
    else if(dm41vsth24) {
      profTh24 = binxcenter;
      profTh34 = 0.5;
      profDm41 = binycenter;
    }
    else if(dm41vsth34) {
      profTh24 = 0.5;
      profTh34 = binxcenter;
      profDm41 = binycenter;
    }

    // Vector of parameters to fit: #theta_{23}, #theta_{34}, #theta_{24}, #deltam_{32}
    std::vector<const IFitVar*> prof_vars;
    prof_vars.push_back(&kFitSinSqTheta23Sterile); // to be constrained to PDG
    prof_vars.push_back(&kFitDmSq32Sterile); // to be constrained to PDG
    prof_vars.push_back(&kFitDelta24InPiUnitsSterile);
    
    if(th24) {
      prof_vars.push_back(&kFitTheta34InDegreesSterile);
    }
    else if(th34) {
      prof_vars.push_back(&kFitTheta24InDegreesSterile);
    }
    else if(dm41vsth24) {
      prof_vars.push_back(&kFitSinSqTheta34Sterile);   
    }
    else if(dm41vsth34) {
      prof_vars.push_back(&kFitSinSqTheta24Sterile);   
    }
    
    std::vector<const IFitVar*> fitvars = prof_vars;

    if(th24) {
      fitvars.push_back(&kFitTheta24InDegreesSterile);
    }
    else if(th34) {
      fitvars.push_back(&kFitTheta34InDegreesSterile);
    }
    else if(th24vsth34) {
      fitvars.push_back(&kFitTheta24InDegreesSterile);
      fitvars.push_back(&kFitTheta34InDegreesSterile);
    }
    else if(dm41vsth24) {
      fitvars.push_back(&kFitSinSqTheta24Sterile);   
      fitvars.push_back(&kFitDmSq41Sterile);   
    }
    else if(dm41vsth34) {
      fitvars.push_back(&kFitSinSqTheta34Sterile);   
      fitvars.push_back(&kFitDmSq41Sterile);   
    }
    
    for(int iTrial = 0; iTrial < NPts; ++iTrial){
      
      std::clock_t tistart = std::clock();
      if (max_time_in_sec > 0.0){
        double elapstime = ( tistart - tstart ) / (double) CLOCKS_PER_SEC;
        if (elapstime >= max_time_in_sec) break;
      }
      
      SetAngles(calc);

      std::cout << "On trial: " << iTrial << std::endl;

      auto trueX = rnd.Uniform(binxlowedge, binxuppedge);
      auto trueY = rnd.Uniform(binylowedge, binyuppedge);
      auto trueDelta24 = rnd.Uniform(0.,2.);

      double trueTh24 = 0.;    
      double trueTh34 = 0.;    
      double trueDm41 = 0.;

      int mock_seed = 0;
      std::cout << "Mock Seed: " << mock_seed << std::endl;

      if(th24 || th34 || th24vsth34) {
        trueTh24 = trueY;
        trueTh34 = trueX;
        trueDm41 = 0.5;
        kFitTheta34InDegreesSterile.SetValue(calc, trueTh34);
        kFitTheta24InDegreesSterile.SetValue(calc, trueTh24);
      }
      else if(dm41vsth24) {
        trueDm41 = trueY;
        trueTh24 = trueX;
        trueTh34 = rnd.Uniform(0.,1.);
        kFitSinSqTheta24Sterile.SetValue(calc, trueTh24);  
        kFitSinSqTheta34Sterile.SetValue(calc, trueTh34);  
      }
      else if(dm41vsth34) {
        trueDm41 = trueY;
        trueTh34 = trueX;
        trueTh24 = rnd.Uniform(0.,1.);
        kFitSinSqTheta24Sterile.SetValue(calc, trueTh24);  
        kFitSinSqTheta34Sterile.SetValue(calc, trueTh34);  
      }
      kFitDelta24InPiUnitsSterile.SetValue(calc, trueDelta24);    
      kFitDmSq41Sterile.SetValue(calc, trueDm41);  

      // Generate Mock Data
      Spectrum* mock_spec = nullptr;
      
      for (iter=data.begin(); iter != data.end(); ++iter){
        delete *iter;
      }
      data.clear();

      // For CovMx Fit
      if (covmx) {

        for (size_t i = 0; i < preds.size(); ++i) {
          data.push_back(new Spectrum(preds[i]->Predict(calc).MockData(12e20)));
          if(samples[i].detector == covmx::Detector::kFarDet)
            data.back()->OverrideLivetime(kNus19SensLivetime);
          if (cosmic[i]) *data[i] += *cosmic[i];
        }

      } // if covmx
      // For Extrap Fit
      else {

        Spectrum prediction = pred->Predict(calc);
        cosmic_ext->Scale(fd_livetime/cosmic_ext->Livetime());
        cosmic_ext->OverridePOT(fd_pot);  
        Spectrum cos_spec = cosmic_ext->MockData(fd_pot);
        mock_spec = new Spectrum(prediction.MockData(fd_pot) + cos_spec);

      } // if extrap
       
      const IExperiment* expt;
      if (covmx) {
        if(corrSysts) expt = new OscCovMxExperiment(preds, mx, data, cosmic);
        else expt = new OscCovMxExperiment(preds, data, cosmic);
      }
      else {
        expt = new SingleSampleExperiment(pred, *mock_spec, *cosmic_ext);
      }
      
      std::vector<const IExperiment*> expts = GetConstraints();
      expts.push_back(expt);
      MultiExperiment multi_expt(expts);

      // -------- TRUE FIT -------- //
      MinuitFitter fittrue(&multi_expt, prof_vars, systs);
      
      //double trueChi = fittrue.Fit(calc)->EvalMetricVal();
      double trueChi = fittrue.Fit(calc, IFitter::kQuiet)->EvalMetricVal();
      double trueDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc);
      double truessTh13 = kFitSinSqTheta13Sterile.GetValue(calc);
      double truessTh23 = kFitSinSqTheta23Sterile.GetValue(calc);
      double trueDmSq = kFitDmSq32Sterile.GetValue(calc);
      trueDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc);
    
      if(th24) {
        trueTh34 = kFitTheta34InDegreesSterile.GetValue(calc);
      }
      else if(th34) {
        trueTh24 = kFitTheta24InDegreesSterile.GetValue(calc);
      }
      else if(dm41vsth24) {
        trueTh34 = kFitSinSqTheta34Sterile.GetValue(calc);
      }
      else if(dm41vsth34) {
        trueTh24 = kFitSinSqTheta24Sterile.GetValue(calc);
      }
    
      // -------- BEST FIT -------- //
      MinuitFitter fit(&multi_expt, fitvars, systs);
    
      double bestChi = 1E20;
      double bestDelta13 = 0;
      double bestssTh13 = 0;
      double bestssTh23 = 0;
      double bestDmSq = 0;
      double bestDelta24 = 0;
      double bestTh24 = 0;
      double bestTh34 = 0;
      double bestDm41 = 0;

      if(th24 || th34 || th24vsth34) {
        for (double th34_seed: {15, 30, 45, 60, 75}) {
          for (double th24_seed: {9, 18, 27, 36}) {
            std::vector<double> th34_seed_vec; th34_seed_vec.push_back(th34_seed);
            std::vector<double> th24_seed_vec; th24_seed_vec.push_back(th24_seed);
            std::map<const IFitVar*, std::vector<double> > seedPts;
            seedPts[&kFitTheta34InDegreesSterile] = th34_seed_vec;
            seedPts[&kFitTheta24InDegreesSterile] = th24_seed_vec;

            SystShifts auxShifts = SystShifts::Nominal();
            std::vector<SystShifts> seedShifts = {};
            
            //auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts)->EvalMetricVal();
            auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts, IFitter::kQuiet)->EvalMetricVal();
            
            if (bestChi_mtplSeed < bestChi) {
              bestChi = bestChi_mtplSeed;
              bestDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc);
              bestssTh13 = kFitSinSqTheta13Sterile.GetValue(calc);
              bestssTh23 = kFitSinSqTheta23Sterile.GetValue(calc);
              bestDmSq = kFitDmSq32Sterile.GetValue(calc);
              bestDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc);
              bestTh24 = kFitTheta24InDegreesSterile.GetValue(calc);
              bestTh34 = kFitTheta34InDegreesSterile.GetValue(calc);
              bestDm41 = kFitDmSq41Sterile.GetValue(calc);
            }
          
            SetAngles(calc);

          }
        }
      }
      
      else if(dm41vsth24 || dm41vsth34) {
        for (double th_seed: {1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1}) {
          for (double dm_seed: {5.0e-2, 5.0e-1, 5.0, 5.0e+1}) {
            std::vector<double> th_seed_vec; th_seed_vec.push_back(th_seed);
            std::vector<double> dm_seed_vec; dm_seed_vec.push_back(dm_seed);
            std::map<const IFitVar*, std::vector<double> > seedPts;
            if(dm41vsth24) {
              seedPts[&kFitSinSqTheta24Sterile] = th_seed_vec;
              seedPts[&kFitDmSq41Sterile] = dm_seed_vec;
              seedPts[&kFitSinSqTheta34Sterile] = {0.5};
            }
            else if(dm41vsth34) {
              seedPts[&kFitSinSqTheta34Sterile] = th_seed_vec;
              seedPts[&kFitDmSq41Sterile] = dm_seed_vec;
              seedPts[&kFitSinSqTheta24Sterile] = {0.5};
            }

            SystShifts auxShifts = SystShifts::Nominal();
            std::vector<SystShifts> seedShifts = {};
      
            //auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts)->EvalMetricVal();
            auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts, IFitter::kQuiet)->EvalMetricVal();
            
            if (bestChi_mtplSeed < bestChi) {
              bestChi = bestChi_mtplSeed;
              bestDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc);
              bestssTh13 = kFitSinSqTheta13Sterile.GetValue(calc);
              bestssTh23 = kFitSinSqTheta23Sterile.GetValue(calc);
              bestDmSq = kFitDmSq32Sterile.GetValue(calc);
              bestDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc);
              bestTh24 = kFitSinSqTheta24Sterile.GetValue(calc);
              bestTh34 = kFitSinSqTheta34Sterile.GetValue(calc);
              bestDm41 = kFitDmSq41Sterile.GetValue(calc);
            }
    
            SetAngles(calc);

          }
        }
      }
 
      std::clock_t tiend = std::clock();
      double calctime = ( tiend - tistart ) / (double) CLOCKS_PER_SEC;

      FCPoint pt(profDelta13_PiUnits, profssTh23, profssTh13, profDmSq, 
                 profDelta24_PiUnits, profTh24, profTh34, profDm41, 
                 0., 0., 
                 trueDelta13, truessTh23, truessTh13, trueDmSq, 
                 trueDelta24, trueTh24, trueTh34, trueDm41,
                 trueChi, 
                 bestDelta13, bestssTh23, bestssTh13, bestDmSq,
                 bestDelta24, bestTh24, bestTh34, bestDm41, 
                 bestChi,
                 trueChi-bestChi, loadtime, calctime);
     
      double bestX = 0;
      double bestY = 0;

      if(th24 || th34 || th24vsth34) {
        bestY = bestTh24;
        bestX = bestTh34;
      }
      else if(dm41vsth24) {
        bestY = bestDm41;
        bestX = bestTh24;
      }
      else if(dm41vsth34) {
        bestY = bestDm41;
        bestX = bestTh34;
      }

      std::cout << "Trial: " << iTrial << ", TrueX: " << trueX << ", Y: " << trueY << ", trueChi: " << trueChi << ", bestX: " << bestX << ", Y: " << bestY << ", bestChi: " << bestChi << ", true-best: " << trueChi-bestChi << ", loadTime: " << loadtime << ", calcTime: " << calctime << std::endl;

      fccol.AddPoint(pt);

      delete expt;
      delete mock_spec;
    }
  }
  fccol.SaveToFile(fcout);
}