*************************************************************************
*     
*     Program reduction.f
*     
*     - reads readraw.f output files: LEVEL0 ntuple, and ped, sig and bad histograms
*     - decodes raw data (DATATRACKER) using DSP ped, sig and bad values
*     - looks for clusters information using ped, sig and bad values from
*     DSP histograms 
*     - fills LEVEL1 ntuple
* 
*
*    
*     11/9/2005 modified by david fedele to include general variables
*    
*************************************************************************


      program reduction

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

c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------


      parameter(ncalibmax=50)

c     local variables
      character*74 data_file    !data file name
      character*74 data_dir     !data file name
      character*74 data_file_calib !calibration list file
      character*40 file_calib(ncalibmax)
      character*74 data_file_level0
      character*74 data_file_level1
      integer minevent          !first event to be analysed
      integer nevent            !number of events to be analysed
      character*24 processing_date

      integer which_calib_last

      parameter (lun_data_file=70) !data file id number
      parameter (lun_data_level0=71) !data file id number
      parameter (lun_data_level1=72) !data file id number
      parameter(lun_calib_list=66) !calibration list file id


      integer ind
      real value(nstrips_view)  !per trovare i cluster
      real clseedcut(nstrips_view)
      real clinclcut(nstrips_view)
      common/searchcluster/value,clseedcut,clinclcut,ind

      integer nshowers
      logical flag_shower
      common/shower/nshowers,flag_shower

      COMMON/QUEST/IQUEST(100)  !permette di ottenere ntuple funzionanti nonostante
                                ! il messaggio dei 64K di RZOUT...!???


c*****************************************************
cccccc 11/9/2005 modified by david fedele
      swcode1=202
c****************************************************
c------------------------------------------------------------------------
c     
c     HBOOK initialization
c     
c------------------------------------------------------------------------

      call HLIMIT(NWPAWC)


c------------------------------------------------------------------------
c     
c     reads input informations (through < go_reduction)
c     
c------------------------------------------------------------------------

      call fdate(processing_date)
      write(*,101)
     $     processing_date
 101  format(/
     $     ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $     ,'*                                                 *',/
     $     ,'* REDUCTION                                       *',/     
     $     ,'*                                                 *',/
     $     ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $     ,a24,/
     $     )

 111  format(a)
      print*,'Data directory:' 
      read(*,111)data_dir
      print*,data_dir
      print*,'File identifier: (DATE_NUM)'    
      read(*,*)data_file
      print*,data_file
c$$$      print*,
c$$$     $     data_dir(1:LNBLNK(data_dir))
c$$$     $     //'DW_'//data_file(1:LNBLNK(data_file))
      minevent=1
      print*,'Maximum number of events to be analized:' !20000
      read(*,*) nevent
      print*,nevent
      print*,'---------------------------------------------------'
 

c-------------------------------------------------------------
c-------------------------------------------------------------
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c     
c     reads calibration file
c     
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c-------------------------------------------------------------
c-------------------------------------------------------------

      print*,' '
      print*,'OPENING CALIBRATION-LIST FILE:'
 501  format(a,'DW_',a,'_calib.txt')
      write(data_file_calib,501)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
      print*,data_file_calib
      open(lun_calib_list,
     $     FILE=data_file_calib(1:LNBLNK(data_file_calib))
     $     ,STATUS='UNKNOWN'
     $     ,IOSTAT=iostat
     $     )

 113  format(i5,'    ',a25)

      do i=1,ncalibmax

         read(lun_calib_list,113,IOSTAT=iostat) n_cal_list
     $        ,file_calib(i)(1:LNBLNK(file_calib(i)))
         if(iostat.ne.0)then
            ncal=i-1
            goto 2000
         endif
         print*,n_cal_list,' - '
     $        ,file_calib(i)(1:LNBLNK(file_calib(i)))

      enddo

 2000 close(lun_calib_list)
      which_calib_last=0


c--------------------------------------------------------------
c--------------------------------------------------------------
c     
c     opening LEVEL 0 ntuple 
c     
c--------------------------------------------------------------
c--------------------------------------------------------------

 504  format(a,'DW_',a,'_level0.rz')
      write(data_file_level0,504)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
      print*,''
      print*,'OPENING LEVEL0 FILE:'
      print*,data_file_level0
      IQUEST(10)=65000          !permette di ottenere ntuple funzionanti nonostante
                                ! il messaggio dei 64K di RZOUT...!???
      call HROPEN(lun_data_level0,
     $     'LEVEL0',data_file_level0,'QP',4096,istat) !opens rz  
      if(istat.ne.0) goto 19      
      call HRIN(0,9999,0)       
      call access_level0
c      call HPRNTU(ntp_level0)  
      call HNOENT(ntp_level0,iemax0)
      print*,'  events',iemax0
     
c------------------------------------------------------------------------
c     
c     LEVEL 1 ntuple booking
c     
c------------------------------------------------------------------------

 503  format(a,'DW_',a,'_level1.rz')
      write(data_file_level1,503)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
      print*,''
      print*,'------------------------------------'
      print*,'  Creating LEVEL1 rz file'
      print*,'  ',data_file_level1
      print*,'------------------------------------'

C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     largest RZ file: IQUEST(10) records x LREC words x 4 byte
C     with the following settings: 65000 x 4096 x 4 = 1G
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      IQUEST(10)=65000          
c     permette di ottenere ntuple funzionanti nonostante
c     il messaggio dei 64K di RZOUT
      call HROPEN(lun_data_level1,
     $     'LEVEL1',data_file_level1,'QNP',4096,istat) !opens rz  
      if(istat.ne.0) goto 19      
      call book_level1

*     -------------------------------------------------------
*     STRIP MASK
*     -------------------------------------------------------
      call stripmask(data_file)


      nshowers=0                
      nocalib=0                
c     counts number of events with more than nclstrmax clusters
      maxevent=minevent+nevent-1
      nevent = 0
      nevent_good = 0
      
      do iev = minevent,MIN(iemax0,maxevent) !loop on events
         
         nevent = nevent + 1
         if(mod(nevent,100).eq.0)print*,nevent

         call init_level1

         call HCDIR('//LEVEL0',' ')
         call HGNT(ntp_level0,iev,ierr) !read event at LEVEL0
         if(ierr.ne.0)then 
            print*,'HGNT ---> error reading event ',iev
            good1 = .false.
            nev1 = nev1 + 1
            goto 200            !jump to fill n-tuple
         endif
c         if(mod(iev,100).eq.0) 
c     $        print*,'              event number: ',iev
C---------------------------------------------------

c*****************************************************
cccccc 11/9/2005 modified by david fedele
ccC     variables in blocks EVENT and CPU are anyway filled
C     variables in blocks GENERAL and CPU are anyway filled
c********************************************************
C     in order to mantain sincronization among 
C     events at different levels
C---------------------------------------------------
c*****************************************************
cccccc 11/9/2005 modified by david fedele
c         good1=good0
c         nev1=nev0 
c         pkt_type1=pkt_type
c         pkt_num1=pkt_num
c         obt1=obt
c         which_calib1=which_calib
         good1=good0
         nev1=nev0
         which_calib1=which_calib
         pkt_type1=pkt_type
         pkt_num1=pkt_num
         obt1=obt
         cpu_crc1=cpu_crc
         do iv=1,12
           crc1(iv)=crc(iv)
         enddo
c*******************************************************

C---------------------------------------------------
C     retrieve info about calib to be used
C---------------------------------------------------

         if(which_calib.ne.which_calib_last.and.
     $        which_calib.ne.0)then
            
            data_file_calib=
     $           data_dir(1:LNBLNK(data_dir))//
     $           file_calib(which_calib)
     $           (1:LNBLNK(file_calib(which_calib)))
c            print*,data_file_calib
            print*,''
            print*,
     $           '@ event ',nev0
     $           ,'   (CPU pkt N.',pkt_num,')'
            print*,'--> ',data_file_calib            
            IQUEST(10)=65000          
            call HROPEN(lun_data_file,
     $           'CALIB',data_file_calib,'QP',4096,istat) !opens
            if(istat.ne.0) goto 19      
            call HRIN(0,9999,0)       
            call fillpedsig
            do iview=1,nviews
               call HDELET(id_hi_bad+iview)
               call HDELET(id_hi_ped+iview)
               call HDELET(id_hi_sig+iview)
            enddo
            call HREND('CALIB')
            close(lun_data_file)            
            which_calib_last=which_calib


         elseif(which_calib.eq.0)then

            nocalib=nocalib+1
            good1=.false.

         endif

         if(.not.good1) goto 200 !jump to fill n-tuple

c--------------------------------------------------
c     read the variable DATATRACKER from LEVEL0
c     and fill the variable ADC (inverting view 11)
c--------------------------------------------------
         call filladc(iflag)
         if(iflag.ne.0)then
            good1=.false.
            print*,'event ',iev,' >>>>>  decode ERROR'
            goto 200
         endif
ccc /////////////////////////////////////////////////// da qui
c--------------------------------------------------
c     computes common noise for each VA1
c     (excluding strips affected by signal, 
c     tagged with the flag CLSTR)
c--------------------------------------------------
         do iv=1,nviews
            do ik=1,nva1_view
               cn(iv,ik)=0      !initializes cn variable
               if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik) 
            enddo
         enddo
c---------------------------------------------
c     loops on views, VA1 and strips, 
c     and computes strips signals using 
c     badstrip, pedestals, and 
c     sigma informations from histograms
c---------------------------------------------
         flag_shower = .false.
         ind=1                  !clsignal array index
         do iv=1,nviews         !loop on views
            do is=1,nstrips_view !loop on strips (1)
               if(mod(iv,2).eq.1) then 
C===  > Y view
                  value(is)= -(DBLE(adc(iv,nvk(is),nst(is))) 
     $                 -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
     $                 *mask(iv,nvk(is),nst(is))
                  clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is)) 
     $                 *mask(iv,nvk(is),nst(is))
                  clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is)) 
     $                 *mask(iv,nvk(is),nst(is))
               else            
C===  > X view
                  value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
     $                 -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
     $                 *mask(iv,nvk(is),nst(is))
                  clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
     $                 *mask(iv,nvk(is),nst(is))
                  clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
     $                 *mask(iv,nvk(is),nst(is))
               endif
            enddo               !end loop on strips (1)
            call search_cluster(iv)
            if(flag_shower.eqv..true.)then
               call init_level1               
               good1=.false.
               nev1 = nev0
               goto 200         !jump to next event
            endif
         enddo                  ! end loop on views
         do iv=1,nviews
            do ik=1,nva1_view
               cnev(iv,ik)=cn(iv,ik) !assigns computed CN to ntuple variables
            enddo
         enddo
         nevent_good = nevent_good + 1

 200     continue               
C---------------------------------------------
C     come here if GOOD1=.false.
C     or the event has too many clusters
C---------------------------------------------
         
         call HCDIR('//LEVEL1',' ')
         call HFNT(ntp_level1)  !fills LEVEL1  ntuple
         
ccc /////////////////////////////////////////////////// a qui

      enddo                     !end loop on events

c$$$      call HPRNT(ntp_level1)

      write(*,108)
     $     nevent
     $     ,nevent_good
     $     ,nshowers
     $     ,nclstrmax
     $     ,nocalib
 108  format(/
     $     ,'|------------------------------------------------|',/
     $     ,'|  Number of processed events                    |',/
     $     ,'|  Total   ',i8,'                              |'/
     $     ,'|  Good    ',i8,'                              |'/
     $     ,'|  (',i5,'  events with > ',i5,'  clusters)        |'/     
     $     ,'|  (',i5,'  events without calibration)           |'/     
     $     ,'|------------------------------------------------|',/
     $     )
      
      goto 9000

c------------------------------------------------------------------------
c     
c     no error exit
c     
c------------------------------------------------------------------------

      print*,' '
      print*,'REDUCTION SUCCESSFULLY COMPLETED'
      print*,' '
      print*,' '
      
      goto 9000                 !happy ending


c------------------------------------------------------------------------
c     
c     data file opening error
c     
c------------------------------------------------------------------------
      
 19   continue
      
      print*,' '
      print*,'reduction: ERROR OPENING DATA FILE: ',data_file
      print*,' '
      print*,' '
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     data file event error
c     
c------------------------------------------------------------------------
      
 20   continue
      
      print*,' '
      print*,'reduction: ERROR IN DATA FILE: ',data_file
      print*,'NO EVENTS FOUND'
      print*,' '
      print*,' '
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     data file eleboration level error
c     
c------------------------------------------------------------------------
      
 21   continue
      
      print*,' '
      print*,'reduction: ERROR IN DATA FILE: ',data_file
      print*,'LEVEL0= ',level0
      print*,' '
      print*,' '
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     level0 ntuple event reading error
c     
c------------------------------------------------------------------------
      
 22   continue
      
      print*,' '
      print*,'reduction: ERROR WHILE READING LEVEL0 NTUPLE, AT EVENT: ',
     $     iev
      print*,ierr
      print*,' '
      print*,' '
      
      goto 9000                 !the end


c------------------------------------------------------------------------
c     
c     closes files and exits
c     
c------------------------------------------------------------------------

 9000 continue
      

      call HCDIR('//LEVEL1',' ')
      call RZPURG(-1)   
      call HROUT(ntp_level1,ICYCLE,'T')
      call HREND('level1')
      close(lun_data_level1)

      call HCDIR('//LEVEL0',' ')
      call HREND('level0')
      close(lun_data_level0)

c$$$  close(lun_data_runinfo)


      stop
      end
***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
*
*
*
*
*
*
*
*
*
***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
      subroutine book_level1

      include '../common/commontracker.f'
      include '../common/level1.f'

      character*35 block1,block2,block3,
     $     block4,block5,block6,block7
     $     ,block8              

*     LEVEL1 ntuple:
      call HBNT(ntp_level1,'LEVEL1',' ')
      
c*****************************************************
cccccc 11/9/2005 modified by david fedele
c      call HBNAME(ntp_level1,'EVENT',good1,'GOOD1:L,NEV1:I')
c      call HBNAME(ntp_level1,'GENERAL',good1,'GOOD1:L,NEV1:I
c     $     ,WHIC_CALIB1:I,SWCODE1:I')
cccccc 06/10/2005 modified by elena vannuccini
      call HBNAME(ntp_level1,'GENERAL',good1,'GOOD1:L,NEV1:I
     $     ,WHIC_CALIB1:I,SWCODE1:I,CRC1(12):L')
c*********************************************************
      call HBNAME(ntp_level1,'CPU',pkt_type1
     $     ,'PKT_TYPE1:I
     $     ,PKT_NUM1:I
     $     ,OBT1:I'//
c********************************************************
cccccc 11/9/2005 modified by david fedele
c     $     ,WHICH_CALIB1:I')
     $     ',CPU_CRC1:L')
c********************************************************
 411  format('NCLSTR1:I::[0,',I4,']') !cluster number
 412  format(',VIEW(NCLSTR1):I::[0,',I2,']') !view the cluster belongs to
 413  format(',LADDER(NCLSTR1):I::[0,',I1,']') !ladder the cluster belongs to
 414  format(',MAXS(NCLSTR1):I::[0,',I5,']') !cluster stip carrying the largest 
                                ! signal value
 415  format(',MULT(NCLSTR1):I::[0,',I2,']') !cluster multiplicity
 416  format(',INDSTART(NCLSTR1):I::[0,',I5,']') !cluster starting point index in 
                                ! clsignal array
 417  format(',INDMAX(NCLSTR1):I::[0,',I5,']') !cluster maximum point index in
     $                          ! clsignal array
 418  format(',TOTCLLENGTH:I::[0,',I5,']') !sum of all clusters length
      write(block1,411) nclstrmax !maximum number of clusters per event
      write(block2,412) nviews  !number of views
      write(block3,413) nladders_view !number of ladders per view
      write(block4,414) nstrips_view !number of strips per view
      write(block5,415) nclstrp !maximum number of strips in a cluster
      write(block6,416) maxlength !maximum number of strip belonging to clusters
      write(block7,417) maxlength ! for the whole event
      write(block8,418) maxlength
      call HBNAME(ntp_level1,'CLUSTER',nclstr1,block1//block2//block3
     $     //block4//block5//',DEDX(NCLSTR1):r'//block6//block7
c*****************************************************
cccccc 11/9/2005 modified by david fedele
c     $     //block8//',CLSIGNAL(TOTCLLENGTH):R')
c     $     //block8//',CLSIGNAL(TOTCLLENGTH):R,CRC1(12):L')
cccccc 06/10/2005 modified by elena vannuccini
     $     //block8//',CLSIGNAL(TOTCLLENGTH):R')
c*****************************************************
 419  format('CNEV(',I4,',',I4,')') !common noise of the event for a certain view
      write(block8,419) nviews,nva1_view ! and VA1
      call HBNAME(ntp_level1,'CNOISE',cnev,block8)


      return
      end

c.............................................................

      subroutine init_level1
      
      include '../common/commontracker.f'
      include '../common/level1.f'
      include '../common/level0.f'

      nclstr1=0
      totCLlength=0
      do ic=1,nclstrmax
         view(ic)=0
         ladder(ic)=0
         indstart(ic)=0
         indmax(ic)=0
         maxs(ic)=0
         mult(ic)=0          
         dedx(ic)=0
      enddo
      do id=1,maxlength         !???
         clsignal(id)=0.
      enddo
      do iv=1,nviews
c*****************************************************
cccccc 11/9/2005 modified by david fedele
        crc1(iv)=.false.
c**************************************************
         do ik=1,nva1_view
            cnev(iv,ik)=0
         enddo
      enddo
      
      return
      end
*---***---***---***---***---***---***---***---***
*
*
*
*
*
*---***---***---***---***---***---***---***---***

      subroutine search_cluster(iv)

      include '../common/commontracker.f'
      include '../common/level1.f'
      include '../common/calib.f'

      integer ind
      real value(nstrips_view)  !per trovare i cluster
      real clseedcut(nstrips_view)
      real clinclcut(nstrips_view)
      common/searchcluster/value,clseedcut,clinclcut,ind

      integer nshowers
      logical flag_shower
      common/shower/nshowers,flag_shower

c     local variables
      integer rmax,lmax         !estremi del cluster
      integer rstop,lstop       !per decidere quali strip includere nel cluster 
                                ! oltre il seed
      integer first,last,diff   !per includere le strip giuste... !???

      integer multtemp          !temporary multiplicity variable

      integer CLlength          !lunghezza in strip del cluster

c      external nst

c------------------------------------------------------------------------
c     looks for clusters on each view
C     : CERCO STRIP SOPRA CLSEEDCUT, POI SCORRO A DX FINCHE' 
c     NON TROVO 
C     STRIP PIU' BASSA (in segnale/rumore) 
C     => L'ULTIMA DELLA SERIE CRESCENTE 
C     (LA PIU' ALTA) E' IL
C     CLUSTER SEED. POI SCORRO A SX E DX INCLUDENDO TUTTE 
C     LE STRIP (FINO A 17 AL 
C     MAX) CHE SUPERANO CLINCLCUT. 
C     QUANDO CERCO IL CLUSTER SEED SUCCESSIVO SALTO LA STRIP 
C     ADIACENTE A DESTRA
C     DELL'ULTIMO CLUSTER SEED (CHE SARA' NECESSARIAMENTE 
C     PIU' BASSA) E PRENDO 
C     COME SEED UNA STRIP SOLO SE IL SUO SEGNALE E' 
C     MAGGIORE DI QUELLO DELLA STRIP
C     PRECEDENTE (PRATICAMENTE PER EVITARE CHE L'ULTIMA 
C     STRIP DI UN GRUPPO DI STRIP
C     TUTTE SOPRA IL CLSEEDCUT VENGA AUTOMATICAMENTE PRESA 
C     COME SEED... DEVE ESSERE
C     PRESA SOLO SE IL CLUSTER E' DOUBLE PEAKED...)
c------------------------------------------------------------------------
c     6 ottobre 2003
c     Elena: CLSEEDCUT = 7 (old value 10)
c     Elena: CLINCLCUT = 4 (old value 5)

      iseed=-999                !cluster seed index initialization

      do jl=1,nladders_view     !1..3 !loops on ladders
         first=1+nstrips_ladder*(jl-1) !1,1025,2049
         last=nstrips_ladder*jl !1024,2048,3072
c     X views have 1018 strips instead of 1024
         if(mod(iv,2).eq.0) then 
            first=first+3
            last=last-3
         endif
         do is=first,last       !loop on strips in each ladder
            if(is.le.iseed+1) goto 220
c-----------------------------------------
c     after a cluster seed as been found,
c     look for next one skipping one strip on the right
c     (i.e. look for double peak cluster)
c-----------------------------------------
            if(is.ne.first) then
               if(value(is).le.value(is-1)) goto 220
            endif
c-----------------------------------------
c     skips cluster seed 
c     finding if strips values are descreasing (a strip
c     can be a cluster seed only if previous strip value
c     is lower)
c-----------------------------------------
            if(value(is).gt.clseedcut(is)) then
c-----------------------------------------
c     possible SEED...
c-----------------------------------------
               itemp=is
               if(itemp.eq.last) goto 230 !estremo...
c$$$               do while(value(itemp)
c$$$     $              /sigma(iv,nvk(itemp),nst(itemp))
c$$$     $              .le.value(itemp+1)
c$$$     $              /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!???
               do while(value(itemp)
     $              .le.value(itemp+1)) !BIAS: aggiustare il caso uguale!???
                  itemp=itemp+1
                  if(itemp.eq.last) goto 230 !stops if reaches last strip 
               enddo            ! of the ladder
 230           continue
c-----------------------------------------
c     fownd SEED!!!
c-----------------------------------------
               iseed=itemp     
c----------------------------------------------------------
c     after finding a cluster seed, checks also adjacent strips, 
C     and marks the ones exceeding clinclcut
c----------------------------------------------------------
               ir=iseed         !indici destro
               il=iseed         ! e sinistro
               
               rmax=ir          !estremo destro del cluster
               lmax=il          ! e sinistro
               
               rstop=0          !initialize flags used to exit from
               lstop=0          ! inclusion loop

               do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
                  ir=ir+1       !position index for strips on right side of 
                                ! cluster seed
                  il=il-1       !and for left side
c------------------------------------------------------------------------
c     checks for last or first strip of the ladder
c------------------------------------------------------------------------
                  if(ir.gt.last) then !when index goes beyond last strip 
                     rstop=1    ! of the ladder, change rstop flag in order
                                ! to "help" exiting from loop
                  endif
                  
                  if(il.lt.first) then !idem when index goes beyond
                     lstop=1    ! first strip of the ladder
                  endif
                  
c------------------------------------------------------------------------
c     check for clusters including more than nclstrp strips
c------------------------------------------------------------------------
                  if((rmax-lmax+1).ge.nclstrp) then
                     goto 210   !exits inclusion loop:
                                ! lmax and rmax maintain last value
                                ! NB .ge.!???
                  endif
c------------------------------------------------------------------------
c     marks strips exceeding inclusion cut
c------------------------------------------------------------------------
                  if(rstop.eq.0) then !if last strip of the ladder or last 
                                ! over-cut strip has not been reached
                     if(value(ir).gt.clinclcut(ir)) then !puts in rmax the 
                        rmax=ir ! last right over-cut strip
                     else
                        rstop=1 !otherwise cluster ends on right and rstop 
                     endif      ! flag=1 signals it
                  endif
                  if(lstop.eq.0) then
                     if(value(il).gt.clinclcut(il)) then
                        lmax=il
                     else
                        lstop=1
                     endif
                  endif

               enddo            !ends strip inclusion loop
 210           continue         !jumps here if more than nclstrp have been included
                     
               multtemp=rmax-lmax+1 !stores multiplicity in temp
                                ! variable. NB rmax and lmax can change later in 
                                ! order to include enough strips to calculate eta3
                                ! and eta4. so mult is not always equal to cllength
c------------------------------------------------------------------------
c     NB per essere sicuro di poter calcolare eta3 e eta4 devo includere
c     sempre e comunque le 2 strip adiacenti al cluster seed e quella 
c     adiacente ulteriore dalla parte della piu' alta fra queste due
c     (vedi oltre...)!???
c------------------------------------------------------------------------

c     nel caso di estremi del ladder...!???

c     ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder 
c     costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi
c     vado oltre (aggiungero' quindi strip a sx e dx in modo da poter calcolare
c     eta3e4)
               if((rmax-lmax+1).lt.4) then 

                  if(iseed.eq.first) then !estremi...
                     rmax=iseed+2 !NB in questo modo puo' anche capitare di
                     lmax=iseed ! includere strip sotto taglio di inclusione
                     goto 250   ! che non serviranno per eta3e4!???
                  endif
                  
                  if(iseed.eq.last) then !estremi...
                     rmax=iseed
                     lmax=iseed-2 !NB 2 e non 3, perche' altrimenti sarei in 
                     goto 250   ! ((rmax-lmax+1).lt.4).eq.false. !???
                  endif         !NMB questo e' l'unico caso di cllength=3!???
                  
                  if(iseed.eq.first+1) then !quasi estremi...
                     rmax=iseed+2
                     lmax=iseed-1
                     goto 250
                  endif
                  if(iseed.eq.last-1) then
                     rmax=iseed+1
                     lmax=iseed-2
                     goto 250
                  endif
c     se ho 4 o piu' strip --> se sono sui bordi esco, se sono sui quasi bordi
c     includo la strip del bordo
               else
                  
                  if(iseed.eq.first) goto 250 !estremi... non includo altro                  
                  if(iseed.eq.last) goto 250
                  if(iseed.eq.first+1) then !quasi estremi... mi assicuro di
                     lmax=first ! avere le strip adiacenti al seed 
                     if((rmax-lmax+1).gt.nclstrp) rmax=rmax-1 !NB effetto 
                     goto 250   ! coperta: se la lunghezza del cluster era gia'
                  endif         ! al limite (nclstrp), per poter aggiungere questa
                                ! strip a sinistra devo toglierne una a destra...!???
                  if(iseed.eq.last-1) then
                     rmax=last
                     if((rmax-lmax+1).gt.nclstrp) lmax=lmax+1
                     goto 250
                  endif                  
               endif
c------------------------------------------------------------------------
c     be sure to include in the cluster the cluster seed with its 2 adjacent 
c     strips, and the one adjacent to the greatest between this two strip, as the 
c     fourth one. if the strips have the same value (!) the fourth one is chosen
c     as the one having the greatest value between the second neighbors
c------------------------------------------------------------------------
               if(value(iseed+1).eq.value(iseed-1)) then
                  if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'
                     diff=(iseed+2)-rmax
                     if(diff.gt.0) then
                        rmax=rmax+diff
                        if((rmax-lmax+1).gt.nclstrp) then
                           lmax=rmax-nclstrp+1
                        endif
                     endif
                     diff=(iseed-1)-lmax
                     if(diff.lt.0) then
                        lmax=lmax+diff
                        if((rmax-lmax+1).gt.nclstrp) then
                           rmax=lmax+nclstrp-1
                        endif
                     endif
                  else
                     diff=(iseed-2)-lmax
                     if(diff.lt.0) then
                        lmax=lmax+diff
                        if((rmax-lmax+1).gt.nclstrp) then
                           rmax=lmax+nclstrp-1
                        endif
                     endif
                     diff=(iseed+1)-rmax
                     if(diff.gt.0) then
                        rmax=rmax+diff
                        if((rmax-lmax+1).gt.nclstrp) then
                           lmax=rmax-nclstrp+1
                        endif
                     endif
                     
                  endif
               elseif(value(iseed+1).gt.value(iseed-1)) then
c     !??? sposto il limite del cluster a destra per includere sempre le strip 
c     necessarie al calcolo di eta-i 
c     se il cluster diventa  troppo lungo lo accorcio a sinistra per avere non piu'
c     di nclstrp (in questo caso sono sicuro di aver gia' incluso le strip 
c     necessarie al calcolo di eta-i a sinistra, quindi se voglio posso uscire)
                  diff=(iseed+2)-rmax
                  if(diff.gt.0) then
                     rmax=rmax+diff
                     if((rmax-lmax+1).gt.nclstrp) then
                        lmax=rmax-nclstrp+1
c     goto 250 
                     endif
                  endif
                  diff=(iseed-1)-lmax
                  if(diff.lt.0) then
                     lmax=lmax+diff
                     if((rmax-lmax+1).gt.nclstrp) then
                        rmax=lmax+nclstrp-1
c     goto 250 !inutile!???
                     endif
                  endif
               else
                  diff=(iseed-2)-lmax
                  if(diff.lt.0) then
                     lmax=lmax+diff
                     if((rmax-lmax+1).gt.nclstrp) then
                        rmax=lmax+nclstrp-1
c     goto 250 
                     endif
                  endif
                  diff=(iseed+1)-rmax
                  if(diff.gt.0) then
                     rmax=rmax+diff
                     if((rmax-lmax+1).gt.nclstrp) then
                        lmax=rmax-nclstrp+1
c     goto 250 !inutile!???
                     endif
                  endif
               endif
 250           continue

c--------------------------------------------------------
c     fills ntuple variables
c--------------------------------------------------------
               nclstr1=nclstr1+1 !cluster number

               if(nclstr1.gt.nclstrmax) then !too many clusters for the event:
                  nshowers=nshowers+1 ! skips variables filling and go to next
                  good1=.false. ! event
                  nclstr1=0
                  totCLlength=0
                  flag_shower = .true.
                  print*,'Event ',nev1,
     $                 ': more than ',nclstrmax,' clusters'
                  goto 2000
               endif
               view(nclstr1)=iv !vista del cluster
               ladder(nclstr1)=nld(iseed,iv) !ladder a cui appartiene il cluster seed
               maxs(nclstr1)=iseed !strip del cluster seed
               mult(nclstr1)=multtemp !molteplicita'
               
               indstart(nclstr1)=ind !posizione dell'inizio del cluster nell'
                                ! array clsignal
               indmax(nclstr1)=indstart(nclstr1)+(iseed-lmax) !posizione del
                                ! cluster seed nell'array clsignal
               
               CLlength=rmax-lmax+1 !numero di strip del cluster
               totCLlength=totCLlength+CLlength
               dedx(nclstr1)=0
               do j=lmax,rmax   !stores sequentially cluter strip values in
                  clsignal(ind)=value(j) ! clsignal array
                  ind=ind+1
c                  if(value(j).gt.0)
                  if(value(j).gt.clinclcut(j))
     $                 dedx(nclstr1)=dedx(nclstr1)+value(j) !cluster charge
               enddo
c--------------------------------------------------------
c     
c--------------------------------------------------------
            endif               !end possible seed conditio
 220        continue            !jumps here to skip strips left of last seed

         enddo                  ! end loop on strips
      enddo                     !end loop on ladders
 2000 continue
      return
      end


*---***---***---***---***---***---***---***---***
*
*
*
*
*
*---***---***---***---***---***---***---***---***


      subroutine stripmask(data_file)

*     this routine set va1 and single-strip masks,
*     on the basis of date and id downlink
*
*     mask(nviews,nva1_view,nstrips_va1) !strip mask
*     mask_vk(nviews,nva1_view)          !VA1 mask
*
      include '../common/commontracker.f'
      include '../common/level1.f'
      include '../common/calib.f'

      character*10 data_file

      character*3 aid
      character*6 adate
      integer id
      integer date

*     ----------------------
*     retrieve date and id
      aid=data_file(8:10)
      adate=data_file(1:6)
      READ (aid, '(I3)'), id
      READ (adate, '(I6)'), date
*     ----------------------

*     init mask
      do iv=1,nviews
         do ivk=1,nva1_view
            mask_vk(iv,ivk)=1
            do is=1,nstrips_va1
               mask(iv,ivk,is)=1
            enddo
         enddo
      enddo

*     ---------------------
*     VIEW 2 - VK 23-24
*     couple of vk damaged during integration
      if(date.ge.50208)then
         print*,'MASK: view 2 - vk 23/24'
         mask_vk(2,23)=0
         mask_vk(2,24)=0
         do is=1,nstrips_va1
            mask(2,23,is)=0
            mask(2,24,is)=0
         enddo
      endif

*     ---------------------
*     VIEW 7 - VK 11-12
      if(date.ge.50209)then
        if(.not.(date.eq.50209.and.id.le.6)) then
          print*,'MASK: view 7 - vk 11/12'
          mask_vk(7,11)=0
          mask_vk(7,12)=0
          do is=1,nstrips_va1
            mask(7,11,is)=0
            mask(7,12,is)=0
          enddo
        endif
      endif

*     ---------------------
*     VIEW 7 - VK 21-22
      if(date.ge.50316)then
         print*,'MASK: view 7 - vk 21/22'
         mask_vk(7,21)=0
         mask_vk(7,22)=0
         do is=1,nstrips_va1
            mask(7,21,is)=0
            mask(7,22,is)=0
         enddo
      endif

*     ---------------------
*     VIEW 12 - VK 1-2-3-4
      if((date.eq.50317).and.(id.le.3))then
        print*,'MASK: view 12 - vk 1/2/3/4'
        mask_vk(12,1)=0
        mask_vk(12,2)=0
        mask_vk(12,3)=0
        mask_vk(12,4)=0
        do is=1,nstrips_va1
          mask(12,1,is)=0
          mask(12,2,is)=0
          mask(12,3,is)=0
          mask(12,4,is)=0
        enddo
      endif

*     ---------------------
*     VIEW 7 - VK 5-6
      if(date.ge.50320)then
        if(.not.(date.eq.50320.and.id.le.3)) then
          print*,'MASK: view 7 - vk 5/6'
          mask_vk(7,5)=0
          mask_vk(7,6)=0
          do is=1,nstrips_va1
            mask(7,5,is)=0
            mask(7,6,is)=0
          enddo
        endif
      endif

*     ---------------------
*     VIEW 7 - VK 13-14
      if(date.ge.50320)then
        if(.not.(date.eq.50320.and.id.le.5)) then
          print*,'MASK: view 7 - vk 13/14'
          mask_vk(7,13)=0
          mask_vk(7,14)=0
          do is=1,nstrips_va1
            mask(7,13,is)=0
            mask(7,14,is)=0
          enddo
        endif
      endif

***   SAMARA
***   SAMARA
***   SAMARA
*     it needs further checks...

*     ---------------------
*     VIEW 7 - VK 9-10
*     VIEW 12 - VK 1-2-3-4
      if((date.eq.50516).and.(id.le.8))then
        print*,'MASK: view 7 - vk 9/10'
        print*,'MASK: view 12 - vk 1/2/3/4'
        mask_vk(7,9)=0
        mask_vk(7,10)=0
        mask_vk(12,1)=0
        mask_vk(12,2)=0
        mask_vk(12,3)=0
        mask_vk(12,4)=0
        do is=1,nstrips_va1
          mask(7,9,is)=0
          mask(7,10,is)=0
          mask(12,1,is)=0
          mask(12,2,is)=0
          mask(12,3,is)=0
          mask(12,4,is)=0
        enddo
      endif

*     ---------------------
*     VIEW 7 - VK 9-10
      if(date.ge.50516)then
        if(.not.(date.eq.50516.and.id.le.8)) then
          print*,'MASK: view 7 - vk 9/10'
          mask_vk(7,9)=0
          mask_vk(7,10)=0
          do is=1,nstrips_va1
            mask(7,9,is)=0
            mask(7,10,is)=0
          enddo
        endif
      endif

*     ---------------------
*     VIEW 12 - VK 7-8
      if(date.ge.50523)then
        if(.not.(date.eq.50523.and.id.le.3)) then
          print*,'MASK: view 12 - vk 7/8'
          mask_vk(12,7)=0
          mask_vk(12,8)=0
          do is=1,nstrips_va1
            mask(12,7,is)=0
            mask(12,8,is)=0
          enddo
        endif
      endif

      return
      end 
