// @(#)root/hist:$Name:  $:$Id: TGraph2D.cxx,v 1.00
// Author: Olivier Couet

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include "Riostream.h"
#include "TROOT.h"
#include "TMath.h"
#include "TGraph2D.h"
#include "TGraphDelaunay.h"
#include "TVirtualPad.h"
#include "TVirtualFitter.h"

ClassImp(TGraph2D)


//______________________________________________________________________________
//
// A Graph2D is a graphics object made of three arrays X, Y and Z with the same
// number of points each.
//
// This class has different constructors:
//
// 1) With an array's dimension and three arrays x, y, and z:
//
//    TGraph2D *g = new TGraph2D(n, x, y, z);
//
//    x, y, z arrays can be doubles, floats, or ints.
//
// 2) With an array's dimension only:
//
//    TGraph2D *g = new TGraph2D(n);
//
//    The internal arrays are then filled with SetPoint. The following line
//    fills the the internal arrays at the position "i" with the values x,y,z.
//
//    g->SetPoint(i, x, y, z);
//
// 3) Without parameters:
//
//    TGraph2D *g = new TGraph2D();
//
//    again SetPoint must be used to fill the internal arrays.
//
// 4) From a file:
//
//    TGraph2D *g = new TGraph2D("graph.dat");
//
//    Arrays are read from the ASCII file "graph.dat" according to a specifies
//    format. The format's default value is "%lg %lg %lg"
//
// Note that in any of these three cases, SetPoint can be used to change a data
// point or add a new one. If the data point index (i) is greater than the
// current size of the internal arrays, they are automatically extended.
//
// Specific drawing options can be used to paint a TGraph2D:
//   "TRI"  : The Delaunay triangles are drawn using filled area.
//            An hidden surface drawing technique is used. The surface is
//            painted with the current fill area color. The edges of each
//            triangles are painted with the current line color.
//   "TRIW" : The Delaunay triangles are drawn as wire frame
//   "TRI1" : The Delaunay triangles are painted with color levels. The edges
//            of each triangles are painted with the current line color.
//   "TRI2" : the Delaunay triangles are painted with color levels.
//   "P"    : Draw a marker at each vertex
//   "P0"   : Draw a circle at each vertex. Each circle background is white.
//
// A TGraph2D can be also drawn with ANY options valid to draw a 2D histogram.
//
// When a TGraph2D is drawn with one of the 2D histogram drawing option,
// a intermediate 2D histogram is filled using the Delaunay triangles
// technique to interpolate the data set.
//
// TGraph2D linearly interpolate a Z value for any (X,Y) point given some
// existing (X,Y,Z) points. The existing (X,Y,Z) points can be randomly
// scattered. The algorithm works by joining the existing points to make
// Delaunay triangles in (X,Y). These are then used to define flat planes
// in (X,Y,Z) over which to interpolate. The interpolated surface thus takes
// the form of tessellating triangles at various angles. Output can take the
// form of a 2D histogram or a vector. The triangles found can be drawn in 3D.
//
// This software cannot be guaranteed to work under all circumstances. They
// were originally written to work with a few hundred points in an XY space
// with similar X and Y ranges.
//
//  The picture below has been generated by the following macro:
//
//{
//   TCanvas *c = new TCanvas("c","Graph2D example",0,0,700,600);
//   Double_t x, y, z, P = 6.;
//   Int_t np = 200;
//   TGraph2D *dt = new TGraph2D();
//   TRandom *r = new TRandom();
//   for (Int_t N=0; N<np; N++) {
//      x = 2*P*(r->Rndm(N))-P;
//      y = 2*P*(r->Rndm(N))-P;
//      z = (sin(x)/x)*(sin(y)/y)+0.2;
//      dt->SetPoint(N,x,y,z);
//   }
//   gStyle->SetPalette(1);
//   dt->Draw("surf1");
//}
//

//
// A more complete example can be find in $ROOTSYS/tutorial/graph2dfit.C. It
// produces the following output:
//
//

//
// Definition of Delaunay triangulation (After B. Delaunay):
// For a set S of points in the Euclidean plane, the unique triangulation DT(S)
// of S such that no point in S is inside the circumcircle of any triangle in
// DT(S). DT(S) is the dual of the Voronoi diagram of S.
// If n is the number of points in S, the Voronoi diagram of S is the partitioning
// of the plane containing S points into n convex polygons such that each polygon
// contains exactly one point and every point in a given polygon is closer to its
// central point than to any other. A Voronoi diagram is sometimes also known as
// a Dirichlet tessellation.
//

//______________________________________________________________________________
 TGraph2D::TGraph2D()
         : TNamed("Graph2D","Graph2D"), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(0)
{
   // Graph2D default constructor

   fSize      = 0;
   fMargin    = 0.;
   fNpx       = 40;
   fNpy       = 40;
   fDirectory = 0;
   fHistogram = 0;
   fMaximum   = -1111;
   fMinimum   = -1111;
   fX         = 0;
   fY         = 0;
   fZ         = 0;
   fZout      = 0;
   fMaxIter   = 100000;
   fPainter   = 0;
   fFunctions = new TList;
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(Int_t n, Int_t *x, Int_t *y, Int_t *z)
         : TNamed("Graph2D","Graph2D"), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(n)
{
   // Graph2D constructor with three vectors of ints as input.

   Build(n);

   // Copy the input vectors into local arrays
   for (Int_t N=0; N<fNpoints; N++) {
      fX[N] = (Double_t)x[N];
      fY[N] = (Double_t)y[N];
      fZ[N] = (Double_t)z[N];
   }
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(Int_t n, Float_t *x, Float_t *y, Float_t *z)
         : TNamed("Graph2D","Graph2D"), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(n)
{
   // Graph2D constructor with three vectors of floats as input.

   Build(n);

   // Copy the input vectors into local arrays
   for (Int_t N=0; N<fNpoints; N++) {
      fX[N] = x[N];
      fY[N] = y[N];
      fZ[N] = z[N];
   }
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(Int_t n, Double_t *x, Double_t *y, Double_t *z)
         : TNamed("Graph2D","Graph2D"), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(n)
{
   // Graph2D constructor with three vectors of doubles as input.

   Build(n);

   // Copy the input vectors into local arrays
   for (Int_t N=0; N<fNpoints; N++) {
      fX[N] = x[N];
      fY[N] = y[N];
      fZ[N] = z[N];
   }
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(TH2 *h2)
         : TNamed("Graph2D","Graph2D"), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(0)
{
   // Graph2D constructor with a TH2 (h2) as input.
   // Only the h2's bins within the X and Y axis ranges are used.

   Build(10);

   SetName(h2->GetName());
   SetTitle(h2->GetTitle());

   TAxis *xaxis = h2->GetXaxis();
   TAxis *yaxis = h2->GetYaxis();
   Int_t Xfirst = xaxis->GetFirst();
   Int_t Xlast  = xaxis->GetLast();
   Int_t Yfirst = yaxis->GetFirst();
   Int_t Ylast  = yaxis->GetLast();

   Double_t hmin = h2->GetMinimum();
   Double_t hmax = h2->GetMaximum();

   Double_t x, y, z;
   Int_t k=0;

   for (Int_t i=Xfirst; i<= Xlast; i++) {
      for (Int_t j=Yfirst; j<= Ylast; j++) {
         x = xaxis->GetBinCenter(i);
         y = yaxis->GetBinCenter(j);
         z = h2->GetBinContent(i,j);
         if (z>hmin && z<hmax) {
            SetPoint(k, x, y, z);
            k++;
         }
      }
   }
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(const char *name,const char *title,
                   Int_t n, Double_t *x, Double_t *y, Double_t *z)
         : TNamed(name,title), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(n)
{
   // Graph2D constructor with name, title and three vectors of doubles as input.
   // name   : name of 2D graph (avoid blanks)
   // title  : 2D graph title
   //          if title is of the form "stringt;stringx;stringy;stringz"
   //          the 2D graph title is set to stringt, the x axis title to stringy,
   //          the y axis title to stringy,etc

   Build(n);

   // Copy the input vectors into local arrays
   for (Int_t N=0; N<fNpoints; N++) {
      fX[N] = x[N];
      fY[N] = y[N];
      fZ[N] = z[N];
   }
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(Int_t n)
         : TNamed("Graph2D","Graph2D"), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(0)
{
   // Graph2D constructor. The arrays fX, fY and fZ should be filled via
   // calls to SetPoint

   Build(n);
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(const char *filename, const char *format, Option_t *)
         : TNamed("Graph2D",filename), TAttLine(1,1,1), TAttFill(0,1001),
           TAttMarker(), fNpoints(0)
{
   // Graph2D constructor reading input from filename
   // filename is assumed to contain at least three columns of numbers

   Build(100);

   Double_t x,y,z;
   FILE *fp = fopen(filename,"r");
   if (!fp) {
      MakeZombie();
      Error("TGraph2D", "Cannot open file: %s, TGraph2D is Zombie",filename);
      return;
   }
   char line[80];
   Int_t np = 0;
   while (fgets(line,80,fp)) {
      sscanf(&line[0],format,&x, &y, &z);
      SetPoint(np,x,y,z);
      np++;
   }

   fclose(fp);
}


//______________________________________________________________________________
 TGraph2D::TGraph2D(const TGraph2D &g)
         : TNamed(g), TAttLine(g), TAttFill(g), TAttMarker(g)
{
   // Graph2D copy constructor.

   fNpoints = g.fNpoints;
   Build(fNpoints);

   for (Int_t N=0; N<fNpoints; N++) {
      fX[N] = g.fX[N];
      fY[N] = g.fY[N];
      fZ[N] = g.fZ[N];
   }
}


//______________________________________________________________________________
 TGraph2D::~TGraph2D()
{
   // TGraph2D destructor.

   if (fX)          delete [] fX;
   if (fY)          delete [] fY;
   if (fZ)          delete [] fZ;
   if (fHistogram)  delete fHistogram;
   if (fFunctions) {
      fFunctions->SetBit(kInvalidObject);
      fFunctions->Delete();
      delete fFunctions;
   }
   if (fDirectory) {
      if (!fDirectory->TestBit(TDirectory::kCloseDirectory))
      fDirectory->GetList()->Remove(this);
   }
   if (fPainter)  delete fPainter;
   fX         = 0;
   fY         = 0;
   fZ         = 0;
   fHistogram = 0;
   fDirectory = 0;
   fFunctions = 0;
}


//______________________________________________________________________________
TGraph2D TGraph2D::operator=(const TGraph2D &g)
{
   // Graph2D operator "="

   if (this == &g) return *this;

   fNpoints = g.fNpoints;
   Build(fNpoints);

   for (Int_t N=0; N<fNpoints; N++) {
      fX[N] = g.fX[N];
      fY[N] = g.fY[N];
      fZ[N] = g.fZ[N];
   }
   return g;
}


//______________________________________________________________________________
 void TGraph2D::Build(Int_t n)
{
   // Creates the 2D graph basic data structure

   if (n <= 0) {
      Error("TGraph2D", "Invalid number of points (%d)", n);
      return;
   }

   fSize      = n;
   fMargin    = 0.;
   fNpx       = 40;
   fNpy       = 40;
   fDirectory = 0;
   fHistogram = 0;
   fMaximum   = -1111;
   fMinimum   = -1111;
   fX         = new Double_t[fSize];
   fY         = new Double_t[fSize];
   fZ         = new Double_t[fSize];
   fZout      = 0;
   fMaxIter   = 100000;
   fFunctions = new TList;
   fPainter   = 0;

   Bool_t add = TH1::AddDirectoryStatus();
   if (add && gDirectory) {
      TObject *old = (TObject*)gDirectory->GetList()->FindObject(GetName());
      if (old) {
         Warning("Build","Replacing existing 2D graph: %s (Potential memory leak).",GetName());
         gDirectory->GetList()->Remove(old);
      }
      gDirectory->Append(this);
      fDirectory = gDirectory;
   }
}


//______________________________________________________________________________
 Int_t TGraph2D::DistancetoPrimitive(Int_t px, Int_t py)
{
   // Computes distance from point px,py to a graph

   Int_t distance = 9999;
   if (fHistogram) distance = fHistogram->DistancetoPrimitive(px,py);
   return distance;
}


//______________________________________________________________________________
 void TGraph2D::Draw(Option_t *option)
{
// Specific drawing options can be used to paint a TGraph2D:
//   "TRI"  : The Delaunay triangles are drawn using filled area.
//            An hidden surface drawing technique is used. The surface is
//            painted with the current fill area color. The edges of each
//            triangles are painted with the current line color.
//   "TRIW" : The Delaunay triangles are drawn as wire frame
//   "TRI1" : The Delaunay triangles are painted with color levels. The edges
//            of each triangles are painted with the current line color.
//   "TRI2" : the Delaunay triangles are painted with color levels.
//   "P"    : Draw a marker at each vertex
//   "P0"   : Draw a circle at each vertex. Each circle background is white.
//
// A TGraph2D can be also drawn with ANY options valid to draw a 2D histogram.
//
// When a TGraph2D is drawn with one of the 2D histogram drawing option,
// a intermediate 2D histogram is filled using the Delaunay triangles
// technique to interpolate the data set.

   TString opt = option;
   opt.ToLower();
   if (gPad) {
      if (!gPad->IsEditable()) (gROOT->GetMakeDefCanvas())();
      if (!opt.Contains("same")) {
         //the following statement is necessary in case one attempts to draw
         //a temporary histogram already in the current pad
         if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
         gPad->Clear();
      }
   }
   AppendPad(opt.Data());
}


//______________________________________________________________________________
 void TGraph2D::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
   // Executes action corresponding to one event

   if (fHistogram) fHistogram->ExecuteEvent(event, px, py);
}


//______________________________________________________________________________
 TObject *TGraph2D::FindObject(const char *name) const
{
   // search object named name in the list of functions

   if (fFunctions) return fFunctions->FindObject(name);
   return 0;
}


//______________________________________________________________________________
 TObject *TGraph2D::FindObject(const TObject *obj) const
{
   // search object obj in the list of functions

   if (fFunctions) return fFunctions->FindObject(obj);
   return 0;
}


//______________________________________________________________________________
 Int_t TGraph2D::Fit(const char *fname, Option_t *option, Option_t *)
{
   // Fits this graph with function with name fname

   TF2 *f2 = (TF2*)gROOT->GetFunction(fname);
   if (!f2) {
      Error("Fit","Unknown function: %s",fname);
      return -1;
   }
   return Fit(f2,option,"");
}


//______________________________________________________________________________
 Int_t TGraph2D::Fit(TF2 *f2, Option_t *option, Option_t *)
{
   // Fits this 2D graph with function f2
   //
   //  f2 is an already predefined function created by TF2.
   //  Predefined functions such as gaus, expo and poln are automatically
   //  created by ROOT.
   //
   //  The list of fit options is given in parameter option.
   //     option = "W"  Set all errors to 1
   //            = "U" Use a User specified fitting algorithm (via SetFCN)
   //            = "Q" Quiet mode (minimum printing)
   //            = "V" Verbose mode (default is between Q and V)
   //            = "R" Use the Range specified in the function range
   //            = "N" Do not store the graphics function, do not draw
   //            = "0" Do not plot the result of the fit. By default the fitted function
   //                  is drawn unless the option"N" above is specified.
   //            = "+" Add this new fitted function to the list of fitted functions
   //                  (by default, any previous function is deleted)
   //
   //  In order to use the Range option, one must first create a function
   //  with the expression to be fitted. For example, if your graph2d
   //  has a defined range between -4 and 4 and you want to fit a gaussian
   //  only in the interval 1 to 3, you can do:
   //       TF2 *f2 = new TF2("f2","gaus",1,3);
   //       graph2d->Fit("f2","R");
   //
   //
   //  Setting initial conditions
   //  ==========================
   //  Parameters must be initialized before invoking the Fit function.
   //  The setting of the parameter initial values is automatic for the
   //  predefined functions : poln, expo, gaus. One can however disable
   //  this automatic computation by specifying the option "B".
   //  You can specify boundary limits for some or all parameters via
   //       f2->SetParLimits(p_number, parmin, parmax);
   //  if parmin>=parmax, the parameter is fixed
   //  Note that you are not forced to fix the limits for all parameters.
   //  For example, if you fit a function with 6 parameters, you can do:
   //    func->SetParameters(0,3.1,1.e-6,0.1,-8,100);
   //    func->SetParLimits(4,-10,-4);
   //    func->SetParLimits(5, 1,1);
   //  With this setup, parameters 0->3 can vary freely
   //  Parameter 4 has boundaries [-10,-4] with initial value -8
   //  Parameter 5 is fixed to 100.
   //
   //  Fit range
   //  =========
   //  The fit range can be specified in two ways:
   //    - specify rxmax > rxmin (default is rxmin=rxmax=0)
   //    - specify the option "R". In this case, the function will be taken
   //      instead of the full graph range.
   //
   //  Changing the fitting function
   //  =============================
   //  By default the fitting function Graph2DFitChisquare is used.
   //  To specify a User defined fitting function, specify option "U" and
   //  call the following functions:
   //    TVirtualFitter::Fitter(mygraph)->SetFCN(MyFittingFunction)
   //  where MyFittingFunction is of type:
   //  extern void MyFittingFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
   //
   //  Associated functions
   //  ====================
   //  One or more object (typically a TF2*) can be added to the list
   //  of functions (fFunctions) associated to each graph.
   //  When TGraph::Fit is invoked, the fitted function is added to this list.
   //  Given a graph gr, one can retrieve an associated function
   //  with:  TF2 *myfunc = gr->GetFunction("myfunc");
   //
   //  Access to the fit results
   //  =========================
   //  If the graph is made persistent, the list of
   //  associated functions is also persistent. Given a pointer (see above)
   //  to an associated function myfunc, one can retrieve the function/fit
   //  parameters with calls such as:
   //    Double_t chi2 = myfunc->GetChisquare();
   //    Double_t par0 = myfunc->GetParameter(0); //value of 1st parameter
   //    Double_t err0 = myfunc->GetParError(0);  //error on first parameter
   //
   //  Fit Statistics
   //  ==============
   //  You can change the statistics box to display the fit parameters with
   //  the TStyle::SetOptFit(mode) method. This mode has four digits.
   //  mode = pcev  (default = 0111)
   //    v = 1;  print name/values of parameters
   //    e = 1;  print errors (if e=1, v must be 1)
   //    c = 1;  print Chisquare/Number of degress of freedom
   //    p = 1;  print Probability
   //
   //  For example: gStyle->SetOptFit(1011);
   //  prints the fit probability, parameter names/values, and errors.
   //  You can change the position of the statistics box with these lines
   //  (where g is a pointer to the TGraph):
   //
   //  Root > TPaveStats *st = (TPaveStats*)g->GetListOfFunctions()->FindObject("stats")
   //  Root > st->SetX1NDC(newx1); //new x start position
   //  Root > st->SetX2NDC(newx2); //new x end position

   Int_t fitResult = 0;
   Double_t xmin=0, xmax=0;
   Int_t i, npar,nvpar,nparx;
   Double_t par, we, al, bl;
   Double_t eplus,eminus,eparab,globcc,amin,edm,errdef,werr;
   TF2 *fnew2;

   // Check validity of function
   if (!f2) {
      Error("Fit", "function may not be null pointer");
      return 0;
   }
   if (f2->IsZombie()) {
      Error("Fit", "function is zombie");
      return 0;
   }
   npar = f2->GetNpar();
   if (npar <= 0) {
      Error("Fit", "function %s has illegal number of parameters = %d", f2->GetName(), npar);
      return 0;
   }

   // Check that function has same dimension as graph
   if (f2->GetNdim() != 2) {
      Error("Fit","function %s is not 2-D", f2->GetName());
      return 0;
   }

   Double_t *arglist = new Double_t[100];

   // Decode string choptin and fill fitOption structure
   Foption_t fitOption;
   fitOption.Quiet   = 0;
   fitOption.Verbose = 0;
   fitOption.Bound   = 0;
   fitOption.Like    = 0;
   fitOption.W1      = 0;
   fitOption.Errors  = 0;
   fitOption.Range   = 0;
   fitOption.Gradient= 0;
   fitOption.Nograph = 0;
   fitOption.Nostore = 0;
   fitOption.Plus    = 0;
   fitOption.User    = 0;

   TString opt = option;
   opt.ToUpper();

   if (opt.Contains("U")) fitOption.User    = 1;
   if (opt.Contains("Q")) fitOption.Quiet   = 1;
   if (opt.Contains("V")){fitOption.Verbose = 1; fitOption.Quiet   = 0;}
   if (opt.Contains("W")) fitOption.W1      = 1;
   if (opt.Contains("E")) fitOption.Errors  = 1;
   if (opt.Contains("R")) fitOption.Range   = 1;
   if (opt.Contains("N")) fitOption.Nostore = 1;
   if (opt.Contains("0")) fitOption.Nograph = 1;
   if (opt.Contains("+")) fitOption.Plus    = 1;
   if (opt.Contains("B")) fitOption.Bound   = 1;

///xmin    = fX[0];
///xmax    = fX[fNpoints-1];
///ymin    = fY[0];
///ymax    = fY[fNpoints-1];
///Double_t err0 = GetErrorX(0);
///Double_t errn = GetErrorX(fNpoints-1);
///if (err0 > 0) xmin -= 2*err0;
///if (errn > 0) xmax += 2*errn;
///for (i=0;i<fNpoints;i++) {
///   if (fX[i] < xmin) xmin = fX[i];
///   if (fX[i] > xmax) xmax = fX[i];
///   if (fY[i] < ymin) ymin = fY[i];
///   if (fY[i] > ymax) ymax = fY[i];
///}
///if (rxmax > rxmin) {
///   xmin = rxmin;
///   xmax = rxmax;
///}

   // Check if Minuit is initialized and create special functions
   TVirtualFitter *grFitter = TVirtualFitter::Fitter(this, f2->GetNpar());
   grFitter->Clear();

   // Get pointer to the function by searching in the list of functions in ROOT
   grFitter->SetUserFunc(f2);
   grFitter->SetFitOption(fitOption);

//*-*- Is a Fit range specified?
   Int_t gxfirst, gxlast;
///if (fitOption.Range) {
///   f2->GetRange(xmin, ymin, xmax, ymax);
///   gxfirst = fNpoints +1;
///   gxlast  = -1;
///   for (i=0;i<fNpoints;i++) {
///      if (fX[i] >= xmin && gxfirst > i) gxfirst = i;
///      if (fX[i] <= xmax  && gxlast < i) gxlast  = i;
///   }
///} else {
///   f2->SetRange(xmin, ymin, xmax, ymax);
      gxfirst = 0;
      gxlast  = fNpoints-1;
///}

   // Some initialisations
   if (!fitOption.Verbose) {
      arglist[0] = -1;
      grFitter->ExecuteCommand("SET PRINT", arglist,1);
      arglist[0] = 0;
      grFitter->ExecuteCommand("SET NOW",   arglist,0);
   }

   // Set error criterion for chisquare
   arglist[0] = 1;
   if (!fitOption.User) grFitter->SetFitMethod("Graph2DFitChisquare");
   fitResult = grFitter->ExecuteCommand("SET ERR",arglist,1);
   if (fitResult != 0) {
      // Abnormal termination, MIGRAD might not have converged on a minimum.
      if (!fitOption.Quiet) {
         Warning("Fit","Abnormal termination of minimization.");
      }
      delete [] arglist;
      return fitResult;
   }

   // Transfer names and initial values of parameters to Minuit
   Int_t nfixed = 0;
   for (i=0;i<npar;i++) {
      par = f2->GetParameter(i);
      f2->GetParLimits(i,al,bl);
      if (al*bl != 0 && al >= bl) {
         al = bl = 0;
         arglist[nfixed] = i+1;
         nfixed++;
      }
      we  = 0.3*TMath::Abs(par);
      if (we <= TMath::Abs(par)*1e-6) we = 1;
      grFitter->SetParameter(i,f2->GetParName(i),par,we,al,bl);
   }
   if(nfixed > 0)grFitter->ExecuteCommand("FIX",arglist,nfixed); // Otto

   // Reset Print level
   if (fitOption.Verbose) {
      arglist[0] = 0; grFitter->ExecuteCommand("SET PRINT", arglist,1);
   }

   // Compute sum of squares of errors in the bin range
   Bool_t hasErrors = kFALSE;
   Double_t ex, ey, ez, sumw2=0;
   for (i=gxfirst;i<=gxlast;i++) {
      ex = GetErrorX(i);
      ey = GetErrorY(i);
      ez = GetErrorZ(i);
      if (ex > 0 || ey > 0 || ez > 0) hasErrors = kTRUE;
      sumw2 += ez*ez;
   }
//*-*- Perform minimization
   if (!InheritsFrom("TGraph2DErrors")) SetBit(kFitInit);
   arglist[0] = TVirtualFitter::GetMaxIterations();
   arglist[1] = sumw2*TVirtualFitter::GetPrecision();
   grFitter->ExecuteCommand("MIGRAD",arglist,2);
   if (fitOption.Errors) {
      grFitter->ExecuteCommand("HESSE",arglist,0);
      grFitter->ExecuteCommand("MINOS",arglist,0);
   }

   grFitter->GetStats(amin,edm,errdef,nvpar,nparx);
   f2->SetChisquare(amin);
   Int_t ndf = f2->GetNumberFitPoints()-npar+nfixed;
   f2->SetNDF(ndf);

   // Get return status
   char parName[50];
   for (i=0;i<npar;i++) {
      grFitter->GetParameter(i,parName, par,we,al,bl);
      if (!fitOption.Errors) werr = we;
      else {
         grFitter->GetErrors(i,eplus,eminus,eparab,globcc);
         if (eplus > 0 && eminus < 0) werr = 0.5*(eplus-eminus);
         else                         werr = we;
      }
      if (!hasErrors && ndf > 1) werr *= TMath::Sqrt(amin/(ndf-1));
      f2->SetParameter(i,par);
      f2->SetParError(i,werr);
   }

   // Print final values of parameters.
   if (!fitOption.Quiet) {
      if (fitOption.Errors) grFitter->PrintResults(4,amin);
      else                  grFitter->PrintResults(3,amin);
   }
   delete [] arglist;

   // Store fitted function in histogram functions list and draw
   if (!fitOption.Nostore) {
      if (!fFunctions) fFunctions = new TList;
      if (!fitOption.Plus) {
         TIter next(fFunctions, kIterBackward);
         TObject *obj;
         while ((obj = next())) {
            if (obj->InheritsFrom(TF1::Class())) delete obj;
         }
      }
      fnew2 = new TF2();
      f2->Copy(*fnew2);
      fFunctions->Add(fnew2);
      fnew2->SetParent(this);
      fnew2->Save(xmin,xmax,0,0,0,0);
      if (fitOption.Nograph) fnew2->SetBit(TF1::kNotDraw);
      fnew2->SetBit(TFormula::kNotGlobal);

      if (TestBit(kCanDelete)) return fitResult;
      if (gPad) gPad->Modified();
   }
   return fitResult;
}

//______________________________________________________________________________
 TAxis *TGraph2D::GetXaxis() const
{
   // Get x axis of the graph.

   //if (!gPad) return 0;
   TH1 *h = ((TGraph2D*)this)->GetHistogram();
   if (!h) return 0;
   return h->GetXaxis();
}

//______________________________________________________________________________
 TAxis *TGraph2D::GetYaxis() const
{
   // Get y axis of the graph.

   //if (!gPad) return 0;
   TH1 *h = ((TGraph2D*)this)->GetHistogram();
   if (!h) return 0;
   return h->GetYaxis();
}


//______________________________________________________________________________
 TAxis *TGraph2D::GetZaxis() const
{
   // Get z axis of the graph.

   //if (!gPad) return 0;
   TH1 *h = ((TGraph2D*)this)->GetHistogram();
   if (!h) return 0;
   return h->GetZaxis();
}


//______________________________________________________________________________
 TList *TGraph2D::GetContourList(Double_t contour)
{
   // Returns the X and Y graphs building a contour. A contour level may
   // consist in several parts not connected to each other. This function
   // returns them in a graphs' list.

   if (fNpoints <= 0) {
      Error("GetContourList", "Empty TGraph2D");
      return 0;
   }

   if(!fHistogram) GetHistogram("empty");

   if (!fPainter) fPainter = fHistogram->GetPainter();

   return fPainter->GetContourList(contour);
}


//______________________________________________________________________________
 Double_t TGraph2D::GetErrorX(Int_t) const
{
   // This function is called by Graph2DFitChisquare.
   // It always returns a negative value. Real implementation in TGraph2DErrors

   return -1;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetErrorY(Int_t) const
{
   // This function is called by Graph2DFitChisquare.
   // It always returns a negative value. Real implementation in TGraph2DErrors

   return -1;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetErrorZ(Int_t) const
{
   // This function is called by Graph2DFitChisquare.
   // It always returns a negative value. Real implementation in TGraph2DErrors

   return -1;
}


//______________________________________________________________________________
 TH2D *TGraph2D::GetHistogram(Option_t *option)
{
   // By default returns a pointer to the Delaunay histogram. If fHistogram
   // doesn't exist, books the 2D histogram fHistogram with a margin around
   // the hull. Calls TGraphDelaunay::Interpolate at each bin centre to build up
   // an interpolated 2D histogram.
   // If the "empty" option is selected, returns an empty histogram booked with
   // the limits of fX, fY and fZ. This option is used when the data set is
   // drawn with markers only. In that particular case there is no need to
   // find the Delaunay triangles.

   if (fNpoints <= 0) {
      Error("GetHistogram", "Empty TGraph2D");
      return 0;
   }

   TString opt = option;
   opt.ToLower();
   Bool_t empty = opt.Contains("empty");

   if (fHistogram) {
      if (!empty && fHistogram->GetEntries()==0) {
         delete fHistogram; fHistogram = 0;
      } else {
         return fHistogram;
      }
   }

   Double_t xmax  = GetXmax();
   Double_t ymax  = GetYmax();
   Double_t xmin  = GetXmin();
   Double_t ymin  = GetYmin();
   Double_t hxmax = xmax+fMargin*(xmax-xmin);
   Double_t hymax = ymax+fMargin*(ymax-ymin);
   Double_t hxmin = xmin-fMargin*(xmax-xmin);
   Double_t hymin = ymin-fMargin*(ymax-ymin);

   // Book fHistogram. It is not added in the current directory
   Bool_t add = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
   fHistogram = new TH2D(GetName(),GetTitle(),
                         fNpx ,hxmin, hxmax,
                         fNpy, hymin, hymax);
   TH1::AddDirectory(add);
   fHistogram->SetBit(TH1::kNoStats);

   // Add a TGraphDelaunay in the list of the fHistogram's functions
   TGraphDelaunay *dt = new TGraphDelaunay(this);
   dt->SetMaxIter(fMaxIter);
   dt->SetMarginBinsContent(fZout);
   TList *hl = fHistogram->GetListOfFunctions();
   hl->Add(dt);

   // Option "empty" is selected. An empty histogram is returned.
   if (empty) {
      if (fMinimum != -1111) {
         fHistogram->SetMinimum(fMinimum);
      } else {
         fHistogram->SetMinimum(GetZmin());
      }
      if (fMaximum != -1111) {
         fHistogram->SetMaximum(fMaximum);
      } else {
         fHistogram->SetMaximum(GetZmax());
      }
      return fHistogram;
   }

   Double_t dx = (hxmax-hxmin)/fNpx;
   Double_t dy = (hymax-hymin)/fNpy;

   Double_t x, y, z;
   for (Int_t ix=1; ix<=fNpx; ix++) {
      x  = hxmin+(ix-0.5)*dx;
      for (Int_t iy=1; iy<=fNpy; iy++) {
         y  = hymin+(iy-0.5)*dy;
         z  = dt->ComputeZ(x, y);
         fHistogram->Fill(x, y, z);
      }
   }

   if (fMinimum != -1111) fHistogram->SetMinimum(fMinimum);
   if (fMaximum != -1111) fHistogram->SetMaximum(fMaximum);

   return fHistogram;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetXmax() const
{
   // Returns the X maximum

   Double_t v = fX[0];
   for (Int_t i=1; i<fNpoints; i++) if (fX[i]>v) v=fX[i];
   return v;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetXmin() const
{
   // Returns the X minimum

   Double_t v = fX[0];
   for (Int_t i=1; i<fNpoints; i++) if (fX[i]<v) v=fX[i];
   return v;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetYmax() const
{
   // Returns the Y maximum

   Double_t v = fY[0];
   for (Int_t i=1; i<fNpoints; i++) if (fY[i]>v) v=fY[i];
   return v;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetYmin() const
{
   // Returns the Y minimum

   Double_t v = fY[0];
   for (Int_t i=1; i<fNpoints; i++) if (fY[i]<v) v=fY[i];
   return v;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetZmax() const
{
   // Returns the Z maximum

   Double_t v = fZ[0];
   for (Int_t i=1; i<fNpoints; i++) if (fZ[i]>v) v=fZ[i];
   return v;
}


//______________________________________________________________________________
 Double_t TGraph2D::GetZmin() const
{
   // Returns the Z minimum

   Double_t v = fZ[0];
   for (Int_t i=1; i<fNpoints; i++) if (fZ[i]<v) v=fZ[i];
   return v;
}


//______________________________________________________________________________
 Double_t TGraph2D::Interpolate(Double_t x, Double_t y)
{
   // Finds the z value at the position (x,y) thanks to
   // the Delaunay interpolation.

   if (fNpoints <= 0) {
      Error("Interpolate", "Empty TGraph2D");
      return 0;
   }

   TGraphDelaunay *dt;

   if(!fHistogram) GetHistogram("empty");

   TList *hl = fHistogram->GetListOfFunctions();
   dt = (TGraphDelaunay*)hl->FindObject("TGraphDelaunay");

   return dt->ComputeZ(x, y);
}


//______________________________________________________________________________
 void TGraph2D::Paint(Option_t *option)
{
   // Paints this 2D graph with its current attributes

   if (fNpoints <= 0) {
      Error("Paint", "Empty TGraph2D");
      return;
   }

   TString opt = option;
   opt.ToLower();
   if (opt.Contains("p") && !opt.Contains("tri")) {
      if (!opt.Contains("pol") &&
          !opt.Contains("sph") &&
          !opt.Contains("psr")) opt.Append("tri0");
   }

   if (opt.Contains("tri0")) {
      GetHistogram("empty");
   } else {
      GetHistogram();
   }

   fHistogram->SetLineColor(GetLineColor());
   fHistogram->SetLineStyle(GetLineStyle());
   fHistogram->SetLineWidth(GetLineWidth());
   fHistogram->SetFillColor(GetFillColor());
   fHistogram->SetFillStyle(GetFillStyle());
   fHistogram->SetMarkerColor(GetMarkerColor());
   fHistogram->SetMarkerStyle(GetMarkerStyle());
   fHistogram->SetMarkerSize(GetMarkerSize());
   fHistogram->Paint(opt.Data());
}


//______________________________________________________________________________
 TH1 *TGraph2D::Project(Option_t *option) const
{
   // Projects a 2-d graph into 1 or 2-d histograms depending on the
   // option parameter
   // option may contain a combination of the characters x,y,z
   // option = "x" return the x projection into a TH1D histogram
   // option = "y" return the y projection into a TH1D histogram
   // option = "xy" return the x versus y projection into a TH2D histogram
   // option = "yx" return the y versus x projection into a TH2D histogram

   if (fNpoints <= 0) {
      Error("Project", "Empty TGraph2D");
      return 0;
   }

   TString opt = option; opt.ToLower();

   Int_t pcase = 0;
   if (opt.Contains("x"))  pcase = 1;
   if (opt.Contains("y"))  pcase = 2;
   if (opt.Contains("xy")) pcase = 3;
   if (opt.Contains("yx")) pcase = 4;

   // Create the projection histogram
   TH1D *h1 = 0;
   TH2D *h2 = 0;
   Int_t nch = strlen(GetName()) +opt.Length() +2;
   char *name = new char[nch];
   sprintf(name,"%s_%s",GetName(),option);
   nch = strlen(GetTitle()) +opt.Length() +2;
   char *title = new char[nch];
   sprintf(title,"%s_%s",GetTitle(),option);

   Double_t hxmin = GetXmin();
   Double_t hxmax = GetXmax();
   Double_t hymin = GetYmin();
   Double_t hymax = GetYmax();

   switch (pcase) {
      case 1:
         // "x"
         h1 = new TH1D(name,title,fNpx,hxmin,hxmax);
         break;
      case 2:
         // "y"
         h1 = new TH1D(name,title,fNpy,hymin,hymax);
         break;
      case 3:
         // "xy"
         h2 = new TH2D(name,title,fNpx,hxmin,hxmax,fNpy,hymin,hymax);
         break;
      case 4:
         // "yx"
         h2 = new TH2D(name,title,fNpy,hymin,hymax,fNpx,hxmin,hxmax);
         break;
   }

   delete [] name;
   delete [] title;
   TH1 *h = h1;
   if (h2) h = h2;
   if (h == 0) return 0;

   // Fill the projected histogram
   Double_t entries = 0;
   for (Int_t N=0; N<fNpoints; N++) {
      switch (pcase) {
         case 1:
            // "x"
            h1->Fill(fX[N], fZ[N]);
            break;
         case 2:
            // "y"
            h1->Fill(fY[N], fZ[N]);
            break;
         case 3:
            // "xy"
            h2->Fill(fX[N], fY[N], fZ[N]);
            break;
         case 4:
            // "yx"
            h2->Fill(fY[N], fX[N], fZ[N]);
            break;
      }
      entries += fZ[N];
   }
   h->SetEntries(entries);
   return h;
}


//______________________________________________________________________________
 Int_t TGraph2D::RemovePoint(Int_t ipoint)
{
   // Deletes point number ipoint

   if (ipoint < 0) return -1;
   if (ipoint >= fNpoints) return -1;

   fNpoints--;
   Double_t *newX = new Double_t[fNpoints];
   Double_t *newY = new Double_t[fNpoints];
   Double_t *newZ = new Double_t[fNpoints];
   Int_t j = -1;
   for (Int_t i=0;i<fNpoints+1;i++) {
      if (i == ipoint) continue;
      j++;
      newX[j] = fX[i];
      newY[j] = fY[i];
      newZ[j] = fZ[i];
   }
   delete [] fX;
   delete [] fY;
   delete [] fZ;
   fX = newX;
   fY = newY;
   fZ = newZ;
   if (fHistogram) {delete fHistogram; fHistogram = 0;}
   return ipoint;
}


//______________________________________________________________________________
 void TGraph2D::SavePrimitive(ofstream &out, Option_t *option)
{
   // Saves primitive as a C++ statement(s) on output stream out

   char quote = '"';
   out<<"   "<<endl;
   if (gROOT->ClassSaved(TGraph2D::Class())) {
      out<<"   ";
   } else {
      out<<"   TGraph2D *";
   }

   out<<"graph = new TGraph2D("<<fNpoints<<");"<<endl;
   out<<"   graph->SetName("<<quote<<GetName()<<quote<<");"<<endl;
   out<<"   graph->SetTitle("<<quote<<GetTitle()<<quote<<");"<<endl;

   if (fDirectory == 0) {
      out<<"   "<<GetName()<<"->SetDirectory(0);"<<endl;
   }

   SaveFillAttributes(out,"graph",0,1001);
   SaveLineAttributes(out,"graph",1,1,1);
   SaveMarkerAttributes(out,"graph",1,1,1);

   for (Int_t i=0;i<fNpoints;i++) {
      out<<"   graph->SetPoint("<<i<<","<<fX[i]<<","<<fY[i]<<","<<fZ[i]<<");"<<endl;
   }

   // save list of functions
   TIter next(fFunctions);
   TObject *obj;
   while ((obj=next())) {
      obj->SavePrimitive(out,"nodraw");
      out<<"   graph->GetListOfFunctions()->Add("<<obj->GetName()<<");"<<endl;
      if (obj->InheritsFrom("TPaveStats")) {
         out<<"   ptstats->SetParent(graph->GetListOfFunctions());"<<endl;
      }
   }

   out<<"   graph->Draw("<<quote<<option<<quote<<");"<<endl;
}


//______________________________________________________________________________
 void TGraph2D::SetDirectory(TDirectory *dir)
{
   // By default when an 2D graph is created, it is added to the list
   // of 2D graph objects in the current directory in memory.
   // Remove reference to this 2D graph from current directory and add
   // reference to new directory dir. dir can be 0 in which case the
   // 2D graph does not belong to any directory.

   if (fDirectory == dir) return;
   if (fDirectory) fDirectory->GetList()->Remove(this);
   fDirectory = dir;
   if (fDirectory) fDirectory->GetList()->Add(this);
}


//______________________________________________________________________________
 void TGraph2D::SetMargin(Double_t m)
{
   // Sets the extra space (in %) around interpolated area for the 2D histogram

   if (m<0 || m>1) {
      Warning("SetMargin","The margin must be >= 0 && <= 1, fMargin set to 0.1");
      fMargin = 0.1;
   } else {
      fMargin = m;
   }
   if (fHistogram) {delete fHistogram; fHistogram = 0;}
}


//______________________________________________________________________________
 void TGraph2D::SetMarginBinsContent(Double_t z)
{
   // Sets the histogram bin height for points lying outside the TGraphDelaunay
   // convex hull ie: the bins in the margin.

   fZout = z;
   if (fHistogram) {delete fHistogram; fHistogram = 0;}
}


//______________________________________________________________________________
 void TGraph2D::SetMaximum(Double_t maximum)
{
   fMaximum = maximum;
   GetHistogram()->SetMaximum(maximum);
}


//______________________________________________________________________________
 void TGraph2D::SetMinimum(Double_t minimum)
{
   fMinimum = minimum;
   GetHistogram()->SetMinimum(minimum);
}


//______________________________________________________________________________
 void TGraph2D::SetName(const char *name)
{
   // Changes the name of this 2D graph

   //  2D graphs are named objects in a THashList.
   //  We must update the hashlist if we change the name
   if (fDirectory) fDirectory->GetList()->Remove(this);
   fName = name;
   if (fDirectory) fDirectory->GetList()->Add(this);
}


//______________________________________________________________________________
 void TGraph2D::SetNpx(Int_t npx)
{
   // Sets the number of bins along X used to draw the function

   if (npx < 4) {
      Warning("SetNpx","Number of points must be >4 && < 500, fNpx set to 4");
      fNpx = 4;
   } else if(npx > 500) {
      Warning("SetNpx","Number of points must be >4 && < 500, fNpx set to 500");
      fNpx = 100000;
   } else {
      fNpx = npx;
   }
   if (fHistogram) {delete fHistogram; fHistogram = 0;}
}


//______________________________________________________________________________
 void TGraph2D::SetNpy(Int_t npy)
{
   // Sets the number of bins along Y used to draw the function

   if (npy < 4) {
      Warning("SetNpy","Number of points must be >4 && < 500, fNpy set to 4");
      fNpy = 4;
   } else if(npy > 500) {
      Warning("SetNpy","Number of points must be >4 && < 500, fNpy set to 500");
      fNpy = 100000;
   } else {
      fNpy = npy;
   }
   if (fHistogram) {delete fHistogram; fHistogram = 0;}
}


//______________________________________________________________________________
 void TGraph2D::SetPoint(Int_t n, Double_t x, Double_t y, Double_t z)
{
   // Sets point number n.
   // If n is greater than the current size, the arrays are automatically
   // extended.

   if (n < 0) return;

   if (!fX || !fY || !fZ || n >= fSize) {
   // re-allocate the object
      Int_t newN = TMath::Max(2*fSize,n+1);
      Double_t *savex = new Double_t [newN];
      Double_t *savey = new Double_t [newN];
      Double_t *savez = new Double_t [newN];
      if (fX && fSize) {
         memcpy(savex,fX,fSize*sizeof(Double_t));
         memset(&savex[fSize],0,(newN-fSize)*sizeof(Double_t));
         delete [] fX;
      }
      if (fY && fSize) {
         memcpy(savey,fY,fSize*sizeof(Double_t));
         memset(&savey[fSize],0,(newN-fSize)*sizeof(Double_t));
         delete [] fY;
      }
      if (fZ && fSize) {
         memcpy(savez,fZ,fSize*sizeof(Double_t));
         memset(&savez[fSize],0,(newN-fSize)*sizeof(Double_t));
         delete [] fZ;
      }
      fX    = savex;
      fY    = savey;
      fZ    = savez;
      fSize = newN;
   }
   fX[n]    = x;
   fY[n]    = y;
   fZ[n]    = z;
   fNpoints = TMath::Max(fNpoints,n+1);
}


//______________________________________________________________________________
 void TGraph2D::SetTitle(const char* title)
{
   // Sets graph title

   fTitle = title;
   if (fHistogram) fHistogram->SetTitle(title);
}


//_______________________________________________________________________
 void TGraph2D::Streamer(TBuffer &b)
{
   // Stream a class object

   if (b.IsReading()) {
      UInt_t R__s, R__c;
      Version_t R__v = b.ReadVersion(&R__s, &R__c);
      TGraph2D::Class()->ReadBuffer(b, this, R__v, R__s, R__c);

      if (!gROOT->ReadingObject()) {
         fDirectory = gDirectory;
         if (!gDirectory->GetList()->FindObject(this)) gDirectory->Append(this);
      }
      ResetBit(kCanDelete);
   } else {
      TGraph2D::Class()->WriteBuffer(b,this);
   }
}