/[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.34 by pam-fi, Wed Mar 5 17:00:20 2008 UTC revision 1.39 by pam-fi, Mon Jan 26 14:01:41 2009 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(planehit.eq.3) goto 881                  if(nclouds_yz.eq.0)goto 880 !go to next event    
186    
187    
188    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  11111    continue               !<<<<<<< come here when performing a new search  11111    continue               !<<<<<<< come here when performing a new search
275    
276             if(nclouds_xz.eq.0)goto 880 !go to next event    
277             if(nclouds_yz.eq.0)goto 880 !go to next event    
278    
279  c         iflag=0  c         iflag=0
280           call clouds_to_ctrack(iflag)           call clouds_to_ctrack(iflag)
281           if(iflag.eq.1)then     !no candidate tracks found           if(iflag.eq.1)then     !no candidate tracks found
282              goto 880            !fill ntp and go to next event                goto 880            !fill ntp and go to next event  
283           endif           endif
284             if(ntracks.eq.0)goto 880 !go to next event    
285    
286           FIMAGE=.false.         !processing best track (not track image)           FIMAGE=.false.         !processing best track (not track image)
287           ibest=0                !best track among candidates           ibest=0                !best track among candidates
# Line 295  c$$$         if(ibest.eq.0)goto 880 !>> Line 325  c$$$         if(ibest.eq.0)goto 880 !>>
325             endif             endif
326           enddo           enddo
327    
 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  
328    
329           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
330  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
# Line 350  c$$$         enddo Line 362  c$$$         enddo
362           do i=1,5           do i=1,5
363              AL_GUESS(i)=AL(i)              AL_GUESS(i)=AL(i)
364           enddo           enddo
 c         print*,'## guess: ',al  
365    
366           do i=1,5           do i=1,5
367              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
# Line 365  c         if(DEBUG.EQ.1)iprint=1 Line 376  c         if(DEBUG.EQ.1)iprint=1
376           if(VERBOSE.EQ.1)iprint=1           if(VERBOSE.EQ.1)iprint=1
377           if(DEBUG.EQ.1)iprint=2           if(DEBUG.EQ.1)iprint=2
378           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
379           if(ifail.ne.0) then  cc         if(ifail.ne.0) then
380             if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
381                if(CHI2.ne.CHI2)CHI2=-9999. !new
382              if(VERBOSE.EQ.1)then              if(VERBOSE.EQ.1)then
383                 print *,                 print *,
384       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
385       $              ,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*,'----------------------------------------------'  
386              endif              endif
 c            chi2=-chi2  
387           endif           endif
388                    
389           if(DEBUG.EQ.1)then           if(DEBUG.EQ.1)then
# Line 424  c         rchi2=chi2/dble(ndof) Line 429  c         rchi2=chi2/dble(ndof)
429              if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'              if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
430              goto 122            !jump              goto 122            !jump
431           endif           endif
 c         print*,'is1 ',is1  
432           do ip=1,NPLANES           do ip=1,NPLANES
433              is2 = SENSOR_STORE(ip,icand) !sensor              is2 = SENSOR_STORE(ip,icand) !sensor
 c            print*,'is2 ',is2,' ip ',ip  
434              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
435              if(              if(
436       $           (is1.ne.is2.and.is2.ne.0)       $           (is1.ne.is2.and.is2.ne.0)
437       $           )goto 122      !jump (this track cannot have an image)       $           )goto 122      !jump (this track cannot have an image)
438           enddo           enddo
439           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  
440  *     ---------------------------------------------------------------  *     ---------------------------------------------------------------
441  *     take the candidate with the greatest number of matching couples  *     take the candidate with the greatest number of matching couples
442  *     if more than one satisfies the requirement,  *     if more than one satisfies the requirement,
# Line 486  c$$$         enddo Line 468  c$$$         enddo
468                                            
469                    endif                    endif
470                 endif                 endif
 c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)  
 c$$$     $              ,CP_STORE(nplanes-ip+1,icand),ncp  
471              enddo              enddo
472   117        continue   117        continue
473              trackimage(i)=ncp   !number of matching couples              trackimage(i)=ncp   !number of matching couples
# Line 501  c$$$     $              ,CP_STORE(nplane Line 481  c$$$     $              ,CP_STORE(nplane
481           do i=1,ntracks           do i=1,ntracks
482              if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1              if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
483           enddo           enddo
484             if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood
485  *     if there are, choose the best one  *     if there are, choose the best one
486           if(ngood.gt.1)then           if(ngood.gt.0)then
487  *     -------------------------------------------------------  *     -------------------------------------------------------
488  *     order track-candidates according to:  *     order track-candidates according to:
489  *     1st) decreasing n.points  *     1st) decreasing n.points
# Line 532  c$$$     $              ,CP_STORE(nplane Line 513  c$$$     $              ,CP_STORE(nplane
513                    endif                    endif
514                 endif                 endif
515              enddo              enddo
516    
517                if(DEBUG.EQ.1)then
518                   print*,'Track candidate ',iimage
519         $              ,' >>> TRACK IMAGE >>> of'
520         $              ,ibest
521                endif
522                            
523           endif           endif
524    
          if(DEBUG.EQ.1)print*,'Track candidate ',iimage  
      $        ,' >>> TRACK IMAGE >>> of'  
      $        ,ibest  
525    
526   122     continue   122     continue
527    
# Line 550  c$$$     $              ,CP_STORE(nplane Line 534  c$$$     $              ,CP_STORE(nplane
534       $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next       $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
535           if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous           if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
536           call fill_level2_tracks(ntrk)     !==> good2=.true.           call fill_level2_tracks(ntrk)     !==> good2=.true.
 c         print*,'++++++++++ iimage,fimage,ntrk,image '  
 c     $        ,iimage,fimage,ntrk,image(ntrk)  
537    
538           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
539              if(verbose.eq.1)              if(verbose.eq.1)
# Line 567  cc            good2=.false. Line 549  cc            good2=.false.
549              goto 1212           !>>> fit image-track              goto 1212           !>>> fit image-track
550           endif           endif
551    
 *     --- 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  
 *     **********************************************  
   
   
552    
553   880     return   880     return
554        end        end
# Line 699  c     $        rchi2best.lt.15..and. Line 649  c     $        rchi2best.lt.15..and.
649        parameter (ndivx=30)        parameter (ndivx=30)
650    
651    
 c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy  
652                
653        resxPAM = 0        resxPAM = 0
654        resyPAM = 0        resyPAM = 0
# Line 713  c$$$      print*,icx,icy,sensor,PFAx,PFA Line 662  c$$$      print*,icx,icy,sensor,PFAx,PFA
662        xPAM_B = 0.D0        xPAM_B = 0.D0
663        yPAM_B = 0.D0        yPAM_B = 0.D0
664        zPAM_B = 0.D0        zPAM_B = 0.D0
 cc      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy  
665    
666        if(sensor.lt.1.or.sensor.gt.2)then        if(sensor.lt.1.or.sensor.gt.2)then
667           print*,'xyz_PAM   ***ERROR*** wrong input '           print*,'xyz_PAM   ***ERROR*** wrong input '
# Line 805  c         resyPAM = RESYAV Line 753  c         resyPAM = RESYAV
753    
754   20   endif   20   endif
755    
 cc      print*,'## stripx,stripy ',stripx,stripy  
756    
757  c===========================================================  c===========================================================
758  C     COUPLE  C     COUPLE
# Line 885  C======================================= Line 832  C=======================================
832              nldx = nldy              nldx = nldy
833              viewx = viewy + 1              viewx = viewy + 1
834    
835              yi   = dcoordsi(stripy,viewy)              xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center
836                yi = dcoordsi(stripy,viewy)
837                zi = 0.D0
838    
839              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
840              yi_A = yi              yi_A = yi
# Line 895  C======================================= Line 844  C=======================================
844              yi_B = yi              yi_B = yi
845              zi_B = 0.              zi_B = 0.
846    
 c            print*,'Y-cl ',icy,stripy,' --> ',yi  
 c            print*,xi_A,' <--> ',xi_B  
847                            
848           elseif(icx.ne.0)then           elseif(icx.ne.0)then
849  c===========================================================  c===========================================================
# Line 907  C======================================= Line 854  C=======================================
854              nldy = nldx              nldy = nldx
855              viewy = viewx - 1              viewy = viewx - 1
856    
 c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx  
 c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  
857              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
858       $           .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...
859                 if(DEBUG.EQ.1) then                 if(DEBUG.EQ.1) then
# Line 916  c            if((stripx.le.3).or.(stripx Line 861  c            if((stripx.le.3).or.(stripx
861       $                 ' WARNING: false X strip: strip ',stripx       $                 ' WARNING: false X strip: strip ',stripx
862                 endif                 endif
863              endif              endif
864              xi   = dcoordsi(stripx,viewx)  
865                xi = dcoordsi(stripx,viewx)
866                yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
867                zi = 0.D0
868    
869              xi_A = xi              xi_A = xi
870              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 932  c            if((stripx.le.3).or.(stripx Line 880  c            if((stripx.le.3).or.(stripx
880                 yi_B = yi                 yi_B = yi
881              endif              endif
882    
 c            print*,'X-cl ',icx,stripx,' --> ',xi  
 c            print*,yi_A,' <--> ',yi_B  
883    
884           else           else
885              if(DEBUG.EQ.1) then              if(DEBUG.EQ.1) then
# Line 979  c     N.B. I convert angles from microra Line 925  c     N.B. I convert angles from microra
925       $        + zi_B       $        + zi_B
926       $        + dz(nplx,nldx,sensor)       $        + dz(nplx,nldx,sensor)
927    
928    
929    
930             xrt = xi
931         $        - omega(nplx,nldx,sensor)*yi
932         $        + gamma(nplx,nldx,sensor)*zi
933         $        + dx(nplx,nldx,sensor)
934            
935             yrt = omega(nplx,nldx,sensor)*xi
936         $        + yi
937         $        - beta(nplx,nldx,sensor)*zi
938         $        + dy(nplx,nldx,sensor)
939            
940             zrt = -gamma(nplx,nldx,sensor)*xi
941         $        + beta(nplx,nldx,sensor)*yi
942         $        + zi
943         $        + dz(nplx,nldx,sensor)
944    
945    
946                    
947  c      xrt = xi  c      xrt = xi
948  c      yrt = yi  c      yrt = yi
# Line 989  c     (xPAM,yPAM,zPAM) = measured coordi Line 953  c     (xPAM,yPAM,zPAM) = measured coordi
953  c                        in PAMELA reference system  c                        in PAMELA reference system
954  c------------------------------------------------------------------------  c------------------------------------------------------------------------
955    
956           xPAM = 0.D0           xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
957           yPAM = 0.D0           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
958           zPAM = 0.D0           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
959    c$$$         xPAM = 0.D0
960    c$$$         yPAM = 0.D0
961    c$$$         zPAM = 0.D0
962    
963           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
964           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 1002  c--------------------------------------- Line 969  c---------------------------------------
969           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
970                    
971    
 c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  
972    
973        else        else
974           if(DEBUG.EQ.1) then           if(DEBUG.EQ.1) then
# Line 1013  c         print*,'A-(',xPAM_A,yPAM_A,') Line 979  c         print*,'A-(',xPAM_A,yPAM_A,')
979        endif        endif
980                    
981    
 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  
982    
983   100  continue   100  continue
984        end        end
# Line 1063  c$$$      PFAy = 'COG4'!PFA Line 1026  c$$$      PFAy = 'COG4'!PFA
1026        call idtoc(pfaid,PFAx)        call idtoc(pfaid,PFAx)
1027        call idtoc(pfaid,PFAy)        call idtoc(pfaid,PFAy)
1028    
 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  
1029                
1030        if(icx.ne.0.and.icy.ne.0)then        if(icx.ne.0.and.icy.ne.0)then
1031    
1032           ipx=npl(VIEW(icx))           ipx=npl(VIEW(icx))
1033           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  
1034                    
1035           if( (nplanes-ipx+1).ne.ip )then           if( (nplanes-ipx+1).ne.ip )then
1036              print*,'xyzpam: ***WARNING*** cluster ',icx              print*,'xyzpam: ***WARNING*** cluster ',icx
# Line 1110  c         zv(ip) = zPAM Line 1065  c         zv(ip) = zPAM
1065        elseif(icx.eq.0.and.icy.ne.0)then        elseif(icx.eq.0.and.icy.ne.0)then
1066    
1067           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  
1068           if( (nplanes-ipy+1).ne.ip )then           if( (nplanes-ipy+1).ne.ip )then
1069              print*,'xyzpam: ***WARNING*** cluster ',icy              print*,'xyzpam: ***WARNING*** cluster ',icy
1070       $           ,' does not belong to plane: ',ip       $           ,' does not belong to plane: ',ip
# Line 1127  c$$$     $        ,' does not belong to Line 1079  c$$$     $        ,' does not belong to
1079           resx(ip)  = 1000.           resx(ip)  = 1000.
1080           resy(ip)  = resyPAM           resy(ip)  = resyPAM
1081    
1082           xm(ip) = -100.  c$$$         xm(ip) = -100.
1083           ym(ip) = -100.  c$$$         ym(ip) = -100.
1084           zm(ip) = (zPAM_A+zPAM_B)/2.  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1085             xm(ip) = xPAM
1086             ym(ip) = yPAM
1087             zm(ip) = zPAM
1088           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1089           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1090           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
# Line 1140  c         zv(ip) = (zPAM_A+zPAM_B)/2. Line 1095  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1095        elseif(icx.ne.0.and.icy.eq.0)then        elseif(icx.ne.0.and.icy.eq.0)then
1096    
1097           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  
1098    
1099           if( (nplanes-ipx+1).ne.ip )then           if( (nplanes-ipx+1).ne.ip )then
1100              print*,'xyzpam: ***WARNING*** cluster ',icx              print*,'xyzpam: ***WARNING*** cluster ',icx
# Line 1158  c$$$     $        ,' does not belong to Line 1110  c$$$     $        ,' does not belong to
1110           resx(ip)  = resxPAM           resx(ip)  = resxPAM
1111           resy(ip)  = 1000.           resy(ip)  = 1000.
1112    
1113           xm(ip) = -100.  c$$$         xm(ip) = -100.
1114           ym(ip) = -100.  c$$$         ym(ip) = -100.
1115           zm(ip) = (zPAM_A+zPAM_B)/2.  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1116             xm(ip) = xPAM
1117             ym(ip) = yPAM
1118             zm(ip) = zPAM
1119           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1120           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1121           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
# Line 1174  c         zv(ip) = (zPAM_A+zPAM_B)/2. Line 1129  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1129           if(lad.ne.0)il=lad           if(lad.ne.0)il=lad
1130           is = 1           is = 1
1131           if(sensor.ne.0)is=sensor           if(sensor.ne.0)is=sensor
 c         print*,nplanes-ip+1,il,is  
1132    
1133           xgood(ip) = 0.           xgood(ip) = 0.
1134           ygood(ip) = 0.           ygood(ip) = 0.
# Line 1194  c         zv(ip) = z_mech_sensor(nplanes Line 1148  c         zv(ip) = z_mech_sensor(nplanes
1148        endif        endif
1149    
1150        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
 c         print*,'----------------------------- track coord'  
1151  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1152           write(*,22222)ip,zm(ip),xm(ip),ym(ip)           write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1153       $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)       $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1154       $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)       $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
 c$$$         print*,'-----------------------------'  
1155        endif        endif
1156        end        end
1157    
# Line 1242  c$$$         print*,'------------------- Line 1194  c$$$         print*,'-------------------
1194    
1195  *     ----------------------  *     ----------------------
1196        if (        if (
1197       +     xPAM.eq.0.and.  c$$$     +     xPAM.eq.0.and.
1198       +     yPAM.eq.0.and.  c$$$     +     yPAM.eq.0.and.
1199       +     zPAM.eq.0.and.  c$$$     +     zPAM.eq.0.and.
1200       +     xPAM_A.ne.0.and.       +     xPAM_A.ne.0.and.
1201       +     yPAM_A.ne.0.and.       +     yPAM_A.ne.0.and.
1202       +     zPAM_A.ne.0.and.       +     zPAM_A.ne.0.and.
# Line 1286  c$$$         print*,'------------------- Line 1238  c$$$         print*,'-------------------
1238  cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2  cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1239           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1240    
 c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',re  
1241    
1242                    
1243  *     ----------------------  *     ----------------------
1244        elseif(        elseif(
1245       +     xPAM.ne.0.and.  c$$$     +     xPAM.ne.0.and.
1246       +     yPAM.ne.0.and.  c$$$     +     yPAM.ne.0.and.
1247       +     zPAM.ne.0.and.  c$$$     +     zPAM.ne.0.and.
1248       +     xPAM_A.eq.0.and.       +     xPAM_A.eq.0.and.
1249       +     yPAM_A.eq.0.and.       +     yPAM_A.eq.0.and.
1250       +     zPAM_A.eq.0.and.       +     zPAM_A.eq.0.and.
# Line 1316  c$$$     $        + Line 1265  c$$$     $        +
1265  c$$$     $        ((yPAM-YPP)/resyPAM)**2  c$$$     $        ((yPAM-YPP)/resyPAM)**2
1266           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1267    
 c$$$         print*,xPAM,yPAM,zPAM  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',resxPAM,resyPAM  
1268                    
1269        else        else
1270                    
 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  
1271        endif          endif  
1272    
1273        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1391  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1332  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1332  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1333  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1334  c------------------------------------------------------------------------  c------------------------------------------------------------------------
 c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)  
 c     $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...  
 c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  
 c                  print*,'whichsensor: ',  
 c     $                ' WARNING: false X strip: strip ',stripx  
 c               endif  
1335                 xi = dcoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1336                 yi = dcoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1337                 zi = 0.D0                 zi = 0.D0
# Line 1424  c--------------------------------------- Line 1359  c---------------------------------------
1359    
1360                 yvv(iv)=sngl(yvvv)                 yvv(iv)=sngl(yvvv)
1361                 xvv(iv)=sngl(xvvv)                 xvv(iv)=sngl(xvvv)
 c               print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '  
 c     $              ,iv,xvv(iv),yvv(iv)  
1362              enddo               !end loop on sensor vertexes              enddo               !end loop on sensor vertexes
1363    
1364              dtot=0.              dtot=0.
# Line 1550  c      include 'common_analysis.f' Line 1483  c      include 'common_analysis.f'
1483        is_cp=0        is_cp=0
1484        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1485        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 !!!'  
1486    
1487        return        return
1488        end        end
# Line 1652  c      include 'level1.f' Line 1584  c      include 'level1.f'
1584            
1585        if(DEBUG.EQ.1)print*,'cl_to_couples:'        if(DEBUG.EQ.1)print*,'cl_to_couples:'
1586    
1587    cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
1588    
1589  *     init variables  *     init variables
       ncp_tot=0  
1590        do ip=1,nplanes        do ip=1,nplanes
1591           do ico=1,ncouplemax           do ico=1,ncouplemax
1592              clx(ip,ico)=0              clx(ip,ico)=0
# Line 1666  c      include 'level1.f' Line 1599  c      include 'level1.f'
1599           ncls(ip)=0           ncls(ip)=0
1600        enddo        enddo
1601        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1602           cl_single(icl) = 1           cl_single(icl) = 1     !all are single
1603           cl_good(icl)   = 0           cl_good(icl)   = 0     !all are bad
1604        enddo        enddo
1605        do iv=1,nviews        do iv=1,nviews
1606           ncl_view(iv)  = 0           ncl_view(iv)  = 0
# Line 1695  c            mask_view(iv) = 1 Line 1628  c            mask_view(iv) = 1
1628        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1629           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1630                    
1631             if(cl_used(icx).ne.0)goto 10
1632    
1633  *     ----------------------------------------------------  *     ----------------------------------------------------
1634  *     jump masked views (X VIEW)  *     jump masked views (X VIEW)
1635  *     ----------------------------------------------------  *     ----------------------------------------------------
# Line 1702  c            mask_view(iv) = 1 Line 1637  c            mask_view(iv) = 1
1637  *     ----------------------------------------------------  *     ----------------------------------------------------
1638  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1639  *     ----------------------------------------------------  *     ----------------------------------------------------
1640           if(sgnl(icx).lt.dedx_x_min)then           if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
1641              cl_single(icx)=0              cl_single(icx)=0
1642              goto 10              goto 10
1643           endif           endif
# Line 1743  c     endif Line 1678  c     endif
1678                    
1679           do icy=1,nclstr1       !loop on cluster (Y)           do icy=1,nclstr1       !loop on cluster (Y)
1680              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1681    
1682                if(cl_used(icx).ne.0)goto 20
1683                            
1684  *     ----------------------------------------------------  *     ----------------------------------------------------
1685  *     jump masked views (Y VIEW)  *     jump masked views (Y VIEW)
# Line 1752  c     endif Line 1689  c     endif
1689  *     ----------------------------------------------------  *     ----------------------------------------------------
1690  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1691  *     ----------------------------------------------------  *     ----------------------------------------------------
1692              if(sgnl(icy).lt.dedx_y_min)then              if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
1693                 cl_single(icy)=0                 cl_single(icy)=0
1694                 goto 20                 goto 20
1695              endif              endif
# Line 1851  c                  mask_view(nviewy(nply Line 1788  c                  mask_view(nviewy(nply
1788   10      continue   10      continue
1789        enddo                     !end loop on clusters(X)        enddo                     !end loop on clusters(X)
1790                
         
1791        do icl=1,nclstr1        do icl=1,nclstr1
1792           if(cl_single(icl).eq.1)then           if(cl_single(icl).eq.1)then
1793              ip=npl(VIEW(icl))              ip=npl(VIEW(icl))
# Line 1859  c                  mask_view(nviewy(nply Line 1795  c                  mask_view(nviewy(nply
1795              cls(ip,ncls(ip))=icl              cls(ip,ncls(ip))=icl
1796           endif           endif
1797        enddo        enddo
1798    
1799     80   continue
1800                
1801                
1802        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
1803           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1804           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1805             print*,'used    ',(cl_used(i),i=1,nclstr1)
1806           print*,'singlets',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
1807           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1808        endif        endif
1809    
1810      
1811          if(.not.RECOVER_SINGLETS)goto 81
1812    
1813    *     ////////////////////////////////////////////////
1814    *     PATCH to recover events with less than 3 couples
1815    *     ////////////////////////////////////////////////    
1816    *     loop over singlet and create "fake" couples
1817    *     (with clx=0 or cly=0)
1818    *    
1819    
1820          if(DEBUG.EQ.1)
1821         $     print*,'>>> Recover singlets '
1822         $     ,'(creates fake couples) <<<'
1823          do icl=1,nclstr1
1824             if(
1825         $        cl_single(icl).eq.1.and.
1826         $        cl_used(icl).eq.0.and.
1827         $        .true.)then
1828    *     ----------------------------------------------------
1829    *     jump masked views
1830    *     ----------------------------------------------------
1831                if( mask_view(VIEW(icl)).ne.0 ) goto 21
1832                if(mod(VIEW(icl),2).eq.0)then !=== X views
1833    *     ----------------------------------------------------
1834    *     cut on charge (X VIEW)
1835    *     ----------------------------------------------------
1836                   if(sgnl(icl).lt.dedx_x_min) goto 21
1837                  
1838                   nplx=npl(VIEW(icl))
1839    *     ------------------> (FAKE) COUPLE <-----------------
1840                   ncp_plane(nplx) = ncp_plane(nplx) + 1
1841                   clx(nplx,ncp_plane(nplx))=icl
1842                   cly(nplx,ncp_plane(nplx))=0
1843    c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1844    *     ----------------------------------------------------
1845                  
1846                else                !=== Y views
1847    *     ----------------------------------------------------
1848    *     cut on charge (Y VIEW)
1849    *     ----------------------------------------------------
1850                   if(sgnl(icl).lt.dedx_y_min) goto 21
1851                  
1852                   nply=npl(VIEW(icl))
1853    *     ------------------> (FAKE) COUPLE <-----------------
1854                   ncp_plane(nply) = ncp_plane(nply) + 1
1855                   clx(nply,ncp_plane(nply))=0
1856                   cly(nply,ncp_plane(nply))=icl
1857    c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1858    *     ----------------------------------------------------
1859                  
1860                endif
1861             endif                  !end "single" condition
1862     21      continue
1863          enddo                     !end loop over clusters
1864    
1865          if(DEBUG.EQ.1)
1866         $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1867    
1868    
1869     81   continue
1870                
1871        do ip=1,6        ncp_tot=0
1872          do ip=1,NPLANES
1873           ncp_tot = ncp_tot + ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1874        enddo        enddo
1875          if(DEBUG.EQ.1)
1876         $     print*,'n.couple tot:      ',ncp_tot
1877                
1878        return        return
1879        end        end
# Line 1933  c      double precision xm3,ym3,zm3 Line 1936  c      double precision xm3,ym3,zm3
1936  *     --------------------------------------------  *     --------------------------------------------
1937        do ip=1,nplanes        do ip=1,nplanes
1938           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  
1939              mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7              mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1940              mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7              mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1941           endif           endif
# Line 1945  c            mask_view(nviewy(ip)) = 8 Line 1946  c            mask_view(nviewy(ip)) = 8
1946        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1947                
1948        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1949    c$$$         print*,'(1) ip ',ip1
1950    c$$$     $        ,mask_view(nviewx(ip1))
1951    c$$$     $        ,mask_view(nviewy(ip1))        
1952           if(  mask_view(nviewx(ip1)).ne.0 .or.           if(  mask_view(nviewx(ip1)).ne.0 .or.
1953       $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane       $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1954           do is1=1,2             !loop on sensors - COPPIA 1                       do is1=1,2             !loop on sensors - COPPIA 1            
1955              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1956                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1957                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1958                  
1959    c$$$               print*,'(1) ip ',ip1,' icp ',icp1
1960    
1961  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1962  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1963                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1964                 xm1=xPAM                 xm1=xPAM
1965                 ym1=yPAM                 ym1=yPAM
1966                 zm1=zPAM                                   zm1=zPAM                  
 c     print*,'***',is1,xm1,ym1,zm1  
1967    
1968                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1969    c$$$                  print*,'(2) ip ',ip2
1970    c$$$     $                 ,mask_view(nviewx(ip2))
1971    c$$$     $                 ,mask_view(nviewy(ip2))
1972                    if(  mask_view(nviewx(ip2)).ne.0 .or.                    if(  mask_view(nviewx(ip2)).ne.0 .or.
1973       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1974                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                                        do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
1975                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1976                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1977                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1978    
1979    c$$$                        print*,'(2) ip ',ip2,' icp ',icp2
1980    
1981  c                        call xyz_PAM  c                        call xyz_PAM
1982  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1983  c                        call xyz_PAM  c                        call xyz_PAM
# Line 1975  c     $                       (icx2,icy2 Line 1987  c     $                       (icx2,icy2
1987                          xm2=xPAM                          xm2=xPAM
1988                          ym2=yPAM                          ym2=yPAM
1989                          zm2=zPAM                          zm2=zPAM
1990                                                    
1991    *                       ---------------------------------------------------
1992    *                       both couples must have a y-cluster
1993    *                       (condition necessary when in RECOVER_SINGLETS mode)
1994    *                       ---------------------------------------------------
1995                            if(icy1.eq.0.or.icy2.eq.0)goto 111
1996    
1997                            if(cl_used(icy1).ne.0)goto 111
1998                            if(cl_used(icy2).ne.0)goto 111
1999    
2000                            
2001  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2002  *     track parameters on Y VIEW  *     track parameters on Y VIEW
2003  *     (2 couples needed)  *     (2 couples needed)
# Line 1985  c     $                       (icx2,icy2 Line 2007  c     $                       (icx2,icy2
2007       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2008       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2009       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2010  c                           good2=.false.  c     good2=.false.
2011  c                           goto 880 !fill ntp and go to next event  c     goto 880 !fill ntp and go to next event
2012                             do iv=1,12                             do iv=1,12
2013  c                              mask_view(iv) = 3  c     mask_view(iv) = 3
2014                                mask_view(iv) = mask_view(iv)+ 2**2                                mask_view(iv) = mask_view(iv)+ 2**2
2015                             enddo                             enddo
2016                             iflag=1                             iflag=1
2017                             return                             return
2018                          endif                          endif
2019                            
2020                            
2021    ccc                        print*,'<doublet> ',icp1,icp2
2022    
2023                          ndblt = ndblt + 1                          ndblt = ndblt + 1
2024  *     store doublet info  *     store doublet info
2025                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
# Line 2002  c                              mask_view Line 2028  c                              mask_view
2028                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2029  *     y0 (cm)  *     y0 (cm)
2030                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
2031                                                      
2032  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2033  ****  reject non phisical couples                    ****  ****  reject non phisical couples                    ****
2034  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2035                            if(SECOND_SEARCH)goto 111
2036                          if(                          if(
2037       $                       abs(alfayz2(ndblt)).gt.alfyz2_max       $                       abs(alfayz2(ndblt)).gt.alfyz2_max
2038       $                       .or.       $                       .or.
2039       $                       abs(alfayz1(ndblt)).gt.alfyz1_max       $                       abs(alfayz1(ndblt)).gt.alfyz1_max
2040       $                       )ndblt = ndblt-1       $                       )ndblt = ndblt-1
2041                                                    
2042  c$$$      if(iev.eq.33)then                          
2043  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$$$  
2044  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2045  *     track parameters on Y VIEW - end  *     track parameters on Y VIEW - end
2046  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2047    
2048    
2049                            if(icx1.ne.0)then
2050                               if(cl_used(icx1).ne.0)goto 31
2051                            endif
2052                            if(icx2.ne.0)then
2053                               if(cl_used(icx2).ne.0)goto 31
2054                            endif
2055    
2056                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2057    
2058                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2059    c$$$                           print*,'(3) ip ',ip3
2060    c$$$     $                          ,mask_view(nviewx(ip3))
2061    c$$$     $                          ,mask_view(nviewy(ip3))                          
2062                             if(  mask_view(nviewx(ip3)).ne.0 .or.                             if(  mask_view(nviewx(ip3)).ne.0 .or.
2063       $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane       $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2064                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
# Line 2034  c$$$ Line 2066  c$$$
2066                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
2067                                   icx3=clx(ip3,icp3)                                   icx3=clx(ip3,icp3)
2068                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2069    
2070    c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
2071    
2072    *     ---------------------------------------------------
2073    *     all three couples must have a x-cluster
2074    *     (condition necessary when in RECOVER_SINGLETS mode)
2075    *     ---------------------------------------------------
2076                                     if(
2077         $                                icx1.eq.0.or.
2078         $                                icx2.eq.0.or.
2079         $                                icx3.eq.0.or.
2080         $                                .false.)goto 29
2081                                    
2082                                     if(cl_used(icx1).ne.0)goto 29
2083                                     if(cl_used(icx2).ne.0)goto 29
2084                                     if(cl_used(icx3).ne.0)goto 29
2085    
2086  c                                 call xyz_PAM  c                                 call xyz_PAM
2087  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2088  c                                 call xyz_PAM  c                                 call xyz_PAM
# Line 2044  c     $                               (i Line 2093  c     $                               (i
2093                                   xm3=xPAM                                   xm3=xPAM
2094                                   ym3=yPAM                                   ym3=yPAM
2095                                   zm3=zPAM                                   zm3=zPAM
2096    
2097    
2098  *     find the circle passing through the three points  *     find the circle passing through the three points
2099  c$$$                                 call tricircle(3,xp,zp,angp,resp,chi                                   iflag_t = DEBUG
 c$$$     $                                ,xc,zc,radius,iflag)  
                                  iflag = DEBUG  
2100                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2101       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag_t)
 c     print*,xc,zc,radius  
2102  *     the circle must intersect the reference plane  *     the circle must intersect the reference plane
2103                                   if(  cc                                 if(iflag.ne.0)goto 29
2104  c     $                                 (xc.le.-1.*xclimit.or.                                   if(iflag_t.ne.0)then
2105  c     $                                 xc.ge.xclimit).and.  *     if tricircle fails, evaluate a straight line
2106       $                                radius**2.ge.(ZINI-zc)**2.and.                                      if(DEBUG.eq.1)
2107       $                                iflag.eq.0.and.       $                                   print*,'TRICIRCLE failure'
2108       $                                .true.)then       $                                   ,' >>> straight line'
2109                                                                        radius=0.
2110                                        xc=0.
2111                                        yc=0.
2112                                        
2113                                        SZZ=0.                  
2114                                        SZX=0.
2115                                        SSX=0.
2116                                        SZ=0.
2117                                        S1=0.
2118                                        X0=0.
2119                                        Ax=0.
2120                                        BX=0.
2121                                        DO I=1,3
2122                                           XX = XP(I)
2123                                           SZZ=SZZ+ZP(I)*ZP(I)
2124                                           SZX=SZX+ZP(I)*XX
2125                                           SSX=SSX+XX
2126                                           SZ=SZ+ZP(I)
2127                                           S1=S1+1.                                    
2128                                        ENDDO
2129                                        DET=SZZ*S1-SZ*SZ
2130                                        AX=(SZX*S1-SZ*SSX)/DET
2131                                        BX=(SZZ*SSX-SZX*SZ)/DET
2132                                        X0  = AX*ZINI+BX
2133                                        
2134                                     endif
2135    
2136                                     if(  .not.SECOND_SEARCH.and.
2137         $                                radius**2.lt.(ZINI-zc)**2)goto 29
2138                                                                      
2139  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2140  *     track parameters on X VIEW  *     track parameters on X VIEW
2141  *     (3 couples needed)  *     (3 couples needed)
2142  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2143                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2144                                      if(verbose.eq.1)print*,                                      if(verbose.eq.1)print*,
2145       $                     '** warning ** number of identified '//       $                                   '** warning **'//
2146       $                     'triplets exceeds vector dimention '       $                                   ' number of identified '//
2147       $                    ,'( ',ntrpt_max,' )'       $                                   'triplets exceeds'//
2148  c                                    good2=.false.       $                                   ' vector dimention '
2149  c                                    goto 880 !fill ntp and go to next event       $                                   ,'( ',ntrpt_max,' )'
2150    c     good2=.false.
2151    c     goto 880 !fill ntp and go to next event
2152                                      do iv=1,nviews                                      do iv=1,nviews
2153  c                                       mask_view(iv) = 4  c     mask_view(iv) = 4
2154                                         mask_view(iv)=mask_view(iv)+ 2**3                                         mask_view(iv) =
2155         $                                      mask_view(iv)+ 2**3
2156                                      enddo                                      enddo
2157                                      iflag=1                                      iflag=1
2158                                      return                                      return
2159                                   endif                                   endif
2160    
2161    ccc                                 print*,'<triplet> ',icp1,icp2,icp3
2162                                    
2163                                   ntrpt = ntrpt +1                                   ntrpt = ntrpt +1
2164  *     store triplet info  *     store triplet info
2165                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
2166                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
2167                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
2168                                                                    
2169                                   if(xc.lt.0)then                                   if(radius.ne.0.and.xc.lt.0)then
2170  *************POSITIVE DEFLECTION  *************POSITIVE DEFLECTION
2171                alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)               alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
2172                alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)               alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2173                alfaxz3(ntrpt) = 1/radius               alfaxz3(ntrpt) = 1/radius
2174                                   else                                  else if(radius.ne.0.and.xc.ge.0)then
2175  *************NEGATIVE DEFLECTION  *************NEGATIVE DEFLECTION
2176                alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)               alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
2177                alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)               alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2178                alfaxz3(ntrpt) = -1/radius               alfaxz3(ntrpt) = -1/radius
2179                                   endif                                  else if(radius.eq.0)then
2180                                    *************straight fit
2181                 alfaxz1(ntrpt) = X0
2182                 alfaxz2(ntrpt) = AX
2183                 alfaxz3(ntrpt) = 0.
2184                                    endif
2185    
2186    c$$$                                print*,'alfaxz1 ', alfaxz1(ntrpt)
2187    c$$$                                print*,'alfaxz2 ', alfaxz2(ntrpt)
2188    c$$$                                print*,'alfaxz3 ', alfaxz3(ntrpt)
2189                                        
2190  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2191  ****  reject non phisical triplets                   ****  ****  reject non phisical triplets                   ****
2192  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2193                                   if(                                  if(SECOND_SEARCH)goto 29
2194       $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max                                  if(
2195       $                                .or.       $                               abs(alfaxz2(ntrpt)).gt.
2196       $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max       $                               alfxz2_max
2197       $                                )ntrpt = ntrpt-1       $                               .or.
2198                                         $                               abs(alfaxz1(ntrpt)).gt.
2199                                         $                               alfxz1_max
2200  c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)       $                               )ntrpt = ntrpt-1
2201                                                                    
2202  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2203  *     track parameters on X VIEW - end  *     track parameters on X VIEW - end
2204  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2205                                endif                                
2206     29                           continue
2207                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2208                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2209   30                     continue   30                     continue
2210                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2211   31                  continue  
2212                         31                  continue                    
2213   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
2214                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2215   20            continue   20            continue
2216              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2217                
2218     11         continue
2219           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2220        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2221   10   continue   10   continue
# Line 2192  c      include 'momanhough_init.f' Line 2287  c      include 'momanhough_init.f'
2287        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2288           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
2289                            
 c     print*,'--------------'  
 c     print*,'** ',idb1,' **'  
               
2290           do icp=1,ncp_tot           do icp=1,ncp_tot
2291              cp_useds1(icp)=0    !init              cp_useds1(icp)=0    !init
2292              cp_useds2(icp)=0    !init              cp_useds2(icp)=0    !init
# Line 2233  ccccc if(db_used(idbref).eq.1)goto 1188 Line 2325  ccccc if(db_used(idbref).eq.1)goto 1188
2325       $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2                     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2              
2326                 distance = sqrt(distance)                 distance = sqrt(distance)
2327                                
 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  
2328                 if(distance.lt.cutdistyz)then                 if(distance.lt.cutdistyz)then
2329    
 c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  
2330                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
2331                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
2332                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
# Line 2258  c     print*,idb1,idb2,distance,' cloud Line 2342  c     print*,idb1,idb2,distance,' cloud
2342    
2343                    temp1 = temp1 + alfayz1(idb2)                    temp1 = temp1 + alfayz1(idb2)
2344                    temp2 = temp2 + alfayz2(idb2)                    temp2 = temp2 + alfayz2(idb2)
 c     print*,'*   idbref,idb2 ',idbref,idb2  
2345                 endif                               endif              
2346                                
2347   1118          continue   1118          continue
# Line 2281  c     print*,'*   idbref,idb2 ',idbref,i Line 2364  c     print*,'*   idbref,idb2 ',idbref,i
2364           enddo           enddo
2365           ncpused=0           ncpused=0
2366           do icp=1,ncp_tot           do icp=1,ncp_tot
2367              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2368         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2369         $           .true.)then
2370                 ncpused=ncpused+1                 ncpused=ncpused+1
2371                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2372                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2291  c     print*,'*   idbref,idb2 ',idbref,i Line 2376  c     print*,'*   idbref,idb2 ',idbref,i
2376           do ip=1,nplanes           do ip=1,nplanes
2377              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2378           enddo           enddo
2379  c     print*,'>>>> ',ncpused,npt,nplused          
 c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  
          if(npt.lt.nptyz_min)goto 2228 !next doublet  
2380           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2381                    
2382  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2325  c               mask_view(iv) = 5 Line 2408  c               mask_view(iv) = 5
2408  c     ptcloud_yz_nt(nclouds_yz)=npt  c     ptcloud_yz_nt(nclouds_yz)=npt
2409           do ipt=1,npt           do ipt=1,npt
2410              db_cloud(npt_tot+ipt) = db_all(ipt)              db_cloud(npt_tot+ipt) = db_all(ipt)
 c     print*,'>> ',ipt,db_all(ipt)  
2411           enddo             enddo  
2412           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2413           if(DEBUG.EQ.1)then           if(DEBUG.EQ.1)then
             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'  
2414              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2415              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2416              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
# Line 2338  c     print*,'>> ',ipt,db_all(ipt) Line 2419  c     print*,'>> ',ipt,db_all(ipt)
2419              print*,'cpcloud_yz '              print*,'cpcloud_yz '
2420       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2421              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)  
2422           endif           endif
2423  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2424  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2356  c$$$     $           ,(db_cloud(iii),iii Line 2433  c$$$     $           ,(db_cloud(iii),iii
2433        endif                            endif                    
2434                
2435        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
          print*,'---------------------- '  
2436           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
          print*,' '  
2437        endif        endif
2438                
2439                
# Line 2421  c      include 'momanhough_init.f' Line 2496  c      include 'momanhough_init.f'
2496   91   continue                     91   continue                  
2497        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2498           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,' **'  
2499                    
2500           do icp=1,ncp_tot           do icp=1,ncp_tot
2501              cp_useds1(icp)=0              cp_useds1(icp)=0
# Line 2456  c         tr_temp(1)=itr1 Line 2529  c         tr_temp(1)=itr1
2529              do itr2=1,ntrpt     !loop (2) on TRIPLETS              do itr2=1,ntrpt     !loop (2) on TRIPLETS
2530                 if(itr2.eq.itr1)goto 11188       !next triplet                 if(itr2.eq.itr1)goto 11188       !next triplet
2531                 if(tr_used(itr2).eq.1)goto 11188 !next triplet                               if(tr_used(itr2).eq.1)goto 11188 !next triplet              
2532    
2533    
2534  *     triplet distance in parameter space  *     triplet distance in parameter space
2535  *     solo i due parametri spaziali per il momemnto  *     solo i due parametri spaziali per il momemnto
2536                 distance=                 distance=
2537       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2538       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2539                 distance = sqrt(distance)                 distance = sqrt(distance)
2540                  
2541    
2542  *     ------------------------------------------------------------------------  *     ------------------------------------------------------------------------
2543  *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE  *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2544  *     ------------------------------------------------------------------------  *     ------------------------------------------------------------------------
# Line 2475  c         tr_temp(1)=itr1 Line 2551  c         tr_temp(1)=itr1
2551       $              .true.)istrimage=1       $              .true.)istrimage=1
2552    
2553                 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  
2554                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2555                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
2556                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
# Line 2494  c     print*,idb1,idb2,distance,' cloud Line 2569  c     print*,idb1,idb2,distance,' cloud
2569                    temp1 = temp1 + alfaxz1(itr2)                    temp1 = temp1 + alfaxz1(itr2)
2570                    temp2 = temp2 + alfaxz2(itr2)                    temp2 = temp2 + alfaxz2(itr2)
2571                    temp3 = temp3 + alfaxz3(itr2)                    temp3 = temp3 + alfaxz3(itr2)
 c     print*,'*   itrref,itr2 ',itrref,itr2,distance  
2572                 endif                               endif              
2573                                
2574  11188          continue  11188          continue
# Line 2515  c     print*,'*   itrref,itr2 ',itrref,i Line 2589  c     print*,'*   itrref,itr2 ',itrref,i
2589  *     1bis)  *     1bis)
2590  *     2) it is not already stored  *     2) it is not already stored
2591  *     ------------------------------------------  *     ------------------------------------------
 c     print*,'check cp_used'  
2592           do ip=1,nplanes           do ip=1,nplanes
2593              hit_plane(ip)=0              hit_plane(ip)=0
2594           enddo           enddo
2595           ncpused=0           ncpused=0
2596           do icp=1,ncp_tot           do icp=1,ncp_tot
2597              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2598         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2599         $           .true.)then
2600                 ncpused=ncpused+1                 ncpused=ncpused+1
2601                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2602                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2531  c     print*,'check cp_used' Line 2606  c     print*,'check cp_used'
2606           do ip=1,nplanes           do ip=1,nplanes
2607              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2608           enddo           enddo
 c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  
          if(npt.lt.nptxz_min)goto 22288     !next triplet  
2609           if(nplused.lt.nplxz_min)goto 22288 !next triplet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2610                    
2611  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2566  c               mask_view(iv) = 6 Line 2639  c               mask_view(iv) = 6
2639           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2640                    
2641           if(DEBUG.EQ.1)then           if(DEBUG.EQ.1)then
             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'  
2642              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2643              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2644              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
# Line 2576  c               mask_view(iv) = 6 Line 2648  c               mask_view(iv) = 6
2648              print*,'cpcloud_xz '              print*,'cpcloud_xz '
2649       $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)       $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2650              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)  
2651           endif           endif
2652  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2653  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2593  c$$$     $           ,(tr_cloud(iii),iii Line 2661  c$$$     $           ,(tr_cloud(iii),iii
2661         endif                             endif                    
2662                
2663        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
          print*,'---------------------- '  
2664           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
          print*,' '  
2665        endif        endif
2666                
2667                
# Line 2651  c$$$     $           ,(tr_cloud(iii),iii Line 2717  c$$$     $           ,(tr_cloud(iii),iii
2717    
2718        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
2719                
2720        do iyz=1,nclouds_yz       !loop on YZ couds        do iyz=1,nclouds_yz       !loop on YZ clouds
2721           do ixz=1,nclouds_xz    !loop on XZ couds           do ixz=1,nclouds_xz    !loop on XZ clouds
2722                            
2723  *     --------------------------------------------------  *     --------------------------------------------------
2724  *     check of consistency of the clouds  *     check of consistency of the clouds
# Line 2661  c$$$     $           ,(tr_cloud(iii),iii Line 2727  c$$$     $           ,(tr_cloud(iii),iii
2727  *     of the two clouds  *     of the two clouds
2728  *     --------------------------------------------------  *     --------------------------------------------------
2729              do ip=1,nplanes              do ip=1,nplanes
2730                 hit_plane(ip)=0                 hit_plane(ip)=0  !n.matching couples (REAL couples, not SINGLETS)
2731                 ncp_match(ip)=0                 ncp_match(ip)=0  !n.matching couples per plane
2732                 do icpp=1,ncouplemax                 do icpp=1,ncouplemax
2733                    cp_match(ip,icpp)=0 !init couple list                    cp_match(ip,icpp)=0 !init couple list
2734                 enddo                 enddo
2735              enddo              enddo
2736              ncp_ok=0              ncp_ok=0            !count n.matching-couples
2737                ncpx_ok=0           !count n.matching-couples with x cluster
2738                ncpy_ok=0           !count n.matching-couples with y cluster
2739    
2740    
2741              do icp=1,ncp_tot    !loop over couples              do icp=1,ncp_tot    !loop over couples
2742  *     get info on  
2743                 cpintersec(icp)=min(                 if(.not.RECOVER_SINGLETS)then
2744       $              cpcloud_yz(iyz,icp),  *     ------------------------------------------------------
2745       $              cpcloud_xz(ixz,icp))  *     if NOT in RECOVER_SINGLETS mode, take the intersection
2746                 if(  *     between xz yz clouds
2747       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.  *     ------------------------------------------------------
2748       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.                    cpintersec(icp)=min(
2749       $              .false.)cpintersec(icp)=0       $                 cpcloud_yz(iyz,icp),
2750         $                 cpcloud_xz(ixz,icp))
2751  *     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
2752    *     ------------------------------------------------------
2753    *     discard the couple if the sensor is in conflict
2754    *     ------------------------------------------------------
2755                      if(
2756         $       (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2757         $       (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2758         $                 .false.)cpintersec(icp)=0
2759                   else
2760    *     ------------------------------------------------------
2761    *     if RECOVER_SINGLETS take the union
2762    *     (otherwise the fake couples formed by singlets would be
2763    *     discarded)    
2764    *     ------------------------------------------------------
2765                      cpintersec(icp)=max(
2766         $                 cpcloud_yz(iyz,icp),
2767         $                 cpcloud_xz(ixz,icp))
2768    c$$$                  if(cpcloud_yz(iyz,icp).gt.0)
2769    c$$$     $                 cpintersec(icp)=cpcloud_yz(iyz,icp)
2770    *     cpintersec is >0 if either yz or xz cloud contains the couple icp
2771                   endif
2772    
2773    c$$$               print*,icp,ip_cp(icp),' -- ',cpintersec(icp)
2774    
2775                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
                   ncp_ok=ncp_ok+1    
2776                                        
2777                    ip=ip_cp(icp)                    ip=ip_cp(icp)
2778                    hit_plane(ip)=1                    hit_plane(ip)=1
2779    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0)
2780    c$$$     $                 ncp_ok=ncp_ok+1  
2781    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0)
2782    c$$$     $                 ncpx_ok=ncpx_ok+1
2783    c$$$                  if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0)
2784    c$$$     $                 ncpy_ok=ncpy_ok+1
2785    
2786                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2787         $                 cpcloud_xz(ixz,icp).gt.0)
2788         $                 ncp_ok=ncp_ok+1  
2789                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2790         $                 cpcloud_xz(ixz,icp).eq.0)
2791         $                 ncpy_ok=ncpy_ok+1  
2792                      if(  cpcloud_yz(iyz,icp).eq.0.and.
2793         $                 cpcloud_xz(ixz,icp).gt.0)
2794         $                 ncpx_ok=ncpx_ok+1  
2795    
2796                    if(cpintersec(icp).eq.1)then                    if(cpintersec(icp).eq.1)then
2797  *     1) only the couple image in sensor 1 matches  *     1) only the couple image in sensor 1 matches
2798                       id=-icp                       id=-icp
# Line 2709  c$$$     $           ,(tr_cloud(iii),iii Line 2819  c$$$     $           ,(tr_cloud(iii),iii
2819              do ip=1,nplanes              do ip=1,nplanes
2820                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2821              enddo              enddo
2822                                        
2823                            if(nplused.lt.3)goto 888 !next combination
2824    ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
2825    ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
2826    *     -----------------------------------------------------------
2827    *     if in RECOVER_SINGLET mode, the two clouds must have
2828    *     at least ONE intersecting real couple
2829    *     -----------------------------------------------------------
2830                if(ncp_ok.lt.1)goto 888 !next combination
2831    
2832              if(DEBUG.EQ.1)then              if(DEBUG.EQ.1)then
2833                 print*,'Combination ',iyz,ixz                 print*,'////////////////////////////'
2834       $              ,' db ',ptcloud_yz(iyz)                 print*,'Cloud combination (Y,X): ',iyz,ixz
2835       $              ,' tr ',ptcloud_xz(ixz)                 print*,' db ',ptcloud_yz(iyz)
2836       $              ,'  -----> # matching couples ',ncp_ok                 print*,' tr ',ptcloud_xz(ixz)
2837              endif                 print*,'  -----> # matching couples ',ncp_ok
2838                   print*,'  -----> # fake couples (X)',ncpx_ok
2839  c            if(nplused.lt.nplxz_min)goto 888 !next combination                 print*,'  -----> # fake couples (Y)',ncpy_ok
2840              if(nplused.lt.nplyz_min)goto 888 !next combination                 do icp=1,ncp_tot
2841              if(ncp_ok.lt.ncpok)goto 888 !next combination                    print*,'cp ',icp,' >'
2842         $                 ,' x ',cpcloud_xz(ixz,icp)
2843  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'       $                 ,' y ',cpcloud_yz(iyz,icp)
2844  c$$$  print*,'Configurazione cluster XZ'       $                 ,' ==> ',cpintersec(icp)
2845  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))                 enddo
2846  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  
2847                                                    
2848              if(DEBUG.EQ.1)then              if(DEBUG.EQ.1)then
                print*,'track candidate', ntracks+1  
2849                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2850                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2851                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2801  c$$$            if(AL_INI(5).gt.defmax)g Line 2878  c$$$            if(AL_INI(5).gt.defmax)g
2878                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2879                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2880                                                                
2881                                  if(DEBUG.eq.1)
2882         $                             print*,'combination: '
2883         $                             ,cp_match(6,icp1)
2884         $                             ,cp_match(5,icp2)
2885         $                             ,cp_match(4,icp3)
2886         $                             ,cp_match(3,icp4)
2887         $                             ,cp_match(2,icp5)
2888         $                             ,cp_match(1,icp6)
2889    
2890    
2891  *                             ---------------------------------------  *                             ---------------------------------------
2892  *                             check if this group of couples has been  *                             check if this group of couples has been
2893  *                             already fitted      *                             already fitted    
# Line 2849  c     $                                 Line 2936  c     $                                
2936       $                                   PFAdef,PFAdef,0.,0.,0.,0.)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2937  *                                   *************************  *                                   *************************
2938  *                                   -----------------------------  *                                   -----------------------------
2939                                      xgood(nplanes-ip+1)=1.                                      if(icx.gt.0.and.icy.gt.0)then
2940                                      ygood(nplanes-ip+1)=1.                                         xgood(nplanes-ip+1)=1.
2941                                      xm(nplanes-ip+1)=xPAM                                         ygood(nplanes-ip+1)=1.
2942                                      ym(nplanes-ip+1)=yPAM                                         xm(nplanes-ip+1)=xPAM
2943                                      zm(nplanes-ip+1)=zPAM                                         ym(nplanes-ip+1)=yPAM
2944                                      resx(nplanes-ip+1)=resxPAM                                         zm(nplanes-ip+1)=zPAM
2945                                      resy(nplanes-ip+1)=resyPAM                                         resx(nplanes-ip+1)=resxPAM                                      
2946                                           resy(nplanes-ip+1)=resyPAM
2947                                           if(DEBUG.EQ.1)print*,'(X,Y)'
2948         $                                      ,nplanes-ip+1,xPAM,yPAM
2949                                        else
2950                                           xm_A(nplanes-ip+1) = xPAM_A
2951                                           ym_A(nplanes-ip+1) = yPAM_A
2952                                           xm_B(nplanes-ip+1) = xPAM_B
2953                                           ym_B(nplanes-ip+1) = yPAM_B
2954                                           zm(nplanes-ip+1)
2955         $                                      = (zPAM_A+zPAM_B)/2.
2956                                           resx(nplanes-ip+1) = resxPAM                                      
2957                                           resy(nplanes-ip+1) = resyPAM
2958                                           if(icx.eq.0.and.icy.gt.0)then
2959                                              xgood(nplanes-ip+1)=0.
2960                                              ygood(nplanes-ip+1)=1.
2961                                              resx(nplanes-ip+1) = 1000.
2962                                              if(DEBUG.EQ.1)print*,'(  Y)'
2963         $                                         ,nplanes-ip+1,xPAM,yPAM
2964                                           elseif(icx.gt.0.and.icy.eq.0)then
2965                                              xgood(nplanes-ip+1)=1.
2966                                              ygood(nplanes-ip+1)=0.
2967                                              if(DEBUG.EQ.1)print*,'(X  )'
2968         $                                         ,nplanes-ip+1,xPAM,yPAM
2969                                              resy(nplanes-ip+1) = 1000.
2970                                           else
2971                                              print*,'both icx=0 and icy=0'
2972         $                                         ,' ==> IMPOSSIBLE!!'
2973                                           endif
2974                                        endif
2975  *                                   -----------------------------  *                                   -----------------------------
2976                                   endif                                   endif
2977                                enddo !end loop on planes                                enddo !end loop on planes
# Line 2896  c                                 chi2=- Line 3012  c                                 chi2=-
3012  *     **********************************************************  *     **********************************************************
3013    
3014                                if(chi2.le.0.)goto 666                                              if(chi2.le.0.)goto 666              
3015                                  if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
3016                                  if(chi2.ne.chi2)goto 666  !OPTIMIZATION
3017    
3018  *     --------------------------  *     --------------------------
3019  *     STORE candidate TRACK INFO  *     STORE candidate TRACK INFO
# Line 2941  c                                    mas Line 3059  c                                    mas
3059       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
3060                                      SENSOR_STORE(nplanes-ip+1,ntracks)                                      SENSOR_STORE(nplanes-ip+1,ntracks)
3061       $                              = is_cp(cp_match(ip,hit_plane(ip)))       $                              = is_cp(cp_match(ip,hit_plane(ip)))
3062                                      LADDER_STORE(nplanes-ip+1,ntracks)                                      
3063       $                                   = LADDER(                                      icl=
3064       $                                   clx(ip,icp_cp(       $                                   clx(ip,icp_cp(
3065       $                                   cp_match(ip,hit_plane(ip)       $                                   cp_match(ip,hit_plane(ip)
3066       $                                   ))));       $                                   )));
3067                                        if(icl.eq.0)
3068         $                                   icl=
3069         $                                   cly(ip,icp_cp(
3070         $                                   cp_match(ip,hit_plane(ip)
3071         $                                   )));
3072    
3073                                        LADDER_STORE(nplanes-ip+1,ntracks)
3074         $                                   = LADDER(icl);
3075                                   else                                   else
3076                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
3077                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0
# Line 2979  c                                    mas Line 3105  c                                    mas
3105                
3106        if(ntracks.eq.0)then        if(ntracks.eq.0)then
3107           iflag=1           iflag=1
3108           return  cc         return
3109        endif        endif
3110                
 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  
3111        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
3112          print*,'****** TRACK CANDIDATES *****************'          print*,'****** TRACK CANDIDATES *****************'
3113          print*,'#         R. chi2        RIG         ndof'          print*,'#         R. chi2        RIG         ndof'
# Line 3070  c$$$         bxyz(3)=0 Line 3187  c$$$         bxyz(3)=0
3187  *     using improved PFAs  *     using improved PFAs
3188  *     -------------------------------------------------  *     -------------------------------------------------
3189  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3190           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.  c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3191    c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3192    c$$$            
3193    c$$$            id=CP_STORE(nplanes-ip+1,ibest)
3194    c$$$            
3195    c$$$            is=is_cp(id)
3196    c$$$            icp=icp_cp(id)
3197    c$$$            if(ip_cp(id).ne.ip)
3198    c$$$     $           print*,'OKKIO!!'
3199    c$$$     $           ,'id ',id,is,icp
3200    c$$$     $           ,ip_cp(id),ip
3201    c$$$            icx=clx(ip,icp)
3202    c$$$            icy=cly(ip,icp)
3203    c$$$c            call xyz_PAM(icx,icy,is,
3204    c$$$c     $           PFA,PFA,
3205    c$$$c     $           AXV_STORE(nplanes-ip+1,ibest),
3206    c$$$c     $           AYV_STORE(nplanes-ip+1,ibest))
3207    c$$$            call xyz_PAM(icx,icy,is,
3208    c$$$     $           PFA,PFA,
3209    c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
3210    c$$$     $           AYV_STORE(nplanes-ip+1,ibest),
3211    c$$$     $           bxyz(1),
3212    c$$$     $           bxyz(2)
3213    c$$$     $           )
3214    c$$$
3215    c$$$            xm(nplanes-ip+1) = xPAM
3216    c$$$            ym(nplanes-ip+1) = yPAM
3217    c$$$            zm(nplanes-ip+1) = zPAM
3218    c$$$            xgood(nplanes-ip+1) = 1
3219    c$$$            ygood(nplanes-ip+1) = 1
3220    c$$$            resx(nplanes-ip+1) = resxPAM
3221    c$$$            resy(nplanes-ip+1) = resyPAM
3222    c$$$
3223    c$$$            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3224    c$$$            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3225             if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
3226       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3227                            
3228              id=CP_STORE(nplanes-ip+1,ibest)              id=CP_STORE(nplanes-ip+1,ibest)
# Line 3095  c     $           AYV_STORE(nplanes-ip+1 Line 3247  c     $           AYV_STORE(nplanes-ip+1
3247       $           bxyz(2)       $           bxyz(2)
3248       $           )       $           )
3249    
3250              xm(nplanes-ip+1) = xPAM              if(icx.gt.0.and.icy.gt.0)then
3251              ym(nplanes-ip+1) = yPAM                 xm(nplanes-ip+1) = xPAM
3252              zm(nplanes-ip+1) = zPAM                 ym(nplanes-ip+1) = yPAM
3253              xgood(nplanes-ip+1) = 1                 zm(nplanes-ip+1) = zPAM
3254              ygood(nplanes-ip+1) = 1                 xm_A(nplanes-ip+1) = 0.
3255              resx(nplanes-ip+1) = resxPAM                 ym_A(nplanes-ip+1) = 0.
3256              resy(nplanes-ip+1) = resyPAM                 xm_B(nplanes-ip+1) = 0.
3257                   ym_B(nplanes-ip+1) = 0.
3258              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))                 xgood(nplanes-ip+1) = 1
3259              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))                 ygood(nplanes-ip+1) = 1
3260                   resx(nplanes-ip+1) = resxPAM
3261                   resy(nplanes-ip+1) = resyPAM
3262                   dedxtrk_x(nplanes-ip+1)=
3263         $              sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3264                   dedxtrk_y(nplanes-ip+1)=
3265         $              sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3266                else
3267                   xm(nplanes-ip+1) = 0.
3268                   ym(nplanes-ip+1) = 0.
3269                   zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
3270                   xm_A(nplanes-ip+1) = xPAM_A
3271                   ym_A(nplanes-ip+1) = yPAM_A
3272                   xm_B(nplanes-ip+1) = xPAM_B
3273                   ym_B(nplanes-ip+1) = yPAM_B
3274                   xgood(nplanes-ip+1) = 0
3275                   ygood(nplanes-ip+1) = 0
3276                   resx(nplanes-ip+1) = 1000.!resxPAM
3277                   resy(nplanes-ip+1) = 1000.!resyPAM
3278                   dedxtrk_x(nplanes-ip+1)= 0
3279                   dedxtrk_y(nplanes-ip+1)= 0
3280                   if(icx.gt.0)then
3281                      xgood(nplanes-ip+1) = 1
3282                      resx(nplanes-ip+1) = resxPAM
3283                      dedxtrk_x(nplanes-ip+1)=
3284         $                 sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3285                   elseif(icy.gt.0)then
3286                      ygood(nplanes-ip+1) = 1
3287                      resy(nplanes-ip+1) = resyPAM
3288                      dedxtrk_y(nplanes-ip+1)=
3289         $                 sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3290                   endif
3291                endif
3292                            
3293  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3294  *     -------------------------------------------------  *     -------------------------------------------------
# Line 3116  c     $           AYV_STORE(nplanes-ip+1 Line 3300  c     $           AYV_STORE(nplanes-ip+1
3300                                
3301              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
3302              ygood(nplanes-ip+1)=0              ygood(nplanes-ip+1)=0
3303    
3304                CP_STORE(nplanes-ip+1,ibest)=0 !re-init
3305                CLS_STORE(nplanes-ip+1,ibest)=0
3306    
3307                                
3308  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3309  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
# Line 3138  c     $           AYV_STORE(nplanes-ip+1 Line 3326  c     $           AYV_STORE(nplanes-ip+1
3326  *     ===========================================  *     ===========================================
3327  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3328  *     ===========================================  *     ===========================================
3329  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'              distmin=100000000.
             distmin=1000000.  
3330              xmm = 0.              xmm = 0.
3331              ymm = 0.              ymm = 0.
3332              zmm = 0.              zmm = 0.
# Line 3152  c            if(DEBUG.EQ.1)print*,'>>>> Line 3339  c            if(DEBUG.EQ.1)print*,'>>>>
3339              do icp=1,ncp_plane(ip) !loop on couples on plane icp              do icp=1,ncp_plane(ip) !loop on couples on plane icp
3340                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3341                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3342                   if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
3343                 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
3344  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
3345  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
# Line 3170  c     $              cl_used(icy).eq.1.o Line 3358  c     $              cl_used(icy).eq.1.o
3358                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3359  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3360                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3361                 if(DEBUG.EQ.1)print*,'( couple ',id                 if(DEBUG.EQ.1)
3362         $              print*,'( couple ',id
3363       $              ,' ) distance ',distance       $              ,' ) distance ',distance
3364                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3365                    xmm = xPAM                    xmm = xPAM
# Line 3182  c               distance = distance / RC Line 3371  c               distance = distance / RC
3371                    idm = id                                      idm = id                  
3372                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3373                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3374  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!                    clincnewc=10*sqrt(rymm**2+rxmm**2
3375                    clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI       $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))
      $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI  
3376                 endif                 endif
3377   1188          continue   1188          continue
3378              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3379  c            if(distmin.le.clinc)then     !QUIQUI                            if(distmin.le.clincnewc)then    
             if(distmin.le.clincnewc)then     !QUIQUI                
3380  *              -----------------------------------  *              -----------------------------------
3381                 xm(nplanes-ip+1) = xmm !<<<                 xm(nplanes-ip+1) = xmm !<<<
3382                 ym(nplanes-ip+1) = ymm !<<<                 ym(nplanes-ip+1) = ymm !<<<
# Line 3203  c            if(distmin.le.clinc)then   Line 3390  c            if(distmin.le.clinc)then  
3390  *              -----------------------------------  *              -----------------------------------
3391                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3392                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3393       $              ,' (dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clincnewc,' )'
3394                 goto 133         !next plane                 goto 133         !next plane
3395              endif              endif
3396  *     ================================================  *     ================================================
3397  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3398  *                     either from a couple or single  *                     either from a couple or single
3399  *     ================================================  *     ================================================
 c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'  
3400              distmin=1000000.              distmin=1000000.
3401              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3402              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3229  c            if(DEBUG.EQ.1)print*,'>>>> Line 3415  c            if(DEBUG.EQ.1)print*,'>>>>
3415              do icp=1,ncp_plane(ip) !loop on cluster inside couples              do icp=1,ncp_plane(ip) !loop on cluster inside couples
3416                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3417                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3418                   if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
3419                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3420                 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
3421  *                                                !jump to the next couple  *                                                !jump to the next couple
# Line 3247  c     $              AXV_STORE(nplanes-i Line 3434  c     $              AXV_STORE(nplanes-i
3434       $              )                     $              )              
3435                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3436  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3437                 if(DEBUG.EQ.1)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)
3438         $              print*,'( cl-X ',icx
3439       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3440                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3441                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3280  c     $              0.,AYV_STORE(nplane Line 3468  c     $              0.,AYV_STORE(nplane
3468       $              )       $              )
3469                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3470  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3471                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)
3472         $              print*,'( cl-Y ',icy
3473       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3474                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3475                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3300  c                 dedxmm = sgnl(icy)  !( Line 3489  c                 dedxmm = sgnl(icy)  !(
3489  11882          continue  11882          continue
3490              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3491  *----- single clusters -----------------------------------------------    *----- single clusters -----------------------------------------------  
 c            print*,'## ncls(',ip,') ',ncls(ip)  
3492              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
 c               print*,'-',ic,'-'  
3493                 icl=cls(ip,ic)                 icl=cls(ip,ic)
 c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  
3494                 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)
3495       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3496       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
# Line 3326  c               if(cl_used(icl).eq.1.or. Line 3512  c               if(cl_used(icl).eq.1.or.
3512    
3513                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3514  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3515                 if(DEBUG.EQ.1)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)
3516         $              print*,'( cl-s ',icl
3517       $              ,' ) distance ',distance       $              ,' ) distance ',distance
3518                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
 c                  if(DEBUG.EQ.1)print*,'YES'  
3519                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3520                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3521                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3350  c                  if(DEBUG.EQ.1)print*, Line 3536  c                  if(DEBUG.EQ.1)print*,
3536                 endif                                   endif                  
3537  18882          continue  18882          continue
3538              enddo               !end loop on single clusters              enddo               !end loop on single clusters
 c            print*,'## distmin ', distmin,' clinc ',clinc  
3539    
 c     QUIQUI------------  
 c     anche qui: non ci vuole clinc???  
3540              if(iclm.ne.0)then              if(iclm.ne.0)then
3541                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3542                    clincnew=                    clincnew=
# Line 3364  c     anche qui: non ci vuole clinc??? Line 3547  c     anche qui: non ci vuole clinc???
3547       $                 10*       $                 10*
3548       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3549                 endif                 endif
3550  c     QUIQUI------------  
3551                                 if(distmin.le.clincnew)then  
                if(distmin.le.clincnew)then   !QUIQUI  
 c     if(distmin.le.clinc)then          !QUIQUI            
3552                                        
3553                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                        CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3554  *     ----------------------------  *     ----------------------------
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3555                    if(mod(VIEW(iclm),2).eq.0)then                    if(mod(VIEW(iclm),2).eq.0)then
3556                       XGOOD(nplanes-ip+1)=1.                       XGOOD(nplanes-ip+1)=1.
3557                       resx(nplanes-ip+1)=rxmm                       resx(nplanes-ip+1)=rxmm
3558                       if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm                       if(DEBUG.EQ.1)
3559         $                    print*,'%%%% included X-cl ',iclm
3560       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3561       $                    ,', dist.= ',distmin       $                    ,', dist.= ',distmin
3562       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clincnew,' )'
3563                    else                    else
3564                       YGOOD(nplanes-ip+1)=1.                       YGOOD(nplanes-ip+1)=1.
3565                       resy(nplanes-ip+1)=rymm                       resy(nplanes-ip+1)=rymm
3566                       if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm                       if(DEBUG.EQ.1)
3567         $                    print*,'%%%% included Y-cl ',iclm
3568       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3569       $                    ,', dist.= ', distmin       $                    ,', dist.= ', distmin
3570       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clincnew,' )'
3571                    endif                    endif
 c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3572  *     ----------------------------  *     ----------------------------
3573                    xm_A(nplanes-ip+1) = xmm_A                    xm_A(nplanes-ip+1) = xmm_A
3574                    ym_A(nplanes-ip+1) = ymm_A                    ym_A(nplanes-ip+1) = ymm_A
# Line 3408  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3589  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3589        return        return
3590        end        end
3591    
3592    
3593  ***************************************************  ***************************************************
3594  *                                                 *  *                                                 *
3595  *                                                 *  *                                                 *
# Line 3417  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3599  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3599  *                                                 *  *                                                 *
3600  **************************************************  **************************************************
3601  *  *
       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  
   
   
   
3602    
3603    
3604    
# Line 3688  c$$$               cl_used(icl)=ntrk   ! Line 3798  c$$$               cl_used(icl)=ntrk   !
3798        integer ssensor,sladder        integer ssensor,sladder
3799        pig=acos(-1.)        pig=acos(-1.)
3800    
3801    
3802    
3803  *     -------------------------------------  *     -------------------------------------
3804        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3805        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
# Line 3735  c$$$               cl_used(icl)=ntrk   ! Line 3847  c$$$               cl_used(icl)=ntrk   !
3847           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3848           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3849        
3850    
3851    ccc         print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3852    
3853           ax   = axv_nt(ip,ntr)           ax   = axv_nt(ip,ntr)
3854           ay   = ayv_nt(ip,ntr)           ay   = ayv_nt(ip,ntr)
3855           bfx  = BX_STORE(ip,IDCAND)           bfx  = BX_STORE(ip,IDCAND)
# Line 3750  c$$$         angy    = 180.*atan(tgtemp) Line 3865  c$$$         angy    = 180.*atan(tgtemp)
3865           angy = effectiveangle(ay,2*ip-1,bfx)           angy = effectiveangle(ay,2*ip-1,bfx)
3866                    
3867                    
 c         print*,'* ',ip,bfx,bfy,angx,angy  
3868    
3869           id  = CP_STORE(ip,IDCAND) ! couple id           id  = CP_STORE(ip,IDCAND) ! couple id
3870           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
# Line 3766  c           >>> is a couple Line 3880  c           >>> is a couple
3880              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3881              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3882    
3883              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            
3884    
3885              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)))  
3886    
3887                   xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))              
3888    
3889              if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)                 if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3890       $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)       $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3891              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)                
3892       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)                 multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3893         $              +10000*mult(cltrx(ip,ntr))
3894                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3895         $              /clsigma(indmax(cltrx(ip,ntr)))
3896                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3897                   xpu(ip,ntr)      = corr
3898    
3899              multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))              endif
3900       $                         +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  
3901    
3902              multmaxy(ip,ntr) = maxs(cltry(ip,ntr))                 cl_used(cltry(ip,ntr)) = 1 !tag used clusters          
      $                         +10000*mult(cltry(ip,ntr))  
             seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))  
      $           /clsigma(indmax(cltry(ip,ntr)))  
             call applypfa(PFA,cltry(ip,ntr),angy,corr,res)  
             ypu(ip,ntr)      = corr  
3903    
3904           elseif(icl.ne.0)then                 ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3905    
3906                   if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3907         $              dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3908                  
3909                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3910         $              +10000*mult(cltry(ip,ntr))
3911                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3912         $              /clsigma(indmax(cltry(ip,ntr)))
3913                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3914                   ypu(ip,ntr)      = corr
3915                endif
3916    
3917             elseif(icl.ne.0)then
3918    
3919              cl_used(icl) = 1    !tag used clusters                        cl_used(icl) = 1    !tag used clusters          
3920    
# Line 3835  c           >>> is a couple Line 3956  c           >>> is a couple
3956           do ip=1,6           do ip=1,6
3957              print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)              print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3958           enddo           enddo
3959             print*,'dedx: '
3960             do ip=1,6
3961                print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3962             enddo
3963        endif        endif
3964    
 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*,'-----------------------'  
   
3965        end        end
3966    
3967        subroutine fill_level2_siglets        subroutine fill_level2_siglets
# Line 3890  c         if( mask_view(iv).ne.0 )good2( Line 4008  c         if( mask_view(iv).ne.0 )good2(
4008                 sxbad(nclsx)  = nbadstrips(1,icl)                 sxbad(nclsx)  = nbadstrips(1,icl)
4009                 multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)                 multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
4010                                
 cc               print*,icl,' >>>> ',sxbad(nclsx)  
4011    
4012                 do is=1,2                 do is=1,2
4013  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
# Line 3898  c                  call xyz_PAM(icl,0,is Line 4015  c                  call xyz_PAM(icl,0,is
4015                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
4016                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
4017                 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)  
4018              else                          !=== Y views              else                          !=== Y views
4019                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4020                 planey(nclsy) = ip                 planey(nclsy) = ip
# Line 3912  c$$$               print*,'xs(2,nclsx)   Line 4024  c$$$               print*,'xs(2,nclsx)  
4024                 sybad(nclsy)  = nbadstrips(1,icl)                 sybad(nclsy)  = nbadstrips(1,icl)
4025                 multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)                 multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
4026    
 cc               print*,icl,' >>>> ',sybad(nclsy)  
4027    
4028                 do is=1,2                 do is=1,2
4029  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
# Line 3920  c                  call xyz_PAM(0,icl,is Line 4031  c                  call xyz_PAM(0,icl,is
4031                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
4032                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
4033                 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)  
4034              endif              endif
4035           endif           endif
4036    
# Line 3999  c$$$               print*,'ys(2,nclsy)   Line 4105  c$$$               print*,'ys(2,nclsy)  
4105                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4106                 nnn=nnn+ptcloud_yz(iyz)                 nnn=nnn+ptcloud_yz(iyz)
4107              enddo              enddo
4108              do ipt=1,nnn              do ipt=1,min(ndblt_max_nt,nnn)
4109                 db_cloud_nt(ipt)=db_cloud(ipt)                 db_cloud_nt(ipt)=db_cloud(ipt)
4110               enddo               enddo
4111           endif           endif
# Line 4012  c$$$               print*,'ys(2,nclsy)   Line 4118  c$$$               print*,'ys(2,nclsy)  
4118                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4119                 nnn=nnn+ptcloud_xz(ixz)                               nnn=nnn+ptcloud_xz(ixz)              
4120              enddo              enddo
4121              do ipt=1,nnn              do ipt=1,min(ntrpt_max_nt,nnn)
4122                tr_cloud_nt(ipt)=tr_cloud(ipt)                tr_cloud_nt(ipt)=tr_cloud(ipt)
4123               enddo               enddo
4124           endif           endif
Legend:
Removed from v.1.34  
changed lines
  Added in v.1.39
  ViewVC Help
Powered by ViewVC 1.1.23