NCExtrapolationFarNear.cxx

Go to the documentation of this file.

//
// NCExtrapolationFarNear.cxx
//
//
// J. Koskinen 6/2007
//
// $Id: NCExtrapolationFarNear.cxx,v 1.48 2009/11/26 12:03:29 bckhouse Exp $

#include "NCUtils/Extrapolation/NCExtrapolationFarNear.h"
#include "NCUtils/Extrapolation/NCEnergyBin.h"

#include "MessageService/MsgService.h"

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

#include "TCanvas.h"
#include "TDirectory.h"
#include "TGraphErrors.h"
#include "TMath.h"

#include <fstream>
#include <cassert>
#include <cmath>

CVSID("$Id: NCExtrapolationFarNear.cxx,v 1.48 2009/11/26 12:03:29 bckhouse Exp $");

#include "NCExtrapolationModule.h"
#include "NCSpectrumInterpolator.h"

REGISTER_NCEXTRAPOLATION(NCExtrapolationFarNear, FarNear)


static const int kNCSig     = 0; //NC spectrum
static const int kCCSig     = 1; //CC spectrum
static const int kNCBkg     = 2; //CC background in NC spectrum
static const int kCCBkg     = 3; //NC background in CC spectrum
static const int kNCTau     = 4; //CC tau background in NC spectrum
static const int kCCTau     = 5; //CC tau background in CC spectrum
static const int kNCBeamNue = 6; //CC nue background in NC spectrum
static const int kCCBeamNue = 7; //CC nue background in CC spectrum
static const int kNCOscNue  = 8; //CC nue background in NC spectrum
static const int kCCOscNue  = 9; //CC nue background in CC spectrum

//......................................................................
NCExtrapolationFarNear::NCExtrapolationFarNear()
{
  fNCalls = 0;
  fDoExtrapolation = true;

  int numBins = kNumEnergyBinsFar;
  double bins[kNumEnergyBinsFar+1];

  //use NCType to fill array of bins as it has the right binning already
  for( int i = 0; i < numBins+1; ++i ){
    bins[i] = kEnergyBinsFar[i];
  }

  // Create a true energy binning scheme for the
  // F/N matrix

  fTrue_bin_width = bins[numBins]/double(kNumTrueEnergyBins);
  vector<double> true_bins(kNumTrueEnergyBins+1, 0);

  for( int i = 0; i < kNumTrueEnergyBins+1; ++i )
    true_bins[i] = i*fTrue_bin_width;


  //--------------------------------------------------
  // Far Detector Histograms

  fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco", "NC True vs. Reco (Nom)",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco", "CC True vs. Reco (Nom)",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_bkg", "NC True vs. Reco (Nom)bkg",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_bkg", "CC True vs. Reco (Nom)bkg",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_tau", "NC True vs. Reco (Nom)tau",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_tau", "CC True vs. Reco (Nom)tau",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_beamnue", "NC True vs. Reco (Nom)beamnue",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_beamnue", "CC True vs. Reco (Nom)beamnue",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_oscnue", "NC True vs. Reco (Nom)oscnue",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
  fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_oscnue", "CC True vs. Reco (Nom)oscnue",
                                         numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));

  fFN.push_back( new TH1D( "FNNC" , "FN NC", numBins, bins ) );
  fFN.push_back( new TH1D( "FNCC" , "FN CC", numBins, bins ) );

  fFD_fit.push_back( new TH1D( "FD_fit_nc" , "FD_fit NC", numBins, bins ) );
  fFD_fit.push_back( new TH1D( "FD_fit_cc" , "FD_fit CC", numBins, bins ) );

  //--------------------------------------------------
  // Near Detector Histograms

  fND_data.push_back( new TH1D( "ND_data_NC" , "ND_data_NC", numBins, bins ) );
  fND_data.push_back( new TH1D( "ND_data_CC" , "ND_data_CC", numBins, bins ) );

  //------------------------------
  // muon neutrinos

  fND_MC.push_back( new TH1D( "ND_MC_NC" , "ND_MC_NC", numBins, bins ) );
  fND_MC.push_back( new TH1D( "ND_MC_CC" , "ND_MC_CC", numBins, bins ) );
  fND_MC.push_back( new TH1D( "ND_MC_NC_ccbkg" , "ND_MC_NC-ccbkg", numBins, bins ) );
  fND_MC.push_back( new TH1D( "ND_MC_CC_ncbkg" , "ND_MC_CC-ncbkg", numBins, bins ) );

  fND_MC.push_back( new TH1D( "ND_MC_NC_tau" , "ND_MC_NC-tau", numBins, bins ) );
  fND_MC.push_back( new TH1D( "ND_MC_CC_tau" , "ND_MC_CC-tau", numBins, bins ) );

  fND_MC.push_back( new TH1D( "ND_MC_NC_beamnue" , "ND_MC_NC-beamnue", numBins, bins ) );
  fND_MC.push_back( new TH1D( "ND_MC_CC_beamnue" , "ND_MC_CC-beamnue", numBins, bins ) );
  fND_MC.push_back( new TH1D( "ND_MC_NC_oscnue" , "ND_MC_NC-oscnue", numBins, bins ) );
  fND_MC.push_back( new TH1D( "ND_MC_CC_oscnue" , "ND_MC_CC-oscnue", numBins, bins ) );

  fFD_data[kNCSig] = new TH1D("FD_data_NC", "FD data NC", kNumEnergyBinsFar, kEnergyBinsFar);
  fFD_data[kCCSig] = new TH1D("FD_data_CC", "FD data CC", kNumEnergyBinsFar, kEnergyBinsFar);

  fSmoothWidth = 0; // Don't do any smoothing
}
//----------------------------------------------------------------------



//......................................................................
NCExtrapolationFarNear::~NCExtrapolationFarNear(){}


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

  r.UnLockValues();

  r.Set("FarNearSmoothingWidth", 0.05);

  r.LockValues();
  return r;
}


//----------------------------------------------------------------------
//this method allows the user to pass in appropriate settings for the
//extrapolation such as fLocations of files, whether to generate the files
//again, etc.  the registry comes from the job module GetConfig method
void NCExtrapolationFarNear::Prepare(const Registry& r)
{
  MSG("NCExtrapolationFarNear", Msg::kInfo) <<  "PREPARING F/N EXTRAPOLATION" << endl;

  // the base class takes care of turning the systematic parameters on/off
  // and setting their adjusted values if any

  NCExtrapolation::Prepare(r);

  double      tmpd;

  //  if(r.Get("UE3SqrVal",             tmpd)) fUE3Sqr                     = tmpd;

  if(r.Get("FarNearSmoothingWidth", tmpd)) fSmoothWidth = tmpd;
}


//----------------------------------------------------------------------
template<class T> void NCExtrapolationFarNear::
HistSmoothHelper(T* h, bool is2D, double weight, double x, double y)
{
  if(fSmoothWidth == 0){
    // Shortcut all the below and just fill the histogram
    if(is2D)
      ((TH2D*)h)->Fill(x, y, weight);
    else
      h->Fill(x, weight);
    return;
  }

  // Paranoia over some scary text in the documentation for TAxis::GetBin
  h->SetBit(TH1::kCanRebin, false);

  TAxis* xaxis = h->GetXaxis();
  const int nbinsx = xaxis->GetNbins();

  // Warning - this number will only work in axis calls, not histogram calls
  const int xbin = xaxis->FindBin(x);
  const double lo = xaxis->GetBinLowEdge(xbin);
  const double hi = xaxis->GetBinUpEdge(xbin);

  const int ybin = is2D ? h->GetYaxis()->FindBin(y) : 0;

  const int histBin = xbin + (nbinsx+2)*ybin;

  const double dist_lo = x-lo;
  const double dist_hi = hi-x;

  // Don't smear events into the underflow bin
  const bool do_lo = xbin > 1 && dist_lo < fSmoothWidth*x;
  // Don't smear events into the overflow bin
  const bool do_hi = xbin < nbinsx && dist_hi < fSmoothWidth*x;

  double weight_here = 1;
  double weight_lo = 0;
  double weight_hi = 0;

  if(do_lo){
    const double transfer = .5*(1-dist_lo/(fSmoothWidth*x));
    weight_here -= transfer;
    weight_lo += transfer;
  }
  if(do_hi){
    const double transfer = .5*(1-dist_hi/(fSmoothWidth*x));
    weight_here -= transfer;
    weight_hi += transfer;
  }

  assert(weight_here >= .499);
  assert(weight_lo <= .501);
  assert(weight_hi <= .501);
  assert(weight_here <= 1);
  assert(weight_lo >= 0);
  assert(weight_hi >= 0);
  assert(fabs(weight_here+weight_lo+weight_hi-1) < .001);

  h->fArray[histBin] += weight*weight_here;

  if(do_lo){
    h->fArray[histBin-1] += weight*weight_lo;
  }
  if(do_hi){
    h->fArray[histBin+1] += weight*weight_hi;
  }
}

template void NCExtrapolationFarNear::
HistSmoothHelper<TH1D>(TH1D*, bool, double, double, double);
template void NCExtrapolationFarNear::
HistSmoothHelper<TH2D>(TH2D*, bool, double, double, double);

void NCExtrapolationFarNear::FillHistogramSmoothed(TH1D* h, double x, double weight)
{
  HistSmoothHelper(h, false, weight, x);
}


void NCExtrapolationFarNear::FillHistogramSmoothed2D(TH2D* h, double x, double y, double weight)
{
  HistSmoothHelper(h, true, weight, x, y);
}


//----------------------------------------------------------------------
void NCExtrapolationFarNear::FillDataMCHistogramsNear(NCBeam* beam,
                                                      NCType::EEventType nccc,
                                                      NC::SystPars systPars)
{
  using Detector::kNear;

  //loop over NCEnergy bin objects in the near detector
  //to get the energies from the data events in each bin
  //and fill the fND_data histogram


  // Alters the systematics 'weights' to reflect the current
  // value returned by the fitter.

  int sig         = kNCSig;
  int bkg         = kNCBkg;
  int tau         = kNCTau;
  int beamnue     = kNCBeamNue;
  int oscnue      = kNCOscNue;
  double bkgShift = systPars.CCBkgScale();
  double cleanShift = 1.;

  if(nccc == NCType::kCC){
    sig      = kCCSig;
    bkg      = kCCBkg;
    tau      = kCCTau;
    beamnue  = kCCBeamNue;
    oscnue   = kCCOscNue;
    bkgShift = systPars.NCBkgScale();
  }

  fND_data[sig]->Reset( "ICE" );
  fND_MC[sig]->Reset( "ICE" );
  fND_MC[bkg]->Reset( "ICE" );
  fND_MC[tau]->Reset( "ICE" );
  fND_MC[beamnue]->Reset( "ICE" );
  fND_MC[oscnue]->Reset( "ICE" );

  for(int b  = 0; b < beam->GetNumberEnergyBins(nccc); ++b){
    NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);

    const int dataSize       = bin->GetDataVectorSize();
    const int mcSize         = bin->GetMCSignalVectorSize();
    const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
    const int mcSize_tau     = bin->GetMCNuTauVectorSize();
    const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
    const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();

    // The variables that get filled with event information
    double trueEnergy    = 0;
    double recoShowerE   = 0;
    double recoMuonE     = 0;
    double energy        = 0;
    double weight        = 0;
    int    flavor        = 0;

    for(int e = 0; e < dataSize; ++e){
      bin->GetDataInformation(energy, weight, e);
      fND_data[sig]->Fill(energy, weight);
    }

    //------------------------------
    // numu
    for(int e = 0; e < mcSize; ++e){
      bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
      FillHistogramSmoothed(fND_MC[sig], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
    }

    for(int e = 0; e < mcSize_bkg; ++e){
      bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
      FillHistogramSmoothed(fND_MC[bkg], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*bkgShift*cleanShift);
    }

    //------------------------------
    // nutau
    for(int e = 0; e < mcSize_tau; ++e){
      bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
      FillHistogramSmoothed(fND_MC[tau], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
    }

    //------------------------------
    // beam nue
    for(int e = 0; e < mcSize_beamnue; ++e){
      bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
      FillHistogramSmoothed(fND_MC[beamnue], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
    }

    //------------------------------
    // oscillated nue
    for(int e = 0; e < mcSize_oscnue; ++e){
      bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
      FillHistogramSmoothed(fND_MC[oscnue], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
    }

  } // end for b
}

//----------------------------------------------------------------------
void NCExtrapolationFarNear::FillDataMCHistogramsFar(NCBeam* beam,
                                                     NCType::EEventType nccc,
                                                     NC::SystPars systPars)
{
  using Detector::kFar;

  //loop over NCEnergy bin objects in the near detector
  //to get the energies from the data events in each bin
  //and fill the fFD_data histogram
  double trueEnergy    = 0.;
  double recoShowerE   = 0.;
  double recoMuonE     = 0.;
  double weight        = 0.;
  int    flavor        = 0;

  int sig         = kNCSig;
  int bkg         = kNCBkg;
  int tau         = kNCTau;
  int beamnue     = kNCBeamNue;
  int oscnue      = kNCOscNue;
  double bkgShift = systPars.CCBkgScale();

  if(nccc == NCType::kCC){
    sig      = kCCSig;
    bkg      = kCCBkg;
    tau      = kCCTau;
    beamnue  = kCCBeamNue;
    oscnue   = kCCOscNue;
    bkgShift = systPars.NCBkgScale();
  }

  fNom_true_vs_reco[sig]->Reset( "ICE" );
  fNom_true_vs_reco[bkg]->Reset( "ICE" );
  fNom_true_vs_reco[tau]->Reset( "ICE" );
  fNom_true_vs_reco[beamnue]->Reset( "ICE" );
  fNom_true_vs_reco[oscnue]->Reset( "ICE" );


  const int B = beam->GetNumberEnergyBins(nccc);
  for(int b = 0; b < B; ++b){
    const NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);

    const int mcSize         = bin->GetMCSignalVectorSize();
    const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
    const int mcSize_tau     = bin->GetMCNuTauVectorSize();
    const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
    const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();

    //numu
    for(int e = 0; e < mcSize; ++e){
      bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed2D(fNom_true_vs_reco[sig], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
    }

    for(int e = 0; e < mcSize_bkg; ++e){
      bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed2D(fNom_true_vs_reco[bkg], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale()*bkgShift);
    }

    //nutau
    for(int e = 0; e < mcSize_tau; ++e){
      bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed2D(fNom_true_vs_reco[tau], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
    }

    //beamnue
    for(int e = 0; e < mcSize_beamnue; ++e){
      bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed2D(fNom_true_vs_reco[beamnue], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
    }

    //oscnue
    for(int e = 0; e < mcSize_oscnue; ++e){
      bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
      if(nccc == NCType::kNC) recoMuonE = 0;
      FillHistogramSmoothed2D(fNom_true_vs_reco[oscnue], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
    }

    const int N = bin->GetDataVectorSize();
    for(int n = 0; n < N; ++n){
      double energy, weight;
      bin->GetDataInformation(energy, weight, n);
      fFD_data[sig]->Fill(energy, weight);
    } // end for n
  } // end for b
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationFarNear::
ConstructFarSpectrum(const NC::OscProb::OscPars* pars)
{
  static vector<double*> survivalProbability;

  static bool first = true;
  if(first){
    first = false;

    fNominal.resize(fNom_true_vs_reco.size());
    fOsc.resize(fNom_true_vs_reco.size());

    for(unsigned int h = 0; h < fNom_true_vs_reco.size(); ++h){
      fNominal[h].resize(kNumEnergyBinsFar);
      fOsc[h].resize(kNumEnergyBinsFar);
    }

    survivalProbability.resize(fNom_true_vs_reco.size());
    for(unsigned int n = 0; n < fNom_true_vs_reco.size(); ++n)
      survivalProbability[n] = new double[kNumTrueEnergyBins];
  }

  static const NC::OscProb::OscPars* prevpars = 0;

  // Transition probability does not depend on reco energy - so we
  // can precalculate it out here.
  // Reuse the calculations from last time round if the oscillations
  // parameters haven't changed since then.
  if(!prevpars || !pars->IsEquiv(prevpars)){
    for(int h = 0; h < int(fNom_true_vs_reco.size()); ++h){
      NCType::EOscMode oscType;
      NCType::EEventType interaction;

      switch(h){
      case kNCSig:
      case kCCBkg:
        interaction = NCType::kNC;
        oscType     = NCType::kNuMuToNuMu;
        break;
      case kCCSig:
      case kNCBkg:
        interaction = NCType::kCC;
        oscType     = NCType::kNuMuToNuMu;
        break;
      case kNCTau:
      case kCCTau:
        interaction = NCType::kCC;
        oscType     = NCType::kNuMuToNuTau;
        break;
      case kNCBeamNue:
      case kCCBeamNue:
        interaction = NCType::kCC;
        oscType     = NCType::kNuEToNuE;
        break;
      case kNCOscNue:
      case kCCOscNue:
        interaction = NCType::kCC;
        oscType     = NCType::kNuMuToNuE;
        break;
      default:
        MSG("NCExtrapolationFarNear", Msg::kError) << "Interaction type or flavor is incorrect" << endl;
        return;
      }

      double* sph = survivalProbability[h];
      for(int j = 0; j < kNumTrueEnergyBins; ++j){
        sph[j] = pars->TransitionProbability(oscType,
                                             interaction,
                                             NCType::kBaseLineFar,
                                             (j+0.5)*fTrue_bin_width);
      } // end for j
    } // end for h
  } // end if


  delete prevpars;
  // Sadly OscPars::Clone isn't labelled const for technical reasons.
  // Need to cast away constness so that we can call it.
  prevpars = const_cast<NC::OscProb::OscPars*>(pars)->Clone();


  const int nx = kNumEnergyBinsFar+2;
  for(unsigned int h = 0; h < fNom_true_vs_reco.size(); ++h){
    const double* t2ra = fNom_true_vs_reco[h]->fArray;
    vector<double>& nh = fNominal[h];
    vector<double>& oh = fOsc[h];
    for(int i = 0; i < kNumEnergyBinsFar; ++i){
      double& nhi = nh[i];
      double& ohi = oh[i];
      nhi = 0;
      ohi = 0;
      const double* sph = survivalProbability[h];
      const int oset = i+1+nx;
      for(int j = 0; j < kNumTrueEnergyBins; ++j){
        //      const double tmpd = t2r->GetCellContent(i+1, j+1);
        const double tmpd = t2ra[nx*j+oset];
        nhi += tmpd;
        ohi += sph[j]*tmpd;
      }//end loop over true bins
    } // end for i
  }//end loop over histograms

  for(int i = 0; i < kNumEnergyBinsFar; ++i){
    const double totalMCNC_fd = fOsc[kNCSig][i]
      + fOsc[kNCBkg][i]
      + fOsc[kNCTau][i]
      + fOsc[kNCBeamNue][i]
      + fOsc[kNCOscNue][i];

    const double totalMCCC_fd = fOsc[kCCSig][i]
      + fOsc[kCCBkg][i]
      + fOsc[kCCTau][i]
      + fOsc[kCCBeamNue][i]
      + fOsc[kCCOscNue][i];

    // TODO - implement this bit in terms of GetNearMCSpectra. Needs some
    // work to keep the cacheing correct.
    const double totalMCNC_nd = fND_MC.at(kNCSig)->GetBinContent(i+1)
      + fND_MC.at(kNCBkg)->GetBinContent(i+1)      // numu
      + fND_MC.at(kNCTau)->GetBinContent(i+1)      // nutau
      + fND_MC.at(kNCBeamNue)->GetBinContent(i+1)  // beamnue
      + fND_MC.at(kNCOscNue)->GetBinContent(i+1);  // oscnue

    const double totalMCCC_nd = fND_MC.at(kCCSig)->GetBinContent(i+1)
      + fND_MC.at(kCCBkg)->GetBinContent(i+1)      // numu
      + fND_MC.at(kCCTau)->GetBinContent(i+1)      // nutau
      + fND_MC.at(kCCBeamNue)->GetBinContent(i+1)  // beamnue
      + fND_MC.at(kCCOscNue)->GetBinContent(i+1);  // oscnue

    if(totalMCNC_nd > 0){
      fFN[kNCSig]->SetBinContent(i+1, totalMCNC_fd/totalMCNC_nd);
      // Set this spectrum up without F/N correction. Will be overwritten if
      // fDoExtrapolation. If not then we need it done this way.
      fFD_fit[kNCSig]->SetBinContent(i+1, totalMCNC_fd);
    }
    else if(totalMCNC_nd == 0 && fUseNC){
      MSG("NCExtrapolationFarNear", Msg::kInfo) << "MC ND NC HISTO "
                                               << "HAS ZERO BIN:"
                                               << i+1 << ". PLEASE FIX IT"
                                               << endl;
    }

    if(totalMCCC_nd > 0){
      fFN[kCCSig]->SetBinContent(i+1, totalMCCC_fd/totalMCCC_nd);
      // Set this spectrum up without F/N correction. Will be overwritten if
      // fDoExtrapolation. If not then we need it done this way.
      fFD_fit[kCCSig]->SetBinContent(i+1, totalMCCC_fd);
    }
    else if(totalMCCC_nd == 0. && fUseCC){
      MSG("NCExtrapolationFarNear", Msg::kInfo) << "MC ND CC HISTO "
                                               << "HAS ZERO BIN:"
                                               << i+1 << ". PLEASE FIX IT"
                                               << endl;
    }
  }//end loop over reco energy bins

  if(fDoExtrapolation){
    fFD_fit[kNCSig]->Multiply(fND_data.at(kNCSig), fFN.at(kNCSig));
    fFD_fit[kCCSig]->Multiply(fND_data.at(kCCSig), fFN.at(kCCSig));
  }
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationFarNear::
FindSpectraForPars(const NC::OscProb::OscPars* oscPars,
                   const NC::SystPars& systPars_org,
                   vector<NCBeam::Info> beamsToUse,
                   vector<TH1*>& expvec,
                   vector<TH1*>& obsvec)
{
  for(int n = 0; n < NCType::kNumSystematicParameters; ++n){
    if(n == NCType::kNormalization) continue;
    if(n == NCType::kRelativeHadronicCalibration) continue;
    if(n == NCType::kAbsoluteHadronicCalibration) continue;
    if(n == NCType::kNCNearClean) continue;
    if(n == NCType::kCCBackground) continue;
    if(n == NCType::kTrackEnergy) continue;
    if(fCoordConv.IsFit(NCType::EFitParam(n))){
      MAXMSG("NCExtrapolationFarNear", Msg::kError, 10)
        << "Trying to fit systematic " << NCType::kParams[n].name
        << " (" << n << ") but FarNear doesn't know how."
        << " Systematic will essentially be ignored" << endl;
    }
  }

  if(fDoSystematicInterpolation) InitializeInterpolator(oscPars);

  /*
  // We definitely aren't smart enough to do more than one beam at a time
  assert(beamsToUse.size() == 1);
  const NCBeam::Info beamInfo = beamsToUse[0];
  */

  // This is a hack to make FarNear ~ work with interpolator. Since we ignore
  // beamsToUse and just do RunAll, then we should only send spectra back for
  // one run - arbitrarily RunI.
  if(fDoSystematicInterpolation &&
     beamsToUse[0].GetRunType() != NC::RunUtil::kRunI) return;

  NCBeam::Info beamInfo = NCBeam::Info(BeamType::kL010z185i,
                                       NC::RunUtil::kRunAll);
  vector<NCType::EFitParam> adj; vector<double> shift;
  beamsToUse[0].GetSystShifts(adj, shift);
  beamInfo.SetSystShifts(adj, shift);

  // Need a mutable copy
  NC::SystPars systPars = systPars_org;

  NCBeam* beam_fd = GetBeam(Detector::kFar, beamInfo);
  NCBeam* beam_nd = GetBeam(Detector::kNear, beamInfo);

  assert(beam_fd);
  assert(beam_nd);


  // These checks don't really belong here, but the fitter code is too stupid
  // to understand one-sided limits (and indeed limits that are different to
  // the one-sigma ranges) Until it's fixed just do it here.

  // TODO - these limits should go somewhere
  /*
  if(fCoordConv.IsFit(NCType::kNormalization)){
    double& norm = coords[fCoordConv.FitterIndex(NCType::kNormalization)];
    if(norm < -1){
      chisqr += 1e3*(norm+1)*(norm+1);
      norm = -1;
    }
  }
  if(fCoordConv.IsFit(NCType::kRelativeHadronicCalibration)){
    double& relhadr = coords[fCoordConv.FitterIndex(NCType::kRelativeHadronicCalibration)];
    if(relhadr < -1){
      chisqr += 1e3*(relhadr+1)*(relhadr+1);
      relhadr = -1;
    }
  }
  if(fCoordConv.IsFit(NCType::kAbsoluteHadronicCalibration)){
    double& abshadr = coords[fCoordConv.FitterIndex(NCType::kAbsoluteHadronicCalibration)];
    if(abshadr < -1){
      chisqr += 1e3*(abshadr+1)*(abshadr+1);
      abshadr = -1;
    }
  }
  */

  // TODO - THIS UNITARITY STUFF SERIOUSLY NEEDS TO BE DONE A LOT BETTER
  /*
  // 4 Flavour models unitarity constraints:
  //     |U_{mu1}|^{2} + |U_{mu2}|^{2} + |U_{mu3}|^{2}  + |U_{mu4}|^{2} > 1.0
  //  or |U_{e3}|^{2}  + |U_{mu3}|^{2} + |U_{tau3}|^{2} + |U_{s3}|^{2}  > 1.0

  // SterileFraction unitarity constraint:
  //     |U_mu3|^2 + f_s*(1-|U_mu3|^2) + |U_e3|^2 < 1

  if(fOscillationModel == NCType::kSterileFraction){
    const double umu3 = coords.at(fCoordConv.FitterIndex(NCType::kUMu3Sqr));
    const double ue3 = fUE3Sqr;
    const double fs = coords.at(fCoordConv.FitterIndex(NCType::kUS3Sqr));

    if(umu3+fs*(1-umu3)+ue3 > 1.001){
      if(!fSoftPenalty){
        MSG("NCExtrapolationFarNear", Msg::kVerbose) << "Unitarity violation - returning 1e6 fs=" << fs << " umu3=" << umu3 << " ue3=" << ue3 << endl;
        return 1.e6;
      }
      else{
        if(fs >= 1) return 1e6; // I don't think this can happen, but don't want to divide by zero if it does
        const double muproj = (1-ue3-fs)/(1-fs);
        const double dist = umu3-muproj;

        coords[NCType::kUMu3Sqr] = muproj;

        chisqr += dist*dist*1e3;

        MSG("NCExtrapolationFarNear", Msg::kVerbose) << "Unitarity violation - penalizing " << chisqr << " fs=" << fs << " umu3=" << umu3 << " ue3=" << ue3 << endl;
      }
    }
  }
  else if(fOscillationModel == NCType::k4Flavor ||
          fOscillationModel == NCType::k4FlavorDelta43Is0 ||
          fOscillationModel == NCType::k4FlavorDelta41Is0 ||
          fOscillationModel == NCType::kSterileFraction){
    double umu3 = coords.at(fCoordConv.FitterIndex(NCType::kUMu3Sqr));
//     if(fCoordConv.fLocUMu4Sqr > 0) umu3 += coords.at(fCoordConv.fLocUMu4Sqr);
    double us3 = coords.at(fCoordConv.FitterIndex(NCType::kUS3Sqr));
    double ue3 = fUE3Sqr;
    if(fOscillationModel == NCType::k4FlavorDelta41Is0
       && (us3 + umu3 + ue3) > 1.001){
      if(!fSoftPenalty){
        MSG("NCExtrapolationFarNear", Msg::kVerbose) << "Unitarity violation - returning 1e6 us3=" << us3 << " umu3=" << umu3 << " ue3=" << ue3 << endl;
        return 1.e6;
      }

      const double root2 = sqrt(2.0);
      const double dist = (us3+umu3+fUE3Sqr-1)/root2;
      const double sproj = us3-dist/root2;
      const double muproj = umu3-dist/root2;

      coords[NCType::kUMu3Sqr] = muproj;
      coords[NCType::kUS3Sqr] = sproj;

      chisqr = dist*dist*1e3;

      MSG("NCExtrapolationFarNear", Msg::kVerbose) << "Unitarity violation - penalizing " << chisqr << " us3=" << us3 << " umu3=" << umu3 << " ue3=" << ue3 << endl;
    }
  }//end checks of unitarity
  */

  // If none of the systematic shifts have changed since last time
  // we were called then there's no need to redo all the
  // FillDataMCHistograms() calls. The actual oscillation happens
  // in ConstructFarSpectrum()

  // TODO - there is no reason to refill the data histograms when
  // systematics change.

  vector<TH1*> nd_ratios, fd_ratios;

  if(fInterpolator){
    const vector<double> shift = fCoordConv.VectorFromSystPars(systPars);

    vector<TH1*> interp = fInterpolator->GetInterpolatedSpectra(shift);

    assert(interp.size()%2 == 0); // Every FD has a corresponding ND
    // The ND ratios are the odd elements of the vector. Copy them out.
    for(unsigned int n = 1; n < interp.size(); n += 2)
      nd_ratios.push_back(interp[n]);
    // The FD ratios are the even elements of the vector. Copy them out.
    for(unsigned int n = 0; n < interp.size(); n += 2)
      fd_ratios.push_back(interp[n]);

    systPars = NC::SystPars(); // Turn off all the systematic shifts
  }





  static NC::SystPars oldsyst;
  static NCBeam::Info oldinfo;
  static bool firstTime = true;

  if(systPars != oldsyst || !(beamsToUse[0] == oldinfo) ||  firstTime){
    oldsyst = systPars;
    oldinfo = beamsToUse[0];
    firstTime = false;

    fFD_fit.at(kNCSig)->Reset("ICE"); //integral contents error
    fFD_fit.at(kCCSig)->Reset("ICE");

    if(fUseNC) FillDataMCHistogramsNear(beam_nd, NCType::kNC, systPars);
    if(fUseCC) FillDataMCHistogramsNear(beam_nd, NCType::kCC, systPars);

    fFD_data[kNCSig]->Reset("ICE");
    fFD_data[kCCSig]->Reset("ICE");

    if(fUseNC) FillDataMCHistogramsFar(beam_fd, NCType::kNC, systPars);
    if(fUseCC) FillDataMCHistogramsFar(beam_fd, NCType::kCC, systPars);
  }

  ConstructFarSpectrum(oscPars);

  if(fUseNC){
    expvec.push_back(fFD_fit[kNCSig]);
    obsvec.push_back(fFD_data[kNCSig]);
  }

  if(fUseCC){
    expvec.push_back(fFD_fit[kCCSig]);
    obsvec.push_back(fFD_data[kCCSig]);
  }

  if(!fd_ratios.empty()){
    assert(nd_ratios.size() == expvec.size());
    assert(fd_ratios.size() == expvec.size());
    for(unsigned int n = 0; n < expvec.size(); ++n){
      // The ND systematic shifts end up on the bottom of the F/N ratio
      // and hence we can just divide by them at the FD.
      if(fDoExtrapolation) expvec[n]->Divide(nd_ratios[n]);
      expvec[n]->Multiply(fd_ratios[n]);
    }
  }

  ++fNCalls;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationFarNear::CleanupSpectra(vector<TH1*> /*exps*/,
                                            vector<TH1*> /*obss*/)
{
  // No need to clear anything. We don't allocate anything new each time
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationFarNear::
WriteResources(const NC::OscProb::OscPars* trueOscPars)
{
  MSG("NCExtrapolationFarNear", Msg::kInfo) << "nCalls: " << fNCalls << endl;

  NCBeam* beam = GetBeam(Detector::kFar,
                         NCBeam::Info(BeamType::kL010z185i,
                                      NC::RunUtil::kRunAll));
  NCBeam* beam_nd = GetBeam(Detector::kNear,
                            NCBeam::Info(BeamType::kL010z185i,
                                         NC::RunUtil::kRunAll));


  NC::SystPars systPars = GetBestFitSysts();
  const NC::OscProb::OscPars* oscPars = GetBestFitOscPars();

  if(fUseNC) FillDataMCHistogramsNear(beam_nd, NCType::kNC, systPars);
  if(fUseCC) FillDataMCHistogramsNear(beam_nd, NCType::kCC, systPars);

  fFD_data[kNCSig]->Reset("ICE");
  fFD_data[kCCSig]->Reset("ICE");

  if(fUseNC) FillDataMCHistogramsFar(beam, NCType::kNC, systPars);
  if(fUseCC) FillDataMCHistogramsFar(beam, NCType::kCC, systPars);

  //resets the histograms - integral contents error
  fFD_fit.at(kNCSig)->Reset("ICE");
  fFD_fit.at(kCCSig)->Reset("ICE");

  if(oscPars) ConstructFarSpectrum(oscPars);

  delete oscPars;

  // We need the near detector spectrum to be set up so we can extrapolate it
  NCExtrapolation::FillResultHistograms(beam_nd);

  // these calls fill the NCBeam histograms using our calculation, which
  // allows for Far/Near ratio to be used. NB that 'beam' is only kRunAll
  // which means that the RunI plots in NCBeam come out without extrapolation.
  if(fUseNC) FillResultHistograms(beam, NCType::kNC, fNominal, fOsc);
  if(fUseCC) FillResultHistograms(beam, NCType::kCC, fNominal, fOsc);

  // Now that the beams should be ready, call up to NCExtrapolation

  NCExtrapolation::WriteResources(trueOscPars);

  // get a pointer to the current directory
  // this is one of the output files

  TDirectory* saveDir = gDirectory;

  for(unsigned int i = 0; i < fFD_fit.size();  ++i)  fFD_fit[i]->Write();
  for(unsigned int i = 0; i < fFN.size();      ++i)      fFN[i]->Write();
  for(unsigned int i = 0; i < fND_data.size(); ++i) fND_data[i]->Write();
  for(unsigned int i = 0; i < fND_MC.size();   ++i)   fND_MC[i]->Write();
  for(unsigned int i = 0; i < fNom_true_vs_reco.size();++i)
    fNom_true_vs_reco[i]->Write();

  saveDir->cd();
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
void NCExtrapolationFarNear::
FillResultHistograms(NCBeam* beam,
                     NCType::EEventType nccc,
                     vector<vector<double> >& nominal,
                     vector<vector<double> >& osc)
{
  int sig     = kNCSig;
  int bkg     = kNCBkg;
  int tau     = kNCTau;
  int beamnue = kNCBeamNue;
  int oscnue  = kNCOscNue;

  if(nccc == NCType::kCC){
    sig     = kCCSig;
    bkg     = kCCBkg;
    tau     = kCCTau;
    beamnue = kCCBeamNue;
    oscnue  = kCCOscNue;
  }

  for(int i = 0; i < kNumEnergyBinsFar; ++i){

    double totalNDMC = fND_MC.at(sig)->GetBinContent(i+1)
      + fND_MC.at(bkg)->GetBinContent(i+1)      // numu
      + fND_MC.at(tau)->GetBinContent(i+1)      // nutau
      + fND_MC.at(beamnue)->GetBinContent(i+1); // nue

    double ndRatio = totalNDMC ? fND_data.at(sig)->GetBinContent(i+1)/totalNDMC : 0;
    double eBeamVal = nominal[beamnue][i];
    double nooscVal = nominal[sig][i] + nominal[bkg][i] + eBeamVal;
    double fitVal   = fFD_fit[sig]->GetBinContent(i+1);
    double bkgVal   = osc[bkg][i];
    double tauVal   = osc[tau][i];
    double eOscVal  = osc[oscnue][i];

    beam->GetMCHistogram(nccc)->SetBinContent(i+1, nooscVal*ndRatio);
    beam->GetMCFitHistogram(nccc)->SetBinContent(i+1, fitVal);//already had ndRatio
    beam->GetMCBackgroundHistogram(nccc)->SetBinContent(i+1, bkgVal*ndRatio);
    beam->GetMCFitNuMuToNuTauHistogram(nccc)->SetBinContent(i+1, tauVal*ndRatio);
    beam->GetMCFitNuEToNuEHistogram(nccc)->SetBinContent(i+1, eBeamVal*ndRatio);
    beam->GetMCFitNuMuToNuEHistogram(nccc)->SetBinContent(i+1, eOscVal*ndRatio);
  }

  // We filled the histograms manually, and don't need NCBeam's automatic stuff
  // to happen. So inform it of that fact.
  beam->MarkHistogramsFilled();
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
bool NCExtrapolationFarNear::EnableNearToFarExtrapolation(bool enable)
{
  const bool ret = fDoExtrapolation;
  fDoExtrapolation = enable;
  return ret;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
// This interface is most convenient for internal use
void NCExtrapolationFarNear::GetNearMCSpectra(NCBeam* beam, NC::SystPars systs,
                                              TH1** nc, TH1** cc)
{
  if(nc){
    FillDataMCHistogramsNear(beam, NCType::kNC, systs);

    *nc = new TH1D("", "", kNumEnergyBinsFar, kEnergyBinsFar);

    for(int i = 1; i <= kNumEnergyBinsFar; ++i){
      const double totalMCNC_nd = fND_MC.at(kNCSig)->GetBinContent(i)
        + fND_MC.at(kNCBkg)->GetBinContent(i)      // numu
        + fND_MC.at(kNCTau)->GetBinContent(i)      // nutau
        + fND_MC.at(kNCBeamNue)->GetBinContent(i)  // beamnue
        + fND_MC.at(kNCOscNue)->GetBinContent(i);  // oscnue

      (*nc)->SetBinContent(i, totalMCNC_nd);
    } // end for i
  }
  if(cc){
    FillDataMCHistogramsNear(beam, NCType::kCC, systs);

    *cc = new TH1D("", "", kNumEnergyBinsFar, kEnergyBinsFar);

    for(int i = 1; i <= kNumEnergyBinsFar; ++i){
      const double totalMCCC_nd = fND_MC.at(kCCSig)->GetBinContent(i)
        + fND_MC.at(kCCBkg)->GetBinContent(i)      // numu
        + fND_MC.at(kCCTau)->GetBinContent(i)      // nutau
        + fND_MC.at(kCCBeamNue)->GetBinContent(i)  // beamnue
        + fND_MC.at(kCCOscNue)->GetBinContent(i);  // oscnue

      (*cc)->SetBinContent(i, totalMCCC_nd);
    } // end for i
  }
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
// This is the interface we need to implement to keep the interpolator happy
vector<const TH1*> NCExtrapolationFarNear::
GetNearMCSpectra(vector<NCBeam::Info> beamsToUse)
{
  // We're too stupid to handle anything else
  //  assert(beamsToUse.size() == 1);

  vector<const TH1*> ret;

  // This is a hack to make FarNear ~ work with interpolator. Since we ignore
  // beamsToUse and just do RunAll, then we should only send spectra back for
  // one run - arbitrarily RunI.
  if(fDoSystematicInterpolation &&
     beamsToUse[0].GetRunType() != NC::RunUtil::kRunI) return ret;

  NCBeam* beam = GetBeam(Detector::kNear, beamsToUse[0]);
  assert(beam);

  NC::SystPars systs; // no systematic shifts to be applied

  TH1 *nc, *cc;
  GetNearMCSpectra(beam, systs, fUseNC ? &nc : 0, fUseCC ? &cc : 0);

  if(fUseNC) ret.push_back(nc);
  if(fUseCC) ret.push_back(cc);
  return ret;
}
//----------------------------------------------------------------------

---