Dear fellow ROOTers,
In performing some physics analysis I needed to compute the
confidence level based on a chi-squared and number of statistical
constraints (i.e. the functionality of the old cernlib routine prob).
Since I couldn't find such a facility in the ROOT framework (in case
there is already one, please let me know), I decided to dig into
my old private fortran utilities.
Below you find attached the function "conlev" which is a C++
translation of the old fortran code (the original code is also
included as comment) which does the job (I tested the reproduction
of the CL values as noted in the particle data book.
Note : the result for 1 constraint is incorrect, but normally one
doesn't encounter this situation).
I know it won't win a prize in a contest for nice programming,
but at least it works.
Feel free to use, modifiy etc.... the code and in case the ROOT
team wants to include it in the standard ROOT distribution, just
go ahead with it.
Curently the easiest way I found to have it always available within
ROOT is to just include the C++ code in the file rootalias.C.
--
Cheers,
_/_/ _/ _/ _/_/_/_/ _/ _/
_/ _/ _/ _/ _/ _/ _/
_/ _/ _/ _/ _/ _/_/
_/ _/_/ _/ _/ _/ _/
_/ _/ _/ _/_/_/_/ _/ _/
*----------------------------------------------------------------------*
Dr. Nick van Eijndhoven Department of Subatomic Physics
email : nick@phys.uu.nl Utrecht University / NIKHEF
tel. +31-30-2532331 (direct) P.O. Box 80.000
tel. +31-30-2531492 (secr.) NL-3508 TA Utrecht
fax. +31-30-2518689 The Netherlands
WWW : http://www.phys.uu.nl/~nick Office : Ornstein lab. 172
----------------------------------------------------------------------
tel. +41-22-7679751 (direct) CERN PPE Division / ALICE exp.
tel. +41-22-7675857 (secr.) CH-1211 Geneva 23
fax. +41-22-7679480 Switzerland
CERN beep : 13+7294 Office : B 160 1-012
*----------------------------------------------------------------------*
///////////////////////////////////////////////////////////////////////////////
// Computation of the confidence level from Chi-squared (chi2)
// and number of constraints (ncons).
//
// Note :
// for even ncons ==> same result as the Cernlib function PROB
// for odd ncons ==> confidence level < result of the Cernlib function PROB
//
//--- Original Fortan routine copied from Lau Gatignon 1980
//--- Modified by NvE 29-sep-1980 K.U. Nijmegen
//--- Modified by NvE 24-apr-1985 NIKHEF-H Amsterdam
//--- Converted to C++ by Nve 06-nov-1998 UU-SAP Utrecht
///////////////////////////////////////////////////////////////////////////////
float conlev(float chi2,int ncons)
{
const float a1=0.705230784e-1, a2=0.422820123e-1,
a3=0.92705272e-2, a4=0.1520143e-3,
a5=0.2765672e-3, a6=0.430638e-4;
// Change to dummy variables
int n=ncons;
float y0=chi2;
// Set CL to zero in case ncons<=0
if (n <= 0) return 0;
if (y0 <= 0.)
{
if (y0 < 0.)
{
return 0;
}
else
{
return 1;
}
}
if (n > 100)
{
float x=sqrt(2.*y0)-sqrt(float(n+n-1));
float t=fabs(x)/1.1412135;
float denom=1.0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6)))));
float vfun=1.0/denom
if (fabs(vfun) < 1.3e-5) vfun=0.; // Prevent underflow
vfun=pow(vfun,16);
float v=0.5*vfun;
if (x < 0.) v=1.-v;
if (v < 0.) v=0.;
return v;
}
else
{
float y1=0.5*y0;
int m=n/2;
if (2*m == n)
{
float sum=1.;
if (m > 1)
{
float term=1.;
for (int i=2; i<=m; i++)
{
if (term > 1.e20) break; // Prevent overflow
float fi=float(i-1);
term=term*y1/fi;
sum+=term;
}
}
float rnick=y1;
if (rnick >= 1.75e2) rnick=1.75e2; // Prevent underflow
float v=sum*exp(-rnick);
if (v < 0.) v=0.;
return v;
}
else
{
float x=sqrt(y0);
float t=fabs(x)/1.1412135;
float denom=1.0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6)))));
float vfun=1.0/denom
if (fabs(vfun) < 1.3e-5) vfun=0.; // Prevent underflow
vfun=pow(vfun,16);
float v=vfun;
if (n < 1) return 0;
if (n == 1)
{
if (v < 0.) v=0.;
return v;
}
float value=v;
float sum=1.;
if (n >= 4)
{
float term=1.;
int k=m-1;
for (int j=1; j<=k; j++)
{
if (term > 1.e20) break; // Prevent overflow
float fj=float(j);
term=term*y0/(fj+fj+1.);
sum+=term;
}
}
if (y1 > 1.75e2) y1=1.75e2; // Prevent underflow
float vi=sum*0.797846*sqrt(y0)*exp(-y1);
v=vi+value;
if (v < 0.) v=0.;
return v;
}
}
}
///////////////////////////////////////////////////////////////////////////////
// SUBROUTINE CONLEV (CL,NCONS,CHI2)
//C *** ROUTINE COPIED FROM LAU GATIGNON ***
//C *** MODIFIED BY NVE 29-SEPT-1980 K.U. NIJMEGEN ***
//C --- LAST MOD. NVE 24-APR-1985 NIKHEF-H AMSTERDAM ---
//C
//C --- THE COMPUTATION OF THE CONFIDENCE LEVEL (CL) FROM CHI-SQAURED (CHI2) ---
//C --- AND NUMBER OF CONSTRAINTS (NCONS) ---
//C NOTE ;
//C FOR EVEN NCONS ==> SAME RESULT AS THE CERN FUNCTION PROB
//C FOR ODD NCONS ==> CL < RESULT OF THE CERN FUNCTION PROB
//C
//* DATA A1 /0.705230784E-1/, A2 /0.422820123E-1/
//* DATA A3 /0.92705272 E-2/, A4 /0.1520143 E-3/
//* DATA A5 /0.2765672 E-3/, A6 /0.430638 E-4/
//*C
//*C --- CHANGE TO DUMMY VARIABLES ---
//* N=NCONS
//* Y0=CHI2
//C
//*C --- SET CL TO ZERO IN CASE OF NCONS .LE. 0 ---
//* IF (N .GT. 0) GO TO 1
//* V=0.0
//* GO TO 130
//*C
//* 1 IF (Y0 .GT. 0.0) GO TO 10
//* V=0.0
//* IF (Y0 .EQ. 0.0) V=1.0
//* GO TO 130
//*C
//* 10 IF (N .LT. 101) GO TO 30
//* X=SQRT(Y0+Y0)-SQRT(FLOAT(N+N-1))
//* ASSIGN 20 TO JUMP
//* GO TO 140
//*C
//* 20 V=0.5*VFUN
//* IF (X .LT. 0.0) V=1.0-V
//* GO TO 110
//*C
//* 30 Y1=0.5*Y0
//* M = N/2
//* IF (M+M .NE. N) GO TO 60
//* SUM=1.0
//* IF (M .LE. 1) GO TO 50
//* TERM=1.0
//* DO 40 I=2,M
//*C
//*C *** PREVENT OVERFLOW ***
//* IF (TERM .GT. 1.0E20) GO TO 50
//* FI=I-1
//* TERM=TERM*Y1/FI
//* 40 SUM=SUM+TERM
//*C *** PREVENT UNDERFLOW ***
//* 50 RNICK=Y1
//* IF (RNICK .GE. 1.75E2) RNICK=1.75E2
//* V=SUM*EXP(-RNICK)
//* GO TO 110
//C
//* 60 X=SQRT(Y0)
//* ASSIGN 70 TO JUMP
//* GO TO 140
//*C
//* 70 V=VFUN
//* IF (N-1) 120,110,80
//* 80 VALUE=V
//* SUM=1.0
//* IF (N .LT. 4) GO TO 100
//* TERM=1.0
//* K=M-1
//* DO 90 I=1,K
//* IF (TERM .GT. 1.0E20) GO TO 100
//* FI=I
//* TERM=TERM*Y0/(FI+FI+1.0)
//* 90 SUM=SUM+TERM
//*C *** PREVENT UNDERFLOW ***
//* 100 IF (Y1 .GT. 1.75E2) Y1=1.75E2
//* VI=SUM*0.797846*SQRT(Y0)*EXP(-Y1)
//* V=VI+VALUE
//* 110 IF (V) 120,130,130
//* 120 V=0.0
//* 130 CONTINUE
//* CL=V
//* RETURN
//*C
//* 140 T=ABS(X)/1.1412135
//* DENOM=1.0+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+T*A6)))))
//* VFUN=1.0/DENOM
//*C
//*C *** PREVENT UNDERFLOW ***
//* IF (ABS(VFUN) .LT. 1.3E-5) VFUN=0.0
//* VFUN=VFUN**16
//* GO TO JUMP,(20,70)
//*C
//* END
//////////////////////////////////////////////////////////////////////////
This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:34:39 MET