/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/analysissubroutines.f

Diff of /DarthVader/TrackerLevel2/src/F77/analysissubroutines.f

Parent Directory | Revision Log | View Patch Patch

-revision 1.27 by pam-fi,
Fri Aug 17 14:36:05 2007 UTC
+revision 1.40 by mocchiut,
Tue Aug  4 14:01:35 2009 UTC
 Line 20
        include 'calib.f'
        include 'level2.f'
  c      print*,'======================================================'
  c$$$      do ic=1,NCLSTR1
  c$$$         if(.false.
-Line 74 
 c$$$      enddo
+Line 72 
 c$$$      enddo
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
        call cl_to_couples(iflag)
        if(iflag.eq.1)then        !bad event
           goto 880               !go to next event
        endif
+       if(ncp_tot.eq.0)goto 880  !go to next event
  *-----------------------------------------------------
  *-----------------------------------------------------
  *     HOUGH TRASFORM
-Line 113 
 c$$$      enddo
+Line 110 
 c$$$      enddo
        if(iflag.eq.1)then        !bad event
           goto 880               !go to next event
        endif
+       if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event
  *-------------------------------------------------------------------------------
-Line 163 
 c$$$      enddo
+Line 161 
 c$$$      enddo
        if(iflag.eq.1)then        !bad event
           goto 880               !fill ntp and go to next event
        endif
+ *     ------------------------------------------------
+ *     first try to release the tolerance
+ *     ------------------------------------------------
        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
-         if(cutdistyz.lt.maxcuty/2)then
+          if(cutdistyz.lt.maxcuty/2)then
-           cutdistyz=cutdistyz+cutystep
+             cutdistyz=cutdistyz+cutystep
-         else
+          else
-           cutdistyz=cutdistyz+(3*cutystep)
+             cutdistyz=cutdistyz+(3*cutystep)
-         endif
+          endif
-         goto 878
+          if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
-       endif
+          goto 878
+       endif
+ *     ------------------------------------------------
+ *     hence reduce the minimum number of plane
+ *     ------------------------------------------------
+       if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then
+          nplyz_min=nplyz_min-1
+          if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min
+          goto 878
+       endif
+       if(nclouds_yz.eq.0)goto 880 !go to next event
-       if(planehit.eq.3) goto 881
+ ccc   if(planehit.eq.3) goto 881
  continue
        call trip_to_XZcloud(iflag)
        if(iflag.eq.1)then        !bad event
           goto 880               !fill ntp and go to next event
        endif
+ *     ------------------------------------------------
+ *     first try to release the tolerance
+ *     ------------------------------------------------
        if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
          cutdistxz=cutdistxz+cutxstep
+          if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz
          goto 879
        endif
+ *     ------------------------------------------------
+ *     hence reduce the minimum number of plane
+ *     ------------------------------------------------
+       if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then
+          nplxz_min=nplxz_min-1
+          if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min
+          goto 879
+       endif
+       if(nclouds_xz.eq.0)goto 880 !go to next event
- continue
+ c$$$ 881  continue
- *     if there is at least three planes on the Y view decreases cuts on X view
+ c$$$*     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.
+ c$$$      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
-      $     nplxz_min.ne.y_min_start)then
+ c$$$     $     nplxz_min.ne.y_min_start)then
-         nptxz_min=x_min_step
+ c$$$        nptxz_min=x_min_step
-         nplxz_min=x_min_start-x_min_step
+ c$$$        nplxz_min=x_min_start-x_min_step
-         goto 879
+ c$$$        goto 879
-       endif
+ c$$$      endif
  return
        end
-Line 215 
 c$$$      enddo
+Line 241 
 c$$$      enddo
  c      include 'momanhough_init.f'
        logical FIMAGE            !
+       real trackimage(NTRACKSMAX)
        real*8 AL_GUESS(5)
  *-------------------------------------------------------------------------------
-Line 242 
 c      include 'momanhough_init.f'
+Line 269 
 c      include 'momanhough_init.f'
  *
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
-          ntrk=0                 !counter of identified physical tracks
+ ccc         ntrk=0                 !counter of identified physical tracks
-   continue               !<<<<<<< come here when performing a new search
+ c11111    continue               !<<<<<<< come here when performing a new search
+          continue                  !<<<<<<< come here when performing a new search
+          if(nclouds_xz.eq.0)goto 880 !go to next event
+          if(nclouds_yz.eq.0)goto 880 !go to next event
  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
+          if(ntracks.eq.0)goto 880 !go to next event
           FIMAGE=.false.         !processing best track (not track image)
           ibest=0                !best track among candidates
-Line 294 
 c$$$         if(ibest.eq.0)goto 880 !>>
+Line 326 
 c$$$         if(ibest.eq.0)goto 880 !>>
             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
  *-------------------------------------------------------------------------------
-Line 331 
 c$$$         enddo
+Line 345 
 c$$$         enddo
              iimage=0
           endif
           if(icand.eq.0)then
-             if(VERBOSE)then
+             if(VERBOSE.EQ.1)then
                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
       $              ,ibest,iimage
              endif
-Line 349 
 c$$$         enddo
+Line 363 
 c$$$         enddo
           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))
-Line 360 
 c         print*,'## guess: ',al
+Line 373 
 c         print*,'## guess: ',al
           jstep=0                !# minimization steps
           iprint=0
- c         if(DEBUG)iprint=1
+ c         if(DEBUG.EQ.1)iprint=1
-          if(VERBOSE)iprint=1
+          if(VERBOSE.EQ.1)iprint=1
-          if(DEBUG)iprint=2
+          if(DEBUG.EQ.1)iprint=2
           call mini2(jstep,ifail,iprint)
-          if(ifail.ne.0) then
+ cc         if(ifail.ne.0) then
-             if(VERBOSE)then
+          if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
+             if(CHI2.ne.CHI2)CHI2=-9999. !new
+             if(VERBOSE.EQ.1)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
+          if(DEBUG.EQ.1)then
              print*,'----------------------------- improved track coord'
       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
              do ip=1,6
-Line 391 
 c            chi2=-chi2
+Line 398 
 c            chi2=-chi2
           endif
  c         rchi2=chi2/dble(ndof)
-          if(DEBUG)then
+          if(DEBUG.EQ.1)then
              print*,' '
              print*,'****** SELECTED TRACK *************'
              print*,'#         R. chi2        RIG'
-Line 407 
 c         rchi2=chi2/dble(ndof)
+Line 414 
 c         rchi2=chi2/dble(ndof)
  *     ------------- 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
+ *     -----------------------------------------------------
+ *     first check if the track is ambiguous
+ *     -----------------------------------------------------
+ *     (modified on august 2007 by ElenaV)
+          is1=0
+          do ip=1,NPLANES
+             if(SENSOR_STORE(ip,icand).ne.0)then
+                is1=SENSOR_STORE(ip,icand)
+                if(ip.eq.6)is1=3-is1 !last plane inverted
+             endif
+          enddo
+          if(is1.eq.0)then
+             if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
+             goto 122            !jump
+          endif
+          do ip=1,NPLANES
+             is2 = SENSOR_STORE(ip,icand) !sensor
+             if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
+             if(
+      $           (is1.ne.is2.and.is2.ne.0)
+      $           )goto 122      !jump (this track cannot have an image)
+          enddo
+          if(DEBUG.eq.1)print*,' >>> ambiguous track! '
+ *     ---------------------------------------------------------------
+ *     take the candidate with the greatest number of matching couples
+ *     if more than one satisfies the requirement,
+ *     choose the one with more points and lower chi2
+ *     ---------------------------------------------------------------
+ *     count the number of matching couples
+          ncpmax = 0
+          iimage   = 0           !id of candidate with better matching
           do i=1,ntracks
-             iimage=i
+             ncp=0
              do ip=1,nplanes
-                if(     CP_STORE(nplanes-ip+1,icand).ne.
+                if(CP_STORE(nplanes-ip+1,icand).ne.0)then
-      $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0
+                   if(
+      $                 CP_STORE(nplanes-ip+1,i).ne.0
+      $                 .and.
+      $                 CP_STORE(nplanes-ip+1,icand).eq.
+      $                 -1*CP_STORE(nplanes-ip+1,i)
+      $                 )then
+                      ncp=ncp+1  !ok
+                   elseif(
+      $                    CP_STORE(nplanes-ip+1,i).ne.0
+      $                    .and.
+      $                    CP_STORE(nplanes-ip+1,icand).ne.
+      $                    -1*CP_STORE(nplanes-ip+1,i)
+      $                    )then
+                      ncp = 0
+                      goto 117   !it is not an image candidate
+                   else
+                   endif
+                endif
+             enddo
+       continue
+             trackimage(i)=ncp   !number of matching couples
+             if(ncp>ncpmax)then
+                ncpmax=ncp
+                iimage=i
+             endif
+          enddo
+ *     check if there are more than one image candidates
+          ngood=0
+          do i=1,ntracks
+             if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
+          enddo
+          if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood
+ *     if there are, choose the best one
+          if(ngood.gt.0)then
+ *     -------------------------------------------------------
+ *     order track-candidates according to:
+ *     1st) decreasing n.points
+ *     2nd) increasing chi**2
+ *     -------------------------------------------------------
+             rchi2best=1000000000.
+             ndofbest=0
+             do i=1,ntracks
+                if( trackimage(i).eq.ncpmax )then
+                   ndof=0
+                   do ii=1,nplanes
+                      ndof=ndof
+      $                    +int(xgood_store(ii,i))
+      $                    +int(ygood_store(ii,i))
+                   enddo
+                   if(ndof.gt.ndofbest)then
+                      iimage=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
+                         iimage=i
+                         rchi2best=RCHI2_STORE(i)
+                         ndofbest=ndof
+                      endif
+                   endif
+                endif
              enddo
-             if(  iimage.ne.0.and.
- c     $           RCHI2_STORE(i).le.CHI2MAX.and.
+             if(DEBUG.EQ.1)then
- c     $           RCHI2_STORE(i).gt.0.and.
+                print*,'Track candidate ',iimage
-      $           .true.)then
-                if(DEBUG)print*,'Track candidate ',iimage
       $              ,' >>> TRACK IMAGE >>> of'
       $              ,ibest
-                goto 122         !image track found
              endif
-          enddo
+          endif
     continue
  *     --- and store the results --------------------------------
           ntrk = ntrk + 1                   !counter of found tracks
           if(.not.FIMAGE
-Line 434 
 c     $           RCHI2_STORE(i).gt.0.an
+Line 535 
 c     $           RCHI2_STORE(i).gt.0.an
       $        .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)
+             if(verbose.eq.1)
       $           print*,
       $           '** warning ** number of identified '//
       $           'tracks exceeds vector dimension '
-Line 451 
 cc            good2=.false.
+Line 550 
 cc            good2=.false.
              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
- *     **********************************************
     return
        end
-Line 572 
 c     $        rchi2best.lt.15..and.
+Line 639 
 c     $        rchi2best.lt.15..and.
        real stripx,stripy
+       double precision xi,yi,zi
+       double precision xi_A,yi_A,zi_A
+       double precision xi_B,yi_B,zi_B
        double precision xrt,yrt,zrt
        double precision xrt_A,yrt_A,zrt_A
        double precision xrt_B,yrt_B,zrt_B
-Line 580 
 c     $        rchi2best.lt.15..and.
+Line 650 
 c     $        rchi2best.lt.15..and.
        parameter (ndivx=30)
- c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
        resxPAM = 0
        resyPAM = 0
-       xPAM = 0.
+       xPAM = 0.D0
-       yPAM = 0.
+       yPAM = 0.D0
-       zPAM = 0.
+       zPAM = 0.D0
-       xPAM_A = 0.
+       xPAM_A = 0.D0
-       yPAM_A = 0.
+       yPAM_A = 0.D0
-       zPAM_A = 0.
+       zPAM_A = 0.D0
-       xPAM_B = 0.
+       xPAM_B = 0.D0
-       yPAM_B = 0.
+       yPAM_B = 0.D0
-       zPAM_B = 0.
+       zPAM_B = 0.D0
- c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
        if(sensor.lt.1.or.sensor.gt.2)then
           print*,'xyz_PAM   ***ERROR*** wrong input '
-Line 611 
 c      print*,'## xyz_PAM: ',icx,icy,sen
+Line 679 
 c      print*,'## xyz_PAM: ',icx,icy,sen
           viewx   = VIEW(icx)
           nldx    = nld(MAXS(icx),VIEW(icx))
           nplx    = npl(VIEW(icx))
-          resxPAM = RESXAV
+ c         resxPAM = RESXAV
           stripx  = float(MAXS(icx))
           if(
-Line 631 
 c      print*,'## xyz_PAM: ',icx,icy,sen
+Line 699 
 c      print*,'## xyz_PAM: ',icx,icy,sen
  *        --------------------------
  *        magnetic-field corrections
  *        --------------------------
-          angtemp  = ax
+          stripx  = stripx +  fieldcorr(viewx,bfy)
-          bfytemp  = bfy
+          angx    = effectiveangle(ax,viewx,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)
+          call applypfa(PFAx,icx,angx,corr,res)
-          if(PFAx.eq.'COG1')then
+          stripx  = stripx + corr
+          resxPAM = res
-             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 = risxeta2(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 = risxeta3(abs(angx))
-             resxPAM = resxPAM*fbad_cog(3,icx)
- c            if(DEBUG.and.fbad_cog(3,icx).ne.1)
- c     $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
-          elseif(PFAx.eq.'ETA4')then
-             stripx  = stripx + pfaeta4(icx,angx)
-             resxPAM = risxeta4(abs(angx))
-             resxPAM = resxPAM*fbad_cog(4,icx)
- c            if(DEBUG.and.fbad_cog(4,icx).ne.1)
- c     $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
-          elseif(PFAx.eq.'ETA')then
-             stripx  = stripx + pfaeta(icx,angx)
- c            resxPAM = riseta(icx,angx)
-             resxPAM = riseta(viewx,angx)
-             resxPAM = resxPAM*fbad_eta(icx,angx)
- c            if(DEBUG.and.fbad_cog(2,icx).ne.1)
- c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
-          elseif(PFAx.eq.'ETAL')then
-             stripx  = stripx + pfaetal(icx,angx)
-             resxPAM = riseta(viewx,angx)
-             resxPAM = resxPAM*fbad_eta(icx,angx)
- c            if(DEBUG.and.fbad_cog(2,icx).ne.1)
- c     $           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
-Line 760 
 c$$$            print*,icx,' *** ',resxP
+Line 717 
 c$$$            print*,icx,' *** ',resxP
           viewy = VIEW(icy)
           nldy = nld(MAXS(icy),VIEW(icy))
           nply = npl(VIEW(icy))
-          resyPAM = RESYAV
+ c         resyPAM = RESYAV
           stripy = float(MAXS(icy))
           if(
-Line 778 
 c$$$            print*,icx,' *** ',resxP
+Line 735 
 c$$$            print*,icx,' *** ',resxP
           endif
           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
-             if(DEBUG) then
+             if(DEBUG.EQ.1) 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
+          stripy  = stripy +  fieldcorr(viewy,bfx)
-          angy    = 180.*atan(tgtemp)/acos(-1.)
+          angy    = effectiveangle(ay,viewy,bfx)
-          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 ---'
+          call applypfa(PFAy,icy,angy,corr,res)
- c$$$         istart = INDSTART(ICY)
+          stripy  = stripy + corr
- c$$$         istop  = TOTCLLENGTH
+          resyPAM = res
- 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 = risyeta2(abs(angy))
-             resyPAM = resyPAM*fbad_cog(2,icy)
- c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
- c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
-          elseif(PFAy.eq.'ETA3')then
-             stripy  = stripy + pfaeta3(icy,angy)
-             resyPAM = resyPAM*fbad_cog(3,icy)
- c            if(DEBUG.and.fbad_cog(3,icy).ne.1)
- c     $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
-          elseif(PFAy.eq.'ETA4')then
-             stripy  = stripy + pfaeta4(icy,angy)
-             resyPAM = resyPAM*fbad_cog(4,icy)
- c            if(DEBUG.and.fbad_cog(4,icy).ne.1)
- c     $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
-          elseif(PFAy.eq.'ETA')then
-             stripy  = stripy + pfaeta(icy,angy)
- c            resyPAM = riseta(icy,angy)
-             resyPAM = riseta(viewy,angy)
-             resyPAM = resyPAM*fbad_eta(icy,angy)
- c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
- c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
-          elseif(PFAy.eq.'ETAL')then
-             stripy  = stripy + pfaetal(icy,angy)
-             resyPAM = riseta(viewy,angy)
-             resyPAM = resyPAM*fbad_eta(icy,angy)
- c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
- c     $           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
-Line 903 
 c     (xi,yi,zi) = mechanical coordinate
+Line 765 
 c     (xi,yi,zi) = mechanical coordinate
  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
+             if(DEBUG.EQ.1) then
                 print*,'xyz_PAM (couple):',
       $              ' WARNING: false X strip: strip ',stripx
              endif
           endif
-          xi = acoordsi(stripx,viewx)
+          xi = dcoordsi(stripx,viewx)
-          yi = acoordsi(stripy,viewy)
+          yi = dcoordsi(stripy,viewy)
-          zi = 0.
+          zi = 0.D0
  c------------------------------------------------------------------------
  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-Line 946 
 c---------------------------------------
+Line 807 
 c---------------------------------------
           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
-          xPAM_A = 0.
+          xPAM_A = 0.D0
-          yPAM_A = 0.
+          yPAM_A = 0.D0
-          zPAM_A = 0.
+          zPAM_A = 0.D0
-          xPAM_B = 0.
+          xPAM_B = 0.D0
-          yPAM_B = 0.
+          yPAM_B = 0.D0
-          zPAM_B = 0.
+          zPAM_B = 0.D0
        elseif(
       $        (icx.ne.0.and.icy.eq.0).or.
-Line 972 
 C=======================================
+Line 833 
 C=======================================
              nldx = nldy
              viewx = viewy + 1
-             yi   = acoordsi(stripy,viewy)
+             xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center
+             yi = dcoordsi(stripy,viewy)
+             zi = 0.D0
              xi_A = edgeY_d - SiDimX/2
              yi_A = yi
-Line 982 
 C=======================================
+Line 845 
 C=======================================
              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===========================================================
-Line 994 
 C=======================================
+Line 855 
 C=======================================
              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
+                if(DEBUG.EQ.1) then
                    print*,'xyz_PAM (X-singlet):',
       $                 ' WARNING: false X strip: strip ',stripx
                 endif
              endif
-             xi   = acoordsi(stripx,viewx)
+             xi = dcoordsi(stripx,viewx)
+             yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
+             zi = 0.D0
              xi_A = xi
              yi_A = edgeX_d - SiDimY/2
-Line 1019 
 c            if((stripx.le.3).or.(stripx
+Line 881 
 c            if((stripx.le.3).or.(stripx
                 yi_B = yi
              endif
- c            print*,'X-cl ',icx,stripx,' --> ',xi
- c            print*,yi_A,' <--> ',yi_B
           else
-             if(DEBUG) then
+             if(DEBUG.EQ.1) then
                 print *,'routine xyz_PAM ---> not properly used !!!'
                 print *,'icx = ',icx
                 print *,'icy = ',icy
-Line 1066 
 c     N.B. I convert angles from microra
+Line 926 
 c     N.B. I convert angles from microra
       $        + zi_B
       $        + dz(nplx,nldx,sensor)
+          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
-Line 1076 
 c     (xPAM,yPAM,zPAM) = measured coordi
+Line 954 
 c     (xPAM,yPAM,zPAM) = measured coordi
  c                        in PAMELA reference system
  c------------------------------------------------------------------------
-          xPAM = 0.
+          xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
-          yPAM = 0.
+          yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
-          zPAM = 0.
+          zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
+ c$$$         xPAM = 0.D0
+ c$$$         yPAM = 0.D0
+ c$$$         zPAM = 0.D0
           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
-Line 1089 
 c---------------------------------------
+Line 970 
 c---------------------------------------
           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
+          if(DEBUG.EQ.1) then
              print *,'routine xyz_PAM ---> not properly used !!!'
              print *,'icx = ',icx
              print *,'icy = ',icy
-Line 1100 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 980 
 c         print*,'A-(',xPAM_A,yPAM_A,')
        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
  continue
        end
-Line 1137 
 c$$$      common/FINALPFA/PFA
+Line 1014 
 c$$$      common/FINALPFA/PFA
  c$$$      PFAx = 'COG4'!PFA
  c$$$      PFAy = 'COG4'!PFA
        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
              print*,'xyzpam: ***WARNING*** clusters ',icx,icy
-      $           ,' does not exists (nclstr1=',nclstr1,')'
+      $           ,' do not exists (n.clusters=',nclstr1,')'
              icx = -1*icx
              icy = -1*icy
              return
-Line 1149 
 c$$$      PFAy = 'COG4'!PFA
+Line 1027 
 c$$$      PFAy = 'COG4'!PFA
        call idtoc(pfaid,PFAx)
        call idtoc(pfaid,PFAy)
- c$$$      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))
- c$$$         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
- c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
- c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
           if( (nplanes-ipx+1).ne.ip )then
              print*,'xyzpam: ***WARNING*** cluster ',icx
-Line 1184 
 c$$$     $        ,' does not belong to
+Line 1056 
 c$$$     $        ,' does not belong to
           xm(ip) = xPAM
           ym(ip) = yPAM
           zm(ip) = zPAM
-          xm_A(ip) = 0.
+          xm_A(ip) = 0.D0
-          ym_A(ip) = 0.
+          ym_A(ip) = 0.D0
-          xm_B(ip) = 0.
+          xm_B(ip) = 0.D0
-          ym_B(ip) = 0.
+          ym_B(ip) = 0.D0
  c         zv(ip) = zPAM
        elseif(icx.eq.0.and.icy.ne.0)then
           ipy=npl(VIEW(icy))
- c$$$         if((nplanes-ipy+1).ne.ip)
- c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
- c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
           if( (nplanes-ipy+1).ne.ip )then
              print*,'xyzpam: ***WARNING*** cluster ',icy
       $           ,' does not belong to plane: ',ip
-Line 1211 
 c$$$     $        ,' does not belong to
+Line 1080 
 c$$$     $        ,' does not belong to
           resx(ip)  = 1000.
           resy(ip)  = resyPAM
-          xm(ip) = -100.
+ c$$$         xm(ip) = -100.
-          ym(ip) = -100.
+ c$$$         ym(ip) = -100.
-          zm(ip) = (zPAM_A+zPAM_B)/2.
+ c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
+          xm(ip) = xPAM
+          ym(ip) = yPAM
+          zm(ip) = zPAM
           xm_A(ip) = xPAM_A
           ym_A(ip) = yPAM_A
           xm_B(ip) = xPAM_B
-Line 1224 
 c         zv(ip) = (zPAM_A+zPAM_B)/2.
+Line 1096 
 c         zv(ip) = (zPAM_A+zPAM_B)/2.
        elseif(icx.ne.0.and.icy.eq.0)then
           ipx=npl(VIEW(icx))
- c$$$         if((nplanes-ipx+1).ne.ip)
- c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
- c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
           if( (nplanes-ipx+1).ne.ip )then
              print*,'xyzpam: ***WARNING*** cluster ',icx
-Line 1242 
 c$$$     $        ,' does not belong to
+Line 1111 
 c$$$     $        ,' does not belong to
           resx(ip)  = resxPAM
           resy(ip)  = 1000.
-          xm(ip) = -100.
+ c$$$         xm(ip) = -100.
-          ym(ip) = -100.
+ c$$$         ym(ip) = -100.
-          zm(ip) = (zPAM_A+zPAM_B)/2.
+ c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
+          xm(ip) = xPAM
+          ym(ip) = yPAM
+          zm(ip) = zPAM
           xm_A(ip) = xPAM_A
           ym_A(ip) = yPAM_A
           xm_B(ip) = xPAM_B
-Line 1258 
 c         zv(ip) = (zPAM_A+zPAM_B)/2.
+Line 1130 
 c         zv(ip) = (zPAM_A+zPAM_B)/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.
-Line 1277 
 c         zv(ip) = z_mech_sensor(nplanes
+Line 1148 
 c         zv(ip) = z_mech_sensor(nplanes
        endif
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
- c         print*,'----------------------------- track coord'
    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
-Line 1307 
 c$$$         print*,'-------------------
+Line 1176 
 c$$$         print*,'-------------------
  *
  ********************************************************************************
-       real function distance_to(XPP,YPP)
+       real function distance_to(rXPP,rYPP)
        include 'common_xyzPAM.f'
-Line 1316 
 c$$$         print*,'-------------------
+Line 1185 
 c$$$         print*,'-------------------
  *     ( i.e. distance/resolution )
  *     -----------------------------------
+       real rXPP,rYPP
+       double precision XPP,YPP
        double precision distance,RE
        double precision BETA,ALFA,xmi,ymi
+       XPP=DBLE(rXPP)
+       YPP=DBLE(rYPP)
  *     ----------------------
        if (
-      +     xPAM.eq.0.and.
+ c$$$     +     xPAM.eq.0.and.
-      +     yPAM.eq.0.and.
+ c$$$     +     yPAM.eq.0.and.
-      +     zPAM.eq.0.and.
+ c$$$     +     zPAM.eq.0.and.
       +     xPAM_A.ne.0.and.
       +     yPAM_A.ne.0.and.
       +     zPAM_A.ne.0.and.
-Line 1365 
 c$$$         print*,'-------------------
+Line 1239 
 c$$$         print*,'-------------------
  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.
+ c$$$     +     xPAM.ne.0.and.
-      +     yPAM.ne.0.and.
+ c$$$     +     yPAM.ne.0.and.
-      +     zPAM.ne.0.and.
+ c$$$     +     zPAM.ne.0.and.
       +     xPAM_A.eq.0.and.
       +     yPAM_A.eq.0.and.
       +     zPAM_A.eq.0.and.
-Line 1395 
 c$$$     $        +
+Line 1266 
 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)
-Line 1445 
 c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
+Line 1308 
 c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
        data c1/1.,0.,0.,1./
        data c2/1.,-1.,-1.,1./
        data c3/1.,1.,0.,0./
-       real*8 yvvv,xvvv
+       double precision yvvv,xvvv
        double precision xi,yi,zi
        double precision xrt,yrt,zrt
        real AA,BB
-       real yvv(4),xvv(4)
+       double precision yvv(4),xvv(4)
  *     tollerance to consider the track inside the sensitive area
        real ptoll
        data ptoll/150./          !um
-       external nviewx,nviewy,acoordsi,dcoord
+       external nviewx,nviewy,dcoordsi,dcoord
        nplpt = nplPAM            !plane
        viewx = nviewx(nplpt)
-Line 1470 
 c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
+Line 1333 
 c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
  c------------------------------------------------------------------------
  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-                if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
+                xi = dcoordsi(stripx,viewx)
-      $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
+                yi = dcoordsi(stripy,viewy)
- c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
+                zi = 0.D0
- 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------------------------------------------------------------------------
-Line 1503 
 c---------------------------------------
+Line 1360 
 c---------------------------------------
                 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.
-Line 1629 
 c      include 'common_analysis.f'
+Line 1484 
 c      include 'common_analysis.f'
        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
-Line 1729 
 c      include 'level1.f'
+Line 1583 
 c      include 'level1.f'
        integer badseed,badclx,badcly
+       iflag = iflag
+       if(DEBUG.EQ.1)print*,'cl_to_couples:'
+ cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
  *     init variables
-       ncp_tot=0
        do ip=1,nplanes
           do ico=1,ncouplemax
              clx(ip,ico)=0
-Line 1743 
 c      include 'level1.f'
+Line 1601 
 c      include 'level1.f'
           ncls(ip)=0
        enddo
        do icl=1,nclstrmax_level2
-          cl_single(icl) = 1
+          cl_single(icl) = 1     !all are single
-          cl_good(icl)   = 0
+          cl_good(icl)   = 0     !all are bad
        enddo
        do iv=1,nviews
           ncl_view(iv)  = 0
-Line 1760 
 c      include 'level1.f'
+Line 1618 
 c      include 'level1.f'
           if( ncl_view(iv).gt. nclusterlimit)then
  c            mask_view(iv) = 1
              mask_view(iv) = mask_view(iv) + 2**0
-             if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
+             if(DEBUG.EQ.1)
-      $           ,nclusterlimit,' on view ', iv,' --> masked!'
+      $        print*,' * WARNING * cl_to_couple: n.clusters > '
+      $        ,nclusterlimit,' on view ', iv,' --> masked!'
           endif
        enddo
-Line 1771 
 c            mask_view(iv) = 1
+Line 1630 
 c            mask_view(iv) = 1
        do icx=1,nclstr1          !loop on cluster (X)
           if(mod(VIEW(icx),2).eq.1)goto 10
+          if(cl_used(icx).ne.0)goto 10
  *     ----------------------------------------------------
  *     jump masked views (X VIEW)
  *     ----------------------------------------------------
-Line 1778 
 c            mask_view(iv) = 1
+Line 1639 
 c            mask_view(iv) = 1
  *     ----------------------------------------------------
  *     cut on charge (X VIEW)
  *     ----------------------------------------------------
-          if(sgnl(icx).lt.dedx_x_min)then
+          if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
              cl_single(icx)=0
              goto 10
           endif
-Line 1819 
 c     endif
+Line 1680 
 c     endif
           do icy=1,nclstr1       !loop on cluster (Y)
              if(mod(VIEW(icy),2).eq.0)goto 20
+             if(cl_used(icx).ne.0)goto 20
  *     ----------------------------------------------------
  *     jump masked views (Y VIEW)
-Line 1828 
 c     endif
+Line 1691 
 c     endif
  *     ----------------------------------------------------
  *     cut on charge (Y VIEW)
  *     ----------------------------------------------------
-             if(sgnl(icy).lt.dedx_y_min)then
+             if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
                 cl_single(icy)=0
                 goto 20
              endif
-Line 1899 
 c                  cut = chcut * sch(npl
+Line 1762 
 c                  cut = chcut * sch(npl
                 endif
                 if(ncp_plane(nplx).gt.ncouplemax)then
-                   if(verbose)print*,
+                   if(verbose.eq.1)print*,
       $                 '** warning ** number of identified '//
       $                 'couples on plane ',nplx,
       $                 'exceeds vector dimention '
       $                 ,'( ',ncouplemax,' ) --> masked!'
  c                  mask_view(nviewx(nplx)) = 2
  c                  mask_view(nviewy(nply)) = 2
-                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
+                  mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
-                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
+                  mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
                    goto 10
                 endif
-Line 1927 
 c                  mask_view(nviewy(nply
+Line 1790 
 c                  mask_view(nviewy(nply
      continue
        enddo                     !end loop on clusters(X)
        do icl=1,nclstr1
           if(cl_single(icl).eq.1)then
              ip=npl(VIEW(icl))
-Line 1935 
 c                  mask_view(nviewy(nply
+Line 1797 
 c                  mask_view(nviewy(nply
              cls(ip,ncls(ip))=icl
           endif
        enddo
+ c 80   continue
+       continue
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
           print*,'clusters  ',nclstr1
           print*,'good    ',(cl_good(i),i=1,nclstr1)
-          print*,'singles ',(cl_single(i),i=1,nclstr1)
+          print*,'used    ',(cl_used(i),i=1,nclstr1)
+          print*,'singlets',(cl_single(i),i=1,nclstr1)
           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
        endif
+       if(.not.RECOVER_SINGLETS)goto 81
+ *     ////////////////////////////////////////////////
+ *     PATCH to recover events with less than 3 couples
+ *     ////////////////////////////////////////////////
+ *     loop over singlet and create "fake" couples
+ *     (with clx=0 or cly=0)
+ *
+       if(DEBUG.EQ.1)
+      $     print*,'>>> Recover singlets '
+      $     ,'(creates fake couples) <<<'
+       do icl=1,nclstr1
+          if(
+      $        cl_single(icl).eq.1.and.
+      $        cl_used(icl).eq.0.and.
+      $        .true.)then
+ *     ----------------------------------------------------
+ *     jump masked views
+ *     ----------------------------------------------------
+             if( mask_view(VIEW(icl)).ne.0 ) goto 21
+             if(mod(VIEW(icl),2).eq.0)then !=== X views
+ *     ----------------------------------------------------
+ *     cut on charge (X VIEW)
+ *     ----------------------------------------------------
+                if(sgnl(icl).lt.dedx_x_min) goto 21
+                nplx=npl(VIEW(icl))
+ *     ------------------> (FAKE) COUPLE <-----------------
+                ncp_plane(nplx) = ncp_plane(nplx) + 1
+                clx(nplx,ncp_plane(nplx))=icl
+                cly(nplx,ncp_plane(nplx))=0
+ c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
+ *     ----------------------------------------------------
+             else                !=== Y views
+ *     ----------------------------------------------------
+ *     cut on charge (Y VIEW)
+ *     ----------------------------------------------------
+                if(sgnl(icl).lt.dedx_y_min) goto 21
+                nply=npl(VIEW(icl))
+ *     ------------------> (FAKE) COUPLE <-----------------
+                ncp_plane(nply) = ncp_plane(nply) + 1
+                clx(nply,ncp_plane(nply))=0
+                cly(nply,ncp_plane(nply))=icl
+ c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
+ *     ----------------------------------------------------
+             endif
+          endif                  !end "single" condition
+     continue
+       enddo                     !end loop over clusters
+       if(DEBUG.EQ.1)
+      $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
+  continue
-       do ip=1,6
+       ncp_tot=0
+       do ip=1,NPLANES
           ncp_tot = ncp_tot + ncp_plane(ip)
        enddo
+       if(DEBUG.EQ.1)
+      $     print*,'n.couple tot:      ',ncp_tot
        return
        end
-Line 2000 
 c      double precision xm3,ym3,zm3
+Line 1930 
 c      double precision xm3,ym3,zm3
        real xc,zc,radius
  *     -----------------------------
+       if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
  *     --------------------------------------------
  *     put a limit to the maximum number of couples
-Line 2008 
 c      double precision xm3,ym3,zm3
+Line 1939 
 c      double precision xm3,ym3,zm3
  *     --------------------------------------------
        do ip=1,nplanes
           if(ncp_plane(ip).gt.ncouplelimit)then
- c            mask_view(nviewx(ip)) = 8
- c            mask_view(nviewy(ip)) = 8
              mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
              mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
           endif
-Line 2020 
 c            mask_view(nviewy(ip)) = 8
+Line 1949 
 c            mask_view(nviewy(ip)) = 8
        ntrpt=0                   !number of triplets
        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
+ c$$$         print*,'(1) ip ',ip1
+ c$$$     $        ,mask_view(nviewx(ip1))
+ c$$$     $        ,mask_view(nviewy(ip1))
           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$$$               print*,'(1) ip ',ip1,' icp ',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
+ c$$$                  print*,'(2) ip ',ip2
+ c$$$     $                 ,mask_view(nviewx(ip2))
+ c$$$     $                 ,mask_view(nviewy(ip2))
                    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$$$                        print*,'(2) ip ',ip2,' icp ',icp2
  c                        call xyz_PAM
  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
  c                        call xyz_PAM
-Line 2051 
 c     $                       (icx2,icy2
+Line 1990 
 c     $                       (icx2,icy2
                          xm2=xPAM
                          ym2=yPAM
                          zm2=zPAM
+ *                       ---------------------------------------------------
+ *                       both couples must have a y-cluster
+ *                       (condition necessary when in RECOVER_SINGLETS mode)
+ *                       ---------------------------------------------------
+                         if(icy1.eq.0.or.icy2.eq.0)goto 111
+                         if(cl_used(icy1).ne.0)goto 111
+                         if(cl_used(icy2).ne.0)goto 111
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  *     track parameters on Y VIEW
  *     (2 couples needed)
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                          if(ndblt.eq.ndblt_max)then
-                            if(verbose)print*,
+                            if(verbose.eq.1)print*,
       $                          '** warning ** number of identified '//
       $                          'doublets exceeds vector dimention '
       $                          ,'( ',ndblt_max,' )'
- c                           good2=.false.
+ c     good2=.false.
- c                           goto 880 !fill ntp and go to next event
+ c     goto 880 !fill ntp and go to next event
                             do iv=1,12
- c                              mask_view(iv) = 3
+ c     mask_view(iv) = 3
                                mask_view(iv) = mask_view(iv)+ 2**2
                             enddo
                             iflag=1
                             return
                          endif
+ ccc                        print*,'<doublet> ',icp1,icp2
                          ndblt = ndblt + 1
  *     store doublet info
                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
-Line 2078 
 c                              mask_view
+Line 2031 
 c                              mask_view
                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
  *     y0 (cm)
                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
  ****  -----------------------------------------------****
  ****  reject non phisical couples                    ****
  ****  -----------------------------------------------****
+                         if(SECOND_SEARCH)goto 111
                          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
+                   continue
- 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(icx1.ne.0)then
+                            if(cl_used(icx1).ne.0)goto 31
+                         endif
+                         if(icx2.ne.0)then
+                            if(cl_used(icx2).ne.0)goto 31
+                         endif
                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
+ c$$$                           print*,'(3) ip ',ip3
+ c$$$     $                          ,mask_view(nviewx(ip3))
+ c$$$     $                          ,mask_view(nviewy(ip3))
                             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
-Line 2110 
 c$$$
+Line 2069 
 c$$$
                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
                                   icx3=clx(ip3,icp3)
                                   icy3=cly(ip3,icp3)
+ c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
+ *     ---------------------------------------------------
+ *     all three couples must have a x-cluster
+ *     (condition necessary when in RECOVER_SINGLETS mode)
+ *     ---------------------------------------------------
+                                  if(
+      $                                icx1.eq.0.or.
+      $                                icx2.eq.0.or.
+      $                                icx3.eq.0.or.
+      $                                .false.)goto 29
+                                  if(cl_used(icx1).ne.0)goto 29
+                                  if(cl_used(icx2).ne.0)goto 29
+                                  if(cl_used(icx3).ne.0)goto 29
  c                                 call xyz_PAM
  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
  c                                 call xyz_PAM
-Line 2120 
 c     $                               (i
+Line 2096 
 c     $                               (i
                                   xm3=xPAM
                                   ym3=yPAM
                                   zm3=zPAM
  *     find the circle passing through the three points
+                                  iflag_t = DEBUG
                                   call tricircle(3,xp,zp,angp,resp,chi
-      $                                ,xc,zc,radius,iflag)
+      $                                ,xc,zc,radius,iflag_t)
- c     print*,xc,zc,radius
  *     the circle must intersect the reference plane
-                                  if(
+ cc                                 if(iflag.ne.0)goto 29
- c     $                                 (xc.le.-1.*xclimit.or.
+                                  if(iflag_t.ne.0)then
- c     $                                 xc.ge.xclimit).and.
+ *     if tricircle fails, evaluate a straight line
-      $                                radius**2.ge.(ZINI-zc)**2.and.
+                                     if(DEBUG.eq.1)
-      $                                iflag.eq.0.and.
+      $                                   print*,'TRICIRCLE failure'
-      $                                .true.)then
+      $                                   ,' >>> straight line'
+                                     radius=0.
+                                     xc=0.
+                                     yc=0.
+                                     SZZ=0.
+                                     SZX=0.
+                                     SSX=0.
+                                     SZ=0.
+                                     S1=0.
+                                     X0=0.
+                                     Ax=0.
+                                     BX=0.
+                                     DO I=1,3
+                                        XX = XP(I)
+                                        SZZ=SZZ+ZP(I)*ZP(I)
+                                        SZX=SZX+ZP(I)*XX
+                                        SSX=SSX+XX
+                                        SZ=SZ+ZP(I)
+                                        S1=S1+1.
+                                     ENDDO
+                                     DET=SZZ*S1-SZ*SZ
+                                     AX=(SZX*S1-SZ*SSX)/DET
+                                     BX=(SZZ*SSX-SZX*SZ)/DET
+                                     X0  = AX*ZINI+BX
+                                  endif
+                                  if(  .not.SECOND_SEARCH.and.
+      $                                radius**2.lt.(ZINI-zc)**2)goto 29
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  *     track parameters on X VIEW
  *     (3 couples needed)
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                                   if(ntrpt.eq.ntrpt_max)then
-                                     if(verbose)print*,
+                                     if(verbose.eq.1)print*,
-      $                     '** warning ** number of identified '//
+      $                                   '** warning **'//
-      $                     'triplets exceeds vector dimention '
+      $                                   ' number of identified '//
-      $                    ,'( ',ntrpt_max,' )'
+      $                                   'triplets exceeds'//
- c                                    good2=.false.
+      $                                   ' vector dimention '
- c                                    goto 880 !fill ntp and go to next event
+      $                                   ,'( ',ntrpt_max,' )'
+ c     good2=.false.
+ c     goto 880 !fill ntp and go to next event
                                      do iv=1,nviews
- c                                       mask_view(iv) = 4
+ c     mask_view(iv) = 4
-                                        mask_view(iv)=mask_view(iv)+ 2**3
+                                        mask_view(iv) =
+      $                                      mask_view(iv)+ 2**3
                                      enddo
                                      iflag=1
                                      return
                                   endif
+ ccc                                 print*,'<triplet> ',icp1,icp2,icp3
                                   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
+                                  if(radius.ne.0.and.xc.lt.0)then
  *************POSITIVE DEFLECTION
-               alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
+              alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
-               alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
+              alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
-               alfaxz3(ntrpt) = 1/radius
+              alfaxz3(ntrpt) = 1/radius
-                                  else
+                                 else if(radius.ne.0.and.xc.ge.0)then
  *************NEGATIVE DEFLECTION
-               alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
+              alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
-               alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
+              alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
-               alfaxz3(ntrpt) = -1/radius
+              alfaxz3(ntrpt) = -1/radius
-                                  endif
+                                 else if(radius.eq.0)then
+ *************straight fit
+              alfaxz1(ntrpt) = X0
+              alfaxz2(ntrpt) = AX
+              alfaxz3(ntrpt) = 0.
+                                 endif
+ c$$$                                print*,'alfaxz1 ', alfaxz1(ntrpt)
+ c$$$                                print*,'alfaxz2 ', alfaxz2(ntrpt)
+ c$$$                                print*,'alfaxz3 ', alfaxz3(ntrpt)
  ****  -----------------------------------------------****
  ****  reject non phisical triplets                   ****
  ****  -----------------------------------------------****
-                                  if(
+                                 if(SECOND_SEARCH)goto 29
-      $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max
+                                 if(
-      $                                .or.
+      $                               abs(alfaxz2(ntrpt)).gt.
-      $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max
+      $                               alfxz2_max
-      $                                )ntrpt = ntrpt-1
+      $                               .or.
+      $                               abs(alfaxz1(ntrpt)).gt.
+      $                               alfxz1_max
- c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)
+      $                               )ntrpt = ntrpt-1
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  *     track parameters on X VIEW - end
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-                               endif
+                          continue
                             enddo !end loop on COPPIA 3
                          enddo   !end loop on sensors - COPPIA 3
                     continue
                       enddo      !end loop on planes  - COPPIA 3
-                 continue
+                 continue
-               enddo         !end loop on COPPIA 2
+ c 1                enddo         !end loop on COPPIA 2
+                  enddo         !end loop on COPPIA 2
                 enddo            !end loop on sensors - COPPIA 2
            continue
              enddo               !end loop on planes  - COPPIA 2
+ c 11         continue
+           continue
           enddo                  !end loop on COPPIA1
        enddo                     !end loop on sensors - COPPIA 1
   continue
        enddo                     !end loop on planes  - COPPIA 1
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
           print*,'--- doublets ',ndblt
           print*,'--- triplets ',ntrpt
        endif
-Line 2246 
 c      include 'momanhough_init.f'
+Line 2273 
 c      include 'momanhough_init.f'
        integer cp_useds1(ncouplemaxtot) ! sensor 1
        integer cp_useds2(ncouplemaxtot) ! sensor 2
+       if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  *     classification of DOUBLETS
-Line 2264 
 c      include 'momanhough_init.f'
+Line 2292 
 c      include 'momanhough_init.f'
        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
-Line 2302 
 ccccc if(db_used(idbref).eq.1)goto 1188
+Line 2327 
 ccccc if(db_used(idbref).eq.1)goto 1188
  *     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
-Line 2330 
 c     print*,idb1,idb2,distance,' cloud
+Line 2347 
 c     print*,idb1,idb2,distance,' cloud
                    temp1 = temp1 + alfayz1(idb2)
                    temp2 = temp2 + alfayz2(idb2)
- c     print*,'*   idbref,idb2 ',idbref,idb2
                 endif
          continue
              enddo               !end loop (2) on DOUBLETS
-      continue
+ c 1188       continue
+             continue
           enddo                  !end loop on... bo?
           nptloop=npv
-Line 2353 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2370 
 c     print*,'*   idbref,idb2 ',idbref,i
           enddo
           ncpused=0
           do icp=1,ncp_tot
-             if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
+             if(
+      $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
+      $           .true.)then
                 ncpused=ncpused+1
                 ip=ip_cp(icp)
                 hit_plane(ip)=1
-Line 2363 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2382 
 c     print*,'*   idbref,idb2 ',idbref,i
           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*,
+             if(verbose.eq.1)print*,
       $           '** warning ** number of identified '//
       $           'YZ clouds exceeds vector dimention '
       $           ,'( ',ncloyz_max,' )'
-Line 2397 
 c               mask_view(iv) = 5
+Line 2414 
 c               mask_view(iv) = 5
  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
+          if(DEBUG.EQ.1)then
-             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
-             print*,'- alfayz1 ',alfayz1_av(nclouds_yz)
+             print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
-             print*,'- alfayz2 ',alfayz2_av(nclouds_yz)
+             print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
-             print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
+             print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
-             print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
+             print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
-             print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
+             print*,'cpcloud_yz '
- c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
+      $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
- c$$$     $           ,ptcloud_yz(nclouds_yz)
+             print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
- c$$$            print*,'nt-uple: db_cloud(...) = '
- c$$$     $           ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)
           endif
  *     >>> NEW CLOUD <<<
  *     ~~~~~~~~~~~~~~~~~
-Line 2425 
 c$$$     $           ,(db_cloud(iii),iii
+Line 2438 
 c$$$     $           ,(db_cloud(iii),iii
          goto 90
        endif
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
-          print*,'---------------------- '
           print*,'Y-Z total clouds ',nclouds_yz
-          print*,' '
        endif
-Line 2474 
 c      include 'momanhough_init.f'
+Line 2485 
 c      include 'momanhough_init.f'
        integer cp_useds1(ncouplemaxtot) ! sensor 1
        integer cp_useds2(ncouplemaxtot) ! sensor 2
+       if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  *     classification of TRIPLETS
  *     according to distance in parameter space
-Line 2489 
 c      include 'momanhough_init.f'
+Line 2502 
 c      include 'momanhough_init.f'
   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
-Line 2524 
 c         tr_temp(1)=itr1
+Line 2535 
 c         tr_temp(1)=itr1
              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
+ *     ------------------------------------------------------------------------
+ *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
+ *     ------------------------------------------------------------------------
+ *     (added in august 2007)
+                istrimage=0
+                if(
+      $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
+      $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
+      $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
+      $              .true.)istrimage=1
+                if(distance.lt.cutdistxz.or.istrimage.eq.1)then
                    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
-Line 2551 
 c     print*,idb1,idb2,distance,' cloud
+Line 2575 
 c     print*,idb1,idb2,distance,' cloud
                    temp1 = temp1 + alfaxz1(itr2)
                    temp2 = temp2 + alfaxz2(itr2)
                    temp3 = temp3 + alfaxz3(itr2)
- c     print*,'*   itrref,itr2 ',itrref,itr2,distance
                 endif
          continue
              enddo               !end loop (2) on TRIPLETS
-      continue
+ c11888       continue
+             continue
           enddo                  !end loop on... bo?
           nptloop=npv
-Line 2572 
 c     print*,'*   itrref,itr2 ',itrref,i
+Line 2596 
 c     print*,'*   itrref,itr2 ',itrref,i
  *     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
+             if(
+      $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
+      $           .true.)then
                 ncpused=ncpused+1
                 ip=ip_cp(icp)
                 hit_plane(ip)=1
-Line 2588 
 c     print*,'check cp_used'
+Line 2613 
 c     print*,'check cp_used'
           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*,
+             if(verbose.eq.1)print*,
       $           '** warning ** number of identified '//
       $           'XZ clouds exceeds vector dimention '
       $           ,'( ',ncloxz_max,' )'
-Line 2622 
 c               mask_view(iv) = 6
+Line 2645 
 c               mask_view(iv) = 6
           enddo
           npt_tot=npt_tot+npt
-          if(DEBUG)then
+          if(DEBUG.EQ.1)then
-             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
+             print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
-             print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
+             print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
-             print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
+             print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
-             print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
+             print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
-             print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
+             print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
-             print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
+             print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
-             print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
+             print*,'cpcloud_xz '
+      $           ,(cpcloud_xz(nclouds_xz,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 <<<
  *     ~~~~~~~~~~~~~~~~~
-Line 2647 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2667 
 c$$$     $           ,(tr_cloud(iii),iii
           goto 91
         endif
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
-          print*,'---------------------- '
           print*,'X-Z total clouds ',nclouds_xz
-          print*,' '
        endif
-Line 2699 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2717 
 c$$$     $           ,(tr_cloud(iii),iii
  *     -----------------------------------------------------------
  *     variables for track fitting
        double precision AL_INI(5)
- c      double precision tath
  *     -----------------------------------------------------------
- c      real fitz(nplanes)        !z coordinates of the planes in cm
+       if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
        ntracks=0                 !counter of track candidates
-       do iyz=1,nclouds_yz       !loop on YZ couds
+       do iyz=1,nclouds_yz       !loop on YZ clouds
-          do ixz=1,nclouds_xz    !loop on XZ couds
+          do ixz=1,nclouds_xz    !loop on XZ clouds
  *     --------------------------------------------------
  *     check of consistency of the clouds
-Line 2717 
 c      real fitz(nplanes)        !z coor
+Line 2734 
 c      real fitz(nplanes)        !z coor
  *     of the two clouds
  *     --------------------------------------------------
              do ip=1,nplanes
-                hit_plane(ip)=0
+                hit_plane(ip)=0  !n.matching couples (REAL couples, not SINGLETS)
-                ncp_match(ip)=0
+                ncp_match(ip)=0  !n.matching couples per plane
                 do icpp=1,ncouplemax
                    cp_match(ip,icpp)=0 !init couple list
                 enddo
              enddo
-             ncp_ok=0
+             ncp_ok=0            !count n.matching-couples
-             do icp=1,ncp_tot    !loop on couples
+             ncpx_ok=0           !count n.matching-couples with x cluster
- *     get info on
+             ncpy_ok=0           !count n.matching-couples with y cluster
-                cpintersec(icp)=min(
-      $              cpcloud_yz(iyz,icp),
-      $              cpcloud_xz(ixz,icp))
+             do icp=1,ncp_tot    !loop over couples
-                if(
-      $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
+                if(.not.RECOVER_SINGLETS)then
-      $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
+ *     ------------------------------------------------------
-      $              .false.)cpintersec(icp)=0
+ *     if NOT in RECOVER_SINGLETS mode, take the intersection
+ *     between xz yz clouds
+ *     ------------------------------------------------------
+                   cpintersec(icp)=min(
+      $                 cpcloud_yz(iyz,icp),
+      $                 cpcloud_xz(ixz,icp))
+ *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
+ *     ------------------------------------------------------
+ *     discard the couple if the sensor is in conflict
+ *     ------------------------------------------------------
+                   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
+                else
+ *     ------------------------------------------------------
+ *     if RECOVER_SINGLETS take the union
+ *     (otherwise the fake couples formed by singlets would be
+ *     discarded)
+ *     ------------------------------------------------------
+                   cpintersec(icp)=max(
+      $                 cpcloud_yz(iyz,icp),
+      $                 cpcloud_xz(ixz,icp))
+ c$$$                  if(cpcloud_yz(iyz,icp).gt.0)
+ c$$$     $                 cpintersec(icp)=cpcloud_yz(iyz,icp)
+ *     cpintersec is >0 if either yz or xz cloud contains the couple icp
+                endif
+ c$$$               print*,icp,ip_cp(icp),' -- ',cpintersec(icp)
                 if(cpintersec(icp).ne.0)then
-                   ncp_ok=ncp_ok+1
                    ip=ip_cp(icp)
                    hit_plane(ip)=1
+ c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0)
+ c$$$     $                 ncp_ok=ncp_ok+1
+ c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0)
+ c$$$     $                 ncpx_ok=ncpx_ok+1
+ c$$$                  if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0)
+ c$$$     $                 ncpy_ok=ncpy_ok+1
+                   if(  cpcloud_yz(iyz,icp).gt.0.and.
+      $                 cpcloud_xz(ixz,icp).gt.0)
+      $                 ncp_ok=ncp_ok+1
+                   if(  cpcloud_yz(iyz,icp).gt.0.and.
+      $                 cpcloud_xz(ixz,icp).eq.0)
+      $                 ncpy_ok=ncpy_ok+1
+                   if(  cpcloud_yz(iyz,icp).eq.0.and.
+      $                 cpcloud_xz(ixz,icp).gt.0)
+      $                 ncpx_ok=ncpx_ok+1
                    if(cpintersec(icp).eq.1)then
  *     1) only the couple image in sensor 1 matches
                       id=-icp
-Line 2764 
 c      real fitz(nplanes)        !z coor
+Line 2826 
 c      real fitz(nplanes)        !z coor
              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
+             if(nplused.lt.3)goto 888 !next combination
+ ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
+ ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
+ *     -----------------------------------------------------------
+ *     if in RECOVER_SINGLET mode, the two clouds must have
+ *     at least ONE intersecting real couple
+ *     -----------------------------------------------------------
+             if(ncp_ok.lt.1)goto 888 !next combination
+             if(DEBUG.EQ.1)then
+                print*,'////////////////////////////'
+                print*,'Cloud combination (Y,X): ',iyz,ixz
+                print*,' db ',ptcloud_yz(iyz)
+                print*,' tr ',ptcloud_xz(ixz)
+                print*,'  -----> # matching couples ',ncp_ok
+                print*,'  -----> # fake couples (X)',ncpx_ok
+                print*,'  -----> # fake couples (Y)',ncpy_ok
+                do icp=1,ncp_tot
+                   print*,'cp ',icp,' >'
+      $                 ,' x ',cpcloud_xz(ixz,icp)
+      $                 ,' y ',cpcloud_yz(iyz,icp)
+      $                 ,' ==> ',cpintersec(icp)
+                enddo
+             endif
+             if(DEBUG.EQ.1)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))
-Line 2853 
 c$$$            if(AL_INI(5).gt.defmax)g
+Line 2885 
 c$$$            if(AL_INI(5).gt.defmax)g
                                hit_plane(6)=icp6
                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
+                               if(DEBUG.eq.1)
+      $                             print*,'combination: '
+      $                             ,cp_match(6,icp1)
+      $                             ,cp_match(5,icp2)
+      $                             ,cp_match(4,icp3)
+      $                             ,cp_match(3,icp4)
+      $                             ,cp_match(2,icp5)
+      $                             ,cp_match(1,icp6)
+ *                             ---------------------------------------
+ *                             check if this group of couples has been
+ *                             already fitted
+ *                             ---------------------------------------
+                               do ica=1,ntracks
+                                  isthesame=1
+                                  do ip=1,NPLANES
+                                     if(hit_plane(ip).ne.0)then
+                                        if(  CP_STORE(nplanes-ip+1,ica)
+      $                                      .ne.
+      $                                      cp_match(ip,hit_plane(ip)) )
+      $                                      isthesame=0
+                                     else
+                                        if(  CP_STORE(nplanes-ip+1,ica)
+      $                                      .ne.
+      $                                      0 )
+      $                                      isthesame=0
+                                     endif
+                                  enddo
+                                  if(isthesame.eq.1)then
+                                     if(DEBUG.eq.1)
+      $                                   print*,'(already fitted)'
+                                     goto 666 !jump to next combination
+                                  endif
+                               enddo
                                call track_init !init TRACK common
-                               do ip=1,nplanes !loop on planes
+                               do ip=1,nplanes !loop on planes (bottom to top)
                                   if(hit_plane(ip).ne.0)then
                                      id=cp_match(ip,hit_plane(ip))
                                      is=is_cp(id)
-Line 2876 
 c     $
+Line 2943 
 c     $
       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
  *                                   *************************
  *                                   -----------------------------
-                                     xgood(nplanes-ip+1)=1.
+                                     if(icx.gt.0.and.icy.gt.0)then
-                                     ygood(nplanes-ip+1)=1.
+                                        xgood(nplanes-ip+1)=1.
-                                     xm(nplanes-ip+1)=xPAM
+                                        ygood(nplanes-ip+1)=1.
-                                     ym(nplanes-ip+1)=yPAM
+                                        xm(nplanes-ip+1)=xPAM
-                                     zm(nplanes-ip+1)=zPAM
+                                        ym(nplanes-ip+1)=yPAM
-                                     resx(nplanes-ip+1)=resxPAM
+                                        zm(nplanes-ip+1)=zPAM
-                                     resy(nplanes-ip+1)=resyPAM
+                                        resx(nplanes-ip+1)=resxPAM
+                                        resy(nplanes-ip+1)=resyPAM
+                                        if(DEBUG.EQ.1)print*,'(X,Y)'
+      $                                      ,nplanes-ip+1,xPAM,yPAM
+                                     else
+                                        xm_A(nplanes-ip+1) = xPAM_A
+                                        ym_A(nplanes-ip+1) = yPAM_A
+                                        xm_B(nplanes-ip+1) = xPAM_B
+                                        ym_B(nplanes-ip+1) = yPAM_B
+                                        zm(nplanes-ip+1)
+      $                                      = (zPAM_A+zPAM_B)/2.
+                                        resx(nplanes-ip+1) = resxPAM
+                                        resy(nplanes-ip+1) = resyPAM
+                                        if(icx.eq.0.and.icy.gt.0)then
+                                           xgood(nplanes-ip+1)=0.
+                                           ygood(nplanes-ip+1)=1.
+                                           resx(nplanes-ip+1) = 1000.
+                                           if(DEBUG.EQ.1)print*,'(  Y)'
+      $                                         ,nplanes-ip+1,xPAM,yPAM
+                                        elseif(icx.gt.0.and.icy.eq.0)then
+                                           xgood(nplanes-ip+1)=1.
+                                           ygood(nplanes-ip+1)=0.
+                                           if(DEBUG.EQ.1)print*,'(X  )'
+      $                                         ,nplanes-ip+1,xPAM,yPAM
+                                           resy(nplanes-ip+1) = 1000.
+                                        else
+                                           print*,'both icx=0 and icy=0'
+      $                                         ,' ==> IMPOSSIBLE!!'
+                                        endif
+                                     endif
  *                                   -----------------------------
                                   endif
                                enddo !end loop on planes
-Line 2900 
 c$$$                              enddo
+Line 2996 
 c$$$                              enddo
                                ifail=0 !error flag in chi^2 computation
                                jstep=0 !number of  minimization steps
                                iprint=0
- c                              if(DEBUG)iprint=1
+ c                              if(DEBUG.EQ.1)iprint=1
-                               if(DEBUG)iprint=2
+                               if(DEBUG.EQ.1)iprint=2
                                call mini2(jstep,ifail,iprint)
                                if(ifail.ne.0) then
-                                  if(DEBUG)then
+                                  if(DEBUG.EQ.1)then
                                      print *,
       $                              '*** MINIMIZATION FAILURE *** '
       $                              //'(clouds_to_ctrack)'
-Line 2923 
 c                                 chi2=-
+Line 3019 
 c                                 chi2=-
  *     **********************************************************
                                if(chi2.le.0.)goto 666
+                               if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
+                               if(chi2.ne.chi2)goto 666  !OPTIMIZATION
  *     --------------------------
  *     STORE candidate TRACK INFO
  *     --------------------------
                                if(ntracks.eq.NTRACKSMAX)then
-                                  if(verbose)print*,
+                                  if(verbose.eq.1)print*,
       $                 '** warning ** number of candidate tracks '//
       $                 ' exceeds vector dimension '
       $                ,'( ',NTRACKSMAX,' )'
-Line 2945 
 c                                    mas
+Line 3043 
 c                                    mas
                                ntracks = ntracks + 1
-                               do ip=1,nplanes
+                               do ip=1,nplanes !top to bottom
                                   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))
-Line 2962 
 c                                    mas
+Line 3060 
 c                                    mas
                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
+ *                                NB! hit_plane is defined from bottom to top
                                   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(
+                                     icl=
       $                                   clx(ip,icp_cp(
       $                                   cp_match(ip,hit_plane(ip)
-      $                                   ))));
+      $                                   )));
+                                     if(icl.eq.0)
+      $                                   icl=
+      $                                   cly(ip,icp_cp(
+      $                                   cp_match(ip,hit_plane(ip)
+      $                                   )));
+                                     LADDER_STORE(nplanes-ip+1,ntracks)
+      $                                   = LADDER(icl);
                                   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
+                                  BX_STORE(ip,ntracks)=0!I dont need it now
-                                  BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
+                                  BY_STORE(ip,ntracks)=0!I dont need it now
-                                  CLS_STORE(nplanes-ip+1,ntracks)=0
+                                  CLS_STORE(ip,ntracks)=0
                                   do i=1,5
                                      AL_STORE(i,ntracks)=sngl(AL(i))
                                   enddo
-Line 3005 
 c                                    mas
+Line 3112 
 c                                    mas
        if(ntracks.eq.0)then
           iflag=1
-          return
+ cc         return
        endif
- c$$$      if(DEBUG)then
+       if(DEBUG.EQ.1)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
-Line 3069 
 c$$$      endif
+Line 3167 
 c$$$      endif
        real xyzp(3),bxyz(3)
        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
+       if(DEBUG.EQ.1)print*,'refine_track:'
  *     =================================================
  *     new estimate of positions using ETA algorithm
  *                          and
-Line 3077 
 c$$$      endif
+Line 3176 
 c$$$      endif
        call track_init
        do ip=1,nplanes           !loop on planes
+          if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
           xP=XV_STORE(nplanes-ip+1,ibest)
           yP=YV_STORE(nplanes-ip+1,ibest)
           zP=ZV_STORE(nplanes-ip+1,ibest)
-Line 3093 
 c$$$         bxyz(3)=0
+Line 3194 
 c$$$         bxyz(3)=0
  *     using improved PFAs
  *     -------------------------------------------------
  *     |||||||||||||||||||||||||||||||||||||||||||||||||
-          if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
+ c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
+ c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
+ c$$$
+ c$$$            id=CP_STORE(nplanes-ip+1,ibest)
+ c$$$
+ c$$$            is=is_cp(id)
+ c$$$            icp=icp_cp(id)
+ c$$$            if(ip_cp(id).ne.ip)
+ c$$$     $           print*,'OKKIO!!'
+ c$$$     $           ,'id ',id,is,icp
+ c$$$     $           ,ip_cp(id),ip
+ c$$$            icx=clx(ip,icp)
+ c$$$            icy=cly(ip,icp)
+ c$$$c            call xyz_PAM(icx,icy,is,
+ c$$$c     $           PFA,PFA,
+ c$$$c     $           AXV_STORE(nplanes-ip+1,ibest),
+ c$$$c     $           AYV_STORE(nplanes-ip+1,ibest))
+ c$$$            call xyz_PAM(icx,icy,is,
+ c$$$     $           PFA,PFA,
+ c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
+ c$$$     $           AYV_STORE(nplanes-ip+1,ibest),
+ c$$$     $           bxyz(1),
+ c$$$     $           bxyz(2)
+ c$$$     $           )
+ c$$$
+ c$$$            xm(nplanes-ip+1) = xPAM
+ c$$$            ym(nplanes-ip+1) = yPAM
+ c$$$            zm(nplanes-ip+1) = zPAM
+ c$$$            xgood(nplanes-ip+1) = 1
+ c$$$            ygood(nplanes-ip+1) = 1
+ c$$$            resx(nplanes-ip+1) = resxPAM
+ c$$$            resy(nplanes-ip+1) = resyPAM
+ c$$$
+ c$$$            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
+ c$$$            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+          if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
              id=CP_STORE(nplanes-ip+1,ibest)
-Line 3118 
 c     $           AYV_STORE(nplanes-ip+1
+Line 3254 
 c     $           AYV_STORE(nplanes-ip+1
       $           bxyz(2)
       $           )
-             xm(nplanes-ip+1) = xPAM
+             if(icx.gt.0.and.icy.gt.0)then
-             ym(nplanes-ip+1) = yPAM
+                xm(nplanes-ip+1) = xPAM
-             zm(nplanes-ip+1) = zPAM
+                ym(nplanes-ip+1) = yPAM
-             xgood(nplanes-ip+1) = 1
+                zm(nplanes-ip+1) = zPAM
-             ygood(nplanes-ip+1) = 1
+                xm_A(nplanes-ip+1) = 0.
-             resx(nplanes-ip+1) = resxPAM
+                ym_A(nplanes-ip+1) = 0.
-             resy(nplanes-ip+1) = resyPAM
+                xm_B(nplanes-ip+1) = 0.
+                ym_B(nplanes-ip+1) = 0.
-             dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
+                xgood(nplanes-ip+1) = 1
-             dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+                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))
+             else
+                xm(nplanes-ip+1) = 0.
+                ym(nplanes-ip+1) = 0.
+                zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
+                xm_A(nplanes-ip+1) = xPAM_A
+                ym_A(nplanes-ip+1) = yPAM_A
+                xm_B(nplanes-ip+1) = xPAM_B
+                ym_B(nplanes-ip+1) = yPAM_B
+                xgood(nplanes-ip+1) = 0
+                ygood(nplanes-ip+1) = 0
+                resx(nplanes-ip+1) = 1000.!resxPAM
+                resy(nplanes-ip+1) = 1000.!resyPAM
+                dedxtrk_x(nplanes-ip+1)= 0
+                dedxtrk_y(nplanes-ip+1)= 0
+                if(icx.gt.0)then
+                   xgood(nplanes-ip+1) = 1
+                   resx(nplanes-ip+1) = resxPAM
+                   dedxtrk_x(nplanes-ip+1)=
+      $                 sgnl(icx)/mip(VIEW(icx),LADDER(icx))
+                elseif(icy.gt.0)then
+                   ygood(nplanes-ip+1) = 1
+                   resy(nplanes-ip+1) = resyPAM
+                   dedxtrk_y(nplanes-ip+1)=
+      $                 sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+                endif
+             endif
  *     |||||||||||||||||||||||||||||||||||||||||||||||||
  *     -------------------------------------------------
-Line 3139 
 c     $           AYV_STORE(nplanes-ip+1
+Line 3307 
 c     $           AYV_STORE(nplanes-ip+1
              xgood(nplanes-ip+1)=0
              ygood(nplanes-ip+1)=0
+             CP_STORE(nplanes-ip+1,ibest)=0 !re-init
+             CLS_STORE(nplanes-ip+1,ibest)=0
  *     --------------------------------------------------------------
  *     determine which ladder and sensor are intersected by the track
-Line 3150 
 c     $           AYV_STORE(nplanes-ip+1
+Line 3322 
 c     $           AYV_STORE(nplanes-ip+1
              LADDER_STORE(nplanes-ip+1,IBEST)=nldt
  *     --------------------------------------------------------------
-             if(DEBUG)then
+             if(DEBUG.EQ.1)then
                 print*,
       $              '------ Plane ',ip,' intersected on LADDER ',nldt
       $              ,' SENSOR ',ist
-Line 3161 
 c     $           AYV_STORE(nplanes-ip+1
+Line 3333 
 c     $           AYV_STORE(nplanes-ip+1
  *     ===========================================
  *     STEP 1 >>>>>>>  try to include a new couple
  *     ===========================================
- c            if(DEBUG)print*,'>>>> try to include a new couple'
+             distmin=100000000.
-             distmin=1000000.
              xmm = 0.
              ymm = 0.
              zmm = 0.
-Line 3175 
 c            if(DEBUG)print*,'>>>> try t
+Line 3346 
 c            if(DEBUG)print*,'>>>> try t
              do icp=1,ncp_plane(ip) !loop on couples on plane icp
                 icx=clx(ip,icp)
                 icy=cly(ip,icp)
+                if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
                 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(icx).ne.0.or. !or the X cluster is already used
-      $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
+      $              cl_used(icy).ne.0.or. !or the Y cluster is already used
       $              .false.)goto 1188 !then jump to next couple.
  *
                 call xyz_PAM(icx,icy,ist,
-Line 3193 
 c     $              cl_used(icy).eq.1.o
+Line 3365 
 c     $              cl_used(icy).eq.1.o
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
                 id=id_cp(ip,icp,ist)
-                if(DEBUG)print*,'( couple ',id
+                if(DEBUG.EQ.1)
+      $              print*,'( couple ',id
       $              ,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm = xPAM
-Line 3205 
 c               distance = distance / RC
+Line 3378 
 c               distance = distance / RC
                    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
-                   clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
+      $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))
-      $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
                 endif
          continue
              enddo               !end loop on couples on plane icp
- c            if(distmin.le.clinc)then     !QUIQUI
+             if(distmin.le.clincnewc)then
-             if(distmin.le.clincnewc)then     !QUIQUI
  *              -----------------------------------
                 xm(nplanes-ip+1) = xmm !<<<
                 ym(nplanes-ip+1) = ymm !<<<
-Line 3225 
 c            if(distmin.le.clinc)then
+Line 3396 
 c            if(distmin.le.clinc)then
                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
  *              -----------------------------------
                 CP_STORE(nplanes-ip+1,ibest)=idm
-                if(DEBUG)print*,'%%%% included couple ',idm
+                if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
-      $              ,' (dist.= ',distmin,', cut ',clinc,' )'
+      $              ,' (dist.= ',distmin,', cut ',clincnewc,' )'
                 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
-Line 3252 
 c            if(DEBUG)print*,'>>>> try t
+Line 3422 
 c            if(DEBUG)print*,'>>>> try t
              do icp=1,ncp_plane(ip) !loop on cluster inside couples
                 icx=clx(ip,icp)
                 icy=cly(ip,icp)
+                if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
                 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
-Line 3270 
 c     $              AXV_STORE(nplanes-i
+Line 3441 
 c     $              AXV_STORE(nplanes-i
       $              )
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG)print*,'( cl-X ',icx
+                if(DEBUG.EQ.1)
+      $              print*,'( cl-X ',icx
       $              ,' in cp ',id,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
-Line 3303 
 c     $              0.,AYV_STORE(nplane
+Line 3475 
 c     $              0.,AYV_STORE(nplane
       $              )
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG)print*,'( cl-Y ',icy
+                if(DEBUG.EQ.1)
+      $              print*,'( cl-Y ',icy
       $              ,' in cp ',id,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
-Line 3323 
 c                 dedxmm = sgnl(icy)  !(
+Line 3496 
 c                 dedxmm = sgnl(icy)  !(
          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
-Line 3348 
 c               if(cl_used(icl).eq.1.or.
+Line 3519 
 c               if(cl_used(icl).eq.1.or.
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG)print*,'( cl-s ',icl
+                if(DEBUG.EQ.1)
-      $              ,' ) distance ',distance,'<',distmin,' ?'
+      $              print*,'( cl-s ',icl
+      $              ,' ) distance ',distance
                 if(distance.lt.distmin)then
-                   if(DEBUG)print*,'YES'
                    xmm_A = xPAM_A
                    ymm_A = yPAM_A
                    zmm_A = zPAM_A
-Line 3372 
 c               distance = distance / RC
+Line 3543 
 c               distance = distance / RC
                 endif
          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=
-Line 3386 
 c     anche qui: non ci vuole clinc???
+Line 3554 
 c     anche qui: non ci vuole clinc???
       $                 10*
       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
                 endif
- c     QUIQUI------------
+                if(distmin.le.clincnew)then
-                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
+                      if(DEBUG.EQ.1)
+      $                    print*,'%%%% included X-cl ',iclm
       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
       $                    ,', dist.= ',distmin
-      $                    ,', cut ',clinc,' )'
+      $                    ,', cut ',clincnew,' )'
                    else
                       YGOOD(nplanes-ip+1)=1.
                       resy(nplanes-ip+1)=rymm
-                      if(DEBUG)print*,'%%%% included Y-cl ',iclm
+                      if(DEBUG.EQ.1)
+      $                    print*,'%%%% included Y-cl ',iclm
       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
       $                    ,', dist.= ', distmin
-      $                    ,', cut ',clinc,' )'
+      $                    ,', cut ',clincnew,' )'
                    endif
- c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
  *     ----------------------------
                    xm_A(nplanes-ip+1) = xmm_A
                    ym_A(nplanes-ip+1) = ymm_A
-Line 3430 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3596 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
        return
        end
  ***************************************************
  *                                                 *
  *                                                 *
-Line 3439 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3606 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
  *                                                 *
  **************************************************
  *
-       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
-Line 3557 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3654 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
              DEDX_Y(IP,IT) = 0
              CLTRX(IP,IT) = 0
              CLTRY(IP,IT) = 0
+             multmaxx(ip,it) = 0
+             seedx(ip,it)    = 0
+             xpu(ip,it)      = 0
+             multmaxy(ip,it) = 0
+             seedy(ip,it)    = 0
+             ypu(ip,it)      = 0
           enddo
           do ipa=1,5
              AL_nt(IPA,IT) = 0
-Line 3576 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3679 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
          ys(1,ip)=0
          ys(2,ip)=0
          sgnlys(ip)=0
+         sxbad(ip)=0
+         sybad(ip)=0
+         multmaxsx(ip)=0
+         multmaxsy(ip)=0
        enddo
        end
-Line 3698 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3805 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
        integer ssensor,sladder
        pig=acos(-1.)
  *     -------------------------------------
        chi2_nt(ntr)        = sngl(chi2)
        nstep_nt(ntr)       = nstep
-Line 3736 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3845 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
           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)
+ ccc         print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
           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)
+ c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
-          if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
+ c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
-          tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
+ c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
-          angx     = 180.*atan(tgtemp)/acos(-1.)
+ c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
-          tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001
+ c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001
-          angy    = 180.*atan(tgtemp)/acos(-1.)
+ c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
+          angx = effectiveangle(ax,2*ip,bfy)
+          angy = effectiveangle(ay,2*ip-1,bfx)
- c         print*,'* ',ip,bfx,bfy,angx,angy
           id  = CP_STORE(ip,IDCAND) ! couple id
           icl = CLS_STORE(ip,IDCAND)
-Line 3771 
 c         print*,'* ',ip,bfx,bfy,angx,an
+Line 3886 
 c         print*,'* ',ip,bfx,bfy,angx,an
  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$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
- c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)
- c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
- c$$$            ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
-             xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
-             ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
+             if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
-             if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
+                cl_used(cltrx(ip,ntr)) = 1 !tag used clusters
-      $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
-             if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
+                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
-      $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
+                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
+      $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
+                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
+      $              +10000*mult(cltrx(ip,ntr))
+                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
+      $              /clsigma(indmax(cltrx(ip,ntr)))
+                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
+                xpu(ip,ntr)      = corr
+             endif
+             if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then
+                cl_used(cltry(ip,ntr)) = 1 !tag used clusters
+                ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
+                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
+      $              dedx_y(ip,ntr)=-dedx_y(ip,ntr)
+                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
+      $              +10000*mult(cltry(ip,ntr))
+                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
+      $              /clsigma(indmax(cltry(ip,ntr)))
+                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
+                ypu(ip,ntr)      = corr
+             endif
           elseif(icl.ne.0)then
+             cl_used(icl) = 1    !tag used clusters
              if(mod(VIEW(icl),2).eq.0)then
                 cltrx(ip,ntr)=icl
- c$$$               nnnnn = npfastrips(icl,PFA,angx)
- c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)
                 xbad(ip,ntr) = nbadstrips(4,icl)
                 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
+                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
+      $                         +10000*mult(cltrx(ip,ntr))
+                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
+      $           /clsigma(indmax(cltrx(ip,ntr)))
+                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
+                xpu(ip,ntr)      = corr
              elseif(mod(VIEW(icl),2).eq.1)then
                 cltry(ip,ntr)=icl
- c$$$               nnnnn = npfastrips(icl,PFA,angy)
- c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)
                 ybad(ip,ntr) = nbadstrips(4,icl)
                 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
+                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
+      $                         +10000*mult(cltry(ip,ntr))
+                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
+      $           /clsigma(indmax(cltry(ip,ntr)))
+                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
+                ypu(ip,ntr)      = corr
              endif
           endif
        enddo
+       if(DEBUG.eq.1)then
- c$$$      print*,(xgood(i),i=1,6)
+          print*,'> STORING TRACK ',ntr
- c$$$      print*,(ygood(i),i=1,6)
+          print*,'clusters: '
- c$$$      print*,(ls(i,ntr),i=1,6)
+          do ip=1,6
- c$$$      print*,(dedx_x(i,ntr),i=1,6)
+             print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
- c$$$      print*,(dedx_y(i,ntr),i=1,6)
+          enddo
- c$$$      print*,'-----------------------'
+          print*,'dedx: '
+          do ip=1,6
+             print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
+          enddo
+       endif
        end
-Line 3840 
 c         if( mask_view(iv).ne.0 )good2(
+Line 3995 
 c         if( mask_view(iv).ne.0 )good2(
           good2(iv) = good2(iv) + mask_view(iv)*2**8
        enddo
+       if(DEBUG.eq.1)then
+          print*,'> STORING SINGLETS '
+       endif
        do icl=1,nclstr1
+          ip=nplanes-npl(VIEW(icl))+1
           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
+                sxbad(nclsx)  = nbadstrips(1,icl)
+                multmaxsx(nclsx) = maxs(icl)+10000*mult(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
+                sybad(nclsy)  = nbadstrips(1,icl)
+                multmaxsy(nclsy) = maxs(icl)+10000*mult(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)
+ *     --------------------------------------------------
+ *     per non perdere la testa...
+ *     whichtrack(icl) e` una variabile del common level1
+ *     che serve solo per sapere quali cluster sono stati
+ *     associati ad una traccia, e permettere di salvare
+ *     solo questi nell'albero di uscita
+ *     --------------------------------------------------
        enddo
        end
-Line 3944 
 c$$$               print*,'ys(2,nclsy)
+Line 4112 
 c$$$               print*,'ys(2,nclsy)
                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)
                 nnn=nnn+ptcloud_yz(iyz)
              enddo
-             do ipt=1,nnn
+             do ipt=1,min(ndblt_max_nt,nnn)
                 db_cloud_nt(ipt)=db_cloud(ipt)
               enddo
           endif
-Line 3957 
 c$$$               print*,'ys(2,nclsy)
+Line 4125 
 c$$$               print*,'ys(2,nclsy)
                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
                 nnn=nnn+ptcloud_xz(ixz)
              enddo
-             do ipt=1,nnn
+             do ipt=1,min(ntrpt_max_nt,nnn)
                tr_cloud_nt(ipt)=tr_cloud(ipt)
               enddo
           endif

 Legend:



Removed from v.1.27
 


changed lines


 
Added in v.1.40
 Legend:



Removed from v.1.27
 


changed lines


 
Added in v.1.40
-Removed from v.1.27
+Added in v.1.40

	ViewVC Help
Powered by ViewVC 1.1.23