Rene Brun
Nick van Eijndhoven wrote:
>
> 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
> //////////////////////////////////////////////////////////////////////////