/[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.16 by pam-fi,
Thu Nov 30 17:04:27 2006 UTC
+revision 1.46 by mocchiut,
Fri Jan 17 10:14:38 2014 UTC
 Line 20
        include 'calib.f'
        include 'level2.f'
- c      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
  *-------------------------------------------------------------------------------
  *     STEP 1
-Line 41 
 c      include 'momanhough_init.f'
+Line 72 
 c      include 'momanhough_init.f'
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
- c      iflag=0
        call cl_to_couples(iflag)
        if(iflag.eq.1)then        !bad event
           goto 880               !go to next event
        endif
+       if(ncp_tot.eq.0)goto 880  !go to next event
  *-----------------------------------------------------
  *-----------------------------------------------------
  *     HOUGH TRASFORM
-Line 75 
 c      iflag=0
+Line 105 
 c      iflag=0
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
- c      iflag=0
        call cp_to_doubtrip(iflag)
        if(iflag.eq.1)then        !bad event
           goto 880               !go to next event
        endif
+       if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event
  *-------------------------------------------------------------------------------
-Line 130 
 c      iflag=0
+Line 161 
 c      iflag=0
        if(iflag.eq.1)then        !bad event
           goto 880               !fill ntp and go to next event
        endif
+ *     ------------------------------------------------
+ *     first try to release the tolerance
+ *     ------------------------------------------------
        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
-         if(cutdistyz.lt.maxcuty/2)then
+          if(cutdistyz.lt.maxcuty/2)then
-           cutdistyz=cutdistyz+cutystep
+             cutdistyz=cutdistyz+cutystep
-         else
+          else
-           cutdistyz=cutdistyz+(3*cutystep)
+             cutdistyz=cutdistyz+(3*cutystep)
-         endif
+          endif
-         goto 878
+          if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
-       endif
+          goto 878
+       endif
+ *     ------------------------------------------------
+ *     hence reduce the minimum number of plane
+ *     ------------------------------------------------
+       if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then
+          nplyz_min=nplyz_min-1
+          if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min
+          goto 878
+       endif
-       if(planehit.eq.3) goto 881
+       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
- continue
+ c$$$ 881  continue
- *     if there is at least three planes on the Y view decreases cuts on X view
+ c$$$*     if there is at least three planes on the Y view decreases cuts on X view
-       if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
+ c$$$      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
-      $     nplxz_min.ne.y_min_start)then
+ c$$$     $     nplxz_min.ne.y_min_start)then
-         nptxz_min=x_min_step
+ c$$$        nptxz_min=x_min_step
-         nplxz_min=x_min_start-x_min_step
+ c$$$        nplxz_min=x_min_start-x_min_step
-         goto 879
+ c$$$        goto 879
-       endif
+ c$$$      endif
  return
        end
-Line 182 
 c      iflag=0
+Line 241 
 c      iflag=0
  c      include 'momanhough_init.f'
        logical FIMAGE            !
+       real trackimage(NTRACKSMAX)
        real*8 AL_GUESS(5)
  *-------------------------------------------------------------------------------
-Line 209 
 c      include 'momanhough_init.f'
+Line 269 
 c      include 'momanhough_init.f'
  *
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
-          ntrk=0                 !counter of identified physical tracks
+ ccc         ntrk=0                 !counter of identified physical tracks
+ c11111    continue               !<<<<<<< come here when performing a new search
+          continue                  !<<<<<<< come here when performing a new search
-   continue               !<<<<<<< come here when performing a new search
+          if(nclouds_xz.eq.0)goto 880 !go to next event
+          if(nclouds_yz.eq.0)goto 880 !go to next event
  c         iflag=0
           call clouds_to_ctrack(iflag)
           if(iflag.eq.1)then     !no candidate tracks found
              goto 880            !fill ntp and go to next event
           endif
+          if(ntracks.eq.0)goto 880 !go to next event
           FIMAGE=.false.         !processing best track (not track image)
           ibest=0                !best track among candidates
-Line 238 
 c$$$         if(ibest.eq.0)goto 880 !>>
+Line 303 
 c$$$         if(ibest.eq.0)goto 880 !>>
  *     1st) decreasing n.points
  *     2nd) increasing chi**2
  *     -------------------------------------------------------
           rchi2best=1000000000.
-          ndofbest=0             !(1)
+          ndofbest=0
           do i=1,ntracks
-            ndof=0               !(1)
+            ndof=0
-            do ii=1,nplanes      !(1)
+            do ii=1,nplanes
-              ndof=ndof          !(1)
+              ndof=ndof
-      $            +int(xgood_store(ii,i)) !(1)
+      $            +int(xgood_store(ii,i))
-      $            +int(ygood_store(ii,i)) !(1)
+      $            +int(ygood_store(ii,i))
-            enddo                !(1)
+            enddo
-            if(ndof.gt.ndofbest)then !(1)
+            if(ndof.gt.ndofbest)then
               ibest=i
               rchi2best=RCHI2_STORE(i)
-              ndofbest=ndof      !(1)
+              ndofbest=ndof
-            elseif(ndof.eq.ndofbest)then !(1)
+            elseif(ndof.eq.ndofbest)then
               if(RCHI2_STORE(i).lt.rchi2best.and.
       $            RCHI2_STORE(i).gt.0)then
                 ibest=i
                 rchi2best=RCHI2_STORE(i)
-                ndofbest=ndof    !(1)
+                ndofbest=ndof
-              endif              !(1)
+              endif
             endif
           enddo
- c$$$         rchi2best=1000000000.
- c$$$         ndofbest=0             !(1)
- c$$$         do i=1,ntracks
- c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
- c$$$     $          RCHI2_STORE(i).gt.0)then
- c$$$             ndof=0             !(1)
- c$$$             do ii=1,nplanes    !(1)
- c$$$               ndof=ndof        !(1)
- c$$$     $              +int(xgood_store(ii,i)) !(1)
- c$$$     $              +int(ygood_store(ii,i)) !(1)
- c$$$             enddo              !(1)
- c$$$             if(ndof.ge.ndofbest)then !(1)
- c$$$               ibest=i
- c$$$               rchi2best=RCHI2_STORE(i)
- c$$$               ndofbest=ndof    !(1)
- c$$$             endif              !(1)
- c$$$           endif
- c$$$         enddo
           if(ibest.eq.0)goto 880 !>> no good candidates
  *-------------------------------------------------------------------------------
-Line 298 
 c$$$         enddo
+Line 345 
 c$$$         enddo
              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 324 
 c$$$         enddo
+Line 373 
 c$$$         enddo
           jstep=0                !# minimization steps
           iprint=0
- c         if(DEBUG)iprint=1
+ c         if(DEBUG.EQ.1)iprint=1
-          if(VERBOSE)iprint=1
+          if(VERBOSE.EQ.1)iprint=1
-          if(DEBUG)iprint=2
+          if(DEBUG.EQ.1)iprint=2
           call mini2(jstep,ifail,iprint)
-          if(ifail.ne.0) then
+ cc         if(ifail.ne.0) then
-             if(VERBOSE)then
+          if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
+             if(CHI2.ne.CHI2)CHI2=-9999. !new
+             if(VERBOSE.EQ.1)then
                 print *,
       $              '*** MINIMIZATION FAILURE *** (after refinement) '
       $              ,iev
- c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
- c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
- c$$$               print*,'result:   ',(al(i),i=1,5)
- c$$$               print*,'xgood ',xgood
- c$$$               print*,'ygood ',ygood
- c$$$               print*,'----------------------------------------------'
              endif
- c            chi2=-chi2
           endif
-          if(DEBUG)then
+          if(DEBUG.EQ.1)then
              print*,'----------------------------- improved track coord'
       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
              do ip=1,6
-Line 355 
 c            chi2=-chi2
+Line 398 
 c            chi2=-chi2
           endif
  c         rchi2=chi2/dble(ndof)
-          if(DEBUG)then
+          if(DEBUG.EQ.1)then
              print*,' '
              print*,'****** SELECTED TRACK *************'
              print*,'#         R. chi2        RIG'
-Line 371 
 c         rchi2=chi2/dble(ndof)
+Line 414 
 c         rchi2=chi2/dble(ndof)
  *     ------------- search if the track has an IMAGE -------------
  *     ------------- (also this is stored )           -------------
           if(FIMAGE)goto 122     !>>> jump! (this is already an image)
- *     now search for track-image, by comparing couples IDs
+ *     -----------------------------------------------------
+ *     first check if the track is ambiguous
+ *     -----------------------------------------------------
+ *     (modified on august 2007 by ElenaV)
+          is1=0
+          do ip=1,NPLANES
+             if(SENSOR_STORE(ip,icand).ne.0)then
+                is1=SENSOR_STORE(ip,icand)
+                if(ip.eq.6)is1=3-is1 !last plane inverted
+             endif
+          enddo
+          if(is1.eq.0)then
+             if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
+             goto 122            !jump
+          endif
+          do ip=1,NPLANES
+             is2 = SENSOR_STORE(ip,icand) !sensor
+             if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
+             if(
+      $           (is1.ne.is2.and.is2.ne.0)
+      $           )goto 122      !jump (this track cannot have an image)
+          enddo
+          if(DEBUG.eq.1)print*,' >>> ambiguous track! '
+ *     ---------------------------------------------------------------
+ *     take the candidate with the greatest number of matching couples
+ *     if more than one satisfies the requirement,
+ *     choose the one with more points and lower chi2
+ *     ---------------------------------------------------------------
+ *     count the number of matching couples
+          ncpmax = 0
+          iimage   = 0           !id of candidate with better matching
           do i=1,ntracks
-             iimage=i
+             ncp=0
              do ip=1,nplanes
-                if(     CP_STORE(nplanes-ip+1,icand).ne.
+                if(CP_STORE(nplanes-ip+1,icand).ne.0)then
-      $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0
+                   if(
+      $                 CP_STORE(nplanes-ip+1,i).ne.0
+      $                 .and.
+      $                 CP_STORE(nplanes-ip+1,icand).eq.
+      $                 -1*CP_STORE(nplanes-ip+1,i)
+      $                 )then
+                      ncp=ncp+1  !ok
+                   elseif(
+      $                    CP_STORE(nplanes-ip+1,i).ne.0
+      $                    .and.
+      $                    CP_STORE(nplanes-ip+1,icand).ne.
+      $                    -1*CP_STORE(nplanes-ip+1,i)
+      $                    )then
+                      ncp = 0
+                      goto 117   !it is not an image candidate
+                   else
+                   endif
+                endif
+             enddo
+       continue
+             trackimage(i)=ncp   !number of matching couples
+             if(ncp>ncpmax)then
+                ncpmax=ncp
+                iimage=i
+             endif
+          enddo
+ *     check if there are more than one image candidates
+          ngood=0
+          do i=1,ntracks
+             if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
+          enddo
+          if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood
+ *     if there are, choose the best one
+          if(ngood.gt.0)then
+ *     -------------------------------------------------------
+ *     order track-candidates according to:
+ *     1st) decreasing n.points
+ *     2nd) increasing chi**2
+ *     -------------------------------------------------------
+             rchi2best=1000000000.
+             ndofbest=0
+             do i=1,ntracks
+                if( trackimage(i).eq.ncpmax )then
+                   ndof=0
+                   do ii=1,nplanes
+                      ndof=ndof
+      $                    +int(xgood_store(ii,i))
+      $                    +int(ygood_store(ii,i))
+                   enddo
+                   if(ndof.gt.ndofbest)then
+                      iimage=i
+                      rchi2best=RCHI2_STORE(i)
+                      ndofbest=ndof
+                   elseif(ndof.eq.ndofbest)then
+                      if(RCHI2_STORE(i).lt.rchi2best.and.
+      $                    RCHI2_STORE(i).gt.0)then
+                         iimage=i
+                         rchi2best=RCHI2_STORE(i)
+                         ndofbest=ndof
+                      endif
+                   endif
+                endif
              enddo
-             if(  iimage.ne.0.and.
- c     $           RCHI2_STORE(i).le.CHI2MAX.and.
+             if(DEBUG.EQ.1)then
- c     $           RCHI2_STORE(i).gt.0.and.
+                print*,'Track candidate ',iimage
-      $           .true.)then
-                if(DEBUG)print*,'Track candidate ',iimage
       $              ,' >>> TRACK IMAGE >>> of'
       $              ,ibest
-                goto 122         !image track found
              endif
-          enddo
+          endif
     continue
- *     --- and store the results --------------------------------
-          ntrk = ntrk + 1                   !counter of found tracks
-          if(.not.FIMAGE
-      $        .and.iimage.eq.0) image(ntrk)= 0
-          if(.not.FIMAGE
-      $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
-          if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
-          call fill_level2_tracks(ntrk)     !==> good2=.true.
- c         print*,'++++++++++ iimage,fimage,ntrk,image '
- c     $        ,iimage,fimage,ntrk,image(ntrk)
+ *     --- and store the results --------------------------------
           if(ntrk.eq.NTRKMAX)then
-             if(verbose)
+             if(verbose.eq.1)
       $           print*,
       $           '** warning ** number of identified '//
       $           'tracks exceeds vector dimension '
-Line 410 
 c     $        ,iimage,fimage,ntrk,image
+Line 537 
 c     $        ,iimage,fimage,ntrk,image
  cc            good2=.false.
              goto 880            !fill ntp and go to next event
           endif
+          ntrk = ntrk + 1                   !counter of found tracks
+          if(.not.FIMAGE
+      $        .and.iimage.eq.0) image(ntrk)= 0
+          if(.not.FIMAGE
+      $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
+          if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
+          call fill_level2_tracks(ntrk)     !==> good2=.true.
+ c$$$         if(ntrk.eq.NTRKMAX)then
+ c$$$            if(verbose.eq.1)
+ c$$$     $           print*,
+ c$$$     $           '** warning ** number of identified '//
+ c$$$     $           'tracks exceeds vector dimension '
+ c$$$     $           ,'( ',NTRKMAX,' )'
+ c$$$cc            good2=.false.
+ c$$$            goto 880            !fill ntp and go to next event
+ c$$$         endif
           if(iimage.ne.0)then
              FIMAGE=.true.       !
              goto 1212           !>>> fit image-track
           endif
- *     --- then remove selected clusters (ibest+iimage) from clouds ----
-          call clean_XYclouds(ibest,iflag)
-          if(iflag.eq.1)then     !bad event
-             goto 880            !fill ntp and go to next event
-          endif
- *     **********************************************************
- *     condition to start a new search
- *     **********************************************************
-          ixznew=0
-          do ixz=1,nclouds_xz
-             if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1
-          enddo
-          iyznew=0
-          do iyz=1,nclouds_yz
-             if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1
-          enddo
-          if(ixznew.ne.0.and.
-      $      iyznew.ne.0.and.
-      $        rchi2best.le.CHI2MAX.and.
- c     $        rchi2best.lt.15..and.
-      $        .true.)then
-             if(DEBUG)then
-                print*,'***** NEW SEARCH ****'
-             endif
-             goto 11111          !try new search
-          endif
- *     **********************************************
     return
        end
-Line 476 
 c     $        rchi2best.lt.15..and.
+Line 589 
 c     $        rchi2best.lt.15..and.
  *     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 514 
 c     $        rchi2best.lt.15..and.
+Line 629 
 c     $        rchi2best.lt.15..and.
  *
  *
-       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 'level1.f'
        include 'calib.f'
- c      include 'level1.f'
        include 'common_align.f'
        include 'common_mech.f'
        include 'common_xyzPAM.f'
- c      include 'common_resxy.f'
- c      logical DEBUG
- c      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 = pfaeta2(cog2,viewx,nldx,angx)
+      $        .false.)then
- c            stripx = stripx + etacorr
+             print*,'xyz_PAM   ***ERROR*** wrong input '
-             stripx = stripx + pfaeta2(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 + pfaeta3(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 + pfaeta4(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 + pfaeta(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
+ *        --------------------------
+ *        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
+  continue
        endif
- c      if(icy.eq.0.and.icx.ne.0)
- c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'
  *     -----------------
  *     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
+ *        --------------------------
+ *        magnetic-field corrections
+ *        --------------------------
+          stripy  = stripy +  fieldcorr(viewy,bfx)
+          angy    = effectiveangle(ay,viewy,bfx)
-          stripy = float(MAXS(icy))
+          call applypfa(PFAy,icy,angy,corr,res)
-          if(PFAy.eq.'COG1')then !(1)
+          stripy  = stripy + corr
-             stripy = stripy     !(1)
+          resyPAM = res
-             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 = pfaeta2(cog2,viewy,nldy,angy)
- c            stripy = stripy + etacorr
-             stripy = stripy + pfaeta2(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 + pfaeta3(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 + pfaeta4(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 + pfaeta(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
+  continue
        endif
  c===========================================================
  C     COUPLE
  C===========================================================
-Line 697 
 c     (xi,yi,zi) = mechanical coordinate
+Line 777 
 c     (xi,yi,zi) = mechanical coordinate
  c------------------------------------------------------------------------
           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
-             print*,'xyz_PAM (couple):',
+             if(DEBUG.EQ.1) then
-      $          ' WARNING: false X strip: strip ',stripx
+                print*,'xyz_PAM (couple):',
+      $              ' WARNING: false X strip: strip ',stripx
+             endif
           endif
-          xi = acoordsi(stripx,viewx)
+          xi = dcoordsi(stripx,viewx)
-          yi = acoordsi(stripy,viewy)
+          yi = dcoordsi(stripy,viewy)
-          zi = 0.
+          zi = 0.D0
  c------------------------------------------------------------------------
  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-Line 738 
 c---------------------------------------
+Line 819 
 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 764 
 C=======================================
+Line 845 
 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 774 
 C=======================================
+Line 857 
 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 786 
 C=======================================
+Line 867 
 C=======================================
              nldy = nldx
              viewy = viewx - 1
- c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
- c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
-                print*,'xyz_PAM (X-singlet):',
+                if(DEBUG.EQ.1) then
-      $             ' WARNING: false X strip: strip ',stripx
+                   print*,'xyz_PAM (X-singlet):',
+      $                 ' WARNING: false X strip: strip ',stripx
+                endif
              endif
-             xi   = acoordsi(stripx,viewx)
+             xi = dcoordsi(stripx,viewx)
+             yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
+             zi = 0.D0
              xi_A = xi
              yi_A = edgeX_d - SiDimY/2
-Line 809 
 c            if((stripx.le.3).or.(stripx
+Line 893 
 c            if((stripx.le.3).or.(stripx
                 yi_B = yi
              endif
- c            print*,'X-cl ',icx,stripx,' --> ',xi
- c            print*,yi_A,' <--> ',yi_B
           else
+             if(DEBUG.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 855 
 c     N.B. I convert angles from microra
+Line 938 
 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 865 
 c     (xPAM,yPAM,zPAM) = measured coordi
+Line 966 
 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 878 
 c---------------------------------------
+Line 982 
 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 912 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 1188 
 c         print*,'A-(',xPAM_A,yPAM_A,')
  *
  ********************************************************************************
-       real function distance_to(XPP,YPP)
+       real function distance_to(rXPP,rYPP)
        include 'common_xyzPAM.f'
-Line 921 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 1197 
 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 966 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 1247 
 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 991 
 c$$$         print*,' resolution ',re
+Line 1270 
 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 1046 
 c$$$         print*,' resolution ',resxP
+Line 1320 
 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 1071 
 c$$$         print*,' resolution ',resxP
+Line 1345 
 c$$$         print*,' resolution ',resxP
  c------------------------------------------------------------------------
  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-                if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
+                xi = dcoordsi(stripx,viewx)
-      $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
+                yi = dcoordsi(stripy,viewy)
- c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
+                zi = 0.D0
-                   print*,'whichsensor: ',
-      $                ' WARNING: false X strip: strip ',stripx
-                endif
-                xi = acoordsi(stripx,viewx)
-                yi = acoordsi(stripy,viewy)
-                zi = 0.
  c------------------------------------------------------------------------
  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
-Line 1104 
 c---------------------------------------
+Line 1372 
 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 1113 
 c     $              ,iv,xvv(iv),yvv(iv)
+Line 1379 
 c     $              ,iv,xvv(iv),yvv(iv)
                 iv1=iside
                 iv2=mod(iside,4)+1
  *     straight line passing trhough two consecutive vertexes
-                AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))
+                AA = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7
-                BB = yvv(iv1) - AA*xvv(iv1)
+                BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7
  *     point along the straight line closer to the track
                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
                 yoo = AA*xoo + BB
-Line 1126 
 c     $              ,iv,xvv(iv),yvv(iv)
+Line 1392 
 c     $              ,iv,xvv(iv),yvv(iv)
                 iv1=iside
                 iv2=mod(iside,4)+1
  *     straight line passing trhough two consecutive vertexes
-                AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))
+                AA = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7
-                BB = xvv(iv1) - AA*yvv(iv1)
+                BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7
  *     point along the straight line closer to the track
                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
                 xoo = AA*yoo + BB
-Line 1230 
 c      include 'common_analysis.f'
+Line 1496 
 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 1330 
 c      include 'level1.f'
+Line 1595 
 c      include 'level1.f'
        integer badseed,badclx,badcly
+       iflag = iflag
+       if(DEBUG.EQ.1)print*,'cl_to_couples:'
+ cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
  *     init variables
-       ncp_tot=0
        do ip=1,nplanes
           do ico=1,ncouplemax
              clx(ip,ico)=0
-Line 1344 
 c      include 'level1.f'
+Line 1613 
 c      include 'level1.f'
           ncls(ip)=0
        enddo
        do icl=1,nclstrmax_level2
-          cl_single(icl) = 1
+          cl_single(icl) = 1     !all are single
-          cl_good(icl)   = 0
+          cl_good(icl)   = 0     !all are bad
        enddo
        do iv=1,nviews
           ncl_view(iv)  = 0
-Line 1359 
 c      include 'level1.f'
+Line 1628 
 c      include 'level1.f'
  *     mask views with too many clusters
        do iv=1,nviews
           if( ncl_view(iv).gt. nclusterlimit)then
-             mask_view(iv) = 1
+ c            mask_view(iv) = 1
-             if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
+             mask_view(iv) = mask_view(iv) + 2**0
-      $           ,nclusterlimit,' on view ', iv,' --> masked!'
+             if(DEBUG.EQ.1)
+      $        print*,' * WARNING * cl_to_couple: n.clusters > '
+      $        ,nclusterlimit,' on view ', iv,' --> masked!'
           endif
        enddo
-Line 1371 
 c      include 'level1.f'
+Line 1642 
 c      include 'level1.f'
        do icx=1,nclstr1          !loop on cluster (X)
           if(mod(VIEW(icx),2).eq.1)goto 10
+          if(cl_used(icx).ne.0)goto 10
  *     ----------------------------------------------------
  *     jump masked views (X VIEW)
  *     ----------------------------------------------------
-Line 1378 
 c      include 'level1.f'
+Line 1651 
 c      include 'level1.f'
  *     ----------------------------------------------------
  *     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
-Line 1419 
 c     endif
+Line 1692 
 c     endif
           do icy=1,nclstr1       !loop on cluster (Y)
              if(mod(VIEW(icy),2).eq.0)goto 20
+             if(cl_used(icx).ne.0)goto 20
  *     ----------------------------------------------------
  *     jump masked views (Y VIEW)
-Line 1428 
 c     endif
+Line 1703 
 c     endif
  *     ----------------------------------------------------
  *     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
-Line 1474 
 c     endif
+Line 1749 
 c     endif
  *     charge correlation
  *     (modified to be applied only below saturation... obviously)
-                if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
+                if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
       $              .and.
-      $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
+      $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
       $              .and.
       $              (badclx.eq.1.and.badcly.eq.1)
       $              .and.
       $              .true.)then
-                   ddd=(dedx(icy)
+                   ddd=(sgnl(icy)
-      $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
+      $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
  c                  cut = chcut * sch(nplx,nldx)
-                   sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
+                   sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
       $                 -kch(nplx,nldx)*cch(nplx,nldx))
                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
                    cut = chcut * (16 + sss/50.)
-Line 1499 
 c                  cut = chcut * sch(npl
+Line 1774 
 c                  cut = chcut * sch(npl
                 endif
                 if(ncp_plane(nplx).gt.ncouplemax)then
-                   if(verbose)print*,
+                   if(verbose.eq.1)print*,
       $                 '** warning ** number of identified '//
       $                 'couples on plane ',nplx,
       $                 'exceeds vector dimention '
       $                 ,'( ',ncouplemax,' ) --> masked!'
-                   mask_view(nviewx(nplx)) = 2
+ c                  mask_view(nviewx(nplx)) = 2
-                   mask_view(nviewy(nply)) = 2
+ c                  mask_view(nviewy(nply)) = 2
+                  mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
+                  mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
                    goto 10
                 endif
-Line 1525 
 c                  cut = chcut * sch(npl
+Line 1802 
 c                  cut = chcut * sch(npl
      continue
        enddo                     !end loop on clusters(X)
        do icl=1,nclstr1
           if(cl_single(icl).eq.1)then
              ip=npl(VIEW(icl))
-Line 1533 
 c                  cut = chcut * sch(npl
+Line 1809 
 c                  cut = chcut * sch(npl
              cls(ip,ncls(ip))=icl
           endif
        enddo
+ c 80   continue
+       continue
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
           print*,'clusters  ',nclstr1
           print*,'good    ',(cl_good(i),i=1,nclstr1)
-          print*,'singles ',(cl_single(i),i=1,nclstr1)
+          print*,'used    ',(cl_used(i),i=1,nclstr1)
+          print*,'singlets',(cl_single(i),i=1,nclstr1)
           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
        endif
+       if(.not.RECOVER_SINGLETS)goto 81
+ *     ////////////////////////////////////////////////
+ *     PATCH to recover events with less than 3 couples
+ *     ////////////////////////////////////////////////
+ *     loop over singlet and create "fake" couples
+ *     (with clx=0 or cly=0)
+ *
+       if(DEBUG.EQ.1)
+      $     print*,'>>> Recover singlets '
+      $     ,'(creates fake couples) <<<'
+       do icl=1,nclstr1
+          if(
+      $        cl_single(icl).eq.1.and.
+      $        cl_used(icl).eq.0.and.
+      $        .true.)then
+ *     ----------------------------------------------------
+ *     jump masked views
+ *     ----------------------------------------------------
+             if( mask_view(VIEW(icl)).ne.0 ) goto 21
+             if(mod(VIEW(icl),2).eq.0)then !=== X views
+ *     ----------------------------------------------------
+ *     cut on charge (X VIEW)
+ *     ----------------------------------------------------
+                if(sgnl(icl).lt.dedx_x_min) goto 21
+                nplx=npl(VIEW(icl))
+ *     ------------------> (FAKE) COUPLE <-----------------
+                ncp_plane(nplx) = ncp_plane(nplx) + 1
+                clx(nplx,ncp_plane(nplx))=icl
+                cly(nplx,ncp_plane(nplx))=0
+ c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
+ *     ----------------------------------------------------
+             else                !=== Y views
+ *     ----------------------------------------------------
+ *     cut on charge (Y VIEW)
+ *     ----------------------------------------------------
+                if(sgnl(icl).lt.dedx_y_min) goto 21
+                nply=npl(VIEW(icl))
+ *     ------------------> (FAKE) COUPLE <-----------------
+                ncp_plane(nply) = ncp_plane(nply) + 1
+                clx(nply,ncp_plane(nply))=0
+                cly(nply,ncp_plane(nply))=icl
+ c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
+ *     ----------------------------------------------------
+             endif
+          endif                  !end "single" condition
+     continue
+       enddo                     !end loop over clusters
+       if(DEBUG.EQ.1)
+      $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
+  continue
-       do ip=1,6
+       ncp_tot=0
+       do ip=1,NPLANES
           ncp_tot = ncp_tot + ncp_plane(ip)
        enddo
+       if(DEBUG.EQ.1)
+      $     print*,'n.couple tot:      ',ncp_tot
        return
        end
-Line 1598 
 c      double precision xm3,ym3,zm3
+Line 1942 
 c      double precision xm3,ym3,zm3
        real xc,zc,radius
  *     -----------------------------
+       if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
  *     --------------------------------------------
  *     put a limit to the maximum number of couples
-Line 1606 
 c      double precision xm3,ym3,zm3
+Line 1951 
 c      double precision xm3,ym3,zm3
  *     --------------------------------------------
        do ip=1,nplanes
           if(ncp_plane(ip).gt.ncouplelimit)then
-             mask_view(nviewx(ip)) = 8
+             mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
-             mask_view(nviewy(ip)) = 8
+             mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
           endif
        enddo
-Line 1616 
 c      double precision xm3,ym3,zm3
+Line 1961 
 c      double precision xm3,ym3,zm3
        ntrpt=0                   !number of triplets
        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
+ c$$$         print*,'(1) ip ',ip1
+ c$$$     $        ,mask_view(nviewx(ip1))
+ c$$$     $        ,mask_view(nviewy(ip1))
           if(  mask_view(nviewx(ip1)).ne.0 .or.
       $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
           do is1=1,2             !loop on sensors - COPPIA 1
              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
                 icx1=clx(ip1,icp1)
                 icy1=cly(ip1,icp1)
+ c$$$               print*,'(1) ip ',ip1,' icp ',icp1
  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
-                call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
+ c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
-                xm1=xPAM
+                call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
-                ym1=yPAM
+                xm1=REAL(xPAM) !EM GCC4.7
-                zm1=zPAM
+                ym1=REAL(yPAM) !EM GCC4.7
- c     print*,'***',is1,xm1,ym1,zm1
+                zm1=REAL(zPAM) !EM GCC4.7
                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
+ c$$$                  print*,'(2) ip ',ip2
+ c$$$     $                 ,mask_view(nviewx(ip2))
+ c$$$     $                 ,mask_view(nviewy(ip2))
                    if(  mask_view(nviewx(ip2)).ne.0 .or.
       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
                          icx2=clx(ip2,icp2)
                          icy2=cly(ip2,icp2)
+ c$$$                        print*,'(2) ip ',ip2,' icp ',icp2
  c                        call xyz_PAM
  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
+ c                        call xyz_PAM
+ 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
+                         xm2=REAL(xPAM) !EM GCC4.7
-                         ym2=yPAM
+                         ym2=REAL(yPAM) !EM GCC4.7
-                         zm2=zPAM
+                         zm2=REAL(zPAM) !EM GCC4.7
+ *                       ---------------------------------------------------
+ *                       both couples must have a y-cluster
+ *                       (condition necessary when in RECOVER_SINGLETS mode)
+ *                       ---------------------------------------------------
+                         if(icy1.eq.0.or.icy2.eq.0)goto 111
+                         if(cl_used(icy1).ne.0)goto 111
+                         if(cl_used(icy2).ne.0)goto 111
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  *     track parameters on Y VIEW
  *     (2 couples needed)
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                          if(ndblt.eq.ndblt_max)then
-                            if(verbose)print*,
+                            if(verbose.eq.1)print*,
       $                          '** warning ** number of identified '//
       $                          'doublets exceeds vector dimention '
       $                          ,'( ',ndblt_max,' )'
- c                           good2=.false.
+ c     good2=.false.
- c                           goto 880 !fill ntp and go to next event
+ c     goto 880 !fill ntp and go to next event
                             do iv=1,12
-                               mask_view(iv) = 3
+ c     mask_view(iv) = 3
+                               mask_view(iv) = mask_view(iv)+ 2**2
                             enddo
                             iflag=1
                             return
                          endif
+ ccc                        print*,'<doublet> ',icp1,icp2
                          ndblt = ndblt + 1
  *     store doublet info
                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
-Line 1669 
 c                           goto 880 !fi
+Line 2042 
 c                           goto 880 !fi
  *     tg(th_yz)
                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
  *     y0 (cm)
-                         alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
+                     alfayz1(ndblt)=alfayz2(ndblt)*(REAL(zini)-zm1)+ym1! EM GCC4.7 zm1, ym1 and alfayz1/2 are REAL
  ****  -----------------------------------------------****
  ****  reject non phisical couples                    ****
  ****  -----------------------------------------------****
+                         if(SECOND_SEARCH)goto 111
                          if(
       $                       abs(alfayz2(ndblt)).gt.alfyz2_max
       $                       .or.
       $                       abs(alfayz1(ndblt)).gt.alfyz1_max
       $                       )ndblt = ndblt-1
- c$$$      if(iev.eq.33)then
- c$$$      print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2
+                   continue
- c$$$     $        ,' || ',icx1,icy1,icx2,icy2
- c$$$     $        ,' || ',xm1,ym1,xm2,ym2
- c$$$     $        ,' || ',alfayz2(ndblt),alfayz1(ndblt)
- c$$$      endif
- c$$$
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  *     track parameters on Y VIEW - end
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+                         if(icx1.ne.0)then
+                            if(cl_used(icx1).ne.0)goto 31
+                         endif
+                         if(icx2.ne.0)then
+                            if(cl_used(icx2).ne.0)goto 31
+                         endif
                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
+ c$$$                           print*,'(3) ip ',ip3
+ c$$$     $                          ,mask_view(nviewx(ip3))
+ c$$$     $                          ,mask_view(nviewy(ip3))
                             if(  mask_view(nviewx(ip3)).ne.0 .or.
       $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
                             do is3=1,2 !loop on sensors - COPPIA 3
-Line 1702 
 c$$$
+Line 2081 
 c$$$
                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
                                   icx3=clx(ip3,icp3)
                                   icy3=cly(ip3,icp3)
+ c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
+ *     ---------------------------------------------------
+ *     all three couples must have a x-cluster
+ *     (condition necessary when in RECOVER_SINGLETS mode)
+ *     ---------------------------------------------------
+                                  if(
+      $                                icx1.eq.0.or.
+      $                                icx2.eq.0.or.
+      $                                icx3.eq.0.or.
+      $                                .false.)goto 29
+                                  if(cl_used(icx1).ne.0)goto 29
+                                  if(cl_used(icx2).ne.0)goto 29
+                                  if(cl_used(icx3).ne.0)goto 29
  c                                 call xyz_PAM
  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
+ c                                 call xyz_PAM
+ 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
-                                  xm3=xPAM
+      $                                ,0.,0.,0.,0.)
-                                  ym3=yPAM
+                                  xm3=REAL(xPAM) !EM GCC4.7
-                                  zm3=zPAM
+                                  ym3=REAL(yPAM) !EM GCC4.7
+                                  zm3=REAL(zPAM) !EM GCC4.7
  *     find the circle passing through the three points
+                                  iflag_t = DEBUG
                                   call tricircle(3,xp,zp,angp,resp,chi
-      $                                ,xc,zc,radius,iflag)
+      $                                ,xc,zc,radius,iflag_t)
- c     print*,xc,zc,radius
  *     the circle must intersect the reference plane
-                                  if(
+ cc                                 if(iflag.ne.0)goto 29
- c     $                                 (xc.le.-1.*xclimit.or.
+                                  if(iflag_t.ne.0)then
- c     $                                 xc.ge.xclimit).and.
+ *     if tricircle fails, evaluate a straight line
-      $                                radius**2.ge.(ZINI-zc)**2.and.
+                                     if(DEBUG.eq.1)
-      $                                iflag.eq.0.and.
+      $                                   print*,'TRICIRCLE failure'
-      $                                .true.)then
+      $                                   ,' >>> straight line'
+                                     radius=0.
+                                     xc=0.
+                                     yc=0.
+                                     SZZ=0.
+                                     SZX=0.
+                                     SSX=0.
+                                     SZ=0.
+                                     S1=0.
+                                     X0=0.
+                                     Ax=0.
+                                     BX=0.
+                                     DO I=1,3
+                                        XX = XP(I)
+                                        SZZ=SZZ+ZP(I)*ZP(I)
+                                        SZX=SZX+ZP(I)*XX
+                                        SSX=SSX+XX
+                                        SZ=SZ+ZP(I)
+                                        S1=S1+1.
+                                     ENDDO
+                                     DET=SZZ*S1-SZ*SZ
+                                     AX=(SZX*S1-SZ*SSX)/DET
+                                     BX=(SZZ*SSX-SZX*SZ)/DET
+                                     X0  = AX*REAL(ZINI)+BX ! EM GCC4.7
+                                  endif
+                                  if(  .not.SECOND_SEARCH.and.
+      $                                radius**2.lt.(ZINI-zc)**2)goto 29
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  *     track parameters on X VIEW
  *     (3 couples needed)
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                                   if(ntrpt.eq.ntrpt_max)then
-                                     if(verbose)print*,
+                                     if(verbose.eq.1)print*,
-      $                     '** warning ** number of identified '//
+      $                                   '** warning **'//
-      $                     'triplets exceeds vector dimention '
+      $                                   ' number of identified '//
-      $                    ,'( ',ntrpt_max,' )'
+      $                                   'triplets exceeds'//
- c                                    good2=.false.
+      $                                   ' vector dimention '
- c                                    goto 880 !fill ntp and go to next event
+      $                                   ,'( ',ntrpt_max,' )'
+ c     good2=.false.
+ c     goto 880 !fill ntp and go to next event
                                      do iv=1,nviews
-                                        mask_view(iv) = 4
+ 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-(REAL(ZINI)-zc)**2) !EM GCC4.7
-               alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
+            alfaxz2(ntrpt) = (REAL(ZINI)-zc)/
-               alfaxz3(ntrpt) = 1/radius
+      $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
-                                  else
+            alfaxz3(ntrpt) = 1/radius
+                                 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-(REAL(ZINI)-zc)**2)
-               alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
+            alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/
-               alfaxz3(ntrpt) = -1/radius
+      $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
-                                  endif
+            alfaxz3(ntrpt) = -1/radius
+                                 else if(radius.eq.0)then
+ *************straight fit
+              alfaxz1(ntrpt) = X0
+              alfaxz2(ntrpt) = AX
+              alfaxz3(ntrpt) = 0.
+                                 endif
+ c$$$                                print*,'alfaxz1 ', alfaxz1(ntrpt)
+ c$$$                                print*,'alfaxz2 ', alfaxz2(ntrpt)
+ c$$$                                print*,'alfaxz3 ', alfaxz3(ntrpt)
  ****  -----------------------------------------------****
  ****  reject non phisical triplets                   ****
  ****  -----------------------------------------------****
-                                  if(
+                                 if(SECOND_SEARCH)goto 29
-      $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max
+                                 if(
-      $                                .or.
+      $                               abs(alfaxz2(ntrpt)).gt.
-      $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max
+      $                               alfxz2_max
-      $                                )ntrpt = ntrpt-1
+      $                               .or.
+      $                               abs(alfaxz1(ntrpt)).gt.
+      $                               alfxz1_max
- c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)
+      $                               )ntrpt = ntrpt-1
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  *     track parameters on X VIEW - end
  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-                               endif
+                          continue
                             enddo !end loop on COPPIA 3
                          enddo   !end loop on sensors - COPPIA 3
                     continue
                       enddo      !end loop on planes  - COPPIA 3
-                 continue
+                 continue
-               enddo         !end loop on COPPIA 2
+ c 1                enddo         !end loop on COPPIA 2
+                  enddo         !end loop on COPPIA 2
                 enddo            !end loop on sensors - COPPIA 2
            continue
              enddo               !end loop on planes  - COPPIA 2
+ c 11         continue
+           continue
           enddo                  !end loop on COPPIA1
        enddo                     !end loop on sensors - COPPIA 1
   continue
        enddo                     !end loop on planes  - COPPIA 1
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
           print*,'--- doublets ',ndblt
           print*,'--- triplets ',ntrpt
        endif
-Line 1834 
 c      include 'momanhough_init.f'
+Line 2287 
 c      include 'momanhough_init.f'
        integer cp_useds1(ncouplemaxtot) ! sensor 1
        integer cp_useds2(ncouplemaxtot) ! sensor 2
+       if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  *     classification of DOUBLETS
-Line 1852 
 c      include 'momanhough_init.f'
+Line 2306 
 c      include 'momanhough_init.f'
        do idb1=1,ndblt           !loop (1) on DOUBLETS
           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
- c     print*,'--------------'
- c     print*,'** ',idb1,' **'
           do icp=1,ncp_tot
              cp_useds1(icp)=0    !init
              cp_useds2(icp)=0    !init
-Line 1890 
 ccccc if(db_used(idbref).eq.1)goto 1188
+Line 2341 
 ccccc if(db_used(idbref).eq.1)goto 1188
  *     doublet distance in parameter space
                 distance=
       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
       $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
                 distance = sqrt(distance)
- c$$$      if(iev.eq.33)then
- c$$$      if(distance.lt.100)
- c$$$     $ print*,'********* ',idb1,idbref,idb2,distance
- c$$$      if(distance.lt.100)
- c$$$     $ print*,'********* ',alfayz1(idbref),alfayz1(idb2)
- c$$$     $                    ,alfayz2(idbref),alfayz2(idb2)
- c$$$      endif
                 if(distance.lt.cutdistyz)then
- c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
-Line 1918 
 c     print*,idb1,idb2,distance,' cloud
+Line 2361 
 c     print*,idb1,idb2,distance,' cloud
                    temp1 = temp1 + alfayz1(idb2)
                    temp2 = temp2 + alfayz2(idb2)
- c     print*,'*   idbref,idb2 ',idbref,idb2
                 endif
          continue
              enddo               !end loop (2) on DOUBLETS
-      continue
+ c 1188       continue
+             continue
           enddo                  !end loop on... bo?
           nptloop=npv
-Line 1941 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2384 
 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 1951 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2396 
 c     print*,'*   idbref,idb2 ',idbref,i
           do ip=1,nplanes
              nplused=nplused+ hit_plane(ip)
           enddo
- c     print*,'>>>> ',ncpused,npt,nplused
- c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
-          if(npt.lt.nptyz_min)goto 2228 !next doublet
           if(nplused.lt.nplyz_min)goto 2228 !next doublet
  *     ~~~~~~~~~~~~~~~~~
  *     >>> NEW CLOUD <<<
           if(nclouds_yz.ge.ncloyz_max)then
-             if(verbose)print*,
+             if(verbose.eq.1)print*,
       $           '** warning ** number of identified '//
       $           'YZ clouds exceeds vector dimention '
       $           ,'( ',ncloyz_max,' )'
  c               good2=.false.
  c     goto 880         !fill ntp and go to next event
              do iv=1,nviews
-                mask_view(iv) = 5
+ c               mask_view(iv) = 5
+                mask_view(iv) = mask_view(iv) + 2**4
              enddo
              iflag=1
              return
-Line 1984 
 c     goto 880         !fill ntp and go
+Line 2428 
 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 2012 
 c$$$     $           ,(db_cloud(iii),iii
+Line 2452 
 c$$$     $           ,(db_cloud(iii),iii
          goto 90
        endif
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
-          print*,'---------------------- '
           print*,'Y-Z total clouds ',nclouds_yz
-          print*,' '
        endif
-Line 2061 
 c      include 'momanhough_init.f'
+Line 2499 
 c      include 'momanhough_init.f'
        integer cp_useds1(ncouplemaxtot) ! sensor 1
        integer cp_useds2(ncouplemaxtot) ! sensor 2
+       if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  *     classification of TRIPLETS
  *     according to distance in parameter space
-Line 2076 
 c      include 'momanhough_init.f'
+Line 2516 
 c      include 'momanhough_init.f'
   continue
        do itr1=1,ntrpt           !loop (1) on TRIPLETS
           if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
- c     print*,'--------------'
- c     print*,'** ',itr1,' **'
           do icp=1,ncp_tot
              cp_useds1(icp)=0
-Line 2111 
 c         tr_temp(1)=itr1
+Line 2549 
 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 2138 
 c     print*,idb1,idb2,distance,' cloud
+Line 2589 
 c     print*,idb1,idb2,distance,' cloud
                    temp1 = temp1 + alfaxz1(itr2)
                    temp2 = temp2 + alfaxz2(itr2)
                    temp3 = temp3 + alfaxz3(itr2)
- c     print*,'*   itrref,itr2 ',itrref,itr2,distance
                 endif
          continue
              enddo               !end loop (2) on TRIPLETS
-      continue
+ c11888       continue
+             continue
           enddo                  !end loop on... bo?
           nptloop=npv
-Line 2159 
 c     print*,'*   itrref,itr2 ',itrref,i
+Line 2610 
 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 2175 
 c     print*,'check cp_used'
+Line 2627 
 c     print*,'check cp_used'
           do ip=1,nplanes
              nplused=nplused+ hit_plane(ip)
           enddo
- c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
+          if(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(verbose)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
-                mask_view(iv) = 6
+ c               mask_view(iv) = 6
+                mask_view(iv) =  mask_view(iv) + 2**5
              enddo
              iflag=1
              return
-Line 2208 
 c     goto 880         !fill ntp and go
+Line 2659 
 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*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
+             print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
-             print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
+             print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
-             print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
+             print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
-             print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
+             print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
-             print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
+             print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
-             print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
+             print*,'cpcloud_xz '
+      $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
- c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
- c$$$     $           ,ptcloud_xz(nclouds_xz)
- c$$$            print*,'nt-uple: tr_cloud(...) = '
- c$$$     $           ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot)
           endif
  *     >>> NEW CLOUD <<<
  *     ~~~~~~~~~~~~~~~~~
-Line 2233 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2681 
 c$$$     $           ,(tr_cloud(iii),iii
           goto 91
         endif
-       if(DEBUG)then
+       if(DEBUG.EQ.1)then
-          print*,'---------------------- '
           print*,'X-Z total clouds ',nclouds_xz
-          print*,' '
        endif
-Line 2254 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2700 
 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'
-Line 2264 
 c***************************************
+Line 2707 
 c***************************************
        include 'common_xyzPAM.f'
        include 'common_mini_2.f'
        include 'common_mech.f'
- c      include 'momanhough_init.f'
  *     output flag
-Line 2288 
 c      include 'momanhough_init.f'
+Line 2731 
 c      include 'momanhough_init.f'
  *     -----------------------------------------------------------
  *     variables for track fitting
        double precision AL_INI(5)
- c      double precision tath
  *     -----------------------------------------------------------
- c      real fitz(nplanes)        !z coordinates of the planes in cm
+       if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
        ntracks=0                 !counter of track candidates
-       do iyz=1,nclouds_yz       !loop on YZ couds
+       do iyz=1,nclouds_yz       !loop on YZ clouds
-          do ixz=1,nclouds_xz    !loop on XZ couds
+          do ixz=1,nclouds_xz    !loop on XZ clouds
  *     --------------------------------------------------
  *     check of consistency of the clouds
-Line 2306 
 c      real fitz(nplanes)        !z coor
+Line 2748 
 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 2353 
 c      real fitz(nplanes)        !z coor
+Line 2840 
 c      real fitz(nplanes)        !z coor
              do ip=1,nplanes
                 nplused=nplused+ hit_plane(ip)
              enddo
- c            if(nplused.lt.nplxz_min)goto 888 !next doublet
-             if(nplused.lt.nplyz_min)goto 888 !next doublet
-             if(ncp_ok.lt.ncpok)goto 888 !next cloud
-             if(DEBUG)then
-                print*,'Combination ',iyz,ixz
-      $              ,' db ',ptcloud_yz(iyz)
-      $              ,' tr ',ptcloud_xz(ixz)
-      $              ,'  -----> # matching couples ',ncp_ok
-             endif
- c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
- c$$$  print*,'Configurazione cluster XZ'
- c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
- c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))
- c$$$  print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))
- c$$$  print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))
- c$$$  print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))
- c$$$  print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))
- c$$$  print*,'Configurazione cluster YZ'
- c$$$  print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))
- c$$$  print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))
- c$$$  print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))
- c$$$  print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))
- c$$$  print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))
- c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))
- c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
- *     -------> INITIAL GUESS <-------
- cccc       SBAGLIATO
- c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
- c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
- c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
- c$$$     $           /dreal(alfaxz2_av(ixz)))
- c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
- c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
- c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
- cccc       GIUSTO (ma si sua guess())
- c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
- c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
- c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
- c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
- c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
- c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
- c$$$     $           /dreal(alfaxz2_av(ixz)))
- c$$$            ELSE
- c$$$               AL_INI(4) = acos(-1.)/2
- c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
- c$$$            ENDIF
- c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
- c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
- c$$$
- c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
- c$$$
- c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
-             if(DEBUG)then
+             if(nplused.lt.3)goto 888 !next combination
+ ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
+ ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
+ *     -----------------------------------------------------------
+ *     if in RECOVER_SINGLET mode, the two clouds must have
+ *     at least ONE intersecting real couple
+ *     -----------------------------------------------------------
+             if(ncp_ok.lt.1)goto 888 !next combination
+             if(DEBUG.EQ.1)then
+                print*,'////////////////////////////'
+                print*,'Cloud combination (Y,X): ',iyz,ixz
+                print*,' db ',ptcloud_yz(iyz)
+                print*,' tr ',ptcloud_xz(ixz)
+                print*,'  -----> # matching couples ',ncp_ok
+                print*,'  -----> # fake couples (X)',ncpx_ok
+                print*,'  -----> # fake couples (Y)',ncpy_ok
+                do icp=1,ncp_tot
+                   print*,'cp ',icp,' >'
+      $                 ,' x ',cpcloud_xz(ixz,icp)
+      $                 ,' y ',cpcloud_yz(iyz,icp)
+      $                 ,' ==> ',cpintersec(icp)
+                enddo
+             endif
+             if(DEBUG.EQ.1)then
                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
-Line 2442 
 c$$$            if(AL_INI(5).gt.defmax)g
+Line 2899 
 c$$$            if(AL_INI(5).gt.defmax)g
                                hit_plane(6)=icp6
                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
+                               if(DEBUG.eq.1)
+      $                             print*,'combination: '
+      $                             ,cp_match(6,icp1)
+      $                             ,cp_match(5,icp2)
+      $                             ,cp_match(4,icp3)
+      $                             ,cp_match(3,icp4)
+      $                             ,cp_match(2,icp5)
+      $                             ,cp_match(1,icp6)
+ *                             ---------------------------------------
+ *                             check if this group of couples has been
+ *                             already fitted
+ *                             ---------------------------------------
+                               do ica=1,ntracks
+                                  isthesame=1
+                                  do ip=1,NPLANES
+                                     if(hit_plane(ip).ne.0)then
+                                        if(  CP_STORE(nplanes-ip+1,ica)
+      $                                      .ne.
+      $                                      cp_match(ip,hit_plane(ip)) )
+      $                                      isthesame=0
+                                     else
+                                        if(  CP_STORE(nplanes-ip+1,ica)
+      $                                      .ne.
+      $                                      0 )
+      $                                      isthesame=0
+                                     endif
+                                  enddo
+                                  if(isthesame.eq.1)then
+                                     if(DEBUG.eq.1)
+      $                                   print*,'(already fitted)'
+                                     goto 666 !jump to next combination
+                                  endif
+                               enddo
                                call track_init !init TRACK common
-                               do ip=1,nplanes !loop on planes
+                               do ip=1,nplanes !loop on planes (bottom to top)
                                   if(hit_plane(ip).ne.0)then
                                      id=cp_match(ip,hit_plane(ip))
                                      is=is_cp(id)
-Line 2459 
 c$$$            if(AL_INI(5).gt.defmax)g
+Line 2951 
 c$$$            if(AL_INI(5).gt.defmax)g
  *                                   *************************
  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
-Line 2487 
 c$$$                              enddo
+Line 3010 
 c$$$                              enddo
                                ifail=0 !error flag in chi^2 computation
                                jstep=0 !number of  minimization steps
                                iprint=0
- c                              if(DEBUG)iprint=1
+ c                              if(DEBUG.EQ.1)iprint=1
-                               if(DEBUG)iprint=2
+                               if(DEBUG.EQ.1)iprint=2
                                call mini2(jstep,ifail,iprint)
                                if(ifail.ne.0) then
-                                  if(DEBUG)then
+                                  if(DEBUG.EQ.1)then
                                      print *,
       $                              '*** MINIMIZATION FAILURE *** '
       $                              //'(clouds_to_ctrack)'
-Line 2510 
 c                                 chi2=-
+Line 3033 
 c                                 chi2=-
  *     **********************************************************
                                if(chi2.le.0.)goto 666
+                               if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
+                               if(chi2.ne.chi2)goto 666  !OPTIMIZATION
  *     --------------------------
  *     STORE candidate TRACK INFO
  *     --------------------------
                                if(ntracks.eq.NTRACKSMAX)then
-                                  if(verbose)print*,
+                                  if(verbose.eq.1)print*,
       $                 '** warning ** number of candidate tracks '//
       $                 ' exceeds vector dimension '
       $                ,'( ',NTRACKSMAX,' )'
  c                                 good2=.false.
  c                                 goto 880 !fill ntp and go to next event
                                   do iv=1,nviews
-                                     mask_view(iv) = 7
+ c                                    mask_view(iv) = 7
+                                     mask_view(iv) = mask_view(iv) + 2**6
                                   enddo
                                   iflag=1
                                   return
-Line 2531 
 c                                 goto 8
+Line 3057 
 c                                 goto 8
                                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))
                                   XM_STORE(ip,ntracks)=sngl(xm(ip))
                                   YM_STORE(ip,ntracks)=sngl(ym(ip))
                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))
-Line 2551 
 c$$$     $
+Line 3074 
 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
+                               RCHI2_STORE(ntracks)=REAL(chi2)
- c$$$  ndof=ndof-5
- c$$$  *                             reduced chi^2
- c$$$  rchi2=chi2/dble(ndof)
-                               RCHI2_STORE(ntracks)=chi2
  *     --------------------------------
  *     STORE candidate TRACK INFO - end
-Line 2587 
 c$$$  rchi2=chi2/dble(ndof)
+Line 3126 
 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 2616 
 c$$$  rchi2=chi2/dble(ndof)
+Line 3161 
 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)
- cccccc 12/08/2006 modified by elena vannucicni ---> (3)
- c******************************************************
        include 'commontracker.f'
        include 'level1.f'
-Line 2628 
 c***************************************
+Line 3168 
 c***************************************
        include 'common_xyzPAM.f'
        include 'common_mini_2.f'
        include 'common_mech.f'
- c      include 'momanhough_init.f'
- c      include 'level1.f'
        include 'calib.f'
  *     flag to chose PFA
        character*10 PFA
        common/FINALPFA/PFA
+       double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7
+       double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7
+       double precision zmm,zmm_A,zmm_B !EM GCC4.7
+       double precision clincnewc !EM GCC4.7
+       double precision clincnew !EM GCC4.7
+       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 2645 
 c      include 'level1.f'
+Line 3196 
 c      include 'level1.f'
        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
-Line 2652 
 c      include 'level1.f'
+Line 3214 
 c      include 'level1.f'
  *     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 2665 
 c      include 'level1.f'
+Line 3262 
 c      include 'level1.f'
       $           ,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
  *     |||||||||||||||||||||||||||||||||||||||||||||||||
  *     -------------------------------------------------
-Line 2696 
 c            dedxtrk(nplanes-ip+1) = (de
+Line 3327 
 c            dedxtrk(nplanes-ip+1) = (de
              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 2718 
 c            dedxtrk(nplanes-ip+1) = (de
+Line 3353 
 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 2732 
 c            if(DEBUG)print*,'>>>> try t
+Line 3366 
 c            if(DEBUG)print*,'>>>> try t
              do icp=1,ncp_plane(ip) !loop on couples on plane icp
                 icx=clx(ip,icp)
                 icy=cly(ip,icp)
+                if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
-      $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)
+      $              cl_used(icx).ne.0.or. !or the X cluster is already used
-      $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
+      $              cl_used(icy).ne.0.or. !or the Y cluster is already used
       $              .false.)goto 1188 !then jump to next couple.
  *
                 call xyz_PAM(icx,icy,ist,
       $              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)
-                distance = distance / RCHI2_STORE(ibest)!<<< MS
+ 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 2758 
 c     $              'ETA2','ETA2',
+Line 3396 
 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.*dsqrt(rymm**2+rxmm**2
+      $            +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))))! EM GCC4.7
                 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 2805 
 c            if(DEBUG)print*,'>>>> try t
+Line 3442 
 c            if(DEBUG)print*,'>>>> try t
              do icp=1,ncp_plane(ip) !loop on cluster inside couples
                 icx=clx(ip,icp)
                 icy=cly(ip,icp)
+                if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
                 id=id_cp(ip,icp,ist)
                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
  *                                                !jump to the next couple
-Line 2812 
 c            if(DEBUG)print*,'>>>> try t
+Line 3450 
 c            if(DEBUG)print*,'>>>> try t
  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)
-                distance = distance / RCHI2_STORE(ibest)!<<< MS
+ 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 2831 
 c     $              'ETA2','ETA2',
+Line 3475 
 c     $              'ETA2','ETA2',
                    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
-Line 2840 
 c                  dedxmm = dedx(icx) !(
+Line 3484 
 c                  dedxmm = dedx(icx) !(
  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)
-                distance = distance / RCHI2_STORE(ibest)!<<< MS
+ c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG)print*,'( cl-Y ',icy
+                if(DEBUG.EQ.1)
-      $              ,' in cp ',id,' ) normalized distance ',distance
+      $              print*,'( cl-Y ',icy
+      $              ,' in cp ',id,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
                    ymm_A = yPAM_A
-Line 2859 
 c     $              'ETA2','ETA2',
+Line 3509 
 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)
- c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
       $              .false.)goto 18882      !jump to the next singlet
                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
                    call xyz_PAM(icl,0,ist,
- 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)
-                distance = distance / RCHI2_STORE(ibest)!<<< MS
+ c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG)print*,'( cl-s ',icl
+                if(DEBUG.EQ.1)
-      $              ,' ) normalized distance ',distance
+      $              print*,'( cl-s ',icl
+      $              ,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
                    ymm_A = yPAM_A
-Line 2899 
 c     $                 'ETA2','ETA2',
+Line 3553 
 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 !<<<<
- *              ----------------------------
- c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
                 if(mod(VIEW(iclm),2).eq.0)then
-                   XGOOD(nplanes-ip+1)=1.
+                   clincnew=
-                   resx(nplanes-ip+1)=rxmm
+      $            20.*     !EM GCC4.7
-                   if(DEBUG)print*,'%%%% included X-cl ',iclm
+      $            dsqrt(rxmm**2 +
- c                  if(.true.)print*,'%%%% included X-cl ',iclm
+      $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1))
-      $                 ,'( chi^2, ',RCHI2_STORE(ibest)
+                else if(mod(VIEW(iclm),2).ne.0)then
-      $                 ,', norm.dist.= ',distmin
+                   clincnew=
-      $                 ,', cut ',clinc,' )'
+      $            10.* !EM GCC4.7
-                else
+      $            dsqrt(rymm**2 +
-                   YGOOD(nplanes-ip+1)=1.
+      $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2))
-                   resy(nplanes-ip+1)=rymm
+                endif
-                   if(DEBUG)print*,'%%%% included Y-cl ',iclm
- c                  if(.true.)print*,'%%%% included Y-cl ',iclm
+                if(distmin.le.clincnew)then
-      $                 ,'( chi^2, ',RCHI2_STORE(ibest)
-      $                 ,', norm.dist.= ', distmin
+                   CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
-      $                 ,', cut ',clinc,' )'
+ *     ----------------------------
+                   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
- c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
- *              ----------------------------
-                xm_A(nplanes-ip+1) = xmm_A
-                ym_A(nplanes-ip+1) = ymm_A
-                xm_B(nplanes-ip+1) = xmm_B
-                ym_B(nplanes-ip+1) = ymm_B
-                zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
- 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 2954 
 c              dedxtrk(nplanes-ip+1) = d
+Line 3618 
 c              dedxtrk(nplanes-ip+1) = d
        return
        end
  ***************************************************
  *                                                 *
  *                                                 *
-Line 2962 
 c              dedxtrk(nplanes-ip+1) = d
+Line 3627 
 c              dedxtrk(nplanes-ip+1) = d
  *                                                 *
  *                                                 *
  **************************************************
- cccccc 12/08/2006 modified by elena ---> (1)
  *
-       subroutine clean_XYclouds(ibest,iflag)
-       include 'commontracker.f'
-       include 'level1.f'
-       include 'common_momanhough.f'
- c      include 'momanhough_init.f'
-       include 'level2.f'        !(1)
- c      include 'calib.f'
- c      include 'level1.f'
-       do ip=1,nplanes           !loop on planes
-          id=CP_STORE(nplanes-ip+1,ibest)
-          icl=CLS_STORE(nplanes-ip+1,ibest)
-          if(id.ne.0.or.icl.ne.0)then
-             if(id.ne.0)then
-                iclx=clx(ip,icp_cp(id))
-                icly=cly(ip,icp_cp(id))
- c               cl_used(iclx)=1  !tag used clusters
- c               cl_used(icly)=1  !tag used clusters
-                cl_used(iclx)=ntrk  !tag used clusters !(1)
-                cl_used(icly)=ntrk  !tag used clusters !(1)
-             elseif(icl.ne.0)then
- c               cl_used(icl)=1   !tag used clusters
-                cl_used(icl)=ntrk   !tag used clusters !1)
-             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
-Line 3056 
 c               endif
+Line 3639 
 c               endif
        include 'level1.f'
        include 'common_momanhough.f'
        include 'level2.f'
- c      include 'level1.f'
+ *     ---------------------------------
+ *     variables initialized from level1
+ *     ---------------------------------
        do i=1,nviews
           good2(i)=good1(i)
+          do j=1,nva1_view
+             vkflag(i,j)=1
+             if(cnnev(i,j).le.0)then
+                vkflag(i,j)=cnnev(i,j)
+             endif
+          enddo
        enddo
+ *     ----------------
+ *     level2 variables
+ *     ----------------
        NTRK = 0
        do it=1,NTRKMAX
           IMAGE(IT)=0
-Line 3073 
 c      include 'level1.f'
+Line 3665 
 c      include 'level1.f'
              ZM_nt(IP,IT) = 0
              RESX_nt(IP,IT) = 0
              RESY_nt(IP,IT) = 0
+             TAILX_nt(IP,IT) = 0
+             TAILY_nt(IP,IT) = 0
+             XBAD(IP,IT) = 0
+             YBAD(IP,IT) = 0
              XGOOD_nt(IP,IT) = 0
              YGOOD_nt(IP,IT) = 0
+             LS(IP,IT) = 0
              DEDX_X(IP,IT) = 0
              DEDX_Y(IP,IT) = 0
              CLTRX(IP,IT) = 0
              CLTRY(IP,IT) = 0
+             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 3098 
 c      include 'level1.f'
+Line 3701 
 c      include 'level1.f'
          ys(1,ip)=0
          ys(2,ip)=0
          sgnlys(ip)=0
+         sxbad(ip)=0
+         sybad(ip)=0
+         multmaxsx(ip)=0
+         multmaxsy(ip)=0
        enddo
        end
-Line 3207 
 c      include 'level1.f'
+Line 3814 
 c      include 'level1.f'
        include 'commontracker.f'
- c      include 'level1.f'
        include 'level1.f'
        include 'common_momanhough.f'
        include 'level2.f'
        include 'common_mini_2.f'
-       real sinth,phi,pig
+       include 'calib.f'
+       character*10 PFA
+       common/FINALPFA/PFA
+       real sinth,phi,pig, npig ! EM GCC4.7
+       integer ssensor,sladder
        pig=acos(-1.)
+ *     -------------------------------------
        chi2_nt(ntr)        = sngl(chi2)
        nstep_nt(ntr)       = nstep
+ *     -------------------------------------
-       phi   = al(4)
+       phi   = REAL(al(4))
-       sinth = al(3)
+       sinth = REAL(al(3))
        if(sinth.lt.0)then
           sinth = -sinth
           phi = phi + pig
-Line 3230 
 c      include 'level1.f'
+Line 3845 
 c      include 'level1.f'
       $     phi = phi + 2*pig
        al(4) = phi
        al(3) = sinth
        do i=1,5
           al_nt(i,ntr)     = sngl(al(i))
           do j=1,5
              coval(i,j,ntr) = sngl(cov(i,j))
           enddo
        enddo
+ *     -------------------------------------
        do ip=1,nplanes           ! loop on planes
           xgood_nt(ip,ntr) = int(xgood(ip))
           ygood_nt(ip,ntr) = int(ygood(ip))
-Line 3246 
 c      include 'level1.f'
+Line 3860 
 c      include 'level1.f'
           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))
-          dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
-          dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)
+          factor = sqrt(
+      $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
+      $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
+      $        1. )
+          dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
+          dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)
-          id  = CP_STORE(ip,IDCAND)
+ 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)
-          if(id.ne.0)then
+          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
+ c         if(id.ne.0)then
+ CCCCCC(10 novembre 2009) PATCH X NUCLEI
+ C     non ho capito perche', ma durante il ritracciamento dei nuclei
+ C     (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
+ C     che non viene reinizializzata correttamente e i cluster esclusi
+ C     dal fit risultano ancora inclusi...
+ C
+          cltrx(ip,ntr)   = 0
+          cltry(ip,ntr)   = 0
+ c$$$         if(
+ c$$$     $        xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
+ c$$$     $        .and.
+ c$$$     $        id.ne.0)then
+          if(id.ne.0)then        !patch 30/12/09
+ c           >>> is a couple
              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
- c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
-          elseif(icl.ne.0)then
+             if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
-             if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl
-             if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl
+                cl_used(cltrx(ip,ntr)) = 1 !tag used clusters
- c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
+                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
+ c$$$         elseif(icl.ne.0)then
+          endif                  !patch 30/12/09
+          if(icl.ne.0)then       !patch 30/12/09
+             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
-Line 3280 
 c            print*,ip,' ',cltrx(ip,ntr)
+Line 4018 
 c            print*,ip,' ',cltrx(ip,ntr)
  *     -------------------------------------------------------
        include 'commontracker.f'
- c      include 'level1.f'
        include 'calib.f'
        include 'level1.f'
        include 'common_momanhough.f'
-Line 3288 
 c      include 'level1.f'
+Line 4025 
 c      include 'level1.f'
        include 'common_xyzPAM.f'
  *     count #cluster per plane not associated to any track
- c      good2=1!.true.
        nclsx = 0
        nclsy = 0
        do iv = 1,nviews
-          if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
+ 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.)
-                   xs(is,nclsx) = (xPAM_A+xPAM_B)/2
+                   call xyz_PAM(icl,0,is,PFAdef,'    ',0.,0.,0.,0.)
+                   xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7
                 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.)
-                   ys(is,nclsy) = (yPAM_A+yPAM_B)/2
+                   call xyz_PAM(0,icl,is,'    ',PFAdef,0.,0.,0.,0.)
+                   ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7
                 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
-Line 3397 
 c      print*,icl,cl_used(icl),cl_good(i
+Line 4150 
 c      print*,icl,cl_used(icl),cl_good(i
                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)
                 nnn=nnn+ptcloud_yz(iyz)
              enddo
-             do ipt=1,nnn
+             do ipt=1,min(ndblt_max_nt,nnn)
                 db_cloud_nt(ipt)=db_cloud(ipt)
               enddo
           endif
-Line 3410 
 c      print*,icl,cl_used(icl),cl_good(i
+Line 4163 
 c      print*,icl,cl_used(icl),cl_good(i
                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
                 nnn=nnn+ptcloud_xz(ixz)
              enddo
-             do ipt=1,nnn
+             do ipt=1,min(ntrpt_max_nt,nnn)
                tr_cloud_nt(ipt)=tr_cloud(ipt)
               enddo
           endif

 Legend:



Removed from v.1.16
 


changed lines


 
Added in v.1.46
 Legend:



Removed from v.1.16
 


changed lines


 
Added in v.1.46
-Removed from v.1.16
+Added in v.1.46

	ViewVC Help
Powered by ViewVC 1.1.23