diff --git a/graf2d/graf/inc/TGraph.h b/graf2d/graf/inc/TGraph.h index 76100feb073c447e4b560819266dbdc0be67e3e4..ed3ddd786d0978bec3d84fe81a15367d0e6ed08b 100644 --- a/graf2d/graf/inc/TGraph.h +++ b/graf2d/graf/inc/TGraph.h @@ -178,6 +178,7 @@ public: virtual void Sort(Bool_t (*greater)(const TGraph*, Int_t, Int_t)=&TGraph::CompareX, Bool_t ascending=kTRUE, Int_t low=0, Int_t high=-1111); virtual void UseCurrentStyle(); + void Zero(Int_t &k,Double_t AZ,Double_t BZ,Double_t E2,Double_t &X,Double_t &Y,Int_t maxiterations); ClassDef(TGraph,4) //Graph graphics class }; diff --git a/graf2d/graf/src/TGraph.cxx b/graf2d/graf/src/TGraph.cxx index 102d0097095bd3ce7eea72ce6f703d3d14087010..e1c292ce058a8e29aed5154ee6fbc2095e9b19bb 100644 --- a/graf2d/graf/src/TGraph.cxx +++ b/graf2d/graf/src/TGraph.cxx @@ -2307,3 +2307,130 @@ Int_t TGraph::Merge(TCollection* li) } return GetN(); } + + +//______________________________________________________________________________ +void TGraph::Zero(Int_t &k,Double_t AZ,Double_t BZ,Double_t E2,Double_t &X,Double_t &Y + ,Int_t maxiterations) +{ + // Find zero of a continuous function. + // This function finds a real zero of the continuous real + // function Y(X) in a given interval (A,B). See accompanying + // notes for details of the argument list and calling sequence + + static Double_t a, b, ya, ytest, y1, x1, h; + static Int_t j1, it, j3, j2; + Double_t yb, x2; + yb = 0; + + // Calculate Y(X) at X=AZ. + if (k <= 0) { + a = AZ; + b = BZ; + X = a; + j1 = 1; + it = 1; + k = j1; + return; + } + + // Test whether Y(X) is sufficiently small. + + if (TMath::Abs(Y) <= E2) { k = 2; return; } + + // Calculate Y(X) at X=BZ. + + if (j1 == 1) { + ya = Y; + X = b; + j1 = 2; + return; + } + // Test whether the signs of Y(AZ) and Y(BZ) are different. + // if not, begin the binary subdivision. + + if (j1 != 2) goto L100; + if (ya*Y < 0) goto L120; + x1 = a; + y1 = ya; + j1 = 3; + h = b - a; + j2 = 1; + x2 = a + 0.5*h; + j3 = 1; + it++; //*-*- Check whether (maxiterations) function values have been calculated. + if (it >= maxiterations) k = j1; + else X = x2; + return; + + // Test whether a bracket has been found . + // If not,continue the search + +L100: + if (j1 > 3) goto L170; + if (ya*Y >= 0) { + if (j3 >= j2) { + h = 0.5*h; j2 = 2*j2; + a = x1; ya = y1; x2 = a + 0.5*h; j3 = 1; + } + else { + a = X; ya = Y; x2 = X + h; j3++; + } + it++; + if (it >= maxiterations) k = j1; + else X = x2; + return; + } + + // The first bracket has been found.calculate the next X by the + // secant method based on the bracket. + +L120: + b = X; + yb = Y; + j1 = 4; +L130: + if (TMath::Abs(ya) > TMath::Abs(yb)) { x1 = a; y1 = ya; X = b; Y = yb; } + else { x1 = b; y1 = yb; X = a; Y = ya; } + + // Use the secant method based on the function values y1 and Y. + // check that x2 is inside the interval (a,b). + +L150: + x2 = X-Y*(X-x1)/(Y-y1); + x1 = X; + y1 = Y; + ytest = 0.5*TMath::Min(TMath::Abs(ya),TMath::Abs(yb)); + if ((x2-a)*(x2-b) < 0) { + it++; + if (it >= maxiterations) k = j1; + else X = x2; + return; + } + + // Calculate the next value of X by bisection . Check whether + // the maximum accuracy has been achieved. + +L160: + x2 = 0.5*(a+b); + ytest = 0; + if ((x2-a)*(x2-b) >= 0) { k = 2; return; } + it++; + if (it >= maxiterations) k = j1; + else X = x2; + return; + + + // Revise the bracket (a,b). + +L170: + if (j1 != 4) return; + if (ya*Y < 0) { b = X; yb = Y; } + else { a = X; ya = Y; } + + // Use ytest to decide the method for the next value of X. + + if (ytest <= 0) goto L130; + if (TMath::Abs(Y)-ytest <= 0) goto L150; + goto L160; +} diff --git a/hist/histpainter/inc/TGraphPainter.h b/hist/histpainter/inc/TGraphPainter.h index 3a05136c37bfb4553d01f8453f39beada0ab4ac4..520dd40749ad39862a479c1d7539d6c8493f0659 100644 --- a/hist/histpainter/inc/TGraphPainter.h +++ b/hist/histpainter/inc/TGraphPainter.h @@ -50,7 +50,6 @@ public: void PaintPolyLineHatches(TGraph *theGraph, Int_t n, const Double_t *x, const Double_t *y); void PaintStats(TGraph *theGraph, TF1 *fit); void Smooth(TGraph *theGraph, Int_t npoints, Double_t *x, Double_t *y, Int_t drawtype); - void Zero(Int_t &k,Double_t AZ,Double_t BZ,Double_t E2,Double_t &X,Double_t &Y,Int_t maxiterations); ClassDef(TGraphPainter,0) // TGraph painter }; diff --git a/hist/histpainter/src/TGraphPainter.cxx b/hist/histpainter/src/TGraphPainter.cxx index e2f26a06ba5ec75226ab2d400001fb3a20781ddd..54fa707efa9791f085b0578c658faeb647cdbbbb 100644 --- a/hist/histpainter/src/TGraphPainter.cxx +++ b/hist/histpainter/src/TGraphPainter.cxx @@ -3042,7 +3042,7 @@ L240: c = z; sb = s; L250: - Zero(kp,0,sb,err,s,z,maxiterations); + theGraph->Zero(kp,0,sb,err,s,z,maxiterations); if (kp == 2) goto L210; if (kp > 2) { Error("Smooth", "Attempt to plot outside plot limits"); @@ -3145,131 +3145,3 @@ L390: delete [] qly; } - -//______________________________________________________________________________ -void TGraphPainter::Zero(Int_t &k,Double_t AZ,Double_t BZ,Double_t E2,Double_t &X,Double_t &Y - ,Int_t maxiterations) -{ - // Find zero of a continuous function. - // - // Underlaying routine for PaintGraph - // This function finds a real zero of the continuous real - // function Y(X) in a given interval (A,B). See accompanying - // notes for details of the argument list and calling sequence - - static Double_t a, b, ya, ytest, y1, x1, h; - static Int_t j1, it, j3, j2; - Double_t yb, x2; - yb = 0; - - // Calculate Y(X) at X=AZ. - if (k <= 0) { - a = AZ; - b = BZ; - X = a; - j1 = 1; - it = 1; - k = j1; - return; - } - - // Test whether Y(X) is sufficiently small. - - if (TMath::Abs(Y) <= E2) { k = 2; return; } - - // Calculate Y(X) at X=BZ. - - if (j1 == 1) { - ya = Y; - X = b; - j1 = 2; - return; - } - // Test whether the signs of Y(AZ) and Y(BZ) are different. - // if not, begin the binary subdivision. - - if (j1 != 2) goto L100; - if (ya*Y < 0) goto L120; - x1 = a; - y1 = ya; - j1 = 3; - h = b - a; - j2 = 1; - x2 = a + 0.5*h; - j3 = 1; - it++; //*-*- Check whether (maxiterations) function values have been calculated. - if (it >= maxiterations) k = j1; - else X = x2; - return; - - // Test whether a bracket has been found . - // If not,continue the search - -L100: - if (j1 > 3) goto L170; - if (ya*Y >= 0) { - if (j3 >= j2) { - h = 0.5*h; j2 = 2*j2; - a = x1; ya = y1; x2 = a + 0.5*h; j3 = 1; - } - else { - a = X; ya = Y; x2 = X + h; j3++; - } - it++; - if (it >= maxiterations) k = j1; - else X = x2; - return; - } - - // The first bracket has been found.calculate the next X by the - // secant method based on the bracket. - -L120: - b = X; - yb = Y; - j1 = 4; -L130: - if (TMath::Abs(ya) > TMath::Abs(yb)) { x1 = a; y1 = ya; X = b; Y = yb; } - else { x1 = b; y1 = yb; X = a; Y = ya; } - - // Use the secant method based on the function values y1 and Y. - // check that x2 is inside the interval (a,b). - -L150: - x2 = X-Y*(X-x1)/(Y-y1); - x1 = X; - y1 = Y; - ytest = 0.5*TMath::Min(TMath::Abs(ya),TMath::Abs(yb)); - if ((x2-a)*(x2-b) < 0) { - it++; - if (it >= maxiterations) k = j1; - else X = x2; - return; - } - - // Calculate the next value of X by bisection . Check whether - // the maximum accuracy has been achieved. - -L160: - x2 = 0.5*(a+b); - ytest = 0; - if ((x2-a)*(x2-b) >= 0) { k = 2; return; } - it++; - if (it >= maxiterations) k = j1; - else X = x2; - return; - - - // Revise the bracket (a,b). - -L170: - if (j1 != 4) return; - if (ya*Y < 0) { b = X; yb = Y; } - else { a = X; ya = Y; } - - // Use ytest to decide the method for the next value of X. - - if (ytest <= 0) goto L130; - if (TMath::Abs(Y)-ytest <= 0) goto L150; - goto L160; -}