*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
*     this file contains all subroutines and functions
*     that are needed for position finding algorithms
*     
*
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***


*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function pfa_eta(ic,angle)
*--------------------------------------------------------------
*     this function returns the position (in strip units)
*     it calls: 
*     - pfa_eta2(ic,angle)
*     - pfa_eta3(ic,angle)
*     - pfa_eta4(ic,angle)
*     according to the angle
*--------------------------------------------------------------
      include 'commontracker.f'
c      include 'calib.f'
      include 'level1.f'
      
      pfa_eta = 0

      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
      
         pfa_eta = pfa_eta2(ic,angle)
   
      else                      !X-view

         if(abs(angle).le.10.)then
            pfa_eta = pfa_eta2(ic,angle)
         elseif(abs(angle).gt.10..and.abs(angle).le.15.)then
            pfa_eta = pfa_eta3(ic,angle)
         elseif(abs(angle).gt.15.)then
            pfa_eta = pfa_eta4(ic,angle)
         endif
            
      endif


 100  return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function ris_eta(ic,angle)
*--------------------------------------------------------------
*     this function returns the average spatial resolution 
*     (in cm) for the ETA algorithm (function pfa_eta(ic,angle))
*     it calls: 
*     - risx_eta2(angle)
*     - risy_eta2(angle)
*     - risx_eta3(angle)
*     - risx_eta4(angle)
*     according to the angle
*--------------------------------------------------------------
      include 'commontracker.f'
c      include 'calib.f'
      include 'level1.f'

c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
      
      ris_eta = 0

      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
      
         ris_eta = risy_eta2(angle)
         if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)

      else                      !X-view

         if(abs(angle).le.10.)then
            ris_eta = risx_eta2(angle)
         elseif(abs(angle).gt.10..and.abs(angle).le.15.)then
            ris_eta = risx_eta3(angle)
         elseif(abs(angle).gt.15..and.abs(angle).le.21.)then
            ris_eta = risx_eta4(angle)
         elseif(abs(angle).gt.21.)then
            ris_eta = risx_eta4(21.)
         endif
            
      endif

c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'
c$$$     $     ,' -->',ris_eta


 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
      
         fbad_eta = fbad_cog(2,ic)

      else                      !X-view

         if(abs(angle).le.10.)then
            fbad_eta = fbad_cog(2,ic)
         elseif(abs(angle).gt.10..and.abs(angle).le.15.)then
            fbad_eta = fbad_cog(3,ic)
         elseif(abs(angle).gt.15.)then
            fbad_eta = fbad_cog(4,ic)
         endif
            
      endif

      return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c*****************************************************
cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
c*****************************************************
c      real function pfa_eta2(cog2,view,lad,angle)
      real function pfa_eta2(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

c      logical DEBUG
c      common/dbg/DEBUG

      iview = VIEW(ic)              !(1)
      lad = nld(MAXS(ic),VIEW(ic)) !(1)
      cog2 = cog(2,ic)             !(1)
      pfa_eta2=cog2

*     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)
     $     print*,'pfa_eta2 *** warning *** angle out of range: ',angle
      if(angle.lt.angL(1))iang=1
      if(angle.gt.angR(nangbin))iang=nangbin
 98   continue                  !jump here if ok


c$$$*     find extremes of interpolation
c$$$      iflag=0
c$$$*     --------------------------------
c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then
c$$$c         print*,'pfa_eta2 *** warning *** argument out of range: ',cog2
c$$$*         goto 100
c$$$*     ----------------------------------------------
c$$$*     non salto piu`, ma scalo di 1 o -1
c$$$*     nel caso si tratti di un cluster 
c$$$*     in cui la strip con il segnale massimo non coincide
c$$$*     con la strip con il rapposto s/n massimo!!!
c$$$*     ----------------------------------------------
c$$$         if(cog2.lt.eta2(1,iang))then !temp
c$$$            cog2=cog2+1.        !temp
c$$$            iflag=1             !temp
c$$$         else                   !temp
c$$$            cog2=cog2-1.        !temp
c$$$            iflag=-1            !temp
c$$$         endif                  !temp
c$$$c         print*,'shifted >>> ',cog2
c$$$      endif

      iadd=0
 10   continue
      if(cog2.lt.eta2(1,iang))then
         cog2 = cog2 + 1
         iadd = iadd + 1
         goto 10
      endif
 20   continue
      if(cog2.gt.eta2(netaval,iang))then
         cog2 = cog2 - 1
         iadd = iadd - 1
         goto 20
      endif

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

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

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

      if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'
     $     ,cog2-iadd,' -->',pfa_eta2


 100  return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c*****************************************************
cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
c*****************************************************
c      real function pfa_eta3(cog3,view,lad,angle)
      real function pfa_eta3(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

c      logical DEBUG
c      common/dbg/DEBUG


      iview = VIEW(ic)              !(1)
      lad = nld(MAXS(ic),VIEW(ic)) !(1)
      cog3 = cog(3,ic)             !(1)
      pfa_eta3=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)
     $     print*,'pfa_eta3 *** warning *** angle out of range: ',angle
      if(angle.lt.angL(1))iang=1
      if(angle.gt.angR(nangbin))iang=nangbin
 98   continue                  !jump here if ok


c$$$*     find extremes of interpolation
c$$$      iflag=0
c$$$*     --------------------------------
c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then
c$$$*     ----------------------------------------------
c$$$*     non salto piu`, ma scalo di 1 o -1
c$$$*     nel caso si tratti di un cluster 
c$$$*     in cui la strip con il segnale massimo non coincide
c$$$*     con la strip con il rapposto s/n massimo!!!
c$$$*     ----------------------------------------------
c$$$         if(cog2.lt.eta2(1,iang))then !temp
c$$$            cog2=cog2+1.        !temp
c$$$            iflag=1             !temp
c$$$         else                   !temp
c$$$            cog2=cog2-1.        !temp
c$$$            iflag=-1            !temp
c$$$         endif                  !temp
c$$$c         print*,'shifted >>> ',cog2
c$$$      endif

      
      iadd=0
 10   continue
      if(cog3.lt.eta3(1,iang))then
         cog3 = cog3 + 1
         iadd = iadd + 1
         goto 10
      endif
 20   continue
      if(cog3.gt.eta3(netaval,iang))then
         cog3 = cog3 - 1
         iadd = iadd - 1
         goto 20
      endif

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

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

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

      if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'
     $     ,cog3-iadd,' -->',pfa_eta3

 100  return
      end

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c*****************************************************
cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
c*****************************************************
c      real function pfa_eta4(cog4,view,lad,angle)
      real function pfa_eta4(ic,angle) !(1)
*--------------------------------------------------------------
*     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

c      logical DEBUG
c      common/dbg/DEBUG

      iview = VIEW(ic)              !(1)
      lad = nld(MAXS(ic),VIEW(ic)) !(1)
      cog4=cog(4,ic)               !(1)
      pfa_eta4=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)
     $     print*,'pfa_eta4 *** warning *** angle out of range: ',angle
      if(angle.lt.angL(1))iang=1
      if(angle.gt.angR(nangbin))iang=nangbin
 98   continue                  !jump here if ok


c$$$*     find extremes of interpolation
c$$$      iflag=0
c$$$*     --------------------------------
c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then
c$$$*     ----------------------------------------------
c$$$*     non salto piu`, ma scalo di 1 o -1
c$$$*     nel caso si tratti di un cluster 
c$$$*     in cui la strip con il segnale massimo non coincide
c$$$*     con la strip con il rapposto s/n massimo!!!
c$$$*     ----------------------------------------------
c$$$         if(cog2.lt.eta2(1,iang))then !temp
c$$$            cog2=cog2+1.        !temp
c$$$            iflag=1             !temp
c$$$         else                   !temp
c$$$            cog2=cog2-1.        !temp
c$$$            iflag=-1            !temp
c$$$         endif                  !temp
c$$$c         print*,'shifted >>> ',cog2
c$$$      endif

      
      iadd=0
 10   continue
      if(cog4.lt.eta4(1,iang))then
         cog4 = cog4 + 1
         iadd = iadd + 1
         goto 10
      endif
 20   continue
      if(cog4.gt.eta4(netaval,iang))then
         cog4 = cog4 - 1
         iadd = iadd - 1
         goto 20
      endif

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

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

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

      if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'
     $     ,cog4-iadd,' -->',pfa_eta4

 100  return
      end



*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function cog0(ncog,ic)
*-------------------------------------------------
*     this function returns 
*
*     - the Center-Of-Gravity of the cluster IC 
*     evaluated using NCOG strips, 
*     calculated relative to MAXS(IC)
*     
*     - zero in case that not  enough strips 
*     have a positive signal
*     
*     NOTE:
*     This is the old definition, used by Straulino.
*     The new routine, according to Landi, 
*     is COG(NCOG,IC)
*-------------------------------------------------


      include 'commontracker.f'
      include 'level1.f'
      
*     --> signal of the central strip
      sc = CLSIGNAL(INDMAX(ic)) !center

*     signal of adjacent strips
*     --> left
      sl1 = 0                  !left 1
      if(
     $     (INDMAX(ic)-1).ge.INDSTART(ic)
     $     )
     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))

      sl2 = 0                  !left 2
      if(
     $     (INDMAX(ic)-2).ge.INDSTART(ic)
     $     )
     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))

*     --> right
      sr1 = 0                  !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 = max(0.,CLSIGNAL(INDMAX(ic)+1))

      sr2 = 0                  !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 = max(0.,CLSIGNAL(INDMAX(ic)+2))
      
************************************************************
*     COG computation
************************************************************

c      print*,sl2,sl1,sc,sr1,sr2

      COG = 0.
         
      if(sl1.gt.sr1.and.sl1.gt.0.)then
         
         if(ncog.eq.2.and.sl1.ne.0)then
            COG = -sl1/(sl1+sc)        
         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
            COG = (sr1-sl1)/(sl1+sc+sr1) 
         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) 
         else
            COG = 0.
         endif
         
      elseif(sl1.le.sr1.and.sr1.gt.0.)then

         if(ncog.eq.2.and.sr1.ne.0)then
            COG = sr1/(sc+sr1)            
         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
            COG = (sr1-sl1)/(sl1+sc+sr1) 
         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) 
         else
            COG = 0.
         endif

      endif

      COG0 = COG

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

      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'
      
c      logical DEBUG
c      common/dbg/DEBUG


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

*     --> signal of the central strip
         sc = CLSIGNAL(INDMAX(ic)) !center
*     signal of adjacent strips
         sl1 = 0                !left 1
         if(
     $        (INDMAX(ic)-1).ge.INDSTART(ic)
     $        )
     $        sl1 = CLSIGNAL(INDMAX(ic)-1)
         
         sl2 = 0                !left 2
         if(
     $        (INDMAX(ic)-2).ge.INDSTART(ic)
     $        )
     $        sl2 = CLSIGNAL(INDMAX(ic)-2)
         
         sr1 = 0                !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 = 0                !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.
         
         if(ncog.eq.2)then
            if(sl1.gt.sr1)then
               COG = -sl1/(sl1+sc)        
            elseif(sl1.le.sr1)then
               COG = sr1/(sc+sr1)            
            endif
         elseif(ncog.eq.3)then
            COG = (sr1-sl1)/(sl1+sc+sr1) 
         elseif(ncog.eq.4)then
            if(sl2.gt.sr2)then
               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) 
            elseif(sl2.le.sr2)then
               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) 
            endif
         else
            print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
            print*,'                 (NCOG must be <= 4)'
            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  
         mu=0
         do i=istart,istop
            ipos=i-INDMAX(ic)
            ivk=nvk(MAXS(ic)+ipos)
            is=nst(MAXS(ic)+ipos)
*            print*,'******************',istart,istop,ipos
*     $           ,MAXS(ic)+ipos,iv,ivk,is
            cut=incut*SIGMA(iv,ivk,is)
            if(CLSIGNAL(i).ge.cut)then
               COG = COG + ipos*CLSIGNAL(i)
               mu = mu + 1
c               print*,ipos,CLSIGNAL(i),incut,cut
            endif
         enddo
         COG=COG/DEDX(ic)
c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'
c     $        ,cog
         
      else
         
         COG=0
         print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
         print*,'                 (NCOG must be >= 0)'
         

      endif

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

      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
         f  = 4.
         si = 8.4
         incut=incuty
      else                      !X-view
         f  = 6.
         si = 3.9
         incut=incutx
      endif
      
      fbad_cog = 1.

      if (ncog.gt.0) then

*     --> signal of the central strip
         sc = CLSIGNAL(INDMAX(ic)) !center
         fsc = 1
         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f
*     --> 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)
            if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f
c         else
c            fsl1 = 0
         endif

         sl2 = 0                !left 2
         fsl2 = 1               !left 2
         if(
     $        (INDMAX(ic)-2).ge.INDSTART(ic)
     $        )then
            sl2 = CLSIGNAL(INDMAX(ic)-2)
            if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f
c         else
c            fsl2 = 0
         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)
            if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f
c         else
c            fsr1 = 0
         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)
            if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f
c         else
c            fsr2 = 0
         endif



************************************************************
*     COG computation
************************************************************

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

         
         iv=VIEW(ic)
         istart=INDSTART(IC)
         istop=TOTCLLENGTH
         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
         COG=0.
         SNU=0.
         SDE=0.
         do i=istart,istop
            ipos=i-INDMAX(ic)
            il=nvk(MAXS(ic)+ipos)
            is=nst(MAXS(ic)+ipos)
            cut=incut*SIGMA(iv,il,is)
            if(CLSIGNAL(i).gt.cut)then
               COG = COG + ipos*CLSIGNAL(i)
            endif
         enddo
         COG=COG/DEDX(ic)
         do i=istart,istop
            ipos=i-INDMAX(ic)
            il=nvk(MAXS(ic)+ipos)
            is=nst(MAXS(ic)+ipos)
            cut=incut*SIGMA(iv,il,is)
            if(CLSIGNAL(i).gt.cut)then
               fs=1
               if(BAD(iv,il,is).eq.0)fs=f
               SNU = SNU + fs*(ipos-COG)**2
               SDE = SDE + (ipos-COG)**2
            endif            
         enddo
         if(SDE.ne.0)FBAD_COG=SNU/SDE

      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


*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      real function fbad_cog0(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.
*
*     NB!!!
*     (this is the old version. It consider only the two 
*     strips with the greatest signal. The new one is 
*     fbad_cog(ncog,ic) )
*     
*-------------------------------------------------------

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

*     --> signal of the central strip
      sc = CLSIGNAL(INDMAX(ic)) !center

*     signal of adjacent strips
*     --> left
      sl1 = 0                  !left 1
      if(
     $     (INDMAX(ic)-1).ge.INDSTART(ic)
     $     )
     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))

      sl2 = 0                  !left 2
      if(
     $     (INDMAX(ic)-2).ge.INDSTART(ic)
     $     )
     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))

*     --> right
      sr1 = 0                  !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 = max(0.,CLSIGNAL(INDMAX(ic)+1))

      sr2 = 0                  !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 = max(0.,CLSIGNAL(INDMAX(ic)+2))


      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
         f  = 4.
         si = 8.4
      else                              !X-view
         f  = 6.
         si = 3.9
      endif

      fbad_cog = 1.
      f0 = 1
      f1 = 1
      f2 = 1
      f3 = 1   
      if(sl1.gt.sr1.and.sl1.gt.0.)then
         
         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f

         if(ncog.eq.2.and.sl1.ne.0)then
            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
            fbad_cog = 1.
         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
            fbad_cog = 1.
         else
            fbad_cog = 1.
         endif
         
      elseif(sl1.le.sr1.and.sr1.gt.0.)then


         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f

         if(ncog.eq.2.and.sr1.ne.0)then
            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
            fbad_cog = 1.
         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
            fbad_cog = 1.
         else
            fbad_cog = 1.
         endif

      endif

      fbad_cog0 = sqrt(fbad_cog)

      return
      end




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

      FUNCTION risx_eta2(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)

      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_eta2=HQUADF* 1e-4

      END

*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      FUNCTION risx_eta3(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)
      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_eta3 = HQUADF* 1e-4

      END
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      FUNCTION risx_eta4(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)

      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_eta4=HQUADF* 1e-4

      END
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
      FUNCTION risy_eta2(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)
      
      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_eta2=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)

      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)

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