src/F77/analysissubroutines.f

************************************************************
*     The following subroutines
*     - track_finding  >> hough transform
*     - track_fitting  >> bob golden fitting
*     all the procedures to create LEVEL2 data, starting from LEVEL1 data.
*
*     
*     
*     (This subroutine and all the dependent subroutines
*      will be included in the flight software)
************************************************************
      subroutine track_finding(iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'common_mech.f'
      include 'common_xyzPAM.f'
      include 'common_mini_2.f'
      include 'calib.f'
      include 'level2.f'


c      print*,'======================================================'
c$$$      do ic=1,NCLSTR1
c$$$         if(.false.
c$$$     $        .or.nsatstrips(ic).gt.0
c$$$c     $        .or.nbadstrips(0,ic).gt.0
c$$$c     $        .or.nbadstrips(4,ic).gt.0
c$$$c     $        .or.nbadstrips(3,ic).gt.0
c$$$     $        .or..false.)then
c$$$            print*,'--- cl-',ic,' ------------------------'
c$$$            istart = INDSTART(IC)
c$$$            istop  = TOTCLLENGTH
c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
c$$$            print*,'view ',VIEW(ic)
c$$$            print*,'maxs ',MAXS(ic)
c$$$            print*,'COG4 ',cog(4,ic)
c$$$            ff = fbad_cog(4,ic)
c$$$            print*,'fbad ',ff
c$$$            print*,(CLBAD(i),i=istart,istop)
c$$$            bb=nbadstrips(0,ic)
c$$$            print*,'#BAD (tot)',bb
c$$$            bb=nbadstrips(4,ic)
c$$$            print*,'#BAD (4)',bb
c$$$            bb=nbadstrips(3,ic)
c$$$            print*,'#BAD (3)',bb
c$$$            ss=nsatstrips(ic)
c$$$            print*,'#saturated ',ss
c$$$         endif
c$$$      enddo
      
*-------------------------------------------------------------------------------
*     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
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------


      call cl_to_couples(iflag)
      if(iflag.eq.1)then        !bad event
         goto 880               !go to next event
      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
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------


      call cp_to_doubtrip(iflag)
      if(iflag.eq.1)then        !bad event
         goto 880               !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     xz2_av,alfaxz3_av 
*     $     ,ptcloud_xz,tr_cloud,cpcloud_xz           
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
*     count number of hit planes
      planehit=0                
      do np=1,nplanes           
        if(ncp_plane(np).ne.0)then  
          planehit=planehit+1   
        endif                   
      enddo                     
      if(planehit.lt.3) goto 880 ! exit              
      
      nptxz_min=x_min_start               
      nplxz_min=x_min_start               
             
      nptyz_min=y_min_start               
      nplyz_min=y_min_start               

      cutdistyz=cutystart       
      cutdistxz=cutxstart       

 878  continue
      call doub_to_YZcloud(iflag)
      if(iflag.eq.1)then        !bad event
         goto 880               !fill ntp and go to next event            
      endif 
      if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then 
        if(cutdistyz.lt.maxcuty/2)then
          cutdistyz=cutdistyz+cutystep 
        else
          cutdistyz=cutdistyz+(3*cutystep) 
        endif
        goto 878                
      endif                     

      if(planehit.eq.3) goto 881          
      
 879  continue  
      call trip_to_XZcloud(iflag)
      if(iflag.eq.1)then        !bad event
         goto 880               !fill ntp and go to next event            
      endif
                          
      if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then 
        cutdistxz=cutdistxz+cutxstep 
        goto 879                
      endif                     

     
 881  continue   
*     if there is at least three planes on the Y view decreases cuts on X view
      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
     $     nplxz_min.ne.y_min_start)then 
        nptxz_min=x_min_step    
        nplxz_min=x_min_start-x_min_step              
        goto 879                
      endif                     
        
 880  return
      end

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

      
      subroutine track_fitting(iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'common_mech.f'
      include 'common_xyzPAM.f'
      include 'common_mini_2.f'
      include 'calib.f'
      include 'level2.f'

c      include 'momanhough_init.f'
      
      logical FIMAGE            !
      real*8 AL_GUESS(5)

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

c         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 -------------
c$$$         rchi2best=1000000000.
c$$$         do i=1,ntracks
c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
c$$$     $         RCHI2_STORE(i).gt.0)then
c$$$               ibest=i
c$$$               rchi2best=RCHI2_STORE(i)
c$$$            endif
c$$$         enddo
c$$$         if(ibest.eq.0)goto 880 !>> no good candidates 

*     -------------------------------------------------------
*     order track-candidates according to:
*     1st) decreasing n.points
*     2nd) increasing chi**2 
*     -------------------------------------------------------
         rchi2best=1000000000.
         ndofbest=0            
         do i=1,ntracks
           ndof=0              
           do ii=1,nplanes     
             ndof=ndof         
     $            +int(xgood_store(ii,i)) 
     $            +int(ygood_store(ii,i))
           enddo               
           if(ndof.gt.ndofbest)then 
             ibest=i
             rchi2best=RCHI2_STORE(i)
             ndofbest=ndof    
           elseif(ndof.eq.ndofbest)then
             if(RCHI2_STORE(i).lt.rchi2best.and.
     $            RCHI2_STORE(i).gt.0)then
               ibest=i
               rchi2best=RCHI2_STORE(i)
               ndofbest=ndof   
             endif            
           endif
         enddo

c$$$         rchi2best=1000000000.
c$$$         ndofbest=0             !(1)
c$$$         do i=1,ntracks
c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
c$$$     $          RCHI2_STORE(i).gt.0)then
c$$$             ndof=0             !(1)
c$$$             do ii=1,nplanes    !(1)
c$$$               ndof=ndof        !(1)
c$$$     $              +int(xgood_store(ii,i)) !(1)
c$$$     $              +int(ygood_store(ii,i)) !(1)
c$$$             enddo              !(1)
c$$$             if(ndof.ge.ndofbest)then !(1)
c$$$               ibest=i
c$$$               rchi2best=RCHI2_STORE(i)
c$$$               ndofbest=ndof    !(1)
c$$$             endif              !(1)
c$$$           endif
c$$$         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 candidate 
            icand=ibest            
         else                   !processing image
            icand=iimage
            iimage=0
         endif
         if(icand.eq.0)then
            if(VERBOSE)then
               print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
     $              ,ibest,iimage
            endif
            return
         endif

*     *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
         call refine_track(icand)
*     *-*-*-*-*-*-*-*-*-*-*-*-*-*-*

*     **********************************************************
*     ************************** FIT *** FIT *** FIT *** FIT ***
*     **********************************************************
         call guess()
         do i=1,5
            AL_GUESS(i)=AL(i)
         enddo
c         print*,'## guess: ',al

         do i=1,5
            AL(i)=dble(AL_STORE(i,icand))            
         enddo
         
         IDCAND = icand         !fitted track-candidate
         ifail=0                !error flag in chi2 computation
         jstep=0                !# minimization steps

         iprint=0
c         if(DEBUG)iprint=1
         if(VERBOSE)iprint=1
         if(DEBUG)iprint=2
         call mini2(jstep,ifail,iprint)
         if(ifail.ne.0) then
            if(VERBOSE)then
               print *,
     $              '*** MINIMIZATION FAILURE *** (after refinement) '
     $              ,iev

c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
c$$$               print*,'result:   ',(al(i),i=1,5)
c$$$               print*,'xgood ',xgood
c$$$               print*,'ygood ',ygood
c$$$               print*,'----------------------------------------------'
            endif
c            chi2=-chi2 
         endif 
         
         if(DEBUG)then
            print*,'----------------------------- improved track coord'
22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.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(verbose)
     $           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,iflag)
         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
*     **********************************************


 880     return
      end


************************************************************
************************************************************
************************************************************
************************************************************
*
*     This routine provides the coordinates (in cm) in the PAMELA reference system:
*       - of the point associated with a COUPLE ---> (xPAM,yPAM,zPAM)
*       - of the extremes of the segment
*         associated with a SINGLET ---------------> (xPAM_A,yPAM_A,zPAM_A) 
*                                               ---> (xPAM_B,yPAM_B,zPAM_B)
*
*     It also assigns the spatial resolution to the evaluated coordinates,
*     as a function (in principle) of the multiplicity, the angle, the PFA etc...
*
*     
*     To call the routine you must pass the arguments:
*     icx    - ID of cluster x
*     icy    - ID of cluster y
*     sensor - sensor (1,2)
*     PFAx   - Position Finding Algorithm in x (COG2,ETA2,...)
*     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
*     angx   - Projected angle in x
*     angy   - Projected angle in y
*     bfx    - x-component of magnetci field
*     bfy    - y-component of magnetic field
*
*     --------- COUPLES -------------------------------------------------------
*     The couple defines a point in the space. 
*     The coordinates of the point are evaluated as follows:
*     1 - the corrected coordinates relative to the sensor are evaluated
*         according to the chosen PFA --> (xi,yi,0)
*     2 - coordinates are rotated and traslated, according to the aligmnet 
*         parameters, and expressed in the reference system of the mechanical 
*         sensor --> (xrt,yrt,zrt)
*     3 - coordinates are finally converted to the PAMELA reference system 
*         --> (xPAM,yPAM,zPAM)
*
*     --------- SINGLETS -------------------------------------------------------
*     Since a coordinate is missing, the singlet defines not a point 
*     in the space but a segment AB (parallel to the strips). 
*     In this case the routine returns the coordinates in the PAMELA reference
*     system of the two extremes A and B of the segment:
*         --> (xPAM_A,yPAM_A,zPAM_A)
*         --> (xPAM_B,yPAM_B,zPAM_B)
*     
*     ==========================================================
*
*     The output of the routine is stored in the commons:
*
*      double precision xPAM,yPAM,zPAM
*      common/coord_xyz_PAM/xPAM,yPAM,zPAM
*     
*      double precision xPAM_A,yPAM_A,zPAM_A
*      double precision xPAM_B,yPAM_B,zPAM_B
*      common/coord_AB_PAM/xPAM_A,yPAM_A,zPAM_A,xPAM_B,yPAM_B,zPAM_B
*
*      double precision resxPAM,resyPAM
*      common/resolution_PAM/resxPAM,resyPAM
*
*     (in file common_xyzPAM.f)
*
*

      subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)

      
      include 'commontracker.f'
      include 'level1.f'
      include 'calib.f'
      include 'common_align.f'
      include 'common_mech.f'
      include 'common_xyzPAM.f'

      integer icx,icy           !X-Y cluster ID
      integer sensor
      integer viewx,viewy
      character*4 PFAx,PFAy     !PFA to be used
      real ax,ay                !X-Y geometric angle
      real angx,angy            !X-Y effective angle
      real bfx,bfy              !X-Y b-field components

      real stripx,stripy

      double precision xrt,yrt,zrt
      double precision xrt_A,yrt_A,zrt_A
      double precision xrt_B,yrt_B,zrt_B
      

      parameter (ndivx=30)


c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
      
      resxPAM = 0
      resyPAM = 0

      xPAM = 0.
      yPAM = 0.
      zPAM = 0.
      xPAM_A = 0.
      yPAM_A = 0.
      zPAM_A = 0.
      xPAM_B = 0.
      yPAM_B = 0.
      zPAM_B = 0.
c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy

*     -----------------
*     CLUSTER X
*     -----------------

      if(icx.ne.0)then

         viewx   = VIEW(icx)
         nldx    = nld(MAXS(icx),VIEW(icx))
         nplx    = npl(VIEW(icx))
         resxPAM = RESXAV 
         stripx  = float(MAXS(icx))
*        --------------------------
*        magnetic-field corrections
*        --------------------------
         angtemp  = ax
         bfytemp  = bfy
*        /////////////////////////////////
*        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!
*        *grvzkkjsdgjhhhgngbn###>:(
*        /////////////////////////////////
c         if(nplx.eq.6) angtemp = -1. * ax
c         if(nplx.eq.6) bfytemp = -1. * bfy
         if(viewx.eq.12) angtemp = -1. * ax
         if(viewx.eq.12) bfytemp = -1. * bfy
         tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
         angx     = 180.*atan(tgtemp)/acos(-1.)
         stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
c$$$         print*,nplx,ax,bfy/10.
c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
c$$$         print*,'========================'
c$$$         if(bfy.ne.0.)print*,viewx,'-x- '
c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ
*        --------------------------

c$$$         print*,'--- x-cl ---'
c$$$         istart = INDSTART(ICX)
c$$$         istop  = TOTCLLENGTH
c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
c$$$         print*,(CLSIGNAL(i),i=istart,istop)
c$$$         print*,INDMAX(icx)
c$$$         print*,cog(4,icx)
c$$$         print*,fbad_cog(4,icx)
         

         if(PFAx.eq.'COG1')then 

            stripx  = stripx
            resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM

         elseif(PFAx.eq.'COG2')then

            stripx  = stripx + cog(2,icx)            
            resxPAM = risx_cog(abs(angx))!TEMPORANEO               
            resxPAM = resxPAM*fbad_cog(2,icx)

         elseif(PFAx.eq.'COG3')then

            stripx  = stripx + cog(3,icx)            
            resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
            resxPAM = resxPAM*fbad_cog(3,icx)

         elseif(PFAx.eq.'COG4')then

            stripx  = stripx + cog(4,icx)            
            resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
            resxPAM = resxPAM*fbad_cog(4,icx)

         elseif(PFAx.eq.'ETA2')then

            stripx  = stripx + pfaeta2(icx,angx)           
            resxPAM = risx_eta2(abs(angx))
            resxPAM = resxPAM*fbad_cog(2,icx)
            if(DEBUG.and.fbad_cog(2,icx).ne.1)
     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)

         elseif(PFAx.eq.'ETA3')then                        

            stripx  = stripx + pfaeta3(icx,angx)           
            resxPAM = risx_eta3(abs(angx))                       
            resxPAM = resxPAM*fbad_cog(3,icx)               
            if(DEBUG.and.fbad_cog(3,icx).ne.1)             
     $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)

         elseif(PFAx.eq.'ETA4')then                         

            stripx  = stripx + pfaeta4(icx,angx)            
            resxPAM = risx_eta4(abs(angx))                       
            resxPAM = resxPAM*fbad_cog(4,icx)               
            if(DEBUG.and.fbad_cog(4,icx).ne.1)              
     $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)

         elseif(PFAx.eq.'ETA')then   

            stripx  = stripx + pfaeta(icx,angx)             
            resxPAM = ris_eta(icx,angx)                     
            resxPAM = resxPAM*fbad_eta(icx,angx)             
            if(DEBUG.and.fbad_cog(2,icx).ne.1)              
     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)

         elseif(PFAx.eq.'COG')then          

            stripx  = stripx + cog(0,icx)             
            resxPAM = risx_cog(abs(angx))                     
            resxPAM = resxPAM*fbad_cog(0,icx)

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


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

      endif
      
*     ----------------- 
*     CLUSTER Y
*     -----------------

      if(icy.ne.0)then

         viewy = VIEW(icy)
         nldy = nld(MAXS(icy),VIEW(icy))
         nply = npl(VIEW(icy))
         resyPAM = RESYAV 
         stripy = float(MAXS(icy))

         if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
            if(DEBUG) then
               print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
     $              ,icx,icy
            endif
            goto 100
         endif
*        --------------------------
*        magnetic-field corrections
*        --------------------------
         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001         
         angy    = 180.*atan(tgtemp)/acos(-1.)
         stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
*        --------------------------
         
c$$$         print*,'--- y-cl ---'
c$$$         istart = INDSTART(ICY)
c$$$         istop  = TOTCLLENGTH
c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
c$$$         print*,(CLSIGNAL(i),i=istart,istop)
c$$$         print*,INDMAX(icy)
c$$$         print*,cog(4,icy)
c$$$         print*,fbad_cog(4,icy)

         if(PFAy.eq.'COG1')then

            stripy  = stripy    
            resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  

         elseif(PFAy.eq.'COG2')then

            stripy  = stripy + cog(2,icy)
            resyPAM = risy_cog(abs(angy))!TEMPORANEO 
            resyPAM = resyPAM*fbad_cog(2,icy)

         elseif(PFAy.eq.'COG3')then

            stripy  = stripy + cog(3,icy)
            resyPAM = risy_cog(abs(angy))!TEMPORANEO 
            resyPAM = resyPAM*fbad_cog(3,icy)

         elseif(PFAy.eq.'COG4')then

            stripy  = stripy + cog(4,icy)
            resyPAM = risy_cog(abs(angy))!TEMPORANEO 
            resyPAM = resyPAM*fbad_cog(4,icy)

         elseif(PFAy.eq.'ETA2')then

            stripy  = stripy + pfaeta2(icy,angy)          
            resyPAM = risy_eta2(abs(angy))              
            resyPAM = resyPAM*fbad_cog(2,icy)
            if(DEBUG.and.fbad_cog(2,icy).ne.1)
     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)

         elseif(PFAy.eq.'ETA3')then                      

            stripy  = stripy + pfaeta3(icy,angy) 
            resyPAM = resyPAM*fbad_cog(3,icy)  
            if(DEBUG.and.fbad_cog(3,icy).ne.1) 
     $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)

         elseif(PFAy.eq.'ETA4')then  

            stripy  = stripy + pfaeta4(icy,angy) 
            resyPAM = resyPAM*fbad_cog(4,icy)
            if(DEBUG.and.fbad_cog(4,icy).ne.1)
     $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)

         elseif(PFAy.eq.'ETA')then 

            stripy  = stripy + pfaeta(icy,angy)
            resyPAM = ris_eta(icy,angy)  
            resyPAM = resyPAM*fbad_eta(icy,angy)
            if(DEBUG.and.fbad_cog(2,icy).ne.1) 
     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)

         elseif(PFAy.eq.'COG')then

            stripy  = stripy + cog(0,icy)            
            resyPAM = risy_cog(abs(angy))
            resyPAM = resyPAM*fbad_cog(0,icy)

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


*     ======================================
*     temporary patch for saturated clusters
*     ======================================
         if( nsatstrips(icy).gt.0 )then
            stripy  = stripy + cog(4,icy)             
            resyPAM = pitchY*1e-4/sqrt(12.)
            cc=cog(4,icy)
c$$$            print*,icy,' *** ',cc
c$$$            print*,icy,' *** ',resyPAM
         endif


      endif

c$$$      print*,'## stripx,stripy ',stripx,stripy

c===========================================================
C     COUPLE
C===========================================================
      if(icx.ne.0.and.icy.ne.0)then 

c------------------------------------------------------------------------
c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
c------------------------------------------------------------------------
         if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
     $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
            if(DEBUG) then
               print*,'xyz_PAM (couple):', 
     $              ' WARNING: false X strip: strip ',stripx
            endif
         endif
         xi = acoordsi(stripx,viewx)
         yi = acoordsi(stripy,viewy)
         zi = 0.
         

c------------------------------------------------------------------------
c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c------------------------------------------------------------------------
c     N.B. I convert angles from microradiants to radiants

         xrt = xi
     $        - omega(nplx,nldx,sensor)*yi 
     $        + gamma(nplx,nldx,sensor)*zi
     $        + dx(nplx,nldx,sensor)
         
         yrt = omega(nplx,nldx,sensor)*xi
     $        + yi
     $        - beta(nplx,nldx,sensor)*zi
     $        + dy(nplx,nldx,sensor)
         
         zrt = -gamma(nplx,nldx,sensor)*xi
     $        + beta(nplx,nldx,sensor)*yi
     $        + zi
     $        + dz(nplx,nldx,sensor)
         
c      xrt = xi
c      yrt = yi
c      zrt = zi
    
c------------------------------------------------------------------------
c     (xPAM,yPAM,zPAM) = measured coordinates (in cm) 
c                        in PAMELA reference system
c------------------------------------------------------------------------

         xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
         yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
         zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4

         xPAM_A = 0.
         yPAM_A = 0.
         zPAM_A = 0.

         xPAM_B = 0.
         yPAM_B = 0.
         zPAM_B = 0.

      elseif(
     $        (icx.ne.0.and.icy.eq.0).or.
     $        (icx.eq.0.and.icy.ne.0).or.
     $        .false.
     $        )then

c------------------------------------------------------------------------
c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
c------------------------------------------------------------------------

         if(icy.ne.0)then
c===========================================================
C     Y-SINGLET
C===========================================================
            nplx = nply
            nldx = nldy
            viewx = viewy + 1

            yi   = acoordsi(stripy,viewy)

            xi_A = edgeY_d - SiDimX/2
            yi_A = yi
            zi_A = 0.
            
            xi_B = SiDimX/2 - edgeY_u
            yi_B = yi
            zi_B = 0.

c            print*,'Y-cl ',icy,stripy,' --> ',yi
c            print*,xi_A,' <--> ',xi_B
            
         elseif(icx.ne.0)then
c===========================================================
C     X-SINGLET
C===========================================================

            nply = nplx
            nldy = nldx
            viewy = viewx - 1

c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
            if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
     $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
               if(DEBUG) then
                  print*,'xyz_PAM (X-singlet):', 
     $                 ' WARNING: false X strip: strip ',stripx
               endif
            endif
            xi   = acoordsi(stripx,viewx)

            xi_A = xi
            yi_A = edgeX_d - SiDimY/2
            zi_A = 0.
            
            xi_B = xi
            yi_B = SiDimY/2 - edgeX_u
            zi_B = 0.
            
            if(viewy.eq.11)then
               yi = yi_A
               yi_A = yi_B
               yi_B = yi
            endif

c            print*,'X-cl ',icx,stripx,' --> ',xi
c            print*,yi_A,' <--> ',yi_B

         else
            if(DEBUG) then
               print *,'routine xyz_PAM ---> not properly used !!!'
               print *,'icx = ',icx
               print *,'icy = ',icy
            endif
            goto 100
            
         endif
c------------------------------------------------------------------------
c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c------------------------------------------------------------------------
c     N.B. I convert angles from microradiants to radiants

         xrt_A = xi_A
     $        - omega(nplx,nldx,sensor)*yi_A
     $        + gamma(nplx,nldx,sensor)*zi_A
     $        + dx(nplx,nldx,sensor)
         
         yrt_A = omega(nplx,nldx,sensor)*xi_A
     $        + yi_A
     $        - beta(nplx,nldx,sensor)*zi_A
     $        + dy(nplx,nldx,sensor)
         
         zrt_A = -gamma(nplx,nldx,sensor)*xi_A
     $        + beta(nplx,nldx,sensor)*yi_A
     $        + zi_A
     $        + dz(nplx,nldx,sensor)

         xrt_B = xi_B
     $        - omega(nplx,nldx,sensor)*yi_B 
     $        + gamma(nplx,nldx,sensor)*zi_B
     $        + dx(nplx,nldx,sensor)
         
         yrt_B = omega(nplx,nldx,sensor)*xi_B
     $        + yi_B
     $        - beta(nplx,nldx,sensor)*zi_B
     $        + dy(nplx,nldx,sensor)
         
         zrt_B = -gamma(nplx,nldx,sensor)*xi_B
     $        + beta(nplx,nldx,sensor)*yi_B
     $        + zi_B
     $        + dz(nplx,nldx,sensor)

         
c      xrt = xi
c      yrt = yi
c      zrt = zi
    
c------------------------------------------------------------------------
c     (xPAM,yPAM,zPAM) = measured coordinates (in cm) 
c                        in PAMELA reference system
c------------------------------------------------------------------------

         xPAM = 0.
         yPAM = 0.
         zPAM = 0.

         xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
         yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
         zPAM_A = ( zrt_A + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4

         xPAM_B = dcoord(xrt_B,viewx,nldx,sensor) / 1.d4
         yPAM_B = dcoord(yrt_B,viewy,nldy,sensor) / 1.d4
         zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
         

c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'

      else
         if(DEBUG) then
            print *,'routine xyz_PAM ---> not properly used !!!'
            print *,'icx = ',icx
            print *,'icy = ',icy
         endif
      endif
         

c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B

 100  continue
      end

************************************************************************
*     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
*     (done to be called from c/c++)
************************************************************************

      subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_mini_2.f'
      include 'common_xyzPAM.f'
      include 'common_mech.f'
      include 'calib.f'
      
*     flag to chose PFA
c$$$      character*10 PFA
c$$$      common/FINALPFA/PFA

      integer icx,icy           !X-Y cluster ID
      integer sensor
      character*4 PFAx,PFAy     !PFA to be used
      real ax,ay                !X-Y geometric angle
      real bfx,bfy              !X-Y b-field components

      ipx=0
      ipy=0      
      
c$$$      PFAx = 'COG4'!PFA
c$$$      PFAy = 'COG4'!PFA
      
      call idtoc(pfaid,PFAx)
      call idtoc(pfaid,PFAy)

      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)

c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
      
      if(icx.ne.0.and.icy.ne.0)then

         ipx=npl(VIEW(icx))
         ipy=npl(VIEW(icy))
         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
     $        ,' does not belong to the correct plane: ',ip,ipx,ipy

         xgood(ip) = 1.
         ygood(ip) = 1.
         resx(ip)  = resxPAM
         resy(ip)  = resyPAM

         xm(ip) = xPAM
         ym(ip) = yPAM
         zm(ip) = zPAM
         xm_A(ip) = 0.
         ym_A(ip) = 0.
         xm_B(ip) = 0.
         ym_B(ip) = 0.

c         zv(ip) = zPAM

      elseif(icx.eq.0.and.icy.ne.0)then

         ipy=npl(VIEW(icy))
         if((nplanes-ipy+1).ne.ip)
     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
     $        ,' does not belong to the correct plane: ',ip,ipx,ipy

         xgood(ip) = 0.
         ygood(ip) = 1.
         resx(ip)  = 1000.
         resy(ip)  = resyPAM

         xm(ip) = -100.
         ym(ip) = -100.
         zm(ip) = (zPAM_A+zPAM_B)/2.
         xm_A(ip) = xPAM_A
         ym_A(ip) = yPAM_A
         xm_B(ip) = xPAM_B
         ym_B(ip) = yPAM_B

c         zv(ip) = (zPAM_A+zPAM_B)/2.
         
      elseif(icx.ne.0.and.icy.eq.0)then

         ipx=npl(VIEW(icx))
         if((nplanes-ipx+1).ne.ip)
     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
     $        ,' does not belong to the correct plane: ',ip,ipx,ipy

         xgood(ip) = 1.
         ygood(ip) = 0.
         resx(ip)  = resxPAM
         resy(ip)  = 1000.

         xm(ip) = -100.
         ym(ip) = -100.
         zm(ip) = (zPAM_A+zPAM_B)/2.
         xm_A(ip) = xPAM_A
         ym_A(ip) = yPAM_A
         xm_B(ip) = xPAM_B
         ym_B(ip) = yPAM_B
         
c         zv(ip) = (zPAM_A+zPAM_B)/2.

      else

         il = 2
         if(lad.ne.0)il=lad
         is = 1
         if(sensor.ne.0)is=sensor
c         print*,nplanes-ip+1,il,is

         xgood(ip) = 0.
         ygood(ip) = 0.
         resx(ip)  = 1000.
         resy(ip)  = 1000.

         xm(ip) = -100.
         ym(ip) = -100.           
         zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
         xm_A(ip) = 0.
         ym_A(ip) = 0.
         xm_B(ip) = 0.
         ym_B(ip) = 0.

c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4

      endif

      if(DEBUG)then
c         print*,'----------------------------- track coord'
22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
         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)
c$$$         print*,'-----------------------------'
      endif
      end

********************************************************************************
********************************************************************************
********************************************************************************
*
*     The function distance_to(XP,YP) should be used after 
*     a call to the xyz_PAM routine and it evaluate the 
*     NORMALIZED distance (PROJECTED on the XY plane) between 
*     the point  (XP,YP), argument of the function, 
*     and:
*     
*     - the point (xPAM,yPAM,zPAM), in the case of a COUPLE
*     or 
*     - the segment (xPAM_A,yPAM_A,zPAM_A)-(xPAM_B,yPAM_B,zPAM_B), 
*       in the case of a SINGLET.
*
*     ( The routine xyz_PAM fills the common defined in "common_xyzPAM.f",
*      which stores the coordinates of the couple/singlet )
*     
********************************************************************************

      real function distance_to(XPP,YPP)

      include 'common_xyzPAM.f'

*     -----------------------------------
*     it computes the normalized distance 
*     ( i.e. distance/resolution )
*     -----------------------------------

      double precision distance,RE
      double precision BETA,ALFA,xmi,ymi

*     ----------------------
      if (
     +     xPAM.eq.0.and.
     +     yPAM.eq.0.and.
     +     zPAM.eq.0.and.
     +     xPAM_A.ne.0.and.
     +     yPAM_A.ne.0.and.
     +     zPAM_A.ne.0.and.
     +     xPAM_B.ne.0.and.
     +     yPAM_B.ne.0.and.
     +     zPAM_B.ne.0.and.
     +     .true.)then
*     -----------------------
*     DISTANCE TO --- SINGLET
*     -----------------------
         if(abs(sngl(xPAM_B-xPAM_A)).lt.abs(sngl(yPAM_B-yPAM_A)))then
*        |||---------- X CLUSTER
            
            BETA = (xPAM_B-xPAM_A)/(yPAM_B-yPAM_A)
            ALFA = xPAM_A - BETA * yPAM_A

            ymi = ( YPP + BETA*XPP - BETA*ALFA )/(1+BETA**2)
            if(ymi.lt.dmin1(yPAM_A,yPAM_B))ymi=dmin1(yPAM_A,yPAM_B)
            if(ymi.gt.dmax1(yPAM_A,yPAM_B))ymi=dmax1(yPAM_A,yPAM_B)
            xmi = ALFA + BETA * ymi
            RE = resxPAM
                        
         else
*        |||---------- Y CLUSTER

            BETA = (yPAM_B-yPAM_A)/(xPAM_B-xPAM_A)
            ALFA = yPAM_A - BETA * xPAM_A

            xmi = ( XPP + BETA*YPP - BETA*ALFA )/(1+BETA**2)
            if(xmi.lt.dmin1(xPAM_A,xPAM_B))xmi=dmin1(xPAM_A,xPAM_B)
            if(xmi.gt.dmax1(xPAM_A,xPAM_B))xmi=dmax1(xPAM_A,xPAM_B)
            ymi = ALFA + BETA * xmi
            RE = resyPAM
            
         endif         

         distance=
     $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
         distance=dsqrt(distance)                     

c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
c$$$     $        ,' --- distance_to --- ',xpp,ypp
c$$$         print*,' resolution ',re

         
*     ----------------------
      elseif(
     +     xPAM.ne.0.and.
     +     yPAM.ne.0.and.
     +     zPAM.ne.0.and.
     +     xPAM_A.eq.0.and.
     +     yPAM_A.eq.0.and.
     +     zPAM_A.eq.0.and.
     +     xPAM_B.eq.0.and.
     +     yPAM_B.eq.0.and.
     +     zPAM_B.eq.0.and.
     +     .true.)then
*     ----------------------
*     DISTANCE TO --- COUPLE
*     ----------------------
         
         distance=
     $       ((xPAM-XPP))**2 !QUIQUI
     $       +
     $       ((yPAM-YPP))**2
c$$$     $        ((xPAM-XPP)/resxPAM)**2
c$$$     $        +
c$$$     $        ((yPAM-YPP)/resyPAM)**2
         distance=dsqrt(distance)                     

c$$$         print*,xPAM,yPAM,zPAM
c$$$     $        ,' --- distance_to --- ',xpp,ypp
c$$$         print*,' resolution ',resxPAM,resyPAM
         
      else
         
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
      endif   

      distance_to = sngl(distance)

      return
      end

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

      subroutine whichsensor(nplPAM,xPAM,yPAM,ladder,sensor)
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *   
*     Given the plane (1-6 from BOTTOM to TOP!!) and the (xPAM,yPAM) 
*     coordinates (in the PAMELA reference system), it returns 
*     the ladder and the sensor which the point belongs to.
*
*     The method to assign a point to a sensor consists in 
*     - calculating the sum of the distances between the point
*     and the sensor edges
*     - requiring that it is less-equal than (SiDimX+SiDimY)
*
*     NB -- SiDimX and SiDimY are not the dimentions of the SENSITIVE volume
*           but of the whole silicon sensor
*     
*     CONVENTION:
*     - sensor 1 is the one closest to the hybrid
*     - ladder 1 is the first to be read out (strips from 1 to 1024)
*
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
      include 'commontracker.f'
      include 'common_align.f'

      integer ladder,sensor,viewx,viewy
      real c1(4),c2(4),c3(4)
      data c1/1.,0.,0.,1./
      data c2/1.,-1.,-1.,1./
      data c3/1.,1.,0.,0./
      real*8 yvvv,xvvv
      double precision xi,yi,zi
      double precision xrt,yrt,zrt
      real AA,BB
      real yvv(4),xvv(4)

*     tollerance to consider the track inside the sensitive area
      real ptoll
      data ptoll/150./          !um

      external nviewx,nviewy,acoordsi,dcoord

      nplpt = nplPAM            !plane
      viewx = nviewx(nplpt)
      viewy = nviewy(nplpt)

      do il=1,nladders_view
         do is=1,2

            do iv=1,4           !loop on sensor vertexes
               stripx = (il-c1(iv))*1024 + c1(iv) + c2(iv)*3
               stripy = (il-c3(iv))*1024 + c3(iv)               
c------------------------------------------------------------------------
c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
c------------------------------------------------------------------------
               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
     $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
c                  print*,'whichsensor: ',
c     $                ' WARNING: false X strip: strip ',stripx
               endif
               xi = acoordsi(stripx,viewx)
               yi = acoordsi(stripy,viewy)
               zi = 0.
c------------------------------------------------------------------------
c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c------------------------------------------------------------------------
c     N.B. I convert angles from microradiants to radiants
               xrt = xi
     $              - omega(nplpt,il,is)*yi 
     $              + gamma(nplpt,il,is)*zi
     $              + dx(nplpt,il,is)               
               yrt = omega(nplpt,il,is)*xi
     $              + yi
     $              - beta(nplpt,il,is)*zi
     $              + dy(nplpt,il,is)               
               zrt = -gamma(nplpt,il,is)*xi
     $              + beta(nplpt,il,is)*yi
     $              + zi
     $              + dz(nplpt,il,is)
c------------------------------------------------------------------------
c     measured coordinates (in cm) in PAMELA reference system
c------------------------------------------------------------------------
               yvvv = dcoord(yrt,viewy,il,is) / 1.d4
               xvvv = dcoord(xrt,viewx,il,is) / 1.d4

               yvv(iv)=sngl(yvvv)
               xvv(iv)=sngl(xvvv)
c               print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '
c     $              ,iv,xvv(iv),yvv(iv)
            enddo               !end loop on sensor vertexes

            dtot=0.
            do iside=1,4,2        !loop on sensor edges X
               iv1=iside
               iv2=mod(iside,4)+1
*     straight line passing trhough two consecutive vertexes
               AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))
               BB = yvv(iv1) - AA*xvv(iv1)
*     point along the straight line closer to the track
               xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
               yoo = AA*xoo + BB
*     sum of the distances
               dtot = dtot + 
     $              sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
            enddo               !end loop on sensor edges
            do iside=2,4,2        !loop on sensor edges Y
               iv1=iside
               iv2=mod(iside,4)+1
*     straight line passing trhough two consecutive vertexes
               AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))
               BB = xvv(iv1) - AA*yvv(iv1)
*     point along the straight line closer to the track
               yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
               xoo = AA*yoo + BB
*     sum of the distances
               dtot = dtot + 
     $              sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
            enddo               !end loop on sensor edges


*     half-perimeter of sensitive area
            Perim = 
     $            SiDimX - edgeX_l - edgeX_r
     $           +SiDimY - edgeY_l - edgeY_r
            Perim = (Perim + ptoll)/1.e4
            if(dtot.le.Perim)goto 100


         enddo
      enddo
      
      ladder = 0 
      sensor = 0 
      goto 200

 100  continue
      ladder = il
      sensor = is
      

 200  return
      end


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

      subroutine reverse(v,n,temp) !invert the order of the components of v(n) vector

      implicit double precision (A-H,O-Z)

      dimension v(*)
      dimension temp(*)
      integer i,n

      do i=1,n
        temp(i)=v(n+1-i)
      enddo

      do i=1,n
        v(i)=temp(i)
      enddo

      return
      end

*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
      integer function ip_cp(id)
*
*     given the couple id, 
*     it returns the plane number
*     
      include 'commontracker.f'
      include 'level1.f'
c      include 'common_analysis.f'
      include 'common_momanhough.f'
      
      ip_cp=0
      ncpp=0
      do ip=1,nplanes
         ncpp=ncpp+ncp_plane(ip)
         if(ncpp.ge.abs(id))then
            ip_cp=ip
            goto 100
         endif
      enddo
 100  continue
      return
      end


      integer function is_cp(id)
*
*     given the couple id, 
*     it returns the sensor number 
*     
      is_cp=0
      if(id.lt.0)is_cp=1
      if(id.gt.0)is_cp=2
c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'

      return
      end


      integer function icp_cp(id)
*
*     given the couple id, 
*     it returns the id number ON THE PLANE
*     
      include 'commontracker.f'
      include 'level1.f'
c      include 'common_analysis.f'
      include 'common_momanhough.f'
      
      icp_cp=0

      ncpp=0
      do ip=1,nplanes
         ncppold=ncpp
         ncpp=ncpp+ncp_plane(ip)
         if(ncpp.ge.abs(id))then
            icp_cp=abs(id)-ncppold
            goto 100
         endif
      enddo
 100  continue
      return
      end
      

      integer function id_cp(ip,icp,is)
*
*     given a plane, a couple and the sensor 
*     it returns the absolute couple id 
*     negative if sensor =1
*     positive if sensor =2
*
      include 'commontracker.f'
      include 'level1.f'
c      include 'calib.f'
c      include 'level1.f'
c      include 'common_analysis.f'
      include 'common_momanhough.f'
      
      id_cp=0

      if(ip.gt.1)then
         do i=1,ip-1
            id_cp = id_cp + ncp_plane(i)
         enddo
      endif

      id_cp = id_cp + icp

      if(is.eq.1) id_cp = -id_cp

      return
      end


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

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

      subroutine cl_to_couples(iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
c      include 'momanhough_init.f'
      include 'calib.f'
c      include 'level1.f'

*     output flag
*     --------------
*     0 = good event
*     1 = bad event
*     --------------
      integer iflag

      integer badseed,badclx,badcly

*     init variables
      ncp_tot=0 
      do ip=1,nplanes 
         do ico=1,ncouplemax
            clx(ip,ico)=0
            cly(ip,ico)=0
         enddo
         ncp_plane(ip)=0
         do icl=1,nclstrmax_level2
            cls(ip,icl)=1
         enddo
         ncls(ip)=0
      enddo
      do icl=1,nclstrmax_level2
         cl_single(icl) = 1
         cl_good(icl)   = 0
      enddo
      do iv=1,nviews
         ncl_view(iv)  = 0
         mask_view(iv) = 0      !all included
      enddo
      
*     count number of cluster per view
      do icl=1,nclstr1
         ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1         
      enddo
*     mask views with too many clusters
      do iv=1,nviews
         if( ncl_view(iv).gt. nclusterlimit)then
            mask_view(iv) = 1 
            if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
     $           ,nclusterlimit,' on view ', iv,' --> masked!'
         endif
      enddo


*     start association
      ncouples=0
      do icx=1,nclstr1          !loop on cluster (X)
         if(mod(VIEW(icx),2).eq.1)goto 10
         
*     ----------------------------------------------------
*     jump masked views (X VIEW)
*     ----------------------------------------------------
         if( mask_view(VIEW(icx)).ne.0 ) goto 10
*     ----------------------------------------------------
*     cut on charge (X VIEW)
*     ----------------------------------------------------
         if(sgnl(icx).lt.dedx_x_min)then
            cl_single(icx)=0
            goto 10
         endif
*     ----------------------------------------------------
*     cut BAD (X VIEW)            
*     ----------------------------------------------------
         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
         ifirst=INDSTART(icx)
         if(icx.ne.nclstr1) then
            ilast=INDSTART(icx+1)-1
         else
            ilast=TOTCLLENGTH
         endif
         badclx=badseed
         do igood=-ngoodstr,ngoodstr
            ibad=1
            if((INDMAX(icx)+igood).gt.ifirst.and.
     $           (INDMAX(icx)+igood).lt.ilast.and.
     $           .true.)then
               ibad=BAD(VIEW(icx),
     $              nvk(MAXS(icx)+igood),
     $              nst(MAXS(icx)+igood))
            endif
            badclx=badclx*ibad
         enddo
*     ----------------------------------------------------
*     >>> eliminato il taglio sulle BAD <<<
*     ----------------------------------------------------
c     if(badcl.eq.0)then 
c     cl_single(icx)=0
c     goto 10
c     endif
*     ----------------------------------------------------
         
         cl_good(icx)=1
         nplx=npl(VIEW(icx)) 
         nldx=nld(MAXS(icx),VIEW(icx))
         
         do icy=1,nclstr1       !loop on cluster (Y)
            if(mod(VIEW(icy),2).eq.0)goto 20
            
*     ----------------------------------------------------
*     jump masked views (Y VIEW)
*     ----------------------------------------------------
            if( mask_view(VIEW(icy)).ne.0 ) goto 20

*     ----------------------------------------------------
*     cut on charge (Y VIEW) 
*     ----------------------------------------------------
            if(sgnl(icy).lt.dedx_y_min)then 
               cl_single(icy)=0
               goto 20
            endif
*     ----------------------------------------------------
*     cut BAD (Y VIEW)            
*     ----------------------------------------------------
            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
            ifirst=INDSTART(icy)
            if(icy.ne.nclstr1) then
               ilast=INDSTART(icy+1)-1 
            else
               ilast=TOTCLLENGTH 
            endif
            badcly=badseed
            do igood=-ngoodstr,ngoodstr
               ibad=1
               if((INDMAX(icy)+igood).gt.ifirst.and.
     $              (INDMAX(icy)+igood).lt.ilast.and.
     $              .true.)
     $              ibad=BAD(VIEW(icy),
     $              nvk(MAXS(icy)+igood),
     $              nst(MAXS(icy)+igood))
               badcly=badcly*ibad
            enddo
*     ----------------------------------------------------
*     >>> eliminato il taglio sulle BAD <<<
*     ----------------------------------------------------
c     if(badcl.eq.0)then 
c     cl_single(icy)=0
c     goto 20
c     endif
*     ----------------------------------------------------
            
            cl_good(icy)=1                  
            nply=npl(VIEW(icy)) 
            nldy=nld(MAXS(icy),VIEW(icy))
            
*     ----------------------------------------------
*     CONDITION TO FORM A COUPLE 
*     ----------------------------------------------
*     geometrical consistency (same plane and ladder)
            if(nply.eq.nplx.and.nldy.eq.nldx)then
*     charge correlation
*     (modified to be applied only below saturation... obviously)

               if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
     $              .and.
     $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
     $              .and. 
     $              (badclx.eq.1.and.badcly.eq.1)
     $              .and.
     $              .true.)then

                  ddd=(sgnl(icy)
     $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
                  ddd=ddd/sqrt(kch(nplx,nldx)**2+1)

c                  cut = chcut * sch(nplx,nldx)

                  sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
     $                 -kch(nplx,nldx)*cch(nplx,nldx))
                  sss=sss/sqrt(kch(nplx,nldx)**2+1)
                  cut = chcut * (16 + sss/50.)

                  if(abs(ddd).gt.cut)then 
                     goto 20    !charge not consistent
                  endif
               endif

               if(ncp_plane(nplx).gt.ncouplemax)then
                  if(verbose)print*,
     $                 '** warning ** number of identified '// 
     $                 'couples on plane ',nplx,
     $                 'exceeds vector dimention '
     $                 ,'( ',ncouplemax,' ) --> masked!'
                  mask_view(nviewx(nplx)) = 2
                  mask_view(nviewy(nply)) = 2
                  goto 10 
               endif
               
*     ------------------> COUPLE <------------------
               ncp_plane(nplx) = ncp_plane(nplx) + 1
               clx(nplx,ncp_plane(nplx))=icx
               cly(nply,ncp_plane(nplx))=icy
               cl_single(icx)=0
               cl_single(icy)=0
*     ----------------------------------------------

            endif                              

 20         continue
         enddo                  !end loop on clusters(Y)
         
 10      continue
      enddo                     !end loop on clusters(X)
      
      
      do icl=1,nclstr1
         if(cl_single(icl).eq.1)then
            ip=npl(VIEW(icl)) 
            ncls(ip)=ncls(ip)+1
            cls(ip,ncls(ip))=icl
         endif
      enddo
      
      
      if(DEBUG)then
         print*,'clusters  ',nclstr1
         print*,'good    ',(cl_good(i),i=1,nclstr1)
         print*,'singles ',(cl_single(i),i=1,nclstr1)
         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
      endif
      
      do ip=1,6
         ncp_tot = ncp_tot + ncp_plane(ip)
      enddo
      
      return
      end
      
***************************************************
*                                                 *
*                                                 *
*                                                 *
*                                                 *
*                                                 *
*                                                 *
**************************************************

      subroutine cp_to_doubtrip(iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'common_xyzPAM.f'
      include 'common_mini_2.f'
      include 'calib.f'


*     output flag
*     --------------
*     0 = good event
*     1 = bad event
*     --------------
      integer iflag


*     -----------------------------
*     DOUBLETS/TRIPLETS coordinates
c      double precision xm1,ym1,zm1
c      double precision xm2,ym2,zm2
c      double precision xm3,ym3,zm3

      real xm1,ym1,zm1
      real xm2,ym2,zm2
      real xm3,ym3,zm3
*     -----------------------------
*     variable needed for tricircle:
      real xp(3),zp(3)!TRIPLETS coordinates, to find a circle
      EQUIVALENCE (xm1,xp(1))
      EQUIVALENCE (xm2,xp(2))
      EQUIVALENCE (xm3,xp(3))
      EQUIVALENCE (zm1,zp(1))
      EQUIVALENCE (zm2,zp(2))
      EQUIVALENCE (zm3,zp(3))
      real angp(3),resp(3),chi
      real xc,zc,radius
*     -----------------------------


*     --------------------------------------------
*     put a limit to the maximum number of couples
*     per plane, in order to apply hough transform
*     (couples recovered during track refinement)
*     --------------------------------------------
      do ip=1,nplanes
         if(ncp_plane(ip).gt.ncouplelimit)then
            mask_view(nviewx(ip)) = 8
            mask_view(nviewy(ip)) = 8
         endif
      enddo


      ndblt=0                   !number of doublets
      ntrpt=0                   !number of triplets
      
      do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
         if(  mask_view(nviewx(ip1)).ne.0 .or. 
     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
         do is1=1,2             !loop on sensors - COPPIA 1            
            do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
               icx1=clx(ip1,icp1)
               icy1=cly(ip1,icp1)
c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.) 
               xm1=xPAM
               ym1=yPAM
               zm1=zPAM                  
c     print*,'***',is1,xm1,ym1,zm1

               do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
                  if(  mask_view(nviewx(ip2)).ne.0 .or. 
     $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
                  do is2=1,2    !loop on sensors -ndblt COPPIA 2
                     
                     do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
                        icx2=clx(ip2,icp2)
                        icy2=cly(ip2,icp2)
c                        call xyz_PAM
c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
c                        call xyz_PAM
c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
                        call xyz_PAM
     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.) 
                        xm2=xPAM
                        ym2=yPAM
                        zm2=zPAM
                                                
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*     track parameters on Y VIEW
*     (2 couples needed)
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                        if(ndblt.eq.ndblt_max)then
                           if(verbose)print*,
     $                          '** warning ** number of identified '// 
     $                          'doublets exceeds vector dimention '
     $                          ,'( ',ndblt_max,' )'
c                           good2=.false.
c                           goto 880 !fill ntp and go to next event
                           do iv=1,12 
                              mask_view(iv) = 3
                           enddo
                           iflag=1
                           return
                        endif
                        ndblt = ndblt + 1
*     store doublet info
                        cpyz1(ndblt)=id_cp(ip1,icp1,is1)
                        cpyz2(ndblt)=id_cp(ip2,icp2,is2)
*     tg(th_yz)
                        alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
*     y0 (cm)
                        alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
                           
****  -----------------------------------------------****
****  reject non phisical couples                    ****
****  -----------------------------------------------****
                        if(
     $                       abs(alfayz2(ndblt)).gt.alfyz2_max
     $                       .or.
     $                       abs(alfayz1(ndblt)).gt.alfyz1_max
     $                       )ndblt = ndblt-1
                        
c$$$      if(iev.eq.33)then
c$$$      print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2
c$$$     $        ,' || ',icx1,icy1,icx2,icy2
c$$$     $        ,' || ',xm1,ym1,xm2,ym2
c$$$     $        ,' || ',alfayz2(ndblt),alfayz1(ndblt)
c$$$      endif
c$$$
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*     track parameters on Y VIEW - end
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 

                        if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples

                        do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
                           if(  mask_view(nviewx(ip3)).ne.0 .or. 
     $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
                           do is3=1,2 !loop on sensors - COPPIA 3
                              
                              do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
                                 icx3=clx(ip3,icp3)
                                 icy3=cly(ip3,icp3)
c                                 call xyz_PAM
c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
c                                 call xyz_PAM
c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1) 
                                 call xyz_PAM
     $                                (icx3,icy3,is3,PFAdef,PFAdef
     $                                ,0.,0.,0.,0.) 
                                 xm3=xPAM
                                 ym3=yPAM
                                 zm3=zPAM
*     find the circle passing through the three points
                                 call tricircle(3,xp,zp,angp,resp,chi
     $                                ,xc,zc,radius,iflag)
c     print*,xc,zc,radius
*     the circle must intersect the reference plane
                                 if(
c     $                                 (xc.le.-1.*xclimit.or.
c     $                                 xc.ge.xclimit).and.
     $                                radius**2.ge.(ZINI-zc)**2.and.
     $                                iflag.eq.0.and.
     $                                .true.)then
                                 
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*     track parameters on X VIEW
*     (3 couples needed)
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                                 if(ntrpt.eq.ntrpt_max)then
                                    if(verbose)print*,
     $                     '** warning ** number of identified '// 
     $                     'triplets exceeds vector dimention '
     $                    ,'( ',ntrpt_max,' )'
c                                    good2=.false.
c                                    goto 880 !fill ntp and go to next event
                                    do iv=1,nviews
                                       mask_view(iv) = 4
                                    enddo
                                    iflag=1
                                    return
                                 endif
                                 ntrpt = ntrpt +1
*     store triplet info
                                 cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
                                 cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
                                 cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
                                 
                                 if(xc.lt.0)then
*************POSITIVE DEFLECTION
              alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
              alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
              alfaxz3(ntrpt) = 1/radius
                                 else
*************NEGATIVE DEFLECTION
              alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
              alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
              alfaxz3(ntrpt) = -1/radius
                                 endif
                                 
****  -----------------------------------------------****
****  reject non phisical triplets                   ****
****  -----------------------------------------------****
                                 if(
     $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max
     $                                .or.
     $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max
     $                                )ntrpt = ntrpt-1
                                 
                                 
c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*     track parameters on X VIEW - end
*     - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                              endif
                           enddo !end loop on COPPIA 3
                        enddo   !end loop on sensors - COPPIA 3
 30                     continue
                     enddo      !end loop on planes  - COPPIA 3
 31                  continue
                     
 1                enddo         !end loop on COPPIA 2
               enddo            !end loop on sensors - COPPIA 2
 20            continue
            enddo               !end loop on planes  - COPPIA 2
            
         enddo                  !end loop on COPPIA1
      enddo                     !end loop on sensors - COPPIA 1
 10   continue
      enddo                     !end loop on planes  - COPPIA 1
      
      if(DEBUG)then
         print*,'--- doublets ',ndblt
         print*,'--- triplets ',ntrpt
      endif
      
c     goto 880               !ntp fill
      
      
      return
      end


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

      subroutine doub_to_YZcloud(iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
c      include 'momanhough_init.f'


*     output flag
*     --------------
*     0 = good event
*     1 = bad event
*     --------------
      integer iflag

      integer db_used(ndblt_max)
      integer db_temp(ndblt_max)
      integer db_all(ndblt_max) !stores db ID in each cloud

      integer hit_plane(nplanes)

*     mask for used couples 
      integer cp_useds1(ncouplemaxtot) ! sensor 1
      integer cp_useds2(ncouplemaxtot) ! sensor 2


*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*     classification of DOUBLETS
*     according to distance in parameter space
*     (cloud = group of points (doublets) in parameter space)
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      do idb=1,ndblt
         db_used(idb)=0
      enddo

      distance=0
      nclouds_yz=0              !number of clouds
      npt_tot=0
      nloop=0                   
 90   continue                  
      do idb1=1,ndblt           !loop (1) on DOUBLETS
         if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
            
c     print*,'--------------'
c     print*,'** ',idb1,' **'
            
         do icp=1,ncp_tot
            cp_useds1(icp)=0    !init
            cp_useds2(icp)=0    !init
         enddo
         do idb=1,ndblt
            db_all(idb)=0
         enddo
         if(cpyz1(idb1).gt.0)cp_useds2(cpyz1(idb1))=1
         if(cpyz1(idb1).lt.0)cp_useds1(-cpyz1(idb1))=1
         if(cpyz2(idb1).gt.0)cp_useds2(cpyz2(idb1))=1
         if(cpyz2(idb1).lt.0)cp_useds1(-cpyz2(idb1))=1
         temp1 = alfayz1(idb1)
         temp2 = alfayz2(idb1)
         npt=1                  !counter of points in the cloud

         db_all(npt) = idb1
         
         nptloop=1
         db_temp(1)=idb1
         
 88      continue
         
         npv=0                  !# new points inlcuded
         do iloop=1,nptloop
            idbref=db_temp(iloop) !local point of reference
ccccc if(db_used(idbref).eq.1)goto 1188 !next
            
            do idb2=1,ndblt     !loop (2) on DOUBLETS
               if(idb2.eq.idbref)goto 1118 !next doublet
               if(db_used(idb2).eq.1)goto 1118
               
                  
*     doublet distance in parameter space
               distance=
     $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2               
               distance = sqrt(distance)
               
c$$$      if(iev.eq.33)then
c$$$      if(distance.lt.100)
c$$$     $ print*,'********* ',idb1,idbref,idb2,distance
c$$$      if(distance.lt.100)
c$$$     $ print*,'********* ',alfayz1(idbref),alfayz1(idb2)
c$$$     $                    ,alfayz2(idbref),alfayz2(idb2)
c$$$      endif
               if(distance.lt.cutdistyz)then

c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
                  if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
                  if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
                  if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
                  if(cpyz2(idb2).lt.0)cp_useds1(-cpyz2(idb2))=1
                  npt = npt + 1 !counter of points in the cloud

                  npv = npv +1
                  db_temp(npv) = idb2
                  db_used(idbref) = 1 
                  db_used(idb2) = 1  
                  
                  db_all(npt) = idb2

                  temp1 = temp1 + alfayz1(idb2)
                  temp2 = temp2 + alfayz2(idb2)
c     print*,'*   idbref,idb2 ',idbref,idb2
               endif               
               
 1118          continue 
            enddo               !end loop (2) on DOUBLETS
            
 1188       continue
         enddo                  !end loop on... bo?
         
         nptloop=npv
         if(nptloop.ne.0)goto 88
         
*     ------------------------------------------
*     stores the cloud only if 
*     1) it includes a minimum number of REAL couples
*     1bis) it inlcudes a minimum number of doublets
*     2) it is not already stored
*     ------------------------------------------
         do ip=1,nplanes
            hit_plane(ip)=0
         enddo
         ncpused=0
         do icp=1,ncp_tot
            if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
               ncpused=ncpused+1
               ip=ip_cp(icp)
               hit_plane(ip)=1
            endif
         enddo
         nplused=0
         do ip=1,nplanes
            nplused=nplused+ hit_plane(ip)
         enddo
c     print*,'>>>> ',ncpused,npt,nplused
c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
         if(npt.lt.nptyz_min)goto 2228 !next doublet
         if(nplused.lt.nplyz_min)goto 2228 !next doublet
         
*     ~~~~~~~~~~~~~~~~~
*     >>> NEW CLOUD <<<

         if(nclouds_yz.ge.ncloyz_max)then
            if(verbose)print*,
     $           '** warning ** number of identified '// 
     $           'YZ clouds exceeds vector dimention '
     $           ,'( ',ncloyz_max,' )'
c               good2=.false.
c     goto 880         !fill ntp and go to next event
            do iv=1,nviews
               mask_view(iv) = 5
            enddo
            iflag=1
            return
         endif
         
         nclouds_yz = nclouds_yz + 1 !increase counter
         alfayz1_av(nclouds_yz) = temp1/npt !store average parameter 
         alfayz2_av(nclouds_yz) = temp2/npt ! "
         do icp=1,ncp_tot
            cpcloud_yz(nclouds_yz,icp)=
     $           cp_useds1(icp)+2*cp_useds2(icp) !store cp info
         enddo 
         ptcloud_yz(nclouds_yz)=npt 
c     ptcloud_yz_nt(nclouds_yz)=npt 
         do ipt=1,npt
            db_cloud(npt_tot+ipt) = db_all(ipt)
c     print*,'>> ',ipt,db_all(ipt)
         enddo  
         npt_tot=npt_tot+npt
         if(DEBUG)then
            print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
            print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
            print*,'- alfayz1 ',alfayz1_av(nclouds_yz)
            print*,'- alfayz2 ',alfayz2_av(nclouds_yz)
            print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
            print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
            print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
c$$$     $           ,ptcloud_yz(nclouds_yz)
c$$$            print*,'nt-uple: db_cloud(...) = '
c$$$     $           ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)
         endif
*     >>> NEW CLOUD <<<
*     ~~~~~~~~~~~~~~~~~
 2228    continue
      enddo                     !end loop (1) on DOUBLETS
      
      
      if(nloop.lt.nstepy)then       
        cutdistyz = cutdistyz+cutystep
        nloop     = nloop+1           
        goto 90                 
      endif                     
      
      if(DEBUG)then
         print*,'---------------------- '
         print*,'Y-Z total clouds ',nclouds_yz
         print*,' '
      endif
      
      
      return
      end
      

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

      subroutine trip_to_XZcloud(iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
c      include 'momanhough_init.f'


*     output flag
*     --------------
*     0 = good event
*     1 = bad event
*     --------------
      integer iflag

      integer tr_used(ntrpt_max)
      integer tr_temp(ntrpt_max)
      integer tr_incl(ntrpt_max)
      integer tr_all(ntrpt_max) !stores tr ID in each cloud

      integer hit_plane(nplanes)

*     mask for used couples 
      integer cp_useds1(ncouplemaxtot) ! sensor 1
      integer cp_useds2(ncouplemaxtot) ! sensor 2

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*     classification of TRIPLETS
*     according to distance in parameter space
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      do itr=1,ntrpt
         tr_used(itr)=0
      enddo

      distance=0
      nclouds_xz=0              !number of clouds         
      npt_tot=0                 !total number of selected triplets
      nloop=0                   
 91   continue                  
      do itr1=1,ntrpt           !loop (1) on TRIPLETS
         if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
c     print*,'--------------'
c     print*,'** ',itr1,' **'
         
         do icp=1,ncp_tot
            cp_useds1(icp)=0
            cp_useds2(icp)=0
         enddo
         do itr=1,ntrpt
            tr_all(itr)=0       !list of included triplets
         enddo
         if(cpxz1(itr1).gt.0)cp_useds2(cpxz1(itr1))=1
         if(cpxz1(itr1).lt.0)cp_useds1(-cpxz1(itr1))=1
         if(cpxz2(itr1).gt.0)cp_useds2(cpxz2(itr1))=1
         if(cpxz2(itr1).lt.0)cp_useds1(-cpxz2(itr1))=1
         if(cpxz3(itr1).gt.0)cp_useds2(cpxz3(itr1))=1
         if(cpxz3(itr1).lt.0)cp_useds1(-cpxz3(itr1))=1
         temp1 = alfaxz1(itr1)
         temp2 = alfaxz2(itr1)
         temp3 = alfaxz3(itr1) 
         npt=1                  !counter of points in the cloud
         
         tr_all(npt) = itr1
         
         nptloop=1
c         tr_temp(1)=itr1
         tr_incl(1)=itr1
         
 8881    continue
         
         npv=0                  !# new points inlcuded
         do iloop=1,nptloop
            itrref=tr_incl(iloop) !local point of reference
            do itr2=1,ntrpt     !loop (2) on TRIPLETS
               if(itr2.eq.itr1)goto 11188       !next triplet
               if(tr_used(itr2).eq.1)goto 11188 !next triplet               
*     triplet distance in parameter space
*     solo i due parametri spaziali per il momemnto
               distance=
     $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2               
               distance = sqrt(distance)
               
               if(distance.lt.cutdistxz)then
c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
                  if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
                  if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
                  if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
                  if(cpxz2(itr2).lt.0)cp_useds1(-cpxz2(itr2))=1
                  if(cpxz3(itr2).gt.0)cp_useds2(cpxz3(itr2))=1
                  if(cpxz3(itr2).lt.0)cp_useds1(-cpxz3(itr2))=1
                  npt = npt + 1 !counter of points in the cloud
                  
                  npv = npv +1
                  tr_temp(npv) = itr2
                  tr_used(itrref) = 1
                  tr_used(itr2) = 1
                  
                  tr_all(npt) = itr2
                  
                  temp1 = temp1 + alfaxz1(itr2)
                  temp2 = temp2 + alfaxz2(itr2) 
                  temp3 = temp3 + alfaxz3(itr2)
c     print*,'*   itrref,itr2 ',itrref,itr2,distance 
               endif               
               
11188          continue 
            enddo               !end loop (2) on TRIPLETS
                       
11888       continue
         enddo                  !end loop on... bo?     
         
         nptloop=npv
         do i=1,npv
            tr_incl(i)=tr_temp(i)
         enddo
         if(nptloop.ne.0)goto 8881 
         
*     ------------------------------------------
*     stores the cloud only if 
*     1) it includes a minimum number of REAL couples
*     1bis)
*     2) it is not already stored
*     ------------------------------------------
c     print*,'check cp_used'
         do ip=1,nplanes
            hit_plane(ip)=0
         enddo
         ncpused=0
         do icp=1,ncp_tot
            if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
               ncpused=ncpused+1
               ip=ip_cp(icp)
               hit_plane(ip)=1
            endif
         enddo
         nplused=0
         do ip=1,nplanes
            nplused=nplused+ hit_plane(ip)
         enddo
c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
         if(npt.lt.nptxz_min)goto 22288     !next triplet
         if(nplused.lt.nplxz_min)goto 22288 !next triplet
         
*     ~~~~~~~~~~~~~~~~~
*     >>> NEW CLOUD <<<
         if(nclouds_xz.ge.ncloxz_max)then
            if(verbose)print*,
     $           '** warning ** number of identified '// 
     $           'XZ clouds exceeds vector dimention '
     $           ,'( ',ncloxz_max,' )'
c     good2=.false.
c     goto 880         !fill ntp and go to next event
            do iv=1,nviews
               mask_view(iv) = 6
            enddo
            iflag=1
            return
         endif
         nclouds_xz = nclouds_xz + 1 !increase counter
         alfaxz1_av(nclouds_xz) = temp1/npt !store average parameter 
         alfaxz2_av(nclouds_xz) = temp2/npt ! "
         alfaxz3_av(nclouds_xz) = temp3/npt ! " 
         do icp=1,ncp_tot
            cpcloud_xz(nclouds_xz,icp)=
     $           cp_useds1(icp)+2*cp_useds2(icp) !store cp info
         enddo
         ptcloud_xz(nclouds_xz)=npt
         do ipt=1,npt
            tr_cloud(npt_tot+ipt) = tr_all(ipt)
         enddo
         npt_tot=npt_tot+npt
         
         if(DEBUG)then
            print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'               
            print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
            print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
            print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
            print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
            print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
            print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
c$$$     $           ,ptcloud_xz(nclouds_xz)
c$$$            print*,'nt-uple: tr_cloud(...) = '
c$$$     $           ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot)
         endif
*     >>> NEW CLOUD <<<
*     ~~~~~~~~~~~~~~~~~
22288    continue
      enddo                     !end loop (1) on DOUBLETS

       if(nloop.lt.nstepx)then       
         cutdistxz=cutdistxz+cutxstep 
         nloop=nloop+1          
         goto 91                
       endif                    
       
      if(DEBUG)then
         print*,'---------------------- '
         print*,'X-Z total clouds ',nclouds_xz
         print*,' '
      endif
      
      
      return
      end
      

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

      subroutine clouds_to_ctrack(iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'common_xyzPAM.f'
      include 'common_mini_2.f'
      include 'common_mech.f'


*     output flag
*     --------------
*     0 = good event
*     1 = bad event
*     --------------
      integer iflag

*     -----------------------------------------------------------
*     mask to store (locally) the couples included
*     in the intersection bewteen a XZ and YZ cloud
      integer cpintersec(ncouplemaxtot)
*     -----------------------------------------------------------
*     list of matching couples in the combination 
*     between a XZ and YZ cloud
      integer cp_match(nplanes,2*ncouplemax)
      integer ncp_match(nplanes)
*     -----------------------------------------------------------
      integer hit_plane(nplanes)
*     -----------------------------------------------------------
*     variables for track fitting
      double precision AL_INI(5)
c      double precision tath
*     -----------------------------------------------------------
c      real fitz(nplanes)        !z coordinates of the planes in cm


      ntracks=0                 !counter of track candidates
      
      do iyz=1,nclouds_yz       !loop on YZ couds
         do ixz=1,nclouds_xz    !loop on XZ couds
            
*     --------------------------------------------------
*     check of consistency of the clouds
*     ---> required a minimum number of matching couples
*     the track fit will be performed on the INTERSECTION
*     of the two clouds
*     --------------------------------------------------
            do ip=1,nplanes
               hit_plane(ip)=0
               ncp_match(ip)=0
               do icpp=1,ncouplemax
                  cp_match(ip,icpp)=0 !init couple list
               enddo
            enddo
            ncp_ok=0
            do icp=1,ncp_tot    !loop on couples
*     get info on 
               cpintersec(icp)=min(
     $              cpcloud_yz(iyz,icp),
     $              cpcloud_xz(ixz,icp))
               if(
     $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
     $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
     $              .false.)cpintersec(icp)=0
               if(cpintersec(icp).ne.0)then
                  ncp_ok=ncp_ok+1   
                  
                  ip=ip_cp(icp)
                  hit_plane(ip)=1
                  if(cpintersec(icp).eq.1)then 
*     1) only the couple image in sensor 1 matches
                     id=-icp
                     ncp_match(ip)=ncp_match(ip)+1
                     cp_match(ip,ncp_match(ip))=id
                  elseif(cpintersec(icp).eq.2)then
*     2) only the couple image in sensor 2 matches
                     id=icp
                     ncp_match(ip)=ncp_match(ip)+1
                     cp_match(ip,ncp_match(ip))=id
                  else
*     3) both couple images match
                     id=icp
                     do is=1,2
                        id=-id
                        ncp_match(ip)=ncp_match(ip)+1
                        cp_match(ip,ncp_match(ip))=id
                     enddo
                  endif                     
               endif            !end matching condition
            enddo               !end loop on couples
            
            nplused=0
            do ip=1,nplanes
               nplused=nplused+ hit_plane(ip)
            enddo
            
c            if(nplused.lt.nplxz_min)goto 888 !next doublet
            if(nplused.lt.nplyz_min)goto 888 !next doublet
            if(ncp_ok.lt.ncpok)goto 888 !next cloud
            
            if(DEBUG)then
               print*,'Combination ',iyz,ixz
     $              ,' db ',ptcloud_yz(iyz)
     $              ,' tr ',ptcloud_xz(ixz)
     $              ,'  -----> # matching couples ',ncp_ok
            endif
c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
c$$$  print*,'Configurazione cluster XZ'
c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))
c$$$  print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))
c$$$  print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))
c$$$  print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))
c$$$  print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))
c$$$  print*,'Configurazione cluster YZ'
c$$$  print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))
c$$$  print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))
c$$$  print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))
c$$$  print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))
c$$$  print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))
c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))
c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
            
*     -------> INITIAL GUESS <-------
cccc       SBAGLIATO
c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
c$$$            AL_INI(4) = PIGR + 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.
cccc       GIUSTO (ma si sua guess())
c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
c$$$     $           /dreal(alfaxz2_av(ixz)))
c$$$            ELSE
c$$$               AL_INI(4) = acos(-1.)/2 
c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
c$$$            ENDIF
c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
c$$$            
c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
c$$$            
c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
                        
            if(DEBUG)then
               print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
               print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
               print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
               print*,'4 >>> ',(cp_match(3,i),i=1,ncp_match(3))
               print*,'5 >>> ',(cp_match(2,i),i=1,ncp_match(2))
               print*,'6 >>> ',(cp_match(1,i),i=1,ncp_match(1))
            endif
            
            do icp1=1,max(1,ncp_match(1))
               hit_plane(1)=icp1
               if(ncp_match(1).eq.0)hit_plane(1)=0 !-icp1
               
               do icp2=1,max(1,ncp_match(2))
                  hit_plane(2)=icp2
                  if(ncp_match(2).eq.0)hit_plane(2)=0 !-icp2
                  
                  do icp3=1,max(1,ncp_match(3))
                     hit_plane(3)=icp3
                     if(ncp_match(3).eq.0)hit_plane(3)=0 !-icp3
                     
                     do icp4=1,max(1,ncp_match(4))
                        hit_plane(4)=icp4
                        if(ncp_match(4).eq.0)hit_plane(4)=0 !-icp4
                        
                        do icp5=1,max(1,ncp_match(5))
                           hit_plane(5)=icp5
                           if(ncp_match(5).eq.0)hit_plane(5)=0 !-icp5
                           
                           do icp6=1,max(1,ncp_match(6))
                              hit_plane(6)=icp6
                              if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
                              
                              
                              call track_init !init TRACK common

                              do ip=1,nplanes !loop on planes
                                 if(hit_plane(ip).ne.0)then 
                                    id=cp_match(ip,hit_plane(ip))
                                    is=is_cp(id)
                                    icp=icp_cp(id)
                                    if(ip_cp(id).ne.ip)
     $                                   print*,'OKKIO!!'
     $                                   ,'id ',id,is,icp
     $                                   ,ip_cp(id),ip
                                    icx=clx(ip,icp)
                                    icy=cly(ip,icp)
*                                   *************************
c                                    call xyz_PAM(icx,icy,is,
c     $                                   'COG2','COG2',0.,0.)
c                                    call xyz_PAM(icx,icy,is, !(1)
c     $                                   PFAdef,PFAdef,0.,0.) !(1)
                                    call xyz_PAM(icx,icy,is, !(1)
     $                                   PFAdef,PFAdef,0.,0.,0.,0.) 
*                                   *************************
*                                   -----------------------------
                                    xgood(nplanes-ip+1)=1.
                                    ygood(nplanes-ip+1)=1.
                                    xm(nplanes-ip+1)=xPAM
                                    ym(nplanes-ip+1)=yPAM
                                    zm(nplanes-ip+1)=zPAM
                                    resx(nplanes-ip+1)=resxPAM
                                    resy(nplanes-ip+1)=resyPAM
*                                   -----------------------------
                                 endif
                              enddo !end loop on planes
*     **********************************************************
*     ************************** FIT *** FIT *** FIT *** FIT ***
*     **********************************************************
cccc  scommentare se si usa al_ini della nuvola
c$$$                              do i=1,5
c$$$                                 AL(i)=AL_INI(i)
c$$$                              enddo
                              call guess()
                              do i=1,5
                                 AL_INI(i)=AL(i)
                              enddo
                              ifail=0 !error flag in chi^2 computation
                              jstep=0 !number of  minimization steps
                              iprint=0
c                              if(DEBUG)iprint=1
                              if(DEBUG)iprint=2
                              call mini2(jstep,ifail,iprint)
                              if(ifail.ne.0) then
                                 if(DEBUG)then
                                    print *,
     $                              '*** MINIMIZATION FAILURE *** ' 
     $                              //'(clouds_to_ctrack)'
                                    print*,'initial guess: '

                                    print*,'AL_INI(1) = ',AL_INI(1)
                                    print*,'AL_INI(2) = ',AL_INI(2)
                                    print*,'AL_INI(3) = ',AL_INI(3)
                                    print*,'AL_INI(4) = ',AL_INI(4)
                                    print*,'AL_INI(5) = ',AL_INI(5)
                                 endif
c                                 chi2=-chi2 
                              endif 
*     **********************************************************
*     ************************** FIT *** FIT *** FIT *** FIT ***
*     **********************************************************

                              if(chi2.le.0.)goto 666              

*     --------------------------
*     STORE candidate TRACK INFO
*     --------------------------
                              if(ntracks.eq.NTRACKSMAX)then
                                 
                                 if(verbose)print*,
     $                 '** warning ** number of candidate tracks '// 
     $                 ' exceeds vector dimension '
     $                ,'( ',NTRACKSMAX,' )'
c                                 good2=.false.
c                                 goto 880 !fill ntp and go to next event                     
                                 do iv=1,nviews
                                    mask_view(iv) = 7
                                 enddo
                                 iflag=1
                                 return
                              endif
                              
                              ntracks = ntracks + 1
                              
                              do ip=1,nplanes

                                 XV_STORE(ip,ntracks)=sngl(xv(ip))
                                 YV_STORE(ip,ntracks)=sngl(yv(ip))
                                 ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
                                 XM_STORE(ip,ntracks)=sngl(xm(ip))
                                 YM_STORE(ip,ntracks)=sngl(ym(ip))
                                 ZM_STORE(ip,ntracks)=sngl(zm(ip))
                                 RESX_STORE(ip,ntracks)=sngl(resx(ip))
                                 RESY_STORE(ip,ntracks)=sngl(resy(ip))
                                 XV_STORE(ip,ntracks)=sngl(xv(ip))
                                 YV_STORE(ip,ntracks)=sngl(yv(ip))
                                 ZV_STORE(ip,ntracks)=sngl(zv(ip))
                                 AXV_STORE(ip,ntracks)=sngl(axv(ip))
                                 AYV_STORE(ip,ntracks)=sngl(ayv(ip))
                                 XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
                                 YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
                                 if(hit_plane(ip).ne.0)then
                                    CP_STORE(nplanes-ip+1,ntracks)=
     $                                   cp_match(ip,hit_plane(ip))
                                    SENSOR_STORE(nplanes-ip+1,ntracks)
     $                              = is_cp(cp_match(ip,hit_plane(ip)))
                                    LADDER_STORE(nplanes-ip+1,ntracks)
     $                                   = LADDER(
     $                                   clx(ip,icp_cp(
     $                                   cp_match(ip,hit_plane(ip)
     $                                   ))));
                                 else
                                    CP_STORE(nplanes-ip+1,ntracks)=0
                                    SENSOR_STORE(nplanes-ip+1,ntracks)=0
                                    LADDER_STORE(nplanes-ip+1,ntracks)=0
                                 endif
                                 BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
                                 BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
                                 CLS_STORE(nplanes-ip+1,ntracks)=0
                                 do i=1,5
                                    AL_STORE(i,ntracks)=sngl(AL(i))
                                 enddo
                              enddo
                              
                              RCHI2_STORE(ntracks)=chi2
                              
*     --------------------------------
*     STORE candidate TRACK INFO - end
*     --------------------------------
                              
 666                          continue
                           enddo !end loop on cp in plane 6
                        enddo   !end loop on cp in plane 5
                     enddo      !end loop on cp in plane 4
                  enddo         !end loop on cp in plane 3
               enddo            !end loop on cp in plane 2
            enddo               !end loop on cp in plane 1
            
 888        continue
         enddo                  !end loop on XZ couds
      enddo                     !end loop on YZ couds
      
      if(ntracks.eq.0)then
         iflag=1
         return
      endif
      
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
      if(DEBUG)then
        print*,'****** TRACK CANDIDATES *****************'
        print*,'#         R. chi2        RIG         ndof'
        do i=1,ntracks
          ndof=0                !(1)
          do ii=1,nplanes       !(1)
            ndof=ndof           !(1)
     $           +int(xgood_store(ii,i)) !(1)
     $           +int(ygood_store(ii,i)) !(1)
          enddo                 !(1)
          print*,i,' --- ',rchi2_store(i),' --- '
     $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
        enddo
        print*,'*****************************************'
      endif
      
      
      return
      end
      

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

      subroutine refine_track(ibest)


      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'common_xyzPAM.f'
      include 'common_mini_2.f'
      include 'common_mech.f'
      include 'calib.f'

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

      real k(6)
      DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/

      real xp,yp,zp
      real xyzp(3),bxyz(3)
      equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))

*     =================================================
*     new estimate of positions using ETA algorithm
*                          and
*     search for new couples and single clusters to add
*     =================================================
      call track_init
      do ip=1,nplanes           !loop on planes

         xP=XV_STORE(nplanes-ip+1,ibest)
         yP=YV_STORE(nplanes-ip+1,ibest)
         zP=ZV_STORE(nplanes-ip+1,ibest)
         call gufld(xyzp,bxyz)
         BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
         BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
c$$$  bxyz(1)=0
c$$$         bxyz(2)=0
c$$$         bxyz(3)=0
*     |||||||||||||||||||||||||||||||||||||||||||||||||
*     -------------------------------------------------
*     If the plane has been already included, it just
*     computes again the coordinates of the x-y couple
*     using improved PFAs
*     -------------------------------------------------
*     |||||||||||||||||||||||||||||||||||||||||||||||||
         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
            
            id=CP_STORE(nplanes-ip+1,ibest)
            
            is=is_cp(id)
            icp=icp_cp(id)
            if(ip_cp(id).ne.ip)
     $           print*,'OKKIO!!'
     $           ,'id ',id,is,icp
     $           ,ip_cp(id),ip
            icx=clx(ip,icp)
            icy=cly(ip,icp)
c            call xyz_PAM(icx,icy,is,
c     $           PFA,PFA,
c     $           AXV_STORE(nplanes-ip+1,ibest),
c     $           AYV_STORE(nplanes-ip+1,ibest))
            call xyz_PAM(icx,icy,is,
     $           PFA,PFA,
     $           AXV_STORE(nplanes-ip+1,ibest),
     $           AYV_STORE(nplanes-ip+1,ibest),
     $           bxyz(1),
     $           bxyz(2)
     $           )

            xm(nplanes-ip+1) = xPAM
            ym(nplanes-ip+1) = yPAM
            zm(nplanes-ip+1) = zPAM
            xgood(nplanes-ip+1) = 1
            ygood(nplanes-ip+1) = 1
            resx(nplanes-ip+1) = resxPAM
            resy(nplanes-ip+1) = resyPAM

            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) 
            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
            
*     |||||||||||||||||||||||||||||||||||||||||||||||||
*     -------------------------------------------------
*     If the plane has NOT  been already included, 
*     it tries to include a COUPLE or a single cluster
*     -------------------------------------------------
*     |||||||||||||||||||||||||||||||||||||||||||||||||
         else                   
               
            xgood(nplanes-ip+1)=0
            ygood(nplanes-ip+1)=0
               
*     --------------------------------------------------------------
*     determine which ladder and sensor are intersected by the track 
            call whichsensor(ip,xP,yP,nldt,ist)
*     if the track hit the plane in a dead area, go to the next plane
            if(nldt.eq.0.or.ist.eq.0)goto 133 

            SENSOR_STORE(nplanes-ip+1,IBEST)=ist
            LADDER_STORE(nplanes-ip+1,IBEST)=nldt
*     --------------------------------------------------------------

            if(DEBUG)then
               print*,
     $              '------ Plane ',ip,' intersected on LADDER ',nldt
     $              ,' SENSOR ',ist
               print*,
     $              '------ coord: ',XP,YP
            endif
            
*     ===========================================
*     STEP 1 >>>>>>>  try to include a new couple
*     ===========================================
c            if(DEBUG)print*,'>>>> try to include a new couple'
            distmin=1000000.
            xmm = 0.
            ymm = 0.
            zmm = 0.
            rxmm = 0.
            rymm = 0.
            dedxmmx = 0.        !(1)
            dedxmmy = 0.        !(1)
            idm = 0             !ID of the closer couple
            distance=0.
            do icp=1,ncp_plane(ip) !loop on couples on plane icp
               icx=clx(ip,icp)
               icy=cly(ip,icp)
               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
     $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)
     $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
     $              .false.)goto 1188 !then jump to next couple.
*          
               call xyz_PAM(icx,icy,ist,
     $              PFA,PFA,
     $              AXV_STORE(nplanes-ip+1,ibest),
     $              AYV_STORE(nplanes-ip+1,ibest),
     $              bxyz(1),
     $              bxyz(2)
     $              )
               
               distance = distance_to(XP,YP)
c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
               id=id_cp(ip,icp,ist)
               if(DEBUG)print*,'( couple ',id
     $              ,' ) distance ',distance
               if(distance.lt.distmin)then
                  xmm = xPAM
                  ymm = yPAM
                  zmm = zPAM
                  rxmm = resxPAM
                  rymm = resyPAM
                  distmin = distance
                  idm = id                  
                  dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
                  dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
                  clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
     $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
               endif
 1188          continue
            enddo               !end loop on couples on plane icp
c            if(distmin.le.clinc)then     !QUIQUI              
            if(distmin.le.clincnewc)then     !QUIQUI              
*              -----------------------------------
               xm(nplanes-ip+1) = xmm !<<<
               ym(nplanes-ip+1) = ymm !<<<
               zm(nplanes-ip+1) = zmm !<<<
               xgood(nplanes-ip+1) = 1 !<<<
               ygood(nplanes-ip+1) = 1 !<<<
               resx(nplanes-ip+1)=rxmm !<<<
               resy(nplanes-ip+1)=rymm !<<<
               dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
               dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
*              -----------------------------------
               CP_STORE(nplanes-ip+1,ibest)=idm      
               if(DEBUG)print*,'%%%% included couple ',idm
     $              ,' (dist.= ',distmin,', cut ',clinc,' )'
               goto 133         !next plane
            endif
*     ================================================
*     STEP 2 >>>>>>>  try to include a single cluster
*                     either from a couple or single
*     ================================================
c            if(DEBUG)print*,'>>>> try to include a new cluster'
            distmin=1000000.
            xmm_A = 0.          !---------------------------
            ymm_A = 0.          ! init variables that 
            zmm_A = 0.          ! define the SINGLET
            xmm_B = 0.          !
            ymm_B = 0.          !
            zmm_B = 0.          !
            rxmm = 0.           !
            rymm = 0.           !
            dedxmmx = 0.        !(1)
            dedxmmy = 0.        !(1)
            iclm=0              !---------------------------
            distance=0.

*----- clusters inside couples -------------------------------------     
            do icp=1,ncp_plane(ip) !loop on cluster inside couples
               icx=clx(ip,icp)
               icy=cly(ip,icp)
               id=id_cp(ip,icp,ist)
               if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
*                                                !jump to the next couple
*----- try cluster x -----------------------------------------------
c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used 
               if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
*                                              !jump to the Y cluster 
c               call xyz_PAM(icx,0,ist,
c     $              PFA,PFA,
c     $              AXV_STORE(nplanes-ip+1,ibest),0.)               
               call xyz_PAM(icx,0,ist,
     $              PFA,PFA,
     $              AXV_STORE(nplanes-ip+1,ibest),0.,
     $              bxyz(1),
     $              bxyz(2)
     $              )               
               distance = distance_to(XP,YP)
c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
               if(DEBUG)print*,'( cl-X ',icx
     $              ,' in cp ',id,' ) distance ',distance
               if(distance.lt.distmin)then
                  xmm_A = xPAM_A
                  ymm_A = yPAM_A
                  zmm_A = zPAM_A
                  xmm_B = xPAM_B
                  ymm_B = yPAM_B
                  zmm_B = zPAM_B
                  rxmm = resxPAM
                  rymm = resyPAM
                  distmin = distance
                  iclm = icx
c                  dedxmm = sgnl(icx) !(1)
                  dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
                  dedxmmy = 0.        !(1)
               endif                  
11881          continue
*----- try cluster y -----------------------------------------------
c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used 
               if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
*                                              !jump to the next couple
c               call xyz_PAM(0,icy,ist,
c     $              PFA,PFA,
c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
               call xyz_PAM(0,icy,ist,
     $              PFA,PFA,
     $              0.,AYV_STORE(nplanes-ip+1,ibest),
     $              bxyz(1),
     $              bxyz(2)
     $              )
               distance = distance_to(XP,YP)
c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
               if(DEBUG)print*,'( cl-Y ',icy
     $              ,' in cp ',id,' ) distance ',distance
               if(distance.lt.distmin)then
                  xmm_A = xPAM_A
                  ymm_A = yPAM_A
                  zmm_A = zPAM_A
                  xmm_B = xPAM_B
                  ymm_B = yPAM_B
                  zmm_B = zPAM_B
                  rxmm = resxPAM
                  rymm = resyPAM
                  distmin = distance
                  iclm = icy
c                 dedxmm = sgnl(icy)  !(1)
                  dedxmmx = 0.        !(1)
                  dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
               endif                  
11882          continue
            enddo               !end loop on cluster inside couples
*----- single clusters -----------------------------------------------   
c            print*,'## ncls(',ip,') ',ncls(ip)
            do ic=1,ncls(ip)    !loop on single clusters
               icl=cls(ip,ic)
c               if(cl_used(icl).eq.1.or.     !if the cluster is already used 
               if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
     $              LADDER(icl).ne.nldt.or. !or the ladder number does not match 
     $              .false.)goto 18882      !jump to the next singlet
               if(mod(VIEW(icl),2).eq.0)then!<---- X view
                  call xyz_PAM(icl,0,ist,
     $                 PFA,PFA,
     $                 AXV_STORE(nplanes-ip+1,ibest),0.,
     $                 bxyz(1),
     $                 bxyz(2)
     $                 )
               else                         !<---- Y view
                  call xyz_PAM(0,icl,ist,
     $                 PFA,PFA,
     $                 0.,AYV_STORE(nplanes-ip+1,ibest),
     $                 bxyz(1),
     $                 bxyz(2)
     $                 )
               endif

               distance = distance_to(XP,YP)
c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
               if(DEBUG)print*,'( cl-s ',icl
     $              ,' ) distance ',distance,'<',distmin,' ?'
               if(distance.lt.distmin)then
                  if(DEBUG)print*,'YES'
                  xmm_A = xPAM_A
                  ymm_A = yPAM_A
                  zmm_A = zPAM_A
                  xmm_B = xPAM_B
                  ymm_B = yPAM_B
                  zmm_B = zPAM_B
                  rxmm = resxPAM
                  rymm = resyPAM
                  distmin = distance   
                  iclm = icl
                  if(mod(VIEW(icl),2).eq.0)then !<---- X view
                     dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) 
                     dedxmmy = 0.                  
                  else          !<---- Y view
                     dedxmmx = 0.                   
                     dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) 
                  endif
               endif                  
18882          continue
            enddo               !end loop on single clusters
c            print*,'## distmin ', distmin,' clinc ',clinc

c     QUIQUI------------
c     anche qui: non ci vuole clinc???
            if(iclm.ne.0)then
               if(mod(VIEW(iclm),2).eq.0)then
                  clincnew=
     $                 20*
     $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
               else if(mod(VIEW(iclm),2).ne.0)then
                  clincnew=
     $                 10*
     $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
               endif
c     QUIQUI------------
               
               if(distmin.le.clincnew)then   !QUIQUI
c     if(distmin.le.clinc)then          !QUIQUI           
                  
                  CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<     
*     ----------------------------
c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
                  if(mod(VIEW(iclm),2).eq.0)then
                     XGOOD(nplanes-ip+1)=1.
                     resx(nplanes-ip+1)=rxmm
                     if(DEBUG)print*,'%%%% included X-cl ',iclm
     $                    ,'( chi^2, ',RCHI2_STORE(ibest)
     $                    ,', dist.= ',distmin
     $                    ,', cut ',clinc,' )'
                  else
                     YGOOD(nplanes-ip+1)=1.
                     resy(nplanes-ip+1)=rymm
                     if(DEBUG)print*,'%%%% included Y-cl ',iclm
     $                    ,'( chi^2, ',RCHI2_STORE(ibest)
     $                    ,', dist.= ', distmin
     $                    ,', cut ',clinc,' )'
                  endif
c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
*     ----------------------------
                  xm_A(nplanes-ip+1) = xmm_A
                  ym_A(nplanes-ip+1) = ymm_A
                  xm_B(nplanes-ip+1) = xmm_B
                  ym_B(nplanes-ip+1) = ymm_B
                  zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
                  dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
                  dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
*     ----------------------------
               endif
            endif
         endif
 133     continue
      enddo                     !end loop on planes
      
      
      return
      end

***************************************************
*                                                 *
*                                                 *
*                                                 *
*                                                 *
*                                                 *
*                                                 *
**************************************************
*
      subroutine clean_XYclouds(ibest,iflag)

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'level2.f'       

      do ip=1,nplanes           !loop on planes

         id=CP_STORE(nplanes-ip+1,ibest)
         icl=CLS_STORE(nplanes-ip+1,ibest)
         if(id.ne.0.or.icl.ne.0)then               
            if(id.ne.0)then
               iclx=clx(ip,icp_cp(id))
               icly=cly(ip,icp_cp(id))
               cl_used(iclx)=ntrk  !tag used clusters 
               cl_used(icly)=ntrk  !tag used clusters 
            elseif(icl.ne.0)then
               cl_used(icl)=ntrk   !tag used clusters 
            endif
            
*     ----------------------------- 
*     remove the couple from clouds
*     remove also vitual couples containing the 
*     selected clusters
*     ----------------------------- 
            do icp=1,ncp_plane(ip)
               if(
     $              clx(ip,icp).eq.iclx
     $              .or.
     $              clx(ip,icp).eq.icl
     $              .or.
     $              cly(ip,icp).eq.icly
     $              .or.
     $              cly(ip,icp).eq.icl
     $              )then
                  id=id_cp(ip,icp,1)
                  if(DEBUG)then
                     print*,ip,' <<< cp ',id
     $                    ,' ( cl-x '
     $                    ,clx(ip,icp)
     $                    ,' cl-y '
     $                    ,cly(ip,icp),' ) --> removed'
                  endif
*     ----------------------------- 
*     remove the couple from clouds
                  do iyz=1,nclouds_yz
                     if(cpcloud_yz(iyz,abs(id)).ne.0)then
                        ptcloud_yz(iyz)=ptcloud_yz(iyz)-1
                        cpcloud_yz(iyz,abs(id))=0
                     endif
                  enddo
                  do ixz=1,nclouds_xz
                     if(cpcloud_xz(ixz,abs(id)).ne.0)then
                        ptcloud_xz(ixz)=ptcloud_xz(ixz)-1
                        cpcloud_xz(ixz,abs(id))=0
                     endif
                  enddo                     
*     ----------------------------- 
               endif
            enddo
            
         endif               
      enddo                     !end loop on planes
      
      return
      end


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

      subroutine init_level2

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'level2.f'

*     ---------------------------------
*     variables initialized from level1
*     ---------------------------------
      do i=1,nviews
         good2(i)=good1(i)
         do j=1,nva1_view
            vkflag(i,j)=1
            if(cnnev(i,j).le.0)then 
               vkflag(i,j)=cnnev(i,j)
            endif
         enddo
      enddo
*     ----------------
*     level2 variables
*     ----------------
      NTRK = 0
      do it=1,NTRKMAX
         IMAGE(IT)=0
         CHI2_nt(IT) = -100000.
         do ip=1,nplanes
            XM_nt(IP,IT) = 0
            YM_nt(IP,IT) = 0
            ZM_nt(IP,IT) = 0
            RESX_nt(IP,IT) = 0
            RESY_nt(IP,IT) = 0
            TAILX_nt(IP,IT) = 0
            TAILY_nt(IP,IT) = 0
            XBAD(IP,IT) = 0
            YBAD(IP,IT) = 0
            XGOOD_nt(IP,IT) = 0
            YGOOD_nt(IP,IT) = 0
            LS(IP,IT) = 0
            DEDX_X(IP,IT) = 0
            DEDX_Y(IP,IT) = 0
            CLTRX(IP,IT) = 0
            CLTRY(IP,IT) = 0
         enddo
         do ipa=1,5
            AL_nt(IPA,IT) = 0
            do ipaa=1,5
               coval(ipa,ipaa,IT)=0
            enddo                  
         enddo                  
      enddo
      nclsx=0
      nclsy=0      
      do ip=1,NSINGMAX
        planex(ip)=0
        xs(1,ip)=0
        xs(2,ip)=0
        sgnlxs(ip)=0
        planey(ip)=0
        ys(1,ip)=0
        ys(2,ip)=0
        sgnlys(ip)=0
      enddo
      end


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


      subroutine init_hough

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f' 
      include 'common_hough.f'
      include 'level2.f'

      ntrpt_nt=0
      ndblt_nt=0
      NCLOUDS_XZ_nt=0
      NCLOUDS_YZ_nt=0
      do idb=1,ndblt_max_nt
         db_cloud_nt(idb)=0
         alfayz1_nt(idb)=0      
         alfayz2_nt(idb)=0      
      enddo
      do itr=1,ntrpt_max_nt
         tr_cloud_nt(itr)=0
         alfaxz1_nt(itr)=0      
         alfaxz2_nt(itr)=0      
         alfaxz3_nt(itr)=0      
      enddo
      do idb=1,ncloyz_max       
        ptcloud_yz_nt(idb)=0    
        alfayz1_av_nt(idb)=0    
        alfayz2_av_nt(idb)=0    
      enddo                     
      do itr=1,ncloxz_max       
        ptcloud_xz_nt(itr)=0    
        alfaxz1_av_nt(itr)=0    
        alfaxz2_av_nt(itr)=0    
        alfaxz3_av_nt(itr)=0    
      enddo                     

      ntrpt=0                   
      ndblt=0                   
      NCLOUDS_XZ=0              
      NCLOUDS_YZ=0              
      do idb=1,ndblt_max        
        db_cloud(idb)=0         
        cpyz1(idb)=0            
        cpyz2(idb)=0            
        alfayz1(idb)=0          
        alfayz2(idb)=0          
      enddo                     
      do itr=1,ntrpt_max        
        tr_cloud(itr)=0         
        cpxz1(itr)=0            
        cpxz2(itr)=0            
        cpxz3(itr)=0            
        alfaxz1(itr)=0          
        alfaxz2(itr)=0          
        alfaxz3(itr)=0          
      enddo                     
      do idb=1,ncloyz_max       
        ptcloud_yz(idb)=0       
        alfayz1_av(idb)=0       
        alfayz2_av(idb)=0       
        do idbb=1,ncouplemaxtot 
          cpcloud_yz(idb,idbb)=0 
        enddo                   
      enddo                     
      do itr=1,ncloxz_max       
        ptcloud_xz(itr)=0       
        alfaxz1_av(itr)=0       
        alfaxz2_av(itr)=0       
        alfaxz3_av(itr)=0       
        do itrr=1,ncouplemaxtot 
          cpcloud_xz(itr,itrr)=0 
        enddo                   
      enddo                     
      end
************************************************************
*
*
*
*
*
*
*
************************************************************


      subroutine fill_level2_tracks(ntr)

*     -------------------------------------------------------
*     This routine fills the ntr-th element of the variables
*     inside the level2_tracks common, which correspond
*     to the ntr-th track info.
*     -------------------------------------------------------

     
      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'level2.f'
      include 'common_mini_2.f'
      include 'calib.f'

      character*10 PFA
      common/FINALPFA/PFA

      real sinth,phi,pig
      integer ssensor,sladder
      pig=acos(-1.)

*     -------------------------------------
      chi2_nt(ntr)        = sngl(chi2)
      nstep_nt(ntr)       = nstep
*     -------------------------------------
      phi   = al(4)           
      sinth = al(3)            
      if(sinth.lt.0)then       
         sinth = -sinth        
         phi = phi + pig       
      endif                    
      npig = aint(phi/(2*pig)) 
      phi = phi - npig*2*pig   
      if(phi.lt.0)             
     $     phi = phi + 2*pig   
      al(4) = phi              
      al(3) = sinth            
      do i=1,5
         al_nt(i,ntr)     = sngl(al(i))
         do j=1,5
            coval(i,j,ntr) = sngl(cov(i,j))
         enddo
      enddo
*     -------------------------------------      
      do ip=1,nplanes           ! loop on planes
         xgood_nt(ip,ntr) = int(xgood(ip))
         ygood_nt(ip,ntr) = int(ygood(ip))
         xm_nt(ip,ntr)    = sngl(xm(ip))
         ym_nt(ip,ntr)    = sngl(ym(ip))
         zm_nt(ip,ntr)    = sngl(zm(ip))
         RESX_nt(IP,ntr)  = sngl(resx(ip))
         RESY_nt(IP,ntr)  = sngl(resy(ip))
         TAILX_nt(IP,ntr) = 0.
         TAILY_nt(IP,ntr) = 0.
         xv_nt(ip,ntr)    = sngl(xv(ip))
         yv_nt(ip,ntr)    = sngl(yv(ip))
         zv_nt(ip,ntr)    = sngl(zv(ip))
         axv_nt(ip,ntr)   = sngl(axv(ip))
         ayv_nt(ip,ntr)   = sngl(ayv(ip))   
c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         factor = sqrt( 
     $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
     $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 + 
     $        1. )
c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor) 
         dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
   
         ax   = axv_nt(ip,ntr)
         ay   = ayv_nt(ip,ntr)
         bfx  = BX_STORE(ip,IDCAND)
         bfy  = BY_STORE(ip,IDCAND)
         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
         angx     = 180.*atan(tgtemp)/acos(-1.)
         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001         
         angy    = 180.*atan(tgtemp)/acos(-1.)
         
c         print*,'* ',ip,bfx,bfy,angx,angy

         id  = CP_STORE(ip,IDCAND) ! couple id
         icl = CLS_STORE(ip,IDCAND)
         ssensor = -1
         sladder = -1
         ssensor = SENSOR_STORE(ip,IDCAND)
         sladder = LADDER_STORE(ip,IDCAND)
         if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
         LS(IP,ntr)      = ssensor+10*sladder 

         if(id.ne.0)then
c           >>> is a couple
            cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
            cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))

c$$$            if(is_cp(id).ne.ssensor)
c$$$     $           print*,'ERROR is sensor assignment (couple)'
c$$$     $           ,is_cp(id),ssensor
c$$$            if(LADDER(clx(nplanes-ip+1,icp_cp(id))).ne.sladder)
c$$$     $           print*,'ERROR is ladder assignment (couple)'
c$$$     $           ,LADDER(clx(nplanes-ip+1,icp_cp(id))),sladder
            
            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
            ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))

            if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
     $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
            if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
     $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)

         elseif(icl.ne.0)then
c           >>> is a singlet
c$$$            if(LADDER(icl).ne.sladder)
c$$$     $           print*,'ERROR is ladder assignment (single)'
c$$$     $           ,LADDER(icl),sladder
            if(mod(VIEW(icl),2).eq.0)then 
               cltrx(ip,ntr)=icl
               nnnnn = npfastrips(icl,PFA,angx)
               xbad(ip,ntr) = nbadstrips(nnnnn,icl) 
               if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
            elseif(mod(VIEW(icl),2).eq.1)then
               cltry(ip,ntr)=icl
               nnnnn = npfastrips(icl,PFA,angy)
               ybad(ip,ntr) = nbadstrips(nnnnn,icl) 
               if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
            endif
         endif          

      enddo


c$$$      print*,(xgood(i),i=1,6)
c$$$      print*,(ygood(i),i=1,6)
c$$$      print*,(ls(i,ntr),i=1,6)
c$$$      print*,(dedx_x(i,ntr),i=1,6)
c$$$      print*,(dedx_y(i,ntr),i=1,6)
c$$$      print*,'-----------------------'

      end

      subroutine fill_level2_siglets

*     -------------------------------------------------------
*     This routine fills the  elements of the variables
*     inside the level2_singletsx and level2_singletsy commons, 
*     which store info on clusters outside the tracks
*     -------------------------------------------------------

      include 'commontracker.f'
      include 'calib.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'level2.f'
      include 'common_xyzPAM.f'

*     count #cluster per plane not associated to any track
      nclsx = 0
      nclsy = 0

      do iv = 1,nviews
         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
      enddo

      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)then !=== X views
               nclsx = nclsx + 1
               planex(nclsx) = ip
               sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
               if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
               clsx(nclsx)   = icl
               do is=1,2
c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
                  xs(is,nclsx) = (xPAM_A+xPAM_B)/2
               enddo
c$$$               print*,'nclsx         ',nclsx
c$$$               print*,'planex(nclsx) ',planex(nclsx)
c$$$               print*,'sgnlxs(nclsx) ',sgnlxs(nclsx)
c$$$               print*,'xs(1,nclsx)   ',xs(1,nclsx)
c$$$               print*,'xs(2,nclsx)   ',xs(2,nclsx)
            else                          !=== Y views
               nclsy = nclsy + 1
               planey(nclsy) = ip
               sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
               if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
               clsy(nclsy)   = icl
               do is=1,2
c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
                  ys(is,nclsy) = (yPAM_A+yPAM_B)/2
               enddo
c$$$               print*,'nclsy         ',nclsy
c$$$               print*,'planey(nclsy) ',planey(nclsy)
c$$$               print*,'sgnlys(nclsy) ',sgnlys(nclsy)
c$$$               print*,'ys(1,nclsy)   ',ys(1,nclsy)
c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
            endif
         endif

***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
         whichtrack(icl) = cl_used(icl)

      enddo
      end

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

      subroutine fill_hough

*     -------------------------------------------------------
*     This routine fills the  variables related to the hough 
*     transform, for the debig n-tuple
*     -------------------------------------------------------

      include 'commontracker.f'
      include 'level1.f'
      include 'common_momanhough.f'
      include 'common_hough.f'
      include 'level2.f'

      if(.false.
     $     .or.ntrpt.gt.ntrpt_max_nt
     $     .or.ndblt.gt.ndblt_max_nt
     $     .or.NCLOUDS_XZ.gt.ncloxz_max
     $     .or.NCLOUDS_yZ.gt.ncloyz_max
     $     )then
         ntrpt_nt=0
         ndblt_nt=0
         NCLOUDS_XZ_nt=0
         NCLOUDS_YZ_nt=0
      else
         ndblt_nt=ndblt
         ntrpt_nt=ntrpt
         if(ndblt.ne.0)then
            do id=1,ndblt
               alfayz1_nt(id)=alfayz1(id) !Y0
               alfayz2_nt(id)=alfayz2(id) !tg theta-yz
            enddo
         endif
         if(ndblt.ne.0)then
            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
            enddo
         endif
         nclouds_yz_nt=nclouds_yz
         nclouds_xz_nt=nclouds_xz
         if(nclouds_yz.ne.0)then
            nnn=0
            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)
               nnn=nnn+ptcloud_yz(iyz)
            enddo
            do ipt=1,nnn
               db_cloud_nt(ipt)=db_cloud(ipt)
             enddo
         endif
         if(nclouds_xz.ne.0)then
            nnn=0
            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)
               nnn=nnn+ptcloud_xz(ixz)               
            enddo
            do ipt=1,nnn
              tr_cloud_nt(ipt)=tr_cloud(ipt)
             enddo
         endif
      endif
      end
      