Diff of /trunk/hist/hist/src/TMultiGraph.cxx

-revision 10847, Mon Dec 27 15:42:36 2004 UTC
+revision 11226, Fri Mar  4 09:06:37 2005 UTC
 Line 1
- // @(#)root/graf:$Name:  $:$Id: TMultiGraph.cxx,v 1.14 2004/06/18 10:46:58 brun Exp $
+ // @(#)root/graf:$Name:  $:$Id: TMultiGraph.cxx,v 1.15 2004/12/27 15:42:36 brun Exp $
  // Author: Rene Brun   12/10/2000
  /*************************************************************************
 Line 15
  #include "TH1.h"
  #include "TVirtualPad.h"
  #include "Riostream.h"
+ #include "TVirtualFitter.h"
  #include <ctype.h>
+ extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);
  ClassImp(TMultiGraph)
-Line 47
+Line 50
  // TMultiGraph default constructor
     fGraphs    = 0;
+    fFunctions = 0;
     fHistogram = 0;
     fMaximum   = -1111;
     fMinimum   = -1111;
-Line 58
+Line 62
  {
  // constructor with name and title
     fGraphs    = 0;
+    fFunctions = 0;
     fHistogram = 0;
     fMaximum   = -1111;
     fMinimum   = -1111;
-Line 68
+Line 73
  {
  // TMultiGraph destructor
     if (!fGraphs) return;
     TGraph *g;
     TIter   next(fGraphs);
-Line 80
+Line 84
     fGraphs = 0;
     delete fHistogram;
     fHistogram = 0;
+    if (fFunctions) {
+       fFunctions->SetBit(kInvalidObject);
+       //special logic to support the case where the same object is
+       //added multiple times in fFunctions.
+       //This case happens when the same object is added with different
+       //drawing modes
+       TObject *obj;
+       while ((obj  = fFunctions->First())) {
+          while(fFunctions->Remove(obj));
+          delete obj;
+       }
+       delete fFunctions;
+    }
  }
  //______________________________________________________________________________
-Line 147
+Line 164
  }
  //______________________________________________________________________________
+ Int_t TMultiGraph::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...
+    char *linear;
+    linear=strstr(fname, "++");
+    TF1 *f1=0;
+    if (linear)
+          f1=new TF1(fname, fname, xmin, xmax);
+    else {
+       f1 = (TF1*)gROOT->GetFunction(fname);
+       if (!f1) { Printf("Unknown function: %s",fname); return -1; }
+    }
+    return Fit(f1,option,"",xmin,xmax);
+ }
+ //______________________________________________________________________________
+ Int_t TMultiGraph::Fit(TF1 *f1, Option_t *option, Option_t *, Axis_t rxmin, Axis_t rxmax)
+ {
+ //*-*-*-*-*-*-*-*-*-*-*Fit this multigraph with function f1*-*-*-*-*-*-*-*-*-*
+ //*-*                  ==================================
+ //
+ //   In this function all graphs of the multigraph are fitted simultaneously
+ //
+ //   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)
+ //             = "B" Use this option when you want to fix one or more parameters
+ //                   and the fitting function is like "gaus","expo","poln","landau".
+ //             = "R" Use the Range specified in the function range
+ //             = "N" Do not store the graphics function, do not draw
+ //             = "0" Do not plot the result of the fit. By default the fitted function
+ //                   is drawn unless the option"N" above is specified.
+ //             = "+" Add this new fitted function to the list of fitted functions
+ //                   (by default, any previous function is deleted)
+ //             = "C" In case of linear fitting, not calculate the chisquare
+ //                    (saves time)
+ //
+ //   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");
+ //
+ //
+ //   who is calling this function
+ //   ============================
+ //   Note that this function is called when calling TGraphErrors::Fit
+ //   or TGraphAsymmErrors::Fit ot TGraphBentErrors::Fit
+ //   see the discussion below on the errors calulation.
+ //
+ //   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, landau. 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);
+ //
+ //  How errors are used in the chisquare function (see TFitter GraphFitChisquare)//   Access to the fit results
+ //   ============================================
+ // 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-exlow and
+ // x+exhigh.
+ //
+ // The chisquare is computed as the sum of the quantity below at each point:
+ //
+ //                     (y - f(x))**2
+ //         -----------------------------------
+ //         ey**2 + ((f(x+exhigh) - f(x-exlow))/2)**2
+ //
+ // where x and y are the point coordinates.
+ //
+ // In case the function lies below (above) the data point, ey is ey_low (ey_high).
+ //
+ //  thanks to Andy Haas (haas@yahoo.com) for adding the case with TGraphasymmerrors
+ //            University of Washington
+ //
+ // a little different approach to approximating the uncertainty in y because of the
+ // errors in x, is to make it equal the error in x times the slope of the line.
+ // The improvement, compared to the first method (f(x+ exhigh) - f(x-exlow))/2
+ // is of (error of x)**2 order. This approach is called "effective variance method".
+ // This improvement has been made in version 4.00/08 by Anna Kreshuk.
+ //
+ //   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");
+ //
+ //  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;
+    Int_t np;
+    TF1 *fnew1;
+    // Check validity of function
+    if (!f1) {
+       Error("Fit", "function may not be null pointer");
+       return 0;
+    }
+    if (f1->IsZombie()) {
+       Error("Fit", "function is zombie");
+       return 0;
+    }
+    npar = f1->GetNpar();
+    if (npar <= 0) {
+       Error("Fit", "function %s has illegal number of parameters = %d", f1->GetName(), npar);
+       return 0;
+    }
+    // Check that function has same dimension as graph
+    if (f1->GetNdim() > 1) {
+       Error("Fit", "function %s is not 1-D", f1->GetName());
+       return 0;
+    }
+    TGraph *g;
+    TIter next(fGraphs);
+    Double_t *arglist = new Double_t[100];
+    // Decode string choptin and fill fitOption structure
+    Foption_t fitOption;
+    fitOption.Quiet   = 0;
+    fitOption.Verbose = 0;
+    fitOption.Bound   = 0;
+    fitOption.Like    = 0;
+    fitOption.W1      = 0;
+    fitOption.Errors  = 0;
+    fitOption.Range   = 0;
+    fitOption.Gradient= 0;
+    fitOption.Nograph = 0;
+    fitOption.Nostore = 0;
+    fitOption.Plus    = 0;
+    fitOption.User    = 0;
+    fitOption.Nochisq = 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;
+    if (opt.Contains("C")) fitOption.Nochisq = 1;
+    if (rxmax > rxmin) {
+       xmin = rxmin;
+       xmax = rxmax;
+    } else {
+       g=(TGraph *)fGraphs->First();
+       if (!g) {
+          Error("Fit", "No graphs in the multigraph");
+          return 0;
+       }
+       Double_t *px, *py;
+       np=g->GetN();
+       px=g->GetX();
+       py=g->GetY();
+       xmin=px[0];
+       xmax=py[np-1];
+       ymin=px[0];
+       ymax=py[np-1];
+       Double_t err0=g->GetErrorX(0);
+       Double_t errn=g->GetErrorX(np-1);
+       if (err0 > 0) xmin -= 2*err0;
+       if (errn > 0) xmax += 2*errn;
+       next.Reset();
+       while ((g = (TGraph*) next())) {
+          np=g->GetN();
+          px=g->GetX();
+          py=g->GetY();
+          for (i=0; i<np; i++) {
+             if (px[i] < xmin) xmin = px[i];
+             if (px[i] > xmax) xmax = px[i];
+             if (py[i] < ymin) ymin = py[i];
+             if (py[i] > ymax) ymax = py[i];
+          }
+       }
+    }
+    ///////////////
+    //set the fitter
+    //////////////
+    //TClass *cl=gROOT->GetClass("TLinearFitter");
+    //
+    Int_t special=f1->GetNumber();
+    Bool_t linear = f1->IsLinear();
+    if (special==299+npar)
+       linear=kTRUE;
+    char l[]="TLinearFitter";
+    Int_t strdiff = 0;
+    Bool_t IsSet = kFALSE;
+    if (TVirtualFitter::GetFitter()){
+       //Is a fitter already set? Is it linear?
+       IsSet = kTRUE;
+       strdiff = strcmp(TVirtualFitter::GetFitter()->IsA()->GetName(), l);
+    }
+    if (linear){
+       //
+       TClass *cl = gROOT->GetClass("TLinearFitter");
+       if (IsSet && strdiff!=0) {
+          delete TVirtualFitter::GetFitter();
+          IsSet=kFALSE;
+       }
+       if (!IsSet) {
+         //TLinearFitter *lf=(TLinearFitter *)cl->New();
+          TVirtualFitter::SetFitter((TVirtualFitter *)cl->New());
+       }
+    } else {
+       if (IsSet && strdiff==0){
+          delete TVirtualFitter::GetFitter();
+          IsSet=kFALSE;
+       }
+       if (!IsSet)
+          TVirtualFitter::SetFitter(0);
+    }
+    TVirtualFitter *grFitter = TVirtualFitter::Fitter(this, f1->GetNpar());
+    grFitter->Clear();
+ //*-*- Get pointer to the function by searching in the list of functions in ROOT
+    grFitter->SetUserFunc(f1);
+    grFitter->SetFitOption(fitOption);
+ //*-*- Is a Fit range specified?
+    if (fitOption.Range) {
+       f1->GetRange(xmin, xmax);
+    } else {
+       f1->SetRange(xmin, xmax);
+    }
+    if (linear && !fitOption.Bound && !fitOption.Like && !fitOption.Errors){
+       grFitter->ExecuteCommand("FitMultiGraph", 0, 0);
+    } else {
+       //Int_t special = f1->GetNumber();
+       if (fitOption.Bound) special = 0;
+       if      (special == 100)      InitGaus(xmin,xmax);
+       else if (special == 400)      InitGaus(xmin,xmax);
+       else if (special == 200)      InitExpo(xmin,xmax);
+       else if (special == 299+npar) InitPolynom(xmin,xmax);
+       //*-*- 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] = TVirtualFitter::GetErrorDef();
+       if (!fitOption.User) grFitter->SetFitMethod("MultiGraphFitChisquare");
+       fitResult = grFitter->ExecuteCommand("SET ERR",arglist,1);
+       if (fitResult != 0) {
+          //   Abnormal termination, MIGRAD might not have converged on a
+          //   minimum.
+          if (!fitOption.Quiet) {
+             Warning("Fit","Abnormal termination of minimization.");
+          }
+          delete [] arglist;
+          return fitResult;
+       }
+       //*-*- Transfer names and initial values of parameters to Minuit
+       Int_t nfixed = 0;
+       for (i=0;i<npar;i++) {
+          par = f1->GetParameter(i);
+          f1->GetParLimits(i,al,bl);
+          if (al*bl != 0 && al >= bl) {
+             al = bl = 0;
+             arglist[nfixed] = i+1;
+             nfixed++;
+          }
+          we  = 0.3*TMath::Abs(par);
+          if (we <= TMath::Abs(par)*1e-6) we = 1;
+          grFitter->SetParameter(i,f1->GetParName(i),par,we,al,bl);
+       }
+       if(nfixed > 0)grFitter->ExecuteCommand("FIX",arglist,nfixed); // Otto
+       //*-*- Reset Print level
+       if (!fitOption.Quiet) {
+          if (fitOption.Verbose) { arglist[0] = 2; grFitter->ExecuteCommand("SET PRINT", arglist,1); }
+          else                   { arglist[0] = 0; grFitter->ExecuteCommand("SET PRINT", arglist,1); }
+       }
+       //*-*- Compute sum of squares of errors in the bin range
+       Bool_t hasErrors = kFALSE;
+       Double_t ex, ey, sumw2=0;
+       next.Reset();
+       while ((g = (TGraph*) next())) {
+          np=g->GetN();
+          for (i=0; i<np; i++){
+             ex=g->GetErrorX(i);
+             ey=g->GetErrorY(i);
+             if (ex > 0 || ey > 0) hasErrors=kTRUE;
+             sumw2+=ey*ey;
+          }
+       }
+       //*-*- Perform minimization
+       arglist[0] = TVirtualFitter::GetMaxIterations();
+       arglist[1] = sumw2*TVirtualFitter::GetPrecision();
+       grFitter->ExecuteCommand("MIGRAD",arglist,2);
+       if (fitOption.Errors) {
+          grFitter->ExecuteCommand("HESSE",arglist,0);
+          grFitter->ExecuteCommand("MINOS",arglist,0);
+       }
+       grFitter->GetStats(amin,edm,errdef,nvpar,nparx);
+       f1->SetChisquare(amin);
+       Int_t ndf = f1->GetNumberFitPoints()-npar+nfixed;
+       f1->SetNDF(ndf);
+       //*-*- Get return status
+       char parName[50];
+       for (i=0;i<npar;i++) {
+          grFitter->GetParameter(i,parName, par,we,al,bl);
+          if (!fitOption.Errors) werr = we;
+          else {
+             grFitter->GetErrors(i,eplus,eminus,eparab,globcc);
+             if (eplus > 0 && eminus < 0) werr = 0.5*(eplus-eminus);
+             else                         werr = we;
+          }
+          if (!hasErrors && ndf > 1) werr *= TMath::Sqrt(amin/(ndf-1));
+          f1->SetParameter(i,par);
+          f1->SetParError(i,werr);
+       }
+    }
+ //*-*- Print final values of parameters.
+    if (!fitOption.Quiet) {
+       if (fitOption.Errors) grFitter->PrintResults(4,amin);
+       else                  grFitter->PrintResults(3,amin);
+    }
+    delete [] arglist;
+ //*-*- Store fitted function in histogram functions list and draw
+    if (!fitOption.Nostore) {
+       if (!fFunctions) fFunctions = new TList;
+       if (!fitOption.Plus) {
+          TIter next2(fFunctions, kIterBackward);
+          TObject *obj;
+          while ((obj = next2())) {
+             if (obj->InheritsFrom(TF1::Class())){
+                obj = fFunctions->Remove(obj);
+                delete obj;
+             }
+          }
+       }
+       fnew1 = new TF1();
+       f1->Copy(*fnew1);
+       fFunctions->Add(fnew1);
+       fnew1->SetParent(this);
+       fnew1->Save(xmin,xmax,0,0,0,0);
+       if (fitOption.Nograph) fnew1->SetBit(TF1::kNotDraw);
+       fnew1->SetBit(TFormula::kNotGlobal);
+       if (TestBit(kCanDelete)) return fitResult;
+       if (gPad) gPad->Modified();
+    }
+    return fitResult;
+ }
+ //______________________________________________________________________________
+ void TMultiGraph::InitGaus(Double_t xmin, Double_t xmax)
+ {
+ //*-*-*-*-*-*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
+    Int_t np = 0;
+    allcha = sumx = sumx2 = 0;
+    TGraph *g;
+    TIter next(fGraphs);
+    Double_t *px, *py;
+    Int_t npp; //number of points in each graph
+    while ((g = (TGraph*) next())) {
+       px=g->GetX();
+       py=g->GetY();
+       npp=g->GetN();
+       for (bin=0; bin<npp; bin++){
+          x=px[bin];
+          if (x<xmin || x>xmax) continue;
+          np++;
+          val=py[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((xmax-xmin)/np);
+    if (rms == 0) rms = 1;
+    TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
+    TF1 *f1 = (TF1*)grFitter->GetUserFunc();
+    f1->SetParameter(0,binwidx*allcha/(sqrtpi*rms));
+    f1->SetParameter(1,mean);
+    f1->SetParameter(2,rms);
+    f1->SetParLimits(2,0,10*rms);
+ }
+ //______________________________________________________________________________
+ void TMultiGraph::InitExpo(Double_t xmin, Double_t xmax)
+ {
+ //*-*-*-*-*-*Compute Initial values of parameters for an exponential*-*-*-*-*
+ //*-*        =======================================================
+    Double_t constant, slope;
+    Int_t ifail;
+    LeastSquareLinearFit(-1, constant, slope, ifail, xmin, xmax);
+    TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
+    TF1 *f1 = (TF1*)grFitter->GetUserFunc();
+    f1->SetParameter(0,constant);
+    f1->SetParameter(1,slope);
+ }
+ //______________________________________________________________________________
+ void TMultiGraph::InitPolynom(Double_t xmin, Double_t xmax)
+ {
+ //*-*-*-*-*-*Compute Initial values of parameters for a polynom*-*-*-*-*-*-*
+ //*-*        ===================================================
+    Double_t fitpar[25];
+    TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
+    TF1 *f1 = (TF1*)grFitter->GetUserFunc();
+    Int_t npar   = f1->GetNpar();
+    LeastSquareFit(npar, fitpar, xmin, xmax);
+    for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
+ }
+ //______________________________________________________________________________
+ void TMultiGraph::LeastSquareFit(Int_t m, Double_t *a, Double_t xmin, Double_t xmax)
+ {
+ //*-*-*-*-*-*-*-*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, bin;
+     Double_t power;
+     Double_t da[20], xk, yk;
+     //count the total number of points to fit
+     TGraph *g;
+     TIter next(fGraphs);
+     Double_t *px, *py;
+     Int_t n=0;
+     Int_t npp;
+     while ((g = (TGraph*) next())) {
+        px=g->GetX();
+        py=g->GetY();
+        npp=g->GetN();
+        for (bin=0; bin<npp; bin++){
+           xk=px[bin];
+           if (xk < xmin || xk > xmax) continue;
+           n++;
+        }
+     }
+     if (m <= 2) {
+        LeastSquareLinearFit(n, a[0], a[1], ifail, xmin, xmax);
+        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;
+     next.Reset();
+     while ((g = (TGraph*) next())) {
+        px=g->GetX();
+        py=g->GetY();
+        npp=g->GetN();
+        for (k = 0; k <= npp; ++k) {
+           xk     = px[k];
+           if (xk < xmin || xk > xmax) continue;
+           np++;
+           yk     = py[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];
+        py=((TGraph *)fGraphs->First())->GetY();
+        a[0]=py[0];
+        for (i=1; i<m; ++i) a[i] = 0;
+        return;
+     }
+     for (i=0; i<m; ++i) a[i] = da[i];
+ }
+ //______________________________________________________________________________
+ void TMultiGraph::LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail, Double_t xmin, Double_t xmax)
+ {
+ //*-*-*-*-*-*-*-*-*-*Least square linear fit without weights*-*-*-*-*-*-*-*-*
+ //*-*                =======================================
+ //
+ //  Fit a straight line (a0 + a1*x) to the data in this graph.
+ //  ndata:  number of points to fit
+ //  first:  first point number to fit
+ //  last:   last point to fit O(ndata should be last-first
+ //  ifail:  return parameter indicating the status of the fit (ifail=0, fit is OK)
+ //
+ //   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;
+     TGraph *g;
+     TIter next(fGraphs);
+     Double_t *px, *py;
+     Int_t npp;
+     while ((g = (TGraph*) next())) {
+        px=g->GetX();
+        py=g->GetY();
+        npp=g->GetN();
+        for (i = 0; i < npp; ++i) {
+           xk = px[i];
+           if (xk < xmin || xk > xmax) continue;
+           np++;
+           yk = py[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;
+ }
+ //______________________________________________________________________________
  TH1F *TMultiGraph::GetHistogram() const
  {
  //    Returns a pointer to the histogram used to draw the axis
-Line 164
+Line 877
  }
  //______________________________________________________________________________
+ TF1 *TMultiGraph::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);
+ }
+ //______________________________________________________________________________
  TAxis *TMultiGraph::GetXaxis() const
  {
     // Get x axis of the graph.
-Line 336
+Line 1062
           lnk = (TObjOptLink*)lnk->Next();
        }
     }
+    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();
+       }
+      }
+    }
  }
  //______________________________________________________________________________

 Legend:



Removed from v.10847
 


changed lines


 
Added in v.11226
 Legend:



Removed from v.10847
 


changed lines


 
Added in v.11226
-Removed from v.10847
+Added in v.11226

Subversion Admin	ViewVC Help
Powered by ViewVC 1.0.9