Hi Nick, I have added a few new functions in the class TMath, including TMath::Probability TMath::Landau where Probability is based on your function below. Thanks for submitting the code. 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 > //////////////////////////////////////////////////////////////////////////
This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:34:39 MET