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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.31 by pam-fi, Fri Aug 31 14:56:52 2007 UTC revision 1.45 by mocchiut, Thu Jan 16 15:29:50 2014 UTC
# Line 20  Line 20 
20        include 'calib.f'        include 'calib.f'
21        include 'level2.f'        include 'level2.f'
22    
   
   
23  c      print*,'======================================================'  c      print*,'======================================================'
24  c$$$      do ic=1,NCLSTR1  c$$$      do ic=1,NCLSTR1
25  c$$$         if(.false.  c$$$         if(.false.
# Line 74  c$$$      enddo Line 72  c$$$      enddo
72  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
73  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
74    
   
75        call cl_to_couples(iflag)        call cl_to_couples(iflag)
76        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
77           goto 880               !go to next event           goto 880               !go to next event
78        endif        endif
79                if(ncp_tot.eq.0)goto 880  !go to next event    
80  *-----------------------------------------------------  *-----------------------------------------------------
81  *-----------------------------------------------------  *-----------------------------------------------------
82  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 113  c$$$      enddo Line 110  c$$$      enddo
110        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
111           goto 880               !go to next event                       goto 880               !go to next event            
112        endif        endif
113          if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event    
114                
115                
116  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 163  c$$$      enddo Line 161  c$$$      enddo
161        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
162           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
163        endif        endif
164    *     ------------------------------------------------
165    *     first try to release the tolerance
166    *     ------------------------------------------------
167        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
168          if(cutdistyz.lt.maxcuty/2)then           if(cutdistyz.lt.maxcuty/2)then
169            cutdistyz=cutdistyz+cutystep              cutdistyz=cutdistyz+cutystep
170          else           else
171            cutdistyz=cutdistyz+(3*cutystep)              cutdistyz=cutdistyz+(3*cutystep)
172          endif           endif
173          goto 878                           if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
174        endif                               goto 878                
175          endif    
176    *     ------------------------------------------------
177    *     hence reduce the minimum number of plane
178    *     ------------------------------------------------
179          if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then
180             nplyz_min=nplyz_min-1
181             if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min
182             goto 878
183          endif
184    
185          if(nclouds_yz.eq.0)goto 880 !go to next event    
186    
187    
188        if(planehit.eq.3) goto 881            ccc   if(planehit.eq.3) goto 881    
189                
190   879  continue     879  continue  
191        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
192        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
193           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
194        endif        endif
195                              *     ------------------------------------------------
196    *     first try to release the tolerance
197    *     ------------------------------------------------                          
198        if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then        if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
199          cutdistxz=cutdistxz+cutxstep          cutdistxz=cutdistxz+cutxstep
200             if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz
201          goto 879                          goto 879                
202        endif                            endif                    
203    *     ------------------------------------------------
204    *     hence reduce the minimum number of plane
205    *     ------------------------------------------------
206          if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then
207             nplxz_min=nplxz_min-1
208             if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min
209             goto 879
210          endif
211    
212          if(nclouds_xz.eq.0)goto 880 !go to next event    
213    
214            
215   881  continue    c$$$ 881  continue  
216  *     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
217        if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.  c$$$      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
218       $     nplxz_min.ne.y_min_start)then  c$$$     $     nplxz_min.ne.y_min_start)then
219          nptxz_min=x_min_step      c$$$        nptxz_min=x_min_step    
220          nplxz_min=x_min_start-x_min_step                c$$$        nplxz_min=x_min_start-x_min_step              
221          goto 879                  c$$$        goto 879                
222        endif                      c$$$      endif                    
223                    
224   880  return   880  return
225        end        end
# Line 243  c      include 'momanhough_init.f' Line 269  c      include 'momanhough_init.f'
269  *  *
270  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
271  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
272           ntrk=0                 !counter of identified physical tracks  ccc         ntrk=0                 !counter of identified physical tracks
273    
274    c11111    continue               !<<<<<<< come here when performing a new search
275             continue                  !<<<<<<< come here when performing a new search
276    
277  11111    continue               !<<<<<<< come here when performing a new search           if(nclouds_xz.eq.0)goto 880 !go to next event    
278             if(nclouds_yz.eq.0)goto 880 !go to next event    
279    
280  c         iflag=0  c         iflag=0
281           call clouds_to_ctrack(iflag)           call clouds_to_ctrack(iflag)
282           if(iflag.eq.1)then     !no candidate tracks found           if(iflag.eq.1)then     !no candidate tracks found
283              goto 880            !fill ntp and go to next event                goto 880            !fill ntp and go to next event  
284           endif           endif
285             if(ntracks.eq.0)goto 880 !go to next event    
286    
287           FIMAGE=.false.         !processing best track (not track image)           FIMAGE=.false.         !processing best track (not track image)
288           ibest=0                !best track among candidates           ibest=0                !best track among candidates
# Line 272  c$$$         if(ibest.eq.0)goto 880 !>> Line 303  c$$$         if(ibest.eq.0)goto 880 !>>
303  *     1st) decreasing n.points  *     1st) decreasing n.points
304  *     2nd) increasing chi**2  *     2nd) increasing chi**2
305  *     -------------------------------------------------------  *     -------------------------------------------------------
306           rchi2best=1000000000.           rchi2best=1000000000.
307           ndofbest=0                       ndofbest=0            
308           do i=1,ntracks           do i=1,ntracks
309             ndof=0                           ndof=0              
# Line 295  c$$$         if(ibest.eq.0)goto 880 !>> Line 326  c$$$         if(ibest.eq.0)goto 880 !>>
326             endif             endif
327           enddo           enddo
328    
 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  
329    
330           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
331  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
# Line 350  c$$$         enddo Line 363  c$$$         enddo
363           do i=1,5           do i=1,5
364              AL_GUESS(i)=AL(i)              AL_GUESS(i)=AL(i)
365           enddo           enddo
 c         print*,'## guess: ',al  
366    
367           do i=1,5           do i=1,5
368              AL(i)=dble(AL_STORE(i,icand))                          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
377           if(VERBOSE.EQ.1)iprint=1           if(VERBOSE.EQ.1)iprint=1
378           if(DEBUG.EQ.1)iprint=2           if(DEBUG.EQ.1)iprint=2
379           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
380           if(ifail.ne.0) then  cc         if(ifail.ne.0) then
381             if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
382                if(CHI2.ne.CHI2)CHI2=-9999. !new
383              if(VERBOSE.EQ.1)then              if(VERBOSE.EQ.1)then
384                 print *,                 print *,
385       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
386       $              ,iev       $              ,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*,'----------------------------------------------'  
387              endif              endif
 c            chi2=-chi2  
388           endif           endif
389                    
390           if(DEBUG.EQ.1)then           if(DEBUG.EQ.1)then
# Line 424  c         rchi2=chi2/dble(ndof) Line 430  c         rchi2=chi2/dble(ndof)
430              if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'              if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
431              goto 122            !jump              goto 122            !jump
432           endif           endif
 c         print*,'is1 ',is1  
433           do ip=1,NPLANES           do ip=1,NPLANES
434              is2 = SENSOR_STORE(ip,icand) !sensor              is2 = SENSOR_STORE(ip,icand) !sensor
 c            print*,'is2 ',is2,' ip ',ip  
435              if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted              if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
436              if(              if(
437       $           (is1.ne.is2.and.is2.ne.0)       $           (is1.ne.is2.and.is2.ne.0)
438       $           )goto 122      !jump (this track cannot have an image)       $           )goto 122      !jump (this track cannot have an image)
439           enddo           enddo
440           if(DEBUG.eq.1)print*,' >>> ambiguous track! '           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  
441  *     ---------------------------------------------------------------  *     ---------------------------------------------------------------
442  *     take the candidate with the greatest number of matching couples  *     take the candidate with the greatest number of matching couples
443  *     if more than one satisfies the requirement,  *     if more than one satisfies the requirement,
# Line 486  c$$$         enddo Line 469  c$$$         enddo
469                                            
470                    endif                    endif
471                 endif                 endif
 c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)  
 c$$$     $              ,CP_STORE(nplanes-ip+1,icand),ncp  
472              enddo              enddo
473   117        continue   117        continue
474              trackimage(i)=ncp   !number of matching couples              trackimage(i)=ncp   !number of matching couples
# Line 501  c$$$     $              ,CP_STORE(nplane Line 482  c$$$     $              ,CP_STORE(nplane
482           do i=1,ntracks           do i=1,ntracks
483              if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1              if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
484           enddo           enddo
485             if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood
486  *     if there are, choose the best one  *     if there are, choose the best one
487           if(ngood.gt.1)then           if(ngood.gt.0)then
488  *     -------------------------------------------------------  *     -------------------------------------------------------
489  *     order track-candidates according to:  *     order track-candidates according to:
490  *     1st) decreasing n.points  *     1st) decreasing n.points
# Line 532  c$$$     $              ,CP_STORE(nplane Line 514  c$$$     $              ,CP_STORE(nplane
514                    endif                    endif
515                 endif                 endif
516              enddo              enddo
517    
518                if(DEBUG.EQ.1)then
519                   print*,'Track candidate ',iimage
520         $              ,' >>> TRACK IMAGE >>> of'
521         $              ,ibest
522                endif
523                            
524           endif           endif
525    
          if(DEBUG.EQ.1)print*,'Track candidate ',iimage  
      $        ,' >>> TRACK IMAGE >>> of'  
      $        ,ibest  
526    
527   122     continue   122     continue
528    
529    
530  *     --- and store the results --------------------------------  *     --- 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)  
   
531           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
532              if(verbose.eq.1)              if(verbose.eq.1)
533       $           print*,       $           print*,
# Line 562  c     $        ,iimage,fimage,ntrk,image Line 537  c     $        ,iimage,fimage,ntrk,image
537  cc            good2=.false.  cc            good2=.false.
538              goto 880            !fill ntp and go to next event              goto 880            !fill ntp and go to next event
539           endif           endif
540    
541             ntrk = ntrk + 1                   !counter of found tracks
542             if(.not.FIMAGE
543         $        .and.iimage.eq.0) image(ntrk)= 0
544             if(.not.FIMAGE
545         $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
546             if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
547             call fill_level2_tracks(ntrk)     !==> good2=.true.
548    
549    c$$$         if(ntrk.eq.NTRKMAX)then
550    c$$$            if(verbose.eq.1)
551    c$$$     $           print*,
552    c$$$     $           '** warning ** number of identified '//
553    c$$$     $           'tracks exceeds vector dimension '
554    c$$$     $           ,'( ',NTRKMAX,' )'
555    c$$$cc            good2=.false.
556    c$$$            goto 880            !fill ntp and go to next event
557    c$$$         endif
558           if(iimage.ne.0)then           if(iimage.ne.0)then
559              FIMAGE=.true.       !              FIMAGE=.true.       !
560              goto 1212           !>>> fit image-track              goto 1212           !>>> fit image-track
561           endif           endif
562    
 *     --- 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  
 *     **********************************************  
   
   
563    
564   880     return   880     return
565        end        end
# Line 688  c     $        rchi2best.lt.15..and. Line 649  c     $        rchi2best.lt.15..and.
649    
650        real stripx,stripy        real stripx,stripy
651    
652          double precision xi,yi,zi
653          double precision xi_A,yi_A,zi_A
654          double precision xi_B,yi_B,zi_B
655        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
656        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
657        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
# Line 696  c     $        rchi2best.lt.15..and. Line 660  c     $        rchi2best.lt.15..and.
660        parameter (ndivx=30)        parameter (ndivx=30)
661    
662    
 c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy  
663                
664        resxPAM = 0        resxPAM = 0
665        resyPAM = 0        resyPAM = 0
666    
667        xPAM = 0.        xPAM = 0.D0
668        yPAM = 0.        yPAM = 0.D0
669        zPAM = 0.        zPAM = 0.D0
670        xPAM_A = 0.        xPAM_A = 0.D0
671        yPAM_A = 0.        yPAM_A = 0.D0
672        zPAM_A = 0.        zPAM_A = 0.D0
673        xPAM_B = 0.        xPAM_B = 0.D0
674        yPAM_B = 0.        yPAM_B = 0.D0
675        zPAM_B = 0.        zPAM_B = 0.D0
 c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy  
676    
677        if(sensor.lt.1.or.sensor.gt.2)then        if(sensor.lt.1.or.sensor.gt.2)then
678           print*,'xyz_PAM   ***ERROR*** wrong input '           print*,'xyz_PAM   ***ERROR*** wrong input '
# Line 754  c         resxPAM = RESXAV Line 716  c         resxPAM = RESXAV
716           stripx  = stripx + corr           stripx  = stripx + corr
717           resxPAM = res           resxPAM = res
718    
719   10   endif   10   continue
720          endif
721            
722  *     -----------------  *     -----------------
723  *     CLUSTER Y  *     CLUSTER Y
# Line 800  c         resyPAM = RESYAV Line 763  c         resyPAM = RESYAV
763           stripy  = stripy + corr           stripy  = stripy + corr
764           resyPAM = res           resyPAM = res
765    
766   20   endif   20   continue
767          endif
768    
 c$$$      print*,'## stripx,stripy ',stripx,stripy  
769    
770  c===========================================================  c===========================================================
771  C     COUPLE  C     COUPLE
# Line 819  c--------------------------------------- Line 782  c---------------------------------------
782       $              ' WARNING: false X strip: strip ',stripx       $              ' WARNING: false X strip: strip ',stripx
783              endif              endif
784           endif           endif
785           xi = acoordsi(stripx,viewx)           xi = dcoordsi(stripx,viewx)
786           yi = acoordsi(stripy,viewy)           yi = dcoordsi(stripy,viewy)
787           zi = 0.           zi = 0.D0
788                    
789  c------------------------------------------------------------------------  c------------------------------------------------------------------------
790  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
# Line 856  c--------------------------------------- Line 819  c---------------------------------------
819           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
820           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
821    
822           xPAM_A = 0.           xPAM_A = 0.D0
823           yPAM_A = 0.           yPAM_A = 0.D0
824           zPAM_A = 0.           zPAM_A = 0.D0
825    
826           xPAM_B = 0.           xPAM_B = 0.D0
827           yPAM_B = 0.           yPAM_B = 0.D0
828           zPAM_B = 0.           zPAM_B = 0.D0
829    
830        elseif(        elseif(
831       $        (icx.ne.0.and.icy.eq.0).or.       $        (icx.ne.0.and.icy.eq.0).or.
# Line 882  C======================================= Line 845  C=======================================
845              nldx = nldy              nldx = nldy
846              viewx = viewy + 1              viewx = viewy + 1
847    
848              yi   = acoordsi(stripy,viewy)              xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center
849                yi = dcoordsi(stripy,viewy)
850                zi = 0.D0
851    
852              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
853              yi_A = yi              yi_A = yi
# Line 892  C======================================= Line 857  C=======================================
857              yi_B = yi              yi_B = yi
858              zi_B = 0.              zi_B = 0.
859    
 c            print*,'Y-cl ',icy,stripy,' --> ',yi  
 c            print*,xi_A,' <--> ',xi_B  
860                            
861           elseif(icx.ne.0)then           elseif(icx.ne.0)then
862  c===========================================================  c===========================================================
# Line 904  C======================================= Line 867  C=======================================
867              nldy = nldx              nldy = nldx
868              viewy = viewx - 1              viewy = viewx - 1
869    
 c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx  
 c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  
870              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
871       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
872                 if(DEBUG.EQ.1) then                 if(DEBUG.EQ.1) then
# Line 913  c            if((stripx.le.3).or.(stripx Line 874  c            if((stripx.le.3).or.(stripx
874       $                 ' WARNING: false X strip: strip ',stripx       $                 ' WARNING: false X strip: strip ',stripx
875                 endif                 endif
876              endif              endif
877              xi   = acoordsi(stripx,viewx)  
878                xi = dcoordsi(stripx,viewx)
879                yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
880                zi = 0.D0
881    
882              xi_A = xi              xi_A = xi
883              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 929  c            if((stripx.le.3).or.(stripx Line 893  c            if((stripx.le.3).or.(stripx
893                 yi_B = yi                 yi_B = yi
894              endif              endif
895    
 c            print*,'X-cl ',icx,stripx,' --> ',xi  
 c            print*,yi_A,' <--> ',yi_B  
896    
897           else           else
898              if(DEBUG.EQ.1) then              if(DEBUG.EQ.1) then
# Line 976  c     N.B. I convert angles from microra Line 938  c     N.B. I convert angles from microra
938       $        + zi_B       $        + zi_B
939       $        + dz(nplx,nldx,sensor)       $        + dz(nplx,nldx,sensor)
940    
941    
942    
943             xrt = xi
944         $        - omega(nplx,nldx,sensor)*yi
945         $        + gamma(nplx,nldx,sensor)*zi
946         $        + dx(nplx,nldx,sensor)
947            
948             yrt = omega(nplx,nldx,sensor)*xi
949         $        + yi
950         $        - beta(nplx,nldx,sensor)*zi
951         $        + dy(nplx,nldx,sensor)
952            
953             zrt = -gamma(nplx,nldx,sensor)*xi
954         $        + beta(nplx,nldx,sensor)*yi
955         $        + zi
956         $        + dz(nplx,nldx,sensor)
957    
958    
959                    
960  c      xrt = xi  c      xrt = xi
961  c      yrt = yi  c      yrt = yi
# Line 986  c     (xPAM,yPAM,zPAM) = measured coordi Line 966  c     (xPAM,yPAM,zPAM) = measured coordi
966  c                        in PAMELA reference system  c                        in PAMELA reference system
967  c------------------------------------------------------------------------  c------------------------------------------------------------------------
968    
969           xPAM = 0.           xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
970           yPAM = 0.           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
971           zPAM = 0.           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
972    c$$$         xPAM = 0.D0
973    c$$$         yPAM = 0.D0
974    c$$$         zPAM = 0.D0
975    
976           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
977           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 999  c--------------------------------------- Line 982  c---------------------------------------
982           zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4           zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
983                    
984    
 c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  
985    
986        else        else
987           if(DEBUG.EQ.1) then           if(DEBUG.EQ.1) then
# Line 1010  c         print*,'A-(',xPAM_A,yPAM_A,') Line 992  c         print*,'A-(',xPAM_A,yPAM_A,')
992        endif        endif
993                    
994    
 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  
995    
996   100  continue   100  continue
997        end        end
# Line 1050  c$$$      PFAy = 'COG4'!PFA Line 1029  c$$$      PFAy = 'COG4'!PFA
1029    
1030        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1031              print*,'xyzpam: ***WARNING*** clusters ',icx,icy              print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1032       $           ,' does not exists (nclstr1=',nclstr1,')'       $           ,' do not exists (n.clusters=',nclstr1,')'
1033              icx = -1*icx              icx = -1*icx
1034              icy = -1*icy              icy = -1*icy
1035              return              return
# Line 1060  c$$$      PFAy = 'COG4'!PFA Line 1039  c$$$      PFAy = 'COG4'!PFA
1039        call idtoc(pfaid,PFAx)        call idtoc(pfaid,PFAx)
1040        call idtoc(pfaid,PFAy)        call idtoc(pfaid,PFAy)
1041    
 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  
1042                
1043        if(icx.ne.0.and.icy.ne.0)then        if(icx.ne.0.and.icy.ne.0)then
1044    
1045           ipx=npl(VIEW(icx))           ipx=npl(VIEW(icx))
1046           ipy=npl(VIEW(icy))           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  
1047                    
1048           if( (nplanes-ipx+1).ne.ip )then           if( (nplanes-ipx+1).ne.ip )then
1049              print*,'xyzpam: ***WARNING*** cluster ',icx              print*,'xyzpam: ***WARNING*** cluster ',icx
# Line 1097  c$$$     $        ,' does not belong to Line 1068  c$$$     $        ,' does not belong to
1068           xm(ip) = xPAM           xm(ip) = xPAM
1069           ym(ip) = yPAM           ym(ip) = yPAM
1070           zm(ip) = zPAM           zm(ip) = zPAM
1071           xm_A(ip) = 0.           xm_A(ip) = 0.D0
1072           ym_A(ip) = 0.           ym_A(ip) = 0.D0
1073           xm_B(ip) = 0.           xm_B(ip) = 0.D0
1074           ym_B(ip) = 0.           ym_B(ip) = 0.D0
1075    
1076  c         zv(ip) = zPAM  c         zv(ip) = zPAM
1077    
1078        elseif(icx.eq.0.and.icy.ne.0)then        elseif(icx.eq.0.and.icy.ne.0)then
1079    
1080           ipy=npl(VIEW(icy))           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  
1081           if( (nplanes-ipy+1).ne.ip )then           if( (nplanes-ipy+1).ne.ip )then
1082              print*,'xyzpam: ***WARNING*** cluster ',icy              print*,'xyzpam: ***WARNING*** cluster ',icy
1083       $           ,' does not belong to plane: ',ip       $           ,' does not belong to plane: ',ip
# Line 1124  c$$$     $        ,' does not belong to Line 1092  c$$$     $        ,' does not belong to
1092           resx(ip)  = 1000.           resx(ip)  = 1000.
1093           resy(ip)  = resyPAM           resy(ip)  = resyPAM
1094    
1095           xm(ip) = -100.  c$$$         xm(ip) = -100.
1096           ym(ip) = -100.  c$$$         ym(ip) = -100.
1097           zm(ip) = (zPAM_A+zPAM_B)/2.  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1098             xm(ip) = xPAM
1099             ym(ip) = yPAM
1100             zm(ip) = zPAM
1101           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1102           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1103           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
# Line 1137  c         zv(ip) = (zPAM_A+zPAM_B)/2. Line 1108  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1108        elseif(icx.ne.0.and.icy.eq.0)then        elseif(icx.ne.0.and.icy.eq.0)then
1109    
1110           ipx=npl(VIEW(icx))           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  
1111    
1112           if( (nplanes-ipx+1).ne.ip )then           if( (nplanes-ipx+1).ne.ip )then
1113              print*,'xyzpam: ***WARNING*** cluster ',icx              print*,'xyzpam: ***WARNING*** cluster ',icx
# Line 1155  c$$$     $        ,' does not belong to Line 1123  c$$$     $        ,' does not belong to
1123           resx(ip)  = resxPAM           resx(ip)  = resxPAM
1124           resy(ip)  = 1000.           resy(ip)  = 1000.
1125    
1126           xm(ip) = -100.  c$$$         xm(ip) = -100.
1127           ym(ip) = -100.  c$$$         ym(ip) = -100.
1128           zm(ip) = (zPAM_A+zPAM_B)/2.  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1129             xm(ip) = xPAM
1130             ym(ip) = yPAM
1131             zm(ip) = zPAM
1132           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1133           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1134           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
# Line 1171  c         zv(ip) = (zPAM_A+zPAM_B)/2. Line 1142  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1142           if(lad.ne.0)il=lad           if(lad.ne.0)il=lad
1143           is = 1           is = 1
1144           if(sensor.ne.0)is=sensor           if(sensor.ne.0)is=sensor
 c         print*,nplanes-ip+1,il,is  
1145    
1146           xgood(ip) = 0.           xgood(ip) = 0.
1147           ygood(ip) = 0.           ygood(ip) = 0.
# Line 1191  c         zv(ip) = z_mech_sensor(nplanes Line 1161  c         zv(ip) = z_mech_sensor(nplanes
1161        endif        endif
1162    
1163        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
 c         print*,'----------------------------- track coord'  
1164  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1165           write(*,22222)ip,zm(ip),xm(ip),ym(ip)           write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1166       $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)       $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1167       $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)       $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
 c$$$         print*,'-----------------------------'  
1168        endif        endif
1169        end        end
1170    
# Line 1220  c$$$         print*,'------------------- Line 1188  c$$$         print*,'-------------------
1188  *      *    
1189  ********************************************************************************  ********************************************************************************
1190    
1191        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1192    
1193        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1194    
# Line 1229  c$$$         print*,'------------------- Line 1197  c$$$         print*,'-------------------
1197  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1198  *     -----------------------------------  *     -----------------------------------
1199    
1200          real rXPP,rYPP
1201          double precision XPP,YPP
1202        double precision distance,RE        double precision distance,RE
1203        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1204    
1205          XPP=DBLE(rXPP)
1206          YPP=DBLE(rYPP)
1207    
1208  *     ----------------------  *     ----------------------
1209        if (        if (
1210       +     xPAM.eq.0.and.  c$$$     +     xPAM.eq.0.and.
1211       +     yPAM.eq.0.and.  c$$$     +     yPAM.eq.0.and.
1212       +     zPAM.eq.0.and.  c$$$     +     zPAM.eq.0.and.
1213       +     xPAM_A.ne.0.and.       +     xPAM_A.ne.0.and.
1214       +     yPAM_A.ne.0.and.       +     yPAM_A.ne.0.and.
1215       +     zPAM_A.ne.0.and.       +     zPAM_A.ne.0.and.
# Line 1278  c$$$         print*,'------------------- Line 1251  c$$$         print*,'-------------------
1251  cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2  cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1252           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1253    
 c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',re  
1254    
1255                    
1256  *     ----------------------  *     ----------------------
1257        elseif(        elseif(
1258       +     xPAM.ne.0.and.  c$$$     +     xPAM.ne.0.and.
1259       +     yPAM.ne.0.and.  c$$$     +     yPAM.ne.0.and.
1260       +     zPAM.ne.0.and.  c$$$     +     zPAM.ne.0.and.
1261       +     xPAM_A.eq.0.and.       +     xPAM_A.eq.0.and.
1262       +     yPAM_A.eq.0.and.       +     yPAM_A.eq.0.and.
1263       +     zPAM_A.eq.0.and.       +     zPAM_A.eq.0.and.
# Line 1308  c$$$     $        + Line 1278  c$$$     $        +
1278  c$$$     $        ((yPAM-YPP)/resyPAM)**2  c$$$     $        ((yPAM-YPP)/resyPAM)**2
1279           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1280    
 c$$$         print*,xPAM,yPAM,zPAM  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',resxPAM,resyPAM  
1281                    
1282        else        else
1283                    
 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  
1284        endif          endif  
1285    
1286        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1358  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1320  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1320        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1321        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1322        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1323        real*8 yvvv,xvvv        double precision yvvv,xvvv
1324        double precision xi,yi,zi        double precision xi,yi,zi
1325        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1326        real AA,BB        real AA,BB
1327        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1328    
1329  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1330        real ptoll        real ptoll
1331        data ptoll/150./          !um        data ptoll/150./          !um
1332    
1333        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1334    
1335        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1336        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1383  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1345  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1345  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1346  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1347  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1348                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 xi = dcoordsi(stripx,viewx)
1349       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...                 yi = dcoordsi(stripy,viewy)
1350  c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...                 zi = 0.D0
 c                  print*,'whichsensor: ',  
 c     $                ' WARNING: false X strip: strip ',stripx  
                endif  
                xi = acoordsi(stripx,viewx)  
                yi = acoordsi(stripy,viewy)  
                zi = 0.  
1351  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1352  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1353  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1416  c--------------------------------------- Line 1372  c---------------------------------------
1372    
1373                 yvv(iv)=sngl(yvvv)                 yvv(iv)=sngl(yvvv)
1374                 xvv(iv)=sngl(xvvv)                 xvv(iv)=sngl(xvvv)
 c               print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '  
 c     $              ,iv,xvv(iv),yvv(iv)  
1375              enddo               !end loop on sensor vertexes              enddo               !end loop on sensor vertexes
1376    
1377              dtot=0.              dtot=0.
# Line 1425  c     $              ,iv,xvv(iv),yvv(iv) Line 1379  c     $              ,iv,xvv(iv),yvv(iv)
1379                 iv1=iside                 iv1=iside
1380                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1381  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1382                 AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))                 AA = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7
1383                 BB = yvv(iv1) - AA*xvv(iv1)                 BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7
1384  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1385                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
1386                 yoo = AA*xoo + BB                 yoo = AA*xoo + BB
# Line 1438  c     $              ,iv,xvv(iv),yvv(iv) Line 1392  c     $              ,iv,xvv(iv),yvv(iv)
1392                 iv1=iside                 iv1=iside
1393                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1394  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1395                 AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))                 AA = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7
1396                 BB = xvv(iv1) - AA*yvv(iv1)                 BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7
1397  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1398                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
1399                 xoo = AA*yoo + BB                 xoo = AA*yoo + BB
# Line 1542  c      include 'common_analysis.f' Line 1496  c      include 'common_analysis.f'
1496        is_cp=0        is_cp=0
1497        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1498        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
 c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  
1499    
1500        return        return
1501        end        end
# Line 1641  c      include 'level1.f' Line 1594  c      include 'level1.f'
1594        integer iflag        integer iflag
1595    
1596        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1597        
1598          iflag = iflag
1599        if(DEBUG.EQ.1)print*,'cl_to_couples:'        if(DEBUG.EQ.1)print*,'cl_to_couples:'
1600    
1601    cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
1602    
1603  *     init variables  *     init variables
       ncp_tot=0  
1604        do ip=1,nplanes        do ip=1,nplanes
1605           do ico=1,ncouplemax           do ico=1,ncouplemax
1606              clx(ip,ico)=0              clx(ip,ico)=0
# Line 1658  c      include 'level1.f' Line 1613  c      include 'level1.f'
1613           ncls(ip)=0           ncls(ip)=0
1614        enddo        enddo
1615        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1616           cl_single(icl) = 1           cl_single(icl) = 1     !all are single
1617           cl_good(icl)   = 0           cl_good(icl)   = 0     !all are bad
1618        enddo        enddo
1619        do iv=1,nviews        do iv=1,nviews
1620           ncl_view(iv)  = 0           ncl_view(iv)  = 0
# Line 1687  c            mask_view(iv) = 1 Line 1642  c            mask_view(iv) = 1
1642        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1643           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1644                    
1645             if(cl_used(icx).ne.0)goto 10
1646    
1647  *     ----------------------------------------------------  *     ----------------------------------------------------
1648  *     jump masked views (X VIEW)  *     jump masked views (X VIEW)
1649  *     ----------------------------------------------------  *     ----------------------------------------------------
# Line 1694  c            mask_view(iv) = 1 Line 1651  c            mask_view(iv) = 1
1651  *     ----------------------------------------------------  *     ----------------------------------------------------
1652  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1653  *     ----------------------------------------------------  *     ----------------------------------------------------
1654           if(sgnl(icx).lt.dedx_x_min)then           if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
1655              cl_single(icx)=0              cl_single(icx)=0
1656              goto 10              goto 10
1657           endif           endif
# Line 1735  c     endif Line 1692  c     endif
1692                    
1693           do icy=1,nclstr1       !loop on cluster (Y)           do icy=1,nclstr1       !loop on cluster (Y)
1694              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1695    
1696                if(cl_used(icx).ne.0)goto 20
1697                            
1698  *     ----------------------------------------------------  *     ----------------------------------------------------
1699  *     jump masked views (Y VIEW)  *     jump masked views (Y VIEW)
# Line 1744  c     endif Line 1703  c     endif
1703  *     ----------------------------------------------------  *     ----------------------------------------------------
1704  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1705  *     ----------------------------------------------------  *     ----------------------------------------------------
1706              if(sgnl(icy).lt.dedx_y_min)then              if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
1707                 cl_single(icy)=0                 cl_single(icy)=0
1708                 goto 20                 goto 20
1709              endif              endif
# Line 1822  c                  cut = chcut * sch(npl Line 1781  c                  cut = chcut * sch(npl
1781       $                 ,'( ',ncouplemax,' ) --> masked!'       $                 ,'( ',ncouplemax,' ) --> masked!'
1782  c                  mask_view(nviewx(nplx)) = 2  c                  mask_view(nviewx(nplx)) = 2
1783  c                  mask_view(nviewy(nply)) = 2  c                  mask_view(nviewy(nply)) = 2
1784                    mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1785                    mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1786                    goto 10                    goto 10
1787                 endif                 endif
1788                                
# Line 1843  c                  mask_view(nviewy(nply Line 1802  c                  mask_view(nviewy(nply
1802   10      continue   10      continue
1803        enddo                     !end loop on clusters(X)        enddo                     !end loop on clusters(X)
1804                
         
1805        do icl=1,nclstr1        do icl=1,nclstr1
1806           if(cl_single(icl).eq.1)then           if(cl_single(icl).eq.1)then
1807              ip=npl(VIEW(icl))              ip=npl(VIEW(icl))
# Line 1851  c                  mask_view(nviewy(nply Line 1809  c                  mask_view(nviewy(nply
1809              cls(ip,ncls(ip))=icl              cls(ip,ncls(ip))=icl
1810           endif           endif
1811        enddo        enddo
1812    
1813    c 80   continue
1814          continue
1815                
1816                
1817        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
1818           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1819           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1820             print*,'used    ',(cl_used(i),i=1,nclstr1)
1821           print*,'singlets',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
1822           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1823        endif        endif
1824    
1825      
1826          if(.not.RECOVER_SINGLETS)goto 81
1827    
1828    *     ////////////////////////////////////////////////
1829    *     PATCH to recover events with less than 3 couples
1830    *     ////////////////////////////////////////////////    
1831    *     loop over singlet and create "fake" couples
1832    *     (with clx=0 or cly=0)
1833    *    
1834    
1835          if(DEBUG.EQ.1)
1836         $     print*,'>>> Recover singlets '
1837         $     ,'(creates fake couples) <<<'
1838          do icl=1,nclstr1
1839             if(
1840         $        cl_single(icl).eq.1.and.
1841         $        cl_used(icl).eq.0.and.
1842         $        .true.)then
1843    *     ----------------------------------------------------
1844    *     jump masked views
1845    *     ----------------------------------------------------
1846                if( mask_view(VIEW(icl)).ne.0 ) goto 21
1847                if(mod(VIEW(icl),2).eq.0)then !=== X views
1848    *     ----------------------------------------------------
1849    *     cut on charge (X VIEW)
1850    *     ----------------------------------------------------
1851                   if(sgnl(icl).lt.dedx_x_min) goto 21
1852                  
1853                   nplx=npl(VIEW(icl))
1854    *     ------------------> (FAKE) COUPLE <-----------------
1855                   ncp_plane(nplx) = ncp_plane(nplx) + 1
1856                   clx(nplx,ncp_plane(nplx))=icl
1857                   cly(nplx,ncp_plane(nplx))=0
1858    c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1859    *     ----------------------------------------------------
1860                  
1861                else                !=== Y views
1862    *     ----------------------------------------------------
1863    *     cut on charge (Y VIEW)
1864    *     ----------------------------------------------------
1865                   if(sgnl(icl).lt.dedx_y_min) goto 21
1866                  
1867                   nply=npl(VIEW(icl))
1868    *     ------------------> (FAKE) COUPLE <-----------------
1869                   ncp_plane(nply) = ncp_plane(nply) + 1
1870                   clx(nply,ncp_plane(nply))=0
1871                   cly(nply,ncp_plane(nply))=icl
1872    c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1873    *     ----------------------------------------------------
1874                  
1875                endif
1876             endif                  !end "single" condition
1877     21      continue
1878          enddo                     !end loop over clusters
1879    
1880          if(DEBUG.EQ.1)
1881         $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1882    
1883    
1884     81   continue
1885                
1886        do ip=1,6        ncp_tot=0
1887          do ip=1,NPLANES
1888           ncp_tot = ncp_tot + ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1889        enddo        enddo
1890          if(DEBUG.EQ.1)
1891         $     print*,'n.couple tot:      ',ncp_tot
1892                
1893        return        return
1894        end        end
# Line 1925  c      double precision xm3,ym3,zm3 Line 1951  c      double precision xm3,ym3,zm3
1951  *     --------------------------------------------  *     --------------------------------------------
1952        do ip=1,nplanes        do ip=1,nplanes
1953           if(ncp_plane(ip).gt.ncouplelimit)then           if(ncp_plane(ip).gt.ncouplelimit)then
 c            mask_view(nviewx(ip)) = 8  
 c            mask_view(nviewy(ip)) = 8  
1954              mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7              mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1955              mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7              mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1956           endif           endif
# Line 1937  c            mask_view(nviewy(ip)) = 8 Line 1961  c            mask_view(nviewy(ip)) = 8
1961        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1962                
1963        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1964    c$$$         print*,'(1) ip ',ip1
1965    c$$$     $        ,mask_view(nviewx(ip1))
1966    c$$$     $        ,mask_view(nviewy(ip1))        
1967           if(  mask_view(nviewx(ip1)).ne.0 .or.           if(  mask_view(nviewx(ip1)).ne.0 .or.
1968       $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane       $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1969           do is1=1,2             !loop on sensors - COPPIA 1                       do is1=1,2             !loop on sensors - COPPIA 1            
1970              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1971                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1972                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1973                  
1974    c$$$               print*,'(1) ip ',ip1,' icp ',icp1
1975    
1976  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1977  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1978                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1979                 xm1=xPAM                 xm1=REAL(xPAM) !EM GCC4.7
1980                 ym1=yPAM                 ym1=REAL(yPAM) !EM GCC4.7
1981                 zm1=zPAM                                   zm1=REAL(zPAM) !EM GCC4.7
 c     print*,'***',is1,xm1,ym1,zm1  
1982    
1983                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1984    c$$$                  print*,'(2) ip ',ip2
1985    c$$$     $                 ,mask_view(nviewx(ip2))
1986    c$$$     $                 ,mask_view(nviewy(ip2))
1987                    if(  mask_view(nviewx(ip2)).ne.0 .or.                    if(  mask_view(nviewx(ip2)).ne.0 .or.
1988       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1989                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                                        do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
1990                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1991                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1992                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1993    
1994    c$$$                        print*,'(2) ip ',ip2,' icp ',icp2
1995    
1996  c                        call xyz_PAM  c                        call xyz_PAM
1997  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1998  c                        call xyz_PAM  c                        call xyz_PAM
1999  c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)  c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2000                          call xyz_PAM                          call xyz_PAM
2001       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2002                          xm2=xPAM                          xm2=REAL(xPAM) !EM GCC4.7
2003                          ym2=yPAM                          ym2=REAL(yPAM) !EM GCC4.7
2004                          zm2=zPAM                          zm2=REAL(zPAM) !EM GCC4.7
2005                                                    
2006    *                       ---------------------------------------------------
2007    *                       both couples must have a y-cluster
2008    *                       (condition necessary when in RECOVER_SINGLETS mode)
2009    *                       ---------------------------------------------------
2010                            if(icy1.eq.0.or.icy2.eq.0)goto 111
2011    
2012                            if(cl_used(icy1).ne.0)goto 111
2013                            if(cl_used(icy2).ne.0)goto 111
2014    
2015                            
2016  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2017  *     track parameters on Y VIEW  *     track parameters on Y VIEW
2018  *     (2 couples needed)  *     (2 couples needed)
# Line 1977  c     $                       (icx2,icy2 Line 2022  c     $                       (icx2,icy2
2022       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2023       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2024       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2025  c                           good2=.false.  c     good2=.false.
2026  c                           goto 880 !fill ntp and go to next event  c     goto 880 !fill ntp and go to next event
2027                             do iv=1,12                             do iv=1,12
2028  c                              mask_view(iv) = 3  c     mask_view(iv) = 3
2029                                mask_view(iv) = mask_view(iv)+ 2**2                                mask_view(iv) = mask_view(iv)+ 2**2
2030                             enddo                             enddo
2031                             iflag=1                             iflag=1
2032                             return                             return
2033                          endif                          endif
2034                            
2035                            
2036    ccc                        print*,'<doublet> ',icp1,icp2
2037    
2038                          ndblt = ndblt + 1                          ndblt = ndblt + 1
2039  *     store doublet info  *     store doublet info
2040                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
# Line 1993  c                              mask_view Line 2042  c                              mask_view
2042  *     tg(th_yz)  *     tg(th_yz)
2043                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2044  *     y0 (cm)  *     y0 (cm)
2045                          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
2046                                                      
2047  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2048  ****  reject non phisical couples                    ****  ****  reject non phisical couples                    ****
2049  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2050                            if(SECOND_SEARCH)goto 111
2051                          if(                          if(
2052       $                       abs(alfayz2(ndblt)).gt.alfyz2_max       $                       abs(alfayz2(ndblt)).gt.alfyz2_max
2053       $                       .or.       $                       .or.
2054       $                       abs(alfayz1(ndblt)).gt.alfyz1_max       $                       abs(alfayz1(ndblt)).gt.alfyz1_max
2055       $                       )ndblt = ndblt-1       $                       )ndblt = ndblt-1
2056                                                    
2057  c$$$      if(iev.eq.33)then                          
2058  c$$$      print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2   111                    continue
 c$$$     $        ,' || ',icx1,icy1,icx2,icy2  
 c$$$     $        ,' || ',xm1,ym1,xm2,ym2  
 c$$$     $        ,' || ',alfayz2(ndblt),alfayz1(ndblt)  
 c$$$      endif  
 c$$$  
2059  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2060  *     track parameters on Y VIEW - end  *     track parameters on Y VIEW - end
2061  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2062    
2063    
2064                            if(icx1.ne.0)then
2065                               if(cl_used(icx1).ne.0)goto 31
2066                            endif
2067                            if(icx2.ne.0)then
2068                               if(cl_used(icx2).ne.0)goto 31
2069                            endif
2070    
2071                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2072    
2073                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2074    c$$$                           print*,'(3) ip ',ip3
2075    c$$$     $                          ,mask_view(nviewx(ip3))
2076    c$$$     $                          ,mask_view(nviewy(ip3))                          
2077                             if(  mask_view(nviewx(ip3)).ne.0 .or.                             if(  mask_view(nviewx(ip3)).ne.0 .or.
2078       $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane       $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2079                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
# Line 2026  c$$$ Line 2081  c$$$
2081                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
2082                                   icx3=clx(ip3,icp3)                                   icx3=clx(ip3,icp3)
2083                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2084    
2085    c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
2086    
2087    *     ---------------------------------------------------
2088    *     all three couples must have a x-cluster
2089    *     (condition necessary when in RECOVER_SINGLETS mode)
2090    *     ---------------------------------------------------
2091                                     if(
2092         $                                icx1.eq.0.or.
2093         $                                icx2.eq.0.or.
2094         $                                icx3.eq.0.or.
2095         $                                .false.)goto 29
2096                                    
2097                                     if(cl_used(icx1).ne.0)goto 29
2098                                     if(cl_used(icx2).ne.0)goto 29
2099                                     if(cl_used(icx3).ne.0)goto 29
2100    
2101  c                                 call xyz_PAM  c                                 call xyz_PAM
2102  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2103  c                                 call xyz_PAM  c                                 call xyz_PAM
# Line 2033  c     $                               (i Line 2105  c     $                               (i
2105                                   call xyz_PAM                                   call xyz_PAM
2106       $                                (icx3,icy3,is3,PFAdef,PFAdef       $                                (icx3,icy3,is3,PFAdef,PFAdef
2107       $                                ,0.,0.,0.,0.)       $                                ,0.,0.,0.,0.)
2108                                   xm3=xPAM                                   xm3=REAL(xPAM) !EM GCC4.7
2109                                   ym3=yPAM                                   ym3=REAL(yPAM) !EM GCC4.7
2110                                   zm3=zPAM                                   zm3=REAL(zPAM) !EM GCC4.7
2111    
2112    
2113  *     find the circle passing through the three points  *     find the circle passing through the three points
2114  c$$$                                 call tricircle(3,xp,zp,angp,resp,chi                                   iflag_t = DEBUG
 c$$$     $                                ,xc,zc,radius,iflag)  
                                  iflag = DEBUG  
2115                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2116       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag_t)
 c     print*,xc,zc,radius  
2117  *     the circle must intersect the reference plane  *     the circle must intersect the reference plane
2118                                   if(  cc                                 if(iflag.ne.0)goto 29
2119  c     $                                 (xc.le.-1.*xclimit.or.                                   if(iflag_t.ne.0)then
2120  c     $                                 xc.ge.xclimit).and.  *     if tricircle fails, evaluate a straight line
2121       $                                radius**2.ge.(ZINI-zc)**2.and.                                      if(DEBUG.eq.1)
2122       $                                iflag.eq.0.and.       $                                   print*,'TRICIRCLE failure'
2123       $                                .true.)then       $                                   ,' >>> straight line'
2124                                                                        radius=0.
2125                                        xc=0.
2126                                        yc=0.
2127                                        
2128                                        SZZ=0.                  
2129                                        SZX=0.
2130                                        SSX=0.
2131                                        SZ=0.
2132                                        S1=0.
2133                                        X0=0.
2134                                        Ax=0.
2135                                        BX=0.
2136                                        DO I=1,3
2137                                           XX = XP(I)
2138                                           SZZ=SZZ+ZP(I)*ZP(I)
2139                                           SZX=SZX+ZP(I)*XX
2140                                           SSX=SSX+XX
2141                                           SZ=SZ+ZP(I)
2142                                           S1=S1+1.
2143                                        ENDDO
2144                                        DET=SZZ*S1-SZ*SZ
2145                                        AX=(SZX*S1-SZ*SSX)/DET
2146                                        BX=(SZZ*SSX-SZX*SZ)/DET
2147                                        X0  = AX*REAL(ZINI)+BX ! EM GCC4.7
2148                                        
2149                                     endif
2150    
2151                                     if(  .not.SECOND_SEARCH.and.
2152         $                                radius**2.lt.(ZINI-zc)**2)goto 29
2153                                                                      
2154  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2155  *     track parameters on X VIEW  *     track parameters on X VIEW
2156  *     (3 couples needed)  *     (3 couples needed)
2157  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2158                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2159                                      if(verbose.eq.1)print*,                                      if(verbose.eq.1)print*,
2160       $                     '** warning ** number of identified '//       $                                   '** warning **'//
2161       $                     'triplets exceeds vector dimention '       $                                   ' number of identified '//
2162       $                    ,'( ',ntrpt_max,' )'       $                                   'triplets exceeds'//
2163  c                                    good2=.false.       $                                   ' vector dimention '
2164  c                                    goto 880 !fill ntp and go to next event       $                                   ,'( ',ntrpt_max,' )'
2165    c     good2=.false.
2166    c     goto 880 !fill ntp and go to next event
2167                                      do iv=1,nviews                                      do iv=1,nviews
2168  c                                       mask_view(iv) = 4  c     mask_view(iv) = 4
2169                                         mask_view(iv)=mask_view(iv)+ 2**3                                         mask_view(iv) =
2170         $                                      mask_view(iv)+ 2**3
2171                                      enddo                                      enddo
2172                                      iflag=1                                      iflag=1
2173                                      return                                      return
2174                                   endif                                   endif
2175    
2176    ccc                                 print*,'<triplet> ',icp1,icp2,icp3
2177                                    
2178                                   ntrpt = ntrpt +1                                   ntrpt = ntrpt +1
2179  *     store triplet info  *     store triplet info
2180                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
2181                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
2182                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
2183                                                                    
2184                                   if(xc.lt.0)then                                   if(radius.ne.0.and.xc.lt.0)then
2185  *************POSITIVE DEFLECTION  *************POSITIVE DEFLECTION
2186                alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc+sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2187                alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = (REAL(ZINI)-zc)/
2188                alfaxz3(ntrpt) = 1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2189                                   else             alfaxz3(ntrpt) = 1/radius
2190                                    else if(radius.ne.0.and.xc.ge.0)then
2191  *************NEGATIVE DEFLECTION  *************NEGATIVE DEFLECTION
2192                alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc-sqrt(radius**2-(REAL(ZINI)-zc)**2)
2193                alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/
2194                alfaxz3(ntrpt) = -1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2195                                   endif             alfaxz3(ntrpt) = -1/radius
2196                                                                    else if(radius.eq.0)then
2197    *************straight fit
2198                 alfaxz1(ntrpt) = X0
2199                 alfaxz2(ntrpt) = AX
2200                 alfaxz3(ntrpt) = 0.
2201                                    endif
2202    
2203    c$$$                                print*,'alfaxz1 ', alfaxz1(ntrpt)
2204    c$$$                                print*,'alfaxz2 ', alfaxz2(ntrpt)
2205    c$$$                                print*,'alfaxz3 ', alfaxz3(ntrpt)
2206                                        
2207  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2208  ****  reject non phisical triplets                   ****  ****  reject non phisical triplets                   ****
2209  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2210                                   if(                                  if(SECOND_SEARCH)goto 29
2211       $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max                                  if(
2212       $                                .or.       $                               abs(alfaxz2(ntrpt)).gt.
2213       $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max       $                               alfxz2_max
2214       $                                )ntrpt = ntrpt-1       $                               .or.
2215                                         $                               abs(alfaxz1(ntrpt)).gt.
2216                                         $                               alfxz1_max
2217  c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)       $                               )ntrpt = ntrpt-1
2218                                                                    
2219  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2220  *     track parameters on X VIEW - end  *     track parameters on X VIEW - end
2221  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2222                                endif                                
2223     29                           continue
2224                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2225                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2226   30                     continue   30                     continue
2227                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2228   31                  continue  
2229                         31                  continue                    
2230   1                enddo         !end loop on COPPIA 2  c 1                enddo         !end loop on COPPIA 2
2231                     enddo         !end loop on COPPIA 2
2232                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2233   20            continue   20            continue
2234              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2235                
2236    c 11         continue
2237              continue
2238           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2239        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2240   10   continue   10   continue
# Line 2184  c      include 'momanhough_init.f' Line 2306  c      include 'momanhough_init.f'
2306        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2307           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
2308                            
 c     print*,'--------------'  
 c     print*,'** ',idb1,' **'  
               
2309           do icp=1,ncp_tot           do icp=1,ncp_tot
2310              cp_useds1(icp)=0    !init              cp_useds1(icp)=0    !init
2311              cp_useds2(icp)=0    !init              cp_useds2(icp)=0    !init
# Line 2222  ccccc if(db_used(idbref).eq.1)goto 1188 Line 2341  ccccc if(db_used(idbref).eq.1)goto 1188
2341  *     doublet distance in parameter space  *     doublet distance in parameter space
2342                 distance=                 distance=
2343       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
2344       $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2                     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
2345                 distance = sqrt(distance)                 distance = sqrt(distance)
2346                                
 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  
2347                 if(distance.lt.cutdistyz)then                 if(distance.lt.cutdistyz)then
2348    
 c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  
2349                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
2350                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
2351                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
# Line 2250  c     print*,idb1,idb2,distance,' cloud Line 2361  c     print*,idb1,idb2,distance,' cloud
2361    
2362                    temp1 = temp1 + alfayz1(idb2)                    temp1 = temp1 + alfayz1(idb2)
2363                    temp2 = temp2 + alfayz2(idb2)                    temp2 = temp2 + alfayz2(idb2)
 c     print*,'*   idbref,idb2 ',idbref,idb2  
2364                 endif                               endif              
2365                                
2366   1118          continue   1118          continue
2367              enddo               !end loop (2) on DOUBLETS              enddo               !end loop (2) on DOUBLETS
2368                            
2369   1188       continue  c 1188       continue
2370                continue
2371           enddo                  !end loop on... bo?           enddo                  !end loop on... bo?
2372                    
2373           nptloop=npv           nptloop=npv
# Line 2273  c     print*,'*   idbref,idb2 ',idbref,i Line 2384  c     print*,'*   idbref,idb2 ',idbref,i
2384           enddo           enddo
2385           ncpused=0           ncpused=0
2386           do icp=1,ncp_tot           do icp=1,ncp_tot
2387              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2388         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2389         $           .true.)then
2390                 ncpused=ncpused+1                 ncpused=ncpused+1
2391                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2392                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2283  c     print*,'*   idbref,idb2 ',idbref,i Line 2396  c     print*,'*   idbref,idb2 ',idbref,i
2396           do ip=1,nplanes           do ip=1,nplanes
2397              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2398           enddo           enddo
2399  c     print*,'>>>> ',ncpused,npt,nplused          
 c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  
          if(npt.lt.nptyz_min)goto 2228 !next doublet  
2400           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2401                    
2402  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2317  c               mask_view(iv) = 5 Line 2428  c               mask_view(iv) = 5
2428  c     ptcloud_yz_nt(nclouds_yz)=npt  c     ptcloud_yz_nt(nclouds_yz)=npt
2429           do ipt=1,npt           do ipt=1,npt
2430              db_cloud(npt_tot+ipt) = db_all(ipt)              db_cloud(npt_tot+ipt) = db_all(ipt)
 c     print*,'>> ',ipt,db_all(ipt)  
2431           enddo             enddo  
2432           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2433           if(DEBUG.EQ.1)then           if(DEBUG.EQ.1)then
             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'  
2434              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2435              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2436              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
# Line 2330  c     print*,'>> ',ipt,db_all(ipt) Line 2439  c     print*,'>> ',ipt,db_all(ipt)
2439              print*,'cpcloud_yz '              print*,'cpcloud_yz '
2440       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2441              print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)              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)  
2442           endif           endif
2443  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2444  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2348  c$$$     $           ,(db_cloud(iii),iii Line 2453  c$$$     $           ,(db_cloud(iii),iii
2453        endif                            endif                    
2454                
2455        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
          print*,'---------------------- '  
2456           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
          print*,' '  
2457        endif        endif
2458                
2459                
# Line 2413  c      include 'momanhough_init.f' Line 2516  c      include 'momanhough_init.f'
2516   91   continue                     91   continue                  
2517        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2518           if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud           if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
 c     print*,'--------------'  
 c     print*,'** ',itr1,' **'  
2519                    
2520           do icp=1,ncp_tot           do icp=1,ncp_tot
2521              cp_useds1(icp)=0              cp_useds1(icp)=0
# Line 2448  c         tr_temp(1)=itr1 Line 2549  c         tr_temp(1)=itr1
2549              do itr2=1,ntrpt     !loop (2) on TRIPLETS              do itr2=1,ntrpt     !loop (2) on TRIPLETS
2550                 if(itr2.eq.itr1)goto 11188       !next triplet                 if(itr2.eq.itr1)goto 11188       !next triplet
2551                 if(tr_used(itr2).eq.1)goto 11188 !next triplet                               if(tr_used(itr2).eq.1)goto 11188 !next triplet              
2552    
2553    
2554  *     triplet distance in parameter space  *     triplet distance in parameter space
2555  *     solo i due parametri spaziali per il momemnto  *     solo i due parametri spaziali per il momemnto
2556                 distance=                 distance=
2557       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2558       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
2559                 distance = sqrt(distance)                 distance = sqrt(distance)
2560                  
2561    
2562  *     ------------------------------------------------------------------------  *     ------------------------------------------------------------------------
2563  *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE  *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2564  *     ------------------------------------------------------------------------  *     ------------------------------------------------------------------------
# Line 2467  c         tr_temp(1)=itr1 Line 2571  c         tr_temp(1)=itr1
2571       $              .true.)istrimage=1       $              .true.)istrimage=1
2572    
2573                 if(distance.lt.cutdistxz.or.istrimage.eq.1)then                 if(distance.lt.cutdistxz.or.istrimage.eq.1)then
 c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  
2574                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2575                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
2576                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
# Line 2486  c     print*,idb1,idb2,distance,' cloud Line 2589  c     print*,idb1,idb2,distance,' cloud
2589                    temp1 = temp1 + alfaxz1(itr2)                    temp1 = temp1 + alfaxz1(itr2)
2590                    temp2 = temp2 + alfaxz2(itr2)                    temp2 = temp2 + alfaxz2(itr2)
2591                    temp3 = temp3 + alfaxz3(itr2)                    temp3 = temp3 + alfaxz3(itr2)
 c     print*,'*   itrref,itr2 ',itrref,itr2,distance  
2592                 endif                               endif              
2593                                
2594  11188          continue  11188          continue
2595              enddo               !end loop (2) on TRIPLETS              enddo               !end loop (2) on TRIPLETS
2596                                                
2597  11888       continue  c11888       continue
2598                continue
2599           enddo                  !end loop on... bo?               enddo                  !end loop on... bo?    
2600                    
2601           nptloop=npv           nptloop=npv
# Line 2507  c     print*,'*   itrref,itr2 ',itrref,i Line 2610  c     print*,'*   itrref,itr2 ',itrref,i
2610  *     1bis)  *     1bis)
2611  *     2) it is not already stored  *     2) it is not already stored
2612  *     ------------------------------------------  *     ------------------------------------------
 c     print*,'check cp_used'  
2613           do ip=1,nplanes           do ip=1,nplanes
2614              hit_plane(ip)=0              hit_plane(ip)=0
2615           enddo           enddo
2616           ncpused=0           ncpused=0
2617           do icp=1,ncp_tot           do icp=1,ncp_tot
2618              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2619         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2620         $           .true.)then
2621                 ncpused=ncpused+1                 ncpused=ncpused+1
2622                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2623                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2523  c     print*,'check cp_used' Line 2627  c     print*,'check cp_used'
2627           do ip=1,nplanes           do ip=1,nplanes
2628              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2629           enddo           enddo
 c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  
          if(npt.lt.nptxz_min)goto 22288     !next triplet  
2630           if(nplused.lt.nplxz_min)goto 22288 !next triplet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2631                    
2632  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2558  c               mask_view(iv) = 6 Line 2660  c               mask_view(iv) = 6
2660           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2661                    
2662           if(DEBUG.EQ.1)then           if(DEBUG.EQ.1)then
2663              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
             print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                
2664              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2665              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2666              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
# Line 2568  c               mask_view(iv) = 6 Line 2669  c               mask_view(iv) = 6
2669              print*,'cpcloud_xz '              print*,'cpcloud_xz '
2670       $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)       $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2671              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              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)  
2672           endif           endif
2673  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2674  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2585  c$$$     $           ,(tr_cloud(iii),iii Line 2682  c$$$     $           ,(tr_cloud(iii),iii
2682         endif                             endif                    
2683                
2684        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
          print*,'---------------------- '  
2685           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
          print*,' '  
2686        endif        endif
2687                
2688                
# Line 2643  c$$$     $           ,(tr_cloud(iii),iii Line 2738  c$$$     $           ,(tr_cloud(iii),iii
2738    
2739        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
2740                
2741        do iyz=1,nclouds_yz       !loop on YZ couds        do iyz=1,nclouds_yz       !loop on YZ clouds
2742           do ixz=1,nclouds_xz    !loop on XZ couds           do ixz=1,nclouds_xz    !loop on XZ clouds
2743                            
2744  *     --------------------------------------------------  *     --------------------------------------------------
2745  *     check of consistency of the clouds  *     check of consistency of the clouds
# Line 2653  c$$$     $           ,(tr_cloud(iii),iii Line 2748  c$$$     $           ,(tr_cloud(iii),iii
2748  *     of the two clouds  *     of the two clouds
2749  *     --------------------------------------------------  *     --------------------------------------------------
2750              do ip=1,nplanes              do ip=1,nplanes
2751                 hit_plane(ip)=0                 hit_plane(ip)=0  !n.matching couples (REAL couples, not SINGLETS)
2752                 ncp_match(ip)=0                 ncp_match(ip)=0  !n.matching couples per plane
2753                 do icpp=1,ncouplemax                 do icpp=1,ncouplemax
2754                    cp_match(ip,icpp)=0 !init couple list                    cp_match(ip,icpp)=0 !init couple list
2755                 enddo                 enddo
2756              enddo              enddo
2757              ncp_ok=0              ncp_ok=0            !count n.matching-couples
2758                ncpx_ok=0           !count n.matching-couples with x cluster
2759                ncpy_ok=0           !count n.matching-couples with y cluster
2760    
2761    
2762              do icp=1,ncp_tot    !loop over couples              do icp=1,ncp_tot    !loop over couples
2763  *     get info on  
2764                 cpintersec(icp)=min(                 if(.not.RECOVER_SINGLETS)then
2765       $              cpcloud_yz(iyz,icp),  *     ------------------------------------------------------
2766       $              cpcloud_xz(ixz,icp))  *     if NOT in RECOVER_SINGLETS mode, take the intersection
2767                 if(  *     between xz yz clouds
2768       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.  *     ------------------------------------------------------
2769       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.                    cpintersec(icp)=min(
2770       $              .false.)cpintersec(icp)=0       $                 cpcloud_yz(iyz,icp),
2771         $                 cpcloud_xz(ixz,icp))
2772  *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp  *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2773    *     ------------------------------------------------------
2774    *     discard the couple if the sensor is in conflict
2775    *     ------------------------------------------------------
2776                      if(
2777         $       (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2778         $       (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2779         $                 .false.)cpintersec(icp)=0
2780                   else
2781    *     ------------------------------------------------------
2782    *     if RECOVER_SINGLETS take the union
2783    *     (otherwise the fake couples formed by singlets would be
2784    *     discarded)    
2785    *     ------------------------------------------------------
2786                      cpintersec(icp)=max(
2787         $                 cpcloud_yz(iyz,icp),
2788         $                 cpcloud_xz(ixz,icp))
2789    c$$$                  if(cpcloud_yz(iyz,icp).gt.0)
2790    c$$$     $                 cpintersec(icp)=cpcloud_yz(iyz,icp)
2791    *     cpintersec is >0 if either yz or xz cloud contains the couple icp
2792                   endif
2793    
2794    c$$$               print*,icp,ip_cp(icp),' -- ',cpintersec(icp)
2795    
2796                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
                   ncp_ok=ncp_ok+1    
2797                                        
2798                    ip=ip_cp(icp)                    ip=ip_cp(icp)
2799                    hit_plane(ip)=1                    hit_plane(ip)=1
2800    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0)
2801    c$$$     $                 ncp_ok=ncp_ok+1  
2802    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0)
2803    c$$$     $                 ncpx_ok=ncpx_ok+1
2804    c$$$                  if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0)
2805    c$$$     $                 ncpy_ok=ncpy_ok+1
2806    
2807                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2808         $                 cpcloud_xz(ixz,icp).gt.0)
2809         $                 ncp_ok=ncp_ok+1  
2810                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2811         $                 cpcloud_xz(ixz,icp).eq.0)
2812         $                 ncpy_ok=ncpy_ok+1  
2813                      if(  cpcloud_yz(iyz,icp).eq.0.and.
2814         $                 cpcloud_xz(ixz,icp).gt.0)
2815         $                 ncpx_ok=ncpx_ok+1  
2816    
2817                    if(cpintersec(icp).eq.1)then                    if(cpintersec(icp).eq.1)then
2818  *     1) only the couple image in sensor 1 matches  *     1) only the couple image in sensor 1 matches
2819                       id=-icp                       id=-icp
# Line 2701  c$$$     $           ,(tr_cloud(iii),iii Line 2840  c$$$     $           ,(tr_cloud(iii),iii
2840              do ip=1,nplanes              do ip=1,nplanes
2841                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2842              enddo              enddo
2843                                        
2844                            if(nplused.lt.3)goto 888 !next combination
2845    ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
2846    ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
2847    *     -----------------------------------------------------------
2848    *     if in RECOVER_SINGLET mode, the two clouds must have
2849    *     at least ONE intersecting real couple
2850    *     -----------------------------------------------------------
2851                if(ncp_ok.lt.1)goto 888 !next combination
2852    
2853              if(DEBUG.EQ.1)then              if(DEBUG.EQ.1)then
2854                 print*,'Combination ',iyz,ixz                 print*,'////////////////////////////'
2855       $              ,' db ',ptcloud_yz(iyz)                 print*,'Cloud combination (Y,X): ',iyz,ixz
2856       $              ,' tr ',ptcloud_xz(ixz)                 print*,' db ',ptcloud_yz(iyz)
2857       $              ,'  -----> # matching couples ',ncp_ok                 print*,' tr ',ptcloud_xz(ixz)
2858              endif                 print*,'  -----> # matching couples ',ncp_ok
2859                   print*,'  -----> # fake couples (X)',ncpx_ok
2860  c            if(nplused.lt.nplxz_min)goto 888 !next combination                 print*,'  -----> # fake couples (Y)',ncpy_ok
2861              if(nplused.lt.nplyz_min)goto 888 !next combination                 do icp=1,ncp_tot
2862              if(ncp_ok.lt.ncpok)goto 888 !next combination                    print*,'cp ',icp,' >'
2863         $                 ,' x ',cpcloud_xz(ixz,icp)
2864  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'       $                 ,' y ',cpcloud_yz(iyz,icp)
2865  c$$$  print*,'Configurazione cluster XZ'       $                 ,' ==> ',cpintersec(icp)
2866  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))                 enddo
2867  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  
2868                                                    
2869              if(DEBUG.EQ.1)then              if(DEBUG.EQ.1)then
                print*,'track candidate', ntracks+1  
2870                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2871                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2872                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2793  c$$$            if(AL_INI(5).gt.defmax)g Line 2899  c$$$            if(AL_INI(5).gt.defmax)g
2899                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2900                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2901                                                                
2902                                  if(DEBUG.eq.1)
2903         $                             print*,'combination: '
2904         $                             ,cp_match(6,icp1)
2905         $                             ,cp_match(5,icp2)
2906         $                             ,cp_match(4,icp3)
2907         $                             ,cp_match(3,icp4)
2908         $                             ,cp_match(2,icp5)
2909         $                             ,cp_match(1,icp6)
2910    
2911    
2912  *                             ---------------------------------------  *                             ---------------------------------------
2913  *                             check if this group of couples has been  *                             check if this group of couples has been
2914  *                             already fitted      *                             already fitted    
# Line 2841  c     $                                 Line 2957  c     $                                
2957       $                                   PFAdef,PFAdef,0.,0.,0.,0.)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2958  *                                   *************************  *                                   *************************
2959  *                                   -----------------------------  *                                   -----------------------------
2960                                      xgood(nplanes-ip+1)=1.                                      if(icx.gt.0.and.icy.gt.0)then
2961                                      ygood(nplanes-ip+1)=1.                                         xgood(nplanes-ip+1)=1.
2962                                      xm(nplanes-ip+1)=xPAM                                         ygood(nplanes-ip+1)=1.
2963                                      ym(nplanes-ip+1)=yPAM                                         xm(nplanes-ip+1)=xPAM
2964                                      zm(nplanes-ip+1)=zPAM                                         ym(nplanes-ip+1)=yPAM
2965                                      resx(nplanes-ip+1)=resxPAM                                         zm(nplanes-ip+1)=zPAM
2966                                      resy(nplanes-ip+1)=resyPAM                                         resx(nplanes-ip+1)=resxPAM
2967                                           resy(nplanes-ip+1)=resyPAM
2968                                           if(DEBUG.EQ.1)print*,'(X,Y)'
2969         $                                      ,nplanes-ip+1,xPAM,yPAM
2970                                        else
2971                                           xm_A(nplanes-ip+1) = xPAM_A
2972                                           ym_A(nplanes-ip+1) = yPAM_A
2973                                           xm_B(nplanes-ip+1) = xPAM_B
2974                                           ym_B(nplanes-ip+1) = yPAM_B
2975                                           zm(nplanes-ip+1)
2976         $                                      = (zPAM_A+zPAM_B)/2.
2977                                           resx(nplanes-ip+1) = resxPAM
2978                                           resy(nplanes-ip+1) = resyPAM
2979                                           if(icx.eq.0.and.icy.gt.0)then
2980                                              xgood(nplanes-ip+1)=0.
2981                                              ygood(nplanes-ip+1)=1.
2982                                              resx(nplanes-ip+1) = 1000.
2983                                              if(DEBUG.EQ.1)print*,'(  Y)'
2984         $                                         ,nplanes-ip+1,xPAM,yPAM
2985                                           elseif(icx.gt.0.and.icy.eq.0)then
2986                                              xgood(nplanes-ip+1)=1.
2987                                              ygood(nplanes-ip+1)=0.
2988                                              if(DEBUG.EQ.1)print*,'(X  )'
2989         $                                         ,nplanes-ip+1,xPAM,yPAM
2990                                              resy(nplanes-ip+1) = 1000.
2991                                           else
2992                                              print*,'both icx=0 and icy=0'
2993         $                                         ,' ==> IMPOSSIBLE!!'
2994                                           endif
2995                                        endif
2996  *                                   -----------------------------  *                                   -----------------------------
2997                                   endif                                   endif
2998                                enddo !end loop on planes                                enddo !end loop on planes
# Line 2888  c                                 chi2=- Line 3033  c                                 chi2=-
3033  *     **********************************************************  *     **********************************************************
3034    
3035                                if(chi2.le.0.)goto 666                                              if(chi2.le.0.)goto 666              
3036                                  if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
3037                                  if(chi2.ne.chi2)goto 666  !OPTIMIZATION
3038    
3039  *     --------------------------  *     --------------------------
3040  *     STORE candidate TRACK INFO  *     STORE candidate TRACK INFO
# Line 2914  c                                    mas Line 3061  c                                    mas
3061    
3062                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
3063                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
3064                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))
3065                                   XM_STORE(ip,ntracks)=sngl(xm(ip))                                   XM_STORE(ip,ntracks)=sngl(xm(ip))
3066                                   YM_STORE(ip,ntracks)=sngl(ym(ip))                                   YM_STORE(ip,ntracks)=sngl(ym(ip))
3067                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))
# Line 2933  c                                    mas Line 3080  c                                    mas
3080       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
3081                                      SENSOR_STORE(nplanes-ip+1,ntracks)                                      SENSOR_STORE(nplanes-ip+1,ntracks)
3082       $                              = is_cp(cp_match(ip,hit_plane(ip)))       $                              = is_cp(cp_match(ip,hit_plane(ip)))
3083                                      LADDER_STORE(nplanes-ip+1,ntracks)                                      
3084       $                                   = LADDER(                                      icl=
3085       $                                   clx(ip,icp_cp(       $                                   clx(ip,icp_cp(
3086       $                                   cp_match(ip,hit_plane(ip)       $                                   cp_match(ip,hit_plane(ip)
3087       $                                   ))));       $                                   )));
3088                                        if(icl.eq.0)
3089         $                                   icl=
3090         $                                   cly(ip,icp_cp(
3091         $                                   cp_match(ip,hit_plane(ip)
3092         $                                   )));
3093    
3094                                        LADDER_STORE(nplanes-ip+1,ntracks)
3095         $                                   = LADDER(icl);
3096                                   else                                   else
3097                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
3098                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0
# Line 2951  c                                    mas Line 3106  c                                    mas
3106                                   enddo                                   enddo
3107                                enddo                                enddo
3108                                                                
3109                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=REAL(chi2)
3110                                                                
3111  *     --------------------------------  *     --------------------------------
3112  *     STORE candidate TRACK INFO - end  *     STORE candidate TRACK INFO - end
# Line 2971  c                                    mas Line 3126  c                                    mas
3126                
3127        if(ntracks.eq.0)then        if(ntracks.eq.0)then
3128           iflag=1           iflag=1
3129           return  cc         return
3130        endif        endif
3131                
 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  
3132        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
3133          print*,'****** TRACK CANDIDATES *****************'          print*,'****** TRACK CANDIDATES *****************'
3134          print*,'#         R. chi2        RIG         ndof'          print*,'#         R. chi2        RIG         ndof'
# Line 3028  c$$$      endif Line 3174  c$$$      endif
3174        character*10 PFA        character*10 PFA
3175        common/FINALPFA/PFA        common/FINALPFA/PFA
3176    
3177          double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7
3178          double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7
3179          double precision zmm,zmm_A,zmm_B !EM GCC4.7
3180          double precision clincnewc !EM GCC4.7
3181          double precision clincnew !EM GCC4.7
3182    
3183        real k(6)        real k(6)
3184        DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/        DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3185    
# Line 3044  c$$$      endif Line 3196  c$$$      endif
3196        call track_init        call track_init
3197        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3198    
3199             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3200    
3201           xP=XV_STORE(nplanes-ip+1,ibest)           xP=XV_STORE(nplanes-ip+1,ibest)
3202           yP=YV_STORE(nplanes-ip+1,ibest)           yP=YV_STORE(nplanes-ip+1,ibest)
3203           zP=ZV_STORE(nplanes-ip+1,ibest)           zP=ZV_STORE(nplanes-ip+1,ibest)
# Line 3060  c$$$         bxyz(3)=0 Line 3214  c$$$         bxyz(3)=0
3214  *     using improved PFAs  *     using improved PFAs
3215  *     -------------------------------------------------  *     -------------------------------------------------
3216  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3217           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.  c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3218    c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3219    c$$$            
3220    c$$$            id=CP_STORE(nplanes-ip+1,ibest)
3221    c$$$            
3222    c$$$            is=is_cp(id)
3223    c$$$            icp=icp_cp(id)
3224    c$$$            if(ip_cp(id).ne.ip)
3225    c$$$     $           print*,'OKKIO!!'
3226    c$$$     $           ,'id ',id,is,icp
3227    c$$$     $           ,ip_cp(id),ip
3228    c$$$            icx=clx(ip,icp)
3229    c$$$            icy=cly(ip,icp)
3230    c$$$c            call xyz_PAM(icx,icy,is,
3231    c$$$c     $           PFA,PFA,
3232    c$$$c     $           AXV_STORE(nplanes-ip+1,ibest),
3233    c$$$c     $           AYV_STORE(nplanes-ip+1,ibest))
3234    c$$$            call xyz_PAM(icx,icy,is,
3235    c$$$     $           PFA,PFA,
3236    c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
3237    c$$$     $           AYV_STORE(nplanes-ip+1,ibest),
3238    c$$$     $           bxyz(1),
3239    c$$$     $           bxyz(2)
3240    c$$$     $           )
3241    c$$$
3242    c$$$            xm(nplanes-ip+1) = xPAM
3243    c$$$            ym(nplanes-ip+1) = yPAM
3244    c$$$            zm(nplanes-ip+1) = zPAM
3245    c$$$            xgood(nplanes-ip+1) = 1
3246    c$$$            ygood(nplanes-ip+1) = 1
3247    c$$$            resx(nplanes-ip+1) = resxPAM
3248    c$$$            resy(nplanes-ip+1) = resyPAM
3249    c$$$
3250    c$$$            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3251    c$$$            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3252             if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
3253       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3254                            
3255              id=CP_STORE(nplanes-ip+1,ibest)              id=CP_STORE(nplanes-ip+1,ibest)
# Line 3085  c     $           AYV_STORE(nplanes-ip+1 Line 3274  c     $           AYV_STORE(nplanes-ip+1
3274       $           bxyz(2)       $           bxyz(2)
3275       $           )       $           )
3276    
3277              xm(nplanes-ip+1) = xPAM              if(icx.gt.0.and.icy.gt.0)then
3278              ym(nplanes-ip+1) = yPAM                 xm(nplanes-ip+1) = xPAM
3279              zm(nplanes-ip+1) = zPAM                 ym(nplanes-ip+1) = yPAM
3280              xgood(nplanes-ip+1) = 1                 zm(nplanes-ip+1) = zPAM
3281              ygood(nplanes-ip+1) = 1                 xm_A(nplanes-ip+1) = 0.
3282              resx(nplanes-ip+1) = resxPAM                 ym_A(nplanes-ip+1) = 0.
3283              resy(nplanes-ip+1) = resyPAM                 xm_B(nplanes-ip+1) = 0.
3284                   ym_B(nplanes-ip+1) = 0.
3285              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))                 xgood(nplanes-ip+1) = 1
3286              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))                 ygood(nplanes-ip+1) = 1
3287                   resx(nplanes-ip+1) = resxPAM
3288                   resy(nplanes-ip+1) = resyPAM
3289                   dedxtrk_x(nplanes-ip+1)=
3290         $              sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3291                   dedxtrk_y(nplanes-ip+1)=
3292         $              sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3293                else
3294                   xm(nplanes-ip+1) = 0.
3295                   ym(nplanes-ip+1) = 0.
3296                   zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
3297                   xm_A(nplanes-ip+1) = xPAM_A
3298                   ym_A(nplanes-ip+1) = yPAM_A
3299                   xm_B(nplanes-ip+1) = xPAM_B
3300                   ym_B(nplanes-ip+1) = yPAM_B
3301                   xgood(nplanes-ip+1) = 0
3302                   ygood(nplanes-ip+1) = 0
3303                   resx(nplanes-ip+1) = 1000.!resxPAM
3304                   resy(nplanes-ip+1) = 1000.!resyPAM
3305                   dedxtrk_x(nplanes-ip+1)= 0
3306                   dedxtrk_y(nplanes-ip+1)= 0
3307                   if(icx.gt.0)then
3308                      xgood(nplanes-ip+1) = 1
3309                      resx(nplanes-ip+1) = resxPAM
3310                      dedxtrk_x(nplanes-ip+1)=
3311         $                 sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3312                   elseif(icy.gt.0)then
3313                      ygood(nplanes-ip+1) = 1
3314                      resy(nplanes-ip+1) = resyPAM
3315                      dedxtrk_y(nplanes-ip+1)=
3316         $                 sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3317                   endif
3318                endif
3319                            
3320  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3321  *     -------------------------------------------------  *     -------------------------------------------------
# Line 3106  c     $           AYV_STORE(nplanes-ip+1 Line 3327  c     $           AYV_STORE(nplanes-ip+1
3327                                
3328              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
3329              ygood(nplanes-ip+1)=0              ygood(nplanes-ip+1)=0
3330    
3331                CP_STORE(nplanes-ip+1,ibest)=0 !re-init
3332                CLS_STORE(nplanes-ip+1,ibest)=0
3333    
3334                                
3335  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3336  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
# Line 3128  c     $           AYV_STORE(nplanes-ip+1 Line 3353  c     $           AYV_STORE(nplanes-ip+1
3353  *     ===========================================  *     ===========================================
3354  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3355  *     ===========================================  *     ===========================================
3356  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'              distmin=100000000.
             distmin=1000000.  
3357              xmm = 0.              xmm = 0.
3358              ymm = 0.              ymm = 0.
3359              zmm = 0.              zmm = 0.
# Line 3142  c            if(DEBUG.EQ.1)print*,'>>>> Line 3366  c            if(DEBUG.EQ.1)print*,'>>>>
3366              do icp=1,ncp_plane(ip) !loop on couples on plane icp              do icp=1,ncp_plane(ip) !loop on couples on plane icp
3367                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3368                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3369                   if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
3370                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3371  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
3372  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
3373       $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)       $              cl_used(icx).ne.0.or. !or the X cluster is already used
3374       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)       $              cl_used(icy).ne.0.or. !or the Y cluster is already used
3375       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3376  *            *          
3377                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3160  c     $              cl_used(icy).eq.1.o Line 3385  c     $              cl_used(icy).eq.1.o
3385                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3386  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3387                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3388                 if(DEBUG.EQ.1)print*,'( couple ',id                 if(DEBUG.EQ.1)
3389         $              print*,'( couple ',id
3390       $              ,' ) distance ',distance       $              ,' ) distance ',distance
3391                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3392                    xmm = xPAM                    xmm = xPAM
# Line 3172  c               distance = distance / RC Line 3398  c               distance = distance / RC
3398                    idm = id                                      idm = id                  
3399                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3400                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3401  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!                    clincnewc=10.*dsqrt(rymm**2+rxmm**2
3402                    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  
3403                 endif                 endif
3404   1188          continue   1188          continue
3405              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3406  c            if(distmin.le.clinc)then     !QUIQUI                            if(distmin.le.clincnewc)then    
             if(distmin.le.clincnewc)then     !QUIQUI                
3407  *              -----------------------------------  *              -----------------------------------
3408                 xm(nplanes-ip+1) = xmm !<<<                 xm(nplanes-ip+1) = xmm !<<<
3409                 ym(nplanes-ip+1) = ymm !<<<                 ym(nplanes-ip+1) = ymm !<<<
# Line 3193  c            if(distmin.le.clinc)then   Line 3417  c            if(distmin.le.clinc)then  
3417  *              -----------------------------------  *              -----------------------------------
3418                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3419                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3420       $              ,' (dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clincnewc,' )'
3421                 goto 133         !next plane                 goto 133         !next plane
3422              endif              endif
3423  *     ================================================  *     ================================================
3424  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3425  *                     either from a couple or single  *                     either from a couple or single
3426  *     ================================================  *     ================================================
 c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'  
3427              distmin=1000000.              distmin=1000000.
3428              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3429              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3219  c            if(DEBUG.EQ.1)print*,'>>>> Line 3442  c            if(DEBUG.EQ.1)print*,'>>>>
3442              do icp=1,ncp_plane(ip) !loop on cluster inside couples              do icp=1,ncp_plane(ip) !loop on cluster inside couples
3443                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3444                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3445                   if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
3446                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3447                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
3448  *                                                !jump to the next couple  *                                                !jump to the next couple
# Line 3237  c     $              AXV_STORE(nplanes-i Line 3461  c     $              AXV_STORE(nplanes-i
3461       $              )                     $              )              
3462                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3463  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3464                 if(DEBUG.EQ.1)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)
3465         $              print*,'( cl-X ',icx
3466       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3467                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3468                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3270  c     $              0.,AYV_STORE(nplane Line 3495  c     $              0.,AYV_STORE(nplane
3495       $              )       $              )
3496                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3497  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3498                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)
3499         $              print*,'( cl-Y ',icy
3500       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3501                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3502                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3290  c                 dedxmm = sgnl(icy)  !( Line 3516  c                 dedxmm = sgnl(icy)  !(
3516  11882          continue  11882          continue
3517              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3518  *----- single clusters -----------------------------------------------    *----- single clusters -----------------------------------------------  
 c            print*,'## ncls(',ip,') ',ncls(ip)  
3519              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3520                 icl=cls(ip,ic)                 icl=cls(ip,ic)
 c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  
3521                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
3522       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3523       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
# Line 3315  c               if(cl_used(icl).eq.1.or. Line 3539  c               if(cl_used(icl).eq.1.or.
3539    
3540                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3541  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3542                 if(DEBUG.EQ.1)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)
3543         $              print*,'( cl-s ',icl
3544       $              ,' ) distance ',distance       $              ,' ) distance ',distance
3545                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
 c                  if(DEBUG.EQ.1)print*,'YES'  
3546                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3547                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3548                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3339  c                  if(DEBUG.EQ.1)print*, Line 3563  c                  if(DEBUG.EQ.1)print*,
3563                 endif                                   endif                  
3564  18882          continue  18882          continue
3565              enddo               !end loop on single clusters              enddo               !end loop on single clusters
 c            print*,'## distmin ', distmin,' clinc ',clinc  
3566    
 c     QUIQUI------------  
 c     anche qui: non ci vuole clinc???  
3567              if(iclm.ne.0)then              if(iclm.ne.0)then
3568                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3569                    clincnew=                    clincnew=
3570       $                 20*       $            20.*     !EM GCC4.7
3571       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))       $            dsqrt(rxmm**2 +
3572         $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1))
3573                 else if(mod(VIEW(iclm),2).ne.0)then                 else if(mod(VIEW(iclm),2).ne.0)then
3574                    clincnew=                    clincnew=
3575       $                 10*       $            10.* !EM GCC4.7
3576       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))       $            dsqrt(rymm**2 +
3577         $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2))
3578                 endif                 endif
3579  c     QUIQUI------------  
3580                                 if(distmin.le.clincnew)then  
                if(distmin.le.clincnew)then   !QUIQUI  
 c     if(distmin.le.clinc)then          !QUIQUI            
3581                                        
3582                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                        CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3583  *     ----------------------------  *     ----------------------------
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3584                    if(mod(VIEW(iclm),2).eq.0)then                    if(mod(VIEW(iclm),2).eq.0)then
3585                       XGOOD(nplanes-ip+1)=1.                       XGOOD(nplanes-ip+1)=1.
3586                       resx(nplanes-ip+1)=rxmm                       resx(nplanes-ip+1)=rxmm
3587                       if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm                       if(DEBUG.EQ.1)
3588         $                    print*,'%%%% included X-cl ',iclm
3589       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3590       $                    ,', dist.= ',distmin       $                    ,', dist.= ',distmin
3591       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clincnew,' )'
3592                    else                    else
3593                       YGOOD(nplanes-ip+1)=1.                       YGOOD(nplanes-ip+1)=1.
3594                       resy(nplanes-ip+1)=rymm                       resy(nplanes-ip+1)=rymm
3595                       if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm                       if(DEBUG.EQ.1)
3596         $                    print*,'%%%% included Y-cl ',iclm
3597       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3598       $                    ,', dist.= ', distmin       $                    ,', dist.= ', distmin
3599       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clincnew,' )'
3600                    endif                    endif
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3601  *     ----------------------------  *     ----------------------------
3602                    xm_A(nplanes-ip+1) = xmm_A                    xm_A(nplanes-ip+1) = xmm_A
3603                    ym_A(nplanes-ip+1) = ymm_A                    ym_A(nplanes-ip+1) = ymm_A
# Line 3397  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3618  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3618        return        return
3619        end        end
3620    
3621    
3622  ***************************************************  ***************************************************
3623  *                                                 *  *                                                 *
3624  *                                                 *  *                                                 *
# Line 3406  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3628  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3628  *                                                 *  *                                                 *
3629  **************************************************  **************************************************
3630  *  *
       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  
   
   
   
3631    
3632    
3633    
# Line 3551  c$$$               cl_used(icl)=ntrk   ! Line 3701  c$$$               cl_used(icl)=ntrk   !
3701          ys(1,ip)=0          ys(1,ip)=0
3702          ys(2,ip)=0          ys(2,ip)=0
3703          sgnlys(ip)=0          sgnlys(ip)=0
3704            sxbad(ip)=0
3705            sybad(ip)=0
3706            multmaxsx(ip)=0
3707            multmaxsy(ip)=0
3708        enddo        enddo
3709        end        end
3710    
# Line 3669  c$$$               cl_used(icl)=ntrk   ! Line 3823  c$$$               cl_used(icl)=ntrk   !
3823        character*10 PFA        character*10 PFA
3824        common/FINALPFA/PFA        common/FINALPFA/PFA
3825    
3826        real sinth,phi,pig        real sinth,phi,pig, npig ! EM GCC4.7
3827        integer ssensor,sladder        integer ssensor,sladder
3828        pig=acos(-1.)        pig=acos(-1.)
3829    
3830    
3831    
3832  *     -------------------------------------  *     -------------------------------------
3833        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3834        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3835  *     -------------------------------------  *     -------------------------------------
3836        phi   = al(4)                  phi   = REAL(al(4))
3837        sinth = al(3)                    sinth = REAL(al(3))
3838        if(sinth.lt.0)then              if(sinth.lt.0)then      
3839           sinth = -sinth                   sinth = -sinth        
3840           phi = phi + pig                 phi = phi + pig      
# Line 3720  c$$$               cl_used(icl)=ntrk   ! Line 3876  c$$$               cl_used(icl)=ntrk   !
3876           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3877           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3878        
3879    
3880    ccc         print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3881    
3882           ax   = axv_nt(ip,ntr)           ax   = axv_nt(ip,ntr)
3883           ay   = ayv_nt(ip,ntr)           ay   = ayv_nt(ip,ntr)
3884           bfx  = BX_STORE(ip,IDCAND)           bfx  = BX_STORE(ip,IDCAND)
# Line 3735  c$$$         angy    = 180.*atan(tgtemp) Line 3894  c$$$         angy    = 180.*atan(tgtemp)
3894           angy = effectiveangle(ay,2*ip-1,bfx)           angy = effectiveangle(ay,2*ip-1,bfx)
3895                    
3896                    
 c         print*,'* ',ip,bfx,bfy,angx,angy  
3897    
3898           id  = CP_STORE(ip,IDCAND) ! couple id           id  = CP_STORE(ip,IDCAND) ! couple id
3899           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
# Line 3746  c         print*,'* ',ip,bfx,bfy,angx,an Line 3904  c         print*,'* ',ip,bfx,bfy,angx,an
3904           if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align           if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3905           LS(IP,ntr)      = ssensor+10*sladder           LS(IP,ntr)      = ssensor+10*sladder
3906    
3907           if(id.ne.0)then  c         if(id.ne.0)then
3908    CCCCCC(10 novembre 2009) PATCH X NUCLEI
3909    C     non ho capito perche', ma durante il ritracciamento dei nuclei
3910    C     (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
3911    C     che non viene reinizializzata correttamente e i cluster esclusi
3912    C     dal fit risultano ancora inclusi...
3913    C
3914             cltrx(ip,ntr)   = 0
3915             cltry(ip,ntr)   = 0
3916    c$$$         if(
3917    c$$$     $        xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
3918    c$$$     $        .and.
3919    c$$$     $        id.ne.0)then
3920             if(id.ne.0)then        !patch 30/12/09
3921    
3922  c           >>> is a couple  c           >>> is a couple
3923              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3924              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3925    
3926              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            
3927    
3928              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)))  
3929    
3930                   xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3931    
3932              if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)                 if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3933       $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)       $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3934              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)                
3935       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)                 multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3936         $              +10000*mult(cltrx(ip,ntr))
3937                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3938         $              /clsigma(indmax(cltrx(ip,ntr)))
3939                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3940                   xpu(ip,ntr)      = corr
3941    
3942              multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))              endif
3943       $                         +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  
3944    
3945              multmaxy(ip,ntr) = maxs(cltry(ip,ntr))                 cl_used(cltry(ip,ntr)) = 1 !tag used clusters          
3946       $                         +10000*mult(cltry(ip,ntr))  
3947              seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))                 ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3948       $           /clsigma(indmax(cltry(ip,ntr)))  
3949              call applypfa(PFA,cltry(ip,ntr),angy,corr,res)                 if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3950              ypu(ip,ntr)      = corr       $              dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3951                  
3952                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3953         $              +10000*mult(cltry(ip,ntr))
3954                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3955         $              /clsigma(indmax(cltry(ip,ntr)))
3956                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3957                   ypu(ip,ntr)      = corr
3958                endif
3959    
3960           elseif(icl.ne.0)then  c$$$         elseif(icl.ne.0)then
3961             endif                  !patch 30/12/09
3962             if(icl.ne.0)then       !patch 30/12/09
3963    
3964              cl_used(icl) = 1    !tag used clusters                        cl_used(icl) = 1    !tag used clusters          
3965    
3966              if(mod(VIEW(icl),2).eq.0)then              if(mod(VIEW(icl),2).eq.0)then
3967                 cltrx(ip,ntr)=icl                 cltrx(ip,ntr)=icl
   
3968                 xbad(ip,ntr) = nbadstrips(4,icl)                 xbad(ip,ntr) = nbadstrips(4,icl)
3969    
3970                 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)                 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
# Line 3797  c           >>> is a couple Line 3978  c           >>> is a couple
3978    
3979              elseif(mod(VIEW(icl),2).eq.1)then              elseif(mod(VIEW(icl),2).eq.1)then
3980                 cltry(ip,ntr)=icl                 cltry(ip,ntr)=icl
   
3981                 ybad(ip,ntr) = nbadstrips(4,icl)                 ybad(ip,ntr) = nbadstrips(4,icl)
3982    
3983                 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)                 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
# Line 3821  c           >>> is a couple Line 4001  c           >>> is a couple
4001           do ip=1,6           do ip=1,6
4002              print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)              print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
4003           enddo           enddo
4004             print*,'dedx: '
4005             do ip=1,6
4006                print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
4007             enddo
4008        endif        endif
4009    
 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*,'-----------------------'  
   
4010        end        end
4011    
4012        subroutine fill_level2_siglets        subroutine fill_level2_siglets
# Line 3865  c         if( mask_view(iv).ne.0 )good2( Line 4042  c         if( mask_view(iv).ne.0 )good2(
4042           ip=nplanes-npl(VIEW(icl))+1                       ip=nplanes-npl(VIEW(icl))+1            
4043                    
4044           if(cl_used(icl).eq.0)then !cluster not included in any track           if(cl_used(icl).eq.0)then !cluster not included in any track
4045    
4046              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
4047    
4048                 nclsx = nclsx + 1                 nclsx = nclsx + 1
4049                 planex(nclsx) = ip                 planex(nclsx) = ip
4050                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4051                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
4052                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
4053                   sxbad(nclsx)  = nbadstrips(1,icl)
4054                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
4055                  
4056    
4057                 do is=1,2                 do is=1,2
4058  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4059  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4060                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)                    call xyz_PAM(icl,0,is,PFAdef,'    ',0.,0.,0.,0.)
4061                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7
4062                 enddo                 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)  
4063              else                          !=== Y views              else                          !=== Y views
4064                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4065                 planey(nclsy) = ip                 planey(nclsy) = ip
4066                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4067                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
4068                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
4069                   sybad(nclsy)  = nbadstrips(1,icl)
4070                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
4071    
4072    
4073                 do is=1,2                 do is=1,2
4074  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4075  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4076                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)                    call xyz_PAM(0,icl,is,'    ',PFAdef,0.,0.,0.,0.)
4077                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7
4078                 enddo                 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)  
4079              endif              endif
4080           endif           endif
4081    
# Line 3911  c$$$               print*,'ys(2,nclsy)   Line 4088  c$$$               print*,'ys(2,nclsy)  
4088  *     associati ad una traccia, e permettere di salvare  *     associati ad una traccia, e permettere di salvare
4089  *     solo questi nell'albero di uscita  *     solo questi nell'albero di uscita
4090  *     --------------------------------------------------  *     --------------------------------------------------
4091                            
   
 c$$$         print*,' cl ',icl,' --> ',cl_used(icl)  
 c$$$  
 c$$$         if( cl_used(icl).ne.0 )then  
 c$$$            if(  
 c$$$     $           mod(VIEW(icl),2).eq.0.and.  
 c$$$     $           cltrx(ip,whichtrack(icl)).ne.icl )  
 c$$$     $           print*,'**WARNING** cltrx(',ip,',',whichtrack(icl)  
 c$$$     $           ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl  
 c$$$            if(  
 c$$$     $           mod(VIEW(icl),2).eq.1.and.  
 c$$$     $           cltry(ip,whichtrack(icl)).ne.icl )  
 c$$$     $           print*,'**WARNING** cltry(',ip,',',whichtrack(icl)  
 c$$$     $           ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl  
 c$$$         endif  
           
   
4092        enddo        enddo
4093        end        end
4094    
# Line 3990  c$$$         endif Line 4150  c$$$         endif
4150                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4151                 nnn=nnn+ptcloud_yz(iyz)                 nnn=nnn+ptcloud_yz(iyz)
4152              enddo              enddo
4153              do ipt=1,nnn              do ipt=1,min(ndblt_max_nt,nnn)
4154                 db_cloud_nt(ipt)=db_cloud(ipt)                 db_cloud_nt(ipt)=db_cloud(ipt)
4155               enddo               enddo
4156           endif           endif
# Line 4003  c$$$         endif Line 4163  c$$$         endif
4163                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4164                 nnn=nnn+ptcloud_xz(ixz)                               nnn=nnn+ptcloud_xz(ixz)              
4165              enddo              enddo
4166              do ipt=1,nnn              do ipt=1,min(ntrpt_max_nt,nnn)
4167                tr_cloud_nt(ipt)=tr_cloud(ipt)                tr_cloud_nt(ipt)=tr_cloud(ipt)
4168               enddo               enddo
4169           endif           endif
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.45
  ViewVC Help
Powered by ViewVC 1.1.23