// @(#)root/graf:$Name: $:$Id: TGraph2D.cxx,v 1.00
// Author: Olivier Couet, Luke Jones (Royal Holloway, University of London)
/*************************************************************************
* 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 "TGraph2D.h"
#include "TMath.h"
#include "TPolyLine.h"
#include "TPolyMarker.h"
#include "TVirtualPad.h"
#include "TVirtualFitter.h"
#include "TView.h"
#include "THLimitsFinder.h"
#include "TStyle.h"
ClassImp(TGraph2D)
//______________________________________________________________________________
//
// A Graph2D is a graphics object made of three arrays X, Y and Z with the same
// number of points each. Graph2D uses Delaunay triangulation to draw the X, Y
// Z arrays. This triangulation code derives from an implementation done by
// Luke Jones (Royal Holloway, University of London) in April 2002 in the PAW
// context.
//
// This class has different constructors:
//
// 1) With an array 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 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
Build(100);
}
//______________________________________________________________________________
TGraph2D::TGraph2D(Int_t n, Int_t *x, Int_t *y, Int_t *z, Option_t *)
: 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, Option_t *)
: 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, Option_t *)
: 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(const char *name,const char *title,
Int_t n, Double_t *x, Double_t *y, Double_t *z, Option_t *)
: 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, Option_t *)
: 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 (fPTried) delete [] fPTried;
if (fNTried) delete [] fNTried;
if (fMTried) delete [] fMTried;
if (fHullPoints) delete [] fHullPoints;
if (fOrder) delete [] fOrder;
if (fDist) delete [] fDist;
if (fXN) delete [] fXN;
if (fYN) delete [] fYN;
if (fFunctions) {
fFunctions->SetBit(kInvalidObject);
fFunctions->Delete();
delete fFunctions;
}
if (fDirectory) {
if (!fDirectory->TestBit(TDirectory::kCloseDirectory))
fDirectory->GetList()->Remove(this);
}
fX = 0;
fY = 0;
fZ = 0;
fHistogram = 0;
fPTried = 0;
fNTried = 0;
fMTried = 0;
fHullPoints = 0;
fOrder = 0;
fDist = 0;
fXN = 0;
fYN = 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,
fTriedSize = 0;
fMargin = 0.;
fNpx = 40;
fNpy = 40;
fZout = 0.;
fNdt = 0;
fNhull = 0;
fDirectory = 0;
fFunctions = 0;
fHistogram = 0;
fMaximum = -1111;
fMinimum = -1111;
fHullPoints = 0;
fXN = 0;
fYN = 0;
fOrder = 0;
fDist = 0;
fPTried = 0;
fNTried = 0;
fMTried = 0;
fGridLevels = 0;
fNbLevels = 0;
fView = 0;
SetMaxIter();
fX = new Double_t[fSize];
fY = new Double_t[fSize];
fZ = new Double_t[fSize];
fFunctions = new TList;
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;
}
}
//______________________________________________________________________________
Double_t TGraph2D::ComputeZ(Double_t xx, Double_t yy)
{
// Finds the Delaunay triangle that the point (xx,yy) sits in (if any) and
// calculate a z-value for it by linearly interpolating the z-values that
// make up that triangle.
Double_t thevalue;
Int_t IT, ntris_tried, P, N, M;
Int_t I,J,K,L,Z,F,D,O1,O2,A,B,T1,T2,T3;
Int_t ndegen=0,degen=0,fdegen=0,o1degen=0,o2degen=0;
Double_t vxN,vyN;
Double_t d1,d2,d3,c1,c2,dko1,dko2,dfo1;
Double_t dfo2,sin_sum,cfo1k,co2o1k,co2o1f;
Bool_t shouldbein;
Double_t dx1,dx2,dx3,dy1,dy2,dy3,U,V,dxz[3],dyz[3];
// create vectors needed for sorting
if (!fOrder) {
fOrder = new Int_t[fNpoints];
fDist = new Double_t[fNpoints];
}
// the input point will be point zero.
fXN[0] = xx;
fYN[0] = yy;
// set the output value to the default value for now
thevalue = fZout;
// some counting
ntris_tried = 0;
// no point in proceeding if xx or yy are silly
if ((xx>fXNmax) || (xx<fXNmin) || (yy>fYNmax) || (yy<fYNmin)) return thevalue;
// check existing Delaunay triangles for a good one
for (IT=1; IT<=fNdt; IT++) {
P = fPTried[IT-1];
N = fNTried[IT-1];
M = fMTried[IT-1];
// P, N and M form a previously found Delaunay triangle, does it
// enclose the point?
if (Enclose(P,N,M,0)) {
// yes, we have the triangle
thevalue = InterpolateOnPlane(P,N,M,0);
return thevalue;
}
}
// is this point inside the convex hull?
shouldbein = InHull(0,-1);
if (!shouldbein) return thevalue;
// it must be in a Delaunay triangle - find it...
// order mass points by distance in mass plane from desired point
for (IT=1; IT<=fNpoints; IT++) {
vxN = fXN[IT];
vyN = fYN[IT];
fDist[IT-1] = TMath::Sqrt((xx-vxN)*(xx-vxN)+(yy-vyN)*(yy-vyN));
}
// sort array 'fDist' to find closest points
TMath::Sort(fNpoints, fDist, fOrder, kFALSE);
for (IT=0; IT<fNpoints; IT++) fOrder[IT]++;
// loop over triplets of close points to try to find a triangle that
// encloses the point.
for (K=3; K<=fNpoints; K++) {
M = fOrder[K-1];
for (J=2; J<=K-1; J++) {
N = fOrder[J-1];
for (I=1; I<=J-1; I++) {
P = fOrder[I-1];
if (ntris_tried > fMaxIter) {
// perhaps this point isn't in the hull after all
/// Warning("ComputeZ",
/// "Abandoning the effort to find a Delaunay triangle (and thus interpolated Z-value) for point %g %g"
/// ,xx,yy);
return thevalue;
}
ntris_tried++;
// check the points aren't colinear
d1 = TMath::Sqrt((fXN[P]-fXN[N])*(fXN[P]-fXN[N])+(fYN[P]-fYN[N])*(fYN[P]-fYN[N]));
d2 = TMath::Sqrt((fXN[P]-fXN[M])*(fXN[P]-fXN[M])+(fYN[P]-fYN[M])*(fYN[P]-fYN[M]));
d3 = TMath::Sqrt((fXN[N]-fXN[M])*(fXN[N]-fXN[M])+(fYN[N]-fYN[M])*(fYN[N]-fYN[M]));
if ((d1+d2<=d3) || (d1+d3<=d2) || (d2+d3<=d1)) goto L90;
// does the triangle enclose the point?
if (!Enclose(P,N,M,0)) goto L90;
// is it a Delaunay triangle? (ie. are there any other points
// inside the circle that is defined by its vertices?)
// test the triangle for Delaunay'ness
// loop over all other points testing each to see if it's
// inside the triangle's circle
ndegen = 0;
for ( Z=1; Z<=fNpoints; Z++) {
if ((Z==P) || (Z==N) || (Z==M)) goto L50;
// An easy first check is to see if point Z is inside the triangle
// (if it's in the triangle it's also in the circle)
// point Z cannot be inside the triangle if it's further from (xx,yy)
// than the furthest pointing making up the triangle - test this
for (L=1; L<=fNpoints; L++) {
if (fOrder[L-1] == Z) {
if ((L<I) || (L<J) || (L<K)) {
// point Z is nearer to (xx,yy) than M, N or P - it could be in the
// triangle so call enclose to find out
// if it is inside the triangle this can't be a Delaunay triangle
if (Enclose(P,N,M,Z)) goto L90;
} else {
// there's no way it could be in the triangle so there's no point
// calling enclose
goto L1;
}
}
}
// is point Z colinear with any pair of the triangle points?
L1:
if (((fXN[P]-fXN[Z])*(fYN[P]-fYN[N])) == ((fYN[P]-fYN[Z])*(fXN[P]-fXN[N]))) {
// Z is colinear with P and N
A = P;
B = N;
} else if (((fXN[P]-fXN[Z])*(fYN[P]-fYN[M])) == ((fYN[P]-fYN[Z])*(fXN[P]-fXN[M]))) {
// Z is colinear with P and M
A = P;
B = M;
} else if (((fXN[N]-fXN[Z])*(fYN[N]-fYN[M])) == ((fYN[N]-fYN[Z])*(fXN[N]-fXN[M]))) {
// Z is colinear with N and M
A = N;
B = M;
} else {
A = 0;
B = 0;
}
if (A != 0) {
// point Z is colinear with 2 of the triangle points, if it lies
// between them it's in the circle otherwise it's outside
if (fXN[A] != fXN[B]) {
if (((fXN[Z]-fXN[A])*(fXN[Z]-fXN[B])) < 0) {
goto L90;
} else if (((fXN[Z]-fXN[A])*(fXN[Z]-fXN[B])) == 0) {
// At least two points are sitting on top of each other, we will
// treat these as one and not consider this a 'multiple points lying
// on a common circle' situation. It is a sign something could be
// wrong though, especially if the two coincident points have
// different fZ's. If they don't then this is harmless.
Warning("ComputeZ", "Two of these three points are coincident %d %d %d",A,B,Z);
}
} else {
if (((fYN[Z]-fYN[A])*(fYN[Z]-fYN[B])) < 0) {
goto L90;
} else if (((fYN[Z]-fYN[A])*(fYN[Z]-fYN[B])) == 0) {
// At least two points are sitting on top of each other - see above.
Warning("ComputeZ", "Two of these three points are coincident %d %d %d",A,B,Z);
}
}
// point is outside the circle, move to next point
goto L50;
}
// if point Z were to look at the triangle, which point would it see
// lying between the other two? (we're going to form a quadrilateral
// from the points, and then demand certain properties of that
// quadrilateral)
dxz[0] = fXN[P]-fXN[Z];
dyz[0] = fYN[P]-fYN[Z];
dxz[1] = fXN[N]-fXN[Z];
dyz[1] = fYN[N]-fYN[Z];
dxz[2] = fXN[M]-fXN[Z];
dyz[2] = fYN[M]-fYN[Z];
for(L=1; L<=3; L++) {
dx1 = dxz[L-1];
dx2 = dxz[L%3];
dx3 = dxz[(L+1)%3];
dy1 = dyz[L-1];
dy2 = dyz[L%3];
dy3 = dyz[(L+1)%3];
U = (dy3*dx2-dx3*dy2)/(dy1*dx2-dx1*dy2);
V = (dy3*dx1-dx3*dy1)/(dy2*dx1-dx2*dy1);
if ((U>=0) && (V>=0)) {
// vector (dx3,dy3) is expressible as a sum of the other two vectors
// with positive coefficents -> i.e. it lies between the other two vectors
if (L == 1) {
F = M;
O1 = P;
O2 = N;
} else if (L == 2) {
F = P;
O1 = N;
O2 = M;
} else {
F = N;
O1 = M;
O2 = P;
}
goto L2;
}
}
Error("ComputeZ", "Should not get to here");
// may as well soldier on
F = M;
O1 = P;
O2 = N;
L2:
// this is not a valid quadrilateral if the diagonals don't cross,
// check that points F and Z lie on opposite side of the line O1-O2,
// this is true if the angle F-O1-Z is greater than O2-O1-Z and O2-O1-F
cfo1k = ((fXN[F]-fXN[O1])*(fXN[Z]-fXN[O1])+(fYN[F]-fYN[O1])*(fYN[Z]-fYN[O1]))/
TMath::Sqrt(((fXN[F]-fXN[O1])*(fXN[F]-fXN[O1])+(fYN[F]-fYN[O1])*(fYN[F]-fYN[O1]))*
((fXN[Z]-fXN[O1])*(fXN[Z]-fXN[O1])+(fYN[Z]-fYN[O1])*(fYN[Z]-fYN[O1])));
co2o1k = ((fXN[O2]-fXN[O1])*(fXN[Z]-fXN[O1])+(fYN[O2]-fYN[O1])*(fYN[Z]-fYN[O1]))/
TMath::Sqrt(((fXN[O2]-fXN[O1])*(fXN[O2]-fXN[O1])+(fYN[O2]-fYN[O1])*(fYN[O2]-fYN[O1]))*
((fXN[Z]-fXN[O1])*(fXN[Z]-fXN[O1]) + (fYN[Z]-fYN[O1])*(fYN[Z]-fYN[O1])));
co2o1f = ((fXN[O2]-fXN[O1])*(fXN[F]-fXN[O1])+(fYN[O2]-fYN[O1])*(fYN[F]-fYN[O1]))/
TMath::Sqrt(((fXN[O2]-fXN[O1])*(fXN[O2]-fXN[O1])+(fYN[O2]-fYN[O1])*(fYN[O2]-fYN[O1]))*
((fXN[F]-fXN[O1])*(fXN[F]-fXN[O1]) + (fYN[F]-fYN[O1])*(fYN[F]-fYN[O1]) ));
if ((cfo1k>co2o1k) || (cfo1k>co2o1f)) {
// not a valid quadrilateral - point Z is definitely outside the circle
goto L50;
}
// calculate the 2 internal angles of the quadrangle formed by joining
// points Z and F to points O1 and O2, at Z and F. If they sum to less
// than 180 degrees then Z lies outside the circle
dko1 = TMath::Sqrt((fXN[Z]-fXN[O1])*(fXN[Z]-fXN[O1])+(fYN[Z]-fYN[O1])*(fYN[Z]-fYN[O1]));
dko2 = TMath::Sqrt((fXN[Z]-fXN[O2])*(fXN[Z]-fXN[O2])+(fYN[Z]-fYN[O2])*(fYN[Z]-fYN[O2]));
dfo1 = TMath::Sqrt((fXN[F]-fXN[O1])*(fXN[F]-fXN[O1])+(fYN[F]-fYN[O1])*(fYN[F]-fYN[O1]));
dfo2 = TMath::Sqrt((fXN[F]-fXN[O2])*(fXN[F]-fXN[O2])+(fYN[F]-fYN[O2])*(fYN[F]-fYN[O2]));
c1 = ((fXN[Z]-fXN[O1])*(fXN[Z]-fXN[O2])+(fYN[Z]-fYN[O1])*(fYN[Z]-fYN[O2]))/dko1/dko2;
c2 = ((fXN[F]-fXN[O1])*(fXN[F]-fXN[O2])+(fYN[F]-fYN[O1])*(fYN[F]-fYN[O2]))/dfo1/dfo2;
sin_sum = c1*TMath::Sqrt(1-c2*c2)+c2*TMath::Sqrt(1-c1*c1);
// sin_sum doesn't always come out as zero when it should do.
if (sin_sum < -1.E-6) {
// Z is inside the circle, this is not a Delaunay triangle
goto L90;
} else if (TMath::Abs(sin_sum) <= 1.E-6) {
// point Z lies on the circumference of the circle (within rounding errors)
// defined by the triangle, so there is potential for degeneracy in the
// triangle set (Delaunay triangulation does not give a unique way to split
// a polygon whose points lie on a circle into constituent triangles). Make
// a note of the additional point number.
ndegen++;
degen = Z;
fdegen = F;
o1degen = O1;
o2degen = O2;
}
L50:
continue;
}
// This is a good triangle
if (ndegen > 0) {
// but is degenerate with at least one other,
// haven't figured out what to do if more than 4 points are involved
if (ndegen > 1) {
Error("ComputeZ",
"More than 4 points lying on a circle. No decision making process formulated for triangulating this region in a non-arbitrary way %d %d %d %d",
P,N,M,degen);
return thevalue;
}
// we have a quadrilateral which can be split down either diagonal
// (D<->F or O1<->O2) to form valid Delaunay triangles. Choose diagonal
// with highest average z-value. Whichever we choose we will have
// verified two triangles as good and two as bad, only note the good ones
D = degen;
F = fdegen;
O1 = o1degen;
O2 = o2degen;
if ((fZ[O1-1]+fZ[O2-1]) > (fZ[D-1]+fZ[F-1])) {
// best diagonalisation of quadrilateral is current one, we have
// the triangle
T1 = P;
T2 = N;
T3 = M;
// file the good triangles
FileIt(P, N, M);
FileIt(D, O1, O2);
} else {
// use other diagonal to split quadrilateral, use triangle formed by
// point F, the degnerate point D and whichever of O1 and O2 create
// an enclosing triangle
T1 = F;
T2 = D;
if (Enclose(F,D,O1,0)) {
T3 = O1;
} else {
T3 = O2;
}
// file the good triangles
FileIt(F, D, O1);
FileIt(F, D, O2);
}
} else {
// this is a Delaunay triangle, file it
FileIt(P, N, M);
T1 = P;
T2 = N;
T3 = M;
}
// do the interpolation
thevalue = InterpolateOnPlane(T1,T2,T3,0);
return thevalue;
L90:
continue;
}
}
}
if (shouldbein) {
Error("ComputeZ",
"Point outside hull when expected inside: this point could be dodgy %g %g %d",
xx, yy, ntris_tried);
}
return thevalue;
}
//______________________________________________________________________________
void TGraph2D::DefineGridLevels()
{
// Define the grid levels drawn on the triangles.
// The grid levels are aligned on the Z axis' main tick marks.
// The function assumes that fView has been defined.
Int_t i, nbins;
Double_t BinLow, BinHigh, BinWidth;
// Find the main tick marks positions.
Double_t *rmin = fView->GetRmin();
Double_t *rmax = fView->GetRmax();
Int_t ndivz = fHistogram->GetZaxis()->GetNdivisions()%100;
if (ndivz > 0) {
THLimitsFinder::Optimize(rmin[2], rmax[2], ndivz,
BinLow, BinHigh, nbins, BinWidth, " ");
} else {
nbins = TMath::Abs(ndivz);
BinLow = rmin[2];
BinHigh = rmax[2];
BinWidth = (BinHigh-BinLow)/nbins;
}
// Define the grid levels
fNbLevels = nbins+1;
fGridLevels = new Double_t[fNbLevels];
for (i = 0; i < fNbLevels; ++i) fGridLevels[i] = BinLow+i*BinWidth;
}
//______________________________________________________________________________
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());
}
//______________________________________________________________________________
Bool_t TGraph2D::Enclose(Int_t T1, Int_t T2, Int_t T3, Int_t Ex) const
{
// Is point E inside the triangle T1-T2-T3 ?
Int_t E=0,A=0,B=0;
Double_t dx1,dx2,dx3,dy1,dy2,dy3,U,V;
Bool_t enclose = kFALSE;
E = TMath::Abs(Ex);
// First ask if point E is colinear with any pair of the triangle points
A = 0;
if (((fXN[T1]-fXN[E])*(fYN[T1]-fYN[T2])) == ((fYN[T1]-fYN[E])*(fXN[T1]-fXN[T2]))) {
// E is colinear with T1 and T2
A = T1;
B = T2;
} else if (((fXN[T1]-fXN[E])*(fYN[T1]-fYN[T3])) == ((fYN[T1]-fYN[E])*(fXN[T1]-fXN[T3]))) {
// E is colinear with T1 and T3
A = T1;
B = T3;
} else if (((fXN[T2]-fXN[E])*(fYN[T2]-fYN[T3])) == ((fYN[T2]-fYN[E])*(fXN[T2]-fXN[T3]))) {
// E is colinear with T2 and T3
A = T2;
B = T3;
}
if (A != 0) {
// point E is colinear with 2 of the triangle points, if it lies
// between them it's in the circle otherwise it's outside
if (fXN[A] != fXN[B]) {
if (((fXN[E]-fXN[A])*(fXN[E]-fXN[B])) <= 0) {
enclose = kTRUE;
return enclose;
}
} else {
if (((fYN[E]-fYN[A])*(fYN[E]-fYN[B])) <= 0) {
enclose = kTRUE;
return enclose;
}
}
// point is outside the triangle
return enclose;
}
// E is not colinear with any pair of triangle points, if it is inside
// the triangle then the vector from E to one of the corners must be
// expressible as a sum with positive coefficients of the vectors from
// the two other corners to E. Say vector3=U*vector1+V*vector2
// vector1==T1->E
dx1 = fXN[E]-fXN[T1];
dy1 = fYN[E]-fYN[T1];
// vector2==T2->E
dx2 = fXN[E]-fXN[T2];
dy2 = fYN[E]-fYN[T2];
// vector3==E->T3
dx3 = fXN[T3]-fXN[E];
dy3 = fYN[T3]-fYN[E];
U = (dx2*dy3-dx3*dy2)/(dx2*dy1-dx1*dy2);
V = (dx1*dy3-dx3*dy1)/(dx1*dy2-dx2*dy1);
if ((U>=0) && (V>=0)) enclose = kTRUE;
return enclose;
}
//______________________________________________________________________________
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);
}
//______________________________________________________________________________
void TGraph2D::FileIt(Int_t P, Int_t N, Int_t M)
{
// Files the triangle defined by the 3 vertices P, N and M into the
// fxTried arrays. If these arrays are to small they are automatically
// expanded.
Bool_t swap;
Int_t tmp, Ps = P, Ns = N, Ms = M;
// order the vertices before storing them
L1:
swap = kFALSE;
if (Ns > Ps) { tmp = Ps; Ps = Ns; Ns = tmp; swap = kTRUE;}
if (Ms > Ns) { tmp = Ns; Ns = Ms; Ms = tmp; swap = kTRUE;}
if (swap) goto L1;
// expand the triangles storage if needed
if (fNdt> fTriedSize) {
Int_t newN = 2*fTriedSize;
Int_t *savep = new Int_t [newN];
Int_t *saven = new Int_t [newN];
Int_t *savem = new Int_t [newN];
memcpy(savep,fPTried,fTriedSize*sizeof(Int_t));
memset(&savep[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
delete [] fPTried;
memcpy(saven,fNTried,fTriedSize*sizeof(Int_t));
memset(&saven[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
delete [] fNTried;
memcpy(savem,fMTried,fTriedSize*sizeof(Int_t));
memset(&savem[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
delete [] fMTried;
fPTried = savep;
fNTried = saven;
fMTried = savem;
fTriedSize = newN;
}
// store a new Delaunay triangle
fNdt++;
fPTried[fNdt-1] = Ps;
fNTried[fNdt-1] = Ns;
fMTried[fNdt-1] = Ms;
}
//______________________________________________________________________________
void TGraph2D::FindAllTriangles()
{
// Attempt to find all the Delaunay triangles of the point set. It is not
// guaranteed that it will fully succeed, and no check is made that it has
// fully succeeded (such a check would be possible by referencing the points
// that make up the convex hull). The method is to check if each triangle
// shares all three of its sides with other triangles. If not, a point is
// generated just outside the triangle on the side(s) not shared, and a new
// triangle is found for that point. If this method is not working properly
// (many triangles are not being found) it's probably because the new points
// are too far beyond or too close to the non-shared sides. Fiddling with the
// size of the `alittlebit' parameter may help.
Double_t xcntr,ycntr,xm,ym,xx,yy;
Double_t sx,sy,nx,ny,mx,my,mdotn,nn,A;
Int_t T1,T2,Pa,Na,Ma,Pb,Nb,Mb,P1=0,P2=0,M,N,P3=0;
Bool_t s[3];
Double_t alittlebit = 0.0001;
// start with a point that is guaranteed to be inside the hull (the
// centre of the hull)
xcntr = 0;
ycntr = 0;
for (N=1; N<=fNhull; N++) {
xcntr = xcntr+fXN[fHullPoints[N-1]];
ycntr = ycntr+fYN[fHullPoints[N-1]];
}
xcntr = xcntr/fNhull;
ycntr = ycntr/fNhull;
// and calculate it's triangle
ComputeZ(xcntr,ycntr);
// loop over all Delaunay triangles (including those constantly being
// produced within the loop) and check to see if their 3 sides also
// correspond to the sides of other Delaunay triangles, i.e. that they
// have all their neighbours.
T1 = 1;
while (T1 <= fNdt) {
// get the three points that make up this triangle
Pa = fPTried[T1-1];
Na = fNTried[T1-1];
Ma = fMTried[T1-1];
// produce three integers which will represent the three sides
s[0] = kFALSE;
s[1] = kFALSE;
s[2] = kFALSE;
// loop over all other Delaunay triangles
for (T2=1; T2<=fNdt; T2++) {
if (T2 != T1) {
// get the points that make up this triangle
Pb = fPTried[T2-1];
Nb = fNTried[T2-1];
Mb = fMTried[T2-1];
// and generate integers which represent its sides
// (note that Mb>Nb>Pb - see routine triencode)
// do triangles T1 and T2 share a side?
if ((Pa==Pb && Na==Nb) || (Pa==Pb && Na==Mb) || (Pa==Nb && Na==Mb)) {
// they share side 1
s[0] = kTRUE;
} else if ((Pa==Pb && Ma==Nb) || (Pa==Pb && Ma==Mb) || (Pa==Nb && Ma==Mb)) {
// they share side 2
s[1] = kTRUE;
} else if ((Na==Pb && Ma==Nb) || (Na==Pb && Ma==Mb) || (Na==Nb && Ma==Mb)) {
// they share side 3
s[2] = kTRUE;
}
}
// if T1 shares all its sides with other Delaunay triangles then
// forget about it
if (s[0] && s[1] && s[2]) continue;
}
// Looks like T1 is missing a neighbour on at least one side.
// For each side, take a point a little bit beyond it and calculate
// the Delaunay triangle for that point, this should be the triangle
// which shares the side.
for (M=1; M<=3; M++) {
if (!s[M-1]) {
// get the two points that make up this side
if (M == 1) {
P1 = Pa;
P2 = Na;
P3 = Ma;
} else if (M == 2) {
P1 = Pa;
P2 = Ma;
P3 = Na;
} else if (M == 3) {
P1 = Na;
P2 = Ma;
P3 = Pa;
}
// get the coordinates of the centre of this side
xm = (fXN[P1]+fXN[P2])/2.;
ym = (fYN[P1]+fYN[P2])/2.;
// we want to add a little to these coordinates to get a point just
// outside the triangle; (sx,sy) will be the vector that represents
// the side
sx = fXN[P1]-fXN[P2];
sy = fYN[P1]-fYN[P2];
// (nx,ny) will be the normal to the side, but don't know if it's
// pointing in or out yet
nx = sy;
ny = -sx;
nn = TMath::Sqrt(nx*nx+ny*ny);
nx = nx/nn;
ny = ny/nn;
mx = fXN[P3]-xm;
my = fYN[P3]-ym;
mdotn = mx*nx+my*ny;
if (mdotn > 0) {
// (nx,ny) is pointing in, we want it pointing out
nx = -nx;
ny = -ny;
}
// increase/decrease xm and ym a little to produce a point
// just outside the triangle (ensuring that the amount added will
// be large enough such that it won't be lost in rounding errors)
A = TMath::Abs(TMath::Max(alittlebit*xm,alittlebit*ym));
xx = xm+nx*A;
yy = ym+ny*A;
// try and find a new Delaunay triangle for this point
ComputeZ(xx,yy);
// this side of T1 should now, hopefully, if it's not part of the
// hull, be shared with a new Delaunay triangle just calculated by ComputeZ
}
}
T1++;
}
}
//______________________________________________________________________________
void TGraph2D::FindHull()
{
// Finds those points which make up the convex hull of the set. If the xy
// plane were a sheet of wood, and the points were nails hammered into it
// at the respective coordinates, then if an elastic band were stretched
// over all the nails it would form the shape of the convex hull. Those
// nails in contact with it are the points that make up the hull.
Int_t N,nhull_tmp;
Bool_t in;
if (!fHullPoints) fHullPoints = new Int_t[fNpoints];
nhull_tmp = 0;
for(N=1; N<=fNpoints; N++) {
// if the point is not inside the hull of the set of all points
// bar it, then it is part of the hull of the set of all points
// including it
in = InHull(N,N);
if (!in) {
// cannot increment fNhull directly - InHull needs to know that
// the hull has not yet been completely found
nhull_tmp++;
fHullPoints[nhull_tmp-1] = N;
}
}
fNhull = nhull_tmp;
}
//______________________________________________________________________________
Bool_t TGraph2D::InHull(Int_t E, Int_t X) const
{
// Is point E inside the hull defined by all points apart from X ?
Int_t n1,n2,N,M,Ntry;
Double_t lastdphi,dd1,dd2,dx1,dx2,dx3,dy1,dy2,dy3;
Double_t U,V,vNv1,vNv2,phi1,phi2,dphi,xx,yy;
Bool_t DTinhull = kFALSE;
xx = fXN[E];
yy = fYN[E];
if (fNhull > 0) {
// The hull has been found - no need to use any points other than
// those that make up the hull
Ntry = fNhull;
} else {
// The hull has not yet been found, will have to try every point
Ntry = fNpoints;
}
// N1 and N2 will represent the two points most separated by angle
// from point E. Initially the angle between them will be <180 degs.
// But subsequent points will increase the N1-E-N2 angle. If it
// increases above 180 degrees then point E must be surrounded by
// points - it is not part of the hull.
n1 = 1;
n2 = 2;
if (n1 == X) {
n1 = n2;
n2++;
} else if (n2 == X) {
n2++;
}
// Get the angle N1-E-N2 and set it to lastdphi
dx1 = xx-fXN[n1];
dy1 = yy-fYN[n1];
dx2 = xx-fXN[n2];
dy2 = yy-fYN[n2];
phi1 = TMath::ATan2(dy1,dx1);
phi2 = TMath::ATan2(dy2,dx2);
dphi = (phi1-phi2)-((Int_t)((phi1-phi2)/TMath::TwoPi())*TMath::TwoPi());
if (dphi < 0) dphi = dphi+TMath::TwoPi();
lastdphi = dphi;
for (N=1; N<=Ntry; N++) {
if (fNhull > 0) {
// Try hull point N
M = fHullPoints[N-1];
} else {
M = N;
}
if ((M!=n1) && (M!=n2) && (M!=X)) {
// Can the vector E->M be represented as a sum with positive
// coefficients of vectors E->N1 and E->N2?
dx1 = xx-fXN[n1];
dy1 = yy-fYN[n1];
dx2 = xx-fXN[n2];
dy2 = yy-fYN[n2];
dx3 = xx-fXN[M];
dy3 = yy-fYN[M];
dd1 = (dx2*dy1-dx1*dy2);
dd2 = (dx1*dy2-dx2*dy1);
if (dd1*dd2!=0) {
U = (dx2*dy3-dx3*dy2)/dd1;
V = (dx1*dy3-dx3*dy1)/dd2;
if ((U<0) || (V<0)) {
// No, it cannot - point M does not lie inbetween N1 and N2 as
// viewed from E. Replace either N1 or N2 to increase the
// N1-E-N2 angle. The one to replace is the one which makes the
// smallest angle with E->M
vNv1 = (dx1*dx3+dy1*dy3)/TMath::Sqrt(dx1*dx1+dy1*dy1);
vNv2 = (dx2*dx3+dy2*dy3)/TMath::Sqrt(dx2*dx2+dy2*dy2);
if (vNv1 > vNv2) {
n1 = M;
phi1 = TMath::ATan2(dy3,dx3);
phi2 = TMath::ATan2(dy2,dx2);
} else {
n2 = M;
phi1 = TMath::ATan2(dy1,dx1);
phi2 = TMath::ATan2(dy3,dx3);
}
dphi = (phi1-phi2)-((Int_t)((phi1-phi2)/TMath::TwoPi())*TMath::TwoPi());
if (dphi < 0) dphi = dphi+TMath::TwoPi();
if (((dphi-TMath::Pi())*(lastdphi-TMath::Pi())) < 0) {
// The addition of point M means the angle N1-E-N2 has risen
// above 180 degs, the point is in the hull.
goto L10;
}
lastdphi = dphi;
}
}
}
}
// Point E is not surrounded by points - it is not in the hull.
goto L999;
L10:
DTinhull = kTRUE;
L999:
return DTinhull;
}
//______________________________________________________________________________
Double_t TGraph2D::Interpolate(Double_t x, Double_t y) const
{
// Finds the z value at the position (x,y) thanks to
// the Delaunay interpolation.
// If needed, offset fX and fY so they average zero, and scale so the average
// of the X and Y ranges is one. The normalized version of fX and fY used
// in ComputeZ.
if (!fXN) {
Double_t xmax = GetXmax();
Double_t ymax = GetYmax();
Double_t xmin = GetXmin();
Double_t ymin = GetYmin();
((TGraph2D*)this)->fXoffset = -(xmax+xmin)/2.;
((TGraph2D*)this)->fYoffset = -(ymax+ymin)/2.;
((TGraph2D*)this)->fScaleFactor = 2./((xmax-xmin)+(ymax-ymin));
((TGraph2D*)this)->fXNmax = (xmax+fXoffset)*fScaleFactor;
((TGraph2D*)this)->fXNmin = (xmin+fXoffset)*fScaleFactor;
((TGraph2D*)this)->fYNmax = (ymax+fYoffset)*fScaleFactor;
((TGraph2D*)this)->fYNmin = (ymin+fYoffset)*fScaleFactor;
((TGraph2D*)this)->fXN = new Double_t[fNpoints+1];
((TGraph2D*)this)->fYN = new Double_t[fNpoints+1];
for (Int_t N=0; N<fNpoints; N++) {
fXN[N+1] = (fX[N]+fXoffset)*fScaleFactor;
fYN[N+1] = (fY[N]+fYoffset)*fScaleFactor;
}
}
// If needed, creates the arrays to hold the Delaunay triangles.
// A maximum number of 2*fNpoints is guessed. If more triangles will be
// find, FillIt will automatically enlarge these arrays.
if (!fPTried) {
((TGraph2D*)this)->fTriedSize = 2*fNpoints;
((TGraph2D*)this)->fPTried = new Int_t[fTriedSize];
((TGraph2D*)this)->fNTried = new Int_t[fTriedSize];
((TGraph2D*)this)->fMTried = new Int_t[fTriedSize];
}
Double_t sx = (x+fXoffset)*fScaleFactor;
Double_t sy = (y+fYoffset)*fScaleFactor;
return ((TGraph2D*)this)->ComputeZ(sx, sy);
}
//______________________________________________________________________________
Double_t TGraph2D::InterpolateOnPlane(Int_t TI1, Int_t TI2, Int_t TI3, Int_t E) const
{
// Finds the z-value at point E given that it lies
// on the plane defined by T1,T2,T3
Int_t tmp;
Bool_t swap;
Double_t x1,x2,x3,y1,y2,y3,f1,f2,f3,U,V,W;
Int_t T1 = TI1;
Int_t T2 = TI2;
Int_t T3 = TI3;
// order the vertices
L1:
swap = kFALSE;
if (T2 > T1) { tmp = T1; T1 = T2; T2 = tmp; swap = kTRUE;}
if (T3 > T2) { tmp = T2; T2 = T3; T3 = tmp; swap = kTRUE;}
if (swap) goto L1;
x1 = fXN[T1];
x2 = fXN[T2];
x3 = fXN[T3];
y1 = fYN[T1];
y2 = fYN[T2];
y3 = fYN[T3];
f1 = fZ[T1-1];
f2 = fZ[T2-1];
f3 = fZ[T3-1];
U = (f1*(y2-y3)+f2*(y3-y1)+f3*(y1-y2))/(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2));
V = (f1*(x2-x3)+f2*(x3-x1)+f3*(x1-x2))/(y1*(x2-x3)+y2*(x3-x1)+y3*(x1-x2));
W = f1-U*x1-V*y1;
return U*fXN[E]+V*fYN[E]+W;
}
//______________________________________________________________________________
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;
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);
grFitter->Clear();
// Get pointer to the function by searching in the list of functions in ROOT
grFitter->SetUserFunc(f2);
grFitter->SetFitOption(fitOption);
if (!f2) { Printf("Function is a null pointer"); return 0; }
npar = f2->GetNpar();
if (npar <=0) { Printf("Illegal number of parameters = %d",npar); return 0; }
// Check that function has same dimension as histogram
if (f2->GetNdim() != 2) {
Error("Fit","Function %s is not 1-D",f2->GetName());
return 0;
}
//*-*- 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.");
}
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;
}
//______________________________________________________________________________
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) const
{
// 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 ComputeZ at each bin centre to build up 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.
TString opt = option;
opt.ToLower();
if (fHistogram) {
if (!opt.Contains("empty") && fHistogram->GetEntries()==0) {
((TGraph2D*)this)->Reset(1);
} else {
return fHistogram;
}
}
Double_t x, y, z, sx, sy;
Int_t N;
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);
((TGraph2D*)this)->fHistogram = new TH2D(GetName(),GetTitle(),
fNpx ,hxmin, hxmax,
fNpy, hymin, hymax);
TH1::AddDirectory(add);
fHistogram->SetBit(TH1::kNoStats);
// Option "empty" is selected. An empty histogram is returned.
if (opt.Contains("empty")) {
if (fMinimum != -1111) {
fHistogram->SetMinimum(fMinimum);
} else {
fHistogram->SetMinimum(GetZmin());
}
if (fMaximum != -1111) {
fHistogram->SetMaximum(fMaximum);
} else {
fHistogram->SetMaximum(GetZmax());
}
return fHistogram;
}
((TGraph2D*)this)->fXoffset = -(xmax+xmin)/2.;
((TGraph2D*)this)->fYoffset = -(ymax+ymin)/2.;
((TGraph2D*)this)->fScaleFactor = 2./((xmax-xmin)+(ymax-ymin));
// If needed, offset fX and fY so they average zero, and scale so the average
// of the X and Y ranges is one. The normalized version of fX and fY used
// in ComputeZ.
if (!fXN) {
((TGraph2D*)this)->fXNmax = (xmax+fXoffset)*fScaleFactor;
((TGraph2D*)this)->fXNmin = (xmin+fXoffset)*fScaleFactor;
((TGraph2D*)this)->fYNmax = (ymax+fYoffset)*fScaleFactor;
((TGraph2D*)this)->fYNmin = (ymin+fYoffset)*fScaleFactor;
((TGraph2D*)this)->fXN = new Double_t[fNpoints+1];
((TGraph2D*)this)->fYN = new Double_t[fNpoints+1];
for (N=0; N<fNpoints; N++) {
fXN[N+1] = (fX[N]+fXoffset)*fScaleFactor;
fYN[N+1] = (fY[N]+fYoffset)*fScaleFactor;
}
}
// If needed, creates the arrays to hold the Delaunay triangles.
// A maximum number of 2*fNpoints is guessed. If more triangles will be
// find, FillIt will automatically enlarge these arrays.
if (!fPTried) {
((TGraph2D*)this)->fTriedSize = 2*fNpoints;
((TGraph2D*)this)->fPTried = new Int_t[fTriedSize];
((TGraph2D*)this)->fNTried = new Int_t[fTriedSize];
((TGraph2D*)this)->fMTried = new Int_t[fTriedSize];
}
((TGraph2D*)this)->FindHull();
Double_t dx = (hxmax-hxmin)/fNpx;
Double_t dy = (hymax-hymin)/fNpy;
for (Int_t ix=1; ix<=fNpx; ix++) {
x = hxmin+(ix-0.5)*dx;
sx = (x+fXoffset)*fScaleFactor;
for (Int_t iy=1; iy<=fNpy; iy++) {
y = hymin+(iy-0.5)*dy;
sy = (y+fYoffset)*fScaleFactor;
z = ((TGraph2D*)this)->ComputeZ(sx, sy);
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;
}
//______________________________________________________________________________
void TGraph2D::Paint(Option_t *option)
{
// Paints this 2D graph with its current attributes
TString opt = option;
opt.ToLower();
if (opt.Contains("tri") || opt.Contains("p")) {
PaintTriangles(option);
} 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(option);
}
}
//______________________________________________________________________________
void TGraph2D::PaintLevels(Int_t *T,Double_t *x, Double_t *y,
Double_t zmin, Double_t zmax, Int_t grid)
{
// Paints one triangle according to the "grid" value.
// grid = 0 : paint the color levels
// grid = 1 : paint the grid
Int_t i, FC, ncolors, theColor0, theColor2;
Int_t P0=T[0]-1;
Int_t P1=T[1]-1;
Int_t P2=T[2]-1;
Double_t xl[2],yl[2];
Double_t Zl, R21, R20, R10;
Double_t X0 = x[0] , X2 = x[0];
Double_t Y0 = y[0] , Y2 = y[0];
Double_t Z0 = fZ[P0], Z2 = fZ[P0];
// Order along Z axis the points (Xi,Yi,Zi) where "i" belongs to {0,1,2}
// After this Z0 < Z1 < Z2
Int_t I0=0, I1=0, I2=0;
if (fZ[P1]<=Z0) {Z0=fZ[P1]; X0=x[1]; Y0=y[1]; I0=1;}
if (fZ[P1]>Z2) {Z2=fZ[P1]; X2=x[1]; Y2=y[1]; I2=1;}
if (fZ[P2]<=Z0) {Z0=fZ[P2]; X0=x[2]; Y0=y[2]; I0=2;}
if (fZ[P2]>Z2) {Z2=fZ[P2]; X2=x[2]; Y2=y[2]; I2=2;}
I1 = 3-I2-I0;
Double_t X1 = x[I1];
Double_t Y1 = y[I1];
Double_t Z1 = fZ[T[I1]-1];
Double_t Zi=0, Zip=0;
switch (grid) {
case 0:
// Paint the colors levels
// Compute the color associated to Z0 (theColor0) and Z2 (theColor2)
ncolors = gStyle->GetNumberOfColors();
theColor0 = (Int_t)( ((Z0-zmin)/(zmax-zmin))*(ncolors-1) );
theColor2 = (Int_t)( ((Z2-zmin)/(zmax-zmin))*(ncolors-1) );
// The stripes drawn to fill the triangles may have up to 5 points
Double_t xp[5], yp[5];
// Rl = Ratio between Z0 and Z2 (long)
// Rs = Ratio between Z0 and Z1 or Z1 and Z2 (short)
Double_t Rl,Rs;
// Zi = Z values of the stripe number i
// Zip = Previous Zi
// Ci = Color of the stripe number i
// npf = number of point needed to draw the current stripe
Int_t Ci,npf;
FC = GetFillColor();
// If the Z0's color and Z2's colors are the same, the whole triangle
// can be painted in one go.
if(theColor0 == theColor2) {
SetFillColor(gStyle->GetColorPalette(theColor0));
TAttFill::Modify();
gPad->PaintFillArea(3,x,y);
// The triangle must be painted with several colors
} else {
for(Ci=theColor0; Ci<=theColor2; Ci++) {
SetFillColor(gStyle->GetColorPalette(Ci));
TAttFill::Modify();
if (Ci==theColor0) {
Zi = (((Ci+1)*(zmax-zmin))/(ncolors-1))+zmin;
xp[0] = X0;
yp[0] = Y0;
Rl = (Zi-Z0)/(Z2-Z0);
xp[1] = Rl*(X2-X0)+X0;
yp[1] = Rl*(Y2-Y0)+Y0;
if (Zi>=Z1 || Z0==Z1) {
Rs = (Zi-Z1)/(Z2-Z1);
xp[2] = Rs*(X2-X1)+X1;
yp[2] = Rs*(Y2-Y1)+Y1;
xp[3] = X1;
yp[3] = Y1;
npf = 4;
} else {
Rs = (Zi-Z0)/(Z1-Z0);
xp[2] = Rs*(X1-X0)+X0;
yp[2] = Rs*(Y1-Y0)+Y0;
npf = 3;
}
} else if (Ci==theColor2) {
xp[0] = xp[1];
yp[0] = yp[1];
xp[1] = X2;
yp[1] = Y2;
if (Zi<Z1 || Z2==Z1) {
xp[3] = xp[2];
yp[3] = yp[2];
xp[2] = X1;
yp[2] = Y1;
npf = 4;
} else {
npf = 3;
}
} else {
Zi = (((Ci+1)*(zmax-zmin))/(ncolors-1))+zmin;
xp[0] = xp[1];
yp[0] = yp[1];
Rl = (Zi-Z0)/(Z2-Z0);
xp[1] = Rl*(X2-X0)+X0;
yp[1] = Rl*(Y2-Y0)+Y0;
if ( Zi>=Z1 && Zip<=Z1) {
xp[3] = X1;
yp[3] = Y1;
xp[4] = xp[2];
yp[4] = yp[2];
npf = 5;
} else {
xp[3] = xp[2];
yp[3] = yp[2];
npf = 4;
}
if (Zi<Z1) {
Rs = (Zi-Z0)/(Z1-Z0);
xp[2] = Rs*(X1-X0)+X0;
yp[2] = Rs*(Y1-Y0)+Y0;
} else {
Rs = (Zi-Z1)/(Z2-Z1);
xp[2] = Rs*(X2-X1)+X1;
yp[2] = Rs*(Y2-Y1)+Y1;
}
}
Zip = Zi;
// Paint a stripe
gPad->PaintFillArea(npf,xp,yp);
}
}
SetFillColor(FC);
TAttFill::Modify();
break;
case 1:
// Paint the grid levels
SetLineStyle(3);
TAttLine::Modify();
for(i=0; i<fNbLevels; i++){
Zl=fGridLevels[i];
if(Zl >= Z0 && Zl <=Z2) {
R21=(Zl-Z1)/(Z2-Z1);
R20=(Zl-Z0)/(Z2-Z0);
R10=(Zl-Z0)/(Z1-Z0);
xl[0]=R20*(X2-X0)+X0;
yl[0]=R20*(Y2-Y0)+Y0;
if(Zl >= Z1 && Zl <=Z2) {
xl[1]=R21*(X2-X1)+X1;
yl[1]=R21*(Y2-Y1)+Y1;
} else {
xl[1]=R10*(X1-X0)+X0;
yl[1]=R10*(Y1-Y0)+Y0;
}
gPad->PaintPolyLine(2,xl,yl);
}
}
SetLineStyle(1);
TAttLine::Modify();
break;
default:
break;
}
}
//______________________________________________________________________________
void TGraph2D::PaintPolyMarker0(Int_t n, Double_t *x, Double_t *y)
{
// Paints a circle at each vertex. Each circle background is white.
SetMarkerSize(GetMarkerSize());
Int_t MC = GetMarkerColor();
for (Int_t i=0; i<n; i++) {
SetMarkerStyle(20);
SetMarkerColor(0);
TAttMarker::Modify();
gPad->PaintPolyMarker(1,&x[i],&y[i]);
SetMarkerStyle(24);
SetMarkerColor(MC);
TAttMarker::Modify();
gPad->PaintPolyMarker(1,&x[i],&y[i]);
}
}
//______________________________________________________________________________
void TGraph2D::PaintTriangles(Option_t *option)
{
// Paints the 2D graph triangles
Double_t x[4], y[4], temp1[3],temp2[3];
Int_t IT,T[3];
TString opt = option;
Bool_t triangles = opt.Contains("tri");
Bool_t tri1 = opt.Contains("tri1");
Bool_t tri2 = opt.Contains("tri2");
Bool_t markers = opt.Contains("p");
Bool_t markers0 = opt.Contains("p0");
Bool_t wire = opt.Contains("w");
Bool_t same = opt.Contains("s");
Bool_t backbox = opt.Contains("bb");
Bool_t frontbox = opt.Contains("fb");
Bool_t axis = opt.Contains("a");
Bool_t logx = gPad->GetLogx();
Bool_t logy = gPad->GetLogy();
Bool_t logz = gPad->GetLogz();
if (markers && !triangles) {
GetHistogram("empty");
} else {
GetHistogram();
}
if (!same) {
if (!backbox) {
fHistogram->Paint("BB");
} else {
fHistogram->Paint("bb");
}
if (!axis) fHistogram->Paint("abb");
}
fView = gPad->GetView();
if (!fView) {
Error("PaintTriangles", "No TView in current pad");
return;
}
if (!tri1 && !tri2 && !wire) DefineGridLevels();
// Compute minimums and maximums
TAxis *xaxis = fHistogram->GetXaxis();
Int_t first = xaxis->GetFirst();
Double_t xmin = xaxis->GetBinLowEdge(first);
if (logx && xmin <= 0) xmin = xaxis->GetBinUpEdge(xaxis->FindFixBin(0.01*xaxis->GetBinWidth(first)));
Double_t xmax = xaxis->GetBinUpEdge(xaxis->GetLast());
TAxis *yaxis = fHistogram->GetYaxis();
first = yaxis->GetFirst();
Double_t ymin = yaxis->GetBinLowEdge(first);
if (logy && ymin <= 0) ymin = yaxis->GetBinUpEdge(yaxis->FindFixBin(0.01*yaxis->GetBinWidth(first)));
Double_t ymax = yaxis->GetBinUpEdge(yaxis->GetLast());
Double_t zmax = fHistogram->GetMaximum();
Double_t zmin = fHistogram->GetMinimum();
if (logz && zmin <= 0) zmin = TMath::Min((Double_t)1, (Double_t)0.001*fHistogram->GetMaximum());
// For each triangle, compute the distance between the triangle centre
// and the back planes. Then these distances are sorted in order to draw
// the triangles from back to front.
if (triangles) {
FindAllTriangles();
Double_t cp = TMath::Cos(fView->GetLongitude()*TMath::Pi()/180.);
Double_t sp = TMath::Sin(fView->GetLongitude()*TMath::Pi()/180.);
if (fOrder) {delete [] fOrder; fOrder = 0;}
if (fDist) {delete [] fDist; fDist = 0;}
fOrder = new Int_t[fNdt];
fDist = new Double_t[fNdt];
Double_t xd,yd;
Int_t P, N, M;
Bool_t o = kFALSE;
for (IT=0; IT<fNdt; IT++) {
P = fPTried[IT];
N = fNTried[IT];
M = fMTried[IT];
xd = (fXN[P]+fXN[N]+fXN[M])/3;
yd = (fYN[P]+fYN[N]+fYN[M])/3;
if ((cp >= 0) && (sp >= 0.)) {
fDist[IT] = -(fXNmax-xd+fYNmax-yd);
} else if ((cp <= 0) && (sp >= 0.)) {
fDist[IT] = -(fXNmax-xd+yd-fYNmin);
o = kTRUE;
} else if ((cp <= 0) && (sp <= 0.)) {
fDist[IT] = -(xd-fXNmin+yd-fYNmin);
} else {
fDist[IT] = -(xd-fXNmin+fYNmax-yd);
o = kTRUE;
}
}
TMath::Sort(fNdt, fDist, fOrder, o);
}
// Draw markers only
if (markers && !triangles) {
Double_t *xm = new Double_t[fNpoints];
Double_t *ym = new Double_t[fNpoints];
Int_t npd = 0;
for (IT=0; IT<fNpoints; IT++) {
if(fX[IT] < xmin || fX[IT] > xmax) continue;
if(fY[IT] < ymin || fY[IT] > ymax) continue;
npd++;
temp1[0] = fX[IT];
temp1[1] = fY[IT];
temp1[2] = fZ[IT];
temp1[0] = TMath::Max(temp1[0],xmin);
temp1[1] = TMath::Max(temp1[1],ymin);
temp1[2] = TMath::Max(temp1[2],zmin);
temp1[2] = TMath::Min(temp1[2],zmax);
if (logx) temp1[0] = TMath::Log10(temp1[0]);
if (logy) temp1[1] = TMath::Log10(temp1[1]);
if (logz) temp1[2] = TMath::Log10(temp1[2]);
fView->WCtoNDC(temp1, &temp2[0]);
xm[IT] = temp2[0];
ym[IT] = temp2[1];
}
if (markers0) {
PaintPolyMarker0(npd,xm,ym);
} else {
SetMarkerStyle(GetMarkerStyle());
SetMarkerSize(GetMarkerSize());
SetMarkerColor(GetMarkerColor());
TAttMarker::Modify();
gPad->PaintPolyMarker(npd,xm,ym);
}
delete [] xm;
delete [] ym;
// Draw the triangles and markers if requested
} else if (triangles) {
SetFillColor(GetFillColor());
Int_t FS = GetFillStyle();
SetFillStyle(1001);
TAttFill::Modify();
SetLineColor(GetLineColor());
TAttLine::Modify();
Int_t LS = GetLineStyle();
for (IT=0; IT<fNdt; IT++) {
T[0] = fPTried[fOrder[IT]];
T[1] = fNTried[fOrder[IT]];
T[2] = fMTried[fOrder[IT]];
for (Int_t t=0; t<3; t++) {
if(fX[T[t]-1] < xmin || fX[T[t]-1] > xmax) goto endloop;
if(fY[T[t]-1] < ymin || fY[T[t]-1] > ymax) goto endloop;
temp1[0] = fX[T[t]-1];
temp1[1] = fY[T[t]-1];
temp1[2] = fZ[T[t]-1];
temp1[0] = TMath::Max(temp1[0],xmin);
temp1[1] = TMath::Max(temp1[1],ymin);
temp1[2] = TMath::Max(temp1[2],zmin);
temp1[2] = TMath::Min(temp1[2],zmax);
if (logx) temp1[0] = TMath::Log10(temp1[0]);
if (logy) temp1[1] = TMath::Log10(temp1[1]);
if (logz) temp1[2] = TMath::Log10(temp1[2]);
fView->WCtoNDC(temp1, &temp2[0]);
x[t] = temp2[0];
y[t] = temp2[1];
}
x[3] = x[0];
y[3] = y[0];
if (tri1 || tri2) PaintLevels(T,x,y,zmin,zmax,0);
if (!tri1 && !tri2 && !wire) {
gPad->PaintFillArea(3,x,y);
PaintLevels(T,x,y,zmin,zmax,1);
}
if (!tri2) gPad->PaintPolyLine(4,x,y);
if (markers) {
if (markers0) {
PaintPolyMarker0(3,x,y);
} else {
SetMarkerStyle(GetMarkerStyle());
SetMarkerSize(GetMarkerSize());
SetMarkerColor(GetMarkerColor());
TAttMarker::Modify();
gPad->PaintPolyMarker(3,x,y);
}
}
endloop:
continue;
}
SetFillStyle(FS);
SetLineStyle(LS);
TAttLine::Modify();
TAttFill::Modify();
delete [] fOrder; fOrder = 0;
delete [] fDist; fDist = 0;
}
if (!same && !frontbox) fHistogram->Paint("fb");
if (fGridLevels) {delete [] fGridLevels; fGridLevels = 0;}
}
//______________________________________________________________________________
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
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;
Reset(1);
return ipoint;
}
//_______________________________________________________________________
void TGraph2D::Reset(Int_t level)
{
// Called each time fHistogram should be recreated.
// level = 0 : it is enough to delete fHistogram
// level = 1 : the data set has changed, the hull and triangles
// must be recomputed.
if (fHistogram) delete fHistogram;
fHistogram = 0;
if (level == 0) return;
if (fHullPoints) delete [] fHullPoints;
if (fPTried) delete [] fPTried;
if (fNTried) delete [] fNTried;
if (fMTried) delete [] fMTried;
if (fOrder) delete [] fOrder;
if (fDist) delete [] fDist;
if (fXN) delete [] fXN;
if (fYN) delete [] fYN;
fHullPoints = 0;
fPTried = 0;
fNTried = 0;
fMTried = 0;
fOrder = 0;
fDist = 0;
fXN = 0;
fYN = 0;
fNhull = 0;
fNdt = 0;
}
//______________________________________________________________________________
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;
}
Reset();
}
//______________________________________________________________________________
void TGraph2D::SetMaximum(Double_t maximum)
{
fMaximum = maximum;
GetHistogram()->SetMaximum(maximum);
}
//______________________________________________________________________________
void TGraph2D::SetMinimum(Double_t minimum)
{
fMinimum = minimum;
GetHistogram()->SetMinimum(minimum);
}
//______________________________________________________________________________
void TGraph2D::SetMaxIter(Int_t n)
{
// Defines the number of triangles tested for a Delaunay triangle
// (number of iterations) before abandoning the search
fMaxIter = n;
}
//______________________________________________________________________________
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;
}
Reset();
}
//______________________________________________________________________________
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;
}
Reset();
}
//______________________________________________________________________________
void TGraph2D::SetMarginBinsContent(Double_t z)
{
// Sets the histogram bin height for points lying outside the convex hull ie:
// the bins in the margin.
fZout = z;
Reset();
}
//______________________________________________________________________________
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);
}
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.