/[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.34 by pam-fi, Wed Mar 5 17:00:20 2008 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 572  c     $        rchi2best.lt.15..and. Line 688  c     $        rchi2best.lt.15..and.
688    
689        real stripx,stripy        real stripx,stripy
690    
691          double precision xi,yi,zi
692          double precision xi_A,yi_A,zi_A
693          double precision xi_B,yi_B,zi_B
694        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
695        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
696        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
# Line 585  c$$$      print*,icx,icy,sensor,PFAx,PFA Line 704  c$$$      print*,icx,icy,sensor,PFAx,PFA
704        resxPAM = 0        resxPAM = 0
705        resyPAM = 0        resyPAM = 0
706    
707        xPAM = 0.        xPAM = 0.D0
708        yPAM = 0.        yPAM = 0.D0
709        zPAM = 0.        zPAM = 0.D0
710        xPAM_A = 0.        xPAM_A = 0.D0
711        yPAM_A = 0.        yPAM_A = 0.D0
712        zPAM_A = 0.        zPAM_A = 0.D0
713        xPAM_B = 0.        xPAM_B = 0.D0
714        yPAM_B = 0.        yPAM_B = 0.D0
715        zPAM_B = 0.        zPAM_B = 0.D0
716  c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy  cc      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
717    
718        if(sensor.lt.1.or.sensor.gt.2)then        if(sensor.lt.1.or.sensor.gt.2)then
719           print*,'xyz_PAM   ***ERROR*** wrong input '           print*,'xyz_PAM   ***ERROR*** wrong input '
# Line 611  c      print*,'## xyz_PAM: ',icx,icy,sen Line 730  c      print*,'## xyz_PAM: ',icx,icy,sen
730           viewx   = VIEW(icx)           viewx   = VIEW(icx)
731           nldx    = nld(MAXS(icx),VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
732           nplx    = npl(VIEW(icx))           nplx    = npl(VIEW(icx))
733           resxPAM = RESXAV  c         resxPAM = RESXAV
734           stripx  = float(MAXS(icx))           stripx  = float(MAXS(icx))
735    
736           if(           if(
# Line 631  c      print*,'## xyz_PAM: ',icx,icy,sen Line 750  c      print*,'## xyz_PAM: ',icx,icy,sen
750  *        --------------------------  *        --------------------------
751  *        magnetic-field corrections  *        magnetic-field corrections
752  *        --------------------------  *        --------------------------
753           angtemp  = ax           stripx  = stripx +  fieldcorr(viewx,bfy)      
754           bfytemp  = bfy           angx    = effectiveangle(ax,viewx,bfy)
 *        /////////////////////////////////  
 *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!  
 *        *grvzkkjsdgjhhhgngbn###>:(  
 *        /////////////////////////////////  
 c         if(nplx.eq.6) angtemp = -1. * ax  
 c         if(nplx.eq.6) bfytemp = -1. * bfy  
          if(viewx.eq.12) angtemp = -1. * ax  
          if(viewx.eq.12) bfytemp = -1. * bfy  
          tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001  
          angx     = 180.*atan(tgtemp)/acos(-1.)  
          stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX  
 c$$$         print*,nplx,ax,bfy/10.  
 c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX  
 c$$$         print*,'========================'  
 c$$$         if(bfy.ne.0.)print*,viewx,'-x- '  
 c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ  
 *        --------------------------  
   
 c$$$         print*,'--- x-cl ---'  
 c$$$         istart = INDSTART(ICX)  
 c$$$         istop  = TOTCLLENGTH  
 c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1  
 c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)  
 c$$$         print*,(CLSIGNAL(i),i=istart,istop)  
 c$$$         print*,INDMAX(icx)  
 c$$$         print*,cog(4,icx)  
 c$$$         print*,fbad_cog(4,icx)  
755                    
756             call applypfa(PFAx,icx,angx,corr,res)
757           if(PFAx.eq.'COG1')then           stripx  = stripx + corr
758             resxPAM = res
             stripx  = stripx  
             resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM  
   
          elseif(PFAx.eq.'COG2')then  
   
             stripx  = stripx + cog(2,icx)              
             resxPAM = risx_cog(abs(angx))!TEMPORANEO                
             resxPAM = resxPAM*fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'COG3')then  
   
             stripx  = stripx + cog(3,icx)              
             resxPAM = risx_cog(abs(angx))!TEMPORANEO                        
             resxPAM = resxPAM*fbad_cog(3,icx)  
   
          elseif(PFAx.eq.'COG4')then  
   
             stripx  = stripx + cog(4,icx)              
             resxPAM = risx_cog(abs(angx))!TEMPORANEO                        
             resxPAM = resxPAM*fbad_cog(4,icx)  
   
          elseif(PFAx.eq.'ETA2')then  
   
             stripx  = stripx + pfaeta2(icx,angx)            
             resxPAM = risxeta2(abs(angx))  
             resxPAM = resxPAM*fbad_cog(2,icx)  
             if(DEBUG.and.fbad_cog(2,icx).ne.1)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'ETA3')then                          
   
             stripx  = stripx + pfaeta3(icx,angx)            
             resxPAM = risxeta3(abs(angx))                        
             resxPAM = resxPAM*fbad_cog(3,icx)                
 c            if(DEBUG.and.fbad_cog(3,icx).ne.1)              
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)  
   
          elseif(PFAx.eq.'ETA4')then                          
   
             stripx  = stripx + pfaeta4(icx,angx)              
             resxPAM = risxeta4(abs(angx))                        
             resxPAM = resxPAM*fbad_cog(4,icx)                
 c            if(DEBUG.and.fbad_cog(4,icx).ne.1)                
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)  
   
          elseif(PFAx.eq.'ETA')then    
   
             stripx  = stripx + pfaeta(icx,angx)              
 c            resxPAM = riseta(icx,angx)                      
             resxPAM = riseta(viewx,angx)                      
             resxPAM = resxPAM*fbad_eta(icx,angx)              
 c            if(DEBUG.and.fbad_cog(2,icx).ne.1)                
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'ETAL')then    
   
             stripx  = stripx + pfaetal(icx,angx)              
             resxPAM = riseta(viewx,angx)                      
             resxPAM = resxPAM*fbad_eta(icx,angx)              
 c            if(DEBUG.and.fbad_cog(2,icx).ne.1)                
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'COG')then            
   
             stripx  = stripx + cog(0,icx)              
             resxPAM = risx_cog(abs(angx))                      
             resxPAM = resxPAM*fbad_cog(0,icx)  
   
          else  
             if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx  
          endif  
   
   
 *     ======================================  
 *     temporary patch for saturated clusters  
 *     ======================================  
          if( nsatstrips(icx).gt.0 )then  
             stripx  = stripx + cog(4,icx)              
             resxPAM = pitchX*1e-4/sqrt(12.)  
             cc=cog(4,icx)  
 c$$$            print*,icx,' *** ',cc  
 c$$$            print*,icx,' *** ',resxPAM  
          endif  
759    
760   10   endif   10   endif
   
761            
762  *     -----------------  *     -----------------
763  *     CLUSTER Y  *     CLUSTER Y
# Line 760  c$$$            print*,icx,' *** ',resxP Line 768  c$$$            print*,icx,' *** ',resxP
768           viewy = VIEW(icy)           viewy = VIEW(icy)
769           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
770           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
771           resyPAM = RESYAV  c         resyPAM = RESYAV
772           stripy = float(MAXS(icy))           stripy = float(MAXS(icy))
773    
774           if(           if(
# Line 778  c$$$            print*,icx,' *** ',resxP Line 786  c$$$            print*,icx,' *** ',resxP
786           endif           endif
787    
788           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
789              if(DEBUG) then              if(DEBUG.EQ.1) then
790                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
791       $              ,icx,icy       $              ,icx,icy
792              endif              endif
793              goto 100              goto 100
794           endif           endif
795    
796  *        --------------------------  *        --------------------------
797  *        magnetic-field corrections  *        magnetic-field corrections
798  *        --------------------------  *        --------------------------
799           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001                   stripy  = stripy +  fieldcorr(viewy,bfx)      
800           angy    = 180.*atan(tgtemp)/acos(-1.)           angy    = effectiveangle(ay,viewy,bfx)
          stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY  
 c$$$         if(bfx.ne.0.)print*,viewy,'-y- '  
 c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ  
 *        --------------------------  
801                    
802  c$$$         print*,'--- y-cl ---'           call applypfa(PFAy,icy,angy,corr,res)
803  c$$$         istart = INDSTART(ICY)           stripy  = stripy + corr
804  c$$$         istop  = TOTCLLENGTH           resyPAM = res
 c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1  
 c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)  
 c$$$         print*,(CLSIGNAL(i),i=istart,istop)  
 c$$$         print*,INDMAX(icy)  
 c$$$         print*,cog(4,icy)  
 c$$$         print*,fbad_cog(4,icy)  
   
          if(PFAy.eq.'COG1')then  
   
             stripy  = stripy      
             resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM    
   
          elseif(PFAy.eq.'COG2')then  
   
             stripy  = stripy + cog(2,icy)  
             resyPAM = risy_cog(abs(angy))!TEMPORANEO  
             resyPAM = resyPAM*fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'COG3')then  
   
             stripy  = stripy + cog(3,icy)  
             resyPAM = risy_cog(abs(angy))!TEMPORANEO  
             resyPAM = resyPAM*fbad_cog(3,icy)  
   
          elseif(PFAy.eq.'COG4')then  
   
             stripy  = stripy + cog(4,icy)  
             resyPAM = risy_cog(abs(angy))!TEMPORANEO  
             resyPAM = resyPAM*fbad_cog(4,icy)  
   
          elseif(PFAy.eq.'ETA2')then  
   
             stripy  = stripy + pfaeta2(icy,angy)            
             resyPAM = risyeta2(abs(angy))                
             resyPAM = resyPAM*fbad_cog(2,icy)  
 c            if(DEBUG.and.fbad_cog(2,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'ETA3')then                        
   
             stripy  = stripy + pfaeta3(icy,angy)  
             resyPAM = resyPAM*fbad_cog(3,icy)    
 c            if(DEBUG.and.fbad_cog(3,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)  
   
          elseif(PFAy.eq.'ETA4')then    
   
             stripy  = stripy + pfaeta4(icy,angy)  
             resyPAM = resyPAM*fbad_cog(4,icy)  
 c            if(DEBUG.and.fbad_cog(4,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)  
   
          elseif(PFAy.eq.'ETA')then  
   
             stripy  = stripy + pfaeta(icy,angy)  
 c            resyPAM = riseta(icy,angy)    
             resyPAM = riseta(viewy,angy)    
             resyPAM = resyPAM*fbad_eta(icy,angy)  
 c            if(DEBUG.and.fbad_cog(2,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'ETAL')then  
   
             stripy  = stripy + pfaetal(icy,angy)  
             resyPAM = riseta(viewy,angy)    
             resyPAM = resyPAM*fbad_eta(icy,angy)  
 c            if(DEBUG.and.fbad_cog(2,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'COG')then  
   
             stripy  = stripy + cog(0,icy)              
             resyPAM = risy_cog(abs(angy))  
             resyPAM = resyPAM*fbad_cog(0,icy)  
   
          else  
             if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx  
          endif  
   
   
 *     ======================================  
 *     temporary patch for saturated clusters  
 *     ======================================  
          if( nsatstrips(icy).gt.0 )then  
             stripy  = stripy + cog(4,icy)              
             resyPAM = pitchY*1e-4/sqrt(12.)  
             cc=cog(4,icy)  
 c$$$            print*,icy,' *** ',cc  
 c$$$            print*,icy,' *** ',resyPAM  
          endif  
   
805    
806   20   endif   20   endif
807    
808  c$$$      print*,'## stripx,stripy ',stripx,stripy  cc      print*,'## stripx,stripy ',stripx,stripy
809    
810  c===========================================================  c===========================================================
811  C     COUPLE  C     COUPLE
# Line 903  c     (xi,yi,zi) = mechanical coordinate Line 817  c     (xi,yi,zi) = mechanical coordinate
817  c------------------------------------------------------------------------  c------------------------------------------------------------------------
818           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
819       $        .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...
820              if(DEBUG) then              if(DEBUG.EQ.1) then
821                 print*,'xyz_PAM (couple):',                 print*,'xyz_PAM (couple):',
822       $              ' WARNING: false X strip: strip ',stripx       $              ' WARNING: false X strip: strip ',stripx
823              endif              endif
824           endif           endif
825           xi = acoordsi(stripx,viewx)           xi = dcoordsi(stripx,viewx)
826           yi = acoordsi(stripy,viewy)           yi = dcoordsi(stripy,viewy)
827           zi = 0.           zi = 0.D0
828                    
   
829  c------------------------------------------------------------------------  c------------------------------------------------------------------------
830  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
831  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 946  c--------------------------------------- Line 859  c---------------------------------------
859           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
860           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
861    
862           xPAM_A = 0.           xPAM_A = 0.D0
863           yPAM_A = 0.           yPAM_A = 0.D0
864           zPAM_A = 0.           zPAM_A = 0.D0
865    
866           xPAM_B = 0.           xPAM_B = 0.D0
867           yPAM_B = 0.           yPAM_B = 0.D0
868           zPAM_B = 0.           zPAM_B = 0.D0
869    
870        elseif(        elseif(
871       $        (icx.ne.0.and.icy.eq.0).or.       $        (icx.ne.0.and.icy.eq.0).or.
# Line 972  C======================================= Line 885  C=======================================
885              nldx = nldy              nldx = nldy
886              viewx = viewy + 1              viewx = viewy + 1
887    
888              yi   = acoordsi(stripy,viewy)              yi   = dcoordsi(stripy,viewy)
889    
890              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
891              yi_A = yi              yi_A = yi
# Line 998  c            print*,'X-singlet ',icx,npl Line 911  c            print*,'X-singlet ',icx,npl
911  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...
912              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
913       $           .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...
914                 if(DEBUG) then                 if(DEBUG.EQ.1) then
915                    print*,'xyz_PAM (X-singlet):',                    print*,'xyz_PAM (X-singlet):',
916       $                 ' WARNING: false X strip: strip ',stripx       $                 ' WARNING: false X strip: strip ',stripx
917                 endif                 endif
918              endif              endif
919              xi   = acoordsi(stripx,viewx)              xi   = dcoordsi(stripx,viewx)
920    
921              xi_A = xi              xi_A = xi
922              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 1023  c            print*,'X-cl ',icx,stripx,' Line 936  c            print*,'X-cl ',icx,stripx,'
936  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
937    
938           else           else
939              if(DEBUG) then              if(DEBUG.EQ.1) then
940                 print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
941                 print *,'icx = ',icx                 print *,'icx = ',icx
942                 print *,'icy = ',icy                 print *,'icy = ',icy
# Line 1076  c     (xPAM,yPAM,zPAM) = measured coordi Line 989  c     (xPAM,yPAM,zPAM) = measured coordi
989  c                        in PAMELA reference system  c                        in PAMELA reference system
990  c------------------------------------------------------------------------  c------------------------------------------------------------------------
991    
992           xPAM = 0.           xPAM = 0.D0
993           yPAM = 0.           yPAM = 0.D0
994           zPAM = 0.           zPAM = 0.D0
995    
996           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
997           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 1092  c--------------------------------------- Line 1005  c---------------------------------------
1005  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1006    
1007        else        else
1008           if(DEBUG) then           if(DEBUG.EQ.1) then
1009              print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1010              print *,'icx = ',icx              print *,'icx = ',icx
1011              print *,'icy = ',icy              print *,'icy = ',icy
# Line 1137  c$$$      common/FINALPFA/PFA Line 1050  c$$$      common/FINALPFA/PFA
1050  c$$$      PFAx = 'COG4'!PFA  c$$$      PFAx = 'COG4'!PFA
1051  c$$$      PFAy = 'COG4'!PFA  c$$$      PFAy = 'COG4'!PFA
1052    
1053    
1054        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1055              print*,'xyzpam: ***WARNING*** clusters ',icx,icy              print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1056       $           ,' does not exists (nclstr1=',nclstr1,')'       $           ,' do not exists (n.clusters=',nclstr1,')'
1057              icx = -1*icx              icx = -1*icx
1058              icy = -1*icy              icy = -1*icy
1059              return              return
# Line 1149  c$$$      PFAy = 'COG4'!PFA Line 1063  c$$$      PFAy = 'COG4'!PFA
1063        call idtoc(pfaid,PFAx)        call idtoc(pfaid,PFAx)
1064        call idtoc(pfaid,PFAy)        call idtoc(pfaid,PFAy)
1065    
1066    cc      print*,PFAx,PFAy
1067    
1068  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)
1069    
1070  c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy  c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
# Line 1184  c$$$     $        ,' does not belong to Line 1100  c$$$     $        ,' does not belong to
1100           xm(ip) = xPAM           xm(ip) = xPAM
1101           ym(ip) = yPAM           ym(ip) = yPAM
1102           zm(ip) = zPAM           zm(ip) = zPAM
1103           xm_A(ip) = 0.           xm_A(ip) = 0.D0
1104           ym_A(ip) = 0.           ym_A(ip) = 0.D0
1105           xm_B(ip) = 0.           xm_B(ip) = 0.D0
1106           ym_B(ip) = 0.           ym_B(ip) = 0.D0
1107    
1108  c         zv(ip) = zPAM  c         zv(ip) = zPAM
1109    
# Line 1277  c         zv(ip) = z_mech_sensor(nplanes Line 1193  c         zv(ip) = z_mech_sensor(nplanes
1193    
1194        endif        endif
1195    
1196        if(DEBUG)then        if(DEBUG.EQ.1)then
1197  c         print*,'----------------------------- track coord'  c         print*,'----------------------------- track coord'
1198  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1199           write(*,22222)ip,zm(ip),xm(ip),ym(ip)           write(*,22222)ip,zm(ip),xm(ip),ym(ip)
# Line 1307  c$$$         print*,'------------------- Line 1223  c$$$         print*,'-------------------
1223  *      *    
1224  ********************************************************************************  ********************************************************************************
1225    
1226        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1227    
1228        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1229    
# Line 1316  c$$$         print*,'------------------- Line 1232  c$$$         print*,'-------------------
1232  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1233  *     -----------------------------------  *     -----------------------------------
1234    
1235          real rXPP,rYPP
1236          double precision XPP,YPP
1237        double precision distance,RE        double precision distance,RE
1238        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1239    
1240          XPP=DBLE(rXPP)
1241          YPP=DBLE(rYPP)
1242    
1243  *     ----------------------  *     ----------------------
1244        if (        if (
1245       +     xPAM.eq.0.and.       +     xPAM.eq.0.and.
# Line 1445  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1366  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1366        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1367        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1368        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1369        real*8 yvvv,xvvv        double precision yvvv,xvvv
1370        double precision xi,yi,zi        double precision xi,yi,zi
1371        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1372        real AA,BB        real AA,BB
1373        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1374    
1375  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1376        real ptoll        real ptoll
1377        data ptoll/150./          !um        data ptoll/150./          !um
1378    
1379        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1380    
1381        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1382        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1470  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1391  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1391  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1392  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1393  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1394                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)  c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1395       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...  c     $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1396  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...
1397  c                  print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1398  c     $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1399                 endif  c               endif
1400                 xi = acoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1401                 yi = acoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1402                 zi = 0.                 zi = 0.D0
1403  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1404  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1405  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1728  c      include 'level1.f' Line 1649  c      include 'level1.f'
1649        integer iflag        integer iflag
1650    
1651        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1652        
1653          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1654    
1655  *     init variables  *     init variables
1656        ncp_tot=0        ncp_tot=0
# Line 1760  c      include 'level1.f' Line 1683  c      include 'level1.f'
1683           if( ncl_view(iv).gt. nclusterlimit)then           if( ncl_view(iv).gt. nclusterlimit)then
1684  c            mask_view(iv) = 1  c            mask_view(iv) = 1
1685              mask_view(iv) = mask_view(iv) + 2**0              mask_view(iv) = mask_view(iv) + 2**0
1686              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG.EQ.1)
1687       $           ,nclusterlimit,' on view ', iv,' --> masked!'       $        print*,' * WARNING * cl_to_couple: n.clusters > '
1688         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1689           endif           endif
1690        enddo        enddo
1691    
# Line 1899  c                  cut = chcut * sch(npl Line 1823  c                  cut = chcut * sch(npl
1823                 endif                 endif
1824    
1825                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1826                    if(verbose)print*,                    if(verbose.eq.1)print*,
1827       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1828       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1829       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
# Line 1937  c                  mask_view(nviewy(nply Line 1861  c                  mask_view(nviewy(nply
1861        enddo        enddo
1862                
1863                
1864        if(DEBUG)then        if(DEBUG.EQ.1)then
1865           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1866           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1867           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
1868           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1869        endif        endif
1870                
# Line 2000  c      double precision xm3,ym3,zm3 Line 1924  c      double precision xm3,ym3,zm3
1924        real xc,zc,radius        real xc,zc,radius
1925  *     -----------------------------  *     -----------------------------
1926    
1927          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1928    
1929  *     --------------------------------------------  *     --------------------------------------------
1930  *     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 1962  c     print*,'***',is1,xm1,ym1,zm1
1962                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1963                    if(  mask_view(nviewx(ip2)).ne.0 .or.                    if(  mask_view(nviewx(ip2)).ne.0 .or.
1964       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1965                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
                       
1966                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1967                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1968                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
# Line 2057  c     $                       (icx2,icy2 Line 1981  c     $                       (icx2,icy2
1981  *     (2 couples needed)  *     (2 couples needed)
1982  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1983                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
1984                             if(verbose)print*,                             if(verbose.eq.1)print*,
1985       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
1986       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
1987       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
# Line 2121  c     $                               (i Line 2045  c     $                               (i
2045                                   ym3=yPAM                                   ym3=yPAM
2046                                   zm3=zPAM                                   zm3=zPAM
2047  *     find the circle passing through the three points  *     find the circle passing through the three points
2048    c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
2049    c$$$     $                                ,xc,zc,radius,iflag)
2050                                     iflag = DEBUG
2051                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2052       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag)
2053  c     print*,xc,zc,radius  c     print*,xc,zc,radius
# Line 2137  c     $                                 Line 2064  c     $                                
2064  *     (3 couples needed)  *     (3 couples needed)
2065  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2066                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2067                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2068       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2069       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2070       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
# Line 2199  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2126  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2126   10   continue   10   continue
2127        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2128                
2129        if(DEBUG)then        if(DEBUG.EQ.1)then
2130           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2131           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2132        endif        endif
# Line 2246  c      include 'momanhough_init.f' Line 2173  c      include 'momanhough_init.f'
2173        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2174        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2175    
2176          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2177    
2178  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2179  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2372  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2300  c         if(ncpused.lt.ncpyz_min)goto 2
2300  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2301    
2302           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2303              if(verbose)print*,              if(verbose.eq.1)print*,
2304       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2305       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2306       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
# Line 2400  c     ptcloud_yz_nt(nclouds_yz)=npt Line 2328  c     ptcloud_yz_nt(nclouds_yz)=npt
2328  c     print*,'>> ',ipt,db_all(ipt)  c     print*,'>> ',ipt,db_all(ipt)
2329           enddo             enddo  
2330           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2331           if(DEBUG)then           if(DEBUG.EQ.1)then
2332              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2333              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2334              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2335              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2336              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2337              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2338              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2339         $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2340                print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
2341  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
2342  c$$$     $           ,ptcloud_yz(nclouds_yz)  c$$$     $           ,ptcloud_yz(nclouds_yz)
2343  c$$$            print*,'nt-uple: db_cloud(...) = '  c$$$            print*,'nt-uple: db_cloud(...) = '
# Line 2425  c$$$     $           ,(db_cloud(iii),iii Line 2355  c$$$     $           ,(db_cloud(iii),iii
2355          goto 90                          goto 90                
2356        endif                            endif                    
2357                
2358        if(DEBUG)then        if(DEBUG.EQ.1)then
2359           print*,'---------------------- '           print*,'---------------------- '
2360           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
2361           print*,' '           print*,' '
# Line 2474  c      include 'momanhough_init.f' Line 2404  c      include 'momanhough_init.f'
2404        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2405        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2406    
2407          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2408    
2409  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2410  *     classification of TRIPLETS  *     classification of TRIPLETS
2411  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2531  c         tr_temp(1)=itr1 Line 2463  c         tr_temp(1)=itr1
2463       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2464                 distance = sqrt(distance)                 distance = sqrt(distance)
2465                                
2466                 if(distance.lt.cutdistxz)then  *     ------------------------------------------------------------------------
2467    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2468    *     ------------------------------------------------------------------------
2469    *     (added in august 2007)
2470                   istrimage=0
2471                   if(
2472         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2473         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2474         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2475         $              .true.)istrimage=1
2476    
2477                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2478  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
2479                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2480                    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 2538  c         if(ncpused.lt.ncpxz_min)goto 2
2538  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2539  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2540           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2541              if(verbose)print*,              if(verbose.eq.1)print*,
2542       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2543       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2544       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
# Line 2622  c               mask_view(iv) = 6 Line 2565  c               mask_view(iv) = 6
2565           enddo           enddo
2566           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2567                    
2568           if(DEBUG)then           if(DEBUG.EQ.1)then
2569              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2570              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2571              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2572              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2573              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2574              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2575              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2576                print*,'cpcloud_xz '
2577         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2578              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2579  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2580  c$$$     $           ,ptcloud_xz(nclouds_xz)  c$$$     $           ,ptcloud_xz(nclouds_xz)
# Line 2647  c$$$     $           ,(tr_cloud(iii),iii Line 2592  c$$$     $           ,(tr_cloud(iii),iii
2592           goto 91                           goto 91                
2593         endif                             endif                    
2594                
2595        if(DEBUG)then        if(DEBUG.EQ.1)then
2596           print*,'---------------------- '           print*,'---------------------- '
2597           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
2598           print*,' '           print*,' '
# Line 2699  c$$$     $           ,(tr_cloud(iii),iii Line 2644  c$$$     $           ,(tr_cloud(iii),iii
2644  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2645  *     variables for track fitting  *     variables for track fitting
2646        double precision AL_INI(5)        double precision AL_INI(5)
 c      double precision tath  
2647  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2648    
2649          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2650    
2651    
2652        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
# Line 2724  c      real fitz(nplanes)        !z coor Line 2668  c      real fitz(nplanes)        !z coor
2668                 enddo                 enddo
2669              enddo              enddo
2670              ncp_ok=0              ncp_ok=0
2671              do icp=1,ncp_tot    !loop on couples              do icp=1,ncp_tot    !loop over couples
2672  *     get info on  *     get info on
2673                 cpintersec(icp)=min(                 cpintersec(icp)=min(
2674       $              cpcloud_yz(iyz,icp),       $              cpcloud_yz(iyz,icp),
# Line 2733  c      real fitz(nplanes)        !z coor Line 2677  c      real fitz(nplanes)        !z coor
2677       $    (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.
2678       $    (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.
2679       $              .false.)cpintersec(icp)=0       $              .false.)cpintersec(icp)=0
2680    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2681                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
2682                    ncp_ok=ncp_ok+1                      ncp_ok=ncp_ok+1  
2683                                        
# Line 2765  c      real fitz(nplanes)        !z coor Line 2710  c      real fitz(nplanes)        !z coor
2710                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2711              enddo              enddo
2712                            
 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  
2713                            
2714              if(DEBUG)then              if(DEBUG.EQ.1)then
2715                 print*,'Combination ',iyz,ixz                 print*,'Combination ',iyz,ixz
2716       $              ,' db ',ptcloud_yz(iyz)       $              ,' db ',ptcloud_yz(iyz)
2717       $              ,' tr ',ptcloud_xz(ixz)       $              ,' tr ',ptcloud_xz(ixz)
2718       $              ,'  -----> # matching couples ',ncp_ok       $              ,'  -----> # matching couples ',ncp_ok
2719              endif              endif
2720    
2721    c            if(nplused.lt.nplxz_min)goto 888 !next combination
2722                if(nplused.lt.nplyz_min)goto 888 !next combination
2723                if(ncp_ok.lt.ncpok)goto 888 !next combination
2724    
2725  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2726  c$$$  print*,'Configurazione cluster XZ'  c$$$  print*,'Configurazione cluster XZ'
2727  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 2767  c$$$            AL_INI(5) = (1.e2*alfaxz
2767  c$$$              c$$$            
2768  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2769                                                    
2770              if(DEBUG)then              if(DEBUG.EQ.1)then
2771                   print*,'track candidate', ntracks+1
2772                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2773                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2774                 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 2801  c$$$            if(AL_INI(5).gt.defmax)g
2801                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2802                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2803                                                                
2804                                  *                             ---------------------------------------
2805    *                             check if this group of couples has been
2806    *                             already fitted    
2807    *                             ---------------------------------------
2808                                  do ica=1,ntracks
2809                                     isthesame=1
2810                                     do ip=1,NPLANES
2811                                        if(hit_plane(ip).ne.0)then
2812                                           if(  CP_STORE(nplanes-ip+1,ica)
2813         $                                      .ne.
2814         $                                      cp_match(ip,hit_plane(ip)) )
2815         $                                      isthesame=0
2816                                        else
2817                                           if(  CP_STORE(nplanes-ip+1,ica)
2818         $                                      .ne.
2819         $                                      0 )
2820         $                                      isthesame=0
2821                                        endif
2822                                     enddo
2823                                     if(isthesame.eq.1)then
2824                                        if(DEBUG.eq.1)
2825         $                                   print*,'(already fitted)'
2826                                        goto 666 !jump to next combination
2827                                     endif
2828                                  enddo
2829    
2830                                call track_init !init TRACK common                                call track_init !init TRACK common
2831    
2832                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2833                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2834                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2835                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2900  c$$$                              enddo Line 2873  c$$$                              enddo
2873                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2874                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2875                                iprint=0                                iprint=0
2876  c                              if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
2877                                if(DEBUG)iprint=2                                if(DEBUG.EQ.1)iprint=2
2878                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2879                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2880                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
2881                                      print *,                                      print *,
2882       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2883       $                              //'(clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
# Line 2929  c                                 chi2=- Line 2902  c                                 chi2=-
2902  *     --------------------------  *     --------------------------
2903                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2904                                                                    
2905                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
2906       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2907       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2908       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
# Line 2945  c                                    mas Line 2918  c                                    mas
2918                                                                
2919                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2920                                                                
2921                                do ip=1,nplanes                                do ip=1,nplanes !top to bottom
2922    
2923                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2924                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
# Line 2962  c                                    mas Line 2935  c                                    mas
2935                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
2936                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
2937                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
2938    *                                NB! hit_plane is defined from bottom to top
2939                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2940                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2941       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
# Line 2977  c                                    mas Line 2951  c                                    mas
2951                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0
2952                                      LADDER_STORE(nplanes-ip+1,ntracks)=0                                      LADDER_STORE(nplanes-ip+1,ntracks)=0
2953                                   endif                                   endif
2954                                   BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now                                   BX_STORE(ip,ntracks)=0!I dont need it now
2955                                   BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now                                   BY_STORE(ip,ntracks)=0!I dont need it now
2956                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(ip,ntracks)=0
2957                                   do i=1,5                                   do i=1,5
2958                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2959                                   enddo                                   enddo
# Line 3008  c                                    mas Line 2982  c                                    mas
2982           return           return
2983        endif        endif
2984                
2985  c$$$      if(DEBUG)then  c$$$      if(DEBUG.EQ.1)then
2986  c$$$         print*,'****** TRACK CANDIDATES ***********'  c$$$         print*,'****** TRACK CANDIDATES ***********'
2987  c$$$         print*,'#         R. chi2        RIG'  c$$$         print*,'#         R. chi2        RIG'
2988  c$$$         do i=1,ntracks  c$$$         do i=1,ntracks
# Line 3017  c$$$     $           ,1./abs(AL_STORE(5, Line 2991  c$$$     $           ,1./abs(AL_STORE(5,
2991  c$$$         enddo  c$$$         enddo
2992  c$$$         print*,'***********************************'  c$$$         print*,'***********************************'
2993  c$$$      endif  c$$$      endif
2994        if(DEBUG)then        if(DEBUG.EQ.1)then
2995          print*,'****** TRACK CANDIDATES *****************'          print*,'****** TRACK CANDIDATES *****************'
2996          print*,'#         R. chi2        RIG         ndof'          print*,'#         R. chi2        RIG         ndof'
2997          do i=1,ntracks          do i=1,ntracks
# Line 3069  c$$$      endif Line 3043  c$$$      endif
3043        real xyzp(3),bxyz(3)        real xyzp(3),bxyz(3)
3044        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3045    
3046          if(DEBUG.EQ.1)print*,'refine_track:'
3047  *     =================================================  *     =================================================
3048  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3049  *                          and  *                          and
# Line 3077  c$$$      endif Line 3052  c$$$      endif
3052        call track_init        call track_init
3053        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3054    
3055             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3056    
3057           xP=XV_STORE(nplanes-ip+1,ibest)           xP=XV_STORE(nplanes-ip+1,ibest)
3058           yP=YV_STORE(nplanes-ip+1,ibest)           yP=YV_STORE(nplanes-ip+1,ibest)
3059           zP=ZV_STORE(nplanes-ip+1,ibest)           zP=ZV_STORE(nplanes-ip+1,ibest)
# Line 3150  c     $           AYV_STORE(nplanes-ip+1 Line 3127  c     $           AYV_STORE(nplanes-ip+1
3127              LADDER_STORE(nplanes-ip+1,IBEST)=nldt              LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3128  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3129    
3130              if(DEBUG)then              if(DEBUG.EQ.1)then
3131                 print*,                 print*,
3132       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3133       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3161  c     $           AYV_STORE(nplanes-ip+1 Line 3138  c     $           AYV_STORE(nplanes-ip+1
3138  *     ===========================================  *     ===========================================
3139  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3140  *     ===========================================  *     ===========================================
3141  c            if(DEBUG)print*,'>>>> try to include a new couple'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
3142              distmin=1000000.              distmin=1000000.
3143              xmm = 0.              xmm = 0.
3144              ymm = 0.              ymm = 0.
# Line 3178  c            if(DEBUG)print*,'>>>> try t Line 3155  c            if(DEBUG)print*,'>>>> try t
3155                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3156  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
3157  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
3158       $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)       $              cl_used(icx).ne.0.or. !or the X cluster is already used
3159       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)       $              cl_used(icy).ne.0.or. !or the Y cluster is already used
3160       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3161  *            *          
3162                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3193  c     $              cl_used(icy).eq.1.o Line 3170  c     $              cl_used(icy).eq.1.o
3170                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3171  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3172                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3173                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)print*,'( couple ',id
3174       $              ,' ) distance ',distance       $              ,' ) distance ',distance
3175                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3176                    xmm = xPAM                    xmm = xPAM
# Line 3225  c            if(distmin.le.clinc)then   Line 3202  c            if(distmin.le.clinc)then  
3202                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3203  *              -----------------------------------  *              -----------------------------------
3204                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3205                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3206       $              ,' (dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3207                 goto 133         !next plane                 goto 133         !next plane
3208              endif              endif
# Line 3233  c            if(distmin.le.clinc)then   Line 3210  c            if(distmin.le.clinc)then  
3210  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3211  *                     either from a couple or single  *                     either from a couple or single
3212  *     ================================================  *     ================================================
3213  c            if(DEBUG)print*,'>>>> try to include a new cluster'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
3214              distmin=1000000.              distmin=1000000.
3215              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3216              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3270  c     $              AXV_STORE(nplanes-i Line 3247  c     $              AXV_STORE(nplanes-i
3247       $              )                     $              )              
3248                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3249  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3250                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)print*,'( cl-X ',icx
3251       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3252                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3253                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3303  c     $              0.,AYV_STORE(nplane Line 3280  c     $              0.,AYV_STORE(nplane
3280       $              )       $              )
3281                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3282  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3283                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy
3284       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3285                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3286                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3325  c                 dedxmm = sgnl(icy)  !( Line 3302  c                 dedxmm = sgnl(icy)  !(
3302  *----- single clusters -----------------------------------------------    *----- single clusters -----------------------------------------------  
3303  c            print*,'## ncls(',ip,') ',ncls(ip)  c            print*,'## ncls(',ip,') ',ncls(ip)
3304              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3305    c               print*,'-',ic,'-'
3306                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3307  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
3308                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
# Line 3348  c               if(cl_used(icl).eq.1.or. Line 3326  c               if(cl_used(icl).eq.1.or.
3326    
3327                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3328  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3329                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)print*,'( cl-s ',icl
3330       $              ,' ) distance ',distance,'<',distmin,' ?'       $              ,' ) distance ',distance
3331                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3332                    if(DEBUG)print*,'YES'  c                  if(DEBUG.EQ.1)print*,'YES'
3333                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3334                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3335                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3397  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3375  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3375                    if(mod(VIEW(iclm),2).eq.0)then                    if(mod(VIEW(iclm),2).eq.0)then
3376                       XGOOD(nplanes-ip+1)=1.                       XGOOD(nplanes-ip+1)=1.
3377                       resx(nplanes-ip+1)=rxmm                       resx(nplanes-ip+1)=rxmm
3378                       if(DEBUG)print*,'%%%% included X-cl ',iclm                       if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
3379       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3380       $                    ,', dist.= ',distmin       $                    ,', dist.= ',distmin
3381       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clinc,' )'
3382                    else                    else
3383                       YGOOD(nplanes-ip+1)=1.                       YGOOD(nplanes-ip+1)=1.
3384                       resy(nplanes-ip+1)=rymm                       resy(nplanes-ip+1)=rymm
3385                       if(DEBUG)print*,'%%%% included Y-cl ',iclm                       if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
3386       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3387       $                    ,', dist.= ', distmin       $                    ,', dist.= ', distmin
3388       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clinc,' )'
# Line 3446  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3424  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3424        include 'common_momanhough.f'        include 'common_momanhough.f'
3425        include 'level2.f'              include 'level2.f'      
3426    
3427          if(DEBUG.EQ.1)print*,'clean_XYclouds:'
3428    
3429        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3430    
3431           id=CP_STORE(nplanes-ip+1,ibest)           id=CP_STORE(nplanes-ip+1,ibest)
# Line 3454  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3434  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3434              if(id.ne.0)then              if(id.ne.0)then
3435                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3436                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3437                 cl_used(iclx)=ntrk  !tag used clusters  c$$$               cl_used(iclx)=ntrk  !tag used clusters
3438                 cl_used(icly)=ntrk  !tag used clusters  c$$$               cl_used(icly)=ntrk  !tag used clusters
3439              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3440                 cl_used(icl)=ntrk   !tag used clusters  c$$$               cl_used(icl)=ntrk   !tag used clusters
3441              endif              endif
3442                            
3443  *     -----------------------------  *     -----------------------------
# Line 3476  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3456  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3456       $              cly(ip,icp).eq.icl       $              cly(ip,icp).eq.icl
3457       $              )then       $              )then
3458                    id=id_cp(ip,icp,1)                    id=id_cp(ip,icp,1)
3459                    if(DEBUG)then                    if(DEBUG.EQ.1)then
3460                       print*,ip,' <<< cp ',id                       print*,ip,' <<< cp ',id
3461       $                    ,' ( cl-x '       $                    ,' ( cl-x '
3462       $                    ,clx(ip,icp)       $                    ,clx(ip,icp)
# Line 3557  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3537  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3537              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3538              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3539              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3540                multmaxx(ip,it) = 0
3541                seedx(ip,it)    = 0
3542                xpu(ip,it)      = 0
3543                multmaxy(ip,it) = 0
3544                seedy(ip,it)    = 0
3545                ypu(ip,it)      = 0
3546           enddo           enddo
3547           do ipa=1,5           do ipa=1,5
3548              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3576  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3562  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3562          ys(1,ip)=0          ys(1,ip)=0
3563          ys(2,ip)=0          ys(2,ip)=0
3564          sgnlys(ip)=0          sgnlys(ip)=0
3565            sxbad(ip)=0
3566            sybad(ip)=0
3567            multmaxsx(ip)=0
3568            multmaxsy(ip)=0
3569        enddo        enddo
3570        end        end
3571    
# Line 3736  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3726  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3726           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3727           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3728           ayv_nt(ip,ntr)   = sngl(ayv(ip))             ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3729  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
3730           factor = sqrt(           factor = sqrt(
3731       $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3732       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3733       $        1. )       $        1. )
3734  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
3735           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3736           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3737        
# Line 3749  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Line 3739  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3739           ay   = ayv_nt(ip,ntr)           ay   = ayv_nt(ip,ntr)
3740           bfx  = BX_STORE(ip,IDCAND)           bfx  = BX_STORE(ip,IDCAND)
3741           bfy  = BY_STORE(ip,IDCAND)           bfy  = BY_STORE(ip,IDCAND)
3742           if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)  c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3743           if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)  c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3744           tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001  c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3745           angx     = 180.*atan(tgtemp)/acos(-1.)  c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3746           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001          c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3747           angy    = 180.*atan(tgtemp)/acos(-1.)  c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3748    
3749             angx = effectiveangle(ax,2*ip,bfy)
3750             angy = effectiveangle(ay,2*ip-1,bfx)
3751            
3752                    
3753  c         print*,'* ',ip,bfx,bfy,angx,angy  c         print*,'* ',ip,bfx,bfy,angx,angy
3754    
# Line 3771  c         print*,'* ',ip,bfx,bfy,angx,an Line 3765  c         print*,'* ',ip,bfx,bfy,angx,an
3765  c           >>> is a couple  c           >>> is a couple
3766              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3767              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3768                
3769  c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)              cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3770  c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)                          cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3771  c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))  
 c$$$            ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))  
3772              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3773              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3774    
# Line 3785  c$$$            ybad(ip,ntr)= nbadstrips Line 3778  c$$$            ybad(ip,ntr)= nbadstrips
3778              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3779       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3780    
3781                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3782         $                         +10000*mult(cltrx(ip,ntr))
3783                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3784         $           /clsigma(indmax(cltrx(ip,ntr)))
3785                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3786                xpu(ip,ntr)      = corr
3787    
3788                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3789         $                         +10000*mult(cltry(ip,ntr))
3790                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3791         $           /clsigma(indmax(cltry(ip,ntr)))
3792                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3793                ypu(ip,ntr)      = corr
3794    
3795           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3796    
3797    
3798                cl_used(icl) = 1    !tag used clusters          
3799    
3800              if(mod(VIEW(icl),2).eq.0)then              if(mod(VIEW(icl),2).eq.0)then
3801                 cltrx(ip,ntr)=icl                 cltrx(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angx)  
 c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)  
3802                 xbad(ip,ntr) = nbadstrips(4,icl)                 xbad(ip,ntr) = nbadstrips(4,icl)
3803    
3804                 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)                 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3805    
3806                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3807         $                         +10000*mult(cltrx(ip,ntr))
3808                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3809         $           /clsigma(indmax(cltrx(ip,ntr)))
3810                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3811                   xpu(ip,ntr)      = corr
3812    
3813              elseif(mod(VIEW(icl),2).eq.1)then              elseif(mod(VIEW(icl),2).eq.1)then
3814                 cltry(ip,ntr)=icl                 cltry(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angy)  
 c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)  
3815                 ybad(ip,ntr) = nbadstrips(4,icl)                 ybad(ip,ntr) = nbadstrips(4,icl)
3816    
3817                 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)                 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3818    
3819                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3820         $                         +10000*mult(cltry(ip,ntr))
3821                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3822         $           /clsigma(indmax(cltry(ip,ntr)))
3823                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3824                   ypu(ip,ntr)      = corr
3825                  
3826              endif              endif
3827    
3828           endif                     endif          
3829    
3830        enddo        enddo
3831    
3832          if(DEBUG.eq.1)then
3833             print*,'> STORING TRACK ',ntr
3834             print*,'clusters: '
3835             do ip=1,6
3836                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3837             enddo
3838          endif
3839    
3840  c$$$      print*,(xgood(i),i=1,6)  c$$$      print*,(xgood(i),i=1,6)
3841  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 3870  c         if( mask_view(iv).ne.0 )good2(
3870           good2(iv) = good2(iv) + mask_view(iv)*2**8           good2(iv) = good2(iv) + mask_view(iv)*2**8
3871        enddo        enddo
3872    
3873          if(DEBUG.eq.1)then
3874             print*,'> STORING SINGLETS '
3875          endif
3876    
3877        do icl=1,nclstr1        do icl=1,nclstr1
3878    
3879             ip=nplanes-npl(VIEW(icl))+1            
3880            
3881           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
3882              ip=nplanes-npl(VIEW(icl))+1              
3883              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3884    
3885                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3886                 planex(nclsx) = ip                 planex(nclsx) = ip
3887                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3888                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3889                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3890                   sxbad(nclsx)  = nbadstrips(1,icl)
3891                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
3892                  
3893    cc               print*,icl,' >>>> ',sxbad(nclsx)
3894    
3895                 do is=1,2                 do is=1,2
3896  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3897  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
# Line 3866  c$$$               print*,'xs(2,nclsx)   Line 3909  c$$$               print*,'xs(2,nclsx)  
3909                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3910                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3911                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3912                   sybad(nclsy)  = nbadstrips(1,icl)
3913                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
3914    
3915    cc               print*,icl,' >>>> ',sybad(nclsy)
3916    
3917                 do is=1,2                 do is=1,2
3918  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3919  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
# Line 3882  c$$$               print*,'ys(2,nclsy)   Line 3930  c$$$               print*,'ys(2,nclsy)  
3930    
3931  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3932           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
3933    *     --------------------------------------------------
3934    *     per non perdere la testa...
3935    *     whichtrack(icl) e` una variabile del common level1
3936    *     che serve solo per sapere quali cluster sono stati
3937    *     associati ad una traccia, e permettere di salvare
3938    *     solo questi nell'albero di uscita
3939    *     --------------------------------------------------
3940                    
3941        enddo        enddo
3942        end        end
3943    

Legend:
Removed from v.1.27  
changed lines
  Added in v.1.34

  ViewVC Help
Powered by ViewVC 1.1.23