* ATTENTION: in this version err_bar_inf and err_bar_sup are double precision; also all other real variables are double precision. * IMPORTANT: care has been taken that factors in expressions (before assignment to a double precision variable) are evaluated in double precision format. * NOTE: in general, double precision variables must be used for all cases where the result of an exponential or factorial is assigned. LOGICAL FUNCTION efficiency(ntot,nhit,eff,err_bar_inf,err_bar_sup) ! err_bar_inf and err_bar_sup are defined .ge.0 IMPLICIT DOUBLE PRECISION(A-H), DOUBLE PRECISION(O-Z) * NOTE: confmin,conmax and conres are ignored when the Gaussian approximation is used (simple formula with nhit and ntot only). DOUBLE PRECISION k_e_val LOGICAL choice, end_r DOUBLE PRECISION fact_r ! real to avoid integer (incorrect) result DOUBLE PRECISION p_r, s_r, r_1, r_2, f_1, f_2, e_val DOUBLE PRECISION eval_min, eval_max, eff, e_res INTEGER k_r_meas DOUBLE PRECISION e(1:2) DOUBLE PRECISION err_bar_inf, err_bar_sup * VERY IMPORTANT NOTE: at present the Gaussian approximation is used also when not justified, near eff=0; in particular, with this solution, for eff=0 the error bar is by construction 0. * NOTE !! the algorithm can be improved introducing the real error bar values (when necessary, i.e. near eff=0) also for eff < 50% (at present the Gaussian approximation is automatically used for this condition, see below) if (ntot.lt.8.OR.nhit.lt.0) then * do not perform any efficiency calculation if ntot < 8 (this, together with condition nhit.le.ntot-4, assures that when efficiency is < 50% the Gaussian approximation is always used) OR nhit < 0 (non-physical values) eff=1.1d+0; err_bar_inf=0.d+0; err_bar_sup=0.d+0; efficiency=.false. else ! perform efficiency calculation eff=dfloat(nhit)/dfloat(ntot) ! double precision c print*, nhit, ntot, eff k_r_meas=nhit if (nhit.le.ntot-4) then ! the Gaussian approximation is possible * NOTE: ntot-4 is the highest number which is less than ntot/ntot*(ntot+4), for ntot.ge.12 (assumed here). This condition assures that the Gaussian approximation does not violate the physical meaning of confidence intervals, but does not guarantee that the Gaussian appproximation well reproduces the binominal confidence intervals. err_bar_sup=sqrt(eff*(1.-eff)/dfloat(ntot)) err_bar_inf=err_bar_sup efficiency=.true. else ! the binomial confidence intervals must be calculated * The needed resolution in the determination of the lower and upper limits of the confidence interval is automatically and safely determined depending on the expected size of the lower confidence error bar. For ntot=10^3 AND nhit=ntot-3 or greater, the lower bar is between 1.2 and 2.3 * 10^-3, meaning that a reasonable resolution is 1*10^-4. Since the size scales as the inverse of ntot, it can be concluded that the needed resolution is given by: 0.1 / ntot e_res=0.1/dfloat(ntot) * The upper and lower limit of the range which is checked to estimate the error bars are automatically and safely determined. * For the lower limit of the checked range, a reasonable estimate is given by: measured efficiency - highest possible value of the lower error bar - 3 times the resolution (0.1 / ntot). * Since the size of the lower error bar scales as the inverse of ntot, it can be concluded that the highest possible value of the lower error bar is given by: 2.3 / ntot * Finally, the lower limit of the range is given by: measured efficiency - (2.3+3*0.1)/ntot = measured efficiency - 2.6/ntot eval_min=eff-2.6/dfloat(ntot) * For the upper limit of the checked range, it can be shown that the size of the lower part of the range, as chosen above, is of the same order of 1 - measured efficiency (precisely, it ~86% or greater). Then for the upper limit of the range, 1 can be chosen without problems. eval_max=1. * NOTE: in general it is better to use k_r integer values and not r real (and hence approximate values) if possible. * for each eff. value e_ival, determine if the eff. estimate (r measured) lays inside the acceptance interval for r * NOTE: it is convenient to start scanning the e_val range from 1-e_res toward 0. * First the e_max for confidence interval will be determined , then the e_min. * Initialization i_end_ph=0 ! index of ended phase * Two phases * 1) search for e_max starting from eval_max (typically few steps) * 2) search for e_min starting from eval_min (typically many steps!!) istep=0 ! initialization of e_val step index for the starting phase do 100 while (i_end_ph .lt. 2) ! loop on e_val (two phases): exit when phase 2 (e_min search) is completed c print*, istep istep=istep+1 * Initialization: new e_val step if (i_end_ph .eq. 0) then ! phase 1 e_val=eval_max-(istep-1)*e_res ! starts from veval_max(1) (or 1) c print*, istep, e_val else ! phase 2 * NOTE: to avoid problems in very high (or low) number management, it is a good idea to start the scan from as near as possible to k_r_meas and with a not too high resolution (to limit the number of steps). * In doing this it is assumed that the confidence interval is continuous, as follows from properties of binomial distribution (to be demonstrated!!!). * NOTE (VERY IMPORTANT!!): the minimum (starting point) eval_min for phase 2 scan must be chosen in such a way that the lower limit results greater than this value (otherwise it means that the true confidence interval extends also below this limit!!!) * NOTE (IMPORTANT): also it is very useful in the praphs to choose the lower limit of the vertical axis to be equal to the minimum (starting point) eval_min for phase 2, and the upper limit to be greater than 1. e_val=eval_min+(istep-1)*e_res c print*, istep, e_val endif if (i_end_ph .eq. 0) then if (e_val .eq. 1) then if (k_r_meas .eq. ntot) then * If this condition is true, then e_max is automatically determined (end first phase). Otherwise e_val=1 simply does not belong to confidence interval. i_end_ph=1 e(i_end_ph)=abs(e_val-eff) endif goto 150 endif if (e_val .eq. 0) then ! last e_val possible for e_max i_end_ph=1 e(i_end_ph)=abs(e_val-eff) goto 150 endif endif if (i_end_ph .eq. 1) then if (e_val .eq. 0) then if (k_r_meas .eq. 0) then * If this condition is true, then e_min is automatically determined (end second phase). Otherwise e_val=0 simply does not belong to confidence interval. i_end_ph=2 e(i_end_ph)=abs(e_val-eff) endif goto 150 endif if (e_val .eq. 1) then ! last e_val possible for e_min i_end_ph=2 e(i_end_ph)=abs(e_val-eff) goto 150 endif endif * NOTE: k_e_val is not integer in general!!! k_e_val=e_val*ntot ! real * construction of acceptance interval (with Feldman-Cousins ordering method). In particular one must find the r value r_max that maximizes f(r,e_val)=b(r,e_val)/b(r,r). * NOTE: k_r is an integer index for possible discrete r values: k_r=r*ntot * NOTE: f(r,e_val)=e_val**(k_r)*(1-e_val)**(ntot-k_r)/(r**(k_r)*r**(ntot-k_r)) * **** determination of r_m, the new discrete value to be added to the acceptance interval * for m=1 one takes the two discrete r values r_1_1 and r_1_2 such that r_1_11, the two values are simply r_m_1 and r_m_2, such that r_m_1