*     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_1<e_val<r_1_2 and evaluates f(r,e_val) to find the maximum (r_1)
*     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

*     initializations for each fixed e_val

            m=1
            s_r=0.
            end_r=.false.
            do 200 while (end_r .eqv. .false.) ! loop on k_r values to construct acceptance interval

	      choice=.false.

              if (m .eq. 1) then

*     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.
*     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

*     NOTE: the difference between the two k_r value is 1

	        k_r_1=int(k_e_val)
	        k_r_2=k_r_1+1
	        choice=.true.
	
	      else 

*     NOTE: the difference between the two k_r value is at least 2
*     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).
                if (k_r_min .eq. 0) then
	          k_r=k_r_max+1
	          k_r_max=k_r_max+1
                elseif (k_r_max .eq. ntot) then
	          k_r=k_r_min-1
	          k_r_min=k_r_min-1
                else ! two allowed values	  
	          k_r_1=k_r_min-1
	          k_r_2=k_r_max+1
	          choice=.true.
	        endif
	      endif

*     ** Choice among the two (if allowed) r values.
	
	      if (choice .eqv. .true.) then

*     useful definitions
	
	        r_1=k_r_1 ! cast into double precision
	        r_1=r_1/ntot ! left side; pay attention to the division by integer
	  
	        r_2=k_r_2 ! cast into double precision
	        r_2=r_2/ntot ! right side; pay attention to the division by integer
	
*     Calculation of the ratio of probabilities for the two possible (and allowed) r values.
*     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).

	        if (k_r_1 .eq. 0) then
	          f_1=e_val**k_r_1*(1-e_val)**(ntot-k_r_1) ! limit
	        else ! usual definition
	          f_1=e_val**k_r_1*(1-e_val)**(ntot-k_r_1)/
     $            (r_1**k_r_1*(1-r_1)**(ntot-k_r_1))
	        endif
	  
	        if (k_r_2 .eq. ntot) then
	          f_2=e_val**k_r_2*(1-e_val)**(ntot-k_r_2)
	        else
	          f_2=e_val**k_r_2*(1-e_val)**(ntot-k_r_2)/
     $	          (r_2**k_r_2*(1-r_2)**(ntot-k_r_2))	
	        endif

c	print*, f_1, f_2

	        if (f_1 .gt. f_2) then
	          k_r=k_r_1
	          k_r_min=k_r_1 ! and k_r_max unchanged
	        elseif (f_1 .lt. f_2) then
	          k_r=k_r_2
	          k_r_max=k_r_2 ! and k_r_min unchanged
*     NOTE: f_1=f_2 case can happen also because of rounding to 6-th significant digit.
	        else ! when equal , the k_r value nearer to k_meas is chosen
	          if (k_e_val .lt. k_r_meas) then
	            k_r=k_r_2
	            k_r_max=k_r_2 ! and k_r_min unchanged
	          else
	            k_r=k_r_1
	            k_r_min=k_r_1 ! and k_r_max unchanged
	          endif
	        endif
	
	        if (m .eq. 1) then ! necessary to initialize the remaining acceptance interval limit
	          if (f_1 .gt. f_2) then
	            k_r_max=k_r_min
	          elseif (f_1 .lt. f_2) then
	            k_r_min=k_r_max
	          else
	            if (k_e_val .lt. k_r_meas) then
	              k_r_min=k_r_max
	            else
	              k_r_max=k_r_min
	            endif
	          endif	  
	        endif
	  
	      endif
	
*     Check if r_meas belongs to the acceptance interval, i.e. if the new value r_m is equal to r_meas

	      if (k_r_meas .eq. k_r) then
	        end_r=.true.
	        i_end_ph=i_end_ph+1
	        if (i_end_ph.eq.1) then ! first phase: search for e_max
	          if (e_val.lt.eff) then ! possibly caused by finite resolution in e_val scan
	            e(i_end_ph)=0
	          else
	            e(i_end_ph)=abs(e_val-eff)
	          endif
	        else ! second phase: search for e_min
	          if (e_val.gt.eff) then ! possibly caused by finite resolution in e_val scan
	            e(i_end_ph)=0
	          else
	            e(i_end_ph)=abs(e_val-eff)
	          endif
	        endif
	        istep=0 ! initializes step counter from next phase
*     Acceptance interval not completed but r_meas belongs to it, hence e_i=e_val.
	      else ! otherwise compare the sum of probabilities against 1-alfa=0.683
	
*     Calculate probability p(r_m,e_val) for the new value r_m.

*     Calculate binomial coefficient: ntot*...*(ntot - k_r + 1)/k_r!
*     NOTE: if k_r=0, then numerator is 1; if k_r=0, then denominator is 1.
*     Numerator: ntot*...*(ntot - k_r + 1): k_r terms
*     Denominator: k_r!: k_r terms

*     Probability: p(r_m,e_val)=(n_r/d_r)*(e_val**k_r)*(1-e_val)**(ntot-k_r)
	
	      fact_r=1D0 ! initialization (and final value for k_r=0)
c	  print*, fact_r
	      fact_r=fact_r*e_val**k_r*(1-e_val)**(ntot-k_r)
c	  print*, fact_r
	      if (k_r .gt. 0) then
	        do i_r=k_r,1,-1
c	    if(i_r.le.ntot-k_r)then
c              fact_r=(fact_r*e_val*(1-e_val)*(ntot-i_r+1))/(k_r-i_r+1)
c	    else
	          fact_r=(fact_r*(ntot-i_r+1))/(k_r-i_r+1) ! pay attention to the division by integer!!
c              if (i_end_ph.eq.1) then
c	       print*, e_val, k_r, i_r, fact_r, p_r, s_r
c              endif		
c	    endif
	        enddo
	      endif

c	  p_r=fact_r*e_val**k_r*(1-e_val)**(ntot-k_r)
	      p_r=fact_r
	  
*     Add the probability to the previous overall sum.

	      s_r=s_r+p_r

c              if (i_end_ph.eq.0) then
c	      print*, i_end_ph, k_r_meas, e_val, k_r, fact_r, p_r
c              endif
	      
*     Compare the new sum with 1-alfa

	      if (s_r .ge. 0.683) then
	        end_r=.true.
*     End scan on r variable but continue scan on e_val variable (acceptance interval defined but does contain r_meas)
	      endif
*     Otherwise continue scan on r variable to complete the acceptance interval.
	
	    endif
	
	    m=m+1

200         continue

150	    continue

100	    continue

c      if (plane_ind.eq.1) then
c      write(*,*) 'strip 7', vnhit_strip_sel(1,7),vngood_strip_sel(1,7)
c      write(*,*) eff(7), 'end'
c      endif

            err_bar_sup=e(1)
            err_bar_inf=e(2)

            efficiency=.true.

        endif
	
      endif

c      print*, e(1), e(2)

      end
