/[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.27 by pam-fi, Fri Aug 17 14:36:05 2007 UTC revision 1.29 by pam-fi, Tue Aug 21 15:21:52 2007 UTC
# Line 215  c$$$      enddo Line 215  c$$$      enddo
215  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
216                
217        logical FIMAGE            !        logical FIMAGE            !
218          real trackimage(NTRACKSMAX)
219        real*8 AL_GUESS(5)        real*8 AL_GUESS(5)
220    
221  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 331  c$$$         enddo Line 332  c$$$         enddo
332              iimage=0              iimage=0
333           endif           endif
334           if(icand.eq.0)then           if(icand.eq.0)then
335              if(VERBOSE)then              if(VERBOSE.EQ.1)then
336                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
337       $              ,ibest,iimage       $              ,ibest,iimage
338              endif              endif
# Line 360  c         print*,'## guess: ',al Line 361  c         print*,'## guess: ',al
361           jstep=0                !# minimization steps           jstep=0                !# minimization steps
362    
363           iprint=0           iprint=0
364  c         if(DEBUG)iprint=1  c         if(DEBUG.EQ.1)iprint=1
365           if(VERBOSE)iprint=1           if(VERBOSE.EQ.1)iprint=1
366           if(DEBUG)iprint=2           if(DEBUG.EQ.1)iprint=2
367           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
368           if(ifail.ne.0) then           if(ifail.ne.0) then
369              if(VERBOSE)then              if(VERBOSE.EQ.1)then
370                 print *,                 print *,
371       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
372       $              ,iev       $              ,iev
# Line 380  c$$$               print*,'------------- Line 381  c$$$               print*,'-------------
381  c            chi2=-chi2  c            chi2=-chi2
382           endif           endif
383                    
384           if(DEBUG)then           if(DEBUG.EQ.1)then
385              print*,'----------------------------- improved track coord'              print*,'----------------------------- improved track coord'
386  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
387              do ip=1,6              do ip=1,6
# Line 391  c            chi2=-chi2 Line 392  c            chi2=-chi2
392           endif           endif
393    
394  c         rchi2=chi2/dble(ndof)  c         rchi2=chi2/dble(ndof)
395           if(DEBUG)then           if(DEBUG.EQ.1)then
396              print*,' '              print*,' '
397              print*,'****** SELECTED TRACK *************'              print*,'****** SELECTED TRACK *************'
398              print*,'#         R. chi2        RIG'              print*,'#         R. chi2        RIG'
# Line 407  c         rchi2=chi2/dble(ndof) Line 408  c         rchi2=chi2/dble(ndof)
408  *     ------------- search if the track has an IMAGE -------------  *     ------------- search if the track has an IMAGE -------------
409  *     ------------- (also this is stored )           -------------  *     ------------- (also this is stored )           -------------
410           if(FIMAGE)goto 122     !>>> jump! (this is already an image)           if(FIMAGE)goto 122     !>>> jump! (this is already an image)
411  *     now search for track-image, by comparing couples IDs  
412    *     -----------------------------------------------------
413    *     first check if the track is ambiguous
414    *     -----------------------------------------------------
415    *     (modified on august 2007 by ElenaV)
416             is1=0
417             do ip=1,NPLANES
418                if(SENSOR_STORE(ip,icand).ne.0)then
419                   is1=SENSOR_STORE(ip,icand)
420                   if(ip.eq.6)is1=3-is1 !last plane inverted
421                endif
422             enddo
423             if(is1.eq.0)then
424                if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
425                goto 122            !jump
426             endif
427    c         print*,'is1 ',is1
428             do ip=1,NPLANES
429                is2 = SENSOR_STORE(ip,icand) !sensor
430    c            print*,'is2 ',is2,' ip ',ip
431                if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
432                if(
433         $           (is1.ne.is2.and.is2.ne.0)
434         $           )goto 122      !jump (this track cannot have an image)
435             enddo
436             if(DEBUG.eq.1)print*,' >>> ambiguous track! '
437    *     now search for track-image among track candidates
438    c$$$         do i=1,ntracks
439    c$$$            iimage=i
440    c$$$            do ip=1,nplanes
441    c$$$               if(     CP_STORE(nplanes-ip+1,icand).ne.
442    c$$$     $              -1*CP_STORE(nplanes-ip+1,i).and.
443    c$$$     $              CP_STORE(nplanes-ip+1,i).ne.0.and.
444    c$$$     $              CP_STORE(nplanes-ip+1,icand).ne.0 )iimage=0
445    c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
446    c$$$     $              ,CP_STORE(nplanes-ip+1,icand),iimage
447    c$$$            enddo
448    c$$$            if(  iimage.ne.0.and.
449    c$$$c     $           RCHI2_STORE(i).le.CHI2MAX.and.
450    c$$$c     $           RCHI2_STORE(i).gt.0.and.
451    c$$$     $           .true.)then
452    c$$$               if(DEBUG.EQ.1)print*,'Track candidate ',iimage
453    c$$$     $              ,' >>> TRACK IMAGE >>> of'
454    c$$$     $              ,ibest
455    c$$$               goto 122         !image track found
456    c$$$            endif
457    c$$$         enddo
458    *     ---------------------------------------------------------------
459    *     take the candidate with the greatest number of matching couples
460    *     if more than one satisfies the requirement,
461    *     choose the one with more points and lower chi2
462    *     ---------------------------------------------------------------
463    *     count the number of matching couples
464             ncpmax = 0
465             iimage   = 0           !id of candidate with better matching
466           do i=1,ntracks           do i=1,ntracks
467              iimage=i              ncp=0
468              do ip=1,nplanes              do ip=1,nplanes
469                 if(     CP_STORE(nplanes-ip+1,icand).ne.                 if(CP_STORE(nplanes-ip+1,icand).ne.0)then
470       $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0                    if(
471         $                 CP_STORE(nplanes-ip+1,i).ne.0
472         $                 .and.
473         $                 CP_STORE(nplanes-ip+1,icand).eq.
474         $                 -1*CP_STORE(nplanes-ip+1,i)
475         $                 )then
476                         ncp=ncp+1  !ok
477                      elseif(
478         $                    CP_STORE(nplanes-ip+1,i).ne.0
479         $                    .and.
480         $                    CP_STORE(nplanes-ip+1,icand).ne.
481         $                    -1*CP_STORE(nplanes-ip+1,i)
482         $                    )then
483                         ncp = 0
484                         goto 117   !it is not an image candidate
485                      else
486                        
487                      endif
488                   endif
489    c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
490    c$$$     $              ,CP_STORE(nplanes-ip+1,icand),ncp
491              enddo              enddo
492              if(  iimage.ne.0.and.   117        continue
493  c     $           RCHI2_STORE(i).le.CHI2MAX.and.              trackimage(i)=ncp   !number of matching couples
494  c     $           RCHI2_STORE(i).gt.0.and.              if(ncp>ncpmax)then
495       $           .true.)then                 ncpmax=ncp
496                 if(DEBUG)print*,'Track candidate ',iimage                 iimage=i
      $              ,' >>> TRACK IMAGE >>> of'  
      $              ,ibest  
                goto 122         !image track found  
497              endif              endif
498           enddo           enddo
499    *     check if there are more than one image candidates
500             ngood=0
501             do i=1,ntracks
502                if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
503             enddo
504    *     if there are, choose the best one
505             if(ngood.gt.1)then
506    *     -------------------------------------------------------
507    *     order track-candidates according to:
508    *     1st) decreasing n.points
509    *     2nd) increasing chi**2
510    *     -------------------------------------------------------
511                rchi2best=1000000000.
512                ndofbest=0            
513                do i=1,ntracks
514                   if( trackimage(i).eq.ncpmax )then
515                      ndof=0              
516                      do ii=1,nplanes    
517                         ndof=ndof        
518         $                    +int(xgood_store(ii,i))
519         $                    +int(ygood_store(ii,i))
520                      enddo              
521                      if(ndof.gt.ndofbest)then
522                         iimage=i
523                         rchi2best=RCHI2_STORE(i)
524                         ndofbest=ndof    
525                      elseif(ndof.eq.ndofbest)then
526                         if(RCHI2_STORE(i).lt.rchi2best.and.
527         $                    RCHI2_STORE(i).gt.0)then
528                            iimage=i
529                            rchi2best=RCHI2_STORE(i)
530                            ndofbest=ndof  
531                         endif            
532                      endif
533                   endif
534                enddo
535                
536             endif
537    
538             if(DEBUG.EQ.1)print*,'Track candidate ',iimage
539         $        ,' >>> TRACK IMAGE >>> of'
540         $        ,ibest
541    
542   122     continue   122     continue
543    
544    
545  *     --- and store the results --------------------------------  *     --- and store the results --------------------------------
546           ntrk = ntrk + 1                   !counter of found tracks           ntrk = ntrk + 1                   !counter of found tracks
547           if(.not.FIMAGE           if(.not.FIMAGE
# Line 438  c         print*,'++++++++++ iimage,fima Line 554  c         print*,'++++++++++ iimage,fima
554  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
555    
556           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
557              if(verbose)              if(verbose.eq.1)
558       $           print*,       $           print*,
559       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
560       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 474  cc            good2=.false. Line 590  cc            good2=.false.
590       $        rchi2best.le.CHI2MAX.and.       $        rchi2best.le.CHI2MAX.and.
591  c     $        rchi2best.lt.15..and.  c     $        rchi2best.lt.15..and.
592       $        .true.)then       $        .true.)then
593              if(DEBUG)then              if(DEBUG.EQ.1)then
594                 print*,'***** NEW SEARCH ****'                 print*,'***** NEW SEARCH ****'
595              endif              endif
596              goto 11111          !try new search              goto 11111          !try new search
# Line 690  c$$$         print*,fbad_cog(4,icx) Line 806  c$$$         print*,fbad_cog(4,icx)
806              stripx  = stripx + pfaeta2(icx,angx)                        stripx  = stripx + pfaeta2(icx,angx)          
807              resxPAM = risxeta2(abs(angx))              resxPAM = risxeta2(abs(angx))
808              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
809              if(DEBUG.and.fbad_cog(2,icx).ne.1)  c$$$            if(DEBUG.EQ.1.and.fbad_cog(2,icx).ne.1)
810       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  c$$$     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
811    
812           elseif(PFAx.eq.'ETA3')then                                   elseif(PFAx.eq.'ETA3')then                        
813    
# Line 733  c     $           print*,'BAD icx >>> ', Line 849  c     $           print*,'BAD icx >>> ',
849              resxPAM = resxPAM*fbad_cog(0,icx)              resxPAM = resxPAM*fbad_cog(0,icx)
850    
851           else           else
852              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx
853           endif           endif
854    
855    
# Line 743  c     $           print*,'BAD icx >>> ', Line 859  c     $           print*,'BAD icx >>> ',
859           if( nsatstrips(icx).gt.0 )then           if( nsatstrips(icx).gt.0 )then
860              stripx  = stripx + cog(4,icx)                          stripx  = stripx + cog(4,icx)            
861              resxPAM = pitchX*1e-4/sqrt(12.)              resxPAM = pitchX*1e-4/sqrt(12.)
862              cc=cog(4,icx)  cc            cc=cog(4,icx)
863  c$$$            print*,icx,' *** ',cc  c$$$            print*,icx,' *** ',cc
864  c$$$            print*,icx,' *** ',resxPAM  c$$$            print*,icx,' *** ',resxPAM
865           endif           endif
# Line 778  c$$$            print*,icx,' *** ',resxP Line 894  c$$$            print*,icx,' *** ',resxP
894           endif           endif
895    
896           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
897              if(DEBUG) then              if(DEBUG.EQ.1) then
898                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
899       $              ,icx,icy       $              ,icx,icy
900              endif              endif
# Line 873  c     $           print*,'BAD icy >>> ', Line 989  c     $           print*,'BAD icy >>> ',
989              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
990    
991           else           else
992              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx
993           endif           endif
994    
995    
# Line 883  c     $           print*,'BAD icy >>> ', Line 999  c     $           print*,'BAD icy >>> ',
999           if( nsatstrips(icy).gt.0 )then           if( nsatstrips(icy).gt.0 )then
1000              stripy  = stripy + cog(4,icy)                          stripy  = stripy + cog(4,icy)            
1001              resyPAM = pitchY*1e-4/sqrt(12.)              resyPAM = pitchY*1e-4/sqrt(12.)
1002              cc=cog(4,icy)  cc            cc=cog(4,icy)
1003  c$$$            print*,icy,' *** ',cc  c$$$            print*,icy,' *** ',cc
1004  c$$$            print*,icy,' *** ',resyPAM  c$$$            print*,icy,' *** ',resyPAM
1005           endif           endif
# Line 903  c     (xi,yi,zi) = mechanical coordinate Line 1019  c     (xi,yi,zi) = mechanical coordinate
1019  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1020           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1021       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1022              if(DEBUG) then              if(DEBUG.EQ.1) then
1023                 print*,'xyz_PAM (couple):',                 print*,'xyz_PAM (couple):',
1024       $              ' WARNING: false X strip: strip ',stripx       $              ' WARNING: false X strip: strip ',stripx
1025              endif              endif
# Line 998  c            print*,'X-singlet ',icx,npl Line 1114  c            print*,'X-singlet ',icx,npl
1114  c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1115              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1116       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1117                 if(DEBUG) then                 if(DEBUG.EQ.1) then
1118                    print*,'xyz_PAM (X-singlet):',                    print*,'xyz_PAM (X-singlet):',
1119       $                 ' WARNING: false X strip: strip ',stripx       $                 ' WARNING: false X strip: strip ',stripx
1120                 endif                 endif
# Line 1023  c            print*,'X-cl ',icx,stripx,' Line 1139  c            print*,'X-cl ',icx,stripx,'
1139  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
1140    
1141           else           else
1142              if(DEBUG) then              if(DEBUG.EQ.1) then
1143                 print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
1144                 print *,'icx = ',icx                 print *,'icx = ',icx
1145                 print *,'icy = ',icy                 print *,'icy = ',icy
# Line 1092  c--------------------------------------- Line 1208  c---------------------------------------
1208  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1209    
1210        else        else
1211           if(DEBUG) then           if(DEBUG.EQ.1) then
1212              print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1213              print *,'icx = ',icx              print *,'icx = ',icx
1214              print *,'icy = ',icy              print *,'icy = ',icy
# Line 1137  c$$$      common/FINALPFA/PFA Line 1253  c$$$      common/FINALPFA/PFA
1253  c$$$      PFAx = 'COG4'!PFA  c$$$      PFAx = 'COG4'!PFA
1254  c$$$      PFAy = 'COG4'!PFA  c$$$      PFAy = 'COG4'!PFA
1255    
1256    
1257        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1258              print*,'xyzpam: ***WARNING*** clusters ',icx,icy              print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1259       $           ,' does not exists (nclstr1=',nclstr1,')'       $           ,' does not exists (nclstr1=',nclstr1,')'
# Line 1149  c$$$      PFAy = 'COG4'!PFA Line 1266  c$$$      PFAy = 'COG4'!PFA
1266        call idtoc(pfaid,PFAx)        call idtoc(pfaid,PFAx)
1267        call idtoc(pfaid,PFAy)        call idtoc(pfaid,PFAy)
1268    
1269    cc      print*,PFAx,PFAy
1270    
1271  c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)  c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1272    
1273  c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy  c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
# Line 1277  c         zv(ip) = z_mech_sensor(nplanes Line 1396  c         zv(ip) = z_mech_sensor(nplanes
1396    
1397        endif        endif
1398    
1399        if(DEBUG)then        if(DEBUG.EQ.1)then
1400  c         print*,'----------------------------- track coord'  c         print*,'----------------------------- track coord'
1401  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1402           write(*,22222)ip,zm(ip),xm(ip),ym(ip)           write(*,22222)ip,zm(ip),xm(ip),ym(ip)
# Line 1728  c      include 'level1.f' Line 1847  c      include 'level1.f'
1847        integer iflag        integer iflag
1848    
1849        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1850        
1851          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1852    
1853  *     init variables  *     init variables
1854        ncp_tot=0        ncp_tot=0
# Line 1760  c      include 'level1.f' Line 1881  c      include 'level1.f'
1881           if( ncl_view(iv).gt. nclusterlimit)then           if( ncl_view(iv).gt. nclusterlimit)then
1882  c            mask_view(iv) = 1  c            mask_view(iv) = 1
1883              mask_view(iv) = mask_view(iv) + 2**0              mask_view(iv) = mask_view(iv) + 2**0
1884              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG.EQ.1)
1885       $           ,nclusterlimit,' on view ', iv,' --> masked!'       $        print*,' * WARNING * cl_to_couple: n.clusters > '
1886         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1887           endif           endif
1888        enddo        enddo
1889    
# Line 1899  c                  cut = chcut * sch(npl Line 2021  c                  cut = chcut * sch(npl
2021                 endif                 endif
2022    
2023                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
2024                    if(verbose)print*,                    if(verbose.eq.1)print*,
2025       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
2026       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
2027       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
# Line 1937  c                  mask_view(nviewy(nply Line 2059  c                  mask_view(nviewy(nply
2059        enddo        enddo
2060                
2061                
2062        if(DEBUG)then        if(DEBUG.EQ.1)then
2063           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
2064           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
2065           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
2066           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
2067        endif        endif
2068                
# Line 2000  c      double precision xm3,ym3,zm3 Line 2122  c      double precision xm3,ym3,zm3
2122        real xc,zc,radius        real xc,zc,radius
2123  *     -----------------------------  *     -----------------------------
2124    
2125          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
2126    
2127  *     --------------------------------------------  *     --------------------------------------------
2128  *     put a limit to the maximum number of couples  *     put a limit to the maximum number of couples
# Line 2037  c     print*,'***',is1,xm1,ym1,zm1 Line 2160  c     print*,'***',is1,xm1,ym1,zm1
2160                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
2161                    if(  mask_view(nviewx(ip2)).ne.0 .or.                    if(  mask_view(nviewx(ip2)).ne.0 .or.
2162       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
2163                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
                       
2164                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
2165                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
2166                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
# Line 2057  c     $                       (icx2,icy2 Line 2179  c     $                       (icx2,icy2
2179  *     (2 couples needed)  *     (2 couples needed)
2180  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2181                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2182                             if(verbose)print*,                             if(verbose.eq.1)print*,
2183       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2184       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2185       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
# Line 2137  c     $                                 Line 2259  c     $                                
2259  *     (3 couples needed)  *     (3 couples needed)
2260  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2261                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2262                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2263       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2264       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2265       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
# Line 2199  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2321  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2321   10   continue   10   continue
2322        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2323                
2324        if(DEBUG)then        if(DEBUG.EQ.1)then
2325           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2326           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2327        endif        endif
# Line 2246  c      include 'momanhough_init.f' Line 2368  c      include 'momanhough_init.f'
2368        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2369        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2370    
2371          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2372    
2373  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2374  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2372  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2495  c         if(ncpused.lt.ncpyz_min)goto 2
2495  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2496    
2497           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2498              if(verbose)print*,              if(verbose.eq.1)print*,
2499       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2500       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2501       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
# Line 2400  c     ptcloud_yz_nt(nclouds_yz)=npt Line 2523  c     ptcloud_yz_nt(nclouds_yz)=npt
2523  c     print*,'>> ',ipt,db_all(ipt)  c     print*,'>> ',ipt,db_all(ipt)
2524           enddo             enddo  
2525           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2526           if(DEBUG)then           if(DEBUG.EQ.1)then
2527              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2528              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2529              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2530              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2531              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2532              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2533              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2534         $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2535                print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
2536  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
2537  c$$$     $           ,ptcloud_yz(nclouds_yz)  c$$$     $           ,ptcloud_yz(nclouds_yz)
2538  c$$$            print*,'nt-uple: db_cloud(...) = '  c$$$            print*,'nt-uple: db_cloud(...) = '
# Line 2425  c$$$     $           ,(db_cloud(iii),iii Line 2550  c$$$     $           ,(db_cloud(iii),iii
2550          goto 90                          goto 90                
2551        endif                            endif                    
2552                
2553        if(DEBUG)then        if(DEBUG.EQ.1)then
2554           print*,'---------------------- '           print*,'---------------------- '
2555           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
2556           print*,' '           print*,' '
# Line 2474  c      include 'momanhough_init.f' Line 2599  c      include 'momanhough_init.f'
2599        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2600        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2601    
2602          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2603    
2604  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2605  *     classification of TRIPLETS  *     classification of TRIPLETS
2606  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2531  c         tr_temp(1)=itr1 Line 2658  c         tr_temp(1)=itr1
2658       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2659                 distance = sqrt(distance)                 distance = sqrt(distance)
2660                                
2661                 if(distance.lt.cutdistxz)then  *     ------------------------------------------------------------------------
2662    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2663    *     ------------------------------------------------------------------------
2664    *     (added in august 2007)
2665                   istrimage=0
2666                   if(
2667         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2668         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2669         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2670         $              .true.)istrimage=1
2671    
2672                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2673  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
2674                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2675                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
# Line 2595  c         if(ncpused.lt.ncpxz_min)goto 2 Line 2733  c         if(ncpused.lt.ncpxz_min)goto 2
2733  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2734  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2735           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2736              if(verbose)print*,              if(verbose.eq.1)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 2622  c               mask_view(iv) = 6 Line 2760  c               mask_view(iv) = 6
2760           enddo           enddo
2761           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2762                    
2763           if(DEBUG)then           if(DEBUG.EQ.1)then
2764              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2765              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2766              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2767              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2768              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2769              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2770              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2771                print*,'cpcloud_xz '
2772         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2773              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2774  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2775  c$$$     $           ,ptcloud_xz(nclouds_xz)  c$$$     $           ,ptcloud_xz(nclouds_xz)
# Line 2647  c$$$     $           ,(tr_cloud(iii),iii Line 2787  c$$$     $           ,(tr_cloud(iii),iii
2787           goto 91                           goto 91                
2788         endif                             endif                    
2789                
2790        if(DEBUG)then        if(DEBUG.EQ.1)then
2791           print*,'---------------------- '           print*,'---------------------- '
2792           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
2793           print*,' '           print*,' '
# Line 2699  c$$$     $           ,(tr_cloud(iii),iii Line 2839  c$$$     $           ,(tr_cloud(iii),iii
2839  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2840  *     variables for track fitting  *     variables for track fitting
2841        double precision AL_INI(5)        double precision AL_INI(5)
 c      double precision tath  
2842  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2843    
2844          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2845    
2846    
2847        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
# Line 2724  c      real fitz(nplanes)        !z coor Line 2863  c      real fitz(nplanes)        !z coor
2863                 enddo                 enddo
2864              enddo              enddo
2865              ncp_ok=0              ncp_ok=0
2866              do icp=1,ncp_tot    !loop on couples              do icp=1,ncp_tot    !loop over couples
2867  *     get info on  *     get info on
2868                 cpintersec(icp)=min(                 cpintersec(icp)=min(
2869       $              cpcloud_yz(iyz,icp),       $              cpcloud_yz(iyz,icp),
# Line 2733  c      real fitz(nplanes)        !z coor Line 2872  c      real fitz(nplanes)        !z coor
2872       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2873       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2874       $              .false.)cpintersec(icp)=0       $              .false.)cpintersec(icp)=0
2875    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2876                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
2877                    ncp_ok=ncp_ok+1                      ncp_ok=ncp_ok+1  
2878                                        
# Line 2765  c      real fitz(nplanes)        !z coor Line 2905  c      real fitz(nplanes)        !z coor
2905                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2906              enddo              enddo
2907                            
 c            if(nplused.lt.nplxz_min)goto 888 !next doublet  
             if(nplused.lt.nplyz_min)goto 888 !next doublet  
             if(ncp_ok.lt.ncpok)goto 888 !next cloud  
2908                            
2909              if(DEBUG)then              if(DEBUG.EQ.1)then
2910                 print*,'Combination ',iyz,ixz                 print*,'Combination ',iyz,ixz
2911       $              ,' db ',ptcloud_yz(iyz)       $              ,' db ',ptcloud_yz(iyz)
2912       $              ,' tr ',ptcloud_xz(ixz)       $              ,' tr ',ptcloud_xz(ixz)
2913       $              ,'  -----> # matching couples ',ncp_ok       $              ,'  -----> # matching couples ',ncp_ok
2914              endif              endif
2915    
2916    c            if(nplused.lt.nplxz_min)goto 888 !next combination
2917                if(nplused.lt.nplyz_min)goto 888 !next combination
2918                if(ncp_ok.lt.ncpok)goto 888 !next combination
2919    
2920  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2921  c$$$  print*,'Configurazione cluster XZ'  c$$$  print*,'Configurazione cluster XZ'
2922  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
# Line 2820  c$$$            AL_INI(5) = (1.e2*alfaxz Line 2962  c$$$            AL_INI(5) = (1.e2*alfaxz
2962  c$$$              c$$$            
2963  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2964                                                    
2965              if(DEBUG)then              if(DEBUG.EQ.1)then
2966                   print*,'track candidate', ntracks+1
2967                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2968                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2969                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2853  c$$$            if(AL_INI(5).gt.defmax)g Line 2996  c$$$            if(AL_INI(5).gt.defmax)g
2996                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2997                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2998                                                                
2999                                  *                             ---------------------------------------
3000    *                             check if this group of couples has been
3001    *                             already fitted    
3002    *                             ---------------------------------------
3003                                  do ica=1,ntracks
3004                                     isthesame=1
3005                                     do ip=1,NPLANES
3006                                        if(hit_plane(ip).ne.0)then
3007                                           if(  CP_STORE(nplanes-ip+1,ica)
3008         $                                      .ne.
3009         $                                      cp_match(ip,hit_plane(ip)) )
3010         $                                      isthesame=0
3011                                        else
3012                                           if(  CP_STORE(nplanes-ip+1,ica)
3013         $                                      .ne.
3014         $                                      0 )
3015         $                                      isthesame=0
3016                                        endif
3017                                     enddo
3018                                     if(isthesame.eq.1)then
3019                                        if(DEBUG.eq.1)
3020         $                                   print*,'(already fitted)'
3021                                        goto 666 !jump to next combination
3022                                     endif
3023                                  enddo
3024    
3025                                call track_init !init TRACK common                                call track_init !init TRACK common
3026    
3027                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
3028                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
3029                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
3030                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2900  c$$$                              enddo Line 3068  c$$$                              enddo
3068                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3069                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3070                                iprint=0                                iprint=0
3071  c                              if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
3072                                if(DEBUG)iprint=2                                if(DEBUG.EQ.1)iprint=2
3073                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
3074                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3075                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
3076                                      print *,                                      print *,
3077       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3078       $                              //'(clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
# Line 2929  c                                 chi2=- Line 3097  c                                 chi2=-
3097  *     --------------------------  *     --------------------------
3098                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3099                                                                    
3100                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
3101       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3102       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3103       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
# Line 2945  c                                    mas Line 3113  c                                    mas
3113                                                                
3114                                ntracks = ntracks + 1                                ntracks = ntracks + 1
3115                                                                
3116                                do ip=1,nplanes                                do ip=1,nplanes !top to bottom
3117    
3118                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
3119                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
# Line 2962  c                                    mas Line 3130  c                                    mas
3130                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
3131                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
3132                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
3133    *                                NB! hit_plane is defined from bottom to top
3134                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
3135                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
3136       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
# Line 2977  c                                    mas Line 3146  c                                    mas
3146                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0
3147                                      LADDER_STORE(nplanes-ip+1,ntracks)=0                                      LADDER_STORE(nplanes-ip+1,ntracks)=0
3148                                   endif                                   endif
3149                                   BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now                                   BX_STORE(ip,ntracks)=0!I dont need it now
3150                                   BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now                                   BY_STORE(ip,ntracks)=0!I dont need it now
3151                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(ip,ntracks)=0
3152                                   do i=1,5                                   do i=1,5
3153                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
3154                                   enddo                                   enddo
# Line 3008  c                                    mas Line 3177  c                                    mas
3177           return           return
3178        endif        endif
3179                
3180  c$$$      if(DEBUG)then  c$$$      if(DEBUG.EQ.1)then
3181  c$$$         print*,'****** TRACK CANDIDATES ***********'  c$$$         print*,'****** TRACK CANDIDATES ***********'
3182  c$$$         print*,'#         R. chi2        RIG'  c$$$         print*,'#         R. chi2        RIG'
3183  c$$$         do i=1,ntracks  c$$$         do i=1,ntracks
# Line 3017  c$$$     $           ,1./abs(AL_STORE(5, Line 3186  c$$$     $           ,1./abs(AL_STORE(5,
3186  c$$$         enddo  c$$$         enddo
3187  c$$$         print*,'***********************************'  c$$$         print*,'***********************************'
3188  c$$$      endif  c$$$      endif
3189        if(DEBUG)then        if(DEBUG.EQ.1)then
3190          print*,'****** TRACK CANDIDATES *****************'          print*,'****** TRACK CANDIDATES *****************'
3191          print*,'#         R. chi2        RIG         ndof'          print*,'#         R. chi2        RIG         ndof'
3192          do i=1,ntracks          do i=1,ntracks
# Line 3069  c$$$      endif Line 3238  c$$$      endif
3238        real xyzp(3),bxyz(3)        real xyzp(3),bxyz(3)
3239        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3240    
3241          if(DEBUG.EQ.1)print*,'refine_track:'
3242  *     =================================================  *     =================================================
3243  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3244  *                          and  *                          and
# Line 3150  c     $           AYV_STORE(nplanes-ip+1 Line 3320  c     $           AYV_STORE(nplanes-ip+1
3320              LADDER_STORE(nplanes-ip+1,IBEST)=nldt              LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3321  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3322    
3323              if(DEBUG)then              if(DEBUG.EQ.1)then
3324                 print*,                 print*,
3325       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3326       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3161  c     $           AYV_STORE(nplanes-ip+1 Line 3331  c     $           AYV_STORE(nplanes-ip+1
3331  *     ===========================================  *     ===========================================
3332  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3333  *     ===========================================  *     ===========================================
3334  c            if(DEBUG)print*,'>>>> try to include a new couple'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
3335              distmin=1000000.              distmin=1000000.
3336              xmm = 0.              xmm = 0.
3337              ymm = 0.              ymm = 0.
# Line 3193  c     $              cl_used(icy).eq.1.o Line 3363  c     $              cl_used(icy).eq.1.o
3363                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3364  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3365                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3366                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)print*,'( couple ',id
3367       $              ,' ) distance ',distance       $              ,' ) distance ',distance
3368                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3369                    xmm = xPAM                    xmm = xPAM
# Line 3225  c            if(distmin.le.clinc)then   Line 3395  c            if(distmin.le.clinc)then  
3395                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3396  *              -----------------------------------  *              -----------------------------------
3397                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3398                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3399       $              ,' (dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3400                 goto 133         !next plane                 goto 133         !next plane
3401              endif              endif
# Line 3233  c            if(distmin.le.clinc)then   Line 3403  c            if(distmin.le.clinc)then  
3403  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3404  *                     either from a couple or single  *                     either from a couple or single
3405  *     ================================================  *     ================================================
3406  c            if(DEBUG)print*,'>>>> try to include a new cluster'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
3407              distmin=1000000.              distmin=1000000.
3408              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3409              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3270  c     $              AXV_STORE(nplanes-i Line 3440  c     $              AXV_STORE(nplanes-i
3440       $              )                     $              )              
3441                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3442  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3443                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)print*,'( cl-X ',icx
3444       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3445                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3446                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3303  c     $              0.,AYV_STORE(nplane Line 3473  c     $              0.,AYV_STORE(nplane
3473       $              )       $              )
3474                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3475  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3476                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy
3477       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3478                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3479                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3348  c               if(cl_used(icl).eq.1.or. Line 3518  c               if(cl_used(icl).eq.1.or.
3518    
3519                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3520  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3521                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)print*,'( cl-s ',icl
3522       $              ,' ) distance ',distance,'<',distmin,' ?'       $              ,' ) distance ',distance
3523                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3524                    if(DEBUG)print*,'YES'  c                  if(DEBUG.EQ.1)print*,'YES'
3525                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3526                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3527                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3397  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3567  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3567                    if(mod(VIEW(iclm),2).eq.0)then                    if(mod(VIEW(iclm),2).eq.0)then
3568                       XGOOD(nplanes-ip+1)=1.                       XGOOD(nplanes-ip+1)=1.
3569                       resx(nplanes-ip+1)=rxmm                       resx(nplanes-ip+1)=rxmm
3570                       if(DEBUG)print*,'%%%% included X-cl ',iclm                       if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
3571       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3572       $                    ,', dist.= ',distmin       $                    ,', dist.= ',distmin
3573       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clinc,' )'
3574                    else                    else
3575                       YGOOD(nplanes-ip+1)=1.                       YGOOD(nplanes-ip+1)=1.
3576                       resy(nplanes-ip+1)=rymm                       resy(nplanes-ip+1)=rymm
3577                       if(DEBUG)print*,'%%%% included Y-cl ',iclm                       if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
3578       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3579       $                    ,', dist.= ', distmin       $                    ,', dist.= ', distmin
3580       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clinc,' )'
# Line 3446  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3616  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3616        include 'common_momanhough.f'        include 'common_momanhough.f'
3617        include 'level2.f'              include 'level2.f'      
3618    
3619          if(DEBUG.EQ.1)print*,'clean_XYclouds:'
3620    
3621        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3622    
3623           id=CP_STORE(nplanes-ip+1,ibest)           id=CP_STORE(nplanes-ip+1,ibest)
# Line 3454  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3626  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3626              if(id.ne.0)then              if(id.ne.0)then
3627                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3628                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3629                 cl_used(iclx)=ntrk  !tag used clusters  c$$$               cl_used(iclx)=ntrk  !tag used clusters
3630                 cl_used(icly)=ntrk  !tag used clusters  c$$$               cl_used(icly)=ntrk  !tag used clusters
3631              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3632                 cl_used(icl)=ntrk   !tag used clusters  c$$$               cl_used(icl)=ntrk   !tag used clusters
3633              endif              endif
3634                            
3635  *     -----------------------------  *     -----------------------------
# Line 3476  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3648  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3648       $              cly(ip,icp).eq.icl       $              cly(ip,icp).eq.icl
3649       $              )then       $              )then
3650                    id=id_cp(ip,icp,1)                    id=id_cp(ip,icp,1)
3651                    if(DEBUG)then                    if(DEBUG.EQ.1)then
3652                       print*,ip,' <<< cp ',id                       print*,ip,' <<< cp ',id
3653       $                    ,' ( cl-x '       $                    ,' ( cl-x '
3654       $                    ,clx(ip,icp)       $                    ,clx(ip,icp)
# Line 3736  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3908  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3908           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3909           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3910           ayv_nt(ip,ntr)   = sngl(ayv(ip))             ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3911  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
3912           factor = sqrt(           factor = sqrt(
3913       $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3914       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3915       $        1. )       $        1. )
3916  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
3917           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3918           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3919        
# Line 3771  c         print*,'* ',ip,bfx,bfy,angx,an Line 3943  c         print*,'* ',ip,bfx,bfy,angx,an
3943  c           >>> is a couple  c           >>> is a couple
3944              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3945              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3946                
3947                cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3948                cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3949    
3950  c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)  c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3951  c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)              c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3952  c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))  c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
# Line 3787  c$$$            ybad(ip,ntr)= nbadstrips Line 3962  c$$$            ybad(ip,ntr)= nbadstrips
3962    
3963           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3964    
3965                cl_used(icl) = 1    !tag used clusters          
3966    
3967              if(mod(VIEW(icl),2).eq.0)then              if(mod(VIEW(icl),2).eq.0)then
3968                 cltrx(ip,ntr)=icl                 cltrx(ip,ntr)=icl
3969  c$$$               nnnnn = npfastrips(icl,PFA,angx)  c$$$               nnnnn = npfastrips(icl,PFA,angx)
# Line 3806  c$$$               ybad(ip,ntr) = nbadst Line 3983  c$$$               ybad(ip,ntr) = nbadst
3983    
3984        enddo        enddo
3985    
3986          if(DEBUG.eq.1)then
3987             print*,'> STORING TRACK ',ntr
3988             print*,'clusters: '
3989             do ip=1,6
3990                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3991             enddo
3992          endif
3993    
3994  c$$$      print*,(xgood(i),i=1,6)  c$$$      print*,(xgood(i),i=1,6)
3995  c$$$      print*,(ygood(i),i=1,6)  c$$$      print*,(ygood(i),i=1,6)
# Line 3840  c         if( mask_view(iv).ne.0 )good2( Line 4024  c         if( mask_view(iv).ne.0 )good2(
4024           good2(iv) = good2(iv) + mask_view(iv)*2**8           good2(iv) = good2(iv) + mask_view(iv)*2**8
4025        enddo        enddo
4026    
4027          if(DEBUG.eq.1)then
4028             print*,'> STORING SINGLETS '
4029          endif
4030    
4031        do icl=1,nclstr1        do icl=1,nclstr1
4032    
4033             ip=nplanes-npl(VIEW(icl))+1            
4034            
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
             ip=nplanes-npl(VIEW(icl))+1              
4036              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
4037                 nclsx = nclsx + 1                 nclsx = nclsx + 1
4038                 planex(nclsx) = ip                 planex(nclsx) = ip
# Line 3882  c$$$               print*,'ys(2,nclsy)   Line 4072  c$$$               print*,'ys(2,nclsy)  
4072    
4073  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4074           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
4075    *     --------------------------------------------------
4076    *     per non perdere la testa...
4077    *     whichtrack(icl) e` una variabile del common level1
4078    *     che serve solo per sapere quali cluster sono stati
4079    *     associati ad una traccia, e permettere di salvare
4080    *     solo questi nell'albero di uscita
4081    *     --------------------------------------------------
4082            
4083    
4084    c$$$         print*,' cl ',icl,' --> ',cl_used(icl)
4085    c$$$
4086    c$$$         if( cl_used(icl).ne.0 )then
4087    c$$$            if(
4088    c$$$     $           mod(VIEW(icl),2).eq.0.and.
4089    c$$$     $           cltrx(ip,whichtrack(icl)).ne.icl )
4090    c$$$     $           print*,'**WARNING** cltrx(',ip,',',whichtrack(icl)
4091    c$$$     $           ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl
4092    c$$$            if(
4093    c$$$     $           mod(VIEW(icl),2).eq.1.and.
4094    c$$$     $           cltry(ip,whichtrack(icl)).ne.icl )
4095    c$$$     $           print*,'**WARNING** cltry(',ip,',',whichtrack(icl)
4096    c$$$     $           ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl
4097    c$$$         endif
4098            
4099    
4100        enddo        enddo
4101        end        end
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.23