NuSystFitter.cxx

Go to the documentation of this file.

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TAxis.h"
#include "TH2.h"

#include "Minuit2/ContoursError.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/FunctionMinimum.h"

#include "MessageService/MsgService.h"
#include "NtupleUtils/NuMMRunPRL.h"
#include "NtupleUtils/NuMMHelperPRL.h"
#include "NtupleUtils/NuMMParameters.h"
#include "NtupleUtils/NuSystFitter.h"
#include "NtupleUtils/NuFCFitter.h"
#include "NtupleUtils/NuUtilities.h"

ClassImp(NuSystFitter)

CVSID("$Id: NuSystFitter.cxx,v 1.47 2010/01/11 10:55:38 nickd Exp $");

//____________________________________________________________________72
NuSystFitter::NuSystFitter()
  : fQuietMode(false),
    fUp(1.0),
    fNumContourPoints(80),
    fsn2PointsForNeutrinoContour(),
    fsn2PointsForNeutrinoContourPhysical(),
    fsn2PointsForNeutrinoContourUnphysical()
{
  fRuns = new vector<NuMMRun*>();

  for (Double_t sn2 = 1.0; sn2 > 0.871; sn2 -= 0.02){
    fsn2PointsForNeutrinoContourPhysical.push_back(sn2);
  }
  fsn2PointsForNeutrinoContourPhysical.push_back(0.87);
  fsn2PointsForNeutrinoContourPhysical.push_back(0.86);
  for (Double_t sn2 = 0.8595; sn2 > 0.7999; sn2 -= 0.0005){
    fsn2PointsForNeutrinoContourPhysical.push_back(sn2);
  }

  for (Double_t sn2 = 1.0; sn2 < 3; sn2 += 0.02){
    fsn2PointsForNeutrinoContourUnphysical.push_back(sn2);
  }

  fsn2PointsForNeutrinoContour = fsn2PointsForNeutrinoContourPhysical;
}

//____________________________________________________________________72
NuSystFitter::~NuSystFitter()
{
}

//____________________________________________________________________72
void NuSystFitter::push_back(NuMMRun* run)
{
  fRuns->push_back(run);
}
//____________________________________________________________________72
void NuSystFitter::Setup(const std::vector<string> helperInputs,
                         const std::vector<NuHistInterpolator*> fdDataInput,
                         const std::vector<NuHistInterpolator*> ndDataInput,
                         const std::string xsecFile,
                         const std::vector<string> /*outFileList*/,
                         const std::string /*combinedOutputFilename*/,
                         Int_t /*model*/)
{
  vector<NuHistInterpolator*>::const_iterator fdDataIt = fdDataInput.begin();
  vector<NuHistInterpolator*>::const_iterator ndDataIt = ndDataInput.begin();
  for (vector<string>::const_iterator helperIt = helperInputs.begin();
       helperIt != helperInputs.end();
       ++helperIt,++ndDataIt,++fdDataIt){
    NuMMHelperPRL* mmHelper = new NuMMHelperPRL(*helperIt,xsecFile);
    NuMMRun* mmRun = new NuMMRunPRL(mmHelper,*ndDataIt,*fdDataIt);
    this->push_back(mmRun);
  }
}

//____________________________________________________________________72
void NuSystFitter::NoOscPrediction()
{

}

//____________________________________________________________________72
const NuMMParameters NuSystFitter::Minimise(const NuMMParameters& mmPars)
{ 
  
  //Set up Minuit parameters
  ROOT::Minuit2::MnUserParameters mnPars = mmPars.MinuitParameters();
  vector<double> param_vec = mnPars.Params();
  //cout << "Fitting with " << param_vec.size() << " parameters." << endl;
  //cout << "Fitting with " << mnPars.VariableParameters() << " variable parameters." << endl;
  if (!fQuietMode) mmPars.PrintStatus();
    
  //Run Minuit fit
  ROOT::Minuit2::MnMigrad mnMigrad(*this,mnPars);
  ROOT::Minuit2::FunctionMinimum funcMin = mnMigrad();
 
  Bool_t didpass = false;
 
  if (!funcMin.IsValid()){
    cout << "Function minimum is not valid!" << endl;
  }
  else{
    didpass = true;
    if (!fQuietMode){
      cout << "Function minimum is: " << endl
      << "Sn2: " << funcMin.UserParameters().Parameter(1).Value()
      << ", dm2: " << funcMin.UserParameters().Parameter(0).Value()
      << "; norm: " << funcMin.UserParameters().Parameter(2).Value()
      << "; shwEn: " << funcMin.UserParameters().Parameter(3).Value()
      << "; NCBack: " << funcMin.UserParameters().Parameter(4).Value()
      << "; Sn2bar: " << funcMin.UserParameters().Parameter(6).Value()
      << "; dm2bar: " << funcMin.UserParameters().Parameter(5).Value()
      << endl;
    }
  }
  
  //cout << "Diagonal covar elements:" << endl;
  //ROOT::Minuit2::MnUserCovariance covar = funcMin.UserCovariance();
  //for (unsigned int i = 0; i < covar.Nrow(); i++) {
  //  cout << "element " << i << " = " << covar(i,i) << endl;
  //}
  
  NuMMParameters bestFitPars(funcMin);
  bestFitPars.MinuitPass(didpass);
    //bestFitPars.PrintStatus();
    return bestFitPars;
  
}

//____________________________________________________________________72
const NuMMParameters NuSystFitter::Minimise(const NuMMParameters& mmPars, bool &success)
{ 

  //Set up Minuit parameters
  ROOT::Minuit2::MnUserParameters mnPars = mmPars.MinuitParameters();
 
  //Run Minuit fit
  ROOT::Minuit2::MnMigrad mnMigrad(*this,mnPars);
  ROOT::Minuit2::FunctionMinimum funcMin = mnMigrad();
  
  if (!funcMin.IsValid()){
    cout << "Function minimum is not valid!" << endl;
  }
  else{
    if (!fQuietMode){
      cout << "Function minimum is: " << endl
           << "Sn2: " << funcMin.UserParameters().Parameter(1).Value()
           << ", dm2: " << funcMin.UserParameters().Parameter(0).Value()
           << "; norm: " << funcMin.UserParameters().Parameter(2).Value()
           << "; shwEn: " << funcMin.UserParameters().Parameter(3).Value()
           << "; NCBack: " << funcMin.UserParameters().Parameter(4).Value()
           << "; Sn2bar: " << funcMin.UserParameters().Parameter(6).Value()
           << "; dm2bar: " << funcMin.UserParameters().Parameter(5).Value()
           << endl;
    }
  }
  
  //cout << "Diagonal covar elements:" << endl;
  //ROOT::Minuit2::MnUserCovariance covar = funcMin.UserCovariance();
  //for (unsigned int i = 0; i < covar.Nrow(); i++) {
  //  cout << "element " << i << " = " << covar(i,i) << endl;
  //}

  success = funcMin.IsValid();
  NuMMParameters bestFitPars(funcMin);
  return bestFitPars;

}

//____________________________________________________________________72
const TGraph* NuSystFitter::PlotContours(const NuMMParameters& mmPars)
{ 

  
  //Set up Minuit parameters
  ROOT::Minuit2::MnUserParameters mnPars = mmPars.MinuitParameters();
  
  //Run Minuit fit
  ROOT::Minuit2::MnMigrad mnMigrad(*this,mnPars);
  ROOT::Minuit2::FunctionMinimum funcMin = mnMigrad();
  
  if (!funcMin.IsValid()){
    cout << "Function minimum is not valid!" << endl;
  }
  else{
    cout << "Contour Minimization: " << endl;
    cout << "Function minimum is: " << endl
         << "Sn2: " << funcMin.UserParameters().Parameter(1).Value()
         << ", dm2: " << funcMin.UserParameters().Parameter(0).Value()
         << "; norm: " << funcMin.UserParameters().Parameter(2).Value()
         << "; shwEn: " << funcMin.UserParameters().Parameter(3).Value()
         << "; NCBack: " << funcMin.UserParameters().Parameter(4).Value()
         << "; Sn2bar: " << funcMin.UserParameters().Parameter(6).Value()
         << "; dm2bar: " << funcMin.UserParameters().Parameter(5).Value()
         << endl;
  }
  
  ROOT::Minuit2::MnContours contours(*this,funcMin);
  ROOT::Minuit2::ContoursError contErr = contours.Contour(mmPars.ContourPars().first,
                                                          mmPars.ContourPars().second,
                                                          fNumContourPoints);
  
  std::vector<std::pair<Double_t,Double_t> > vContPoints = contErr();
  //std::vector<std::pair<Double_t,Double_t> > vContPoints(10);// = contErr();
  
  Double_t* x = new Double_t[vContPoints.size()];
  for (UInt_t pointCount = 0; pointCount<vContPoints.size(); ++pointCount){
    x[pointCount] = vContPoints[pointCount].first;
  }
  
  Double_t* y = new Double_t[vContPoints.size()];
  for (UInt_t pointCount = 0; pointCount<vContPoints.size(); ++pointCount){
    y[pointCount] = vContPoints[pointCount].second;  
  }
  TGraph* gCont = new TGraph(vContPoints.size(),y,x);
//  cout << "drawing contour" << endl;
//  gCont->SetName("gCont");
//  gDirectory->Print();
//  TFile* file=TFile::Open("contourplot.root","RECREATE");
//  gDirectory->Print();
//  gCont->Write("gCont");
//  file->Close();
  return gCont;
}

//____________________________________________________________________72
const TGraph* NuSystFitter
::OneSidedNeutrinoContourFinder(const NuMMParameters& inputPars)
{
  std::vector<std::pair<Double_t,Double_t> > contPoints =
    this->NeutrinoContourFinder(inputPars);

  //Put the contour points in a TGraph
  TGraph* contour = new TGraph(contPoints.size());
  Int_t point=0;
  for (vector<pair<Double_t, Double_t> >::const_iterator it =
         contPoints.begin();
       it != contPoints.end();
       ++it){
    contour->SetPoint(point,(*it).first,(*it).second);
    cout << "Set contour point " << point
         << ", " << (*it).first
         << ", " << (*it).second
         << endl;
    ++point;
  }
  return contour;
}

//____________________________________________________________________72
const TGraph* NuSystFitter
::TwoSidedNeutrinoContourFinder(const NuMMParameters& inputPars)
{
  fsn2PointsForNeutrinoContour = fsn2PointsForNeutrinoContourPhysical;
  std::vector<std::pair<Double_t,Double_t> > contPointsPhys =
    this->NeutrinoContourFinder(inputPars);

  fsn2PointsForNeutrinoContour = fsn2PointsForNeutrinoContourUnphysical;
  std::vector<std::pair<Double_t,Double_t> > contPointsUnphys =
    this->NeutrinoContourFinder(inputPars);

  //Make a single vector of points
  std::vector<pair<Double_t,Double_t> > contourPoints;
  for (vector<pair<Double_t, Double_t> >::const_iterator it =
         contPointsPhys.begin();
       it != contPointsPhys.end();
       ++it){
    pair<Double_t,Double_t> currentPoint((*it).first,(*it).second);
    contourPoints.push_back(currentPoint);
  }
  for (vector<pair<Double_t, Double_t> >::reverse_iterator
         itr = contPointsUnphys.rbegin();
       itr != contPointsUnphys.rend();
       ++itr){
    pair<Double_t,Double_t> currentPoint((*itr).first,(*itr).second);
    contourPoints.push_back(currentPoint);
  }


  //Put the contour points in a TGraph
  TGraph* contour = new TGraph(contourPoints.size());
  Int_t point=0;
  for (vector<pair<Double_t, Double_t> >::const_iterator it =
         contourPoints.begin();
       it != contourPoints.end();
       ++it){
    contour->SetPoint(point,(*it).first,(*it).second);
    cout << "Set contour point " << point
         << ", " << (*it).first
         << ", " << (*it).second
         << endl;
    ++point;
  }
  return contour;
}

//____________________________________________________________________72
const std::vector<std::pair<Double_t,Double_t> > NuSystFitter
::NeutrinoContourFinder(const NuMMParameters& inputPars)
{ 
  //Find the target chi2 for the contour
  NuMMParameters bestFitPars = this->Minimise(inputPars);
  Double_t bestFitChi2 = (*this)(bestFitPars.VectorParameters());
  Double_t targetChi2 = bestFitChi2 + this->Up();
  cout << "Target chi2 = " << targetChi2 << endl;

  //Fix sn2 ready for the dm2 scans:
  NuMMParameters mmPars = inputPars;
  mmPars.FixSn2();
  //Make somewhere to put the contour points
  std::vector<pair<Double_t,Double_t> > pointsUp;
  std::vector<pair<Double_t,Double_t> > pointsDown;

  //Run the contour in quiet mode
  this->QuietModeOn();

  //Get the points for the contour
  for (vector<Double_t>::iterator snIt = fsn2PointsForNeutrinoContour.begin();
       snIt != fsn2PointsForNeutrinoContour.end();
       ++snIt){
    Double_t sn2 = *snIt;
    mmPars.Sn2(sn2);
    const std::pair<Double_t, Double_t> values =
      this->DmScanForContour(mmPars,targetChi2);
    cout << "Contour points: sn2 = " << sn2
         << "; dm2 = " << values.first << " and " << values.second
         << endl;
    if ((values.first < 0.0) || (values.second < 0.0)){
      cout << "End of contour at sn2 = " << mmPars.Sn2() << endl;
      break;
    }
    else{
      pair<Double_t, Double_t> pointDown(sn2,values.first);
      pair<Double_t, Double_t> pointUp(sn2,values.second);
      pointsDown.push_back(pointDown);
      pointsUp.push_back(pointUp);
    }
  }

  //Make a single vector of points
  std::vector<pair<Double_t,Double_t> > contour;
  for (vector<pair<Double_t, Double_t> >::const_iterator it = pointsUp.begin();
       it != pointsUp.end();
       ++it){
    pair<Double_t,Double_t> currentPoint((*it).first,(*it).second);
    contour.push_back(currentPoint);
  }
  for (vector<pair<Double_t, Double_t> >::reverse_iterator
         itr = pointsDown.rbegin();
       itr != pointsDown.rend();
       ++itr){
    pair<Double_t,Double_t> currentPoint((*itr).first,(*itr).second);
    contour.push_back(currentPoint);
  }
  return contour;
}

//____________________________________________________________________72
const std::pair<Double_t,Double_t> NuSystFitter
::DmScanForContour(const NuMMParameters& mmPars,
                   Double_t targetChi2)
{
  std::pair<Double_t,Double_t> errors;

  NuMMParameters bestFitPoint = this->Minimise(mmPars);
  Double_t bestFitChi2 = (*this)(bestFitPoint.VectorParameters());
  Double_t bestFitDm2 = bestFitPoint.Dm2();

  if (bestFitChi2 > targetChi2){
    errors.first = -1.0;
    errors.second = -1.0;
    return errors;
  }

  cout << "Best chi2: " << bestFitChi2
       << " at dm2 = " << bestFitDm2
       << ", sn2 = " << bestFitPoint.Sn2()
       << ", norm = " << bestFitPoint.Normalisation()
       << ", NC = " << bestFitPoint.NCBackgroundScale()
       << ", shwEn = " << bestFitPoint.ShwEnScale()
       << endl
       << ". Target chi2: " << targetChi2
       << endl;

  for (Int_t direction = -1; direction < 2; direction += 2){

    cout << "Direction = " << direction << endl;
    
    NuMMParameters scanMMPars = mmPars;
    scanMMPars.Dm2(bestFitDm2);
    scanMMPars.FixDm2();
    
    Double_t incrementArray[9] = {0.1e-3,
                                  0.03e-3,
                                  0.01e-3,
                                  0.003e-3,
                                  0.001e-3,
                                  0.0003e-3,
                                  0.0001e-3,
                                  0.00003e-3,
                                  0.00001e-3};
    Double_t firstChi2 = bestFitChi2;
    Double_t firstDm2 = bestFitDm2;
    Double_t currentChi2 = -1.0;
    Double_t lastChi2 = -1.0;
    Double_t lastDm2 = -1.0;
    for (Int_t incCount = 0; incCount<9; ++incCount){
      Double_t increment = incrementArray[incCount];
      
      currentChi2 = firstChi2;
      scanMMPars.Dm2(firstDm2);   
      
      cout << "Starting while with firstChi2 = " << firstChi2
           << ", firstDm2 = " << firstDm2
           << ", increment = " << increment
           << endl;
      
      while (currentChi2 < targetChi2){
        lastChi2 = currentChi2;
        lastDm2 = scanMMPars.Dm2();

        scanMMPars.Dm2(scanMMPars.Dm2() + direction*increment);
        cout << "Trying dm2 = " << scanMMPars.Dm2() << "... ";
        flush(cout);
        if (scanMMPars.AreAllParametersFixed()){
          currentChi2 = (*this)(scanMMPars.VectorParameters());
        }
        else{
          currentChi2 = (*this)(this->Minimise(scanMMPars).VectorParameters());
        }
        cout << "chi2 = " << currentChi2 << endl;
      }
      
      firstDm2 = lastDm2;
      firstChi2 = lastChi2;
      
    }

    Double_t chi2Difference = lastChi2 - currentChi2;
    Double_t dm2Difference = lastDm2 - scanMMPars.Dm2();
    Double_t value = -1.0;
    cout << "chi2Difference = " << chi2Difference << endl;
    cout << "dm2Difference = " << dm2Difference << endl;
    if ((fabs(chi2Difference)<1e-12) || (fabs(dm2Difference)<1e-12)){
      cout << "Chi2 or Dm2 difference is zero. Supplying currentChi2."
           << endl;
      value = scanMMPars.Dm2();
    }
    else{
      Double_t fractionBetween = (targetChi2 - currentChi2)/chi2Difference;
      value = fractionBetween*dm2Difference + scanMMPars.Dm2();
    }
    cout << "Error is between " << lastDm2
         << " and " << scanMMPars.Dm2()
         << ", with chi2 between " << lastChi2
         << " and " << currentChi2
         << endl;
    cout << "Value is " << value << endl;
    if (-1 == direction){errors.first = value;}
    if (1 == direction){errors.second = value;}
  }    

  cout << "Final answer is dm2 = " << bestFitDm2
       << " + " << errors.second - bestFitDm2
       << " - " << errors.first - bestFitDm2
       << endl;
  return errors;
}

//____________________________________________________________________72
const std::pair<Double_t,Double_t> NuSystFitter
::SingleParDm2Errors(const NuMMParameters& mmPars)
{
  std::pair<Double_t,Double_t> errors;

  NuMMParameters bestFitPoint = this->Minimise(mmPars);
  Double_t bestFitChi2 = (*this)(bestFitPoint.VectorParameters());
  Double_t targetChi2 = bestFitChi2 + this->Up();
  Double_t bestFitDm2 = bestFitPoint.Dm2();

  cout << "Best chi2: " << bestFitChi2
       << " at dm2 = " << bestFitDm2
       << ". Target chi2: " << targetChi2
       << endl;

  for (Int_t direction = -1; direction < 2; direction += 2){

    cout << "Direction = " << direction << endl;
    
    NuMMParameters scanMMPars = mmPars;
    scanMMPars.Dm2(bestFitDm2);
    scanMMPars.FixDm2();
    
    Double_t incrementArray[5] = {0.1e-3,
                                  0.01e-3,
                                  0.001e-3,
                                  0.0001e-3,
                                  0.00001e-3};
    Double_t firstChi2 = bestFitChi2;
    Double_t firstDm2 = bestFitDm2;
    Double_t currentChi2 = -1.0;
    Double_t lastChi2 = -1.0;
    Double_t lastDm2 = -1.0;
    for (Int_t incCount = 0; incCount<5; ++incCount){
      Double_t increment = incrementArray[incCount];
      
      currentChi2 = firstChi2;
      scanMMPars.Dm2(firstDm2);   
      
      cout << "Starting while with firstChi2 = " << firstChi2
           << ", firstDm2 = " << firstDm2
           << ", increment = " << increment
           << endl;
      
      while (currentChi2 < targetChi2){
        lastChi2 = currentChi2;
        lastDm2 = scanMMPars.Dm2();

        scanMMPars.Dm2(scanMMPars.Dm2() + direction*increment);
        cout << "Trying dm2 = " << scanMMPars.Dm2() << "... ";
        if (scanMMPars.AreAllParametersFixed()){
          currentChi2 = (*this)(scanMMPars.VectorParameters());
        }
        else{
          currentChi2 = (*this)(this->Minimise(scanMMPars).VectorParameters());
        }
        cout << "chi2 = " << currentChi2 << endl;
      }
      
      firstDm2 = lastDm2;
      firstChi2 = lastChi2;
      
    }

    Double_t chi2Difference = lastChi2 - currentChi2;
    Double_t dm2Difference = lastDm2 - scanMMPars.Dm2();
    Double_t value = -1.0;
    cout << "chi2Difference = " << chi2Difference << endl;
    cout << "dm2Difference = " << dm2Difference << endl;
    if ((fabs(chi2Difference)<1e-12) || (fabs(dm2Difference)<1e-12)){
      cout << "Chi2 or Dm2 difference is zero. Supplying currentChi2."
           << endl;
      value = scanMMPars.Dm2();
    }
    else{
      Double_t fractionBetween = (targetChi2 - currentChi2)/chi2Difference;
      value = fractionBetween*dm2Difference + scanMMPars.Dm2();
    }
    cout << "Error is between " << lastDm2
         << " and " << scanMMPars.Dm2()
         << ", with chi2 between " << lastChi2
         << " and " << currentChi2
         << endl;
    cout << "Value is " << value << endl;
    if (-1 == direction){errors.first = value;}
    if (1 == direction){errors.second = value;}
  }    

  cout << "Final answer is dm2 = " << bestFitDm2
       << " + " << errors.second - bestFitDm2
       << " - " << errors.first - bestFitDm2
       << endl;
  return errors;
}

//____________________________________________________________________72
double NuSystFitter::operator()(const std::vector<double>& pars) const
{
  NuMMParameters mmPars(pars);
  //cout << "Called operator()" << endl;
  //mmPars.PrintStatus();
  return Likelihood(mmPars);
}

//____________________________________________________________________72
double NuSystFitter::Likelihood(const NuMMParameters& mmPars) const
{
    Double_t likelihood = 0.0;
    for (std::vector<NuMMRun*>::iterator runIt = fRuns->begin();
         fRuns->end() != runIt;
         ++runIt){
        likelihood += (*runIt)->ComparePredWithData(mmPars);
    }
    if (!fQuietMode){
        cout << "Stats likelihood: " << likelihood;
    }
    likelihood += mmPars.PenaltyTerm();
    if (!fQuietMode){
        cout << " + penalty = " << likelihood << endl;
    }
    
    return likelihood;
}

//____________________________________________________________________72
void NuSystFitter::QuietModeOn()
{
  fQuietMode = true;
  for (std::vector<NuMMRun*>::iterator runIt = fRuns->begin();
       fRuns->end() != runIt;
       ++runIt){
    (*runIt)->QuietModeOn();
  }
}

//____________________________________________________________________72
void NuSystFitter::QuietModeOff()
{
  fQuietMode = true;
  for (std::vector<NuMMRun*>::iterator runIt = fRuns->begin();
       fRuns->end() != runIt;
       ++runIt){
    (*runIt)->QuietModeOff();
  }
}

//_______________________________________________________________________
const TGraph *NuSystFitter::Chi2SliceSin2bar(const NuMMParameters & mmParsC, Int_t numpoints, const Double_t minval, const Double_t maxval)
{
        Double_t increment = (maxval - minval) / (Double_t)numpoints;

        if (!fQuietMode) cout << "Generating Likelyhood surface for sin2bar" << endl;

        // Set the parameters object
        NuMMParameters mmPars(mmParsC);

        Double_t sn2value = 0, like = 0, minimum = 0;

        // Make the TGraph to fill
        TGraph *chisurf = new TGraph(numpoints);
        chisurf->SetNameTitle("Chi2Sin2Bar", "#Delta#chi^{2} Surface for Sin^{2}2#bar#theta");
        chisurf->GetXaxis()->SetTitle("Sin^{2}2#bar#theta");

        // We want to scan over the sin2bar surface and generate a likelyhood plot
        for (Int_t i = 0; i < numpoints; i++)
        {
                sn2value = minval + increment*i;
                mmPars.Sn2Bar(sn2value);

                // Generate the likelyhood....
                like = (*this)(mmPars.VectorParameters());
        
                // Fill the array
                //sin2s[i] = sn2value;
                //likelyhoods[i] = like;
                chisurf->SetPoint(i, sn2value, like);
        
                // Is this the lowest value?    
                if (i == 0)
                        minimum = like;
                else
                        minimum = like < minimum ? like : minimum;
        }       

        // Offset the data points
        for (Int_t i = 0; i < numpoints; i++)
        {
                Double_t x, y;
                chisurf->GetPoint(i, x, y);
                chisurf->SetPoint(i, x, y-minimum);
        }
        
        // Now
        if (!fQuietMode) cout << "Surface Generation Done." << endl;

        // Return the completed TGraph
        return chisurf;
}


//_______________________________________________________________________
const TGraph *NuSystFitter::Chi2SliceDm2bar(const NuMMParameters & mmParsC, const Int_t numpoints, const Double_t minval, const Double_t maxval)
{
        Double_t increment = (maxval - minval) / (Double_t)numpoints;


        if (!fQuietMode) cout << "Generating Likelyhood surface for dm2bar" << endl;

        // Set the parameters object
        NuMMParameters mmPars(mmParsC);

        // Make the TGraph to fill
        TGraph *chisurf = new TGraph(numpoints);
        chisurf->SetNameTitle("Chi2Dm2Bar", "#Delta#chi^{2} Surface for #Delta#barm^{2}");
        chisurf->GetXaxis()->SetTitle("#Delta#barm^{2}");
        
        Double_t dm2value = 0, like = 0, minimum = 0;

        // We want to scan over the sin2bar surface and generate a likelyhood plot
        for (Int_t i = 0; i < numpoints; i++)
        {
                dm2value = minval + increment*i;
                mmPars.Dm2Bar(dm2value);

                // Generate the likelyhood....
                like = (*this)(mmPars.VectorParameters());
        
                // Fill the array
//              dm2s[i] = dm2value;
//              likelyhoods[i] = like;
                
                chisurf->SetPoint(i, dm2value, like);
        
                // Is this the lowest value?    
                if (i == 0)
                        minimum = like;
                else
                        minimum = like < minimum ? like : minimum;
        }       

        // Offset the data points
        for (Int_t i = 0; i < numpoints; i++)
        {
                Double_t x, y;
                chisurf->GetPoint(i, x, y);
                chisurf->SetPoint(i, x, y-minimum);
        }
                
        // Now
        if (!fQuietMode) cout << "Surface Generation Done." << endl;

        // Return the completed TGraph
        return chisurf;
}

const TGraph *NuSystFitter::Chi2ValleyBar(const NuMMParameters &mmParsC, const Int_t points, const Double_t min, const Double_t max )
{
        NuMMParameters mmPars(mmParsC);

        //  mmPars.Dm2Bar(2.5e-3);
         // mmPars.Sn2Bar(1.0);
          mmPars.FixNorm();
          mmPars.FixShwScale();
          mmPars.FixNCBackground();
          mmPars.FixDm2();
          mmPars.FixSn2();
          mmPars.FixDm2Bar(false);
          mmPars.FixSn2Bar();

//      Double_t sin2[points], likelyhood[points], mass[points];
        Double_t sinval=0, likeval=0, minimum=0;
        Double_t interval =  (max - min) / (points -1);;

        if(!fQuietMode) cout << "Generating Valley trace plot..." << endl;

        TGraph *sin2like = new TGraph(points);

        for (Int_t i = 0; i < points; i++)
        {
//        cout << endl << "Sin2 " << i+1 << " / " << points << endl;
          // Set up the parameter structure
          sinval = min + interval*i;
          mmPars.Sn2Bar(sinval);

          // Now, minimize the likelyhood for this value
          NuMMParameters mmFit = this->Minimise(mmPars);
          //NuMMParameters mmFit(mmPars);

    likeval = (*this)(mmFit.VectorParameters());
    sin2like->SetPoint(i, sinval, likeval);
        
          // Calculate Minimum
          if (i == 0)
            minimum = likeval;
          else
            minimum = likeval < minimum ? likeval : minimum;
//              cout << "Minimum is now " << minimum << endl;
        }

        if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;
        // Offset the data points
        for (Int_t i = 0; i < points; i++)
        {
                Double_t x, y;
                sin2like->GetPoint(i, x, y);
//              cout << "From: " << x << ", " << y << " to ";
                sin2like->SetPoint(i, x, y-minimum);
//              cout << x << ", " << y-minimum << endl;
        }

        if (!fQuietMode) {
    cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
  }

        // Make the graphs
        return sin2like;
}

const TGraph *NuSystFitter::Chi2ValleyBarDm2Bar( const NuMMParameters &mmParsC,
                                        const Int_t points,
                                        const Double_t min,
                                        const Double_t max )
{
    NuMMParameters mmPars(mmParsC);

    mmPars.FixNorm();
    mmPars.FixShwScale();
    mmPars.FixNCBackground();
    mmPars.FixDm2();
    mmPars.FixSn2();
    mmPars.FixDm2Bar();
    mmPars.FixSn2Bar(false);
    mmPars.ConstrainPhysical();
    

    //  Double_t sin2[points], likelyhood[points], mass[points];
    Double_t dmval=0, likeval=0, minimum=0;
    Double_t interval =  (max - min) / (points -1);;

    if(!fQuietMode) cout << "Generating Valley trace plot for dm2..." << endl;

    TGraph *dm2like = new TGraph(points);
    dm2like->SetName("valleyplot_mass");
    
    for (Int_t i = 0; i < points; i++)
    {
        if (i % 10 == 0) {
          cout << i << "/" << points << endl;
        }
        // Set up the parameter structure
        dmval = min + interval*i;
        mmPars.Dm2Bar(dmval);

        // Now, minimize the likelyhood for this value
        NuMMParameters mmFit = this->Minimise(mmPars);
        //NuMMParameters mmFit(mmPars);

        likeval = (*this)(mmFit.VectorParameters());
        dm2like->SetPoint(i, dmval, likeval);

        // Calculate Minimum
        if (i == 0)
            minimum = likeval;
        else
            minimum = likeval < minimum ? likeval : minimum;
    }

    if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;

    // Offset the data points
    for (Int_t i = 0; i < points; i++)
    {
        Double_t x, y;
        dm2like->GetPoint(i, x, y);
        dm2like->SetPoint(i, x, y-minimum);
    }

    if (!fQuietMode) {
        cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
    }

    return dm2like;
}


const TGraph *NuSystFitter::Chi2ValleyDm2(  const NuMMParameters &mmParsC,
                                      const Int_t points,
                                      const Double_t min,
                                      const Double_t max)
{
    NuMMParameters mmPars(mmParsC);

    mmPars.FixNorm();
    mmPars.FixShwScale();
    mmPars.FixNCBackground();
    mmPars.FixDm2();
    mmPars.ReleaseSn2();
    mmPars.FixDm2Bar();
    mmPars.FixSn2Bar();
    mmPars.ConstrainPhysical();


    //  Double_t sin2[points], likelyhood[points], mass[points];
    Double_t dmval=0, likeval=0, minimum=0;
    Double_t interval =  (max - min) / (points -1);;

    if(!fQuietMode) cout << "Generating Valley trace plot for dm2..." << endl;

    TGraph *dm2like = new TGraph(points);
    dm2like->SetName("valleyplot_mass");

    for (Int_t i = 0; i < points; i++)
    {
        // Print a progress bar
        NuUtilities::ProgressBar(i, points-1,5);
      
        // Set up the parameter structure
        dmval = min + interval*i;
        mmPars.Dm2(dmval);

        // Now, minimize the likelyhood for this value
        NuMMParameters mmFit = this->Minimise(mmPars);
        //NuMMParameters mmFit(mmPars);

        likeval = (*this)(mmFit.VectorParameters());
        dm2like->SetPoint(i, dmval, likeval);

        // Calculate Minimum
        if (i == 0)
            minimum = likeval;
        else
            minimum = likeval < minimum ? likeval : minimum;
    }

    if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;

    // Offset the data points
    for (Int_t i = 0; i < points; i++)
    {
        Double_t x, y;
        dm2like->GetPoint(i, x, y);
        dm2like->SetPoint(i, x, y-minimum);
    }

    if (!fQuietMode) {
        cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
    }

    return dm2like;
}


const TGraph *NuSystFitter::Chi2ValleySn2(  const NuMMParameters &mmParsC,
                                      const Int_t points,
                                      const Double_t min,
                                      const Double_t max)
{
    NuMMParameters mmPars(mmParsC);

    mmPars.FixNorm();
    mmPars.FixShwScale();
    mmPars.FixNCBackground();
    mmPars.ReleaseDm2();
    mmPars.FixSn2();
    mmPars.FixDm2Bar();
    mmPars.FixSn2Bar();
    mmPars.ConstrainPhysical();


    //  Double_t sin2[points], likelyhood[points], mass[points];
    Double_t snval=0, likeval=0, minimum=0;
    Double_t interval =  (max - min) / (points -1);;

    if(!fQuietMode) cout << "Generating Valley trace plot for sn2..." << endl;

    TGraph *sn2like = new TGraph(points);
    sn2like->SetName("valleyplot_sin");

    for (Int_t i = 0; i < points; i++)
    {
        // Print a progress bar
        NuUtilities::ProgressBar(i, points-1,5);

        // Set up the parameter structure
        snval = min + interval*i;
        mmPars.Sn2(snval);

        // Now, minimize the likelyhood for this value
        NuMMParameters mmFit = this->Minimise(mmPars);
        //NuMMParameters mmFit(mmPars);

        likeval = (*this)(mmFit.VectorParameters());
        sn2like->SetPoint(i, snval, likeval);

        // Calculate Minimum
        if (i == 0)
            minimum = likeval;
        else
            minimum = likeval < minimum ? likeval : minimum;
    }

    if (!fQuietMode) cout << endl << "Adjusting Graph Minimum" << endl;

    // Offset the data points
    for (Int_t i = 0; i < points; i++)
    {
        Double_t x, y;
        sn2like->GetPoint(i, x, y);
        sn2like->SetPoint(i, x, y-minimum);
    }

    if (!fQuietMode) {
        cout << endl << "Valley trace plot generation completed" << endl;
        cout << "   Minimum found was " << minimum << endl;
    }

    return sn2like;
}


TH2D *NuSystFitter::LikelihoodSurfaceBar (const NuMMParameters &mmParsIn,
                                        Int_t sinPoints,
                                        Int_t massPoints,
                                        Double_t sn2BarMin,
                                        Double_t sn2BarMax,
                                        Double_t dm2BarMin,
                                        Double_t dm2BarMax,
                                        Bool_t True_Minimum )
{
    // Calculate the bin spacing. This is -1 because we want to 'add'
    // one bin outside of the specified range, so that the bin contents
    // are representative of the center of the bin.
    Double_t sininterval =  (sn2BarMax - sn2BarMin) / (sinPoints -1);
    Double_t massinterval = (dm2BarMax - dm2BarMin) / (massPoints - 1);

    // Offsets. This is how much to shift the bins min and max by
    Double_t sinoffset = sininterval / 2.0;
    Double_t massoffset = massinterval / 2.0;
    // Create the destination histogram
    TH2D *surface = new TH2D(   "chisurf", "#chi^{2} Likelihood Surface",
                                sinPoints, sn2BarMin-sinoffset, sn2BarMax+sinoffset,
                                massPoints, dm2BarMin-massoffset, dm2BarMax+massoffset );
    // Set it up
    surface->GetXaxis()->SetTitle("sin^{2}2#bar#theta");
    surface->GetYaxis()->SetTitle("#Delta#barm^{2}");
    
    // A variable to store the minimum likelihood encountered
    // Double_t minimum = -1;
   
    // Make a copy of the parameters
    NuMMParameters mmPars(mmParsIn);
    
    // Variables to store the current value. These are declared here
    // so that we can verify we have counted them correctly after the
    // looping ends.
    // Double_t sinVal = 0., dmVal = 0.;
    
    if (!fQuietMode) cout << "Generating Likelihood Surface..." << endl;
    
    Double_t minimum = this->FillLikelihoodSurfaceBar(*surface, mmPars);
    
    // // Now, start looping over the surface
    // for (Int_t sinBin = 1; sinBin <= sinPoints; sinBin++) 
    // {
    //     // Calculate and set this sin (bin 0 starts at -sinoffset but
    //     // we want it to contain the value at 0)
    //     sinVal = sn2BarMin + sininterval*(sinBin-1);
    //     mmPars.Sn2Bar(sinVal);
    //     
    //     for (Int_t massBin = 1; massBin <= massPoints; massBin++)
    //     {
    //         // Calculate and set the dmass value
    //         dmVal = dm2BarMin + massinterval*(massBin-1);
    //         mmPars.Dm2Bar(dmVal);
    //         
    //         // Now, calculate the likelihood using this fitter
    //         Double_t likelihood = (*this)(mmPars.VectorParameters());
    //         
    //         // Minimim calculations
    //         if (minimum < 0 || likelihood < minimum)
    //             minimum = likelihood;
    //         
    //         // Set the bin content!
    //         surface->SetBinContent(sinBin, massBin, likelihood);
    //     }
    // }
    // Print a message to say we have finished
    if (!fQuietMode) {
        cout << "Surface Generation Done. Minimum Likelihood = "
             << minimum << endl;
    }
             
    
    // Get the minimum bin on this histogram
    Int_t minbinX, minbinY, minbinZ;
    surface->GetMinimumBin(minbinX, minbinY, minbinZ);
    // Set the fit parameters to this point as a start point
    mmPars.Sn2Bar(surface->GetXaxis()->GetBinCenter(minbinX));
    mmPars.Dm2Bar(surface->GetYaxis()->GetBinCenter(minbinY));

  if (True_Minimum) {
    // Now minimise at this point, and grab the likelihood
    if (!fQuietMode) cout << "Searching for true minimum" << endl;

    NuMMParameters mmFit = this->Minimise(mmPars);
    Double_t minlike = (*this)(mmFit.VectorParameters());

    // Check if the fit failed, if it did, then don't use it
    if (mmFit.MinuitPass())
    {
        // If this is smaller than this minimum we have found (it should be)
        if (minlike <= minimum)
          minimum = minlike;
        else
          cout << "WARNING: Fit Minimisation found minimum higher than grid search" << endl;
    }
  }
    // Now, we have generated likelihoods for all grid points.
    // Offset each one by the minimum, so that it is a delta chi
    // squared surface.
    for (Int_t sinBin = 1; sinBin <= sinPoints; sinBin++) {
        for (Int_t massBin = 1; massBin <= massPoints; massBin++) {
            Double_t content = surface->GetBinContent(sinBin, massBin);
            surface->SetBinContent(sinBin, massBin, content - minimum);
        }
    }
    
    // Now we have generated and offset the histogram. Return it!
    return surface;
}


Double_t NuSystFitter::FillLikelihoodSurface (TH2D &surface, NuMMParameters mmPars)
{
  Double_t minimum = -1;
  Int_t bincount = surface.GetNbinsX() * surface.GetNbinsY();
  
  for (Int_t x = 1; x <= surface.GetNbinsX(); x++)
  {
      for (Int_t y = 1; y <= surface.GetNbinsY(); y++)
      {
          // Draw a progress bar, but only if we have quiet mode on
          if (fQuietMode) NuUtilities::ProgressBar(x*y, bincount, 5);
          
          mmPars.Sn2(surface.GetXaxis()->GetBinCenter(x));
          mmPars.Dm2(surface.GetYaxis()->GetBinCenter(y));
          Double_t likelihood = (*this)(mmPars.VectorParameters());
          
          surface.SetBinContent(x, y, likelihood);
          
          // Work out if this is the minimum
          if (minimum < 0 || minimum > likelihood) minimum = likelihood;
      }
  }
  // Now, return the minimum!
  return minimum;
}


Double_t NuSystFitter::FillLikelihoodSurfaceBar (TH2D &surface, NuMMParameters mmPars)
{
  Double_t minimum = -1;
  Int_t bincount = surface.GetNbinsX() * surface.GetNbinsY();
  
  for (Int_t x = 1; x <= surface.GetNbinsX(); x++)
  {
      for (Int_t y = 1; y <= surface.GetNbinsY(); y++)
      {
          if (fQuietMode) NuUtilities::ProgressBar(x*y, bincount, 5);
          
          mmPars.Sn2Bar(surface.GetXaxis()->GetBinCenter(x));
          mmPars.Dm2Bar(surface.GetYaxis()->GetBinCenter(y));
          Double_t likelihood = (*this)(mmPars.VectorParameters());
          
          surface.SetBinContent(x, y, likelihood);
          
          // Work out if this is the minimum
          if (minimum < 0 || minimum > likelihood) minimum = likelihood;
      }
  }
  // Now, return the minimum!
  return minimum;
}

Double_t NuSystFitter::FillLikelihoodSurfaceCPT (TH2D &surface, NuMMParameters mmPars)
{
  Double_t minimum = -1;
  for (Int_t x = 1; x <= surface.GetNbinsX(); x++)
  {
      for (Int_t y = 1; y <= surface.GetNbinsY(); y++)
      {
          mmPars.Sn2(surface.GetXaxis()->GetBinCenter(x));
          mmPars.Dm2(surface.GetYaxis()->GetBinCenter(y));
          mmPars.Sn2Bar(surface.GetXaxis()->GetBinCenter(x));
          mmPars.Dm2Bar(surface.GetYaxis()->GetBinCenter(y));
          Double_t likelihood = (*this)(mmPars.VectorParameters());
          
          surface.SetBinContent(x, y, likelihood);
          
          // Work out if this is the minimum
          if (minimum < 0 || minimum > likelihood) minimum = likelihood;
      }
  }
  // Now, return the minimum!
  return minimum;
}


const NuMMParameters NuSystFitter::FCMinimise(const NuMMParameters& mmPars)
{
  // Start with a coarse grid search, to seed minuit
  const NuMMParameters mmStart = CoarseGridSearch(mmPars);    
  
  // Now run the normal minuit procedure with this as the starting point
  NuMMParameters mmFit = Minimise(mmStart);
  
  // If minuit failed, do a finer grid search around the coarse starting point
  // Did minuit succeed?
  // if (mmFit.MinuitPass())
  // {

  // } else {
  if (!mmFit.MinuitPass())
  {
      // Minuit failed to fit. Use a more intensive grid search,
      // starting on the minuit starting point (that we get from
      // a grid search)
      MSG("NuSystFitter", Msg::kInfo) << "Minuit fit failed. Doing Finer Grid search" << endl;
      mmFit = FineGridSearch(mmStart);
  } else {
      // Tell them what minuit found
      MSG("NuSystFitter", Msg::kInfo)
        << "Minuit found: dm2: " << mmFit.Dm2()
        << " sn2: " << mmFit.Sn2()
        << " dm2bar: " << mmFit.Dm2Bar()
        << " sn2bar: " << mmFit.Sn2Bar()
        << '\n';
  }
  return mmFit;
}

NuMMParameters NuSystFitter::CoarseGridSearch(NuMMParameters mmGrid)
{
  // The width and height of the grid
  const Int_t gridWidth = 5;
  const Int_t gridHeight = 15;

  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, 0.2, 1, gridHeight, 2e-3, 30e-3);
  // TH2D likesurf("likesurf", "likesurf", gridWidth, 0.2, 1, 998*10, 2e-3, 1000e-3);
  mmGrid = GridSearch(likesurf, mmGrid);
  // cout << "Found dm2: " << mmGrid.Dm2Bar() << ", sn2: " << mmGrid.Sn2Bar() << endl;

  // Now do a log search
  vector<Double_t> bins;
  for (Double_t i = -1.5; i <= -0.; i += 0.1) {
    bins.push_back(pow(10, i));
  }
  Int_t bincount = bins.size()-1;
  TH2D logsurf("logsurf", "logsurf", gridWidth, 0.2, 1, bincount, &(bins.front()));
  mmGrid = GridSearch(logsurf, mmGrid, &mmGrid);

  // This was code to do a much sparser grid search instead of the log search
  // TH2D sparser("sparser", "sparser", gridWidth, 0.2, 1, 400, 30e-3, 1000e-3);
  //  mmGrid = GridSearch(sparser, mmGrid, &mmGrid);

  // Output the minimum
  MSG("NuSystFitter", Msg::kInfo) << "Coarse grid search found sn2: "
    << mmGrid.Sn2Bar() << ", dm: " << mmGrid.Dm2Bar() << '\n';
  
  return mmGrid;
}

NuMMParameters NuSystFitter::FineGridSearch(NuMMParameters mmStart)
{
  // Calculate the location of the fine grid search
  Double_t gridL, gridR, gridT, gridB;
  gridL = gridR = gridT = gridB = -1.0;
  
  // Go +- 0.2 in sin space
  gridL = mmStart.Sn2Bar() - 0.2;
  gridR = gridL + 0.4;
  // Now limit it to the physical realm
  if (gridL <= mmStart.LowerSn2BarLimit()) gridL = mmStart.LowerSn2BarLimit();
  if (gridR > 1.0) gridR = 1.0;

  // And go +-2e-3 in mass space
  gridB = mmStart.Dm2Bar() - 2e-3;
  gridT = gridB + 4e-3;
  // We only really need to worry about lower limit here
  if (gridB <= mmStart.LowerDm2BarLimit()) gridB = mmStart.LowerSn2BarLimit();
  

  MSG("NuSystFitter", Msg::kInfo) 
    << "Doing fine grid search on sin = [ "
    << gridL << ", " << gridR << " ]\n"
    << "                              dm  = [ "
    << gridB << ", " << gridT << " ]\n";

  // Ensure that we haven't done anything stupid like set low above high
  if (gridL >= gridR || gridT <= gridB) {
//      cout << "Error: Negative range in fine grid search" << endl;
      throw runtime_error("Negative range in fine grid search");
  }
  
  // Make the surface
  
  // The width and height of the grid
  const Int_t gridWidth = 19;
  const Int_t gridHeight = 19;
  
  // Generate the surface
  TH2D likesurf("likesurf", "likesurf", gridWidth, gridL, gridR,
    gridHeight, gridB, gridT);
    
  // And now do a grid search on this
  NuMMParameters mmGrid = GridSearch(likesurf, mmStart);
  
  // Output the minimum
  MSG("NusystFitter", Msg::kInfo) << "Fine grid search found sn2bar: "
    << mmGrid.Sn2Bar() << ", dm2bar: " << mmGrid.Dm2Bar() << '\n';
    
  return mmGrid;
}

NuMMParameters NuSystFitter::GridSearch(TH2D &likesurf, NuMMParameters mmGrid, NuMMParameters *Ref)
{
    Double_t likelihood = FillLikelihoodSurfaceBar(likesurf, mmGrid);
    
    // Find the point on this surface that is minimum
    Int_t minbinX = 0, minbinY = 0, minbinZ = 0;
    likesurf.GetMinimumBin(minbinX, minbinY, minbinZ);
    Double_t minsin = likesurf.GetXaxis()->GetBinCenter(minbinX);
    Double_t minmass = likesurf.GetYaxis()->GetBinCenter(minbinY);

    // Check the zero point, to see if it is lower than anything else on the grid
    mmGrid.Sn2Bar(0);
    mmGrid.Dm2Bar(0);
    Double_t zeropoint = (*this)(mmGrid.VectorParameters());
    if (zeropoint < likelihood) {
        // Then the no-oscillation case is the minimum
        minsin = 0.0;
        minmass = 0.0;
    }
    
    // If we have been given a reference point, check to see if this is lower
    if (Ref) {
      Double_t refpoint = (*this)(Ref->VectorParameters());
      if (refpoint < likelihood) {
        MAXMSG("NuSystFitter", 5, Msg::kDebug) << "Reference point is lower than grid search" << endl;
        minsin = Ref->Sn2Bar();
        minmass = Ref->Dm2Bar();
      }
    }
    
    // Set mmPars to the newly found point, then return it
    mmGrid.Sn2Bar(minsin);
    mmGrid.Dm2Bar(minmass);
    
    return mmGrid;
}


//_______________________________________________________________________
const TGraph *NuSystFitter::Chi2SliceDm2(const NuMMParameters & mmParsC, const Int_t numpoints, const Double_t minval, const Double_t maxval)
{
        Double_t increment = (maxval - minval) / (Double_t)numpoints;


        if (!fQuietMode) cout << "Generating Likelyhood surface for dm2" << endl;

        // Set the parameters object
        NuMMParameters mmPars(mmParsC);

        // Make the TGraph to fill
        TGraph *chisurf = new TGraph(numpoints);
        chisurf->SetNameTitle("Chi2Dm2", "#Delta#chi^{2} Surface for #Deltam^{2}");
        chisurf->GetXaxis()->SetTitle("#Deltam^{2}");
        
        Double_t dm2value = 0, like = 0, minimum = 0;

        // We want to scan over the sin2bar surface and generate a likelyhood plot
        for (Int_t i = 0; i < numpoints; i++)
        {
                dm2value = minval + increment*i;
                mmPars.Dm2(dm2value);

                // Generate the likelyhood....
                like = (*this)(mmPars.VectorParameters());
        
                // Fill the array
//              dm2s[i] = dm2value;
//              likelyhoods[i] = like;
                
                chisurf->SetPoint(i, dm2value, like);
        
                // Is this the lowest value?    
                if (i == 0)
                        minimum = like;
                else
                        minimum = like < minimum ? like : minimum;
        }       

        // Offset the data points
        for (Int_t i = 0; i < numpoints; i++)
        {
                Double_t x, y;
                chisurf->GetPoint(i, x, y);
                chisurf->SetPoint(i, x, y-minimum);
        }
                
        // Now
        if (!fQuietMode) cout << "Surface Generation Done." << endl;

        // Return the completed TGraph
        return chisurf;
}

---