/[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.3 by pam-fi, Tue Sep 5 12:52:20 2006 UTC revision 1.8 by pam-fi, Wed Oct 25 16:18:41 2006 UTC
# Line 242  c         iflag=0 Line 242  c         iflag=0
242           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
243           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
244           jstep=0                !# minimization steps           jstep=0                !# minimization steps
245             iprint=0
246           call mini_2(jstep,ifail)           if(DEBUG)iprint=1
247             call mini2(jstep,ifail,iprint)
248           if(ifail.ne.0) then           if(ifail.ne.0) then
249              if(DEBUG)then              if(DEBUG)then
250                 print *,                 print *,
251       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (mini2) '
252       $              ,iev       $              ,iev
253              endif              endif
254              chi2=-chi2              chi2=-chi2
# Line 311  c         print*,'++++++++++ iimage,fima Line 312  c         print*,'++++++++++ iimage,fima
312  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
313    
314           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
315              if(DEBUG)              if(verbose)
316       $           print*,       $           print*,
317       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
318       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 667  c      double precision xi_B,yi_B,zi_B Line 668  c      double precision xi_B,yi_B,zi_B
668              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
669           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
670  c            cog2 = cog(2,icx)  c            cog2 = cog(2,icx)
671  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)              c            etacorr = pfaeta2(cog2,viewx,nldx,angx)            
672  c            stripx = stripx + etacorr  c            stripx = stripx + etacorr
673              stripx = stripx + pfa_eta2(icx,angx)            !(3)              stripx = stripx + pfaeta2(icx,angx)            !(3)
674              resxPAM = risx_eta2(angx)                       !   (4)              resxPAM = risx_eta2(angx)                       !   (4)
675              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
676       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
677              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
678           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                         !(3)
679              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)            !(3)
680              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(angx)                       !   (4)
681              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
682       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
683              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
684           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                         !(3)
685              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            !(3)
686              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(angx)                       !   (4)
687              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
688       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
689              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
690           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then                          !(3)
691              stripx = stripx + pfa_eta(icx,angx)             !(3)              stripx = stripx + pfaeta(icx,angx)             !(3)
692              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = ris_eta(icx,angx)                     !   (4)
693              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
694       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
# Line 702  c            resxPAM = resxPAM*fbad_cog( Line 703  c            resxPAM = resxPAM*fbad_cog(
703           endif           endif
704    
705        endif        endif
706    c      if(icy.eq.0.and.icx.ne.0)
707    c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'
708                
709  *     -----------------  *     -----------------
710  *     CLUSTER Y  *     CLUSTER Y
# Line 729  c            resxPAM = resxPAM*fbad_cog( Line 732  c            resxPAM = resxPAM*fbad_cog(
732              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
733           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
734  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
735  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
736  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
737              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
738              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(angy)                       !   (4)
739              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
740              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
741       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
742           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
743              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
744              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
745              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
746       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
747           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
748              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
749              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
750              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
751       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
752           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
753              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
754              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
755  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
756              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 773  C======================================= Line 776  C=======================================
776  c------------------------------------------------------------------------  c------------------------------------------------------------------------
777  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
778  c------------------------------------------------------------------------  c------------------------------------------------------------------------
779             if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
780         $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
781                print*,'xyz_PAM (couple):',
782         $          ' WARNING: false X strip: strip ',stripx
783             endif
784           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
785           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
786           zi = 0.           zi = 0.
# Line 859  C======================================= Line 867  C=======================================
867              nldy = nldx              nldy = nldx
868              viewy = viewx - 1              viewy = viewx - 1
869    
870    c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
871    c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
872                if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
873         $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
874                   print*,'xyz_PAM (X-singlet):',
875         $             ' WARNING: false X strip: strip ',stripx
876                endif
877              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
878    
879              xi_A = xi              xi_A = xi
# Line 1137  c$$$         print*,' resolution ',resxP Line 1152  c$$$         print*,' resolution ',resxP
1152  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1153  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1154  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1155                   if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1156         $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1157    c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1158                      print*,'whichsensor: ',
1159         $                ' WARNING: false X strip: strip ',stripx
1160                   endif
1161                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1162                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
1163                 zi = 0.                 zi = 0.
# Line 1655  c      common/dbg/DEBUG Line 1676  c      common/dbg/DEBUG
1676  *     --------------  *     --------------
1677        integer iflag        integer iflag
1678    
1679        integer badseed,badcl        integer badseed,badclx,badcly
1680    
1681  *     init variables  *     init variables
1682        ncp_tot=0        ncp_tot=0
# Line 1682  c      common/dbg/DEBUG Line 1703  c      common/dbg/DEBUG
1703                    
1704  *     ----------------------------------------------------  *     ----------------------------------------------------
1705  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1706    *     ----------------------------------------------------
1707           if(dedx(icx).lt.dedx_x_min)then           if(dedx(icx).lt.dedx_x_min)then
1708              cl_single(icx)=0              cl_single(icx)=0
1709              goto 10              goto 10
1710           endif           endif
1711    *     ----------------------------------------------------
1712    *     cut on multiplicity (X VIEW)
1713    *     ----------------------------------------------------
1714             if(mult(icx).ge.mult_x_max)then
1715                cl_single(icx)=0
1716                goto 10
1717             endif
1718    *     ----------------------------------------------------
1719  *     cut BAD (X VIEW)              *     cut BAD (X VIEW)            
1720    *     ----------------------------------------------------
1721           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1722           ifirst=INDSTART(icx)           ifirst=INDSTART(icx)
1723           if(icx.ne.nclstr1) then           if(icx.ne.nclstr1) then
# Line 1694  c      common/dbg/DEBUG Line 1725  c      common/dbg/DEBUG
1725           else           else
1726              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1727           endif           endif
1728           badcl=badseed           badclx=badseed
1729           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1730              ibad=1              ibad=1
1731              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1704  c      common/dbg/DEBUG Line 1735  c      common/dbg/DEBUG
1735       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1736       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1737              endif              endif
1738              badcl=badcl*ibad              badclx=badclx*ibad
1739           enddo           enddo
1740    *     ----------------------------------------------------
1741    *     >>> eliminato il taglio sulle BAD <<<
1742    *     ----------------------------------------------------
1743  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1744  c     cl_single(icx)=0  c     cl_single(icx)=0
1745  c     goto 10  c     goto 10
# Line 1721  c     endif Line 1755  c     endif
1755                            
1756  *     ----------------------------------------------------  *     ----------------------------------------------------
1757  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1758    *     ----------------------------------------------------
1759              if(dedx(icy).lt.dedx_y_min)then              if(dedx(icy).lt.dedx_y_min)then
1760                 cl_single(icy)=0                 cl_single(icy)=0
1761                 goto 20                 goto 20
1762              endif              endif
1763    *     ----------------------------------------------------
1764    *     cut on multiplicity (X VIEW)
1765    *     ----------------------------------------------------
1766                if(mult(icy).ge.mult_y_max)then
1767                   cl_single(icy)=0
1768                   goto 20
1769                endif
1770    *     ----------------------------------------------------
1771  *     cut BAD (Y VIEW)              *     cut BAD (Y VIEW)            
1772    *     ----------------------------------------------------
1773              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1774              ifirst=INDSTART(icy)              ifirst=INDSTART(icy)
1775              if(icy.ne.nclstr1) then              if(icy.ne.nclstr1) then
# Line 1733  c     endif Line 1777  c     endif
1777              else              else
1778                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1779              endif              endif
1780              badcl=badseed              badcly=badseed
1781              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1782                 ibad=1                 ibad=1
1783                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1742  c     endif Line 1786  c     endif
1786       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1787       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1788       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1789                 badcl=badcl*ibad                 badcly=badcly*ibad
1790              enddo              enddo
1791    *     ----------------------------------------------------
1792    *     >>> eliminato il taglio sulle BAD <<<
1793    *     ----------------------------------------------------
1794  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1795  c     cl_single(icy)=0  c     cl_single(icy)=0
1796  c     goto 20  c     goto 20
1797  c     endif  c     endif
1798  *     ----------------------------------------------------  *     ----------------------------------------------------
1799                            
               
1800              cl_good(icy)=1                                cl_good(icy)=1                  
1801              nply=npl(VIEW(icy))              nply=npl(VIEW(icy))
1802              nldy=nld(MAXS(icy),VIEW(icy))              nldy=nld(MAXS(icy),VIEW(icy))
# Line 1761  c     endif Line 1807  c     endif
1807  *     geometrical consistency (same plane and ladder)  *     geometrical consistency (same plane and ladder)
1808              if(nply.eq.nplx.and.nldy.eq.nldx)then              if(nply.eq.nplx.and.nldy.eq.nldx)then
1809  *     charge correlation  *     charge correlation
1810                 ddd=(dedx(icy)  *     (modified to be applied only below saturation... obviously)
1811       $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
1812                 ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  *     -------------------------------------------------------------
1813                 cut=chcut*sch(nplx,nldx)  *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<
1814                 if(abs(ddd).gt.cut)goto 20 !charge not consistent  *     -------------------------------------------------------------
1815                                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1816         $              .and.
1817         $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1818         $              .and.
1819         $              (badclx.eq.1.and.badcly.eq.1)
1820         $              .and.
1821         $              .true.)then
1822    
1823                      ddd=(dedx(icy)
1824         $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1825                      ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1826    
1827    c                  cut = chcut * sch(nplx,nldx)
1828    
1829                      sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1830         $                 -kch(nplx,nldx)*cch(nplx,nldx))
1831                      sss=sss/sqrt(kch(nplx,nldx)**2+1)
1832                      cut = chcut * (16 + sss/50.)
1833    
1834                      if(abs(ddd).gt.cut)then
1835                         goto 20    !charge not consistent
1836                      endif
1837                   endif
1838                                
1839  *     ------------------> COUPLE <------------------  *     ------------------> COUPLE <------------------
1840  *     check to do not overflow vector dimentions  *     check to do not overflow vector dimentions
1841                 if(ncp_plane(nplx).gt.ncouplemax)then  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
1842                    if(DEBUG)print*,  c$$$                  if(DEBUG)print*,
1843       $                    ' ** warning ** number of identified'//  c$$$     $                    ' ** warning ** number of identified'//
1844       $                    ' couples on plane ',nplx,  c$$$     $                    ' couples on plane ',nplx,
1845       $                    ' exceeds vector dimention'//  c$$$     $                    ' exceeds vector dimention'//
1846       $                    ' ( ',ncouplemax,' )'  c$$$     $                    ' ( ',ncouplemax,' )'
1847  c     good2=.false.  c$$$c     good2=.false.
1848  c     goto 880   !fill ntp and go to next event  c$$$c     goto 880   !fill ntp and go to next event
1849                    iflag=1  c$$$                  iflag=1
1850                    return  c$$$                  return
1851                 endif  c$$$               endif
1852                                
1853                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).eq.ncouplemax)then
1854                    if(DEBUG)print*,                    if(verbose)print*,
1855       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1856       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1857       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
# Line 1831  c     goto 880   !fill ntp and go to nex Line 1899  c     goto 880   !fill ntp and go to nex
1899  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
1900                
1901        if(ncp_tot.gt.ncp_max)then        if(ncp_tot.gt.ncp_max)then
1902           if(DEBUG)print*,           if(verbose)print*,
1903       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
1904       $           'couples exceeds upper limit for Hough tr. '       $           'couples exceeds upper limit for Hough tr. '
1905       $           ,'( ',ncp_max,' )'                   $           ,'( ',ncp_max,' )'            
# Line 1990  c$$$               if(abs(ddd).gt.cut)go Line 2058  c$$$               if(abs(ddd).gt.cut)go
2058                                
2059  *     ------------------> COUPLE <------------------  *     ------------------> COUPLE <------------------
2060  *     check to do not overflow vector dimentions  *     check to do not overflow vector dimentions
2061                 if(ncp_plane(nplx).gt.ncouplemax)then  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
2062                    if(DEBUG)print*,  c$$$                  if(DEBUG)print*,
2063       $                    ' ** warning ** number of identified'//  c$$$     $                    ' ** warning ** number of identified'//
2064       $                    ' couples on plane ',nplx,  c$$$     $                    ' couples on plane ',nplx,
2065       $                    ' exceeds vector dimention'//  c$$$     $                    ' exceeds vector dimention'//
2066       $                    ' ( ',ncouplemax,' )'  c$$$     $                    ' ( ',ncouplemax,' )'
2067  c     good2=.false.  c$$$c     good2=.false.
2068  c     goto 880   !fill ntp and go to next event  c$$$c     goto 880   !fill ntp and go to next event
2069                    iflag=1  c$$$                  iflag=1
2070                    return  c$$$                  return
2071                 endif  c$$$               endif
2072                                
2073                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).eq.ncouplemax)then
2074                    if(DEBUG)print*,                    if(verbose)print*,
2075       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
2076       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
2077       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
# Line 2051  c     goto 880   !fill ntp and go to nex Line 2119  c     goto 880   !fill ntp and go to nex
2119  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
2120                
2121        if(ncp_tot.gt.ncp_max)then        if(ncp_tot.gt.ncp_max)then
2122           if(DEBUG)print*,           if(verbose)print*,
2123       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2124       $           'couples exceeds upper limit for Hough tr. '       $           'couples exceeds upper limit for Hough tr. '
2125       $           ,'( ',ncp_max,' )'                   $           ,'( ',ncp_max,' )'            
# Line 2064  c     goto 880       !fill ntp and go to Line 2132  c     goto 880       !fill ntp and go to
2132        return        return
2133        end        end
2134    
 c$$$      subroutine cl_to_couples_2(iflag)  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'momanhough_init.f'  
 c$$$      include 'calib.f'  
 c$$$      include 'level1.f'  
 c$$$  
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
 c$$$  
 c$$$*     output flag  
 c$$$*     --------------  
 c$$$*     0 = good event  
 c$$$*     1 = bad event  
 c$$$*     --------------  
 c$$$      integer iflag  
 c$$$  
 c$$$      integer badseed,badcl  
 c$$$  
 c$$$*     init variables  
 c$$$      ncp_tot=0  
 c$$$      do ip=1,nplanes  
 c$$$         do ico=1,ncouplemax  
 c$$$            clx(ip,ico)=0  
 c$$$            cly(ip,ico)=0  
 c$$$         enddo  
 c$$$         ncp_plane(ip)=0  
 c$$$         do icl=1,nclstrmax_level2  
 c$$$            cls(ip,icl)=1  
 c$$$         enddo  
 c$$$         ncls(ip)=0  
 c$$$      enddo  
 c$$$      do icl=1,nclstrmax_level2  
 c$$$         cl_single(icl)=1  
 c$$$         cl_good(icl)=0  
 c$$$      enddo  
 c$$$        
 c$$$*     start association  
 c$$$      ncouples=0  
 c$$$      do icx=1,nclstr1          !loop on cluster (X)  
 c$$$         if(mod(VIEW(icx),2).eq.1)goto 10  
 c$$$          
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (X VIEW)  
 c$$$         if(dedx(icx).lt.dedx_x_min)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     cut BAD (X VIEW)              
 c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
 c$$$         ifirst=INDSTART(icx)  
 c$$$         if(icx.ne.nclstr1) then  
 c$$$            ilast=INDSTART(icx+1)-1  
 c$$$         else  
 c$$$            ilast=TOTCLLENGTH  
 c$$$         endif  
 c$$$         badcl=badseed  
 c$$$         do igood=-ngoodstr,ngoodstr  
 c$$$            ibad=1  
 c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.  
 c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.  
 c$$$     $           .true.)then  
 c$$$               ibad=BAD(VIEW(icx),  
 c$$$     $              nvk(MAXS(icx)+igood),  
 c$$$     $              nst(MAXS(icx)+igood))  
 c$$$            endif  
 c$$$            badcl=badcl*ibad  
 c$$$         enddo  
 c$$$*         print*,'icx ',icx,badcl  
 c$$$         if(badcl.eq.0)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     ----------------------------------------------------  
 c$$$          
 c$$$         cl_good(icx)=1  
 c$$$         nplx=npl(VIEW(icx))  
 c$$$         nldx=nld(MAXS(icx),VIEW(icx))  
 c$$$          
 c$$$         do icy=1,nclstr1       !loop on cluster (Y)  
 c$$$            if(mod(VIEW(icy),2).eq.0)goto 20  
 c$$$              
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (Y VIEW)  
 c$$$            if(dedx(icy).lt.dedx_y_min)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     cut BAD (Y VIEW)              
 c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
 c$$$            ifirst=INDSTART(icy)  
 c$$$            if(icy.ne.nclstr1) then  
 c$$$               ilast=INDSTART(icy+1)-1  
 c$$$            else  
 c$$$               ilast=TOTCLLENGTH  
 c$$$            endif  
 c$$$            badcl=badseed  
 c$$$            do igood=-ngoodstr,ngoodstr  
 c$$$               ibad=1  
 c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.  
 c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.  
 c$$$     $              .true.)  
 c$$$     $              ibad=BAD(VIEW(icy),  
 c$$$     $              nvk(MAXS(icy)+igood),  
 c$$$     $              nst(MAXS(icy)+igood))  
 c$$$               badcl=badcl*ibad  
 c$$$            enddo  
 c$$$*            print*,'icy ',icy,badcl  
 c$$$            if(badcl.eq.0)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     ----------------------------------------------------  
 c$$$              
 c$$$              
 c$$$            cl_good(icy)=1                    
 c$$$            nply=npl(VIEW(icy))  
 c$$$            nldy=nld(MAXS(icy),VIEW(icy))  
 c$$$              
 c$$$*     ----------------------------------------------  
 c$$$*     CONDITION TO FORM A COUPLE  
 c$$$*     ----------------------------------------------  
 c$$$*     geometrical consistency (same plane and ladder)  
 c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  
 c$$$  
 c$$$c$$$*     charge correlation  
 c$$$c$$$               ddd=(dedx(icy)  
 c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 c$$$                
 c$$$*     ------------------> COUPLE <------------------  
 c$$$*     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.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  
 c$$$                
 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  
 c$$$                
 c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1  
 c$$$               clx(nplx,ncp_plane(nplx))=icx  
 c$$$               cly(nply,ncp_plane(nplx))=icy  
 c$$$               cl_single(icx)=0  
 c$$$               cl_single(icy)=0  
 c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy  
 c$$$            endif                                
 c$$$*     ----------------------------------------------  
 c$$$  
 c$$$ 20         continue  
 c$$$         enddo                  !end loop on clusters(Y)  
 c$$$          
 c$$$ 10      continue  
 c$$$      enddo                     !end loop on clusters(X)  
 c$$$        
 c$$$        
 c$$$      do icl=1,nclstr1  
 c$$$         if(cl_single(icl).eq.1)then  
 c$$$            ip=npl(VIEW(icl))  
 c$$$            ncls(ip)=ncls(ip)+1  
 c$$$            cls(ip,ncls(ip))=icl  
 c$$$         endif  
 c$$$      enddo  
 c$$$        
 c$$$        
 c$$$      if(DEBUG)then  
 c$$$         print*,'clusters  ',nclstr1  
 c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)  
 c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)  
 c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
 c$$$      endif  
 c$$$        
 c$$$      do ip=1,6  
 c$$$         ncp_tot=ncp_tot+ncp_plane(ip)  
 c$$$      enddo  
 c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
 c$$$        
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(DEBUG)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$c            good2=.false.  
 c$$$c     goto 880       !fill ntp and go to next event  
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
 c$$$        
 c$$$      return  
 c$$$      end  
2135                
2136  ***************************************************  ***************************************************
2137  *                                                 *  *                                                 *
# Line 2365  c     $                       (icx2,icy2 Line 2224  c     $                       (icx2,icy2
2224  *     (2 couples needed)  *     (2 couples needed)
2225  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2226                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2227                             if(DEBUG)print*,                             if(verbose)print*,
2228       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2229       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2230       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
# Line 2435  c     $                                 Line 2294  c     $                                
2294  *     (3 couples needed)  *     (3 couples needed)
2295  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2296                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2297                                      if(DEBUG)print*,                                      if(verbose)print*,
2298       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2299       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2300       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
# Line 2662  c     print*,'>>>> ',ncpused,npt,nplused Line 2521  c     print*,'>>>> ',ncpused,npt,nplused
2521  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2522    
2523           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2524              if(DEBUG)print*,              if(verbose)print*,
2525       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2526       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2527       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
# Line 2874  c     print*,'check cp_used' Line 2733  c     print*,'check cp_used'
2733  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2734  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2735           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2736              if(DEBUG)print*,              if(verbose)print*,
2737       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2738       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2739       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
# Line 3154  c     $                                 Line 3013  c     $                                
3013                                enddo                                enddo
3014                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3015                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3016                                call mini_2(jstep,ifail)                                iprint=0
3017                                  if(DEBUG)iprint=1
3018                                  call mini2(jstep,ifail,iprint)
3019                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3020                                   if(DEBUG)then                                   if(DEBUG)then
3021                                      print *,                                      print *,
3022       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3023       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(mini2 in clouds_to_ctrack)'
3024                                   endif                                   endif
3025                                   chi2=-chi2                                   chi2=-chi2
3026                                endif                                endif
# Line 3174  c     $                                 Line 3035  c     $                                
3035  *     --------------------------  *     --------------------------
3036                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3037                                                                    
3038                                   if(DEBUG)print*,                                   if(verbose)print*,
3039       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3040       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3041       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
# Line 3800  c*************************************** Line 3661  c***************************************
3661        include 'level2.f'        include 'level2.f'
3662        include 'level1.f'        include 'level1.f'
3663    
3664          do i=1,nviews
3665             good2(i)=good1(i)
3666          enddo
3667    
3668    c      good2 = 0!.false.
       good2 = 0!.false.  
3669  c$$$      nev2 = nev1  c$$$      nev2 = nev1
3670    
3671  c$$$# ifndef TEST2003  c$$$# ifndef TEST2003
# Line 3904  c*************************************** Line 3767  c***************************************
3767        real sinth,phi,pig        !(4)        real sinth,phi,pig        !(4)
3768        pig=acos(-1.)        pig=acos(-1.)
3769    
3770        good2=1!.true.  c      good2=1!.true.
3771        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3772          nstep_nt(ntr)       = nstep
3773    
3774        phi   = al(4)             !(4)        phi   = al(4)             !(4)
3775        sinth = al(3)             !(4)        sinth = al(3)             !(4)
# Line 3989  c*************************************** Line 3853  c***************************************
3853        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3854    
3855  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
3856        good2=1!.true.  c      good2=1!.true.
3857        nclsx = 0        nclsx = 0
3858        nclsy = 0        nclsy = 0
3859    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23