NCFitMaster.cxx

Go to the documentation of this file.

//$Id: NCFitMaster.cxx,v 1.21 2009/06/18 18:12:22 bckhouse Exp $
//
//C. Backhouse 11/2008
//

#include "NCFitMaster.h"

#include "NCUtils/Extrapolation/NCExtrapolation.h"

#include "MessageService/MsgService.h"

#include "TCanvas.h"
#include "TDirectory.h"
#include "TLegend.h"
#include "TMarker.h"
#include "TStopwatch.h"

CVSID("$Id: NCFitMaster.cxx,v 1.21 2009/06/18 18:12:22 bckhouse Exp $");


using namespace std;
using namespace NC::Fitter;


namespace NC{
  class GetChiSqrFromDerived: public ICallableND
  {
  public:
    GetChiSqrFromDerived(){}
    GetChiSqrFromDerived(NCExtrapolation* P,
                         std::vector<NCParameter> pars,
                         bool penalty,
                         bool analyticNorm)
      :p(P), params(pars), usePenalty(penalty),
       fAnalyticNormalization(analyticNorm){}
    virtual ~GetChiSqrFromDerived(){} // shut gcc up
    double EvalAtEx(const CoordNDim& r, CoordNDim* ret) const;
    void UsePenalty(bool p){usePenalty = p;}
    void DrawAndWriteGraph() const;
  private:
    NCExtrapolation* p;
    // Cache must be static, otherwise each instance will have a seperate copy.
    static std::map<CoordNDim, double> cache;
    std::vector<NCParameter> params;
    bool usePenalty;
    bool fAnalyticNormalization;
    static std::vector<CoordNDim> evaledAt;
    static std::vector<double> chisqRets;
  };
} // end namespace


map<CoordNDim, double> NC::GetChiSqrFromDerived::cache;
vector<CoordNDim> NC::GetChiSqrFromDerived::evaledAt;
vector<double> NC::GetChiSqrFromDerived::chisqRets;

//----------------------------------------------------------------------
double NC::GetChiSqrFromDerived::EvalAtEx(const CoordNDim& r,
                                          CoordNDim* retCoord) const
{
  if(retCoord) *retCoord = r;

  evaledAt.push_back(r);

  const map<CoordNDim, double>::const_iterator it = cache.find(r);
  if(it != cache.end()){
//    cout << "chisq from cache " << r << endl;
    chisqRets.push_back(it->second);
    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();
  //cout << "UNCACHED CHISQ REQUESTED " << r << endl;

  CoordNDim evalAt = r;

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

  const CoordinateConverter* cc = p->GetCoordConv();
  NC::OscProb::OscPars* oscPars = cc->OscParsFromCoordNDim(evalAt);

  double norm = 1;
  double ret = p->FindChiSqrForPars(oscPars,
                                    cc->SystParsFromCoordNDim(evalAt),
                                    fAnalyticNormalization ? &norm : 0);

  // TODO WE NEED TO STORE NORM INTO RETCOORD, BUT DON'T KNOW WHICH ONE!

  delete oscPars;

  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.
    const double amplitude = 1e3;
    for(unsigned int n = 0; n < params.size(); ++n){
      if(params[n].UseLimits()){
        const double wsqr = (params[n].Min()-params[n].Max())*(params[n].Min()-params[n].Max());
        assert(wsqr);
        if(r[n] < params[n].Min())
          ret += amplitude*(params[n].Min()-r[n])*(params[n].Min()-r[n])/wsqr;
        if(r[n] > params[n].Max())
          ret += amplitude*(r[n]-params[n].Max())*(r[n]-params[n].Max())/wsqr;
      }
    }
  }

  //  for(unsigned int n = 0; n < evalAt.size(); ++n)
  //    cout << names[n] << " = " << evalAt[n] << "  ";
  //  cout << "-> " << ret << endl;

  // TODO - this cache is broken when using more than one extrapolation
  // disable it for now
  //  cache[r] = ret;
  chisqRets.push_back(ret);

  return ret;
}

//----------------------------------------------------------------------
void NC::GetChiSqrFromDerived::DrawAndWriteGraph() const
{
  TCanvas c("minimizer_eval_positions");
  TH2D axes("", "Evaluation Positions;time;", 100, 0, evaledAt.size(), 100, -1, 5);
  axes.Draw();
  for(unsigned int m = 0; m < evaledAt[0].size(); ++m){
    TGraph* g = new TGraph;
    g->Set(evaledAt.size());
    g->SetMarkerColor(m+1);
    for(unsigned int n = 0; n < evaledAt.size(); ++n){
      g->SetPoint(n, n, evaledAt[n][m]);
      g->Draw("p same");
    }
  }

  c.Write();
}


//----------------------------------------------------------------------
NC::FitMaster::FitMaster(NCExtrapolation* child)
  : fExtrap(child), fBestFitOscPars(0), fFixedValsForMarg(0)
{
  // We need to teach MinFinderSimple about where our variables are degenerate
  vector<string> names;
  names.push_back("theta23");
  names.push_back("theta");
  names.push_back("umu3");
  vector<double> poss;
  poss.push_back(TMath::Pi()/4);
  poss.push_back(TMath::Pi()/4);
  poss.push_back(0.5);
  NC::Fitter::MinFinderSimple::SetSpacePartitions(names, poss);
}


//----------------------------------------------------------------------
const Registry& NC::FitMaster::DefaultConfig()
{
  static Registry r;

  r.UnLockValues();

  // Whether the contour finder should use Minuit and if so how much (0 to 2).
  r.Set("LazyContoursUseMinuit", 0);

  r.Set("LazyContoursWant68",        true);
  r.Set("LazyContoursWant90",        true);
  r.Set("LazyContoursWant99",        true);
  r.Set("AllowBestFitOutsideRange",  false);
  r.Set("AllowContoursOutsideRange", false);
  r.Set("WantMinimum",               true);
  r.Set("Want2DProjections",         false);
  r.Set("WantGridSearch",            false);
  r.Set("WantSlices",                false);
  r.Set("Want2DSlices",              false);
  r.Set("WantCovarianceMatrix",      false);
  r.Set("TryPreviousMinimum",        false);
  r.Set("ForceBestFitToTruth",       false);
  r.Set("PrecScale",                 1.0);
  r.Set("AnalyticNormalization",     false);
  r.Set("ContoursList",              "");
  r.Set("ProjectionsList",           "");

  r.Merge(NC::CoordinateConverter::DefaultConfig());

  r.LockValues();

  return r;
}


//----------------------------------------------------------------------
void NC::FitMaster::Prepare(const Registry& r)
{
  MSG("NCFitMaster", Msg::kInfo) << "NC::FitMaster::Prepare(const Registry& r)" << endl;

  int         tmpb;  // a temp bool.
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string

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

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

  /*?
  MSG("NCExtrapolation", Msg::kInfo) << "PARAMETERS WE ARE USING:" << endl;
  for(int i = 0; i < NCType::kNumSystematicParameters; ++i){
    MSG("NCExtrapolation", Msg::kInfo) << i << ") "
                                       << NCType::kParams[i].name
                                       << " = " << fSystematicsToUse[i]
                                       << " with default = "
                                       << NCType::kParams[i].defaultVal
                                       << " adjust = "
                                       << fSystematicsAdjust[i] << endl;
  }
  */

  if(r.Get("LazyContoursUseMinuit",     tmpi)) fContourMethod = tmpi;
  if(r.Get("LazyContoursWant68",        tmpb)) fWant68 = tmpb;
  if(r.Get("LazyContoursWant90",        tmpb)) fWant90 = tmpb;
  if(r.Get("LazyContoursWant99",        tmpb)) fWant99 = tmpb;
  if(r.Get("WantMinimum",               tmpb)) fWantMinimum = tmpb;
  if(r.Get("Want2DProjections",         tmpb)) fWant2DProjections = tmpb;
  if(r.Get("WantGridSearch",            tmpb)) fWantGridSearch = tmpb;
  if(r.Get("WantSlices",                tmpb)) fWantSlices = tmpb;
  if(r.Get("Want2DSlices",              tmpb)) fWant2DSlices = tmpb;
  if(r.Get("WantCovarianceMatrix",      tmpb)) fWantCovariance = tmpb;
  if(r.Get("AllowBestFitOutsideRange",  tmpb)) fAllowBestFitOutsideRange = tmpb;
  if(r.Get("AllowContoursOutsideRange", tmpb)) fAllowContoursOutsideRange = tmpb;
  if(r.Get("TryPreviousMinimum",        tmpb)) fTryPreviousMinimum = tmpb;
  if(r.Get("ForceBestFitToTruth",       tmpb)) fForceBestFitToTruth = tmpb;
  if(r.Get("PrecScale",                 tmpd)) fPrecScale = tmpd;

  if(r.Get("AnalyticNormalization",     tmpb)) fAnalyticNormalization = tmpb;

  fCoordConv.Prepare(r);

  if(r.Get("ContoursList", tmps)){
    vector<int> xs, ys;
    ParseContourList(tmps, xs, ys);
    fXaxes.clear(); fYaxes.clear();
    for(unsigned int n = 0; n < xs.size(); ++n){
      fXaxes.push_back(NCType::EFitParam(xs[n]));
      fYaxes.push_back(NCType::EFitParam(ys[n]));
    }
  }

  if(r.Get("ProjectionsList", tmps)){
    const vector<int> pl = NC::Utility::ParseNumberList(tmps);
    for(unsigned int n = 0; n < pl.size(); ++n)
      fProjectionsList.push_back(NCType::EFitParam(pl[n]));
  }
}


//----------------------------------------------------------------------
void NC::FitMaster::Run(const NC::OscProb::OscPars* trueOscPars)
{
  assert(fXaxes.size() == fYaxes.size());

  assert(fContourMethod == 0 || fContourMethod == 1 || fContourMethod == 2);

  //?!? TODO  fSoftPenalty = true;

  // Create functionoid that redirects to derived-class's implementation.
  GetChiSqrFromDerived chiSqrFunc(fExtrap, fCoordConv.fFitParams, !fAllowBestFitOutsideRange, fAnalyticNormalization);

  TMinuit* minu; // If we use a minuit object for minimization,
                 // keep hold of it in case we want minuit contours.

  // A minimum is required to find contours or to do the projections.
  if((fWantMinimum ||
      fWant68 || fWant90 || fWant99 ||
      !fProjectionsList.empty() ||
      fWantSlices || fWant2DSlices) &&
     !fForceBestFitToTruth){
    MSG("NCFitMaster", Msg::kInfo) << "Begin fit..." << endl;

    TStopwatch fit_sw;

    if(fContourMethod == 0){ // simple
      vector<NCParameter> preciseParams = fCoordConv.fFitParams;
      for(unsigned int n = 0; n < preciseParams.size(); ++n)
        preciseParams[n].SetPrecision(preciseParams[n].Precision()/100);

      MinFinderSimple mf(preciseParams,
                         !fAllowBestFitOutsideRange,
                         fTryPreviousMinimum,
                         &chiSqrFunc);
      fMinChiSqr = mf.FindMin(fMinCoord, true);

      if(fWantCovariance){
        MSG("NCFitMaster", Msg::kInfo)
          << "Finding covariance matrix using finite differences...\n";
        FindCovarianceMatrix(&chiSqrFunc, fMinCoord, preciseParams).Print();
      }
    }
    if(fContourMethod == 1 || fContourMethod == 2){ // hybrid or minuit
      fMinChiSqr = MinFinderMinuit(&chiSqrFunc, fCoordConv.fFitParams, fMinCoord, &minu);
      if(fWantCovariance){
        MSG("NCFitMaster", Msg::kInfo)
          << "Printing Minuit covariance matrix...\n";
        minu->mnmatu(1);
      }
    }

    MSG("NCFitMaster", Msg::kInfo) << endl << "Best fit  - chisq = "
                                   << fMinChiSqr << endl;

    for(unsigned int n = 0; n < fCoordConv.fFitParams.size(); ++n){
      MSG("NCFitMaster", Msg::kInfo) << "  "
                                     << fCoordConv.fFitParams[n].ShortName()
                                     << " = " << fMinCoord[n] << endl;
    }
    MSG("NCFitMaster", Msg::kInfo) << endl;

    // This populates the syst and oscpars translations
    SetBestFitCoordNDim(fMinCoord);
    fExtrap->SetBestFitCoordNDim(fMinCoord); // TODO - make so we don't need this

    MSG("NCFitMaster", Msg::kInfo) << "Fit took ";
    fit_sw.Print("m");
  } // end if need minimum


  if(fForceBestFitToTruth){
    MSG("NCFitMaster", Msg::kInfo) << "Forcing best fit point to be "
                                   << "input truth point" << endl;
    assert(trueOscPars);
    fMinCoord = fCoordConv.CoordNDimFromOscPars(trueOscPars);
    // This populates the syst and oscpars translations
    SetBestFitCoordNDim(fMinCoord);
    fExtrap->SetBestFitCoordNDim(fMinCoord); // TODO - make so we don't need this
  }


  if(fFixedValsForMarg) MargWithFixedVals(&chiSqrFunc);


  // TODO - This doesn't actually work??
  chiSqrFunc.UsePenalty(!fAllowContoursOutsideRange);


  //? fSoftPenalty = false;

  // Want any contours?
  if(fWant68 || fWant90 || fWant99){
    TStopwatch cont_sw;

    // These represent the 3 confidence levels 68%, 90%, 99%
    const double contVals[] = {2.30, 4.61, 9.21};

    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 NCType::EFitParam x = fXaxes[axisPair];
        const NCType::EFitParam y = fYaxes[axisPair];

        const NCParameter xPar = fCoordConv.ParameterForFitParam(x);
        const NCParameter yPar = fCoordConv.ParameterForFitParam(y);

        MSG("NCFitMaster", Msg::kInfo) << "Building contour in "
                                           << xPar.ShortName()
                                           << " and " << yPar.ShortName()
                                           << " (up = " << contVals[conf] << ")..."
                                           << endl;

        if(fContourMethod == 2){ // minuit
          fContourGraphs.push_back(GetContourMinuit(minu, x, y, contVals[conf]));
          delete minu;
        }
        else{
          vector<TGraph> marg_graphs;
          vector<Coord> contour = GetContour(&chiSqrFunc,
                                             fCoordConv.fFitParams,
                                             x, y, contVals[conf],
                                             marg_graphs);
          /* TODO - hacked out because it's broken. FIX!
          assert(marg_graphs.size() == 4*fCoordConv.fFitParams.size());

          // This gets picked up later and made the title of the canvas
          marg_graphs[0].SetName(TString::Format("%s_%s",
                                                 xPar.ShortName().c_str(),
                                                 yPar.ShortName().c_str()));
          marg_graphs[0].SetTitle(xPar.LatexName().c_str());
          marg_graphs[4].SetTitle(yPar.LatexName().c_str());

          int m = 0;
          for(unsigned int n = 8; n < marg_graphs.size(); n+=4, ++m){
            if(m == fCoordConv.FitterIndex(x)) ++m;
            if(m == fCoordConv.FitterIndex(y)) ++m;
            marg_graphs[n].SetTitle(fCoordConv.ParameterForFitterIndex(m).LatexName().c_str());
          }
          fMarginalGraphs.push_back(marg_graphs);
          */

          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("NCFitMaster", Msg::kInfo) << "\nthere are " << contour.size()
                                             << " points on this contour" << endl;
        }
      } // end for axisPair
    } // end for conf

    MSG("NCFitMaster", Msg::kInfo) << "Contours took ";
    cont_sw.Print("m");
  } // end if fWant


  //? fSoftPenalty = true;

  if(fWant2DProjections){
    for(unsigned int axisPair = 0; axisPair < fXaxes.size(); ++axisPair){
      // Look up which axes we should be using.
      const NCType::EFitParam x = fXaxes[axisPair];
      const NCType::EFitParam y = fYaxes[axisPair];

      vector<NCParameter> lowResParams = fCoordConv.fFitParams;
      for(unsigned int n = 0; n < lowResParams.size(); ++n)
        lowResParams[n].SetPrecision(lowResParams[n].Precision()*10);

      const NCParameter xPar = fCoordConv.ParameterForFitParam(x);
      const NCParameter yPar = fCoordConv.ParameterForFitParam(y);

      TH2D* grid = Get2DGrid(&chiSqrFunc, lowResParams, x, y);

      grid->SetName(("grid_"+xPar.ShortName()+"-"+yPar.ShortName()).c_str());
      grid->SetTitle(("Grid;"+xPar.LatexName()+";"+yPar.LatexName()).c_str());
      fLowResGrids.push_back(grid);
    } // end for axisPair
  } // end if fWant2D


  if(!fProjectionsList.empty()){
    TStopwatch proj_sw;

    for(unsigned int nAxis = 0; nAxis < fProjectionsList.size(); ++nAxis){
      const NCType::EFitParam x = fProjectionsList[nAxis];

      const NCParameter xPar = fCoordConv.ParameterForFitParam(x);

      MSG("NCFitMaster", Msg::kInfo) << "Building 1D projection in "
                                     << xPar.ShortName() << "..." << endl;

      double leftErr, rightErr;
      vector<TGraph> marg_graphs;

      // Seemingly for the projections we need more precision in the marginalised variables.
      vector<NCParameter> pars = fCoordConv.fFitParams;
      for(unsigned int n = 0; n < pars.size(); ++n){
        if(n != UInt_t(x)) pars[n].SetPrecision(pars[n].Precision()/10);
      }

      TGraph* proj = Get1DProjection(&chiSqrFunc, pars, x,
                                     leftErr, rightErr, marg_graphs);

      const double bestFit1D = fMinCoord[fCoordConv.FitterIndex(x)];

      MSG("NCFitMaster", Msg::kInfo)
        << "Asymmetric errors in " << xPar.ShortName()
        << " " << leftErr-bestFit1D << " "
        << rightErr-bestFit1D << endl;

      fOneSigmaErrors[xPar.ShortName()] =
        pair<double, double>(bestFit1D-leftErr, rightErr-bestFit1D);

      proj->SetName(("projection_"+xPar.ShortName()).c_str());
      proj->GetHistogram()->SetTitle((";"+xPar.LatexName()+";#chi^{2}").c_str());
      fChiSqGraphs.push_back(proj);

      // This gets picked up later and made the title of the canvas
      marg_graphs[0].SetName(TString::Format("%s", xPar.ShortName().c_str()));
      unsigned int m = 0;
      for(unsigned int n = 0; n < marg_graphs.size(); ++n, ++m){
        if(m == UInt_t(x)) ++m;
        marg_graphs[n].SetTitle(fCoordConv.ParameterForFitterIndex(m).LatexName().c_str());
      }
      fMarginalGraphs1D.push_back(marg_graphs);

    } // end for nAxis

    MSG("NCFitMaster", Msg::kInfo) << "Projections took ";
    proj_sw.Print("m");
  } // end if projections


  if(fWantGridSearch){
    // Arbitrary rescaling to make plots look nice.
    //    if(fCoordConv.IsFit(NCType::kUMu4Sqr)){
      //      fCoordConv.fFitParams[fCoordConv.FitterIndex(NCType::kUMu4Sqr)].Max() = .2;
      //      fCoordConv.fFitParams[fCoordConv.FitterIndex(NCType::kUMu4Sqr)].prec *= .2;
      //      fCoordConv.fFitParams[2].prec = .02; // dmsq
    //    }
    PrintGridSearch(&chiSqrFunc, fCoordConv.fFitParams, cout);
  }

  if(fWantSlices){
    MSG("NCFitMaster", Msg::kInfo) << "Starting slices" << endl;
    fSlices = MakeSlices(&chiSqrFunc, fCoordConv.fFitParams, fMinCoord);
  }

  if(fWant2DSlices) 
    f2DSlices = Make2DSlices(&chiSqrFunc, fCoordConv.fFitParams, fMinCoord);

//  TFile f("~/eval_pos.root", "RECREATE");
//  chiSqrFunc.DrawAndWriteGraph();
//  f.Close();
}


//----------------------------------------------------------------------
void NC::FitMaster::WriteResources()
{
  MSG("NCFitMaster", Msg::kInfo)
    << "NC::FitMaster::WriteResources() - "
    << fExtrap->GetShortName() << endl;

  /* ???
  DrawBestFitSpectra(Detector::kNear);
  DrawBestFitSpectra(Detector::kFar);

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

  gDirectory->mkdir("fit", "fit")->cd();

  GetBestFitPointAsRegistry()->Write();
  if(!fContourGraphs.empty()) DrawAndWriteContourGraphs();
  if(!fChiSqGraphs.empty()) DrawAndWriteChiSqGraphs();
  if(!fLowResGrids.empty()) DrawAndWriteLowResGrids();
  if(!fSlices.empty()) DrawAndWriteSlices();
  if(!f2DSlices.empty()) DrawAndWrite2DSlices();
}


//----------------------------------------------------------------------
void NC::FitMaster::ParseContourList(TString contList,
                                     vector<int>& xs, vector<int>& ys) const
{
  const int p1 = contList.First('(');
  const int p2 = contList.First(',');
  const int p3 = contList.First(')');

  // Assuming failure is signalled by < 0
  // Hopefully this means we reached the end of the string
  if(p1 < 0 || p2 < 0 || p3 < 0) return;

  // Get the bits that should contain the numbers
  // Slicing operator takes start and length
  const TString strX = contList(p1+1, p2-p1-1);
  const TString strY = contList(p2+1, p3-p2-1);

  // If they're not numbers we're well and truly stuck
  assert(strX.IsDigit());
  assert(strY.IsDigit());

  xs.push_back(strX.Atoi());
  ys.push_back(strY.Atoi());

  // Repeat the same operation on the rest of the string
  const TString remain = contList(p3+1, contList.Length()-p3);
  ParseContourList(remain, xs, ys);
}


//----------------------------------------------------------------------
vector<NCType::EFitParam> NC::FitMaster::GetUniqueAxes() const
{
  vector<NCType::EFitParam> ret;

  for(unsigned int n = 0; n < fXaxes.size(); ++n){
    const NCType::EFitParam ax = fXaxes[n];
    bool add = true;
    for(unsigned int m = 0; m < ret.size(); ++m){
      if(ret[m] == ax) add = false;
    } // end for m
    if(add) ret.push_back(ax);
  } // end for n
  for(unsigned int n = 0; n < fYaxes.size(); ++n){
    const NCType::EFitParam ax = fYaxes[n];
    bool add = true;
    for(unsigned int m = 0; m < ret.size(); ++m){
      if(ret[m] == ax) add = false;
    } // end for m
    if(add) ret.push_back(ax);
  } // end for n

  return ret;
}


//----------------------------------------------------------------------
vector<Coord> NC::FitMaster::GetContour(ICallableND* chiSqrFunc,
                                        const vector<NCParameter>& params,
                                        NCType::EFitParam x,
                                        NCType::EFitParam y,
                                        double up,
                                        vector<TGraph>& marg_graphs)
{
  const int xi = fCoordConv.FitterIndex(x);
  const int yi = fCoordConv.FitterIndex(y);

  const NCParameter xp = fCoordConv.ParameterForFitParam(x);
  const NCParameter yp = fCoordConv.ParameterForFitParam(y);

  if(fContourMethod == 0){ // simple
    MarginalizeSimple func(chiSqrFunc, params, fTryPreviousMinimum,
                           xi, yi);
    return FindContour(&func,
                       Coord(fMinCoord[xi], fMinCoord[yi]),
                       fMinChiSqr+up,
                       Coord(xp.Precision(), yp.Precision()),
                       &marg_graphs);
  }

  if(fContourMethod == 1){ // hybrid
    MarginalizeHybrid func(chiSqrFunc, params, fTryPreviousMinimum, x, y);
    return FindContour(&func,
                       Coord(fMinCoord[xi], fMinCoord[yi]),
                       fMinChiSqr+up,
                       Coord(xp.Precision(), yp.Precision()),
                       &marg_graphs);
  }

  assert(0 && "This should never happen");
}


//----------------------------------------------------------------------
TGraph* NC::FitMaster::GetContourMinuit(TMinuit* minu,
                                        NCType::EFitParam x,
                                        NCType::EFitParam y, double up)
{
  const int numContourPoints = 300;
  minu->SetErrorDef(up);
  TGraph* g = (TGraph*)minu->Contour(numContourPoints,
                                     fCoordConv.FitterIndex(x),
                                     fCoordConv.FitterIndex(y));
  if(minu->GetStatus()) MSG("NCExtrapolation", Msg::kWarning) << "Couldn't find contour.\n";
  return g;
}


//----------------------------------------------------------------------
TH2D* NC::FitMaster::Get2DGrid(ICallableND* chiSqrFunc,
                               const vector<NCParameter>& params,
                               NCType::EFitParam x, NCType::EFitParam y)
{
  const int xi = fCoordConv.FitterIndex(x);
  const int yi = fCoordConv.FitterIndex(y);

  const NCParameter xp = fCoordConv.ParameterForFitParam(x);
  const NCParameter yp = fCoordConv.ParameterForFitParam(y);

  if(fContourMethod == 0){ // simple
    MarginalizeSimple func(chiSqrFunc, params, fTryPreviousMinimum, xi, yi);
    return GetLowResGrid(&func,
                         Coord(xp.Min(), yp.Min()),
                         Coord(xp.Max(), yp.Max()),
                         Coord(xp.Precision(), yp.Precision()));
  }

  if(fContourMethod == 1 || fContourMethod == 2){ // hybrid or minuit
    MarginalizeHybrid func(chiSqrFunc, params, fTryPreviousMinimum, xi, yi);
    return GetLowResGrid(&func,
                         Coord(xp.Min(), yp.Min()),
                         Coord(xp.Max(), yp.Max()),
                         Coord(xp.Precision(), yp.Precision()));
  }

  assert(0 && "This should never happen");
}


//----------------------------------------------------------------------
TGraph* NC::FitMaster::Get1DProjection(ICallableND* chiSqrFunc,
                                       const vector<NCParameter>& params,
                                       NCType::EFitParam x,
                                       double& leftErr,
                                       double& rightErr,
                                       vector<TGraph>& marg_graphs)
{
  const int xi = fCoordConv.FitterIndex(x);
  const NCParameter xp = fCoordConv.ParameterForFitParam(x);

  if(fContourMethod == 0){ // simple
    MarginalizeSimple func(chiSqrFunc, params, fTryPreviousMinimum, xi, -1);
    return CreateOneDimProjection(&func, xp, fMinChiSqr, leftErr, rightErr, marg_graphs);
  }

  if(fContourMethod == 1 || fContourMethod == 2){ // hybrid or minuit
    MarginalizeHybrid func(chiSqrFunc, params, fTryPreviousMinimum, xi, -1);
    return CreateOneDimProjection(&func, xp, fMinChiSqr, leftErr, rightErr, marg_graphs);
  }

  assert(0 && "This should never happen");
}


//----------------------------------------------------------------------
Registry* NC::FitMaster::GetBestFitPointAsRegistry() const
{
  Registry* rBestFit = new Registry;
  rBestFit->SetName("bestFit_registry");
  rBestFit->UnLockKeys();
  rBestFit->UnLockValues();

  for(unsigned int n = 0; n < fMinCoord.size(); ++n){
    rBestFit->Set(fCoordConv.ParameterForFitterIndex(n).ShortName().c_str(),
                  fMinCoord[n]);
  }

  rBestFit->Set("best_chisq", fMinChiSqr);

  for(OneSigmaErrors::const_iterator it = fOneSigmaErrors.begin();
      it != fOneSigmaErrors.end(); ++it){
    const string name = it->first;
    const double lo = it->second.first;
    const double hi = it->second.second;
    rBestFit->Set((name+"_one_sigma_up").c_str(), lo);
    rBestFit->Set((name+"_one_sigma_down").c_str(), hi);
  }

  for(ChisqMargWithFixed::const_iterator it = fChisqsMargWithFixed.begin();
      it != fChisqsMargWithFixed.end(); ++it){
    const string name = it->first;
    const double chisq = it->second;
    rBestFit->Set(("chisq_marg_"+name+"_fixed").c_str(), chisq);
  }

  return rBestFit;
}


//----------------------------------------------------------------------
void NC::FitMaster::
SetFixedValuesForMarginalization(const NC::OscProb::OscPars* pars)
{
  fFixedValsForMarg = pars;
}


//----------------------------------------------------------------------
void NC::FitMaster::MargWithFixedVals(const ICallableND* chiSqrFunc)
{
  TStopwatch sw;

  const vector<NCType::EFitParam> axes = GetUniqueAxes();

  for(unsigned int nAxis = 0; nAxis < axes.size(); ++nAxis){
    const NCType::EFitParam x = axes[nAxis];

    const NCParameter xPar = fCoordConv.ParameterForFitParam(x);

    const int fixIndex = fCoordConv.FitterIndex(x);
    CoordNDim fix = fCoordConv.CoordNDimFromOscPars(fFixedValsForMarg);
    vector<int> vary;
    for(unsigned int n = 0; n < fix.size(); ++n){
      if(int(n) != fixIndex) vary.push_back(n);
    }
    FixVars fv(chiSqrFunc, fix, vary);

    vector<NCParameter> pars;
    for(unsigned int n = 0; n < fCoordConv.fFitParams.size(); ++n){
      if(int(n) != fixIndex) pars.push_back(fCoordConv.fFitParams[n]);
    }
    // Copy the increased precision used in projections, for fairness.
    for(unsigned int n = 0; n < pars.size(); ++n){
      pars[n].SetPrecision(pars[n].Precision()/10);
    }

    MinFinderSimple mf(pars, false, false, &fv); // ???

    MSG("NCFitMaster", Msg::kInfo) << "Holding " << xPar.ShortName()
                                   << " fixed at " << fix[fixIndex]
                                   << " while marginalizing over other"
                                   << " parameters..." << endl;

    CoordNDim junk;
    const double chisq = mf.FindMin(junk);

    fChisqsMargWithFixed[xPar.ShortName()] = chisq;
  } // end for nAxis

  MSG("NCFitMaster", Msg::kInfo) << "Marginalizations took ";
  sw.Print("m");
}


//----------------------------------------------------------------------
void NC::FitMaster::DrawAndWriteContourGraphs()
{
  assert(fXaxes.size() == fYaxes.size());
  for(unsigned int n = 0; n < fXaxes.size(); ++n){
    const NCParameter xPar = fCoordConv.ParameterForFitParam(fXaxes[n]);
    const NCParameter yPar = fCoordConv.ParameterForFitParam(fYaxes[n]);

    const TString canvasName = "contours_"+xPar.ShortName()+"_"+yPar.ShortName();

    const double xBest = fMinCoord[fCoordConv.FitterIndex(fXaxes[n])];
    const double yBest = fMinCoord[fCoordConv.FitterIndex(fYaxes[n])];

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

    TCanvas canv(canvasName);
    TH2* forceAxes = fCoordConv.AxesForParameters("", "", fXaxes[n], fYaxes[n]);
    forceAxes->SetDirectory(0); // Don't write out useless axes to disk.
    forceAxes->GetXaxis()->CenterTitle();
    forceAxes->GetYaxis()->CenterTitle();
    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();
  }

  for(unsigned int n = 0; n < fMarginalGraphs.size(); ++n){
    vector<TGraph>& gs = fMarginalGraphs[n];
    if(gs.empty()) continue;
    TCanvas c(TString::Format("marg_vars_%s_canv", gs[0].GetName()).Data());

    double maxx = 0;
    for(unsigned int i = 0; i < gs.size(); ++i){
      double x, y;
      gs[i].GetPoint(gs[i].GetN()-1, x, y);
      if(x > maxx) maxx = x;
    }

    TH2D h("h", "Marginalized variables",
           100, 0, maxx, 100, -2, 4.5);
    h.Draw();
    TLegend l(.5, .6, .85, .85);
    l.SetBorderSize(0);
    for(unsigned int m = 0; m < gs.size(); ++m){
      TGraph& g = gs[m];
      g.SetName(TString::Format("marg_graph_%d_%d", n, m).Data());
      g.SetMarkerStyle(kFullCircle);
      g.SetLineColor(1+m/4);
      g.SetMarkerColor(1+m/4);
      g.Draw("same lp");
      if(m%4==0) l.AddEntry(&g, "", "lp");
      //      g.Write();
    }
    l.Draw();
    c.Write();
  }
}


//----------------------------------------------------------------------
void NC::FitMaster::DrawAndWriteChiSqGraphs()
{
  for(unsigned int n = 0; n < fChiSqGraphs.size(); ++n){
    TGraph* graph = fChiSqGraphs[n];

    TCanvas* canv = new TCanvas(TString::Format("canv_%s", graph->GetName()));

    graph->GetXaxis()->CenterTitle();

    TAxis* yax = graph->GetYaxis();
    yax->CenterTitle();
    graph->Draw("al");
    const double ymax = yax->GetXmax();
    yax->SetRangeUser(0, ymax < 20 ? ymax : 20);

    canv->Write();
    graph->Write();
  }

  for(unsigned int n = 0; n < fMarginalGraphs1D.size(); ++n){
    vector<TGraph>& gs = fMarginalGraphs1D[n];
    if(gs.empty()) continue;
    TCanvas c(TString::Format("marg_vars_%s_canv", gs[0].GetName()).Data());

    double minx, maxx, y;
    gs[0].GetPoint(0, minx, y);
    gs[0].GetPoint(gs[0].GetN()-1, maxx, y);

    TH2D h("h", "Marginalized variables",
           100, minx, maxx, 100, -2, 4.5);
    h.Draw();
    TLegend l(.5, .6, .85, .85);
    l.SetBorderSize(0);
    for(unsigned int m = 0; m < gs.size(); ++m){
      TGraph& g = gs[m];
      g.SetName(TString::Format("marg_graph_%d_%d", n, m).Data());
      g.SetMarkerStyle(kFullCircle);
      g.SetLineColor(1+m);
      g.SetMarkerColor(1+m);
      g.Draw("same lp");
      l.AddEntry(&g, "", "lp");
    }
    l.Draw();
    c.Write();
  }
}


//---------------------------------------------------------------------
void NC::FitMaster::DrawAndWriteLowResGrids()
{
  for(unsigned int n = 0; n < fLowResGrids.size(); ++n){
    TString name = "canv";
    name += fLowResGrids[n]->GetName();
    TCanvas* canv = new TCanvas(name);
    fLowResGrids[n]->Draw("colz");
    canv->Write();
    fLowResGrids[n]->Write();
  }
}


//---------------------------------------------------------------------
void NC::FitMaster::DrawAndWriteSlices()
{
  for(unsigned int n = 0; n < fSlices.size(); ++n){
    TGraph& graph = fSlices[n];

    TCanvas* canv = new TCanvas(TString::Format("canv_%s", graph.GetName()));

    graph.GetXaxis()->CenterTitle();

    TAxis* yax = graph.GetYaxis();
    yax->CenterTitle();
    graph.Draw("al");
    //    const double ymax = yax->GetXmax();
    //    yax->SetRangeUser(0, ymax < 20 ? ymax : 20);
    yax->SetRangeUser(fMinChiSqr, fMinChiSqr+20);

    canv->Write();
    graph.Write();
  }
}

//---------------------------------------------------------------------
void NC::FitMaster::DrawAndWrite2DSlices()
{
  for(unsigned int n = 0; n < f2DSlices.size(); ++n){
    TH2D* h = f2DSlices[n];

    TCanvas* canv = new TCanvas(TString::Format("canv_%s", h->GetName()));

    h->GetXaxis()->CenterTitle();
    h->GetYaxis()->CenterTitle();

    h->Draw("colz");

    canv->Write();
    h->Write();
  }
}

---