#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 "TBrowser.h"
#include "TColor.h"
#include "TClass.h"
Bool_t TF1::fgAbsValue = kFALSE;
Bool_t TF1::fgRejectPoint = kFALSE;
static TF1 *gHelper = 0;
static Double_t gErrorTF1 = 0;
ClassImp(TF1)
//Begin_Html
/*
<img src="gif/function1.gif">
*/
//End_Html
TF1 *TF1::fgCurrent = 0;
TF1::TF1(): TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
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()
{
if (xmin < xmax ) {
fXmin = xmin;
fXmax = xmax;
} else {
fXmin = xmax;
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()
{
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;
fMethodCall = new TMethodCall();
fMethodCall->InitWithPrototype(name,"Double_t*,Double_t*");
fNumber = -1;
gROOT->GetListOfFunctions()->Add(this);
if (! fMethodCall->IsValid() ) {
Error("TF1","No function found with the signature %s(Double_t*,Double_t*)",name);
}
} else {
Error("TF1","requires a proper function name!");
}
}
TF1::TF1(const char *name,void *fcn, Double_t xmin, Double_t xmax, Int_t npar)
:TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
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);
if (! fMethodCall->IsValid() ) {
Error("TF1","No function found with the signature %s(Double_t*,Double_t*)",funcname);
}
} else {
Error("TF1","can not find any function at the address 0x%x. This function requested for %s",fcn,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()
{
fXmin = xmin;
fXmax = xmax;
fNpx = 100;
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;
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(const char *name,Double_t (*fcn)(const Double_t *, const Double_t *), Double_t xmin, Double_t xmax, Int_t npar)
:TFormula(), TAttLine(), TAttFill(), TAttMarker()
{
fXmin = xmin;
fXmax = xmax;
fNpx = 100;
fType = 1;
fMethodCall = 0;
typedef Double_t (*Function_t) (Double_t *, Double_t *);
fFunction = (Function_t)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;
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()
{
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) fParent->RecursiveRemove(this);
}
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)
{
fgAbsValue = flag;
}
void TF1::Browse(TBrowser *b)
{
Draw(b ? b->GetDrawOption() : "");
gPad->Update();
}
void TF1::Copy(TObject &obj) const
{
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
{
const Double_t kC1 = 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);
Double_t h2 = 1/(2.*h);
Double_t d0 = f1 - f2;
Double_t d2 = 2*(g1 - g2);
gErrorTF1 = kC1*h2*fx;
Double_t deriv = h2*(4*d2 - d0)/3.;
return deriv;
}
Double_t TF1::Derivative2(Double_t x, Double_t *params, Double_t eps) const
{
const Double_t kC1 = 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);
Double_t hh = 1/(h*h);
Double_t d0 = f3 - 2*f2 + f1;
Double_t d2 = 4*g3 - 8*f2 +4*g1;
gErrorTF1 = kC1*hh*f2;
Double_t deriv = hh*(4*d2 - d0)/3.;
return deriv;
}
Double_t TF1::Derivative3(Double_t x, Double_t *params, Double_t eps) const
{
const Double_t kC1 = 1e-15;
if(eps< 1e-4 || eps > 1e-2) {
Warning("Derivative3","parameter esp=%g out of allowed range[1e-4,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+2*h; Double_t f1 = func->EvalPar(xx, params);
xx[0] = x+h; Double_t f2 = func->EvalPar(xx, params);
xx[0] = x-h; Double_t f3 = func->EvalPar(xx, params);
xx[0] = x-2*h; Double_t f4 = func->EvalPar(xx, params);
xx[0] = x; Double_t fx = func->EvalPar(xx, params);
xx[0] = x+h/2; Double_t g2 = func->EvalPar(xx, params);
xx[0] = x-h/2; Double_t g3 = func->EvalPar(xx, params);
Double_t hhh = 1/(h*h*h);
Double_t d0 = 0.5*f1 - f2 +f3 - 0.5*f4;
Double_t d2 = 4*f2 - 8*g2 +8*g3 - 4*f3;
gErrorTF1 = kC1*hhh*fx;
Double_t deriv = hhh*(4*d2 - d0)/3.;
return deriv;
}
Double_t TF1::DerivativeError()
{
return gErrorTF1;
}
Int_t TF1::DistancetoPrimitive(Int_t px, Int_t py)
{
if (!fHistogram) return 9999;
Int_t distance = 9999;
if (px >= 0) {
distance = fHistogram->DistancetoPrimitive(px,py);
if (distance <= 1) return distance;
} else {
px = -px;
}
Double_t xx[1];
Double_t x = gPad->AbsPixeltoX(px);
xx[0] = gPad->PadtoX(x);
if (xx[0] < fXmin || xx[0] > fXmax) return distance;
Double_t fval = Eval(xx[0]);
Double_t y = gPad->YtoPad(fval);
Int_t pybin = gPad->YtoAbsPixel(y);
return TMath::Abs(py - pybin);
}
void TF1::Draw(Option_t *option)
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
AppendPad(option);
}
TF1 *TF1::DrawCopy(Option_t *option) const
{
TF1 *newf1 = new TF1();
Copy(*newf1);
newf1->AppendPad(option);
newf1->SetBit(kCanDelete);
return newf1;
}
void TF1::DrawDerivative(Option_t *option)
{
TVirtualPad *pad = gROOT->GetSelectedPad();
TVirtualPad *padsav = gPad;
if (pad) pad->cd();
char cmd[512];
sprintf(cmd,"{TGraph *R__%s_Derivative = new TGraph((TF1*)0x%lx,\"d\");R__%s_Derivative->Draw(\"%s\");}",GetName(),(Long_t)this,GetName(),option);
gROOT->ProcessLine(cmd);
if (padsav) padsav->cd();
}
void TF1::DrawIntegral(Option_t *option)
{
TVirtualPad *pad = gROOT->GetSelectedPad();
TVirtualPad *padsav = gPad;
if (pad) pad->cd();
char cmd[512];
sprintf(cmd,"{TGraph *R__%s_Integral = new TGraph((TF1*)0x%lx,\"i\");R__%s_Integral->Draw(\"%s\");}",GetName(),(Long_t)this,GetName(),option);
gROOT->ProcessLine(cmd);
if (padsav) padsav->cd();
}
void TF1::DrawF1(const char *formula, Double_t xmin, Double_t xmax, Option_t *option)
{
if (Compile(formula)) return;
SetRange(xmin, xmax);
Draw(option);
}
Double_t TF1::Eval(Double_t x, Double_t y, Double_t z, Double_t t) const
{
Double_t xx[4];
xx[0] = x;
xx[1] = y;
xx[2] = z;
xx[3] = t;
((TF1*)this)->InitArgs(xx,fParams);
return ((TF1*)this)->EvalPar(xx,fParams);
}
Double_t TF1::EvalPar(const Double_t *x, const Double_t *params)
{
fgCurrent = this;
if (fType == 0) return TFormula::EvalPar(x,params);
Double_t result = 0;
if (fType == 1) {
if (fFunction) {
if (params) result = (*fFunction)((Double_t*)x,(Double_t*)params);
else result = (*fFunction)((Double_t*)x,fParams);
}else result = GetSave(x);
}
if (fType == 2) {
if (fMethodCall) fMethodCall->Execute(result);
else result = GetSave(x);
}
return result;
}
void TF1::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
fHistogram->ExecuteEvent(event,px,py);
if (!gPad->GetView()) {
if (event == kMouseMotion) gPad->SetCursor(kHand);
}
}
void TF1::FixParameter(Int_t ipar, Double_t value)
{
if (ipar < 0 || ipar > fNpar-1) return;
SetParameter(ipar,value);
if (value != 0) SetParLimits(ipar,value,value);
else SetParLimits(ipar,1,1);
}
TF1 *TF1::GetCurrent()
{
return fgCurrent;
}
TH1 *TF1::GetHistogram() const
{
if (fHistogram) return fHistogram;
((TF1*)this)->Paint();
return fHistogram;
}
Double_t TF1::GetMaximum(Double_t xmin, Double_t xmax) const
{
if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
Double_t x;
Int_t niter=0;
x = MinimStep(3, xmin, xmax, 0);
Bool_t ok = kTRUE;
x = MinimBrent(3, xmin, xmax, x, 0, ok);
while (!ok){
if (niter>10){
Error("GetMaximum", "maximum search didn't converge");
break;
}
x=MinimStep(3, xmin, xmax,0);
x = MinimBrent(3, xmin, xmax, x, 0, ok);
niter++;
}
return x;
}
Double_t TF1::GetMaximumX(Double_t xmin, Double_t xmax) const
{
if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
Double_t x;
Int_t niter=0;
x = MinimStep(2, xmin, xmax, 0);
Bool_t ok = kTRUE;
x = MinimBrent(2, xmin, xmax, x, 0, ok);
while (!ok){
if (niter>10){
Error("GetMaximumX", "maximum search didn't converge");
break;
}
x=MinimStep(2, xmin, xmax, 0);
x = MinimBrent(2, xmin, xmax, x, 0, ok);
niter++;
}
return x;
}
Double_t TF1::GetMinimum(Double_t xmin, Double_t xmax) const
{
if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
Double_t x;
Int_t niter=0;
x = MinimStep(1, xmin, xmax, 0);
Bool_t ok = kTRUE;
x = MinimBrent(1, xmin, xmax, x, 0, ok);
while (!ok){
if (niter>10){
Error("GetMinimum", "minimum search didn't converge");
break;
}
x=MinimStep(1, xmin, xmax,0);
x = MinimBrent(1, xmin, xmax, x, 0, ok);
niter++;
}
return x;
}
Double_t TF1::GetMinimumX(Double_t xmin, Double_t xmax) const
{
if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
Int_t niter=0;
Double_t x;
x = MinimStep(0, xmin, xmax, 0);
Bool_t ok = kTRUE;
x = MinimBrent(0, xmin, xmax, x, 0, ok);
while (!ok){
if (niter>10){
Error("GetMinimumX", "minimum search didn't converge");
break;
}
x=MinimStep(0, xmin, xmax,0);
x = MinimBrent(0, xmin, xmax, x, 0, ok);
niter++;
}
return x;
}
Double_t TF1::GetX(Double_t fy, Double_t xmin, Double_t xmax) const
{
if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
Int_t niter=0;
Double_t x;
x = MinimStep(4, xmin, xmax, fy);
Bool_t ok = kTRUE;
x = MinimBrent(4, xmin, xmax, x, fy, ok);
while (!ok){
if (niter>10){
Error("GetX", "Search didn't converge");
break;
}
x=MinimStep(4, xmin, xmax, fy);
x = MinimBrent(4, xmin, xmax, x, fy, ok);
niter++;
}
return x;
}
Double_t TF1::MinimStep(Int_t type, Double_t &xmin, Double_t &xmax, Double_t fy) const
{
Double_t x,y, dx;
dx = (xmax-xmin)/(fNpx-1);
Double_t xxmin = xmin;
Double_t yymin;
if (type < 2)
yymin = Eval(xmin);
else if (type < 4)
yymin = -Eval(xmin);
else
yymin = TMath::Abs(Eval(xmin)-fy);
for (Int_t i=1; i<=fNpx-1; i++) {
x = xmin + i*dx;
if (type < 2)
y = Eval(x);
else if (type < 4)
y = -Eval(x);
else
y = TMath::Abs(Eval(x)-fy);
if (y < yymin) {xxmin = x; yymin = y;}
}
xmin = TMath::Max(xmin,xxmin-dx);
xmax = TMath::Min(xmax,xxmin+dx);
return TMath::Min(xxmin, xmax);
}
Double_t TF1::MinimBrent(Int_t type, Double_t &xmin, Double_t &xmax, Double_t xmiddle, Double_t fy, Bool_t &ok) const
{
Double_t eps = 1e-10;
Double_t t = 1e-8;
Int_t itermax = 100;
Double_t c = (3.-TMath::Sqrt(5.))/2.;
Double_t u, v, w, x, fv, fu, fw, fx, e, p, q, r, t2, d=0, m, tol;
v = w = x = xmiddle;
e=0;
Double_t a=xmin;
Double_t b=xmax;
if (type < 2)
fv = fw = fx = Eval(x);
else if (type < 4)
fv = fw = fx = -Eval(x);
else
fv = fw = fx = TMath::Abs(Eval(x)-fy);
for (Int_t i=0; i<itermax; i++){
m=0.5*(a + b);
tol = eps*(TMath::Abs(x))+t;
t2 = 2*tol;
if (TMath::Abs(x-m) <= (t2-0.5*(b-a))) {
ok=kTRUE;
if (type==1)
return fx;
else if (type==3)
return -fx;
else
return x;
}
if (TMath::Abs(e)>tol){
r = (x-w)*(fx-fv);
q = (x-v)*(fx-fw);
p = (x-v)*q - (x-w)*r;
q = 2*(q-r);
if (q>0) p=-p;
else q=-q;
r=e;
e=d;
if (TMath::Abs(p) < TMath::Abs(0.5*q*r) || p < q*(a-x) || p < q*(b-x)) {
d = p/q;
u = x+d;
if (u-a < t2 || b-u < t2)
d=TMath::Sign(tol, m-x);
} else {
e=(x>=m ? a-x : b-x);
d = c*e;
}
} else {
e=(x>=m ? a-x : b-x);
d = c*e;
}
u = (TMath::Abs(d)>=tol ? x+d : x+TMath::Sign(tol, d));
if (type < 2)
fu = Eval(u);
else if (type < 4)
fu = -Eval(u);
else
fu = TMath::Abs(Eval(u)-fy);
if (fu<=fx){
if (u<x) b=x;
else a=x;
v=w; fv=fw; w=x; fw=fx; x=u; fx=fu;
} else {
if (u<x) a=u;
else b=u;
if (fu<=fw || w==x){
v=w; fv=fw; w=u; fw=fu;
}
else if (fu<=fv || v==x || v==w){
v=u; fv=fu;
}
}
}
ok = kFALSE;
xmin = a;
xmax = b;
return x;
}
Int_t TF1::GetNDF() const
{
if (fNDF == 0) return fNpfits-fNpar;
return fNDF;
}
Int_t TF1::GetNumberFreeParameters() const
{
Int_t nfree = fNpar;
Double_t al,bl;
for (Int_t i=0;i<fNpar;i++) {
((TF1*)this)->GetParLimits(i,al,bl);
if (al*bl != 0 && al >= bl) nfree--;
}
return nfree;
}
char *TF1::GetObjectInfo(Int_t px, Int_t ) const
{
static char info[64];
Double_t x = gPad->PadtoX(gPad->AbsPixeltoX(px));
sprintf(info,"(x=%g, f=%g)",x,((TF1*)this)->Eval(x));
return info;
}
Double_t TF1::GetParError(Int_t ipar) const
{
if (ipar < 0 || ipar > fNpar-1) return 0;
return fParErrors[ipar];
}
void TF1::GetParLimits(Int_t ipar, Double_t &parmin, Double_t &parmax) const
{
parmin = 0;
parmax = 0;
if (ipar < 0 || ipar > fNpar-1) return;
if (fParMin) parmin = fParMin[ipar];
if (fParMax) parmax = fParMax[ipar];
}
Double_t TF1::GetProb() const
{
if (fNDF <= 0) return 0;
return TMath::Prob(fChisquare,fNDF);
}
Int_t TF1::GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum)
{
const Int_t npx = TMath::Min(250,TMath::Max(50,2*nprobSum));
const Double_t xMin = GetXmin();
const Double_t xMax = GetXmax();
const Double_t dx = (xMax-xMin)/npx;
TArrayD integral(npx+1);
TArrayD alpha(npx);
TArrayD beta(npx);
TArrayD gamma(npx);
integral[0] = 0;
Int_t intNegative = 0;
Int_t i;
for (i = 0; i < npx; i++) {
const Double_t *params = 0;
Double_t integ = Integral(Double_t(xMin+i*dx),Double_t(xMin+i*dx+dx),params);
if (integ < 0) {intNegative++; integ = -integ;}
integral[i+1] = integral[i] + integ;
}
if (intNegative > 0)
Warning("GetQuantiles","function:%s has %d negative values: abs assumed",
GetName(),intNegative);
if (integral[npx] == 0) {
Error("GetQuantiles","Integral of function is zero");
return 0;
}
const Double_t total = integral[npx];
for (i = 1; i <= npx; i++) integral[i] /= total;
for (i = 0; i < npx; i++) {
const Double_t x0 = xMin+dx*i;
const Double_t r2 = integral[i+1]-integral[i];
const Double_t r1 = Integral(x0,x0+0.5*dx)/total;
gamma[i] = (2*r2-4*r1)/(dx*dx);
beta[i] = r2/dx-gamma[i]*dx;
alpha[i] = x0;
gamma[i] *= 2;
}
for (i = 0; i < nprobSum; i++) {
const Double_t r = probSum[i];
Int_t bin = TMath::Max(TMath::BinarySearch(npx+1,integral.GetArray(),r)-1,(Long64_t)0);
while (bin < npx-1 && integral[bin+1] == r) {
if (integral[bin+2] == r) bin++;
else break;
}
const Double_t rr = r-integral[bin];
if (rr != 0.0) {
Double_t xx = 0.0;
const Double_t fac = -2.*gamma[bin]*rr/beta[bin]/beta[bin];
if (fac != 0 && fac <= 1)
xx = (-beta[bin]+TMath::Sqrt(beta[bin]*beta[bin]+2*gamma[bin]*rr))/gamma[bin];
else if (beta[bin] != 0.)
xx = rr/beta[bin];
q[i] = alpha[bin]+xx;
} else {
q[i] = alpha[bin];
if (integral[bin+1] == r) q[i] += dx;
}
}
return nprobSum;
}
Double_t TF1::GetRandom()
{
if (fIntegral == 0) {
Double_t dx = (fXmax-fXmin)/fNpx;
fIntegral = new Double_t[fNpx+1];
fAlpha = new Double_t[fNpx];
fBeta = new Double_t[fNpx];
fGamma = new Double_t[fNpx];
fIntegral[0] = 0;
Double_t integ;
Int_t intNegative = 0;
Int_t i;
for (i=0;i<fNpx;i++) {
integ = Integral(Double_t(fXmin+i*dx), Double_t(fXmin+i*dx+dx));
if (integ < 0) {intNegative++; integ = -integ;}
fIntegral[i+1] = fIntegral[i] + integ;
}
if (intNegative > 0) {
Warning("GetRandom","function:%s has %d negative values: abs assumed",GetName(),intNegative);
}
if (fIntegral[fNpx] == 0) {
Error("GetRandom","Integral of function is zero");
return 0;
}
Double_t total = fIntegral[fNpx];
for (i=1;i<=fNpx;i++) {
fIntegral[i] /= total;
}
Double_t x0,r1,r2,r3;
for (i=0;i<fNpx;i++) {
x0 = fXmin+i*dx;
r2 = fIntegral[i+1] - fIntegral[i];
r1 = Integral(x0,x0+0.5*dx)/total;
r3 = 2*r2 - 4*r1;
if (TMath::Abs(r3) > 1e-8) fGamma[i] = r3/(dx*dx);
else fGamma[i] = 0;
fBeta[i] = r2/dx - fGamma[i]*dx;
fAlpha[i] = x0;
fGamma[i] *= 2;
}
}
Double_t r = gRandom->Rndm();
Int_t bin = TMath::BinarySearch(fNpx,fIntegral,r);
Double_t rr = r - fIntegral[bin];
Double_t xx;
if(fGamma[bin] != 0)
xx = (-fBeta[bin] + TMath::Sqrt(fBeta[bin]*fBeta[bin]+2*fGamma[bin]*rr))/fGamma[bin];
else
xx = rr/fBeta[bin];
Double_t x = fAlpha[bin] + xx;
return x;
}
Double_t TF1::GetRandom(Double_t xmin, Double_t xmax)
{
if (fIntegral == 0) {
Double_t dx = (fXmax-fXmin)/fNpx;
fIntegral = new Double_t[fNpx+1];
fAlpha = new Double_t[fNpx];
fBeta = new Double_t[fNpx];
fGamma = new Double_t[fNpx];
fIntegral[0] = 0;
Double_t integ;
Int_t intNegative = 0;
Int_t i;
for (i=0;i<fNpx;i++) {
integ = Integral(Double_t(fXmin+i*dx), Double_t(fXmin+i*dx+dx));
if (integ < 0) {intNegative++; integ = -integ;}
fIntegral[i+1] = fIntegral[i] + integ;
}
if (intNegative > 0) {
Warning("GetRandom","function:%s has %d negative values: abs assumed",GetName(),intNegative);
}
if (fIntegral[fNpx] == 0) {
Error("GetRandom","Integral of function is zero");
return 0;
}
Double_t total = fIntegral[fNpx];
for (i=1;i<=fNpx;i++) {
fIntegral[i] /= total;
}
Double_t x0,r1,r2,r3;
for (i=0;i<fNpx;i++) {
x0 = fXmin+i*dx;
r2 = fIntegral[i+1] - fIntegral[i];
r1 = Integral(x0,x0+0.5*dx)/total;
r3 = 2*r2 - 4*r1;
if (TMath::Abs(r3) > 1e-8) fGamma[i] = r3/(dx*dx);
else fGamma[i] = 0;
fBeta[i] = r2/dx - fGamma[i]*dx;
fAlpha[i] = x0;
fGamma[i] *= 2;
}
}
Double_t dx = (fXmax-fXmin)/fNpx;
Int_t nbinmin = (Int_t)((xmin-fXmin)/dx);
Int_t nbinmax = (Int_t)((xmax-fXmin)/dx)+2;
if(nbinmax>fNpx) nbinmax=fNpx;
Double_t pmin=fIntegral[nbinmin];
Double_t pmax=fIntegral[nbinmax];
Double_t r,x,xx,rr;
do {
r = gRandom->Uniform(pmin,pmax);
Int_t bin = TMath::BinarySearch(fNpx,fIntegral,r);
rr = r - fIntegral[bin];
if(fGamma[bin] != 0)
xx = (-fBeta[bin] + TMath::Sqrt(fBeta[bin]*fBeta[bin]+2*fGamma[bin]*rr))/fGamma[bin];
else
xx = rr/fBeta[bin];
x = fAlpha[bin] + xx;
} while(x<xmin || x>xmax);
return x;
}
void TF1::GetRange(Double_t &xmin, Double_t &xmax) const
{
xmin = fXmin;
xmax = fXmax;
}
void TF1::GetRange(Double_t &xmin, Double_t &ymin, Double_t &xmax, Double_t &ymax) const
{
xmin = fXmin;
xmax = fXmax;
ymin = 0;
ymax = 0;
}
void TF1::GetRange(Double_t &xmin, Double_t &ymin, Double_t &zmin, Double_t &xmax, Double_t &ymax, Double_t &zmax) const
{
xmin = fXmin;
xmax = fXmax;
ymin = 0;
ymax = 0;
zmin = 0;
zmax = 0;
}
Double_t TF1::GetSave(const Double_t *xx)
{
if (fNsave <= 0) return 0;
if (fSave == 0) return 0;
Int_t np = fNsave - 3;
Double_t xmin = Double_t(fSave[np+1]);
Double_t xmax = Double_t(fSave[np+2]);
Double_t x = Double_t(xx[0]);
Double_t dx = (xmax-xmin)/np;
if (x < xmin || x > xmax) return 0;
if (dx <= 0) return 0;
Int_t bin = Int_t((x-xmin)/dx);
Double_t xlow = xmin + bin*dx;
Double_t xup = xlow + dx;
Double_t ylow = fSave[bin];
Double_t yup = fSave[bin+1];
Double_t y = ((xup*ylow-xlow*yup) + x*(yup-ylow))/dx;
return y;
}
TAxis *TF1::GetXaxis() const
{
TH1 *h = GetHistogram();
if (!h) return 0;
return h->GetXaxis();
}
TAxis *TF1::GetYaxis() const
{
TH1 *h = GetHistogram();
if (!h) return 0;
return h->GetYaxis();
}
TAxis *TF1::GetZaxis() const
{
TH1 *h = GetHistogram();
if (!h) return 0;
return h->GetZaxis();
}
void TF1::GradientPar(const Double_t *x, Double_t *grad, Double_t eps)
{
if(eps< 1e-10 || eps > 1) {
Warning("Derivative","parameter esp=%g out of allowed range[1e-10,1], reset to 0.01",eps);
eps = 0.01;
}
Double_t h;
TF1 *func = (TF1*)this;
Double_t *params=0;
Double_t par_local[20];
Bool_t isAllocated=kFALSE;
if (fNpar > 20){
params = new Double_t[fNpar];
isAllocated = kTRUE;
} else
params = par_local;
Bool_t errorsComputed=kFALSE;
for (Int_t ipar=0; ipar<fNpar; ipar++){
params[ipar]=fParams[ipar];
if (func->GetParError(ipar)!=0)
errorsComputed=kTRUE;
}
Double_t al, bl;
Double_t f1, f2, g1, g2, h2, d0, d2;
for (Int_t ipar=0; ipar<fNpar; ipar++){
func->InitArgs(x, params);
((TF1*)this)->GetParLimits(ipar,al,bl);
if (al*bl != 0 && al >= bl) {
grad[ipar]=0;
continue;
}
if (errorsComputed)
h = eps*func->GetParError(ipar);
else
h=eps;
params[ipar] = fParams[ipar]+h; f1 = func->EvalPar(x,params);
params[ipar] = fParams[ipar]-h; f2 = func->EvalPar(x,params);
params[ipar] = fParams[ipar]+h/2; g1 = func->EvalPar(x,params);
params[ipar] = fParams[ipar]-h/2; g2 = func->EvalPar(x,params);
h2 = 1/(2.*h);
d0 = f1 - f2;
d2 = 2*(g1 - g2);
grad[ipar] = h2*(4*d2 - d0)/3.;
params[ipar]=fParams[ipar];
}
if (isAllocated)
delete [] params;
}
void TF1::InitArgs(const Double_t *x, const Double_t *params)
{
if (fMethodCall) {
Long_t args[2];
args[0] = (Long_t)x;
if (params) args[1] = (Long_t)params;
else args[1] = (Long_t)fParams;
fMethodCall->SetParamPtrs(args);
}
}
void TF1::InitStandardFunctions()
{
TF1 *f1;
if (!gROOT->GetListOfFunctions()->FindObject("gaus")) {
f1 = new TF1("gaus","gaus",-1,1); f1->SetParameters(1,0,1);
f1 = new TF1("gausn","gausn",-1,1); f1->SetParameters(1,0,1);
f1 = new TF1("landau","landau",-1,1); f1->SetParameters(1,0,1);
f1 = new TF1("landaun","landaun",-1,1); f1->SetParameters(1,0,1);
f1 = new TF1("expo","expo",-1,1); f1->SetParameters(1,1);
for (Int_t i=0;i<10;i++) {
f1 = new TF1(Form("pol%d",i),Form("pol%d",i),-1,1);
f1->SetParameters(1,1,1,1,1,1,1,1,1,1);
}
}
}
Double_t TF1::Integral(Double_t a, Double_t b, const Double_t *params, Double_t epsilon)
{
//Begin_Html
/*
<BR>
<P>
This function computes,
to an attempted specified accuracy, the value of the integral
<P> <IMG WIDTH=140 HEIGHT=45 ALIGN=BOTTOM ALT="displaymath120" SRC="gif/gaus_img1.gif" > <P>
<p><b>Usage:</b><p>
In any arithmetic expression, this function
has the approximate value of the integral I.
<DL COMPACT>
<DT>a,b
<DD>End-points of integration interval. Note that <TT>B</TT>
may be less than <TT>A</TT>.
<DT>params
<DD>Array of function parameters. If 0, use current parameters.
<DT>epsilon
<DD>Accuracy parameter (see <B>Accuracy</B>).
</DL>
<P>
<p><b>Method:</b><p>
For any interval [<I>a</I>,<I>b</I>] we define <IMG WIDTH=49 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline128" SRC="gif/gaus_img3.gif" > and <IMG WIDTH=55 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline130" SRC="gif/gaus_img4.gif" > to be the
8-point and 16-point Gaussian quadrature approximations to
<P> <IMG WIDTH=120 HEIGHT=40 ALIGN=BOTTOM ALT="displaymath132" SRC="gif/gaus_img5.gif" > <P>
and define
<P> <IMG WIDTH=220 HEIGHT=57 ALIGN=BOTTOM ALT="displaymath134" SRC="gif/gaus_img6.gif" > <P>
Then,
<P> <IMG WIDTH=180 HEIGHT=57 ALIGN=BOTTOM ALT="displaymath138" SRC="gif/gaus_img7.gif" > <P>
<P>
where, starting with <IMG WIDTH=48 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline140" SRC="gif/gaus_img8.gif" > and finishing with <IMG WIDTH=49 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline142" SRC="gif/gaus_img9.gif" > ,
the subdivision points <IMG WIDTH=107 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline144" SRC="gif/gaus_img10.gif" > are given by
<P> <IMG WIDTH=180 HEIGHT=40 ALIGN=BOTTOM ALT="displaymath146" SRC="gif/gaus_img11.gif" > <P>
with <IMG WIDTH=8 HEIGHT=11 ALIGN=BOTTOM ALT="tex2html_wrap_inline148" SRC="gif/gaus_img12.gif" > equal to the first member of the sequence
<IMG WIDTH=66 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline150" SRC="gif/gaus_img13.gif" > for which
<IMG WIDTH=117 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline152" SRC="gif/gaus_img14.gif" > .
If, at any stage in the process of subdivision, the ratio
<P> <IMG WIDTH=180 HEIGHT=42 ALIGN=BOTTOM ALT="displaymath154" SRC="gif/gaus_img15.gif" > <P>
is so small that 1+0.005<I>q</I> is indistinguishable from 1 to
machine accuracy, an error exit occurs with the function value
set equal to zero.
<P>
<p><b>Accuracy:</b><p>
Unless there is severe cancellation of positive and negative
values of <I>f</I>(<I>x</I>) over the interval [<I>A</I>,<I>B</I>], the argument <TT>EPS</TT>
may be considered as specifying a bound on the <I>relative</I> error of
<I>I</I> in the case |<I>I</I>|>1, and a bound on the <I>absolute</I> error in
the case |<I>I</I>|<1. More precisely, if <I>k</I> is the number of sub-intervals
contributing to the approximation (see Method), and if
<P> <IMG WIDTH=180 HEIGHT=50 ALIGN=BOTTOM ALT="displaymath170" SRC="gif/gaus_img16.gif" > <P>
then the relation
<P> <IMG WIDTH=140 HEIGHT=50 ALIGN=BOTTOM ALT="displaymath172" SRC="gif/gaus_img17.gif" > <P>
will nearly always be true, provided the routine terminates
without printing an error message. For functions
<I>f</I> having no singularities in the closed interval [<I>A</I>,<I>B</I>]
the accuracy will usually be much higher than this.
<P>
<p><b>Error handling:</b><p>
The requested accuracy cannot be
obtained (see <B>Method</B>).
The function value is set equal to zero.
<P>
<p><b>Note1:</b><p>
Values of the function <I>f</I>(<I>x</I>) at the interval end-points
<I>A</I> and <I>B</I> are not required. The subprogram may therefore
be used when these values are undefined.
<BR><HR>
<p><b>Note2:</b><p>
Instead of TF1::Integral, you may want to use the combination of
<b>TF1::CalcGaussLegendreSamplingPoints</b> and <b>TF1::IntegralFast</b>.
See an example with the following script:
<pre>
void gint() {
TF1 *g = new TF1("g","gaus",-5,5);
g->SetParameters(1,0,1);
//default gaus integration method uses 6 points
//not suitable to integrate on a large domain
double r1 = g->Integral(0,5);
double r2 = g->Integral(0,1000);
//try with user directives computing more points
Int_t np = 1000;
double *x=new double[np];
double *w=new double[np];
g->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
double r3 = g->IntegralFast(np,x,w,0,5);
double r4 = g->IntegralFast(np,x,w,0,1000);
double r5 = g->IntegralFast(np,x,w,0,10000);
double r6 = g->IntegralFast(np,x,w,0,100000);
printf("g->Integral(0,5) = %g\n",r1);
printf("g->Integral(0,1000) = %g\n",r2);
printf("g->IntegralFast(n,x,w,0,5) = %g\n",r3);
printf("g->IntegralFast(n,x,w,0,1000) = %g\n",r4);
printf("g->IntegralFast(n,x,w,0,10000) = %g\n",r5);
printf("g->IntegralFast(n,x,w,0,100000)= %g\n",r6);
delete [] x;
delete [] w;
}
</pre>
<p>This example produces the following results:
<pre>
g->Integral(0,5) = 1.25331
g->Integral(0,1000) = 1.25319
g->IntegralFast(n,x,w,0,5) = 1.25331
g->IntegralFast(n,x,w,0,1000) = 1.25331
g->IntegralFast(n,x,w,0,10000) = 1.25331
g->IntegralFast(n,x,w,0,100000)= 1.253
</pre>
<BR><HR>
*/
//End_Html
const Double_t kHF = 0.5;
const Double_t kCST = 5./1000;
Double_t x[12] = { 0.96028985649753623, 0.79666647741362674,
0.52553240991632899, 0.18343464249564980,
0.98940093499164993, 0.94457502307323258,
0.86563120238783174, 0.75540440835500303,
0.61787624440264375, 0.45801677765722739,
0.28160355077925891, 0.09501250983763744};
Double_t w[12] = { 0.10122853629037626, 0.22238103445337447,
0.31370664587788729, 0.36268378337836198,
0.02715245941175409, 0.06225352393864789,
0.09515851168249278, 0.12462897125553387,
0.14959598881657673, 0.16915651939500254,
0.18260341504492359, 0.18945061045506850};
Double_t h, aconst, bb, aa, c1, c2, u, s8, s16, f1, f2;
Double_t xx[1];
Int_t i;
InitArgs(xx,params);
h = 0;
if (b == a) return h;
aconst = kCST/TMath::Abs(b-a);
bb = a;
CASE1:
aa = bb;
bb = b;
CASE2:
c1 = kHF*(bb+aa);
c2 = kHF*(bb-aa);
s8 = 0;
for (i=0;i<4;i++) {
u = c2*x[i];
xx[0] = c1+u;
f1 = EvalPar(xx,params);
if (fgAbsValue) f1 = TMath::Abs(f1);
xx[0] = c1-u;
f2 = EvalPar(xx,params);
if (fgAbsValue) f2 = TMath::Abs(f2);
s8 += w[i]*(f1 + f2);
}
s16 = 0;
for (i=4;i<12;i++) {
u = c2*x[i];
xx[0] = c1+u;
f1 = EvalPar(xx,params);
if (fgAbsValue) f1 = TMath::Abs(f1);
xx[0] = c1-u;
f2 = EvalPar(xx,params);
if (fgAbsValue) f2 = TMath::Abs(f2);
s16 += w[i]*(f1 + f2);
}
s16 = c2*s16;
if (TMath::Abs(s16-c2*s8) <= epsilon*(1. + TMath::Abs(s16))) {
h += s16;
if(bb != b) goto CASE1;
} else {
bb = c1;
if(1. + aconst*TMath::Abs(c2) != 1) goto CASE2;
h = s8;
}
return h;
}
Double_t TF1::Integral(Double_t, Double_t, Double_t, Double_t, Double_t)
{
Error("Integral","Must be called with a TF2 only");
return 0;
}
Double_t TF1::Integral(Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t)
{
Error("Integral","Must be called with a TF3 only");
return 0;
}
#ifdef INTHEFUTURE
Double_t TF1::IntegralFast(const TGraph *g, Double_t a, Double_t b, Double_t *params)
{
if (!g) return 0;
return IntegralFast(g->GetN(), g->GetX(), g->GetY(), a, b, params);
}
#endif
Double_t TF1::IntegralFast(Int_t num, Double_t *x, Double_t *w, Double_t a, Double_t b, Double_t *params)
{
if (num<=0 || x == 0 || w == 0)
return 0;
const Double_t a0 = (b + a)/2;
const Double_t b0 = (b - a)/2;
Double_t xx[1];
InitArgs(xx, params);
Double_t result = 0.0;
for (int i=0; i<num; i++)
{
xx[0] = a0 + b0*x[i];
result += w[i] * EvalPar(xx, params);
}
return result*b0;
}
Double_t TF1::IntegralMultiple(Int_t n, const Double_t *a, const Double_t *b, Double_t eps, Double_t &relerr)
{
Int_t nfnevl,ifail;
Int_t minpts = 2+2*n*(n+1)+1;
Int_t maxpts = 1000;
Double_t result = IntegralMultiple(n,a,b,minpts, maxpts,eps,relerr,nfnevl,ifail);
if (ifail > 0) {
Warning("IntegralMultiple","failed code=%d, ",ifail);
}
return result;
}
Double_t TF1::IntegralMultiple(Int_t n, const Double_t *a, const Double_t *b, Int_t minpts, Int_t maxpts, Double_t eps, Double_t &relerr,Int_t &nfnevl, Int_t &ifail)
{
//Begin_Html
/*
<img src="gif/integralmultiple.gif">
*/
//End_Html
Double_t ctr[15], wth[15], wthl[15], z[15];
const Double_t xl2 = 0.358568582800318073;
const Double_t xl4 = 0.948683298050513796;
const Double_t xl5 = 0.688247201611685289;
const Double_t w2 = 980./6561;
const Double_t w4 = 200./19683;
const Double_t wp2 = 245./486;
const Double_t wp4 = 25./729;
Double_t wn1[14] = { -0.193872885230909911, -0.555606360818980835,
-0.876695625666819078, -1.15714067977442459, -1.39694152314179743,
-1.59609815576893754, -1.75461057765584494, -1.87247878880251983,
-1.94970278920896201, -1.98628257887517146, -1.98221815780114818,
-1.93750952598689219, -1.85215668343240347, -1.72615963013768225};
Double_t wn3[14] = { 0.0518213686937966768, 0.0314992633236803330,
0.0111771579535639891,-0.00914494741655235473,-0.0294670527866686986,
-0.0497891581567850424,-0.0701112635269013768, -0.0904333688970177241,
-0.110755474267134071, -0.131077579637250419, -0.151399685007366752,
-0.171721790377483099, -0.192043895747599447, -0.212366001117715794};
Double_t wn5[14] = { 0.871183254585174982e-01, 0.435591627292587508e-01,
0.217795813646293754e-01, 0.108897906823146873e-01, 0.544489534115734364e-02,
0.272244767057867193e-02, 0.136122383528933596e-02, 0.680611917644667955e-03,
0.340305958822333977e-03, 0.170152979411166995e-03, 0.850764897055834977e-04,
0.425382448527917472e-04, 0.212691224263958736e-04, 0.106345612131979372e-04};
Double_t wpn1[14] = { -1.33196159122085045, -2.29218106995884763,
-3.11522633744855959, -3.80109739368998611, -4.34979423868312742,
-4.76131687242798352, -5.03566529492455417, -5.17283950617283939,
-5.17283950617283939, -5.03566529492455417, -4.76131687242798352,
-4.34979423868312742, -3.80109739368998611, -3.11522633744855959};
Double_t wpn3[14] = { 0.0445816186556927292, -0.0240054869684499309,
-0.0925925925925925875, -0.161179698216735251, -0.229766803840877915,
-0.298353909465020564, -0.366941015089163228, -0.435528120713305891,
-0.504115226337448555, -0.572702331961591218, -0.641289437585733882,
-0.709876543209876532, -0.778463648834019195, -0.847050754458161859};
Double_t result = 0;
Double_t abserr = 0;
ifail = 3;
nfnevl = 0;
relerr = 0;
if (n < 2 || n > 15) return 0;
Double_t twondm = TMath::Power(2,n);
Int_t ifncls = 0;
Bool_t ldv = kFALSE;
Int_t irgnst = 2*n+3;
Int_t irlcls = Int_t(twondm) +2*n*(n+1)+1;
Int_t isbrgn = irgnst;
Int_t isbrgs = irgnst;
if (minpts < 1) minpts = irlcls;
if (maxpts < minpts) maxpts = 10*minpts;
Int_t iwk = irgnst*(1 +maxpts/irlcls)/2;
Double_t *wk = new Double_t[iwk+10];
Int_t j;
for (j=0;j<n;j++) {
ctr[j] = (b[j] + a[j])*0.5;
wth[j] = (b[j] - a[j])*0.5;
}
Double_t rgnvol, sum1, sum2, sum3, sum4, sum5, difmax, f2, f3, dif, aresult;
Double_t rgncmp=0, rgnval, rgnerr;
Int_t j1, k, l, m, idvaxn=0, idvax0=0, isbtmp, isbtpp;
InitArgs(z,fParams);
L20:
rgnvol = twondm;
for (j=0;j<n;j++) {
rgnvol *= wth[j];
z[j] = ctr[j];
}
sum1 = EvalPar(z,fParams);
difmax = 0;
sum2 = 0;
sum3 = 0;
for (j=0;j<n;j++) {
z[j] = ctr[j] - xl2*wth[j];
if (fgAbsValue) f2 = TMath::Abs(EvalPar(z,fParams));
else f2 = EvalPar(z,fParams);
z[j] = ctr[j] + xl2*wth[j];
if (fgAbsValue) f2 += TMath::Abs(EvalPar(z,fParams));
else f2 += EvalPar(z,fParams);
wthl[j] = xl4*wth[j];
z[j] = ctr[j] - wthl[j];
if (fgAbsValue) f3 = TMath::Abs(EvalPar(z,fParams));
else f3 = EvalPar(z,fParams);
z[j] = ctr[j] + wthl[j];
if (fgAbsValue) f3 += TMath::Abs(EvalPar(z,fParams));
else f3 += EvalPar(z,fParams);
sum2 += f2;
sum3 += f3;
dif = TMath::Abs(7*f2-f3-12*sum1);
if (dif >= difmax) {
difmax=dif;
idvaxn=j+1;
}
z[j] = ctr[j];
}
sum4 = 0;
for (j=1;j<n;j++) {
j1 = j-1;
for (k=j;k<n;k++) {
for (l=0;l<2;l++) {
wthl[j1] = -wthl[j1];
z[j1] = ctr[j1] + wthl[j1];
for (m=0;m<2;m++) {
wthl[k] = -wthl[k];
z[k] = ctr[k] + wthl[k];
if (fgAbsValue) sum4 += TMath::Abs(EvalPar(z,fParams));
else sum4 += EvalPar(z,fParams);
}
}
z[k] = ctr[k];
}
z[j1] = ctr[j1];
}
sum5 = 0;
for (j=0;j<n;j++) {
wthl[j] = -xl5*wth[j];
z[j] = ctr[j] + wthl[j];
}
L90:
if (fgAbsValue) sum5 += TMath::Abs(EvalPar(z,fParams));
else sum5 += EvalPar(z,fParams);
for (j=0;j<n;j++) {
wthl[j] = -wthl[j];
z[j] = ctr[j] + wthl[j];
if (wthl[j] > 0) goto L90;
}
rgncmp = rgnvol*(wpn1[n-2]*sum1+wp2*sum2+wpn3[n-2]*sum3+wp4*sum4);
rgnval = wn1[n-2]*sum1+w2*sum2+wn3[n-2]*sum3+w4*sum4+wn5[n-2]*sum5;
rgnval *= rgnvol;
rgnerr = TMath::Abs(rgnval-rgncmp);
result += rgnval;
abserr += rgnerr;
ifncls += irlcls;
aresult = TMath::Abs(result);
if (ldv) {
L110:
isbtmp = 2*isbrgn;
if (isbtmp > isbrgs) goto L160;
if (isbtmp < isbrgs) {
isbtpp = isbtmp + irgnst;
if (wk[isbtmp-1] < wk[isbtpp-1]) isbtmp = isbtpp;
}
if (rgnerr >= wk[isbtmp-1]) goto L160;
for (k=0;k<irgnst;k++) {
wk[isbrgn-k-1] = wk[isbtmp-k-1];
}
isbrgn = isbtmp;
goto L110;
}
L140:
isbtmp = (isbrgn/(2*irgnst))*irgnst;
if (isbtmp >= irgnst && rgnerr > wk[isbtmp-1]) {
for (k=0;k<irgnst;k++) {
wk[isbrgn-k-1] = wk[isbtmp-k-1];
}
isbrgn = isbtmp;
goto L140;
}
L160:
wk[isbrgn-1] = rgnerr;
wk[isbrgn-2] = rgnval;
wk[isbrgn-3] = Double_t(idvaxn);
for (j=0;j<n;j++) {
isbtmp = isbrgn-2*j-4;
wk[isbtmp] = ctr[j];
wk[isbtmp-1] = wth[j];
}
if (ldv) {
ldv = kFALSE;
ctr[idvax0-1] += 2*wth[idvax0-1];
isbrgs += irgnst;
isbrgn = isbrgs;
goto L20;
}
relerr = abserr/aresult;
if (relerr < 1e-1 && aresult < 1e-20) ifail = 0;
if (relerr < 1e-3 && aresult < 1e-10) ifail = 0;
if (relerr < 1e-5 && aresult < 1e-5) ifail = 0;
if (isbrgs+irgnst > iwk) ifail = 2;
if (ifncls+2*irlcls > maxpts) {
if (sum1==0 && sum2==0 && sum3==0 && sum4==0 && sum5==0){
ifail = 0;
result = 0;
}
else
ifail = 1;
}
if (relerr < eps && ifncls >= minpts) ifail = 0;
if (ifail == 3) {
ldv = kTRUE;
isbrgn = irgnst;
abserr -= wk[isbrgn-1];
result -= wk[isbrgn-2];
idvax0 = Int_t(wk[isbrgn-3]);
for (j=0;j<n;j++) {
isbtmp = isbrgn-2*j-4;
ctr[j] = wk[isbtmp];
wth[j] = wk[isbtmp-1];
}
wth[idvax0-1] = 0.5*wth[idvax0-1];
ctr[idvax0-1] -= wth[idvax0-1];
goto L20;
}
delete [] wk;
nfnevl = ifncls;
return result;
}
Bool_t TF1::IsInside(const Double_t *x) const
{
if (x[0] < fXmin || x[0] > fXmax) return kFALSE;
return kTRUE;
}
void TF1::Paint(Option_t *option)
{
Int_t i;
Double_t xv[1];
fgCurrent = this;
TString opt = option;
opt.ToLower();
Double_t xmin=fXmin, xmax=fXmax, pmin=fXmin, pmax=fXmax;
if (gPad) {
pmin = gPad->PadtoX(gPad->GetUxmin());
pmax = gPad->PadtoX(gPad->GetUxmax());
}
if (opt.Contains("same")) {
if (xmax < pmin) return;
if (xmin > pmax) return;
if (xmin < pmin) xmin = pmin;
if (xmax > pmax) xmax = pmax;
}
TString xtitle = "";
TString ytitle = "";
char *semicol = (char*)strstr(GetTitle(),";");
if (semicol) {
Int_t nxt = strlen(semicol);
char *ctemp = new char[nxt];
strcpy(ctemp,semicol+1);
semicol = (char*)strstr(ctemp,";");
if (semicol) {
*semicol = 0;
ytitle = semicol+1;
}
xtitle = ctemp;
delete [] ctemp;
}
if (fHistogram) {
xtitle = fHistogram->GetXaxis()->GetTitle();
ytitle = fHistogram->GetYaxis()->GetTitle();
if (!gPad->GetLogx() && fHistogram->TestBit(TH1::kLogX)) { delete fHistogram; fHistogram = 0;}
if ( gPad->GetLogx() && !fHistogram->TestBit(TH1::kLogX)) { delete fHistogram; fHistogram = 0;}
}
if (fHistogram) {
fHistogram->GetXaxis()->SetLimits(xmin,xmax);
} else {
if (xmin > 0 && gPad && gPad->GetLogx()) {
Double_t *xbins = new Double_t[fNpx+1];
Double_t xlogmin = TMath::Log10(xmin);
Double_t xlogmax = TMath::Log10(xmax);
Double_t dlogx = (xlogmax-xlogmin)/((Double_t)fNpx);
for (i=0;i<=fNpx;i++) {
xbins[i] = gPad->PadtoX(xlogmin+ i*dlogx);
}
fHistogram = new TH1D("Func",GetTitle(),fNpx,xbins);
fHistogram->SetBit(TH1::kLogX);
delete [] xbins;
} else {
fHistogram = new TH1D("Func",GetTitle(),fNpx,xmin,xmax);
}
if (!fHistogram) return;
if (fMinimum != -1111) fHistogram->SetMinimum(fMinimum);
if (fMaximum != -1111) fHistogram->SetMaximum(fMaximum);
fHistogram->SetDirectory(0);
}
fHistogram->GetXaxis()->SetTitle(xtitle.Data());
fHistogram->GetYaxis()->SetTitle(ytitle.Data());
InitArgs(xv,fParams);
for (i=1;i<=fNpx;i++) {
xv[0] = fHistogram->GetBinCenter(i);
fHistogram->SetBinContent(i,EvalPar(xv,fParams));
}
Double_t minimum = fHistogram->GetMinimumStored();
Double_t maximum = fHistogram->GetMaximumStored();
if (minimum <= 0 && gPad && gPad->GetLogy()) minimum = -1111;
if (minimum == -1111) {
if (fHistogram->TestBit(TH1::kIsZoomed)) {
minimum = fHistogram->GetYaxis()->GetXmin();
} else {
minimum = fMinimum;
if (minimum == -1111) {
Double_t hmin = fHistogram->GetMinimum();
if (hmin > 0) {
Double_t hmax = fHistogram->GetMaximum();
hmin -= 0.05*(hmax-hmin);
if (hmin < 0) hmin = 0;
if (hmin <= 0 && gPad && gPad->GetLogy()) hmin = 0.001*hmax;
minimum = hmin;
}
}
}
fHistogram->SetMinimum(minimum);
}
if (maximum == -1111) {
if (fHistogram->TestBit(TH1::kIsZoomed)) {
maximum = fHistogram->GetYaxis()->GetXmax();
} else {
maximum = fMaximum;
}
fHistogram->SetMaximum(maximum);
}
fHistogram->SetBit(TH1::kNoStats);
fHistogram->SetLineColor(GetLineColor());
fHistogram->SetLineStyle(GetLineStyle());
fHistogram->SetLineWidth(GetLineWidth());
fHistogram->SetFillColor(GetFillColor());
fHistogram->SetFillStyle(GetFillStyle());
fHistogram->SetMarkerColor(GetMarkerColor());
fHistogram->SetMarkerStyle(GetMarkerStyle());
fHistogram->SetMarkerSize(GetMarkerSize());
if (!gPad) return;
if (opt.Length() == 0) fHistogram->Paint("lf");
else if (opt == "same") fHistogram->Paint("lfsame");
else fHistogram->Paint(option);
}
void TF1::Print(Option_t *option) const
{
TFormula::Print(option);
if (fHistogram) fHistogram->Print(option);
}
void TF1::ReleaseParameter(Int_t ipar)
{
if (ipar < 0 || ipar > fNpar-1) return;
SetParLimits(ipar,0,0);
}
void TF1::Save(Double_t xmin, Double_t xmax, Double_t, Double_t, Double_t, Double_t)
{
if (fSave != 0) {delete [] fSave; fSave = 0;}
fNsave = fNpx+3;
if (fNsave <= 3) {fNsave=0; return;}
fSave = new Double_t[fNsave];
Int_t i;
Double_t dx = (xmax-xmin)/fNpx;
if (dx <= 0) {
dx = (fXmax-fXmin)/fNpx;
fNsave--;
xmin = fXmin +0.5*dx;
xmax = fXmax -0.5*dx;
}
Double_t xv[1];
InitArgs(xv,fParams);
for (i=0;i<=fNpx;i++) {
xv[0] = xmin + dx*i;
fSave[i] = EvalPar(xv,fParams);
}
fSave[fNpx+1] = xmin;
fSave[fNpx+2] = xmax;
}
void TF1::SavePrimitive(ostream &out, Option_t *option )
{
Int_t i;
char quote = '"';
out<<" "<<endl;
if (!fMethodCall) {
out<<" TF1 *"<<GetName()<<" = new TF1("<<quote<<GetName()<<quote<<","<<quote<<GetTitle()<<quote<<","<<fXmin<<","<<fXmax<<");"<<endl;
if (fNpx != 100) {
out<<" "<<GetName()<<"->SetNpx("<<fNpx<<");"<<endl;
}
} else {
out<<" TF1 *"<<GetName()<<" = new TF1("<<quote<<"*"<<GetName()<<quote<<","<<fXmin<<","<<fXmax<<","<<GetNpar()<<");"<<endl;
out<<" //The original function : "<<GetTitle()<<" had originally been created by:" <<endl;
out<<" //TF1 *"<<GetName()<<" = new TF1("<<quote<<GetName()<<quote<<","<<GetTitle()<<","<<fXmin<<","<<fXmax<<","<<GetNpar()<<");"<<endl;
out<<" "<<GetName()<<"->SetRange("<<fXmin<<","<<fXmax<<");"<<endl;
out<<" "<<GetName()<<"->SetName("<<quote<<GetName()<<quote<<");"<<endl;
out<<" "<<GetName()<<"->SetTitle("<<quote<<GetTitle()<<quote<<");"<<endl;
if (fNpx != 100) {
out<<" "<<GetName()<<"->SetNpx("<<fNpx<<");"<<endl;
}
Double_t dx = (fXmax-fXmin)/fNpx;
Double_t xv[1];
InitArgs(xv,fParams);
for (i=0;i<=fNpx;i++) {
xv[0] = fXmin + dx*i;
Double_t save = EvalPar(xv,fParams);
out<<" "<<GetName()<<"->SetSavedPoint("<<i<<","<<save<<");"<<endl;
}
out<<" "<<GetName()<<"->SetSavedPoint("<<fNpx+1<<","<<fXmin<<");"<<endl;
out<<" "<<GetName()<<"->SetSavedPoint("<<fNpx+2<<","<<fXmax<<");"<<endl;
}
if (TestBit(kNotDraw)) {
out<<" "<<GetName()<<"->SetBit(TF1::kNotDraw);"<<endl;
}
if (GetFillColor() != 0) {
if (GetFillColor() > 228) {
TColor::SaveColor(out, GetFillColor());
out<<" "<<GetName()<<"->SetFillColor(ci);" << endl;
} else
out<<" "<<GetName()<<"->SetFillColor("<<GetFillColor()<<");"<<endl;
}
if (GetFillStyle() != 1001) {
out<<" "<<GetName()<<"->SetFillStyle("<<GetFillStyle()<<");"<<endl;
}
if (GetMarkerColor() != 1) {
if (GetMarkerColor() > 228) {
TColor::SaveColor(out, GetMarkerColor());
out<<" "<<GetName()<<"->SetMarkerColor(ci);" << endl;
} else
out<<" "<<GetName()<<"->SetMarkerColor("<<GetMarkerColor()<<");"<<endl;
}
if (GetMarkerStyle() != 1) {
out<<" "<<GetName()<<"->SetMarkerStyle("<<GetMarkerStyle()<<");"<<endl;
}
if (GetMarkerSize() != 1) {
out<<" "<<GetName()<<"->SetMarkerSize("<<GetMarkerSize()<<");"<<endl;
}
if (GetLineColor() != 1) {
if (GetLineColor() > 228) {
TColor::SaveColor(out, GetLineColor());
out<<" "<<GetName()<<"->SetLineColor(ci);" << endl;
} else
out<<" "<<GetName()<<"->SetLineColor("<<GetLineColor()<<");"<<endl;
}
if (GetLineWidth() != 4) {
out<<" "<<GetName()<<"->SetLineWidth("<<GetLineWidth()<<");"<<endl;
}
if (GetLineStyle() != 1) {
out<<" "<<GetName()<<"->SetLineStyle("<<GetLineStyle()<<");"<<endl;
}
if (GetChisquare() != 0) {
out<<" "<<GetName()<<"->SetChisquare("<<GetChisquare()<<");"<<endl;
}
GetXaxis()->SaveAttributes(out,GetName(),"->GetXaxis()");
GetYaxis()->SaveAttributes(out,GetName(),"->GetYaxis()");
Double_t parmin, parmax;
for (i=0;i<fNpar;i++) {
out<<" "<<GetName()<<"->SetParameter("<<i<<","<<GetParameter(i)<<");"<<endl;
out<<" "<<GetName()<<"->SetParError("<<i<<","<<GetParError(i)<<");"<<endl;
GetParLimits(i,parmin,parmax);
out<<" "<<GetName()<<"->SetParLimits("<<i<<","<<parmin<<","<<parmax<<");"<<endl;
}
if (!strstr(option,"nodraw")) {
out<<" "<<GetName()<<"->Draw("
<<quote<<option<<quote<<");"<<endl;
}
}
void TF1::SetCurrent(TF1 *f1)
{
fgCurrent = f1;
}
void TF1::SetMaximum(Double_t maximum)
{
fMaximum = maximum;
if (fHistogram) fHistogram->SetMaximum(maximum);
if (gPad) gPad->Modified();
}
void TF1::SetMinimum(Double_t minimum)
{
fMinimum = minimum;
if (fHistogram) fHistogram->SetMinimum(minimum);
if (gPad) gPad->Modified();
}
void TF1::SetNDF(Int_t ndf)
{
fNDF = ndf;
}
void TF1::SetNpx(Int_t npx)
{
if (npx < 4) {
Warning("SetNpx","Number of points must be >4 && < 100000, fNpx set to 4");
fNpx = 4;
} else if(npx > 100000) {
Warning("SetNpx","Number of points must be >4 && < 100000, fNpx set to 100000");
fNpx = 100000;
} else {
fNpx = npx;
}
Update();
}
void TF1::SetParError(Int_t ipar, Double_t error)
{
if (ipar < 0 || ipar > fNpar-1) return;
fParErrors[ipar] = error;
}
void TF1::SetParErrors(const Double_t *errors)
{
if (!errors) return;
for (Int_t i=0;i<fNpar;i++) fParErrors[i] = errors[i];
}
void TF1::SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
{
if (ipar < 0 || ipar > fNpar-1) return;
Int_t i;
if (!fParMin) {fParMin = new Double_t[fNpar]; for (i=0;i<fNpar;i++) fParMin[i]=0;}
if (!fParMax) {fParMax = new Double_t[fNpar]; for (i=0;i<fNpar;i++) fParMax[i]=0;}
fParMin[ipar] = parmin;
fParMax[ipar] = parmax;
}
void TF1::SetRange(Double_t xmin, Double_t xmax)
{
fXmin = xmin;
fXmax = xmax;
Update();
}
void TF1::SetSavedPoint(Int_t point, Double_t value)
{
if (!fSave) {
fNsave = fNpx+3;
fSave = new Double_t[fNsave];
}
if (point < 0 || point >= fNsave) return;
fSave[point] = value;
}
void TF1::SetTitle(const char *title)
{
if (!title) return;
fTitle = title;
if (!fHistogram) return;
char *semicol = (char*)strstr(title,";");
if (semicol) {
Int_t nxt = strlen(semicol);
char *ctemp = new char[nxt];
strcpy(ctemp,semicol+1);
semicol = (char*)strstr(ctemp,";");
if (semicol) {
*semicol = 0;
fHistogram->GetYaxis()->SetTitle(semicol+1);
}
fHistogram->GetXaxis()->SetTitle(ctemp);
delete [] ctemp;
}
if (gPad) gPad->Modified();
}
void TF1::Streamer(TBuffer &b)
{
if (b.IsReading()) {
UInt_t R__s, R__c;
Version_t v = b.ReadVersion(&R__s, &R__c);
if (v > 4) {
TF1::Class()->ReadBuffer(b, this, v, R__s, R__c);
if (v == 5 && fNsave > 0) {
Int_t np = fNsave - 3;
fSave[np] = fSave[np-1];
fSave[np+1] = fXmin;
fSave[np+2] = fXmax;
}
return;
}
TFormula::Streamer(b);
TAttLine::Streamer(b);
TAttFill::Streamer(b);
TAttMarker::Streamer(b);
if (v < 4) {
Float_t xmin,xmax;
b >> xmin; fXmin = xmin;
b >> xmax; fXmax = xmax;
} else {
b >> fXmin;
b >> fXmax;
}
b >> fNpx;
b >> fType;
b >> fChisquare;
b.ReadArray(fParErrors);
if (v > 1) {
b.ReadArray(fParMin);
b.ReadArray(fParMax);
} else {
fParMin = new Double_t[fNpar+1];
fParMax = new Double_t[fNpar+1];
}
b >> fNpfits;
if (v == 1) {
b >> fHistogram;
delete fHistogram; fHistogram = 0;
}
if (v > 1) {
if (v < 4) {
Float_t minimum,maximum;
b >> minimum; fMinimum =minimum;
b >> maximum; fMaximum =maximum;
} else {
b >> fMinimum;
b >> fMaximum;
}
}
if (v > 2) {
b >> fNsave;
if (fNsave > 0) {
fSave = new Double_t[fNsave+10];
b.ReadArray(fSave);
fSave[fNsave] = fSave[fNsave-1];
fSave[fNsave+1] = fSave[fNsave+2];
fSave[fNsave+2] = fSave[fNsave+3];
fNsave += 3;
} else fSave = 0;
}
b.CheckByteCount(R__s, R__c, TF1::IsA());
} else {
Int_t saved = 0;
if (fType > 0 && fNsave <= 0) { saved = 1; Save(fXmin,fXmax,0,0,0,0);}
TF1::Class()->WriteBuffer(b,this);
if (saved) {delete [] fSave; fSave = 0; fNsave = 0;}
}
}
void TF1::Update()
{
delete fHistogram;
fHistogram = 0;
if (fIntegral) {
delete [] fIntegral; fIntegral = 0;
delete [] fAlpha; fAlpha = 0;
delete [] fBeta; fBeta = 0;
delete [] fGamma; fGamma = 0;
}
}
void TF1::RejectPoint(Bool_t reject)
{
fgRejectPoint = reject;
}
Bool_t TF1::RejectedPoint()
{
return fgRejectPoint;
}
Double_t TF1_ExpValHelperx(Double_t *x, Double_t *par) {
return x[0]*gHelper->EvalPar(x,par);
}
Double_t TF1_ExpValHelper(Double_t *x, Double_t *par) {
Int_t npar = gHelper->GetNpar();
Double_t xbar = par[npar];
Double_t n = par[npar+1];
return TMath::Power(x[0]-xbar,n)*gHelper->EvalPar(x,par);
}
Double_t TF1::Moment(Double_t n, Double_t a, Double_t b, const Double_t *params, Double_t epsilon)
{
fgAbsValue = kTRUE;
Double_t norm = Integral(a,b,params,epsilon);
if (norm == 0) {
fgAbsValue = kFALSE;
Error("Moment", "Integral zero over range");
return 0;
}
gHelper = this;
TF1 fnc("TF1_ExpValHelper",TF1_ExpValHelper,fXmin,fXmax,fNpar+2);
for (Int_t i=0;i<fNpar;i++) {
if(params) fnc.SetParameter(i,params[i]);
else fnc.SetParameter(i,fParams[i]);
}
fnc.SetParameter(fNpar,0);
fnc.SetParameter(fNpar+1,n);
fgAbsValue = kFALSE;
Double_t res = fnc.Integral(a,b,params,epsilon)/norm;
return res;
}
Double_t TF1::CentralMoment(Double_t n, Double_t a, Double_t b, const Double_t *params, Double_t epsilon)
{
fgAbsValue = kTRUE;
Double_t norm = Integral(a,b,params,epsilon);
if (norm == 0) {
fgAbsValue = kFALSE;
Error("CentralMoment", "Integral zero over range");
return 0;
}
gHelper = this;
TF1 fncx("TF1_ExpValHelperx",TF1_ExpValHelperx,fXmin,fXmax,fNpar);
Int_t i;
for (i=0;i<fNpar;i++) {
if(params) fncx.SetParameter(i,params[i]);
else fncx.SetParameter(i,fParams[i]);
}
fgAbsValue = kFALSE;
Double_t xbar = fncx.Integral(a,b,params,epsilon)/norm;
TF1 fnc("TF1_ExpValHelper",TF1_ExpValHelper,fXmin,fXmax,fNpar+2);
for (i=0;i<fNpar;i++) {
if(params) fnc.SetParameter(i,params[i]);
else fnc.SetParameter(i,fParams[i]);
}
fnc.SetParameter(fNpar,0);
fnc.SetParameter(fNpar+1,n);
fnc.SetParameter(fNpar,xbar);
fnc.SetParameter(fNpar+1,n);
Double_t res = fnc.Integral(a,b,params,epsilon)/norm;
return res;
}
#ifdef INTHEFUTURE
void TF1::CalcGaussLegendreSamplingPoints(TGraph *g, Double_t eps)
{
if (!g) return;
CalcGaussLegendreSamplingPoints(g->GetN(), g->GetX(), g->GetY(), eps);
}
TGraph *TF1::CalcGaussLegendreSamplingPoints(Int_t num, Double_t eps)
{
if (num<=0)
return 0;
TGraph *g = new TGraph(num);
CalcGaussLegendreSamplingPoints(g->GetN(), g->GetX(), g->GetY(), eps);
return g;
}
#endif
void TF1::CalcGaussLegendreSamplingPoints(Int_t num, Double_t *x, Double_t *w, Double_t eps)
{
if (num<=0 || eps<=0)
return;
const UInt_t m = (num+1)/2;
Double_t z, pp, p1,p2, p3;
for (UInt_t i=0; i<m; i++) {
z = TMath::Cos(TMath::Pi()*(i+0.75)/(num+0.5));
do {
p1=1.0;
p2=0.0;
for (int j=0; j<num; j++)
{
p3 = p2;
p2 = p1;
p1 = ((2.0*j+1.0)*z*p2-j*p3)/(j+1.0);
}
pp = num*(z*p1-p2)/(z*z-1.0);
z -= p1/pp;
} while (TMath::Abs(p1/pp) > eps);
x[i] = -z;
x[num-i-1] = z;
w[i] = 2.0/((1.0-z*z)*pp*pp);
w[num-i-1] = w[i];
}
}
ROOT page - Class index - Class Hierarchy - 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.