/[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.4 by pam-fi, Thu Sep 28 14:04:40 2006 UTC revision 1.13 by pam-fi, Wed Nov 8 16:42:28 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          real*8 AL_GUESS(5)
187    
188  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
189  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 198  c         iflag=0 Line 224  c         iflag=0
224           ibest=0                !best track among candidates           ibest=0                !best track among candidates
225           iimage=0               !track image           iimage=0               !track image
226  *     ------------- select the best track -------------  *     ------------- select the best track -------------
227    c$$$         rchi2best=1000000000.
228    c$$$         do i=1,ntracks
229    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
230    c$$$     $         RCHI2_STORE(i).gt.0)then
231    c$$$               ibest=i
232    c$$$               rchi2best=RCHI2_STORE(i)
233    c$$$            endif
234    c$$$         enddo
235    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
236    
237           rchi2best=1000000000.           rchi2best=1000000000.
238             ndofbest=0             !(1)
239           do i=1,ntracks           do i=1,ntracks
240              if(RCHI2_STORE(i).lt.rchi2best.and.             if(RCHI2_STORE(i).lt.rchi2best.and.
241       $         RCHI2_STORE(i).gt.0)then       $          RCHI2_STORE(i).gt.0)then
242                 ndof=0             !(1)
243                 do ii=1,nplanes    !(1)
244                   ndof=ndof        !(1)
245         $              +int(xgood_store(ii,i)) !(1)
246         $              +int(ygood_store(ii,i)) !(1)
247                 enddo              !(1)
248                 if(ndof.ge.ndofbest)then !(1)
249                 ibest=i                 ibest=i
250                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
251              endif                 ndofbest=ndof    !(1)
252                 endif              !(1)
253               endif
254           enddo           enddo
255           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
256  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
# Line 236  c         iflag=0 Line 282  c         iflag=0
282  *     **********************************************************  *     **********************************************************
283  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
284  *     **********************************************************  *     **********************************************************
285             call guess()
286             do i=1,5
287                AL_GUESS(i)=AL(i)
288             enddo
289    
290           do i=1,5           do i=1,5
291              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
292           enddo           enddo
293            
294           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
295           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
296           jstep=0                !# minimization steps           jstep=0                !# minimization steps
297    
298           call mini_2(jstep,ifail)           iprint=0
299    c         if(DEBUG)iprint=1
300             if(VERBOSE)iprint=1
301             call mini2(jstep,ifail,iprint)
302           if(ifail.ne.0) then           if(ifail.ne.0) then
303              if(DEBUG)then              if(.true.)then
304                 print *,                 print *,
305       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
306       $              ,iev       $              ,iev
307    
308    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
309    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
310    c$$$               print*,'result:   ',(al(i),i=1,5)
311    c$$$               print*,'xgood ',xgood
312    c$$$               print*,'ygood ',ygood
313    c$$$               print*,'----------------------------------------------'
314              endif              endif
315              chi2=-chi2  c            chi2=-chi2
316           endif           endif
317                    
318           if(DEBUG)then           if(DEBUG)then
# Line 311  c         print*,'++++++++++ iimage,fima Line 373  c         print*,'++++++++++ iimage,fima
373  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
374    
375           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
376              if(DEBUG)              if(verbose)
377       $           print*,       $           print*,
378       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
379       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 607  c                (implemented variable r Line 669  c                (implemented variable r
669  c*****************************************************  c*****************************************************
670                
671        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
672        include 'level1.f'        include 'level1.f'
673          include 'calib.f'
674    c      include 'level1.f'
675        include 'common_align.f'        include 'common_align.f'
676        include 'common_mech.f'        include 'common_mech.f'
677        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
678        include 'common_resxy.f'  c      include 'common_resxy.f'
679    
680  c      logical DEBUG  c      logical DEBUG
681  c      common/dbg/DEBUG  c      common/dbg/DEBUG
# Line 667  c      double precision xi_B,yi_B,zi_B Line 730  c      double precision xi_B,yi_B,zi_B
730              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
731           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
732  c            cog2 = cog(2,icx)  c            cog2 = cog(2,icx)
733  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)              c            etacorr = pfaeta2(cog2,viewx,nldx,angx)            
734  c            stripx = stripx + etacorr  c            stripx = stripx + etacorr
735              stripx = stripx + pfa_eta2(icx,angx)            !(3)              stripx = stripx + pfaeta2(icx,angx)            !(3)
736              resxPAM = risx_eta2(angx)                       !   (4)              resxPAM = risx_eta2(angx)                       !   (4)
737              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
738       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
739              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
740           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                         !(3)
741              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)            !(3)
742              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(angx)                       !   (4)
743              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
744       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
745              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
746           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                         !(3)
747              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            !(3)
748              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(angx)                       !   (4)
749              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
750       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
751              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
752           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then                          !(3)
753              stripx = stripx + pfa_eta(icx,angx)             !(3)              stripx = stripx + pfaeta(icx,angx)             !(3)
754              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = ris_eta(icx,angx)                     !   (4)
755              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
756       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
# Line 731  c     $     print*,PFAx,icx,angx,stripx, Line 794  c     $     print*,PFAx,icx,angx,stripx,
794              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
795           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
796  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
797  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
798  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
799              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
800              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(angy)                       !   (4)
801              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
802              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
803       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
804           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
805              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
806              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
807              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
808       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
809           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
810              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
811              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
812              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
813       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
814           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
815              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
816              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
817  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
818              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 1284  c     $              ,iv,xvv(iv),yvv(iv) Line 1347  c     $              ,iv,xvv(iv),yvv(iv)
1347  *     it returns the plane number  *     it returns the plane number
1348  *      *    
1349        include 'commontracker.f'        include 'commontracker.f'
1350          include 'level1.f'
1351  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1352        include 'common_momanhough.f'        include 'common_momanhough.f'
1353                
# Line 1321  c      include 'common_analysis.f' Line 1385  c      include 'common_analysis.f'
1385  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1386  *      *    
1387        include 'commontracker.f'        include 'commontracker.f'
1388          include 'level1.f'
1389  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1390        include 'common_momanhough.f'        include 'common_momanhough.f'
1391                
# Line 1349  c      include 'common_analysis.f' Line 1414  c      include 'common_analysis.f'
1414  *     positive if sensor =2  *     positive if sensor =2
1415  *  *
1416        include 'commontracker.f'        include 'commontracker.f'
1417          include 'level1.f'
1418  c      include 'calib.f'  c      include 'calib.f'
1419  c      include 'level1.f'  c      include 'level1.f'
1420  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1660  c$$$      end Line 1726  c$$$      end
1726        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1727    
1728        include 'commontracker.f'        include 'commontracker.f'
1729          include 'level1.f'
1730        include 'common_momanhough.f'        include 'common_momanhough.f'
1731        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1732        include 'calib.f'        include 'calib.f'
1733        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1734    
1735  *     output flag  *     output flag
1736  *     --------------  *     --------------
# Line 1675  c      common/dbg/DEBUG Line 1739  c      common/dbg/DEBUG
1739  *     --------------  *     --------------
1740        integer iflag        integer iflag
1741    
1742        integer badseed,badcl        integer badseed,badclx,badcly
1743    
1744  *     init variables  *     init variables
1745        ncp_tot=0        ncp_tot=0
# Line 1691  c      common/dbg/DEBUG Line 1755  c      common/dbg/DEBUG
1755           ncls(ip)=0           ncls(ip)=0
1756        enddo        enddo
1757        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1758           cl_single(icl)=1           cl_single(icl) = 1
1759           cl_good(icl)=0           cl_good(icl)   = 0
1760          enddo
1761          do iv=1,nviews
1762             ncl_view(iv)  = 0
1763             mask_view(iv) = 0      !all included
1764        enddo        enddo
1765                
1766    *     count number of cluster per view
1767          do icl=1,nclstr1
1768             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1769          enddo
1770    *     mask views with too many clusters
1771          do iv=1,nviews
1772             if( ncl_view(iv).gt. nclustermax)then
1773                mask_view(iv) = 1
1774                if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1775         $           ,nclustermax,' on view ', iv,' --> masked!'
1776             endif
1777          enddo
1778    
1779    
1780  *     start association  *     start association
1781        ncouples=0        ncouples=0
1782        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1783           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1784                    
1785  *     ----------------------------------------------------  *     ----------------------------------------------------
1786    *     jump masked views (X VIEW)
1787    *     ----------------------------------------------------
1788             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1789    *     ----------------------------------------------------
1790  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1791  *     ----------------------------------------------------  *     ----------------------------------------------------
1792           if(dedx(icx).lt.dedx_x_min)then           if(dedx(icx).lt.dedx_x_min)then
# Line 1717  c      common/dbg/DEBUG Line 1803  c      common/dbg/DEBUG
1803           else           else
1804              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1805           endif           endif
1806           badcl=badseed           badclx=badseed
1807           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1808              ibad=1              ibad=1
1809              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1727  c      common/dbg/DEBUG Line 1813  c      common/dbg/DEBUG
1813       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1814       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1815              endif              endif
1816              badcl=badcl*ibad              badclx=badclx*ibad
1817           enddo           enddo
1818  *     ----------------------------------------------------  *     ----------------------------------------------------
1819  *     >>> eliminato il taglio sulle BAD <<<  *     >>> eliminato il taglio sulle BAD <<<
# Line 1746  c     endif Line 1832  c     endif
1832              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1833                            
1834  *     ----------------------------------------------------  *     ----------------------------------------------------
1835    *     jump masked views (Y VIEW)
1836    *     ----------------------------------------------------
1837                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1838    
1839    *     ----------------------------------------------------
1840  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1841  *     ----------------------------------------------------  *     ----------------------------------------------------
1842              if(dedx(icy).lt.dedx_y_min)then              if(dedx(icy).lt.dedx_y_min)then
# Line 1762  c     endif Line 1853  c     endif
1853              else              else
1854                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1855              endif              endif
1856              badcl=badseed              badcly=badseed
1857              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1858                 ibad=1                 ibad=1
1859                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1771  c     endif Line 1862  c     endif
1862       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1863       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1864       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1865                 badcl=badcl*ibad                 badcly=badcly*ibad
1866              enddo              enddo
1867  *     ----------------------------------------------------  *     ----------------------------------------------------
1868  *     >>> eliminato il taglio sulle BAD <<<  *     >>> eliminato il taglio sulle BAD <<<
# Line 1794  c     endif Line 1885  c     endif
1885  *     charge correlation  *     charge correlation
1886  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1887    
1888  *     -------------------------------------------------------------                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1889  *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<       $              .and.
1890  *     -------------------------------------------------------------       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1891  c$$$               if(dedx(icy).lt.chsaty.or.dedx(icx).lt.chsatx)then       $              .and.
1892  c$$$                  ddd=(dedx(icy)       $              (badclx.eq.1.and.badcly.eq.1)
1893  c$$$     $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $              .and.
1894  c$$$                  ddd=ddd/sqrt(kch(nplx,nldx)**2+1)       $              .true.)then
1895  c$$$                  cut=chcut*sch(nplx,nldx)  
1896  c$$$                  if(abs(ddd).gt.cut)goto 20 !charge not consistent                    ddd=(dedx(icy)
1897  c$$$               endif       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1898                      ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1899    
1900    c                  cut = chcut * sch(nplx,nldx)
1901    
1902                      sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1903         $                 -kch(nplx,nldx)*cch(nplx,nldx))
1904                      sss=sss/sqrt(kch(nplx,nldx)**2+1)
1905                      cut = chcut * (16 + sss/50.)
1906    
1907                      if(abs(ddd).gt.cut)then
1908                         goto 20    !charge not consistent
1909                      endif
1910                   endif
1911                                
1912  *     ------------------> COUPLE <------------------  *     ------------------> COUPLE <------------------
1913  *     check to do not overflow vector dimentions  *     check to do not overflow vector dimentions
1914                 if(ncp_plane(nplx).gt.ncouplemax)then  c$$$               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  
                 
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
1915  c$$$                  if(DEBUG)print*,  c$$$                  if(DEBUG)print*,
1916  c$$$     $                 '** warning ** number of identified '//  c$$$     $                    ' ** warning ** number of identified'//
1917  c$$$     $                 'couples on plane ',nplx,  c$$$     $                    ' couples on plane ',nplx,
1918  c$$$     $                 'exceeds vector dimention '  c$$$     $                    ' exceeds vector dimention'//
1919  c$$$     $                 ,'( ',ncouplemax,' )'  c$$$     $                    ' ( ',ncouplemax,' )'
1920  c$$$c     good2=.false.  c$$$c     good2=.false.
1921  c$$$c     goto 880   !fill ntp and go to next event                      c$$$c     goto 880   !fill ntp and go to next event
1922  c$$$                  iflag=1  c$$$                  iflag=1
1923  c$$$                  return  c$$$                  return
1924  c$$$               endif  c$$$               endif
1925                                
1926                   if(ncp_plane(nplx).eq.ncouplemax)then
1927                      if(verbose)print*,
1928         $                 '** warning ** number of identified '//
1929         $                 'couples on plane ',nplx,
1930         $                 'exceeds vector dimention '
1931         $                 ,'( ',ncouplemax,' )'
1932    c     good2=.false.
1933    c     goto 880   !fill ntp and go to next event
1934                      mask_view(nviewx(nplx)) = 2
1935                      mask_view(nviewy(nply)) = 2
1936    c                  iflag=1
1937    c                  return
1938                   endif
1939                  
1940                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1941                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1942                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
# Line 1863  c$$$               endif Line 1969  c$$$               endif
1969        endif        endif
1970                
1971        do ip=1,6        do ip=1,6
1972           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1973        enddo        enddo
1974  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
1975                
1976        if(ncp_tot.gt.ncp_max)then  c$$$      if(ncp_tot.gt.ncp_max)then
1977           if(DEBUG)print*,  c$$$         if(verbose)print*,
1978       $           '** warning ** number of identified '//  c$$$     $           '** warning ** number of identified '//
1979       $           'couples exceeds upper limit for Hough tr. '  c$$$     $           'couples exceeds upper limit for Hough tr. '
1980       $           ,'( ',ncp_max,' )'              c$$$     $           ,'( ',ncp_max,' )'            
1981  c            good2=.false.  c$$$         iflag=1
1982  c     goto 880       !fill ntp and go to next event  c$$$         return
1983           iflag=1  c$$$      endif
          return  
       endif  
1984                
1985        return        return
1986        end        end
# Line 1889  c     goto 880       !fill ntp and go to Line 1993  c     goto 880       !fill ntp and go to
1993  *                                                 *  *                                                 *
1994  *                                                 *  *                                                 *
1995  **************************************************  **************************************************
1996        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  
                 
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               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)  
1997  c$$$  c$$$
1998  c$$$      include 'commontracker.f'  c$$$      include 'commontracker.f'
1999    c$$$      include 'level1.f'
2000  c$$$      include 'common_momanhough.f'  c$$$      include 'common_momanhough.f'
2001  c$$$      include 'momanhough_init.f'  c$$$c      include 'momanhough_init.f'
2002  c$$$      include 'calib.f'  c$$$      include 'calib.f'
2003  c$$$      include 'level1.f'  c$$$c      include 'level1.f'
2004  c$$$  c$$$
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
2005  c$$$  c$$$
2006  c$$$*     output flag  c$$$*     output flag
2007  c$$$*     --------------  c$$$*     --------------
# Line 2170  c$$$     $              nst(MAXS(icx)+ig Line 2061  c$$$     $              nst(MAXS(icx)+ig
2061  c$$$            endif  c$$$            endif
2062  c$$$            badcl=badcl*ibad  c$$$            badcl=badcl*ibad
2063  c$$$         enddo  c$$$         enddo
2064  c$$$*         print*,'icx ',icx,badcl  c$$$         if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut
2065  c$$$         if(badcl.eq.0)then  c$$$            cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut
2066  c$$$            cl_single(icx)=0  c$$$            goto 10             !<<<<<<<<<<<<<< BAD cut
2067  c$$$            goto 10  c$$$         endif                  !<<<<<<<<<<<<<< BAD cut
 c$$$         endif  
2068  c$$$*     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2069  c$$$          c$$$        
2070  c$$$         cl_good(icx)=1  c$$$         cl_good(icx)=1
# Line 2209  c$$$     $              nvk(MAXS(icy)+ig Line 2099  c$$$     $              nvk(MAXS(icy)+ig
2099  c$$$     $              nst(MAXS(icy)+igood))  c$$$     $              nst(MAXS(icy)+igood))
2100  c$$$               badcl=badcl*ibad  c$$$               badcl=badcl*ibad
2101  c$$$            enddo  c$$$            enddo
2102  c$$$*            print*,'icy ',icy,badcl  c$$$            if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut
2103  c$$$            if(badcl.eq.0)then  c$$$               cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut
2104  c$$$               cl_single(icy)=0  c$$$               goto 20          !<<<<<<<<<<<<<< BAD cut
2105  c$$$               goto 20  c$$$            endif               !<<<<<<<<<<<<<< BAD cut
 c$$$            endif  
2106  c$$$*     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2107  c$$$              c$$$            
2108  c$$$              c$$$            
# Line 2226  c$$$*     CONDITION TO FORM A COUPLE Line 2115  c$$$*     CONDITION TO FORM A COUPLE
2115  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
2116  c$$$*     geometrical consistency (same plane and ladder)  c$$$*     geometrical consistency (same plane and ladder)
2117  c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then
2118  c$$$  c$$$*     charge correlation
2119  c$$$c$$$*     charge correlation  c$$$*     ===========================================================
2120    c$$$*     this version of the subroutine is used for the calibration
2121    c$$$*     thus charge-correlation selection is obviously removed
2122    c$$$*     ===========================================================
2123  c$$$c$$$               ddd=(dedx(icy)  c$$$c$$$               ddd=(dedx(icy)
2124  c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
2125  c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
2126  c$$$c$$$               cut=chcut*sch(nplx,nldx)  c$$$c$$$               cut=chcut*sch(nplx,nldx)
2127  c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent
2128    c$$$*     ===========================================================
2129    c$$$              
2130  c$$$                c$$$              
2131  c$$$*     ------------------> COUPLE <------------------  c$$$*     ------------------> COUPLE <------------------
2132  c$$$*     check to do not overflow vector dimentions  c$$$*     check to do not overflow vector dimentions
2133  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  c$$$c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
2134  c$$$                  if(DEBUG)print*,  c$$$c$$$                  if(DEBUG)print*,
2135  c$$$     $                    ' ** warning ** number of identified'//  c$$$c$$$     $                    ' ** warning ** number of identified'//
2136  c$$$     $                    ' couples on plane ',nplx,  c$$$c$$$     $                    ' couples on plane ',nplx,
2137  c$$$     $                    ' exceeds vector dimention'//  c$$$c$$$     $                    ' exceeds vector dimention'//
2138  c$$$     $                    ' ( ',ncouplemax,' )'  c$$$c$$$     $                    ' ( ',ncouplemax,' )'
2139  c$$$c     good2=.false.  c$$$c$$$c     good2=.false.
2140  c$$$c     goto 880   !fill ntp and go to next event  c$$$c$$$c     goto 880   !fill ntp and go to next event
2141  c$$$                  iflag=1  c$$$c$$$                  iflag=1
2142  c$$$                  return  c$$$c$$$                  return
2143  c$$$               endif  c$$$c$$$               endif
2144  c$$$                c$$$              
2145  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then
2146  c$$$                  if(DEBUG)print*,  c$$$                  if(verbose)print*,
2147  c$$$     $                 '** warning ** number of identified '//  c$$$     $                 '** warning ** number of identified '//
2148  c$$$     $                 'couples on plane ',nplx,  c$$$     $                 'couples on plane ',nplx,
2149  c$$$     $                 'exceeds vector dimention '  c$$$     $                 'exceeds vector dimention '
# Line 2265  c$$$               clx(nplx,ncp_plane(np Line 2159  c$$$               clx(nplx,ncp_plane(np
2159  c$$$               cly(nply,ncp_plane(nplx))=icy  c$$$               cly(nply,ncp_plane(nplx))=icy
2160  c$$$               cl_single(icx)=0  c$$$               cl_single(icx)=0
2161  c$$$               cl_single(icy)=0  c$$$               cl_single(icy)=0
 c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy  
2162  c$$$            endif                                c$$$            endif                              
2163  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
2164  c$$$  c$$$
# Line 2298  c$$$      enddo Line 2191  c$$$      enddo
2191  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!!!)
2192  c$$$        c$$$      
2193  c$$$      if(ncp_tot.gt.ncp_max)then  c$$$      if(ncp_tot.gt.ncp_max)then
2194  c$$$         if(DEBUG)print*,  c$$$         if(verbose)print*,
2195  c$$$     $           '** warning ** number of identified '//  c$$$     $           '** warning ** number of identified '//
2196  c$$$     $           'couples exceeds upper limit for Hough tr. '  c$$$     $           'couples exceeds upper limit for Hough tr. '
2197  c$$$     $           ,'( ',ncp_max,' )'              c$$$     $           ,'( ',ncp_max,' )'            
# Line 2310  c$$$      endif Line 2203  c$$$      endif
2203  c$$$        c$$$      
2204  c$$$      return  c$$$      return
2205  c$$$      end  c$$$      end
2206    c$$$
2207                
2208  ***************************************************  ***************************************************
2209  *                                                 *  *                                                 *
# Line 2326  c     02/02/2006 modified by Elena Vannu Line 2220  c     02/02/2006 modified by Elena Vannu
2220  c*****************************************************  c*****************************************************
2221    
2222        include 'commontracker.f'        include 'commontracker.f'
2223          include 'level1.f'
2224        include 'common_momanhough.f'        include 'common_momanhough.f'
2225        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2226        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2227        include 'common_mini_2.f'        include 'common_mini_2.f'
2228        include 'calib.f'        include 'calib.f'
2229        include 'level1.f'  c      include 'level1.f'
2230    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2231    
2232  *     output flag  *     output flag
2233  *     --------------  *     --------------
# Line 2402  c     $                       (icx2,icy2 Line 2295  c     $                       (icx2,icy2
2295  *     (2 couples needed)  *     (2 couples needed)
2296  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2297                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2298                             if(DEBUG)print*,                             if(verbose)print*,
2299       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2300       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2301       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2302  c                           good2=.false.  c                           good2=.false.
2303  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
2304                               do iv=1,12
2305                                  mask_view(iv) = 3
2306                               enddo
2307                             iflag=1                             iflag=1
2308                             return                             return
2309                          endif                          endif
# Line 2472  c     $                                 Line 2368  c     $                                
2368  *     (3 couples needed)  *     (3 couples needed)
2369  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2370                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2371                                      if(DEBUG)print*,                                      if(verbose)print*,
2372       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2373       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2374       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
2375  c                                    good2=.false.  c                                    good2=.false.
2376  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2377                                        do iv=1,nviews
2378                                           mask_view(iv) = 4
2379                                        enddo
2380                                      iflag=1                                      iflag=1
2381                                      return                                      return
2382                                   endif                                   endif
# Line 2552  c     goto 880               !ntp fill Line 2451  c     goto 880               !ntp fill
2451        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
2452    
2453        include 'commontracker.f'        include 'commontracker.f'
2454          include 'level1.f'
2455        include 'common_momanhough.f'        include 'common_momanhough.f'
2456        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2457    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2458    
2459  *     output flag  *     output flag
2460  *     --------------  *     --------------
# Line 2588  c      common/dbg/DEBUG Line 2486  c      common/dbg/DEBUG
2486        distance=0        distance=0
2487        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
2488        npt_tot=0        npt_tot=0
2489          nloop=0                  
2490     90   continue                  
2491        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2492           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
2493                            
# Line 2691  c     print*,'*   idbref,idb2 ',idbref,i Line 2591  c     print*,'*   idbref,idb2 ',idbref,i
2591              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2592           enddo           enddo
2593  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
2594           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
2595           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
2596           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2597                    
# Line 2699  c     print*,'>>>> ',ncpused,npt,nplused Line 2599  c     print*,'>>>> ',ncpused,npt,nplused
2599  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2600    
2601           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2602              if(DEBUG)print*,              if(verbose)print*,
2603       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2604       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2605       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2606  c               good2=.false.  c               good2=.false.
2607  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2608                do iv=1,nviews
2609                   mask_view(iv) = 5
2610                enddo
2611              iflag=1              iflag=1
2612              return              return
2613           endif           endif
# Line 2742  c$$$     $           ,(db_cloud(iii),iii Line 2645  c$$$     $           ,(db_cloud(iii),iii
2645        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2646                
2647                
2648          if(nloop.lt.nstepy)then      
2649            cutdistyz = cutdistyz+cutystep
2650            nloop     = nloop+1          
2651            goto 90                
2652          endif                    
2653          
2654        if(DEBUG)then        if(DEBUG)then
2655           print*,'---------------------- '           print*,'---------------------- '
2656           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2768  c$$$     $           ,(db_cloud(iii),iii Line 2677  c$$$     $           ,(db_cloud(iii),iii
2677        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2678    
2679        include 'commontracker.f'        include 'commontracker.f'
2680          include 'level1.f'
2681        include 'common_momanhough.f'        include 'common_momanhough.f'
2682        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2683    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2684    
2685  *     output flag  *     output flag
2686  *     --------------  *     --------------
# Line 2803  c      common/dbg/DEBUG Line 2711  c      common/dbg/DEBUG
2711        distance=0        distance=0
2712        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2713        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2714          nloop=0                  
2715     91   continue                  
2716        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2717           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
2718  c     print*,'--------------'  c     print*,'--------------'
# Line 2904  c     print*,'check cp_used' Line 2814  c     print*,'check cp_used'
2814           do ip=1,nplanes           do ip=1,nplanes
2815              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2816           enddo           enddo
2817           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2818           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2819           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next doublet
2820                    
2821  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2822  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2823           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2824              if(DEBUG)print*,              if(verbose)print*,
2825       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2826       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2827       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2828  c     good2=.false.  c     good2=.false.
2829  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2830                do iv=1,nviews
2831                   mask_view(iv) = 6
2832                enddo
2833              iflag=1              iflag=1
2834              return              return
2835           endif           endif
# Line 2952  c$$$     $           ,(tr_cloud(iii),iii Line 2865  c$$$     $           ,(tr_cloud(iii),iii
2865  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2866  22288    continue  22288    continue
2867        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2868          
2869           if(nloop.lt.nstepx)then      
2870             cutdistxz=cutdistxz+cutxstep
2871             nloop=nloop+1          
2872             goto 91                
2873           endif                    
2874          
2875        if(DEBUG)then        if(DEBUG)then
2876           print*,'---------------------- '           print*,'---------------------- '
2877           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2979  c     02/02/2006 modified by Elena Vannu Line 2898  c     02/02/2006 modified by Elena Vannu
2898  c*****************************************************  c*****************************************************
2899    
2900        include 'commontracker.f'        include 'commontracker.f'
2901          include 'level1.f'
2902        include 'common_momanhough.f'        include 'common_momanhough.f'
2903        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2904        include 'common_mini_2.f'        include 'common_mini_2.f'
2905        include 'common_mech.f'        include 'common_mech.f'
2906        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2907    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2908    
2909  *     output flag  *     output flag
2910  *     --------------  *     --------------
# Line 3002  c      common/dbg/DEBUG Line 2920  c      common/dbg/DEBUG
2920  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2921  *     list of matching couples in the combination  *     list of matching couples in the combination
2922  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2923        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2924        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2925  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2926        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2927  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2928  *     variables for track fitting  *     variables for track fitting
2929        double precision AL_INI(5)        double precision AL_INI(5)
2930        double precision tath  c      double precision tath
2931  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2932  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2933    
# Line 3102  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 3020  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
3020  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
3021                            
3022  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
3023              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
3024              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
3025              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
3026       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
3027              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
3028              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
3029              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
3030                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
3031  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
3032              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
3033                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
3034  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
3035  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
3036                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
3037  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
3038  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
3039                c$$$            ELSE
3040    c$$$               AL_INI(4) = acos(-1.)/2
3041    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
3042    c$$$            ENDIF
3043    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
3044    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
3045    c$$$            
3046    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
3047    c$$$            
3048    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
3049                            
3050              if(DEBUG)then              if(DEBUG)then
3051                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
3052                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3186  c     $                                 Line 3114  c     $                                
3114  *     **********************************************************  *     **********************************************************
3115  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
3116  *     **********************************************************  *     **********************************************************
3117    cccc  scommentare se si usa al_ini della nuvola
3118    c$$$                              do i=1,5
3119    c$$$                                 AL(i)=AL_INI(i)
3120    c$$$                              enddo
3121                                  call guess()
3122                                do i=1,5                                do i=1,5
3123                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
3124                                enddo                                enddo
3125                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3126                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3127                                call mini_2(jstep,ifail)                                iprint=0
3128    c                              if(DEBUG)iprint=1
3129                                  if(DEBUG)iprint=1
3130                                  call mini2(jstep,ifail,iprint)
3131                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3132                                   if(DEBUG)then                                   if(DEBUG)then
3133                                      print *,                                      print *,
3134       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3135       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
3136                                        print*,'initial guess: '
3137    
3138                                        print*,'AL_INI(1) = ',AL_INI(1)
3139                                        print*,'AL_INI(2) = ',AL_INI(2)
3140                                        print*,'AL_INI(3) = ',AL_INI(3)
3141                                        print*,'AL_INI(4) = ',AL_INI(4)
3142                                        print*,'AL_INI(5) = ',AL_INI(5)
3143                                   endif                                   endif
3144                                   chi2=-chi2  c                                 chi2=-chi2
3145                                endif                                endif
3146  *     **********************************************************  *     **********************************************************
3147  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3211  c     $                                 Line 3154  c     $                                
3154  *     --------------------------  *     --------------------------
3155                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3156                                                                    
3157                                   if(DEBUG)print*,                                   if(verbose)print*,
3158       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3159       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3160       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
3161  c                                 good2=.false.  c                                 good2=.false.
3162  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
3163                                     do iv=1,nviews
3164                                        mask_view(iv) = 7
3165                                     enddo
3166                                   iflag=1                                   iflag=1
3167                                   return                                   return
3168                                endif                                endif
# Line 3315  cccccc 12/08/2006 modified by elena vann Line 3261  cccccc 12/08/2006 modified by elena vann
3261  c******************************************************  c******************************************************
3262    
3263        include 'commontracker.f'        include 'commontracker.f'
3264          include 'level1.f'
3265        include 'common_momanhough.f'        include 'common_momanhough.f'
3266        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3267        include 'common_mini_2.f'        include 'common_mini_2.f'
3268        include 'common_mech.f'        include 'common_mech.f'
3269        include 'momanhough_init.f'  c      include 'momanhough_init.f'
3270        include 'level1.f'  c      include 'level1.f'
3271        include 'calib.f'        include 'calib.f'
3272    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
3273    
3274  *     flag to chose PFA  *     flag to chose PFA
3275        character*10 PFA        character*10 PFA
# Line 3338  c      common/dbg/DEBUG Line 3283  c      common/dbg/DEBUG
3283        call track_init        call track_init
3284        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3285    
3286    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3287  *     -------------------------------------------------  *     -------------------------------------------------
3288  *     If the plane has been already included, it just  *     If the plane has been already included, it just
3289  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
3290  *     using improved PFAs  *     using improved PFAs
3291  *     -------------------------------------------------  *     -------------------------------------------------
3292    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3293           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3294       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3295                            
# Line 3377  c            dedxtrk(nplanes-ip+1) = (de Line 3324  c            dedxtrk(nplanes-ip+1) = (de
3324              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3325              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3326                            
3327    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3328  *     -------------------------------------------------  *     -------------------------------------------------
3329  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
3330  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
3331  *     -------------------------------------------------  *     -------------------------------------------------
3332    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3333           else                             else                  
3334                                
3335              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3435  c     $              'ETA2','ETA2', Line 3384  c     $              'ETA2','ETA2',
3384       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest))
3385                                
3386                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3387                   distance = distance / RCHI2_STORE(ibest)!<<< MS
3388                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3389                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
3390       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
# Line 3505  c     $              'ETA2','ETA2', Line 3455  c     $              'ETA2','ETA2',
3455       $              PFA,PFA,       $              PFA,PFA,
3456       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3457                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3458  c     if(DEBUG)print*,'normalized distance ',distance                 distance = distance / RCHI2_STORE(ibest)!<<< MS
3459                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3460       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
3461                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3533  c     $              'ETA2','ETA2', Line 3483  c     $              'ETA2','ETA2',
3483       $              PFA,PFA,       $              PFA,PFA,
3484       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest))
3485                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3486                   distance = distance / RCHI2_STORE(ibest)!<<< MS
3487                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3488       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
3489                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3572  c     $                 'ETA2','ETA2', Line 3523  c     $                 'ETA2','ETA2',
3523                 endif                 endif
3524    
3525                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3526                   distance = distance / RCHI2_STORE(ibest)!<<< MS
3527                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3528       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
3529                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3601  c                  dedxmm = dedx(icl)   Line 3553  c                  dedxmm = dedx(icl)  
3553                                
3554                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3555  *              ----------------------------  *              ----------------------------
3556    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3557                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3558                    XGOOD(nplanes-ip+1)=1.                    XGOOD(nplanes-ip+1)=1.
3559                    resx(nplanes-ip+1)=rxmm                    resx(nplanes-ip+1)=rxmm
3560                    if(DEBUG)print*,'%%%% included X-cl ',iclm                    if(DEBUG)print*,'%%%% included X-cl ',iclm
3561       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included X-cl ',iclm
3562         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
3563         $                 ,', norm.dist.= ',distmin
3564         $                 ,', cut ',clinc,' )'
3565                 else                 else
3566                    YGOOD(nplanes-ip+1)=1.                    YGOOD(nplanes-ip+1)=1.
3567                    resy(nplanes-ip+1)=rymm                    resy(nplanes-ip+1)=rymm
3568                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                    if(DEBUG)print*,'%%%% included Y-cl ',iclm
3569       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included Y-cl ',iclm
3570         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
3571         $                 ,', norm.dist.= ', distmin
3572         $                 ,', cut ',clinc,' )'
3573                 endif                 endif
3574    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3575  *              ----------------------------  *              ----------------------------
3576                 xm_A(nplanes-ip+1) = xmm_A                 xm_A(nplanes-ip+1) = xmm_A
3577                 ym_A(nplanes-ip+1) = ymm_A                 ym_A(nplanes-ip+1) = ymm_A
# Line 3645  cccccc 12/08/2006 modified by elena ---> Line 3605  cccccc 12/08/2006 modified by elena --->
3605        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3606    
3607        include 'commontracker.f'        include 'commontracker.f'
3608          include 'level1.f'
3609        include 'common_momanhough.f'        include 'common_momanhough.f'
3610        include 'momanhough_init.f'  c      include 'momanhough_init.f'
3611        include 'level2.f'        !(1)        include 'level2.f'        !(1)
3612  c      include 'calib.f'  c      include 'calib.f'
3613  c      include 'level1.f'  c      include 'level1.f'
3614    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
3615    
3616    
3617        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3828  c$$$ Line 3787  c$$$
3787    
3788        subroutine init_level2        subroutine init_level2
3789    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3790        include 'commontracker.f'        include 'commontracker.f'
3791          include 'level1.f'
3792        include 'common_momanhough.f'        include 'common_momanhough.f'
3793        include 'level2.f'        include 'level2.f'
3794        include 'level1.f'  c      include 'level1.f'
3795    
3796        do i=1,nviews        do i=1,nviews
3797           good2(i)=good1(i)           good2(i)=good1(i)
3798        enddo        enddo
3799    
 c      good2 = 0!.false.  
 c$$$      nev2 = nev1  
   
 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*****************************************************  
3800    
3801        NTRK = 0        NTRK = 0
3802        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3803           IMAGE(IT)=0           IMAGE(IT)=0
3804           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3805           do ip=1,nplanes           do ip=1,nplanes
3806              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3807              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
# Line 3877  c         BdL(IT) = 0. Line 3810  c         BdL(IT) = 0.
3810              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3811              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3812              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3813              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3814              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3815              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3816              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3817           enddo           enddo
# Line 3893  cccccc 17/8/2006 modified by elena Line 3822  cccccc 17/8/2006 modified by elena
3822              enddo                                enddo                  
3823           enddo                             enddo                  
3824        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3825        nclsx=0        nclsx=0
3826        nclsy=0              nclsy=0      
3827        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3828          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3829          xs(1,ip)=0          xs(1,ip)=0
3830          xs(2,ip)=0          xs(2,ip)=0
3831          sgnlxs(ip)=0          sgnlxs(ip)=0
3832          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3833          ys(1,ip)=0          ys(1,ip)=0
3834          ys(2,ip)=0          ys(2,ip)=0
3835          sgnlys(ip)=0          sgnlys(ip)=0
3836        enddo        enddo
 c*******************************************************  
3837        end        end
3838    
3839    
# Line 3926  c*************************************** Line 3848  c***************************************
3848  ************************************************************  ************************************************************
3849    
3850    
3851          subroutine init_hough
3852    
3853          include 'commontracker.f'
3854          include 'level1.f'
3855          include 'common_momanhough.f'
3856          include 'common_hough.f'
3857          include 'level2.f'
3858    
3859          ntrpt_nt=0
3860          ndblt_nt=0
3861          NCLOUDS_XZ_nt=0
3862          NCLOUDS_YZ_nt=0
3863          do idb=1,ndblt_max_nt
3864             db_cloud_nt(idb)=0
3865             alfayz1_nt(idb)=0      
3866             alfayz2_nt(idb)=0      
3867          enddo
3868          do itr=1,ntrpl_max_nt
3869             tr_cloud_nt(itr)=0
3870             alfaxz1_nt(itr)=0      
3871             alfaxz2_nt(itr)=0      
3872             alfaxz3_nt(itr)=0      
3873          enddo
3874          do idb=1,ncloyz_max      
3875            ptcloud_yz_nt(idb)=0    
3876            alfayz1_av_nt(idb)=0    
3877            alfayz2_av_nt(idb)=0    
3878          enddo                    
3879          do itr=1,ncloxz_max      
3880            ptcloud_xz_nt(itr)=0    
3881            alfaxz1_av_nt(itr)=0    
3882            alfaxz2_av_nt(itr)=0    
3883            alfaxz3_av_nt(itr)=0    
3884          enddo                    
3885    
3886          ntrpt=0                  
3887          ndblt=0                  
3888          NCLOUDS_XZ=0              
3889          NCLOUDS_YZ=0              
3890          do idb=1,ndblt_max        
3891            db_cloud(idb)=0        
3892            cpyz1(idb)=0            
3893            cpyz2(idb)=0            
3894            alfayz1(idb)=0          
3895            alfayz2(idb)=0          
3896          enddo                    
3897          do itr=1,ntrpl_max        
3898            tr_cloud(itr)=0        
3899            cpxz1(itr)=0            
3900            cpxz2(itr)=0            
3901            cpxz3(itr)=0            
3902            alfaxz1(itr)=0          
3903            alfaxz2(itr)=0          
3904            alfaxz3(itr)=0          
3905          enddo                    
3906          do idb=1,ncloyz_max      
3907            ptcloud_yz(idb)=0      
3908            alfayz1_av(idb)=0      
3909            alfayz2_av(idb)=0      
3910            do idbb=1,ncouplemaxtot
3911              cpcloud_yz(idb,idbb)=0
3912            enddo                  
3913          enddo                    
3914          do itr=1,ncloxz_max      
3915            ptcloud_xz(itr)=0      
3916            alfaxz1_av(itr)=0      
3917            alfaxz2_av(itr)=0      
3918            alfaxz3_av(itr)=0      
3919            do itrr=1,ncouplemaxtot
3920              cpcloud_xz(itr,itrr)=0
3921            enddo                  
3922          enddo                    
3923          end
3924    ************************************************************
3925    *
3926    *
3927    *
3928    *
3929    *
3930    *
3931    *
3932    ************************************************************
3933    
3934    
3935        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3936    
3937  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3936  c*************************************** Line 3942  c***************************************
3942    
3943            
3944        include 'commontracker.f'        include 'commontracker.f'
3945    c      include 'level1.f'
3946        include 'level1.f'        include 'level1.f'
3947          include 'common_momanhough.f'
3948        include 'level2.f'        include 'level2.f'
3949        include 'common_mini_2.f'        include 'common_mini_2.f'
3950        include 'common_momanhough.f'        real sinth,phi,pig      
       real sinth,phi,pig        !(4)  
3951        pig=acos(-1.)        pig=acos(-1.)
3952    
 c      good2=1!.true.  
3953        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3954        nstep_nt(ntr)       = 0!nstep        nstep_nt(ntr)       = nstep
3955    
3956          phi   = al(4)          
3957          sinth = al(3)            
3958          if(sinth.lt.0)then      
3959             sinth = -sinth        
3960             phi = phi + pig      
3961          endif                    
3962          npig = aint(phi/(2*pig))
3963          phi = phi - npig*2*pig  
3964          if(phi.lt.0)            
3965         $     phi = phi + 2*pig  
3966          al(4) = phi              
3967          al(3) = sinth            
3968    
       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)  
 *****************************************************  
3969        do i=1,5        do i=1,5
3970           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3971           do j=1,5           do j=1,5
3972              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3973           enddo           enddo
 c     print*,al_nt(i,ntr)  
3974        enddo        enddo
3975                
3976        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
# Line 3981  c     print*,al_nt(i,ntr) Line 3986  c     print*,al_nt(i,ntr)
3986           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3987           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3988           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))
 c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3989           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
3990           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)  
3991        
# Line 3998  c            print*,ip,' ',cltrx(ip,ntr) Line 4002  c            print*,ip,' ',cltrx(ip,ntr)
4002           endif                     endif          
4003    
4004        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)  
4005    
4006    
4007        end        end
4008    
4009        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*****************************************************  
4010    
4011  *     -------------------------------------------------------  *     -------------------------------------------------------
4012  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 4022  c*************************************** Line 4015  c***************************************
4015  *     -------------------------------------------------------  *     -------------------------------------------------------
4016    
4017        include 'commontracker.f'        include 'commontracker.f'
4018        include 'level1.f'  c      include 'level1.f'
       include 'level2.f'  
4019        include 'calib.f'        include 'calib.f'
4020          include 'level1.f'
4021        include 'common_momanhough.f'        include 'common_momanhough.f'
4022          include 'level2.f'
4023        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
4024    
4025  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
# Line 4033  c      good2=1!.true. Line 4027  c      good2=1!.true.
4027        nclsx = 0        nclsx = 0
4028        nclsy = 0        nclsy = 0
4029    
4030          do iv = 1,nviews
4031             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
4032          enddo
4033    
4034        do icl=1,nclstr1        do icl=1,nclstr1
4035           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
4036              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
# Line 4076  c      print*,icl,cl_used(icl),cl_good(i Line 4074  c      print*,icl,cl_used(icl),cl_good(i
4074        enddo        enddo
4075        end        end
4076    
4077    ***************************************************
4078    *                                                 *
4079    *                                                 *
4080    *                                                 *
4081    *                                                 *
4082    *                                                 *
4083    *                                                 *
4084    **************************************************
4085    
4086          subroutine fill_hough
4087    
4088    *     -------------------------------------------------------
4089    *     This routine fills the  variables related to the hough
4090    *     transform, for the debig n-tuple
4091    *     -------------------------------------------------------
4092    
4093          include 'commontracker.f'
4094          include 'level1.f'
4095          include 'common_momanhough.f'
4096          include 'common_hough.f'
4097          include 'level2.f'
4098    
4099          if(.false.
4100         $     .or.ntrpt.gt.ntrpt_max_nt
4101         $     .or.ndblt.gt.ndblt_max_nt
4102         $     .or.NCLOUDS_XZ.gt.ncloxz_max
4103         $     .or.NCLOUDS_yZ.gt.ncloyz_max
4104         $     )then
4105             ntrpt_nt=0
4106             ndblt_nt=0
4107             NCLOUDS_XZ_nt=0
4108             NCLOUDS_YZ_nt=0
4109          else
4110             ndblt_nt=ndblt
4111             ntrpt_nt=ntrpt
4112             if(ndblt.ne.0)then
4113                do id=1,ndblt
4114                   alfayz1_nt(id)=alfayz1(id) !Y0
4115                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
4116                enddo
4117             endif
4118             if(ndblt.ne.0)then
4119                do it=1,ntrpt
4120                   alfaxz1_nt(it)=alfaxz1(it) !X0
4121                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
4122                   alfaxz3_nt(it)=alfaxz3(it) !1/r
4123                enddo
4124             endif
4125             nclouds_yz_nt=nclouds_yz
4126             nclouds_xz_nt=nclouds_xz
4127             if(nclouds_yz.ne.0)then
4128                nnn=0
4129                do iyz=1,nclouds_yz
4130                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
4131                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
4132                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4133                   nnn=nnn+ptcloud_yz(iyz)
4134                enddo
4135                do ipt=1,nnn
4136                   db_cloud_nt(ipt)=db_cloud(ipt)
4137                 enddo
4138             endif
4139             if(nclouds_xz.ne.0)then
4140                nnn=0
4141                do ixz=1,nclouds_xz
4142                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4143                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4144                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4145                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4146                   nnn=nnn+ptcloud_xz(ixz)              
4147                enddo
4148                do ipt=1,nnn
4149                  tr_cloud_nt(ipt)=tr_cloud(ipt)
4150                 enddo
4151             endif
4152          endif
4153          end
4154          

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.13

  ViewVC Help
Powered by ViewVC 1.1.23