// @(#)root/graf:$Name: $:$Id: TGraph.cxx,v 1.83 2002/09/15 19:45:15 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"
TVirtualFitter *grFitter;
TF1 *grF1 = 0;
Foption_t fitOption;
Double_t *gxwork, *gywork, *gxworkl, *gyworkl;
extern void GraphFitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
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());
SetMaximum(h->GetMaximumStored());
SetMinimum(h->GetMinimumStored());
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 *option)
: 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();
}
//______________________________________________________________________________
void TGraph::ComputeRange(Double_t &, Double_t &, Double_t &, Double_t &)
{
// 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
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
grFitter = TVirtualFitter::Fitter(this);
grFitter->Clear();
//*-*- Get pointer to the function by searching in the list of functions in ROOT
grF1 = f1;
if (!grF1) { Printf("Function is a null pointer"); return 0; }
npar = grF1->GetNpar();
if (npar <=0) { Printf("Illegal number of parameters = %d",npar); return 0; }
//*-*- Check that function has same dimension as histogram
if (grF1->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) {
grF1->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 {
grF1->SetRange(xmin, xmax);
gxfirst = 0;
gxlast = fNpoints-1;
}
//*-*- If case of a predefined function, then compute initial values of parameters
Int_t special = grF1->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->SetFCN(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 = grF1->GetParameter(i);
grF1->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,grF1->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);
grF1->SetChisquare(amin);
Int_t ndf = grF1->GetNumberFitPoints()-npar+nfixed;
grF1->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));
grF1->SetParameter(i,par);
grF1->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();
grF1->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;
if (!gPad) return 0;
gPad->Modified();
gPad->Update();
if (fHistogram) return fHistogram;
TH1F *h1 = (TH1F*)gPad->FindObject("hframe");
return h1;
}
//______________________________________________________________________________
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 GraphFitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
//*-*-*-*-*-*Minimization function for Graphs using a Chisquare method*-*-*-*-*
//*-* =========================================================
//
// In case of a TGraphErrors object, ex, the error along x, is projected
// along the y-direction by calculating the function at the points x-ex and
// x+ex.
//
// The chisquare is computed as the sum of the quantity below at each point:
//
// (y - f(x))**2
// -----------------------------------
// ey**2 + ((f(x+ex) - f(x-ex))/2)**2
//
// where x and y are the point coordinates
Double_t cu,eu,ex,ey,eux,fu,fsum,fm,fp;
Double_t x[1], xx[1];
Double_t xm,xp;
Int_t bin, npfits=0;
TGraph *gr = (TGraph*)grFitter->GetObjectFit();
Int_t n = gr->GetN();
Double_t *gx = gr->GetX();
Double_t *gy = gr->GetY();
Double_t fxmin = grF1->GetXmin();
Double_t fxmax = grF1->GetXmax();
npar = grF1->GetNpar();
grF1->InitArgs(x,u);
f = 0;
for (bin=0;bin<n;bin++) {
x[0] = gx[bin];
if (!grF1->IsInside(x)) continue;
cu = gy[bin];
TF1::RejectPoint(kFALSE);
fu = grF1->EvalPar(x,u);
if (TF1::RejectedPoint()) continue;
fsum = (cu-fu);
npfits++;
if (fitOption.W1) {
f += fsum*fsum;
continue;
}
ex = gr->GetErrorX(bin);
ey = gr->GetErrorY(bin);
if (ex < 0) ex = 0;
if (ey < 0) ey = 0;
if (ex >= 0) {
xm = x[0] - ex; if (xm < fxmin) xm = fxmin;
xp = x[0] + ex; if (xp > fxmax) xp = fxmax;
xx[0] = xm; fm = grF1->EvalPar(xx,u);
xx[0] = xp; fp = grF1->EvalPar(xx,u);
eux = 0.5*(fp-fm);
} else
eux = 0.;
eu = ey*ey+eux*eux;
if (eu <= 0) eu = 1;
f += fsum*fsum/eu;
}
grF1->SetNumberFitPoints(npfits);
}
//______________________________________________________________________________
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;
grF1->SetParameter(0,binwidx*allcha/(sqrtpi*rms));
grF1->SetParameter(1,mean);
grF1->SetParameter(2,rms);
grF1->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);
grF1->SetParameter(0,constant);
grF1->SetParameter(1,slope);
}
//______________________________________________________________________________
void TGraph::InitPolynom(Int_t first, Int_t last)
{
//*-*-*-*-*-*Compute Initial values of parameters for a polynom*-*-*-*-*-*-*
//*-* ===================================================
Double_t fitpar[25];
Int_t npar = grF1->GetNpar();
LeastSquareFit(npar, fitpar, first, last);
for (Int_t i=0;i<npar;i++) grF1->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 beetwen 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.
//
//
const char *where = "PaintGraphHist";
Int_t OptionLine , OptionAxis , OptionCurve, OptionStar , OptionMark;
Int_t OptionBar , OptionRot , OptionOne;
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("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
gPad->PaintPolyLine(npt-2,&gxworkl[1],&gyworkl[1]);
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;
else xw = 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]);
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]))-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] = wminstep + (i-first)*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];
}
if (gxwork[npt-1] < uxmin || gxwork[npt-1] > uxmax) { npt--; continue;}
if ((OptionMark != 10) && (OptionLine == 0)) {
if (gPad->GetLogy()) {
if (y[i-1] < 0) {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;
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;
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;
if (fHistogram) fHistogram->SetMaximum(maximum);
}
//______________________________________________________________________________
void TGraph::SetMinimum(Double_t minimum)
{
fMinimum = minimum;
if (fHistogram) fHistogram->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-10;
if (jtype > 0) { ktype = jtype; loptx = kTRUE; }
else ktype = drawtype;
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::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::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.