*************************************************************************
*     
*     Subroutine cncomp.f 
*     
*     iterates common noise computation subroutine (./cnoise.f) and cluster
*     cutting subroutine (./cutcn.f) till no more clusters are found
*          
*************************************************************************

      subroutine cncomp(i,j,errflag)    !(view, VA1)

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

      integer errflag           !error flag to mark no signal free VA1
      integer clstr_old(nstrips_va1) !flag storage vector
      real signal(nstrips_va1)  !"signal" (=adc-ped) value storage vector
      real smean, ssigma        !"signal" mean and sigma
      real cut                  !"strange" strip exclusion cut
      integer newclstr          !flag to warn about new found clusters to be
c                               ! excluded from common noise computation 

c------------------------------------------------------------------------
c
c     variables initialization
c
c------------------------------------------------------------------------
      do k=1,nstrips_va1        !loops on strips
         clstr(i,j,k)=1         !initializes signal affected strips flag
         clstr_old(k)=1         !initializes signal affected strips storage 
         strange(i,j,k)=1       !initializes unusually high or low signal 
      enddo                     ! affected strips flag

      newclstr=1                !flag to warn about new found signal 
c                               ! affected strips
c------------------------------------------------------------------------
c
c     high or low signal affected strips exclusion: computes "signal" (=adc-ped)
c     mean value and sigma, and cuts from common noise computation strips
c     whose ABS(signal) exceeds scut*sigma
c
c------------------------------------------------------------------------
      countme=0                 !???
 666  continue                  !???

      smean=0.                  !initialization
      ssigma=0.
      nstr=0
      
      do k=1,nstrips_va1
         nstr = nstr + strange(i,j,k) !uses only
         if(mod(i,2).eq.1) then ! ---> Y view
            signal(k) = - (DBLE(adc(i,j,k)) - pedestal(i,j,k)) !negative signal
         else                   ! ---> X view
            signal(k) =    DBLE(adc(i,j,k)) - pedestal(i,j,k) !positive signal
         endif
         smean = smean + signal(k)*strange(i,j,k)
         ssigma = ssigma + (signal(k)**2)*strange(i,j,k)
      enddo
      
      smean=smean/nstr          !strips value distribution mean      
      ssigma=SQRT((ssigma/nstr)-smean**2) !strips value distribution sigma
      
      cut=scut*ssigma           !exclusion cut
      
      nco=0
      nbo=0
      do k=1,nstrips_va1
         if(ABS(signal(k)-smean).gt.cut) then 
            strange(i,j,k)=0    !marks strips exceeding cut
c            print*,i,j,k,signal(k),smean
         endif
         nco=nco+strange(i,j,k)
         nbo=nbo+bad(i,j,k)
      enddo                     ! in order not to use them in CN computation

c$$$      if(i.eq.12.and.(j.eq.2.or.j.eq.3))then
c$$$         print*,'view ',i,' vk ',j
c$$$         print*,'ADC (1-51-128) = ',adc(i,j,1),adc(i,j,52),adc(i,j,128)
c$$$         print*,'<ADC-PED> = ',smean
c$$$         print*,'s         = ',ssigma
c$$$         print*,'nstrange  = ',128-nco
c$$$         print*,'nbad      = ',128-nbo
c$$$      endif

      countme = countme + 1         !???
      if (countme.le.3) goto 666 !???

c------------------------------------------------------------------------
c     
c     common noise computation
c     
c-----------------------------------------------------------------------
*     loops on this VA1 till no new signal affected strips are found
      do while(newclstr.eq.1)  

         newclstr=0             !to exit from loop if no new cluster is found
         
         errflag=0
         call cnoise(i,j,errflag) !(view, VA1, error flag) computes cn
         if(errflag.eq.1) goto 10 !goes to next VA1: this one has no signal-free strips...
         
         call cutcn(i,j)        !(view, VA1) excludes clusters from cn computation
         
         ncs=0                  !initializes number of strips not excluded by cncut
         do k=1,nstrips_va1     !loops on strips
*           checks if there are new found clusters, and if so sets 
            if(clstr(i,j,k).ne.clstr_old(k)) then 
               newclstr=1                      
               clstr_old(k)=clstr(i,j,k)  
            endif               
            iok=strange(i,j,k)*bad(i,j,k)*clstr(i,j,k)
            ncs=ncs+iok         !counts number of good strips for cn computation
         enddo

      enddo                     !ends do while 

 10   continue

      return
      end




*************************************************************************
*     
*     Subroutine cnoise.f!DA COMMENTARE!???
*     
*     uses adc(nviews,nva1_view,nstrips_va1) and 
*     pedestal(nviews,nva1_view,nstrips_va1) variables to compute common noise,
*     and fills cn(nviews,nva1_view) variable. in the computation only 
*     not-bad and not-signal-affected strips are used
*     (bad(nviews,nva1_view,nstrips_va1) and 
*     clstr(nviews,nva1_view,nstrips_va1) flags)
*     
*     needs:
*     - ./common_calib.f
*     
*     to be called inside ./cncomp.f
*     
*************************************************************************

      subroutine cnoise(i,j,gulp) !(view, VA1)

      include 'commontracker.f'
      include 'level0.f'
      include 'level1.f'
      include 'common_reduction.f'
      include 'calib.f'
      

      integer gulp              !error flag
      
      ncn=0                     !number of strips in cn computation
      cn(i,j)=0                 !initializes cn variable
      cnrms(i,j)=0              !initializes cn rms
      cnn(i,j)=0                !initialize cn flag 

      do k=1,nstrips_va1        !loops on strips
*        tags strange, bad or signal-affected strips
         iok = strange(i,j,k)*bad(i,j,k)*clstr(i,j,k) 
         cn(i,j) = cn(i,j) + (DBLE(adc(i,j,k)) - pedestal(i,j,k))*iok 
         cnrms(i,j) = cnrms(i,j) 
     $        + (DBLE(adc(i,j,k)) - pedestal(i,j,k))
     $        *(DBLE(adc(i,j,k)) - pedestal(i,j,k))*iok 
         ncn = ncn + iok            !counts number of strips in cn computation
      enddo
      
      if(ncn.lt.NSTRIPMIN) then         !no signal free strips on this VA1...
         if(ncn.eq.0)then
            if(debug)print*,' WARNING - cnoise: ',
     $        'no strips for CN computation on VA1 ',j,
     $        ', VIEW ',i,'  >>> FAILED '
         else
            if(debug)print*,' WARNING - cnoise: ',
     $        'less than ',NSTRIPMIN
     $           ,' strips for CN computation on VA1 ',j,
     $        ', VIEW ',i,'  >>> FAILED '
         endif
         gulp=1
         cnn(i,j) = 0
      else
         cn(i,j)=cn(i,j)/DBLE(ncn) !<<<< computes common noise
         cnrms(i,j)= SQRT( cnrms(i,j)/DBLE(ncn) - cn(i,j)**2 ) 
         cnn(i,j) = ncn
         gulp=0                 
         if(debug.and.ABS(cn(i,j)).gt.1000)
     $        print*,'Event ',eventn(1)
     $        ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn
      endif

      return
      end


*************************************************************************
*     
*     Subroutine cutcn.f!DA COMMENTARE!???
*     
*     excludes strips with particle signals and/or noisy strips from common
*     noise calculation, marking their clstr(nviews,nva1_view,nstrips_va1)
*     flag:
*     clstr=0 ---> not to be used in CN computation
*     clstr=1 ---> to be used in CN computation
*     
*     needs:
*     - ./common_calib.f
*     
*     to be called inside ./cncomp.f
*     
*************************************************************************
      
      subroutine cutcn(i,j)     !(view, VA1)

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


      integer skip              !used to skip strips (see later...)

      integer kr, kl            !position indexes to check signal affected
                                ! strips on right and left side of cluster
                                ! seed
      integer ir, il            !flags to exit loop on reaching VA1 extremes

      real valuec                !cluster seed signal
      real cut,stripcut         !cluster seed cut

      real valuel, valuer       !left and right strips signal
      real stripcnincut         !strip include cut

      skip = 0                  !initializes skip

      do k=1,nstrips_va1        !loops on strips searching for cluster seeds

         if(k.le.skip) goto 20  !continues only if k strip has not been
                                ! checked yet

         clstr(i,j,k)=1         !reinitializes strip to be used in CN!???
                                ! computation, in order to be able to exclude
                                ! different strips at every CN computation loop

c------------------------------------------------------------------------
c     
c     selects cut according to view
c     
c------------------------------------------------------------------------
         if(mod(i,2).eq.1) then !odd strip ---> Y view
            valuec= - (DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k)) !negative signal
            cut=clcuty          !sets Y cut to find cluster seeds
         else                   !even strip ---> X view
            valuec= DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k) !positive signal
            cut=clcutx          !sets X cut to find cluster seeds
         endif


c------------------------------------------------------------------------
c     
c     seeks clusters
c     
c------------------------------------------------------------------------
         stripcut=cut*sigma(i,j,k) !cluster seed cut

c     if(ABS(valuec).gt.stripcut) then !checks if signal exceeds threshold!???
         if(valuec.gt.stripcut) then !checks if signal exceeds threshold

c$$$            print*,'cut',i,j,k,valuec,stripcut,adc(i,j,k),cn(i,j)
c$$$     $           ,pedestal(i,j,k) !???

            clstr(i,j,k)=0      !if so, marks this strip as part of a 
                                ! cluster

c------------------------------------------------------------------------
c     after finding a cluster seed, checks also adiacent strips, and marks 
c     the ones exceeding cnincut
c------------------------------------------------------------------------
            kr=k                !initializes position indexes to be equal to
            kl=k                ! cluster seed position

            ir=0                !initialize flags used to exit from
            il=0                ! inclusion loop

            do while (il.eq.0.or.ir.eq.0) !shifts left and right from
                                ! cluster seed till it finds a strip below
                                ! the threshold, or till it reaches first or 
                                ! last VA1 strip
               kr=kr+1          !position index for strips on right side of 
                                ! cluster seed
               kl=kl-1          !and for left side

c------------------------------------------------------------------------
c     checks for last or first strip
c------------------------------------------------------------------------
               if(kr.gt.nstrips_va1.and.ir.eq.0) then !when index goes 
                  ir=1          ! beyond last VA1 strip, change ir flag in 
                                ! order to "help" exiting from loop
                  skip=nstrips_va1+1 !sets skip beyond last strip: all 
                                ! strips on the right have been included in
                                ! the cluster, so skips all next strips 
                                ! (goto 20 condition is now always true)
               endif

               if(kl.lt.1.and.il.eq.0) then !idem when index goes beyond
                  il=1          ! first strip
               endif

c     P.S.: the "....and.i#.eq.0" term in above conditions is needed. In
c     fact, even if I reach a under-cut strip on the right (so I get ir=1),
c     the "do while loop" continues till such strip will be found on the 
c     left too.
c     Thus kl and kr (! too) keep increasing, and it can happen kr gets
c     greater than nstrips_va1 before kl reaches a under-cut strip. In this
c     case it would pass this "if condition", so setting skip=nstrips_va1+1
c     and skipping right strips never checked, if the "....and.i#.eq.0" term
c     weren't the: instead, including this part it won't pass it
c     because when I found reach the last VA1 strip on the right I set ir=1.
c     (AAAAAAHHHHHHHHH!!!!!!!!!!!)

c------------------------------------------------------------------------
c     marks strips exceeding inclusion cut
c------------------------------------------------------------------------
c     for right strips (kr index)
               if(ir.eq.0) then !if last strip or last over-cut strip has
                                ! not been reached

                  if(mod(i,2).eq.1) then !Y view
                     valuer= - (DBLE(adc(i,j,kr))-cn(i,j) !puts in valuer
     $                    -pedestal(i,j,kr)) ! right strip value
                  else          !X view
                     valuer=DBLE(adc(i,j,kr))-cn(i,j)-pedestal(i,j,kr)
                  endif

                  stripcnincut=cnincut*sigma(i,j,kr) !defines include cut
c     if(ABS(valuer).gt.stripcnincut) then !marks right strip if it !???
                  if(valuer.gt.stripcnincut) then !marks right strip if it 
                     clstr(i,j,kr)=0 !exceedes include cut
c$$$      print*,'inclcut_r',i,j,kr,valuer,stripcnincut
c$$$     $                    ,adc(i,j,kr),cn(i,j),pedestal(i,j,kr) !???
                  else
                     ir=1       !otherwise cluster ends and ir flag =1
                                ! signals it
                     skip=kr    !putting skip=kr, next k to be checked is
                                ! k=kr
                  endif 

               endif

c     for left strips (kl index)
               if(il.eq.0) then !if first strip or last over-cut strip has
                                ! not been reached

                  if (mod(i,2).eq.1) then !Y view
                     valuel= - (DBLE(adc(i,j,kl))-cn(i,j) !puts in valuel
     $                    -pedestal(i,j,kl)) ! left strip value
                  else          !X view
                     valuel=DBLE(adc(i,j,kl))-cn(i,j)-pedestal(i,j,kl)
                  endif

                  stripcnincut=cnincut*sigma(i,j,kl) !defines include cut
c     if(ABS(valuel).gt.stripcnincut) then !marks left strip if it!???
                  if(valuel.gt.stripcnincut) then !marks left strip if it
                     clstr(i,j,kl)=0 !exceedes include cut
c$$$      print*,'inclcut_l',i,j,kl,valuel,stripcnincut
c$$$     $                    ,adc(i,j,kl),cn(i,j),pedestal(i,j,kl) !???
                  else
                     il=1       !otherwise cluster ends and il flag =1
                                ! signals it
                  endif

               endif

            enddo               !ends lateral strips loop

         endif                  !ends cluster seed condition

 20      continue               !comes here if next strip on the right has 
                                ! already been included in a cluster

      enddo                     !ends principal strip loop

      return
      end
