--------------85C135B5E88
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
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 *----------------------------------------------------------------------*
--------------85C135B5E88 Content-Type: text/plain; charset=us-ascii; name="conlev.cc" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="conlev.cc"
/////////////////////////////////////////////////////////////////////////////// // 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 //////////////////////////////////////////////////////////////////////////
--------------85C135B5E88--