source/analysis/momanhough-efficiency.F

*************************************************************************
*     Program analysis.f
*
*     - reads cluster information (LEVEL1, reduction.f output ntuple)
*     - perform track identification and fit
*     - create LEVEL2 ntuple
*
*************************************************************************
      program momanhough

      include '../common/commontracker.f'
      include '../common/common_momanhough.f'
      include '../common/momanhough_init.f'
      include '../common/common_level2debug.f'
      include '../common/common_mech.f'
      include '../common/common_xyzPAM.f'
      include '../common/common_mini_2.f'
      include '../common/calib.f'
      include '../common/level1.f'
      include '../common/level2.f'

*     flag to set debug mode
      logical DEBUG
      common/dbg/DEBUG

*     flag to chose PFA
      character*10 PFA
      common/FINALPFA/PFA


c      parameter (inf=1.e8)      !just a huge number...
        
*     external functions
      external npl
      external acoordsi,coordsi,nld,coord,dcoord
c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------

      character*24 processing_date

      parameter (lun_data_level1=71) !data file id number
      parameter (lun_data_level2=72) !data file id number
      parameter (lun_data_calib=74) !data file id number

      character*74 data_file    !data file name
      character*74 data_dir    !data file name
      character*74 data_file_calib
      character*74 data_file_level1
      character*74 data_file_level2
# ifndef TEST2003
      parameter(ncalibmax=50)
      character*40 file_calib(ncalibmax)
      parameter(lun_calib_list=66) !calibration list file id
      integer which_calib_last
# endif

      integer minevent          !first event to be analysed
      logical FIMAGE            !


*     ******************************************
*     PLANE TO EXCLUDE (for efficiency studies)
*     --> provided as input parameter
*     (set according to "tracking" notation): 
*            1-6 from TOP to BOTTOM
*     ******************************************
      integer ipexclude
      common/plane_exclude/ipexclude


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

c------------------------------------------------------------------------
c     
c     HBOOK initialization
c     
c------------------------------------------------------------------------

      call HLIMIT(NWPAWC)


c------------------------------------------------------------------------
c     
c     reads input informations
c     
c------------------------------------------------------------------------
      call fdate(processing_date)
      write(*,101)
     $     processing_date
 101  format(/
     $      ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $      ,'*                                                 *',/
     $      ,'* ANALYSIS                                        *',/     
     $      ,'*                                                 *',/
     $      ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $     ,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
      minevent=1
      print*,'Maximum number of events to be analized:' !20000
      read(*,*) ntotev
      print*,ntotev    
      print*,'Position-finding algorythm: (GOG2,ETA2)'
      read(*,*)PFA
      PFA=PFA(1:lnblnk(PFA))
      print*,PFA
      print*,'DEBUG mode: (T, F)'
 11   format(l1)
      read(*,11)DEBUG
      print*,DEBUG
      print*,'Plane to exclude from fitting (1-6 from TOP to BOTTOM):'
      read(*,*) ipexclude
      print*, ipexclude    
      
      print*,'---------------------------------------------------'

****** INITIALIZATIONS *************************************

*     1) read charge-correlation parameters
      print*,' '
      print*,'- read charge-correlation parameters'
      print*,' '
      call readchargeparam

*     2) read z coordinates of the planes
      print*,' '
      print*,'- read z coordinates of the planes'
      print*,' '
      call mech_sensor          !reads sensors centres coordinates
      do ip=1,nplanes
        fitz(ip)=z_mech_sensor(ip,1,1)*0.1 !cm
*     gets planes mechanical z positions
*     (in mm) and sets them in micrometers
      enddo                     


*     3) read eta PFA parameters
      print*,' '
      print*,'- read PFA parameters'
      print*,' '
      call readetaparam

*     4) read magnetic field map      
      print*,' '
      print*,'- read magnetic field map'
      print*,' '
c      call read_B(5)            !legge il nome da STI!!.. temporaneo
c      call read_B_2maps(5)            !legge il nome da STI!!.. temporaneo
      call read_B

*     5) read allignment parameters
      print*,' '
      print*,'- read aligment parameters'
      print*,' '
      call readalignparam
 
************************************************************


c------------------------------------------------------------------------
c     
c     LEVEL 2 ntuple booking
c     
c------------------------------------------------------------------------


 503  format(a,'DW_',a,'_level2_no-p',i1,'.rz')
      write(data_file_level2,503)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
     $     ,ipexclude

      print*,''
      print*,'------------------------------------'
      print*,'  Creating LEVEL2 rz file'
      print*,'  ',data_file_level2
      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_level2,
     $     'LEVEL2',data_file_level2,'QNP',4096,istat) !opens rz  
      if(istat.ne.0) goto 19      
      call book_level2

      if(DEBUG)then
         print*,'(creating DEBUG nt-uple and histos)'
         call book_debug
      endif
      

c------------------------------------------------------------------------
c     
c     fills bad variables with online DSP bad strips from calib histograms
c     for bad strip exclusion
c     
c------------------------------------------------------------------------
        
# ifdef TEST2003
c------------------------------------------------------------------------
c     
c     opens calib file
c     
c------------------------------------------------------------------------

 505  format(a,'output_',a,'_calib.rz')
      write(data_file_calib,505)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
         
      print*,' '
      print*,'OPENING CALIB FILE...', data_file_calib
      
      IQUEST(10)=65000
      call HROPEN(lun_data_calib,'IN',data_file_calib,'QP',4096,istat)
      if(istat.ne.0) goto 17
      call HCDIR('//IN',' ')
      
      call HRIN(0,9999,0)       !puts histograms in memory
      
      print*,' '

      print*,' '
      print*,'READING PEDESTAL, SIGMA AND BADSTRIP HISTOGRAMS...'
      print*,' '
      
      call fillpedsig
      
      call HREND('IN')
      close (lun_data_calib)  
# else

      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

# endif 
      
c------------------------------------------------------------------------
c     
c     opens level1 file
c     
c------------------------------------------------------------------------

# ifdef TEST2003
 504  format(a,'output_',a,'_level1.rz')
# else
 504  format(a,'DW_',a,'_level1.rz')
# endif
      write(data_file_level1,504)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
      print*,''
      print*,'OPENING LEVEL1 FILE:'
      print*,data_file_level1
      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,'QP',4096,istat) !opens rz  
      if(istat.ne.0) goto 19      
      call HRIN(ntp_level1,9999,20)       
      call access_level1
c      call HPRNTU(ntp_level1+20)  
      call HNOENT(ntp_level1+20,iemax0)
      print*,'  events',iemax0


************************************************************
************************************************************
************************************************************
*
*     start track analysis
*
************************************************************
************************************************************
************************************************************
      maxevent=minevent+ntotev-1
      do iev = minevent,MIN(iemax0,maxevent) !loop on events

         call HCDIR('//LEVEL1',' ')
         call HGNT(ntp_level1+20,iev,ierr) !reads an event
         if(ierr.ne.0) goto 21
*------------------------------------------------------
*        LEVEL2 N-TUPLE INITIALIZATIONS
         call init_level2
         if(.not.good1)then
            goto 8800           !fill nt-uple and go to next event
         endif
*------------------------------------------------------


# ifndef TEST2003
      
         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 ',nev2
     $           ,'   (CPU pkt N.',pkt_num1,')'
            print*,'--> ',data_file_calib            
            IQUEST(10)=65000          
            call HROPEN(lun_data_calib,
     $           '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_calib)            
            which_calib_last=which_calib           
            
         elseif(which_calib.eq.0)then
            
            nocalib=nocalib+1
            good2=.false.
            goto 8800 !fill nt-uple and go to next event
            
         endif

# endif

      
*------------------------------------------------------
*     cut on maximum number of clusters
*------------------------------------------------------
         if(nclstr1.gt.nclstrmax_level2)then
            goto 8800           !fill nt-uple and go to next event
         endif

         do i=1,nclstr1
            cl_used(i)=0        !init mask of clusters associated to a track
         enddo

         if(DEBUG)then
            print*,'----------------------------------'
            print*,iev,'   ** ',nev2
         endif

*-------------------------------------------------------------------------------
*     STEP 1
*-------------------------------------------------------------------------------
*     X-Y cluster association
*     
*     Clusters are associated to form COUPLES
*     Clusters not associated in any couple are called SINGLETS
*
*     Track identification (Hough transform) and fitting is first done on couples. 
*     Hence singlets are possibly added to the track.
*    
*     Variables assigned by the routine "cl_to_couples" are those in the 
*     common blocks:
*     - common/clusters/cl_good
*     - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2
*     - common/singlets/ncls,cls,cl_single
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

         iflag=0
         call cl_to_couples(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*     selezione di tracce pulite per diagnostica
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c$$$         if(DEBUG)then
c$$$            do ip=1,nplanes 
c$$$               if(ncp_plane(ip).ne.1)good2=.false.
c$$$            enddo 
c$$$c            if(good2.eq.0)goto 100!next event 
c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event 
c$$$         endif
         
 
*-----------------------------------------------------
*-----------------------------------------------------
*     HOUGH TRASFORM
*-----------------------------------------------------
*-----------------------------------------------------


*-------------------------------------------------------------------------------
*     STEP 2
*-------------------------------------------------------------------------------
*     
*     Association of couples to form
*     - DOUBLETS in YZ view
*     - TRIPLETS in XZ view
*    
*     Variables assigned by the routine "cp_to_doubtrip" are those in the 
*     common blocks:
*     - common/hough_param/
*     $     alfayz1,  !Y0
*     $     alfayz2,  !tg theta-yz
*     $     alfaxz1,  !X0
*     $     alfaxz2,  !tg theta-xz
*     $     alfaxz3   !1/r
*     - common/doublets/ndblt,cpyz1,cpyz2
*     - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

         iflag=0
         call cp_to_doubtrip(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif


*-------------------------------------------------------------------------------
*     STEP 3
*-------------------------------------------------------------------------------
*     
*     Classification of doublets and triplets to form CLOUDS,  
*     according to distance in parameter space. 
*     
*     cloud = cluster of points (doublets/triplets) in parameter space
*
*     
*    
*     Variables assigned by the routine "doub_to_YZcloud" are those in the 
*     common blocks:
*     - common/clouds_yz/                   
*     $     nclouds_yz                        
*     $     ,alfayz1_av,alfayz2_av           
*     $     ,ptcloud_yz,db_cloud,cpcloud_yz
*
*     Variables assigned by the routine "trip_to_XZcloud" are those in the 
*     common blocks:
*      common/clouds_xz/                   
*     $      nclouds_xz                        
*     $     ,alfaxz1_av,alfaxz2_av,alfaxz3_av 
*     $     ,ptcloud_xz,tr_cloud,cpcloud_xz           
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

         iflag=0
         call doub_to_YZcloud(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif
         iflag=0
         call trip_to_XZcloud(iflag)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif


*-------------------------------------------------------------------------------
*     STEP 4   (ITERATED until any other physical track isn't found)
*-------------------------------------------------------------------------------
*     
*     YZ and XZ clouds are combined in order to obtain the initial guess
*     of the candidate-track parameters.
*     A minimum number of matching couples between YZ and XZ clouds is required.
*
*     A TRACK CANDIDATE is defined by 
*     - the couples resulting from the INTERSECTION of the two clouds, and 
*     - the associated track parameters (evaluated by performing a zero-order  
*       track fitting) 
*
*     The NTRACKS candidate-track parameters are stored in common block:
*    
*     - common/track_candidates/NTRACKS,AL_STORE
*     $     ,XV_STORE,YV_STORE,ZV_STORE
*     $     ,XM_STORE,YM_STORE,ZM_STORE
*     $     ,RESX_STORE,RESY_STORE
*     $     ,AXV_STORE,AYV_STORE
*     $     ,XGOOD_STORE,YGOOD_STORE
*     $     ,CP_STORE,RCHI2_STORE
*
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
         ntrk=0                 !counter of identified physical tracks

11111    continue               !<<<<<<< come here when performing a new search

         iflag=0
         call clouds_to_ctrack(iflag)
         if(iflag.eq.1)then     !no candidate tracks found
            goto 880            !fill ntp and go to next event   
         endif

         FIMAGE=.false.         !processing best track (not track image)
         ibest=0                !best track among candidates
         iimage=0               !track image
*     ------------- select the best track -------------
         rchi2best=100000.
         do i=1,ntracks
            if(RCHI2_STORE(i).lt.rchi2best.and.
     $         RCHI2_STORE(i).gt.0)then
               ibest=i
               rchi2best=RCHI2_STORE(i)
            endif
         enddo
         if(ibest.eq.0)goto 880 !>> no good candidates 
*-------------------------------------------------------------------------------     
*     The best track candidate (ibest) is selected and a new fitting is performed. 
*     Previous to this, the track is refined by:
*     - possibly adding new COUPLES or SINGLETS from the missing planes
*     - evaluating the coordinates with improved PFAs 
*       ( angle-dependent ETA algorithms )
*-------------------------------------------------------------------------------

 1212    continue               !<<<<< come here to fit track-image

         if(.not.FIMAGE)then    !processing best track 
            icand=ibest            
         else                   !processing image
            icand=iimage
            iimage=0
         endif
         if(icand.eq.0)then
            print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
     $           ,ibest,iimage
            return
         endif
         call refine_track(icand)


*     ***************************************
*     Checks if a plane should be excluded 
*     (for efficiency study)
*     ***************************************
         xgood(ipexclude) = 0 !<<<
         ygood(ipexclude) = 0 !<<<
*     ***************************************
         

*     **********************************************************
*     ************************** FIT *** FIT *** FIT *** FIT ***
*     **********************************************************
         do i=1,5
            AL(i)=dble(AL_STORE(i,icand))
         enddo
         ifail=0                !error flag in chi2 computation
         jstep=0                !# minimization steps

         call mini_2(jstep,ifail)
         if(ifail.ne.0) then
            if(DEBUG)then
               print *,
     $              '*** MINIMIZATION FAILURE *** (mini_2) '
     $              ,iev
            endif
            chi2=-chi2 
         endif 
         
         if(DEBUG)then
            print*,'----------------------------- improved track coord'
22222       format(i2,' * ',3f9.4,' --- ',4f9.4,' --- ',2f4.0,2f8.5)
            do ip=1,6
               write(*,22222)ip,zm(ip),xm(ip),ym(ip)
     $              ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
     $              ,xgood(ip),ygood(ip),resx(ip),resy(ip)
            enddo
         endif

c         rchi2=chi2/dble(ndof)
         if(DEBUG)then
            print*,' ' 
            print*,'****** SELECTED TRACK *************'
            print*,'#         R. chi2        RIG'
            print*,' --- ',chi2,' --- '
     $           ,1./abs(AL(5)) 
            print*,'***********************************'
         endif
*     **********************************************************
*     ************************** FIT *** FIT *** FIT *** FIT ***
*     **********************************************************


*     ------------- search if the track has an IMAGE -------------
*     ------------- (also this is stored )           -------------
         if(FIMAGE)goto 122     !>>> jump! (this is already an image)
*     now search for track-image, by comparing couples IDs
         do i=1,ntracks
            iimage=i
            do ip=1,nplanes
               if(     CP_STORE(nplanes-ip+1,icand).ne.
     $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0
            enddo
            if(  iimage.ne.0.and.
c     $           RCHI2_STORE(i).le.CHI2MAX.and.
c     $           RCHI2_STORE(i).gt.0.and.
     $           .true.)then 
               if(DEBUG)print*,'Track candidate ',iimage
     $              ,' >>> TRACK IMAGE >>> of'
     $              ,ibest
               goto 122         !image track found
            endif
         enddo
 122     continue

*     --- and store the results --------------------------------
         ntrk = ntrk + 1                   !counter of found tracks
         if(.not.FIMAGE
     $        .and.iimage.eq.0) image(ntrk)= 0
         if(.not.FIMAGE
     $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
         if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
         call fill_level2_tracks(ntrk)     !==> good2=.true.

c         print*,'++++++++++ iimage,fimage,ntrk,image '
c     $        ,iimage,fimage,ntrk,image(ntrk)

         if(ntrk.eq.NTRKMAX)then
            if(DEBUG)
     $           print*,
     $           '** warning ** number of identified '// 
     $           'tracks exceeds vector dimension '
     $           ,'( ',NTRKMAX,' )'
cc            good2=.false.
            goto 880            !fill ntp and go to next event
         endif
         if(iimage.ne.0)then
            FIMAGE=.true.       !
            goto 1212           !>>> fit image-track
         endif

*     --- then remove selected clusters (ibest+iimage) from clouds ----
         call clean_XYclouds(ibest)
         if(iflag.eq.1)then     !bad event
            goto 880            !fill ntp and go to next event            
         endif

*     **********************************************************
*     condition to start a new search
*     **********************************************************
         ixznew=0
         do ixz=1,nclouds_xz
            if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1
         enddo
         iyznew=0
         do iyz=1,nclouds_yz
            if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1
         enddo
         
         if(ixznew.ne.0.and.
     $      iyznew.ne.0.and.
     $        rchi2best.le.CHI2MAX.and.
c     $        rchi2best.lt.15..and.
     $        .true.)then
            if(DEBUG)then
               print*,'***** NEW SEARCH ****'
            endif
            goto 11111          !try new search
            
         endif
*     **********************************************

*     >>>>>>>>>>>>>>>>>>> NT-UPLE filling <<<<<<<<<<<<<<<<<<<<
*       +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +
*      / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ /
*     +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   

 880     continue

*     count #cluster per plane not associated to any track
         do icl=1,nclstr1
            if(cl_used(icl).eq.0)then !cluster not included in any track
               ip=nplanes-npl(VIEW(icl))+1
               if(mod(VIEW(icl),2).eq.0)nclsx(ip)=nclsx(ip)+1
               if(mod(VIEW(icl),2).eq.1)nclsy(ip)=nclsy(ip)+1
            endif
c            print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
         enddo
         if(DEBUG)then

            print*,''
            print*,'DONE!'
            print*,''
            print*,'* summary *'
            print*,'tracks    ',ntrk
            print*,'cl used   ',(cl_used(i),i=1,nclstr1)
            print*,'cl unused (x-y)'
            do ip=1,nplanes
               print*,ip,' << ',nclsx(ip),nclsy(ip)
            enddo
            print*,''
            print*,''
         endif


 8800    continue

         call HCDIR('//LEVEL2',' ')

         call HFNT(ntp_level2)  !fill LEVEL2 nt-uple         

         if(DEBUG)then
            good2_nt=good2
            nev2_nt=nev2
            
            if(ntrpt.gt.ntrpt_max_nt.or.ndblt.gt.ndblt_max_nt)then
               good2_nt=.false.
               ntrpt_nt=0
               ndblt_nt=0
               NCLOUDS_XZ=0
               NCLOUDS_YZ=0
            else
               ndblt_nt=ndblt
               ntrpt_nt=ntrpt
               do id=1,ndblt
                  alfayz1_nt(id)=alfayz1(id) !Y0
                  alfayz2_nt(id)=alfayz2(id) !tg theta-yz
                  db_cloud_nt(id)=db_cloud(id)
               enddo
               do it=1,ntrpt
                  alfaxz1_nt(it)=alfaxz1(it) !X0
                  alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
                  alfaxz3_nt(it)=alfaxz3(it) !1/r
                  tr_cloud_nt(it)=tr_cloud(it)
               enddo
               nclouds_yz_nt=nclouds_yz
               nclouds_xz_nt=nclouds_xz
               do iyz=1,nclouds_yz
                  ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
                  alfayz1_av_nt(iyz)=alfayz1_av(iyz)
                  alfayz2_av_nt(iyz)=alfayz2_av(iyz)
               enddo
               do ixz=1,nclouds_xz
                  ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
                  alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
                  alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
                  alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
               enddo
            endif
c            print*,ntrpt,ntrpt_max_nt,ndblt,ndblt_max_nt
c            do iii=1,ndblt
c               print*,iii,alfayz1(iii),alfayz2(iii)
c            enddo
c            print*,good2_nt,good2
            call HFNT(ntp_level2+1) !fill DEBUG nt-uple
         endif


 100     continue
      enddo                     !end loop on events

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

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


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

c------------------------------------------------------------------------
c     
c     level1 ntuple event reading error
c     
c------------------------------------------------------------------------
      
 21   continue
      
      print*,' '
      print*,'ERROR WHILE READING LEVEL1 NTUPLE, AT EVENT
     $     : ',iev
      print*,' '
      print*,' '
      
      goto 9000                 !the end

c------------------------------------------------------------------------
c     
c     calib file opening error
c     
c------------------------------------------------------------------------
      
 17   continue
      
      print*,' '
      print*,'preanalysis: ERROR OPENING INPUT FILE: ',data_file_calib
      print*,' '
      print*,' '
      
      goto 9000                 !the end

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

 9000 continue
      
      if(DEBUG)call HPRNTU(ntp_level2+1)
      call HPRNTU(ntp_level2)
      print*,''
      call HCDIR('//LEVEL2',' ')
      call HROUT(ntp_level2,ICYCLE,'T')
      print *,'- Stored LEVEL2 nt-uple (',nev2,' entries )'
      if(DEBUG)call HROUT(ntp_level2+1,ICYCLE,'T')
      if(DEBUG)print *,'- Stored DEBUG nt-uple '
      call HREND('level2')
      close(lun_data_level2)

      call HCDIR('//LEVEL1',' ')
      call HREND('level1')
      close(lun_data_level1)

      stop
      end


# include "momanhough-subroutines.F"


c$$$
c$$$
c$$$
c$$$************************************************************
c$$$
c$$$      subroutine readchargeparam
c$$$      
c$$$      
c$$$      include '../common/commontracker.f'
c$$$      include '../common/calib.f'
c$$$
c$$$c      character*40 fname_binning
c$$$      character*60 fname_param 
c$$$c      character*120 cmd1
c$$$c      character*120 cmd2
c$$$
c$$$c 201  format('../common/charge_param/charge-l',i1,'.dat')      
c$$$c 201  format('$TRK_GRND/source/common/charge_param/charge-l',i1,'.dat')
c$$$ 201  format('charge-l',i1,'.dat')
c$$$      do ilad=1,nladders_view         
c$$$         write(fname_param,201)ilad
c$$$c$$$         cmd1='cp $TRK_GRND/source/common/charge_param/'
c$$$c$$$     $        //fname_param(1:LNBLNK(fname_param))//' .'         
c$$$c$$$         call system(cmd1)
c$$$         print *,'Opening file: ',fname_param
c$$$         open(10,
c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))
c$$$     $        ,STATUS='UNKNOWN'
c$$$     $        ,IOSTAT=iostat
c$$$     $        )
c$$$         if(iostat.ne.0)then
c$$$            print*,'READCHARGEPARAM: *** Error in opening file ***'
c$$$            return
c$$$         endif
c$$$         do ip=1,nplanes
c$$$            read(10,*
c$$$     $           ,IOSTAT=iostat
c$$$     $           )pip,
c$$$     $            kch(ip,ilad),cch(ip,ilad),sch(ip,ilad)         
c$$$c            print*,ilad,ip,pip,kch(ip,ilad),
c$$$c     $           cch(ip,ilad),sch(ip,ilad) 
c$$$         enddo
c$$$         close(10)
c$$$c$$$         cmd2='rm -f '
c$$$c$$$     $        //fname_param(1:LNBLNK(fname_param))
c$$$c$$$         call system(cmd2)
c$$$      enddo
c$$$
c$$$      return
c$$$      end
c$$$      
c$$$************************************************************
c$$$************************************************************
c$$$************************************************************
c$$$************************************************************
c$$$*
c$$$*     This routine provides the coordinates (in cm) in the PAMELA reference system:
c$$$*       - of the point associated with a COUPLE ---> (xPAM,yPAM,zPAM)
c$$$*       - of the extremes of the segment
c$$$*         associated with a SINGLET ---------------> (xPAM_A,yPAM_A,zPAM_A) 
c$$$*                                               ---> (xPAM_B,yPAM_B,zPAM_B)
c$$$*
c$$$*     It also assigns the spatial resolution to the evaluated coordinates,
c$$$*     as a function (in principle) of the multiplicity, the angle, the PFA etc...
c$$$*
c$$$*     
c$$$*     To call the routine you must pass the arguments:
c$$$*     icx    - ID of cluster x
c$$$*     icy    - ID of cluster y
c$$$*     sensor - sensor (1,2)
c$$$*     PFAx   - Position Finding Algorithm in x (COG2,ETA2,...)
c$$$*     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
c$$$*     angx   - Projected angle in x
c$$$*     angy   - Projected angle in y
c$$$*
c$$$*     --------- COUPLES -------------------------------------------------------
c$$$*     The couple defines a point in the space. 
c$$$*     The coordinates of the point are evaluated as follows:
c$$$*     1 - the corrected coordinates relative to the sensor are evaluated
c$$$*         according to the chosen PFA --> (xi,yi,0)
c$$$*     2 - coordinates are rotated and traslated, according to the aligmnet 
c$$$*         parameters, and expressed in the reference system of the mechanical 
c$$$*         sensor --> (xrt,yrt,zrt)
c$$$*     3 - coordinates are finally converted to the PAMELA reference system 
c$$$*         --> (xPAM,yPAM,zPAM)
c$$$*
c$$$*     --------- SINGLETS -------------------------------------------------------
c$$$*     Since a coordinate is missing, the singlet defines not a point 
c$$$*     in the space but a segment AB (parallel to the strips). 
c$$$*     In this case the routine returns the coordinates in the PAMELA reference
c$$$*     system of the two extremes A and B of the segment:
c$$$*         --> (xPAM_A,yPAM_A,zPAM_A)
c$$$*         --> (xPAM_B,yPAM_B,zPAM_B)
c$$$*     
c$$$*     ==========================================================
c$$$*
c$$$*     The output of the routine is stored in the commons:
c$$$*
c$$$*      double precision xPAM,yPAM,zPAM
c$$$*      common/coord_xyz_PAM/xPAM,yPAM,zPAM
c$$$*     
c$$$*      double precision xPAM_A,yPAM_A,zPAM_A
c$$$*      double precision xPAM_B,yPAM_B,zPAM_B
c$$$*      common/coord_AB_PAM/xPAM_A,yPAM_A,zPAM_A,xPAM_B,yPAM_B,zPAM_B
c$$$*
c$$$*      double precision resxPAM,resyPAM
c$$$*      common/resolution_PAM/resxPAM,resyPAM
c$$$*
c$$$*     (in file ../common/common_xyzPAM.f)
c$$$*
c$$$*
c$$$
c$$$      subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)
c$$$      
c$$$      include '../common/commontracker.f'
c$$$      include '../common/calib.f'
c$$$      include '../common/level1.f'
c$$$      include '../common/common_align.f'
c$$$      include '../common/common_mech.f'
c$$$      include '../common/common_xyzPAM.f'
c$$$      include '../common/common_resxy.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$      integer icx,icy           !X-Y cluster ID
c$$$      integer sensor
c$$$      integer viewx,viewy
c$$$      character*4 PFAx,PFAy     !PFA to be used
c$$$      real angx,angy            !X-Y angle
c$$$
c$$$      real stripx,stripy
c$$$
c$$$      double precision xrt,yrt,zrt
c$$$      double precision xrt_A,yrt_A,zrt_A
c$$$      double precision xrt_B,yrt_B,zrt_B
c$$$c      double precision xi,yi,zi
c$$$c      double precision xi_A,yi_A,zi_A
c$$$c      double precision xi_B,yi_B,zi_B
c$$$      
c$$$
c$$$      parameter (ndivx=30)
c$$$      
c$$$
c$$$      resxPAM = 0
c$$$      resyPAM = 0
c$$$
c$$$      xPAM = 0.
c$$$      yPAM = 0.
c$$$      zPAM = 0.
c$$$      xPAM_A = 0.
c$$$      yPAM_A = 0.
c$$$      zPAM_A = 0.
c$$$      xPAM_B = 0.
c$$$      yPAM_B = 0.
c$$$      zPAM_B = 0.
c$$$
c$$$*     -----------------
c$$$*     CLUSTER X
c$$$*     -----------------
c$$$
c$$$      if(icx.ne.0)then
c$$$         viewx = VIEW(icx)
c$$$         nldx = nld(MAXS(icx),VIEW(icx))
c$$$         nplx = npl(VIEW(icx))
c$$$         resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
c$$$         
c$$$         stripx = float(MAXS(icx))
c$$$         if(PFAx.eq.'COG2')then
c$$$            stripx = stripx + cog(2,icx)            
c$$$            resxPAM = resxPAM*fbad_cog(2,icx)
c$$$         elseif(PFAx.eq.'ETA2')then
c$$$            cog2 = cog(2,icx)
c$$$            etacorr = pfa_eta2(cog2,viewx,nldx,angx)
c$$$c            print*,'---- ',viewx,cog2,etacorr
c$$$            stripx = stripx + etacorr
c$$$c     print*,'COG2 ',cog2,' >>> ETA2 ',etacorr
c$$$            if(DEBUG.and.fbad_cog(2,icx).ne.1)
c$$$     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
c$$$            resxPAM = resxPAM*fbad_cog(2,icx)
c$$$         else
c$$$            print*,'Accontentati di MAXS !!!!'
c$$$         endif
c$$$
c$$$
c$$$      endif
c$$$      
c$$$*     ----------------- 
c$$$*     CLUSTER Y
c$$$*     -----------------
c$$$
c$$$      if(icy.ne.0)then
c$$$         viewy = VIEW(icy)
c$$$         nldy = nld(MAXS(icy),VIEW(icy))
c$$$         nply = npl(VIEW(icy))
c$$$         resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
c$$$
c$$$
c$$$         if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
c$$$            print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
c$$$     $           ,icx,icy
c$$$            goto 100
c$$$         endif
c$$$         
c$$$         stripy = float(MAXS(icy))
c$$$         if(PFAy.eq.'COG2')then
c$$$            stripy = stripy + cog(2,icy)
c$$$            resyPAM = resyPAM*fbad_cog(2,icy)
c$$$         elseif(PFAy.eq.'ETA2')then
c$$$            cog2 = cog(2,icy)
c$$$            etacorr = pfa_eta2(cog2,viewy,nldy,angy)
c$$$            stripy = stripy + etacorr
c$$$c     print*,'COG2 ',cog2,' >>> ETA2 ',etacorr
c$$$            resyPAM = resyPAM*fbad_cog(2,icy)
c$$$            if(DEBUG.and.fbad_cog(2,icy).ne.1)
c$$$     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
c$$$         else
c$$$            print*,'Accontentati di MAXS !!!!'
c$$$         endif
c$$$
c$$$      endif
c$$$
c$$$      
c$$$c===========================================================
c$$$C     COUPLE
c$$$C===========================================================
c$$$      if(icx.ne.0.and.icy.ne.0)then 
c$$$
c$$$c------------------------------------------------------------------------
c$$$c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
c$$$c------------------------------------------------------------------------
c$$$         xi = acoordsi(stripx,viewx)
c$$$         yi = acoordsi(stripy,viewy)
c$$$         zi = 0.
c$$$         
c$$$
c$$$c------------------------------------------------------------------------
c$$$c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c$$$c------------------------------------------------------------------------
c$$$c     N.B. I convert angles from microradiants to radiants
c$$$
c$$$         xrt = xi
c$$$     $        - omega(nplx,nldx,sensor)*yi 
c$$$     $        + gamma(nplx,nldx,sensor)*zi
c$$$     $        + dx(nplx,nldx,sensor)
c$$$         
c$$$         yrt = omega(nplx,nldx,sensor)*xi
c$$$     $        + yi
c$$$     $        - beta(nplx,nldx,sensor)*zi
c$$$     $        + dy(nplx,nldx,sensor)
c$$$         
c$$$         zrt = -gamma(nplx,nldx,sensor)*xi
c$$$     $        + beta(nplx,nldx,sensor)*yi
c$$$     $        + zi
c$$$     $        + dz(nplx,nldx,sensor)
c$$$         
c$$$c      xrt = xi
c$$$c      yrt = yi
c$$$c      zrt = zi
c$$$    
c$$$c------------------------------------------------------------------------
c$$$c     (xPAM,yPAM,zPAM) = measured coordinates (in cm) 
c$$$c                        in PAMELA reference system
c$$$c------------------------------------------------------------------------
c$$$
c$$$         xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
c$$$         yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
c$$$         zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
c$$$
c$$$         xPAM_A = 0.
c$$$         yPAM_A = 0.
c$$$         zPAM_A = 0.
c$$$
c$$$         xPAM_B = 0.
c$$$         yPAM_B = 0.
c$$$         zPAM_B = 0.
c$$$
c$$$      elseif(
c$$$     $        (icx.ne.0.and.icy.eq.0).or.
c$$$     $        (icx.eq.0.and.icy.ne.0).or.
c$$$     $        .false.
c$$$     $        )then
c$$$
c$$$c------------------------------------------------------------------------
c$$$c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
c$$$c------------------------------------------------------------------------
c$$$
c$$$         if(icy.ne.0)then
c$$$c===========================================================
c$$$C     Y-SINGLET
c$$$C===========================================================
c$$$            nplx = nply
c$$$            nldx = nldy
c$$$            viewx = viewy + 1
c$$$
c$$$            yi   = acoordsi(stripy,viewy)
c$$$
c$$$            xi_A = edgeY_d - SiDimX/2
c$$$            yi_A = yi
c$$$            zi_A = 0.
c$$$            
c$$$            xi_B = SiDimX/2 - edgeY_u
c$$$            yi_B = yi
c$$$            zi_B = 0.
c$$$
c$$$c            print*,'Y-cl ',icy,stripy,' --> ',yi
c$$$c            print*,xi_A,' <--> ',xi_B
c$$$            
c$$$         elseif(icx.ne.0)then
c$$$c===========================================================
c$$$C     X-SINGLET
c$$$C===========================================================
c$$$
c$$$            nply = nplx
c$$$            nldy = nldx
c$$$            viewy = viewx - 1
c$$$
c$$$            xi   = acoordsi(stripx,viewx)
c$$$
c$$$            xi_A = xi
c$$$            yi_A = edgeX_d - SiDimY/2
c$$$            zi_A = 0.
c$$$            
c$$$            xi_B = xi
c$$$            yi_B = SiDimY/2 - edgeX_u
c$$$            zi_B = 0.
c$$$            
c$$$            if(viewy.eq.11)then
c$$$               yi = yi_A
c$$$               yi_A = yi_B
c$$$               yi_B = yi
c$$$            endif
c$$$
c$$$c            print*,'X-cl ',icx,stripx,' --> ',xi
c$$$c            print*,yi_A,' <--> ',yi_B
c$$$
c$$$         else
c$$$
c$$$            print *,'routine xyz_PAM ---> not properly used !!!'
c$$$            print *,'icx = ',icx
c$$$            print *,'icy = ',icy
c$$$            goto 100
c$$$            
c$$$         endif
c$$$c------------------------------------------------------------------------
c$$$c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c$$$c------------------------------------------------------------------------
c$$$c     N.B. I convert angles from microradiants to radiants
c$$$
c$$$         xrt_A = xi_A
c$$$     $        - omega(nplx,nldx,sensor)*yi_A
c$$$     $        + gamma(nplx,nldx,sensor)*zi_A
c$$$     $        + dx(nplx,nldx,sensor)
c$$$         
c$$$         yrt_A = omega(nplx,nldx,sensor)*xi_A
c$$$     $        + yi_A
c$$$     $        - beta(nplx,nldx,sensor)*zi_A
c$$$     $        + dy(nplx,nldx,sensor)
c$$$         
c$$$         zrt_A = -gamma(nplx,nldx,sensor)*xi_A
c$$$     $        + beta(nplx,nldx,sensor)*yi_A
c$$$     $        + zi_A
c$$$     $        + dz(nplx,nldx,sensor)
c$$$
c$$$         xrt_B = xi_B
c$$$     $        - omega(nplx,nldx,sensor)*yi_B 
c$$$     $        + gamma(nplx,nldx,sensor)*zi_B
c$$$     $        + dx(nplx,nldx,sensor)
c$$$         
c$$$         yrt_B = omega(nplx,nldx,sensor)*xi_B
c$$$     $        + yi_B
c$$$     $        - beta(nplx,nldx,sensor)*zi_B
c$$$     $        + dy(nplx,nldx,sensor)
c$$$         
c$$$         zrt_B = -gamma(nplx,nldx,sensor)*xi_B
c$$$     $        + beta(nplx,nldx,sensor)*yi_B
c$$$     $        + zi_B
c$$$     $        + dz(nplx,nldx,sensor)
c$$$
c$$$         
c$$$c      xrt = xi
c$$$c      yrt = yi
c$$$c      zrt = zi
c$$$    
c$$$c------------------------------------------------------------------------
c$$$c     (xPAM,yPAM,zPAM) = measured coordinates (in cm) 
c$$$c                        in PAMELA reference system
c$$$c------------------------------------------------------------------------
c$$$
c$$$         xPAM = 0.
c$$$         yPAM = 0.
c$$$         zPAM = 0.
c$$$
c$$$         xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
c$$$         yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
c$$$         zPAM_A = ( zrt_A + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
c$$$
c$$$         xPAM_B = dcoord(xrt_B,viewx,nldx,sensor) / 1.d4
c$$$         yPAM_B = dcoord(yrt_B,viewy,nldy,sensor) / 1.d4
c$$$         zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
c$$$         
c$$$
c$$$c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
c$$$
c$$$      else
c$$$          
c$$$         print *,'routine xyz_PAM ---> not properly used !!!'
c$$$         print *,'icx = ',icx
c$$$         print *,'icy = ',icy
c$$$            
c$$$      endif
c$$$         
c$$$ 100  continue
c$$$      end
c$$$
c$$$
c$$$********************************************************************************
c$$$********************************************************************************
c$$$********************************************************************************
c$$$*
c$$$*     The function distance_to(XP,YP) should be used after 
c$$$*     a call to the xyz_PAM routine and it evaluate the 
c$$$*     NORMALIZED distance (PROJECTED on the XY plane) between 
c$$$*     the point  (XP,YP), argument of the function, 
c$$$*     and:
c$$$*     
c$$$*     - the point (xPAM,yPAM,zPAM), in the case of a COUPLE
c$$$*     or 
c$$$*     - the segment (xPAM_A,yPAM_A,zPAM_A)-(xPAM_B,yPAM_B,zPAM_B), 
c$$$*       in the case of a SINGLET.
c$$$*
c$$$*     ( The routine xyz_PAM fills the common defined in "common_xyzPAM.f",
c$$$*      which stores the coordinates of the couple/singlet )
c$$$*     
c$$$********************************************************************************
c$$$
c$$$      real function distance_to(XPP,YPP)
c$$$
c$$$      include '../common/common_xyzPAM.f'
c$$$
c$$$*     -----------------------------------
c$$$*     it computes the normalized distance 
c$$$*     ( i.e. distance/resolution )
c$$$*     -----------------------------------
c$$$
c$$$      double precision distance,RE
c$$$      double precision BETA,ALFA,xmi,ymi
c$$$
c$$$*     ----------------------
c$$$      if (
c$$$     +     xPAM.eq.0.and.
c$$$     +     yPAM.eq.0.and.
c$$$     +     zPAM.eq.0.and.
c$$$     +     xPAM_A.ne.0.and.
c$$$     +     yPAM_A.ne.0.and.
c$$$     +     zPAM_A.ne.0.and.
c$$$     +     xPAM_B.ne.0.and.
c$$$     +     yPAM_B.ne.0.and.
c$$$     +     zPAM_B.ne.0.and.
c$$$     +     .true.)then
c$$$*     -----------------------
c$$$*     DISTANCE TO --- SINGLET
c$$$*     -----------------------
c$$$         if(abs(sngl(xPAM_B-xPAM_A)).lt.abs(sngl(yPAM_B-yPAM_A)))then
c$$$*        |||---------- X CLUSTER
c$$$            
c$$$            BETA = (xPAM_B-xPAM_A)/(yPAM_B-yPAM_A)
c$$$            ALFA = xPAM_A - BETA * yPAM_A
c$$$
c$$$            ymi = ( YPP + BETA*XPP - BETA*ALFA )/(1+BETA**2)
c$$$            if(ymi.lt.dmin1(yPAM_A,yPAM_B))ymi=dmin1(yPAM_A,yPAM_B)
c$$$            if(ymi.gt.dmax1(yPAM_A,yPAM_B))ymi=dmax1(yPAM_A,yPAM_B)
c$$$            xmi = ALFA + BETA * ymi
c$$$            RE = resxPAM
c$$$                        
c$$$         else
c$$$*        |||---------- Y CLUSTER
c$$$
c$$$            BETA = (yPAM_B-yPAM_A)/(xPAM_B-xPAM_A)
c$$$            ALFA = yPAM_A - BETA * xPAM_A
c$$$
c$$$            xmi = ( XPP + BETA*YPP - BETA*ALFA )/(1+BETA**2)
c$$$            if(xmi.lt.dmin1(xPAM_A,xPAM_B))xmi=dmin1(xPAM_A,xPAM_B)
c$$$            if(xmi.gt.dmax1(xPAM_A,xPAM_B))xmi=dmax1(xPAM_A,xPAM_B)
c$$$            ymi = ALFA + BETA * xmi
c$$$            RE = resyPAM
c$$$            
c$$$         endif         
c$$$
c$$$         distance=
c$$$     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
c$$$         distance=dsqrt(distance)                     
c$$$
c$$$c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
c$$$c$$$     $        ,' --- distance_to --- ',xpp,ypp
c$$$c$$$         print*,' resolution ',re
c$$$
c$$$         
c$$$*     ----------------------
c$$$      elseif(
c$$$     +     xPAM.ne.0.and.
c$$$     +     yPAM.ne.0.and.
c$$$     +     zPAM.ne.0.and.
c$$$     +     xPAM_A.eq.0.and.
c$$$     +     yPAM_A.eq.0.and.
c$$$     +     zPAM_A.eq.0.and.
c$$$     +     xPAM_B.eq.0.and.
c$$$     +     yPAM_B.eq.0.and.
c$$$     +     zPAM_B.eq.0.and.
c$$$     +     .true.)then
c$$$*     ----------------------
c$$$*     DISTANCE TO --- COUPLE
c$$$*     ----------------------
c$$$         
c$$$         distance=
c$$$     $        ((xPAM-XPP)/resxPAM)**2
c$$$     $        +
c$$$     $        ((yPAM-YPP)/resyPAM)**2
c$$$         distance=dsqrt(distance)                     
c$$$
c$$$c$$$         print*,xPAM,yPAM,zPAM
c$$$c$$$     $        ,' --- distance_to --- ',xpp,ypp
c$$$c$$$         print*,' resolution ',resxPAM,resyPAM
c$$$         
c$$$      else
c$$$         
c$$$         print*
c$$$     $        ,' function distance_to ---> wrong usage!!!'
c$$$         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM 
c$$$         print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
c$$$     $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
c$$$      endif   
c$$$
c$$$      distance_to = sngl(distance)
c$$$
c$$$      return
c$$$      end
c$$$
c$$$********************************************************************************
c$$$********************************************************************************
c$$$********************************************************************************
c$$$********************************************************************************
c$$$
c$$$      subroutine whichsensor(nplPAM,xPAM,yPAM,ladder,sensor)
c$$$* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *   
c$$$*     Given the view and the (xPAM,yPAM) coordinates (in the PAMELA
c$$$*     reference system), it returns the ladder and the sensor 
c$$$*     which the point belongs to.
c$$$*
c$$$*     The method to assign a point to a sensor consists in 
c$$$*     - calculating the sum of the distances between the point
c$$$*     and the sensor edges
c$$$*     - requiring that it is less-equal than (SiDimX+SiDimY)
c$$$*     (should be strictly equal actually...)
c$$$*
c$$$*     NB -- SiDimX and SiDimY are not the dimentions of the SENSITIVE volume
c$$$*           but of the whole silicon sensor
c$$$* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_align.f'
c$$$
c$$$      integer ladder,sensor,viewx,viewy
c$$$      real c1(4),c2(4),c3(4)
c$$$      data c1/1.,0.,0.,1./
c$$$      data c2/1.,-1.,-1.,1./
c$$$      data c3/1.,1.,0.,0./
c$$$      real*8 yvvv,xvvv
c$$$      double precision xi,yi,zi
c$$$      double precision xrt,yrt,zrt
c$$$      real AA,BB
c$$$      real yvv(4),xvv(4)
c$$$
c$$$*     tollerance to consider the track inside the sensitive area
c$$$      real ptoll
c$$$      data ptoll/150./          !um
c$$$
c$$$      external nviewx,nviewy,acoordsi,dcoord
c$$$
c$$$
c$$$
c$$$      nplpt = nplPAM            !plane
c$$$      viewx = nviewx(nplpt)
c$$$      viewy = nviewy(nplpt)
c$$$
c$$$      do il=1,nladders_view
c$$$         do is=1,2
c$$$
c$$$            do iv=1,4           !loop on sensor vertexes
c$$$               stripx = (il-c1(iv))*1024 + c1(iv) + c2(iv)*3
c$$$               stripy = (il-c3(iv))*1024 + c3(iv)               
c$$$c------------------------------------------------------------------------
c$$$c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
c$$$c------------------------------------------------------------------------
c$$$               xi = acoordsi(stripx,viewx)
c$$$               yi = acoordsi(stripy,viewy)
c$$$               zi = 0.
c$$$c------------------------------------------------------------------------
c$$$c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c$$$c------------------------------------------------------------------------
c$$$c     N.B. I convert angles from microradiants to radiants
c$$$               xrt = xi
c$$$     $              - omega(nplpt,il,is)*yi 
c$$$     $              + gamma(nplpt,il,is)*zi
c$$$     $              + dx(nplpt,il,is)               
c$$$               yrt = omega(nplpt,il,is)*xi
c$$$     $              + yi
c$$$     $              - beta(nplpt,il,is)*zi
c$$$     $              + dy(nplpt,il,is)               
c$$$               zrt = -gamma(nplpt,il,is)*xi
c$$$     $              + beta(nplpt,il,is)*yi
c$$$     $              + zi
c$$$     $              + dz(nplpt,il,is)
c$$$c------------------------------------------------------------------------
c$$$c     measured coordinates (in cm) in PAMELA reference system
c$$$c------------------------------------------------------------------------
c$$$               yvvv = dcoord(yrt,viewy,il,is) / 1.d4
c$$$               xvvv = dcoord(xrt,viewx,il,is) / 1.d4
c$$$
c$$$               yvv(iv)=sngl(yvvv)
c$$$               xvv(iv)=sngl(xvvv)
c$$$c               print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '
c$$$c     $              ,iv,xvv(iv),yvv(iv)
c$$$            enddo               !end loop on sensor vertexes
c$$$
c$$$            dtot=0.
c$$$            do iside=1,4,2        !loop on sensor edges X
c$$$               iv1=iside
c$$$               iv2=mod(iside,4)+1
c$$$*     straight line passing trhough two consecutive vertexes
c$$$               AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))
c$$$               BB = yvv(iv1) - AA*xvv(iv1)
c$$$*     point along the straight line closer to the track
c$$$               xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
c$$$               yoo = AA*xoo + BB
c$$$*     sum of the distances
c$$$               dtot = dtot + 
c$$$     $              sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
c$$$            enddo               !end loop on sensor edges
c$$$            do iside=2,4,2        !loop on sensor edges Y
c$$$               iv1=iside
c$$$               iv2=mod(iside,4)+1
c$$$*     straight line passing trhough two consecutive vertexes
c$$$               AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))
c$$$               BB = xvv(iv1) - AA*yvv(iv1)
c$$$*     point along the straight line closer to the track
c$$$               yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
c$$$               xoo = AA*yoo + BB
c$$$*     sum of the distances
c$$$               dtot = dtot + 
c$$$     $              sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
c$$$            enddo               !end loop on sensor edges
c$$$
c$$$
c$$$*     half-perimeter of sensitive area
c$$$            Perim = 
c$$$     $            SiDimX - edgeX_l - edgeX_r
c$$$     $           +SiDimY - edgeY_l - edgeY_r
c$$$            Perim = (Perim + ptoll)/1.e4
c$$$            if(dtot.le.Perim)goto 100
c$$$
c$$$
c$$$         enddo
c$$$      enddo
c$$$      
c$$$      ladder = 0 
c$$$      sensor = 0 
c$$$      goto 200
c$$$
c$$$ 100  continue
c$$$      ladder = il
c$$$      sensor = is
c$$$      
c$$$
c$$$ 200  return
c$$$      end
c$$$
c$$$
c$$$
c$$$*************************************************************************
c$$$
c$$$      subroutine reverse(v,n,temp) !invert the order of the components of v(n) vector
c$$$
c$$$      implicit double precision (A-H,O-Z)
c$$$
c$$$      dimension v(*)
c$$$      dimension temp(*)
c$$$      integer i,n
c$$$
c$$$      do i=1,n
c$$$        temp(i)=v(n+1-i)
c$$$      enddo
c$$$
c$$$      do i=1,n
c$$$        v(i)=temp(i)
c$$$      enddo
c$$$
c$$$      return
c$$$      end
c$$$
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$      integer function ip_cp(id)
c$$$*
c$$$*     given the couple id, 
c$$$*     it returns the plane number
c$$$*     
c$$$      include '../common/commontracker.f'
c$$$c      include '../common/common_analysis.f'
c$$$      include '../common/common_momanhough.f'
c$$$      
c$$$      ip_cp=0
c$$$      ncpp=0
c$$$      do ip=1,nplanes
c$$$         ncpp=ncpp+ncp_plane(ip)
c$$$         if(ncpp.ge.abs(id))then
c$$$            ip_cp=ip
c$$$            goto 100
c$$$         endif
c$$$      enddo
c$$$ 100  continue
c$$$      return
c$$$      end
c$$$
c$$$
c$$$      integer function is_cp(id)
c$$$*
c$$$*     given the couple id, 
c$$$*     it returns the sensor number 
c$$$*     
c$$$      is_cp=0
c$$$      if(id.lt.0)is_cp=1
c$$$      if(id.gt.0)is_cp=2
c$$$      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
c$$$
c$$$      return
c$$$      end
c$$$
c$$$
c$$$      integer function icp_cp(id)
c$$$*
c$$$*     given the couple id, 
c$$$*     it returns the id number ON THE PLANE
c$$$*     
c$$$      include '../common/commontracker.f'
c$$$c      include '../common/common_analysis.f'
c$$$      include '../common/common_momanhough.f'
c$$$      
c$$$      icp_cp=0
c$$$
c$$$      ncpp=0
c$$$      do ip=1,nplanes
c$$$         ncppold=ncpp
c$$$         ncpp=ncpp+ncp_plane(ip)
c$$$         if(ncpp.ge.abs(id))then
c$$$            icp_cp=abs(id)-ncppold
c$$$            goto 100
c$$$         endif
c$$$      enddo
c$$$ 100  continue
c$$$      return
c$$$      end
c$$$      
c$$$
c$$$
c$$$      integer function id_cp(ip,icp,is)
c$$$*
c$$$*     given a plane, a couple and the sensor 
c$$$*     it returns the absolute couple id 
c$$$*     negative if sensor =1
c$$$*     positive if sensor =2
c$$$*
c$$$      include '../common/commontracker.f'
c$$$c      include '../common/calib.f'
c$$$c      include '../common/level1.f'
c$$$c      include '../common/common_analysis.f'
c$$$      include '../common/common_momanhough.f'
c$$$      
c$$$      id_cp=0
c$$$
c$$$      if(ip.gt.1)then
c$$$         do i=1,ip-1
c$$$            id_cp = id_cp + ncp_plane(i)
c$$$         enddo
c$$$      endif
c$$$
c$$$      id_cp = id_cp + icp
c$$$
c$$$      if(is.eq.1) id_cp = -id_cp
c$$$
c$$$      return
c$$$      end
c$$$
c$$$
c$$$
c$$$
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$*************************************************************************
c$$$      subroutine book_debug
c$$$
c$$$      include '../common/commontracker.f'
c$$$c      include '../common/calib.f'
c$$$c      include '../common/level1.f'
c$$$c      include '../common/common_analysis.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/common_level2debug.f'
c$$$
c$$$c      include '../common/common_mini.f'
c$$$
c$$$      character*35 block1,block2,block3,block4
c$$$     $     ,block5,block6
c$$$      
c$$$      
c$$$* * * * * * * * * * * * * * * * * * * * * * * *
c$$$*     HOUGH TRANSFORM PARAMETERS
c$$$c      call HBOOK1(1001
c$$$c     $     ,'y'
c$$$c     $     ,3000,-70.,70.,0.)
c$$$c      call HBOOK1(1002
c$$$c     $     ,'tg thyz'
c$$$c     $     ,300,-1.,1.,0.)         
c$$$      
c$$$      call HBOOK2(1003
c$$$     $     ,'y vs tg thyz'
c$$$     $     ,300,-1.,1.         !x         
c$$$     $     ,3000,-70.,70.,0.)   !y
c$$$
c$$$      call HBOOK1(1004
c$$$     $     ,'Dy'
c$$$     $     ,100,0.,2.,0.)  
c$$$
c$$$      call HBOOK1(1005
c$$$     $     ,'D thyz'
c$$$     $     ,100,0.,.05,0.)   
c$$$
c$$$
c$$$
c$$$*     DEBUG ntuple:
c$$$      call HBNT(ntp_level2+1,'LEVEL2',' ')
c$$$      
c$$$      call HBNAME(ntp_level2+1,'EVENT',good2_nt,
c$$$     $     'GOOD2:L,NEV2:I')
c$$$
c$$$ 411  format('NDBLT:I::[0,',I5,']')
c$$$      write(block1,411) ndblt_max_nt
c$$$      call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt,
c$$$     $     block1//'
c$$$     $     ,ALFAYZ1(NDBLT):R
c$$$     $     ,ALFAYZ2(NDBLT):R
c$$$     $     ')   
c$$$
c$$$ 412  format('NTRPT:I::[0,',I5,']')
c$$$      write(block2,412) ntrpt_max_nt
c$$$      call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt,
c$$$     $     block2//'
c$$$     $     ,ALFAXZ1(NTRPT):R 
c$$$     $     ,ALFAXZ2(NTRPT):R  
c$$$     $     ,ALFAXZ3(NTRPT):R
c$$$     $     ')
c$$$       
c$$$  
c$$$ 413  format('NCLOUDS_YZ:I::[0,',I4,']')
c$$$ 414  format('DB_CLOUD(',I4,'):I')
c$$$      write(block3,413) ncloyz_max
c$$$      write(block4,414) ndblt_max_nt
c$$$      call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ,
c$$$     $     block3//'
c$$$     $     ,ALFAYZ1_AV(NCLOUDS_YZ):R
c$$$     $     ,ALFAYZ2_AV(NCLOUDS_YZ):R
c$$$     $     ,PTCLOUD_YZ(NCLOUDS_YZ):I
c$$$     $     ,'//block4
c$$$     $     ) 
c$$$ 
c$$$ 415  format('NCLOUDS_XZ:I::[0,',I4,']')
c$$$ 416  format('TR_CLOUD(',I5,'):I')
c$$$      write(block5,415) ncloxz_max
c$$$      write(block6,416) ntrpt_max_nt
c$$$      call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ,
c$$$     $     block5//'
c$$$     $     ,ALFAXZ1_AV(NCLOUDS_XZ):R
c$$$     $     ,ALFAXZ2_AV(NCLOUDS_XZ):R
c$$$     $     ,ALFAXZ3_AV(NCLOUDS_XZ):R
c$$$     $     ,PTCLOUD_XZ(NCLOUDS_XZ):I
c$$$     $     ,'//block6
c$$$     $     ) 
c$$$
c$$$      
c$$$      return
c$$$      end
c$$$***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
c$$$*
c$$$*
c$$$*
c$$$*
c$$$*
c$$$*
c$$$*
c$$$*
c$$$*
c$$$***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
c$$$      subroutine book_level2
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/level2.f'
c$$$
c$$$      character*35 block1,block2
c$$$
c$$$c      print*,'__________ booking LEVEL2 n-tuple __________'
c$$$
c$$$*     LEVEL1 ntuple:
c$$$      call HBNT(ntp_level2,'LEVEL2',' ')
c$$$      
c$$$      call HBNAME(ntp_level2,'EVENT',good2,
c$$$     $     'GOOD2:L,NEV2:I')
c$$$
c$$$# ifndef TEST2003
c$$$
c$$$      call HBNAME(ntp_level2,'CPU',pkt_type
c$$$     $     ,'PKT_TYPE:I::[0,50]
c$$$     $     ,PKT_NUM:I
c$$$     $     ,OBT:I
c$$$     $     ,WHICH_CALIB:I::[0,50]')
c$$$
c$$$# endif
c$$$
c$$$ 417  format('NTRK:I::[0,',I4,']')
c$$$ 418  format(',IMAGE(NTRK):I::[0,',I4,']')
c$$$      write(block1,417)NTRKMAX
c$$$      write(block2,418)NTRKMAX
c$$$      call HBNAME(ntp_level2,'TRACKS',NTRK,
c$$$     $     block1//
c$$$     $     block2//'
c$$$     $     ,XM(6,NTRK):R 
c$$$     $     ,YM(6,NTRK):R
c$$$     $     ,ZM(6,NTRK):R
c$$$     $     ,RESX(6,NTRK):R
c$$$     $     ,RESY(6,NTRK):R
c$$$     $     ,AL(5,NTRK):R
c$$$     $     ,COVAL(5,5,NTRK):R
c$$$     $     ,CHI2(NTRK):R
c$$$     $     ,XGOOD(6,NTRK):I::[0,1]
c$$$     $     ,YGOOD(6,NTRK):I::[0,1]
c$$$     $     ,XV(6,NTRK):R 
c$$$     $     ,YV(6,NTRK):R
c$$$     $     ,ZV(6,NTRK):R
c$$$     $     ,AXV(6,NTRK):R 
c$$$     $     ,AYV(6,NTRK):R
c$$$     $     ,DEDXP(6,NTRK):R
c$$$     $     ')
c$$$
c$$$  
c$$$      call HBNAME(ntp_level2,'SINGLETS',nclsx,
c$$$     $     'NCLSX(6):I,NCLSY(6):I')
c$$$
c$$$      return
c$$$      end
c$$$
c$$$*     ****************************************************
c$$$
c$$$      subroutine init_level2
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/level2.f'
c$$$      include '../common/level1.f'
c$$$
c$$$
c$$$      good2 = .false.
c$$$      nev2 = nev1
c$$$
c$$$# ifndef TEST2003
c$$$
c$$$      pkt_type = pkt_type1
c$$$      pkt_num = pkt_num1
c$$$      obt = obt1
c$$$      which_calib = which_calib1
c$$$
c$$$# endif
c$$$
c$$$      NTRK = 0
c$$$      do it=1,NTRACKSMAX
c$$$         IMAGE(IT)=0
c$$$         CHI2_NT(IT) = -100000.
c$$$         do ip=1,nplanes
c$$$            XM_nt(IP,IT) = 0
c$$$            YM_nt(IP,IT) = 0
c$$$            ZM_nt(IP,IT) = 0
c$$$            RESX_nt(IP,IT) = 0
c$$$            RESY_nt(IP,IT) = 0
c$$$            XGOOD_nt(IP,IT) = 0
c$$$            YGOOD_nt(IP,IT) = 0
c$$$         enddo
c$$$         do ipa=1,5
c$$$            AL_nt(IPA,IT) = 0
c$$$            do ipaa=1,5
c$$$               coval(ipa,ipaa,IT)=0
c$$$            enddo                  
c$$$         enddo                  
c$$$      enddo
c$$$      
c$$$      do ip=1,nplanes
c$$$         nclsx(ip)=0
c$$$         nclsy(ip)=0
c$$$      enddo
c$$$
c$$$
c$$$      end
c$$$
c$$$*     -------------------------------------------------------
c$$$*     This routine fills the ntr-th element of the variables
c$$$*     inside the level2_tracks common, which correspond
c$$$*     to the ntr-th track info.
c$$$*     -------------------------------------------------------
c$$$
c$$$      subroutine fill_level2_tracks(ntr)
c$$$     
c$$$      include '../common/commontracker.f'
c$$$      include '../common/level2.f'
c$$$      include '../common/common_mini_2.f'
c$$$
c$$$c      print*,'TRACK ',ntr
c$$$      
c$$$      good2=.true.
c$$$      chi2_nt(ntr)        = sngl(chi2)
c$$$c      print*,chi2_nt(ntr)
c$$$      do i=1,5
c$$$         al_nt(i,ntr)     = sngl(al(i))
c$$$         do j=1,5
c$$$            coval(i,j,ntr) = sngl(cov(i,j))
c$$$         enddo
c$$$c         print*,al_nt(i,ntr)
c$$$      enddo
c$$$      do ip=1,nplanes           ! loop on planes
c$$$         xgood_nt(ip,ntr) = int(xgood(ip))
c$$$         ygood_nt(ip,ntr) = int(ygood(ip))
c$$$         xm_nt(ip,ntr)    = sngl(xm(ip))
c$$$         ym_nt(ip,ntr)    = sngl(ym(ip))
c$$$         zm_nt(ip,ntr)    = sngl(zm(ip))
c$$$         RESX_nt(IP,ntr)   = sngl(resx(ip))
c$$$         RESY_nt(IP,ntr)   = sngl(resy(ip))
c$$$         xv_nt(ip,ntr)    = sngl(xv(ip))
c$$$         yv_nt(ip,ntr)    = sngl(yv(ip))
c$$$         zv_nt(ip,ntr)    = sngl(zv(ip))
c$$$         axv_nt(ip,ntr)   = sngl(axv(ip))
c$$$         ayv_nt(ip,ntr)   = sngl(ayv(ip))
c$$$         dedxp(ip,ntr)    = sngl(dedxtrk(ip))
c$$$c         print*,ip,'|',xm_nt(ip,ntr),ym_nt(ip,ntr),zm_nt(ip,ntr)
c$$$      enddo
c$$$      
c$$$      end
c$$$
c$$$***************************************************
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$**************************************************
c$$$
c$$$      subroutine cl_to_couples(iflag)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/momanhough_init.f'
c$$$      include '../common/calib.f'
c$$$      include '../common/level1.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$*     output flag
c$$$*     --------------
c$$$*     0 = good event
c$$$*     1 = bad event
c$$$*     --------------
c$$$      integer iflag
c$$$
c$$$      integer badseed,badcl
c$$$
c$$$*     init variables
c$$$      ncp_tot=0 
c$$$      do ip=1,nplanes 
c$$$         do ico=1,ncouplemax
c$$$            clx(ip,ico)=0
c$$$            cly(ip,ico)=0
c$$$         enddo
c$$$         ncp_plane(ip)=0
c$$$         do icl=1,nclstrmax_level2
c$$$            cls(ip,icl)=1
c$$$         enddo
c$$$         ncls(ip)=0
c$$$      enddo
c$$$      do icl=1,nclstrmax_level2
c$$$         cl_single(icl)=1
c$$$         cl_good(icl)=0
c$$$      enddo
c$$$      
c$$$*     start association
c$$$      ncouples=0
c$$$      do icx=1,nclstr1          !loop on cluster (X)
c$$$         if(mod(VIEW(icx),2).eq.1)goto 10
c$$$         
c$$$*     ----------------------------------------------------
c$$$*     cut on charge (X VIEW)
c$$$         if(dedx(icx).lt.dedx_x_min)then
c$$$            cl_single(icx)=0
c$$$            goto 10
c$$$         endif
c$$$*     cut BAD (X VIEW)            
c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
c$$$         ifirst=INDSTART(icx)
c$$$         if(icx.ne.nclstr1) then
c$$$            ilast=INDSTART(icx+1)-1
c$$$         else
c$$$            ilast=TOTCLLENGTH
c$$$         endif
c$$$         badcl=badseed
c$$$         do igood=-ngoodstr,ngoodstr
c$$$            ibad=1
c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.
c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.
c$$$     $           .true.)then
c$$$               ibad=BAD(VIEW(icx),
c$$$     $              nvk(MAXS(icx)+igood),
c$$$     $              nst(MAXS(icx)+igood))
c$$$            endif
c$$$            badcl=badcl*ibad
c$$$         enddo
c$$$c     if(badcl.eq.0)then 
c$$$c     cl_single(icx)=0
c$$$c     goto 10
c$$$c     endif
c$$$*     ----------------------------------------------------
c$$$         
c$$$         cl_good(icx)=1
c$$$         nplx=npl(VIEW(icx)) 
c$$$         nldx=nld(MAXS(icx),VIEW(icx))
c$$$         
c$$$         do icy=1,nclstr1       !loop on cluster (Y)
c$$$            if(mod(VIEW(icy),2).eq.0)goto 20
c$$$            
c$$$*     ----------------------------------------------------
c$$$*     cut on charge (Y VIEW) 
c$$$            if(dedx(icy).lt.dedx_y_min)then 
c$$$               cl_single(icy)=0
c$$$               goto 20
c$$$            endif
c$$$*     cut BAD (Y VIEW)            
c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
c$$$            ifirst=INDSTART(icy)
c$$$            if(icy.ne.nclstr1) then
c$$$               ilast=INDSTART(icy+1)-1 
c$$$            else
c$$$               ilast=TOTCLLENGTH 
c$$$            endif
c$$$            badcl=badseed
c$$$            do igood=-ngoodstr,ngoodstr
c$$$               ibad=1
c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.
c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.
c$$$     $              .true.)
c$$$     $              ibad=BAD(VIEW(icy),
c$$$     $              nvk(MAXS(icy)+igood),
c$$$     $              nst(MAXS(icy)+igood))
c$$$               badcl=badcl*ibad
c$$$            enddo
c$$$c     if(badcl.eq.0)then 
c$$$c     cl_single(icy)=0
c$$$c     goto 20
c$$$c     endif
c$$$*     ----------------------------------------------------
c$$$            
c$$$            
c$$$            cl_good(icy)=1                  
c$$$            nply=npl(VIEW(icy)) 
c$$$            nldy=nld(MAXS(icy),VIEW(icy))
c$$$            
c$$$*     ----------------------------------------------
c$$$*     CONDITION TO FORM A COUPLE 
c$$$*     ----------------------------------------------
c$$$*     geometrical consistency (same plane and ladder)
c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then
c$$$*     charge correlation
c$$$               ddd=(dedx(icy)
c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
c$$$               cut=chcut*sch(nplx,nldx)
c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent
c$$$               
c$$$               
c$$$*     ------------------> COUPLE <------------------
c$$$*     check to do not overflow vector dimentions
c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
c$$$                  if(DEBUG)print*,
c$$$     $                    ' ** warning ** number of identified'// 
c$$$     $                    ' couples on plane ',nplx,
c$$$     $                    ' exceeds vector dimention'//
c$$$     $                    ' ( ',ncouplemax,' )'
c$$$c     good2=.false.
c$$$c     goto 880   !fill ntp and go to next event
c$$$                  iflag=1
c$$$                  return
c$$$               endif
c$$$               
c$$$               if(ncp_plane(nplx).eq.ncouplemax)then
c$$$                  if(DEBUG)print*,
c$$$     $                 '** warning ** number of identified '// 
c$$$     $                 'couples on plane ',nplx,
c$$$     $                 'exceeds vector dimention '
c$$$     $                 ,'( ',ncouplemax,' )'
c$$$c     good2=.false.
c$$$c     goto 880   !fill ntp and go to next event                     
c$$$                  iflag=1
c$$$                  return
c$$$               endif
c$$$               
c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1
c$$$               clx(nplx,ncp_plane(nplx))=icx
c$$$               cly(nply,ncp_plane(nplx))=icy
c$$$               cl_single(icx)=0
c$$$               cl_single(icy)=0
c$$$            endif                              
c$$$*     ----------------------------------------------
c$$$
c$$$ 20         continue
c$$$         enddo                  !end loop on clusters(Y)
c$$$         
c$$$ 10      continue
c$$$      enddo                     !end loop on clusters(X)
c$$$      
c$$$      
c$$$      do icl=1,nclstr1
c$$$         if(cl_single(icl).eq.1)then
c$$$            ip=npl(VIEW(icl)) 
c$$$            ncls(ip)=ncls(ip)+1
c$$$            cls(ip,ncls(ip))=icl
c$$$         endif
c$$$      enddo
c$$$      
c$$$      
c$$$      if(DEBUG)then
c$$$         print*,'clusters  ',nclstr1
c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)
c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)
c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
c$$$      endif
c$$$      
c$$$      do ip=1,6
c$$$         ncp_tot=ncp_tot+ncp_plane(ip)
c$$$      enddo
c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
c$$$      
c$$$      if(ncp_tot.gt.ncp_max)then
c$$$         if(DEBUG)print*,
c$$$     $           '** warning ** number of identified '// 
c$$$     $           'couples exceeds upper limit for Hough tr. '
c$$$     $           ,'( ',ncp_max,' )'            
c$$$c            good2=.false.
c$$$c     goto 880       !fill ntp and go to next event
c$$$         iflag=1
c$$$         return
c$$$      endif
c$$$      
c$$$      return
c$$$      end
c$$$      
c$$$      
c$$$***************************************************
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$**************************************************
c$$$
c$$$      subroutine cp_to_doubtrip(iflag)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/momanhough_init.f'
c$$$      include '../common/common_xyzPAM.f'
c$$$      include '../common/calib.f'
c$$$      include '../common/level1.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$*     output flag
c$$$*     --------------
c$$$*     0 = good event
c$$$*     1 = bad event
c$$$*     --------------
c$$$      integer iflag
c$$$
c$$$
c$$$*     -----------------------------
c$$$*     DOUBLETS/TRIPLETS coordinates
c$$$c      double precision xm1,ym1,zm1
c$$$c      double precision xm2,ym2,zm2
c$$$c      double precision xm3,ym3,zm3
c$$$
c$$$      real xm1,ym1,zm1
c$$$      real xm2,ym2,zm2
c$$$      real xm3,ym3,zm3
c$$$*     -----------------------------
c$$$*     variable needed for tricircle:
c$$$      real xp(3),zp(3)!TRIPLETS coordinates, to find a circle
c$$$      EQUIVALENCE (xm1,xp(1))
c$$$      EQUIVALENCE (xm2,xp(2))
c$$$      EQUIVALENCE (xm3,xp(3))
c$$$      EQUIVALENCE (zm1,zp(1))
c$$$      EQUIVALENCE (zm2,zp(2))
c$$$      EQUIVALENCE (zm3,zp(3))
c$$$      real angp(3),resp(3),chi
c$$$      real xc,zc,radius
c$$$*     -----------------------------
c$$$
c$$$
c$$$
c$$$      ndblt=0                   !number of doublets
c$$$      ntrpt=0                   !number of triplets
c$$$      
c$$$      do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
c$$$         do is1=1,2             !loop on sensors - COPPIA 1
c$$$            
c$$$            do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
c$$$               icx1=clx(ip1,icp1)
c$$$               icy1=cly(ip1,icp1)
c$$$               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)
c$$$               xm1=xPAM
c$$$               ym1=yPAM
c$$$               zm1=zPAM                  
c$$$c     print*,'***',is1,xm1,ym1,zm1
c$$$               do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
c$$$                  do is2=1,2    !loop on sensors -ndblt COPPIA 2
c$$$                     
c$$$                     do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
c$$$                        icx2=clx(ip2,icp2)
c$$$                        icy2=cly(ip2,icp2)
c$$$                        call xyz_PAM
c$$$     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)
c$$$                        xm2=xPAM
c$$$                        ym2=yPAM
c$$$                        zm2=zPAM
c$$$                                                
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$*     track parameters on Y VIEW
c$$$*     (2 couples needed)
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$                        if(ndblt.eq.ndblt_max)then
c$$$                           if(DEBUG)print*,
c$$$     $                          '** warning ** number of identified '// 
c$$$     $                          'doublets exceeds vector dimention '
c$$$     $                          ,'( ',ndblt_max,' )'
c$$$c                           good2=.false.
c$$$c                           goto 880 !fill ntp and go to next event
c$$$                           iflag=1
c$$$                           return
c$$$                        endif
c$$$                        ndblt = ndblt + 1
c$$$*     store doublet info
c$$$                        cpyz1(ndblt)=id_cp(ip1,icp1,is1)
c$$$                        cpyz2(ndblt)=id_cp(ip2,icp2,is2)
c$$$*     tg(th_yz)
c$$$                        alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
c$$$*     y0 (cm)
c$$$                        alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
c$$$                           
c$$$****  -----------------------------------------------****
c$$$****  reject non phisical couples                    ****
c$$$****  -----------------------------------------------****
c$$$                        if(
c$$$     $                       abs(alfayz2(ndblt)).gt.alfyz2_max
c$$$     $                       .or.
c$$$     $                       abs(alfayz1(ndblt)).gt.alfyz1_max
c$$$     $                       )ndblt = ndblt-1
c$$$                        
c$$$c$$$      if(iev.eq.33)then
c$$$c$$$      print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2
c$$$c$$$     $        ,' || ',icx1,icy1,icx2,icy2
c$$$c$$$     $        ,' || ',xm1,ym1,xm2,ym2
c$$$c$$$     $        ,' || ',alfayz2(ndblt),alfayz1(ndblt)
c$$$c$$$      endif
c$$$c$$$
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$*     track parameters on Y VIEW - end
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$ 
c$$$
c$$$                        if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples
c$$$                        do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
c$$$                           do is3=1,2 !loop on sensors - COPPIA 3
c$$$                              
c$$$                              do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
c$$$                                 icx3=clx(ip3,icp3)
c$$$                                 icy3=cly(ip3,icp3)
c$$$                                 call xyz_PAM
c$$$     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)
c$$$                                 xm3=xPAM
c$$$                                 ym3=yPAM
c$$$                                 zm3=zPAM
c$$$*     find the circle passing through the three points
c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
c$$$     $                                ,xc,zc,radius,iflag)
c$$$c     print*,xc,zc,radius
c$$$*     the circle must intersect the reference plane
c$$$                                 if(
c$$$c     $                                 (xc.le.-1.*xclimit.or.
c$$$c     $                                 xc.ge.xclimit).and.
c$$$     $                                radius**2.ge.(ZINI-zc)**2.and.
c$$$     $                                iflag.eq.0.and.
c$$$     $                                .true.)then
c$$$                                 
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$*     track parameters on X VIEW
c$$$*     (3 couples needed)
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$                                 if(ntrpt.eq.ntrpt_max)then
c$$$                                    if(DEBUG)print*,
c$$$     $                     '** warning ** number of identified '// 
c$$$     $                     'triplets exceeds vector dimention '
c$$$     $                    ,'( ',ntrpt_max,' )'
c$$$c                                    good2=.false.
c$$$c                                    goto 880 !fill ntp and go to next event
c$$$                                    iflag=1
c$$$                                    return
c$$$                                 endif
c$$$                                 ntrpt = ntrpt +1
c$$$*     store triplet info
c$$$                                 cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
c$$$                                 cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
c$$$                                 cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
c$$$                                 
c$$$                                 if(xc.lt.0)then
c$$$*************POSITIVE DEFLECTION
c$$$              alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
c$$$              alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
c$$$              alfaxz3(ntrpt) = 1/radius
c$$$                                 else
c$$$*************NEGATIVE DEFLECTION
c$$$              alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
c$$$              alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
c$$$              alfaxz3(ntrpt) = -1/radius
c$$$                                 endif
c$$$                                 
c$$$****  -----------------------------------------------****
c$$$****  reject non phisical triplets                   ****
c$$$****  -----------------------------------------------****
c$$$                                 if(
c$$$     $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max
c$$$     $                                .or.
c$$$     $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max
c$$$     $                                )ntrpt = ntrpt-1
c$$$                                 
c$$$                                 
c$$$c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$*     track parameters on X VIEW - end
c$$$*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c$$$                              endif
c$$$                           enddo !end loop on COPPIA 3
c$$$                        enddo   !end loop on sensors - COPPIA 3
c$$$                     enddo      !end loop on planes  - COPPIA 3
c$$$ 30                  continue
c$$$                     
c$$$ 1                enddo         !end loop on COPPIA 2
c$$$               enddo            !end loop on sensors - COPPIA 2
c$$$            enddo               !end loop on planes  - COPPIA 2
c$$$            
c$$$         enddo                  !end loop on COPPIA1
c$$$      enddo                     !end loop on sensors - COPPIA 1
c$$$      enddo                     !end loop on planes  - COPPIA 1
c$$$      
c$$$      if(DEBUG)then
c$$$         print*,'--- doublets ',ndblt
c$$$         print*,'--- triplets ',ntrpt
c$$$      endif
c$$$      
c$$$c     goto 880               !ntp fill
c$$$      
c$$$      
c$$$      return
c$$$      end
c$$$
c$$$
c$$$
c$$$***************************************************
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$**************************************************
c$$$
c$$$      subroutine doub_to_YZcloud(iflag)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/momanhough_init.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$*     output flag
c$$$*     --------------
c$$$*     0 = good event
c$$$*     1 = bad event
c$$$*     --------------
c$$$      integer iflag
c$$$
c$$$      integer db_used(ndblt_max)
c$$$      integer db_temp(ndblt_max)
c$$$      integer db_all(ndblt_max) !stores db ID in each cloud
c$$$
c$$$      integer hit_plane(nplanes)
c$$$
c$$$*     mask for used couples 
c$$$      integer cp_useds1(ncouplemaxtot) ! sensor 1
c$$$      integer cp_useds2(ncouplemaxtot) ! sensor 2
c$$$
c$$$*     ******************************************
c$$$*     PLANE TO EXCLUDE (for efficiency studies)
c$$$*     --> provided as input parameter
c$$$*     (set according to "tracking" notation): 
c$$$*            1-6 from TOP to BOTTOM
c$$$*     ******************************************
c$$$      integer ipexclude
c$$$      common/plane_exclude/ipexclude
c$$$
c$$$
c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c$$$*     classification of DOUBLETS
c$$$*     according to distance in parameter space
c$$$*     (cloud = group of points (doublets) in parameter space)
c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c$$$      do idb=1,ndblt
c$$$         db_used(idb)=0
c$$$      enddo
c$$$
c$$$      distance=0
c$$$      nclouds_yz=0              !number of clouds
c$$$      npt_tot=0
c$$$      do idb1=1,ndblt           !loop (1) on DOUBLETS
c$$$         if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
c$$$*     ***************************************
c$$$*     Checks if a plane should be excluded 
c$$$*     (for efficiency study)
c$$$*     ***************************************
c$$$         ip1=ip_cp(cpyz1(idb1))
c$$$         ip2=ip_cp(cpyz2(idb1))
c$$$         if(
c$$$     $        ip1.eq.(nplanes-ipexclude+1).or.
c$$$     $        ip2.eq.(nplanes-ipexclude+1).or.
c$$$     $        .false.)goto 2228        !exclude DOUBLET from cloud
c$$$*     ***************************************
c$$$
c$$$c     print*,'--------------'
c$$$c     print*,'** ',idb1,' **'
c$$$            
c$$$         do icp=1,ncp_tot
c$$$            cp_useds1(icp)=0    !init
c$$$            cp_useds2(icp)=0    !init
c$$$         enddo
c$$$         do idb=1,ndblt
c$$$            db_all(idb)=0
c$$$         enddo
c$$$         if(cpyz1(idb1).gt.0)cp_useds2(cpyz1(idb1))=1
c$$$         if(cpyz1(idb1).lt.0)cp_useds1(-cpyz1(idb1))=1
c$$$         if(cpyz2(idb1).gt.0)cp_useds2(cpyz2(idb1))=1
c$$$         if(cpyz2(idb1).lt.0)cp_useds1(-cpyz2(idb1))=1
c$$$         temp1 = alfayz1(idb1)
c$$$         temp2 = alfayz2(idb1)
c$$$         npt=1                  !counter of points in the cloud
c$$$
c$$$         db_all(npt) = idb1
c$$$         
c$$$         nptloop=1
c$$$         db_temp(1)=idb1
c$$$         
c$$$ 88      continue
c$$$         
c$$$         npv=0                  !# new points inlcuded
c$$$         do iloop=1,nptloop
c$$$            idbref=db_temp(iloop) !local point of reference
c$$$ccccc if(db_used(idbref).eq.1)goto 1188 !next
c$$$            
c$$$            do idb2=1,ndblt     !loop (2) on DOUBLETS
c$$$               if(idb2.eq.idbref)goto 1118 !next doublet
c$$$               if(db_used(idb2).eq.1)goto 1118
c$$$*     ***************************************
c$$$*     Checks if a plane should be excluded 
c$$$*     (for efficiency study)
c$$$*     ***************************************
c$$$               ip1=ip_cp(cpyz1(idb2))
c$$$               ip2=ip_cp(cpyz2(idb2))
c$$$               if(
c$$$     $              ip1.eq.(nplanes-ipexclude+1).or.
c$$$     $              ip2.eq.(nplanes-ipexclude+1).or.
c$$$     $              .false.)goto 1118 !exclude DOUBLET from cloud
c$$$*     ***************************************
c$$$               
c$$$                  
c$$$*     doublet distance in parameter space
c$$$               distance=
c$$$     $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
c$$$     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2               
c$$$               distance = sqrt(distance)
c$$$               
c$$$c$$$      if(iev.eq.33)then
c$$$c$$$      if(distance.lt.100)
c$$$c$$$     $ print*,'********* ',idb1,idbref,idb2,distance
c$$$c$$$      if(distance.lt.100)
c$$$c$$$     $ print*,'********* ',alfayz1(idbref),alfayz1(idb2)
c$$$c$$$     $                    ,alfayz2(idbref),alfayz2(idb2)
c$$$c$$$      endif
c$$$               if(distance.lt.cutdistyz)then
c$$$
c$$$c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
c$$$                  if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
c$$$                  if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
c$$$                  if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
c$$$                  if(cpyz2(idb2).lt.0)cp_useds1(-cpyz2(idb2))=1
c$$$                  npt = npt + 1 !counter of points in the cloud
c$$$
c$$$                  npv = npv +1
c$$$                  db_temp(npv) = idb2
c$$$                  db_used(idbref) = 1 
c$$$                  db_used(idb2) = 1  
c$$$                  
c$$$                  db_all(npt) = idb2
c$$$
c$$$                  temp1 = temp1 + alfayz1(idb2)
c$$$                  temp2 = temp2 + alfayz2(idb2)
c$$$c     print*,'*   idbref,idb2 ',idbref,idb2
c$$$               endif               
c$$$               
c$$$ 1118          continue 
c$$$            enddo               !end loop (2) on DOUBLETS
c$$$            
c$$$ 1188       continue
c$$$         enddo                  !end loop on... bo?
c$$$         
c$$$         nptloop=npv
c$$$         if(nptloop.ne.0)goto 88
c$$$         
c$$$*     ------------------------------------------
c$$$*     stores the cloud only if 
c$$$*     1) it includes a minimum number of REAL couples
c$$$*     1bis) it inlcudes a minimum number of doublets
c$$$*     2) it is not already stored
c$$$*     ------------------------------------------
c$$$         do ip=1,nplanes
c$$$            hit_plane(ip)=0
c$$$         enddo
c$$$         ncpused=0
c$$$         do icp=1,ncp_tot
c$$$            if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
c$$$               ncpused=ncpused+1
c$$$               ip=ip_cp(icp)
c$$$               hit_plane(ip)=1
c$$$            endif
c$$$         enddo
c$$$         nplused=0
c$$$         do ip=1,nplanes
c$$$            nplused=nplused+ hit_plane(ip)
c$$$         enddo
c$$$c     print*,'>>>> ',ncpused,npt,nplused
c$$$         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
c$$$         if(npt.lt.nptyz_min)goto 2228 !next doublet
c$$$         if(nplused.lt.nplyz_min)goto 2228 !next doublet
c$$$         
c$$$*     ~~~~~~~~~~~~~~~~~
c$$$*     >>> NEW CLOUD <<<
c$$$
c$$$         if(nclouds_yz.ge.ncloyz_max)then
c$$$            if(DEBUG)print*,
c$$$     $           '** warning ** number of identified '// 
c$$$     $           'YZ clouds exceeds vector dimention '
c$$$     $           ,'( ',ncloyz_max,' )'
c$$$c               good2=.false.
c$$$c     goto 880         !fill ntp and go to next event
c$$$            iflag=1
c$$$            return
c$$$         endif
c$$$         
c$$$         nclouds_yz = nclouds_yz + 1 !increase counter
c$$$         alfayz1_av(nclouds_yz) = temp1/npt !store average parameter 
c$$$         alfayz2_av(nclouds_yz) = temp2/npt ! "
c$$$         do icp=1,ncp_tot
c$$$            cpcloud_yz(nclouds_yz,icp)=
c$$$     $           cp_useds1(icp)+2*cp_useds2(icp) !store cp info
c$$$         enddo 
c$$$         ptcloud_yz(nclouds_yz)=npt 
c$$$c     ptcloud_yz_nt(nclouds_yz)=npt 
c$$$         do ipt=1,npt
c$$$            db_cloud(npt_tot+ipt) = db_all(ipt)
c$$$c     print*,'>> ',ipt,db_all(ipt)
c$$$         enddo  
c$$$         npt_tot=npt_tot+npt
c$$$         if(DEBUG)then
c$$$            print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
c$$$            print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
c$$$            print*,'- alfayz1 ',alfayz1_av(nclouds_yz)
c$$$            print*,'- alfayz2 ',alfayz2_av(nclouds_yz)
c$$$            print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
c$$$            print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
c$$$            print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
c$$$         endif
c$$$*     >>> NEW CLOUD <<<
c$$$*     ~~~~~~~~~~~~~~~~~
c$$$ 2228    continue
c$$$      enddo                     !end loop (1) on DOUBLETS
c$$$      
c$$$      
c$$$      if(DEBUG)then
c$$$         print*,'---------------------- '
c$$$         print*,'Y-Z total clouds ',nclouds_yz
c$$$         print*,' '
c$$$      endif
c$$$      
c$$$      
c$$$      return
c$$$      end
c$$$      
c$$$
c$$$
c$$$
c$$$
c$$$***************************************************
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$**************************************************
c$$$
c$$$      subroutine trip_to_XZcloud(iflag)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/momanhough_init.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$*     output flag
c$$$*     --------------
c$$$*     0 = good event
c$$$*     1 = bad event
c$$$*     --------------
c$$$      integer iflag
c$$$
c$$$      integer tr_used(ntrpt_max)
c$$$      integer tr_temp(ntrpt_max)
c$$$      integer tr_incl(ntrpt_max)
c$$$      integer tr_all(ntrpt_max) !stores tr ID in each cloud
c$$$
c$$$      integer hit_plane(nplanes)
c$$$
c$$$*     mask for used couples 
c$$$      integer cp_useds1(ncouplemaxtot) ! sensor 1
c$$$      integer cp_useds2(ncouplemaxtot) ! sensor 2
c$$$
c$$$*     ******************************************
c$$$*     PLANE TO EXCLUDE (for efficiency studies)
c$$$*     --> provided as input parameter
c$$$*     (set according to "tracking" notation): 
c$$$*            1-6 from TOP to BOTTOM
c$$$*     ******************************************
c$$$      integer ipexclude
c$$$      common/plane_exclude/ipexclude
c$$$
c$$$
c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c$$$*     classification of TRIPLETS
c$$$*     according to distance in parameter space
c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c$$$      do itr=1,ntrpt
c$$$         tr_used(itr)=0
c$$$      enddo
c$$$
c$$$      distance=0
c$$$      nclouds_xz=0              !number of clouds         
c$$$      npt_tot=0                 !total number of selected triplets
c$$$      do itr1=1,ntrpt           !loop (1) on TRIPLETS
c$$$         if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
c$$$*     ***************************************
c$$$*     Checks if a plane should be excluded 
c$$$*     (for efficiency study)
c$$$*     ***************************************
c$$$         ip1=ip_cp(cpxz1(itr1))
c$$$         ip2=ip_cp(cpxz2(itr1))
c$$$         ip3=ip_cp(cpxz3(itr1))
c$$$         if(
c$$$     $        ip1.eq.(nplanes-ipexclude+1).or.
c$$$     $        ip2.eq.(nplanes-ipexclude+1).or.
c$$$     $        ip3.eq.(nplanes-ipexclude+1).or.
c$$$     $        .false.)goto 22288        !exclude TRIPLET from cloud
c$$$*     ***************************************
c$$$c     print*,'--------------'
c$$$c     print*,'** ',itr1,' **'
c$$$         
c$$$         do icp=1,ncp_tot
c$$$            cp_useds1(icp)=0
c$$$            cp_useds2(icp)=0
c$$$         enddo
c$$$         do itr=1,ntrpt
c$$$            tr_all(itr)=0       !list of included triplets
c$$$         enddo
c$$$         if(cpxz1(itr1).gt.0)cp_useds2(cpxz1(itr1))=1
c$$$         if(cpxz1(itr1).lt.0)cp_useds1(-cpxz1(itr1))=1
c$$$         if(cpxz2(itr1).gt.0)cp_useds2(cpxz2(itr1))=1
c$$$         if(cpxz2(itr1).lt.0)cp_useds1(-cpxz2(itr1))=1
c$$$         if(cpxz3(itr1).gt.0)cp_useds2(cpxz3(itr1))=1
c$$$         if(cpxz3(itr1).lt.0)cp_useds1(-cpxz3(itr1))=1
c$$$         temp1 = alfaxz1(itr1)
c$$$         temp2 = alfaxz2(itr1)
c$$$         temp3 = alfaxz3(itr1) 
c$$$         npt=1                  !counter of points in the cloud
c$$$         
c$$$         tr_all(npt) = itr1
c$$$         
c$$$         nptloop=1
c$$$c         tr_temp(1)=itr1
c$$$         tr_incl(1)=itr1
c$$$         
c$$$ 8881    continue
c$$$         
c$$$c     print*,'start search pv ',nptloop
c$$$         npv=0                  !# new points inlcuded
c$$$         do iloop=1,nptloop
c$$$c            itrref=tr_temp(iloop) !local point of reference
c$$$            itrref=tr_incl(iloop) !local point of reference
c$$$c     print*,'ref ',itrref
c$$$ccccc if(tr_used(itrref).eq.1)goto 11888 !next
c$$$            do itr2=1,ntrpt     !loop (2) on TRIPLETS
c$$$               if(itr2.eq.itr1)goto 11188       !next triplet
c$$$               if(tr_used(itr2).eq.1)goto 11188 !next triplet
c$$$*     ***************************************
c$$$*     Checks if a plane should be excluded 
c$$$*     (for efficiency study)
c$$$*     ***************************************
c$$$               ip1=ip_cp(cpxz1(itr2))
c$$$               ip2=ip_cp(cpxz2(itr2))
c$$$               ip3=ip_cp(cpxz3(itr2))
c$$$               if(
c$$$     $              ip1.eq.(nplanes-ipexclude+1).or.
c$$$     $              ip2.eq.(nplanes-ipexclude+1).or.
c$$$     $              ip3.eq.(nplanes-ipexclude+1).or.
c$$$     $              .false.)goto 11188 !exclude TRIPLET from cloud
c$$$*     ***************************************
c$$$               
c$$$*     triplet distance in parameter space
c$$$*     solo i due parametri spaziali per il momemnto
c$$$               distance=
c$$$     $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
c$$$     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2               
c$$$               distance = sqrt(distance)
c$$$               
c$$$               if(distance.lt.cutdistxz)then
c$$$c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
c$$$                  if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
c$$$                  if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
c$$$                  if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
c$$$                  if(cpxz2(itr2).lt.0)cp_useds1(-cpxz2(itr2))=1
c$$$                  if(cpxz3(itr2).gt.0)cp_useds2(cpxz3(itr2))=1
c$$$                  if(cpxz3(itr2).lt.0)cp_useds1(-cpxz3(itr2))=1
c$$$                  npt = npt + 1 !counter of points in the cloud
c$$$                  
c$$$                  npv = npv +1
c$$$                  tr_temp(npv) = itr2
c$$$                  tr_used(itrref) = 1
c$$$                  tr_used(itr2) = 1
c$$$                  
c$$$                  tr_all(npt) = itr2
c$$$                  
c$$$                  temp1 = temp1 + alfaxz1(itr2)
c$$$                  temp2 = temp2 + alfaxz2(itr2) 
c$$$                  temp3 = temp3 + alfaxz3(itr2)
c$$$c     print*,'*   itrref,itr2 ',itrref,itr2,distance 
c$$$               endif               
c$$$               
c$$$11188          continue 
c$$$            enddo               !end loop (2) on TRIPLETS
c$$$                       
c$$$11888       continue
c$$$         enddo                  !end loop on... bo?     
c$$$         
c$$$         nptloop=npv
c$$$         do i=1,npv
c$$$            tr_incl(i)=tr_temp(i)
c$$$         enddo
c$$$         if(nptloop.ne.0)goto 8881 
c$$$         
c$$$*     ------------------------------------------
c$$$*     stores the cloud only if 
c$$$*     1) it includes a minimum number of REAL couples
c$$$*     1bis)
c$$$*     2) it is not already stored
c$$$*     ------------------------------------------
c$$$c     print*,'check cp_used'
c$$$         do ip=1,nplanes
c$$$            hit_plane(ip)=0
c$$$         enddo
c$$$         ncpused=0
c$$$         do icp=1,ncp_tot
c$$$            if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
c$$$               ncpused=ncpused+1
c$$$               ip=ip_cp(icp)
c$$$               hit_plane(ip)=1
c$$$            endif
c$$$         enddo
c$$$         nplused=0
c$$$         do ip=1,nplanes
c$$$            nplused=nplused+ hit_plane(ip)
c$$$         enddo
c$$$         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
c$$$         if(npt.lt.nptxz_min)goto 22288 !next triplet
c$$$         if(nplused.lt.nplxz_min)goto 22288 !next doublet
c$$$         
c$$$*     ~~~~~~~~~~~~~~~~~
c$$$*     >>> NEW CLOUD <<<
c$$$         if(nclouds_xz.ge.ncloxz_max)then
c$$$            if(DEBUG)print*,
c$$$     $           '** warning ** number of identified '// 
c$$$     $           'XZ clouds exceeds vector dimention '
c$$$     $           ,'( ',ncloxz_max,' )'
c$$$c     good2=.false.
c$$$c     goto 880         !fill ntp and go to next event
c$$$            iflag=1
c$$$            return
c$$$         endif
c$$$         nclouds_xz = nclouds_xz + 1 !increase counter
c$$$         alfaxz1_av(nclouds_xz) = temp1/npt !store average parameter 
c$$$         alfaxz2_av(nclouds_xz) = temp2/npt ! "
c$$$         alfaxz3_av(nclouds_xz) = temp3/npt ! " 
c$$$         do icp=1,ncp_tot
c$$$            cpcloud_xz(nclouds_xz,icp)=
c$$$     $           cp_useds1(icp)+2*cp_useds2(icp) !store cp info
c$$$         enddo
c$$$         ptcloud_xz(nclouds_xz)=npt
c$$$c     ptcloud_xz_nt(nclouds_xz)=npt
c$$$         do ipt=1,npt
c$$$            tr_cloud(npt_tot+ipt) = tr_all(ipt)
c$$$c     print *,nclouds_xz,ipt,tr_all(ipt),npt_tot+ipt,
c$$$c     $              tr_cloud(npt_tot+ipt) 
c$$$         enddo
c$$$         npt_tot=npt_tot+npt
c$$$         
c$$$         if(DEBUG)then
c$$$            print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
c$$$            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'               
c$$$            print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
c$$$            print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
c$$$            print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
c$$$            print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
c$$$            print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
c$$$            print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
c$$$         endif
c$$$*     >>> NEW CLOUD <<<
c$$$*     ~~~~~~~~~~~~~~~~~
c$$$22288    continue
c$$$      enddo                     !end loop (1) on DOUBLETS
c$$$      
c$$$      if(DEBUG)then
c$$$         print*,'---------------------- '
c$$$         print*,'X-Z total clouds ',nclouds_xz
c$$$         print*,' '
c$$$      endif
c$$$      
c$$$      
c$$$      return
c$$$      end
c$$$      
c$$$
c$$$***************************************************
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$**************************************************
c$$$
c$$$      subroutine clouds_to_ctrack(iflag)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/common_xyzPAM.f'
c$$$      include '../common/common_mini_2.f'
c$$$      include '../common/common_mech.f'
c$$$      include '../common/momanhough_init.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$*     output flag
c$$$*     --------------
c$$$*     0 = good event
c$$$*     1 = bad event
c$$$*     --------------
c$$$      integer iflag
c$$$
c$$$*     -----------------------------------------------------------
c$$$*     mask to store (locally) the couples included
c$$$*     in the intersection bewteen a XZ and YZ cloud
c$$$      integer cpintersec(ncouplemaxtot)
c$$$*     -----------------------------------------------------------
c$$$*     list of matching couples in the combination 
c$$$*     between a XZ and YZ cloud
c$$$      integer cp_match(nplanes,ncouplemax)
c$$$      integer ncp_match(nplanes)
c$$$*     -----------------------------------------------------------
c$$$      integer hit_plane(nplanes)
c$$$*     -----------------------------------------------------------
c$$$*     variables for track fitting
c$$$      double precision AL_INI(5)
c$$$      double precision tath
c$$$*     -----------------------------------------------------------
c$$$c      real fitz(nplanes)        !z coordinates of the planes in cm
c$$$
c$$$
c$$$      ntracks=0                 !counter of track candidates
c$$$      
c$$$      do iyz=1,nclouds_yz       !loop on YZ couds
c$$$         do ixz=1,nclouds_xz    !loop on XZ couds
c$$$            
c$$$*     --------------------------------------------------
c$$$*     check of consistency of the clouds
c$$$*     ---> required a minimum number of matching couples
c$$$*     the track fit will be performed on the INTERSECTION
c$$$*     of the two clouds
c$$$*     --------------------------------------------------
c$$$            do ip=1,nplanes
c$$$               hit_plane(ip)=0
c$$$               ncp_match(ip)=0
c$$$               do icpp=1,ncouplemax
c$$$                  cp_match(ip,icpp)=0 !init couple list
c$$$               enddo
c$$$            enddo
c$$$            ncp_ok=0
c$$$            do icp=1,ncp_tot    !loop on couples
c$$$*     get info on 
c$$$               cpintersec(icp)=min(
c$$$     $              cpcloud_yz(iyz,icp),
c$$$     $              cpcloud_xz(ixz,icp))
c$$$               if(
c$$$     $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
c$$$     $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
c$$$     $              .false.)cpintersec(icp)=0
c$$$               if(cpintersec(icp).ne.0)then
c$$$                  ncp_ok=ncp_ok+1   
c$$$                  
c$$$                  ip=ip_cp(icp)
c$$$                  hit_plane(ip)=1
c$$$                  if(cpintersec(icp).eq.1)then 
c$$$*     1) only the couple image in sensor 1 matches
c$$$                     id=-icp
c$$$                     ncp_match(ip)=ncp_match(ip)+1
c$$$                     cp_match(ip,ncp_match(ip))=id
c$$$                  elseif(cpintersec(icp).eq.2)then
c$$$*     2) only the couple image in sensor 2 matches
c$$$                     id=icp
c$$$                     ncp_match(ip)=ncp_match(ip)+1
c$$$                     cp_match(ip,ncp_match(ip))=id
c$$$                  else
c$$$*     3) both couple images match
c$$$                     id=icp
c$$$                     do is=1,2
c$$$                        id=-id
c$$$                        ncp_match(ip)=ncp_match(ip)+1
c$$$                        cp_match(ip,ncp_match(ip))=id
c$$$                     enddo
c$$$                  endif                     
c$$$               endif            !end matching condition
c$$$            enddo               !end loop on couples
c$$$            
c$$$            nplused=0
c$$$            do ip=1,nplanes
c$$$               nplused=nplused+ hit_plane(ip)
c$$$            enddo
c$$$            
c$$$            if(nplused.lt.nplxz_min)goto 888 !next doublet
c$$$            if(ncp_ok.lt.ncpok)goto 888 !next cloud
c$$$            
c$$$            if(DEBUG)then
c$$$               print*,'Combination ',iyz,ixz
c$$$     $              ,' db ',ptcloud_yz(iyz)
c$$$     $              ,' tr ',ptcloud_xz(ixz)
c$$$     $              ,'  -----> # matching couples ',ncp_ok
c$$$            endif
c$$$c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
c$$$c$$$  print*,'Configurazione cluster XZ'
c$$$c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'Configurazione cluster YZ'
c$$$c$$$  print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))
c$$$c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
c$$$            
c$$$*     -------> INITIAL GUESS <-------
c$$$            AL_INI(1)=dreal(alfaxz1_av(ixz))
c$$$            AL_INI(2)=dreal(alfayz1_av(iyz))
c$$$            AL_INI(4)=datan(dreal(alfayz2_av(iyz))
c$$$     $           /dreal(alfaxz2_av(ixz)))
c$$$            tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
c$$$            AL_INI(3)=tath/sqrt(1+tath**2)
c$$$            AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
c$$$            
c$$$c     print*,'*******',AL_INI(5)
c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
c$$$            
c$$$c     print*,'alfaxz2, alfayz2 '
c$$$c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)
c$$$            
c$$$*     -------> INITIAL GUESS <------- 
c$$$c     print*,'AL_INI ',(al_ini(i),i=1,5)
c$$$            
c$$$            if(DEBUG)then
c$$$               print*,'1 >>> ',(cp_match(1,i),i=1,ncp_match(1))
c$$$               print*,'2 >>> ',(cp_match(2,i),i=1,ncp_match(2))
c$$$               print*,'3 >>> ',(cp_match(3,i),i=1,ncp_match(3))
c$$$               print*,'4 >>> ',(cp_match(4,i),i=1,ncp_match(4))
c$$$               print*,'5 >>> ',(cp_match(5,i),i=1,ncp_match(5))
c$$$               print*,'6 >>> ',(cp_match(6,i),i=1,ncp_match(6))
c$$$            endif
c$$$            
c$$$            do icp1=1,max(1,ncp_match(1))
c$$$               hit_plane(1)=icp1
c$$$               if(ncp_match(1).eq.0)hit_plane(1)=0 !-icp1
c$$$               
c$$$               do icp2=1,max(1,ncp_match(2))
c$$$                  hit_plane(2)=icp2
c$$$                  if(ncp_match(2).eq.0)hit_plane(2)=0 !-icp2
c$$$                  
c$$$                  do icp3=1,max(1,ncp_match(3))
c$$$                     hit_plane(3)=icp3
c$$$                     if(ncp_match(3).eq.0)hit_plane(3)=0 !-icp3
c$$$                     
c$$$                     do icp4=1,max(1,ncp_match(4))
c$$$                        hit_plane(4)=icp4
c$$$                        if(ncp_match(4).eq.0)hit_plane(4)=0 !-icp4
c$$$                        
c$$$                        do icp5=1,max(1,ncp_match(5))
c$$$                           hit_plane(5)=icp5
c$$$                           if(ncp_match(5).eq.0)hit_plane(5)=0 !-icp5
c$$$                           
c$$$                           do icp6=1,max(1,ncp_match(6))
c$$$                              hit_plane(6)=icp6
c$$$                              if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
c$$$                              
c$$$                              
c$$$                              call track_init !init TRACK common
c$$$
c$$$                              do ip=1,nplanes !loop on planes
c$$$                                 if(hit_plane(ip).ne.0)then 
c$$$                                    id=cp_match(ip,hit_plane(ip))
c$$$                                    is=is_cp(id)
c$$$                                    icp=icp_cp(id)
c$$$                                    if(ip_cp(id).ne.ip)
c$$$     $                                   print*,'OKKIO!!'
c$$$     $                                   ,'id ',id,is,icp
c$$$     $                                   ,ip_cp(id),ip
c$$$                                    icx=clx(ip,icp)
c$$$                                    icy=cly(ip,icp)
c$$$*                                   *************************
c$$$                                    call xyz_PAM(icx,icy,is,
c$$$     $                                   'COG2','COG2',0.,0.)
c$$$*                                   *************************
c$$$*                                   -----------------------------
c$$$                                    xgood(nplanes-ip+1)=1.
c$$$                                    ygood(nplanes-ip+1)=1.
c$$$                                    xm(nplanes-ip+1)=xPAM
c$$$                                    ym(nplanes-ip+1)=yPAM
c$$$                                    zm(nplanes-ip+1)=zPAM
c$$$                                    resx(nplanes-ip+1)=resxPAM
c$$$                                    resy(nplanes-ip+1)=resyPAM
c$$$*                                   -----------------------------
c$$$                                 endif
c$$$                              enddo !end loop on planes
c$$$*     **********************************************************
c$$$*     ************************** FIT *** FIT *** FIT *** FIT ***
c$$$*     **********************************************************
c$$$                              do i=1,5
c$$$                                 AL(i)=AL_INI(i)
c$$$                              enddo
c$$$                              ifail=0 !error flag in chi^2 computation
c$$$                              jstep=0 !number of  minimization steps
c$$$                              call mini_2(jstep,ifail)
c$$$                              if(ifail.ne.0) then
c$$$                                 if(DEBUG)then
c$$$                                    print *,
c$$$     $                              '*** MINIMIZATION FAILURE *** ' 
c$$$     $                              //'(mini_2 in clouds_to_ctrack)'
c$$$                                 endif
c$$$                                 chi2=-chi2 
c$$$                              endif 
c$$$*     **********************************************************
c$$$*     ************************** FIT *** FIT *** FIT *** FIT ***
c$$$*     **********************************************************
c$$$
c$$$                              if(chi2.le.0.)goto 666              
c$$$
c$$$*     --------------------------
c$$$*     STORE candidate TRACK INFO
c$$$*     --------------------------
c$$$                              if(ntracks.eq.NTRACKSMAX)then
c$$$                                 
c$$$                                 if(DEBUG)print*,
c$$$     $                 '** warning ** number of candidate tracks '// 
c$$$     $                 ' exceeds vector dimention '
c$$$     $                ,'( ',NTRACKSMAX,' )'
c$$$c                                 good2=.false.
c$$$c                                 goto 880 !fill ntp and go to next event                     
c$$$                                 iflag=1
c$$$                                 return
c$$$                              endif
c$$$                              
c$$$                              ntracks = ntracks + 1
c$$$                              
c$$$c$$$                              ndof=0                                
c$$$                              do ip=1,nplanes
c$$$c$$$                                 ndof=ndof
c$$$c$$$     $                                +int(xgood(ip))
c$$$c$$$     $                                +int(ygood(ip))
c$$$                                 XV_STORE(ip,ntracks)=sngl(xv(ip))
c$$$                                 YV_STORE(ip,ntracks)=sngl(yv(ip))
c$$$                                 ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
c$$$                                 XM_STORE(ip,ntracks)=sngl(xm(ip))
c$$$                                 YM_STORE(ip,ntracks)=sngl(ym(ip))
c$$$                                 ZM_STORE(ip,ntracks)=sngl(zm(ip))
c$$$                                 RESX_STORE(ip,ntracks)=sngl(resx(ip))
c$$$                                 RESY_STORE(ip,ntracks)=sngl(resy(ip))
c$$$                                 XV_STORE(ip,ntracks)=sngl(xv(ip))
c$$$                                 YV_STORE(ip,ntracks)=sngl(yv(ip))
c$$$                                 ZV_STORE(ip,ntracks)=sngl(zv(ip))
c$$$                                 AXV_STORE(ip,ntracks)=sngl(axv(ip))
c$$$                                 AYV_STORE(ip,ntracks)=sngl(ayv(ip))
c$$$                                 XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
c$$$                                 YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
c$$$                                 if(hit_plane(ip).ne.0)then
c$$$                                    CP_STORE(nplanes-ip+1,ntracks)=
c$$$     $                                   cp_match(ip,hit_plane(ip))
c$$$                                 else
c$$$                                    CP_STORE(nplanes-ip+1,ntracks)=0
c$$$                                 endif
c$$$                                 CLS_STORE(nplanes-ip+1,ntracks)=0
c$$$                                 do i=1,5
c$$$                                    AL_STORE(i,ntracks)=sngl(AL(i))
c$$$                                 enddo
c$$$                              enddo
c$$$                              
c$$$c$$$  *                             Number of Degree Of Freedom
c$$$c$$$  ndof=ndof-5                          
c$$$c$$$  *                             reduced chi^2
c$$$c$$$  rchi2=chi2/dble(ndof)
c$$$                              RCHI2_STORE(ntracks)=chi2
c$$$                              
c$$$*     --------------------------------
c$$$*     STORE candidate TRACK INFO - end
c$$$*     --------------------------------
c$$$                              
c$$$ 666                          continue
c$$$                           enddo !end loop on cp in plane 6
c$$$                        enddo   !end loop on cp in plane 5
c$$$                     enddo      !end loop on cp in plane 4
c$$$                  enddo         !end loop on cp in plane 3
c$$$               enddo            !end loop on cp in plane 2
c$$$            enddo               !end loop on cp in plane 1
c$$$            
c$$$ 888        continue
c$$$         enddo                  !end loop on XZ couds
c$$$      enddo                     !end loop on YZ couds
c$$$      
c$$$      if(ntracks.eq.0)then
c$$$         iflag=1
c$$$         return
c$$$      endif
c$$$      
c$$$      if(DEBUG)then
c$$$         print*,'****** TRACK CANDIDATES ***********'
c$$$         print*,'#         R. chi2        RIG'
c$$$         do i=1,ntracks
c$$$            print*,i,' --- ',rchi2_store(i),' --- '
c$$$     $           ,1./abs(AL_STORE(5,i)) 
c$$$         enddo
c$$$         print*,'***********************************'
c$$$      endif
c$$$      
c$$$      
c$$$      return
c$$$      end
c$$$      
c$$$
c$$$***************************************************
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$**************************************************
c$$$
c$$$      subroutine refine_track(ibest)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/common_xyzPAM.f'
c$$$      include '../common/common_mini_2.f'
c$$$      include '../common/common_mech.f'
c$$$      include '../common/momanhough_init.f'
c$$$      include '../common/level1.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$*     flag to chose PFA
c$$$      character*10 PFA
c$$$      common/FINALPFA/PFA
c$$$
c$$$*     =================================================
c$$$*     new estimate of positions using ETA algorithm
c$$$*                          and
c$$$*     search for new couples and single clusters to add
c$$$*     =================================================
c$$$      call track_init
c$$$      do ip=1,nplanes           !loop on planes
c$$$
c$$$*     -------------------------------------------------
c$$$*     If the plane has been already included, it just
c$$$*     computes again the coordinates of the x-y couple
c$$$*     using improved PFAs
c$$$*     -------------------------------------------------
c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
c$$$            
c$$$            id=CP_STORE(nplanes-ip+1,ibest)
c$$$            
c$$$            is=is_cp(id)
c$$$            icp=icp_cp(id)
c$$$            if(ip_cp(id).ne.ip)
c$$$     $           print*,'OKKIO!!'
c$$$     $           ,'id ',id,is,icp
c$$$     $           ,ip_cp(id),ip
c$$$            icx=clx(ip,icp)
c$$$            icy=cly(ip,icp)
c$$$            call xyz_PAM(icx,icy,is,
c$$$c     $           'ETA2','ETA2',
c$$$     $           PFA,PFA,
c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
c$$$     $           AYV_STORE(nplanes-ip+1,ibest))
c$$$c$$$  call xyz_PAM(icx,icy,is,
c$$$c$$$  $              'COG2','COG2',
c$$$c$$$  $              0.,
c$$$c$$$  $              0.)
c$$$            xm(nplanes-ip+1) = xPAM
c$$$            ym(nplanes-ip+1) = yPAM
c$$$            zm(nplanes-ip+1) = zPAM
c$$$            xgood(nplanes-ip+1) = 1
c$$$            ygood(nplanes-ip+1) = 1
c$$$            resx(nplanes-ip+1) = resxPAM
c$$$            resy(nplanes-ip+1) = resyPAM
c$$$
c$$$            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2.
c$$$            
c$$$*     -------------------------------------------------
c$$$*     If the plane has NOT  been already included, 
c$$$*     it tries to include a COUPLE or a single cluster
c$$$*     -------------------------------------------------
c$$$         else                   
c$$$               
c$$$            xgood(nplanes-ip+1)=0
c$$$            ygood(nplanes-ip+1)=0
c$$$               
c$$$*     --------------------------------------------------------------
c$$$*     determine which ladder and sensor are intersected by the track 
c$$$            xP=XV_STORE(nplanes-ip+1,ibest)
c$$$            yP=YV_STORE(nplanes-ip+1,ibest)
c$$$            zP=ZV_STORE(nplanes-ip+1,ibest)
c$$$            call whichsensor(ip,xP,yP,nldt,ist)
c$$$*     if the track hit the plane in a dead area, go to the next plane
c$$$            if(nldt.eq.0.or.ist.eq.0)goto 133 
c$$$*     --------------------------------------------------------------
c$$$
c$$$            if(DEBUG)then
c$$$               print*,
c$$$     $              '------ Plane ',ip,' intersected on LADDER ',nldt
c$$$     $              ,' SENSOR ',ist
c$$$               print*,
c$$$     $              '------ coord: ',XP,YP
c$$$            endif
c$$$            
c$$$*     ===========================================
c$$$*     STEP 1 >>>>>>>  try to include a new couple
c$$$*     ===========================================
c$$$c            if(DEBUG)print*,'>>>> try to include a new couple'
c$$$            distmin=1000000.
c$$$            xmm = 0.
c$$$            ymm = 0.
c$$$            zmm = 0.
c$$$            rxmm = 0.
c$$$            rymm = 0.
c$$$            idm = 0               !ID of the closer couple
c$$$            distance=0.
c$$$            do icp=1,ncp_plane(ip) !loop on couples on plane icp
c$$$               icx=clx(ip,icp)
c$$$               icy=cly(ip,icp)
c$$$               if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
c$$$     $              cl_used(icx).eq.1.or. !or the X cluster is already used 
c$$$     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
c$$$     $              .false.)goto 1188 !then jump to next couple.
c$$$*                                     !Else compute the coordinates
c$$$               call xyz_PAM(icx,icy,ist,
c$$$     $              PFA,PFA,
c$$$c     $              'ETA2','ETA2',
c$$$     $              AXV_STORE(nplanes-ip+1,ibest),
c$$$     $              AYV_STORE(nplanes-ip+1,ibest))
c$$$               
c$$$               distance = distance_to(XP,YP)
c$$$               id=id_cp(ip,icp,ist)
c$$$               if(DEBUG)print*,'( couple ',id
c$$$     $              ,' ) normalized distance ',distance
c$$$               if(distance.lt.distmin)then
c$$$                  xmm = xPAM
c$$$                  ymm = yPAM
c$$$                  zmm = zPAM
c$$$                  rxmm = resxPAM
c$$$                  rymm = resyPAM
c$$$                  distmin = distance
c$$$                  idm = id                  
c$$$                  dedxmm = (dedx(icx)+dedx(icy))/2.
c$$$               endif
c$$$ 1188          continue
c$$$            enddo               !end loop on couples on plane icp
c$$$            if(distmin.le.clinc)then                  
c$$$*              -----------------------------------
c$$$               xm(nplanes-ip+1) = xmm         !<<<
c$$$               ym(nplanes-ip+1) = ymm         !<<<
c$$$               zm(nplanes-ip+1) = zmm         !<<<
c$$$               xgood(nplanes-ip+1) = 1        !<<<
c$$$               ygood(nplanes-ip+1) = 1        !<<<
c$$$               resx(nplanes-ip+1)=rxmm        !<<<
c$$$               resy(nplanes-ip+1)=rymm        !<<<
c$$$               dedxtrk(nplanes-ip+1) = dedxmm !<<<
c$$$*              -----------------------------------
c$$$               CP_STORE(nplanes-ip+1,ibest)=idm      
c$$$               if(DEBUG)print*,'%%%% included couple ',idm
c$$$     $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'
c$$$               goto 133         !next plane
c$$$            endif
c$$$*     ================================================
c$$$*     STEP 2 >>>>>>>  try to include a single cluster
c$$$*                     either from a couple or single
c$$$*     ================================================
c$$$c            if(DEBUG)print*,'>>>> try to include a new cluster'
c$$$            distmin=1000000.
c$$$            xmm_A = 0.          !---------------------------
c$$$            ymm_A = 0.          ! init variables that 
c$$$            zmm_A = 0.          ! define the SINGLET
c$$$            xmm_B = 0.          !
c$$$            ymm_B = 0.          !
c$$$            zmm_B = 0.          !
c$$$            rxmm = 0.           !
c$$$            rymm = 0.           !
c$$$            iclm=0              !---------------------------
c$$$            distance=0.
c$$$
c$$$*----- clusters inside couples -------------------------------------     
c$$$            do icp=1,ncp_plane(ip) !loop on cluster inside couples
c$$$               icx=clx(ip,icp)
c$$$               icy=cly(ip,icp)
c$$$               id=id_cp(ip,icp,ist)
c$$$               if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
c$$$*                                                !jump to the next couple
c$$$*----- try cluster x -----------------------------------------------
c$$$               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used 
c$$$*                                              !jump to the Y cluster 
c$$$               call xyz_PAM(icx,0,ist,
c$$$c     $              'ETA2','ETA2',
c$$$     $              PFA,PFA,
c$$$     $              AXV_STORE(nplanes-ip+1,ibest),0.)               
c$$$               distance = distance_to(XP,YP)
c$$$c     if(DEBUG)print*,'normalized distance ',distance
c$$$               if(DEBUG)print*,'( cl-X ',icx
c$$$     $              ,' in cp ',id,' ) normalized distance ',distance
c$$$               if(distance.lt.distmin)then
c$$$                  xmm_A = xPAM_A
c$$$                  ymm_A = yPAM_A
c$$$                  zmm_A = zPAM_A
c$$$                  xmm_B = xPAM_B
c$$$                  ymm_B = yPAM_B
c$$$                  zmm_B = zPAM_B
c$$$                  rxmm = resxPAM
c$$$                  rymm = resyPAM
c$$$                  distmin = distance
c$$$                  iclm = icx
c$$$                  dedxmm = dedx(icx)
c$$$               endif                  
c$$$11881          continue
c$$$*----- try cluster y -----------------------------------------------
c$$$               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used 
c$$$*                                              !jump to the next couple
c$$$               call xyz_PAM(0,icy,ist,
c$$$c     $              'ETA2','ETA2',
c$$$     $              PFA,PFA,
c$$$     $              0.,AYV_STORE(nplanes-ip+1,ibest))
c$$$               distance = distance_to(XP,YP)
c$$$               if(DEBUG)print*,'( cl-Y ',icy
c$$$     $              ,' in cp ',id,' ) normalized distance ',distance
c$$$               if(distance.lt.distmin)then
c$$$                  xmm_A = xPAM_A
c$$$                  ymm_A = yPAM_A
c$$$                  zmm_A = zPAM_A
c$$$                  xmm_B = xPAM_B
c$$$                  ymm_B = yPAM_B
c$$$                  zmm_B = zPAM_B
c$$$                  rxmm = resxPAM
c$$$                  rymm = resyPAM
c$$$                  distmin = distance
c$$$                  iclm = icy
c$$$                  dedxmm = dedx(icy)
c$$$               endif                  
c$$$11882          continue
c$$$            enddo               !end loop on cluster inside couples
c$$$*----- single clusters -----------------------------------------------     
c$$$            do ic=1,ncls(ip)    !loop on single clusters
c$$$               icl=cls(ip,ic)
c$$$               if(cl_used(icl).eq.1.or.     !if the cluster is already used 
c$$$     $              LADDER(icl).ne.nldt.or. !or the ladder number does not match 
c$$$     $              .false.)goto 18882      !jump to the next singlet
c$$$               if(mod(VIEW(icl),2).eq.0)then!<---- X view
c$$$                  call xyz_PAM(icl,0,ist,
c$$$c     $                 'ETA2','ETA2',
c$$$     $                 PFA,PFA,
c$$$     $                 AXV_STORE(nplanes-ip+1,ibest),0.)
c$$$               else                         !<---- Y view
c$$$                  call xyz_PAM(0,icl,ist,
c$$$c     $                 'ETA2','ETA2',
c$$$     $                 PFA,PFA,
c$$$     $                 0.,AYV_STORE(nplanes-ip+1,ibest))
c$$$               endif
c$$$
c$$$               distance = distance_to(XP,YP)
c$$$               if(DEBUG)print*,'( cl-s ',icl
c$$$     $              ,' ) normalized distance ',distance
c$$$               if(distance.lt.distmin)then
c$$$                  xmm_A = xPAM_A
c$$$                  ymm_A = yPAM_A
c$$$                  zmm_A = zPAM_A
c$$$                  xmm_B = xPAM_B
c$$$                  ymm_B = yPAM_B
c$$$                  zmm_B = zPAM_B
c$$$                  rxmm = resxPAM
c$$$                  rymm = resyPAM
c$$$                  distmin = distance
c$$$                  iclm = icl
c$$$                  dedxmm = dedx(icl)
c$$$               endif                  
c$$$18882          continue
c$$$            enddo               !end loop on single clusters
c$$$
c$$$            if(distmin.le.clinc)then                  
c$$$               
c$$$               CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<     
c$$$*              ----------------------------
c$$$               if(mod(VIEW(iclm),2).eq.0)then
c$$$                  XGOOD(nplanes-ip+1)=1.
c$$$                  resx(nplanes-ip+1)=rxmm
c$$$                  if(DEBUG)print*,'%%%% included X-cl ',iclm
c$$$     $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'
c$$$               else
c$$$                  YGOOD(nplanes-ip+1)=1.
c$$$                  resy(nplanes-ip+1)=rymm
c$$$                  if(DEBUG)print*,'%%%% included Y-cl ',iclm
c$$$     $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'
c$$$               endif
c$$$*              ----------------------------
c$$$               xm_A(nplanes-ip+1) = xmm_A
c$$$               ym_A(nplanes-ip+1) = ymm_A
c$$$               xm_B(nplanes-ip+1) = xmm_B
c$$$               ym_B(nplanes-ip+1) = ymm_B
c$$$               zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
c$$$               dedxtrk(nplanes-ip+1) = dedxmm !<<<
c$$$*              ----------------------------
c$$$            endif
c$$$         endif
c$$$ 133     continue
c$$$      enddo                     !end loop on planes
c$$$      
c$$$      
c$$$      
c$$$      return
c$$$      end
c$$$
c$$$***************************************************
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$*                                                 *
c$$$**************************************************
c$$$
c$$$      subroutine clean_XYclouds(ibest,iflag)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_momanhough.f'
c$$$      include '../common/momanhough_init.f'
c$$$c      include '../common/calib.f'
c$$$c      include '../common/level1.f'
c$$$
c$$$      logical DEBUG
c$$$      common/dbg/DEBUG
c$$$
c$$$
c$$$      do ip=1,nplanes           !loop on planes
c$$$
c$$$         id=CP_STORE(nplanes-ip+1,ibest)
c$$$         icl=CLS_STORE(nplanes-ip+1,ibest)
c$$$         if(id.ne.0.or.icl.ne.0)then               
c$$$            if(id.ne.0)then
c$$$               iclx=clx(ip,icp_cp(id))
c$$$               icly=cly(ip,icp_cp(id))
c$$$               cl_used(iclx)=1  !tag used clusters
c$$$               cl_used(icly)=1  !tag used clusters
c$$$            elseif(icl.ne.0)then
c$$$               cl_used(icl)=1   !tag used clusters
c$$$            endif
c$$$            
c$$$c               if(DEBUG)then
c$$$c                  print*,ip,' <<< ',id
c$$$c               endif
c$$$*     ----------------------------- 
c$$$*     remove the couple from clouds
c$$$*     remove also vitual couples containing the 
c$$$*     selected clusters
c$$$*     ----------------------------- 
c$$$            do icp=1,ncp_plane(ip)
c$$$               if(
c$$$     $              clx(ip,icp).eq.iclx
c$$$     $              .or.
c$$$     $              clx(ip,icp).eq.icl
c$$$     $              .or.
c$$$     $              cly(ip,icp).eq.icly
c$$$     $              .or.
c$$$     $              cly(ip,icp).eq.icl
c$$$     $              )then
c$$$                  id=id_cp(ip,icp,1)
c$$$                  if(DEBUG)then
c$$$                     print*,ip,' <<< cp ',id
c$$$     $                    ,' ( cl-x '
c$$$     $                    ,clx(ip,icp)
c$$$     $                    ,' cl-y '
c$$$     $                    ,cly(ip,icp),' ) --> removed'
c$$$                  endif
c$$$*     ----------------------------- 
c$$$*     remove the couple from clouds
c$$$                  do iyz=1,nclouds_yz
c$$$                     if(cpcloud_yz(iyz,abs(id)).ne.0)then
c$$$                        ptcloud_yz(iyz)=ptcloud_yz(iyz)-1
c$$$                        cpcloud_yz(iyz,abs(id))=0
c$$$                     endif
c$$$                  enddo
c$$$                  do ixz=1,nclouds_xz
c$$$                     if(cpcloud_xz(ixz,abs(id)).ne.0)then
c$$$                        ptcloud_xz(ixz)=ptcloud_xz(ixz)-1
c$$$                        cpcloud_xz(ixz,abs(id))=0
c$$$                     endif
c$$$                  enddo                     
c$$$*     ----------------------------- 
c$$$               endif
c$$$            enddo
c$$$            
c$$$         endif               
c$$$      enddo                     !end loop on planes
c$$$      
c$$$      return
c$$$      end
c$$$
c$$$
c$$$
c$$$
c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
c$$$      real function fbad_cog(ncog,ic)
c$$$
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/level1.f'
c$$$      include '../common/calib.f'
c$$$
c$$$*     --> signal of the central strip
c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
c$$$
c$$$*     signal of adjacent strips
c$$$*     --> left
c$$$      sl1 = 0                  !left 1
c$$$      if(
c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
c$$$     $     )
c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
c$$$
c$$$      sl2 = 0                  !left 2
c$$$      if(
c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
c$$$     $     )
c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
c$$$
c$$$*     --> right
c$$$      sr1 = 0                  !right 1
c$$$      if(
c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
c$$$     $     .or.
c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
c$$$     $     )
c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
c$$$
c$$$      sr2 = 0                  !right 2
c$$$      if(
c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
c$$$     $     .or.
c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
c$$$     $     )
c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
c$$$
c$$$
c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
c$$$         f  = 4.
c$$$         si = 8.4
c$$$      else                              !X-view
c$$$         f  = 6.
c$$$         si = 3.9
c$$$      endif
c$$$
c$$$      fbad_cog = 1.
c$$$      f0 = 1
c$$$      f1 = 1
c$$$      f2 = 1
c$$$      f3 = 1   
c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
c$$$         
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
c$$$
c$$$         if(ncog.eq.2.and.sl1.ne.0)then
c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
c$$$            fbad_cog = 1.
c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
c$$$            fbad_cog = 1.
c$$$         else
c$$$            fbad_cog = 1.
c$$$         endif
c$$$         
c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
c$$$
c$$$
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
c$$$
c$$$         if(ncog.eq.2.and.sr1.ne.0)then
c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
c$$$            fbad_cog = 1.
c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
c$$$            fbad_cog = 1.
c$$$         else
c$$$            fbad_cog = 1.
c$$$         endif
c$$$
c$$$      endif
c$$$
c$$$      fbad_cog = sqrt(fbad_cog)
c$$$
c$$$      return
c$$$      end