*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
*     this file contains all subroutines and functions
*     that are needed for position finding algorithms:
*          
*     subroutine idtoc(ipfa,cpfa)
*
*     subroutine applypfa(PFAtt,ic,ang,corr,res)
*
*     integer function npfastrips(ic,angle)
*
*     real function pfaeta(ic,angle)
*     real function pfaetal(ic,angle)
*     real function pfaeta2(ic,angle) 
*     real function pfaeta3(ic,angle) 
*     real function pfaeta4(ic,angle) 
*     real function cog(ncog,ic)
*
*     real function fbad_cog(ncog,ic)
*     real function fbad_eta(ic,angle)
*
*     real function riseta(iview,angle)
*     FUNCTION risxeta2(x)
*     FUNCTION risxeta3(x)
*     FUNCTION risxeta4(x)
*     FUNCTION risyeta2(x)
*     FUNCTION risy_cog(x)
*     FUNCTION risx_cog(x)
*
*     real function pfacorr(ic,angle) 
*
*     real function effectiveangle(ang,iview,bbb)
*     real function fieldcorr(iview,bbb)
*
*     NB - The angle is the "effective angle", which is relative 
*          to the sensor and it takes into account the magnetic field
*
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***

      subroutine idtoc(ipfa,cpfa)
      
      integer ipfa
      character*10 cpfa

      CPFA='COG4'
      if(ipfa.eq.0)CPFA='ETA'
      if(ipfa.eq.2)CPFA='ETA2'
      if(ipfa.eq.3)CPFA='ETA3'
      if(ipfa.eq.4)CPFA='ETA4'
      if(ipfa.eq.5)CPFA='ETAL'
      if(ipfa.eq.10)CPFA='COG'
      if(ipfa.eq.11)CPFA='COG1'
      if(ipfa.eq.12)CPFA='COG2'
      if(ipfa.eq.13)CPFA='COG3'
      if(ipfa.eq.14)CPFA='COG4'
      
      end
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function effectiveangle(ang,iview,bbb)

      include 'commontracker.f'

      effectiveangle = 0.

      if(mod(iview,2).eq.0)then
c     =================================================
c     X view
c     =================================================
c     here bbb is the y component of the m.field
         angx = ang
         by   = bbb
         if(iview.eq.12) angx = -1. * ang
         if(iview.eq.12) by   = -1. * bbb
cc         tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!!
         tgtemp  = tan(angx*acos(-1.)/180.) + pmuH_h*by*0.00001

      elseif(mod(iview,2).eq.1)then
c     =================================================
c     Y view
c     =================================================         
c     here bbb is the x component of the m.filed
         angy = ang
         bx   = bbb
         tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001         

      endif      
      effectiveangle = 180.*atan(tgtemp)/acos(-1.)

      return
      end
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function fieldcorr(iview,bbb)

      include 'commontracker.f'

      fieldcorr = 0.

      if(mod(iview,2).eq.0)then

c     =================================================
c     X view
c     =================================================
c     here bbb is the y component of the m.field
         by   = bbb
         if(iview.eq.12) by = -1. * bbb
         fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX

      elseif(mod(iview,2).eq.1)then
c     =================================================
c     Y view
c     =================================================         
c     here bbb is the x component of the m.filed
         bx   = bbb
         fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY

      endif      
      
      return
      end
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***

      subroutine applypfa(PFAtt,ic,ang,corr,res)
*---------------------------------------------------------------
*     this subroutine calculate the coordinate of cluster ic (in
*     strip units), relative to the strip with the maximum signal,
*     and its spatial resolution (in cm), applying PFAtt.
*     ang is the effective angle, relative to the sensor
*---------------------------------------------------------------

      character*4 PFAtt
      include 'commontracker.f'
      include 'level1.f'

      corr = 0
      res  = 0

      if(ic.le.0)return

      iview   = VIEW(ic)

      if(mod(iview,2).eq.0)then
c     =================================================
c     X view
c     =================================================

         res = RESXAV 

         if(PFAtt.eq.'COG1')then 

            corr   = 0
            res = 1e-4*pitchX/sqrt(12.)!!res

         elseif(PFAtt.eq.'COG2')then

            corr    = cog(2,ic)            
            res = risx_cog(abs(ang))!TEMPORANEO               
            res = res*fbad_cog(2,ic)

         elseif(PFAtt.eq.'COG3')then

            corr    = cog(3,ic)            
            res = risx_cog(abs(ang))!TEMPORANEO                      
            res = res*fbad_cog(3,ic)

         elseif(PFAtt.eq.'COG4')then

            corr    = cog(4,ic)            
            res = risx_cog(abs(ang))!TEMPORANEO                      
            res = res*fbad_cog(4,ic)

         elseif(PFAtt.eq.'ETA2')then

            corr  = pfaeta2(ic,ang)           
            res = risxeta2(abs(ang))
            res = res*fbad_cog(2,ic)

         elseif(PFAtt.eq.'ETA3')then                        

            corr  = pfaeta3(ic,ang)           
            res = risxeta3(abs(ang))                       
            res = res*fbad_cog(3,ic)               

         elseif(PFAtt.eq.'ETA4')then                         

            corr  = pfaeta4(ic,ang)            
            res = risxeta4(abs(ang))                       
            res = res*fbad_cog(4,ic)               

         elseif(PFAtt.eq.'ETA')then   

            corr  = pfaeta(ic,ang)             
c            res = riseta(ic,ang)                     
            res = riseta(iview,ang)                     
            res = res*fbad_eta(ic,ang)             

         elseif(PFAtt.eq.'ETAL')then   

            corr  = pfaetal(ic,ang)             
            res = riseta(iview,ang)                     
            res = res*fbad_eta(ic,ang)             

         elseif(PFAtt.eq.'COG')then          

            corr  = cog(0,ic)             
            res = risx_cog(abs(ang))                     
            res = res*fbad_cog(0,ic)

         else
            if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
         endif


*     ======================================
*     temporary patch for saturated clusters
*     ======================================
         if( nsatstrips(ic).gt.0 )then
            corr  = cog(4,ic)             
            res = pitchX*1e-4/sqrt(12.)
cc            cc=cog(4,ic)
c$$$            print*,ic,' *** ',cc
c$$$            print*,ic,' *** ',res
         endif


      elseif(mod(iview,2).eq.1)then
c     =================================================
c     Y view
c     =================================================

         res = RESYAV 

         if(PFAtt.eq.'COG1')then

            corr  = 0  
            res = 1e-4*pitchY/sqrt(12.)!res  

         elseif(PFAtt.eq.'COG2')then

            corr    = cog(2,ic)
            res = risy_cog(abs(ang))!TEMPORANEO 
            res = res*fbad_cog(2,ic)

         elseif(PFAtt.eq.'COG3')then

            corr    = cog(3,ic)
            res = risy_cog(abs(ang))!TEMPORANEO 
            res = res*fbad_cog(3,ic)

         elseif(PFAtt.eq.'COG4')then

            corr    = cog(4,ic)
            res = risy_cog(abs(ang))!TEMPORANEO 
            res = res*fbad_cog(4,ic)

         elseif(PFAtt.eq.'ETA2')then

            corr    = pfaeta2(ic,ang)          
            res = risyeta2(abs(ang))              
            res = res*fbad_cog(2,ic)

         elseif(PFAtt.eq.'ETA3')then                      

            corr    = pfaeta3(ic,ang) 
            res = res*fbad_cog(3,ic)  

         elseif(PFAtt.eq.'ETA4')then  

            corr    = pfaeta4(ic,ang) 
            res = res*fbad_cog(4,ic)

         elseif(PFAtt.eq.'ETA')then 

            corr    = pfaeta(ic,ang)
c            res = riseta(ic,ang)  
            res = riseta(iview,ang)  
            res = res*fbad_eta(ic,ang)

         elseif(PFAtt.eq.'ETAL')then 

            corr    = pfaetal(ic,ang)
            res = riseta(iview,ang)  
            res = res*fbad_eta(ic,ang)

         elseif(PFAtt.eq.'COG')then

            corr    = cog(0,ic)            
            res = risy_cog(abs(ang))
            res = res*fbad_cog(0,ic)

         else
            if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
         endif


*     ======================================
*     temporary patch for saturated clusters
*     ======================================
         if( nsatstrips(ic).gt.0 )then
            corr    = cog(4,ic)             
            res = pitchY*1e-4/sqrt(12.)
cc            cc=cog(4,ic)
c$$$            print*,ic,' *** ',cc
c$$$            print*,ic,' *** ',res
         endif
        
      endif
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      integer function npfastrips(ic,angle)
*--------------------------------------------------------------
*     thid function returns the number of strips used
*     to evaluate the position of a cluster, according to the p.f.a.
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'level1.f'
      include 'calib.f'

      character*4 usedPFA
      


      call idtoc(pfaid,usedPFA)

      npfastrips=-1

      if(usedPFA.eq.'COG1')npfastrips=1
      if(usedPFA.eq.'COG2')npfastrips=2
      if(usedPFA.eq.'COG3')npfastrips=3
      if(usedPFA.eq.'COG4')npfastrips=4
      if(usedPFA.eq.'ETA2')npfastrips=2
      if(usedPFA.eq.'ETA3')npfastrips=3
      if(usedPFA.eq.'ETA4')npfastrips=4
*     ----------------------------------------------------------------
      if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
c         print*,VIEW(ic),angle
         if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
            if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
               npfastrips=2
            elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
               npfastrips=3
            elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
               npfastrips=4
            else
               npfastrips=4     !COG4
            endif                        
         else                   !X-view
            if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
               npfastrips=2
            elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
               npfastrips=3
            elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
               npfastrips=4
            else
               npfastrips=4     !COG4
            endif                        
         endif
      endif
*     ----------------------------------------------------------------
      if(usedPFA.eq.'COG')then

         npfastrips=0

c$$$         iv=VIEW(ic)
c$$$         if(mod(iv,2).eq.1)incut=incuty
c$$$         if(mod(iv,2).eq.0)incut=incutx
c$$$         istart = INDSTART(IC)
c$$$         istop  = TOTCLLENGTH
c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
c$$$         mu  = 0
c$$$         do i = INDMAX(IC),istart,-1
c$$$            ipos = i-INDMAX(ic)
c$$$            cut  = incut*CLSIGMA(i)
c$$$            if(CLSIGNAL(i).ge.cut)then
c$$$               mu = mu + 1
c$$$               print*,i,mu
c$$$            else
c$$$               goto 10
c$$$            endif
c$$$         enddo
c$$$ 10      continue
c$$$         do i = INDMAX(IC)+1,istop
c$$$            ipos = i-INDMAX(ic)
c$$$            cut  = incut*CLSIGMA(i)
c$$$            if(CLSIGNAL(i).ge.cut)then
c$$$               mu = mu + 1
c$$$               print*,i,mu
c$$$            else
c$$$               goto 20
c$$$            endif
c$$$         enddo
c$$$ 20      continue
c$$$         npfastrips=mu

      endif
*     ----------------------------------------------------------------

c      print*,pfaid,usedPFA,angle,npfastrips

      return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function pfaeta(ic,angle)
*--------------------------------------------------------------
*     this function returns the position (in strip units)
*     it calls: 
*     - pfaeta2(ic,angle)
*     - pfaeta3(ic,angle)
*     - pfaeta4(ic,angle)
*     according to the angle
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'level1.f'
      include 'calib.f'
      
      pfaeta = 0

      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
      
         if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
            pfaeta = pfaeta2(ic,angle)
cc            print*,pfaeta2(ic,angle)
         elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
            pfaeta = pfaeta3(ic,angle)
         elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
            pfaeta = pfaeta4(ic,angle)
         else
            pfaeta = cog(4,ic)
         endif            

      else                      !X-view

         if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
            pfaeta = pfaeta2(ic,angle)
         elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
            pfaeta = pfaeta3(ic,angle)
         elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
            pfaeta = pfaeta4(ic,angle)
         else
            pfaeta = cog(4,ic)
         endif            
            
      endif
      
 100  return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function pfaetal(ic,angle)
*--------------------------------------------------------------
*     this function returns the position (in strip units)
*     it calls: 
*     - pfaeta2(ic,angle)+pfcorr(ic,angle)
*     - pfaeta3(ic,angle)+pfcorr(ic,angle)
*     - pfaeta4(ic,angle)+pfcorr(ic,angle)
*     according to the angle
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'level1.f'
      include 'calib.f'
      
      pfaetal = 0

      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
      
         if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
            pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
         elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
            pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
         elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
            pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
         else
            pfaetal = cog(4,ic)
         endif            

      else                      !X-view

         if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
            pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
         elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
            pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
         elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
            pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
         else
            pfaetal = cog(4,ic)
         endif            
            
      endif
      
 100  return
      end
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c      real function riseta(ic,angle)
      real function riseta(iview,angle)
*--------------------------------------------------------------
*     this function returns the average spatial resolution 
*     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
*     it calls: 
*     - risxeta2(angle)
*     - risyeta2(angle)
*     - risxeta3(angle)
*     - risxeta4(angle)
*     according to the angle
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'level1.f'
      include 'calib.f'

      riseta = 0

c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
      if(mod(iview,2).eq.1)then !Y-view
      

         if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
            riseta = risyeta2(angle)
         elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
            riseta = risy_cog(angle) !ATTENZIONE!!
         elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
            riseta = risy_cog(angle) !ATTENZIONE!!
         else
            riseta = risy_cog(angle)
         endif            

      else                      !X-view

         if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
            riseta = risxeta2(angle)
         elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
            riseta = risxeta3(angle) 
         elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
            riseta = risxeta4(angle)
         else
            riseta = risx_cog(angle)
         endif            
            
      endif


 100  return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function fbad_eta(ic,angle)
*-------------------------------------------------------
*     this function returns a factor that takes into 
*     account deterioration of the spatial resolution
*     in the case BAD strips are included in the cluster.
*     This factor should multiply the nominal spatial 
*     resolution.
*     It calls the function FBAD_COG(NCOG,IC), 
*     accordingto the angle
*-------------------------------------------------------

      include 'commontracker.f'
      include 'level1.f'
      include 'calib.f'
      fbad_eta = 0

      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
      
         if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
            fbad_eta = fbad_cog(2,ic)
         elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
            fbad_eta = fbad_cog(3,ic)
         elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
            fbad_eta = fbad_cog(4,ic)
         else
            fbad_eta = fbad_cog(4,ic)
         endif            

      else                      !X-view

         if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
            fbad_eta = fbad_cog(2,ic)
         elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
            fbad_eta = fbad_cog(3,ic)
         elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
            fbad_eta = fbad_cog(4,ic)
         else
            fbad_eta = fbad_cog(4,ic)
         endif            
            
      endif

      return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function pfaeta2(ic,angle) !(1)
*--------------------------------------------------------------
*     this function returns 
*
*     - the position (in strip units) 
*     corrected according to the ETA2 Position Finding Algorithm.
*     The function performs an interpolation of FETA2%ETA2
*     
*     - if the angle is out of range, the calibration parameters
*     of the lowest or higher bin are used
*     
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'calib.f'
      include 'level1.f'

      real cog2,angle
      integer iview,lad

      iview = VIEW(ic)            
      lad = nld(MAXS(ic),VIEW(ic))
      cog2 = cog(2,ic)           
      pfaeta2=cog2

*     ----------------
*     find angular bin
*     ----------------
*     (in futuro possiamo pensare di interpolare anche sull'angolo)
      do iang=1,nangbin
         if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
            iangle=iang
            goto 98
         endif
      enddo
      if(DEBUG.EQ.1)
     $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
      if(angle.le.angL(1))iang=1
      if(angle.ge.angR(nangbin))iang=nangbin
 98   continue                  !jump here if ok

*     -------------
*     within +/-0.5
*     -------------

      iaddmax=10
      iadd=0
 10   continue
      if(cog2.lt.eta2(1,iang))then
         cog2 = cog2 + 1
         iadd = iadd + 1
         if(iadd>iaddmax)goto 111
         goto 10
      endif
 20   continue
      if(cog2.gt.eta2(netaval,iang))then
         cog2 = cog2 - 1
         iadd = iadd - 1
         if(iadd<-1*iaddmax)goto 111
         goto 20
      endif
      goto 1111
 111  continue
      if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
      if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
      cog2=0
 1111 continue

*     --------------------------------
c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
      do i=2,netaval
         if(eta2(i,iang).gt.cog2)then

            x1 = eta2(i-1,iang)
            x2 = eta2(i,iang)
            y1 = feta2(i-1,iview,lad,iang)
            y2 = feta2(i,iview,lad,iang)

c            print*,'*****',i,view,lad,iang
c            print*,'-----',x1,x2,y1,y2
            goto 99
         endif
      enddo
 99   continue


      AA=(y2-y1)/(x2-x1)
      BB=y1-AA*x1

      pfaeta2 = AA*cog2+BB
      pfaeta2 = pfaeta2 - iadd

c$$$      if(iflag.eq.1)then
c$$$         pfaeta2=pfaeta2-1.   !temp
c$$$         cog2=cog2-1.           !temp
c$$$      endif
c$$$      if(iflag.eq.-1)then
c$$$         pfaeta2=pfaeta2+1.   !temp
c$$$         cog2=cog2+1.           !temp
c$$$      endif

      if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
     $     ,cog2-iadd,' -->',pfaeta2


 100  return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function pfaeta3(ic,angle) !(1)
*--------------------------------------------------------------
*     this function returns 
*
*     - the position (in strip units) 
*     corrected according to the ETA3 Position Finding Algorithm.
*     The function performs an interpolation of FETA3%ETA3
*     
*     - if the angle is out of range, the calibration parameters
*     of the lowest or higher bin are used
*     
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'calib.f'
      include 'level1.f'

      real cog3,angle
      integer iview,lad


      iview = VIEW(ic)            
      lad = nld(MAXS(ic),VIEW(ic))
      cog3 = cog(3,ic)
      cc = cog3
      cog3 = cc
      pfaeta3=cog3

*     ----------------
*     find angular bin
*     ----------------
*     (in futuro possiamo pensare di interpolare anche sull'angolo)
      do iang=1,nangbin
c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
         if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
            iangle=iang
            goto 98
         endif
      enddo
      if(DEBUG.EQ.1)
     $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
      if(angle.le.angL(1))iang=1
      if(angle.ge.angR(nangbin))iang=nangbin
 98   continue                  !jump here if ok

*     -------------
*     within +/-0.5
*     -------------

      iaddmax=10
      iadd=0
 10   continue
      if(cog3.lt.eta3(1,iang))then
         cog3   = cog3 + 1.
         iadd = iadd + 1
         if(iadd>iaddmax) goto 111
         goto 10
      endif
 20   continue
      if(cog3.gt.eta3(netaval,iang))then
         cog3 = cog3 - 1.
         iadd = iadd - 1
         if(iadd<-1*iaddmax) goto 111
         goto 20
      endif
      goto 1111
 111  continue
      if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
      if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
      cog3=0      
 1111 continue

*     --------------------------------
c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
      do i=2,netaval
         if(eta3(i,iang).gt.cog3)then

            x1 = eta3(i-1,iang)
            x2 = eta3(i,iang)
            y1 = feta3(i-1,iview,lad,iang)
            y2 = feta3(i,iview,lad,iang)

c            print*,'*****',i,view,lad,iang
c            print*,'-----',x1,x2,y1,y2
            goto 99
         endif
      enddo
 99   continue


      AA=(y2-y1)/(x2-x1)
      BB=y1-AA*x1

      pfaeta3 = AA*cog3+BB
      pfaeta3 = pfaeta3 - iadd


      if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
     $     ,cog3-iadd,' -->',pfaeta3

 100  return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function pfaeta4(ic,angle) 
*--------------------------------------------------------------
*     this function returns 
*
*     - the position (in strip units) 
*     corrected according to the ETA4 Position Finding Algorithm.
*     The function performs an interpolation of FETA3%ETA3
*     
*     - if the angle is out of range, the calibration parameters
*     of the lowest or higher bin are used
*     
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'calib.f'
      include 'level1.f'

      real cog4,angle
      integer iview,lad


      iview = VIEW(ic)            
      lad = nld(MAXS(ic),VIEW(ic))
      cog4=cog(4,ic)
      pfaeta4=cog4

*     ----------------
*     find angular bin
*     ----------------
*     (in futuro possiamo pensare di interpolare anche sull'angolo)
      do iang=1,nangbin
c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
         if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
            iangle=iang
            goto 98
         endif
      enddo
      if(DEBUG.EQ.1)
     $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
      if(angle.le.angL(1))iang=1
      if(angle.ge.angR(nangbin))iang=nangbin
 98   continue                  !jump here if ok

*     -------------
*     within +/-0.5
*     -------------

      iaddmax=10
      iadd=0
 10   continue
      if(cog4.lt.eta4(1,iang))then
         cog4 = cog4 + 1
         iadd = iadd + 1
         if(iadd>iaddmax)goto 111
         goto 10
      endif
 20   continue
      if(cog4.gt.eta4(netaval,iang))then
         cog4 = cog4 - 1
         iadd = iadd - 1
         if(iadd<-1*iaddmax)goto 111
         goto 20
      endif
      goto 1111
 111  continue
      if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
      if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
      cog4=0
 1111 continue

*     --------------------------------
c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
      do i=2,netaval
         if(eta4(i,iang).gt.cog4)then

            x1 = eta4(i-1,iang)
            x2 = eta4(i,iang)
            y1 = feta4(i-1,iview,lad,iang)
            y2 = feta4(i,iview,lad,iang)

c            print*,'*****',i,view,lad,iang
c            print*,'-----',x1,x2,y1,y2
            goto 99
         endif
      enddo
 99   continue


      AA=(y2-y1)/(x2-x1)
      BB=y1-AA*x1

      pfaeta4 = AA*cog4+BB
      pfaeta4 = pfaeta4 - iadd

c$$$      if(iflag.eq.1)then
c$$$         pfaeta2=pfaeta2-1.   !temp
c$$$         cog2=cog2-1.           !temp
c$$$      endif
c$$$      if(iflag.eq.-1)then
c$$$         pfaeta2=pfaeta2+1.   !temp
c$$$         cog2=cog2+1.           !temp
c$$$      endif

      if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
     $     ,cog4-iadd,' -->',pfaeta4

 100  return
      end



*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function cog(ncog,ic)
*-------------------------------------------------
*     this function returns 
*
*     - if NCOG=0, the Center-Of-Gravity of the 
*     cluster IC, relative to MAXS(IC), according to
*     the cluster multiplicity
*     
*     - if NCOG>0, the Center-Of-Gravity of the cluster IC 
*     evaluated using NCOG strips, even if they have a 
*     negative signal (according to Landi)
*          
*-------------------------------------------------


      include 'commontracker.f'
      include 'calib.f'
      include 'level1.f'
      


      if (ncog.gt.0) then
*     ===========================
*     ETA2 ETA3 ETA4 computation
*     ===========================

*     --> signal of the central strip
         sc = CLSIGNAL(INDMAX(ic)) !center
*     signal of adjacent strips
         sl1 = -9999.           !left 1
         if(
     $        (INDMAX(ic)-1).ge.INDSTART(ic)
     $        )
     $        sl1 = CLSIGNAL(INDMAX(ic)-1)
         
         sl2 = -9999.           !left 2
         if(
     $        (INDMAX(ic)-2).ge.INDSTART(ic)
     $        )
     $        sl2 = CLSIGNAL(INDMAX(ic)-2)
         
         sr1 = -9999.           !right 1
         if(
     $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
     $        .or.
     $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
     $        )
     $        sr1 = CLSIGNAL(INDMAX(ic)+1)
         
         sr2 = -9999.           !right 2
         if(
     $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
     $        .or.
     $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
     $        )
     $        sr2 = CLSIGNAL(INDMAX(ic)+2)
         
         COG = 0.
         
c      print *,'## ',sl2,sl1,sc,sr1,sr2

c     ==============================================================
         if(ncog.eq.1)then
            COG = 0.
	    if(sr1.gt.sc)cog=1.
	    if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
c     ==============================================================
         elseif(ncog.eq.2)then
            COG = 0.
            if(sl1.gt.sr1)then
               if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
            elseif(sl1.lt.sr1)then
               if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
	    elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
               if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
     $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) 
               if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
     $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) 
	    endif 
c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
c     $           ,' : ',sl2,sl1,sc,sr1,sr2
c     ==============================================================
         elseif(ncog.eq.3)then
            COG = 0
            sss = sc
            if( sl1.ne.-9999. )COG = COG-sl1
            if( sl1.ne.-9999. )sss = sss+sl1
            if( sr1.ne.-9999. )COG = COG+sr1
            if( sr1.ne.-9999. )sss = sss+sr1
            if(sss.ne.0)COG=COG/sss

c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1) 
c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
c     $            ,' : ',sl2,sl1,sc,sr1,sr2
c     ==============================================================
         elseif(ncog.eq.4)then

            COG = 0
            sss = sc
            if( sl1.ne.-9999. )COG = COG-sl1
            if( sl1.ne.-9999. )sss = sss+sl1
            if( sr1.ne.-9999. )COG = COG+sr1
            if( sr1.ne.-9999. )sss = sss+sr1
            if(sl2.gt.sr2)then
               if((sl2+sss).ne.0) 
     $              COG = (COG-2*sl2)/(sl2+sss) 
            elseif(sl2.lt.sr2)then
               if((sr2+sss).ne.0)
     $              COG = (2*sr2+COG)/(sr2+sss) 
            elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then 
               if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
     $              .and.(sl2+sss).ne.0 )
     $              cog = (cog-2*sl2)/(sl2+sss) 
               if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
     $              .and.(sr2+sss).ne.0 )
     $              cog = (2*sr2+cog)/(sr2+sss)               
            endif
c     ==============================================================
         elseif(ncog.eq.5)then
            COG = 0
            sss = sc
            if( sl1.ne.-9999. )COG = COG-sl1
            if( sl1.ne.-9999. )sss = sss+sl1
            if( sr1.ne.-9999. )COG = COG+sr1
            if( sr1.ne.-9999. )sss = sss+sr1
            if( sl2.ne.-9999. )COG = COG-2*sl2
            if( sl2.ne.-9999. )sss = sss+sl2
            if( sr2.ne.-9999. )COG = COG+2*sr2
            if( sr2.ne.-9999. )sss = sss+sr2
            if(sss.ne.0)COG=COG/sss
         else
            print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
     $           ,' not implemented'
            COG = 0.
         endif

c         print*,'NCOG ',ncog,ic,' @@@ ',sl1,sc,sr1,' @@@ ',cog
            
      elseif(ncog.eq.0)then
*     =========================
*     COG computation
*     =========================

         iv=VIEW(ic)
         if(mod(iv,2).eq.1)incut=incuty
         if(mod(iv,2).eq.0)incut=incutx
         istart = INDSTART(IC)
         istop  = TOTCLLENGTH
         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
         COG = 0  
         SGN = 0.
         mu  = 0
c         print*,'-------'
         do i = INDMAX(IC),istart,-1
            ipos = i-INDMAX(ic)
            cut  = incut*CLSIGMA(i)
            if(CLSIGNAL(i).ge.cut)then               
               COG = COG + ipos*CLSIGNAL(i)
               SGN = SGN + CLSIGNAL(i)
               mu = mu + 1
c               print*,ipos,CLSIGNAL(i)
            else
               goto 10
            endif
         enddo
 10      continue
         do i = INDMAX(IC)+1,istop
            ipos = i-INDMAX(ic)
            cut  = incut*CLSIGMA(i)
            if(CLSIGNAL(i).ge.cut)then
               COG = COG + ipos*CLSIGNAL(i)
               SGN = SGN + CLSIGNAL(i)
               mu = mu + 1
c               print*,ipos,CLSIGNAL(i)
            else
               goto 20
            endif
         enddo
 20      continue
         if(SGN.le.0)then
            print*,'cog(0,ic) --> ic, dedx ',ic,SGN
            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
            print*,(CLSIGNAL(i),i=istart,istop)
c            print*,'cog(0,ic) --> NOT EVALUATED '
         else
            COG=COG/SGN
         endif
c         print*,'-------'
         
      else
         
         COG=0
         print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
         print*,'                 (NCOG must be >= 0)'
         

      endif

c      print *,'## cog ',ncog,ic,cog,'/////////////'

      if(COG.lt.-0.75.or.COG.gt.+0.75)then
         if(DEBUG.eq.1)
     $        print*,'cog *** warning *** anomalous cluster ??? --> '
         if(DEBUG.eq.1)
     $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG 
      endif


      return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***

      real function fbad_cog(ncog,ic)
*-------------------------------------------------------
*     this function returns a factor that takes into 
*     account deterioration of the spatial resolution
*     in the case BAD strips are included in the cluster.
*     This factor should multiply the nominal spatial 
*     resolution.
*
*-------------------------------------------------------

      include 'commontracker.f'
      include 'level1.f'
      include 'calib.f'

      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
         si = 8.4  !average good-strip noise
         f  = 4.   !average bad-strip noise: f*si
         incut=incuty
      else                      !X-view
         si = 3.9  !average good-strip noise
         f  = 6.   !average bad-strip noise: f*si
         incut=incutx
      endif
      
      fbad_cog = 1.

      if (ncog.gt.0) then

*     --> signal of the central strip
         sc = CLSIGNAL(INDMAX(ic)) !center
         fsc = 1
c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
         fsc = clsigma(INDMAX(ic))/si
*     --> signal of adjacent strips
         sl1 = 0                !left 1
         fsl1 = 1               !left 1
         if(
     $        (INDMAX(ic)-1).ge.INDSTART(ic)
     $        )then
            sl1 = CLSIGNAL(INDMAX(ic)-1)
c            if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f
            fsl1 = clsigma(INDMAX(ic)-1)/si
         endif

         sl2 = 0                !left 2
         fsl2 = 1               !left 2
         if(
     $        (INDMAX(ic)-2).ge.INDSTART(ic)
     $        )then
            sl2 = CLSIGNAL(INDMAX(ic)-2)
c            if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f
            fsl2 = clsigma(INDMAX(ic)-2)/si
         endif
         sr1 = 0                !right 1
         fsr1 = 1               !right 1
         if(
     $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
     $        .or.
     $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
     $        )then
            sr1 = CLSIGNAL(INDMAX(ic)+1)
c            if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f
            fsr1 = clsigma(INDMAX(ic)+1)/si
         endif    
         sr2 = 0                !right 2
         fsr2 = 1               !right 2
         if(
     $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
     $        .or.
     $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
     $        )then
            sr2 = CLSIGNAL(INDMAX(ic)+2)
c            if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f
            fsr2 = clsigma(INDMAX(ic)+2)/si
         endif



************************************************************
*     COG2-3-4 computation
************************************************************

c      print*,sl2,sl1,sc,sr1,sr2
         
         vCOG = cog(ncog,ic)!0.
         
         if(ncog.eq.2)then
            if(sl1.gt.sr1)then
c               COG = -sl1/(sl1+sc)        
               fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
               fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
            elseif(sl1.le.sr1)then
c               COG = sr1/(sc+sr1)            
               fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
               fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
            endif
         elseif(ncog.eq.3)then
c            COG = (sr1-sl1)/(sl1+sc+sr1) 
            fbad_cog = 
     $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
            fbad_cog = 
     $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
         elseif(ncog.eq.4)then
            if(sl2.gt.sr2)then
c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) 
               fbad_cog = 
     $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
     $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
               fbad_cog = 
     $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
     $              +(-vCOG)**2+(1-vCOG)**2)
            elseif(sl2.le.sr2)then
c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) 
               fbad_cog = 
     $              (fsl1*(-1-vCOG)**2
     $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
               fbad_cog = 
     $              fbad_cog / ((-1-vCOG)**2
     $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
            endif
         else
            print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
            print*,'                               (NCOG must be <= 4)'
c            COG = 0.
         endif
         
      elseif(ncog.eq.0)then
*     =========================
*     COG computation
*     =========================
         
         vCOG = cog(0,ic)

         iv     = VIEW(ic)
         istart = INDSTART(IC)
         istop  = TOTCLLENGTH
         if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
         SGN = 0.
         SNU = 0.
         SDE = 0.

         do i=INDMAX(IC),istart,-1
            ipos = i-INDMAX(ic)
            cut  = incut*CLSIGMA(i)
            if(CLSIGNAL(i).gt.cut)then
               fs = clsigma(i)/si
               SNU = SNU + fs*(ipos-vCOG)**2
               SDE = SDE + (ipos-vCOG)**2
            else
               goto 10
            endif            
         enddo
 10      continue
         do i=INDMAX(IC)+1,istop
            ipos = i-INDMAX(ic)
            cut  = incut*CLSIGMA(i)
            if(CLSIGNAL(i).gt.cut)then
               fs = clsigma(i)/si
               SNU = SNU + fs*(ipos-vCOG)**2
               SDE = SDE + (ipos-vCOG)**2
            else
               goto 20
            endif            
         enddo
 20      continue
         if(SDE.ne.0)then
            FBAD_COG=SNU/SDE
         else
            
         endif

      else
                  
         FBAD_COG=0
         print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
         print*,'                               (NCOG must be >= 0)'
         

      endif


      fbad_cog = sqrt(fbad_cog)

      return
      end


c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c$$$      real function fbad_cog0(ncog,ic)
c$$$*-------------------------------------------------------
c$$$*     this function returns a factor that takes into 
c$$$*     account deterioration of the spatial resolution
c$$$*     in the case BAD strips are included in the cluster.
c$$$*     This factor should multiply the nominal spatial 
c$$$*     resolution.
c$$$*
c$$$*     NB!!!
c$$$*     (this is the old version. It consider only the two 
c$$$*     strips with the greatest signal. The new one is 
c$$$*     fbad_cog(ncog,ic) )
c$$$*     
c$$$*-------------------------------------------------------
c$$$
c$$$      include 'commontracker.f'
c$$$      include 'level1.f'
c$$$      include 'calib.f'
c$$$
c$$$*     --> signal of the central strip
c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
c$$$
c$$$*     signal of adjacent strips
c$$$*     --> left
c$$$      sl1 = 0                  !left 1
c$$$      if(
c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
c$$$     $     )
c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
c$$$
c$$$      sl2 = 0                  !left 2
c$$$      if(
c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
c$$$     $     )
c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
c$$$
c$$$*     --> right
c$$$      sr1 = 0                  !right 1
c$$$      if(
c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
c$$$     $     .or.
c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
c$$$     $     )
c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
c$$$
c$$$      sr2 = 0                  !right 2
c$$$      if(
c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
c$$$     $     .or.
c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
c$$$     $     )
c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
c$$$
c$$$
c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
c$$$         f  = 4.
c$$$         si = 8.4
c$$$      else                              !X-view
c$$$         f  = 6.
c$$$         si = 3.9
c$$$      endif
c$$$
c$$$      fbad_cog = 1.
c$$$      f0 = 1
c$$$      f1 = 1
c$$$      f2 = 1
c$$$      f3 = 1   
c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
c$$$         
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
c$$$
c$$$         if(ncog.eq.2.and.sl1.ne.0)then
c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
c$$$            fbad_cog = 1.
c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
c$$$            fbad_cog = 1.
c$$$         else
c$$$            fbad_cog = 1.
c$$$         endif
c$$$         
c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
c$$$
c$$$
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
c$$$
c$$$         if(ncog.eq.2.and.sr1.ne.0)then
c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
c$$$            fbad_cog = 1.
c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
c$$$            fbad_cog = 1.
c$$$         else
c$$$            fbad_cog = 1.
c$$$         endif
c$$$
c$$$      endif
c$$$
c$$$      fbad_cog0 = sqrt(fbad_cog)
c$$$
c$$$      return
c$$$      end
c$$$
c$$$
c$$$

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***

      FUNCTION risxeta2(x)

      DOUBLE PRECISION V( 1)
      INTEGER NPAR, NDIM, IMQFUN, I, J
      DOUBLE PRECISION HQDJ, VV, VCONST
      DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
      DOUBLE PRECISION SIGV(  18, 1)
      DOUBLE PRECISION SIGDEL(  18)
      DOUBLE PRECISION SIGA(  18)
      DATA NPAR, NDIM, IMQFUN /   18,    1,    1/
      DATA VCONST /  0.000000000000    /
      DATA SIGVMI / -20.50000000000    /
      DATA SIGVT /  41.00000000000    /
      DATA SIGV / 0.6097560748458E-01
     +, 0.1097560971975    
     +, 0.1341463327408    
     +, 0.1829268187284    
     +, 0.2317073047161    
     +, 0.4268292486668    
     +, 0.4756097495556    
     +, 0.4999999701977    
     +, 0.5243902206421    
     +, 0.5731707215309    
     +, 0.7682926654816    
     +, 0.8170731663704    
     +, 0.8658536076546    
     +, 0.8902438879013    
     +, 0.9390243291855    
     +,  0.000000000000    
     +,  1.000000000000    
     +, 0.3658536374569    
     +/
      DATA SIGDEL / 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.1999999994950E-05
     +, 0.1999999994950E-05
     +, 0.9756097197533E-01
     +/
      DATA SIGA /  51.65899502118    
     +, -150.4733247841    
     +,  143.0468613786    
     +, -16.56096738997    
     +,  5.149319798083    
     +,  21.57149712673    
     +, -39.46652322782    
     +,  47.13181632948    
     +, -32.93197883680    
     +,  16.38645317092    
     +,  1.453688482992    
     +, -10.00547244421    
     +,  131.3517670587    
     +, -140.6351538257    
     +,  49.05515749582    
     +, -23.00028974788    
     +, -22.58470403729    
     +, -3.824682486418    
     +/

      V(1)= abs(x)
      if(V(1).gt.20.)V(1)=20.

      HQUADF = 0.
      DO 20 J = 1, NPAR
         HQDJ = 0.
         DO 10 I = 1, NDIM
            VV = (V (I) - SIGVMI (I)) / SIGVT (I)
            HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
   10    CONTINUE
         HQDJ = HQDJ + SIGDEL (J) ** 2
         HQDJ = SQRT (HQDJ)
         HQUADF = HQUADF + SIGA (J) * HQDJ
   20 CONTINUE
      IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
      
      risxeta2=HQUADF* 1e-4

      END

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      FUNCTION risxeta3(x)
      DOUBLE PRECISION V( 1)
      INTEGER NPAR, NDIM, IMQFUN, I, J
      DOUBLE PRECISION HQDJ, VV, VCONST
      DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
      DOUBLE PRECISION SIGV(  18, 1)
      DOUBLE PRECISION SIGDEL(  18)
      DOUBLE PRECISION SIGA(  18)
      DATA NPAR, NDIM, IMQFUN /   18,    1,    1/
      DATA VCONST /  0.000000000000    /
      DATA SIGVMI / -20.50000000000    /
      DATA SIGVT /  41.00000000000    /
      DATA SIGV / 0.6097560748458E-01
     +, 0.1097560971975    
     +, 0.1341463327408    
     +, 0.1829268187284    
     +, 0.2317073047161    
     +, 0.4756097495556    
     +, 0.4999999701977    
     +, 0.5243902206421    
     +, 0.7682926654816    
     +, 0.8170731663704    
     +, 0.8658536076546    
     +, 0.8902438879013    
     +, 0.9390243291855    
     +,  0.000000000000    
     +,  1.000000000000    
     +, 0.3658536374569    
     +, 0.4146341383457    
     +, 0.6097560524940    
     +/
      DATA SIGDEL / 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.1999999994950E-05
     +, 0.1999999994950E-05
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +/
      DATA SIGA /  55.18284054458    
     +, -160.3358431242    
     +,  144.6939185763    
     +, -20.45200854118    
     +,  5.223570087108    
     +,-0.4171476953945    
     +, -27.67911907462    
     +,  17.70327157495    
     +, -1.867165491707    
     +, -8.884458169181    
     +,  124.3526608791    
     +, -143.3309398345    
     +,  50.80345027122    
     +, -16.44454904415    
     +, -15.73785568450    
     +, -22.71810502561    
     +,  36.86170101430    
     +,  2.437918198452    
     +/

      V(1) =  abs(x)
      if(V(1).gt.20.)V(1)=20.

      HQUADF = 0.
      DO 20 J = 1, NPAR
         HQDJ = 0.
         DO 10 I = 1, NDIM
            VV = (V (I) - SIGVMI (I)) / SIGVT (I)
            HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
   10    CONTINUE
         HQDJ = HQDJ + SIGDEL (J) ** 2
         HQDJ = SQRT (HQDJ)
         HQUADF = HQUADF + SIGA (J) * HQDJ
   20 CONTINUE
      IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
      
      risxeta3 = HQUADF* 1e-4

      END
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      FUNCTION risxeta4(x)
      DOUBLE PRECISION V( 1)
      INTEGER NPAR, NDIM, IMQFUN, I, J
      DOUBLE PRECISION HQDJ, VV, VCONST
      DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
      DOUBLE PRECISION SIGV(  18, 1)
      DOUBLE PRECISION SIGDEL(  18)
      DOUBLE PRECISION SIGA(  18)
      DATA NPAR, NDIM, IMQFUN /   18,    1,    1/
      DATA VCONST /  0.000000000000    /
      DATA SIGVMI / -20.50000000000    /
      DATA SIGVT /  41.00000000000    /
      DATA SIGV / 0.3658536449075E-01
     +, 0.6097560748458E-01
     +, 0.1097560971975    
     +, 0.1341463327408    
     +, 0.4756097495556    
     +, 0.5243902206421    
     +, 0.8658536076546    
     +, 0.8902438879013    
     +, 0.9390243291855    
     +, 0.9634146094322    
     +,  0.000000000000    
     +,  1.000000000000    
     +, 0.3658536374569    
     +, 0.4146341383457    
     +, 0.6097560524940    
     +, 0.6585365533829    
     +, 0.7560975551605    
     +, 0.2439024299383    
     +/
      DATA SIGDEL / 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.1999999994950E-05
     +, 0.1999999994950E-05
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.1951219439507    
     +/
      DATA SIGA / -43.61551887895    
     +,  57.88466995373    
     +, -92.04113299504    
     +,  74.08166649890    
     +, -9.768686062558    
     +, -4.304496875334    
     +,  72.62237333937    
     +, -91.21920840618    
     +,  56.75519978630    
     +, -43.21115751243    
     +,  12.79984505413    
     +,  12.10074868595    
     +, -6.238587250860    
     +,  23.43447356326    
     +,  17.98221401495    
     +, -7.980332610975    
     +, -3.426733307051    
     +, -8.683439558751    
     +/

      V(1)=abs(x)
      if(V(1).gt.20.)V(1)=20.

      HQUADF = 0.
      DO 20 J = 1, NPAR
         HQDJ = 0.
         DO 10 I = 1, NDIM
            VV = (V (I) - SIGVMI (I)) / SIGVT (I)
            HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
   10    CONTINUE
         HQDJ = HQDJ + SIGDEL (J) ** 2
         HQDJ = SQRT (HQDJ)
         HQUADF = HQUADF + SIGA (J) * HQDJ
   20 CONTINUE
      IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
      
      risxeta4=HQUADF* 1e-4

      END
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      FUNCTION risyeta2(x)
      DOUBLE PRECISION V( 1)
      INTEGER NPAR, NDIM, IMQFUN, I, J
      DOUBLE PRECISION HQDJ, VV, VCONST
      DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
      DOUBLE PRECISION SIGV(  12, 1)
      DOUBLE PRECISION SIGDEL(  12)
      DOUBLE PRECISION SIGA(  12)
      DATA NPAR, NDIM, IMQFUN /   12,    1,    1/
      DATA VCONST /  0.000000000000    /
      DATA SIGVMI / -20.50000000000    /
      DATA SIGVT /  41.00000000000    /
      DATA SIGV / 0.1585365831852    
     +, 0.4024389982224    
     +, 0.4756097495556    
     +, 0.5243902206421    
     +, 0.5975609421730    
     +, 0.8414633870125    
     +,  0.000000000000    
     +,  1.000000000000    
     +, 0.2682926654816    
     +, 0.3170731663704    
     +, 0.7073170542717    
     +, 0.7560975551605    
     +/
      DATA SIGDEL / 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.1999999994950E-05
     +, 0.1999999994950E-05
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +/
      DATA SIGA /  14.57433603529    
     +, -15.93532436156    
     +, -13.24628335221    
     +, -14.31193855410    
     +, -12.67339684488    
     +,  18.19876051780    
     +, -5.270493486725    
     +, -5.107670990828    
     +, -9.553262933901    
     +,  43.34150727448    
     +,  55.91366786432    
     +, -29.38037318563    
     +/

      v(1)= abs(x)
      if(V(1).gt.20.)V(1)=20.
      
      HQUADF = 0.
      DO 20 J = 1, NPAR
         HQDJ = 0.
         DO 10 I = 1, NDIM
            VV = (V (I) - SIGVMI (I)) / SIGVT (I)
            HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
   10    CONTINUE
         HQDJ = HQDJ + SIGDEL (J) ** 2
         HQDJ = SQRT (HQDJ)
         HQUADF = HQUADF + SIGA (J) * HQDJ
   20 CONTINUE
      IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)

      risyeta2=HQUADF* 1e-4

      END
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***

      FUNCTION risy_cog(x)
      DOUBLE PRECISION V( 1)
      INTEGER NPAR, NDIM, IMQFUN, I, J
      DOUBLE PRECISION HQDJ, VV, VCONST
      DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
      DOUBLE PRECISION SIGV(  10, 1)
      DOUBLE PRECISION SIGDEL(  10)
      DOUBLE PRECISION SIGA(  10)
      DATA NPAR, NDIM, IMQFUN /   10,    1,    1/
      DATA VCONST /  0.000000000000    /
      DATA SIGVMI / -20.50000000000    /
      DATA SIGVT /  41.00000000000    /
      DATA SIGV / 0.1585365831852    
     +, 0.8414633870125    
     +,  0.000000000000    
     +,  1.000000000000    
     +, 0.4634146094322    
     +, 0.5121951103210    
     +, 0.5609756112099    
     +, 0.6585365533829    
     +, 0.7073170542717    
     +, 0.3414633870125    
     +/
      DATA SIGDEL / 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.1999999994950E-05
     +, 0.1999999994950E-05
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.9756097197533E-01
     +, 0.1951219439507    
     +/
      DATA SIGA /  23.73833445988    
     +,  24.10182100013    
     +,  1.865894323190    
     +,  1.706006262931    
     +, -1.075607857202    
     +, -22.11489493403    
     +,  1.663100707801    
     +,  4.089852595440    
     +, -4.314993873697    
     +, -2.174479487744    
     +/
      
      V(1)=abs(x)
      if(V(1).gt.20.)V(1)=20.

      HQUADF = 0.
      DO 20 J = 1, NPAR
         HQDJ = 0.
         DO 10 I = 1, NDIM
            VV = (V (I) - SIGVMI (I)) / SIGVT (I)
            HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
   10    CONTINUE
         HQDJ = HQDJ + SIGDEL (J) ** 2
         HQDJ = SQRT (HQDJ)
         HQUADF = HQUADF + SIGA (J) * HQDJ
   20 CONTINUE
      IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)

      risy_cog=HQUADF* 1e-4
      
      END
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      FUNCTION risx_cog(x)
      DOUBLE PRECISION V( 1)
      INTEGER NPAR, NDIM, IMQFUN, I, J
      DOUBLE PRECISION HQDJ, VV, VCONST
      DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
      DOUBLE PRECISION SIGV(  15, 1)
      DOUBLE PRECISION SIGDEL(  15)
      DOUBLE PRECISION SIGA(  15)
      DATA NPAR, NDIM, IMQFUN /   15,    1,    1/
      DATA VCONST /  0.000000000000    /
      DATA SIGVMI / -20.50000000000    /
      DATA SIGVT /  41.00000000000    /
      DATA SIGV / 0.6097560748458E-01
     +, 0.8536584675312E-01
     +, 0.1341463327408    
     +, 0.2317073047161    
     +, 0.2804878056049    
     +, 0.3780487775803    
     +, 0.6219512224197    
     +, 0.7195121645927    
     +, 0.7682926654816    
     +, 0.8658536076546    
     +, 0.9146341085434    
     +, 0.9390243291855    
     +,  0.000000000000    
     +,  1.000000000000    
     +, 0.5121951103210    
     +/
      DATA SIGDEL / 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.4878048598766E-01
     +, 0.1999999994950E-05
     +, 0.1999999994950E-05
     +, 0.9756097197533E-01
     +/
      DATA SIGA /  31.95672945139    
     +, -34.23286209245    
     +, -6.298459168211    
     +,  10.98847700545    
     +,-0.3052213535054    
     +,  13.10517991464    
     +,  15.60290821679    
     +, -1.956118448507    
     +,  12.41453816720    
     +, -7.354056408553    
     +, -32.32512668778    
     +,  30.61116178966    
     +,  1.418505329236    
     +,  1.583492573619    
     +, -18.48799977042    
     +/

      V(1)=abs(x)
      if(V(1).gt.20.)V(1)=20.

      HQUADF = 0.
      DO 20 J = 1, NPAR
         HQDJ = 0.
         DO 10 I = 1, NDIM
            VV = (V (I) - SIGVMI (I)) / SIGVT (I)
            HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
   10    CONTINUE
         HQDJ = HQDJ + SIGDEL (J) ** 2
         HQDJ = SQRT (HQDJ)
         HQUADF = HQUADF + SIGA (J) * HQDJ
   20 CONTINUE
      IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)

      risx_cog = HQUADF * 1e-4

      END


*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function pfacorr(ic,angle) 
*--------------------------------------------------------------
*     this function returns the landi correction for this cluster
*--------------------------------------------------------------
      include 'commontracker.f'
      include 'calib.f'
      include 'level1.f'

      real angle
      integer iview,lad

      iview = VIEW(ic)            
      lad = nld(MAXS(ic),VIEW(ic))

*     find angular bin
*     (in futuro possiamo pensare di interpolare anche sull'angolo)
      do iang=1,nangbin
         if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
            iangle=iang
            goto 98
         endif
      enddo
      if(DEBUG.eq.1)
     $     print*,'pfacorr *** warning *** angle out of range: ',angle
      if(angle.le.angL(1))iang=1
      if(angle.ge.angR(nangbin))iang=nangbin
 98   continue                  !jump here if ok

      pfacorr = fcorr(iview,lad,iang)

      if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr


 100  return
      end
