// @(#)root/graf:$Name:  $:$Id: TGraph.cxx,v 1.101 2003/05/05 16:38:34 brun Exp $
// Author: Rene Brun, Olivier Couet   12/12/94

/*************************************************************************
 * 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 <string.h>

#include "Riostream.h"
#include "TROOT.h"
#include "TGraph.h"
#include "TGaxis.h"
#include "TH1.h"
#include "TF1.h"
#include "TStyle.h"
#include "TMath.h"
#include "TVector.h"
#include "TVectorD.h"
#include "Foption.h"
#include "TRandom.h"
#include "TPaveStats.h"
#include "TPluginManager.h"
#include "TVirtualFitter.h"
#include "TVirtualPad.h"
#include "TVirtualUtilPad.h"
#include "TVirtualHistPainter.h"


Double_t *gxwork, *gywork, *gxworkl, *gyworkl;

extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);


ClassImp(TGraph)

//______________________________________________________________________________
//
//   A Graph is a graphics object made of two arrays X and Y
//   with npoints each.
//   This class supports essentially two graph categories:
//     - General case with non equidistant points
//     - Special case with equidistant points
//   The various format options to draw a Graph are explained in
//     TGraph::PaintGraph  and TGraph::PaintGrapHist
//   These two functions are derived from the HIGZ routines IGRAPH and IGHIST
//   and many modifications.
//
//  The picture below has been generated by the following macro:
//------------------------------------------------------------------
//{
//   TCanvas *c1 = new TCanvas("c1","A Simple Graph Example",200,10,700,500);
//   Double_t x[100], y[100];
//   Int_t n = 20;
//   for (Int_t i=0;i<n;i++) {
//     x[i] = i*0.1;
//     y[i] = 10*sin(x[i]+0.2);
//   }
//   gr = new TGraph(n,x,y);
//   gr->Draw("AC*");
//}
//
/*

*/
//
//

//______________________________________________________________________________
 TGraph::TGraph(): TNamed(), TAttLine(), TAttFill(1,1001), TAttMarker()
{
//*-*-*-*-*-*-*-*-*-*-*Graph default constructor-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =========================
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   fFunctions = 0;
   fHistogram = 0;
   fMaximum = -1111;
   fMinimum = -1111;
   SetBit(kClipFrame);
}

//______________________________________________________________________________
 TGraph::TGraph(Int_t n)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
// constructor with only the number of points set
// the arrsys x and y will be set later

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   if (n <= 0) {
      Error("TGraph", "illegal number of points (%d)", n);
      return;
   }

   fFunctions = new TList;
   fNpoints   = n;
   fX         = new Double_t[n];
   fY         = new Double_t[n];
   fMaximum   = -1111;
   fMinimum   = -1111;
   for (Int_t i=0;i<n;i++) {
      fX[i] = 0;
      fY[i] = 0;
   }
   SetBit(kClipFrame);
}

//______________________________________________________________________________
 TGraph::TGraph(Int_t n, const Int_t *x, const Int_t *y)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
//*-*-*-*-*-*-*-*-*-*-*Graph normal constructor with ints-*-*-*-*-*-*-*-*-*
//*-*                  ==================================

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   if (n <= 0) {
      Error("TGraph", "illegal number of points (%d)", n);
      return;
   }

   fFunctions = new TList;
   fNpoints   = n;
   fX         = new Double_t[n];
   fY         = new Double_t[n];
   fMaximum   = -1111;
   fMinimum   = -1111;
   if (!x || !y) return;
   for (Int_t i=0;i<n;i++) {
      fX[i] = (Double_t)x[i];
      fY[i] = (Double_t)y[i];
   }
   SetBit(kClipFrame);
}

//______________________________________________________________________________
 TGraph::TGraph(Int_t n, const Float_t *x, const Float_t *y)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
//*-*-*-*-*-*-*-*-*-*-*Graph normal constructor with floats-*-*-*-*-*-*-*-*-*
//*-*                  ====================================

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   if (n <= 0) {
      Error("TGraph", "illegal number of points (%d)", n);
      return;
   }

   fFunctions = new TList;
   fNpoints   = n;
   fX         = new Double_t[n];
   fY         = new Double_t[n];
   fMaximum   = -1111;
   fMinimum   = -1111;
   if (!x || !y) return;
   for (Int_t i=0;i<n;i++) {
      fX[i] = x[i];
      fY[i] = y[i];
   }
   SetBit(kClipFrame);
}

//______________________________________________________________________________
 TGraph::TGraph(Int_t n, const Double_t *x, const Double_t *y)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
//*-*-*-*-*-*-*-*-*-*-*Graph normal constructor with doubles-*-*-*-*-*-*-*-*
//*-*                  ========================
//

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   if (n <= 0) {
      Error("TGraph", "illegal number of points (%d)", n);
      return;
   }

   fFunctions = new TList;
   fNpoints   = n;
   fX         = new Double_t[n];
   fY         = new Double_t[n];
   fMaximum   = -1111;
   fMinimum   = -1111;
   SetBit(kClipFrame);
   if (!x || !y) return;
   for (Int_t i=0;i<n;i++) {
      fX[i] = x[i];
      fY[i] = y[i];
   }
}

//______________________________________________________________________________
 TGraph::TGraph(const TVector &vx, const TVector &vy)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
// Graph constructor with two vectors of floats in input
// A graph is build with the X coordinates taken from vx and Y coord from vy
// The number of points in the graph is the minimum of number of points
// in vx and vy.

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   Int_t n  = vx.GetNrows();
   Int_t ny = vy.GetNrows();
   if (ny < n) n = ny;
   if (n <= 0) {
      Error("TGraph", "illegal number of points (%d)", n);
      return;
   }

   fFunctions = new TList;
   fNpoints   = n;
   fX         = new Double_t[n];
   fY         = new Double_t[n];
   fMaximum   = -1111;
   fMinimum   = -1111;
   SetBit(kClipFrame);
   for (Int_t i=0;i<n;i++) {
      fX[i] = vx(i);
      fY[i] = vy(i);
   }
}

//______________________________________________________________________________
 TGraph::TGraph(const TVectorD &vx, const TVectorD &vy)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
// Graph constructor with two vectors of doubles in input
// A graph is build with the X coordinates taken from vx and Y coord from vy
// The number of points in the graph is the minimum of number of points
// in vx and vy.

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   Int_t n  = vx.GetNrows();
   Int_t ny = vy.GetNrows();
   if (ny < n) n = ny;
   if (n <= 0) {
      Error("TGraph", "illegal number of points (%d)", n);
      return;
   }

   fFunctions = new TList;
   fNpoints   = n;
   fX         = new Double_t[n];
   fY         = new Double_t[n];
   fMaximum   = -1111;
   fMinimum   = -1111;
   SetBit(kClipFrame);
   for (Int_t i=0;i<n;i++) {
      fX[i] = vx(i);
      fY[i] = vy(i);
   }
}

//______________________________________________________________________________
 TGraph::TGraph(const TH1 *h)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
// Graph constructor importing its parameters from the TH1 object passed as argument

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   if (!h) {
      Error("TGraph", "Pointer to histogram is null");
      return;
   }
   if (h->GetDimension() != 1) {
      Error("TGraph", "Histogram must be 1-D; h %s is %d-D",h->GetName(),h->GetDimension());
      return;
   }
   TAxis *xaxis = ((TH1*)h)->GetXaxis();
   fFunctions = new TList;
   fNpoints   = xaxis->GetNbins();
   fX         = new Double_t[fNpoints];
   fY         = new Double_t[fNpoints];
   fMaximum   = -1111;
   fMinimum   = -1111;
   SetBit(kClipFrame);
   for (Int_t i=0;i<fNpoints;i++) {
      fX[i] = xaxis->GetBinCenter(i+1);
      fY[i] = h->GetBinContent(i+1);
   }
   SetLineColor(h->GetLineColor());;
   SetLineWidth(h->GetLineWidth());
   SetLineStyle(h->GetLineStyle());
   SetFillColor(h->GetFillColor());
   SetFillStyle(h->GetFillStyle());
   SetMarkerStyle(h->GetMarkerStyle());
   SetMarkerColor(h->GetMarkerColor());
   SetMarkerSize(h->GetMarkerSize());

   SetName(h->GetName());
   SetTitle(h->GetTitle());
}


//______________________________________________________________________________
 TGraph::TGraph(const TF1 *f, Option_t *option)
       : TNamed("Graph","Graph"), TAttLine(), TAttFill(1,1001), TAttMarker()
{
// Graph constructor importing its parameters from the TF1 object passed as argument
// if option =="" (default), a TGraph is created with points computed
//                at the fNpx points of f.
// if option =="d", a TGraph is created with points computed with the derivatives
//                at the fNpx points of f.
// if option =="i", a TGraph is created with points computed with the integral
//                at the fNpx points of f.
// if option =="I", a TGraph is created with points computed with the integral
//                at the fNpx+1 points of f and the integral is normalized to 1.

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   if (!f) {
      Error("TGraph", "Pointer to function is null");
      return;
   }
   char coption = ' ';
   if (option) coption = *option;
   fFunctions = new TList;
   fNpoints   = f->GetNpx();
   Double_t xmin = f->GetXmin();
   Double_t xmax = f->GetXmax();
   Double_t dx   = (xmax-xmin)/fNpoints;
   if (coption == 'i' || coption == 'I') fNpoints++;
   fX         = new Double_t[fNpoints];
   fY         = new Double_t[fNpoints];
   fMaximum   = -1111;
   fMinimum   = -1111;
   SetBit(kClipFrame);
   Double_t integ = 0;
   Int_t i;
   for (i=0;i<fNpoints;i++) {
      if (coption == 'i' || coption == 'I') {
         fX[i] = xmin +i*dx;
         if (i == 0) fY[i] = 0;
         else        fY[i] = integ + ((TF1*)f)->Integral(fX[i]-dx,fX[i]);
         integ = fY[i];
      } else if (coption == 'd' || coption == 'D') {
         fX[i] = xmin + (i+0.5)*dx;
         fY[i] = ((TF1*)f)->Derivative(fX[i]);
      } else {
         fX[i] = xmin + (i+0.5)*dx;
         fY[i] = ((TF1*)f)->Eval(fX[i]);
      }
   }
   if (integ != 0 && coption == 'I') {
      for (i=1;i<fNpoints;i++) fY[i] /= integ;
   }

   SetLineColor(f->GetLineColor());;
   SetLineWidth(f->GetLineWidth());
   SetLineStyle(f->GetLineStyle());
   SetFillColor(f->GetFillColor());
   SetFillStyle(f->GetFillStyle());

   SetName(f->GetName());
   SetTitle(f->GetTitle());
}

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

   fFunctions = 0;
   fHistogram = 0;
   fNpoints   = 0;
   fX         = 0;
   fY         = 0;
   fMaximum   = -1111;
   fMinimum   = -1111;
   fFunctions = new TList;
   SetBit(kClipFrame);


   Double_t x,y;
   FILE *fp = fopen(filename,"r");
   if (!fp) {
      MakeZombie();
      Error("TGraph", "Cannot open file: %s, TGraph is Zombie",filename);
      return;
   }
   Set(100);  //initial number of points
   char line[80];
   Int_t np = 0;
   while (fgets(line,80,fp)) {
      sscanf(&line[0],format,&x, &y);
      if (np >= fNpoints) Set(2*fNpoints);
      SetPoint(np,x,y);
      np++;
   }
   Set(np);

   fclose(fp);
}

//______________________________________________________________________________
 TGraph::~TGraph()
{
//*-*-*-*-*-*-*-*-*-*-*Graph default destructor-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =======================

   delete [] fX;
   delete [] fY;
   if (fFunctions) {
      fFunctions->SetBit(kInvalidObject);
      fFunctions->Delete();
      delete fFunctions;
   }
   fFunctions = 0;
   delete fHistogram;
   fHistogram = 0;
}

//______________________________________________________________________________
 void TGraph::Apply(TF1 *f)
{
  // apply function f to all the data points
  // f may be a 1-D function TF1 or 2-d function TF2
  // The Y values of the graph are replaced by the new values computed
  // using the function

  for (Int_t i=0;i<fNpoints;i++) {
    fY[i] = f->Eval(fX[i],fY[i]);
  }
}


//______________________________________________________________________________
 void TGraph::Browse(TBrowser *)
{
    Draw("alp");
    gPad->Update();
}

//______________________________________________________________________________
 Bool_t TGraph::CompareX(const TGraph* gr, Int_t left, Int_t right) {
// return kTRUE if fX[left] > fX[right]. Can be used by Sort.
   return gr->fX[left]>gr->fX[right];
}

//______________________________________________________________________________
 Bool_t TGraph::CompareY(const TGraph* gr, Int_t left, Int_t right) {
// return kTRUE if fY[left] > fY[right]. Can be used by Sort.
   return gr->fY[left]>gr->fY[right];
}

//______________________________________________________________________________
 Bool_t TGraph::CompareRadius(const TGraph* gr, Int_t left, Int_t right) {
// return kTRUE if point number "left"'s distance to origin is bigger than 
// that of point number "right". Can be used by Sort.
   return gr->fX[left]*gr->fX[left]+gr->fY[left]*gr->fY[left]
      >gr->fX[right]*gr->fX[right]+gr->fY[right]*gr->fY[right];
}

//______________________________________________________________________________
 void TGraph::ComputeRange(Double_t &, Double_t &, Double_t &, Double_t &) const
{
// this function is dummy in TGraph, but redefined by TGraphErrors
}

//______________________________________________________________________________
 void TGraph::Draw(Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this graph with its current attributes*-*-*-*-*-*-*
//*-*                  ==========================================
//
//   Options to draw a graph are described in TGraph::PaintGraph

   TString opt = option;
   opt.ToLower();

   if (opt.Contains("same"))
      Error("Draw", "option "same" not supported,n"
            "see TGraph::PaintGraph() for options");

   // in case of option *, set marker style to 3 (star) and replace
   // * option by option P.
   Ssiz_t pos;
   if ((pos = opt.Index("*")) != kNPOS) {
      SetMarkerStyle(3);
      opt.Replace(pos, 1, "P");
   }
   if (gPad) {
      if (!gPad->IsEditable()) (gROOT->GetMakeDefCanvas())();
   }
   AppendPad(opt);
}

//______________________________________________________________________________
 Int_t TGraph::DistancetoPrimitive(Int_t px, Int_t py)
{
//*-*-*-*-*-*-*-*-*-*-*Compute distance from point px,py to a graph*-*-*-*-*-*
//*-*                  ===========================================
//  Compute the closest distance of approach from point px,py to this line.
//  The distance is computed in pixels units.
//

//*-*- Are we on the axis?
   Int_t distance;
   if (fHistogram) {
      distance = fHistogram->DistancetoPrimitive(px,py);
      if (distance <= 5) return distance;
   }

//*-*- Somewhere on the graph points?
   const Int_t big = 9999;
   const Int_t kMaxDiff = 10;
   Int_t puxmin = gPad->XtoAbsPixel(gPad->GetUxmin());
   Int_t puymin = gPad->YtoAbsPixel(gPad->GetUymin());
   Int_t puxmax = gPad->XtoAbsPixel(gPad->GetUxmax());
   Int_t puymax = gPad->YtoAbsPixel(gPad->GetUymax());

//*-*- return if point is not in the histogram area
   if (px <= puxmin) return big;
   if (py >= puymin) return big;
   if (px >= puxmax) return big;
   if (py <= puymax) return big;

//*-*- check if point is near one of the graph points
   Int_t i, pxp, pyp, d;
   distance = big;

   for (i=0;i<fNpoints;i++) {
      pxp = gPad->XtoAbsPixel(gPad->XtoPad(fX[i]));
      pyp = gPad->YtoAbsPixel(gPad->YtoPad(fY[i]));
      d   = TMath::Abs(pxp-px) + TMath::Abs(pyp-py);
      if (d < distance) distance = d;
   }
   if (distance < kMaxDiff) return distance;

   for (i=0;i<fNpoints-1;i++) {
      d = DistancetoLine(px, py, gPad->XtoPad(fX[i]), gPad->YtoPad(fY[i]), gPad->XtoPad(fX[i+1]), gPad->YtoPad(fY[i+1]));
      if (d < distance) distance = d;
   }

//*-*- Loop on the list of associated functions and user objects
   TObject *f;
   TIter   next(fFunctions);
   while ((f = (TObject*) next())) {
      Int_t dist = f->DistancetoPrimitive(px,py);
      if (dist < kMaxDiff) {
         gPad->SetSelected(f);
         return 0; //must be o and not dist in case of TMultiGraph
      }
   }
   return distance;
}


//______________________________________________________________________________
 void TGraph::DrawGraph(Int_t n, const Int_t *x, const Int_t *y, Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this graph with new attributes*-*-*-*-*-*-*-*-*-*
//*-*                  ===================================

   TGraph *newgraph = new TGraph(n, x, y);
   TAttLine::Copy(*newgraph);
   TAttFill::Copy(*newgraph);
   TAttMarker::Copy(*newgraph);
   newgraph->SetBit(kCanDelete);
   newgraph->AppendPad(option);
}

//______________________________________________________________________________
 void TGraph::DrawGraph(Int_t n, const Float_t *x, const Float_t *y, Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this graph with new attributes*-*-*-*-*-*-*-*-*-*
//*-*                  ===================================

   TGraph *newgraph = new TGraph(n, x, y);
   TAttLine::Copy(*newgraph);
   TAttFill::Copy(*newgraph);
   TAttMarker::Copy(*newgraph);
   newgraph->SetBit(kCanDelete);
   newgraph->AppendPad(option);
}


//______________________________________________________________________________
 void TGraph::DrawGraph(Int_t n, const Double_t *x, const Double_t *y, Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this graph with new attributes*-*-*-*-*-*-*-*-*-*
//*-*                  ===================================

   TGraph *newgraph = new TGraph(n, x, y);
   TAttLine::Copy(*newgraph);
   TAttFill::Copy(*newgraph);
   TAttMarker::Copy(*newgraph);
   newgraph->SetBit(kCanDelete);
   newgraph->AppendPad(option);
}

//______________________________________________________________________________
 void TGraph::DrawPanel()
{
//*-*-*-*-*-*-*Display a panel with all graph drawing options*-*-*-*-*-*
//*-*          ==============================================
//

   printf("TGraph::DrawPanel: not yet implementedn");
}

//______________________________________________________________________________
 void TGraph::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
//*-*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*
//*-*                  =========================================
//  This member function is called when a graph is clicked with the locator
//
//  If Left button clicked on one of the line end points, this point
//     follows the cursor until button is released.
//
//  if Middle button clicked, the line is moved parallel to itself
//     until the button is released.
//
   Int_t i, d;
   Double_t xmin, xmax, ymin, ymax, dx, dy, dxr, dyr;
   const Int_t kMaxDiff = 10;
   static Bool_t MIDDLE, BADCASE;
   static Int_t ipoint, pxp, pyp;
   static Int_t px1,px2,py1,py2;
   static Int_t pxold, pyold, px1old, py1old, px2old, py2old;
   static Int_t dpx, dpy;
   static Int_t *x=0, *y=0;

   if (!IsEditable()) {gPad->SetCursor(kHand); return;}
   if (!gPad->IsEditable()) return;

   switch (event) {

   case kButton1Down:
      BADCASE = kFALSE;
      gVirtualX->SetLineColor(-1);
      TAttLine::Modify();  //Change line attributes only if necessary
      px1 = gPad->XtoAbsPixel(gPad->GetX1());
      py1 = gPad->YtoAbsPixel(gPad->GetY1());
      px2 = gPad->XtoAbsPixel(gPad->GetX2());
      py2 = gPad->YtoAbsPixel(gPad->GetY2());
      ipoint = -1;


      if (x || y) break;
      x = new Int_t[fNpoints+1];
      y = new Int_t[fNpoints+1];
      for (i=0;i<fNpoints;i++) {
         pxp = gPad->XtoAbsPixel(gPad->XtoPad(fX[i]));
         pyp = gPad->YtoAbsPixel(gPad->YtoPad(fY[i]));
         if (pxp < -kMaxPixel || pxp >= kMaxPixel ||
             pyp < -kMaxPixel || pyp >= kMaxPixel) {
            BADCASE = kTRUE;
            continue;
         }
         gVirtualX->DrawLine(pxp-4, pyp-4, pxp+4,  pyp-4);
         gVirtualX->DrawLine(pxp+4, pyp-4, pxp+4,  pyp+4);
         gVirtualX->DrawLine(pxp+4, pyp+4, pxp-4,  pyp+4);
         gVirtualX->DrawLine(pxp-4, pyp+4, pxp-4,  pyp-4);
         x[i] = pxp;
         y[i] = pyp;
         d   = TMath::Abs(pxp-px) + TMath::Abs(pyp-py);
         if (d < kMaxDiff) ipoint =i;
      }
      dpx = 0;
      dpy = 0;
      pxold = px;
      pyold = py;
      if (ipoint < 0) return;
      if (ipoint == 0) {
         px1old = 0;
         py1old = 0;
         px2old = gPad->XtoAbsPixel(fX[1]);
         py2old = gPad->YtoAbsPixel(fY[1]);
      } else if (ipoint == fNpoints-1) {
         px1old = gPad->XtoAbsPixel(gPad->XtoPad(fX[fNpoints-2]));
         py1old = gPad->YtoAbsPixel(gPad->YtoPad(fY[fNpoints-2]));
         px2old = 0;
         py2old = 0;
      } else {
         px1old = gPad->XtoAbsPixel(gPad->XtoPad(fX[ipoint-1]));
         py1old = gPad->YtoAbsPixel(gPad->YtoPad(fY[ipoint-1]));
         px2old = gPad->XtoAbsPixel(gPad->XtoPad(fX[ipoint+1]));
         py2old = gPad->YtoAbsPixel(gPad->YtoPad(fY[ipoint+1]));
      }
      pxold = gPad->XtoAbsPixel(gPad->XtoPad(fX[ipoint]));
      pyold = gPad->YtoAbsPixel(gPad->YtoPad(fY[ipoint]));

      break;


   case kMouseMotion:

      MIDDLE = kTRUE;
      for (i=0;i<fNpoints;i++) {
         pxp = gPad->XtoAbsPixel(gPad->XtoPad(fX[i]));
         pyp = gPad->YtoAbsPixel(gPad->YtoPad(fY[i]));
         d   = TMath::Abs(pxp-px) + TMath::Abs(pyp-py);
         if (d < kMaxDiff) MIDDLE = kFALSE;
      }


//*-*- check if point is close to an axis
      if (MIDDLE) gPad->SetCursor(kMove);
      else gPad->SetCursor(kHand);
      break;

   case kButton1Motion:
      if (MIDDLE) {
         for(i=0;i<fNpoints-1;i++) {
            gVirtualX->DrawLine(x[i]+dpx, y[i]+dpy, x[i+1]+dpx, y[i+1]+dpy);
            pxp = x[i]+dpx;
            pyp = y[i]+dpy;
            if (pxp < -kMaxPixel || pxp >= kMaxPixel ||
                pyp < -kMaxPixel || pyp >= kMaxPixel) continue;
            gVirtualX->DrawLine(pxp-4, pyp-4, pxp+4,  pyp-4);
            gVirtualX->DrawLine(pxp+4, pyp-4, pxp+4,  pyp+4);
            gVirtualX->DrawLine(pxp+4, pyp+4, pxp-4,  pyp+4);
            gVirtualX->DrawLine(pxp-4, pyp+4, pxp-4,  pyp-4);
         }
         pxp = x[fNpoints-1]+dpx;
         pyp = y[fNpoints-1]+dpy;
         gVirtualX->DrawLine(pxp-4, pyp-4, pxp+4,  pyp-4);
         gVirtualX->DrawLine(pxp+4, pyp-4, pxp+4,  pyp+4);
         gVirtualX->DrawLine(pxp+4, pyp+4, pxp-4,  pyp+4);
         gVirtualX->DrawLine(pxp-4, pyp+4, pxp-4,  pyp-4);
         dpx += px - pxold;
         dpy += py - pyold;
         pxold = px;
         pyold = py;
         for(i=0;i<fNpoints-1;i++) {
            gVirtualX->DrawLine(x[i]+dpx, y[i]+dpy, x[i+1]+dpx, y[i+1]+dpy);
            pxp = x[i]+dpx;
            pyp = y[i]+dpy;
            if (pxp < -kMaxPixel || pxp >= kMaxPixel ||
                pyp < -kMaxPixel || pyp >= kMaxPixel) continue;
            gVirtualX->DrawLine(pxp-4, pyp-4, pxp+4,  pyp-4);
            gVirtualX->DrawLine(pxp+4, pyp-4, pxp+4,  pyp+4);
            gVirtualX->DrawLine(pxp+4, pyp+4, pxp-4,  pyp+4);
            gVirtualX->DrawLine(pxp-4, pyp+4, pxp-4,  pyp-4);
         }
         pxp = x[fNpoints-1]+dpx;
         pyp = y[fNpoints-1]+dpy;
         gVirtualX->DrawLine(pxp-4, pyp-4, pxp+4,  pyp-4);
         gVirtualX->DrawLine(pxp+4, pyp-4, pxp+4,  pyp+4);
         gVirtualX->DrawLine(pxp+4, pyp+4, pxp-4,  pyp+4);
         gVirtualX->DrawLine(pxp-4, pyp+4, pxp-4,  pyp-4);
      } else {
         if (px1old) gVirtualX->DrawLine(px1old, py1old, pxold,  pyold);
         if (px2old) gVirtualX->DrawLine(pxold,  pyold,  px2old, py2old);
         gVirtualX->DrawLine(pxold-4, pyold-4, pxold+4,  pyold-4);
         gVirtualX->DrawLine(pxold+4, pyold-4, pxold+4,  pyold+4);
         gVirtualX->DrawLine(pxold+4, pyold+4, pxold-4,  pyold+4);
         gVirtualX->DrawLine(pxold-4, pyold+4, pxold-4,  pyold-4);
         pxold = px;
         pxold = TMath::Max(pxold, px1);
         pxold = TMath::Min(pxold, px2);
         pyold = py;
         pyold = TMath::Max(pyold, py2);
         pyold = TMath::Min(pyold, py1);
         if (px1old) gVirtualX->DrawLine(px1old, py1old, pxold,  pyold);
         if (px2old) gVirtualX->DrawLine(pxold,  pyold,  px2old, py2old);
         gVirtualX->DrawLine(pxold-4, pyold-4, pxold+4,  pyold-4);
         gVirtualX->DrawLine(pxold+4, pyold-4, pxold+4,  pyold+4);
         gVirtualX->DrawLine(pxold+4, pyold+4, pxold-4,  pyold+4);
         gVirtualX->DrawLine(pxold-4, pyold+4, pxold-4,  pyold-4);
      }
      break;

   case kButton1Up:

//*-*- Compute x,y range
      xmin = gPad->GetUxmin();
      xmax = gPad->GetUxmax();
      ymin = gPad->GetUymin();
      ymax = gPad->GetUymax();
      dx   = xmax-xmin;
      dy   = ymax-ymin;
      dxr  = dx/(1 - gPad->GetLeftMargin() - gPad->GetRightMargin());
      dyr  = dy/(1 - gPad->GetBottomMargin() - gPad->GetTopMargin());

   // Range() could change the size of the pad pixmap and therefore should
   // be called before the other paint routines
         gPad->Range(xmin - dxr*gPad->GetLeftMargin(),
                      ymin - dyr*gPad->GetBottomMargin(),
                      xmax + dxr*gPad->GetRightMargin(),
                      ymax + dyr*gPad->GetTopMargin());
         gPad->RangeAxis(xmin, ymin, xmax, ymax);

      if (MIDDLE) {
         for(i=0;i<fNpoints;i++) {
            if (BADCASE) continue;  //do not update if big zoom and points moved
            fX[i] = gPad->PadtoX(gPad->AbsPixeltoX(x[i]+dpx));
            fY[i] = gPad->PadtoY(gPad->AbsPixeltoY(y[i]+dpy));
         }
      } else {
         fX[ipoint] = gPad->PadtoX(gPad->AbsPixeltoX(pxold));
         fY[ipoint] = gPad->PadtoY(gPad->AbsPixeltoY(pyold));
         if (InheritsFrom("TCutG")) {
            //make sure first and last point are the same
            if (ipoint == 0) {
               fX[fNpoints-1] = fX[0];
               fY[fNpoints-1] = fY[0];
            }
            if (ipoint == fNpoints-1) {
               fX[0] = fX[fNpoints-1];
               fY[0] = fY[fNpoints-1];
            }
         }
      }
      BADCASE = kFALSE;
      delete [] x; x = 0;
      delete [] y; y = 0;
      gPad->Modified(kTRUE);
      gVirtualX->SetLineColor(-1);

   }

}

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

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

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

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

//______________________________________________________________________________
 Int_t TGraph::Fit(const char *fname, Option_t *option, Option_t *, Axis_t xmin, Axis_t xmax)
{
//*-*-*-*-*-*Fit this graph with function with name fname*-*-*-*-*-*-*-*-*-*
//*-*        ============================================
//  interface to TF1::Fit(TF1 *f1...

   TF1 *f1 = (TF1*)gROOT->GetFunction(fname);
   if (!f1) { Printf("Unknown function: %s",fname); return -1; }
   return Fit(f1,option,"",xmin,xmax);
}

//______________________________________________________________________________
 Int_t TGraph::Fit(TF1 *f1, Option_t *option, Option_t *, Axis_t rxmin, Axis_t rxmax)
{
//*-*-*-*-*-*-*-*-*-*-*Fit this graph with function f1*-*-*-*-*-*-*-*-*-*
//*-*                  ==================================
//
//   f1 is an already predefined function created by TF1.
//   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)
//
//   When the fit is drawn (by default), the parameter goption may be used
//   to specify a list of graphics options. See TGraph::Paint for a complete
//   list of these options.
//
//   In order to use the Range option, one must first create a function
//   with the expression to be fitted. For example, if your graph
//   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:
//        TF1 *f1 = new TF1("f1","gaus",1,3);
//        graph->Fit("f1","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
//        f1->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 GraphFitChisquare 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 TF1*) 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:  TF1 *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, xmax, ymin, ymax;
   Int_t i, npar,nvpar,nparx;
   Double_t par, we, al, bl;
   Double_t eplus,eminus,eparab,globcc,amin,edm,errdef,werr;
   TF1 *fnew1;

   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(f1);
   grFitter->SetFitOption(fitOption);
   
   if (!f1) { Printf("Function is a null pointer"); return 0; }
   npar = f1->GetNpar();
   if (npar <=0) { Printf("Illegal number of parameters = %d",npar); return 0; }

//*-*- Check that function has same dimension as histogram
   if (f1->GetNdim() > 1) {
      Printf("Error function %s is not 1-D",f1->GetName()); return 0; }

//*-*- Is a Fit range specified?
   Int_t gxfirst, gxlast;
   if (fitOption.Range) {
      f1->GetRange(xmin, xmax);
      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 {
      f1->SetRange(xmin, xmax);
      gxfirst = 0;
      gxlast  = fNpoints-1;
   }

//*-*- If case of a predefined function, then compute initial values of parameters
   Int_t special = f1->GetNumber();
   if (fitOption.Bound) special = 0;
   if      (special == 100)      InitGaus(gxfirst,gxlast);
   else if (special == 200)      InitExpo(gxfirst,gxlast);
   else if (special == 299+npar) InitPolynom(gxfirst,gxlast);

//*-*- 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("GraphFitChisquare");
   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 = f1->GetParameter(i);
      f1->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,f1->GetParName(i),par,we,al,bl);
   }
   if(nfixed > 0)grFitter->ExecuteCommand("FIX",arglist,nfixed); // Otto

//*-*- Reset Print level
   if (!fitOption.Quiet) {
      if (fitOption.Verbose) { arglist[0] = 2; grFitter->ExecuteCommand("SET PRINT", arglist,1); }
      else                   { 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, sumw2=0;
   for (i=gxfirst;i<=gxlast;i++) {
      ex = GetErrorX(i);
      ey = GetErrorY(i);
      if (ex > 0 || ey > 0) hasErrors = kTRUE;
      sumw2 += ey*ey;
   }
//*-*- Perform minimization
   if (!InheritsFrom("TGraphErrors")) 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);
   f1->SetChisquare(amin);
   Int_t ndf = f1->GetNumberFitPoints()-npar+nfixed;
   f1->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));
      f1->SetParameter(i,par);
      f1->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;
         }
      }
      fnew1 = new TF1();
      f1->Copy(*fnew1);
      fFunctions->Add(fnew1);
      fnew1->SetParent(this);
      fnew1->Save(xmin,xmax,0,0,0,0);
      if (fitOption.Nograph) fnew1->SetBit(TF1::kNotDraw);
      fnew1->SetBit(TFormula::kNotGlobal);

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

//______________________________________________________________________________
 void TGraph::FitPanel()
{
//*-*-*-*-*-*-*Display a panel with all graph fit options*-*-*-*-*-*
//*-*          ==========================================
//
//   See class TFitPanel for example

   if (gPad) {
      gROOT->SetSelectedPrimitive(this);
      gROOT->SetSelectedPad(gPad);
   } else {
      Error("FitPanelGraph", "need to draw graph first");
      return;
   }

   //The pad utility manager is required (a plugin)
   TVirtualUtilPad *util = (TVirtualUtilPad*)gROOT->GetListOfSpecials()->FindObject("R__TVirtualUtilPad");
   if (!util) {
      TPluginHandler *h;
      if ((h = gROOT->GetPluginManager()->FindHandler("TVirtualUtilPad"))) {
          if (h->LoadPlugin() == -1)
            return;
          h->ExecPlugin(0);
          util = (TVirtualUtilPad*)gROOT->GetListOfSpecials()->FindObject("R__TVirtualUtilPad");
      }
   }
   util->FitPanelGraph();
}
//______________________________________________________________________________
 Double_t TGraph::GetCorrelationFactor() const
{
// Return graph correlation factor

  Double_t rms1 = GetRMS(1);
  if (rms1 == 0) return 0;
  Double_t rms2 = GetRMS(2);
  if (rms2 == 0) return 0;
  return GetCovariance()/rms1/rms2;
}

//______________________________________________________________________________
 Double_t TGraph::GetCovariance() const
{
// Return covariance of vectors x,y

  if (fNpoints <= 0) return 0;
  Double_t sum = fNpoints, sumx = 0, sumy = 0, sumxy = 0;

  for (Int_t i=0;i<fNpoints;i++) {
     sumx  += fX[i];
     sumy  += fY[i];
     sumxy += fX[i]*fY[i];
  }

  return sumxy/sum - sumx/sum*sumy/sum;
}

//______________________________________________________________________________
 Double_t TGraph::GetMean(Int_t axis) const
{
// Return mean value of X (axis=1)  or Y (axis=2)

  if (axis < 1 || axis > 2) return 0;
  if (fNpoints <= 0) return 0;
  Double_t sumx = 0;
  for (Int_t i=0;i<fNpoints;i++) {
     if (axis == 1) sumx += fX[i];
     else           sumx += fY[i];
  }
  return sumx/fNpoints;
}

//______________________________________________________________________________
 Double_t TGraph::GetRMS(Int_t axis) const
{
// Return RMS of X (axis=1)  or Y (axis=2)

  if (axis < 1 || axis > 2) return 0;
  if (fNpoints <= 0) return 0;
  Double_t sumx = 0, sumx2 = 0;
  for (Int_t i=0;i<fNpoints;i++) {
     if (axis == 1) {sumx += fX[i]; sumx2 += fX[i]*fX[i];}
     else           {sumx += fY[i]; sumx2 += fY[i]*fY[i];}
  }
  Double_t x = sumx/fNpoints;
  Double_t rms2 = TMath::Abs(sumx2/fNpoints -x*x);
  return TMath::Sqrt(rms2);
}

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

   return -1;

}

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

   return -1;

}

//______________________________________________________________________________
 TF1 *TGraph::GetFunction(const char *name) const
{
//*-*-*-*-*Return pointer to function with name*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*      ===================================
//
// Functions such as TGraph::Fit store the fitted function in the list of
// functions of this graph.

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

//______________________________________________________________________________
 TH1F *TGraph::GetHistogram() const
{
//    Returns a pointer to the histogram used to draw the axis
//    Takes into account the two following cases.
//       1- option 'A' was specified in TGraph::Draw. Return fHistogram
//       2- user had called TPad::DrawFrame. return pointer to hframe histogram

  if (fHistogram) return fHistogram;
  Double_t rwxmin,rwxmax, rwymin, rwymax, maximum, minimum, dx, dy;
  Double_t uxmin, uxmax;
  rwxmin = rwxmax = fX[0];
  rwymin = rwymax = fY[0];
  for (Int_t i=1;i<fNpoints;i++) {
     if (fX[i] < rwxmin) rwxmin = fX[i];
     if (fX[i] > rwxmax) rwxmax = fX[i];
     if (fY[i] < rwymin) rwymin = fY[i];
     if (fY[i] > rwymax) rwymax = fY[i];
  }

  ComputeRange(rwxmin, rwymin, rwxmax, rwymax);  //this is redefined in TGraphErrors

  if (rwxmin == rwxmax) rwxmax += 1.;
  if (rwymin == rwymax) rwymax += 1.;
  dx = 0.1*(rwxmax-rwxmin);
  dy = 0.1*(rwymax-rwymin);
  uxmin    = rwxmin - dx;
  uxmax    = rwxmax + dx;
  minimum  = rwymin - dy;
  maximum  = rwymax + dy;
  if (fMinimum != -1111) minimum = fMinimum;
  if (fMaximum != -1111) maximum = fMaximum;
  
      // the graph is created with at least as many channels as there are points
      // to permit zooming on the full range
  if (uxmin < 0 && rwxmin >= 0) {
     if (gPad && gPad->GetLogx()) uxmin = 0.9*rwxmin;
     else                 uxmin = 0;
  }
  if (uxmax > 0 && rwxmax <= 0) {
     if (gPad && gPad->GetLogx()) uxmax = 1.1*rwxmax;
     else                 uxmax = 0;
  }
  if (minimum < 0 && rwymin >= 0) {
     if(gPad && gPad->GetLogy()) minimum = 0.9*rwymin;
     else                minimum = 0;
  }
  if (minimum <= 0 && gPad && gPad->GetLogy()) minimum = 0.001*maximum;
  if (uxmin <= 0 && gPad && gPad->GetLogx()) {
     if (uxmax > 1000) uxmin = 1;
     else              uxmin = 0.001*uxmax;
  }

  rwxmin = uxmin;
  rwxmax = uxmax;
  Int_t npt = 100;
  if (fNpoints > npt) npt = fNpoints;
  ((TGraph*)this)->fHistogram = new TH1F(GetName(),GetTitle(),npt,rwxmin,rwxmax);
  if (!fHistogram) return 0;
  fHistogram->SetMinimum(minimum);
  fHistogram->SetBit(TH1::kNoStats);
  fHistogram->SetMaximum(maximum);
  fHistogram->GetYaxis()->SetLimits(minimum,maximum);
  fHistogram->SetDirectory(0);
  return fHistogram;
}

//______________________________________________________________________________
 void TGraph::GetPoint(Int_t i, Double_t &x, Double_t &y)
{
//*-*-*-*-*-*-*-*-*-*-*Get x and y values for point number i*-*-*-*-*-*-*-*-*
//*-*                  =====================================

   if (i < 0 || i >= fNpoints) return;
   if (!fX || !fY) return;
   x = fX[i];
   y = fY[i];
}

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

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

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

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

//______________________________________________________________________________
 void TGraph::InitGaus(Int_t first, Int_t last)
{
//*-*-*-*-*-*Compute Initial values of parameters for a gaussian*-*-*-*-*-*-*
//*-*        ===================================================

   Double_t allcha, sumx, sumx2, x, val, rms, mean;
   Int_t bin;
   const Double_t sqrtpi = 2.506628;

//*-*- Compute mean value and RMS of the graph in the given range
   if (last <= first) {
      first = 0;
      last  = fNpoints-1;
   }
   Int_t np = 0;
   allcha = sumx = sumx2 = 0;
   for (bin=first;bin<=last;bin++) {
      x       = fX[bin];
      if (x < fX[first] || x > fX[last]) continue;
      np++;
      val     = fY[bin];
      sumx   += val*x;
      sumx2  += val*x*x;
      allcha += val;
   }
   if (np == 0 || allcha == 0) return;
   mean = sumx/allcha;
   rms  = TMath::Sqrt(sumx2/allcha - mean*mean);
   Double_t binwidx = TMath::Abs((fX[last] - fX[first])/np);
   if (rms == 0) rms = 1;
   TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
   TF1 *f1 = (TF1*)grFitter->GetUserFunc();
   f1->SetParameter(0,binwidx*allcha/(sqrtpi*rms));
   f1->SetParameter(1,mean);
   f1->SetParameter(2,rms);
   f1->SetParLimits(2,0,10*rms);
}

//______________________________________________________________________________
 void TGraph::InitExpo(Int_t first, Int_t last)
{
//*-*-*-*-*-*Compute Initial values of parameters for an exponential*-*-*-*-*
//*-*        =======================================================

   Double_t constant, slope;
   Int_t ifail;
   if (last <= first) {
      first = 0;
      last  = fNpoints-1;
   }
   Int_t nchanx = last - first + 1;

   LeastSquareLinearFit(-nchanx, constant, slope, ifail, first, last);

   TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
   TF1 *f1 = (TF1*)grFitter->GetUserFunc();
   f1->SetParameter(0,constant);
   f1->SetParameter(1,slope);

}

//______________________________________________________________________________
 void TGraph::InitPolynom(Int_t first, Int_t last)
{
//*-*-*-*-*-*Compute Initial values of parameters for a polynom*-*-*-*-*-*-*
//*-*        ===================================================

   Double_t fitpar[25];

   TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
   TF1 *f1 = (TF1*)grFitter->GetUserFunc();
   Int_t npar   = f1->GetNpar();

   LeastSquareFit(npar, fitpar, first, last);

   for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
}

//______________________________________________________________________________
 Int_t TGraph::InsertPoint()
{
// Insert a new point at the mouse position

   Int_t px = gPad->GetEventX();
   Int_t py = gPad->GetEventY();

   //localize point where to insert
   Int_t ipoint = -2;
   Int_t i,d=0;
   // start with a small window (in case the mouse is very close to one point)
   for (i=0;i<fNpoints-1;i++) {
      d = DistancetoLine(px, py, gPad->XtoPad(fX[i]), gPad->YtoPad(fY[i]), gPad->XtoPad(fX[i+1]), gPad->YtoPad(fY[i+1]));
      if (d < 5) {ipoint = i+1; break;}
   }
   if (ipoint == -2) {
      //may be we are far from one point, try again with a larger window
      for (i=0;i<fNpoints-1;i++) {
         d = DistancetoLine(px, py, gPad->XtoPad(fX[i]), gPad->YtoPad(fY[i]), gPad->XtoPad(fX[i+1]), gPad->YtoPad(fY[i+1]));
         if (d < 10) {ipoint = i+1; break;}
      }
   }
   if (ipoint == -2) {
      //distinguish between first and last point
      Int_t dpx = px - gPad->XtoAbsPixel(gPad->XtoPad(fX[0]));
      Int_t dpy = py - gPad->YtoAbsPixel(gPad->XtoPad(fY[0]));
      if (dpx*dpx+dpy*dpy < 25) ipoint = 0;
      else                      ipoint = fNpoints;
   }
   fNpoints++;
   Double_t *newX = new Double_t[fNpoints];
   Double_t *newY = new Double_t[fNpoints];
   for (i=0;i<ipoint;i++) {
      newX[i] = fX[i];
      newY[i] = fY[i];
   }
   newX[ipoint] = gPad->PadtoX(gPad->AbsPixeltoX(px));
   newY[ipoint] = gPad->PadtoY(gPad->AbsPixeltoY(py));
   for (i=ipoint+1;i<fNpoints;i++) {
      newX[i] = fX[i-1];
      newY[i] = fY[i-1];
   }
   delete [] fX;
   delete [] fY;
   fX = newX;
   fY = newY;
   gPad->Modified();
   return ipoint;
}

//______________________________________________________________________________
 void TGraph::LeastSquareFit(Int_t m, Double_t *a, Int_t first, Int_t last)
{
//*-*-*-*-*-*-*-*Least squares lpolynomial fitting without weights*-*-*-*-*-*-*
//*-*            =================================================
//
//  m     number of parameters
//  a     array of parameters
//  first 1st point number to fit (default =0)
//  last  last point number to fit (default=fNpoints-1)
//
//   based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
//
//
    const Double_t zero = 0.;
    const Double_t one = 1.;
    const Int_t idim = 20;

    Double_t  b[400]	/* was [20][20] */;
    Int_t i, k, l, ifail;
    Double_t power;
    Double_t da[20], xk, yk;

    if (last <= first) {
       first = 0;
       last  = fNpoints-1;
    }
    Int_t n = last-first+1;

    if (m <= 2) {
       LeastSquareLinearFit(n, a[0], a[1], ifail, first, last);
       return;
    }
    if (m > idim || m > n) return;
    da[0] = zero;
    for (l = 2; l <= m; ++l) {
	b[l-1]           = zero;
	b[m + l*20 - 21] = zero;
	da[l-1]          = zero;
    }
    Int_t np = 0;
    for (k = first; k <= last; ++k) {
	xk     = fX[k];
	if (xk < fX[first] || xk > fX[last]) continue;
        np++;
        yk     = fY[k];
	power  = one;
	da[0] += yk;
	for (l = 2; l <= m; ++l) {
	    power   *= xk;
	    b[l-1]  += power;
	    da[l-1] += power*yk;
	}
	for (l = 2; l <= m; ++l) {
	    power            *= xk;
	    b[m + l*20 - 21] += power;
	}
    }
    b[0]  = Double_t(np);
    for (i = 3; i <= m; ++i) {
	for (k = i; k <= m; ++k) {
	    b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
	}
    }
    H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);

    if (ifail < 0) {
       a[0] = fY[0];
       for (i=1; i<m; ++i) a[i] = 0;
       return;
    }
    for (i=0; i<m; ++i) a[i] = da[i];

}

//______________________________________________________________________________
 void TGraph::LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail, Int_t first, Int_t last)
{
//*-*-*-*-*-*-*-*-*-*Least square linear fit without weights*-*-*-*-*-*-*-*-*
//*-*                =======================================
//
//   extracted from CERNLIB LLSQ: Translated to C++ by Rene Brun
//

    Double_t xbar, ybar, x2bar;
    Int_t i;
    Double_t xybar;
    Double_t fn, xk, yk;
    Double_t det;

    ifail = -2;
    xbar  = ybar = x2bar = xybar = 0;
    Int_t np = 0;
    for (i = first; i <= last; ++i) {
	xk = fX[i];
	if (xk < fX[first] || xk > fX[last]) continue;
        np++;
	yk = fY[i];
	if (ndata < 0) {
	    if (yk <= 0) yk = 1e-9;
	    yk = TMath::Log(yk);
	}
	xbar  += xk;
	ybar  += yk;
	x2bar += xk*xk;
	xybar += xk*yk;
    }
    fn    = Double_t(np);
    det   = fn*x2bar - xbar*xbar;
    ifail = -1;
    if (det <= 0) {
       if (fn > 0) a0 = ybar/fn;
       else        a0 = 0;
       a1 = 0;
       return;
    }
    ifail = 0;
    a0 = (x2bar*ybar - xbar*xybar) / det;
    a1 = (fn*xybar - xbar*ybar) / det;

}

//______________________________________________________________________________
 void TGraph::Paint(Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this graph with its current attributes*-*-*-*-*-*-*
//*-*                  ===========================================

   if (strstr(option,"H") || strstr(option,"h")) {
      PaintGrapHist(fNpoints, fX, fY, option);
   } else {
      PaintGraph(fNpoints, fX, fY, option);
   }

   TObject *f;
   if (fFunctions) {
     TIter   next(fFunctions);
     while ((f = (TObject*) next())) {
      if (f->InheritsFrom(TF1::Class())) {
         if (f->TestBit(TF1::kNotDraw) == 0) f->Paint("lsame");
      } else  {
         f->Paint();
      }
     }
   }
}

//______________________________________________________________________________
 void TGraph::PaintFit(TF1 *fit)
{
//  Paint "stats" box with the fit info

   Int_t dofit;
   TPaveStats *stats  = 0;
   TIter next(fFunctions);
   TObject *obj;
   while ((obj = next())) {
      if (obj->InheritsFrom(TPaveStats::Class())) {
         stats = (TPaveStats*)obj;
         break;
      }
   }

   if (stats) dofit  = stats->GetOptFit();
   else       dofit  = gStyle->GetOptFit();

   if (!dofit) fit = 0;
   if (!fit) return;
   if (dofit  == 1) dofit  =  111;
   Int_t nlines = 0;
   Int_t print_fval    = dofit%10;
   Int_t print_ferrors = (dofit/10)%10;
   Int_t print_fchi2   = (dofit/100)%10;
   Int_t print_fprob   = (dofit/1000)%10;
   Int_t nlinesf = print_fval + print_fchi2 + print_fprob;
   if (fit) nlinesf += fit->GetNpar();
   Bool_t done = kFALSE;
   Double_t  statw  = 1.8*gStyle->GetStatW();
   Double_t  stath  = 0.25*(nlines+nlinesf)*gStyle->GetStatH();
   if (stats) {
      stats->Clear();
      done = kTRUE;
   } else {
      stats  = new TPaveStats(
               gStyle->GetStatX()-statw,
               gStyle->GetStatY()-stath,
               gStyle->GetStatX(),
               gStyle->GetStatY(),"brNDC");

      stats->SetParent(fFunctions);
      stats->SetOptFit(dofit);
      stats->SetOptStat(0);
      stats->SetFillColor(gStyle->GetStatColor());
      stats->SetFillStyle(gStyle->GetStatStyle());
      stats->SetBorderSize(gStyle->GetStatBorderSize());
      stats->SetTextFont(gStyle->GetStatFont());
      if (gStyle->GetStatFont()%10 > 2)
         stats->SetTextSize(gStyle->GetStatFontSize());
      stats->SetFitFormat(gStyle->GetFitFormat());
      stats->SetStatFormat(gStyle->GetStatFormat());
      stats->SetName("stats");

      stats->SetTextColor(gStyle->GetStatTextColor());
      stats->SetTextAlign(12);
      stats->SetBit(kCanDelete);
      stats->SetBit(kMustCleanup);
   }

   char t[64];
   char textstats[50];
   Int_t ndf = fit->GetNDF();
   sprintf(textstats,"#chi^{2} / ndf = %s%s / %d","%",stats->GetFitFormat(),ndf);
   sprintf(t,textstats,(Float_t)fit->GetChisquare());
   if (print_fchi2) stats->AddText(t);
   if (print_fprob) {
      sprintf(textstats,"Prob  = %s%s","%",stats->GetFitFormat());
      sprintf(t,textstats,(Float_t)TMath::Prob(fit->GetChisquare(),ndf));
      stats->AddText(t);
   }
   if (print_fval || print_ferrors) {
      for (Int_t ipar=0;ipar<fit->GetNpar();ipar++) {
         if (print_ferrors) {
            sprintf(textstats,"%-8s = %s%s #pm %s%s ",fit->GetParName(ipar),"%",stats->GetFitFormat(),"%",stats->GetFitFormat());
            sprintf(t,textstats,(Float_t)fit->GetParameter(ipar)
                            ,(Float_t)fit->GetParError(ipar));
         } else {
            sprintf(textstats,"%-8s = %s%s ",fit->GetParName(ipar),"%",stats->GetFitFormat());
            sprintf(t,textstats,(Float_t)fit->GetParameter(ipar));
         }
         t[63] = 0;
         stats->AddText(t);
      }
   }

   if (!done) fFunctions->Add(stats);
   stats->Paint();
}

//______________________________________________________________________________
 void TGraph::PaintGraph(Int_t npoints, const Double_t *x, const Double_t *y, Option_t *chopt)
{
//*-*-*-*-*-*-*-*-*-*-*-*Control function to draw a graph*-*-*-*-*-*-*-*-*-*-*
//*-*                    ================================
//
//
//  Draws one dimensional graphs. The aspect of the graph is done
//  according to the value of the chopt.
//
//  Input parameters:
//
//  npoints : Number of points in X or in Y.
//  x[npoints] or x[2] : X coordinates or (XMIN,XMAX) (WC space).
//  y[npoints] or y[2] : Y coordinates or (YMIN,YMAX) (WC space).
//  chopt : Option.
//
//  chopt='L' :  A simple polyline between every points is drawn
//
//  chopt='F' :  A fill area is drawn ('CF' draw a smooth fill area)
//
//  chopt='A' :  Axis are drawn around the graph
//
//  chopt='C' :  A smooth Curve is drawn
//
//  chopt='*' :  A Star is plotted at each point
//
//  chopt='P' :  Idem with the current marker
//
//  chopt='B' :  A Bar chart is drawn at each point
//
//  chopt='1' :  ylow=rwymin
//
//

   Int_t OptionLine , OptionAxis , OptionCurve,OptionStar ,OptionMark;
   Int_t OptionBar  , OptionR    , OptionOne;
   Int_t OptionFill , OptionZ    , OptionCurveFill;
   Int_t i, npt, nloop;
   Int_t drawtype=0;
   Double_t xlow, xhigh, ylow, yhigh;
   Double_t barxmin, barxmax, barymin, barymax;
   Double_t uxmin, uxmax;
   Double_t x1, xn, y1, yn;
   Double_t dbar, bdelta;

//*-* ______________________________________

  if (npoints <= 0) {
     Error("PaintGraph", "illegal number of points (%d)", npoints);
     return;
  }
   TString opt = chopt;
   opt.ToUpper();

   if(opt.Contains("L")) OptionLine = 1;  else OptionLine = 0;
   if(opt.Contains("A")) OptionAxis = 1;  else OptionAxis = 0;
   if(opt.Contains("C")) OptionCurve= 1;  else OptionCurve= 0;
   if(opt.Contains("*")) OptionStar = 1;  else OptionStar = 0;
   if(opt.Contains("P")) OptionMark = 1;  else OptionMark = 0;
   if(opt.Contains("B")) OptionBar  = 1;  else OptionBar  = 0;
   if(opt.Contains("R")) OptionR    = 1;  else OptionR    = 0;
   if(opt.Contains("1")) OptionOne  = 1;  else OptionOne  = 0;
   if(opt.Contains("F")) OptionFill = 1;  else OptionFill = 0;
   OptionZ    = 0;
//*-*           If no "drawing" option is selected and if chopt<>' '
//*-*           nothing is done.

  if (OptionLine+OptionFill+OptionCurve+OptionStar+OptionMark+OptionBar == 0) {
     if (strlen(chopt) == 0)  OptionLine=1;
     else   return;
  }

  if (OptionStar) SetMarkerStyle(3);

  OptionCurveFill = 0;
  if (OptionCurve && OptionFill) {
     OptionCurveFill = 1;
     OptionFill      = 0;
  }


//*-*-           Draw the Axis with a fixed number of division: 510

  Double_t rwxmin,rwxmax, rwymin, rwymax, maximum, minimum, dx, dy;

  if (OptionAxis) {

     if (fHistogram) {
        rwxmin    = gPad->GetUxmin();
        rwxmax    = gPad->GetUxmax();
        rwymin    = gPad->GetUymin();
        rwymax    = gPad->GetUymax();
        minimum   = fHistogram->GetMinimumStored();
        maximum   = fHistogram->GetMaximumStored();
        if (minimum == -1111) { //this can happen after unzooming
           minimum = fHistogram->GetYaxis()->GetXmin();
           fHistogram->SetMinimum(minimum);
        }
        if (maximum == -1111) {
           maximum = fHistogram->GetYaxis()->GetXmax();
           fHistogram->SetMaximum(maximum);
        }
        uxmin     = gPad->PadtoX(rwxmin);
        uxmax     = gPad->PadtoX(rwxmax);
     } else {
        rwxmin = rwxmax = x[0];
        rwymin = rwymax = y[0];
        for (i=1;i<npoints;i++) {
           if (x[i] < rwxmin) rwxmin = x[i];
           if (x[i] > rwxmax) rwxmax = x[i];
           if (y[i] < rwymin) rwymin = y[i];
           if (y[i] > rwymax) rwymax = y[i];
        }

        ComputeRange(rwxmin, rwymin, rwxmax, rwymax);  //this is redefined in TGraphErrors

        if (rwxmin == rwxmax) rwxmax += 1.;
        if (rwymin == rwymax) rwymax += 1.;
        dx = 0.1*(rwxmax-rwxmin);
        dy = 0.1*(rwymax-rwymin);
        uxmin    = rwxmin - dx;
        uxmax    = rwxmax + dx;
        minimum  = rwymin - dy;
        maximum  = rwymax + dy;
     }
     if (fMinimum != -1111) rwymin = minimum = fMinimum;
     if (fMaximum != -1111) rwymax = maximum = fMaximum;
     if (uxmin < 0 && rwxmin >= 0) {
        if (gPad->GetLogx()) uxmin = 0.9*rwxmin;
        else                 uxmin = 0;
     }
     if (uxmax > 0 && rwxmax <= 0) {
        if (gPad->GetLogx()) uxmax = 1.1*rwxmax;
        else                 uxmax = 0;
     }
     if (minimum < 0 && rwymin >= 0) {
        if(gPad->GetLogy()) minimum = 0.9*rwymin;
        else                minimum = 0;
     }
     if (maximum > 0 && rwymax <= 0) {
        //if(gPad->GetLogy()) maximum = 1.1*rwymax;
        //else                maximum = 0;
     }
     if (minimum <= 0 && gPad->GetLogy()) minimum = 0.001*maximum;
     if (uxmin <= 0 && gPad->GetLogx()) {
        if (uxmax > 1000) uxmin = 1;
        else              uxmin = 0.001*uxmax;
     }
     rwymin = minimum;
     rwymax = maximum;

//*-*-  Create a temporary histogram and fill each channel with the function value
   if (!fHistogram) {
      // the graph is created with at least as many channels as there are points
      // to permit zooming on the full range
      rwxmin = uxmin;
      rwxmax = uxmax;
      npt = 100;
      if (fNpoints > npt) npt = fNpoints;
      fHistogram = new TH1F(GetName(),GetTitle(),npt,rwxmin,rwxmax);
      if (!fHistogram) return;
      fHistogram->SetMinimum(rwymin);
      fHistogram->SetBit(TH1::kNoStats);
      fHistogram->SetMaximum(rwymax);
      fHistogram->GetYaxis()->SetLimits(rwymin,rwymax);
      fHistogram->SetDirectory(0);
   }
   fHistogram->Paint("a");
  }

  //*-*  Set Clipping option
  gPad->SetBit(kClipFrame, TestBit(kClipFrame));

  TF1 *fit = 0;
  if (fFunctions) fit = (TF1*)fFunctions->First();
  TObject *f;
  if (fFunctions) {
     TIter   next(fFunctions);
     while ((f = (TObject*) next())) {
        if (f->InheritsFrom(TF1::Class())) {
           fit = (TF1*)f;
           break;
        }
     }
  }
  if (fit) PaintFit(fit);

  rwxmin   = gPad->GetUxmin();
  rwxmax   = gPad->GetUxmax();
  rwymin   = gPad->GetUymin();
  rwymax   = gPad->GetUymax();
  uxmin    = gPad->PadtoX(rwxmin);
  uxmax    = gPad->PadtoX(rwxmax);
  if (fHistogram) {
     maximum = fHistogram->GetMaximum();
     minimum = fHistogram->GetMinimum();
  } else {
     maximum = gPad->PadtoY(rwymax);
     minimum = gPad->PadtoY(rwymin);
  }

//*-*-           Set attributes
  TAttLine::Modify();
  TAttFill::Modify();
  TAttMarker::Modify();

//*-*-           Draw the graph with a polyline or a fill area

  gxwork  = new Double_t[2*npoints+10];
  gywork  = new Double_t[2*npoints+10];
  gxworkl = new Double_t[2*npoints+10];
  gyworkl = new Double_t[2*npoints+10];

  if (OptionLine || OptionFill) {
     x1       = x[0];
     xn       = x[npoints-1];
     y1       = y[0];
     yn       = y[npoints-1];
     nloop = npoints;
     if (OptionFill && (xn != x1 || yn != y1)) nloop++;
     npt = 0;
     for (i=1;i<=nloop;i++) {
        if (i > npoints) {
           gxwork[npt] = gxwork[0];  gywork[npt] = gywork[0];
        } else {
           gxwork[npt] = x[i-1];      gywork[npt] = y[i-1];
           npt++;
        }
        if (i == nloop) {
           ComputeLogs(npt, OptionZ);
           Int_t bord = gStyle->GetDrawBorder();
           if (OptionR) {
              if (OptionFill) {
                 gPad->PaintFillArea(npt,gyworkl,gxworkl);
                 if (bord) gPad->PaintPolyLine(npt,gyworkl,gxworkl);
              }
              else         gPad->PaintPolyLine(npt,gyworkl,gxworkl);
           }
           else {
              if (OptionFill) {
                 gPad->PaintFillArea(npt,gxworkl,gyworkl);
                 if (bord) gPad->PaintPolyLine(npt,gxworkl,gyworkl);
              }
              else         gPad->PaintPolyLine(npt,gxworkl,gyworkl);
           }
           gxwork[0] = gxwork[npt-1];  gywork[0] = gywork[npt-1];
           npt      = 1;
        }
     }
  }

//*-*-           Draw the graph with a smooth Curve. Smoothing via Smooth

  if (OptionCurve) {
     x1 = x[0];
     xn = x[npoints-1];
     y1 = y[0];
     yn = y[npoints-1];
     drawtype = 1;
     nloop = npoints;
     if (OptionCurveFill) {
        drawtype += 1000;
        if (xn != x1 || yn != y1) nloop++;
     }
     if (!OptionR) {
        npt = 0;
        for (i=1;i<=nloop;i++) {
           if (i > npoints) {
              gxwork[npt] = gxwork[0];  gywork[npt] = gywork[0];
           } else {
              if (y[i-1] < minimum || y[i-1] > maximum) continue;
              if (x[i-1] < uxmin    || x[i-1] > uxmax)  continue;
              gxwork[npt] = x[i-1];      gywork[npt] = y[i-1];
              npt++;
           }
           ComputeLogs(npt, OptionZ);
           if (gyworkl[npt-1] < rwymin || gyworkl[npt-1] > rwymax) {
              if (npt > 2) {
                 ComputeLogs(npt, OptionZ);
                 Smooth(npt,gxworkl,gyworkl,drawtype);
              }
              gxwork[0] = gxwork[npt-1]; gywork[0] = gywork[npt-1];
              npt=1;
              continue;
           }
        }  //endfor (i=0;i<nloop;i++)
        if (npt > 1) {
           ComputeLogs(npt, OptionZ);
           Smooth(npt,gxworkl,gyworkl,drawtype);
        }
     }
     else {
        drawtype += 10;
        npt    = 0;
        for (i=1;i<=nloop;i++) {
           if (i > npoints) {
              gxwork[npt] = gxwork[0];  gywork[npt] = gywork[0];
           } else {
              if (y[i-1] < minimum || y[i-1] > maximum) continue;
              if (x[i-1] < uxmin    || x[i-1] > uxmax)  continue;
              gxwork[npt] = x[i-1];      gywork[npt] = y[i-1];
              npt++;
           }
           ComputeLogs(npt, OptionZ);
           if (gxworkl[npt-1] < rwxmin || gxworkl[npt-1] > rwxmax) {
              if (npt > 2) {
                 ComputeLogs(npt, OptionZ);
                 Smooth(npt,gxworkl,gyworkl,drawtype);
              }
              gxwork[0] = gxwork[npt-1]; gywork[0] = gywork[npt-1];
              npt=1;
              continue;
           }
        } //endfor (i=1;i<=nloop;i++)
        if (npt > 1) {
           ComputeLogs(npt, OptionZ);
           Smooth(npt,gxworkl,gyworkl,drawtype);
        }
     }
  }

//*-*-           Draw the graph with a '*' on every points

  if (OptionStar) {
     SetMarkerStyle(3);
     npt = 0;
     for (i=1;i<=npoints;i++) {
        if (y[i-1] >= minimum && y[i-1] <= maximum && x[i-1] >= uxmin  && x[i-1] <= uxmax) {
           gxwork[npt] = x[i-1];      gywork[npt] = y[i-1];
           npt++;
        }
        if (i == npoints) {
           ComputeLogs(npt, OptionZ);
           if (OptionR)  gPad->PaintPolyMarker(npt,gyworkl,gxworkl);
           else          gPad->PaintPolyMarker(npt,gxworkl,gyworkl);
           npt = 0;
        }
     }
  }

//*-*-           Draw the graph with the current polymarker on
//*-*-           every points

  if (OptionMark) {
     npt = 0;
     for (i=1;i<=npoints;i++) {
        if (y[i-1] >= minimum && y[i-1] <= maximum && x[i-1] >= uxmin  && x[i-1] <= uxmax) {
           gxwork[npt] = x[i-1];      gywork[npt] = y[i-1];
           npt++;
        }
        if (i == npoints) {
           ComputeLogs(npt, OptionZ);
           if (OptionR) gPad->PaintPolyMarker(npt,gyworkl,gxworkl);
           else         gPad->PaintPolyMarker(npt,gxworkl,gyworkl);
           npt = 0;
        }
     }
  }

//*-*-           Draw the graph as a bar chart

  if (OptionBar) {
     if (!OptionR) {
        barxmin = x[0];
        barxmax = x[0];
        for (i=1;i<npoints;i++) {
           if (x[i] < barxmin) barxmin = x[i];
           if (x[i] > barxmax) barxmax = x[i];
        }
        bdelta = (barxmax-barxmin)/Double_t(npoints);
     }
     else {
        barymin = y[0];
        barymax = y[0];
        for (i=1;i<npoints;i++) {
           if (y[i] < barymin) barymin = y[i];
           if (y[i] > barymax) barymax = y[i];
        }
        bdelta = (barymax-barymin)/Double_t(npoints);
     }
     dbar  = 0.5*bdelta*gStyle->GetBarWidth();
     if (!OptionR) {
        for (i=1;i<=npoints;i++) {
           xlow  = x[i-1] - dbar;
           xhigh = x[i-1] + dbar;
           yhigh = y[i-1];
           if (xlow  < uxmin) continue;
           if (xhigh > uxmax) continue;
           if (!OptionOne) ylow = TMath::Max((Double_t)0,gPad->GetUymin());
           else            ylow = gPad->GetUymin();
           gxwork[0] = xlow;
           gywork[0] = ylow;
           gxwork[1] = xhigh;
           gywork[1] = yhigh;
           ComputeLogs(2, OptionZ);
           if (gyworkl[0] < gPad->GetUymin()) gyworkl[0] = gPad->GetUymin();
           if (gyworkl[1] < gPad->GetUymin()) continue;
           if (gyworkl[1] > gPad->GetUymax()) gyworkl[1] = gPad->GetUymax();
           if (gyworkl[0] > gPad->GetUymax()) continue;

           gPad->PaintBox(gxworkl[0],gyworkl[0],gxworkl[1],gyworkl[1], "B");
        }
     }
     else {
        for (i=1;i<=npoints;i++) {
           xhigh = x[i-1];
           ylow  = y[i-1] - dbar;
           yhigh = y[i-1] + dbar;
           xlow     = TMath::Max((Double_t)0, gPad->GetUxmin());
           gxwork[0] = xlow;
           gywork[0] = ylow;
           gxwork[1] = xhigh;
           gywork[1] = yhigh;
           ComputeLogs(2, OptionZ);
           gPad->PaintBox(gxworkl[0],gyworkl[0],gxworkl[1],gyworkl[1], "B");
        }
     }
  }
   gPad->ResetBit(kClipFrame);

   delete [] gxwork;
   delete [] gywork;
   delete [] gxworkl;
   delete [] gyworkl;
}

//______________________________________________________________________________
 void TGraph::PaintGrapHist(Int_t npoints, const Double_t *x, const Double_t *y, Option_t *chopt)
{
//*-*-*-*-*-*-*-*-*Control function to draw a graphistogram*-*-*-*-*-*-*-*-*-*
//*-*              ========================================
//
//
//   Draws one dimensional graphs. The aspect of the graph is done
// according to the value of the chopt.
//
// Input parameters:
//
//  npoints : Number of points in X or in Y.
//  X(N) or x[1] : X coordinates or (XMIN,XMAX) (WC space).
//  Y(N) or y[1] : Y coordinates or (YMIN,YMAX) (WC space).
//  chopt : Option.
//
//  chopt='R' :  Graph is drawn horizontaly, parallel to X axis.
//               (default is vertically, parallel to Y axis)
//               If option R is selected the user must give:
//                 2 values for Y (y[0]=YMIN and y[1]=YMAX)
//                 N values for X, one for each channel.
//               Otherwise the user must give:
//                 N values for Y, one for each channel.
//                 2 values for X (x[0]=XMIN and x[1]=XMAX)
//
//  chopt='L' :  A simple polyline beetwen every points is drawn
//
//  chopt='H' :  An Histogram with equidistant bins is drawn
//               as a polyline.
//
//  chopt='F' :  An histogram with equidistant bins is drawn
//               as a fill area. Contour is not drawn unless
//               chopt='H' is also selected..
//
//  chopt='N' :  Non equidistant bins (default is equidistant)
//               If N is the number of channels array X and Y
//               must be dimensionned as follow:
//               If option R is not selected (default) then
//               the user must give:
//                 (N+1) values for X (limits of channels).
//                  N values for Y, one for each channel.
//               Otherwise the user must give:
//                 (N+1) values for Y (limits of channels).
//                  N values for X, one for each channel.
//
//  chopt='F1':  Idem as 'F' except that fill area is no more
//               reparted arround axis X=0 or Y=0 .
//
//  chopt='F2':  Draw a Fill area polyline connecting the center of bins
//
//  chopt='C' :  A smooth Curve is drawn.
//
//  chopt='*' :  A Star is plotted at the center of each bin.
//
//  chopt='P' :  Idem with the current marker
//  chopt='P0':  Idem with the current marker. Empty bins also drawn
//
//  chopt='B' :  A Bar chart with equidistant bins is drawn as fill
//               areas (Contours are drawn).
//
//  chopt='9' :  Force graph to be drawn in high resolution mode.
//               By default, the graph is drawn in low resolution
//               in case the number of points is greater than the number of pixels
//               in the current pad.
//
//  chopt='][' : "Cutoff" style. When this option is selected together with
//               H option, the first and last vertical lines of the histogram
//               are not drawn.
//

   const char *where = "PaintGraphHist";

   Int_t OptionLine , OptionAxis , OptionCurve, OptionStar , OptionMark;
   Int_t OptionBar  , OptionRot  , OptionOne  , OptionOff;
   Int_t OptionFill , OptionZ;
   Int_t OptionHist , OptionBins , OptionMarker;
   Int_t i, j, npt;
   Int_t drawtype=0, drawborder, drawbordersav;
   Double_t xlow, xhigh, ylow, yhigh;
   Double_t wmin, wmax;
   Double_t dbar, offset, wminstep;
   Double_t delta = 0;
   Double_t ylast = 0;
   Double_t xi, xi1, xj, xj1, yi1, yi, yj, yj1, xwmin, ywmin;
   Int_t first, last, nbins;
   Int_t fillarea;

   char choptaxis[10] = " ";
//*-* ______________________________________

   if (npoints <= 0) {
      Error(where, "illegal number of points (%d)", npoints);
      return;
   }
   TString opt = chopt;
   opt.ToUpper();
   if(opt.Contains("H"))  OptionHist = 1;  else OptionHist = 0;
   if(opt.Contains("F"))  OptionFill = 1;  else OptionFill = 0;
   if(opt.Contains("C"))  OptionCurve= 1;  else OptionCurve= 0;
   if(opt.Contains("*"))  OptionStar = 1;  else OptionStar = 0;
   if(opt.Contains("R"))  OptionRot  = 1;  else OptionRot  = 0;
   if(opt.Contains("1"))  OptionOne  = 1;  else OptionOne  = 0;
   if(opt.Contains("B"))  OptionBar  = 1;  else OptionBar  = 0;
   if(opt.Contains("N"))  OptionBins = 1;  else OptionBins = 0;
   if(opt.Contains("L"))  OptionLine = 1;  else OptionLine = 0;
   if(opt.Contains("P"))  OptionMark = 1;  else OptionMark = 0;
   if(opt.Contains("A"))  OptionAxis = 1;  else OptionAxis = 0;
   if(opt.Contains("][")) OptionOff  = 1;  else OptionOff  = 0;
   if(opt.Contains("P0")) OptionMark = 10;

   Int_t OptionFill2 = 0;
   if(opt.Contains("F") && opt.Contains("2")) {
      OptionFill = 0; OptionFill2 = 1;
   }

//*-*  Set Clipping option
  gPad->SetBit(kClipFrame, TestBit(kClipFrame));

  OptionZ = 1;

  if (OptionStar) SetMarkerStyle(3);

  first = 1;
  last  = npoints;
  nbins = last - first + 1;

//*-*-           Draw the Axis with a fixed number of division: 510

  Double_t baroffset = gStyle->GetBarOffset();
  Double_t barwidth  = gStyle->GetBarWidth();
  Double_t rwxmin    = gPad->GetUxmin();
  Double_t rwxmax    = gPad->GetUxmax();
  Double_t rwymin    = gPad->GetUymin();
  Double_t rwymax    = gPad->GetUymax();
  Double_t uxmin     = gPad->PadtoX(rwxmin);
  Double_t uxmax     = gPad->PadtoX(rwxmax);
  Double_t rounding  = (uxmax-uxmin)*1.e-5;
  drawborder         = gStyle->GetDrawBorder();
  if (OptionAxis) {
     Int_t nx1, nx2, ndivx, ndivy, ndiv;
     choptaxis[0]  = 0;
     Double_t rwmin  = rwxmin;
     Double_t rwmax  = rwxmax;
     ndivx = gStyle->GetNdivisions("X");
     ndivy = gStyle->GetNdivisions("Y");
     if (ndivx > 1000) {
        nx2   = ndivx/100;
        nx1   = TMath::Max(1, ndivx%100);
        ndivx = 100*nx2 + Int_t(Double_t(nx1)*gPad->GetAbsWNDC());
     }
     ndiv  =TMath::Abs(ndivx);
     if (ndivx < 0) strcat(choptaxis, "N");
     if (gPad->GetGridx()) {
        strcat(choptaxis, "W");
     }
     if (gPad->GetLogx()) {
        rwmin = TMath::Power(10,rwxmin);
        rwmax = TMath::Power(10,rwxmax);
        strcat(choptaxis, "G");
     }
     TGaxis *axis = new TGaxis();
     axis->SetLineColor(gStyle->GetAxisColor("X"));
     axis->SetTextColor(gStyle->GetLabelColor("X"));
     axis->SetTextFont(gStyle->GetLabelFont("X"));
     axis->SetLabelSize(gStyle->GetLabelSize("X"));
     axis->SetLabelOffset(gStyle->GetLabelOffset("X"));
     axis->SetTickSize(gStyle->GetTickLength("X"));

     axis->PaintAxis(rwxmin,rwymin,rwxmax,rwymin,rwmin,rwmax,ndiv,choptaxis);

     choptaxis[0]  = 0;
     rwmin  = rwymin;
     rwmax  = rwymax;
     if (ndivy < 0) {
       nx2   = ndivy/100;
       nx1   = TMath::Max(1, ndivy%100);
       ndivy = 100*nx2 + Int_t(Double_t(nx1)*gPad->GetAbsHNDC());
       strcat(choptaxis, "N");
     }
     ndiv  =TMath::Abs(ndivy);
     if (gPad->GetGridy()) {
        strcat(choptaxis, "W");
     }
     if (gPad->GetLogy()) {
        rwmin = TMath::Power(10,rwymin);
        rwmax = TMath::Power(10,rwymax);
        strcat(choptaxis,"G");
     }
     axis->SetLineColor(gStyle->GetAxisColor("Y"));
     axis->SetTextColor(gStyle->GetLabelColor("Y"));
     axis->SetTextFont(gStyle->GetLabelFont("Y"));
     axis->SetLabelSize(gStyle->GetLabelSize("Y"));
     axis->SetLabelOffset(gStyle->GetLabelOffset("Y"));
     axis->SetTickSize(gStyle->GetTickLength("Y"));

     axis->PaintAxis(rwxmin,rwymin,rwxmin,rwymax,rwmin,rwmax,ndiv,choptaxis);
     delete axis;
  }


//*-*-           Set attributes
  TAttLine::Modify();
  TAttFill::Modify();
  TAttMarker::Modify();

//*-*-       Min-Max scope

  if (!OptionRot) {wmin = x[0];   wmax = x[1];}
  else            {wmin = y[0];   wmax = y[1];}

  if (!OptionBins) delta = (wmax - wmin)/ Double_t(nbins);

//*-*-           Draw the histogram with a fill area

  gxwork  = new Double_t[2*npoints+10];
  gywork  = new Double_t[2*npoints+10];
  gxworkl = new Double_t[2*npoints+10];
  gyworkl = new Double_t[2*npoints+10];

  if (OptionFill && !OptionCurve) {
     fillarea = kTRUE;
     if (!OptionRot) {
        gxwork[0] = wmin;
        if (!OptionOne) gywork[0] = TMath::Max((Double_t)0,gPad->GetUymin());
        else            gywork[0] = gPad->GetUymin();
        npt = 2;
        for (j=first; j<=last;j++) {
           if (!OptionBins) {
              gxwork[npt-1]   = gxwork[npt-2];
              gxwork[npt]     = wmin+((j-first+1)*delta);
           }
           else {
              xj1 = x[j];      xj  = x[j-1];
              if (xj1 < xj) {
                 if (j != last) Error(where, "X must be in increasing order");
                 else           Error(where, "X must have N+1 values with option N");
                 return;
              }
              gxwork[npt-1] = x[j-1];       gxwork[npt] = x[j];
           }
           gywork[npt-1] = y[j-1];
           gywork[npt]   = y[j-1];
           if (gxwork[npt-1] >= uxmin-rounding && gxwork[npt] <= uxmax+rounding) npt += 2;
           else gxwork[npt-2] = TMath::Min(gxwork[npt], uxmax);
           if (j == last) {
              gxwork[npt-1] = gxwork[npt-2];
              gywork[npt-1] = gywork[0];
              ComputeLogs(npt, OptionZ);
              gPad->PaintFillArea(npt,gxworkl,gyworkl);
              if (drawborder) {
                 if (!fillarea) gyworkl[0] = ylast;
                 gPad->PaintPolyLine(npt-1,gxworkl,gyworkl);
              }
              continue;
           }
        }  //endfor (j=first; j<=last;j++) {
     }
     else {
        gywork[0] = wmin;
        if (!OptionOne) gxwork[0] = TMath::Max((Double_t)0,gPad->GetUxmin());
        else            gxwork[0] = gPad->GetUxmin();
        npt = 2;
        for (j=first; j<=last;j++) {
           if (!OptionBins) {
              gywork[npt-1] = gywork[npt-2];
              gywork[npt]   = wmin+((j-first+1)*delta);
           }
           else {
              yj1 = y[j];      yj  = y[j-1];
              if (yj1 < yj) {
                 if (j != last) Error(where, "Y must be in increasing order");
                 else           Error(where, "Y must have N+1 values with option N");
                 return;
              }
              gywork[npt-1] = y[j-1];       gywork[npt] = y[j];
           }
           gxwork[npt-1] = x[j-1];      gxwork[npt] = x[j-1];
           if (gxwork[npt-1] >= uxmin-rounding && gxwork[npt] <= uxmax+rounding) npt += 2;
           if (j == last) {
              gywork[npt-1] = gywork[npt-2];
              gxwork[npt-1] = gxwork[0];
              ComputeLogs(npt, OptionZ);
              gPad->PaintFillArea(npt,gxworkl,gyworkl);
              if (drawborder) {
                 if (!fillarea) gyworkl[0] = ylast;
                 gPad->PaintPolyLine(npt-1,gxworkl,gyworkl);
              }
              continue;
           }
        }  //endfor (j=first; j<=last;j++)
     }
     TAttLine::Modify();
     TAttFill::Modify();
  }

//*-*-      Draw a standard Histogram (default)

  if ((OptionHist) || strlen(chopt) == 0) {
     if (!OptionRot) {
        gxwork[0] = wmin;
//      gywork[0] = TMath::Max((Double_t)0,gPad->GetUymin());
        gywork[0] = gPad->GetUymin();
        ywmin    = gywork[0];
        npt      = 2;
        for (i=first; i<=last;i++) {
           if (!OptionBins) {
              gxwork[npt-1] = gxwork[npt-2];
              gxwork[npt]   = wmin+((i-first+1)*delta);
           }
           else {
              xi1 = x[i];      xi  = x[i-1];
              if (xi1 < xi) {
                 if (i != last) Error(where, "X must be in increasing order");
                 else           Error(where, "X must have N+1 values with option N");
                 return;
              }
              gxwork[npt-1] = x[i-1];      gxwork[npt] = x[i];
           }
           gywork[npt-1] = y[i-1];
           gywork[npt]   = y[i-1];
           if (gxwork[npt-1] >= uxmin-rounding && gxwork[npt] <= uxmax+rounding) npt += 2;
           else gxwork[npt-2] = TMath::Min(gxwork[npt], uxmax);
           if (i == last) {
              gxwork[npt-1] = gxwork[npt-2];
              gywork[npt-1] = ywmin;
              ComputeLogs(npt, OptionZ);
              //gPad->PaintPolyLine(npt,gxworkl,gyworkl);
              //do not draw the two vertical lines on the edges
              Int_t npoints = npt-2;
              Int_t point1  = 1;
              if (OptionOff) {
                 // remove points before the low cutoff
                 Int_t Ip;
                 for (Ip=point1; Ip<=npoints; Ip++) {
                    if (gyworkl[Ip] != ywmin) {
                       point1 = Ip;
                       break;
                    }
                 }
                 // remove points after the high cutoff
		 Int_t point2 = npoints;
                 for (Ip=point2; Ip>=point1; Ip--) {
                    if (gyworkl[Ip] != ywmin) {
                       point2 = Ip;
                       break;
                    }
                 }
			  npoints = point2-point1+1;
		    }
              gPad->PaintPolyLine(npoints,&gxworkl[point1],&gyworkl[point1]);
              continue;
           }
        }  //endfor (i=first; i<=last;i++)
     }
     else {
        gywork[0] = wmin;
        gxwork[0] = TMath::Max((Double_t)0,gPad->GetUxmin());
        xwmin    = gxwork[0];
        npt      = 2;
        for (i=first; i<=last;i++) {
           if (!OptionBins) {
              gywork[npt-1]   = gywork[npt-2];
              gywork[npt] = wmin+((i-first+1)*delta);
           }
           else {
              yi1 = y[i];      yi  = y[i-1];
              if (yi1 < yi) {
                 if (i != last) Error(where, "Y must be in increasing order");
                 else           Error(where, "Y must have N+1 values with option N");
                 return;
              }
              gywork[npt-1] = y[i-1];      gywork[npt] = y[i];
           }
           gxwork[npt-1] = x[i-1];      gxwork[npt] = x[i-1];
           if (gxwork[npt-1] >= uxmin-rounding && gxwork[npt] <= uxmax+rounding) npt += 2;
           if (i == last) {
              gywork[npt-1] = gywork[npt-2];
              gxwork[npt-1] = xwmin;
              ComputeLogs(npt, OptionZ);
              gPad->PaintPolyLine(npt,gxworkl,gyworkl);
              continue;
           }
        }  //endfor (i=first; i<=last;i++)
     }
  }

//*-*-              Draw the histogram with a smooth Curve. The computing
//*-*-              of the smoothing is done by the routine IGRAP1

  if (OptionCurve) {
     if (!OptionFill) drawtype = 1;
     else {
        if (!OptionOne) drawtype = 2;
        else            drawtype = 3;
     }
     if (!OptionRot) {
        npt = 0;
        for (i=first; i<=last;i++) {
           npt++;
           if (!OptionBins) gxwork[npt-1] = wmin+(i-first)*delta+0.5*delta;
           else {
              xi1 = x[i];      xi  = x[i-1];
              if (xi1 < xi) {
                 if (i != last) Error(where, "X must be in increasing order");
                 else           Error(where, "X must have N+1 values with option N");
                 return;
              }
              gxwork[npt-1] = x[i-1] + 0.5*(x[i]-x[i-1]);
           }
           if (gxwork[npt-1] < uxmin || gxwork[npt-1] > uxmax) {
              npt--;
              continue;
           }
           gywork[npt-1] = y[i-1];
           ComputeLogs(npt, OptionZ);
           if ((gyworkl[npt-1] < rwymin) || (gyworkl[npt-1] > rwymax)) {
              if (npt > 2) {
                 ComputeLogs(npt, OptionZ);
                 Smooth(npt,gxworkl,gyworkl,drawtype);
              }
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
              continue;
           }
           if (npt >= 50) {
              ComputeLogs(50, OptionZ);
              Smooth(50,gxworkl,gyworkl,drawtype);
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
           }
        }  //endfor (i=first; i<=last;i++)
        if (npt > 1) {
           ComputeLogs(npt, OptionZ);
           Smooth(npt,gxworkl,gyworkl,drawtype);
        }
     }
     else {
        drawtype = drawtype+10;
        npt   = 0;
        for (i=first; i<=last;i++) {
           npt++;
           if (!OptionBins) gywork[npt-1] = wmin+(i-first)*delta+0.5*delta;
           else {
              yi1 = y[i];      yi = y[i-1];
              if (yi1 < yi) {
                 if (i != last) Error(where, "Y must be in increasing order");
                 else           Error(where, "Y must have N+1 values with option N");
                 return;
              }
              gywork[npt-1] = y[i-1] + 0.5*(y[i]-y[i-1]);
           }
           gxwork[npt-1] = x[i-1];
           ComputeLogs(npt, OptionZ);
           if ((gxworkl[npt] < uxmin) || (gxworkl[npt] > uxmax)) {
              if (npt > 2) {
                 ComputeLogs(npt, OptionZ);
                 Smooth(npt,gxworkl,gyworkl,drawtype);
              }
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
              continue;
           }
           if (npt >= 50) {
              ComputeLogs(50, OptionZ);
              Smooth(50,gxworkl,gyworkl,drawtype);
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
           }
        }  //endfor (i=first; i<=last;i++)
        if (npt > 1) {
           ComputeLogs(npt, OptionZ);
           Smooth(npt,gxworkl,gyworkl,drawtype);
        }
     }
  }

//*-*-    Draw the histogram with a simple line or/and a marker

  OptionMarker = 0;
  if ((OptionStar) || (OptionMark))OptionMarker=1;
  if ((OptionMarker) || (OptionLine)) {
    wminstep = wmin + 0.5*delta;
    Axis_t ax1,ax2,ay1,ay2;
    gPad->GetRangeAxis(ax1,ay1,ax2,ay2);

    Int_t ax1Pix = gPad->XtoAbsPixel(ax1);
    Int_t ax2Pix = gPad->XtoAbsPixel(ax2);
    Int_t ay1Pix = gPad->YtoAbsPixel(ay1);
    Int_t ay2Pix = gPad->YtoAbsPixel(ay2);

    Int_t nrPix;
    if (!OptionRot)
      nrPix = ax2Pix-ax1Pix+1;
    else
      nrPix = ay2Pix-ay1Pix+1;

    // Make here decision whether it should be painted in high or low resolution
    Int_t ip, ipix, lowRes = 0;
    if (3*nrPix < last-first+1) {
      lowRes = 1;
    }
    if (OptionFill2)       lowRes = 0;
    if (opt.Contains("9")) lowRes = 0;
    if (lowRes) {
      Double_t *minPix   = new Double_t[nrPix];
      Double_t *maxPix   = new Double_t[nrPix];
      Double_t *centrPix = new Double_t[nrPix];
      Int_t *nrEntries   = new Int_t[nrPix];

      for (ip = 0; ip < nrPix; ip++) {
        minPix[ip]    =  1e100;
        maxPix[ip]    = -1e100;
        nrEntries[ip] = 0;
      }

      for (ip = first; ip < last; ip++) {
        Double_t xw;
        if (!OptionBins) xw = wminstep + (ip-first)*delta+0.5*delta;
        else             xw = x[ip-1] + 0.5*(x[ip]-x[ip-1]);;

        if (!OptionRot) {
          Int_t ix = gPad->XtoAbsPixel(gPad->XtoPad(xw))-ax1Pix;
          if (ix < 0) ix = 0;
          if (ix >= nrPix) ix = nrPix-1;
          Int_t yPixel = gPad->YtoAbsPixel(y[ip-1]);
          if (yPixel >= ay1Pix) continue;
          if (minPix[ix] > yPixel) minPix[ix] = yPixel;
          if (maxPix[ix] < yPixel) maxPix[ix] = yPixel;
          (nrEntries[ix])++;
        } else {
          Int_t iy = gPad->YtoAbsPixel(gPad->YtoPad(y[ip-1]))-ay1Pix;
          if (iy < 0) iy = 0;
          if (iy >= nrPix) iy = nrPix-1;;
          Int_t xPixel = gPad->XtoAbsPixel(gPad->XtoPad(xw));
          if (minPix[iy] > xPixel) minPix[iy] = xPixel;
          if (maxPix[iy] < xPixel) maxPix[iy] = xPixel;
          (nrEntries[iy])++;
        }
      }

      for (ipix = 0; ipix < nrPix; ipix++) {
        if (nrEntries[ipix] > 0)
          centrPix[ipix] = (minPix[ipix]+maxPix[ipix])/2.0;
        else
          centrPix[ipix] = 2*TMath::Max(TMath::Abs(minPix[ipix]),
                                         TMath::Abs(maxPix[ipix]));
      }

      Double_t *xc = new Double_t[nrPix];
      Double_t *yc = new Double_t[nrPix];

      Double_t xcadjust = 0.3*(gPad->AbsPixeltoX(ax1Pix+1) - gPad->AbsPixeltoX(ax1Pix));
      Double_t ycadjust = 0.3*(gPad->AbsPixeltoY(ay1Pix)   - gPad->AbsPixeltoY(ay1Pix+1));
      Int_t nrLine = 0;
      for (ipix = 0; ipix < nrPix; ipix++) {
        if (minPix[ipix] <= maxPix[ipix]) {
          Double_t xl[2]; Double_t yl[2];
          if (!OptionRot) {
            xc[nrLine] = gPad->AbsPixeltoX(ax1Pix+ipix) + xcadjust;
            yc[nrLine] = gPad->AbsPixeltoY((Int_t)centrPix[ipix]);

            xl[0]      = xc[nrLine];
            yl[0]      = gPad->AbsPixeltoY((Int_t)minPix[ipix]);
            xl[1]      = xc[nrLine];
            yl[1]      = gPad->AbsPixeltoY((Int_t)maxPix[ipix]);
          } else {
            yc[nrLine] = gPad->AbsPixeltoY(ay1Pix+ipix) + ycadjust;
            xc[nrLine] = gPad->AbsPixeltoX((Int_t)centrPix[ipix]);

            xl[0]      = gPad->AbsPixeltoX((Int_t)minPix[ipix]);
            yl[0]      = yc[nrLine];
            xl[1]      = gPad->AbsPixeltoX((Int_t)maxPix[ipix]);
            yl[1]      = yc[nrLine];
          }
          if (!OptionZ && gPad->GetLogx()) {
            if (xc[nrLine] > 0) xc[nrLine] = TMath::Log10(xc[nrLine]);
            else                xc[nrLine] = gPad->GetX1();
            for (Int_t il = 0; il < 2; il++) {
              if (xl[il] > 0) xl[il] = TMath::Log10(xl[il]);
              else            xl[il] = gPad->GetX1();
            }
          }
          if (!OptionZ && gPad->GetLogy()) {
            if (yc[nrLine] > 0) yc[nrLine] = TMath::Log10(yc[nrLine]);
            else                yc[nrLine] = gPad->GetY1();
            for (Int_t il = 0; il < 2; il++) {
              if (yl[il] > 0) yl[il] = TMath::Log10(yl[il]);
              else            yl[il] = gPad->GetY1();
            }
          }

          gPad->PaintPolyLine(2,xl,yl);
          nrLine++;
        }
      }

      gPad->PaintPolyLine(nrLine,xc,yc);

      delete [] xc;
      delete [] yc;

      delete [] minPix;
      delete [] maxPix;
      delete [] centrPix;
      delete [] nrEntries;
    } else {
     if (!OptionRot) {
        npt = 0;
        for (i=first; i<=last;i++) {
           npt++;
           if (!OptionBins) gxwork[npt-1] = wmin+(i-first)*delta+0.5*delta;
           else {
              xi1 = x[i];      xi = x[i-1];
              if (xi1 < xi) {
                 if (i != last) Error(where, "X must be in increasing order");
                 else           Error(where, "X must have N+1 values with option N");
                 return;
              }
              gxwork[npt-1] = x[i-1] + 0.5*(x[i]-x[i-1]);
           }
           if (gxwork[npt-1] < uxmin || gxwork[npt-1] > uxmax) { npt--; continue;}
           if ((OptionMark != 10) && (OptionLine == 0)) {
              if (gPad->GetLogy()) {
                 if (y[i-1] <= rwymin)  {npt--; continue;}
              } else {
                 if (y[i-1] == 0)  {npt--; continue;}
              }
           }
           gywork[npt-1] = y[i-1];
           gywork[npt]   = y[i-1]; //new
           if ((gywork[npt-1] < rwymin) || (gywork[npt-1] > rwymax)) {
              if ((gywork[npt-1] < rwymin)) gywork[npt-1] = rwymin;
              if ((gywork[npt-1] > rwymax)) gywork[npt-1] = rwymax;
              if (npt > 2) {
                 if (OptionMarker) {
                    ComputeLogs(npt, OptionZ);
                    gPad->PaintPolyMarker(npt,gxworkl,gyworkl);
                 }
                 if (OptionLine) {
                    if (!OptionMarker) ComputeLogs(npt, OptionZ);
                    if (OptionFill2) {
                       gxworkl[npt]   = gxworkl[npt-1]; gyworkl[npt]   = rwymin;
                       gxworkl[npt+1] = gxworkl[0];     gyworkl[npt+1] = rwymin;
                       gPad->PaintFillArea(npt+2,gxworkl,gyworkl);
                    }
                    gPad->PaintPolyLine(npt,gxworkl,gyworkl);
                 }
              }
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
              continue;
           }
           if (npt >= 50) {
              if (OptionMarker) {
                 ComputeLogs(50, OptionZ);
                 gPad->PaintPolyMarker(50,gxworkl,gyworkl);
              }
              if (OptionLine) {
                 if (!OptionMarker) ComputeLogs(50, OptionZ);
                 if (OptionFill2) {
                    gxworkl[npt]   = gxworkl[npt-1]; gyworkl[npt]   = rwymin;
                    gxworkl[npt+1] = gxworkl[0];     gyworkl[npt+1] = rwymin;
                    gPad->PaintFillArea(52,gxworkl,gyworkl);
                 }
                 gPad->PaintPolyLine(50,gxworkl,gyworkl);
              }
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
           }
        }  //endfor (i=first; i<=last;i++)
        if (OptionMarker && npt > 0) {
           ComputeLogs(npt, OptionZ);
           gPad->PaintPolyMarker(npt,gxworkl,gyworkl);
        }
        if (OptionLine && npt > 1) {
           if (!OptionMarker) ComputeLogs(npt, OptionZ);
           if (OptionFill2) {
              gxworkl[npt]   = gxworkl[npt-1]; gyworkl[npt]   = rwymin;
              gxworkl[npt+1] = gxworkl[0];     gyworkl[npt+1] = rwymin;
              gPad->PaintFillArea(npt+2,gxworkl,gyworkl);
           }
           gPad->PaintPolyLine(npt,gxworkl,gyworkl);
        }
     }
     else {
        npt = 0;
        for (i=first; i<=last;i++) {
           npt++;
           if (!OptionBins) gywork[npt-1] = wminstep+(i-first)*delta+0.5*delta;
           else {
              yi1 = y[i];      yi = y[i-1];
              if (yi1 < yi) {
                 if (i != last) Error(where, "Y must be in increasing order");
                 else           Error(where, "Y must have N+1 values with option N");
                 return;
              }
              gywork[npt-1] = y[i-1] + 0.5*(y[i]-y[i-1]);
           }
           gxwork[npt-1] = x[i-1];
           if ((gxwork[npt-1] < uxmin) || (gxwork[npt-1] > uxmax)) {
              if (npt > 2) {
                if (OptionMarker) {
                    ComputeLogs(npt, OptionZ);
                    gPad->PaintPolyMarker(npt,gxworkl,gyworkl);
                 }
                 if (OptionLine) {
                    if (!OptionMarker) ComputeLogs(npt, OptionZ);
                    gPad->PaintPolyLine(npt,gxworkl,gyworkl);
                 }
              }
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
              continue;
           }
           if (npt >= 50) {
              if (OptionMarker) {
                 ComputeLogs(50, OptionZ);
                 gPad->PaintPolyMarker(50,gxworkl,gyworkl);
              }
              if (OptionLine) {
                 if (!OptionMarker) ComputeLogs(50, OptionZ);
                 gPad->PaintPolyLine(50,gxworkl,gyworkl);
              }
              gxwork[0] = gxwork[npt-1];
              gywork[0] = gywork[npt-1];
              npt      = 1;
           }
        }  //endfor (i=first; i<=last;i++)
        if (OptionMarker && npt > 0) {
           ComputeLogs(npt, OptionZ);
           gPad->PaintPolyMarker(npt,gxworkl,gyworkl);
        }
        if (OptionLine != 0 && npt > 1) {
           if (!OptionMarker) ComputeLogs(npt, OptionZ);
           gPad->PaintPolyLine(npt,gxworkl,gyworkl);
        }
     }
   }
  }

//*-*-              Draw the histogram as a bar chart

  if (OptionBar) {
     if (!OptionBins) { offset = delta*baroffset; dbar = delta*barwidth; }
     else {
        if (!OptionRot) {
           offset = (x[1]-x[0])*baroffset;
           dbar   = (x[1]-x[0])*barwidth;
        } else {
           offset = (y[1]-y[0])*baroffset;
           dbar   = (y[1]-y[0])*barwidth;
        }
     }
     drawbordersav = drawborder;
     gStyle->SetDrawBorder(1);
     if (!OptionRot) {
        xlow  = wmin+offset;
        xhigh = wmin+offset+dbar;
        if (!OptionOne) ylow = TMath::Max((Double_t)0,gPad->GetUymin());
        else            ylow = gPad->GetUymin();

        for (i=first; i<=last;i++) {
           yhigh    = y[i-1];
           gxwork[0] = xlow;
           gywork[0] = ylow;
           gxwork[1] = xhigh;
           gywork[1] = yhigh;
           ComputeLogs(2, OptionZ);
           gPad->PaintBox(gxworkl[0],gyworkl[0],gxworkl[1],gyworkl[1]);
           if (!OptionBins) {
              xlow  = xlow+delta;
              xhigh = xhigh+delta;
           }
           else {
              if (i < last) {
                 xi1 = x[i];      xi = x[i-1];
                 if (xi1 < xi) {
                    Error(where, "X must be in increasing order");
                    return;
                 }
                 offset  = (x[i+1]-x[i])*baroffset;
                 dbar    = (x[i+1]-x[i])*barwidth;
                 xlow    = x[i] + offset;
                 xhigh   = x[i] + offset + dbar;
              }
           }
        }  //endfor (i=first; i<=last;i++)
     }
     else {
        ylow  = wmin + offset;
        yhigh = wmin + offset + dbar;
        if (!OptionOne) xlow = TMath::Max((Double_t)0,gPad->GetUxmin());
        else            xlow = gPad->GetUxmin();

        for (i=first; i<=last;i++) {
           xhigh    = x[i-1];
           gxwork[0] = xlow;
           gywork[0] = ylow;
           gxwork[1] = xhigh;
           gywork[1] = yhigh;
           ComputeLogs(2, OptionZ);
           gPad->PaintBox(gxworkl[0],gyworkl[0],gxworkl[1],gyworkl[1]);
           gPad->PaintBox(xlow,ylow,xhigh,yhigh);
           if (!OptionBins) {
              ylow  = ylow  + delta;
              yhigh = yhigh + delta;
           }
           else {
              if (i < last) {
                 yi1 = y[i];      yi = y[i-1];
                 if (yi1 < yi) {
                    Error(where, "Y must be in increasing order");
                    return;
                 }
                 offset  = (y[i+1]-y[i])*baroffset;
                 dbar    = (y[i+1]-y[i])*barwidth;
                 ylow    = y[i] + offset;
                 yhigh   = y[i] + offset + dbar;
              }
           }
        }  //endfor (i=first; i<=last;i++)
     }
     gStyle->SetDrawBorder(drawbordersav);
  }
  gPad->ResetBit(kClipFrame);

  delete [] gxwork;
  delete [] gywork;
  delete [] gxworkl;
  delete [] gyworkl;
}


//______________________________________________________________________________
 void TGraph::ComputeLogs(Int_t npoints, Int_t opt)
{
//*-*-*-*-*-*-*-*-*-*-*-*Convert WC from Log scales*-*-*-*-*-*-*-*-*-*-*-*
//*-*                    ==========================
//
//   Take the LOG10 of gxwork and gywork according to the value of Options
//   and put it in gxworkl and gyworkl.
//
//  npoints : Number of points in gxwork and in gywork.
//

  for (Int_t i=0;i<npoints;i++) {
     gxworkl[i] = gxwork[i];
     gyworkl[i] = gywork[i];
     if (gPad->GetLogx()) {
        if (gxworkl[i] > 0) gxworkl[i] = TMath::Log10(gxworkl[i]);
        else               gxworkl[i] = gPad->GetX1();
     }
     if (!opt && gPad->GetLogy()) {
        if (gyworkl[i] > 0) gyworkl[i] = TMath::Log10(gyworkl[i]);
        else               gyworkl[i] = gPad->GetY1();
     }
//     if (gyworkl[i] > gPad->GetUymax()) gyworkl[i] = gPad->GetUymax();
  }
}


//______________________________________________________________________________
 void TGraph::Print(Option_t *) const
{
//*-*-*-*-*-*-*-*-*-*-*Print graph values*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ==================

   for (Int_t i=0;i<fNpoints;i++) {
      printf("x[%d]=%g, y[%d]=%gn",i,fX[i],i,fY[i]);
   }
}

//______________________________________________________________________________
 Int_t TGraph::RemovePoint()
{
// Delete point close to the mouse position

   Int_t px = gPad->GetEventX();
   Int_t py = gPad->GetEventY();

   //localize point to be deleted
   Int_t ipoint = -2;
   Int_t i;
   // start with a small window (in case the mouse is very close to one point)
   for (i=0;i<fNpoints;i++) {
      Int_t dpx = px - gPad->XtoAbsPixel(gPad->XtoPad(fX[i]));
      Int_t dpy = py - gPad->YtoAbsPixel(gPad->YtoPad(fY[i]));
      if (dpx*dpx+dpy*dpy < 25) {ipoint = i; break;}
   }
   if (ipoint == -2) return -1;
   fNpoints--;
   Double_t *newX = new Double_t[fNpoints];
   Double_t *newY = new Double_t[fNpoints];
   Int_t j = -1;
   for (i=0;i<fNpoints+1;i++) {
      if (i == ipoint) continue;
      j++;
      newX[j] = fX[i];
      newY[j] = fY[i];
   }
   delete [] fX;
   delete [] fY;
   fX = newX;
   fY = newY;
   gPad->Modified();
   return ipoint;
}

//______________________________________________________________________________
 Int_t TGraph::RemovePoint(Int_t ipoint)
{
// Delete 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];
   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];
   }
   delete [] fX;
   delete [] fY;
   fX = newX;
   fY = newY;
   if (gPad) gPad->Modified();
   return ipoint;
}

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

   char quote = '"';
   out<<"   "<<endl;
   if (gROOT->ClassSaved(TGraph::Class())) {
       out<<"   ";
   } else {
       out<<"   TGraph *";
   }
   out<<"graph = new TGraph("<<fNpoints<<");"<<endl;
   out<<"   graph->SetName("<<quote<<GetName()<<quote<<");"<<endl;
   out<<"   graph->SetTitle("<<quote<<GetTitle()<<quote<<");"<<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]<<");"<<endl;
   }
   if (strstr(option,"multigraph")) {
      out<<"   multigraph->Add(graph);"<<endl;
      return;
   }
   static Int_t frameNumber = 0;
   if (fHistogram) {
      frameNumber++;
      TString hname = fHistogram->GetName();
      hname += frameNumber;
      fHistogram->SetName(hname.Data());
      fHistogram->SavePrimitive(out,"nodraw");
      out<<"   graph->SetHistogram("<<fHistogram->GetName()<<");"<<endl;
      out<<"   "<<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;
      }
   }

   if (!strstr(option,"multigraph")) {
      out<<"   graph->Draw("
         <<quote<<option<<quote<<");"<<endl;
   }
}

//______________________________________________________________________________
 void TGraph::Set(Int_t n)
{
// Set number of points in the graph
// Existing coordinates are preserved
// New coordinates above fNpoints are preset to 0.

   if (n < 0) n = 0;
   if (n == fNpoints) return;
   Double_t *xx=0, *yy=0;
   if (n > 0) {
      xx = new Double_t[n];
      yy = new Double_t[n];
   }
   Int_t i;
   for (i=0; i<fNpoints && i<n;i++) {
      if (fX) xx[i] = fX[i];
      if (fY) yy[i] = fY[i];
   }
   for (i=fNpoints; i<n;i++) {
      xx[i] = 0;
      yy[i] = 0;
   }
   delete [] fX;
   delete [] fY;
   fNpoints =n;
   fX = xx;
   fY = yy;
}

//______________________________________________________________________________
 void TGraph::SetEditable(Bool_t editable)
{
// if editable=kFALSE, the graph cannot be modified with the mouse
//  by default a TGraph is editable

   if (editable) ResetBit(kNotEditable);
   else          SetBit(kNotEditable);
}

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

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

//______________________________________________________________________________
 void TGraph::SetPoint(Int_t i, Double_t x, Double_t y)
{
//*-*-*-*-*-*-*-*-*-*-*Set x and y values for point number i*-*-*-*-*-*-*-*-*
//*-*                  =====================================

   if (i < 0) return;
   if (i >= fNpoints) {
   // re-allocate the object
      Double_t *savex = new Double_t[i+1];
      Double_t *savey = new Double_t[i+1];
      if (fNpoints > 0) {
         memcpy(savex,fX,fNpoints*sizeof(Double_t));
         memcpy(savey,fY,fNpoints*sizeof(Double_t));
      }
      if (fX) delete [] fX;
      if (fY) delete [] fY;
      fX = savex;
      fY = savey;
      fNpoints = i+1;
   }
   fX[i] = x;
   fY[i] = y;
   if (fHistogram) {
      delete fHistogram;
      fHistogram = 0;
   }
}

//______________________________________________________________________________
 void TGraph::SetTitle(const char* title)
{
  fTitle = title;
  if (fHistogram) fHistogram->SetTitle(title);
}

//______________________________________________________________________________
 void TGraph::Smooth(Int_t npoints, Double_t *x, Double_t *y, Int_t drawtype)
{
//*-*-*-*-*-*-*-*-*-*-*-*Smooth a curve given by N points*-*-*-*-*-*-*-*-*-*
//*-*                    ================================
//
//   Underlaying routine for Draw based on the CERN GD3 routine TVIPTE
//
//     Author - Marlow etc.   Modified by - P. Ward     Date -  3.10.1973
//
//   This routine draws a smooth tangentially continuous curve through
// the sequence of data points P(I) I=1,N where P(I)=(X(I),Y(I))
// the curve is approximated by a polygonal arc of short vectors .
// the data points can represent open curves, P(1) != P(N) or closed
// curves P(2) == P(N) . If a tangential discontinuity at P(I) is
// required , then set P(I)=P(I+1) . loops are also allowed .
//
// Reference Marlow and Powell,Harwell report No.R.7092.1972
// MCCONALOGUE,Computer Journal VOL.13,NO4,NOV1970PP392 6
//
// _Input parameters:
//
//  npoints   : Number of data points.
//  x         : Abscissa
//  y         : Ordinate
//
//
// delta is the accuracy required in constructing the curve.
// if it is zero then the routine calculates a value other-
// wise it uses this value. (default is 0.0)
//
//

   Int_t i, k, kp, km, NpointsMax, banksize, n2, npt;
   Int_t maxiterations, finished;
   Int_t jtype, ktype, closed;
   Double_t sxmin, sxmax, symin, symax;
   Double_t delta;
   Double_t xorg, yorg;
   Double_t ratio_signs, xratio, yratio;
   Int_t flgic, flgis;
   Int_t iw, loptx;
   Double_t P1, P2, P3, P4, P5, P6;
   Double_t W1, W2, W3;
   Double_t A, B, C, R, S, T, Z;
   Double_t CO, SO, CT, ST, CTU, STU, XNT;
   Double_t DX1, DY1, DX2, DY2, DK1, DK2;
   Double_t XO, YO, DX, DY, XT, YT;
   Double_t XA, XB, YA, YB;
   Double_t U1, U2, U3, TJ;
   Double_t CC, ERR;
   Double_t SB, STH;
   Double_t wsign, tsquare, tcube;
   C = T = CO = SO = CT = ST = CTU = STU = DX1 = DY1 = DX2 = DY2 = 0;
   XT = YT = XA = XB = YA = YB = U1 = U2 = U3 = TJ = SB = 0;
//*-* ______________________________________

   NpointsMax  = npoints*10;
   n2          = NpointsMax-2;
   banksize    = n2;

   Double_t *qlx = new Double_t[NpointsMax];
   Double_t *qly = new Double_t[NpointsMax];
   if (!qlx || !qly) {
      Error("Smooth", "not enough space in memory");
      return;
   }

//*-*-  Decode the type of curve according to
//*-*-  chopt of IGHIST.
//*-*-  ('S', 'SA', 'SA1' ,'XS', 'XSA', or 'XSA1')

   //if (drawtype >= 1000) drawtype -= 1000;
   loptx = kFALSE;
   jtype  = (drawtype%1000)-10;
   if (jtype > 0) { ktype = jtype; loptx = kTRUE; }
   else             ktype = drawtype%1000;

   Double_t ruxmin = gPad->GetUxmin();
   Double_t ruymin = gPad->GetUymin();
   if (ktype == 3) {
      xorg = ruxmin;
      yorg = ruymin;
   } else {
     xorg = TMath::Max((Double_t)0,ruxmin);
     yorg = TMath::Max((Double_t)0,ruymin);
   }

   maxiterations = 20;
   delta         = 0.00055;

//*-*-       Scale data to the range 0-ratio_signs in X, 0-1 in Y
//*-*-       where ratio_signs is the ratio between the number of changes
//*-*-       of sign in Y divided by the number of changes of sign in X

   sxmin = x[0];
   sxmax = x[0];
   symin = y[0];
   symax = y[0];
   Double_t six   = 1;
   Double_t siy   = 1;
   for (i=1;i<npoints;i++) {
      if (i > 1) {
         if ((x[i]-x[i-1])*(x[i-1]-x[i-2]) < 0) six++;
         if ((y[i]-y[i-1])*(y[i-1]-y[i-2]) < 0) siy++;
      }
      if (x[i] < sxmin) sxmin = x[i];
      if (x[i] > sxmax) sxmax = x[i];
      if (y[i] < symin) symin = y[i];
      if (y[i] > symax) symax = y[i];
   }
   closed = 0;
   Double_t dx1n   = TMath::Abs(x[npoints-1]-x[0]);
   Double_t dy1n   = TMath::Abs(y[npoints-1]-y[0]);
   if (dx1n < 0.01*(sxmax-sxmin) && dy1n < 0.01*(symax-symin))  closed = 1;
   if (sxmin == sxmax) xratio = 1;
   else {
      if (six > 1) ratio_signs = siy/six;
      else         ratio_signs = 20;
      xratio = ratio_signs/(sxmax-sxmin);
   }
   if (symin == symax) yratio = 1;
   else                yratio = 1/(symax-symin);

   qlx[0] = x[0];
   qly[0] = y[0];
   for (i=0;i<npoints;i++) {
      x[i] = (x[i]-sxmin)*xratio;
      y[i] = (y[i]-symin)*yratio;
   }

//*-*-           finished is minus one if we must draw a straight line from P(k-1)
//*-*-           to P(k). finished is one if the last call to IPL has  <  N2
//*-*-           points. finished is zero otherwise. npt counts the X and Y
//*-*-           coordinates in work . When npt=N2 a call to IPL is made.

   finished = 0;
   npt      = 1;
   k        = 1;

//*-*-           Convert coordinates back to original system

//*-*-           Separate the set of data points into arcs P(k-1),P(k).
//*-*-           Calculate the direction cosines. first consider whether
//*-*-           there is a continuous tangent at the endpoints.

   if (!closed) {
      if (x[0] != x[npoints-1] || y[0] != y[npoints-1]) goto L40;
      if (x[npoints-2] == x[npoints-1] && y[npoints-2] == y[npoints-1]) goto L40;
      if (x[0] == x[1] && y[0] == y[1]) goto L40;
   }
   flgic = kFALSE;
   flgis = kTRUE;

//*-*-           flgic is true if the curve is open and false if it is closed.
//*-*-           flgis is true in the main loop, but is false if there is
//*-*-           a deviation from the main loop.

   km = npoints - 1;

//*-*-           Calculate direction cosines at P(1) using P(N-1),P(1),P(2).

    goto L100;
L40:
   flgic = kTRUE;
   flgis = kFALSE;

//*-*-           Skip excessive consecutive equal points.

L50:
   if (k >= npoints) {
      finished = 1;  //*-*-  Prepare to clear out remaining short vectors before returning
      if (npt > 1) goto L310;
      goto L390;
   }
   k++;
   if (x[k-1] == x[k-2] && y[k-1] == y[k-2])  goto L50;
L60:
   km = k-1;
   if (k > npoints) {
      finished = 1;  //*-*-  Prepare to clear out remaining short vectors before returning
      if (npt > 1) goto L310;
      goto L390;
   }
   if (k < npoints) goto L90;
   if (!flgic) { kp = 2; goto L130;}

L80:
   if (flgis) goto L150;

//*-*-           Draw a straight line from P(k-1) to P(k).

   finished = -1;
   goto L170;

//*-*-           Test whether P(k) is a cusp.

L90:
   if (x[k-1] == x[k] && y[k-1] == y[k]) goto L80;
L100:
   kp = k+1;
   goto L130;

//*-*-           Branch if the next section of the curve begins at a cusp.

L110:
   if (!flgis) goto L50;

//*-*-           Carry forward the direction cosines from the previous arc.

L120:
   CO = CT;
   SO = ST;
   k++;
   goto L60;

//*-*-           Calculate the direction cosines at P(k).  If k=1 then
//*-*-           N-1 is used for k-1. If k=N then 2 is used for k+1.
//*-*-           direction cosines at P(k) obtained from P(k-1),P(k),P(k+1).

L130:
   DX1 = x[k-1]  - x[km-1];
   DY1 = y[k-1]  - y[km-1];
   DK1 = DX1*DX1 + DY1*DY1;
   DX2 = x[kp-1] - x[k-1];
   DY2 = y[kp-1] - y[k-1];
   DK2 = DX2*DX2 + DY2*DY2;
   CTU = DX1*DK2 + DX2*DK1;
   STU = DY1*DK2 + DY2*DK1;
   XNT = CTU*CTU + STU*STU;

//*-*-           If both ctu and stu are zero,then default.This can
//*-*-           occur when P(k)=P(k+1). I.E. A loop.

   if (XNT < 1.E-25) {
      CTU = DY1;
      STU =-DX1;
      XNT = DK1;
   }
//*-*-           Normalise direction cosines.

   CT = CTU/TMath::Sqrt(XNT);
   ST = STU/TMath::Sqrt(XNT);
   if (flgis) goto L160;

//*-*-           Direction cosines at P(k-1) obtained from P(k-1),P(k),P(k+1).

   W3    = 2*(DX1*DY2-DX2*DY1);
   CO    = CTU+W3*DY1;
   SO    = STU-W3*DX1;
   XNT   = 1/TMath::Sqrt(CO*CO+SO*SO);
   CO    = CO*XNT;
   SO    = SO*XNT;
   flgis = kTRUE;
   goto L170;

//*-*-           Direction cosines at P(k) obtained from P(k-2),P(k-1),P(k).

L150:
   W3    = 2*(DX1*DY2-DX2*DY1);
   CT    = CTU-W3*DY2;
   ST    = STU+W3*DX2;
   XNT   = 1/TMath::Sqrt(CT*CT+ST*ST);
   CT    = CT*XNT;
   ST    = ST*XNT;
   flgis = kFALSE;
   goto L170;
L160:
   if (k <= 1) goto L120;

//*-*-           For the arc between P(k-1) and P(k) with direction cosines CO,
//*-*-           SO and CT,ST respectively, calculate the coefficients of the
//*-*-           parametric cubic represented by X(T) and Y(T) where
//*-*-           X(T)=XA*T**3 + XB*T**2 + CO*T + XO
//*-*-           Y(T)=YA*T**3 + YB*T**2 + SO*T + YO

L170:
   XO = x[k-2];
   YO = y[k-2];
   DX = x[k-1] - XO;
   DY = y[k-1] - YO;

//*-*-           Initialise the values of X(TI),Y(TI) in XT and YT respectively.

   XT = XO;
   YT = YO;
   if (finished < 0) {  //*-*- Draw a straight line between (XO,YO) and (XT,YT)
      XT += DX;
      YT += DY;
      goto L300;
   }
   C  = DX*DX+DY*DY;
   A  = CO+CT;
   B  = SO+ST;
   R  = DX*A+DY*B;
   T  = C*6/(TMath::Sqrt(R*R+2*(7-CO*CT-SO*ST)*C)+R);
   tsquare = T*T;
   tcube   = T*tsquare;
   XA = (A*T-2*DX)/tcube;
   XB = (3*DX-(CO+A)*T)/tsquare;
   YA = (B*T-2*DY)/tcube;
   YB = (3*DY-(SO+B)*T)/tsquare;

//*-*-           If the curve is close to a straight line then use a straight
//*-*-           line between (XO,YO) and (XT,YT).

   if (.75*TMath::Max(TMath::Abs(DX*SO-DY*CO),TMath::Abs(DX*ST-DY*CT)) <= delta) {
      finished = -1;
      XT += DX;
      YT += DY;
      goto L300;
   }

//*-*-           Calculate a set of values 0 == T(0).LTCT(1) <  ...  < T(M)=TC
//*-*-           such that polygonal arc joining X(T(J)),Y(T(J)) (J=0,1,..M)
//*-*-           is within the required accuracy of the curve

   TJ = 0;
   U1 = YA*XB-YB*XA;
   U2 = YB*CO-XB*SO;
   U3 = SO*XA-YA*CO;

//*-*-           Given T(J), calculate T(J+1). The values of X(T(J)),
//*-*-           Y(T(J)) T(J) are contained in XT,YT and TJ respectively.

L180:
   S  = T - TJ;
   iw = -2;

//*-*-           Define iw here later.

   P1 = (2*U1)*TJ-U3;
   P2 = (U1*TJ-U3)*3*TJ+U2;
   P3 = 3*TJ*YA+YB;
   P4 = (P3+YB)*TJ+SO;
   P5 = 3*TJ*XA+XB;
   P6 = (P5+XB)*TJ+CO;

//*-*-           Test D(TJ,THETA). A is set to (Y(TJ+S)-Y(TJ))/S.B is
//*-*-           set to (X(TJ+S)-X(TJ))/S.

   CC  = 0.8209285;
   ERR = 0.1209835;
L190:
   iw -= 2;
L200:
   A   = (S*YA+P3)*S+P4;
   B   = (S*XA+P5)*S+P6;

//*-*-           Set Z to PSI(D/delta)-CC.

   W1 = -S*(S*U1+P1);
   W2 = S*S*U1-P2;
   W3 = 1.5*W1+W2;

//*-*-           Set the estimate of (THETA-TJ)/S.Then set the numerator
//*-*-           of the expression (EQUATION 4.4)/S. Then set the square
//*-*-           of D(TJ,TJ+S)/delta. Then replace Z by PSI(D/delta)-CC.

   if (W3 > 0) wsign = TMath::Abs(W1);
   else        wsign = -TMath::Abs(W1);
   STH = 0.5+wsign/(3.4*TMath::Abs(W1)+5.2*TMath::Abs(W3));
   Z   = S*STH*(S-S*STH)*(W1*STH+W1+W2);
   Z   = Z*Z/((A*A+B*B)*(delta*delta));
   Z   = (Z+2.642937)*Z/((.3715652*Z+3.063444)*Z+.2441889)-CC;

//*-*-           Branch if Z has been calculated

   if (iw > 0) goto L250;
   if (Z > ERR) goto L240;
   goto L220;
L210:
   iw -= 2;
L220:
   if (iw+2 == 0) goto L190;
   if (iw+2 >  0) goto L290;

//*-*-           Last part of arc.

L230:
   XT = x[k-1];
   YT = y[k-1];
   S  = 0;
   goto L300;

//*-*-           Z(S). find a value of S where 0 <= S <= SB such that
//*-*-           TMath::Abs(Z(S)) < ERR

L240:
   kp = 0;
   C  = Z;
   SB = S;
L250:
   Zero(kp,0,SB,ERR,S,Z,maxiterations);
   if (kp == 2) goto L210;
   if (kp > 2) {
      Error("Smooth", "Attempt to plot outside plot limits");
      goto L230;
   }
   if (iw > 0) goto L200;

//*-*-           Set Z=Z(S) for S=0.

  if (iw < 0) {
     Z  = -CC;
     iw = 0;
     goto L250;
  }

//*-*-           Set Z=Z(S) for S=SB.

   Z  = C;
   iw = 1;
   goto L250;

//*-*-           Update TJ,XT and YT.

L290:
   XT = XT + S*B;
   YT = YT + S*A;
   TJ = S  + TJ;

//*-*-           Convert coordinates to original system

L300:
   qlx[npt] = sxmin + XT/xratio;
   qly[npt] = symin + YT/yratio;
   npt++;

//*-*-           If a fill area must be drawn and if the banks LX and
//*-*-           LY are too small they are enlarged in order to draw
//*-*-           the filled area in one go.

   if (npt < banksize)  goto L320;
   if (drawtype >= 1000 || ktype > 1) {
      Int_t newsize = banksize + n2;
      Double_t *qtemp = new Double_t[banksize];
      for (i=0;i<banksize;i++) qtemp[i] = qlx[i];
      delete [] qlx;
      qlx = new Double_t[newsize];
      for (i=0;i<banksize;i++) qlx[i]   = qtemp[i];
      for (i=0;i<banksize;i++) qtemp[i] = qly[i];
      delete [] qly;
      qly = new Double_t[newsize];
      for (i=0;i<banksize;i++) qly[i] = qtemp[i];
      delete [] qtemp;
      banksize = newsize;
      goto L320;
   }

//*-*-           Draw the graph

L310:
   if (drawtype >= 1000) {
      gPad->PaintFillArea(npt,qlx,qly, "B");
   }
   else {
      if (ktype > 1) {
         if (!loptx) {
            qlx[npt]   = qlx[npt-1];
            qlx[npt+1] = qlx[0];
            qly[npt]   = yorg;
            qly[npt+1] = yorg;
         }
         else {
            qlx[npt]   = xorg;
            qlx[npt+1] = xorg;
            qly[npt]   = qly[npt-1];
            qly[npt+1] = qly[0];
         }
         gPad->PaintFillArea(npt+2,qlx,qly);
      }
      gPad->PaintPolyLine(npt,qlx,qly);
   }
   npt = 1;
   qlx[0] = sxmin + XT/xratio;
   qly[0] = symin + YT/yratio;
L320:
   if (finished > 0) goto L390;
   if (finished < 0) { finished = 0; goto L110;}
   if (S > 0) goto L180;
   goto L110;



//*-*-           Convert coordinates back to original system

L390:
    for (i=0;i<npoints;i++) {
      x[i] = sxmin + x[i]/xratio;
      y[i] = symin + y[i]/yratio;
    }

   delete [] qlx;
   delete [] qly;
}

//______________________________________________________________________________
 void TGraph::Sort(Bool_t (*greater)(const TGraph*, Int_t, Int_t) /*=TGraph::CompareX()*/, 
		  Bool_t ascending /*=kTRUE*/, Int_t low /* =0 */, Int_t high /* =-1111 */) {
// Sorts the points of this TGraph using in-place quicksort (see e.g. older glibc).
// To compare two points the function parameter greater is used (see TGraph::CompareX for an 
// example of such a method, which is also the default comparison function for Sort). After 
// the sort, greater(this, i, j) will return kTRUE for all i>j if ascending == kTRUE, and 
// kFALSE otherwise.
//
// The last two parameters are used for the recursive quick sort, stating the range to be sorted
//
// Examples:
//   // sort points along x axis
//   graph->Sort(); 
//   // sort points along their distance to origin
//   graph->Sort(&TGraph::CompareRadius); 
//   
//   Bool_t CompareErrors(const TGraph* gr, Int_t i, Int_t j) {
//     const TGraphErrors* ge=(const TGraphErrors*)gr;
//     return (ge->GetEY()[i]>ge->GetEY()[j]); }
//   // sort using the above comparison function, largest errors first
//   graph->Sort(&CompareErrors, kFALSE);

   if (high == -1111) high = GetN()-1;
   //  Termination condition
   if (high <= low) return;

   int left, right;
   left = low; // low is the pivot element
   right = high;
   while (left < right) {
      // move left while item < pivot 
      while(left <= high && greater(this, left, low) != ascending) 
	 left++;
      // move right while item > pivot
      while(right > low && greater(this, right, low) == ascending) 
	 right--;
      if (left < right && left < high && right > low) 
	 SwapPoints(left, right);
   }
   // right is final position for the pivot
   if (right > low) 
      SwapPoints(low, right);
   Sort( greater, ascending, low, right-1 );
   Sort( greater, ascending, right+1, high );
}

//______________________________________________________________________________
 void TGraph::Streamer(TBuffer &b)
{
   // Stream an object of class TGraph.

   if (b.IsReading()) {
      UInt_t R__s, R__c;
      Version_t R__v = b.ReadVersion(&R__s, &R__c);
      if (R__v > 2) {
         TGraph::Class()->ReadBuffer(b, this, R__v, R__s, R__c);
         if (fHistogram) fHistogram->SetDirectory(0);
         return;
      }
      //====process old versions before automatic schema evolution
      TNamed::Streamer(b);
      TAttLine::Streamer(b);
      TAttFill::Streamer(b);
      TAttMarker::Streamer(b);
      b >> fNpoints;
      fX = new Double_t[fNpoints];
      fY = new Double_t[fNpoints];
      if (R__v < 2) {
         Float_t *x = new Float_t[fNpoints];
         Float_t *y = new Float_t[fNpoints];
         b.ReadFastArray(x,fNpoints);
         b.ReadFastArray(y,fNpoints);
         for (Int_t i=0;i<fNpoints;i++) {
            fX[i] = x[i];
            fY[i] = y[i];
         }
         delete [] y;
         delete [] x;
      } else {
         b.ReadFastArray(fX,fNpoints);
         b.ReadFastArray(fY,fNpoints);
      }
      b >> fFunctions;
      b >> fHistogram;
      if (fHistogram) fHistogram->SetDirectory(0);
      if (R__v < 2) {
         Float_t mi,ma;
         b >> mi;
         b >> ma;
         fMinimum = mi;
         fMaximum = ma;
      } else {
         b >> fMinimum;
         b >> fMaximum;
      }
      b.CheckByteCount(R__s, R__c, TGraph::IsA());
      //====end of old versions

   } else {
      TGraph::Class()->WriteBuffer(b,this);
   }
}

//______________________________________________________________________________
 void TGraph::SwapPoints(Int_t pos1, Int_t pos2) {
   SwapValues(fX, pos1, pos2);
   SwapValues(fY, pos1, pos2);
}

//______________________________________________________________________________
 void TGraph::SwapValues(Double_t* arr, Int_t pos1, Int_t pos2) {
   Double_t tmp=arr[pos1];
   arr[pos1]=arr[pos2];
   arr[pos2]=tmp;
}

//______________________________________________________________________________
 void TGraph::UseCurrentStyle()
{
   // Set current style settings in this graph
   // This function is called when either TCanvas::UseCurrentStyle
   // or TROOT::ForceStyle have been invoked.

   SetFillColor(gStyle->GetHistFillColor());
   SetFillStyle(gStyle->GetHistFillStyle());
   SetLineColor(gStyle->GetHistLineColor());
   SetLineStyle(gStyle->GetHistLineStyle());
   SetLineWidth(gStyle->GetHistLineWidth());
   SetMarkerColor(gStyle->GetMarkerColor());
   SetMarkerStyle(gStyle->GetMarkerStyle());
   SetMarkerSize(gStyle->GetMarkerSize());
   if (fHistogram) fHistogram->UseCurrentStyle();
}

//______________________________________________________________________________
 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*-*-*-*-*-*-*-*-*-*
//*-*                    ==================================
//
//  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;
}


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.