NCExtrapolation.cxx

Go to the documentation of this file.

//$Id: NCExtrapolation.cxx,v 1.273 2009/12/18 01:30:21 rodriges Exp $
//
//NCExtrapolation.cxx
//
// Base class for the various extrapolations
//
//B. Rebel 2/2007

#include "NCUtils/Extrapolation/NCExtrapolation.h"
#include "NCUtils/Extrapolation/NCSpectrumInterpolator.h"
#include "NCUtils/NCOscProb.h"
#include "NCUtils/NCUtility.h"

#include "MessageService/MsgService.h"

#include "AnalysisNtuples/ANtpAnalysisInfo.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"

#include "Conventions/Detector.h"

#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TGraph.h"
#include "TObjArray.h"
#include "TLine.h"
#include "TList.h"
#include "TLegend.h"
#include "TMath.h"
#include "TMarker.h"
#include "TROOT.h"
#include "TStopwatch.h"

#include <cassert>
#include <fstream>

using namespace NCType;
using NC::Utility::SQR;

CVSID("$Id: NCExtrapolation.cxx,v 1.273 2009/12/18 01:30:21 rodriges Exp $");

// ---------------------------------------------------------------------
NCExtrapolation::~NCExtrapolation()
{
  if(fBestFitOscPars) delete fBestFitOscPars;

  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it)
    delete it->second;

  // Includes fInterpolator somewhere
  for(InterpCache_t::iterator it = fInterpCache.begin();
      it != fInterpCache.end(); ++it){
    delete it->first;
    delete it->second;
  } // end for it
}


const Registry& NCExtrapolation::DefaultConfig()
{
  static Registry r;

  r.UnLockValues();

  r.Set("EnergyThreshold",           0);

  r.Set("UseCCEvents",               true);
  r.Set("UseNCEvents",               true);
  r.Set("Prediction",                false);
  r.Set("PredictionWriteData",       false);

  r.Set("PrintFDEvents",             false);
  r.Set("HideFDData",                false);

  r.Set("DoSystematicInterpolation", false);
  r.Set("DisableSystPenalty",        false);
  r.Set("AddEventsToEnergyBins",     true);

  r.LockValues();
  return r;
}


//----------------------------------------------------------------------
//this method allows the user to pass in appropriate settings for the
//extrapolation such as locations of files, whether to generate the files
//again, etc.  the registry comes from the job module GetConfig method
void NCExtrapolation::Prepare(const Registry& r)
{
  MSG("NCExtrapolation", Msg::kInfo) << "NCExtrapolation::Prepare(const Registry& r)"
                                     << endl;
  int         tmpb;  // a temp bool.
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.

  if(r.Get("OscillationModel",    tmpi)) fOscillationModel = NCType::EOscModel(tmpi);

  if(r.Get("EnergyThreshold",       tmpd)) fEnergyThreshold = tmpd;

  if(r.Get("UseCCEvents",           tmpb)) fUseCC           = tmpb;
  if(r.Get("UseNCEvents",           tmpb)) fUseNC           = tmpb;

  fCoordConv.Prepare(r);

  if(r.Get("PrintFDEvents", tmpb)) fPrintFDEvents = tmpb;
  if(r.Get("HideFDData",    tmpb)) fHideFDData    = tmpb;

  if(r.Get("DoSystematicInterpolation", tmpb))
    fDoSystematicInterpolation = tmpb;

  if(r.Get("DisableSystPenalty",    tmpb)) fUseSystPenalty = !tmpb;

  if(r.Get("AddEventsToEnergyBins", tmpb)) fAddEventsToEnergyBins = tmpb;
}


//----------------------------------------------------------------------
void NCExtrapolation::AddEvent(NCEventInfo info,
                               bool useMCAsData,
                               NCType::EFileType fileType,
                               NCBeam::Info beamInfo)
{
  //method to add events to near detector beams.  all methods should see
  //the same near detector spectra, so there is no loss of generalization.

  //dont bother with events outside the fiducial volume
  if(info.reco->inFiducialVolume < 1 ||
     (info.header->detector == int(Detector::kNear) &&
      info.reco->isSimpleCutsClean < 1)
     )
     return;


  if(info.reco->nuEnergy < fEnergyThreshold) return;


  //need to check if this type of beam already exists in the vector.
  //if so add the event to that beam. if not, make it the beam then
  //the add event to it
  Detector::Detector_t detector = Detector::Detector_t(info.header->detector);

  NCBeam* beam = GetBeam(detector, beamInfo, true);

  beam->AddEvent(info.header, info.beam, info.reco,
                 info.analysis, info.truth, useMCAsData, fileType);

  // TODO - I don't know what this was achieving, but whatever it was it isn't anymore
  /*
  //if this is MC and not useMCAsData, look to see if there are
  //other run types for the same beam.  default run type for MC is
  //NCType::kRunI
  if(info.header->dataType == int(SimFlag::kMC) && !useMCAsData){
    for(int i = NCType::kRunII; i < NCType::kMaxRun+1; ++i){
      beam = FindBeamIndex(detector, beamType, NCType::ERunType(i));
      if(beam < 100)
        fBeams[beam]->AddEvent(info.header, info.beam, info.reco,
                              info.analysis, info.truth, useMCAsData, fileType);
    }//end loop over possible runs
  }//end if mc and needs to be added to other runs for this beam
  */

  //print out the far detector events if requested
  if(fPrintFDEvents && detector == Detector::kFar && 
     info.header->dataType == int(SimFlag::kData)){

    static bool once=true;
    if(once){
      once=false;
      MsgService::Instance()->GetStream("DataEventListAll")->AttachOStream(Msg::kInfo,"ncutils_events_all.txt");
      
      MsgService::Instance()->GetStream("DataEventListCC")->AttachOStream(Msg::kInfo,"ncutils_events_cc.txt");

      MsgService::Instance()->GetStream("DataEventListNC")->AttachOStream(Msg::kInfo,"ncutils_events_nc.txt");
    }
    TString stream = "";
    TString date = TString::Format("%d-%02d ", info.header->year, info.header->month);
    TString run = TString::Format("F000%d", info.header->run);

    if(info.analysis->isNC < 1 && info.analysis->isCC > 0) stream = "DataEventListCC";
    if(info.analysis->isNC > 0 && info.analysis->isCC < 1) stream = "DataEventListNC";

    if(stream!="")
      MSG(stream, Msg::kInfo)
        << date
        << run << " " << info.header->subRun
        << " " << info.header->snarl
        << " " << info.reco->nuEnergy
        << " " << info.analysis->isCC
        << " " << info.analysis->isNC << endl;

      MSG("DataEventListAll", Msg::kInfo)
        << date
        << run << " " << info.header->subRun
        << " " << info.header->snarl
        << " " << info.reco->nuEnergy
        << " " << info.analysis->isCC
        << " " << info.analysis->isNC << endl;

  }
}


//----------------------------------------------------------------------
void NCExtrapolation::DrawBestFitSpectra(Detector::Detector_t det)
{
  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
    NCBeam* beam = it->second;
    if(beam->GetDetector() != det) continue;

    TString name = "bestFitSpectraCanv"+GetShortName();
    name += Detector::AsString(det);
    name += beam->GetInfo().GetDescription();
    TCanvas *canv = new TCanvas(name, name, 150, 10, 1200, 300);
    TPad *pad1 = new TPad(name+"pad1", "", 0.00, 0.0, 0.50, 1.0, 10);
    TPad *pad2 = new TPad(name+"pad2", "", 0.50, 0.0, 1.00, 1.0, 10);
    pad1->Draw();
    pad2->Draw();
    pad1->cd();

    /*
    TString bestDeltaM;
    bestDeltaM.Format("#Deltam^{2} = %5.4f", BestDeltaMSqr());
    TString bestUMu3Sqr;
    bestUMu3Sqr.Format("sin^{2}(2#theta) = %3.2f", BestUMu3Sqr());
    TString bestUS3Sqr;
    bestUS3Sqr.Format("f_{s} = %3.2f", BestUS3Sqr());
    TString minChiSqr;
    minChiSqr.Format("#chi^{2} = %3.2f", fMinChiSqr);
    */

    for(NCType::EEventType i = NCType::kNC;
        i <= NCType::kCC;
        i = NCType::EEventType(int(i)+1)){
      if(i == NCType::kCC) pad2->cd();

      double maxY = 1.3*TMath::Max(beam->GetMCHistogram(i)->GetMaximum(),
                                   beam->GetDataGraph(i)->GetHistogram()->GetMaximum());

      beam->GetMCHistogram(i)->SetMaximum(maxY);
      TString title = Detector::AsString(det);
      if(i == NCType::kNC)
        title += " Detector Neutral Current Selected";
      else
        title += " Detector Charged Current Selected";
      beam->GetMCHistogram(i)->SetTitle(title);
      beam->GetMCHistogram(i)->Draw();
      //if the fit histogram was filled, draw the fits
      if(beam->GetMCFitHistogram(i)->Integral() > 0.){
        beam->GetMCFitHistogram(i)->Draw("same");
        beam->GetMCBackgroundHistogram(i)->Draw("same");
        beam->GetMCFitNuMuToNuTauHistogram(i)->Draw("same");
        beam->GetMCFitNuEToNuEHistogram(i)->Draw("same");
        beam->GetMCFitNuMuToNuEHistogram(i)->Draw("same");
      }
      else{
        //no fit, no point drawing NuMuToNuTau or NuMuToNuE
        beam->GetMCNoFitBackgroundHistogram(i)->Draw("same");
        beam->GetMCNoFitNuEToNuEHistogram(i)->Draw("same");
      }

      // Don't show the annoying 20-120 GeV bin by default. It's still
      // there in case anyone wants it
      beam->GetMCHistogram(i)->GetXaxis()->SetRangeUser(0, 19.9);

      // This option plots only the ND data when making a prediction.
      // This effectively preserves blinding while being able to make
      // a FD prediction using the actual data.

      if(!fHideFDData || beam->GetDetector() != Detector::kFar){
        beam->GetDataGraph(i)->Draw("pe1same");
      }

      TLegend *leg = new TLegend(0.45, 0.65, 0.85, 0.85);
      leg->SetBorderSize(0);
      leg->AddEntry(beam->GetDataGraph(i), "Data", "lep");
      leg->AddEntry(beam->GetMCHistogram(i),
                    "Monte Carlo (Nominal)", "l");
      if(beam->GetMCFitHistogram(i)->Integral() > 0.){
        leg->AddEntry(beam->GetMCFitHistogram(i),
                      "Monte Carlo (Best Fit)", "l");
        leg->AddEntry(beam->GetMCBackgroundHistogram(i),
                      "Background", "bf");
        leg->AddEntry(beam->GetMCFitNuMuToNuTauHistogram(i),
                      "Background - #nu_{#tau} CC", "bf");
        leg->AddEntry(beam->GetMCFitNuEToNuEHistogram(i),
                      "Background - Beam #nu_{e} CC", "bf");
        leg->AddEntry(beam->GetMCFitNuMuToNuEHistogram(i),
                      "Background - #nu_{e} CC", "bf");
      }
      else{
        leg->AddEntry(beam->GetMCNoFitBackgroundHistogram(i),
                      "Background", "bf");
        leg->AddEntry(beam->GetMCNoFitNuEToNuEHistogram(i),
                      "Background - Beam #nu_{e} CC", "bf");
      }
//       leg->AddEntry(beam->GetDataGraph(i), bestDeltaM, "");
//       leg->AddEntry(beam->GetDataGraph(i), bestUS3Sqr, "");
//       leg->AddEntry(beam->GetDataGraph(i), bestUMu3Sqr, "");
//       leg->AddEntry(beam->GetDataGraph(i), minChiSqr, "");
      leg->Draw();

    }//end loop over NC/CC


    canv->Update();

    // Append the canvas object to the current directory
    // otherwise it can't be found and saved later
    gDirectory->Append(canv);

  }//end loop over beams
}


//----------------------------------------------------------------------
void NCExtrapolation::DrawBestFitRatios(Detector::Detector_t det)
{
  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
    NCBeam* beam = it->second;
    if(beam->GetDetector() != det) continue;

    TString name = "bestFitRatiosCanv"+GetShortName();
    name += Detector::AsString(det);
    name += beam->GetInfo().GetDescription();
    TCanvas *canv = new TCanvas(name, name, 150, 10, 1200, 300);
    TPad *pad1 = new TPad(name+"pad1", "", 0.00, 0.0, 0.50, 1.0, 10);
    TPad *pad2 = new TPad(name+"pad2", "", 0.50, 0.0, 1.00, 1.0, 10);
    pad1->Draw();
    pad2->Draw();
    pad1->cd();

    for(NCType::EEventType i = NCType::kNC;
        i <= NCType::kCC;
        i = NCType::EEventType(int(i)+1)){
      if(i == NCType::kCC) pad2->cd();

      TString title = Detector::AsString(det);
      if(i == NCType::kNC)
        title += " Detector Neutral Current Selected";
      else
        title += " Detector Charged Current Selected";

      TH1* mc_ratio = (TH1*)beam->GetMCFitHistogram(i)->Clone();
      TH1* mc_hist = (TH1*)beam->GetMCHistogram(i)->Clone();
      mc_ratio->Divide(mc_hist);

      mc_ratio->SetTitle(title);
      mc_ratio->GetYaxis()->SetRangeUser(0, 1.5);
      mc_ratio->GetYaxis()->SetTitle("Ratio to no oscillations");
      mc_ratio->Draw("][");

      // This option plots only the ND data when making a prediction.
      // This effectively preserves blinding while being able to make
      // a FD prediction using the actual data.

      TGraphAsymmErrors* data_ratio = 0;

      if(!fHideFDData || beam->GetDetector() != Detector::kFar){
        data_ratio = (TGraphAsymmErrors*)beam->GetDataGraph(i)->Clone();

        for(int n = 0; n < data_ratio->GetN(); ++n){
          double x, y;
          data_ratio->GetPoint(n, x, y);
          double hi = y+data_ratio->GetErrorYhigh(n);
          double lo = y-data_ratio->GetErrorYlow(n);
          double div = mc_hist->GetBinContent(mc_hist->FindBin(x));
          if(div == 0){ // Don't want to see this point...
            div = 1;
            x = -1;
          }
          y /= div; hi /= div; lo /= div;
          data_ratio->SetPoint(n, x, y);
          data_ratio->SetPointEYhigh(n, hi-y);
          data_ratio->SetPointEYlow(n, y-lo);
        }

        data_ratio->Draw("pe1same");
      }

      TLine* one = new TLine(0, 1, mc_ratio->GetXaxis()->GetXmax(), 1);
      one->SetLineStyle(7); // medium dashed
      one->Draw("same");


      TLegend* leg = new TLegend(0.45, 0.35, 0.85, 0.15);
      leg->SetBorderSize(0);
      if(data_ratio) leg->AddEntry(data_ratio, "Data", "l");
      leg->AddEntry(mc_ratio, "Best Fit", "l");

      leg->Draw();

    }//end loop over NC/CC

    canv->Update();

    // Append the canvas object to the current directory
    // otherwise it can't be found and saved later
    gDirectory->Append(canv);

  }//end loop over beams
}


//----------------------------------------------------------------------
void IgnoreErrors(int, Bool_t, const char*, const char*){}


//----------------------------------------------------------------------
void NCExtrapolation::WriteResources(const NC::OscProb::OscPars* /*truePars*/)
{
  MSG("NCExtrapolation", Msg::kInfo)
    << "NCExtrapolation::WriteResources() - "
    << GetShortName() << endl;

  // Make sure all those pretty NCBeam plots get filled
  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it)
    FillResultHistograms(it->second);

  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* saveDir = gDirectory;

  TDirectory* beamsDir = gDirectory->mkdir("beams", "beams");
  beamsDir->cd();

  //write out the spectra etc from the beam objects, while
  //checking on the prediction only flag
  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
    NCBeam* beam = it->second;
    TString dirName = beam->GetInfo().GetDescription();

    // Make the directory dirName, or if it already exists then cd into it.
    // There seems to be no way to do this without trigerring an error, so
    // tell ROOT to ignore them, and then restore the old state afterwards.
    ErrorHandlerFunc_t handler = GetErrorHandler();
    SetErrorHandler(IgnoreErrors);
    if(!beamsDir->cd(dirName)) beamsDir->mkdir(dirName)->cd();
    SetErrorHandler(handler);

    beam->MakeResultPlots();
    if(fHideFDData && beam->GetDetector() == Detector::kFar)
      beam->WritePredictionResources();
    else
      beam->WriteResources();

  }//end loop over beams

  saveDir->cd();

  DrawBestFitSpectra(Detector::kNear);
  DrawBestFitSpectra(Detector::kFar);

  DrawBestFitRatios(Detector::kNear);
  DrawBestFitRatios(Detector::kFar);


  if(fInterpolator){
    TDirectory* saveDir = gDirectory;
    gDirectory->mkdir("interp")->cd();
    fInterpolator->WriteResources();
    saveDir->cd();
  }


  saveDir->cd();
}


//----------------------------------------------------------------------
//add says whether to add a new beam if this one isn't found
NCBeam* NCExtrapolation::GetBeam(Detector::Detector_t det,
                                 NCBeam::Info beamInfo,
                                 bool add)

{
  pair<NCBeam::Info, Detector::Detector_t> key(beamInfo, det);
  fBeams_t::iterator it = fBeams.find(key);
  if(it != fBeams.end()) return it->second;

  if(!add) return 0;

  NCBeam* beam = new NCBeam(det, beamInfo, fUseCC, fUseNC, fAddEventsToEnergyBins);
  fBeams[key] = beam;

  MSG("NCExtrapolation", Msg::kInfo) << "adding new beam "
                                     << beamInfo.GetDescription()
                                     << " to "
                                     << Detector::AsString(det)
                                     << endl;

  return beam;
}


//----------------------------------------------------------------------
//method to fill the result histograms based on the best fit to the
//oscillation parameters and penalty terms
void NCExtrapolation::FillResultHistograms(NCBeam* beam)
{
  NC::OscProb::OscPars* oscPars = GetBestFitOscPars();
  if(!oscPars) return;
  const NC::SystPars systPars = GetBestFitSysts();
  beam->FillResultHistograms(oscPars, &systPars);
}

//----------------------------------------------------------------------
double NCExtrapolation::FindChiSqrForPars(const NC::OscProb::OscPars* oscPars,
                                          const NC::SystPars& systPars,
                                          double* analyticNorm)
{
  // Because profiling reveals we spend all our time in ~TH1 removing ourselves
  // from the current directory. We restore this at the end of the function,
  // need to remember not to take any early returns
  const bool addDirectoryRestore = TH1::AddDirectoryStatus();
  TH1::AddDirectory(false);


  vector<TH1*> exps, obss;

  // For the moment just use all the beams
  std::vector<NCBeam::Info> beamsToUse;
  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
    NCBeam* beam = it->second;
    // Only want to ask for each physical beam once (arbitrarily FD)
    // Don't take the overall combined "kRunAll" beam
    if(beam->GetDetector() == Detector::kFar &&
       beam->GetRunType() != NC::RunUtil::kRunAll &&
       beam->GetInfo().IsNominal())
      beamsToUse.push_back(beam->GetInfo());
  }

  FindSpectraForPars(oscPars, systPars, beamsToUse, exps, obss);


  if(analyticNorm){
    assert(exps.size() == obss.size());
    assert(exps.size() == 1); // TODO, this shouldn't be necessary
    const double M = exps[0]->Integral();
    const double D = obss[0]->Integral();
    //    cout << M << " " << D;
    *analyticNorm = CalculateBestNormalization(M, D);
    //    cout << " -> " << *analyticNorm << endl;
    exps[0]->Scale(*analyticNorm);
  }

  const double ret = LogLikelihood(exps, obss)+
    (fUseSystPenalty ? systPars.PenaltyForShifts() : 0);

  CleanupSpectra(exps, obss);


  TH1::AddDirectory(addDirectoryRestore);

  return ret;
}


//----------------------------------------------------------------------
void NCExtrapolation::SetBestFitCoordNDim(const NC::Fitter::CoordNDim minCoord)
{
  fMinCoord = minCoord;

  if(fBestFitOscPars) delete fBestFitOscPars;

  fBestFitSysts = fCoordConv.SystParsFromCoordNDim(minCoord);
  fBestFitOscPars = fCoordConv.OscParsFromCoordNDim(minCoord);
}


//----------------------------------------------------------------------
void NCExtrapolation::Reset(bool nearData, bool farData,
                            bool nearMC, bool farMC)
{
  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
    NCBeam* beam = it->second;
    if(beam->GetDetector() == Detector::kFar){
      if(farData || farMC) beam->Reset(farData, farMC);
    }
    else{
      if(nearData || nearMC) beam->Reset(nearData, nearMC);
    }
  }
}


//----------------------------------------------------------------------
void NCExtrapolation::CleanupSpectra(std::vector<TH1*> /*exps*/,
                                     std::vector<TH1*> /*obss*/)
{
  // Do nothing. This is for derived classes to, optionally, override
}


//----------------------------------------------------------------------
double NCExtrapolation::LogLikelihood(const vector<TH1*>& exps,
                                      const vector<TH1*>& obss) const
{
  assert(exps.size() == obss.size());

  // With this value, negative expected events and one observed
  // events gives a chisq from this one bin of 182.

  const double minexp = 1e-40; // Don't let expectation go lower than this

  double chisqr = 0;

  for(unsigned int n = 0; n < exps.size(); ++n){
    const int expssize = NC::Utility::GetArraySize(exps[n]);
    const int obsssize = NC::Utility::GetArraySize(obss[n]);

    assert(expssize == obsssize);

    for(int i = 0; i < expssize; ++i){
      double exp = exps[n]->GetBinContent(i);
      if(exp < minexp) exp = minexp;

      const double obs = obss[n]->GetBinContent(i);
      assert(obs >= 0);

      chisqr += 2*(exp-obs);
      if(obs) chisqr += 2*obs*log(obs/exp);
    }
  }

  return chisqr;
}


//----------------------------------------------------------------------
double NCExtrapolation::CalculateBestNormalization(double M, double D) const
{
  // sigma -> infinity
  if(!fUseSystPenalty) return D/M;

  const double sigmaSqr = SQR(NCType::kParams[kNormalization].sigma);

  const double x = (1-sigmaSqr*M)/2;

  return x + sqrt(SQR(x)+D*sigmaSqr);
}


//----------------------------------------------------------------------
vector<double> NCExtrapolation::SystsAsSigmas(const NCBeam::Info& info,
                                              bool* simpleScale) const
{
  // Work out which systematics ever have shifts. Go through all beams and take
  // the combination of all shifts we see.
  bool shiftedVars[NCType::kNumSystematicParameters] = {false,};
  for(fBeams_t::const_iterator it = fBeams.begin(); it != fBeams.end(); ++it){
    NCBeam* beam = it->second;
    for(int n = 0; n < NCType::kNumSystematicParameters; ++n){
      if(beam->GetInfo().IsShifted(NCType::EFitParam(n)))
        shiftedVars[n] = true;
    }
  }

  // Does `shift` represent a systematic that is just a scale factor?
  // This is true for combinations of normalization and cleaning, at least
  // in PID and FarNear. It is possible that for some future extrapolation
  // more or less systematics can be treated this way.
  if(simpleScale) *simpleScale = true;

  // Translation of 'info' into number of sigma shifts in each parameter
  vector<double> shift;
  shift.reserve(NCType::kNumSystematicParameters);
  for(int n = 0; n < NCType::kNumSystematicParameters; ++n){
    if(shiftedVars[n]){
      double val;
      if(!info.IsShifted(NCType::EFitParam(n), &val)) val = 0;
      const double sigmas = val/NCType::kParams[n].sigma;
      // Currently the interpolator can't understand shifts that aren't +/-1
      assert(sigmas == 0 || sigmas == 1 || sigmas == -1);
      shift.push_back(sigmas);

      if(simpleScale &&
         sigmas != 0 &&
         n != NCType::kNormalization) *simpleScale = false; //&&
         //n != NCType::kNCNearClean) *simpleScale = false;
    } // end if shiftedVars
  } // end for n

  return shift;
}


//----------------------------------------------------------------------
bool NCExtrapolation::EnableNearToFarExtrapolation(bool /*enable*/)
{
  assert(0 && "EnableNearToFarExtrapolation needs to be overridden");
}


//----------------------------------------------------------------------
vector<const TH1*> NCExtrapolation::
GetNearMCSpectra(vector<NCBeam::Info> /*beamsToUse*/)
{
  assert(0 && "GetNearMCSpectra needs to be overridden");
}


//----------------------------------------------------------------------
void NCExtrapolation::InitializeInterpolator(const NC::OscProb::OscPars* osc)
{
  assert(fDoSystematicInterpolation);

  using namespace NC::Utility;

  // Ignore calls that come as a result of running this function
  static bool already = false;
  if(already) return;
  already = true;

  // If we already have an interpolator cached for these oscillation
  // parameters then use just that.
  if(fInterpCache.find(osc) != fInterpCache.end()){
    fInterpolator = fInterpCache[osc];
    already = false;
    return;
  }

  // Stop the cache getting too huge (interpolators hold lots of histograms)
  if(fInterpCache.size() > 100){
    const NC::OscProb::OscPars* oldKey = 0;
    for(InterpCache_t::iterator it = fInterpCache.begin();
        it != fInterpCache.end(); ++it){
      if(fInterpolator && it->second == fInterpolator) oldKey = it->first;
      if(it->second != fInterpolator){
        delete it->first;
        delete it->second;
      }
    } // end for it
    fInterpCache.clear();
    if(oldKey) fInterpCache[oldKey] = fInterpolator;
  }

  // We don't want any of the FindSpectraForPars to use the interpolator
  NC::ISpectrumInterpolator* oldInterp = fInterpolator;
  fInterpolator = 0;

  // Disable extrapolation so that we get independently shifted spectrum in FD
  const bool farNearRestore = EnableNearToFarExtrapolation(false);

  // Build a list of beam infos describing each separate running condition
  vector<NCBeam::Info> allBeams;
  for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
    const pair<NCBeam::Info, Detector::Detector_t> key = it->first;
    const NCBeam::Info info = key.first;
    const Detector::Detector_t det = key.second;

    // RunAll isn't a "real" beam
    if(info.GetRunType() == NC::RunUtil::kRunAll) continue;
    // Arbitrary - just need to get each beam description once
    if(det != Detector::kFar) continue;

    allBeams.push_back(info);
  }

  // beamsToUseNom = [b for b in allBeams if b.IsNominal()]
  vector<NCBeam::Info> beamsToUseNom;
  for(unsigned int n = 0; n < allBeams.size(); ++n){
    if(allBeams[n].IsNominal()) beamsToUseNom.push_back(allBeams[n]);
  }

  vector<TH1*> junk;    // We don't need to get data spectra
  vector<TH1*> noms_fd; // Nominal FD spectra
  FindSpectraForPars(osc, NC::SystPars(), beamsToUseNom, noms_fd, junk);

  vector<const TH1*> noms_nd = GetNearMCSpectra(beamsToUseNom);

  // The index order of the spectra is:
  // (Run I FD) (Run I ND) (Run II FD) (Run II ND) etc.
  vector<const TH1*> noms;
  assert(noms_nd.size() == noms_fd.size());
  noms.reserve(2*noms_fd.size());
  for(unsigned int n = 0; n < noms_fd.size(); ++n){
    noms.push_back(CloneFast(noms_fd[n]));
    noms.push_back(CloneFast(noms_nd[n]));
  }

  CleanupSpectra(noms_fd, junk);

  // Don't assign to fInterpolator yet, as we still have some more calls
  // to FindSpectraForPars to make that don't require interpolation.
  NC::ISpectrumInterpolator* interp = new NC::SpectrumInterpolatorSimplest;

  // If an entry in this map is filled, then this shift has already been
  // dealt with for all runs by the simpleScale optimization below.
  map<vector<double>, bool> alreadyFilled;

  // We accumulate expected spectra binned by what shifts they have.
  typedef map<vector<double>, vector<const TH1*> > ExpsMap;
  ExpsMap expsMap;

  for(unsigned int beamIdx = 0; beamIdx < allBeams.size(); ++beamIdx){
    const NCBeam::Info info = allBeams[beamIdx];

    bool simpleScale;
    const vector<double> shift = SystsAsSigmas(info, &simpleScale);

    // The simpleScale optimization already fired and there's no need to add
    // any more ratios for this shift setting
    if(alreadyFilled[shift]) continue;

    if(simpleScale && oldInterp){
      // These are the ratios for this shift for _all_ beams
      vector<TH1*> ratios = oldInterp->GetInputSpectra(shift);
      // So set the flag so that we don't do this again and double-count
      alreadyFilled[shift] = true;
      interp->AddInputSpectra(shift, ratios);
      continue;
    } // end if

    vector<NCBeam::Info> beamsToUse(1, info);
    vector<TH1*> exps_fd;
    FindSpectraForPars(osc, NC::SystPars(), beamsToUse, exps_fd, junk);

    vector<const TH1*> exps_nd = GetNearMCSpectra(beamsToUse);

    vector<const TH1*> exps; // shifted expectations
    assert(exps_fd.size() == exps_nd.size());
    exps.reserve(2*exps_fd.size());
    // The spectra go in the vector matching noms, indexed like
    // (Run I FD) (Run I ND) (Run II FD) (Run II ND)
    for(unsigned int n = 0; n < exps_fd.size(); ++n){
      TH1* exp_fd = CloneFast(exps_fd[n]);
      exp_fd->SetName((info.GetDescription()+" FD").Data());
      exps.push_back(exp_fd);

      TH1* exp_nd = CloneFast(exps_nd[n]);
      exp_nd->SetName((info.GetDescription()+" ND").Data());
      exps.push_back(exp_nd);
    }

    CleanupSpectra(exps_fd, junk);

    // expsMap[shift].append(exps)
    vector<const TH1*>& already = expsMap[shift];
    already.insert(already.end(), exps.begin(), exps.end());
  } // end for beamIdx

  for(ExpsMap::iterator it = expsMap.begin(); it != expsMap.end(); ++it){
    vector<double> shift = it->first;
    vector<const TH1*> exps = it->second;
    // The shifted spectra are stored in the interpolator as ratios over the
    // nominal spectra.
    vector<TH1*> ratios;
    ratios.reserve(exps.size());
    assert(exps.size() == noms.size());
    for(unsigned int n = 0; n < exps.size(); ++n){
      ratios.push_back(CloneFast(exps[n]));
      DivideFast(ratios[n], noms[n]);
    }
    interp->AddInputSpectra(shift, ratios);
  } // end for it

  // All of these were Clone'd and need deleting
  for(unsigned int n = 0; n < noms.size(); ++n) delete noms[n];

  // Put extrapolation setting back to however it was.
  EnableNearToFarExtrapolation(farNearRestore);

  // Have to wait till the last moment, because if fSpectrumInterpolator exists
  // then FindSpectraForPars tries to use it -> chicken and egg problem.
  fInterpolator = interp;

  assert(fInterpCache.find(osc) == fInterpCache.end());
  fInterpCache[const_cast<NC::OscProb::OscPars*>(osc)->Clone()] = fInterpolator;

  already = false;
}

---