/[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.33 by pam-fi,
Tue Nov 27 15:28:57 2007 UTC
+revision 1.47 by pam-ts,
Wed Jun  4 07:57:04 2014 UTC
 Line 20
        include 'calib.f'
        include 'level2.f'
  c      print*,'======================================================'
  c$$$      do ic=1,NCLSTR1
  c$$$         if(.false.
-Line 74 
 c$$$      enddo
+Line 72 
 c$$$      enddo
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
        call cl_to_couples(iflag)
        if(iflag.eq.1)then        !bad event
           goto 880               !go to next event
        endif
+       if(ncp_tot.eq.0)goto 880  !go to next event
  *-----------------------------------------------------
  *-----------------------------------------------------
  *     HOUGH TRASFORM
-Line 113 
 c$$$      enddo
+Line 110 
 c$$$      enddo
        if(iflag.eq.1)then        !bad event
           goto 880               !go to next event
        endif
+       if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event
  *-------------------------------------------------------------------------------
-Line 163 
 c$$$      enddo
+Line 161 
 c$$$      enddo
        if(iflag.eq.1)then        !bad event
           goto 880               !fill ntp and go to next event
        endif
+ *     ------------------------------------------------
+ *     first try to release the tolerance
+ *     ------------------------------------------------
        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
-         if(cutdistyz.lt.maxcuty/2)then
+          if(cutdistyz.lt.maxcuty/2)then
-           cutdistyz=cutdistyz+cutystep
+             cutdistyz=cutdistyz+cutystep
-         else
+          else
-           cutdistyz=cutdistyz+(3*cutystep)
+             cutdistyz=cutdistyz+(3*cutystep)
-         endif
+          endif
-         goto 878
+          if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
-       endif
+          goto 878
+       endif
+ *     ------------------------------------------------
+ *     hence reduce the minimum number of plane
+ *     ------------------------------------------------
+       if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then
+          nplyz_min=nplyz_min-1
+          if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min
+          goto 878
+       endif
+       if(nclouds_yz.eq.0)goto 880 !go to next event
-       if(planehit.eq.3) goto 881
+ ccc   if(planehit.eq.3) goto 881
  continue
        call trip_to_XZcloud(iflag)
        if(iflag.eq.1)then        !bad event
           goto 880               !fill ntp and go to next event
        endif
+ *     ------------------------------------------------
+ *     first try to release the tolerance
+ *     ------------------------------------------------
        if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
          cutdistxz=cutdistxz+cutxstep
+          if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz
          goto 879
        endif
+ *     ------------------------------------------------
+ *     hence reduce the minimum number of plane
+ *     ------------------------------------------------
+       if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then
+          nplxz_min=nplxz_min-1
+          if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min
+          goto 879
+       endif
+       if(nclouds_xz.eq.0)goto 880 !go to next event
- continue
+ c$$$ 881  continue
- *     if there is at least three planes on the Y view decreases cuts on X view
+ c$$$*     if there is at least three planes on the Y view decreases cuts on X view
-       if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
+ c$$$      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
-      $     nplxz_min.ne.y_min_start)then
+ c$$$     $     nplxz_min.ne.y_min_start)then
-         nptxz_min=x_min_step
+ c$$$        nptxz_min=x_min_step
-         nplxz_min=x_min_start-x_min_step
+ c$$$        nplxz_min=x_min_start-x_min_step
-         goto 879
+ c$$$        goto 879
-       endif
+ c$$$      endif
  return
        end
-Line 243 
 c      include 'momanhough_init.f'
+Line 269 
 c      include 'momanhough_init.f'
  *
  *-------------------------------------------------------------------------------
  *-------------------------------------------------------------------------------
-          ntrk=0                 !counter of identified physical tracks
+ ccc         ntrk=0                 !counter of identified physical tracks
-   continue               !<<<<<<< come here when performing a new search
+ c11111    continue               !<<<<<<< come here when performing a new search
+          continue                  !<<<<<<< come here when performing a new search
+          if(nclouds_xz.eq.0)goto 880 !go to next event
+          if(nclouds_yz.eq.0)goto 880 !go to next event
  c         iflag=0
           call clouds_to_ctrack(iflag)
           if(iflag.eq.1)then     !no candidate tracks found
              goto 880            !fill ntp and go to next event
           endif
+          if(ntracks.eq.0)goto 880 !go to next event
           FIMAGE=.false.         !processing best track (not track image)
           ibest=0                !best track among candidates
-Line 272 
 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
           do i=1,ntracks
             ndof=0
-Line 295 
 c$$$         if(ibest.eq.0)goto 880 !>>
+Line 326 
 c$$$         if(ibest.eq.0)goto 880 !>>
             endif
           enddo
- c$$$         rchi2best=1000000000.
- c$$$         ndofbest=0             !(1)
- c$$$         do i=1,ntracks
- c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
- c$$$     $          RCHI2_STORE(i).gt.0)then
- c$$$             ndof=0             !(1)
- c$$$             do ii=1,nplanes    !(1)
- c$$$               ndof=ndof        !(1)
- c$$$     $              +int(xgood_store(ii,i)) !(1)
- c$$$     $              +int(ygood_store(ii,i)) !(1)
- c$$$             enddo              !(1)
- c$$$             if(ndof.ge.ndofbest)then !(1)
- c$$$               ibest=i
- c$$$               rchi2best=RCHI2_STORE(i)
- c$$$               ndofbest=ndof    !(1)
- c$$$             endif              !(1)
- c$$$           endif
- c$$$         enddo
           if(ibest.eq.0)goto 880 !>> no good candidates
  *-------------------------------------------------------------------------------
-Line 350 
 c$$$         enddo
+Line 363 
 c$$$         enddo
           do i=1,5
              AL_GUESS(i)=AL(i)
           enddo
- c         print*,'## guess: ',al
           do i=1,5
              AL(i)=dble(AL_STORE(i,icand))
-Line 365 
 c         if(DEBUG.EQ.1)iprint=1
+Line 377 
 c         if(DEBUG.EQ.1)iprint=1
           if(VERBOSE.EQ.1)iprint=1
           if(DEBUG.EQ.1)iprint=2
           call mini2(jstep,ifail,iprint)
-          if(ifail.ne.0) then
+ cc         if(ifail.ne.0) then
+          if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
+             if(CHI2.ne.CHI2)CHI2=-9999. !new
              if(VERBOSE.EQ.1)then
                 print *,
       $              '*** MINIMIZATION FAILURE *** (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.EQ.1)then
-Line 424 
 c         rchi2=chi2/dble(ndof)
+Line 430 
 c         rchi2=chi2/dble(ndof)
              if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
              goto 122            !jump
           endif
- c         print*,'is1 ',is1
           do ip=1,NPLANES
              is2 = SENSOR_STORE(ip,icand) !sensor
- c            print*,'is2 ',is2,' ip ',ip
              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! '
- *     now search for track-image among track candidates
- c$$$         do i=1,ntracks
- c$$$            iimage=i
- c$$$            do ip=1,nplanes
- c$$$               if(     CP_STORE(nplanes-ip+1,icand).ne.
- c$$$     $              -1*CP_STORE(nplanes-ip+1,i).and.
- c$$$     $              CP_STORE(nplanes-ip+1,i).ne.0.and.
- c$$$     $              CP_STORE(nplanes-ip+1,icand).ne.0 )iimage=0
- c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
- c$$$     $              ,CP_STORE(nplanes-ip+1,icand),iimage
- c$$$            enddo
- c$$$            if(  iimage.ne.0.and.
- c$$$c     $           RCHI2_STORE(i).le.CHI2MAX.and.
- c$$$c     $           RCHI2_STORE(i).gt.0.and.
- c$$$     $           .true.)then
- c$$$               if(DEBUG.EQ.1)print*,'Track candidate ',iimage
- c$$$     $              ,' >>> TRACK IMAGE >>> of'
- c$$$     $              ,ibest
- c$$$               goto 122         !image track found
- c$$$            endif
- c$$$         enddo
  *     ---------------------------------------------------------------
  *     take the candidate with the greatest number of matching couples
  *     if more than one satisfies the requirement,
-Line 486 
 c$$$         enddo
+Line 469 
 c$$$         enddo
                    endif
                 endif
- c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
- c$$$     $              ,CP_STORE(nplanes-ip+1,icand),ncp
              enddo
        continue
              trackimage(i)=ncp   !number of matching couples
-Line 501 
 c$$$     $              ,CP_STORE(nplane
+Line 482 
 c$$$     $              ,CP_STORE(nplane
           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.1)then
+          if(ngood.gt.0)then
  *     -------------------------------------------------------
  *     order track-candidates according to:
  *     1st) decreasing n.points
-Line 532 
 c$$$     $              ,CP_STORE(nplane
+Line 514 
 c$$$     $              ,CP_STORE(nplane
                    endif
                 endif
              enddo
+             if(DEBUG.EQ.1)then
+                print*,'Track candidate ',iimage
+      $              ,' >>> TRACK IMAGE >>> of'
+      $              ,ibest
+             endif
           endif
-          if(DEBUG.EQ.1)print*,'Track candidate ',iimage
-      $        ,' >>> TRACK IMAGE >>> of'
-      $        ,ibest
     continue
  *     --- and store the results --------------------------------
-          ntrk = ntrk + 1                   !counter of found tracks
-          if(.not.FIMAGE
-      $        .and.iimage.eq.0) image(ntrk)= 0
-          if(.not.FIMAGE
-      $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
-          if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
-          call fill_level2_tracks(ntrk)     !==> good2=.true.
- c         print*,'++++++++++ iimage,fimage,ntrk,image '
- c     $        ,iimage,fimage,ntrk,image(ntrk)
           if(ntrk.eq.NTRKMAX)then
              if(verbose.eq.1)
       $           print*,
-Line 562 
 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.EQ.1)then
-                print*,'***** NEW SEARCH ****'
-             endif
-             goto 11111          !try new search
-          endif
- *     **********************************************
     return
        end
-Line 699 
 c     $        rchi2best.lt.15..and.
+Line 660 
 c     $        rchi2best.lt.15..and.
        parameter (ndivx=30)
- c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
        resxPAM = 0
        resyPAM = 0
-Line 713 
 c$$$      print*,icx,icy,sensor,PFAx,PFA
+Line 673 
 c$$$      print*,icx,icy,sensor,PFAx,PFA
        xPAM_B = 0.D0
        yPAM_B = 0.D0
        zPAM_B = 0.D0
- c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
        if(sensor.lt.1.or.sensor.gt.2)then
           print*,'xyz_PAM   ***ERROR*** wrong input '
-Line 757 
 c         resxPAM = RESXAV
+Line 716 
 c         resxPAM = RESXAV
           stripx  = stripx + corr
           resxPAM = res
-  endif
+  continue
+       endif
  *     -----------------
  *     CLUSTER Y
-Line 803 
 c         resyPAM = RESYAV
+Line 763 
 c         resyPAM = RESYAV
           stripy  = stripy + corr
           resyPAM = res
-  endif
+  continue
+       endif
- c$$$      print*,'## stripx,stripy ',stripx,stripy
  c===========================================================
  C     COUPLE
-Line 885 
 C=======================================
+Line 845 
 C=======================================
              nldx = nldy
              viewx = viewy + 1
-             yi   = dcoordsi(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 895 
 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 907 
 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...
                 if(DEBUG.EQ.1) then
-Line 916 
 c            if((stripx.le.3).or.(stripx
+Line 874 
 c            if((stripx.le.3).or.(stripx
       $                 ' WARNING: false X strip: strip ',stripx
                 endif
              endif
-             xi   = dcoordsi(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 932 
 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
-Line 979 
 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 989 
 c     (xPAM,yPAM,zPAM) = measured coordi
+Line 966 
 c     (xPAM,yPAM,zPAM) = measured coordi
  c                        in PAMELA reference system
  c------------------------------------------------------------------------
-          xPAM = 0.D0
+          xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
-          yPAM = 0.D0
+          yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
-          zPAM = 0.D0
+          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 1002 
 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
-Line 1013 
 c         print*,'A-(',xPAM_A,yPAM_A,')
+Line 992 
 c         print*,'A-(',xPAM_A,yPAM_A,')
        endif
- c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
- c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
- c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
  continue
        end
-Line 1063 
 c$$$      PFAy = 'COG4'!PFA
+Line 1039 
 c$$$      PFAy = 'COG4'!PFA
        call idtoc(pfaid,PFAx)
        call idtoc(pfaid,PFAy)
- cc      print*,PFAx,PFAy
- c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
- c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
        if(icx.ne.0.and.icy.ne.0)then
           ipx=npl(VIEW(icx))
           ipy=npl(VIEW(icy))
- c$$$         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
- c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
- c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
           if( (nplanes-ipx+1).ne.ip )then
-             print*,'xyzpam: ***WARNING*** cluster ',icx
+             print*,'xyzpam: ***WARNING*** cluster icx=',icx
-      $           ,' does not belong to plane: ',ip
+      $           ,' belongs to plane ',(nplanes-ipx+1)
+      $           ,' and not ',ip
              icx = -1*icx
              return
           endif
           if( (nplanes-ipy+1).ne.ip )then
-             print*,'xyzpam: ***WARNING*** cluster ',icy
+             print*,'xyzpam: ***WARNING*** cluster icy=',icy
-      $           ,' does not belong to plane: ',ip
+      $           ,' belongs to plane ',(nplanes-ipy+1)
-             icy = -1*icy
+      $           ,' and not ',ip
+              icy = -1*icy
              return
           endif
-Line 1102 
 c$$$     $        ,' does not belong to
+Line 1072 
 c$$$     $        ,' does not belong to
           zm(ip) = zPAM
           xm_A(ip) = 0.D0
           ym_A(ip) = 0.D0
+          zm_A(ip) = 0.D0
           xm_B(ip) = 0.D0
           ym_B(ip) = 0.D0
+          zm_B(ip) = 0.D0
  c         zv(ip) = zPAM
        elseif(icx.eq.0.and.icy.ne.0)then
           ipy=npl(VIEW(icy))
- c$$$         if((nplanes-ipy+1).ne.ip)
- c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
- c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
           if( (nplanes-ipy+1).ne.ip )then
-             print*,'xyzpam: ***WARNING*** cluster ',icy
+             print*,'xyzpam: ***WARNING*** cluster icy=',icy
-      $           ,' does not belong to plane: ',ip
+      $           ,' belongs to plane ',(nplanes-ipy+1)
+      $           ,' and not ',ip
              icy = -1*icy
              return
           endif
-Line 1127 
 c$$$     $        ,' does not belong to
+Line 1097 
 c$$$     $        ,' does not belong to
           resx(ip)  = 1000.
           resy(ip)  = resyPAM
-          xm(ip) = -100.
+ c$$$         xm(ip) = -100.
-          ym(ip) = -100.
+ c$$$         ym(ip) = -100.
-          zm(ip) = (zPAM_A+zPAM_B)/2.
+ c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
+          xm(ip) = xPAM
+          ym(ip) = yPAM
+          zm(ip) = zPAM
           xm_A(ip) = xPAM_A
           ym_A(ip) = yPAM_A
+          zm_A(ip) = zPAM_A
           xm_B(ip) = xPAM_B
           ym_B(ip) = yPAM_B
+          zm_B(ip) = zPAM_B
  c         zv(ip) = (zPAM_A+zPAM_B)/2.
        elseif(icx.ne.0.and.icy.eq.0)then
           ipx=npl(VIEW(icx))
- c$$$         if((nplanes-ipx+1).ne.ip)
- c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
- c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
           if( (nplanes-ipx+1).ne.ip )then
-             print*,'xyzpam: ***WARNING*** cluster ',icx
+             print*,'xyzpam: ***WARNING*** cluster icx=',icx
-      $           ,' does not belong to plane: ',ip
+      $           ,' belongs to plane ',(nplanes-ipx+1)
+      $           ,' and not ',ip
              icx = -1*icx
              return
           endif
-Line 1158 
 c$$$     $        ,' does not belong to
+Line 1131 
 c$$$     $        ,' does not belong to
           resx(ip)  = resxPAM
           resy(ip)  = 1000.
-          xm(ip) = -100.
+ c$$$         xm(ip) = -100.
-          ym(ip) = -100.
+ c$$$         ym(ip) = -100.
-          zm(ip) = (zPAM_A+zPAM_B)/2.
+ c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
+          xm(ip) = xPAM
+          ym(ip) = yPAM
+          zm(ip) = zPAM
           xm_A(ip) = xPAM_A
           ym_A(ip) = yPAM_A
+          zm_A(ip) = zPAM_A
           xm_B(ip) = xPAM_B
           ym_B(ip) = yPAM_B
+          zm_B(ip) = zPAM_B
  c         zv(ip) = (zPAM_A+zPAM_B)/2.
-Line 1174 
 c         zv(ip) = (zPAM_A+zPAM_B)/2.
+Line 1152 
 c         zv(ip) = (zPAM_A+zPAM_B)/2.
           if(lad.ne.0)il=lad
           is = 1
           if(sensor.ne.0)is=sensor
- c         print*,nplanes-ip+1,il,is
           xgood(ip) = 0.
           ygood(ip) = 0.
-Line 1186 
 c         print*,nplanes-ip+1,il,is
+Line 1163 
 c         print*,nplanes-ip+1,il,is
           zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
           xm_A(ip) = 0.
           ym_A(ip) = 0.
+          zm_A(ip) = 0.
           xm_B(ip) = 0.
           ym_B(ip) = 0.
+          zm_B(ip) = 0.
  c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
        endif
        if(DEBUG.EQ.1)then
- c         print*,'----------------------------- track coord'
    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
           write(*,22222)ip,zm(ip),xm(ip),ym(ip)
       $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
       $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
- c$$$         print*,'-----------------------------'
        endif
        end
-Line 1242 
 c$$$         print*,'-------------------
+Line 1219 
 c$$$         print*,'-------------------
  *     ----------------------
        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 1286 
 c$$$         print*,'-------------------
+Line 1263 
 c$$$         print*,'-------------------
  cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
           distance=dsqrt(distance)
- c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
- c$$$     $        ,' --- distance_to --- ',xpp,ypp
- c$$$         print*,' resolution ',re
  *     ----------------------
        elseif(
-      +     xPAM.ne.0.and.
+ c$$$     +     xPAM.ne.0.and.
-      +     yPAM.ne.0.and.
+ c$$$     +     yPAM.ne.0.and.
-      +     zPAM.ne.0.and.
+ c$$$     +     zPAM.ne.0.and.
       +     xPAM_A.eq.0.and.
       +     yPAM_A.eq.0.and.
       +     zPAM_A.eq.0.and.
-Line 1316 
 c$$$     $        +
+Line 1290 
 c$$$     $        +
  c$$$     $        ((yPAM-YPP)/resyPAM)**2
           distance=dsqrt(distance)
- c$$$         print*,xPAM,yPAM,zPAM
- c$$$     $        ,' --- distance_to --- ',xpp,ypp
- c$$$         print*,' resolution ',resxPAM,resyPAM
        else
- c         print*
- c     $        ,' function distance_to ---> wrong usage!!!'
- c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
- c         print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
- c     $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
        endif
        distance_to = sngl(distance)
-Line 1391 
 c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
+Line 1357 
 c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
  c------------------------------------------------------------------------
  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
  c------------------------------------------------------------------------
- c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
- c     $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
- c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
- c                  print*,'whichsensor: ',
- c     $                ' WARNING: false X strip: strip ',stripx
- c               endif
                 xi = dcoordsi(stripx,viewx)
                 yi = dcoordsi(stripy,viewy)
                 zi = 0.D0
-Line 1424 
 c---------------------------------------
+Line 1384 
 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 1433 
 c     $              ,iv,xvv(iv),yvv(iv)
+Line 1391 
 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 1446 
 c     $              ,iv,xvv(iv),yvv(iv)
+Line 1404 
 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 1550 
 c      include 'common_analysis.f'
+Line 1508 
 c      include 'common_analysis.f'
        is_cp=0
        if(id.lt.0)is_cp=1
        if(id.gt.0)is_cp=2
- c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
        return
        end
-Line 1649 
 c      include 'level1.f'
+Line 1606 
 c      include 'level1.f'
        integer iflag
        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 1666 
 c      include 'level1.f'
+Line 1625 
 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 1695 
 c            mask_view(iv) = 1
+Line 1654 
 c            mask_view(iv) = 1
        do icx=1,nclstr1          !loop on cluster (X)
           if(mod(VIEW(icx),2).eq.1)goto 10
+          if(cl_used(icx).ne.0)goto 10
  *     ----------------------------------------------------
  *     jump masked views (X VIEW)
  *     ----------------------------------------------------
-Line 1702 
 c            mask_view(iv) = 1
+Line 1663 
 c            mask_view(iv) = 1
  *     ----------------------------------------------------
  *     cut on charge (X VIEW)
  *     ----------------------------------------------------
-          if(sgnl(icx).lt.dedx_x_min)then
+          if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
              cl_single(icx)=0
              goto 10
           endif
-Line 1743 
 c     endif
+Line 1704 
 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 1752 
 c     endif
+Line 1715 
 c     endif
  *     ----------------------------------------------------
  *     cut on charge (Y VIEW)
  *     ----------------------------------------------------
-             if(sgnl(icy).lt.dedx_y_min)then
+             if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
                 cl_single(icy)=0
                 goto 20
              endif
-Line 1830 
 c                  cut = chcut * sch(npl
+Line 1793 
 c                  cut = chcut * sch(npl
       $                 ,'( ',ncouplemax,' ) --> masked!'
  c                  mask_view(nviewx(nplx)) = 2
  c                  mask_view(nviewy(nply)) = 2
-                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
+                  mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
-                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
+                  mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
                    goto 10
                 endif
-Line 1851 
 c                  mask_view(nviewy(nply
+Line 1814 
 c                  mask_view(nviewy(nply
      continue
        enddo                     !end loop on clusters(X)
        do icl=1,nclstr1
           if(cl_single(icl).eq.1)then
              ip=npl(VIEW(icl))
-Line 1859 
 c                  mask_view(nviewy(nply
+Line 1821 
 c                  mask_view(nviewy(nply
              cls(ip,ncls(ip))=icl
           endif
        enddo
+ c 80   continue
+       continue
        if(DEBUG.EQ.1)then
           print*,'clusters  ',nclstr1
           print*,'good    ',(cl_good(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 1933 
 c      double precision xm3,ym3,zm3
+Line 1963 
 c      double precision xm3,ym3,zm3
  *     --------------------------------------------
        do ip=1,nplanes
           if(ncp_plane(ip).gt.ncouplelimit)then
- c            mask_view(nviewx(ip)) = 8
- c            mask_view(nviewy(ip)) = 8
              mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
              mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
           endif
-Line 1945 
 c            mask_view(nviewy(ip)) = 8
+Line 1973 
 c            mask_view(nviewy(ip)) = 8
        ntrpt=0                   !number of triplets
        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
+ c$$$         print*,'(1) ip ',ip1
+ c$$$     $        ,mask_view(nviewx(ip1))
+ c$$$     $        ,mask_view(nviewy(ip1))
           if(  mask_view(nviewx(ip1)).ne.0 .or.
       $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
           do is1=1,2             !loop on sensors - COPPIA 1
              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
                 icx1=clx(ip1,icp1)
                 icy1=cly(ip1,icp1)
+ c$$$               print*,'(1) ip ',ip1,' icp ',icp1
  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
-                xm1=xPAM
+                xm1=REAL(xPAM) !EM GCC4.7
-                ym1=yPAM
+                ym1=REAL(yPAM) !EM GCC4.7
-                zm1=zPAM
+                zm1=REAL(zPAM) !EM GCC4.7
- c     print*,'***',is1,xm1,ym1,zm1
                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
+ c$$$                  print*,'(2) ip ',ip2
+ c$$$     $                 ,mask_view(nviewx(ip2))
+ c$$$     $                 ,mask_view(nviewy(ip2))
                    if(  mask_view(nviewx(ip2)).ne.0 .or.
       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
                          icx2=clx(ip2,icp2)
                          icy2=cly(ip2,icp2)
+ c$$$                        print*,'(2) ip ',ip2,' icp ',icp2
  c                        call xyz_PAM
  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
  c                        call xyz_PAM
  c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
                          call xyz_PAM
       $                       (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)
-Line 1985 
 c     $                       (icx2,icy2
+Line 2034 
 c     $                       (icx2,icy2
       $                          '** warning ** number of identified '//
       $                          'doublets exceeds vector dimention '
       $                          ,'( ',ndblt_max,' )'
- c                           good2=.false.
+ c     good2=.false.
- c                           goto 880 !fill ntp and go to next event
+ c     goto 880 !fill ntp and go to next event
                             do iv=1,12
- c                              mask_view(iv) = 3
+ c     mask_view(iv) = 3
                                mask_view(iv) = mask_view(iv)+ 2**2
                             enddo
                             iflag=1
                             return
                          endif
+ ccc                        print*,'<doublet> ',icp1,icp2
                          ndblt = ndblt + 1
  *     store doublet info
                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
-Line 2001 
 c                              mask_view
+Line 2054 
 c                              mask_view
  *     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 2034 
 c$$$
+Line 2093 
 c$$$
                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
                                   icx3=clx(ip3,icp3)
                                   icy3=cly(ip3,icp3)
+ c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
+ *     ---------------------------------------------------
+ *     all three couples must have a x-cluster
+ *     (condition necessary when in RECOVER_SINGLETS mode)
+ *     ---------------------------------------------------
+                                  if(
+      $                                icx1.eq.0.or.
+      $                                icx2.eq.0.or.
+      $                                icx3.eq.0.or.
+      $                                .false.)goto 29
+                                  if(cl_used(icx1).ne.0)goto 29
+                                  if(cl_used(icx2).ne.0)goto 29
+                                  if(cl_used(icx3).ne.0)goto 29
  c                                 call xyz_PAM
  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
  c                                 call xyz_PAM
-Line 2041 
 c     $                               (i
+Line 2117 
 c     $                               (i
                                   call xyz_PAM
       $                                (icx3,icy3,is3,PFAdef,PFAdef
       $                                ,0.,0.,0.,0.)
-                                  xm3=xPAM
+                                  xm3=REAL(xPAM) !EM GCC4.7
-                                  ym3=yPAM
+                                  ym3=REAL(yPAM) !EM GCC4.7
-                                  zm3=zPAM
+                                  zm3=REAL(zPAM) !EM GCC4.7
  *     find the circle passing through the three points
- c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
+                                  iflag_t = DEBUG
- c$$$     $                                ,xc,zc,radius,iflag)
-                                  iflag = 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.eq.1)print*,
-      $                     '** warning ** number of identified '//
+      $                                   '** warning **'//
-      $                     'triplets exceeds vector dimention '
+      $                                   ' number of identified '//
-      $                    ,'( ',ntrpt_max,' )'
+      $                                   'triplets exceeds'//
- c                                    good2=.false.
+      $                                   ' vector dimention '
- c                                    goto 880 !fill ntp and go to next event
+      $                                   ,'( ',ntrpt_max,' )'
+ c     good2=.false.
+ c     goto 880 !fill ntp and go to next event
                                      do iv=1,nviews
- c                                       mask_view(iv) = 4
+ c     mask_view(iv) = 4
-                                        mask_view(iv)=mask_view(iv)+ 2**3
+                                        mask_view(iv) =
+      $                                      mask_view(iv)+ 2**3
                                      enddo
                                      iflag=1
                                      return
                                   endif
+ ccc                                 print*,'<triplet> ',icp1,icp2,icp3
                                   ntrpt = ntrpt +1
  *     store triplet info
                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
-                                  if(xc.lt.0)then
+                                  if(radius.ne.0.and.xc.lt.0)then
  *************POSITIVE DEFLECTION
-               alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
+            alfaxz1(ntrpt) = xc+sqrt(radius**2-(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
-Line 2192 
 c      include 'momanhough_init.f'
+Line 2318 
 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 2230 
 ccccc if(db_used(idbref).eq.1)goto 1188
+Line 2353 
 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 2258 
 c     print*,idb1,idb2,distance,' cloud
+Line 2373 
 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 2281 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2396 
 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 2291 
 c     print*,'*   idbref,idb2 ',idbref,i
+Line 2408 
 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
  *     ~~~~~~~~~~~~~~~~~
-Line 2325 
 c               mask_view(iv) = 5
+Line 2440 
 c               mask_view(iv) = 5
  c     ptcloud_yz_nt(nclouds_yz)=npt
           do ipt=1,npt
              db_cloud(npt_tot+ipt) = db_all(ipt)
- c     print*,'>> ',ipt,db_all(ipt)
           enddo
           npt_tot=npt_tot+npt
           if(DEBUG.EQ.1)then
-             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
-Line 2338 
 c     print*,'>> ',ipt,db_all(ipt)
+Line 2451 
 c     print*,'>> ',ipt,db_all(ipt)
              print*,'cpcloud_yz '
       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
              print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
- c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
- c$$$     $           ,ptcloud_yz(nclouds_yz)
- c$$$            print*,'nt-uple: db_cloud(...) = '
- c$$$     $           ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)
           endif
  *     >>> NEW CLOUD <<<
  *     ~~~~~~~~~~~~~~~~~
-Line 2356 
 c$$$     $           ,(db_cloud(iii),iii
+Line 2465 
 c$$$     $           ,(db_cloud(iii),iii
        endif
        if(DEBUG.EQ.1)then
-          print*,'---------------------- '
           print*,'Y-Z total clouds ',nclouds_yz
-          print*,' '
        endif
-Line 2421 
 c      include 'momanhough_init.f'
+Line 2528 
 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 2456 
 c         tr_temp(1)=itr1
+Line 2561 
 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)
  *     ------------------------------------------------------------------------
  *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
  *     ------------------------------------------------------------------------
-Line 2475 
 c         tr_temp(1)=itr1
+Line 2583 
 c         tr_temp(1)=itr1
       $              .true.)istrimage=1
                 if(distance.lt.cutdistxz.or.istrimage.eq.1)then
- c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
-Line 2494 
 c     print*,idb1,idb2,distance,' cloud
+Line 2601 
 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 2515 
 c     print*,'*   itrref,itr2 ',itrref,i
+Line 2622 
 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 2531 
 c     print*,'check cp_used'
+Line 2639 
 c     print*,'check cp_used'
           do ip=1,nplanes
              nplused=nplused+ hit_plane(ip)
           enddo
- c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
-          if(npt.lt.nptxz_min)goto 22288     !next triplet
           if(nplused.lt.nplxz_min)goto 22288 !next triplet
  *     ~~~~~~~~~~~~~~~~~
-Line 2566 
 c               mask_view(iv) = 6
+Line 2672 
 c               mask_view(iv) = 6
           npt_tot=npt_tot+npt
           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*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
-Line 2576 
 c               mask_view(iv) = 6
+Line 2681 
 c               mask_view(iv) = 6
              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 2593 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2694 
 c$$$     $           ,(tr_cloud(iii),iii
         endif
        if(DEBUG.EQ.1)then
-          print*,'---------------------- '
           print*,'X-Z total clouds ',nclouds_xz
-          print*,' '
        endif
-Line 2651 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2750 
 c$$$     $           ,(tr_cloud(iii),iii
        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 2661 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2760 
 c$$$     $           ,(tr_cloud(iii),iii
  *     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
+             ncpx_ok=0           !count n.matching-couples with x cluster
+             ncpy_ok=0           !count n.matching-couples with y cluster
              do icp=1,ncp_tot    !loop over couples
- *     get info on
-                cpintersec(icp)=min(
+                if(.not.RECOVER_SINGLETS)then
-      $              cpcloud_yz(iyz,icp),
+ *     ------------------------------------------------------
-      $              cpcloud_xz(ixz,icp))
+ *     if NOT in RECOVER_SINGLETS mode, take the intersection
-                if(
+ *     between xz yz clouds
-      $    (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.
+                   cpintersec(icp)=min(
-      $              .false.)cpintersec(icp)=0
+      $                 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 2709 
 c$$$     $           ,(tr_cloud(iii),iii
+Line 2852 
 c$$$     $           ,(tr_cloud(iii),iii
              do ip=1,nplanes
                 nplused=nplused+ hit_plane(ip)
              enddo
+             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*,'Combination ',iyz,ixz
+                print*,'////////////////////////////'
-      $              ,' db ',ptcloud_yz(iyz)
+                print*,'Cloud combination (Y,X): ',iyz,ixz
-      $              ,' tr ',ptcloud_xz(ixz)
+                print*,' db ',ptcloud_yz(iyz)
-      $              ,'  -----> # matching couples ',ncp_ok
+                print*,' tr ',ptcloud_xz(ixz)
-             endif
+                print*,'  -----> # matching couples ',ncp_ok
+                print*,'  -----> # fake couples (X)',ncpx_ok
- c            if(nplused.lt.nplxz_min)goto 888 !next combination
+                print*,'  -----> # fake couples (Y)',ncpy_ok
-             if(nplused.lt.nplyz_min)goto 888 !next combination
+                do icp=1,ncp_tot
-             if(ncp_ok.lt.ncpok)goto 888 !next combination
+                   print*,'cp ',icp,' >'
+      $                 ,' x ',cpcloud_xz(ixz,icp)
- c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
+      $                 ,' y ',cpcloud_yz(iyz,icp)
- c$$$  print*,'Configurazione cluster XZ'
+      $                 ,' ==> ',cpintersec(icp)
- c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
+                enddo
- c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))
+             endif
- 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.EQ.1)then
-                print*,'track candidate', ntracks+1
                 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 2801 
 c$$$            if(AL_INI(5).gt.defmax)g
+Line 2911 
 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
-Line 2849 
 c     $
+Line 2969 
 c     $
       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
  *                                   *************************
  *                                   -----------------------------
-                                     xgood(nplanes-ip+1)=1.
+                                     if(icx.gt.0.and.icy.gt.0)then
-                                     ygood(nplanes-ip+1)=1.
+                                        xgood(nplanes-ip+1)=1.
-                                     xm(nplanes-ip+1)=xPAM
+                                        ygood(nplanes-ip+1)=1.
-                                     ym(nplanes-ip+1)=yPAM
+                                        xm(nplanes-ip+1)=xPAM
-                                     zm(nplanes-ip+1)=zPAM
+                                        ym(nplanes-ip+1)=yPAM
-                                     resx(nplanes-ip+1)=resxPAM
+                                        zm(nplanes-ip+1)=zPAM
-                                     resy(nplanes-ip+1)=resyPAM
+                                        resx(nplanes-ip+1)=resxPAM
+                                        resy(nplanes-ip+1)=resyPAM
+                                        if(DEBUG.EQ.1)print*,'(X,Y)'
+      $                                      ,nplanes-ip+1,xPAM,yPAM
+                                     else
+                                        xm_A(nplanes-ip+1) = xPAM_A
+                                        ym_A(nplanes-ip+1) = yPAM_A
+                                        xm_B(nplanes-ip+1) = xPAM_B
+                                        ym_B(nplanes-ip+1) = yPAM_B
+                                        zm(nplanes-ip+1)
+      $                                      = (zPAM_A+zPAM_B)/2.
+                                        resx(nplanes-ip+1) = resxPAM
+                                        resy(nplanes-ip+1) = resyPAM
+                                        if(icx.eq.0.and.icy.gt.0)then
+                                           xgood(nplanes-ip+1)=0.
+                                           ygood(nplanes-ip+1)=1.
+                                           resx(nplanes-ip+1) = 1000.
+                                           if(DEBUG.EQ.1)print*,'(  Y)'
+      $                                         ,nplanes-ip+1,xPAM,yPAM
+                                        elseif(icx.gt.0.and.icy.eq.0)then
+                                           xgood(nplanes-ip+1)=1.
+                                           ygood(nplanes-ip+1)=0.
+                                           if(DEBUG.EQ.1)print*,'(X  )'
+      $                                         ,nplanes-ip+1,xPAM,yPAM
+                                           resy(nplanes-ip+1) = 1000.
+                                        else
+                                           print*,'both icx=0 and icy=0'
+      $                                         ,' ==> IMPOSSIBLE!!'
+                                        endif
+                                     endif
  *                                   -----------------------------
                                   endif
                                enddo !end loop on planes
-Line 2896 
 c                                 chi2=-
+Line 3045 
 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
-Line 2922 
 c                                    mas
+Line 3073 
 c                                    mas
                                   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 2941 
 c                                    mas
+Line 3092 
 c                                    mas
       $                                   cp_match(ip,hit_plane(ip))
                                      SENSOR_STORE(nplanes-ip+1,ntracks)
       $                              = is_cp(cp_match(ip,hit_plane(ip)))
-                                     LADDER_STORE(nplanes-ip+1,ntracks)
-      $                                   = LADDER(
+                                     icl=
       $                                   clx(ip,icp_cp(
       $                                   cp_match(ip,hit_plane(ip)
-      $                                   ))));
+      $                                   )));
+                                     if(icl.eq.0)
+      $                                   icl=
+      $                                   cly(ip,icp_cp(
+      $                                   cp_match(ip,hit_plane(ip)
+      $                                   )));
+                                     LADDER_STORE(nplanes-ip+1,ntracks)
+      $                                   = LADDER(icl);
                                   else
                                      CP_STORE(nplanes-ip+1,ntracks)=0
                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0
-Line 2959 
 c                                    mas
+Line 3118 
 c                                    mas
                                   enddo
                                enddo
-                               RCHI2_STORE(ntracks)=chi2
+                               RCHI2_STORE(ntracks)=REAL(chi2)
  *     --------------------------------
  *     STORE candidate TRACK INFO - end
-Line 2979 
 c                                    mas
+Line 3138 
 c                                    mas
        if(ntracks.eq.0)then
           iflag=1
-          return
+ cc         return
        endif
- c$$$      if(DEBUG.EQ.1)then
- c$$$         print*,'****** TRACK CANDIDATES ***********'
- c$$$         print*,'#         R. chi2        RIG'
- c$$$         do i=1,ntracks
- c$$$            print*,i,' --- ',rchi2_store(i),' --- '
- c$$$     $           ,1./abs(AL_STORE(5,i))
- c$$$         enddo
- c$$$         print*,'***********************************'
- c$$$      endif
        if(DEBUG.EQ.1)then
          print*,'****** TRACK CANDIDATES *****************'
          print*,'#         R. chi2        RIG         ndof'
-Line 3036 
 c$$$      endif
+Line 3186 
 c$$$      endif
        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/
-Line 3052 
 c$$$      endif
+Line 3208 
 c$$$      endif
        call track_init
        do ip=1,nplanes           !loop on planes
+          if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
           xP=XV_STORE(nplanes-ip+1,ibest)
           yP=YV_STORE(nplanes-ip+1,ibest)
           zP=ZV_STORE(nplanes-ip+1,ibest)
-Line 3068 
 c$$$         bxyz(3)=0
+Line 3226 
 c$$$         bxyz(3)=0
  *     using improved PFAs
  *     -------------------------------------------------
  *     |||||||||||||||||||||||||||||||||||||||||||||||||
-          if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
+ c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
+ c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
+ c$$$
+ c$$$            id=CP_STORE(nplanes-ip+1,ibest)
+ c$$$
+ c$$$            is=is_cp(id)
+ c$$$            icp=icp_cp(id)
+ c$$$            if(ip_cp(id).ne.ip)
+ c$$$     $           print*,'OKKIO!!'
+ c$$$     $           ,'id ',id,is,icp
+ c$$$     $           ,ip_cp(id),ip
+ c$$$            icx=clx(ip,icp)
+ c$$$            icy=cly(ip,icp)
+ c$$$c            call xyz_PAM(icx,icy,is,
+ c$$$c     $           PFA,PFA,
+ c$$$c     $           AXV_STORE(nplanes-ip+1,ibest),
+ c$$$c     $           AYV_STORE(nplanes-ip+1,ibest))
+ c$$$            call xyz_PAM(icx,icy,is,
+ c$$$     $           PFA,PFA,
+ c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
+ c$$$     $           AYV_STORE(nplanes-ip+1,ibest),
+ c$$$     $           bxyz(1),
+ c$$$     $           bxyz(2)
+ c$$$     $           )
+ c$$$
+ c$$$            xm(nplanes-ip+1) = xPAM
+ c$$$            ym(nplanes-ip+1) = yPAM
+ c$$$            zm(nplanes-ip+1) = zPAM
+ c$$$            xgood(nplanes-ip+1) = 1
+ c$$$            ygood(nplanes-ip+1) = 1
+ c$$$            resx(nplanes-ip+1) = resxPAM
+ c$$$            resy(nplanes-ip+1) = resyPAM
+ c$$$
+ c$$$            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
+ c$$$            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+          if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
              id=CP_STORE(nplanes-ip+1,ibest)
-Line 3093 
 c     $           AYV_STORE(nplanes-ip+1
+Line 3286 
 c     $           AYV_STORE(nplanes-ip+1
       $           bxyz(2)
       $           )
-             xm(nplanes-ip+1) = xPAM
+             if(icx.gt.0.and.icy.gt.0)then
-             ym(nplanes-ip+1) = yPAM
+                xm(nplanes-ip+1) = xPAM
-             zm(nplanes-ip+1) = zPAM
+                ym(nplanes-ip+1) = yPAM
-             xgood(nplanes-ip+1) = 1
+                zm(nplanes-ip+1) = zPAM
-             ygood(nplanes-ip+1) = 1
+                xm_A(nplanes-ip+1) = 0.
-             resx(nplanes-ip+1) = resxPAM
+                ym_A(nplanes-ip+1) = 0.
-             resy(nplanes-ip+1) = resyPAM
+                xm_B(nplanes-ip+1) = 0.
+                ym_B(nplanes-ip+1) = 0.
-             dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
+                xgood(nplanes-ip+1) = 1
-             dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+                ygood(nplanes-ip+1) = 1
+                resx(nplanes-ip+1) = resxPAM
+                resy(nplanes-ip+1) = resyPAM
+                dedxtrk_x(nplanes-ip+1)=
+      $              sgnl(icx)/mip(VIEW(icx),LADDER(icx))
+                dedxtrk_y(nplanes-ip+1)=
+      $              sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+             else
+                xm(nplanes-ip+1) = 0.
+                ym(nplanes-ip+1) = 0.
+                zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
+                xm_A(nplanes-ip+1) = xPAM_A
+                ym_A(nplanes-ip+1) = yPAM_A
+                xm_B(nplanes-ip+1) = xPAM_B
+                ym_B(nplanes-ip+1) = yPAM_B
+                xgood(nplanes-ip+1) = 0
+                ygood(nplanes-ip+1) = 0
+                resx(nplanes-ip+1) = 1000.!resxPAM
+                resy(nplanes-ip+1) = 1000.!resyPAM
+                dedxtrk_x(nplanes-ip+1)= 0
+                dedxtrk_y(nplanes-ip+1)= 0
+                if(icx.gt.0)then
+                   xgood(nplanes-ip+1) = 1
+                   resx(nplanes-ip+1) = resxPAM
+                   dedxtrk_x(nplanes-ip+1)=
+      $                 sgnl(icx)/mip(VIEW(icx),LADDER(icx))
+                elseif(icy.gt.0)then
+                   ygood(nplanes-ip+1) = 1
+                   resy(nplanes-ip+1) = resyPAM
+                   dedxtrk_y(nplanes-ip+1)=
+      $                 sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+                endif
+             endif
  *     |||||||||||||||||||||||||||||||||||||||||||||||||
  *     -------------------------------------------------
-Line 3114 
 c     $           AYV_STORE(nplanes-ip+1
+Line 3339 
 c     $           AYV_STORE(nplanes-ip+1
              xgood(nplanes-ip+1)=0
              ygood(nplanes-ip+1)=0
+             CP_STORE(nplanes-ip+1,ibest)=0 !re-init
+             CLS_STORE(nplanes-ip+1,ibest)=0
  *     --------------------------------------------------------------
  *     determine which ladder and sensor are intersected by the track
-Line 3136 
 c     $           AYV_STORE(nplanes-ip+1
+Line 3365 
 c     $           AYV_STORE(nplanes-ip+1
  *     ===========================================
  *     STEP 1 >>>>>>>  try to include a new couple
  *     ===========================================
- c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
+             distmin=100000000.
-             distmin=1000000.
              xmm = 0.
              ymm = 0.
              zmm = 0.
-Line 3150 
 c            if(DEBUG.EQ.1)print*,'>>>>
+Line 3378 
 c            if(DEBUG.EQ.1)print*,'>>>>
              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
-Line 3168 
 c     $              cl_used(icy).eq.1.o
+Line 3397 
 c     $              cl_used(icy).eq.1.o
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
                 id=id_cp(ip,icp,ist)
-                if(DEBUG.EQ.1)print*,'( couple ',id
+                if(DEBUG.EQ.1)
+      $              print*,'( couple ',id
       $              ,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm = xPAM
-Line 3180 
 c               distance = distance / RC
+Line 3410 
 c               distance = distance / RC
                    idm = id
                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
- c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
+                   clincnewc=10.*dsqrt(rymm**2+rxmm**2
-                   clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
+      $            +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))))! EM GCC4.7
-      $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
                 endif
          continue
              enddo               !end loop on couples on plane icp
- c            if(distmin.le.clinc)then     !QUIQUI
+             if(distmin.le.clincnewc)then
-             if(distmin.le.clincnewc)then     !QUIQUI
  *              -----------------------------------
                 xm(nplanes-ip+1) = xmm !<<<
                 ym(nplanes-ip+1) = ymm !<<<
-Line 3201 
 c            if(distmin.le.clinc)then
+Line 3429 
 c            if(distmin.le.clinc)then
  *              -----------------------------------
                 CP_STORE(nplanes-ip+1,ibest)=idm
                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
-      $              ,' (dist.= ',distmin,', cut ',clinc,' )'
+      $              ,' (dist.= ',distmin,', cut ',clincnewc,' )'
                 goto 133         !next plane
              endif
  *     ================================================
  *     STEP 2 >>>>>>>  try to include a single cluster
  *                     either from a couple or single
  *     ================================================
- c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
              distmin=1000000.
              xmm_A = 0.          !---------------------------
              ymm_A = 0.          ! init variables that
-Line 3227 
 c            if(DEBUG.EQ.1)print*,'>>>>
+Line 3454 
 c            if(DEBUG.EQ.1)print*,'>>>>
              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 3245 
 c     $              AXV_STORE(nplanes-i
+Line 3473 
 c     $              AXV_STORE(nplanes-i
       $              )
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG.EQ.1)print*,'( cl-X ',icx
+                if(DEBUG.EQ.1)
+      $              print*,'( cl-X ',icx
       $              ,' in cp ',id,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
-Line 3278 
 c     $              0.,AYV_STORE(nplane
+Line 3507 
 c     $              0.,AYV_STORE(nplane
       $              )
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG.EQ.1)print*,'( cl-Y ',icy
+                if(DEBUG.EQ.1)
+      $              print*,'( cl-Y ',icy
       $              ,' in cp ',id,' ) distance ',distance
                 if(distance.lt.distmin)then
                    xmm_A = xPAM_A
-Line 3298 
 c                 dedxmm = sgnl(icy)  !(
+Line 3528 
 c                 dedxmm = sgnl(icy)  !(
          continue
              enddo               !end loop on cluster inside couples
  *----- single clusters -----------------------------------------------
- c            print*,'## ncls(',ip,') ',ncls(ip)
              do ic=1,ncls(ip)    !loop on single clusters
                 icl=cls(ip,ic)
- c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
       $              .false.)goto 18882      !jump to the next singlet
-Line 3323 
 c               if(cl_used(icl).eq.1.or.
+Line 3551 
 c               if(cl_used(icl).eq.1.or.
                 distance = distance_to(XP,YP)
  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
-                if(DEBUG.EQ.1)print*,'( cl-s ',icl
+                if(DEBUG.EQ.1)
+      $              print*,'( cl-s ',icl
       $              ,' ) distance ',distance
                 if(distance.lt.distmin)then
- c                  if(DEBUG.EQ.1)print*,'YES'
                    xmm_A = xPAM_A
                    ymm_A = yPAM_A
                    zmm_A = zPAM_A
-Line 3347 
 c                  if(DEBUG.EQ.1)print*,
+Line 3575 
 c                  if(DEBUG.EQ.1)print*,
                 endif
          continue
              enddo               !end loop on single clusters
- c            print*,'## distmin ', distmin,' clinc ',clinc
- c     QUIQUI------------
- c     anche qui: non ci vuole clinc???
              if(iclm.ne.0)then
                 if(mod(VIEW(iclm),2).eq.0)then
                    clincnew=
-      $                 20*
+      $            20.*     !EM GCC4.7
-      $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
+      $            dsqrt(rxmm**2 +
+      $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1))
                 else if(mod(VIEW(iclm),2).ne.0)then
                    clincnew=
-      $                 10*
+      $            10.* !EM GCC4.7
-      $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
+      $            dsqrt(rymm**2 +
+      $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2))
                 endif
- c     QUIQUI------------
+                if(distmin.le.clincnew)then
-                if(distmin.le.clincnew)then   !QUIQUI
- c     if(distmin.le.clinc)then          !QUIQUI
                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
  *     ----------------------------
- c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
                    if(mod(VIEW(iclm),2).eq.0)then
                       XGOOD(nplanes-ip+1)=1.
                       resx(nplanes-ip+1)=rxmm
-                      if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
+                      if(DEBUG.EQ.1)
+      $                    print*,'%%%% included X-cl ',iclm
       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
       $                    ,', dist.= ',distmin
-      $                    ,', cut ',clinc,' )'
+      $                    ,', cut ',clincnew,' )'
                    else
                       YGOOD(nplanes-ip+1)=1.
                       resy(nplanes-ip+1)=rymm
-                      if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
+                      if(DEBUG.EQ.1)
+      $                    print*,'%%%% included Y-cl ',iclm
       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
       $                    ,', dist.= ', distmin
-      $                    ,', cut ',clinc,' )'
+      $                    ,', cut ',clincnew,' )'
                    endif
- c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
  *     ----------------------------
                    xm_A(nplanes-ip+1) = xmm_A
                    ym_A(nplanes-ip+1) = ymm_A
-Line 3405 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3630 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
        return
        end
  ***************************************************
  *                                                 *
  *                                                 *
-Line 3414 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
+Line 3640 
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
  *                                                 *
  **************************************************
  *
-       subroutine clean_XYclouds(ibest,iflag)
-       include 'commontracker.f'
-       include 'level1.f'
-       include 'common_momanhough.f'
-       include 'level2.f'
-       if(DEBUG.EQ.1)print*,'clean_XYclouds:'
-       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)=ntrk  !tag used clusters
- c$$$               cl_used(icly)=ntrk  !tag used clusters
-             elseif(icl.ne.0)then
- c$$$               cl_used(icl)=ntrk   !tag used clusters
-             endif
- *     -----------------------------
- *     remove the couple from clouds
- *     remove also vitual couples containing the
- *     selected clusters
- *     -----------------------------
-             do icp=1,ncp_plane(ip)
-                if(
-      $              clx(ip,icp).eq.iclx
-      $              .or.
-      $              clx(ip,icp).eq.icl
-      $              .or.
-      $              cly(ip,icp).eq.icly
-      $              .or.
-      $              cly(ip,icp).eq.icl
-      $              )then
-                   id=id_cp(ip,icp,1)
-                   if(DEBUG.EQ.1)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 3559 
 c$$$               cl_used(icl)=ntrk   !
+Line 3713 
 c$$$               cl_used(icl)=ntrk   !
          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 3677 
 c$$$               cl_used(icl)=ntrk   !
+Line 3835 
 c$$$               cl_used(icl)=ntrk   !
        character*10 PFA
        common/FINALPFA/PFA
-       real sinth,phi,pig
+       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 3728 
 c$$$               cl_used(icl)=ntrk   !
+Line 3888 
 c$$$               cl_used(icl)=ntrk   !
           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)
+ ccc         print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
           ax   = axv_nt(ip,ntr)
           ay   = ayv_nt(ip,ntr)
           bfx  = BX_STORE(ip,IDCAND)
-Line 3743 
 c$$$         angy    = 180.*atan(tgtemp)
+Line 3906 
 c$$$         angy    = 180.*atan(tgtemp)
           angy = effectiveangle(ay,2*ip-1,bfx)
- c         print*,'* ',ip,bfx,bfy,angx,angy
           id  = CP_STORE(ip,IDCAND) ! couple id
           icl = CLS_STORE(ip,IDCAND)
-Line 3754 
 c         print*,'* ',ip,bfx,bfy,angx,an
+Line 3916 
 c         print*,'* ',ip,bfx,bfy,angx,an
           if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
           LS(IP,ntr)      = ssensor+10*sladder
-          if(id.ne.0)then
+ c         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))
-             cl_used(cltrx(ip,ntr)) = 1     !tag used clusters
+             if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
-             cl_used(cltry(ip,ntr)) = 1     !tag used clusters
-             xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
+                cl_used(cltrx(ip,ntr)) = 1 !tag used clusters
-             ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
+                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
-             if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
+                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
-      $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
+      $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
-             if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
-      $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
+                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
-             multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
+             endif
-      $                         +10000*mult(cltrx(ip,ntr))
+             if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then
-             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
-             multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
+                cl_used(cltry(ip,ntr)) = 1 !tag used clusters
-      $                         +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
-          elseif(icl.ne.0)then
+                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
-Line 3828 
 c           >>> is a couple
+Line 4013 
 c           >>> is a couple
           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
- c$$$      print*,(xgood(i),i=1,6)
- c$$$      print*,(ygood(i),i=1,6)
- c$$$      print*,(ls(i,ntr),i=1,6)
- c$$$      print*,(dedx_x(i,ntr),i=1,6)
- c$$$      print*,(dedx_y(i,ntr),i=1,6)
- c$$$      print*,'-----------------------'
        end
        subroutine fill_level2_siglets
-Line 3872 
 c         if( mask_view(iv).ne.0 )good2(
+Line 4054 
 c         if( mask_view(iv).ne.0 )good2(
           ip=nplanes-npl(VIEW(icl))+1
           if(cl_used(icl).eq.0)then !cluster not included in any track
              if(mod(VIEW(icl),2).eq.0)then !=== X views
                 nclsx = nclsx + 1
                 planex(nclsx) = ip
                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
                 clsx(nclsx)   = icl
+                sxbad(nclsx)  = nbadstrips(1,icl)
+                multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
                 do is=1,2
  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
-                   call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
+                   call xyz_PAM(icl,0,is,PFAdef,'    ',0.,0.,0.,0.)
-                   xs(is,nclsx) = (xPAM_A+xPAM_B)/2
+                   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) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
                 clsy(nclsy)   = icl
+                sybad(nclsy)  = nbadstrips(1,icl)
+                multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
                 do is=1,2
  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
-                   call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
+                   call xyz_PAM(0,icl,is,'    ',PFAdef,0.,0.,0.,0.)
-                   ys(is,nclsy) = (yPAM_A+yPAM_B)/2
+                   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
-Line 3980 
 c$$$               print*,'ys(2,nclsy)
+Line 4162 
 c$$$               print*,'ys(2,nclsy)
                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)
                 nnn=nnn+ptcloud_yz(iyz)
              enddo
-             do ipt=1,nnn
+             do ipt=1,min(ndblt_max_nt,nnn)
                 db_cloud_nt(ipt)=db_cloud(ipt)
               enddo
           endif
-Line 3993 
 c$$$               print*,'ys(2,nclsy)
+Line 4175 
 c$$$               print*,'ys(2,nclsy)
                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
                 nnn=nnn+ptcloud_xz(ixz)
              enddo
-             do ipt=1,nnn
+             do ipt=1,min(ntrpt_max_nt,nnn)
                tr_cloud_nt(ipt)=tr_cloud(ipt)
               enddo
           endif

 Legend:



Removed from v.1.33
 


changed lines


 
Added in v.1.47
 Legend:



Removed from v.1.33
 


changed lines


 
Added in v.1.47
-Removed from v.1.33
+Added in v.1.47

	ViewVC Help
Powered by ViewVC 1.1.23