<commit_before>// @(#)root/hist:$Name:  $:$Id: TF1.cxx,v 1.84 2004/07/23 14:04:19 brun Exp $
// Author: <NAME>   18/08/95

/*************************************************************************
 * Copyright (C) 1995-2000, <NAME> and <NAME>.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include "Riostream.h"
#include "TROOT.h"
#include "TMath.h"
#include "TF1.h"
#include "TH1.h"
#include "TVirtualPad.h"
#include "TStyle.h"
#include "TRandom.h"
#include "Api.h"
#include "TPluginManager.h"
#include "TVirtualUtilPad.h"
#include "TBrowser.h"

Bool_t TF1::fgAbsValue    = kFALSE;
Bool_t TF1::fgRejectPoint = kFALSE;
static TF1 *gHelper = 0;
static Double_t ErrorTF1 = 0;

ClassImp(TF1)

//______________________________________________________________________________
//
// a TF1 object is a 1-Dim function defined between a lower and upper limit.
// The function may be a simple function (see TFormula) or a precompiled
// user function.
// The function may have associated parameters.
// TF1 graphics function is via the TH1/TGraph drawing functions.
//
//  The following types of functions can be created:
//    A- Expression using variable x and no parameters
//    B- Expression using variable x with parameters
//    C- A general C function with parameters
//
//      Example of a function of type A
//
//   TF1 *f1 = new TF1("f1","sin(x)/x",0,10);
//   f1->Draw();
//Begin_Html
/*
<img src="gif/function1.gif">
*/
//End_Html
//
//      Example of a function of type B
//   TF1 *f1 = new TF1("f1","[0]*x*sin([1]*x)",-3,3);
//    This creates a function of variable x with 2 parameters.
//    The parameters must be initialized via:
//      f1->SetParameter(0,value_first_parameter);
//      f1->SetParameter(1,value_second_parameter);
//    Parameters may be given a name:
//      f1->SetParName(0,"Constant");
//
//     Example of function of type C
//   Consider the macro myfunc.C below
//-------------macro myfunc.C-----------------------------
//Double_t myfunction(Double_t *x, Double_t *par)
//{
//   Float_t xx =x[0];
//   Double_t f = TMath::Abs(par[0]*sin(par[1]*xx)/xx);
//   return f;
//}
//void myfunc()
//{
//   TF1 *f1 = new TF1("myfunc",myfunction,0,10,2);
//   f1->SetParameters(2,1);
//   f1->SetParNames("constant","coefficient");
//   f1->Draw();
//}
//void myfit()
//{
//   TH1F *h1=new TH1F("h1","test",100,0,10);
//   h1->FillRandom("myfunc",20000);
//   TF1 *f1=gROOT->GetFunction("myfunc");
//   f1->SetParameters(800,1);
//   h1.Fit("myfunc");
//}
//--------end of macro myfunc.C---------------------------------
//
// In an interactive session you can do:
//   Root > .L myfunc.C
//   Root > myfunc();
//   Root > myfit();
//
//
//  TF1 objects can reference other TF1 objects (thanks <NAME>)
//  of type A or B defined above.This excludes CINT interpreted functions
//  and compiled functions.
//  However, there is a restriction. A function cannot reference a basic
//  function if the basic function is a polynomial polN.
//  Example:
//{
//  TF1 *fcos = new TF1 ("fcos", "[0]*cos(x)", 0., 10.);
//  fcos->SetParNames( "cos");
//  fcos->SetParameter( 0, 1.1);
//
//  TF1 *fsin = new TF1 ("fsin", "[0]*sin(x)", 0., 10.);
//  fsin->SetParNames( "sin");
//  fsin->SetParameter( 0, 2.1);
//
//  TF1 *fsincos = new TF1 ("fsc", "fcos+fsin");
//
//  TF1 *fs2 = new TF1 ("fs2", "fsc+fsc");
//}
//
//     WHY TF1 CANNOT ACCEPT A CLASS MEMBER FUNCTION ?
//     ===============================================
// This is a frequently asked question.
// C++ is a strongly typed language. There is no way for TF1 (without
// recompiling this class) to know about all possible user defined data types.
// This also apply to the case of a static class function.
//
//------------------------------------------------------------------------

TF1 *TF1::fgCurrent = 0;

//______________________________________________________________________________
TF1::TF1(): TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*-*-*-*-*F1 default constructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ======================

   fXmin      = 0;
   fXmax      = 0;
   fNpx       = 100;
   fType      = 0;
   fNpfits    = 0;
   fNDF       = 0;
   fNsave     = 0;
   fChisquare = 0;
   fIntegral  = 0;
   fFunction  = 0;
   fParErrors = 0;
   fParMin    = 0;
   fParMax    = 0;
   fAlpha     = 0;
   fBeta      = 0;
   fGamma     = 0;
   fParent    = 0;
   fSave      = 0;
   fHistogram = 0;
   fMinimum   = -1111;
   fMaximum   = -1111;
   fMethodCall = 0;
   SetFillStyle(0);
}


//______________________________________________________________________________
TF1::TF1(const char *name,const char *formula, Double_t xmin, Double_t xmax)
      :TFormula(name,formula), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using a formula definition*-*-*-*-*-*-*-*-*-*-*
//*-*          =========================================
//*-*
//*-*  See TFormula constructor for explanation of the formula syntax.
//*-*
//*-*  See tutorials: fillrandom, first, fit1, formula1, multifit
//*-*  for real examples.
//*-*
//*-*  Creates a function of type A or B between xmin and xmax
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   if (xmin < xmax ) {
      fXmin      = xmin;
      fXmax      = xmax;
   } else {
      fXmin = xmax; //when called from TF2,TF3
      fXmax = xmin;
   }
   fNpx       = 100;
   fType      = 0;
   fFunction  = 0;
   if (fNpar) {
      fParErrors = new Double_t[fNpar];
      fParMin    = new Double_t[fNpar];
      fParMax    = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fNsave      = 0;
   fSave       = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fMethodCall = 0;
   
   if (fNdim != 1 && xmin < xmax) {
      Error("TF1","function: %s/%s has %d parameters instead of 1",name,formula,fNdim);
      MakeZombie();
   }

   if (!gStyle) return;
   SetLineColor(gStyle->GetFuncColor());
   SetLineWidth(gStyle->GetFuncWidth());
   SetLineStyle(gStyle->GetFuncStyle());
   SetFillStyle(0);
}


//______________________________________________________________________________
TF1::TF1(const char *name, Double_t xmin, Double_t xmax, Int_t npar)
      :TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using name of an interpreted function*-*-*-*
//*-*          =======================================================
//*-*
//*-*  Creates a function of type C between xmin and xmax.
//*-*  name is the name of an interpreted CINT cunction.
//*-*  The function is defined with npar parameters
//*-*  fcn must be a function of type:
//*-*     Double_t fcn(Double_t *x, Double_t *params)
//*-*
//*-*  This constructor is called for functions of type C by CINT.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   fXmin       = xmin;
   fXmax       = xmax;
   fNpx        = 100;
   fType       = 2;
   fFunction   = 0;
   if (npar > 0 ) fNpar = npar;
   if (fNpar) {
      fNames      = new TString[fNpar];
      fParams     = new Double_t[fNpar];
      fParErrors  = new Double_t[fNpar];
      fParMin     = new Double_t[fNpar];
      fParMax     = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParams[i]     = 0;
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fNsave      = 0;
   fSave       = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fMethodCall = 0;
   fNdim       = 1;

   TF1 *f1old = (TF1*)gROOT->GetListOfFunctions()->FindObject(name);
   if (f1old) delete f1old;
   SetName(name);

   if (gStyle) {
      SetLineColor(gStyle->GetFuncColor());
      SetLineWidth(gStyle->GetFuncWidth());
      SetLineStyle(gStyle->GetFuncStyle());
   }
   SetFillStyle(0);

   SetTitle(name);
   if (name) {
      if (*name == '*') return; //case happens via SavePrimitive
      fMethodCall = new TMethodCall();
      fMethodCall->InitWithPrototype(name,"Double_t*,Double_t*");
      fNumber = -1;
      gROOT->GetListOfFunctions()->Add(this);
   } else {
      Printf("Function:%s cannot be compiled",name);
   }
}

//______________________________________________________________________________
TF1::TF1(const char *name,void *fcn, Double_t xmin, Double_t xmax, Int_t npar)
      :TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using pointer to an interpreted function*-*-*-*
//*-*          =======================================================
//*-*
//*-*  See TFormula constructor for explanation of the formula syntax.
//*-*
//*-*  Creates a function of type C between xmin and xmax.
//*-*  The function is defined with npar parameters
//*-*  fcn must be a function of type:
//*-*     Double_t fcn(Double_t *x, Double_t *params)
//*-*
//*-*  see tutorial; myfit for an example of use
//*-*  also test/stress.cxx (see function stress1)
//*-*
//*-*
//*-*  This constructor is called for functions of type C by CINT.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   fXmin       = xmin;
   fXmax       = xmax;
   fNpx        = 100;
   fType       = 2;
   fFunction   = 0;
   if (npar > 0 ) fNpar = npar;
   if (fNpar) {
      fNames      = new TString[fNpar];
      fParams     = new Double_t[fNpar];
      fParErrors  = new Double_t[fNpar];
      fParMin     = new Double_t[fNpar];
      fParMax     = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParams[i]     = 0;
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fNsave      = 0;
   fSave       = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fMethodCall = 0;
   fNdim       = 1;

   TF1 *f1old = (TF1*)gROOT->GetListOfFunctions()->FindObject(name);
   if (f1old) delete f1old;
   SetName(name);

   if (gStyle) {
      SetLineColor(gStyle->GetFuncColor());
      SetLineWidth(gStyle->GetFuncWidth());
      SetLineStyle(gStyle->GetFuncStyle());
   }
   SetFillStyle(0);

   if (!fcn) return;
   char *funcname = G__p2f2funcname(fcn);
   SetTitle(funcname);
   if (funcname) {
      fMethodCall = new TMethodCall();
      fMethodCall->InitWithPrototype(funcname,"Double_t*,Double_t*");
      fNumber = -1;
      gROOT->GetListOfFunctions()->Add(this);
   } else {
      Printf("Function:%s cannot be compiled",name);
   }
}

//______________________________________________________________________________
TF1::TF1(const char *name,Double_t (*fcn)(Double_t *, Double_t *), Double_t xmin, Double_t xmax, Int_t npar)
      :TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using a pointer to real function*-*-*-*-*-*-*-*
//*-*          ===============================================
//*-*
//*-*   npar is the number of free parameters used by the function
//*-*
//*-*   This constructor creates a function of type C when invoked
//*-*   with the normal C++ compiler.
//*-*
//*-*   see test program test/stress.cxx (function stress1) for an example.
//*-*   note the interface with an intermediate pointer.
//*-*
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   fXmin       = xmin;
   fXmax       = xmax;
   fNpx        = 100;
   char *funcname = G__p2f2funcname((void*)fcn);
   if (funcname) {
      fType       = 2;
      SetTitle(funcname);
      fMethodCall = new TMethodCall();
      fMethodCall->InitWithPrototype(funcname,"Double_t*,Double_t*");
      fNumber = -1;
      fFunction   = 0;
   } else {
      fType       = 1;
      fMethodCall = 0;
      fFunction   = fcn;
   }
   if (npar > 0 ) fNpar = npar;
   if (fNpar) {
      fNames      = new TString[fNpar];
      fParams     = new Double_t[fNpar];
      fParErrors  = new Double_t[fNpar];
      fParMin     = new Double_t[fNpar];
      fParMax     = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParams[i]     = 0;
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fNsave      = 0;
   fSave       = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fNdim       = 1;
//*-*- Store formula in linked list of formula in ROOT

   TF1 *f1old = (TF1*)gROOT->GetListOfFunctions()->FindObject(name);
   if (f1old) delete f1old;
   SetName(name);
   gROOT->GetListOfFunctions()->Add(this);

   if (!gStyle) return;
   SetLineColor(gStyle->GetFuncColor());
   SetLineWidth(gStyle->GetFuncWidth());
   SetLineStyle(gStyle->GetFuncStyle());
   SetFillStyle(0);

}

//______________________________________________________________________________
TF1& TF1::operator=(const TF1 &rhs) 
{
   if (this != &rhs) {
      rhs.Copy(*this);
   }
   return *this;
}

//______________________________________________________________________________
TF1::~TF1()
{
//*-*-*-*-*-*-*-*-*-*-*F1 default destructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =====================

   if (fParMin)    delete [] fParMin;
   if (fParMax)    delete [] fParMax;
   if (fParErrors) delete [] fParErrors;
   if (fIntegral)  delete [] fIntegral;
   if (fAlpha)     delete [] fAlpha;
   if (fBeta)      delete [] fBeta;
   if (fGamma)     delete [] fGamma;
   if (fSave)      delete [] fSave;
   delete fHistogram;
   delete fMethodCall;

   if (fParent) {
      if (fParent->InheritsFrom(TH1::Class())) {
         ((TH1*)fParent)->GetListOfFunctions()->Remove(this);
         return;
      }
      //parent may be a graph, or ??
      //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->RemoveObject(fParent,this);
      fParent = 0;
   }
}

//______________________________________________________________________________
TF1::TF1(const TF1 &f1) : TFormula(), TAttLine(f1), TAttFill(f1), TAttMarker(f1)
{

   fXmin      = 0;
   fXmax      = 0;
   fNpx       = 100;
   fType      = 0;
   fNpfits    = 0;
   fNDF       = 0;
   fNsave     = 0;
   fChisquare = 0;
   fIntegral  = 0;
   fFunction  = 0;
   fParErrors = 0;
   fParMin    = 0;
   fParMax    = 0;
   fAlpha     = 0;
   fBeta      = 0;
   fGamma     = 0;
   fParent    = 0;
   fSave      = 0;
   fHistogram = 0;
   fMinimum   = -1111;
   fMaximum   = -1111;
   fMethodCall = 0;
   SetFillStyle(0);

   ((TF1&)f1).Copy(*this);
}

//______________________________________________________________________________
void TF1::AbsValue(Bool_t flag)
{
   // static function: set the fgAbsValue flag.
   // By default TF1::Integral uses the original function value to compute the integral
   // However, TF1::Moment, CentralMoment require to compute the integral
   // using the absolute value of the function.
   
   fgAbsValue = flag;
}

//______________________________________________________________________________
void TF1::Browse(TBrowser *b)
{
    Draw(b ? b->GetDrawOption() : "");
    gPad->Update();
}


//______________________________________________________________________________
void TF1::Copy(TObject &obj) const
{
//*-*-*-*-*-*-*-*-*-*-*Copy this F1 to a new F1*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ========================

   if (((TF1&)obj).fParMin)    delete [] ((TF1&)obj).fParMin;
   if (((TF1&)obj).fParMax)    delete [] ((TF1&)obj).fParMax;
   if (((TF1&)obj).fParErrors) delete [] ((TF1&)obj).fParErrors;
   if (((TF1&)obj).fIntegral)  delete [] ((TF1&)obj).fIntegral;
   if (((TF1&)obj).fAlpha)     delete [] ((TF1&)obj).fAlpha;
   if (((TF1&)obj).fBeta)      delete [] ((TF1&)obj).fBeta;
   if (((TF1&)obj).fGamma)     delete [] ((TF1&)obj).fGamma;
   if (((TF1&)obj).fSave)      delete [] ((TF1&)obj).fSave;
   delete ((TF1&)obj).fHistogram;
   delete ((TF1&)obj).fMethodCall;

   TFormula::Copy(obj);
   TAttLine::Copy((TF1&)obj);
   TAttFill::Copy((TF1&)obj);
   TAttMarker::Copy((TF1&)obj);
   ((TF1&)obj).fXmin = fXmin;
   ((TF1&)obj).fXmax = fXmax;
   ((TF1&)obj).fNpx  = fNpx;
   ((TF1&)obj).fType = fType;
   ((TF1&)obj).fFunction  = fFunction;
   ((TF1&)obj).fChisquare = fChisquare;
   ((TF1&)obj).fNpfits  = fNpfits;
   ((TF1&)obj).fNDF     = fNDF;
   ((TF1&)obj).fMinimum = fMinimum;
   ((TF1&)obj).fMaximum = fMaximum;

   ((TF1&)obj).fParErrors = 0;
   ((TF1&)obj).fParMin    = 0;
   ((TF1&)obj).fParMax    = 0;
   ((TF1&)obj).fIntegral  = 0;
   ((TF1&)obj).fAlpha     = 0;
   ((TF1&)obj).fBeta      = 0;
   ((TF1&)obj).fGamma     = 0;
   ((TF1&)obj).fParent    = fParent;
   ((TF1&)obj).fNsave     = fNsave;
   ((TF1&)obj).fSave      = 0;
   ((TF1&)obj).fHistogram = 0;
   ((TF1&)obj).fMethodCall = 0;
   if (fNsave) {
      ((TF1&)obj).fSave = new Double_t[fNsave];
      for (Int_t j=0;j<fNsave;j++) ((TF1&)obj).fSave[j] = fSave[j];
   }
   if (fNpar) {
      ((TF1&)obj).fParErrors = new Double_t[fNpar];
      ((TF1&)obj).fParMin    = new Double_t[fNpar];
      ((TF1&)obj).fParMax    = new Double_t[fNpar];
      Int_t i;
      for (i=0;i<fNpar;i++)   ((TF1&)obj).fParErrors[i] = fParErrors[i];
      for (i=0;i<fNpar;i++)   ((TF1&)obj).fParMin[i]    = fParMin[i];
      for (i=0;i<fNpar;i++)   ((TF1&)obj).fParMax[i]    = fParMax[i];
   }
   if (fMethodCall) {
      TMethodCall *m = new TMethodCall();
      m->InitWithPrototype(fMethodCall->GetMethodName(),fMethodCall->GetProto());
      ((TF1&)obj).fMethodCall  = m;
   }
}

//______________________________________________________________________________
Double_t TF1::Derivative(Double_t x, Double_t *params, Double_t eps) const
{
  // returns the first derivative of the function at point x, 
  // computed by Richardson's extrapolation method (use 2 derivative estimates 
  // to compute a third, more accurate estimation)
  // first, derivatives with steps h and h/2 are computed by central difference formulas
  // D(h) = (f(x+h) - f(x-h))/2h
  // the final estimate D = (4*D(h/2) - D(h))/3
  //
  // if the argument params is null, the current function parameters are used,
  // otherwise the parameters in params are used.
  //
  // the argument eps may be specified to control the step size (precision).
  // the step size is taken as eps*(xmax-xmin).
  // the default value (0.001) should be good enough for the vast majority
  // of functions. Give a smaller value if your function has many changes
  // of the second derivative in the function range.
  //
  // Getting the error via TF1::DerivativeError
  // -----------------
  //   (total error = roundoff error + interpolation error)
  // the estimate of the roundoff error is taken as follows:
  //    err = k*Sqrt(f(x)*f(x) + x*x*deriv*deriv)*Sqrt(Sum(ai)*(ai)),
  // where k is the double precision, ai are coefficients used in
  // central difference formulas
  // interpolation error is decreased by making the step size h smaller.
  //
  // Author: <NAME>
  
  const Double_t C1 = 1e-15;
   
  if(eps< 1e-10 || eps > 1e-2) {
     Warning("Derivative","parameter esp=%g out of allowed range[1e-10,1e-2], reset to 0.001",eps);
     eps = 0.001;
  }
  Double_t xmin, xmax;
  GetRange(xmin, xmax);
  Double_t h = eps*(xmax-xmin);

  Double_t xx[1];
  TF1 *func = (TF1*)this;
  func->InitArgs(xx, params);
  xx[0] = x+h;     Double_t f1 = func->EvalPar(xx, params);
  xx[0] = x;       Double_t fx = func->EvalPar(xx, params);
  xx[0] = x-h;     Double_t f2 = func->EvalPar(xx, params);

  xx[0] = x+h/2;   Double_t g1 = func->EvalPar(xx, params);
  xx[0] = x-h/2;   Double_t g2 = func->EvalPar(xx, params);

  //compute the central differences
  Double_t h2    = 1/(2.*h);
  Double_t D0    = f1 - f2;
  Double_t D2    = 2*(g1 - g2);
  ErrorTF1       = C1*h2*fx;  //compute the error
  Double_t deriv = h2*(4*D2 - D0)/3.;  
  return deriv;
}


//______________________________________________________________________________
Double_t TF1::Derivative2(Double_t x, Double_t *params, Double_t eps) const
{
  // returns the first derivative of the function at point x, 
  // computed by Richardson's extrapolation method (use 2 derivative estimates 
  // to compute a third, more accurate estimation)
  // first, derivatives with steps h and h/2 are computed by central difference formulas
  //    D(h) = (f(x+h) - 2*f(x) + f(x-h))/(h*h)
  //    the final estimate D = (4*D(h/2) - D(h))/3
  //
  // if the argument params is null, the current function parameters are used,
  // otherwise the parameters in params are used.
  //
  // the argument eps may be specified to control the step size (precision).
  // the step size is taken as eps*(xmax-xmin).
  // the default value (0.001) should be good enough for the vast majority
  // of functions. Give a smaller value if your function has many changes
  // of the second derivative in the function range.
  //
  // Getting the error via TF1::DerivativeError
  // -----------------
  //   (total error = roundoff error + interpolation error)
  // the estimate of the roundoff error is taken as follows:
  //    err = k*Sqrt(f(x)*f(x) + x*x*deriv*deriv)*Sqrt(Sum(ai)*(ai)),
  // where k is the double precision, ai are coefficients used in
  // central difference formulas
  // interpolation error is decreased by making the step size h smaller.
  //
  // Author: <NAME>

  const Double_t C1 = 2*1e-15;
   
  if(eps< 1e-6 || eps > 1e-2) {
     Warning("Derivative2","parameter esp=%g out of allowed range[1e-6,1e-2], reset to 0.001",eps);
     eps = 0.001;
  }
  Double_t xmin, xmax;
  GetRange(xmin, xmax);
  Double_t h = eps*(xmax-xmin);

  Double_t xx[1];
  TF1 *func = (TF1*)this;
  func->InitArgs(xx, params);
  xx[0] = x+h;     Double_t f1 = func->EvalPar(xx, params);
  xx[0] = x;       Double_t f2 = func->EvalPar(xx, params);
  xx[0] = x-h;     Double_t f3 = func->EvalPar(xx, params);

  xx[0] = x+h/2;   Double_t g1 = func->EvalPar(xx, params);
  xx[0] = x-h/2;   Double_t g3 = func->EvalPar(xx, params);

  //compute the central differences
  Double_t hh    = 1/(h*h);
  Double_t D0    = f3 - 2*f2 + f1;
  Double_t D2    = 4*g3 - 8*f2 +4*g1;
  ErrorTF1       = C1*hh*f2;  //compute the error
  Double_t deriv = hh*(4*D2 - D0)/3.;
  return deriv;
}


//______________________________________________________________________________
Double_t TF1::Derivative3(Double_t x, Double_t *params, Double_t eps) const
{
  // returns the first derivative of the function at point x, 
  // computed by Richardson's extrapolation method (use 2 derivative estimates 
  // to compute a third, more accurate estimation)
  // first, derivatives with steps h and h/2 are computed by central difference formulas
  //    D(h) = (f(x+2h) - 2*f(x+h) + 2*f(x-h) - f(x-2h))/(2*h*h*h)
  //    the final estimate D = (4*D(h/2) - D(h))/3
  //
  // if the argument params is null, the current function parameters are used,
  // otherwise the parameters in params are used.
  //
  // the argument eps may be specified to control the step size (precision).
  // the step size is taken as eps*(xmax-xmin).
  // the default value (0.001) should be good enough for the vast majority
  // of functions. Give a smaller value if your function has many changes
  // of the second derivative in the function range.
  //
  // Getting the error via TF1::DerivativeError
  // -----------------
  //   (total error = roundoff error + interpolation error)
  // the estimate of the roundoff error is taken as follows:
  //    err = k*Sqrt(f(x)*f(x) + x*x*deriv*deriv)*Sqrt(Sum(ai)*(ai)),
  // where k is the double precision, ai are coefficients used in
  // central difference formulas
  // interpolation error is decreased by making the step size h smaller.
  //
  // Author: <NAME>
<commit_msg>add refrence to Richardson's paper for the Derivative functions<commit_after>// @(#)root/hist:$Name:  $:$Id: TF1.cxx,v 1.85 2004/08/02 11:22:52 brun Exp $
// Author: <NAME>   18/08/95

/*************************************************************************
 * Copyright (C) 1995-2000, <NAME> and <NAME>.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include "Riostream.h"
#include "TROOT.h"
#include "TMath.h"
#include "TF1.h"
#include "TH1.h"
#include "TVirtualPad.h"
#include "TStyle.h"
#include "TRandom.h"
#include "Api.h"
#include "TPluginManager.h"
#include "TVirtualUtilPad.h"
#include "TBrowser.h"

Bool_t TF1::fgAbsValue    = kFALSE;
Bool_t TF1::fgRejectPoint = kFALSE;
static TF1 *gHelper = 0;
static Double_t ErrorTF1 = 0;

ClassImp(TF1)

//______________________________________________________________________________
//
// a TF1 object is a 1-Dim function defined between a lower and upper limit.
// The function may be a simple function (see TFormula) or a precompiled
// user function.
// The function may have associated parameters.
// TF1 graphics function is via the TH1/TGraph drawing functions.
//
//  The following types of functions can be created:
//    A- Expression using variable x and no parameters
//    B- Expression using variable x with parameters
//    C- A general C function with parameters
//
//      Example of a function of type A
//
//   TF1 *f1 = new TF1("f1","sin(x)/x",0,10);
//   f1->Draw();
//Begin_Html
/*
<img src="gif/function1.gif">
*/
//End_Html
//
//      Example of a function of type B
//   TF1 *f1 = new TF1("f1","[0]*x*sin([1]*x)",-3,3);
//    This creates a function of variable x with 2 parameters.
//    The parameters must be initialized via:
//      f1->SetParameter(0,value_first_parameter);
//      f1->SetParameter(1,value_second_parameter);
//    Parameters may be given a name:
//      f1->SetParName(0,"Constant");
//
//     Example of function of type C
//   Consider the macro myfunc.C below
//-------------macro myfunc.C-----------------------------
//Double_t myfunction(Double_t *x, Double_t *par)
//{
//   Float_t xx =x[0];
//   Double_t f = TMath::Abs(par[0]*sin(par[1]*xx)/xx);
//   return f;
//}
//void myfunc()
//{
//   TF1 *f1 = new TF1("myfunc",myfunction,0,10,2);
//   f1->SetParameters(2,1);
//   f1->SetParNames("constant","coefficient");
//   f1->Draw();
//}
//void myfit()
//{
//   TH1F *h1=new TH1F("h1","test",100,0,10);
//   h1->FillRandom("myfunc",20000);
//   TF1 *f1=gROOT->GetFunction("myfunc");
//   f1->SetParameters(800,1);
//   h1.Fit("myfunc");
//}
//--------end of macro myfunc.C---------------------------------
//
// In an interactive session you can do:
//   Root > .L myfunc.C
//   Root > myfunc();
//   Root > myfit();
//
//
//  TF1 objects can reference other TF1 objects (thanks <NAME>)
//  of type A or B defined above.This excludes CINT interpreted functions
//  and compiled functions.
//  However, there is a restriction. A function cannot reference a basic
//  function if the basic function is a polynomial polN.
//  Example:
//{
//  TF1 *fcos = new TF1 ("fcos", "[0]*cos(x)", 0., 10.);
//  fcos->SetParNames( "cos");
//  fcos->SetParameter( 0, 1.1);
//
//  TF1 *fsin = new TF1 ("fsin", "[0]*sin(x)", 0., 10.);
//  fsin->SetParNames( "sin");
//  fsin->SetParameter( 0, 2.1);
//
//  TF1 *fsincos = new TF1 ("fsc", "fcos+fsin");
//
//  TF1 *fs2 = new TF1 ("fs2", "fsc+fsc");
//}
//
//     WHY TF1 CANNOT ACCEPT A CLASS MEMBER FUNCTION ?
//     ===============================================
// This is a frequently asked question.
// C++ is a strongly typed language. There is no way for TF1 (without
// recompiling this class) to know about all possible user defined data types.
// This also apply to the case of a static class function.
//
//------------------------------------------------------------------------

TF1 *TF1::fgCurrent = 0;

//______________________________________________________________________________
TF1::TF1(): TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*-*-*-*-*F1 default constructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ======================

   fXmin      = 0;
   fXmax      = 0;
   fNpx       = 100;
   fType      = 0;
   fNpfits    = 0;
   fNDF       = 0;
   fNsave     = 0;
   fChisquare = 0;
   fIntegral  = 0;
   fFunction  = 0;
   fParErrors = 0;
   fParMin    = 0;
   fParMax    = 0;
   fAlpha     = 0;
   fBeta      = 0;
   fGamma     = 0;
   fParent    = 0;
   fSave      = 0;
   fHistogram = 0;
   fMinimum   = -1111;
   fMaximum   = -1111;
   fMethodCall = 0;
   SetFillStyle(0);
}


//______________________________________________________________________________
TF1::TF1(const char *name,const char *formula, Double_t xmin, Double_t xmax)
      :TFormula(name,formula), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using a formula definition*-*-*-*-*-*-*-*-*-*-*
//*-*          =========================================
//*-*
//*-*  See TFormula constructor for explanation of the formula syntax.
//*-*
//*-*  See tutorials: fillrandom, first, fit1, formula1, multifit
//*-*  for real examples.
//*-*
//*-*  Creates a function of type A or B between xmin and xmax
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   if (xmin < xmax ) {
      fXmin      = xmin;
      fXmax      = xmax;
   } else {
      fXmin = xmax; //when called from TF2,TF3
      fXmax = xmin;
   }
   fNpx       = 100;
   fType      = 0;
   fFunction  = 0;
   if (fNpar) {
      fParErrors = new Double_t[fNpar];
      fParMin    = new Double_t[fNpar];
      fParMax    = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fNsave      = 0;
   fSave       = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fMethodCall = 0;
   
   if (fNdim != 1 && xmin < xmax) {
      Error("TF1","function: %s/%s has %d parameters instead of 1",name,formula,fNdim);
      MakeZombie();
   }

   if (!gStyle) return;
   SetLineColor(gStyle->GetFuncColor());
   SetLineWidth(gStyle->GetFuncWidth());
   SetLineStyle(gStyle->GetFuncStyle());
   SetFillStyle(0);
}


//______________________________________________________________________________
TF1::TF1(const char *name, Double_t xmin, Double_t xmax, Int_t npar)
      :TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using name of an interpreted function*-*-*-*
//*-*          =======================================================
//*-*
//*-*  Creates a function of type C between xmin and xmax.
//*-*  name is the name of an interpreted CINT cunction.
//*-*  The function is defined with npar parameters
//*-*  fcn must be a function of type:
//*-*     Double_t fcn(Double_t *x, Double_t *params)
//*-*
//*-*  This constructor is called for functions of type C by CINT.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   fXmin       = xmin;
   fXmax       = xmax;
   fNpx        = 100;
   fType       = 2;
   fFunction   = 0;
   if (npar > 0 ) fNpar = npar;
   if (fNpar) {
      fNames      = new TString[fNpar];
      fParams     = new Double_t[fNpar];
      fParErrors  = new Double_t[fNpar];
      fParMin     = new Double_t[fNpar];
      fParMax     = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParams[i]     = 0;
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fNsave      = 0;
   fSave       = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fMethodCall = 0;
   fNdim       = 1;

   TF1 *f1old = (TF1*)gROOT->GetListOfFunctions()->FindObject(name);
   if (f1old) delete f1old;
   SetName(name);

   if (gStyle) {
      SetLineColor(gStyle->GetFuncColor());
      SetLineWidth(gStyle->GetFuncWidth());
      SetLineStyle(gStyle->GetFuncStyle());
   }
   SetFillStyle(0);

   SetTitle(name);
   if (name) {
      if (*name == '*') return; //case happens via SavePrimitive
      fMethodCall = new TMethodCall();
      fMethodCall->InitWithPrototype(name,"Double_t*,Double_t*");
      fNumber = -1;
      gROOT->GetListOfFunctions()->Add(this);
   } else {
      Printf("Function:%s cannot be compiled",name);
   }
}

//______________________________________________________________________________
TF1::TF1(const char *name,void *fcn, Double_t xmin, Double_t xmax, Int_t npar)
      :TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using pointer to an interpreted function*-*-*-*
//*-*          =======================================================
//*-*
//*-*  See TFormula constructor for explanation of the formula syntax.
//*-*
//*-*  Creates a function of type C between xmin and xmax.
//*-*  The function is defined with npar parameters
//*-*  fcn must be a function of type:
//*-*     Double_t fcn(Double_t *x, Double_t *params)
//*-*
//*-*  see tutorial; myfit for an example of use
//*-*  also test/stress.cxx (see function stress1)
//*-*
//*-*
//*-*  This constructor is called for functions of type C by CINT.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   fXmin       = xmin;
   fXmax       = xmax;
   fNpx        = 100;
   fType       = 2;
   fFunction   = 0;
   if (npar > 0 ) fNpar = npar;
   if (fNpar) {
      fNames      = new TString[fNpar];
      fParams     = new Double_t[fNpar];
      fParErrors  = new Double_t[fNpar];
      fParMin     = new Double_t[fNpar];
      fParMax     = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParams[i]     = 0;
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fNsave      = 0;
   fSave       = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fMethodCall = 0;
   fNdim       = 1;

   TF1 *f1old = (TF1*)gROOT->GetListOfFunctions()->FindObject(name);
   if (f1old) delete f1old;
   SetName(name);

   if (gStyle) {
      SetLineColor(gStyle->GetFuncColor());
      SetLineWidth(gStyle->GetFuncWidth());
      SetLineStyle(gStyle->GetFuncStyle());
   }
   SetFillStyle(0);

   if (!fcn) return;
   char *funcname = G__p2f2funcname(fcn);
   SetTitle(funcname);
   if (funcname) {
      fMethodCall = new TMethodCall();
      fMethodCall->InitWithPrototype(funcname,"Double_t*,Double_t*");
      fNumber = -1;
      gROOT->GetListOfFunctions()->Add(this);
   } else {
      Printf("Function:%s cannot be compiled",name);
   }
}

//______________________________________________________________________________
TF1::TF1(const char *name,Double_t (*fcn)(Double_t *, Double_t *), Double_t xmin, Double_t xmax, Int_t npar)
      :TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*F1 constructor using a pointer to real function*-*-*-*-*-*-*-*
//*-*          ===============================================
//*-*
//*-*   npar is the number of free parameters used by the function
//*-*
//*-*   This constructor creates a function of type C when invoked
//*-*   with the normal C++ compiler.
//*-*
//*-*   see test program test/stress.cxx (function stress1) for an example.
//*-*   note the interface with an intermediate pointer.
//*-*
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   fXmin       = xmin;
   fXmax       = xmax;
   fNpx        = 100;
   char *funcname = G__p2f2funcname((void*)fcn);
   if (funcname) {
      fType       = 2;
      SetTitle(funcname);
      fMethodCall = new TMethodCall();
      fMethodCall->InitWithPrototype(funcname,"Double_t*,Double_t*");
      fNumber = -1;
      fFunction   = 0;
   } else {
      fType       = 1;
      fMethodCall = 0;
      fFunction   = fcn;
   }
   if (npar > 0 ) fNpar = npar;
   if (fNpar) {
      fNames      = new TString[fNpar];
      fParams     = new Double_t[fNpar];
      fParErrors  = new Double_t[fNpar];
      fParMin     = new Double_t[fNpar];
      fParMax     = new Double_t[fNpar];
      for (int i = 0; i < fNpar; i++) {
         fParams[i]     = 0;
         fParErrors[i]  = 0;
         fParMin[i]     = 0;
         fParMax[i]     = 0;
      }
   } else {
      fParErrors = 0;
      fParMin    = 0;
      fParMax    = 0;
   }
   fChisquare  = 0;
   fIntegral   = 0;
   fAlpha      = 0;
   fBeta       = 0;
   fGamma      = 0;
   fNsave      = 0;
   fSave       = 0;
   fParent     = 0;
   fNpfits     = 0;
   fNDF        = 0;
   fHistogram  = 0;
   fMinimum    = -1111;
   fMaximum    = -1111;
   fNdim       = 1;
//*-*- Store formula in linked list of formula in ROOT

   TF1 *f1old = (TF1*)gROOT->GetListOfFunctions()->FindObject(name);
   if (f1old) delete f1old;
   SetName(name);
   gROOT->GetListOfFunctions()->Add(this);

   if (!gStyle) return;
   SetLineColor(gStyle->GetFuncColor());
   SetLineWidth(gStyle->GetFuncWidth());
   SetLineStyle(gStyle->GetFuncStyle());
   SetFillStyle(0);

}

//______________________________________________________________________________
TF1& TF1::operator=(const TF1 &rhs) 
{
   if (this != &rhs) {
      rhs.Copy(*this);
   }
   return *this;
}

//______________________________________________________________________________
TF1::~TF1()
{
//*-*-*-*-*-*-*-*-*-*-*F1 default destructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =====================

   if (fParMin)    delete [] fParMin;
   if (fParMax)    delete [] fParMax;
   if (fParErrors) delete [] fParErrors;
   if (fIntegral)  delete [] fIntegral;
   if (fAlpha)     delete [] fAlpha;
   if (fBeta)      delete [] fBeta;
   if (fGamma)     delete [] fGamma;
   if (fSave)      delete [] fSave;
   delete fHistogram;
   delete fMethodCall;

   if (fParent) {
      if (fParent->InheritsFrom(TH1::Class())) {
         ((TH1*)fParent)->GetListOfFunctions()->Remove(this);
         return;
      }
      //parent may be a graph, or ??
      //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->RemoveObject(fParent,this);
      fParent = 0;
   }
}

//______________________________________________________________________________
TF1::TF1(const TF1 &f1) : TFormula(), TAttLine(f1), TAttFill(f1), TAttMarker(f1)
{

   fXmin      = 0;
   fXmax      = 0;
   fNpx       = 100;
   fType      = 0;
   fNpfits    = 0;
   fNDF       = 0;
   fNsave     = 0;
   fChisquare = 0;
   fIntegral  = 0;
   fFunction  = 0;
   fParErrors = 0;
   fParMin    = 0;
   fParMax    = 0;
   fAlpha     = 0;
   fBeta      = 0;
   fGamma     = 0;
   fParent    = 0;
   fSave      = 0;
   fHistogram = 0;
   fMinimum   = -1111;
   fMaximum   = -1111;
   fMethodCall = 0;
   SetFillStyle(0);

   ((TF1&)f1).Copy(*this);
}

//______________________________________________________________________________
void TF1::AbsValue(Bool_t flag)
{
   // static function: set the fgAbsValue flag.
   // By default TF1::Integral uses the original function value to compute the integral
   // However, TF1::Moment, CentralMoment require to compute the integral
   // using the absolute value of the function.
   
   fgAbsValue = flag;
}

//______________________________________________________________________________
void TF1::Browse(TBrowser *b)
{
    Draw(b ? b->GetDrawOption() : "");
    gPad->Update();
}


//______________________________________________________________________________
void TF1::Copy(TObject &obj) const
{
//*-*-*-*-*-*-*-*-*-*-*Copy this F1 to a new F1*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ========================

   if (((TF1&)obj).fParMin)    delete [] ((TF1&)obj).fParMin;
   if (((TF1&)obj).fParMax)    delete [] ((TF1&)obj).fParMax;
   if (((TF1&)obj).fParErrors) delete [] ((TF1&)obj).fParErrors;
   if (((TF1&)obj).fIntegral)  delete [] ((TF1&)obj).fIntegral;
   if (((TF1&)obj).fAlpha)     delete [] ((TF1&)obj).fAlpha;
   if (((TF1&)obj).fBeta)      delete [] ((TF1&)obj).fBeta;
   if (((TF1&)obj).fGamma)     delete [] ((TF1&)obj).fGamma;
   if (((TF1&)obj).fSave)      delete [] ((TF1&)obj).fSave;
   delete ((TF1&)obj).fHistogram;
   delete ((TF1&)obj).fMethodCall;

   TFormula::Copy(obj);
   TAttLine::Copy((TF1&)obj);
   TAttFill::Copy((TF1&)obj);
   TAttMarker::Copy((TF1&)obj);
   ((TF1&)obj).fXmin = fXmin;
   ((TF1&)obj).fXmax = fXmax;
   ((TF1&)obj).fNpx  = fNpx;
   ((TF1&)obj).fType = fType;
   ((TF1&)obj).fFunction  = fFunction;
   ((TF1&)obj).fChisquare = fChisquare;
   ((TF1&)obj).fNpfits  = fNpfits;
   ((TF1&)obj).fNDF     = fNDF;
   ((TF1&)obj).fMinimum = fMinimum;
   ((TF1&)obj).fMaximum = fMaximum;

   ((TF1&)obj).fParErrors = 0;
   ((TF1&)obj).fParMin    = 0;
   ((TF1&)obj).fParMax    = 0;
   ((TF1&)obj).fIntegral  = 0;
   ((TF1&)obj).fAlpha     = 0;
   ((TF1&)obj).fBeta      = 0;
   ((TF1&)obj).fGamma     = 0;
   ((TF1&)obj).fParent    = fParent;
   ((TF1&)obj).fNsave     = fNsave;
   ((TF1&)obj).fSave      = 0;
   ((TF1&)obj).fHistogram = 0;
   ((TF1&)obj).fMethodCall = 0;
   if (fNsave) {
      ((TF1&)obj).fSave = new Double_t[fNsave];
      for (Int_t j=0;j<fNsave;j++) ((TF1&)obj).fSave[j] = fSave[j];
   }
   if (fNpar) {
      ((TF1&)obj).fParErrors = new Double_t[fNpar];
      ((TF1&)obj).fParMin    = new Double_t[fNpar];
      ((TF1&)obj).fParMax    = new Double_t[fNpar];
      Int_t i;
      for (i=0;i<fNpar;i++)   ((TF1&)obj).fParErrors[i] = fParErrors[i];
      for (i=0;i<fNpar;i++)   ((TF1&)obj).fParMin[i]    = fParMin[i];
      for (i=0;i<fNpar;i++)   ((TF1&)obj).fParMax[i]    = fParMax[i];
   }
   if (fMethodCall) {
      TMethodCall *m = new TMethodCall();
      m->InitWithPrototype(fMethodCall->GetMethodName(),fMethodCall->GetProto());
      ((TF1&)obj).fMethodCall  = m;
   }
}

//______________________________________________________________________________
Double_t TF1::Derivative(Double_t x, Double_t *params, Double_t eps) const
{
  // returns the first derivative of the function at point x, 
  // computed by Richardson's extrapolation method (use 2 derivative estimates 
  // to compute a third, more accurate estimation)
  // first, derivatives with steps h and h/2 are computed by central difference formulas
  // D(h) = (f(x+h) - f(x-h))/2h
  // the final estimate D = (4*D(h/2) - D(h))/3
  //  "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition"
  //
  // if the argument params is null, the current function parameters are used,
  // otherwise the parameters in params are used.
  //
  // the argument eps may be specified to control the step size (precision).
  // the step size is taken as eps*(xmax-xmin).
  // the default value (0.001) should be good enough for the vast majority
  // of functions. Give a smaller value if your function has many changes
  // of the second derivative in the function range.
  //
  // Getting the error via TF1::DerivativeError
  // -----------------
  //   (total error = roundoff error + interpolation error)
  // the estimate of the roundoff error is taken as follows:
  //    err = k*Sqrt(f(x)*f(x) + x*x*deriv*deriv)*Sqrt(Sum(ai)*(ai)),
  // where k is the double precision, ai are coefficients used in
  // central difference formulas
  // interpolation error is decreased by making the step size h smaller.
  //
  // Author: <NAME>
  
  const Double_t C1 = 1e-15;
   
  if(eps< 1e-10 || eps > 1e-2) {
     Warning("Derivative","parameter esp=%g out of allowed range[1e-10,1e-2], reset to 0.001",eps);
     eps = 0.001;
  }
  Double_t xmin, xmax;
  GetRange(xmin, xmax);
  Double_t h = eps*(xmax-xmin);

  Double_t xx[1];
  TF1 *func = (TF1*)this;
  func->InitArgs(xx, params);
  xx[0] = x+h;     Double_t f1 = func->EvalPar(xx, params);
  xx[0] = x;       Double_t fx = func->EvalPar(xx, params);
  xx[0] = x-h;     Double_t f2 = func->EvalPar(xx, params);

  xx[0] = x+h/2;   Double_t g1 = func->EvalPar(xx, params);
  xx[0] = x-h/2;   Double_t g2 = func->EvalPar(xx, params);

  //compute the central differences
  Double_t h2    = 1/(2.*h);
  Double_t D0    = f1 - f2;
  Double_t D2    = 2*(g1 - g2);
  ErrorTF1       = C1*h2*fx;  //compute the error
  Double_t deriv = h2*(4*D2 - D0)/3.;  
  return deriv;
}


//______________________________________________________________________________
Double_t TF1::Derivative2(Double_t x, Double_t *params, Double_t eps) const
{
  // returns the first derivative of the function at point x, 
  // computed by Richardson's extrapolation method (use 2 derivative estimates 
  // to compute a third, more accurate estimation)
  // first, derivatives with steps h and h/2 are computed by central difference formulas
  //    D(h) = (f(x+h) - 2*f(x) + f(x-h))/(h*h)
  // the final estimate D = (4*D(h/2) - D(h))/3
  //  "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition"
  //
  // if the argument params is null, the current function parameters are used,
  // otherwise the parameters in params are used.
  //
  // the argument eps may be specified to control the step size (precision).
  // the step size is taken as eps*(xmax-xmin).
  // the default value (0.001) should be good enough for the vast majority
  // of functions. Give a smaller value if your function has many changes
  // of the second derivative in the function range.
  //
  // Getting the error via TF1::DerivativeError
  // -----------------
  //   (total error = roundoff error + interpolation error)
  // the estimate of the roundoff error is taken as follows:
  //    err = k*Sqrt(f(x)*f(x) + x*x*deriv*deriv)*Sqrt(Sum(ai)*(ai)),
  // where k is the double precision, ai are coefficients used in
  // central difference formulas
  // interpolation error is decreased by making the step size h smaller.
  //
  // Author: <NAME>

  const Double_t C1 = 2*1e-15;
   
  if(eps< 1e-6 || eps > 1e-2) {
     Warning("Derivative2","parameter esp=%g out of allowed range[1e-6,1e-2], reset to 0.001",eps);
     eps = 0.001;
  }
  Double_t xmin, xmax;
  GetRange(xmin, xmax);
  Double_t h = eps*(xmax-xmin);

  Double_t xx[1];
  TF1 *func = (TF1*)this;
  func->InitArgs(xx, params);
  xx[0] = x+h;     Double_t f1 = func->EvalPar(xx, params);
  xx[0] = x;       Double_t f2 = func->EvalPar(xx, params);
  xx[0] = x-h;     Double_t f3 = func->EvalPar(xx, params);

  xx[0] = x+h/2;   Double_t g1 = func->EvalPar(xx, params);
  xx[0] = x-h/2;   Double_t g3 = func->EvalPar(xx, params);

  //compute the central differences
  Double_t hh    = 1/(h*h);
  Double_t D0    = f3 - 2*f2 + f1;
  Double_t D2    = 4*g3 - 8*f2 +4*g1;
  ErrorTF1       = C1*hh*f2;  //compute the error
  Double_t deriv = hh*(4*D2 - D0)/3.;
  return deriv;
}


//______________________________________________________________________________
Double_t TF1::Derivative3(Double_t x, Double_t *params, Double_t eps) const
{
  // returns the first derivative of the function at point x, 
  // computed by Richardson's extrapolation method (use 2 derivative estimates 
  // to compute a third, more accurate estimation)
  // first, derivatives with steps h and h/2 are computed by central difference formulas
  //    D(h) = (f(x+2h) - 2*f(x+h) + 2*f(x-h) - f(x-2h))/(2*h*h*h)
  // the final estimate D = (4*D(h/2) - D(h))/3
  //  "Numerical Methods for Scientists and Engineers", H.M.Antia, 2nd edition"
  //
  // if the argument params is null, the current function parameters are used,
  // otherwise the parameters in params are used.
  //
  // the argument eps may be specified to control the step size (precision).
  // the step size is taken as eps*(xmax-xmin).
  // the default value (0.001) should be good enough for the vast majority
  // of functions. Give a smaller value if your function has many changes
  // of the second derivative in the function range.
  //
  // Getting the error via TF1::DerivativeError
  // -----------------
  //   (total error = roundoff error + interpolation error)
  // the estimate of the roundoff error is taken as follows:
  //    err = k*Sqrt(f(x)*f(x) + x*x*deriv*deriv)*Sqrt(Sum(ai)*(ai)),
  // where k is the double precision, ai are coefficients used in
  // central difference formulas
  // interpolation error is decreased by making the step size h smaller.
  //
  // Author: <NAME>

{'fraction_non_alphanumeric': 0.18657538558388417, 'fraction_numerical': 0.02504327982373308, 'mean_word_length': 2.4291014570966003, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 175, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4392249', 'n_tokens_mistral': 18616, 'n_tokens_neox': 17564, 'n_words': 4497}