| 1 |
* ATTENTION: in this version err_bar_inf and err_bar_sup are double precision; also all other real variables are double precision. |
| 2 |
* IMPORTANT: care has been taken that factors in expressions (before assignment to a double precision variable) are evaluated in double precision format. |
| 3 |
* NOTE: in general, double precision variables must be used for all cases where the result of an exponential or factorial is assigned. |
| 4 |
|
| 5 |
LOGICAL FUNCTION efficiency(ntot,nhit,eff,err_bar_inf,err_bar_sup) ! err_bar_inf and err_bar_sup are defined .ge.0 |
| 6 |
|
| 7 |
IMPLICIT DOUBLE PRECISION(A-H), DOUBLE PRECISION(O-Z) |
| 8 |
|
| 9 |
* NOTE: confmin,conmax and conres are ignored when the Gaussian approximation is used (simple formula with nhit and ntot only). |
| 10 |
|
| 11 |
DOUBLE PRECISION k_e_val |
| 12 |
|
| 13 |
LOGICAL choice, end_r |
| 14 |
|
| 15 |
DOUBLE PRECISION fact_r ! real to avoid integer (incorrect) result |
| 16 |
DOUBLE PRECISION p_r, s_r, r_1, r_2, f_1, f_2, e_val |
| 17 |
DOUBLE PRECISION eval_min, eval_max, eff, e_res |
| 18 |
|
| 19 |
INTEGER k_r_meas |
| 20 |
|
| 21 |
DOUBLE PRECISION e(1:2) |
| 22 |
|
| 23 |
DOUBLE PRECISION err_bar_inf, err_bar_sup |
| 24 |
|
| 25 |
* 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. |
| 26 |
* 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) |
| 27 |
|
| 28 |
if (ntot.lt.8.OR.nhit.lt.0) then |
| 29 |
* 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) |
| 30 |
|
| 31 |
eff=1.1d+0; |
| 32 |
err_bar_inf=0.d+0; |
| 33 |
err_bar_sup=0.d+0; |
| 34 |
|
| 35 |
efficiency=.false. |
| 36 |
|
| 37 |
else ! perform efficiency calculation |
| 38 |
|
| 39 |
eff=dfloat(nhit)/dfloat(ntot) ! double precision |
| 40 |
|
| 41 |
c print*, nhit, ntot, eff |
| 42 |
|
| 43 |
k_r_meas=nhit |
| 44 |
|
| 45 |
if (nhit.le.ntot-4) then ! the Gaussian approximation is possible |
| 46 |
* 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. |
| 47 |
|
| 48 |
err_bar_sup=sqrt(eff*(1.-eff)/dfloat(ntot)) |
| 49 |
|
| 50 |
err_bar_inf=err_bar_sup |
| 51 |
|
| 52 |
efficiency=.true. |
| 53 |
|
| 54 |
else ! the binomial confidence intervals must be calculated |
| 55 |
|
| 56 |
* 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 |
| 57 |
|
| 58 |
e_res=0.1/dfloat(ntot) |
| 59 |
|
| 60 |
* The upper and lower limit of the range which is checked to estimate the error bars are automatically and safely determined. |
| 61 |
|
| 62 |
* 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). |
| 63 |
* 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 |
| 64 |
* Finally, the lower limit of the range is given by: measured efficiency - (2.3+3*0.1)/ntot = measured efficiency - 2.6/ntot |
| 65 |
|
| 66 |
eval_min=eff-2.6/dfloat(ntot) |
| 67 |
|
| 68 |
* 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. |
| 69 |
|
| 70 |
eval_max=1. |
| 71 |
|
| 72 |
* NOTE: in general it is better to use k_r integer values and not r real (and hence approximate values) if possible. |
| 73 |
|
| 74 |
* for each eff. value e_ival, determine if the eff. estimate (r measured) lays inside the acceptance interval for r |
| 75 |
* NOTE: it is convenient to start scanning the e_val range from 1-e_res toward 0. |
| 76 |
* First the e_max for confidence interval will be determined , then the e_min. |
| 77 |
|
| 78 |
* Initialization |
| 79 |
|
| 80 |
i_end_ph=0 ! index of ended phase |
| 81 |
|
| 82 |
* Two phases |
| 83 |
* 1) search for e_max starting from eval_max (typically few steps) |
| 84 |
* 2) search for e_min starting from eval_min (typically many steps!!) |
| 85 |
|
| 86 |
istep=0 ! initialization of e_val step index for the starting phase |
| 87 |
|
| 88 |
do 100 while (i_end_ph .lt. 2) ! loop on e_val (two phases): exit when phase 2 (e_min search) is completed |
| 89 |
|
| 90 |
c print*, istep |
| 91 |
|
| 92 |
istep=istep+1 |
| 93 |
|
| 94 |
* Initialization: new e_val step |
| 95 |
if (i_end_ph .eq. 0) then ! phase 1 |
| 96 |
e_val=eval_max-(istep-1)*e_res ! starts from veval_max(1) (or 1) |
| 97 |
c print*, istep, e_val |
| 98 |
|
| 99 |
else ! phase 2 |
| 100 |
* 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). |
| 101 |
* In doing this it is assumed that the confidence interval is continuous, as follows from properties of binomial distribution (to be demonstrated!!!). |
| 102 |
* 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!!!) |
| 103 |
* 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. |
| 104 |
e_val=eval_min+(istep-1)*e_res |
| 105 |
c print*, istep, e_val |
| 106 |
endif |
| 107 |
|
| 108 |
if (i_end_ph .eq. 0) then |
| 109 |
if (e_val .eq. 1) then |
| 110 |
if (k_r_meas .eq. ntot) then |
| 111 |
* 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. |
| 112 |
i_end_ph=1 |
| 113 |
e(i_end_ph)=abs(e_val-eff) |
| 114 |
endif |
| 115 |
goto 150 |
| 116 |
endif |
| 117 |
if (e_val .eq. 0) then ! last e_val possible for e_max |
| 118 |
i_end_ph=1 |
| 119 |
e(i_end_ph)=abs(e_val-eff) |
| 120 |
goto 150 |
| 121 |
endif |
| 122 |
endif |
| 123 |
|
| 124 |
if (i_end_ph .eq. 1) then |
| 125 |
if (e_val .eq. 0) then |
| 126 |
if (k_r_meas .eq. 0) then |
| 127 |
* 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. |
| 128 |
i_end_ph=2 |
| 129 |
e(i_end_ph)=abs(e_val-eff) |
| 130 |
endif |
| 131 |
goto 150 |
| 132 |
endif |
| 133 |
if (e_val .eq. 1) then ! last e_val possible for e_min |
| 134 |
i_end_ph=2 |
| 135 |
e(i_end_ph)=abs(e_val-eff) |
| 136 |
goto 150 |
| 137 |
endif |
| 138 |
endif |
| 139 |
|
| 140 |
* NOTE: k_e_val is not integer in general!!! |
| 141 |
|
| 142 |
k_e_val=e_val*ntot ! real |
| 143 |
|
| 144 |
* 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). |
| 145 |
|
| 146 |
* NOTE: k_r is an integer index for possible discrete r values: k_r=r*ntot |
| 147 |
* NOTE: f(r,e_val)=e_val**(k_r)*(1-e_val)**(ntot-k_r)/(r**(k_r)*r**(ntot-k_r)) |
| 148 |
|
| 149 |
* **** determination of r_m, the new discrete value to be added to the acceptance interval |
| 150 |
* for m=1 one takes the two discrete r values r_1_1 and r_1_2 such that r_1_1<e_val<r_1_2 and evaluates f(r,e_val) to find the maximum (r_1) |
| 151 |
* for m>1, the two values are simply r_m_1 and r_m_2, such that r_m_1<r_(m-1)<r_m_2 |
| 152 |
|
| 153 |
* initializations for each fixed e_val |
| 154 |
|
| 155 |
m=1 |
| 156 |
s_r=0. |
| 157 |
end_r=.false. |
| 158 |
do 200 while (end_r .eqv. .false.) ! loop on k_r values to construct acceptance interval |
| 159 |
|
| 160 |
choice=.false. |
| 161 |
|
| 162 |
if (m .eq. 1) then |
| 163 |
|
| 164 |
* From the binomial distribution properties, it follows that r_max is one of the two discrete values of r which lay on the two sides of e_val. r_max then becomes r_1. |
| 165 |
* NOTE: by construction if the following code is executed, 0<e_val<1, so that the two r values are always such that 0 =< r =< 1 |
| 166 |
|
| 167 |
* NOTE: the difference between the two k_r value is 1 |
| 168 |
|
| 169 |
k_r_1=int(k_e_val) |
| 170 |
k_r_2=k_r_1+1 |
| 171 |
choice=.true. |
| 172 |
|
| 173 |
else |
| 174 |
|
| 175 |
* NOTE: the difference between the two k_r value is at least 2 |
| 176 |
* NOTE: for what follows it cannot be possible to be here and have both k_r_min=0 and k_r_max=1 (loop is exited after s_r=1 check). |
| 177 |
if (k_r_min .eq. 0) then |
| 178 |
k_r=k_r_max+1 |
| 179 |
k_r_max=k_r_max+1 |
| 180 |
elseif (k_r_max .eq. ntot) then |
| 181 |
k_r=k_r_min-1 |
| 182 |
k_r_min=k_r_min-1 |
| 183 |
else ! two allowed values |
| 184 |
k_r_1=k_r_min-1 |
| 185 |
k_r_2=k_r_max+1 |
| 186 |
choice=.true. |
| 187 |
endif |
| 188 |
endif |
| 189 |
|
| 190 |
* ** Choice among the two (if allowed) r values. |
| 191 |
|
| 192 |
if (choice .eqv. .true.) then |
| 193 |
|
| 194 |
* useful definitions |
| 195 |
|
| 196 |
r_1=k_r_1 ! cast into double precision |
| 197 |
r_1=r_1/ntot ! left side; pay attention to the division by integer |
| 198 |
|
| 199 |
r_2=k_r_2 ! cast into double precision |
| 200 |
r_2=r_2/ntot ! right side; pay attention to the division by integer |
| 201 |
|
| 202 |
* Calculation of the ratio of probabilities for the two possible (and allowed) r values. |
| 203 |
* NOTE: in case r_1=0 or r_2=1, then b(r_0,r_0) or b(r_1,r_1) is undefined but the limit is 1 (binomial distribution is degenerated into a step function). |
| 204 |
|
| 205 |
if (k_r_1 .eq. 0) then |
| 206 |
f_1=e_val**k_r_1*(1-e_val)**(ntot-k_r_1) ! limit |
| 207 |
else ! usual definition |
| 208 |
f_1=e_val**k_r_1*(1-e_val)**(ntot-k_r_1)/ |
| 209 |
$ (r_1**k_r_1*(1-r_1)**(ntot-k_r_1)) |
| 210 |
endif |
| 211 |
|
| 212 |
if (k_r_2 .eq. ntot) then |
| 213 |
f_2=e_val**k_r_2*(1-e_val)**(ntot-k_r_2) |
| 214 |
else |
| 215 |
f_2=e_val**k_r_2*(1-e_val)**(ntot-k_r_2)/ |
| 216 |
$ (r_2**k_r_2*(1-r_2)**(ntot-k_r_2)) |
| 217 |
endif |
| 218 |
|
| 219 |
c print*, f_1, f_2 |
| 220 |
|
| 221 |
if (f_1 .gt. f_2) then |
| 222 |
k_r=k_r_1 |
| 223 |
k_r_min=k_r_1 ! and k_r_max unchanged |
| 224 |
elseif (f_1 .lt. f_2) then |
| 225 |
k_r=k_r_2 |
| 226 |
k_r_max=k_r_2 ! and k_r_min unchanged |
| 227 |
* NOTE: f_1=f_2 case can happen also because of rounding to 6-th significant digit. |
| 228 |
else ! when equal , the k_r value nearer to k_meas is chosen |
| 229 |
if (k_e_val .lt. k_r_meas) then |
| 230 |
k_r=k_r_2 |
| 231 |
k_r_max=k_r_2 ! and k_r_min unchanged |
| 232 |
else |
| 233 |
k_r=k_r_1 |
| 234 |
k_r_min=k_r_1 ! and k_r_max unchanged |
| 235 |
endif |
| 236 |
endif |
| 237 |
|
| 238 |
if (m .eq. 1) then ! necessary to initialize the remaining acceptance interval limit |
| 239 |
if (f_1 .gt. f_2) then |
| 240 |
k_r_max=k_r_min |
| 241 |
elseif (f_1 .lt. f_2) then |
| 242 |
k_r_min=k_r_max |
| 243 |
else |
| 244 |
if (k_e_val .lt. k_r_meas) then |
| 245 |
k_r_min=k_r_max |
| 246 |
else |
| 247 |
k_r_max=k_r_min |
| 248 |
endif |
| 249 |
endif |
| 250 |
endif |
| 251 |
|
| 252 |
endif |
| 253 |
|
| 254 |
* Check if r_meas belongs to the acceptance interval, i.e. if the new value r_m is equal to r_meas |
| 255 |
|
| 256 |
if (k_r_meas .eq. k_r) then |
| 257 |
end_r=.true. |
| 258 |
i_end_ph=i_end_ph+1 |
| 259 |
if (i_end_ph.eq.1) then ! first phase: search for e_max |
| 260 |
if (e_val.lt.eff) then ! possibly caused by finite resolution in e_val scan |
| 261 |
e(i_end_ph)=0 |
| 262 |
else |
| 263 |
e(i_end_ph)=abs(e_val-eff) |
| 264 |
endif |
| 265 |
else ! second phase: search for e_min |
| 266 |
if (e_val.gt.eff) then ! possibly caused by finite resolution in e_val scan |
| 267 |
e(i_end_ph)=0 |
| 268 |
else |
| 269 |
e(i_end_ph)=abs(e_val-eff) |
| 270 |
endif |
| 271 |
endif |
| 272 |
istep=0 ! initializes step counter from next phase |
| 273 |
* Acceptance interval not completed but r_meas belongs to it, hence e_i=e_val. |
| 274 |
else ! otherwise compare the sum of probabilities against 1-alfa=0.683 |
| 275 |
|
| 276 |
* Calculate probability p(r_m,e_val) for the new value r_m. |
| 277 |
|
| 278 |
* Calculate binomial coefficient: ntot*...*(ntot - k_r + 1)/k_r! |
| 279 |
* NOTE: if k_r=0, then numerator is 1; if k_r=0, then denominator is 1. |
| 280 |
* Numerator: ntot*...*(ntot - k_r + 1): k_r terms |
| 281 |
* Denominator: k_r!: k_r terms |
| 282 |
|
| 283 |
* Probability: p(r_m,e_val)=(n_r/d_r)*(e_val**k_r)*(1-e_val)**(ntot-k_r) |
| 284 |
|
| 285 |
fact_r=1D0 ! initialization (and final value for k_r=0) |
| 286 |
c print*, fact_r |
| 287 |
fact_r=fact_r*e_val**k_r*(1-e_val)**(ntot-k_r) |
| 288 |
c print*, fact_r |
| 289 |
if (k_r .gt. 0) then |
| 290 |
do i_r=k_r,1,-1 |
| 291 |
c if(i_r.le.ntot-k_r)then |
| 292 |
c fact_r=(fact_r*e_val*(1-e_val)*(ntot-i_r+1))/(k_r-i_r+1) |
| 293 |
c else |
| 294 |
fact_r=(fact_r*(ntot-i_r+1))/(k_r-i_r+1) ! pay attention to the division by integer!! |
| 295 |
c if (i_end_ph.eq.1) then |
| 296 |
c print*, e_val, k_r, i_r, fact_r, p_r, s_r |
| 297 |
c endif |
| 298 |
c endif |
| 299 |
enddo |
| 300 |
endif |
| 301 |
|
| 302 |
c p_r=fact_r*e_val**k_r*(1-e_val)**(ntot-k_r) |
| 303 |
p_r=fact_r |
| 304 |
|
| 305 |
* Add the probability to the previous overall sum. |
| 306 |
|
| 307 |
s_r=s_r+p_r |
| 308 |
|
| 309 |
c if (i_end_ph.eq.0) then |
| 310 |
c print*, i_end_ph, k_r_meas, e_val, k_r, fact_r, p_r |
| 311 |
c endif |
| 312 |
|
| 313 |
* Compare the new sum with 1-alfa |
| 314 |
|
| 315 |
if (s_r .ge. 0.683) then |
| 316 |
end_r=.true. |
| 317 |
* End scan on r variable but continue scan on e_val variable (acceptance interval defined but does contain r_meas) |
| 318 |
endif |
| 319 |
* Otherwise continue scan on r variable to complete the acceptance interval. |
| 320 |
|
| 321 |
endif |
| 322 |
|
| 323 |
m=m+1 |
| 324 |
|
| 325 |
200 continue |
| 326 |
|
| 327 |
150 continue |
| 328 |
|
| 329 |
100 continue |
| 330 |
|
| 331 |
c if (plane_ind.eq.1) then |
| 332 |
c write(*,*) 'strip 7', vnhit_strip_sel(1,7),vngood_strip_sel(1,7) |
| 333 |
c write(*,*) eff(7), 'end' |
| 334 |
c endif |
| 335 |
|
| 336 |
err_bar_sup=e(1) |
| 337 |
err_bar_inf=e(2) |
| 338 |
|
| 339 |
efficiency=.true. |
| 340 |
|
| 341 |
endif |
| 342 |
|
| 343 |
endif |
| 344 |
|
| 345 |
c print*, e(1), e(2) |
| 346 |
|
| 347 |
end |
Parent Directory
| 
Revision Log