    // @(#)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. // /*
This applet gives a nice practical view of Delaunay triangulation and Voronoi diagram. */ //

//______________________________________________________________________________ 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; kUserHisto = kFALSE; } //______________________________________________________________________________ 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; kUserHisto = kFALSE; 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 // Predefined functions such as gaus, expo and poln are automatically // created by ROOT. // fname can also be a formula, accepted by the linear fitter (linear parts divided // by "++" sign), for example "x++sin(x)" for fitting "[0]*x+[1]*sin(x)" char *linear; linear= (char*)strstr(fname, "++"); TF2 *f2=0; if (linear) f2=new TF2(fname, fname); else{ f2 = (TF2*)gROOT->GetFunction(fname); if (!f2) { Printf("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) // = "C" In case of linear fitting, not calculate the chisquare // (saves time) // = "ROB" In case of linear fitting, compute the LTS regression // coefficients (robust(resistant) regression), using // the default fraction of good points // "ROB=0.x" - compute the LTS regression coefficients, using // 0.x as a fraction of good points // // 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; TString opt = option; opt.ToUpper(); Double_t h=0; opt.ReplaceAll("ROB", "H"); //for robust fitting, see if # of good points is defined if (opt.Contains("H=0.")) { int start = opt.Index("H=0."); int numpos = start + strlen("H=0."); int numlen = 0; int len = opt.Length(); while( (numpos+numlen<len) && isdigit(opt[numpos+numlen]) ) numlen++; TString num = opt(numpos,numlen); opt.Remove(start+strlen("H"),strlen("=0.")+numlen); h = atof(num.Data()); h*=TMath::Power(10, -numlen); } 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; if (opt.Contains("C")) fitOption.Nochisq = 1; if (opt.Contains("H")) fitOption.Robust = 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 Bool_t linear = f2->IsLinear(); if (fitOption.Bound || fitOption.User || fitOption.Errors) linear = kFALSE; char l[]="TLinearFitter"; Int_t strdiff = 0; Bool_t IsSet = kFALSE; if (TVirtualFitter::GetFitter()){ //Is a fitter already set? Is it linear? IsSet=kTRUE; strdiff = strcmp(TVirtualFitter::GetFitter()->IsA()->GetName(), l); } if (linear){ // TClass *cl = gROOT->GetClass("TLinearFitter"); if (IsSet && strdiff!=0) { delete TVirtualFitter::GetFitter(); IsSet=kFALSE; } if (!IsSet) { //TLinearFitter *lf=(TLinearFitter *)cl->New(); TVirtualFitter::SetFitter((TVirtualFitter *)cl->New()); } } else { if (IsSet && strdiff==0){ delete TVirtualFitter::GetFitter(); IsSet=kFALSE; } if (!IsSet) TVirtualFitter::SetFitter(0); } 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; ///} if (linear){ if (fitOption.Robust) grFitter->ExecuteCommand("FitGraph2D", &h, 0); else grFitter->ExecuteCommand("FitGraph2D", 0, 0); } else { // 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 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) { if (!kUserHisto) { delete fHistogram; fHistogram = 0; } } else { return fHistogram; } } Double_t hxmax, hymax, hxmin, hymin; // Book fHistogram if needed. It is not added in the current directory if (!kUserHisto) { Bool_t add = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); Double_t xmax = GetXmax(); Double_t ymax = GetYmax(); Double_t xmin = GetXmin(); Double_t ymin = GetYmin(); hxmin = xmin-fMargin*(xmax-xmin); hymin = ymin-fMargin*(ymax-ymin); hxmax = xmax+fMargin*(xmax-xmin); hymax = ymax+fMargin*(ymax-ymin); fHistogram = new TH2D(GetName(),GetTitle(), fNpx ,hxmin, hxmax, fNpy, hymin, hymax); TH1::AddDirectory(add); fHistogram->SetBit(TH1::kNoStats); } else { hxmin = fHistogram->GetXaxis()->GetXmin(); hymin = fHistogram->GetYaxis()->GetXmin(); hxmax = fHistogram->GetXaxis()->GetXmax(); hymax = fHistogram->GetYaxis()->GetXmax(); } // 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::SetHistogram(TH2 *h) { // Sets the histogram to be filled kUserHisto = kTRUE; fHistogram = (TH2D*)h; fNpx = h->GetNbinsX(); fNpy = h->GetNbinsY(); } //______________________________________________________________________________ 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 = 500; } 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 = 500; } 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); } }