NCExtrapolation.cxx

Go to the documentation of this file.

//$Id: NCExtrapolation.cxx,v 1.75 2007/10/26 20:21:43 brebel Exp $
//
//NCExtrapolation.cxx
//
//class to hold pdfs for doing nccc separation
//
//B. Rebel 2/2007

#include "NCUtils/Extrapolation/NCExtrapolation.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 "Conventions/DetectorType.h"

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

#include <cassert>

using namespace NCType;

CVSID("$Id: NCExtrapolation.cxx,v 1.75 2007/10/26 20:21:43 brebel Exp $");


// ---------------------------------------------------------------------
// Minimize func within the hypercube r1-r2. Returns the minimum value found.
// The position of the minimum is returned in minCoord, which has precision
// in each direction given by precision.
template<class Callable>
double NCExtrapolation::FindMinNDim(const Callable& func,
                                    const CoordNDim& r1,
                                    const CoordNDim& r2,
                                    CoordNDim& minCoord,
                                    const CoordNDim& precision)
{
  // How many points to put check in a box, along each axis.
  const int searchGridSize = 2;

  for(unsigned int n = 0; n < precision.size(); ++n) assert(precision[n] > 0);
  // Check that the arguments have consistent dimensionality.
  assert(r1.size() == r2.size() && r2.size() == precision.size());
  const int numDims = r1.size();
  minCoord.resize(numDims);

  // Algorithm assumes that all r1 values are less than the corresponding r2.
  for(int n = 0; n < numDims; ++n) assert(r1[n] <= r2[n]);

  // Check if we are precise enough in all directions.
  bool isPrecise = true;
  for(int n = 0; n < numDims; ++n)
    if(r2[n]-r1[n] > precision[n])
      isPrecise = false;
  // If so, average the corner coordinates and return.
  if(isPrecise){
    for(int n = 0; n < numDims; ++n) minCoord[n] = (r1[n]+r2[n])/2;
    return func(minCoord);
  }

  // How far to go in any direction to reach the next grid point.
  CoordNDim step(numDims);
  for(int n = 0; n < numDims; ++n) step[n] = (r2[n]-r1[n])/(searchGridSize+1);

  // This part is confusing: We need a grid of points inside the search
  // volume. toCheck will eventually hold all the points we want to look at.
  std::vector<CoordNDim> toCheck;
  // To start with just insert the grid point nearest to r1.
  CoordNDim firstPoint(numDims);
  for(int n = 0; n < numDims; ++n) firstPoint[n] = r1[n]+step[n];
  toCheck.push_back(firstPoint);

  // Now, for every dimension...
  for(int n = 0; n < numDims; ++n){
    // ...we need to loop over all the points in the list so far...
    const int iMax = toCheck.size();
    for(int i = 0; i < iMax; ++i){
      for(int d = 2; d <= searchGridSize; ++d){
        //...copy that point...
        CoordNDim toAdd = toCheck[i];
        // ...and modify the value of the nth coordinate to be each of the
        // searchGridSize points along this dimension...
        toAdd[n] = r1[n]+d*step[n];
        // ...before appending this new point into the list of places to check.
        toCheck.push_back(toAdd);
      }
    }
  }


  // The minimimum value of func() found so far. FLT_MAX = infinity.
  double funcMin = FLT_MAX;
  for(unsigned int n = 0; n < toCheck.size(); ++n){
    const double funcVal = func(toCheck[n]);
    if(funcVal < funcMin){
      funcMin = funcVal;
      minCoord = toCheck[n];
    }
  }

  CoordNDim r1New = minCoord;
  CoordNDim r2New = minCoord;

  // The new hypercube to search in is the lowest point we found above
  // plus or minus one grid size in every direction.
  for(int n = 0; n < numDims; ++n){
    r1New[n] -= step[n];
    r2New[n] += step[n];
  }

  return FindMinNDim(func, r1New, r2New, minCoord, precision);
}

// ---------------------------------------------------------------------
// Find contour of height 'height' in the function 'func'. Search will be
// to precision 'precision'. You must pass the previously located minimum
// of the function in 'minCoord'. The contour must be continuous, finite
// and closed.
template<class Callable>
std::vector<Coord> NCExtrapolation::FindContour(const Callable& func,
                                                const Coord& minCoord,
                                                const double& height,
                                                const Coord& precision)
{
  const double dx = precision.x;
  const double dy = precision.y;

  typedef pair<int, int> Box; // first == left, second == top

  std::map<Box, bool> results; // Is this box already in the result set?
  std::vector<Box> todo;

  // Top left corner of the box containing the best fit.
  const int bfxi = int(minCoord.x/dx);
  const int bfyi = int(minCoord.y/dy);

  // Walk to the left
  for(int x = bfxi; ; --x){
    // As soon as we exceed the contour
    if(func(Coord(x*dx, bfyi*dy)) > height){
      // make it the first element in the todo list
      todo.push_back(Box(x, bfyi));
      // and stop going left.
      break;
    }
  }

  std::vector<Coord> contour;

  while(!todo.empty()){
    Box now = todo.back();
    todo.pop_back();

    // We end up evaluating the function multiple times at each point, so are
    // relying heavily on the cacheing done in GetChiSqrFromDerived()

    // Evaluate the function at all the corners.
    const double topleft =     func(Coord( now.first   *dx,  now.second   *dy));
    const double topright =    func(Coord((now.first+1)*dx,  now.second   *dy));
    const double bottomleft =  func(Coord( now.first   *dx, (now.second+1)*dy));
    const double bottomright = func(Coord((now.first+1)*dx, (now.second+1)*dy));

    // For each side of the box: see if one corner is high and the other low.
    // If so, and if the box adjacent on this side isn't already part of the
    // results, then append the adjacent box to the todo list and linearly
    // interpolate to find where to put the contour point on that edge.
    if( (topleft-height)*(topright-height) <= 0 ){
      Box up(now.first, now.second-1);
      if(!results[up]){
        todo.push_back(up);
        results[up] = true;
        contour.push_back(Coord((now.first+(height-topleft)/(topright-topleft))*dx, now.second*dy));
      }
    }
    if( (bottomleft-height)*(bottomright-height) <= 0){
      Box down(now.first, now.second+1);
      if(!results[down]){
        todo.push_back(down);
        results[down] = true;
        contour.push_back(Coord((now.first+(height-bottomleft)/(bottomright-bottomleft))*dx, (now.second+1)*dy));
      }
    }
    if( (topleft-height)*(bottomleft-height) <= 0){
      Box left(now.first-1, now.second);
      if(!results[left]){
        todo.push_back(left);
        results[left] = true;
        contour.push_back(Coord(now.first*dx, (now.second+(height-topleft)/(bottomleft-topleft))*dy));
      }
    }
    if( (topright-height)*(bottomright-height) <=0 ){
      Box right(now.first+1, now.second);
      if(!results[right]){
        todo.push_back(right);
        results[right] = true;
        contour.push_back(Coord((now.first+1)*dx, (now.second+(height-topright)/(bottomright-topright))*dy));
      }
    }
  }

  if(contour.size()) contour.push_back(contour[0]); // Close the loop.
  return contour;
}

// ---------------------------------------------------------------------
// Evaluate callable fun for a single varying value and some fixed values.
// Initialise the fixed values by passing a CoordNDim fix with the
// fixed values required, and the indices of the parameter to be varied.
// Then, calling this object as a function will evaluate fun() with
// these fixed values, and the values passed substituted in.
template<class Callable>
NCExtrapolation::GetFuncWithConstVars<Callable>::
GetFuncWithConstVars(const Callable& fun,
                     const CoordNDim& fix,
                     const std::vector<int>& vary)
  :func(fun), fixed(fix), tovary(vary){}

template<class Callable>
double NCExtrapolation::GetFuncWithConstVars<Callable>::
operator()(const CoordNDim& r) const
{
  assert(tovary.size() == r.size());
  // For every parameter in the list of what we are to vary - substitute
  // in, in order, from the parameters we were passed.
  for(unsigned int n = 0; n < tovary.size(); ++n) fixed[tovary[n]] = r[n];
  return func(fixed);
}


// ---------------------------------------------------------------------
// Minimize callable f with respect to one of its arguments. toMin is the
// position of the argument that will be minimized min, max and prec give
// its range and precision. Call the object like a function with a Coord2D,
// the minimized parameter will be inserted in the correct place in the
// argument sequence and the minimum value returned.
template<class Callable>
NCExtrapolation::GetFuncMinWRTSomeVars<Callable>::
GetFuncMinWRTSomeVars(const Callable& f,
                      const std::vector<int>& toMin,
                      const CoordNDim& min,
                      const CoordNDim& max,
                      const CoordNDim& prec)
  :func(f), tomin(toMin), r1(min), r2(max), precision(prec)
{
  assert(r1.size() == r2.size() && r2.size() == precision.size());

  // Where we are going to substiture in r.x and r.y
  xpos = ypos = -1;
  // Keep increasing the index tried until ypos is set.
  for(int n = 0; ypos < 0 ;++n){
    // Is n in the list of variables to minimize?
    bool inmin = false;
    for(unsigned int m = 0; m < tomin.size(); ++m)
      if(tomin[m] == n) inmin = true;
    // If not, then fill in xpos, and then ypos next time we find one.
    if(!inmin){
      if(xpos < 0)
        xpos = n;
      else
        ypos = n;
    }
  }
}

template<class Callable>
double NCExtrapolation::GetFuncMinWRTSomeVars<Callable>::
operator()(const Coord& r) const
{
  CoordNDim fix(r1.size()+2);
  fix[xpos] = r.x;
  fix[ypos] = r.y;
  GetFuncWithConstVars<Callable> f(func, fix, tomin);
  CoordNDim minCoord;
  return NCExtrapolation::FindMinNDim(f, r1, r2, minCoord, precision);
}


// Override TMinuit's Eval function - so that this gets called as the
// minimization function. Forward to the functoid passed in the constructor.
template<class Callable>
Int_t NCExtrapolation::NCExtrapolationMinuit<Callable>::Eval(Int_t /*npar*/,
                                                             Double_t* /*grad*/,
                                                             Double_t& fval,
                                                             Double_t* par,
                                                             Int_t flag)
{
  // If we are being asked to do something that isn't just
  // calculate the function then we don't know how...
  assert(flag == 4);
  CoordNDim r(numpar);
  for(int n = 0; n < numpar; ++n) r[n] = par[n];
  fval = func(r);
  return 0;
}


// ---------------------------------------------------------------------
// Minimize callable f with respect to some of its arguments. xcoord and
// ycoord are positions of the arguments that shouldn't be minimized over.
// min and max give the range of all parameters (values for the parameters
//  not being minimized over are ignored).
template<class Callable>
NCExtrapolation::GetFuncMinWRTSomeVarsHybrid<Callable>::
GetFuncMinWRTSomeVarsHybrid(const Callable& f,
                            const int xcoord,
                            const int ycoord,
                            const CoordNDim& min,
                            const CoordNDim& max)
  :func(f), xpos(xcoord), ypos(ycoord), r1(min), r2(max)
{
  assert(r1.size() == r2.size());
  assert(xpos < int(r1.size()) && ypos < int(r1.size()));
  assert(ypos > xpos);

  pminu = new NCExtrapolationMinuit<Callable>(func, r1.size());

  pminu->SetPrintLevel(-1);
  //  minu.Command("set strategy 2"); // Be careful, instead of fast.

  for(int n = 0; n < int(r1.size()); ++n){
    const TString title = (n == xpos) ? "x" : (n == ypos) ? "y" : "";
    pminu->DefineParameter(n, title, (r1[n]+r2[n])/2, (r2[n]-r1[n])/4, 0, 0);
  }
}

template<class Callable>
double NCExtrapolation::GetFuncMinWRTSomeVarsHybrid<Callable>::
operator()(const Coord& r) const
{
  assert(pminu);

  // Minuit tries to return parameters we don't want to keep.
  double junkd; int junki;

  // Add one to the index because of ridiculous fortran numbering scheme.
  double args1[] = {xpos+1, r.x};
  pminu->mnexcm("set param", args1, 2, junki);
  double args2[] = {ypos+1, r.y};
  pminu->mnexcm("set param", args2, 2, junki);

  pminu->mnfixp(xpos, junki);
  pminu->mnfixp(ypos-1, junki); // The parameters all got renumbered when we fixed x.

  // Simplex seems to be much faster than Migrad.
  // I think Migrad can fail when put in a bad place, and spends
  // a long time doing it.
  // However - I have had bad results using Simplex, so stick with Migrad.
  //  minu.mnsimp();
  if(pminu->Migrad())
    ;//return FLT_MAX; // this messes up the contours...

  double minval;
  pminu->mnstat(minval, junkd, junkd, junki, junki, junki);
  pminu->mnfree(0); // Free all the parameters, we will refix them next time round.
  return minval;
}


// ---------------------------------------------------------------------
// Wrap NCExtrapolation::GetChiSqrForMap so that it can be called more simply,
// additionally keep a cache to speed up repeated queries.
// Evaluates points outside the physical region using a penalty term.

std::map<CoordNDim, double> NCExtrapolation::GetChiSqrFromDerived::cache;

double NCExtrapolation::GetChiSqrFromDerived::operator()(const CoordNDim& r) const
{
  const std::map<CoordNDim, double>::const_iterator it = cache.find(r);
  if(it != cache.end()) return it->second;
  // Maximum cache size is arbitrary, tune for performance.
  const unsigned int maxCacheSize = int(1e7); // on the order of 10meg
  if(cache.size() > maxCacheSize) cache.clear();

  CoordNDim evalAt = r;

  if(usePenalty){
    // If the point is outside the physical region - move it to the boundary
    for(unsigned int n = 0; n < start.size(); ++n){
      if(r[n] < start[n]) evalAt[n] = start[n];
      if(r[n] > end[n]) evalAt[n] = end[n];
    }
  }

  double ret = p->GetChiSqrForMap(evalAt, names);

  if(usePenalty){
    // If we had to move the point then apply a quadratic penalty term.
    // The factor of 1000 was arrived at through trial and error.
    for(unsigned int n = 0; n < start.size(); ++n){
      const double wsqr = (start[n]-end[n])*(start[n]-end[n]);
      if(r[n] < start[n]) ret+=1e3*(start[n]-r[n])*(start[n]-r[n])/wsqr;
      if(r[n] > end[n]) ret+=1e3*(r[n]-end[n])*(r[n]-end[n])/wsqr;
    }
  }

  cache[r] = ret;
  return ret;
}


//......................................................................
NCExtrapolation::NCExtrapolation() :
  fLocNormalization(0),
  fLocRelativeHadronicCalibration(0),
  fLocAbsoluteHadronicCalibration(0),
  fLocTrackEnergy(0),
  fLocNCBkg(0),
  fLocCCBkg(0),
  fLocPhi43(0),
  fLocPE(0),
  fLocPMu(0),
  fLocPS(0),
  fShowerScale(1.),
  fTrackScale(1.),
  fWeightNormalization(0.),
  fWeightRelativeHadronicCalibration(0.),
  fWeightAbsoluteHadronicCalibration(0.),
  fWeightTrackEnergy(0.),
  fWeightNCBkg(0.),
  fWeightCCBkg(0.),
  fUseWSDeltaMContour(false),
  fUseWSU3Contour(false),
  fUseWSWMuContour(false),
  fUseWSWEContour(false),
  fUseU3DeltaMContour(false),
  fUseU3WMuContour(false),
  fUseU3WEContour(false),
  fUseDeltaMWMuContour(false),
  fUseDeltaWEContour(false),
  fUseWEWMuContour(false)
{
  fUtils = new NCAnalysisUtils();
  fUseCC = false;
  fUseNC = false;

  for(int i = 0; i < NCType::kNumParameters; ++i)
    fSystematicsToUse.push_back(false);


  fDeltaChiSqrDeltaMVsFSterile = new TH2F("DeltaChiSqrDeltaMVsFSterile",
                                          "",
                                          NCType::kNumFSterileBins,
                                          NCType::kFSterileStart,
                                          NCType::kFSterileEnd,
                                          NCType::kNumDeltaMSqrBins,
                                          NCType::kDeltaMSqrStart,
                                          NCType::kDeltaMSqrEnd);
  fDeltaChiSqrDeltaMVsFSterile->SetXTitle("f_{s}");
  fDeltaChiSqrDeltaMVsFSterile->SetYTitle("#Delta m^{2} (eV^{2})");

  fDeltaChiSqrSinSqrVsFSterile = new TH2F("DeltaChiSqrSinSqrVsFSterile",
                                          "",
                                          NCType::kNumFSterileBins,
                                          NCType::kFSterileStart,
                                          NCType::kFSterileEnd,
                                          NCType::kNumSinSqrBins,
                                          NCType::kSinSqrStart,
                                          NCType::kSinSqrEnd);
  fDeltaChiSqrSinSqrVsFSterile->SetXTitle("f_{s}");
  fDeltaChiSqrSinSqrVsFSterile->SetYTitle("sin^{2}(2#theta)");

  fDeltaChiSqrDeltaMVsSinSqr = new TH2F("DeltaChiSqrDeltaMVsSinSqr",
                                        "",
                                        NCType::kNumSinSqrBins,
                                        NCType::kSinSqrStart,
                                        NCType::kSinSqrEnd,
                                        NCType::kNumDeltaMSqrBins,
                                        NCType::kDeltaMSqrStart,
                                        NCType::kDeltaMSqrEnd);
  fDeltaChiSqrDeltaMVsSinSqr->SetXTitle("sin^{2}(2#theta)");
  fDeltaChiSqrDeltaMVsSinSqr->SetYTitle("#Delta m^{2} (eV^{2})");

}

//.......................................................................
NCExtrapolation::NCExtrapolation(NCAnalysisUtils *utils,
                                 std::map<Int_t,double>& dataNear,
                                 std::map<Int_t,double>& mcNear,
                                 std::map<Int_t,double>& dataFar,
                                 std::map<Int_t,double>& mcFar,
                                 bool useCC, bool useNC) :
  fDataPOTNear(dataNear),
  fMCPOTNear(mcNear),
  fDataPOTFar(dataFar),
  fMCPOTFar(mcFar),
  fUseCC(useCC),
  fUseNC(useNC),
  fUtils(utils),
  fLocNormalization(0),
  fLocRelativeHadronicCalibration(0),
  fLocAbsoluteHadronicCalibration(0),
  fLocTrackEnergy(0),
  fLocNCBkg(0),
  fLocCCBkg(0),
  fLocPhi43(0),
  fLocPE(0),
  fLocPMu(0),
  fLocPS(0),
  fShowerScale(1.),
  fTrackScale(1.),
  fWeightNormalization(0.),
  fWeightRelativeHadronicCalibration(0.),
  fWeightAbsoluteHadronicCalibration(0.),
  fWeightTrackEnergy(0.),
  fWeightNCBkg(0.),
  fWeightCCBkg(0.)
{
  BeamType::BeamType_t beam = BeamType::kL010z185i;
  int run = NCType::kRunI;

  for(std::map<Int_t,double>::const_iterator i = fDataPOTNear.begin();
      i != fDataPOTNear.end(); ++i){

    NCType::ConvertIndexToBeamRun(i->first, beam, run);

    MSG("NCExtrapolation", Msg::kInfo)
      << "fDataPOTNear["
      << BeamType::AsString(beam) << ", " << run << "] = "
      << i->second << endl;
  }
  for(std::map<Int_t,double>::const_iterator i = fMCPOTNear.begin();
      i != fMCPOTNear.end(); ++i){

    NCType::ConvertIndexToBeamRun(i->first, beam, run);

    MSG("NCExtrapolation", Msg::kInfo)
      << "fMCPOTNear["
      << BeamType::AsString(beam) << ", " << run << "] = "
      << i->second << endl;
  }

  for(std::map<Int_t,double>::const_iterator i = fDataPOTFar.begin();
      i != fDataPOTFar.end(); ++i){

    NCType::ConvertIndexToBeamRun(i->first, beam, run);

    MSG("NCExtrapolation", Msg::kInfo)
      << "fDataPOTFar["
      << BeamType::AsString(beam) << ", " << run << "] = "
      << i->second << endl;
  }

  for(std::map<Int_t,double>::const_iterator i = fMCPOTFar.begin();
      i != fMCPOTFar.end(); ++i){

    NCType::ConvertIndexToBeamRun(i->first, beam, run);

    MSG("NCExtrapolation", Msg::kInfo)
      << "fMCPOTFar["
      << BeamType::AsString(beam) << ", " << run << "] = "
      << i->second << endl;
  }

  for(int i = 0; i < NCType::kNumParameters; ++i){
    fSystematicsToUse.push_back(false);
    fSystematicsAdjust.push_back(NCType::kParameterDefaults[i]);
  }


  fDeltaChiSqrDeltaMVsFSterile = new TH2F("DeltaChiSqrDeltaMVsFSterile",
                                          "DeltaChiSqrDeltaMVsFSterile",
                                          NCType::kNumFSterileBins,
                                          NCType::kFSterileStart,
                                          NCType::kFSterileEnd,
                                          NCType::kNumDeltaMSqrBins,
                                          NCType::kDeltaMSqrStart,
                                          NCType::kDeltaMSqrEnd);

  fDeltaChiSqrSinSqrVsFSterile = new TH2F("DeltaChiSqrSinSqrVsFSterile",
                                          "DeltaChiSqrSinSqrVsFSterile",
                                          NCType::kNumFSterileBins,
                                          NCType::kFSterileStart,
                                          NCType::kFSterileEnd,
                                          NCType::kNumSinSqrBins,
                                          NCType::kSinSqrStart,
                                          NCType::kSinSqrEnd);

  fDeltaChiSqrDeltaMVsSinSqr = new TH2F("DeltaChiSqrDeltaMVsSinSqr",
                                        "DeltaChiSqrDeltaMVsSinSqr",
                                        NCType::kNumSinSqrBins,
                                        NCType::kSinSqrStart,
                                        NCType::kSinSqrEnd,
                                        NCType::kNumDeltaMSqrBins,
                                        NCType::kDeltaMSqrStart,
                                        NCType::kDeltaMSqrEnd);
}

//----------------------------------------------------------------------
NCExtrapolation::~NCExtrapolation()
{
}

//----------------------------------------------------------------------
//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. See comment under DefaultConfig...
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string

  //get the value of theta_13
  if( r.Get("FixedTheta13", tmpd) ) fTheta13 = tmpd;

  //set the defaults for the systematic parameters
  //set them all off by default
  TString adjust = "Adjust";
  for(int i = 0; i < NCType::kNumParameters; ++i){
    if( r.Get(NCType::kParameterNames[i], tmpb)) fSystematicsToUse[i] = tmpb;
    if( r.Get(adjust + NCType::kParameterNames[i], tmpd)) fSystematicsAdjust[i] = tmpd;
  }

  MSG("NCExtrapolation", Msg::kInfo) << "PARAMETERS WE ARE USING:" << endl;
  for(int i = 0; i < NCType::kNumParameters; ++i){
    MSG("NCExtrapolation", Msg::kInfo) << i << ") "
                                       << NCType::kParameterNames[i]
                                       << " = " << fSystematicsToUse[i]
                                       << " with default = "
                                       << NCType::kParameterDefaults[i]
                                       << " adjust = "
                                       << fSystematicsAdjust[i] << endl;
  }

  tmpi = 0;
  tmps = "s";

  r.Get("LazyContoursUseMinuit", fLazyContoursUseMinuit);
  r.Get("LazyContoursWant68", fWant68);
  r.Get("LazyContoursWant90", fWant90);
  r.Get("LazyContoursWant99", fWant99);
  r.Get("LazyContoursWantFsterSinsq", fWantFsterSinsq);
  r.Get("LazyContoursWantFsterDmsq", fWantFsterDmsq);
  r.Get("LazyContoursWantSinsqDmsq", fWantSinsqDmsq);
  r.Get("AllowBestFitOutsideRange", fAllowBestFitOutsideRange);
  r.Get("AllowContoursOutsideRange", fAllowContoursOutsideRange);

  if(r.Get("UseWSDeltaMContour",  tmpb)) fUseWSDeltaMContour  = tmpb;
  if(r.Get("UseWSU3Contour",      tmpb)) fUseWSU3Contour      = tmpb;
  if(r.Get("UseWSWMuContour",     tmpb)) fUseWSWMuContour     = tmpb;
  if(r.Get("UseWSWEContour",      tmpb)) fUseWSWEContour      = tmpb;
  if(r.Get("UseU3DeltaMContour",  tmpb)) fUseU3DeltaMContour  = tmpb;
  if(r.Get("UseU3WMuContour",     tmpb)) fUseU3WMuContour     = tmpb;
  if(r.Get("UseU3WEContour",      tmpb)) fUseU3WEContour      = tmpb;
  if(r.Get("UseDeltaMWMuContour", tmpb)) fUseDeltaMWMuContour = tmpb;
  if(r.Get("UseDeltaWEContour",   tmpb)) fUseDeltaWEContour   = tmpb;
  if(r.Get("UseWEWMuContour",     tmpb)) fUseWEWMuContour     = tmpb;
  
  fXAxes.clear();
  fYAxes.clear();
    
  //set which axes to use in the contours
  if(fUseWSDeltaMContour)
    fXAxes.push_back(NCType::kFSterile);  fYAxes.push_back(NCType::kDeltaMSqr);
  if(fUseWSU3Contour)
    fXAxes.push_back(NCType::kFSterile);  fYAxes.push_back(NCType::kSinSqr);
  if(fUseWSWMuContour)
    fXAxes.push_back(NCType::kFSterile);  fYAxes.push_back(fLocPMu);
  if(fUseWSWEContour)
    fXAxes.push_back(NCType::kFSterile);  fYAxes.push_back(fLocPE);
  if(fUseU3DeltaMContour)
    fXAxes.push_back(NCType::kSinSqr);    fYAxes.push_back(NCType::kDeltaMSqr);
  if(fUseU3WMuContour)
    fXAxes.push_back(NCType::kSinSqr);    fYAxes.push_back(fLocPMu);
  if(fUseU3WEContour)
    fXAxes.push_back(NCType::kSinSqr);    fYAxes.push_back(fLocPE);
  if(fUseDeltaMWMuContour)
    fXAxes.push_back(NCType::kDeltaMSqr); fYAxes.push_back(fLocPMu);
  if(fUseDeltaWEContour)
    fXAxes.push_back(NCType::kDeltaMSqr); fYAxes.push_back(fLocPE);
  if(fUseWEWMuContour)
    fXAxes.push_back(fLocPE);             fYAxes.push_back(fLocPMu);

  int locCounter   = 2; // first three elements are osc parameters
  double max = 0.;
  double min = 0.;

  if (r.Get("Normalization", tmpb)){
    if ( tmpb == true ){

      min = -NCType::kParameterSigmas[NCType::kNormalization];
      max =  NCType::kParameterSigmas[NCType::kNormalization];

      fExtraVarsNames.push_back("Normalization");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocNormalization = ++locCounter;
    }
  }
  if (r.Get("RelativeHadronicCalibration", tmpb)){
    if ( tmpb == true ){

      min = -NCType::kParameterSigmas[NCType::kRelativeHadronicCalibration];
      max =  NCType::kParameterSigmas[NCType::kRelativeHadronicCalibration];

      fExtraVarsNames.push_back("Relative Hadronic Calibration");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocRelativeHadronicCalibration = ++locCounter;
    }
  }
  if (r.Get("AbsoluteHadronicCalibration", tmpb)){
    if ( tmpb == true ){

      min = -NCType::kParameterSigmas[NCType::kAbsoluteHadronicCalibration];
      max =  NCType::kParameterSigmas[NCType::kAbsoluteHadronicCalibration];

      fExtraVarsNames.push_back("Absolute Hadronic Calibration");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocAbsoluteHadronicCalibration = ++locCounter;
    }
  }
  if (r.Get("TrackEnergy", tmpb)){
    if ( tmpb == true ){

      min = -NCType::kParameterSigmas[NCType::kTrackEnergy];
      max =  NCType::kParameterSigmas[NCType::kTrackEnergy];

      fExtraVarsNames.push_back("Track Energy");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocTrackEnergy = ++locCounter;
    }
  }
  if (r.Get("NCBackground", tmpb)){
    if ( tmpb == true ){

      min = -NCType::kParameterSigmas[NCType::kNCBackground];
      max =  NCType::kParameterSigmas[NCType::kNCBackground];

      fExtraVarsNames.push_back("NC Background");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocNCBkg = ++locCounter;
    }
  }
  if (r.Get("CCBackground", tmpb)){
    if ( tmpb == true ){

      min = -NCType::kParameterSigmas[NCType::kCCBackground];
      max =  NCType::kParameterSigmas[NCType::kCCBackground];

      fExtraVarsNames.push_back("CC Background");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocCCBkg = ++locCounter;
    }
  }
  if (r.Get("WE", tmpb)){
    if ( tmpb == true ){

      min = 0.;
      max = 0.2;

      fExtraVarsNames.push_back("w_{e}");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocPE = ++locCounter;

      cout << "w_e location = " << fLocPE << endl;
    }
  }
  if (r.Get("WMu", tmpb)){
    if ( tmpb == true ){

      min = 0.;
      max = 0.5;

      fExtraVarsNames.push_back("w_{#mub}");
      fExtraVarsMin.push_back( min );
      fExtraVarsMax.push_back( max );
      fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
      fLocPMu = ++locCounter;

      cout << "w_mub location = " << fLocPMu << endl;
    }
  }
//   if (r.Get("Phi43", tmpb)){
//     if ( tmpb == true ){

//       min = 0.;
//       max = 2.*TMath::Pi();

//       fExtraVarsNames.push_back("phi_43");
//       fExtraVarsMin.push_back( min );
//       fExtraVarsMax.push_back( max );
//       fExtraVarsPrec.push_back( TMath::Abs((max - min)*.01) );
//       fLocPhi43 = ++locCounter;

//       cout << "phi_43 location = " << fLocPhi43 << endl;
//     }
//   }

  return;
}

//----------------------------------------------------------------------
void NCExtrapolation::AddEvent(ANtpHeaderInfo *headerInfo,
                               ANtpBeamInfo *beamInfo,
                               ANtpRecoInfo *recoInfo,
                               ANtpAnalysisInfo *analysisInfo,
                               ANtpTruthInfoBeam *truthInfo,
                               int runType,
                               bool useMCAsData,
                               int fileType)
{

  //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(recoInfo->inFiducialVolume < 1 ||
     (headerInfo->detector == (int)Detector::kNear
      && recoInfo->isCleanHighMultSnarl < 1)
     )
     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::kNear;
  BeamType::BeamType_t beamType = BeamType::FromZarko(beamInfo->beamType);
  int beamRunIndex = NCType::ConvertBeamRunToIndex(beamType, runType);
  //MC events only have run type = NCType::kRunI
  int beamRunIndexMC = NCType::ConvertBeamRunToIndex(beamType, NCType::kRunI);
  double dataPOT = fDataPOTNear[beamRunIndex];
  double mcPOT = fMCPOTNear[beamRunIndexMC];
  if(headerInfo->detector == (int)Detector::kFar){
    detector = Detector::kFar;
    dataPOT = fDataPOTFar[beamRunIndex];
    mcPOT = fMCPOTFar[beamRunIndexMC];
  }
  int beam = FindBeamIndex(detector, beamType, runType, true, dataPOT, mcPOT);

  fBeams[beam].AddEvent(headerInfo, beamInfo, recoInfo,
                        analysisInfo, truthInfo, useMCAsData, fileType);

  //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(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){
    for(int i = NCType::kRunII; i < NCType::kMaxRun+1; ++i){
      beam = FindBeamIndex(detector, beamType, i);
      if(beam < 100)
        fBeams[beam].AddEvent(headerInfo, beamInfo, recoInfo,
                              analysisInfo, truthInfo, useMCAsData, fileType);
    }//end loop over possible runs
  }//end if mc and needs to be added to other runs for this beam

  return;
}

//----------------------------------------------------------------------
void NCExtrapolation::PrepareNearDetector()
{
  return;
}

//----------------------------------------------------------------------
int NCExtrapolation::GetExtrapolationType() const
{
  return fExtrapolationType;
}

//----------------------------------------------------------------------
TH1F* NCExtrapolation::GetPredictedSpectrum(Detector::Detector_t detector,
                                            BeamType::BeamType_t beamType,
                                            int nccc)
{
  TH1F *spectrum = 0;

  for(uint i = 0; i < fBeams.size(); ++i){
    if(fBeams[i].GetBeamType() == beamType
       && fBeams[i].GetDetector() == detector)
      spectrum = fBeams[i].GetMCHistogram(nccc);
  }

  return spectrum;
}

//----------------------------------------------------------------------
TH1F* NCExtrapolation::GetFitSpectrum(Detector::Detector_t detector,
                                      BeamType::BeamType_t beamType,
                                      int nccc)
{
  TH1F *spectrum = 0;

  for(uint i = 0; i < fBeams.size(); ++i){
    if(fBeams[i].GetBeamType() == beamType
       && fBeams[i].GetDetector() == detector)
      spectrum = fBeams[i].GetMCFitHistogram(nccc);
  }

  return spectrum;
}

//----------------------------------------------------------------------
TH1F* NCExtrapolation::GetBackgroundSpectrum(Detector::Detector_t detector,
                                             BeamType::BeamType_t beamType,
                                             int nccc)
{
  TH1F *spectrum = 0;

  for(uint i = 0; i < fBeams.size(); ++i){
    if(fBeams[i].GetBeamType() == beamType
       && fBeams[i].GetDetector() == detector)
      spectrum = fBeams[i].GetMCBackgroundHistogram(nccc);
  }

  return spectrum;
}

//----------------------------------------------------------------------
TGraphAsymmErrors* NCExtrapolation::GetDataSpectrum(Detector::Detector_t detector,
                                                    BeamType::BeamType_t beamType,
                                                    int nccc)
{
  TGraphAsymmErrors *spectrum = 0;

  for(uint i = 0; i < fBeams.size(); ++i){
    if(fBeams[i].GetBeamType() == beamType
       && fBeams[i].GetDetector() == detector)
      spectrum = fBeams[i].GetDataGraph(nccc);
  }

  return spectrum;
}

//----------------------------------------------------------------------
void NCExtrapolation::ResetBeams()
{

  for(uint i = 0; i < fBeams.size(); ++i) fBeams[i].Reset();

  return;
}

//----------------------------------------------------------------------
NCBeam NCExtrapolation::GetNCBeam(Detector::Detector_t det,
                                  BeamType::BeamType_t beam,
                                  int runType)
{
  int beamPos = FindBeamIndex(det, beam, runType);

  if(beamPos > 100){

    MSG("NCExtrapolation", Msg::kWarning) << "could not find chosen beam "
                                          << BeamType::AsString(beam)
                                          << " in " << Detector::AsString(det)
                                          << " detector, returning empty beam "
                                          << beamPos
                                          << endl;
    NCBeam empty;
    return empty;
  }

  return fBeams[beamPos];

}

//----------------------------------------------------------------------
void NCExtrapolation::DrawContours(TH2F *deltaChiSqr,
                                   double bestX,
                                   double bestY)
{

  TString name = deltaChiSqr->GetName();
  TString clName[] = {name+"68CL", name+"90CL", name+"99CL"};
  TString drawOpt[] = {"cont1", "cont1same", "cont1same"};
  TCanvas *canv = new TCanvas("contoursCanv"+name,
                              "contours", 150, 10, 900, 600);

  //make histograms with the same dimensions as the one passed
  //into the method
  vector<TH2F *> contours;
  double contVals[] = {2.3, 4.6, 9.2};
  double contVals1[] = {0., 0., 1.e6};
  int lineStyles[] = {3, 2, 1};
  int colors[] = {2, 4, 1};
  int colors1[] = {10, 4, 10};
  for(int i = 0; i < 3; ++i){
    contours.push_back(new TH2F(clName[i], "",
                                deltaChiSqr->GetXaxis()->GetNbins(),
                                deltaChiSqr->GetXaxis()->GetXmin(),
                                deltaChiSqr->GetXaxis()->GetXmax(),
                                deltaChiSqr->GetYaxis()->GetNbins(),
                                deltaChiSqr->GetYaxis()->GetXmin(),
                                deltaChiSqr->GetYaxis()->GetXmax()));
    contVals1[2] = contVals[i];
    colors1[2] = colors[i];
    contours[i]->Add(deltaChiSqr);
    contours[i]->SetContour(3, contVals1);
    contours[i]->SetLineStyle(lineStyles[i]);
    contours[i]->SetLineColor(colors[i]);
    contours[i]->Draw(drawOpt[i]);
  }

  TMarker *bestFit = new TMarker(bestX, bestY, 29);
  bestFit->SetMarkerSize(3.5);
  bestFit->SetMarkerColor(1);
  bestFit->Draw("p");
  TString bestFitLabel = "Best Fit ";
  bestFitLabel += deltaChiSqr->GetXaxis()->GetTitle() + TString::Format(" = %3.2e,", bestX);
  bestFitLabel += deltaChiSqr->GetYaxis()->GetTitle() + TString::Format(" = %3.2e", bestY);

  TLegend *leg = new TLegend(0.75, 0.75, 0.95, 0.95);
  leg->SetBorderSize(0);
  leg->AddEntry(contours[0], "68% Confidence Level", "l");
  leg->AddEntry(contours[1], "90% Confidence Level", "l");
  leg->AddEntry(contours[2], "99% Confidence Level", "l");
  leg->AddEntry(bestFit, bestFitLabel, "p");
  leg->Draw();

  canv->Update();

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

  return;

}

//----------------------------------------------------------------------
//loop over all the beams.  draw pairs of plots for each beam in the
//chosen detector - NC and CC
void NCExtrapolation::DrawBestFitSpectra(Detector::Detector_t det)
{

  for(uint b = 0; b < fBeams.size(); ++b){
    if(fBeams[b].GetDetector() != det) continue;

    TString name = "bestFitSpectraCanv"+NCType::kExtrapolationNames[fExtrapolationType];
    name += DetectorType::AsString(det);
    name += BeamType::AsString(fBeams[b].GetBeamType());
    name += NCType::kRunNames[fBeams[b].GetRunType()];
    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", fBestDeltaMSqr);
    TString bestSinSqr;
    bestSinSqr.Format("sin^{2}(2#theta) = %3.2f", fBestSinSqr);
    TString bestFSterile;
    bestFSterile.Format("f_{s} = %3.2f", fBestFSterile);
    TString minChiSqr;
    minChiSqr.Format("#chi^{2} = %3.2f", fMinChiSqr);

    double maxY = 0.;

    for(int i = 0; i < 2; ++i){
      if(i > 0){
        pad2->cd();
      }

      maxY = 1.3*TMath::Max(fBeams[b].GetMCHistogram(i)->GetMaximum(),
                            fBeams[b].GetDataGraph(i)->GetHistogram()->GetMaximum());

      fBeams[b].GetMCHistogram(i)->SetMaximum(maxY);
      fBeams[b].GetMCHistogram(i)->Draw();
      //if the fit histogram was filled, draw the fits
      if(fBeams[b].GetMCFitHistogram(i)->Integral() > 0.){
        fBeams[b].GetMCFitHistogram(i)->Draw("same");
        fBeams[b].GetMCBackgroundHistogram(i)->Draw("same");
        fBeams[b].GetMCFitNuMuToNuTauHistogram(i)->Draw("same");
        fBeams[b].GetMCFitNuEToNuEHistogram(i)->Draw("same");
        fBeams[b].GetMCFitNuMuToNuEHistogram(i)->Draw("same");
      }
      else{
        //no fit, no point drawing NuMuToNuTau or NuMuToNuE
        fBeams[b].GetMCNoFitBackgroundHistogram(i)->Draw("same");
        fBeams[b].GetMCNoFitNuEToNuEHistogram(i)->Draw("same");
      }

      fBeams[b].GetDataGraph(i)->Draw("pe1same");

      TLegend *leg = new TLegend(0.75, 0.75, 0.95, 0.95);
      leg->SetBorderSize(0);
      leg->AddEntry(fBeams[b].GetDataGraph(i), "Data", "lp");
      leg->AddEntry(fBeams[b].GetMCHistogram(i), "Monte Carlo (Nominal)", "f");
      if(fBeams[b].GetMCFitHistogram(i)->Integral() > 0.){
        leg->AddEntry(fBeams[b].GetMCFitHistogram(i), "Monte Carlo (Best Fit)", "l");
        leg->AddEntry(fBeams[b].GetMCBackgroundHistogram(i), "Background", "lf");
        leg->AddEntry(fBeams[b].GetMCFitNuMuToNuTauHistogram(i), "Background - #nu_{#tau} CC", "lf");
        leg->AddEntry(fBeams[b].GetMCFitNuEToNuEHistogram(i), "Background - Beam #nu_{e} CC", "lf");
        leg->AddEntry(fBeams[b].GetMCFitNuMuToNuEHistogram(i), "Background - #nu_{e} CC", "lf");
      }
      else{
        leg->AddEntry(fBeams[b].GetMCNoFitBackgroundHistogram(i), "Background", "lf");
        leg->AddEntry(fBeams[b].GetMCNoFitNuEToNuEHistogram(i), "Background - Beam #nu_{e} CC", "lf");
      }
//       leg->AddEntry(fBeams[b].GetDataGraph(i), bestDeltaM, "");
//       leg->AddEntry(fBeams[b].GetDataGraph(i), bestFSterile, "");
//       leg->AddEntry(fBeams[b].GetDataGraph(i), bestSinSqr, "");
//       leg->AddEntry(fBeams[b].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

  return;
}

// Only in case of NoMinuit
//----------------------------------------------------------------------
void NCExtrapolation::DrawAndWriteLazyContours()
{
  assert(fLazyContours.size() == 9);

  for(int n = 0; n < 3; ++n){
    TString canvasName, axisLabels;
    double xStart = 0.;
    double xEnd = 0.;
    double yStart = 0.;
    double yEnd = 0.;
    double xBest = 0.;
    double yBest = 0.;

    switch(n){
    case(0):
      canvasName = "contours-sinsq-dmsq";
      axisLabels = ";sin^{2}2\\theta;\\Delta m^{2}";
      xStart = NCType::kSinSqrStart; yStart = NCType::kDeltaMSqrStart*1e-3;
      xEnd   = NCType::kSinSqrEnd;   yEnd   = NCType::kDeltaMSqrEnd*1e-3;
      xBest  = fBestSinSqr;          yBest  = fBestDeltaMSqr;
      break;
    case(1):
      canvasName = "contours-fster-dmsq";
      axisLabels = ";f_{s};\\Delta m^{2}";
      xStart = NCType::kFSterileStart; yStart = NCType::kDeltaMSqrStart*1e-3;
      xEnd   = NCType::kFSterileEnd;   yEnd   = NCType::kDeltaMSqrEnd*1e-3;
      xBest  = fBestFSterile;          yBest  = fBestDeltaMSqr;
      break;
    case(2):
      canvasName = "contours-fster-sinsq";
      axisLabels = ";f_{s};sin^{2}\\theta";
      xStart = NCType::kFSterileStart; yStart = NCType::kSinSqrStart;
      xEnd   = NCType::kFSterileEnd;   yEnd   = NCType::kSinSqrEnd;
      xBest  = fBestFSterile;          yBest  = fBestSinSqr;
      break;
    }

    TLegend leg(0.75, 0.75, 0.95, 0.95);
    leg.SetBorderSize(0);

    TCanvas canv(canvasName);
    TH2F forceAxes ("", axisLabels, 10, xStart, xEnd, 10, yStart, yEnd);
    forceAxes.SetDirectory(0); // Don't write out useless axes to disk.
    forceAxes.Draw();
    TGraph* contGraph[3];
    for(int conf = 0; conf < 3; ++conf){
      const std::vector<Coord>& contVec = fLazyContours[conf*3+n];
      if(contVec.size() > 0){
        contGraph[conf] = new TGraph();
        for(unsigned int i = 0; i < contVec.size(); ++i)
          contGraph[conf]->SetPoint(i, contVec[i].x, contVec[i].y);
        const TString confText = (conf == 0) ? "68" : ( (conf == 1) ? "90" : "99");
        leg.AddEntry(contGraph[conf], confText+"% Confidence Level", "l");
        contGraph[conf]->SetName(canvasName+confText+"Graph");
        contGraph[conf]->Draw("l");
        //      contGraph[conf]->Write();
        gDirectory->Append(contGraph[conf]);
      }
    }

    TMarker bestFit(xBest, yBest, kFullStar);
    bestFit.SetMarkerSize(2);
    bestFit.SetMarkerColor(kBlack);
    bestFit.Draw("p same");

    leg.AddEntry(&bestFit, "Best Fit Point", "p");
    leg.Draw();

    canv.Write();
  }
}

//----------------------------------------------------------------------
void NCExtrapolation::WriteResources()
{

  MSG("NCExtrapolation", Msg::kInfo)
    << "NCExtrapolation::WriteResources() - "
    << NCType::kExtrapolationNames[fExtrapolationType]
    << endl;

  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* saveDir = gDirectory;
  TString plotName = "plots"+NCType::kExtrapolationNames[fExtrapolationType];
  TDirectory* plotsDir = saveDir->mkdir(plotName, plotName);
  plotsDir->cd();

  //write out the spectra etc from the beam objects
  for(uint i = 0; i < fBeams.size(); ++i){
    fBeams[i].MakeResultPlots();
    fBeams[i].WriteResources();
  }//end loop over beams

  if(fDeltaChiSqrDeltaMVsFSterile->GetEntries() > 0){
    fDeltaChiSqrDeltaMVsFSterile->Write();
    this->DrawContours(fDeltaChiSqrDeltaMVsFSterile, fBestFSterile, fBestDeltaMSqr);
  }

  if(fDeltaChiSqrSinSqrVsFSterile->GetEntries() > 0){
    fDeltaChiSqrSinSqrVsFSterile->Write();
    this->DrawContours(fDeltaChiSqrSinSqrVsFSterile, fBestFSterile, fBestSinSqr);
  }

  if(fDeltaChiSqrDeltaMVsSinSqr->GetEntries() > 0){
    fDeltaChiSqrDeltaMVsSinSqr->Write();
    this->DrawContours(fDeltaChiSqrDeltaMVsSinSqr, fBestSinSqr, fBestDeltaMSqr);
  }

  this->DrawBestFitSpectra(Detector::kNear);
  this->DrawBestFitSpectra(Detector::kFar);

  if(fLazyContours.size() > 0) DrawAndWriteLazyContours(); // Only in case of NoMinuit
  if(fContourGraphs.size() > 0) DrawAndWriteLazyContoursMinuit();

  saveDir->cd();

  return;
}

//----------------------------------------------------------------------
void NCExtrapolation::DrawAndWriteLazyContoursMinuit()
{
  // TODO - this function is now hideous:
  //   should eventually stop special-casing fster,sinsq and dmsq
  //   should not rely on the fact that we know fExtraVars* will be filled, need to make this an
  //     explicit part of the contract with the base class.
  for(unsigned int n = 0; n < fXaxes.size(); ++n){
    TString canvasName, axisLabels;
    double xStart = 0.;
    double xEnd = 0.;
    double yStart = 0.;
    double yEnd = 0.;
    double xBest = 0.;
    double yBest = 0.;

    switch(fXaxes[n]){
    case(0):
      canvasName = "contours-fster-";
      axisLabels = ";f_{s};";
      xStart = NCType::kFSterileStart;
      xEnd   = NCType::kFSterileEnd;
      xBest  = fBestFSterile;
      break;
    case(1):
      canvasName = "contours-sinsq-";
      axisLabels = ";sin^{2}2\\theta;";
      xStart = NCType::kSinSqrStart;
      xEnd   = NCType::kSinSqrEnd;
      xBest  = fBestSinSqr;
      break;
    case(2):
      canvasName += "contours-dmsq-";
      axisLabels += ";\\Delta m^{2};";
      xStart = NCType::kDeltaMSqrStart*1e-3;
      xEnd   = NCType::kDeltaMSqrEnd*1e-3;
      xBest  = fBestDeltaMSqr;
      break;
    default:
      canvasName = "contours-"+fExtraVarsNames[fXaxes[n]-3]+"-";
      axisLabels = ";"+fExtraVarsNames[fXaxes[n]-3]+";";
      xStart = fExtraVarsMin[fXaxes[n]-3];
      xEnd = fExtraVarsMax[fXaxes[n]-3];
      xBest = fMinCoord[fXaxes[n]];
    }
    switch(fYaxes[n]){
    case(0):
      canvasName += "fster";
      axisLabels += ";f_{s};";
      yStart = NCType::kFSterileStart;
      yEnd   = NCType::kFSterileEnd;
      yBest  = fBestFSterile;
      break;
    case(1):
      canvasName += "sinsq";
      axisLabels += "sin^{2}2\\theta";
      yStart = NCType::kSinSqrStart;
      yEnd   = NCType::kSinSqrEnd;
      yBest  = fBestSinSqr;
      break;
    case(2):
      canvasName += "dmsq";
      axisLabels += "\\Delta m^{2}";
      yStart = NCType::kDeltaMSqrStart*1e-3;
      yEnd   = NCType::kDeltaMSqrEnd*1e-3;
      yBest  = fBestDeltaMSqr;
      break;
    default:
      canvasName += fExtraVarsNames[fYaxes[n]-3];
      axisLabels += fExtraVarsNames[fYaxes[n]-3];
      yStart = fExtraVarsMin[fYaxes[n]-3];
      yEnd = fExtraVarsMax[fYaxes[n]-3];
      yBest = fMinCoord[fYaxes[n]];
    }

    TLegend leg(0.75, 0.75, 0.95, 0.95);
    leg.SetBorderSize(0);

    TCanvas canv(canvasName);
    TH2F forceAxes ("", axisLabels, 10, xStart, xEnd, 10, yStart, yEnd);
    forceAxes.SetDirectory(0); // Don't write out useless axes to disk.
    forceAxes.Draw();
    int trueconf = 0; // because some confidences were never written, need to keep track of the effective index.
    for(int conf = 0; conf < 3; ++conf){
      if(conf==0 && !fWant68) continue;
      if(conf==1 && !fWant90) continue;
      if(conf==2 && !fWant99) continue;
      TGraph* contGraph = fContourGraphs[trueconf*fXaxes.size()+n];
      if(contGraph && contGraph->GetN() > 0){
        const TString confText = (conf == 0) ? "68" : ( (conf == 1) ? "90" : "99");
        leg.AddEntry(contGraph, confText+"% Confidence Level", "l");
        contGraph->SetName(canvasName+confText+"Graph");
        contGraph->Draw("l");
        gDirectory->Append(contGraph);
      }
      ++trueconf;
    }

    TMarker bestFit(xBest, yBest, kFullStar);
    bestFit.SetMarkerSize(2);
    bestFit.SetMarkerColor(kBlack);
    bestFit.Draw("p same");

    leg.AddEntry(&bestFit, "Best Fit Point", "p");
    leg.Draw();

    canv.Write();
  }
}

//----------------------------------------------------------------------
//add says whether to add a new beam if this one isnt found, scalefactor is
//the pot normalization for data to mc
int NCExtrapolation::FindBeamIndex(Detector::Detector_t det,
                                   BeamType::BeamType_t bt,
                                   int rt,
                                   bool add,
                                   double dataPOT,
                                   double mcPOT)
{
  int index = 999;
  for(uint i = 0; i < fBeams.size(); ++i){
    if(fBeams[i].GetDetector() == det
       && fBeams[i].GetBeamType() == bt
       && fBeams[i].GetRunType() == rt) index = (int)i;
  }
  if(index > 100){
    if(add){
      fBeams.push_back(NCBeam(bt, det, rt, dataPOT, mcPOT, fUseCC, fUseNC));
      index = fBeams.size()-1;
      MSG("NCExtrapolation", Msg::kInfo) << "adding beam type "
                                         << BeamType::AsString(bt)
                                         << " run type "
                                         << rt
                                         << " to "
                                         << Detector::AsString(det)
                                         << " data pot = " << dataPOT
                                         << " mc pot = " << mcPOT
                                         << " index = " << index << endl;
    }//end if add beams when not found
//     else{
//       MAXMSG("NCExtrapolation", Msg::kWarning,1) << "could not find "
//                                               << Detector::AsString(det)
//                                               << " " << BeamType::AsString(bt)
//                                               << " " << rt
//                                               << " beam in vector - return index = "
//                                               << index << endl;
//     }//end if cant find beam

  }//end if beam not found

  return index;
}

//----------------------------------------------------------------------
//method to fill the result histograms based on the best fit to the
//oscillation parameters and penalty terms
void NCExtrapolation::FillResultHistograms(int beam)
{

  NCEnergyBin *bin       = 0;
  int mcSize             = 0;
  int mcSize_bkg         = 0;
  int mcSize_NuMuToNuTau = 0;
  int mcSize_NuEToNuE    = 0;
  int mcSize_NuMuToNuE   = 0;
  double trueEnergy      = 0.;
  double recoShowerE     = 0.;
  double recoMuonE       = 0.;
  double weight          = 0.;
  int    flavor          = 0;
  double survivalProb    = 1.;
  int sigType  = NCType::kNC;
  int bkgType  = NCType::kCC;
  double baseLine = NCType::kBaseLineFar;
  Detector::Detector_t det = fBeams[beam].GetDetector();
  if(det == Detector::kNear) baseLine = NCType::kBaseLineNear;

  CoordNDim coords = GetLazyBestFitPoint();

  // Alters the systematics 'weights' to reflect the current
  // value returned by the fitter.
  if( fExtraVarsNames.size() > 0 ){
    if( fLocNormalization > 0 ){
      fWeightNormalization = coords.at( fLocNormalization );
    }
    if( fLocRelativeHadronicCalibration > 0 ){
      fWeightRelativeHadronicCalibration = coords.at( fLocRelativeHadronicCalibration );
    }
    if( fLocAbsoluteHadronicCalibration > 0 ){
      fWeightAbsoluteHadronicCalibration = coords.at( fLocAbsoluteHadronicCalibration );
    }
    if( fLocTrackEnergy > 0 ){
      fWeightTrackEnergy = coords.at( fLocTrackEnergy );
    }
    if( fLocNCBkg > 0 ){
      fWeightNCBkg = coords.at( fLocNCBkg );
    }
    if( fLocCCBkg > 0 ){
      fWeightCCBkg = coords.at( fLocCCBkg );
    }
  }//end if fitting systematic nuisance parameters

  fShowerScale = 1.+fWeightRelativeHadronicCalibration;
  fShowerScale *= 1.+fWeightAbsoluteHadronicCalibration;
  fTrackScale  = 1.+fWeightTrackEnergy;

  //set oscillation parameters
  vector<double> pars;
  SetParametersForOscillationFit(coords, pars);

  //loop over NC/CC energy bins in beam
  for(int j = 0; j < 2; ++j){

    if(j == NCType::kNC && det == Detector::kFar && !fUseNC) continue;
    if(j == NCType::kCC && det == Detector::kFar && !fUseCC) continue;

    sigType = j;
    bkgType = NCType::kNC;

    if( j == NCType::kNC) bkgType = NCType::kCC;

    for(int i = 0; i < fBeams[beam].GetNumberEnergyBins(j); ++i){

      bin = fBeams[beam].GetEnergyBin(i, j);

      mcSize             = bin->GetMCSignalVectorSize();
      mcSize_bkg         = bin->GetMCBackgroundVectorSize();
      mcSize_NuMuToNuTau = bin->GetMCNuMuToNuTauVectorSize();
      mcSize_NuEToNuE    = bin->GetMCNuEToNuEVectorSize();
      mcSize_NuMuToNuE   = bin->GetMCNuMuToNuEVectorSize();

      for( int e = 0; e < mcSize; ++e ){
        bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
        survivalProb = FindSurvivalProbability(NCType::kNuMuToNuMu,
                                               sigType,
                                               NCType::kBaseLineFar,
                                               trueEnergy,
                                               pars);
        fBeams[beam].GetMCHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                   + recoMuonE*fTrackScale,
                                                   weight*(1.+fWeightNormalization));
        fBeams[beam].GetMCTrueHistogram(sigType)->Fill(trueEnergy,weight);

        fBeams[beam].GetMCFitHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                      + recoMuonE*fTrackScale,
                                                      weight*survivalProb*(1.+fWeightNormalization));
      }//end loop over signal

      for( int e = 0; e < mcSize_bkg; ++e ){
        bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
        survivalProb = FindSurvivalProbability(NCType::kNuMuToNuMu,
                                               bkgType,
                                               NCType::kBaseLineFar,
                                               trueEnergy,
                                               pars);
        fBeams[beam].GetMCBackgroundHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                             + recoMuonE*fTrackScale,
                                                             weight*survivalProb*(1.+fWeightNormalization));
        fBeams[beam].GetMCHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                   + recoMuonE*fTrackScale,
                                                   weight*(1.+fWeightNormalization));
        fBeams[beam].GetMCFitHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                      + recoMuonE*fTrackScale,
                                                      weight*survivalProb*(1.+fWeightNormalization));
      }//end loop over background

      //nutau
      for( int e = 0; e < mcSize_NuMuToNuTau; ++e ){
        bin->GetMCNuMuToNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
        survivalProb = FindSurvivalProbability(NCType::kNuMuToNuTau,
                                               NCType::kCC,
                                               NCType::kBaseLineFar,
                                               trueEnergy,
                                               pars);
        fBeams[beam].GetMCFitNuMuToNuTauHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                                 + recoMuonE*fTrackScale,
                                                                 weight*survivalProb*(1.+fWeightNormalization));
        fBeams[beam].GetMCFitHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                      + recoMuonE*fTrackScale,
                                                      weight*survivalProb*(1.+fWeightNormalization));
        fBeams[beam].GetMCNoFitNuMuToNuTauHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                                   + recoMuonE*fTrackScale,
                                                                   weight*(1.+fWeightNormalization));
      }//end loop over taus

      //nue
      for( int e = 0; e < mcSize_NuEToNuE; ++e ){
        bin->GetMCNuEToNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
        survivalProb = FindSurvivalProbability(NCType::kNuEToNuE,
                                               NCType::kCC,
                                               NCType::kBaseLineFar,
                                               trueEnergy,
                                               pars);
        fBeams[beam].GetMCFitHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                      + recoMuonE*fTrackScale,
                                                      weight*(1.+fWeightNormalization));
        fBeams[beam].GetMCNoFitNuEToNuEHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                                + recoMuonE*fTrackScale,
                                                                weight*(1.+fWeightNormalization));
        fBeams[beam].GetMCFitNuEToNuEHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                              + recoMuonE*fTrackScale,
                                                              weight*(1.+fWeightNormalization));
      }//end loop over electrons

      //nue appearance
      for( int e = 0; e < mcSize_NuMuToNuE; ++e ){
        bin->GetMCNuMuToNuEInformation(trueEnergy, recoShowerE,
                                       recoMuonE, weight, flavor, e);
        survivalProb = FindSurvivalProbability(NCType::kNuMuToNuE,
                                               NCType::kCC,
                                               NCType::kBaseLineFar,
                                               trueEnergy,
                                               pars);
        fBeams[beam].GetMCFitHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                      + recoMuonE*fTrackScale,
                                                      weight*survivalProb*(1.+fWeightNormalization));
        fBeams[beam].GetMCNoFitNuMuToNuEHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                                 + recoMuonE*fTrackScale,
                                                                 weight*(1.+fWeightNormalization));
        fBeams[beam].GetMCFitNuMuToNuEHistogram(sigType)->Fill(recoShowerE*fShowerScale
                                                               + recoMuonE*fTrackScale,
                                                               weight*survivalProb*(1.+fWeightNormalization));
      }//end loop over electron appearance

    }//end loop over energy bins

  }//end loop over nc/cc

  return;
}

//----------------------------------------------------------------------
void NCExtrapolation::SetParametersForOscillationFit(CoordNDim coords,
                                                     std::vector<double>& pars)
{
  pars.clear();
  pars.push_back(coords.at(NCType::kSinSqr));
  pars.push_back(coords.at(NCType::kDeltaMSqr));
  if(fLocPhi43 > 0)
    pars.push_back(coords.at(fLocPhi43));
  else pars.push_back(0.);
  if(fLocPMu > 0)
    pars.push_back(coords.at(fLocPMu));
  else pars.push_back(0.);
  pars.push_back(coords.at(NCType::kFSterile));
  if(fLocPE > 0)
    pars.push_back(coords.at(fLocPE));
  else pars.push_back(0.0);

  for(uint i = 0; i < pars.size(); ++i){
    if(pars[i] < 0. && TMath::Abs(pars[i]) < 1.e-3) pars[i] = 0.;
  }

  //3 Flavor mixing
//   pars.push_back(0.5*TMath::ASin(TMath::Sqrt(coords.at(NCType::kSinSqr))));
//   pars.push_back(coords.at(NCType::kDeltaMSqr));
//   pars.push_back(0.);
//   pars.push_back(0.61);//SNO+KamLAND
//   pars.push_back(7.59e-5);//SNO+KamLAND
//   pars.push_back(2.7);//g/cc standard rock
//   pars.push_back(0.);//why not?
//   pars.push_back(1.);//again, why not?

  return;
}

//----------------------------------------------------------------------
//method to find survival probability based on passed in parameters
//in pars vector.  the order of the parameters in pars depends on
//the function used to calculate the probability.
double NCExtrapolation::FindSurvivalProbability(int oscillationMode,
                                                int interactionType,
                                                double baseline,
                                                double trueEnergy,
                                                std::vector<double>& pars)
{
  double probability = 1.;

//   probability = fUtils->Find3FlavorProbability(oscillationMode,
//                                             interactionType,
//                                             baseline,
//                                             trueEnergy,
//                                             pars);

//   probability = fUtils->FindExtendedModelIProbability(oscillationMode,
//                                                    interactionType,
//                                                    baseline,
//                                                    trueEnergy,
//                                                    pars);

  probability = fUtils->FindParkeModelProbability(oscillationMode,
                                                  interactionType,
                                                  baseline,
                                                  trueEnergy,
                                                  pars);

  return probability;
}

//----------------------------------------------------------------------
void NCExtrapolation::CombineRuns(Detector::Detector_t det)
{
  //make a map of the number of run types for each beam condition
  map<BeamType::BeamType_t, int> numRuns;
  map<BeamType::BeamType_t, int>::iterator itr = numRuns.begin();

  int totIndex = 0;
  int beamIndex = 0;

  for(uint i = 0; i < fBeams.size(); ++i){
    if(fBeams[i].GetDetector() == det){
      itr = numRuns.find(fBeams[i].GetBeamType());
      if(itr != numRuns.end()) itr->second++;
      else numRuns[fBeams[i].GetBeamType()] = 1;
      MSG("NCExtrapolation", Msg::kInfo) << "beam "
                                         << BeamType::AsString(fBeams[i].GetBeamType())
                                         << " in " << Detector::AsString(det)
                                         << " has " << numRuns[fBeams[i].GetBeamType()]
                                         << " runs " << endl;

    }//end if the same detector as what i want
  }//end loop over beams

  //now figure out how many beam objects to add
  for(itr = numRuns.begin(); itr != numRuns.end(); ++itr){

    //see if there is more than one run for this beam type,
    //and if so create a total beam for it, then get the individual
    //runs and add them to the total
    if(itr->second > 1){
      totIndex = FindBeamIndex(det, itr->first, NCType::kRunAll, true, 0., 0.);
      for(int i = NCType::kRunI; i < NCType::kMaxRun+1; ++i){
        beamIndex = FindBeamIndex(det, itr->first, i);
        if(beamIndex < 100){
          fBeams[totIndex].Add(fBeams[beamIndex]);
          MSG("NCExtrapolation", Msg::kInfo) << "combining runs for beam "
                                             << BeamType::AsString(itr->first)
                                             << " in " << Detector::AsString(det)
                                             << endl;
        }//end if run type exists
      }//end loop over run types
    }//end if more than one run for the beam
    else{
      for(int i = NCType::kRunI; i < NCType::kMaxRun+1; ++i){
        beamIndex = FindBeamIndex(det, itr->first, i);
        if(beamIndex < 100){
          fBeams[beamIndex].SetRunType(NCType::kRunAll);
          MSG("NCExtrapolation", Msg::kInfo) << "setting run " << i
                                             << " to "  << NCType::kRunAll
                                             << " for beam "
                                             << BeamType::AsString(itr->first)
                                             << " in " << Detector::AsString(det)
                                             << endl;
        }//end if the run for this beam exists
      }//end loop over runs
    }//end if only one run for the beam
  }//end loop over map of beam type to # of runs

  return;
}

//----------------------------------------------------------------------
//make sure the fMinChiSqr variable is set before calling this method
//assumes the vector of vector of vectors of doubles has the indicies
//[Delta m^2][f_sterile][sin^2(2theta)]
void NCExtrapolation::FillDeltaChiSqrMaps(std::vector< std::vector< std::vector<double> > > &chiSqrPnts)
{

  //loop over the parameter space values of chi^2 to save the delta chi^2 to
  //the histograms
  double minSinSqrChiSqr = 1.e10;
  for (int dm = 0; dm < NCType::kNumDeltaMSqrBins; dm++) {
    float dmsq = (1.*dm + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;

    // fsterile loop
    for (int fs = 0; fs < NCType::kNumFSterileBins; fs++) {
      float fster = (1.*fs + 0.5)*NCType::kDeltaFSterile+NCType::kFSterileStart;
      minSinSqrChiSqr = 1.e10;

      // sin2theta loop
      for (int st = 0; st < NCType::kNumSinSqrBins; st++) {

        if(chiSqrPnts[dm][fs][st] < minSinSqrChiSqr)
          minSinSqrChiSqr = chiSqrPnts[dm][fs][st];

      }//end s2t

      //found the minimum in sin^2 2theta for this point in delta m^2 vs f_s
      fDeltaChiSqrDeltaMVsFSterile->Fill(fster, dmsq, minSinSqrChiSqr-fMinChiSqr);

    }//end fs
  }//end dmsq

  double minDeltaMSqrChiSqr = 1.e10;
  // sin2theta loop
  for (int st=0;st<NCType::kNumSinSqrBins;st++) {
    float s2t=(1.*st + 0.5)*NCType::kDeltaSinSqr+NCType::kSinSqrStart;

    // fsterile loop
    for (int fs=0;fs<NCType::kNumFSterileBins;fs++) {
      float fster=(1.*fs + 0.5)*NCType::kDeltaFSterile+NCType::kFSterileStart;

      minDeltaMSqrChiSqr = 1.e10;
      for (int dm=0;dm<NCType::kNumDeltaMSqrBins;dm++) {

        if(chiSqrPnts[dm][fs][st] < minDeltaMSqrChiSqr)
          minDeltaMSqrChiSqr = chiSqrPnts[dm][fs][st];

      }//end s2t

      //found the minimum in sin^2 2theta for this point in delta m^2 vs f_s
      fDeltaChiSqrSinSqrVsFSterile->Fill(fster, s2t, minDeltaMSqrChiSqr-fMinChiSqr);

    }//end fs
  }//end dmsq

  double minFSterileChiSqr = 1.e10;
  for (int dm=0;dm<NCType::kNumDeltaMSqrBins;dm++) {
    float dmsq=(1.*dm + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;

    // sin2theta loop
    for (int st=0;st<NCType::kNumSinSqrBins;st++) {
      float s2t=(1.*st + 0.5)*NCType::kDeltaSinSqr+NCType::kSinSqrStart;
      minFSterileChiSqr = 1.e10;

      // fsterile loop
      for (int fs=0;fs<NCType::kNumFSterileBins;fs++) {

        if(chiSqrPnts[dm][fs][st] < minFSterileChiSqr)
          minFSterileChiSqr = chiSqrPnts[dm][fs][st];

      }//end s2t

      //found the minimum in sin^2 2theta for this point in delta m^2 vs f_s
      fDeltaChiSqrDeltaMVsSinSqr->Fill(s2t, dmsq, minFSterileChiSqr-fMinChiSqr);

    }//end fs
  }//end dmsq

  return;
}


// ---------------------------------------------------------------------
// Make sure your derived class has overridden GetChiSqrForMap() and that
// everything is ready, and then call this function. It finds the best-fit
// point, and contours at 3 levels of significance, projected onto 3 different
// axes. Sets fMinChiSqr, fBestSinSqr, fBestDeltaMSqr, fBestFSterile.
void NCExtrapolation::MakeDeltaChiSqrMaps(CoordNDim extraVarsMin,
                                          CoordNDim extraVarsMax,
                                          CoordNDim extraVarsPrecision,
                                          std::vector<std::string> extraVarsNames,
                                          std::vector<int> xAxes,
                                          std::vector<int> yAxes)
{
  assert(xAxes.size() == yAxes.size());
  fXaxes = xAxes;
  fYaxes = yAxes;

  switch(fLazyContoursUseMinuit){
  default:
    MSG("NCExtrapolation", Msg::kWarning) << "LazyContoursUseMinuit not set. Assuming zero.\n";
    // fall through
  case(0):
    MakeDeltaChiSqrMapsNoMinuit(extraVarsMin, extraVarsMax, extraVarsPrecision, extraVarsNames, fWant68, fWant90, fWant99);
    break;
  case(1):
    MakeDeltaChiSqrMapsHybrid(extraVarsMin,
                              extraVarsMax,
                              extraVarsPrecision,
                              extraVarsNames);
    break;
  case(2):
    MakeDeltaChiSqrMapsMinuit(extraVarsMin, extraVarsMax, extraVarsNames, fWant68, fWant90, fWant99);
    break;
  }
}

// ---------------------------------------------------------------------
// Make sure your derived class has overridden GetChiSqrForMap() and that
// everything is ready, and then call this function. It finds the best-fit
// point, and contours at 3 levels of significance, projected onto 3 different
// axes. Sets fMinChiSqr, fBestSinSqr, fBestDeltaMSqr, fBestFSterile.
void NCExtrapolation::MakeDeltaChiSqrMapsNoMinuit(CoordNDim extraVarsMin,
                                                  CoordNDim extraVarsMax,
                                                  CoordNDim extraVarsPrecision,
                                                  std::vector<std::string> extraVarsNames,
                                                  bool want68,
                                                  bool want90,
                                                  bool want99)
{
  assert(extraVarsMin.size() == extraVarsMax.size() &&
         extraVarsMax.size() == extraVarsPrecision.size());

  enum{FSTER, SINSQ, DMSQ}; // Is this already defined somewhere else?

  // Use NCType to get the range of the oscillation parameters.
  CoordNDim start(3);
  start[FSTER] = NCType::kFSterileStart;
  start[SINSQ] = NCType::kSinSqrStart;
  start[DMSQ]  = NCType::kDeltaMSqrStart*1e-3;
  // Add the ranges of the additional parameters we were passed.
  start.insert(start.end(), extraVarsMin.begin(), extraVarsMin.end());

  CoordNDim end(3);
  end[FSTER] = NCType::kFSterileEnd;
  end[SINSQ] = NCType::kSinSqrEnd;
  end[DMSQ]  = NCType::kDeltaMSqrEnd*1e-3;
  end.insert(end.end(), extraVarsMax.begin(), extraVarsMax.end());

  CoordNDim minimPrec(3);
  minimPrec[FSTER] = NCType::kDeltaFSterile;
  minimPrec[SINSQ] = NCType::kDeltaSinSqr;
  minimPrec[DMSQ]  = NCType::kDeltaDeltaMSqr*1e-3;
  minimPrec.insert(minimPrec.end(),
                   extraVarsPrecision.begin(), extraVarsPrecision.end());

  std::vector<std::string> allNames(3);
  allNames[FSTER] = "fster";
  allNames[SINSQ] = "sinsq";
  allNames[DMSQ] = "dmsq";
  allNames.insert(allNames.end(), extraVarsNames.begin(), extraVarsNames.end());
  while(allNames.size() < start.size()) allNames.push_back("[unknown]");

  // Create functionoid that redirects to derived-class's implementation.
  // FIXME - currently always using unphysical penalties.
  GetChiSqrFromDerived chiSqrFunc(this, start, end, allNames, true);

  // Store the best-fit point in here.
  CoordNDim minCoord;

  // Find the best fit point.
  fMinChiSqr = FindMinNDim(chiSqrFunc, start, end, minCoord, minimPrec);

  MSG("NCExtrapolation", Msg::kInfo) << "Best fit : f_s = " << minCoord[FSTER]
                                             << " sin^2 = " << minCoord[SINSQ]
                                              << " dm^2 = " << minCoord[DMSQ];
  if(minCoord.size() > 3){
    MSG("NCExtrapolation", Msg::kInfo) << endl
                                       << " Additional parameters:";
    for(unsigned int n = 3; n < minCoord.size(); ++n)
      MSG("NCExtrapolation", Msg::kInfo) << allNames.at(n) << " = "
                                         << minCoord[n] << " ";
  }
  MSG("NCExtrapolation", Msg::kInfo) << endl;

  // The three functions we will be plotting the contours of.
  GetFuncMinWRTSomeVars<GetChiSqrFromDerived> minWRT[3];

  // For each of the oscillation parameters we will minimize over:
  for(int minVar = 0; minVar < 3; ++minVar){
    // Store its index, range and precision.
    std::vector<int> toMin(1);
    toMin[0] = minVar;
    CoordNDim min(1), max(1), prec(1);
    min[0] = start[minVar];
    max[0] = end[minVar];
    prec[0] = minimPrec[minVar];

    // And then append the indices, ranges and precisions of all
    // additional variables we will minimize over.
    for(unsigned int n = 0; n < extraVarsMin.size(); ++n){
      toMin.push_back(3+n);
      min.push_back(extraVarsMin[n]);
      max.push_back(extraVarsMax[n]);
      prec.push_back(extraVarsPrecision[n]);
    }
    // Create the functionoid that minimizes over all of these.
    minWRT[minVar] = GetFuncMinWRTSomeVars<GetChiSqrFromDerived>
                       (chiSqrFunc, toMin, min, max, prec);
  }

  // QUESTION: are these the right values for a fit with more variables?
  const double contVals[] = {2.3, 4.6, 9.2};

  // The axes that should be used when minimizing over each of the three variables.
  const int xAxis[] = {SINSQ, FSTER, FSTER};
  const int yAxis[] = {DMSQ, DMSQ, SINSQ};
  // And human-readable form.
  const TString text[] = {"sin^2 and dm^2", "f_s and dm^2", "f_s and sin^2"};

  // For each of the 3 confidence levels (68, 90, 99)
  for(int conf = 0; conf < 3; ++conf){
    // If we don't want to draw contours at this confidence, add empty dummy contours.
    if((conf == 0 && !want68) || (conf == 1 && !want90) || (conf == 2 && !want99)){
      for(int i = 0; i < 3; ++i) fLazyContours.push_back(std::vector<Coord>());
      continue;
    }

    // For each of the 3 oscillation parameters we are minimizing.
    for(int minVar = 0; minVar < 3; ++minVar){
      MSG("NCExtrapolation", Msg::kInfo) << "Building contour in " << text[minVar]
                                         << " (up = " << contVals[conf] << ")...\n";
      // Look up which axes we should be using.
      const int x = xAxis[minVar];
      const int y = yAxis[minVar];
      const std::vector<Coord> contour = FindContour(minWRT[minVar],
                                                     //Coord(start[x], start[y]),
                                                     //Coord(end[x], end[y]),
                                                     Coord(minCoord[x], minCoord[y]),
                                                     fMinChiSqr+contVals[conf],
                                                     Coord(minimPrec[x], minimPrec[y]));
      fLazyContours.push_back(contour);
      MSG("NCExtrapolation", Msg::kInfo) << "there are " << contour.size()
                                         << " points on this contour" << endl;
    }
  }

  // The best-fit parameters are stored as member variables for later use.
  fBestFSterile = minCoord[FSTER];
  fBestSinSqr = minCoord[SINSQ];
  fBestDeltaMSqr = minCoord[DMSQ];
}

// ---------------------------------------------------------------------
// Make sure your derived class has overridden GetChiSqrForMap() and that
// everything is ready, and then call this function. It finds the best-fit
// point, and contours at 3 levels of significance, projected onto 3 different
// axes. Sets fMinChiSqr, fBestSinSqr, fBestDeltaMSqr, fBestFSterile.
void NCExtrapolation::MakeDeltaChiSqrMapsHybrid(CoordNDim extraVarsMin,
                                                CoordNDim extraVarsMax,
                                                CoordNDim extraVarsPrecision,
                                                std::vector<std::string> extraVarsNames)
{
  assert(extraVarsMin.size() == extraVarsMax.size() &&
         extraVarsMax.size() == extraVarsPrecision.size());

  enum{FSTER, SINSQ, DMSQ}; // Is this already defined somewhere else?

  // Use NCType to get the range of the oscillation parameters.
  CoordNDim start(3);
  start[FSTER] = NCType::kFSterileStart;
  start[SINSQ] = NCType::kSinSqrStart;
  start[DMSQ]  = NCType::kDeltaMSqrStart*1e-3;
  // Add the ranges of the additional parameters we were passed.
  start.insert(start.end(), extraVarsMin.begin(), extraVarsMin.end());

  CoordNDim end(3);
  end[FSTER] = NCType::kFSterileEnd;
  end[SINSQ] = NCType::kSinSqrEnd;
  end[DMSQ]  = NCType::kDeltaMSqrEnd*1e-3;
  end.insert(end.end(), extraVarsMax.begin(), extraVarsMax.end());

  CoordNDim prec(3);
  prec[FSTER] = NCType::kDeltaFSterile;
  prec[SINSQ] = NCType::kDeltaSinSqr;
  prec[DMSQ]  = NCType::kDeltaDeltaMSqr*1e-3;
  prec.insert(prec.end(), extraVarsPrecision.begin(), extraVarsPrecision.end());

  std::vector<std::string> allNames(3);
  allNames[FSTER] = "fster";
  allNames[SINSQ] = "sinsq";
  allNames[DMSQ] = "dmsq";
  allNames.insert(allNames.end(), extraVarsNames.begin(), extraVarsNames.end());
  while(allNames.size() < start.size()) allNames.push_back("[unknown]");

  // Create functionoid that redirects to derived-class's implementation.
  GetChiSqrFromDerived chiSqrFunc(this, start, end, allNames, !fAllowBestFitOutsideRange);


  NCExtrapolationMinuit<GetChiSqrFromDerived> minu(chiSqrFunc, 3+extraVarsMin.size());
  minu.SetPrintLevel(-1); // Shut Minuit up.
  minu.Command("set strategy 2"); // Be careful, instead of fast.

  for(int n = 0; n < int(start.size()); ++n)
    // NB - no parameter limits imposed.
    minu.DefineParameter(n, "", (start[n]+end[n])/2, (end[n]-start[n])/4, 0, 0);

  // Minuit tries to return a lot of parameters we don't want to keep.
  double junkd; int junki; TString junks;

  if(minu.Migrad()){
    MSG("NCExtrapolation", Msg::kWarning) << "Migrad failed, falling back to Simplex.\n";
    minu.SetPrintLevel(0); // Things are going wrong, we should let the user see.
    minu.mnsimp();
    minu.Migrad();
  }

  minu.mnstat(fMinChiSqr, junkd, junkd, junki, junki, junki);

  fMinCoord.resize(start.size());
  for(unsigned int n = 0; n < start.size(); ++n)
    minu.GetParameter(n, fMinCoord[n], junkd);

  // The best-fit parameters are stored as member variables for later use.
  fBestFSterile = fMinCoord[FSTER];
  fBestSinSqr = fMinCoord[SINSQ];
  fBestDeltaMSqr = fMinCoord[DMSQ];

  MSG("NCExtrapolation", Msg::kInfo) << "Best fit :" << endl;
  for(unsigned int n = 0; n < allNames.size(); ++n)
    MSG("NCExtrapolation", Msg::kInfo) << "  " << allNames[n] << " = " << fMinCoord[n] << endl;


  // The axes that should be used when minimizing over each of the three variables.
  //  const int xAxis[] = {SINSQ, FSTER, FSTER};
  //  const int yAxis[] = {DMSQ, DMSQ, SINSQ};

  chiSqrFunc.UsePenalty(!fAllowContoursOutsideRange);

  // These represent the 3 confidence levels 68%, 90%, 99%
  const double contVals[] = {2.3, 4.6, 9.2};

  for(int conf = 0; conf < 3; ++conf){
    for(unsigned int axisPair = 0; axisPair < fXaxes.size(); ++axisPair){
      if(conf == 0 && !fWant68) continue;
      if(conf == 1 && !fWant90) continue;
      if(conf == 2 && !fWant99) continue;

      // Look up which axes we should be using.
      const int x = fXaxes[axisPair];
      const int y = fYaxes[axisPair];

      MSG("NCExtrapolation", Msg::kInfo) << "Building contour in "
                                         << allNames[x]
                                         << " and " << allNames[y]
                                         << " (up = " << contVals[conf] << ")...\n";

      GetFuncMinWRTSomeVarsHybrid<GetChiSqrFromDerived> func(chiSqrFunc,
                                                             x, y,
                                                             start, end);
      const std::vector<Coord> contour = FindContour(func,
                                                     Coord(fMinCoord[x], fMinCoord[y]),
                                                     fMinChiSqr+contVals[conf],
                                                     Coord(prec[x], prec[y]));
      TGraph* g = new TGraph(contour.size());
      for(unsigned int n = 0; n < contour.size(); ++n)
        g->SetPoint(n, contour[n].x, contour[n].y);
      fContourGraphs.push_back(g);

      MSG("NCExtrapolation", Msg::kInfo) << "there are " << contour.size()
                                         << " points on this contour" << endl;
    } // end for minVar
  } // end for conf
}


// ---------------------------------------------------------------------
// Make sure your derived class has overridden GetChiSqrForMap() and that
// everything is ready, and then call this function. It finds the best-fit
// point, and contours at 3 levels of significance, projected onto 3 different
// axes. Sets fMinChiSqr, fBestSinSqr, fBestDeltaMSqr, fBestFSterile.
void NCExtrapolation::MakeDeltaChiSqrMapsMinuit(CoordNDim extraVarsMin,
                                                CoordNDim extraVarsMax,
                                                std::vector<std::string> extraVarsNames,
                                                bool want68,
                                                bool want90,
                                                bool want99)
{
  class NCExtrapolationMinuit: public TMinuit
  {
  public:
    NCExtrapolationMinuit(NCExtrapolation* pDeriv, CoordNDim r1, CoordNDim r2, std::vector<std::string> n)
      :pDerived(pDeriv), start(r1), end(r2), penalty(true), names(n)
    {
      assert(start.size() == end.size());
    }
    // Override TMinuit's Eval function - so that this gets called as the
    // minimization function. Redirects to GetChiSqrForMap in an NCExtrapolation
    // derived class - and applies quadratic penalties around the physical
    // region if asked to.
    virtual Int_t Eval(Int_t /*npar*/, Double_t* /*grad*/, Double_t& fval, Double_t* par, Int_t flag)
    {
      // If we are being asked to do something that isn't just
      // calculate the function then we don't know how...
      assert(flag == 4);

      CoordNDim r(start.size());
      for(unsigned int n = 0; n < start.size(); ++n) r[n] = par[n];

      // If the point is outside the physical region - move it to the boundary
      if(penalty)
        for(unsigned int n = 0; n < start.size(); ++n){
          if(par[n] < start[n]) r[n] = start[n];
          if(par[n] > end[n]) r[n] = end[n];
        }

      fval = pDerived->GetChiSqrForMap(r, names);

      // If we had to move the point then apply a quadratic penalty term.
      // The factor of 1000 was arrived at through trial and error.
      if(penalty)
        for(unsigned int n = 0; n < start.size(); ++n){
          const double wsqr = (start[n]-end[n])*(start[n]-end[n]);
          if(par[n] < start[n]) fval+=1e3*(start[n]-par[n])*(start[n]-par[n])/wsqr;
          if(par[n] > end[n]) fval+=1e3*(par[n]-end[n])*(par[n]-end[n])/wsqr;
        }

      return 0;
    }
    // The penalty terms are required to keep the best-fit point within physical limits.
    // But they interact badly with the contour finder - so turn them off before that.
    // CURRENTLY UNUSED
    void SetPenalty(bool p){penalty = p;}
  protected:
    NCExtrapolation* pDerived;
    int numberOfParameters;
    CoordNDim start, end;
    bool penalty;
    std::vector<std::string> names;
  };


  assert(extraVarsMin.size() == extraVarsMax.size());

  enum{FSTER, SINSQ, DMSQ}; // Is this already defined somewhere else?

  // Use NCType to get the range of the oscillation parameters.
  // We shouldn't be using the dm^2 limits as they aren't physical
  // limits - keep them in because we need them to aim the  fitter
  // in the right general direction.
  CoordNDim start(3);
  start[FSTER] = NCType::kFSterileStart;
  start[SINSQ] = NCType::kSinSqrStart;
  start[DMSQ]  = NCType::kDeltaMSqrStart*1e-3;
  // Add the ranges of the additional parameters we were passed.
  start.insert(start.end(), extraVarsMin.begin(), extraVarsMin.end());

  CoordNDim end(3);
  end[FSTER] = NCType::kFSterileEnd;
  end[SINSQ] = NCType::kSinSqrEnd;
  end[DMSQ]  = NCType::kDeltaMSqrEnd*1e-3;
  end.insert(end.end(), extraVarsMax.begin(), extraVarsMax.end());

  std::vector<std::string> allNames(3);
  allNames[FSTER] = "fster";
  allNames[SINSQ] = "sinsq";
  allNames[DMSQ] = "dmsq";
  allNames.insert(allNames.end(), extraVarsNames.begin(), extraVarsNames.end());
  while(allNames.size() < start.size()) allNames.push_back("[unknown]");

  NCExtrapolationMinuit minu(this, start, end, allNames);
  minu.SetPrintLevel(-1); // Shut Minuit up.
  minu.Command("set strategy 2"); // Be careful, instead of fast.

  for(int n = 0; n < int(3+extraVarsMin.size()); ++n)
    // NB - no parameter limits imposed.
    minu.DefineParameter(n, "", (start[n]+end[n])/2, (end[n]-start[n])/4, 0, 0);

  // Minuit tries to return a lot of parameters we don't want to keep.
  double junkd; int junki; TString junks;

  if(minu.Migrad()){
    MSG("NCExtrapolation", Msg::kWarning) << "Migrad failed, falling back to Simplex.\n";
    minu.SetPrintLevel(0); // Things are going wrong, we should let the user see.
    minu.mnsimp();
    minu.Migrad();
  }

  minu.mnstat(fMinChiSqr, junkd, junkd, junki, junki, junki);

  minu.GetParameter(0, fBestFSterile, junkd);
  minu.GetParameter(1, fBestSinSqr, junkd);
  minu.GetParameter(2, fBestDeltaMSqr, junkd);

  MSG("NCExtrapolation", Msg::kInfo) << "Best fit : f_s = " << fBestFSterile
                                             << " sin^2 = " << fBestSinSqr
                                              << " dm^2 = " << fBestDeltaMSqr;
  if(extraVarsMin.size() > 0){
    MSG("NCExtrapolation", Msg::kInfo) << endl
                                       << " Additional parameters:";
    for(unsigned int n = 3; n < 3+extraVarsMin.size(); ++n){
      double val;
      minu.GetParameter(n, val, junkd);
      MSG("NCExtrapolation", Msg::kInfo) << allNames.at(n) << " "
                                         << val << " ";
    }
  }
  MSG("NCExtrapolation", Msg::kInfo) << endl;

  // The edge penalties break the contour finder so turn them off.
  // It causes wobbles in the contours - but seemingly not too bad.
  // And it allows us to get the right results re physical limits,
  // none of the other methods I tried work.
  //  minu.SetPenalty(false);

  // QUESTION: are these the right values for a fit with more variables?
  const double contVals[] = {2.3, 4.6, 9.2};

  // The axes that should be used when minimizing over each of the three variables.
  const int xAxis[] = {SINSQ, FSTER, FSTER};
  const int yAxis[] = {DMSQ, DMSQ, SINSQ};
  // And human-readable form.
  const TString text[] = {"sin^2 and dm^2", "f_s and dm^2", "f_s and sin^2"};

  const int numContourPoints = 300;

  for(int c = 0; c < 3; ++c){
    for(int n = 0; n < 3; ++n){
      if( (c == 0 && !want68) || (c == 1 && !want90) || (c == 2 && !want99) ){
        fContourGraphs.push_back(new TGraph);
        continue;
      }
      MSG("NCExtrapolation", Msg::kInfo) << "Building contour in " << text[n]
                                         << " (up = " << contVals[c] << ")...\n";

      // Apply physical limits to those parameters being minimized over.
      // for(int p = 0; p < int(3+extraVarsMin.size()); ++p){
      //   if(n != xAxis[n] && n != yAxis[n]){
      //     // Passing the axis number as a double is bad - but unavoidable.
      //     double paramList[] = {p, start[p], end[p]};
      //     minu.mnexcm("set limits", paramList, 3, junki);
      //   }
      // }

      minu.SetErrorDef(contVals[c]);
      TGraph* g = (TGraph*)minu.Contour(numContourPoints, xAxis[n], yAxis[n]);
      fContourGraphs.push_back(g);
      if(minu.GetStatus()) MSG("NCExtrapolation", Msg::kWarning) << "Couldn't find contour.\n";

      //      minu.Command("set limits"); // Remove all parameter limits.
    }
  }
}

---