/[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.2 by pam-fi, Tue May 30 16:30:37 2006 UTC revision 1.10 by pam-fi, Thu Nov 2 11:19:47 2006 UTC
# Line 12  Line 12 
12        subroutine track_finding(iflag)        subroutine track_finding(iflag)
13    
14        include 'commontracker.f'        include 'commontracker.f'
15          include 'level1.f'
16        include 'common_momanhough.f'        include 'common_momanhough.f'
17        include 'common_mech.f'        include 'common_mech.f'
18        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
19        include 'common_mini_2.f'        include 'common_mini_2.f'
20        include 'calib.f'        include 'calib.f'
21        include 'level1.f'  c      include 'level1.f'
22        include 'level2.f'        include 'level2.f'
23    
24        include 'momanhough_init.f'  c      include 'momanhough_init.f'
25                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
26  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
27  *     STEP 1  *     STEP 1
28  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 47  c      common/dbg/DEBUG Line 45  c      common/dbg/DEBUG
45  c      iflag=0  c      iflag=0
46        call cl_to_couples(iflag)        call cl_to_couples(iflag)
47        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
48           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
49        endif        endif
50                
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     selezione di tracce pulite per diagnostica  
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 c$$$         if(DEBUG)then  
 c$$$            do ip=1,nplanes  
 c$$$               if(ncp_plane(ip).ne.1)good2=.false.  
 c$$$            enddo  
 c$$$c            if(good2.eq.0)goto 100!next event  
 c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event  
 c$$$         endif  
           
   
   
51  *-----------------------------------------------------  *-----------------------------------------------------
52  *-----------------------------------------------------  *-----------------------------------------------------
53  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 94  c$$$         endif Line 79  c$$$         endif
79  c      iflag=0  c      iflag=0
80        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
81        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
82           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
83        endif        endif
84                
85                
# Line 123  c      iflag=0 Line 108  c      iflag=0
108  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
111    *     count number of hit planes
112          planehit=0                
113          do np=1,nplanes          
114            if(ncp_plane(np).ne.0)then  
115              planehit=planehit+1  
116            endif                  
117          enddo                    
118          if(planehit.lt.3) goto 880 ! exit              
119          
120          nptxz_min=x_min_start              
121          nplxz_min=x_min_start              
122                
123          nptyz_min=y_min_start              
124          nplyz_min=y_min_start              
125    
126  c      iflag=0        cutdistyz=cutystart      
127          cutdistxz=cutxstart      
128    
129     878  continue
130        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
131        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
132           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
133        endif        endif
134  c      iflag=0        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
135            if(cutdistyz.lt.maxcuty/2)then
136              cutdistyz=cutdistyz+cutystep
137            else
138              cutdistyz=cutdistyz+(3*cutystep)
139            endif
140            goto 878                
141          endif                    
142    
143          if(planehit.eq.3) goto 881          
144          
145     879  continue  
146        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
147        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
148           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
149        endif        endif
150                                    
151          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
152            cutdistxz=cutdistxz+cutxstep
153            goto 879                
154          endif                    
155    
156        
157     881  continue  
158    *     if there is at least three planes on the Y view decreases cuts on X view
159          if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
160         $     nplxz_min.ne.y_min_start)then
161            nptxz_min=x_min_step    
162            nplxz_min=x_min_start-x_min_step              
163            goto 879                
164          endif                    
165            
166   880  return   880  return
167        end        end
168    
# Line 144  c      iflag=0 Line 172  c      iflag=0
172        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
173    
174        include 'commontracker.f'        include 'commontracker.f'
175          include 'level1.f'
176        include 'common_momanhough.f'        include 'common_momanhough.f'
177        include 'common_mech.f'        include 'common_mech.f'
178        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
179        include 'common_mini_2.f'        include 'common_mini_2.f'
180        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
181        include 'level2.f'        include 'level2.f'
182    
183        include 'momanhough_init.f'  c      include 'momanhough_init.f'
184                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
185        logical FIMAGE            !        logical FIMAGE            !
186    
187  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 198  c         iflag=0 Line 223  c         iflag=0
223           ibest=0                !best track among candidates           ibest=0                !best track among candidates
224           iimage=0               !track image           iimage=0               !track image
225  *     ------------- select the best track -------------  *     ------------- select the best track -------------
226    c$$$         rchi2best=1000000000.
227    c$$$         do i=1,ntracks
228    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
229    c$$$     $         RCHI2_STORE(i).gt.0)then
230    c$$$               ibest=i
231    c$$$               rchi2best=RCHI2_STORE(i)
232    c$$$            endif
233    c$$$         enddo
234    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
235    
236           rchi2best=1000000000.           rchi2best=1000000000.
237             ndofbest=0             !(1)
238           do i=1,ntracks           do i=1,ntracks
239              if(RCHI2_STORE(i).lt.rchi2best.and.             if(RCHI2_STORE(i).lt.rchi2best.and.
240       $         RCHI2_STORE(i).gt.0)then       $          RCHI2_STORE(i).gt.0)then
241                 ndof=0             !(1)
242                 do ii=1,nplanes    !(1)
243                   ndof=ndof        !(1)
244         $              +int(xgood_store(ii,i)) !(1)
245         $              +int(ygood_store(ii,i)) !(1)
246                 enddo              !(1)
247                 if(ndof.ge.ndofbest)then !(1)
248                 ibest=i                 ibest=i
249                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
250              endif                 ndofbest=ndof    !(1)
251                 endif              !(1)
252               endif
253           enddo           enddo
254           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
255  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
# Line 237  c         iflag=0 Line 282  c         iflag=0
282  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
283  *     **********************************************************  *     **********************************************************
284           do i=1,5           do i=1,5
285              AL(i)=dble(AL_STORE(i,icand))              AL(i)=dble(AL_STORE(i,icand))            
286           enddo           enddo
287             IDCAND = icand         !fitted track-candidate
288           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
289           jstep=0                !# minimization steps           jstep=0                !# minimization steps
290    
291           call mini_2(jstep,ifail)           iprint=0
292             if(DEBUG)iprint=1
293             call mini2(jstep,ifail,iprint)
294           if(ifail.ne.0) then           if(ifail.ne.0) then
295              if(DEBUG)then              if(DEBUG)then
296                 print *,                 print *,
297       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (mini2) '
298       $              ,iev       $              ,iev
299              endif              endif
300              chi2=-chi2              chi2=-chi2
# Line 310  c         print*,'++++++++++ iimage,fima Line 358  c         print*,'++++++++++ iimage,fima
358  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
359    
360           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
361              if(DEBUG)              if(verbose)
362       $           print*,       $           print*,
363       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
364       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 606  c                (implemented variable r Line 654  c                (implemented variable r
654  c*****************************************************  c*****************************************************
655                
656        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
657        include 'level1.f'        include 'level1.f'
658          include 'calib.f'
659    c      include 'level1.f'
660        include 'common_align.f'        include 'common_align.f'
661        include 'common_mech.f'        include 'common_mech.f'
662        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
# Line 666  c      double precision xi_B,yi_B,zi_B Line 715  c      double precision xi_B,yi_B,zi_B
715              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
716           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
717  c            cog2 = cog(2,icx)  c            cog2 = cog(2,icx)
718  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)              c            etacorr = pfaeta2(cog2,viewx,nldx,angx)            
719  c            stripx = stripx + etacorr  c            stripx = stripx + etacorr
720              stripx = stripx + pfa_eta2(icx,angx)            !(3)              stripx = stripx + pfaeta2(icx,angx)            !(3)
721              resxPAM = risx_eta2(angx)                       !   (4)              resxPAM = risx_eta2(angx)                       !   (4)
722              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
723       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
724              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
725           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                         !(3)
726              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)            !(3)
727              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(angx)                       !   (4)
728              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
729       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
730              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
731           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                         !(3)
732              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            !(3)
733              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(angx)                       !   (4)
734              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
735       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
736              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
737           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then                          !(3)
738              stripx = stripx + pfa_eta(icx,angx)             !(3)              stripx = stripx + pfaeta(icx,angx)             !(3)
739              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = ris_eta(icx,angx)                     !   (4)
740              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
741       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
# Line 701  c            resxPAM = resxPAM*fbad_cog( Line 750  c            resxPAM = resxPAM*fbad_cog(
750           endif           endif
751    
752        endif        endif
753    c      if(icy.eq.0.and.icx.ne.0)
754    c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'
755                
756  *     -----------------  *     -----------------
757  *     CLUSTER Y  *     CLUSTER Y
# Line 728  c            resxPAM = resxPAM*fbad_cog( Line 779  c            resxPAM = resxPAM*fbad_cog(
779              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
780           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
781  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
782  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
783  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
784              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
785              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(angy)                       !   (4)
786              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
787              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
788       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
789           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
790              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
791              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
792              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
793       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
794           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
795              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
796              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
797              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
798       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
799           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
800              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
801              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
802  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
803              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 772  C======================================= Line 823  C=======================================
823  c------------------------------------------------------------------------  c------------------------------------------------------------------------
824  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
825  c------------------------------------------------------------------------  c------------------------------------------------------------------------
826             if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
827         $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
828                print*,'xyz_PAM (couple):',
829         $          ' WARNING: false X strip: strip ',stripx
830             endif
831           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
832           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
833           zi = 0.           zi = 0.
# Line 858  C======================================= Line 914  C=======================================
914              nldy = nldx              nldy = nldx
915              viewy = viewx - 1              viewy = viewx - 1
916    
917    c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
918    c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
919                if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
920         $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
921                   print*,'xyz_PAM (X-singlet):',
922         $             ' WARNING: false X strip: strip ',stripx
923                endif
924              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
925    
926              xi_A = xi              xi_A = xi
# Line 1136  c$$$         print*,' resolution ',resxP Line 1199  c$$$         print*,' resolution ',resxP
1199  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1200  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1201  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1202                   if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1203         $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1204    c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1205                      print*,'whichsensor: ',
1206         $                ' WARNING: false X strip: strip ',stripx
1207                   endif
1208                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1209                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
1210                 zi = 0.                 zi = 0.
# Line 1263  c     $              ,iv,xvv(iv),yvv(iv) Line 1332  c     $              ,iv,xvv(iv),yvv(iv)
1332  *     it returns the plane number  *     it returns the plane number
1333  *      *    
1334        include 'commontracker.f'        include 'commontracker.f'
1335          include 'level1.f'
1336  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1337        include 'common_momanhough.f'        include 'common_momanhough.f'
1338                
# Line 1300  c      include 'common_analysis.f' Line 1370  c      include 'common_analysis.f'
1370  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1371  *      *    
1372        include 'commontracker.f'        include 'commontracker.f'
1373          include 'level1.f'
1374  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1375        include 'common_momanhough.f'        include 'common_momanhough.f'
1376                
# Line 1328  c      include 'common_analysis.f' Line 1399  c      include 'common_analysis.f'
1399  *     positive if sensor =2  *     positive if sensor =2
1400  *  *
1401        include 'commontracker.f'        include 'commontracker.f'
1402          include 'level1.f'
1403  c      include 'calib.f'  c      include 'calib.f'
1404  c      include 'level1.f'  c      include 'level1.f'
1405  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1639  c$$$      end Line 1711  c$$$      end
1711        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1712    
1713        include 'commontracker.f'        include 'commontracker.f'
1714          include 'level1.f'
1715        include 'common_momanhough.f'        include 'common_momanhough.f'
1716        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1717        include 'calib.f'        include 'calib.f'
1718        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1719    
1720  *     output flag  *     output flag
1721  *     --------------  *     --------------
# Line 1654  c      common/dbg/DEBUG Line 1724  c      common/dbg/DEBUG
1724  *     --------------  *     --------------
1725        integer iflag        integer iflag
1726    
1727        integer badseed,badcl        integer badseed,badclx,badcly
1728    
1729  *     init variables  *     init variables
1730        ncp_tot=0        ncp_tot=0
# Line 1670  c      common/dbg/DEBUG Line 1740  c      common/dbg/DEBUG
1740           ncls(ip)=0           ncls(ip)=0
1741        enddo        enddo
1742        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1743           cl_single(icl)=1           cl_single(icl) = 1
1744           cl_good(icl)=0           cl_good(icl)   = 0
1745          enddo
1746          do iv=1,nviews
1747             ncl_view(iv)  = 0
1748             mask_view(iv) = 0      !all included
1749        enddo        enddo
1750                
1751    *     count number of cluster per view
1752          do icl=1,nclstr1
1753             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1754          enddo
1755    *     mask views with too many clusters
1756          do iv=1,nviews
1757             if( ncl_view(iv).gt. nclustermax)then
1758                mask_view(iv) = 1
1759                if(VERBOSE)print*,' * WARNING * cl_to_couple: n.clusters > '
1760         $           ,nclustermax,' on view ', iv,' --> masked!'
1761             endif
1762          enddo
1763    
1764    
1765  *     start association  *     start association
1766        ncouples=0        ncouples=0
1767        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1768           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1769                    
1770  *     ----------------------------------------------------  *     ----------------------------------------------------
1771    *     jump masked views (X VIEW)
1772    *     ----------------------------------------------------
1773             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1774    *     ----------------------------------------------------
1775  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1776    *     ----------------------------------------------------
1777           if(dedx(icx).lt.dedx_x_min)then           if(dedx(icx).lt.dedx_x_min)then
1778              cl_single(icx)=0              cl_single(icx)=0
1779              goto 10              goto 10
1780           endif           endif
1781    *     ----------------------------------------------------
1782  *     cut BAD (X VIEW)              *     cut BAD (X VIEW)            
1783    *     ----------------------------------------------------
1784           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1785           ifirst=INDSTART(icx)           ifirst=INDSTART(icx)
1786           if(icx.ne.nclstr1) then           if(icx.ne.nclstr1) then
# Line 1693  c      common/dbg/DEBUG Line 1788  c      common/dbg/DEBUG
1788           else           else
1789              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1790           endif           endif
1791           badcl=badseed           badclx=badseed
1792           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1793              ibad=1              ibad=1
1794              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1703  c      common/dbg/DEBUG Line 1798  c      common/dbg/DEBUG
1798       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1799       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1800              endif              endif
1801              badcl=badcl*ibad              badclx=badclx*ibad
1802           enddo           enddo
1803    *     ----------------------------------------------------
1804    *     >>> eliminato il taglio sulle BAD <<<
1805    *     ----------------------------------------------------
1806  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1807  c     cl_single(icx)=0  c     cl_single(icx)=0
1808  c     goto 10  c     goto 10
# Line 1719  c     endif Line 1817  c     endif
1817              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1818                            
1819  *     ----------------------------------------------------  *     ----------------------------------------------------
1820    *     jump masked views (Y VIEW)
1821    *     ----------------------------------------------------
1822                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1823    
1824    *     ----------------------------------------------------
1825  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1826    *     ----------------------------------------------------
1827              if(dedx(icy).lt.dedx_y_min)then              if(dedx(icy).lt.dedx_y_min)then
1828                 cl_single(icy)=0                 cl_single(icy)=0
1829                 goto 20                 goto 20
1830              endif              endif
1831    *     ----------------------------------------------------
1832  *     cut BAD (Y VIEW)              *     cut BAD (Y VIEW)            
1833    *     ----------------------------------------------------
1834              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1835              ifirst=INDSTART(icy)              ifirst=INDSTART(icy)
1836              if(icy.ne.nclstr1) then              if(icy.ne.nclstr1) then
# Line 1732  c     endif Line 1838  c     endif
1838              else              else
1839                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1840              endif              endif
1841              badcl=badseed              badcly=badseed
1842              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1843                 ibad=1                 ibad=1
1844                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1741  c     endif Line 1847  c     endif
1847       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1848       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1849       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1850                 badcl=badcl*ibad                 badcly=badcly*ibad
1851              enddo              enddo
1852    *     ----------------------------------------------------
1853    *     >>> eliminato il taglio sulle BAD <<<
1854    *     ----------------------------------------------------
1855  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1856  c     cl_single(icy)=0  c     cl_single(icy)=0
1857  c     goto 20  c     goto 20
1858  c     endif  c     endif
1859  *     ----------------------------------------------------  *     ----------------------------------------------------
1860                            
               
1861              cl_good(icy)=1                                cl_good(icy)=1                  
1862              nply=npl(VIEW(icy))              nply=npl(VIEW(icy))
1863              nldy=nld(MAXS(icy),VIEW(icy))              nldy=nld(MAXS(icy),VIEW(icy))
# Line 1760  c     endif Line 1868  c     endif
1868  *     geometrical consistency (same plane and ladder)  *     geometrical consistency (same plane and ladder)
1869              if(nply.eq.nplx.and.nldy.eq.nldx)then              if(nply.eq.nplx.and.nldy.eq.nldx)then
1870  *     charge correlation  *     charge correlation
1871                 ddd=(dedx(icy)  *     (modified to be applied only below saturation... obviously)
1872       $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
1873                 ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1874                 cut=chcut*sch(nplx,nldx)       $              .and.
1875                 if(abs(ddd).gt.cut)goto 20 !charge not consistent       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1876                       $              .and.
1877         $              (badclx.eq.1.and.badcly.eq.1)
1878         $              .and.
1879         $              .true.)then
1880    
1881                      ddd=(dedx(icy)
1882         $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1883                      ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1884    
1885    c                  cut = chcut * sch(nplx,nldx)
1886    
1887                      sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1888         $                 -kch(nplx,nldx)*cch(nplx,nldx))
1889                      sss=sss/sqrt(kch(nplx,nldx)**2+1)
1890                      cut = chcut * (16 + sss/50.)
1891    
1892                      if(abs(ddd).gt.cut)then
1893                         goto 20    !charge not consistent
1894                      endif
1895                   endif
1896                                
1897  *     ------------------> COUPLE <------------------  *     ------------------> COUPLE <------------------
1898  *     check to do not overflow vector dimentions  *     check to do not overflow vector dimentions
1899                 if(ncp_plane(nplx).gt.ncouplemax)then  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
1900                    if(DEBUG)print*,  c$$$                  if(DEBUG)print*,
1901       $                    ' ** warning ** number of identified'//  c$$$     $                    ' ** warning ** number of identified'//
1902       $                    ' couples on plane ',nplx,  c$$$     $                    ' couples on plane ',nplx,
1903       $                    ' exceeds vector dimention'//  c$$$     $                    ' exceeds vector dimention'//
1904       $                    ' ( ',ncouplemax,' )'  c$$$     $                    ' ( ',ncouplemax,' )'
1905  c     good2=.false.  c$$$c     good2=.false.
1906  c     goto 880   !fill ntp and go to next event  c$$$c     goto 880   !fill ntp and go to next event
1907                    iflag=1  c$$$                  iflag=1
1908                    return  c$$$                  return
1909                 endif  c$$$               endif
1910                                
1911                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).eq.ncouplemax)then
1912                    if(DEBUG)print*,                    if(verbose)print*,
1913       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1914       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1915       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1916       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' )'
1917  c     good2=.false.  c     good2=.false.
1918  c     goto 880   !fill ntp and go to next event                      c     goto 880   !fill ntp and go to next event
1919                    iflag=1                    mask_view(nviewx(nplx)) = 2
1920                    return                    mask_view(nviewy(nply)) = 2
1921    c                  iflag=1
1922    c                  return
1923                 endif                 endif
1924                                
1925                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
# Line 1825  c     goto 880   !fill ntp and go to nex Line 1954  c     goto 880   !fill ntp and go to nex
1954        endif        endif
1955                
1956        do ip=1,6        do ip=1,6
1957           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1958        enddo        enddo
1959  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
1960                
1961        if(ncp_tot.gt.ncp_max)then  c$$$      if(ncp_tot.gt.ncp_max)then
1962           if(DEBUG)print*,  c$$$         if(verbose)print*,
1963       $           '** warning ** number of identified '//  c$$$     $           '** warning ** number of identified '//
1964       $           'couples exceeds upper limit for Hough tr. '  c$$$     $           'couples exceeds upper limit for Hough tr. '
1965       $           ,'( ',ncp_max,' )'              c$$$     $           ,'( ',ncp_max,' )'            
1966  c            good2=.false.  c$$$         iflag=1
1967  c     goto 880       !fill ntp and go to next event  c$$$         return
1968           iflag=1  c$$$      endif
          return  
       endif  
1969                
1970        return        return
1971        end        end
# Line 1851  c     goto 880       !fill ntp and go to Line 1978  c     goto 880       !fill ntp and go to
1978  *                                                 *  *                                                 *
1979  *                                                 *  *                                                 *
1980  **************************************************  **************************************************
1981        subroutine cl_to_couples_nocharge(iflag)  c$$$      subroutine cl_to_couples_nocharge(iflag)
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
   
 *     init variables  
       ncp_tot=0  
       do ip=1,nplanes  
          do ico=1,ncouplemax  
             clx(ip,ico)=0  
             cly(ip,ico)=0  
          enddo  
          ncp_plane(ip)=0  
          do icl=1,nclstrmax_level2  
             cls(ip,icl)=1  
          enddo  
          ncls(ip)=0  
       enddo  
       do icl=1,nclstrmax_level2  
          cl_single(icl)=1  
          cl_good(icl)=0  
       enddo  
         
 *     start association  
       ncouples=0  
       do icx=1,nclstr1          !loop on cluster (X)  
          if(mod(VIEW(icx),2).eq.1)goto 10  
           
 *     ----------------------------------------------------  
 *     cut on charge (X VIEW)  
          if(dedx(icx).lt.dedx_x_min)then  
             cl_single(icx)=0  
             goto 10  
          endif  
 *     cut BAD (X VIEW)              
          badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
          ifirst=INDSTART(icx)  
          if(icx.ne.nclstr1) then  
             ilast=INDSTART(icx+1)-1  
          else  
             ilast=TOTCLLENGTH  
          endif  
          badcl=badseed  
          do igood=-ngoodstr,ngoodstr  
             ibad=1  
             if((INDMAX(icx)+igood).gt.ifirst.and.  
      $           (INDMAX(icx)+igood).lt.ilast.and.  
      $           .true.)then  
                ibad=BAD(VIEW(icx),  
      $              nvk(MAXS(icx)+igood),  
      $              nst(MAXS(icx)+igood))  
             endif  
             badcl=badcl*ibad  
          enddo  
          if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
             cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
             goto 10             !<<<<<<<<<<<<<< BAD cut  
          endif                  !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
           
          cl_good(icx)=1  
          nplx=npl(VIEW(icx))  
          nldx=nld(MAXS(icx),VIEW(icx))  
           
          do icy=1,nclstr1       !loop on cluster (Y)  
             if(mod(VIEW(icy),2).eq.0)goto 20  
               
 *     ----------------------------------------------------  
 *     cut on charge (Y VIEW)  
             if(dedx(icy).lt.dedx_y_min)then  
                cl_single(icy)=0  
                goto 20  
             endif  
 *     cut BAD (Y VIEW)              
             badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
             ifirst=INDSTART(icy)  
             if(icy.ne.nclstr1) then  
                ilast=INDSTART(icy+1)-1  
             else  
                ilast=TOTCLLENGTH  
             endif  
             badcl=badseed  
             do igood=-ngoodstr,ngoodstr  
                ibad=1  
                if((INDMAX(icy)+igood).gt.ifirst.and.  
      $              (INDMAX(icy)+igood).lt.ilast.and.  
      $              .true.)  
      $              ibad=BAD(VIEW(icy),  
      $              nvk(MAXS(icy)+igood),  
      $              nst(MAXS(icy)+igood))  
                badcl=badcl*ibad  
             enddo  
             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
                goto 20          !<<<<<<<<<<<<<< BAD cut  
             endif               !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
               
               
             cl_good(icy)=1                    
             nply=npl(VIEW(icy))  
             nldy=nld(MAXS(icy),VIEW(icy))  
               
 *     ----------------------------------------------  
 *     CONDITION TO FORM A COUPLE  
 *     ----------------------------------------------  
 *     geometrical consistency (same plane and ladder)  
             if(nply.eq.nplx.and.nldy.eq.nldx)then  
 *     charge correlation  
 *     ===========================================================  
 *     this version of the subroutine is used for the calibration  
 *     thus charge-correlation selection is obviously removed  
 *     ===========================================================  
 c$$$               ddd=(dedx(icy)  
 c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 *     ===========================================================  
                 
                 
 *     ------------------> COUPLE <------------------  
 *     check to do not overflow vector dimentions  
                if(ncp_plane(nplx).gt.ncouplemax)then  
                   if(DEBUG)print*,  
      $                    ' ** warning ** number of identified'//  
      $                    ' couples on plane ',nplx,  
      $                    ' exceeds vector dimention'//  
      $                    ' ( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
                   iflag=1  
                   return  
                endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(DEBUG)print*,  
      $                 '** warning ** number of identified '//  
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
                 
                ncp_plane(nplx) = ncp_plane(nplx) + 1  
                clx(nplx,ncp_plane(nplx))=icx  
                cly(nply,ncp_plane(nplx))=icy  
                cl_single(icx)=0  
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
   
  20         continue  
          enddo                  !end loop on clusters(Y)  
           
  10      continue  
       enddo                     !end loop on clusters(X)  
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
   
 c$$$      subroutine cl_to_couples_2(iflag)  
1982  c$$$  c$$$
1983  c$$$      include 'commontracker.f'  c$$$      include 'commontracker.f'
1984    c$$$      include 'level1.f'
1985  c$$$      include 'common_momanhough.f'  c$$$      include 'common_momanhough.f'
1986  c$$$      include 'momanhough_init.f'  c$$$c      include 'momanhough_init.f'
1987  c$$$      include 'calib.f'  c$$$      include 'calib.f'
1988  c$$$      include 'level1.f'  c$$$c      include 'level1.f'
1989  c$$$  c$$$
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
1990  c$$$  c$$$
1991  c$$$*     output flag  c$$$*     output flag
1992  c$$$*     --------------  c$$$*     --------------
# Line 2132  c$$$     $              nst(MAXS(icx)+ig Line 2046  c$$$     $              nst(MAXS(icx)+ig
2046  c$$$            endif  c$$$            endif
2047  c$$$            badcl=badcl*ibad  c$$$            badcl=badcl*ibad
2048  c$$$         enddo  c$$$         enddo
2049  c$$$*         print*,'icx ',icx,badcl  c$$$         if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut
2050  c$$$         if(badcl.eq.0)then  c$$$            cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut
2051  c$$$            cl_single(icx)=0  c$$$            goto 10             !<<<<<<<<<<<<<< BAD cut
2052  c$$$            goto 10  c$$$         endif                  !<<<<<<<<<<<<<< BAD cut
 c$$$         endif  
2053  c$$$*     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2054  c$$$          c$$$        
2055  c$$$         cl_good(icx)=1  c$$$         cl_good(icx)=1
# Line 2171  c$$$     $              nvk(MAXS(icy)+ig Line 2084  c$$$     $              nvk(MAXS(icy)+ig
2084  c$$$     $              nst(MAXS(icy)+igood))  c$$$     $              nst(MAXS(icy)+igood))
2085  c$$$               badcl=badcl*ibad  c$$$               badcl=badcl*ibad
2086  c$$$            enddo  c$$$            enddo
2087  c$$$*            print*,'icy ',icy,badcl  c$$$            if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut
2088  c$$$            if(badcl.eq.0)then  c$$$               cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut
2089  c$$$               cl_single(icy)=0  c$$$               goto 20          !<<<<<<<<<<<<<< BAD cut
2090  c$$$               goto 20  c$$$            endif               !<<<<<<<<<<<<<< BAD cut
 c$$$            endif  
2091  c$$$*     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2092  c$$$              c$$$            
2093  c$$$              c$$$            
# Line 2188  c$$$*     CONDITION TO FORM A COUPLE Line 2100  c$$$*     CONDITION TO FORM A COUPLE
2100  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
2101  c$$$*     geometrical consistency (same plane and ladder)  c$$$*     geometrical consistency (same plane and ladder)
2102  c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then
2103  c$$$  c$$$*     charge correlation
2104  c$$$c$$$*     charge correlation  c$$$*     ===========================================================
2105    c$$$*     this version of the subroutine is used for the calibration
2106    c$$$*     thus charge-correlation selection is obviously removed
2107    c$$$*     ===========================================================
2108  c$$$c$$$               ddd=(dedx(icy)  c$$$c$$$               ddd=(dedx(icy)
2109  c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
2110  c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
2111  c$$$c$$$               cut=chcut*sch(nplx,nldx)  c$$$c$$$               cut=chcut*sch(nplx,nldx)
2112  c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent
2113    c$$$*     ===========================================================
2114    c$$$              
2115  c$$$                c$$$              
2116  c$$$*     ------------------> COUPLE <------------------  c$$$*     ------------------> COUPLE <------------------
2117  c$$$*     check to do not overflow vector dimentions  c$$$*     check to do not overflow vector dimentions
2118  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  c$$$c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
2119  c$$$                  if(DEBUG)print*,  c$$$c$$$                  if(DEBUG)print*,
2120  c$$$     $                    ' ** warning ** number of identified'//  c$$$c$$$     $                    ' ** warning ** number of identified'//
2121  c$$$     $                    ' couples on plane ',nplx,  c$$$c$$$     $                    ' couples on plane ',nplx,
2122  c$$$     $                    ' exceeds vector dimention'//  c$$$c$$$     $                    ' exceeds vector dimention'//
2123  c$$$     $                    ' ( ',ncouplemax,' )'  c$$$c$$$     $                    ' ( ',ncouplemax,' )'
2124  c$$$c     good2=.false.  c$$$c$$$c     good2=.false.
2125  c$$$c     goto 880   !fill ntp and go to next event  c$$$c$$$c     goto 880   !fill ntp and go to next event
2126  c$$$                  iflag=1  c$$$c$$$                  iflag=1
2127  c$$$                  return  c$$$c$$$                  return
2128  c$$$               endif  c$$$c$$$               endif
2129  c$$$                c$$$              
2130  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then
2131  c$$$                  if(DEBUG)print*,  c$$$                  if(verbose)print*,
2132  c$$$     $                 '** warning ** number of identified '//  c$$$     $                 '** warning ** number of identified '//
2133  c$$$     $                 'couples on plane ',nplx,  c$$$     $                 'couples on plane ',nplx,
2134  c$$$     $                 'exceeds vector dimention '  c$$$     $                 'exceeds vector dimention '
# Line 2227  c$$$               clx(nplx,ncp_plane(np Line 2144  c$$$               clx(nplx,ncp_plane(np
2144  c$$$               cly(nply,ncp_plane(nplx))=icy  c$$$               cly(nply,ncp_plane(nplx))=icy
2145  c$$$               cl_single(icx)=0  c$$$               cl_single(icx)=0
2146  c$$$               cl_single(icy)=0  c$$$               cl_single(icy)=0
 c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy  
2147  c$$$            endif                                c$$$            endif                              
2148  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
2149  c$$$  c$$$
# Line 2260  c$$$      enddo Line 2176  c$$$      enddo
2176  c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
2177  c$$$        c$$$      
2178  c$$$      if(ncp_tot.gt.ncp_max)then  c$$$      if(ncp_tot.gt.ncp_max)then
2179  c$$$         if(DEBUG)print*,  c$$$         if(verbose)print*,
2180  c$$$     $           '** warning ** number of identified '//  c$$$     $           '** warning ** number of identified '//
2181  c$$$     $           'couples exceeds upper limit for Hough tr. '  c$$$     $           'couples exceeds upper limit for Hough tr. '
2182  c$$$     $           ,'( ',ncp_max,' )'              c$$$     $           ,'( ',ncp_max,' )'            
# Line 2272  c$$$      endif Line 2188  c$$$      endif
2188  c$$$        c$$$      
2189  c$$$      return  c$$$      return
2190  c$$$      end  c$$$      end
2191    c$$$
2192                
2193  ***************************************************  ***************************************************
2194  *                                                 *  *                                                 *
# Line 2288  c     02/02/2006 modified by Elena Vannu Line 2205  c     02/02/2006 modified by Elena Vannu
2205  c*****************************************************  c*****************************************************
2206    
2207        include 'commontracker.f'        include 'commontracker.f'
2208          include 'level1.f'
2209        include 'common_momanhough.f'        include 'common_momanhough.f'
2210        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2211        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2212        include 'common_mini_2.f'        include 'common_mini_2.f'
2213        include 'calib.f'        include 'calib.f'
2214        include 'level1.f'  c      include 'level1.f'
2215    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2216    
2217  *     output flag  *     output flag
2218  *     --------------  *     --------------
# Line 2364  c     $                       (icx2,icy2 Line 2280  c     $                       (icx2,icy2
2280  *     (2 couples needed)  *     (2 couples needed)
2281  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2282                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2283                             if(DEBUG)print*,                             if(verbose)print*,
2284       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2285       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2286       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2287  c                           good2=.false.  c                           good2=.false.
2288  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
2289                               do iv=1,12
2290                                  mask_view(iv) = 3
2291                               enddo
2292                             iflag=1                             iflag=1
2293                             return                             return
2294                          endif                          endif
# Line 2434  c     $                                 Line 2353  c     $                                
2353  *     (3 couples needed)  *     (3 couples needed)
2354  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2355                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2356                                      if(DEBUG)print*,                                      if(verbose)print*,
2357       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2358       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2359       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
2360  c                                    good2=.false.  c                                    good2=.false.
2361  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2362                                        do iv=1,nviews
2363                                           mask_view(iv) = 4
2364                                        enddo
2365                                      iflag=1                                      iflag=1
2366                                      return                                      return
2367                                   endif                                   endif
# Line 2514  c     goto 880               !ntp fill Line 2436  c     goto 880               !ntp fill
2436        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
2437    
2438        include 'commontracker.f'        include 'commontracker.f'
2439          include 'level1.f'
2440        include 'common_momanhough.f'        include 'common_momanhough.f'
2441        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2442    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2443    
2444  *     output flag  *     output flag
2445  *     --------------  *     --------------
# Line 2550  c      common/dbg/DEBUG Line 2471  c      common/dbg/DEBUG
2471        distance=0        distance=0
2472        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
2473        npt_tot=0        npt_tot=0
2474          nloop=0                  
2475     90   continue                  
2476        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2477           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
2478                            
# Line 2653  c     print*,'*   idbref,idb2 ',idbref,i Line 2576  c     print*,'*   idbref,idb2 ',idbref,i
2576              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2577           enddo           enddo
2578  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
2579           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
2580           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
2581           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2582                    
# Line 2661  c     print*,'>>>> ',ncpused,npt,nplused Line 2584  c     print*,'>>>> ',ncpused,npt,nplused
2584  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2585    
2586           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2587              if(DEBUG)print*,              if(verbose)print*,
2588       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2589       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2590       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2591  c               good2=.false.  c               good2=.false.
2592  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2593                do iv=1,nviews
2594                   mask_view(iv) = 5
2595                enddo
2596              iflag=1              iflag=1
2597              return              return
2598           endif           endif
# Line 2704  c$$$     $           ,(db_cloud(iii),iii Line 2630  c$$$     $           ,(db_cloud(iii),iii
2630        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2631                
2632                
2633          if(nloop.lt.nstepy)then      
2634            cutdistyz = cutdistyz+cutystep
2635            nloop     = nloop+1          
2636            goto 90                
2637          endif                    
2638          
2639        if(DEBUG)then        if(DEBUG)then
2640           print*,'---------------------- '           print*,'---------------------- '
2641           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2730  c$$$     $           ,(db_cloud(iii),iii Line 2662  c$$$     $           ,(db_cloud(iii),iii
2662        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2663    
2664        include 'commontracker.f'        include 'commontracker.f'
2665          include 'level1.f'
2666        include 'common_momanhough.f'        include 'common_momanhough.f'
2667        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2668    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2669    
2670  *     output flag  *     output flag
2671  *     --------------  *     --------------
# Line 2765  c      common/dbg/DEBUG Line 2696  c      common/dbg/DEBUG
2696        distance=0        distance=0
2697        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2698        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2699          nloop=0                  
2700     91   continue                  
2701        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2702           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
2703  c     print*,'--------------'  c     print*,'--------------'
# Line 2866  c     print*,'check cp_used' Line 2799  c     print*,'check cp_used'
2799           do ip=1,nplanes           do ip=1,nplanes
2800              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2801           enddo           enddo
2802           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2803           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2804           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next doublet
2805                    
2806  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2807  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2808           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2809              if(DEBUG)print*,              if(verbose)print*,
2810       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2811       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2812       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2813  c     good2=.false.  c     good2=.false.
2814  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2815                do iv=1,nviews
2816                   mask_view(iv) = 6
2817                enddo
2818              iflag=1              iflag=1
2819              return              return
2820           endif           endif
# Line 2914  c$$$     $           ,(tr_cloud(iii),iii Line 2850  c$$$     $           ,(tr_cloud(iii),iii
2850  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2851  22288    continue  22288    continue
2852        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2853          
2854           if(nloop.lt.nstepx)then      
2855             cutdistxz=cutdistxz+cutxstep
2856             nloop=nloop+1          
2857             goto 91                
2858           endif                    
2859          
2860        if(DEBUG)then        if(DEBUG)then
2861           print*,'---------------------- '           print*,'---------------------- '
2862           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2941  c     02/02/2006 modified by Elena Vannu Line 2883  c     02/02/2006 modified by Elena Vannu
2883  c*****************************************************  c*****************************************************
2884    
2885        include 'commontracker.f'        include 'commontracker.f'
2886          include 'level1.f'
2887        include 'common_momanhough.f'        include 'common_momanhough.f'
2888        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2889        include 'common_mini_2.f'        include 'common_mini_2.f'
2890        include 'common_mech.f'        include 'common_mech.f'
2891        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2892    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2893    
2894  *     output flag  *     output flag
2895  *     --------------  *     --------------
# Line 2964  c      common/dbg/DEBUG Line 2905  c      common/dbg/DEBUG
2905  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2906  *     list of matching couples in the combination  *     list of matching couples in the combination
2907  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2908        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2909        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2910  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2911        integer hit_plane(nplanes)        integer hit_plane(nplanes)
# Line 3064  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 3005  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
3005  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
3006                            
3007  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
3008              AL_INI(1)=dreal(alfaxz1_av(ixz))              AL_INI(1) = dreal(alfaxz1_av(ixz))
3009              AL_INI(2)=dreal(alfayz1_av(iyz))              AL_INI(2) = dreal(alfayz1_av(iyz))
3010              AL_INI(4)=datan(dreal(alfayz2_av(iyz))              AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
3011       $           /dreal(alfaxz2_av(ixz)))       $           /dreal(alfaxz2_av(ixz)))
3012              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))              tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
3013              AL_INI(3)=tath/sqrt(1+tath**2)              AL_INI(3) = tath/sqrt(1+tath**2)
3014              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.              AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
3015                            
3016  c     print*,'*******',AL_INI(5)  c     print*,'*******',AL_INI(5)
3017              if(AL_INI(5).gt.defmax)goto 888 !next cloud              if(AL_INI(5).gt.defmax)goto 888 !next cloud
# Line 3153  c     $                                 Line 3094  c     $                                
3094                                enddo                                enddo
3095                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3096                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3097                                call mini_2(jstep,ifail)                                iprint=0
3098                                  if(DEBUG)iprint=1
3099                                  call mini2(jstep,ifail,iprint)
3100                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3101                                   if(DEBUG)then                                   if(DEBUG)then
3102                                      print *,                                      print *,
3103       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3104       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(mini2 in clouds_to_ctrack)'
3105                                   endif                                   endif
3106                                   chi2=-chi2                                   chi2=-chi2
3107                                endif                                endif
# Line 3173  c     $                                 Line 3116  c     $                                
3116  *     --------------------------  *     --------------------------
3117                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3118                                                                    
3119                                   if(DEBUG)print*,                                   if(verbose)print*,
3120       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3121       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3122       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
3123  c                                 good2=.false.  c                                 good2=.false.
3124  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
3125                                     do iv=1,nviews
3126                                        mask_view(iv) = 7
3127                                     enddo
3128                                   iflag=1                                   iflag=1
3129                                   return                                   return
3130                                endif                                endif
# Line 3273  c$$$  rchi2=chi2/dble(ndof) Line 3219  c$$$  rchi2=chi2/dble(ndof)
3219  c******************************************************  c******************************************************
3220  cccccc 06/10/2005 modified by elena vannuccini ---> (1)  cccccc 06/10/2005 modified by elena vannuccini ---> (1)
3221  cccccc 31/01/2006 modified by elena vannuccini ---> (2)  cccccc 31/01/2006 modified by elena vannuccini ---> (2)
3222    cccccc 12/08/2006 modified by elena vannucicni ---> (3)
3223  c******************************************************  c******************************************************
3224    
3225        include 'commontracker.f'        include 'commontracker.f'
3226          include 'level1.f'
3227        include 'common_momanhough.f'        include 'common_momanhough.f'
3228        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3229        include 'common_mini_2.f'        include 'common_mini_2.f'
3230        include 'common_mech.f'        include 'common_mech.f'
3231        include 'momanhough_init.f'  c      include 'momanhough_init.f'
3232        include 'level1.f'  c      include 'level1.f'
3233        include 'calib.f'        include 'calib.f'
3234    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
3235    
3236  *     flag to chose PFA  *     flag to chose PFA
3237        character*10 PFA        character*10 PFA
# Line 3383  c            if(DEBUG)print*,'>>>> try t Line 3329  c            if(DEBUG)print*,'>>>> try t
3329                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3330                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3331                 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
3332       $              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
3333       $              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
3334         $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)
3335         $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
3336       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3337  *            *          
3338                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3456  c            if(DEBUG)print*,'>>>> try t Line 3404  c            if(DEBUG)print*,'>>>> try t
3404                 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
3405  *                                                !jump to the next couple  *                                                !jump to the next couple
3406  *----- try cluster x -----------------------------------------------  *----- try cluster x -----------------------------------------------
3407                 if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used  c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
3408                   if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
3409  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3410                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
3411  c     $              'ETA2','ETA2',  c     $              'ETA2','ETA2',
# Line 3483  c                  dedxmm = dedx(icx) !( Line 3432  c                  dedxmm = dedx(icx) !(
3432                 endif                                   endif                  
3433  11881          continue  11881          continue
3434  *----- try cluster y -----------------------------------------------  *----- try cluster y -----------------------------------------------
3435                 if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used  c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
3436                   if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
3437  *                                              !jump to the next couple  *                                              !jump to the next couple
3438                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
3439  c     $              'ETA2','ETA2',  c     $              'ETA2','ETA2',
# Line 3512  c                 dedxmm = dedx(icy)  !( Line 3462  c                 dedxmm = dedx(icy)  !(
3462  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------    
3463              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3464                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3465                 if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
3466                   if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
3467       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3468       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3469                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
# Line 3596  c              dedxtrk(nplanes-ip+1) = d Line 3547  c              dedxtrk(nplanes-ip+1) = d
3547  *                                                 *  *                                                 *
3548  *                                                 *  *                                                 *
3549  **************************************************  **************************************************
3550    cccccc 12/08/2006 modified by elena ---> (1)
3551    *
3552        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3553    
3554        include 'commontracker.f'        include 'commontracker.f'
3555          include 'level1.f'
3556        include 'common_momanhough.f'        include 'common_momanhough.f'
3557        include 'momanhough_init.f'  c      include 'momanhough_init.f'
3558          include 'level2.f'        !(1)
3559  c      include 'calib.f'  c      include 'calib.f'
3560  c      include 'level1.f'  c      include 'level1.f'
3561    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
3562    
3563    
3564        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3617  c      common/dbg/DEBUG Line 3569  c      common/dbg/DEBUG
3569              if(id.ne.0)then              if(id.ne.0)then
3570                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3571                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3572                 cl_used(iclx)=1  !tag used clusters  c               cl_used(iclx)=1  !tag used clusters
3573                 cl_used(icly)=1  !tag used clusters  c               cl_used(icly)=1  !tag used clusters
3574                   cl_used(iclx)=ntrk  !tag used clusters !(1)
3575                   cl_used(icly)=ntrk  !tag used clusters !(1)
3576              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3577                 cl_used(icl)=1   !tag used clusters  c               cl_used(icl)=1   !tag used clusters
3578                   cl_used(icl)=ntrk   !tag used clusters !1)
3579              endif              endif
3580                            
3581  c               if(DEBUG)then  c               if(DEBUG)then
# Line 3779  c$$$ Line 3734  c$$$
3734    
3735        subroutine init_level2        subroutine init_level2
3736    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3737        include 'commontracker.f'        include 'commontracker.f'
3738          include 'level1.f'
3739        include 'common_momanhough.f'        include 'common_momanhough.f'
3740        include 'level2.f'        include 'level2.f'
3741        include 'level1.f'  c      include 'level1.f'
   
   
3742    
3743        good2 = 0!.false.        do i=1,nviews
3744  c$$$      nev2 = nev1           good2(i)=good1(i)
3745          enddo
3746    
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
3747    
3748        NTRK = 0        NTRK = 0
3749        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3750           IMAGE(IT)=0           IMAGE(IT)=0
3751           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
          BdL(IT) = 0.  
3752           do ip=1,nplanes           do ip=1,nplanes
3753              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3754              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
# Line 3826  c*************************************** Line 3757  c***************************************
3757              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3758              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3759              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3760              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3761              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3762  c******************************************************              CLTRX(IP,IT) = 0
3763                CLTRY(IP,IT) = 0
3764           enddo           enddo
3765           do ipa=1,5           do ipa=1,5
3766              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3839  c*************************************** Line 3769  c***************************************
3769              enddo                                enddo                  
3770           enddo                             enddo                  
3771        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3772        nclsx=0        nclsx=0
3773        nclsy=0              nclsy=0      
3774        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3775          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3776          xs(1,ip)=0          xs(1,ip)=0
3777          xs(2,ip)=0          xs(2,ip)=0
3778          sgnlxs(ip)=0          sgnlxs(ip)=0
3779          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3780          ys(1,ip)=0          ys(1,ip)=0
3781          ys(2,ip)=0          ys(2,ip)=0
3782          sgnlys(ip)=0          sgnlys(ip)=0
3783        enddo        enddo
 c*******************************************************  
3784        end        end
3785    
3786    
# Line 3872  c*************************************** Line 3795  c***************************************
3795  ************************************************************  ************************************************************
3796    
3797    
3798          subroutine init_hough
3799    
3800          include 'commontracker.f'
3801          include 'level1.f'
3802          include 'common_momanhough.f'
3803          include 'common_hough.f'
3804          include 'level2.f'
3805    
3806          ntrpt_nt=0
3807          ndblt_nt=0
3808          NCLOUDS_XZ_nt=0
3809          NCLOUDS_YZ_nt=0
3810          do idb=1,ndblt_max_nt
3811             db_cloud_nt(idb)=0
3812             alfayz1_nt(idb)=0      
3813             alfayz2_nt(idb)=0      
3814          enddo
3815          do itr=1,ntrpl_max_nt
3816             tr_cloud_nt(itr)=0
3817             alfaxz1_nt(itr)=0      
3818             alfaxz2_nt(itr)=0      
3819             alfaxz3_nt(itr)=0      
3820          enddo
3821          do idb=1,ncloyz_max      
3822            ptcloud_yz_nt(idb)=0    
3823            alfayz1_av_nt(idb)=0    
3824            alfayz2_av_nt(idb)=0    
3825          enddo                    
3826          do itr=1,ncloxz_max      
3827            ptcloud_xz_nt(itr)=0    
3828            alfaxz1_av_nt(itr)=0    
3829            alfaxz2_av_nt(itr)=0    
3830            alfaxz3_av_nt(itr)=0    
3831          enddo                    
3832    
3833          ntrpt=0                  
3834          ndblt=0                  
3835          NCLOUDS_XZ=0              
3836          NCLOUDS_YZ=0              
3837          do idb=1,ndblt_max        
3838            db_cloud(idb)=0        
3839            cpyz1(idb)=0            
3840            cpyz2(idb)=0            
3841            alfayz1(idb)=0          
3842            alfayz2(idb)=0          
3843          enddo                    
3844          do itr=1,ntrpl_max        
3845            tr_cloud(itr)=0        
3846            cpxz1(itr)=0            
3847            cpxz2(itr)=0            
3848            cpxz3(itr)=0            
3849            alfaxz1(itr)=0          
3850            alfaxz2(itr)=0          
3851            alfaxz3(itr)=0          
3852          enddo                    
3853          do idb=1,ncloyz_max      
3854            ptcloud_yz(idb)=0      
3855            alfayz1_av(idb)=0      
3856            alfayz2_av(idb)=0      
3857            do idbb=1,ncouplemaxtot
3858              cpcloud_yz(idb,idbb)=0
3859            enddo                  
3860          enddo                    
3861          do itr=1,ncloxz_max      
3862            ptcloud_xz(itr)=0      
3863            alfaxz1_av(itr)=0      
3864            alfaxz2_av(itr)=0      
3865            alfaxz3_av(itr)=0      
3866            do itrr=1,ncouplemaxtot
3867              cpcloud_xz(itr,itrr)=0
3868            enddo                  
3869          enddo                    
3870          end
3871    ************************************************************
3872    *
3873    *
3874    *
3875    *
3876    *
3877    *
3878    *
3879    ************************************************************
3880    
3881    
3882        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3883    
3884  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3882  c*************************************** Line 3889  c***************************************
3889    
3890            
3891        include 'commontracker.f'        include 'commontracker.f'
3892    c      include 'level1.f'
3893          include 'level1.f'
3894          include 'common_momanhough.f'
3895        include 'level2.f'        include 'level2.f'
3896        include 'common_mini_2.f'        include 'common_mini_2.f'
3897        real sinth,phi,pig        !(4)        real sinth,phi,pig      
3898        pig=acos(-1.)        pig=acos(-1.)
3899    
       good2=1!.true.  
3900        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3901          nstep_nt(ntr)       = nstep
3902    
3903          phi   = al(4)          
3904          sinth = al(3)            
3905          if(sinth.lt.0)then      
3906             sinth = -sinth        
3907             phi = phi + pig      
3908          endif                    
3909          npig = aint(phi/(2*pig))
3910          phi = phi - npig*2*pig  
3911          if(phi.lt.0)            
3912         $     phi = phi + 2*pig  
3913          al(4) = phi              
3914          al(3) = sinth            
3915    
       phi   = al(4)             !(4)  
       sinth = al(3)             !(4)  
       if(sinth.lt.0)then        !(4)  
          sinth = -sinth         !(4)  
          phi = phi + pig        !(4)  
       endif                     !(4)  
       npig = aint(phi/(2*pig))  !(4)  
       phi = phi - npig*2*pig    !(4)  
       if(phi.lt.0)              !(4)  
      $     phi = phi + 2*pig    !(4)  
       al(4) = phi               !(4)  
       al(3) = sinth             !(4)  
 *****************************************************  
3916        do i=1,5        do i=1,5
3917           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3918           do j=1,5           do j=1,5
3919              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3920           enddo           enddo
 c     print*,al_nt(i,ntr)  
3921        enddo        enddo
3922                
3923        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
# Line 3924  c     print*,al_nt(i,ntr) Line 3933  c     print*,al_nt(i,ntr)
3933           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3934           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3935           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))
 c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3936           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
3937           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)                     dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)  
3938      
3939             id  = CP_STORE(ip,IDCAND)
3940             icl = CLS_STORE(ip,IDCAND)
3941             if(id.ne.0)then
3942                cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3943                cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3944    c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
3945             elseif(icl.ne.0)then
3946                if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl
3947                if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl
3948    c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
3949             endif          
3950    
3951        enddo        enddo
 c      call CalcBdL(100,xxxx,IFAIL)  
 c      if(ifps(xxxx).eq.1)BdL(ntr) = xxxx  
 c$$$      print*,'xgood(ip,ntr) ',(xgood_nt(ip,ntr),ip=1,6)  
 c$$$      print*,'ygood(ip,ntr) ',(ygood_nt(ip,ntr),ip=1,6)  
 c$$$      print*,'dedx_x(ip,ntr) ',(dedx_x(ip,ntr),ip=1,6)  
 c$$$      print*,'dedx_y(ip,ntr) ',(dedx_y(ip,ntr),ip=1,6)  
3952    
3953    
3954        end        end
3955    
3956        subroutine fill_level2_siglets        subroutine fill_level2_siglets
 c*****************************************************  
 c     07/10/2005 created by elena vannuccini  
 c     31/01/2006 modified by elena vannuccini  
 *     to convert adc to mip  --> (2)  
 c*****************************************************  
3957    
3958  *     -------------------------------------------------------  *     -------------------------------------------------------
3959  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3952  c*************************************** Line 3962  c***************************************
3962  *     -------------------------------------------------------  *     -------------------------------------------------------
3963    
3964        include 'commontracker.f'        include 'commontracker.f'
3965        include 'level1.f'  c      include 'level1.f'
       include 'level2.f'  
3966        include 'calib.f'        include 'calib.f'
3967          include 'level1.f'
3968        include 'common_momanhough.f'        include 'common_momanhough.f'
3969          include 'level2.f'
3970        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3971    
3972  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
3973        good2=1!.true.  c      good2=1!.true.
3974        nclsx = 0        nclsx = 0
3975        nclsy = 0        nclsy = 0
3976    
3977          do iv = 1,nviews
3978             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3979          enddo
3980    
3981        do icl=1,nclstr1        do icl=1,nclstr1
3982           if(cl_used(icl).eq.0)then !cluster not included in any track           if(cl_used(icl).eq.0)then !cluster not included in any track
3983              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
# Line 3970  c*************************************** Line 3985  c***************************************
3985                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3986                 planex(nclsx) = ip                 planex(nclsx) = ip
3987                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3988                   clsx(nclsx)   = icl
3989                 do is=1,2                 do is=1,2
3990  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3991                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
# Line 3984  c$$$               print*,'xs(2,nclsx)   Line 4000  c$$$               print*,'xs(2,nclsx)  
4000                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4001                 planey(nclsy) = ip                 planey(nclsy) = ip
4002                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
4003                   clsy(nclsy)   = icl
4004                 do is=1,2                 do is=1,2
4005  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4006                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
# Line 3997  c$$$               print*,'ys(2,nclsy)   Line 4014  c$$$               print*,'ys(2,nclsy)  
4014              endif              endif
4015           endif           endif
4016  c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
4017    
4018    ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4019             whichtrack(icl) = cl_used(icl)
4020    
4021        enddo        enddo
4022        end        end
4023    
4024    ***************************************************
4025    *                                                 *
4026    *                                                 *
4027    *                                                 *
4028    *                                                 *
4029    *                                                 *
4030    *                                                 *
4031    **************************************************
4032    
4033          subroutine fill_hough
4034    
4035    *     -------------------------------------------------------
4036    *     This routine fills the  variables related to the hough
4037    *     transform, for the debig n-tuple
4038    *     -------------------------------------------------------
4039    
4040          include 'commontracker.f'
4041          include 'level1.f'
4042          include 'common_momanhough.f'
4043          include 'common_hough.f'
4044          include 'level2.f'
4045    
4046          if(.false.
4047         $     .or.ntrpt.gt.ntrpt_max_nt
4048         $     .or.ndblt.gt.ndblt_max_nt
4049         $     .or.NCLOUDS_XZ.gt.ncloxz_max
4050         $     .or.NCLOUDS_yZ.gt.ncloyz_max
4051         $     )then
4052             ntrpt_nt=0
4053             ndblt_nt=0
4054             NCLOUDS_XZ_nt=0
4055             NCLOUDS_YZ_nt=0
4056          else
4057             ndblt_nt=ndblt
4058             ntrpt_nt=ntrpt
4059             if(ndblt.ne.0)then
4060                do id=1,ndblt
4061                   alfayz1_nt(id)=alfayz1(id) !Y0
4062                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
4063                enddo
4064             endif
4065             if(ndblt.ne.0)then
4066                do it=1,ntrpt
4067                   alfaxz1_nt(it)=alfaxz1(it) !X0
4068                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
4069                   alfaxz3_nt(it)=alfaxz3(it) !1/r
4070                enddo
4071             endif
4072             nclouds_yz_nt=nclouds_yz
4073             nclouds_xz_nt=nclouds_xz
4074             if(nclouds_yz.ne.0)then
4075                nnn=0
4076                do iyz=1,nclouds_yz
4077                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
4078                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
4079                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4080                   nnn=nnn+ptcloud_yz(iyz)
4081                enddo
4082                do ipt=1,nnn
4083                   db_cloud_nt(ipt)=db_cloud(ipt)
4084                 enddo
4085             endif
4086             if(nclouds_xz.ne.0)then
4087                nnn=0
4088                do ixz=1,nclouds_xz
4089                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4090                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4091                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4092                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4093                   nnn=nnn+ptcloud_xz(ixz)              
4094                enddo
4095                do ipt=1,nnn
4096                  tr_cloud_nt(ipt)=tr_cloud(ipt)
4097                 enddo
4098             endif
4099          endif
4100          end
4101          
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23