/[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.9 by pam-fi, Fri Oct 27 16:08:19 2006 UTC revision 1.13 by pam-fi, Wed Nov 8 16:42:28 2006 UTC
# Line 183  c      iflag=0 Line 183  c      iflag=0
183  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
184                
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 281  c$$$         if(ibest.eq.0)goto 880 !>> Line 282  c$$$         if(ibest.eq.0)goto 880 !>>
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           iprint=0           iprint=0
299           if(DEBUG)iprint=1  c         if(DEBUG)iprint=1
300             if(VERBOSE)iprint=1
301           call mini2(jstep,ifail,iprint)           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 *** (mini2) '       $              '*** 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 660  c      include 'level1.f' Line 675  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 1756  c      include 'level1.f' Line 1771  c      include 'level1.f'
1771        do iv=1,nviews        do iv=1,nviews
1772           if( ncl_view(iv).gt. nclustermax)then           if( ncl_view(iv).gt. nclustermax)then
1773              mask_view(iv) = 1              mask_view(iv) = 1
1774              print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1775       $           ,nclustermax,' on view ', iv,' --> masked!'       $           ,nclustermax,' on view ', iv,' --> masked!'
1776           endif           endif
1777        enddo        enddo
# Line 1913  c$$$               endif Line 1928  c$$$               endif
1928       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1929       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1930       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1931       $                 ,'( ',ncouplemax,' ) NB - THIS SHOULD NOT HAPPEN'       $                 ,'( ',ncouplemax,' )'
1932  c     good2=.false.  c     good2=.false.
1933  c     goto 880   !fill ntp and go to next event                      c     goto 880   !fill ntp and go to next event
1934                    iflag=1                    mask_view(nviewx(nplx)) = 2
1935                    return                    mask_view(nviewy(nply)) = 2
1936    c                  iflag=1
1937    c                  return
1938                 endif                 endif
1939                                
1940                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
# Line 1952  c     goto 880   !fill ntp and go to nex Line 1969  c     goto 880   !fill ntp and go to nex
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(verbose)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 1978  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)
1997    c$$$
1998        include 'commontracker.f'  c$$$      include 'commontracker.f'
1999        include 'level1.f'  c$$$      include 'level1.f'
2000        include 'common_momanhough.f'  c$$$      include 'common_momanhough.f'
2001  c      include 'momanhough_init.f'  c$$$c      include 'momanhough_init.f'
2002        include 'calib.f'  c$$$      include 'calib.f'
2003  c      include 'level1.f'  c$$$c      include 'level1.f'
2004    c$$$
2005    c$$$
2006  *     output flag  c$$$*     output flag
2007  *     --------------  c$$$*     --------------
2008  *     0 = good event  c$$$*     0 = good event
2009  *     1 = bad event  c$$$*     1 = bad event
2010  *     --------------  c$$$*     --------------
2011        integer iflag  c$$$      integer iflag
2012    c$$$
2013        integer badseed,badcl  c$$$      integer badseed,badcl
2014    c$$$
2015  *     init variables  c$$$*     init variables
2016        ncp_tot=0  c$$$      ncp_tot=0
2017        do ip=1,nplanes  c$$$      do ip=1,nplanes
2018           do ico=1,ncouplemax  c$$$         do ico=1,ncouplemax
2019              clx(ip,ico)=0  c$$$            clx(ip,ico)=0
2020              cly(ip,ico)=0  c$$$            cly(ip,ico)=0
2021           enddo  c$$$         enddo
2022           ncp_plane(ip)=0  c$$$         ncp_plane(ip)=0
2023           do icl=1,nclstrmax_level2  c$$$         do icl=1,nclstrmax_level2
2024              cls(ip,icl)=1  c$$$            cls(ip,icl)=1
2025           enddo  c$$$         enddo
2026           ncls(ip)=0  c$$$         ncls(ip)=0
2027        enddo  c$$$      enddo
2028        do icl=1,nclstrmax_level2  c$$$      do icl=1,nclstrmax_level2
2029           cl_single(icl)=1  c$$$         cl_single(icl)=1
2030           cl_good(icl)=0  c$$$         cl_good(icl)=0
2031        enddo  c$$$      enddo
2032          c$$$      
2033  *     start association  c$$$*     start association
2034        ncouples=0  c$$$      ncouples=0
2035        do icx=1,nclstr1          !loop on cluster (X)  c$$$      do icx=1,nclstr1          !loop on cluster (X)
2036           if(mod(VIEW(icx),2).eq.1)goto 10  c$$$         if(mod(VIEW(icx),2).eq.1)goto 10
2037            c$$$        
2038  *     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2039  *     cut on charge (X VIEW)  c$$$*     cut on charge (X VIEW)
2040           if(dedx(icx).lt.dedx_x_min)then  c$$$         if(dedx(icx).lt.dedx_x_min)then
2041              cl_single(icx)=0  c$$$            cl_single(icx)=0
2042              goto 10  c$$$            goto 10
2043           endif  c$$$         endif
2044  *     cut BAD (X VIEW)              c$$$*     cut BAD (X VIEW)            
2045           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
2046           ifirst=INDSTART(icx)  c$$$         ifirst=INDSTART(icx)
2047           if(icx.ne.nclstr1) then  c$$$         if(icx.ne.nclstr1) then
2048              ilast=INDSTART(icx+1)-1  c$$$            ilast=INDSTART(icx+1)-1
2049           else  c$$$         else
2050              ilast=TOTCLLENGTH  c$$$            ilast=TOTCLLENGTH
2051           endif  c$$$         endif
2052           badcl=badseed  c$$$         badcl=badseed
2053           do igood=-ngoodstr,ngoodstr  c$$$         do igood=-ngoodstr,ngoodstr
2054              ibad=1  c$$$            ibad=1
2055              if((INDMAX(icx)+igood).gt.ifirst.and.  c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.
2056       $           (INDMAX(icx)+igood).lt.ilast.and.  c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.
2057       $           .true.)then  c$$$     $           .true.)then
2058                 ibad=BAD(VIEW(icx),  c$$$               ibad=BAD(VIEW(icx),
2059       $              nvk(MAXS(icx)+igood),  c$$$     $              nvk(MAXS(icx)+igood),
2060       $              nst(MAXS(icx)+igood))  c$$$     $              nst(MAXS(icx)+igood))
2061              endif  c$$$            endif
2062              badcl=badcl*ibad  c$$$            badcl=badcl*ibad
2063           enddo  c$$$         enddo
2064           if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  c$$$         if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut
2065              cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  c$$$            cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut
2066              goto 10             !<<<<<<<<<<<<<< BAD cut  c$$$            goto 10             !<<<<<<<<<<<<<< BAD cut
2067           endif                  !<<<<<<<<<<<<<< BAD cut  c$$$         endif                  !<<<<<<<<<<<<<< BAD cut
2068  *     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2069            c$$$        
2070           cl_good(icx)=1  c$$$         cl_good(icx)=1
2071           nplx=npl(VIEW(icx))  c$$$         nplx=npl(VIEW(icx))
2072           nldx=nld(MAXS(icx),VIEW(icx))  c$$$         nldx=nld(MAXS(icx),VIEW(icx))
2073            c$$$        
2074           do icy=1,nclstr1       !loop on cluster (Y)  c$$$         do icy=1,nclstr1       !loop on cluster (Y)
2075              if(mod(VIEW(icy),2).eq.0)goto 20  c$$$            if(mod(VIEW(icy),2).eq.0)goto 20
2076                c$$$            
2077  *     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2078  *     cut on charge (Y VIEW)  c$$$*     cut on charge (Y VIEW)
2079              if(dedx(icy).lt.dedx_y_min)then  c$$$            if(dedx(icy).lt.dedx_y_min)then
2080                 cl_single(icy)=0  c$$$               cl_single(icy)=0
2081                 goto 20  c$$$               goto 20
2082              endif  c$$$            endif
2083  *     cut BAD (Y VIEW)              c$$$*     cut BAD (Y VIEW)            
2084              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
2085              ifirst=INDSTART(icy)  c$$$            ifirst=INDSTART(icy)
2086              if(icy.ne.nclstr1) then  c$$$            if(icy.ne.nclstr1) then
2087                 ilast=INDSTART(icy+1)-1  c$$$               ilast=INDSTART(icy+1)-1
2088              else  c$$$            else
2089                 ilast=TOTCLLENGTH  c$$$               ilast=TOTCLLENGTH
2090              endif  c$$$            endif
2091              badcl=badseed  c$$$            badcl=badseed
2092              do igood=-ngoodstr,ngoodstr  c$$$            do igood=-ngoodstr,ngoodstr
2093                 ibad=1  c$$$               ibad=1
2094                 if((INDMAX(icy)+igood).gt.ifirst.and.  c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.
2095       $              (INDMAX(icy)+igood).lt.ilast.and.  c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.
2096       $              .true.)  c$$$     $              .true.)
2097       $              ibad=BAD(VIEW(icy),  c$$$     $              ibad=BAD(VIEW(icy),
2098       $              nvk(MAXS(icy)+igood),  c$$$     $              nvk(MAXS(icy)+igood),
2099       $              nst(MAXS(icy)+igood))  c$$$     $              nst(MAXS(icy)+igood))
2100                 badcl=badcl*ibad  c$$$               badcl=badcl*ibad
2101              enddo  c$$$            enddo
2102              if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  c$$$            if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut
2103                 cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  c$$$               cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut
2104                 goto 20          !<<<<<<<<<<<<<< BAD cut  c$$$               goto 20          !<<<<<<<<<<<<<< BAD cut
2105              endif               !<<<<<<<<<<<<<< BAD cut  c$$$            endif               !<<<<<<<<<<<<<< BAD cut
2106  *     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2107                c$$$            
2108                c$$$            
2109              cl_good(icy)=1                    c$$$            cl_good(icy)=1                  
2110              nply=npl(VIEW(icy))  c$$$            nply=npl(VIEW(icy))
2111              nldy=nld(MAXS(icy),VIEW(icy))  c$$$            nldy=nld(MAXS(icy),VIEW(icy))
2112                c$$$            
2113  *     ----------------------------------------------  c$$$*     ----------------------------------------------
2114  *     CONDITION TO FORM A COUPLE  c$$$*     CONDITION TO FORM A COUPLE
2115  *     ----------------------------------------------  c$$$*     ----------------------------------------------
2116  *     geometrical consistency (same plane and ladder)  c$$$*     geometrical consistency (same plane and ladder)
2117              if(nply.eq.nplx.and.nldy.eq.nldx)then  c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then
2118  *     charge correlation  c$$$*     charge correlation
2119  *     ===========================================================  c$$$*     ===========================================================
2120  *     this version of the subroutine is used for the calibration  c$$$*     this version of the subroutine is used for the calibration
2121  *     thus charge-correlation selection is obviously removed  c$$$*     thus charge-correlation selection is obviously removed
2122  *     ===========================================================  c$$$*     ===========================================================
2123  c$$$               ddd=(dedx(icy)  c$$$c$$$               ddd=(dedx(icy)
2124  c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
2125  c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
2126  c$$$               cut=chcut*sch(nplx,nldx)  c$$$c$$$               cut=chcut*sch(nplx,nldx)
2127  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$$$              
2131  *     ------------------> COUPLE <------------------  c$$$*     ------------------> COUPLE <------------------
2132  *     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$$$c     good2=.false.
2140    c$$$c$$$c     goto 880   !fill ntp and go to next event
2141    c$$$c$$$                  iflag=1
2142    c$$$c$$$                  return
2143    c$$$c$$$               endif
2144    c$$$              
2145    c$$$               if(ncp_plane(nplx).eq.ncouplemax)then
2146    c$$$                  if(verbose)print*,
2147    c$$$     $                 '** warning ** number of identified '//
2148    c$$$     $                 'couples on plane ',nplx,
2149    c$$$     $                 'exceeds vector dimention '
2150    c$$$     $                 ,'( ',ncouplemax,' )'
2151  c$$$c     good2=.false.  c$$$c     good2=.false.
2152  c$$$c     goto 880   !fill ntp and go to next event  c$$$c     goto 880   !fill ntp and go to next event                    
2153  c$$$                  iflag=1  c$$$                  iflag=1
2154  c$$$                  return  c$$$                  return
2155  c$$$               endif  c$$$               endif
2156                  c$$$              
2157                 if(ncp_plane(nplx).eq.ncouplemax)then  c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1
2158                    if(verbose)print*,  c$$$               clx(nplx,ncp_plane(nplx))=icx
2159       $                 '** warning ** number of identified '//  c$$$               cly(nply,ncp_plane(nplx))=icy
2160       $                 'couples on plane ',nplx,  c$$$               cl_single(icx)=0
2161       $                 'exceeds vector dimention '  c$$$               cl_single(icy)=0
2162       $                 ,'( ',ncouplemax,' )'  c$$$            endif                              
2163  c     good2=.false.  c$$$*     ----------------------------------------------
2164  c     goto 880   !fill ntp and go to next event                      c$$$
2165                    iflag=1  c$$$ 20         continue
2166                    return  c$$$         enddo                  !end loop on clusters(Y)
2167                 endif  c$$$        
2168                  c$$$ 10      continue
2169                 ncp_plane(nplx) = ncp_plane(nplx) + 1  c$$$      enddo                     !end loop on clusters(X)
2170                 clx(nplx,ncp_plane(nplx))=icx  c$$$      
2171                 cly(nply,ncp_plane(nplx))=icy  c$$$      
2172                 cl_single(icx)=0  c$$$      do icl=1,nclstr1
2173                 cl_single(icy)=0  c$$$         if(cl_single(icl).eq.1)then
2174              endif                                c$$$            ip=npl(VIEW(icl))
2175  *     ----------------------------------------------  c$$$            ncls(ip)=ncls(ip)+1
2176    c$$$            cls(ip,ncls(ip))=icl
2177   20         continue  c$$$         endif
2178           enddo                  !end loop on clusters(Y)  c$$$      enddo
2179            c$$$      
2180   10      continue  c$$$      
2181        enddo                     !end loop on clusters(X)  c$$$      if(DEBUG)then
2182          c$$$         print*,'clusters  ',nclstr1
2183          c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)
2184        do icl=1,nclstr1  c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)
2185           if(cl_single(icl).eq.1)then  c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
2186              ip=npl(VIEW(icl))  c$$$      endif
2187              ncls(ip)=ncls(ip)+1  c$$$      
2188              cls(ip,ncls(ip))=icl  c$$$      do ip=1,6
2189           endif  c$$$         ncp_tot=ncp_tot+ncp_plane(ip)
2190        enddo  c$$$      enddo
2191          c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
2192          c$$$      
2193        if(DEBUG)then  c$$$      if(ncp_tot.gt.ncp_max)then
2194           print*,'clusters  ',nclstr1  c$$$         if(verbose)print*,
2195           print*,'good    ',(cl_good(i),i=1,nclstr1)  c$$$     $           '** warning ** number of identified '//
2196           print*,'singles ',(cl_single(i),i=1,nclstr1)  c$$$     $           'couples exceeds upper limit for Hough tr. '
2197           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  c$$$     $           ,'( ',ncp_max,' )'            
2198        endif  c$$$c            good2=.false.
2199          c$$$c     goto 880       !fill ntp and go to next event
2200        do ip=1,6  c$$$         iflag=1
2201           ncp_tot=ncp_tot+ncp_plane(ip)  c$$$         return
2202        enddo  c$$$      endif
2203  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c$$$      
2204          c$$$      return
2205        if(ncp_tot.gt.ncp_max)then  c$$$      end
2206           if(verbose)print*,  c$$$
      $           '** 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  
   
2207                
2208  ***************************************************  ***************************************************
2209  *                                                 *  *                                                 *
# Line 2286  c     $                       (icx2,icy2 Line 2301  c     $                       (icx2,icy2
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 2356  c     $                                 Line 2374  c     $                                
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 2584  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2605  c         if(ncpused.lt.ncpyz_min)goto 2
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 2803  c         if(ncpused.lt.ncpxz_min)goto 2 Line 2827  c         if(ncpused.lt.ncpxz_min)goto 2
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 2900  c      include 'momanhough_init.f' Line 2927  c      include 'momanhough_init.f'
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 2993  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) = PIGR + 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 3077  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                                iprint=0                                iprint=0
3128    c                              if(DEBUG)iprint=1
3129                                if(DEBUG)iprint=1                                if(DEBUG)iprint=1
3130                                call mini2(jstep,ifail,iprint)                                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       $                              //'(mini2 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 3110  c     $                                 Line 3160  c     $                                
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 3230  c      include 'level1.f' Line 3283  c      include 'level1.f'
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 3269  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 3327  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 3397  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 3425  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 3464  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 3493  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
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.23