RooRarFit: rarNLL Class Reference

rarNLL::rarNLL ( )

Trivial ctor.

Usually the objects should be created using other ctors.

Definition at line 40 of file rarNLL.cc.

References init().

00041   : TNamed("",""),
00042     _xs(0), _ys(0), _mIdx(-1), _x0s(0), _x1s(0), _x2s(0),
00043     _iLIntegrals(0), _tLIntegrals(0),
00044     _verbose(kFALSE)
00045 {
00046   init();
00047 }

rarNLL::rarNLL	(	RooCurve *	curve,
		const char *	name = `"theRARNLL"`,
		const char *	title = `"the RARNLL"`,
		const Bool_t	verbose = `kFALSE`
	)

Default ctor.

Parameters:

	curve	NLL curve
	name	The name
	title	The title
	verbose	Boolean to control debug output

Default ctor to set several common data members, and then call init.

Definition at line 58 of file rarNLL.cc.

References init().

00060   : TNamed(name, title),
00061     _xs(0), _ys(0), _mIdx(-1), _x0s(0), _x1s(0), _x2s(0),
00062     _iLIntegrals(0), _tLIntegrals(0),
00063     _verbose(verbose)
00064 {
00065   init(curve);
00066 }

rarNLL::~rarNLL ( ) [virtual]

Trivial dtor.

Definition at line 69 of file rarNLL.cc.

00070 {
00071 }

rarNLL::rarNLL ( const rarNLL & ) [private]

rarNLL::ClassDef	(	rarNLL	,
		0
	)			`[private]`

Double_t rarNLL::getLIntegral ( )

Return the likelihood curve integral.

Returns:: The likelihood curve integral

It returns likelihood curve integral

Definition at line 545 of file rarNLL.cc.

Referenced by getLIntegral(), and getLIntegralInverse().

00546 {
00547   const Double_t bignumber = 1e10;
00548   return getLIntegral(0., bignumber);
00549 }

Double_t rarNLL::getLIntegral	(	Double_t	xl,
		Double_t	xh
	)

Return the likelihood curve integral.

Parameters:

	xl	The low limit for integration
	xh	The high limit for integration

Returns:: The likelihood curve integral

It returns likelihood curve integral between the limits

Definition at line 536 of file rarNLL.cc.

References getLIntegral().

00537 {
00538   return getLIntegral(xh)-getLIntegral(xl);
00539 }

Double_t rarNLL::getLIntegral ( Double_t x )

Return the likelihood curve integral.

Parameters:

x

The limit for integration

Returns:: The likelihood curve integral

It returns likelihood curve integral from -inf to the point

Definition at line 495 of file rarNLL.cc.

References _coeffMList, _iLIntegrals, _lCoefMList, _nSteps, _tLIntegrals, _verbose, _x0s, and lIntegralFunc().

Referenced by rarMLFitter::getUL().

00496 {
00497   Double_t retVal(0);
00498   const Double_t bignumber = 1e10;
00499   if (_nSteps<=0) return retVal;
00500   
00501   // first find the bin
00502   Int_t i(0);
00503   for (i=0; i<_nSteps; i++) {
00504     Double_t x0=_x0s[i];
00505     if (0==i) x0 = -bignumber;
00506     if (x==x0) {
00507       if (_verbose)
00508         cout<<" i="<<i<<" x="<<x<<" integral="<<_tLIntegrals[i]<<endl;
00509       return _tLIntegrals[i];
00510     }
00511     if (x<x0) break;
00512   }
00513   // the point falls into the previous bin
00514   i--;
00515   if (i<0) return 0; // below the lowest point
00516   // get x0, a, b, c
00517   Double_t x0=_x0s[i];
00518   TMatrixD *A=(TMatrixD*)_coeffMList.At(i);
00519   TMatrixD *lA=(TMatrixD*)_lCoefMList.At(i);
00520   if (_verbose) cout<<" i="<<i<<" x0="<<x0<<endl;
00521   // get integral for this x, a b c within the bin
00522   Double_t thisIntegral=lIntegralFunc(x, *A, *lA)-_iLIntegrals[i];
00523   retVal=_tLIntegrals[i]+thisIntegral;
00524   
00525   if (_verbose)
00526     cout<<" integral="<<retVal<<"\t(base+"<<thisIntegral<<")"<<endl;
00527   return retVal;
00528 }

Double_t rarNLL::getLIntegralInverse ( Double_t iVal )

Calculate x for given integral.

Parameters:

iVal

The integral

Returns:: The x for the given integral

It calculates x for a given integral value

Definition at line 634 of file rarNLL.cc.

References getLIntegralInverse().

00635 {
00636   return getLIntegralInverse(0., iVal);
00637 }

Double_t rarNLL::getLIntegralInverse	(	Double_t	xl,
		Double_t	iVal
	)

Calculate x for given integral.

Parameters:

	xl	The intergral starting point
	iVal	The integral

Returns:: The x for the given integral

It calculates x for a given integral value

Definition at line 557 of file rarNLL.cc.

References _coeffMList, _iLIntegrals, _lCoefMList, _nSteps, _tLIntegrals, _verbose, _x0s, _x2s, getLIntegral(), and lIntegralFuncInverse().

Referenced by getLIntegralInverse(), and rarMLFitter::getUL().

00558 {
00559   const Double_t bignumber = 1e10;
00560   if (_nSteps<=0) return 0;
00561   // increase iVal by the integral at xl
00562   iVal+=getLIntegral(xl);
00563   
00564   Double_t x=_x0s[0];
00565   // Make sure iVal is valid
00566   if (iVal<=0) return (-bignumber);
00567   if (iVal>=_tLIntegrals[_nSteps]) return (bignumber);
00568   
00569   // first identify which step iVal belongs to
00570   Int_t i(0);
00571   for (i=0; i<_nSteps; i++) {
00572     if (iVal<_tLIntegrals[i+1]) break;
00573   }
00574   
00575   // get x0, a, b, c
00576   Double_t x0=_x0s[i];
00577   TMatrixD *A=(TMatrixD*)_coeffMList.At(i);
00578   Double_t a=A->operator()(0,0);
00579   Double_t b=A->operator()(1,0);
00580   Double_t c=A->operator()(2,0);
00581   // get iVal's residual within the bin
00582   Double_t residual=iVal-_tLIntegrals[i];
00583   if (_verbose) 
00584     cout<<" i="<<i<<" x0="<<x0
00585         <<" a="<<a
00586         <<" b="<<b
00587         <<" c="<<c
00588         <<" residual="<<residual
00589         <<endl;
00590   // get integral from this bin for this x
00591   Double_t thisXI=_iLIntegrals[i]+residual;
00592   if (_verbose) cout<<" thisXI="<<thisXI<<endl;
00593   // identify integral method as in rarNLL::lIntegralFunc().
00594   // for a=b=0
00595   if ((0==a)&&(0==b)) {
00596     cout<<" Using constant fit"<<endl;
00597     x=thisXI*exp(c/2.);
00598   } else if (a>0) { // normal fit
00599     // get erf value
00600     Double_t thisXerf=thisXI*sqrt(2*a/TMath::Pi())*exp(-(b*b-4*a*c)/8/a);
00601     if ((thisXerf<-1)||(1<thisXerf)) {
00602       cout<<" "<<thisXerf<<" should be between -1 and 1"<<endl;
00603     } else {
00604       Double_t y=TMath::ErfInverse(thisXerf);
00605       x=(y*sqrt(8*a)-b)/2/a;
00606     }
00607   } else if ((0!=b)&&(a>-1e-5)&&(TMath::Abs(a/b)<100)) { // a ~ 0
00608     cout<<" Using linear fit"<<endl;
00609     Double_t thisXExp=-thisXI*b/2.;
00610     if (thisXExp<=0) {
00611       cout<<" "<<thisXExp<<" should be >0"<<endl;
00612     } else {
00613       x=-(2*TMath::Log(thisXExp)+c)/b;
00614     }
00615   } else { // numerical integral
00616     TMatrixD *lA=(TMatrixD*)_lCoefMList.At(i);
00617     Double_t la=lA->operator()(0,0);
00618     Double_t lb=lA->operator()(1,0);
00619     Double_t lc=lA->operator()(2,0);
00620     x=lIntegralFuncInverse(_x0s[i], _x2s[i], 10, thisXI, la, lb, lc);
00621     cout<<" Using numerical integral: la="<<la<<" lb="<<lb<<" lc="<<lc<<endl
00622         <<" x="<<x<<endl;
00623   }
00624   
00625   if (_verbose) cout<<" x="<<x<<endl;
00626   return x;
00627 }

TArrayD rarNLL::getMin	(	Double_t	x0,
		Double_t	x1,
		Double_t	x2,
		Double_t	a,
		Double_t	b,
		Double_t	c
	)			`[protected]`

Return the local min point for a given bin.

Parameters:

	x0	low x boundary
	x1	middle point
	x2	high x boundary
	a	2nd order coeff from fit
	b	1st order coeff from fit
	c	0th order coeff from fit

Returns:: Array of min point (x, y)

It calculates the min point on NLL curve for a given bin

Definition at line 404 of file rarNLL.cc.

References _verbose, and getMin().

00406 {
00407   const Double_t bignumber = 1e10;
00408   TArrayD xy(2);
00409   xy.Reset(bignumber);
00410   // get min from boundaries
00411   getMin(xy, x0, a, b, c);
00412   getMin(xy, x2, a, b, c);
00413   // for normal fit
00414   if (a>0) {
00415     const Double_t limit = 1e-5;
00416     Double_t xmin=-b/2/a;
00417     // check if xmin is actually x1
00418     Double_t dLimit=TMath::Abs(x2-x0) * limit;
00419     if (dLimit>limit) dLimit = limit;
00420     if (TMath::Abs(xmin-x1)<dLimit) {
00421       if (_verbose) cout<<" Set xmin from "<<xmin<<" to "<<x1<<endl;
00422       xmin=x1;
00423     }
00424     if ((x0<xmin)&&(xmin<x2)) getMin(xy, xmin, a, b, c);
00425     else if (_verbose) {
00426       Double_t ymin=-b*b/4/a+c;
00427       cout<<"  parabolic min ("<<xmin<<", "<<ymin<<")"
00428           <<" not within x ("<<x0<<", "<<x2<<")"<<endl;
00429     }
00430   }
00431   
00432   if (_verbose) cout<<"  local min: ("<<xy[0]<<", "<<xy[1]<<")"<<endl;
00433   return xy;
00434 }

void rarNLL::getMin	(	TArrayD &	xy,
		Double_t	x,
		Double_t	a,
		Double_t	b,
		Double_t	c
	)			`[protected]`

Replace array contents if y calculated is smaller.

Parameters:

	xy	Array for min (x,y)
	x	x to be evaluated
	a	2nd order coeff from fit
	b	1st order coeff from fit
	c	0th order coeff from fit

It replaces the array contents if y calculated is smaller

Definition at line 379 of file rarNLL.cc.

00380 {
00381   const Double_t bignumber = 1e10;
00382   if (xy.GetSize()<2) {
00383     xy.Set(2);
00384     xy.Reset(bignumber);
00385   }
00386   Double_t y=a*x*x+b*x+c;
00387   if (y<xy[1]) {
00388     xy[0]=x;
00389     xy[1]=y;
00390   }
00391   return;
00392 }

TArrayD rarNLL::getMin ( )

Return the min point (x,y).

Returns:: The array of min point (x,y)

It return the min point in array

Definition at line 483 of file rarNLL.cc.

Referenced by getMin().

00484 {
00485   TArrayD xy(2);
00486   getMin(xy[0], xy[1]);
00487   return xy;
00488 }

void rarNLL::getMin	(	Double_t &	x,
		Double_t &	y
	)

Return the min point.

Parameters:

	x	x var
	y	y var

It calculates the min point on NLL curve

Definition at line 441 of file rarNLL.cc.

References _coeffMList, _nSteps, _verbose, _x0s, _x1s, _x2s, and getMin().

Referenced by rarMLFitter::addErrToCurve(), rarMLFitter::combine(), rarMLFitter::combineNLLCurves(), and rarMLFitter::getMeanErrs().

00442 {
00443   if (_nSteps<=0) return;
00444   const Double_t bignumber = 1e10;
00445   x=y=bignumber;
00446   // go through bins to find min bin and min
00447   for (Int_t i=0; i<_nSteps; i++) {
00448     Double_t x0=_x0s[i];
00449     if (0==i) x0 = -bignumber;
00450     Double_t x1=_x1s[i];
00451     Double_t x2=_x2s[i];
00452     if (_nSteps-1==i) x2 = bignumber;
00453     TMatrixD *A=(TMatrixD*)_coeffMList.At(i);
00454     Double_t a=A->operator()(0,0);
00455     Double_t b=A->operator()(1,0);
00456     Double_t c=A->operator()(2,0);
00457     if (_verbose)
00458       cout<<" i="<<i
00459           <<" x0="<<x0
00460           <<" x1="<<x1
00461           <<" x2="<<x2
00462           <<" a="<<a
00463           <<" b="<<b
00464           <<" c="<<c
00465           <<endl;
00466     TArrayD xy(getMin(x0, x1, x2, a, b, c));
00467     if (2==xy.GetSize()) { // local min
00468       if (xy[1]<y) {
00469         x=xy[0];
00470         y=xy[1];
00471       }
00472     }
00473   }
00474   
00475   if (_verbose) cout<<" global min: ("<<x<<", "<<y<<")"<<endl;
00476   return;
00477 }

Double_t rarNLL::getNLL ( Double_t x )

Return NLL for a given value.

Parameters:

x

The point for NLL

Returns:: The NLL value for a given point

It returns an NLL value for a given point

Definition at line 245 of file rarNLL.cc.

References getY().

Referenced by rarMLFitter::getSignf().

00246 {
00247   return getY(x);
00248 }

TArrayD rarNLL::getX ( Double_t y )

Return x values for a given y.

Parameters:

y

NLL value

Returns:: Array having x's

It returns array of x's for a given y value

Todo:: should use weights

Definition at line 296 of file rarNLL.cc.

References _coeffMList, _nSteps, _verbose, _x0s, and _x2s.

Referenced by rarMLFitter::getMeanErrs().

00297 {
00298   TArrayD X(0); // return array
00299   if (_nSteps<=0) return X;
00300   const Double_t bignumber = 1e10;
00301   // go through bins to find x's
00302   for (Int_t i=0; i<_nSteps; i++) {
00303     TArrayD s(0); // solutions for the bin
00304     Double_t x0=_x0s[i];
00305     if (0==i) x0 = -bignumber;
00306     Double_t x2=_x2s[i];
00307     if (_nSteps-1==i) x2 = bignumber;
00308     TMatrixD *A=(TMatrixD*)_coeffMList.At(i);
00309     Double_t a=A->operator()(0,0);
00310     Double_t b=A->operator()(1,0);
00311     Double_t c=A->operator()(2,0);
00312     if ((0==a)&&(0==b)) continue;
00313     Double_t b24ac=b*b-4*a*(c-y);
00314     if (0==a) { // for a==0
00315       s.Set(1);
00316       s[0]=(y-c)/b;
00317     } else if (b24ac>=0) {
00318       s.Set(2);
00319       s[0]=(-b-sqrt(b24ac))/2/a;
00320       s[1]=(-b+sqrt(b24ac))/2/a;
00321     } else {
00322       continue;
00323     }
00324     // save solutions within the bin
00325     for (Int_t j=0; j<s.GetSize(); j++) {
00326       if ((x0<=s[j])&&(s[j]<x2)) {
00327         X.Set(X.GetSize()+1);
00328         X[X.GetSize()-1]=s[j];
00329       }
00330     }
00331   }
00332   // sort x's
00333   TArrayI I(X.GetSize());
00334   TMath::Sort(X.GetSize(), X.GetArray(), I.GetArray(), kFALSE);
00335   // combine points too close
00336   const Double_t limit = 1e-6;
00337   Double_t dLimit=(_x2s[_nSteps-1]-_x0s[0]) * limit;
00338   if (0==dLimit) dLimit = limit;
00339   if (dLimit>limit) dLimit = limit;
00340   TArrayD XX(0);
00341   if (I.GetSize()>0) {
00342     XX.Set(1);
00343     XX[0]=X[I[0]];
00344   }
00345   for (Int_t i=1; i<I.GetSize(); i++) {
00346     Double_t x=X[I[i]];
00347     Int_t iXX=XX.GetSize()-1;
00348     if (TMath::Abs(XX[iXX]-x)<dLimit) { // combine them
00349       if (_verbose) {
00350         cout<<" Merging "<<X[iXX]<<" and "<<x
00351             <<" w/ diff: "<<x-X[iXX]<<endl;
00352       }
00353       XX[iXX]=(X[iXX]+x)/2; 
00354     } else { // add x to XX
00355       XX.Set(iXX+2);
00356       XX[iXX+1]=x;
00357     }
00358   }
00359   
00360   if (_verbose) {
00361     cout<<" WRT y="<<y<<" x's are:"<<endl;
00362     for (Int_t i=0; i<XX.GetSize(); i++) {
00363       cout<<"  #"<<i<<"="<<XX[i];
00364     }
00365     cout<<endl;
00366   }
00367   
00368   return XX;
00369 }

Double_t rarNLL::getY ( Double_t x )

Return NLL (Y) for a given value.

Parameters:

x

The point for NLL

Returns:: The NLL value for a given point

It returns an NLL value for a given point

Definition at line 255 of file rarNLL.cc.

References _coeffMList, _nPoints, _nSteps, _verbose, _x2s, _xs, and _ys.

Referenced by getNLL().

00256 {
00257   Double_t retVal(0);
00258   if (_nSteps<=0) return retVal;
00259   // if exact match, return cached value
00260   for (Int_t i=0; i<_nPoints; i++) {
00261     if (x==_xs[i]) {
00262       if (_verbose) cout<<" curve point #"<<i<<" = "<<_ys[i]<<endl;
00263       return _ys[i];
00264     }
00265   }
00266   
00267   // first find the bin
00268   Int_t i;
00269   for (i=0; i<_nSteps; i++) {
00270     if (x<=_x2s[i]) break;
00271   }
00272   if (i<0) i=0;
00273   if (i>=_nSteps) i=_nSteps-1;
00274   TMatrixD *A=(TMatrixD*)_coeffMList.At(i);
00275   Double_t a=A->operator()(0,0);
00276   Double_t b=A->operator()(1,0);
00277   Double_t c=A->operator()(2,0);
00278   retVal=a*x*x+b*x+c;
00279   
00280   if (_verbose) {
00281     cout<<" i="<<i
00282         <<" a="<<a
00283         <<" b="<<b
00284         <<" c="<<c
00285         <<endl;
00286     cout<<" NLL("<<x<<")="<<retVal<<endl;
00287   }
00288   return retVal;
00289 }

void rarNLL::init ( RooCurve * curve = 0 ) [virtual]

Initial function called by ctor.

Parameters:

curve

NLL curve to use

It will initialize NLL data members based on NLL curve

Definition at line 77 of file rarNLL.cc.

References _coeffMList, _iLIntegrals, _lCoefMList, _mIdx, _nPoints, _nSteps, _tLIntegrals, _verbose, _x0s, _x1s, _x2s, _xs, _ys, and lIntegralFunc().

Referenced by rarNLL().

00078 {
00079   _nPoints=0;
00080   _nSteps=0;
00081   _coeffMList.Delete();
00082   _lCoefMList.Delete();
00083   
00084   if (!curve) {
00085     cout<<" Can not find any NLL curve"<<endl;
00086     return;
00087   }
00088   // now clone the curve
00089   curve=(RooCurve*)curve->Clone();
00090   // sort it
00091   curve->Sort();
00092   {
00093     // make sure two adjacent points are not too close
00094     Double_t dLimit=curve->GetXaxis()->GetXmax()-curve->GetXaxis()->GetXmin();
00095     dLimit=TMath::Abs(dLimit/2000.);
00096     if (0==dLimit) dLimit=1./2000.;
00097     // check points
00098     for (Int_t i=0; i<curve->GetN()-1; i++) {
00099       Double_t xlo, ylo, xhi, yhi;
00100       curve->GetPoint(i, xlo, ylo);
00101       curve->GetPoint(i+1, xhi, yhi);
00102       if (TMath::Abs(xhi-xlo)<dLimit) {
00103         Int_t rIdx=ylo<yhi?i+1:i;
00104         if (0==i) rIdx=i+1; // do not remove end point
00105         if (curve->GetN()-2==i) rIdx=i; // do not remove end point
00106         curve->RemovePoint(rIdx);
00107         if (_verbose)
00108           cout<<" Point #"<<i<<" ("<<xlo<<", "<<ylo<<") and "<<endl
00109               <<" point #"<<i+1<<" ("<<xhi<<", "<<yhi<<")"
00110               <<" are too close,"<<endl
00111               <<" remove point #"<<rIdx<<" for NLL curve calculation"<<endl;
00112         i--; // check that spot again
00113       }
00114     }
00115   }
00116   // first make sure the curve has at least 3 points
00117   _nPoints=curve->GetN();
00118   if (_nPoints<3) {
00119     cout<<" At least 3 points in the NLL curve ("<<curve->GetName()
00120         <<")are needed"<<endl;
00121     // do not need the curve any longer
00122     delete curve;
00123     return;
00124   }
00125   _nSteps=(_nPoints-1)/2;
00126   // set right size for arrays
00127   _xs.Set(_nPoints);
00128   _ys.Set(_nPoints);
00129   _x0s.Set(_nSteps);
00130   _x1s.Set(_nSteps);
00131   _x2s.Set(_nSteps);
00132   _iLIntegrals.Set(_nSteps);
00133   _tLIntegrals.Set(_nSteps+1); // The last one is the total integral
00134   
00135   // save points, and minY index
00136   _mIdx=0;
00137   for (Int_t i=0; i<_nPoints; i++) {
00138     curve->GetPoint(i, _xs[i], _ys[i]);
00139     if (_ys[i]<_ys[_mIdx]) _mIdx=i;
00140   }
00141   
00142   // calculate coeff matrics, integrals, etc.
00143   for (Int_t i=0; i<_nSteps; i++) {
00144     TMatrixD X(3,3), Y(3, 1, &_ys[2*i]), lY(3,1);
00145     // lY is y's in Likelihood space
00146     lY(0,0)=exp(-.5*Y(0,0));
00147     lY(1,0)=exp(-.5*Y(1,0));
00148     lY(2,0)=exp(-.5*Y(2,0));
00149     
00150     _x0s[i]=_xs[2*i+0];
00151     _x1s[i]=_xs[2*i+1];
00152     _x2s[i]=_xs[2*i+2];
00153     
00154     X(0,0)=_x0s[i]*_x0s[i];
00155     X(0,1)=_x0s[i];
00156     X(0,2)=1;
00157     X(1,0)=_x1s[i]*_x1s[i];
00158     X(1,1)=_x1s[i];
00159     X(1,2)=1;
00160     X(2,0)=_x2s[i]*_x2s[i];
00161     X(2,1)=_x2s[i];
00162     X(2,2)=1;
00163     
00164     // invert X
00165     TMatrixD invX=TMatrixD(TMatrixD::kInverted, X);
00166     // get coeff matrix
00167     TMatrixD *A=new TMatrixD(invX*Y);
00168     _coeffMList.Add(A);
00169     // coeff matrix in likelihood space
00170     TMatrixD *lA=new TMatrixD(invX*lY);
00171     _lCoefMList.Add(lA);
00172     
00173     // calculate integral
00174     Double_t a=A->operator()(0,0);
00175     Double_t b=A->operator()(1,0);
00176     Double_t c=A->operator()(2,0);
00177     // check if a <= 0
00178     if (a<=0) {
00179       cout<<" W A R N I N G !"<<endl
00180           <<" The 3-point fit on (X, Y)"<<endl
00181           <<" ("<<X(0,1)<<", "<<Y(0,0)<<")"<<endl
00182           <<" ("<<X(1,1)<<", "<<Y(1,0)<<")"<<endl
00183           <<" ("<<X(2,1)<<", "<<Y(2,0)<<")"<<endl
00184           <<" returns a="<<a<<" (<=0)"
00185           <<" b="<<b<<" c="<<c<<endl
00186           <<" Please make sure the scanning fits around there are normal"
00187           <<endl;
00188       if (i<=0) { // do not allow
00189         curve->RemovePoint(0);
00190         init(curve);
00191         delete curve;
00192         return;
00193       }
00194       if (i>=_nSteps-1) { // do not allow
00195         curve->RemovePoint(_nPoints-1);
00196         init (curve);
00197         delete curve;
00198         return;
00199       }
00200     }
00201     
00202     Double_t x=_x0s[i];
00203     const Double_t bignumber = 1e10;
00204     if (i<=0) x = -bignumber;
00205     // get lower integral for this step
00206     _iLIntegrals[i]=lIntegralFunc(x, *A, *lA);
00207     x=_x2s[i];
00208     if (i>=_nSteps-1) x = bignumber;
00209     // get integral for this step
00210     Double_t thisIntegral=lIntegralFunc(x, *A, *lA)-_iLIntegrals[i];
00211     if (thisIntegral<0) {
00212       cout<<" Negative integral="<<thisIntegral
00213           <<" for i="<<i<<" x0="<<_x0s[i]
00214           <<" x1="<<_x1s[i]<<" x2="<<_x2s[i]<<endl;
00215       thisIntegral=0;
00216     }
00217     // get total integral for next step
00218     if (i<=0) _tLIntegrals[0]=0;
00219     _tLIntegrals[i+1]=_tLIntegrals[i]+thisIntegral;
00220   }
00221   // printout
00222   if (_verbose) {
00223     cout<<" rarNLL based on:"<<endl;
00224     curve->Print("v");
00225     _coeffMList.Print();
00226     cout<<" "<<_nSteps<<" steps for integral calculation:"<<endl;
00227     for (Int_t i=0; i<_nSteps; i++) {
00228       cout<<" i="<<i<<"\t x0="<<_x0s[i]
00229           <<"\t I="<<_iLIntegrals[i]
00230           <<"\t T="<<_tLIntegrals[i]
00231           <<endl;
00232     }
00233   }
00234   
00235   // do not need the curve any longer
00236   delete curve;
00237   return;
00238 }

Double_t rarNLL::lIntegralFunc	(	Double_t	x,
		TMatrixD &	A,
		TMatrixD &	lA
	)

The likelihood curve integral function.

Parameters:

	x	x value
	A	Coeff matrix
	lA	Coeff matrix in likelihood space

It is the likelihood integral function based on 2NLL curve parameters The formula is based on exp(-y/2) where y is from 2NLL curve

/// Integrate[Exp[(-a x x - b x - c)/2],x]
/// 
///   2
///  b /(8 a) - c/2      Pi          b + 2 a x
/// E               Sqrt[--] Erf[-----------------]
///                      2       2 Sqrt[2] Sqrt[a]
/// -----------------------------------------------
///                     Sqrt[a]
///

If a<0, the analytical integral is a function of imaginary erf

/// Integrate[Exp[(a x x - b x - c)/2],x]
///
///    2
///  -b /(8 a) - c/2      Pi          -b + 2 a x
/// E                Sqrt[--] Erfi[-----------------]
///                       2        2 Sqrt[2] Sqrt[a]
/// -------------------------------------------------
///                      Sqrt[a]
///

We prefer to use numerical integral in likelihood space for a<0 because of the lack of Erfi and InverseErfi in ROOT/RooFit and the complicated nature of this function. If a==0, we use analytical integral

/// Integrate[Exp[(- b x - c)/2], x]
///
///     -c/2 - (b x)/2
/// -2 E
/// ------------------
///         b
///

If a==0 and b==0, we use

/// Integrate[Exp[(- c)/2], x]
///
///  x
/// ----
///  c/2
/// E
///

http://mathworld.wolfram.com/Erf.html
http://mathworld.wolfram.com/Erfi.html

Definition at line 691 of file rarNLL.cc.

References _verbose.

Referenced by getLIntegral(), and init().

00692 {
00693   Double_t retVal(0);
00694   Double_t a=A(0,0);
00695   Double_t b=A(1,0);
00696   Double_t c=A(2,0);
00697   // for a=b=0
00698   if ((0==a)&&(0==b)) {
00699     cout<<" Using constant fit"<<endl;
00700     retVal=x/exp(c/2.);
00701   } else if (a>0) { // normal fit
00702     retVal=sqrt(TMath::Pi()/2/a)*exp((b*b-4*a*c)/8/a)*
00703       TMath::Erf((b+2*a*x)/sqrt(8*a));
00704   } else if ((0!=b)&&(a>-1e-5)&&(TMath::Abs(a/b)<100)) { // a ~ 0
00705     cout<<" Using linear fit"<<endl;
00706     retVal=-2/b*exp((-b*x-c)/2.);
00707   } else { // numerical integral
00708     Double_t la=lA(0,0);
00709     Double_t lb=lA(1,0);
00710     Double_t lc=lA(2,0);
00711     retVal=la/3*x*x*x+lb/2*x*x+lc*x;
00712     cout<<" Using numerical integral: la="<<la<<" lb="<<lb<<" lc="<<lc<<endl
00713         <<" x="<<x<<" retV:"<<retVal<<endl;
00714   }
00715   const Double_t bignumber = 1e200;
00716   if (retVal>bignumber)  { retVal = bignumber; }
00717   if (retVal<-bignumber) { retVal = -bignumber; }
00718   if (_verbose)
00719     cout<<"x:"<<x<<" a:"<<a<<" b:"<<b<<" c:"<<c<<" retV LIFunc:"<<retVal<<endl;
00720   return retVal;
00721 }

Double_t rarNLL::lIntegralFuncInverse	(	Double_t	x0,
		Double_t	x2,
		Int_t	iter,
		Double_t &	thisXI,
		Double_t &	la,
		Double_t &	lb,
		Double_t &	lc
	)

Calculate x for given integral using numerical evaluation.

Parameters:

	x0	Lower limit
	x2	Upper limit
	iter	Number of iterations
	thisXI	The integral
	la	Coeff la
	lb	Coeff lb
	lc	Coeff lc

Returns:: The x for the given integral

It calculates x for a given integral value using numerical evaluation (divide and conquer)

Definition at line 735 of file rarNLL.cc.

Referenced by getLIntegralInverse().

00738 {
00739   Double_t x=(x0+x2)/2.;
00740   if (iter--<=0) return x;
00741   Double_t theIntegral=la/3*x*x*x+lb/2*x*x+lc*x;
00742   if (theIntegral==thisXI) return x;
00743   if (theIntegral<thisXI)
00744     return lIntegralFuncInverse(x, x2, iter, thisXI, la, lb, lc);
00745   return lIntegralFuncInverse(x0, x, iter, thisXI, la, lb, lc);
00746 }

virtual void rarNLL::setVerbose ( Bool_t verbose = kTRUE ) [inline, virtual]

Set verbose mode.

Parameters:

verbose

Verbose mode

Definition at line 35 of file rarNLL.hh.

References _verbose.

00035 {_verbose=verbose;}

TList rarNLL::_coeffMList [protected]

List of coeff matrices.

Definition at line 63 of file rarNLL.hh.

Referenced by getLIntegral(), getLIntegralInverse(), getMin(), getX(), getY(), and init().

TArrayD rarNLL::_iLIntegrals [protected]

Lower likelihood integrals within each step.

Definition at line 68 of file rarNLL.hh.

Referenced by getLIntegral(), getLIntegralInverse(), and init().

TList rarNLL::_lCoefMList [protected]

List of coeff matrices for likelihood curve fit.

Definition at line 64 of file rarNLL.hh.

Referenced by getLIntegral(), getLIntegralInverse(), and init().

Int_t rarNLL::_mIdx [protected]

NLL point index for minY.

Definition at line 61 of file rarNLL.hh.

Referenced by init().

Int_t rarNLL::_nPoints [protected]

Number of points in NLL curve.

Definition at line 58 of file rarNLL.hh.

Referenced by getY(), and init().

Int_t rarNLL::_nSteps [protected]

Number of calculation steps.

Definition at line 62 of file rarNLL.hh.

Referenced by getLIntegral(), getLIntegralInverse(), getMin(), getX(), getY(), and init().

TArrayD rarNLL::_tLIntegrals [protected]

Total likelihood integrals before each step.

Definition at line 69 of file rarNLL.hh.

Referenced by getLIntegral(), getLIntegralInverse(), and init().

Bool_t rarNLL::_verbose [protected]

Boolean to control debug output.

Definition at line 71 of file rarNLL.hh.

Referenced by getLIntegral(), getLIntegralInverse(), getMin(), getX(), getY(), init(), lIntegralFunc(), and setVerbose().

TArrayD rarNLL::_x0s [protected]

Starting point for each step.

Definition at line 65 of file rarNLL.hh.

Referenced by getLIntegral(), getLIntegralInverse(), getMin(), getX(), and init().

TArrayD rarNLL::_x1s [protected]

Middle point for each step.

Definition at line 66 of file rarNLL.hh.

Referenced by getMin(), and init().

TArrayD rarNLL::_x2s [protected]

Ending point for each step.

Definition at line 67 of file rarNLL.hh.

Referenced by getLIntegralInverse(), getMin(), getX(), getY(), and init().

TArrayD rarNLL::_xs [protected]

NLL curve x values.

Definition at line 59 of file rarNLL.hh.

Referenced by getY(), and init().

TArrayD rarNLL::_ys [protected]

NLL curve y values.

Definition at line 60 of file rarNLL.hh.

Referenced by getY(), and init().

rarNLL Class Reference

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	rarNLL ()
	Trivial ctor.
	rarNLL (RooCurve curve, const char name="theRARNLL", const char *title="the RARNLL", const Bool_t verbose=kFALSE)
	Default ctor.
virtual	~rarNLL ()
	Trivial dtor.
virtual void	setVerbose (Bool_t verbose=kTRUE)
	Set verbose mode.
virtual void	init (RooCurve *curve=0)
	Initial function called by ctor.
Double_t	getNLL (Double_t x)
	Return NLL for a given value.
Double_t	getY (Double_t x)
	Return NLL (Y) for a given value.
TArrayD	getX (Double_t y)
	Return x values for a given y.
void	getMin (Double_t &x, Double_t &y)
	Return the min point.
TArrayD	getMin ()
	Return the min point (x,y).
Double_t	getLIntegral (Double_t x)
	Return the likelihood curve integral.
Double_t	getLIntegral (Double_t xl, Double_t xh)
	Return the likelihood curve integral.
Double_t	getLIntegral ()
	Return the likelihood curve integral.
Double_t	getLIntegralInverse (Double_t xl, Double_t iVal)
	Calculate x for given integral.
Double_t	getLIntegralInverse (Double_t iVal)
	Calculate x for given integral.
Double_t	lIntegralFunc (Double_t x, TMatrixD &A, TMatrixD &lA)
	The likelihood curve integral function.
Double_t	lIntegralFuncInverse (Double_t x0, Double_t x2, Int_t iter, Double_t &thisXI, Double_t &la, Double_t &lb, Double_t &lc)
	Calculate x for given integral using numerical evaluation.
Protected Member Functions
void	getMin (TArrayD &xy, Double_t x, Double_t a, Double_t b, Double_t c)
	Replace array contents if y calculated is smaller.
TArrayD	getMin (Double_t x0, Double_t x1, Double_t x2, Double_t a, Double_t b, Double_t c)
	Return the local min point for a given bin.
Protected Attributes
Int_t	_nPoints
	Number of points in NLL curve.
TArrayD	_xs
	NLL curve x values.
TArrayD	_ys
	NLL curve y values.
Int_t	_mIdx
	NLL point index for minY.
Int_t	_nSteps
	Number of calculation steps.
TList	_coeffMList
	List of coeff matrices.
TList	_lCoefMList
	List of coeff matrices for likelihood curve fit.
TArrayD	_x0s
	Starting point for each step.
TArrayD	_x1s
	Middle point for each step.
TArrayD	_x2s
	Ending point for each step.
TArrayD	_iLIntegrals
	Lower likelihood integrals within each step.
TArrayD	_tLIntegrals
	Total likelihood integrals before each step.
Bool_t	_verbose
	Boolean to control debug output.
Private Member Functions
	rarNLL (const rarNLL &)
	ClassDef (rarNLL, 0)