/[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.1 by mocchiut,
Fri May 19 13:15:54 2006 UTC
+revision 1.38 by pam-fi,
Tue Dec 23 11:28:36 2008 UTC
 Line 12
        subroutine track_finding(iflag)
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
        include 'common_mech.f'
        include 'common_xyzPAM.f'
        include 'common_mini_2.f'
        include 'calib.f'
-       include 'level1.f'
        include 'level2.f'
-       include 'momanhough_init.f'
+ c      print*,'======================================================'
+ c$$$      do ic=1,NCLSTR1
+ c$$$         if(.false.
+ c$$$     $        .or.nsatstrips(ic).gt.0
+ c$$$c     $        .or.nbadstrips(0,ic).gt.0
+ c$$$c     $        .or.nbadstrips(4,ic).gt.0
+ c$$$c     $        .or.nbadstrips(3,ic).gt.0
+ c$$$     $        .or..false.)then
+ c$$$            print*,'--- cl-',ic,' ------------------------'
+ c$$$            istart = INDSTART(IC)
+ c$$$            istop  = TOTCLLENGTH
+ c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
+ c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
+ c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
+ c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
+ c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
+ c$$$            print*,'view ',VIEW(ic)
+ c$$$            print*,'maxs ',MAXS(ic)
+ c$$$            print*,'COG4 ',cog(4,ic)
+ c$$$            ff = fbad_cog(4,ic)
+ c$$$            print*,'fbad ',ff
+ c$$$            print*,(CLBAD(i),i=istart,istop)
+ c$$$            bb=nbadstrips(0,ic)
+ c$$$            print*,'#BAD (tot)',bb
+ c$$$            bb=nbadstrips(4,ic)
+ c$$$            print*,'#BAD (4)',bb
+ c$$$            bb=nbadstrips(3,ic)
+ c$$$            print*,'#BAD (3)',bb
+ c$$$            ss=nsatstrips(ic)
+ c$$$            print*,'#saturated ',ss
+ c$$$         endif
+ c$$$      enddo
-       logical DEBUG
-       common/dbg/DEBUG
  *-------------------------------------------------------------------------------
  *     STEP 1
  *-------------------------------------------------------------------------------
-Line 44
+Line 72
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
- c      iflag=0
        call cl_to_couples(iflag)
        if(iflag.eq.1)then        !bad event
-          goto 880               !fill ntp and go to next event
+          goto 880               !go to next event
        endif
+       if(ncp_tot.eq.0)goto 880  !go to next event
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *     selezione di tracce pulite per diagnostica
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- c$$$         if(DEBUG)then
- c$$$            do ip=1,nplanes
- c$$$               if(ncp_plane(ip).ne.1)good2=.false.
- c$$$            enddo
- c$$$c            if(good2.eq.0)goto 100!next event
- c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event
- c$$$         endif
  *-----------------------------------------------------
  *-----------------------------------------------------
  *     HOUGH TRASFORM
-Line 91 
 c$$$         endif
+Line 105 
 c$$$         endif
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
- c      iflag=0
        call cp_to_doubtrip(iflag)
        if(iflag.eq.1)then        !bad event
-          goto 880               !fill ntp and go to next event
+          goto 880               !go to next event
        endif
+       if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event
  *-------------------------------------------------------------------------------
-Line 123 
 c      iflag=0
+Line 138 
 c      iflag=0
  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
+ *     count number of hit planes
+       planehit=0
+       do np=1,nplanes
+         if(ncp_plane(np).ne.0)then
+           planehit=planehit+1
+         endif
+       enddo
+       if(planehit.lt.3) goto 880 ! exit
+       nptxz_min=x_min_start
+       nplxz_min=x_min_start
+       nptyz_min=y_min_start
+       nplyz_min=y_min_start
+       cutdistyz=cutystart
+       cutdistxz=cutxstart
- c      iflag=0
+ continue
        call doub_to_YZcloud(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_yz.eq.0.and.cutdistyz.lt.maxcuty)then
+          if(cutdistyz.lt.maxcuty/2)then
+             cutdistyz=cutdistyz+cutystep
+          else
+             cutdistyz=cutdistyz+(3*cutystep)
+          endif
+          if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
+          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
- c      iflag=0
+       if(nclouds_yz.eq.0)goto 880 !go to next event
+ 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
+ c$$$ 881  continue
+ c$$$*     if there is at least three planes on the Y view decreases cuts on X view
+ c$$$      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
+ c$$$     $     nplxz_min.ne.y_min_start)then
+ c$$$        nptxz_min=x_min_step
+ c$$$        nplxz_min=x_min_start-x_min_step
+ c$$$        goto 879
+ c$$$      endif
  return
        end
-Line 144 
 c      iflag=0
+Line 230 
 c      iflag=0
        subroutine track_fitting(iflag)
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
        include 'common_mech.f'
        include 'common_xyzPAM.f'
        include 'common_mini_2.f'
        include 'calib.f'
-       include 'level1.f'
        include 'level2.f'
-       include 'momanhough_init.f'
+ c      include 'momanhough_init.f'
-       logical DEBUG
-       common/dbg/DEBUG
        logical FIMAGE            !
+       real trackimage(NTRACKSMAX)
+       real*8 AL_GUESS(5)
  *-------------------------------------------------------------------------------
  *     STEP 4   (ITERATED until any other physical track isn't found)
-Line 184 
 c      iflag=0
+Line 269 
 c      iflag=0
  *
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
-          ntrk=0                 !counter of identified physical tracks
+ ccc         ntrk=0                 !counter of identified physical tracks
    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
           iimage=0               !track image
  *     ------------- select the best track -------------
-          rchi2best=1000000000.
+ c$$$         rchi2best=1000000000.
+ c$$$         do i=1,ntracks
+ c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
+ c$$$     $         RCHI2_STORE(i).gt.0)then
+ c$$$               ibest=i
+ c$$$               rchi2best=RCHI2_STORE(i)
+ c$$$            endif
+ c$$$         enddo
+ c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
+ *     -------------------------------------------------------
+ *     order track-candidates according to:
+ *     1st) decreasing n.points
+ *     2nd) increasing chi**2
+ *     -------------------------------------------------------
+          rchi2best=1000000000.
+          ndofbest=0
           do i=1,ntracks
-             if(RCHI2_STORE(i).lt.rchi2best.and.
+            ndof=0
-      $         RCHI2_STORE(i).gt.0)then
+            do ii=1,nplanes
+              ndof=ndof
+      $            +int(xgood_store(ii,i))
+      $            +int(ygood_store(ii,i))
+            enddo
+            if(ndof.gt.ndofbest)then
+              ibest=i
+              rchi2best=RCHI2_STORE(i)
+              ndofbest=ndof
+            elseif(ndof.eq.ndofbest)then
+              if(RCHI2_STORE(i).lt.rchi2best.and.
+      $            RCHI2_STORE(i).gt.0)then
                 ibest=i
                 rchi2best=RCHI2_STORE(i)
-             endif
+                ndofbest=ndof
+              endif
+            endif
           enddo
           if(ibest.eq.0)goto 880 !>> no good candidates
  *-------------------------------------------------------------------------------
  *     The best track candidate (ibest) is selected and a new fitting is performed.
-Line 224 
 c         iflag=0
+Line 344 
 c         iflag=0
              iimage=0
           endif
           if(icand.eq.0)then
-             print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
+             if(VERBOSE.EQ.1)then
-      $           ,ibest,iimage
+                print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
+      $              ,ibest,iimage
+             endif
              return
           endif
-Line 236 
 c         iflag=0
+Line 358 
 c         iflag=0
  *     **********************************************************
  *     ************************** FIT *** FIT *** FIT *** FIT ***
  *     **********************************************************
+          call guess()
+          do i=1,5
+             AL_GUESS(i)=AL(i)
+          enddo
           do i=1,5
              AL(i)=dble(AL_STORE(i,icand))
           enddo
+          IDCAND = icand         !fitted track-candidate
           ifail=0                !error flag in chi2 computation
           jstep=0                !# minimization steps
-          call mini_2(jstep,ifail)
+          iprint=0
-          if(ifail.ne.0) then
+ c         if(DEBUG.EQ.1)iprint=1
-             if(DEBUG)then
+          if(VERBOSE.EQ.1)iprint=1
+          if(DEBUG.EQ.1)iprint=2
+          call mini2(jstep,ifail,iprint)
+ cc         if(ifail.ne.0) 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 *** (mini_2) '
+      $              '*** MINIMIZATION FAILURE *** (after refinement) '
       $              ,iev
              endif
-             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 263 
 c         iflag=0
+Line 397 
 c         iflag=0
           endif
  c         rchi2=chi2/dble(ndof)
-          if(DEBUG)then
+          if(DEBUG.EQ.1)then
              print*,' '
              print*,'****** SELECTED TRACK *************'
              print*,'#         R. chi2        RIG'
-Line 279 
 c         rchi2=chi2/dble(ndof)
+Line 413 
 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 306 
 c     $           RCHI2_STORE(i).gt.0.an
+Line 534 
 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(DEBUG)
+             if(verbose.eq.1)
       $           print*,
       $           '** warning ** number of identified '//
       $           'tracks exceeds vector dimension '
-Line 323 
 cc            good2=.false.
+Line 549 
 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
- c$$$************************************************************
- c$$$
- c$$$      subroutine readmipparam
- c$$$
- c$$$      include 'commontracker.f'
- c$$$      include 'calib.f'
- c$$$
- c$$$      character*60 fname_param
- c$$$ 201  format('trk-LADDER',i1,'-mip.dat')
- c$$$      do ilad=1,nladders_view
- c$$$         write(fname_param,201)ilad
- c$$$         print *,'Opening file: ',fname_param
- c$$$         open(10,
- c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))
- c$$$     $        ,STATUS='UNKNOWN'
- c$$$     $        ,IOSTAT=iostat
- c$$$     $        )
- c$$$         if(iostat.ne.0)then
- c$$$            print*,'READMIPPARAM: *** Error in opening file ***'
- c$$$            return
- c$$$         endif
- c$$$         do iv=1,nviews
- c$$$            read(10,*
- c$$$     $           ,IOSTAT=iostat
- c$$$     $           )pip,
- c$$$     $            mip(int(pip),ilad)
- c$$$c            print*,ilad,iv,pip,mip(int(pip),ilad)
- c$$$         enddo
- c$$$         close(10)
- c$$$      enddo
- c$$$
- c$$$      return
- c$$$      end
- c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
- c$$$      subroutine readchargeparam
- c$$$
- c$$$
- c$$$      include 'commontracker.f'
- c$$$      include 'calib.f'
- c$$$
- c$$$      character*60 fname_param
- c$$$ 201  format('charge-l',i1,'.dat')
- c$$$      do ilad=1,nladders_view
- c$$$         write(fname_param,201)ilad
- c$$$         print *,'Opening file: ',fname_param
- c$$$         open(10,
- c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))
- c$$$     $        ,STATUS='UNKNOWN'
- c$$$     $        ,IOSTAT=iostat
- c$$$     $        )
- c$$$         if(iostat.ne.0)then
- c$$$            print*,'READCHARGEPARAM: *** Error in opening file ***'
- c$$$            return
- c$$$         endif
- c$$$         do ip=1,nplanes
- c$$$            read(10,*
- c$$$     $           ,IOSTAT=iostat
- c$$$     $           )pip,
- c$$$     $            kch(ip,ilad),cch(ip,ilad),sch(ip,ilad)
- c$$$c            print*,ilad,ip,pip,kch(ip,ilad),
- c$$$c     $           cch(ip,ilad),sch(ip,ilad)
- c$$$         enddo
- c$$$         close(10)
- c$$$      enddo
- c$$$
- c$$$      return
- c$$$      end
- c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
- c$$$      subroutine readetaparam
- c$$$*     -----------------------------------------
- c$$$*     read eta2,3,4 calibration parameters
- c$$$*     and fill variables:
- c$$$*
- c$$$*     eta2(netabin,nladders_view,nviews)
- c$$$*     eta3(2*netabin,nladders_view,nviews)
- c$$$*     eta4(2*netabin,nladders_view,nviews)
- c$$$*
- c$$$      include 'commontracker.f'
- c$$$      include 'calib.f'
- c$$$
- c$$$      character*40 fname_binning
- c$$$      character*40 fname_param
- c$$$c      character*120 cmd1
- c$$$c      character*120 cmd2
- c$$$
- c$$$
- c$$$******retrieve ANGULAR BINNING info
- c$$$      fname_binning='binning.dat'
- c$$$      print *,'Opening file: ',fname_binning
- c$$$      open(10,
- c$$$     $     FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning))
- c$$$     $     ,STATUS='UNKNOWN'
- c$$$     $     ,IOSTAT=iostat
- c$$$     $     )
- c$$$      if(iostat.ne.0)then
- c$$$         print*,'READETAPARAM: *** Error in opening file ***'
- c$$$         return
- c$$$      endif
- c$$$      print*,'---- ANGULAR BINNING ----'
- c$$$      print*,'Bin   -   angL   -   angR'
- c$$$ 101  format(i2,'       ',f6.2,'     ',f6.2)
- c$$$      do ibin=1,nangmax
- c$$$         read(10,*
- c$$$     $        ,IOSTAT=iostat
- c$$$     $        )xnn,angL(ibin),angR(ibin)
- c$$$         if(iostat.ne.0)goto 1000
- c$$$         write(*,101)int(xnn),angL(ibin),angR(ibin)
- c$$$      enddo
- c$$$ 1000 nangbin=int(xnn)
- c$$$      close(10)
- c$$$      print*,'-------------------------'
- c$$$
- c$$$
- c$$$
- c$$$      do ieta=2,4               !loop on eta 2,3,4
- c$$$******retrieve correction parameters
- c$$$ 200     format(' Opening eta',i1,' files...')
- c$$$         write(*,200)ieta
- c$$$
- c$$$ 201     format('eta',i1,'-bin',i1,'-l',i1,'.dat')
- c$$$ 202     format('eta',i1,'-bin',i2,'-l',i1,'.dat')
- c$$$         do iang=1,nangbin
- c$$$            do ilad=1,nladders_view
- c$$$               if(iang.lt.10)write(fname_param,201)ieta,iang,ilad
- c$$$               if(iang.ge.10)write(fname_param,202)ieta,iang,ilad
- c$$$c               print *,'Opening file: ',fname_param
- c$$$               open(10,
- c$$$     $             FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))
- c$$$     $              ,STATUS='UNKNOWN'
- c$$$     $              ,IOSTAT=iostat
- c$$$     $              )
- c$$$               if(iostat.ne.0)then
- c$$$                  print*,'READETAPARAM: *** Error in opening file ***'
- c$$$                  return
- c$$$               endif
- c$$$               do ival=1,netavalmax
- c$$$                  if(ieta.eq.2)read(10,*
- c$$$     $                 ,IOSTAT=iostat
- c$$$     $                 )
- c$$$     $                 eta2(ival,iang),
- c$$$     $                 (feta2(ival,iv,ilad,iang),iv=1,nviews)
- c$$$                  if(ieta.eq.3)read(10,*
- c$$$     $                 ,IOSTAT=iostat
- c$$$     $                 )
- c$$$     $                 eta3(ival,iang),
- c$$$     $                 (feta3(ival,iv,ilad,iang),iv=1,nviews)
- c$$$                  if(ieta.eq.4)read(10,*
- c$$$     $                 ,IOSTAT=iostat
- c$$$     $                 )
- c$$$     $                 eta4(ival,iang),
- c$$$     $                 (feta4(ival,iv,ilad,iang),iv=1,nviews)
- c$$$                  if(iostat.ne.0)then
- c$$$                     netaval=ival-1
- c$$$c$$$                     if(eta2(1,iang).ne.-eta2(netaval,iang))
- c$$$c$$$     $                    print*,'**** ERROR on parameters !!! ****'
- c$$$                     goto 2000
- c$$$                  endif
- c$$$               enddo
- c$$$ 2000          close(10)
- c$$$*               print*,'... done'
- c$$$            enddo
- c$$$         enddo
- c$$$
- c$$$      enddo                     !end loop on eta 2,3,4
- c$$$
- c$$$
- c$$$      return
- c$$$      end
- c$$$
  ************************************************************
  ************************************************************
-Line 556 
 c$$$
+Line 578 
 c$$$
  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
  *     angx   - Projected angle in x
  *     angy   - Projected angle in y
+ *     bfx    - x-component of magnetci field
+ *     bfy    - y-component of magnetic field
  *
  *     --------- COUPLES -------------------------------------------------------
  *     The couple defines a point in the space.
-Line 594 
 c$$$
+Line 618 
 c$$$
  *
  *
-       subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)
+       subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
- c*****************************************************
- c     07/10/2005 modified by elena vannuccini --> (1)
- c     01/02/2006 modified by elena vannuccini --> (2)
- c     02/02/2006 modified by Elena Vannuccini --> (3)
- c                (implemented new p.f.a.)
- c     03/02/2006 modified by Elena Vannuccini --> (4)
- c                (implemented variable resolution)
- c*****************************************************
        include 'commontracker.f'
-       include 'calib.f'
        include 'level1.f'
+       include 'calib.f'
        include 'common_align.f'
        include 'common_mech.f'
        include 'common_xyzPAM.f'
-       include 'common_resxy.f'
-       logical DEBUG
-       common/dbg/DEBUG
        integer icx,icy           !X-Y cluster ID
        integer sensor
        integer viewx,viewy
        character*4 PFAx,PFAy     !PFA to be used
-       real angx,angy            !X-Y angle
+       real ax,ay                !X-Y geometric angle
+       real angx,angy            !X-Y effective angle
+       real bfx,bfy              !X-Y b-field components
        real stripx,stripy
+       double precision 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
- c      double precision xi,yi,zi
- c      double precision xi_A,yi_A,zi_A
- c      double precision xi_B,yi_B,zi_B
        parameter (ndivx=30)
        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
+       if(sensor.lt.1.or.sensor.gt.2)then
+          print*,'xyz_PAM   ***ERROR*** wrong input '
+          print*,'sensor ',sensor
+          icx=0
+          icy=0
+       endif
  *     -----------------
  *     CLUSTER X
  *     -----------------
        if(icx.ne.0)then
-          viewx = VIEW(icx)
-          nldx = nld(MAXS(icx),VIEW(icx))
+          viewx   = VIEW(icx)
-          nplx = npl(VIEW(icx))
+          nldx    = nld(MAXS(icx),VIEW(icx))
-          resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
+          nplx    = npl(VIEW(icx))
+ c         resxPAM = RESXAV
-          stripx = float(MAXS(icx))
+          stripx  = float(MAXS(icx))
-          if(PFAx.eq.'COG1')then  !(1)
-             stripx = stripx      !(1)
+          if(
-             resxPAM = resxPAM    !(1)
+      $        viewx.lt.1.or.
-          elseif(PFAx.eq.'COG2')then
+      $        viewx.gt.12.or.
-             stripx = stripx + cog(2,icx)
+      $        nldx.lt.1.or.
-             resxPAM = resxPAM*fbad_cog(2,icx)
+      $        nldx.gt.3.or.
-          elseif(PFAx.eq.'ETA2')then
+      $        stripx.lt.1.or.
- c            cog2 = cog(2,icx)
+      $        stripx.gt.3072.or.
- c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)
+      $        .false.)then
- c            stripx = stripx + etacorr
+             print*,'xyz_PAM   ***ERROR*** wrong input '
-             stripx = stripx + pfa_eta2(icx,angx)            !(3)
+             print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
-             resxPAM = risx_eta2(angx)                       !   (4)
+             icx = 0
-             if(DEBUG.and.fbad_cog(2,icx).ne.1)
+             goto 10
-      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
-             resxPAM = resxPAM*fbad_cog(2,icx)
-          elseif(PFAx.eq.'ETA3')then                         !(3)
-             stripx = stripx + pfa_eta3(icx,angx)            !(3)
-             resxPAM = risx_eta3(angx)                       !   (4)
-             if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
-      $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
-             resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
-          elseif(PFAx.eq.'ETA4')then                         !(3)
-             stripx = stripx + pfa_eta4(icx,angx)            !(3)
-             resxPAM = risx_eta4(angx)                       !   (4)
-             if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
-      $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
-             resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
-          elseif(PFAx.eq.'ETA')then                          !(3)
-             stripx = stripx + pfa_eta(icx,angx)             !(3)
-             resxPAM = ris_eta(icx,angx)                     !   (4)
-             if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
-      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
- c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO
-             resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)
-          elseif(PFAx.eq.'COG')then           !(2)
-             stripx = stripx + cog(0,icx)     !(2)
-             resxPAM = risx_cog(angx)                        !   (4)
-             resxPAM = resxPAM*fbad_cog(0,icx)!(2)
-          else
-             print*,'*** Non valid p.f.a. (x) --> ',PFAx
           endif
-       endif
+ *        --------------------------
+ *        magnetic-field corrections
+ *        --------------------------
+          stripx  = stripx +  fieldcorr(viewx,bfy)
+          angx    = effectiveangle(ax,viewx,bfy)
+          call applypfa(PFAx,icx,angx,corr,res)
+          stripx  = stripx + corr
+          resxPAM = res
+  endif
  *     -----------------
  *     CLUSTER Y
  *     -----------------
        if(icy.ne.0)then
           viewy = VIEW(icy)
           nldy = nld(MAXS(icy),VIEW(icy))
           nply = npl(VIEW(icy))
-          resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
+ c         resyPAM = RESYAV
+          stripy = float(MAXS(icy))
+          if(
+      $        viewy.lt.1.or.
+      $        viewy.gt.12.or.
+      $        nldy.lt.1.or.
+      $        nldy.gt.3.or.
+      $        stripy.lt.1.or.
+      $        stripy.gt.3072.or.
+      $        .false.)then
+             print*,'xyz_PAM   ***ERROR*** wrong input '
+             print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
+             icy = 0
+             goto 20
+          endif
           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
-             print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
+             if(DEBUG.EQ.1) then
-      $           ,icx,icy
+                print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
+      $              ,icx,icy
+             endif
              goto 100
           endif
-          stripy = float(MAXS(icy))
-          if(PFAy.eq.'COG1')then !(1)
-             stripy = stripy     !(1)
-             resyPAM = resyPAM   !(1)
-          elseif(PFAy.eq.'COG2')then
-             stripy = stripy + cog(2,icy)
-             resyPAM = resyPAM*fbad_cog(2,icy)
-          elseif(PFAy.eq.'ETA2')then
- c            cog2 = cog(2,icy)
- c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)
- c            stripy = stripy + etacorr
-             stripy = stripy + pfa_eta2(icy,angy)            !(3)
-             resyPAM = risy_eta2(angy)                       !   (4)
-             resyPAM = resyPAM*fbad_cog(2,icy)
-             if(DEBUG.and.fbad_cog(2,icy).ne.1)
-      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
-          elseif(PFAy.eq.'ETA3')then                         !(3)
-             stripy = stripy + pfa_eta3(icy,angy)            !(3)
-             resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
-             if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
-      $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
-          elseif(PFAy.eq.'ETA4')then                         !(3)
-             stripy = stripy + pfa_eta4(icy,angy)            !(3)
-             resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
-             if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
-      $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
-          elseif(PFAy.eq.'ETA')then                          !(3)
-             stripy = stripy + pfa_eta(icy,angy)             !(3)
-             resyPAM = ris_eta(icy,angy)                     !   (4)
- c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
-             resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
-             if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)
-      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)
-          elseif(PFAy.eq.'COG')then
-             stripy = stripy + cog(0,icy)
-             resyPAM = risy_cog(angy)                        !   (4)
- c            resyPAM = ris_eta(icy,angy)                    !   (4)
-             resyPAM = resyPAM*fbad_cog(0,icy)
-          else
-             print*,'*** Non valid p.f.a. (x) --> ',PFAx
-          endif
-       endif
+ *        --------------------------
+ *        magnetic-field corrections
+ *        --------------------------
+          stripy  = stripy +  fieldcorr(viewy,bfx)
+          angy    = effectiveangle(ay,viewy,bfx)
+          call applypfa(PFAy,icy,angy,corr,res)
+          stripy  = stripy + corr
+          resyPAM = res
+  endif
  c===========================================================
  C     COUPLE
  C===========================================================
-Line 772 
 C=======================================
+Line 762 
 C=======================================
  c------------------------------------------------------------------------
  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-          xi = acoordsi(stripx,viewx)
+          if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
-          yi = acoordsi(stripy,viewy)
+      $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
-          zi = 0.
+             if(DEBUG.EQ.1) then
+                print*,'xyz_PAM (couple):',
+      $              ' WARNING: false X strip: strip ',stripx
+             endif
+          endif
+          xi = dcoordsi(stripx,viewx)
+          yi = dcoordsi(stripy,viewy)
+          zi = 0.D0
  c------------------------------------------------------------------------
  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-Line 810 
 c---------------------------------------
+Line 806 
 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 836 
 C=======================================
+Line 832 
 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 846 
 C=======================================
+Line 844 
 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 858 
 C=======================================
+Line 854 
 C=======================================
              nldy = nldx
              viewy = viewx - 1
-             xi   = acoordsi(stripx,viewx)
+             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.EQ.1) then
+                   print*,'xyz_PAM (X-singlet):',
+      $                 ' WARNING: false X strip: strip ',stripx
+                endif
+             endif
+             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 874 
 C=======================================
+Line 880 
 C=======================================
                 yi_B = yi
              endif
- c            print*,'X-cl ',icx,stripx,' --> ',xi
- c            print*,yi_A,' <--> ',yi_B
           else
+             if(DEBUG.EQ.1) then
-             print *,'routine xyz_PAM ---> not properly used !!!'
+                print *,'routine xyz_PAM ---> not properly used !!!'
-             print *,'icx = ',icx
+                print *,'icx = ',icx
-             print *,'icy = ',icy
+                print *,'icy = ',icy
+             endif
              goto 100
           endif
-Line 920 
 c     N.B. I convert angles from microra
+Line 925 
 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 930 
 c     (xPAM,yPAM,zPAM) = measured coordi
+Line 953 
 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 943 
 c---------------------------------------
+Line 969 
 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.EQ.1) then
-          print *,'routine xyz_PAM ---> not properly used !!!'
+             print *,'routine xyz_PAM ---> not properly used !!!'
-          print *,'icx = ',icx
+             print *,'icx = ',icx
-          print *,'icy = ',icy
+             print *,'icy = ',icy
+          endif
        endif
  continue
        end
+ ************************************************************************
+ *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
+ *     (done to be called from c/c++)
+ ************************************************************************
+       subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
+       include 'commontracker.f'
+       include 'level1.f'
+       include 'common_mini_2.f'
+       include 'common_xyzPAM.f'
+       include 'common_mech.f'
+       include 'calib.f'
+ *     flag to chose PFA
+ c$$$      character*10 PFA
+ c$$$      common/FINALPFA/PFA
+       integer icx,icy           !X-Y cluster ID
+       integer sensor
+       character*4 PFAx,PFAy     !PFA to be used
+       real ax,ay                !X-Y geometric angle
+       real bfx,bfy              !X-Y b-field components
+       ipx=0
+       ipy=0
+ c$$$      PFAx = 'COG4'!PFA
+ c$$$      PFAy = 'COG4'!PFA
+       if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
+             print*,'xyzpam: ***WARNING*** clusters ',icx,icy
+      $           ,' do not exists (n.clusters=',nclstr1,')'
+             icx = -1*icx
+             icy = -1*icy
+             return
+       endif
+       call idtoc(pfaid,PFAx)
+       call idtoc(pfaid,PFAy)
+       if(icx.ne.0.and.icy.ne.0)then
+          ipx=npl(VIEW(icx))
+          ipy=npl(VIEW(icy))
+          if( (nplanes-ipx+1).ne.ip )then
+             print*,'xyzpam: ***WARNING*** cluster ',icx
+      $           ,' does not belong to plane: ',ip
+             icx = -1*icx
+             return
+          endif
+          if( (nplanes-ipy+1).ne.ip )then
+             print*,'xyzpam: ***WARNING*** cluster ',icy
+      $           ,' does not belong to plane: ',ip
+             icy = -1*icy
+             return
+          endif
+          call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
+          xgood(ip) = 1.
+          ygood(ip) = 1.
+          resx(ip)  = resxPAM
+          resy(ip)  = resyPAM
+          xm(ip) = xPAM
+          ym(ip) = yPAM
+          zm(ip) = zPAM
+          xm_A(ip) = 0.D0
+          ym_A(ip) = 0.D0
+          xm_B(ip) = 0.D0
+          ym_B(ip) = 0.D0
+ c         zv(ip) = zPAM
+       elseif(icx.eq.0.and.icy.ne.0)then
+          ipy=npl(VIEW(icy))
+          if( (nplanes-ipy+1).ne.ip )then
+             print*,'xyzpam: ***WARNING*** cluster ',icy
+      $           ,' does not belong to plane: ',ip
+             icy = -1*icy
+             return
+          endif
+          call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
+          xgood(ip) = 0.
+          ygood(ip) = 1.
+          resx(ip)  = 1000.
+          resy(ip)  = resyPAM
+ c$$$         xm(ip) = -100.
+ c$$$         ym(ip) = -100.
+ 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
+          ym_B(ip) = yPAM_B
+ c         zv(ip) = (zPAM_A+zPAM_B)/2.
+       elseif(icx.ne.0.and.icy.eq.0)then
+          ipx=npl(VIEW(icx))
+          if( (nplanes-ipx+1).ne.ip )then
+             print*,'xyzpam: ***WARNING*** cluster ',icx
+      $           ,' does not belong to plane: ',ip
+             icx = -1*icx
+             return
+          endif
+          call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
+          xgood(ip) = 1.
+          ygood(ip) = 0.
+          resx(ip)  = resxPAM
+          resy(ip)  = 1000.
+ c$$$         xm(ip) = -100.
+ c$$$         ym(ip) = -100.
+ 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
+          ym_B(ip) = yPAM_B
+ c         zv(ip) = (zPAM_A+zPAM_B)/2.
+       else
+          il = 2
+          if(lad.ne.0)il=lad
+          is = 1
+          if(sensor.ne.0)is=sensor
+          xgood(ip) = 0.
+          ygood(ip) = 0.
+          resx(ip)  = 1000.
+          resy(ip)  = 1000.
+          xm(ip) = -100.
+          ym(ip) = -100.
+          zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
+          xm_A(ip) = 0.
+          ym_A(ip) = 0.
+          xm_B(ip) = 0.
+          ym_B(ip) = 0.
+ c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
+       endif
+       if(DEBUG.EQ.1)then
+   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)
+       endif
+       end
  ********************************************************************************
  ********************************************************************************
-Line 977 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 1175 
 c         print*,'A-(',xPAM_A,yPAM_A,')
  *
  ********************************************************************************
-       real function distance_to(XPP,YPP)
+       real function distance_to(rXPP,rYPP)
        include 'common_xyzPAM.f'
-Line 986 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 1184 
 c         print*,'A-(',xPAM_A,yPAM_A,')
  *     ( 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 1031 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 1234 
 c         print*,'A-(',xPAM_A,yPAM_A,')
           endif
           distance=
-      $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
+      $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
+ cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
           distance=dsqrt(distance)
- c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
- c$$$     $        ,' --- distance_to --- ',xpp,ypp
- c$$$         print*,' resolution ',re
  *     ----------------------
        elseif(
-      +     xPAM.ne.0.and.
+ 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 1056 
 c$$$         print*,' resolution ',re
+Line 1257 
 c$$$         print*,' resolution ',re
  *     ----------------------
           distance=
-      $        ((xPAM-XPP)/resxPAM)**2
+      $       ((xPAM-XPP))**2 !QUIQUI
-      $        +
+      $       +
-      $        ((yPAM-YPP)/resyPAM)**2
+      $       ((yPAM-YPP))**2
+ c$$$     $        ((xPAM-XPP)/resxPAM)**2
+ c$$$     $        +
+ c$$$     $        ((yPAM-YPP)/resyPAM)**2
           distance=dsqrt(distance)
- c$$$         print*,xPAM,yPAM,zPAM
- c$$$     $        ,' --- distance_to --- ',xpp,ypp
- c$$$         print*,' resolution ',resxPAM,resyPAM
        else
-          print*
-      $        ,' function distance_to ---> wrong usage!!!'
-          print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
-          print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
-      $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
        endif
        distance_to = sngl(distance)
-Line 1111 
 c$$$         print*,' resolution ',resxP
+Line 1307 
 c$$$         print*,' resolution ',resxP
        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 1136 
 c$$$         print*,' resolution ',resxP
+Line 1332 
 c$$$         print*,' resolution ',resxP
  c------------------------------------------------------------------------
  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-                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 1163 
 c---------------------------------------
+Line 1359 
 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 1263 
 c     $              ,iv,xvv(iv),yvv(iv)
+Line 1457 
 c     $              ,iv,xvv(iv),yvv(iv)
  *     it returns the plane number
  *
        include 'commontracker.f'
+       include 'level1.f'
  c      include 'common_analysis.f'
        include 'common_momanhough.f'
-Line 1288 
 c      include 'common_analysis.f'
+Line 1483 
 c      include 'common_analysis.f'
        is_cp=0
        if(id.lt.0)is_cp=1
        if(id.gt.0)is_cp=2
-       if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
        return
        end
-Line 1300 
 c      include 'common_analysis.f'
+Line 1494 
 c      include 'common_analysis.f'
  *     it returns the id number ON THE PLANE
  *
        include 'commontracker.f'
+       include 'level1.f'
  c      include 'common_analysis.f'
        include 'common_momanhough.f'
-Line 1328 
 c      include 'common_analysis.f'
+Line 1523 
 c      include 'common_analysis.f'
  *     positive if sensor =2
  *
        include 'commontracker.f'
+       include 'level1.f'
  c      include 'calib.f'
  c      include 'level1.f'
  c      include 'common_analysis.f'
-Line 1357 
 c      include 'common_analysis.f'
+Line 1553 
 c      include 'common_analysis.f'
  *************************************************************************
  *************************************************************************
  *************************************************************************
- c$$$      subroutine book_debug
- c$$$
- c$$$      include 'commontracker.f'
- c$$$      include 'common_momanhough.f'
- c$$$      include 'common_level2debug.f'
- c$$$
- c$$$      character*35 block1,block2,block3!,block4
- c$$$     $     ,block5!,block6
- c$$$
- c$$$* * * * * * * * * * * * * * * * * * * * * * * *
- c$$$*     HOUGH TRANSFORM PARAMETERS
- c$$$
- c$$$      call HBOOK2(1003
- c$$$     $     ,'y vs tg thyz'
- c$$$     $     ,300,-1.,1.         !x
- c$$$     $     ,3000,-70.,70.,0.)   !y
- c$$$
- c$$$      call HBOOK1(1004
- c$$$     $     ,'Dy'
- c$$$     $     ,100,0.,2.,0.)
- c$$$
- c$$$      call HBOOK1(1005
- c$$$     $     ,'D thyz'
- c$$$     $     ,100,0.,.05,0.)
- c$$$
- c$$$
- c$$$
- c$$$*     DEBUG ntuple:
- c$$$      call HBNT(ntp_level2+1,'LEVEL2',' ')
- c$$$
- c$$$      call HBNAME(ntp_level2+1,'EVENT',good2_nt,
- c$$$     $     'GOOD2:L,NEV2:I')
- c$$$
- c$$$ 411  format('NDBLT:I::[0,',I5,']')
- c$$$      write(block1,411) ndblt_max_nt
- c$$$      call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt,
- c$$$     $     block1//'
- c$$$     $     ,ALFAYZ1(NDBLT):R
- c$$$     $     ,ALFAYZ2(NDBLT):R
- c$$$     $     ,DB_CLOUD(NDBLT):I
- c$$$     $     ')
- c$$$
- c$$$ 412  format('NTRPT:I::[0,',I5,']')
- c$$$      write(block2,412) ntrpt_max_nt
- c$$$      call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt,
- c$$$     $     block2//'
- c$$$     $     ,ALFAXZ1(NTRPT):R
- c$$$     $     ,ALFAXZ2(NTRPT):R
- c$$$     $     ,ALFAXZ3(NTRPT):R
- c$$$     $     ,TR_CLOUD(NTRPT):I
- c$$$     $     ')
- c$$$
- c$$$
- c$$$ 413  format('NCLOUDS_YZ:I::[0,',I4,']')
- c$$$c$$$ 414  format('DB_CLOUD(',I4,'):I')
- c$$$      write(block3,413) ncloyz_max
- c$$$c$$$      write(block4,414) ndblt_max_nt
- c$$$      call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ,
- c$$$     $     block3//'
- c$$$     $     ,ALFAYZ1_AV(NCLOUDS_YZ):R
- c$$$     $     ,ALFAYZ2_AV(NCLOUDS_YZ):R
- c$$$     $     ,PTCLOUD_YZ(NCLOUDS_YZ):I'
- c$$$c$$$     $     ,'//block4
- c$$$     $     )
- c$$$
- c$$$ 415  format('NCLOUDS_XZ:I::[0,',I4,']')
- c$$$c$$$ 416  format('TR_CLOUD(',I5,'):I')
- c$$$      write(block5,415) ncloxz_max
- c$$$c$$$      write(block6,416) ntrpt_max_nt
- c$$$      call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ,
- c$$$     $     block5//'
- c$$$     $     ,ALFAXZ1_AV(NCLOUDS_XZ):R
- c$$$     $     ,ALFAXZ2_AV(NCLOUDS_XZ):R
- c$$$     $     ,ALFAXZ3_AV(NCLOUDS_XZ):R
- c$$$     $     ,PTCLOUD_XZ(NCLOUDS_XZ):I'
- c$$$c$$$     $     ,'//block6
- c$$$     $     )
- c$$$
- c$$$
- c$$$      return
- c$$$      end
- ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
- *
- *
- *
- *
- *
- *
- *
- *
- *
- ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
- c$$$      subroutine book_level2
- c$$$c*****************************************************
- c$$$cccccc 11/9/2005 modified by david fedele
- c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2)
- c$$$c*****************************************************
- c$$$
- c$$$      include 'commontracker.f'
- c$$$      include 'common_momanhough.f'
- c$$$      include 'level2.f'
- c$$$
- c$$$      character*35 block1,block2
- c$$$
- c$$$c      print*,'__________ booking LEVEL2 n-tuple __________'
- c$$$
- c$$$*     LEVEL1 ntuple:
- c$$$      call HBNT(ntp_level2,'LEVEL2',' ')
- c$$$
- c$$$c*****************************************************
- c$$$cccccc 11/9/2005 modified by david fedele
- c$$$c      call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I')
- c$$$cccccc 06/10/2005 modified by elena vannuccini
- c$$$c      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I
- c$$$c     $     ,WHIC_CALIB:I,SWCODE:I')
- c$$$      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I
- c$$$     $     ,WHICH_CALIB:I,SWCODE:I,CRC(12):L')
- c$$$c*********************************************************
- c$$$
- c$$$
- c$$$c$$$# ifndef TEST2003
- c$$$c$$$
- c$$$c$$$      call HBNAME(ntp_level2,'CPU',pkt_type
- c$$$c$$$     $     ,'PKT_TYPE:I::[0,50]
- c$$$c$$$     $     ,PKT_NUM:I
- c$$$c$$$     $     ,OBT:I'//
- c$$$c$$$c********************************************************
- c$$$c$$$cccccc 11/9/2005 modified by david fedele
- c$$$c$$$c     $     ,WHICH_CALIB:I::[0,50]')
- c$$$c$$$     $     ',CPU_CRC:L')
- c$$$c$$$c********************************************************
- c$$$c$$$
- c$$$c$$$# endif
- c$$$
- c$$$ 417  format('NTRK:I::[0,',I4,']')
- c$$$ 418  format(',IMAGE(NTRK):I::[0,',I4,']')
- c$$$      write(block1,417)NTRKMAX
- c$$$      write(block2,418)NTRKMAX
- c$$$      call HBNAME(ntp_level2,'TRACKS',NTRK,
- c$$$     $     block1//
- c$$$     $     block2//'
- c$$$     $     ,XM(6,NTRK):R
- c$$$     $     ,YM(6,NTRK):R
- c$$$     $     ,ZM(6,NTRK):R
- c$$$     $     ,RESX(6,NTRK):R
- c$$$     $     ,RESY(6,NTRK):R
- c$$$     $     ,AL(5,NTRK):R
- c$$$     $     ,COVAL(5,5,NTRK):R
- c$$$     $     ,CHI2(NTRK):R
- c$$$     $     ,XGOOD(6,NTRK):I::[0,1]
- c$$$     $     ,YGOOD(6,NTRK):I::[0,1]
- c$$$     $     ,XV(6,NTRK):R
- c$$$     $     ,YV(6,NTRK):R
- c$$$     $     ,ZV(6,NTRK):R
- c$$$     $     ,AXV(6,NTRK):R
- c$$$     $     ,AYV(6,NTRK):R'//
- c$$$c*****************************************************
- c$$$cccccc 11/9/2005 modified by david fedele
- c$$$c     $     ,DEDXP(6,NTRK):R'//
- c$$$c     $     ')
- c$$$     $     ',DEDX_X(6,NTRK):R
- c$$$     $     ,DEDX_Y(6,NTRK):R'//
- c$$$c****************************************************
- c$$$cccccc 06/10/2005 modified by elena vannuccini
- c$$$c     $     ,CRC(12):L
- c$$$     $     ',BdL(NTRK):R'
- c$$$     $     )
- c$$$c****************************************************
- c$$$
- c$$$
- c$$$      call HBNAME(ntp_level2,'SINGLETX',nclsx,
- c$$$c*****************************************************
- c$$$cccccc 11/9/2005 modified by david fedele
- c$$$c     $     'NCLSX(6):I,NCLSY(6):I')
- c$$$     $     'NCLSX:I::[0,500],PLANEX(NCLSX):I
- c$$$     $     ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2)
- c$$$c    $     ,XS(NCLSX):R,SGNLXS(NCLSX):R')   !(2)
- c$$$      call HBNAME(ntp_level2,'SINGLETY',nclsy,
- c$$$     $     'NCLSY:I::[0,500],PLANEY(NCLSY):I
- c$$$     $     ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2)
- c$$$c    $     ,YS(NCLSY):R,SGNLYS(NCLSY):R')   !(2)
- c$$$      return
- c$$$      end
- c$$$      subroutine fill_level2_clouds
- c$$$c*****************************************************
- c$$$c     29/11/2005 created by elena vannuccini
- c$$$c*****************************************************
- c$$$
- c$$$*     -------------------------------------------------------
- c$$$*     This routine fills the  variables related to the hough
- c$$$*     transform, for the debig n-tuple
- c$$$*     -------------------------------------------------------
- c$$$
- c$$$      include 'commontracker.f'
- c$$$      include 'common_momanhough.f'
- c$$$      include 'common_level2debug.f'
- c$$$      include 'level2.f'
- c$$$
- c$$$      good2_nt=.true.!good2
- c$$$c      nev2_nt=nev2
- c$$$
- c$$$      if(.false.
- c$$$     $     .or.ntrpt.gt.ntrpt_max_nt
- c$$$     $     .or.ndblt.gt.ndblt_max_nt
- c$$$     $     .or.NCLOUDS_XZ.gt.ncloxz_max
- c$$$     $     .or.NCLOUDS_yZ.gt.ncloyz_max
- c$$$     $     )then
- c$$$         good2_nt=.false.
- c$$$         ntrpt_nt=0
- c$$$         ndblt_nt=0
- c$$$         NCLOUDS_XZ_nt=0
- c$$$         NCLOUDS_YZ_nt=0
- c$$$      else
- c$$$         ndblt_nt=ndblt
- c$$$         ntrpt_nt=ntrpt
- c$$$         if(ndblt.ne.0)then
- c$$$            do id=1,ndblt
- c$$$               alfayz1_nt(id)=alfayz1(id) !Y0
- c$$$               alfayz2_nt(id)=alfayz2(id) !tg theta-yz
- c$$$c               db_cloud_nt(id)=db_cloud(id)
- c$$$            enddo
- c$$$         endif
- c$$$         if(ndblt.ne.0)then
- c$$$            do it=1,ntrpt
- c$$$               alfaxz1_nt(it)=alfaxz1(it) !X0
- c$$$               alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
- c$$$               alfaxz3_nt(it)=alfaxz3(it) !1/r
- c$$$c               tr_cloud_nt(it)=tr_cloud(it)
- c$$$            enddo
- c$$$         endif
- c$$$         nclouds_yz_nt=nclouds_yz
- c$$$         nclouds_xz_nt=nclouds_xz
- c$$$         if(nclouds_yz.ne.0)then
- c$$$            nnn=0
- c$$$            do iyz=1,nclouds_yz
- c$$$               ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
- c$$$               alfayz1_av_nt(iyz)=alfayz1_av(iyz)
- c$$$               alfayz2_av_nt(iyz)=alfayz2_av(iyz)
- c$$$               nnn=nnn+ptcloud_yz(iyz)
- c$$$            enddo
- c$$$            do ipt=1,nnn
- c$$$               db_cloud_nt(ipt)=db_cloud(ipt)
- c$$$            enddo
- c$$$c            print*,'#### ntupla #### ptcloud_yz '
- c$$$c     $           ,(ptcloud_yz(i),i=1,nclouds_yz)
- c$$$c            print*,'#### ntupla #### db_cloud '
- c$$$c     $           ,(db_cloud(i),i=1,nnn)
- c$$$         endif
- c$$$         if(nclouds_xz.ne.0)then
- c$$$            nnn=0
- c$$$            do ixz=1,nclouds_xz
- c$$$               ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
- c$$$               alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
- c$$$               alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
- c$$$               alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
- c$$$               nnn=nnn+ptcloud_xz(ixz)
- c$$$            enddo
- c$$$            do ipt=1,nnn
- c$$$               tr_cloud_nt(ipt)=tr_cloud(ipt)
- c$$$            enddo
- c$$$c            print*,'#### ntupla #### ptcloud_xz '
- c$$$c     $           ,(ptcloud_xz(i),i=1,nclouds_xz)
- c$$$c            print*,'#### ntupla #### tr_cloud '
- c$$$c     $           ,(tr_cloud(i),i=1,nnn)
- c$$$         endif
- c$$$      endif
- c$$$      end
  ***************************************************
-Line 1639 
 c$$$      end
+Line 1567 
 c$$$      end
        subroutine cl_to_couples(iflag)
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
-       include 'momanhough_init.f'
+ c      include 'momanhough_init.f'
        include 'calib.f'
-       include 'level1.f'
+ c      include 'level1.f'
-       logical DEBUG
-       common/dbg/DEBUG
  *     output flag
  *     --------------
-Line 1654 
 c$$$      end
+Line 1580 
 c$$$      end
  *     --------------
        integer iflag
-       integer badseed,badcl
+       integer badseed,badclx,badcly
+       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 1670 
 c$$$      end
+Line 1599 
 c$$$      end
           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
+          mask_view(iv) = 0      !all included
        enddo
+ *     count number of cluster per view
+       do icl=1,nclstr1
+          ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1
+       enddo
+ *     mask views with too many clusters
+       do iv=1,nviews
+          if( ncl_view(iv).gt. nclusterlimit)then
+ c            mask_view(iv) = 1
+             mask_view(iv) = mask_view(iv) + 2**0
+             if(DEBUG.EQ.1)
+      $        print*,' * WARNING * cl_to_couple: n.clusters > '
+      $        ,nclusterlimit,' on view ', iv,' --> masked!'
+          endif
+       enddo
  *     start association
        ncouples=0
        do icx=1,nclstr1          !loop on cluster (X)
           if(mod(VIEW(icx),2).eq.1)goto 10
+          if(cl_used(icx).ne.0)goto 10
+ *     ----------------------------------------------------
+ *     jump masked views (X VIEW)
+ *     ----------------------------------------------------
+          if( mask_view(VIEW(icx)).ne.0 ) goto 10
  *     ----------------------------------------------------
  *     cut on charge (X VIEW)
-          if(dedx(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
+ *     ----------------------------------------------------
  *     cut BAD (X VIEW)
+ *     ----------------------------------------------------
           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
           ifirst=INDSTART(icx)
           if(icx.ne.nclstr1) then
-Line 1693 
 c$$$      end
+Line 1651 
 c$$$      end
           else
              ilast=TOTCLLENGTH
           endif
-          badcl=badseed
+          badclx=badseed
           do igood=-ngoodstr,ngoodstr
              ibad=1
              if((INDMAX(icx)+igood).gt.ifirst.and.
-Line 1703 
 c$$$      end
+Line 1661 
 c$$$      end
       $              nvk(MAXS(icx)+igood),
       $              nst(MAXS(icx)+igood))
              endif
-             badcl=badcl*ibad
+             badclx=badclx*ibad
           enddo
+ *     ----------------------------------------------------
+ *     >>> eliminato il taglio sulle BAD <<<
+ *     ----------------------------------------------------
  c     if(badcl.eq.0)then
  c     cl_single(icx)=0
  c     goto 10
-Line 1717 
 c     endif
+Line 1678 
 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)
+ *     ----------------------------------------------------
+             if( mask_view(VIEW(icy)).ne.0 ) goto 20
+ *     ----------------------------------------------------
  *     cut on charge (Y VIEW)
-             if(dedx(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
+ *     ----------------------------------------------------
  *     cut BAD (Y VIEW)
+ *     ----------------------------------------------------
              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
              ifirst=INDSTART(icy)
              if(icy.ne.nclstr1) then
-Line 1732 
 c     endif
+Line 1703 
 c     endif
              else
                 ilast=TOTCLLENGTH
              endif
-             badcl=badseed
+             badcly=badseed
              do igood=-ngoodstr,ngoodstr
                 ibad=1
                 if((INDMAX(icy)+igood).gt.ifirst.and.
-Line 1741 
 c     endif
+Line 1712 
 c     endif
       $              ibad=BAD(VIEW(icy),
       $              nvk(MAXS(icy)+igood),
       $              nst(MAXS(icy)+igood))
-                badcl=badcl*ibad
+                badcly=badcly*ibad
              enddo
+ *     ----------------------------------------------------
+ *     >>> eliminato il taglio sulle BAD <<<
+ *     ----------------------------------------------------
  c     if(badcl.eq.0)then
  c     cl_single(icy)=0
  c     goto 20
  c     endif
  *     ----------------------------------------------------
              cl_good(icy)=1
              nply=npl(VIEW(icy))
              nldy=nld(MAXS(icy),VIEW(icy))
-Line 1760 
 c     endif
+Line 1733 
 c     endif
  *     geometrical consistency (same plane and ladder)
              if(nply.eq.nplx.and.nldy.eq.nldx)then
  *     charge correlation
-                ddd=(dedx(icy)
+ *     (modified to be applied only below saturation... obviously)
-      $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
-                ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
+                if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
-                cut=chcut*sch(nplx,nldx)
+      $              .and.
-                if(abs(ddd).gt.cut)goto 20 !charge not consistent
+      $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
+      $              .and.
+      $              (badclx.eq.1.and.badcly.eq.1)
- *     ------------------> COUPLE <------------------
+      $              .and.
- *     check to do not overflow vector dimentions
+      $              .true.)then
-                if(ncp_plane(nplx).gt.ncouplemax)then
-                   if(DEBUG)print*,
+                   ddd=(sgnl(icy)
-      $                    ' ** warning ** number of identified'//
+      $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
-      $                    ' couples on plane ',nplx,
+                   ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
-      $                    ' exceeds vector dimention'//
-      $                    ' ( ',ncouplemax,' )'
+ c                  cut = chcut * sch(nplx,nldx)
- c     good2=.false.
- c     goto 880   !fill ntp and go to next event
+                   sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
-                   iflag=1
+      $                 -kch(nplx,nldx)*cch(nplx,nldx))
-                   return
+                   sss=sss/sqrt(kch(nplx,nldx)**2+1)
+                   cut = chcut * (16 + sss/50.)
+                   if(abs(ddd).gt.cut)then
+                      goto 20    !charge not consistent
+                   endif
                 endif
-                if(ncp_plane(nplx).eq.ncouplemax)then
+                if(ncp_plane(nplx).gt.ncouplemax)then
-                   if(DEBUG)print*,
+                   if(verbose.eq.1)print*,
       $                 '** warning ** number of identified '//
       $                 'couples on plane ',nplx,
       $                 'exceeds vector dimention '
-      $                 ,'( ',ncouplemax,' )'
+      $                 ,'( ',ncouplemax,' ) --> masked!'
- c     good2=.false.
+ c                  mask_view(nviewx(nplx)) = 2
- c     goto 880   !fill ntp and go to next event
+ c                  mask_view(nviewy(nply)) = 2
-                   iflag=1
+                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
-                   return
+                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
+                   goto 10
                 endif
+ *     ------------------> COUPLE <------------------
                 ncp_plane(nplx) = ncp_plane(nplx) + 1
                 clx(nplx,ncp_plane(nplx))=icx
                 cly(nply,ncp_plane(nplx))=icy
                 cl_single(icx)=0
                 cl_single(icy)=0
-             endif
  *     ----------------------------------------------
+             endif
         continue
           enddo                  !end loop on clusters(Y)
      continue
        enddo                     !end loop on clusters(X)
        do icl=1,nclstr1
           if(cl_single(icl).eq.1)then
              ip=npl(VIEW(icl))
-Line 1815 
 c     goto 880   !fill ntp and go to nex
+Line 1795 
 c     goto 880   !fill ntp and go to nex
              cls(ip,ncls(ip))=icl
           endif
        enddo
+  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
-       do ip=1,6
-          ncp_tot=ncp_tot+ncp_plane(ip)
-       enddo
- c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
-       if(ncp_tot.gt.ncp_max)then
-          if(DEBUG)print*,
-      $           '** warning ** number of identified '//
-      $           'couples exceeds upper limit for Hough tr. '
-      $           ,'( ',ncp_max,' )'
- c            good2=.false.
- c     goto 880       !fill ntp and go to next event
-          iflag=1
-          return
-       endif
-       return
-       end
- ***************************************************
- *                                                 *
- *                                                 *
- *                                                 *
- *                                                 *
- *                                                 *
- *                                                 *
- **************************************************
-       subroutine cl_to_couples_nocharge(iflag)
-       include 'commontracker.f'
-       include 'common_momanhough.f'
-       include 'momanhough_init.f'
-       include 'calib.f'
-       include 'level1.f'
-       logical DEBUG
-       common/dbg/DEBUG
- *     output flag
- *     --------------
+       if(.not.RECOVER_SINGLETS)goto 81
- *     0 = good event
- *     1 = bad event
- *     --------------
-       integer iflag
-       integer badseed,badcl
+ *     ////////////////////////////////////////////////
+ *     PATCH to recover events with less than 3 couples
+ *     ////////////////////////////////////////////////
+ *     loop over singlet and create "fake" couples
+ *     (with clx=0 or cly=0)
+ *
- *     init variables
+       if(DEBUG.EQ.1)
-       ncp_tot=0
+      $     print*,'>>> Recover singlets '
-       do ip=1,nplanes
+      $     ,'(creates fake couples) <<<'
-          do ico=1,ncouplemax
+       do icl=1,nclstr1
-             clx(ip,ico)=0
+          if(
-             cly(ip,ico)=0
+      $        cl_single(icl).eq.1.and.
-          enddo
+      $        cl_used(icl).eq.0.and.
-          ncp_plane(ip)=0
+      $        .true.)then
-          do icl=1,nclstrmax_level2
-             cls(ip,icl)=1
-          enddo
-          ncls(ip)=0
-       enddo
-       do icl=1,nclstrmax_level2
-          cl_single(icl)=1
-          cl_good(icl)=0
-       enddo
- *     start association
-       ncouples=0
-       do icx=1,nclstr1          !loop on cluster (X)
-          if(mod(VIEW(icx),2).eq.1)goto 10
  *     ----------------------------------------------------
- *     cut on charge (X VIEW)
+ *     jump masked views
-          if(dedx(icx).lt.dedx_x_min)then
-             cl_single(icx)=0
-             goto 10
-          endif
- *     cut BAD (X VIEW)
-          badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
-          ifirst=INDSTART(icx)
-          if(icx.ne.nclstr1) then
-             ilast=INDSTART(icx+1)-1
-          else
-             ilast=TOTCLLENGTH
-          endif
-          badcl=badseed
-          do igood=-ngoodstr,ngoodstr
-             ibad=1
-             if((INDMAX(icx)+igood).gt.ifirst.and.
-      $           (INDMAX(icx)+igood).lt.ilast.and.
-      $           .true.)then
-                ibad=BAD(VIEW(icx),
-      $              nvk(MAXS(icx)+igood),
-      $              nst(MAXS(icx)+igood))
-             endif
-             badcl=badcl*ibad
-          enddo
-          if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut
-             cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut
-             goto 10             !<<<<<<<<<<<<<< BAD cut
-          endif                  !<<<<<<<<<<<<<< BAD cut
  *     ----------------------------------------------------
+             if( mask_view(VIEW(icl)).ne.0 ) goto 21
-          cl_good(icx)=1
+             if(mod(VIEW(icl),2).eq.0)then !=== X views
-          nplx=npl(VIEW(icx))
-          nldx=nld(MAXS(icx),VIEW(icx))
-          do icy=1,nclstr1       !loop on cluster (Y)
-             if(mod(VIEW(icy),2).eq.0)goto 20
  *     ----------------------------------------------------
- *     cut on charge (Y VIEW)
+ *     cut on charge (X VIEW)
-             if(dedx(icy).lt.dedx_y_min)then
-                cl_single(icy)=0
-                goto 20
-             endif
- *     cut BAD (Y VIEW)
-             badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
-             ifirst=INDSTART(icy)
-             if(icy.ne.nclstr1) then
-                ilast=INDSTART(icy+1)-1
-             else
-                ilast=TOTCLLENGTH
-             endif
-             badcl=badseed
-             do igood=-ngoodstr,ngoodstr
-                ibad=1
-                if((INDMAX(icy)+igood).gt.ifirst.and.
-      $              (INDMAX(icy)+igood).lt.ilast.and.
-      $              .true.)
-      $              ibad=BAD(VIEW(icy),
-      $              nvk(MAXS(icy)+igood),
-      $              nst(MAXS(icy)+igood))
-                badcl=badcl*ibad
-             enddo
-             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut
-                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut
-                goto 20          !<<<<<<<<<<<<<< BAD cut
-             endif               !<<<<<<<<<<<<<< BAD cut
  *     ----------------------------------------------------
+                if(sgnl(icl).lt.dedx_x_min) goto 21
-             cl_good(icy)=1
-             nply=npl(VIEW(icy))
-             nldy=nld(MAXS(icy),VIEW(icy))
- *     ----------------------------------------------
- *     CONDITION TO FORM A COUPLE
- *     ----------------------------------------------
- *     geometrical consistency (same plane and ladder)
-             if(nply.eq.nplx.and.nldy.eq.nldx)then
- *     charge correlation
- *     ===========================================================
- *     this version of the subroutine is used for the calibration
- *     thus charge-correlation selection is obviously removed
- *     ===========================================================
- c$$$               ddd=(dedx(icy)
- c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
- c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
- c$$$               cut=chcut*sch(nplx,nldx)
- c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent
- *     ===========================================================
+                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!!!
+ *     ----------------------------------------------------
- *     ------------------> COUPLE <------------------
+             else                !=== Y views
- *     check to do not overflow vector dimentions
+ *     ----------------------------------------------------
-                if(ncp_plane(nplx).gt.ncouplemax)then
+ *     cut on charge (Y VIEW)
-                   if(DEBUG)print*,
+ *     ----------------------------------------------------
-      $                    ' ** warning ** number of identified'//
+                if(sgnl(icl).lt.dedx_y_min) goto 21
-      $                    ' couples on plane ',nplx,
-      $                    ' exceeds vector dimention'//
-      $                    ' ( ',ncouplemax,' )'
- c     good2=.false.
- c     goto 880   !fill ntp and go to next event
-                   iflag=1
-                   return
-                endif
-                if(ncp_plane(nplx).eq.ncouplemax)then
+                nply=npl(VIEW(icl))
-                   if(DEBUG)print*,
+ *     ------------------> (FAKE) COUPLE <-----------------
-      $                 '** warning ** number of identified '//
+                ncp_plane(nply) = ncp_plane(nply) + 1
-      $                 'couples on plane ',nplx,
+                clx(nply,ncp_plane(nply))=0
-      $                 'exceeds vector dimention '
+                cly(nply,ncp_plane(nply))=icl
-      $                 ,'( ',ncouplemax,' )'
+ c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
- c     good2=.false.
+ *     ----------------------------------------------------
- c     goto 880   !fill ntp and go to next event
-                   iflag=1
-                   return
-                endif
-                ncp_plane(nplx) = ncp_plane(nplx) + 1
+             endif
-                clx(nplx,ncp_plane(nplx))=icx
+          endif                  !end "single" condition
-                cly(nply,ncp_plane(nplx))=icy
+     continue
-                cl_single(icx)=0
+       enddo                     !end loop over clusters
-                cl_single(icy)=0
-             endif
- *     ----------------------------------------------
-        continue
+       if(DEBUG.EQ.1)
-          enddo                  !end loop on clusters(Y)
+      $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
-     continue
-       enddo                     !end loop on clusters(X)
+  continue
-       do icl=1,nclstr1
-          if(cl_single(icl).eq.1)then
-             ip=npl(VIEW(icl))
-             ncls(ip)=ncls(ip)+1
-             cls(ip,ncls(ip))=icl
-          endif
-       enddo
-       if(DEBUG)then
-          print*,'clusters  ',nclstr1
-          print*,'good    ',(cl_good(i),i=1,nclstr1)
-          print*,'singles ',(cl_single(i),i=1,nclstr1)
-          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
-       endif
-       do ip=1,6
+       ncp_tot=0
-          ncp_tot=ncp_tot+ncp_plane(ip)
+       do ip=1,NPLANES
+          ncp_tot = ncp_tot + ncp_plane(ip)
        enddo
- c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
+       if(DEBUG.EQ.1)
+      $     print*,'n.couple tot:      ',ncp_tot
-       if(ncp_tot.gt.ncp_max)then
-          if(DEBUG)print*,
-      $           '** warning ** number of identified '//
-      $           'couples exceeds upper limit for Hough tr. '
-      $           ,'( ',ncp_max,' )'
- c            good2=.false.
- c     goto 880       !fill ntp and go to next event
-          iflag=1
-          return
-       endif
        return
        end
- c$$$      subroutine cl_to_couples_2(iflag)
- c$$$
- c$$$      include 'commontracker.f'
- c$$$      include 'common_momanhough.f'
- c$$$      include 'momanhough_init.f'
- c$$$      include 'calib.f'
- c$$$      include 'level1.f'
- c$$$
- c$$$      logical DEBUG
- c$$$      common/dbg/DEBUG
- c$$$
- c$$$*     output flag
- c$$$*     --------------
- c$$$*     0 = good event
- c$$$*     1 = bad event
- c$$$*     --------------
- c$$$      integer iflag
- c$$$
- c$$$      integer badseed,badcl
- c$$$
- c$$$*     init variables
- c$$$      ncp_tot=0
- c$$$      do ip=1,nplanes
- c$$$         do ico=1,ncouplemax
- c$$$            clx(ip,ico)=0
- c$$$            cly(ip,ico)=0
- c$$$         enddo
- c$$$         ncp_plane(ip)=0
- c$$$         do icl=1,nclstrmax_level2
- c$$$            cls(ip,icl)=1
- c$$$         enddo
- c$$$         ncls(ip)=0
- c$$$      enddo
- c$$$      do icl=1,nclstrmax_level2
- c$$$         cl_single(icl)=1
- c$$$         cl_good(icl)=0
- c$$$      enddo
- c$$$
- c$$$*     start association
- c$$$      ncouples=0
- c$$$      do icx=1,nclstr1          !loop on cluster (X)
- c$$$         if(mod(VIEW(icx),2).eq.1)goto 10
- c$$$
- c$$$*     ----------------------------------------------------
- c$$$*     cut on charge (X VIEW)
- c$$$         if(dedx(icx).lt.dedx_x_min)then
- c$$$            cl_single(icx)=0
- c$$$            goto 10
- c$$$         endif
- c$$$*     cut BAD (X VIEW)
- c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
- c$$$         ifirst=INDSTART(icx)
- c$$$         if(icx.ne.nclstr1) then
- c$$$            ilast=INDSTART(icx+1)-1
- c$$$         else
- c$$$            ilast=TOTCLLENGTH
- c$$$         endif
- c$$$         badcl=badseed
- c$$$         do igood=-ngoodstr,ngoodstr
- c$$$            ibad=1
- c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.
- c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.
- c$$$     $           .true.)then
- c$$$               ibad=BAD(VIEW(icx),
- c$$$     $              nvk(MAXS(icx)+igood),
- c$$$     $              nst(MAXS(icx)+igood))
- c$$$            endif
- c$$$            badcl=badcl*ibad
- c$$$         enddo
- c$$$*         print*,'icx ',icx,badcl
- c$$$         if(badcl.eq.0)then
- c$$$            cl_single(icx)=0
- c$$$            goto 10
- c$$$         endif
- c$$$*     ----------------------------------------------------
- c$$$
- c$$$         cl_good(icx)=1
- c$$$         nplx=npl(VIEW(icx))
- c$$$         nldx=nld(MAXS(icx),VIEW(icx))
- c$$$
- c$$$         do icy=1,nclstr1       !loop on cluster (Y)
- c$$$            if(mod(VIEW(icy),2).eq.0)goto 20
- c$$$
- c$$$*     ----------------------------------------------------
- c$$$*     cut on charge (Y VIEW)
- c$$$            if(dedx(icy).lt.dedx_y_min)then
- c$$$               cl_single(icy)=0
- c$$$               goto 20
- c$$$            endif
- c$$$*     cut BAD (Y VIEW)
- c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
- c$$$            ifirst=INDSTART(icy)
- c$$$            if(icy.ne.nclstr1) then
- c$$$               ilast=INDSTART(icy+1)-1
- c$$$            else
- c$$$               ilast=TOTCLLENGTH
- c$$$            endif
- c$$$            badcl=badseed
- c$$$            do igood=-ngoodstr,ngoodstr
- c$$$               ibad=1
- c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.
- c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.
- c$$$     $              .true.)
- c$$$     $              ibad=BAD(VIEW(icy),
- c$$$     $              nvk(MAXS(icy)+igood),
- c$$$     $              nst(MAXS(icy)+igood))
- c$$$               badcl=badcl*ibad
- c$$$            enddo
- c$$$*            print*,'icy ',icy,badcl
- c$$$            if(badcl.eq.0)then
- c$$$               cl_single(icy)=0
- c$$$               goto 20
- c$$$            endif
- c$$$*     ----------------------------------------------------
- c$$$
- c$$$
- c$$$            cl_good(icy)=1
- c$$$            nply=npl(VIEW(icy))
- c$$$            nldy=nld(MAXS(icy),VIEW(icy))
- c$$$
- c$$$*     ----------------------------------------------
- c$$$*     CONDITION TO FORM A COUPLE
- c$$$*     ----------------------------------------------
- c$$$*     geometrical consistency (same plane and ladder)
- c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then
- c$$$
- c$$$c$$$*     charge correlation
- c$$$c$$$               ddd=(dedx(icy)
- c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
- c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
- c$$$c$$$               cut=chcut*sch(nplx,nldx)
- c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent
- c$$$
- c$$$*     ------------------> COUPLE <------------------
- c$$$*     check to do not overflow vector dimentions
- c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
- c$$$                  if(DEBUG)print*,
- c$$$     $                    ' ** warning ** number of identified'//
- c$$$     $                    ' couples on plane ',nplx,
- c$$$     $                    ' exceeds vector dimention'//
- c$$$     $                    ' ( ',ncouplemax,' )'
- c$$$c     good2=.false.
- c$$$c     goto 880   !fill ntp and go to next event
- c$$$                  iflag=1
- c$$$                  return
- c$$$               endif
- c$$$
- c$$$               if(ncp_plane(nplx).eq.ncouplemax)then
- c$$$                  if(DEBUG)print*,
- c$$$     $                 '** warning ** number of identified '//
- c$$$     $                 'couples on plane ',nplx,
- c$$$     $                 'exceeds vector dimention '
- c$$$     $                 ,'( ',ncouplemax,' )'
- c$$$c     good2=.false.
- c$$$c     goto 880   !fill ntp and go to next event
- c$$$                  iflag=1
- c$$$                  return
- c$$$               endif
- c$$$
- c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1
- c$$$               clx(nplx,ncp_plane(nplx))=icx
- c$$$               cly(nply,ncp_plane(nplx))=icy
- c$$$               cl_single(icx)=0
- c$$$               cl_single(icy)=0
- c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy
- c$$$            endif
- c$$$*     ----------------------------------------------
- c$$$
- c$$$ 20         continue
- c$$$         enddo                  !end loop on clusters(Y)
- c$$$
- c$$$ 10      continue
- c$$$      enddo                     !end loop on clusters(X)
- c$$$
- c$$$
- c$$$      do icl=1,nclstr1
- c$$$         if(cl_single(icl).eq.1)then
- c$$$            ip=npl(VIEW(icl))
- c$$$            ncls(ip)=ncls(ip)+1
- c$$$            cls(ip,ncls(ip))=icl
- c$$$         endif
- c$$$      enddo
- c$$$
- c$$$
- c$$$      if(DEBUG)then
- c$$$         print*,'clusters  ',nclstr1
- c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)
- c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)
- c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
- c$$$      endif
- c$$$
- c$$$      do ip=1,6
- c$$$         ncp_tot=ncp_tot+ncp_plane(ip)
- c$$$      enddo
- c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
- c$$$
- c$$$      if(ncp_tot.gt.ncp_max)then
- c$$$         if(DEBUG)print*,
- c$$$     $           '** warning ** number of identified '//
- c$$$     $           'couples exceeds upper limit for Hough tr. '
- c$$$     $           ,'( ',ncp_max,' )'
- c$$$c            good2=.false.
- c$$$c     goto 880       !fill ntp and go to next event
- c$$$         iflag=1
- c$$$         return
- c$$$      endif
- c$$$
- c$$$      return
- c$$$      end
  ***************************************************
  *                                                 *
-Line 2283 
 c$$$      end
+Line 1888 
 c$$$      end
  **************************************************
        subroutine cp_to_doubtrip(iflag)
- c*****************************************************
- c     02/02/2006 modified by Elena Vannuccini --> (1)
- c*****************************************************
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
-       include 'momanhough_init.f'
        include 'common_xyzPAM.f'
        include 'common_mini_2.f'
        include 'calib.f'
-       include 'level1.f'
-       logical DEBUG
-       common/dbg/DEBUG
  *     output flag
  *     --------------
-Line 2328 
 c      double precision xm3,ym3,zm3
+Line 1927 
 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
+ *     per plane, in order to apply hough transform
+ *     (couples recovered during track refinement)
+ *     --------------------------------------------
+       do ip=1,nplanes
+          if(ncp_plane(ip).gt.ncouplelimit)then
+             mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
+             mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
+          endif
+       enddo
        ndblt=0                   !number of doublets
        ntrpt=0                   !number of triplets
        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
-          do is1=1,2             !loop on sensors - 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)
-                call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,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
-                   do is2=1,2    !loop on sensors -ndblt 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
+ c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
                          call xyz_PAM
-      $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
+      $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
                          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(DEBUG)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
+                               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 2381 
 c                           goto 880 !fi
+Line 2028 
 c                           goto 880 !fi
                          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(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples
+                         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
                                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
+ c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
                                   call xyz_PAM
-      $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
+      $                                (icx3,icy3,is3,PFAdef,PFAdef
+      $                                ,0.,0.,0.,0.)
                                   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
- *     the circle must intersect the reference plane
+ cc                                 if(iflag.ne.0)goto 29
-                                  if(
+                                  if(iflag_t.ne.0)then
- c     $                                 (xc.le.-1.*xclimit.or.
+ *     if tricircle fails, evaluate a straight line
- c     $                                 xc.ge.xclimit).and.
+                                     if(DEBUG.eq.1)
-      $                                radius**2.ge.(ZINI-zc)**2.and.
+      $                                   print*,'TRICIRCLE failure'
-      $                                iflag.eq.0.and.
+      $                                   ,' >>> straight line'
-      $                                .true.)then
+                                     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(DEBUG)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
+                                        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
                 enddo            !end loop on sensors - COPPIA 2
+           continue
              enddo               !end loop on planes  - COPPIA 2
+        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 2514 
 c     goto 880               !ntp fill
+Line 2246 
 c     goto 880               !ntp fill
        subroutine doub_to_YZcloud(iflag)
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
-       include 'momanhough_init.f'
+ c      include 'momanhough_init.f'
-       logical DEBUG
-       common/dbg/DEBUG
  *     output flag
  *     --------------
-Line 2537 
 c     goto 880               !ntp fill
+Line 2268 
 c     goto 880               !ntp fill
        integer cp_useds1(ncouplemaxtot) ! sensor 1
        integer cp_useds2(ncouplemaxtot) ! sensor 2
+       if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  *     classification of DOUBLETS
-Line 2550 
 c     goto 880               !ntp fill
+Line 2282 
 c     goto 880               !ntp fill
        distance=0
        nclouds_yz=0              !number of clouds
        npt_tot=0
+       nloop=0
+  continue
        do idb1=1,ndblt           !loop (1) on DOUBLETS
           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
- c     print*,'--------------'
- c     print*,'** ',idb1,' **'
           do icp=1,ncp_tot
              cp_useds1(icp)=0    !init
              cp_useds2(icp)=0    !init
-Line 2594 
 ccccc if(db_used(idbref).eq.1)goto 1188
+Line 2325 
 ccccc if(db_used(idbref).eq.1)goto 1188
       $              +((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 2619 
 c     print*,idb1,idb2,distance,' cloud
+Line 2342 
 c     print*,idb1,idb2,distance,' cloud
                    temp1 = temp1 + alfayz1(idb2)
                    temp2 = temp2 + alfayz2(idb2)
- c     print*,'*   idbref,idb2 ',idbref,idb2
                 endif
          continue
-Line 2642 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2364 
 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 2652 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2376 
 c     print*,'*   idbref,idb2 ',idbref,i
           do ip=1,nplanes
              nplused=nplused+ hit_plane(ip)
           enddo
- c     print*,'>>>> ',ncpused,npt,nplused
-          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(DEBUG)print*,
+             if(verbose.eq.1)print*,
       $           '** warning ** number of identified '//
       $           'YZ clouds exceeds vector dimention '
       $           ,'( ',ncloyz_max,' )'
  c               good2=.false.
  c     goto 880         !fill ntp and go to next event
+             do iv=1,nviews
+ c               mask_view(iv) = 5
+                mask_view(iv) = mask_view(iv) + 2**4
+             enddo
              iflag=1
              return
           endif
-Line 2682 
 c     goto 880         !fill ntp and go
+Line 2408 
 c     goto 880         !fill ntp and go
  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 2704 
 c$$$     $           ,(db_cloud(iii),iii
+Line 2426 
 c$$$     $           ,(db_cloud(iii),iii
        enddo                     !end loop (1) on DOUBLETS
-       if(DEBUG)then
+       if(nloop.lt.nstepy)then
-          print*,'---------------------- '
+         cutdistyz = cutdistyz+cutystep
+         nloop     = nloop+1
+         goto 90
+       endif
+       if(DEBUG.EQ.1)then
           print*,'Y-Z total clouds ',nclouds_yz
-          print*,' '
        endif
-Line 2730 
 c$$$     $           ,(db_cloud(iii),iii
+Line 2456 
 c$$$     $           ,(db_cloud(iii),iii
        subroutine trip_to_XZcloud(iflag)
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
-       include 'momanhough_init.f'
+ c      include 'momanhough_init.f'
-       logical DEBUG
-       common/dbg/DEBUG
  *     output flag
  *     --------------
-Line 2754 
 c$$$     $           ,(db_cloud(iii),iii
+Line 2479 
 c$$$     $           ,(db_cloud(iii),iii
        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 2765 
 c$$$     $           ,(db_cloud(iii),iii
+Line 2492 
 c$$$     $           ,(db_cloud(iii),iii
        distance=0
        nclouds_xz=0              !number of clouds
        npt_tot=0                 !total number of selected triplets
+       nloop=0
+  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 2802 
 c         tr_temp(1)=itr1
+Line 2529 
 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 2829 
 c     print*,idb1,idb2,distance,' cloud
+Line 2569 
 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
-Line 2850 
 c     print*,'*   itrref,itr2 ',itrref,i
+Line 2589 
 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 2866 
 c     print*,'check cp_used'
+Line 2606 
 c     print*,'check cp_used'
           do ip=1,nplanes
              nplused=nplused+ hit_plane(ip)
           enddo
-          if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
+          if(nplused.lt.nplxz_min)goto 22288 !next triplet
-          if(npt.lt.nptxz_min)goto 22288     !next triplet
-          if(nplused.lt.nplxz_min)goto 22288 !next doublet
  *     ~~~~~~~~~~~~~~~~~
  *     >>> NEW CLOUD <<<
           if(nclouds_xz.ge.ncloxz_max)then
-             if(DEBUG)print*,
+             if(verbose.eq.1)print*,
       $           '** warning ** number of identified '//
       $           'XZ clouds exceeds vector dimention '
       $           ,'( ',ncloxz_max,' )'
  c     good2=.false.
  c     goto 880         !fill ntp and go to next event
+             do iv=1,nviews
+ c               mask_view(iv) = 6
+                mask_view(iv) =  mask_view(iv) + 2**5
+             enddo
              iflag=1
              return
           endif
-Line 2896 
 c     goto 880         !fill ntp and go
+Line 2638 
 c     goto 880         !fill ntp and go
           enddo
           npt_tot=npt_tot+npt
-          if(DEBUG)then
+          if(DEBUG.EQ.1)then
-             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
              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 <<<
  *     ~~~~~~~~~~~~~~~~~
    continue
        enddo                     !end loop (1) on DOUBLETS
-       if(DEBUG)then
+        if(nloop.lt.nstepx)then
-          print*,'---------------------- '
+          cutdistxz=cutdistxz+cutxstep
+          nloop=nloop+1
+          goto 91
+        endif
+       if(DEBUG.EQ.1)then
           print*,'X-Z total clouds ',nclouds_xz
-          print*,' '
        endif
-Line 2936 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2679 
 c$$$     $           ,(tr_cloud(iii),iii
  **************************************************
        subroutine clouds_to_ctrack(iflag)
- c*****************************************************
- c     02/02/2006 modified by Elena Vannuccini --> (1)
- c*****************************************************
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
        include 'common_xyzPAM.f'
        include 'common_mini_2.f'
        include 'common_mech.f'
-       include 'momanhough_init.f'
-       logical DEBUG
-       common/dbg/DEBUG
  *     output flag
  *     --------------
-Line 2964 
 c***************************************
+Line 2703 
 c***************************************
  *     -----------------------------------------------------------
  *     list of matching couples in the combination
  *     between a XZ and YZ cloud
-       integer cp_match(nplanes,ncouplemax)
+       integer cp_match(nplanes,2*ncouplemax)
        integer ncp_match(nplanes)
  *     -----------------------------------------------------------
        integer hit_plane(nplanes)
  *     -----------------------------------------------------------
  *     variables for track fitting
        double precision AL_INI(5)
-       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 2989 
 c      real fitz(nplanes)        !z coor
+Line 2727 
 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 3036 
 c      real fitz(nplanes)        !z coor
+Line 2819 
 c      real fitz(nplanes)        !z coor
              do ip=1,nplanes
                 nplused=nplused+ hit_plane(ip)
              enddo
-             if(nplused.lt.nplxz_min)goto 888 !next doublet
+             if(nplused.lt.3)goto 888 !next combination
-             if(ncp_ok.lt.ncpok)goto 888 !next cloud
+ ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
+ ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
-             if(DEBUG)then
+ *     -----------------------------------------------------------
-                print*,'Combination ',iyz,ixz
+ *     if in RECOVER_SINGLET mode, the two clouds must have
-      $              ,' db ',ptcloud_yz(iyz)
+ *     at least ONE intersecting real couple
-      $              ,' tr ',ptcloud_xz(ixz)
+ *     -----------------------------------------------------------
-      $              ,'  -----> # matching couples ',ncp_ok
+             if(ncp_ok.lt.1)goto 888 !next combination
-             endif
- c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
+             if(DEBUG.EQ.1)then
- c$$$  print*,'Configurazione cluster XZ'
+                print*,'////////////////////////////'
- c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
+                print*,'Cloud combination (Y,X): ',iyz,ixz
- c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))
+                print*,' db ',ptcloud_yz(iyz)
- c$$$  print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))
+                print*,' tr ',ptcloud_xz(ixz)
- c$$$  print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))
+                print*,'  -----> # matching couples ',ncp_ok
- c$$$  print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))
+                print*,'  -----> # fake couples (X)',ncpx_ok
- c$$$  print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))
+                print*,'  -----> # fake couples (Y)',ncpy_ok
- c$$$  print*,'Configurazione cluster YZ'
+                do icp=1,ncp_tot
- c$$$  print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))
+                   print*,'cp ',icp,' >'
- c$$$  print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))
+      $                 ,' x ',cpcloud_xz(ixz,icp)
- c$$$  print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))
+      $                 ,' y ',cpcloud_yz(iyz,icp)
- c$$$  print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))
+      $                 ,' ==> ',cpintersec(icp)
- c$$$  print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))
+                enddo
- c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))
+             endif
- c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
+             if(DEBUG.EQ.1)then
- *     -------> INITIAL GUESS <-------
-             AL_INI(1)=dreal(alfaxz1_av(ixz))
-             AL_INI(2)=dreal(alfayz1_av(iyz))
-             AL_INI(4)=datan(dreal(alfayz2_av(iyz))
-      $           /dreal(alfaxz2_av(ixz)))
-             tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
-             AL_INI(3)=tath/sqrt(1+tath**2)
-             AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
- c     print*,'*******',AL_INI(5)
-             if(AL_INI(5).gt.defmax)goto 888 !next cloud
- c     print*,'alfaxz2, alfayz2 '
- c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)
- *     -------> INITIAL GUESS <-------
- c     print*,'AL_INI ',(al_ini(i),i=1,5)
-             if(DEBUG)then
                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
-Line 3114 
 c     print*,'AL_INI ',(al_ini(i),i=1,5)
+Line 2878 
 c     print*,'AL_INI ',(al_ini(i),i=1,5)
                                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 3131 
 c     print*,'AL_INI ',(al_ini(i),i=1,5)
+Line 2930 
 c     print*,'AL_INI ',(al_ini(i),i=1,5)
  *                                   *************************
  c                                    call xyz_PAM(icx,icy,is,
  c     $                                   'COG2','COG2',0.,0.)
+ c                                    call xyz_PAM(icx,icy,is, !(1)
+ c     $                                   PFAdef,PFAdef,0.,0.) !(1)
                                      call xyz_PAM(icx,icy,is, !(1)
-      $                                   PFAdef,PFAdef,0.,0.) !(1)
+      $                                   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
  *     **********************************************************
  *     ************************** FIT *** FIT *** FIT *** FIT ***
  *     **********************************************************
+ cccc  scommentare se si usa al_ini della nuvola
+ c$$$                              do i=1,5
+ c$$$                                 AL(i)=AL_INI(i)
+ c$$$                              enddo
+                               call guess()
                                do i=1,5
-                                  AL(i)=AL_INI(i)
+                                  AL_INI(i)=AL(i)
                                enddo
                                ifail=0 !error flag in chi^2 computation
                                jstep=0 !number of  minimization steps
-                               call mini_2(jstep,ifail)
+                               iprint=0
+ c                              if(DEBUG.EQ.1)iprint=1
+                               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 *** '
-      $                              //'(mini_2 in clouds_to_ctrack)'
+      $                              //'(clouds_to_ctrack)'
+                                     print*,'initial guess: '
+                                     print*,'AL_INI(1) = ',AL_INI(1)
+                                     print*,'AL_INI(2) = ',AL_INI(2)
+                                     print*,'AL_INI(3) = ',AL_INI(3)
+                                     print*,'AL_INI(4) = ',AL_INI(4)
+                                     print*,'AL_INI(5) = ',AL_INI(5)
                                   endif
-                                  chi2=-chi2
+ c                                 chi2=-chi2
                                endif
  *     **********************************************************
  *     ************************** FIT *** FIT *** FIT *** FIT ***
  *     **********************************************************
                                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(DEBUG)print*,
+                                  if(verbose.eq.1)print*,
       $                 '** warning ** number of candidate tracks '//
       $                 ' exceeds vector dimension '
       $                ,'( ',NTRACKSMAX,' )'
  c                                 good2=.false.
  c                                 goto 880 !fill ntp and go to next event
+                                  do iv=1,nviews
+ c                                    mask_view(iv) = 7
+                                     mask_view(iv) = mask_view(iv) + 2**6
+                                  enddo
                                   iflag=1
                                   return
                                endif
                                ntracks = ntracks + 1
- c$$$                              ndof=0
+                               do ip=1,nplanes !top to bottom
-                               do ip=1,nplanes
- c$$$                                 ndof=ndof
- c$$$     $                                +int(xgood(ip))
- c$$$     $                                +int(ygood(ip))
                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))
-Line 3205 
 c$$$     $
+Line 3053 
 c$$$     $
                                   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)))
+                                     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
-                                  CLS_STORE(nplanes-ip+1,ntracks)=0
+                                  BX_STORE(ip,ntracks)=0!I dont need it now
+                                  BY_STORE(ip,ntracks)=0!I dont need it now
+                                  CLS_STORE(ip,ntracks)=0
                                   do i=1,5
                                      AL_STORE(i,ntracks)=sngl(AL(i))
                                   enddo
                                enddo
- c$$$  *                             Number of Degree Of Freedom
- c$$$  ndof=ndof-5
- c$$$  *                             reduced chi^2
- c$$$  rchi2=chi2/dble(ndof)
                                RCHI2_STORE(ntracks)=chi2
  *     --------------------------------
-Line 3241 
 c$$$  rchi2=chi2/dble(ndof)
+Line 3105 
 c$$$  rchi2=chi2/dble(ndof)
        if(ntracks.eq.0)then
           iflag=1
-          return
+ cc         return
        endif
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
-          print*,'****** TRACK CANDIDATES ***********'
+         print*,'****** TRACK CANDIDATES *****************'
-          print*,'#         R. chi2        RIG'
+         print*,'#         R. chi2        RIG         ndof'
-          do i=1,ntracks
+         do i=1,ntracks
-             print*,i,' --- ',rchi2_store(i),' --- '
+           ndof=0                !(1)
-      $           ,1./abs(AL_STORE(5,i))
+           do ii=1,nplanes       !(1)
-          enddo
+             ndof=ndof           !(1)
-          print*,'***********************************'
+      $           +int(xgood_store(ii,i)) !(1)
+      $           +int(ygood_store(ii,i)) !(1)
+           enddo                 !(1)
+           print*,i,' --- ',rchi2_store(i),' --- '
+      $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
+         enddo
+         print*,'*****************************************'
        endif
-Line 3270 
 c$$$  rchi2=chi2/dble(ndof)
+Line 3140 
 c$$$  rchi2=chi2/dble(ndof)
        subroutine refine_track(ibest)
- c******************************************************
- cccccc 06/10/2005 modified by elena vannuccini ---> (1)
- cccccc 31/01/2006 modified by elena vannuccini ---> (2)
- c******************************************************
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
        include 'common_xyzPAM.f'
        include 'common_mini_2.f'
        include 'common_mech.f'
-       include 'momanhough_init.f'
-       include 'level1.f'
        include 'calib.f'
-       logical DEBUG
-       common/dbg/DEBUG
  *     flag to chose PFA
        character*10 PFA
        common/FINALPFA/PFA
+       real k(6)
+       DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
+       real xp,yp,zp
+       real xyzp(3),bxyz(3)
+       equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
+       if(DEBUG.EQ.1)print*,'refine_track:'
  *     =================================================
  *     new estimate of positions using ETA algorithm
  *                          and
-Line 3299 
 c***************************************
+Line 3169 
 c***************************************
        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)
+          call gufld(xyzp,bxyz)
+          BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
+          BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
+ c$$$  bxyz(1)=0
+ c$$$         bxyz(2)=0
+ c$$$         bxyz(3)=0
+ *     |||||||||||||||||||||||||||||||||||||||||||||||||
  *     -------------------------------------------------
  *     If the plane has been already included, it just
  *     computes again the coordinates of the x-y couple
  *     using improved PFAs
  *     -------------------------------------------------
-          if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
+ *     |||||||||||||||||||||||||||||||||||||||||||||||||
+ 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 3317 
 c***************************************
+Line 3235 
 c***************************************
       $           ,ip_cp(id),ip
              icx=clx(ip,icp)
              icy=cly(ip,icp)
+ c            call xyz_PAM(icx,icy,is,
+ c     $           PFA,PFA,
+ c     $           AXV_STORE(nplanes-ip+1,ibest),
+ c     $           AYV_STORE(nplanes-ip+1,ibest))
              call xyz_PAM(icx,icy,is,
- c     $           'ETA2','ETA2',
       $           PFA,PFA,
       $           AXV_STORE(nplanes-ip+1,ibest),
-      $           AYV_STORE(nplanes-ip+1,ibest))
+      $           AYV_STORE(nplanes-ip+1,ibest),
- c$$$  call xyz_PAM(icx,icy,is,
+      $           bxyz(1),
- c$$$  $              'COG2','COG2',
+      $           bxyz(2)
- c$$$  $              0.,
+      $           )
- c$$$  $              0.)
-             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.
- c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)
+                xgood(nplanes-ip+1) = 1
-             dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
+                ygood(nplanes-ip+1) = 1
-             dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
+                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
+ *     |||||||||||||||||||||||||||||||||||||||||||||||||
  *     -------------------------------------------------
  *     If the plane has NOT  been already included,
  *     it tries to include a COUPLE or a single cluster
  *     -------------------------------------------------
+ *     |||||||||||||||||||||||||||||||||||||||||||||||||
           else
              xgood(nplanes-ip+1)=0
              ygood(nplanes-ip+1)=0
+             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
-             xP=XV_STORE(nplanes-ip+1,ibest)
-             yP=YV_STORE(nplanes-ip+1,ibest)
-             zP=ZV_STORE(nplanes-ip+1,ibest)
              call whichsensor(ip,xP,yP,nldt,ist)
  *     if the track hit the plane in a dead area, go to the next plane
              if(nldt.eq.0.or.ist.eq.0)goto 133
+             SENSOR_STORE(nplanes-ip+1,IBEST)=ist
+             LADDER_STORE(nplanes-ip+1,IBEST)=nldt
  *     --------------------------------------------------------------
-             if(DEBUG)then
+             if(DEBUG.EQ.1)then
                 print*,
       $              '------ Plane ',ip,' intersected on LADDER ',nldt
       $              ,' SENSOR ',ist
-Line 3368 
 c            dedxtrk(nplanes-ip+1) = (de
+Line 3326 
 c            dedxtrk(nplanes-ip+1) = (de
  *     ===========================================
  *     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 3382 
 c            if(DEBUG)print*,'>>>> try t
+Line 3339 
 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
-      $              cl_used(icx).eq.1.or. !or the X cluster is already used
+ c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
-      $              cl_used(icy).eq.1.or. !or the Y 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
+      $              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,
       $              PFA,PFA,
- c     $              'ETA2','ETA2',
       $              AXV_STORE(nplanes-ip+1,ibest),
-      $              AYV_STORE(nplanes-ip+1,ibest))
+      $              AYV_STORE(nplanes-ip+1,ibest),
+      $              bxyz(1),
+      $              bxyz(2)
+      $              )
                 distance = distance_to(XP,YP)
+ c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
                 id=id_cp(ip,icp,ist)
-                if(DEBUG)print*,'( couple ',id
+                if(DEBUG.EQ.1)
-      $              ,' ) normalized distance ',distance
+      $              print*,'( couple ',id
+      $              ,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm = xPAM
                    ymm = yPAM
-Line 3405 
 c     $              'ETA2','ETA2',
+Line 3369 
 c     $              'ETA2','ETA2',
                    rymm = resyPAM
                    distmin = distance
                    idm = id
- c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)
+                   dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
-                   dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
+                   dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
-                   dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
+                   clincnewc=10*sqrt(rymm**2+rxmm**2
+      $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))
                 endif
          continue
              enddo               !end loop on couples on plane icp
-             if(distmin.le.clinc)then
+             if(distmin.le.clincnewc)then
  *              -----------------------------------
-                xm(nplanes-ip+1) = xmm         !<<<
+                xm(nplanes-ip+1) = xmm !<<<
-                ym(nplanes-ip+1) = ymm         !<<<
+                ym(nplanes-ip+1) = ymm !<<<
-                zm(nplanes-ip+1) = zmm         !<<<
+                zm(nplanes-ip+1) = zmm !<<<
-                xgood(nplanes-ip+1) = 1        !<<<
+                xgood(nplanes-ip+1) = 1 !<<<
-                ygood(nplanes-ip+1) = 1        !<<<
+                ygood(nplanes-ip+1) = 1 !<<<
-                resx(nplanes-ip+1)=rxmm        !<<<
+                resx(nplanes-ip+1)=rxmm !<<<
-                resy(nplanes-ip+1)=rymm        !<<<
+                resy(nplanes-ip+1)=rymm !<<<
- c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)
+                dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
-                dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)
+                dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
-                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)
  *              -----------------------------------
                 CP_STORE(nplanes-ip+1,ibest)=idm
-                if(DEBUG)print*,'%%%% included couple ',idm
+                if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
-      $              ,' (norm.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 3452 
 c            if(DEBUG)print*,'>>>> try t
+Line 3415 
 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
  *----- try cluster x -----------------------------------------------
-                if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
+ c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
+                if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
  *                                              !jump to the Y cluster
+ c               call xyz_PAM(icx,0,ist,
+ c     $              PFA,PFA,
+ c     $              AXV_STORE(nplanes-ip+1,ibest),0.)
                 call xyz_PAM(icx,0,ist,
- c     $              'ETA2','ETA2',
       $              PFA,PFA,
-      $              AXV_STORE(nplanes-ip+1,ibest),0.)
+      $              AXV_STORE(nplanes-ip+1,ibest),0.,
+      $              bxyz(1),
+      $              bxyz(2)
+      $              )
                 distance = distance_to(XP,YP)
- c     if(DEBUG)print*,'normalized distance ',distance
+ c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG)print*,'( cl-X ',icx
+                if(DEBUG.EQ.1)
-      $              ,' in cp ',id,' ) normalized distance ',distance
+      $              print*,'( cl-X ',icx
+      $              ,' in cp ',id,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
                    ymm_A = yPAM_A
-Line 3477 
 c     if(DEBUG)print*,'normalized distan
+Line 3448 
 c     if(DEBUG)print*,'normalized distan
                    rymm = resyPAM
                    distmin = distance
                    iclm = icx
- c                  dedxmm = dedx(icx) !(1)
+ c                  dedxmm = sgnl(icx) !(1)
-                   dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
+                   dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
                    dedxmmy = 0.        !(1)
                 endif
          continue
  *----- try cluster y -----------------------------------------------
-                if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
+ c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
+                if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
  *                                              !jump to the next couple
+ c               call xyz_PAM(0,icy,ist,
+ c     $              PFA,PFA,
+ c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
                 call xyz_PAM(0,icy,ist,
- c     $              'ETA2','ETA2',
       $              PFA,PFA,
-      $              0.,AYV_STORE(nplanes-ip+1,ibest))
+      $              0.,AYV_STORE(nplanes-ip+1,ibest),
+      $              bxyz(1),
+      $              bxyz(2)
+      $              )
                 distance = distance_to(XP,YP)
-                if(DEBUG)print*,'( cl-Y ',icy
+ c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-      $              ,' in cp ',id,' ) normalized distance ',distance
+                if(DEBUG.EQ.1)
+      $              print*,'( cl-Y ',icy
+      $              ,' in cp ',id,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
                    ymm_A = yPAM_A
-Line 3503 
 c     $              'ETA2','ETA2',
+Line 3482 
 c     $              'ETA2','ETA2',
                    rymm = resyPAM
                    distmin = distance
                    iclm = icy
- c                 dedxmm = dedx(icy)  !(1)
+ c                 dedxmm = sgnl(icy)  !(1)
                    dedxmmx = 0.        !(1)
-                   dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
+                   dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
                 endif
          continue
              enddo               !end loop on cluster inside couples
  *----- single clusters -----------------------------------------------
              do ic=1,ncls(ip)    !loop on single clusters
                 icl=cls(ip,ic)
-                if(cl_used(icl).eq.1.or.     !if the cluster is already used
+                if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
       $              .false.)goto 18882      !jump to the next singlet
                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
                    call xyz_PAM(icl,0,ist,
- c     $                 'ETA2','ETA2',
       $                 PFA,PFA,
-      $                 AXV_STORE(nplanes-ip+1,ibest),0.)
+      $                 AXV_STORE(nplanes-ip+1,ibest),0.,
+      $                 bxyz(1),
+      $                 bxyz(2)
+      $                 )
                 else                         !<---- Y view
                    call xyz_PAM(0,icl,ist,
- c     $                 'ETA2','ETA2',
       $                 PFA,PFA,
-      $                 0.,AYV_STORE(nplanes-ip+1,ibest))
+      $                 0.,AYV_STORE(nplanes-ip+1,ibest),
+      $                 bxyz(1),
+      $                 bxyz(2)
+      $                 )
                 endif
                 distance = distance_to(XP,YP)
-                if(DEBUG)print*,'( cl-s ',icl
+ c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-      $              ,' ) normalized distance ',distance
+                if(DEBUG.EQ.1)
+      $              print*,'( cl-s ',icl
+      $              ,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
                    ymm_A = yPAM_A
-Line 3541 
 c     $                 'ETA2','ETA2',
+Line 3526 
 c     $                 'ETA2','ETA2',
                    rymm = resyPAM
                    distmin = distance
                    iclm = icl
- c                  dedxmm = dedx(icl)                   !(1)
                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
-                      dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
+                      dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
-                      dedxmmy = 0.                       !(1)
+                      dedxmmy = 0.
                    else          !<---- Y view
-                      dedxmmx = 0.                       !(1)
+                      dedxmmx = 0.
-                      dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
+                      dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
                    endif
                 endif
          continue
              enddo               !end loop on single clusters
-             if(distmin.le.clinc)then
+             if(iclm.ne.0)then
-                CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
- *              ----------------------------
                 if(mod(VIEW(iclm),2).eq.0)then
-                   XGOOD(nplanes-ip+1)=1.
+                   clincnew=
-                   resx(nplanes-ip+1)=rxmm
+      $                 20*
-                   if(DEBUG)print*,'%%%% included X-cl ',iclm
+      $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
-      $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'
+                else if(mod(VIEW(iclm),2).ne.0)then
-                else
+                   clincnew=
-                   YGOOD(nplanes-ip+1)=1.
+      $                 10*
-                   resy(nplanes-ip+1)=rymm
+      $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
-                   if(DEBUG)print*,'%%%% included Y-cl ',iclm
+                endif
-      $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'
+                if(distmin.le.clincnew)then
+                   CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
+ *     ----------------------------
+                   if(mod(VIEW(iclm),2).eq.0)then
+                      XGOOD(nplanes-ip+1)=1.
+                      resx(nplanes-ip+1)=rxmm
+                      if(DEBUG.EQ.1)
+      $                    print*,'%%%% included X-cl ',iclm
+      $                    ,'( chi^2, ',RCHI2_STORE(ibest)
+      $                    ,', dist.= ',distmin
+      $                    ,', cut ',clincnew,' )'
+                   else
+                      YGOOD(nplanes-ip+1)=1.
+                      resy(nplanes-ip+1)=rymm
+                      if(DEBUG.EQ.1)
+      $                    print*,'%%%% included Y-cl ',iclm
+      $                    ,'( chi^2, ',RCHI2_STORE(ibest)
+      $                    ,', dist.= ', distmin
+      $                    ,', cut ',clincnew,' )'
+                   endif
+ *     ----------------------------
+                   xm_A(nplanes-ip+1) = xmm_A
+                   ym_A(nplanes-ip+1) = ymm_A
+                   xm_B(nplanes-ip+1) = xmm_B
+                   ym_B(nplanes-ip+1) = ymm_B
+                   zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
+                   dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
+                   dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
+ *     ----------------------------
                 endif
- *              ----------------------------
-                xm_A(nplanes-ip+1) = xmm_A
-                ym_A(nplanes-ip+1) = ymm_A
-                xm_B(nplanes-ip+1) = xmm_B
-                ym_B(nplanes-ip+1) = ymm_B
-                zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
- c              dedxtrk(nplanes-ip+1) = dedxmm !<<<    !(1)
-                dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1)
-                dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1)
- *              ----------------------------
              endif
           endif
     continue
-Line 3588 
 c              dedxtrk(nplanes-ip+1) = d
+Line 3589 
 c              dedxtrk(nplanes-ip+1) = d
        return
        end
  ***************************************************
  *                                                 *
  *                                                 *
-Line 3596 
 c              dedxtrk(nplanes-ip+1) = d
+Line 3598 
 c              dedxtrk(nplanes-ip+1) = d
  *                                                 *
  *                                                 *
  **************************************************
+ *
-       subroutine clean_XYclouds(ibest,iflag)
-       include 'commontracker.f'
-       include 'common_momanhough.f'
-       include 'momanhough_init.f'
- c      include 'calib.f'
- c      include 'level1.f'
-       logical DEBUG
-       common/dbg/DEBUG
-       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)=1  !tag used clusters
-                cl_used(icly)=1  !tag used clusters
-             elseif(icl.ne.0)then
-                cl_used(icl)=1   !tag used clusters
-             endif
- c               if(DEBUG)then
- c                  print*,ip,' <<< ',id
- c               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
- c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
- c$$$      real function fbad_cog(ncog,ic)
- c$$$
- c$$$
- c$$$      include 'commontracker.f'
- c$$$      include 'level1.f'
- c$$$      include 'calib.f'
- c$$$
- c$$$*     --> signal of the central strip
- c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
- c$$$
- c$$$*     signal of adjacent strips
- c$$$*     --> left
- c$$$      sl1 = 0                  !left 1
- c$$$      if(
- c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
- c$$$     $     )
- c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
- c$$$
- c$$$      sl2 = 0                  !left 2
- c$$$      if(
- c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
- c$$$     $     )
- c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
- c$$$
- c$$$*     --> right
- c$$$      sr1 = 0                  !right 1
- c$$$      if(
- c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
- c$$$     $     .or.
- c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
- c$$$     $     )
- c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
- c$$$
- c$$$      sr2 = 0                  !right 2
- c$$$      if(
- c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
- c$$$     $     .or.
- c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
- c$$$     $     )
- c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
- c$$$
- c$$$
- c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
- c$$$         f  = 4.
- c$$$         si = 8.4
- c$$$      else                              !X-view
- c$$$         f  = 6.
- c$$$         si = 3.9
- c$$$      endif
- c$$$
- c$$$      fbad_cog = 1.
- c$$$      f0 = 1
- c$$$      f1 = 1
- c$$$      f2 = 1
- c$$$      f3 = 1
- c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
- c$$$
- c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
- c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
- c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
- c$$$
- c$$$         if(ncog.eq.2.and.sl1.ne.0)then
- c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
- c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
- c$$$            fbad_cog = 1.
- c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
- c$$$            fbad_cog = 1.
- c$$$         else
- c$$$            fbad_cog = 1.
- c$$$         endif
- c$$$
- c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
- c$$$
- c$$$
- c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
- c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
- c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
- c$$$
- c$$$         if(ncog.eq.2.and.sr1.ne.0)then
- c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
- c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
- c$$$            fbad_cog = 1.
- c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
- c$$$            fbad_cog = 1.
- c$$$         else
- c$$$            fbad_cog = 1.
- c$$$         endif
- c$$$
- c$$$      endif
- c$$$
- c$$$      fbad_cog = sqrt(fbad_cog)
- c$$$
- c$$$      return
- c$$$      end
- c$$$
-Line 3779 
 c$$$
+Line 3606 
 c$$$
        subroutine init_level2
- c*****************************************************
- c     07/10/2005 modified by elena vannuccini --> (1)
- c*****************************************************
        include 'commontracker.f'
+       include 'level1.f'
        include 'common_momanhough.f'
        include 'level2.f'
-       include 'level1.f'
-       good2 = 0!.false.
- c$$$      nev2 = nev1
- c$$$# ifndef TEST2003
- c$$$c*****************************************************
- c$$$cccccc 11/9/2005 modified by david fedele
- c$$$c      pkt_type = pkt_type1
- c$$$c      pkt_num = pkt_num1
- c$$$c      obt = obt1
- c$$$c      which_calib = which_calib1
- c$$$      swcode = 302
- c$$$
- c$$$      which_calib = which_calib1
- c$$$      pkt_type = pkt_type1
- c$$$      pkt_num = pkt_num1
- c$$$      obt = obt1
- c$$$      cpu_crc = cpu_crc1
- c$$$      do iv=1,12
- c$$$         crc(iv)=crc1(iv)
- c$$$      enddo
- c$$$# endif
- c*****************************************************
+ *     ---------------------------------
+ *     variables initialized from level1
+ *     ---------------------------------
+       do i=1,nviews
+          good2(i)=good1(i)
+          do j=1,nva1_view
+             vkflag(i,j)=1
+             if(cnnev(i,j).le.0)then
+                vkflag(i,j)=cnnev(i,j)
+             endif
+          enddo
+       enddo
+ *     ----------------
+ *     level2 variables
+ *     ----------------
        NTRK = 0
-       do it=1,NTRKMAX!NTRACKSMAX
+       do it=1,NTRKMAX
           IMAGE(IT)=0
           CHI2_nt(IT) = -100000.
-          BdL(IT) = 0.
           do ip=1,nplanes
              XM_nt(IP,IT) = 0
              YM_nt(IP,IT) = 0
              ZM_nt(IP,IT) = 0
              RESX_nt(IP,IT) = 0
              RESY_nt(IP,IT) = 0
+             TAILX_nt(IP,IT) = 0
+             TAILY_nt(IP,IT) = 0
+             XBAD(IP,IT) = 0
+             YBAD(IP,IT) = 0
              XGOOD_nt(IP,IT) = 0
              YGOOD_nt(IP,IT) = 0
- c*****************************************************
+             LS(IP,IT) = 0
- cccccc 11/9/2005 modified by david fedele
              DEDX_X(IP,IT) = 0
              DEDX_Y(IP,IT) = 0
- c******************************************************
+             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 3839 
 c***************************************
+Line 3661 
 c***************************************
              enddo
           enddo
        enddo
- c*****************************************************
- cccccc 11/9/2005 modified by david fedele
        nclsx=0
        nclsy=0
        do ip=1,NSINGMAX
          planex(ip)=0
- c        xs(ip)=0
          xs(1,ip)=0
          xs(2,ip)=0
          sgnlxs(ip)=0
          planey(ip)=0
- c        ys(ip)=0
          ys(1,ip)=0
          ys(2,ip)=0
          sgnlys(ip)=0
+         sxbad(ip)=0
+         sybad(ip)=0
+         multmaxsx(ip)=0
+         multmaxsy(ip)=0
        enddo
- c*******************************************************
        end
-Line 3872 
 c***************************************
+Line 3691 
 c***************************************
  ************************************************************
+       subroutine init_hough
+       include 'commontracker.f'
+       include 'level1.f'
+       include 'common_momanhough.f'
+       include 'common_hough.f'
+       include 'level2.f'
+       ntrpt_nt=0
+       ndblt_nt=0
+       NCLOUDS_XZ_nt=0
+       NCLOUDS_YZ_nt=0
+       do idb=1,ndblt_max_nt
+          db_cloud_nt(idb)=0
+          alfayz1_nt(idb)=0
+          alfayz2_nt(idb)=0
+       enddo
+       do itr=1,ntrpt_max_nt
+          tr_cloud_nt(itr)=0
+          alfaxz1_nt(itr)=0
+          alfaxz2_nt(itr)=0
+          alfaxz3_nt(itr)=0
+       enddo
+       do idb=1,ncloyz_max
+         ptcloud_yz_nt(idb)=0
+         alfayz1_av_nt(idb)=0
+         alfayz2_av_nt(idb)=0
+       enddo
+       do itr=1,ncloxz_max
+         ptcloud_xz_nt(itr)=0
+         alfaxz1_av_nt(itr)=0
+         alfaxz2_av_nt(itr)=0
+         alfaxz3_av_nt(itr)=0
+       enddo
+       ntrpt=0
+       ndblt=0
+       NCLOUDS_XZ=0
+       NCLOUDS_YZ=0
+       do idb=1,ndblt_max
+         db_cloud(idb)=0
+         cpyz1(idb)=0
+         cpyz2(idb)=0
+         alfayz1(idb)=0
+         alfayz2(idb)=0
+       enddo
+       do itr=1,ntrpt_max
+         tr_cloud(itr)=0
+         cpxz1(itr)=0
+         cpxz2(itr)=0
+         cpxz3(itr)=0
+         alfaxz1(itr)=0
+         alfaxz2(itr)=0
+         alfaxz3(itr)=0
+       enddo
+       do idb=1,ncloyz_max
+         ptcloud_yz(idb)=0
+         alfayz1_av(idb)=0
+         alfayz2_av(idb)=0
+         do idbb=1,ncouplemaxtot
+           cpcloud_yz(idb,idbb)=0
+         enddo
+       enddo
+       do itr=1,ncloxz_max
+         ptcloud_xz(itr)=0
+         alfaxz1_av(itr)=0
+         alfaxz2_av(itr)=0
+         alfaxz3_av(itr)=0
+         do itrr=1,ncouplemaxtot
+           cpcloud_xz(itr,itrr)=0
+         enddo
+       enddo
+       end
+ ************************************************************
+ *
+ *
+ *
+ *
+ *
+ *
+ *
+ ************************************************************
        subroutine fill_level2_tracks(ntr)
  *     -------------------------------------------------------
-Line 3882 
 c***************************************
+Line 3785 
 c***************************************
        include 'commontracker.f'
+       include 'level1.f'
+       include 'common_momanhough.f'
        include 'level2.f'
        include 'common_mini_2.f'
-       real sinth,phi,pig        !(4)
+       include 'calib.f'
+       character*10 PFA
+       common/FINALPFA/PFA
+       real sinth,phi,pig
+       integer ssensor,sladder
        pig=acos(-1.)
-       good2=1!.true.
-       chi2_nt(ntr)        = sngl(chi2)
-       phi   = al(4)             !(4)
-       sinth = al(3)             !(4)
+ *     -------------------------------------
-       if(sinth.lt.0)then        !(4)
+       chi2_nt(ntr)        = sngl(chi2)
-          sinth = -sinth         !(4)
+       nstep_nt(ntr)       = nstep
-          phi = phi + pig        !(4)
+ *     -------------------------------------
-       endif                     !(4)
+       phi   = al(4)
-       npig = aint(phi/(2*pig))  !(4)
+       sinth = al(3)
-       phi = phi - npig*2*pig    !(4)
+       if(sinth.lt.0)then
-       if(phi.lt.0)              !(4)
+          sinth = -sinth
-      $     phi = phi + 2*pig    !(4)
+          phi = phi + pig
-       al(4) = phi               !(4)
+       endif
-       al(3) = sinth             !(4)
+       npig = aint(phi/(2*pig))
- *****************************************************
+       phi = phi - npig*2*pig
+       if(phi.lt.0)
+      $     phi = phi + 2*pig
+       al(4) = phi
+       al(3) = sinth
        do i=1,5
           al_nt(i,ntr)     = sngl(al(i))
           do j=1,5
              coval(i,j,ntr) = sngl(cov(i,j))
           enddo
- c     print*,al_nt(i,ntr)
        enddo
+ *     -------------------------------------
        do ip=1,nplanes           ! loop on planes
           xgood_nt(ip,ntr) = int(xgood(ip))
           ygood_nt(ip,ntr) = int(ygood(ip))
-Line 3919 
 c     print*,al_nt(i,ntr)
+Line 3831 
 c     print*,al_nt(i,ntr)
           zm_nt(ip,ntr)    = sngl(zm(ip))
           RESX_nt(IP,ntr)  = sngl(resx(ip))
           RESY_nt(IP,ntr)  = sngl(resy(ip))
+          TAILX_nt(IP,ntr) = 0.
+          TAILY_nt(IP,ntr) = 0.
           xv_nt(ip,ntr)    = sngl(xv(ip))
           yv_nt(ip,ntr)    = sngl(yv(ip))
           zv_nt(ip,ntr)    = sngl(zv(ip))
           axv_nt(ip,ntr)   = sngl(axv(ip))
           ayv_nt(ip,ntr)   = sngl(ayv(ip))
- c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)
-          dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
+          factor = sqrt(
-          dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)
+      $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
-       enddo
+      $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
- c      call CalcBdL(100,xxxx,IFAIL)
+      $        1. )
- c      if(ifps(xxxx).eq.1)BdL(ntr) = xxxx
- c$$$      print*,'xgood(ip,ntr) ',(xgood_nt(ip,ntr),ip=1,6)
+          dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
- c$$$      print*,'ygood(ip,ntr) ',(ygood_nt(ip,ntr),ip=1,6)
+          dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)
- c$$$      print*,'dedx_x(ip,ntr) ',(dedx_x(ip,ntr),ip=1,6)
- c$$$      print*,'dedx_y(ip,ntr) ',(dedx_y(ip,ntr),ip=1,6)
+ 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)
+ c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
+ c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
+ c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
+ c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
+ c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001
+ c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
+          angx = effectiveangle(ax,2*ip,bfy)
+          angy = effectiveangle(ay,2*ip-1,bfx)
+          id  = CP_STORE(ip,IDCAND) ! couple id
+          icl = CLS_STORE(ip,IDCAND)
+          ssensor = -1
+          sladder = -1
+          ssensor = SENSOR_STORE(ip,IDCAND)
+          sladder = LADDER_STORE(ip,IDCAND)
+          if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
+          LS(IP,ntr)      = ssensor+10*sladder
+          if(id.ne.0)then
+ c           >>> is a couple
+             cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
+             cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
+             if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
+                cl_used(cltrx(ip,ntr)) = 1 !tag used clusters
+                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
+                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
+      $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
+                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
+                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
+                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
+          print*,'> STORING TRACK ',ntr
+          print*,'clusters: '
+          do ip=1,6
+             print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
+          enddo
+          print*,'dedx: '
+          do ip=1,6
+             print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
+          enddo
+       endif
        end
        subroutine fill_level2_siglets
- c*****************************************************
- c     07/10/2005 created by elena vannuccini
- c     31/01/2006 modified by elena vannuccini
- *     to convert adc to mip  --> (2)
- c*****************************************************
  *     -------------------------------------------------------
  *     This routine fills the  elements of the variables
-Line 3952 
 c***************************************
+Line 3973 
 c***************************************
  *     -------------------------------------------------------
        include 'commontracker.f'
-       include 'level1.f'
-       include 'level2.f'
        include 'calib.f'
+       include 'level1.f'
        include 'common_momanhough.f'
+       include 'level2.f'
        include 'common_xyzPAM.f'
  *     count #cluster per plane not associated to any track
-       good2=1!.true.
        nclsx = 0
        nclsy = 0
+       do iv = 1,nviews
+ c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
+          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) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
+                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.)
-                   call xyz_PAM(icl,0,is,PFAdef,' ',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) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
+                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.)
-                   call xyz_PAM(0,icl,is,' ',PFAdef,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
- c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
+ ***** 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
+ ***************************************************
+ *                                                 *
+ *                                                 *
+ *                                                 *
+ *                                                 *
+ *                                                 *
+ *                                                 *
+ **************************************************
+       subroutine fill_hough
+ *     -------------------------------------------------------
+ *     This routine fills the  variables related to the hough
+ *     transform, for the debig n-tuple
+ *     -------------------------------------------------------
+       include 'commontracker.f'
+       include 'level1.f'
+       include 'common_momanhough.f'
+       include 'common_hough.f'
+       include 'level2.f'
+       if(.false.
+      $     .or.ntrpt.gt.ntrpt_max_nt
+      $     .or.ndblt.gt.ndblt_max_nt
+      $     .or.NCLOUDS_XZ.gt.ncloxz_max
+      $     .or.NCLOUDS_yZ.gt.ncloyz_max
+      $     )then
+          ntrpt_nt=0
+          ndblt_nt=0
+          NCLOUDS_XZ_nt=0
+          NCLOUDS_YZ_nt=0
+       else
+          ndblt_nt=ndblt
+          ntrpt_nt=ntrpt
+          if(ndblt.ne.0)then
+             do id=1,ndblt
+                alfayz1_nt(id)=alfayz1(id) !Y0
+                alfayz2_nt(id)=alfayz2(id) !tg theta-yz
+             enddo
+          endif
+          if(ndblt.ne.0)then
+             do it=1,ntrpt
+                alfaxz1_nt(it)=alfaxz1(it) !X0
+                alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
+                alfaxz3_nt(it)=alfaxz3(it) !1/r
+             enddo
+          endif
+          nclouds_yz_nt=nclouds_yz
+          nclouds_xz_nt=nclouds_xz
+          if(nclouds_yz.ne.0)then
+             nnn=0
+             do iyz=1,nclouds_yz
+                ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
+                alfayz1_av_nt(iyz)=alfayz1_av(iyz)
+                alfayz2_av_nt(iyz)=alfayz2_av(iyz)
+                nnn=nnn+ptcloud_yz(iyz)
+             enddo
+             do ipt=1,nnn
+                db_cloud_nt(ipt)=db_cloud(ipt)
+              enddo
+          endif
+          if(nclouds_xz.ne.0)then
+             nnn=0
+             do ixz=1,nclouds_xz
+                ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
+                alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
+                alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
+                alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
+                nnn=nnn+ptcloud_xz(ixz)
+             enddo
+             do ipt=1,nnn
+               tr_cloud_nt(ipt)=tr_cloud(ipt)
+              enddo
+          endif
+       endif
+       end

 Legend:



Removed from v.1.1
 


changed lines


 
Added in v.1.38
 Legend:



Removed from v.1.1
 


changed lines


 
Added in v.1.38
-Removed from v.1.1
+Added in v.1.38

	ViewVC Help
Powered by ViewVC 1.1.23