/[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.36 by pam-fi, Tue Nov 25 14:41:38 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 1211  c$$$     $        ,' does not belong to Line 1127  c$$$     $        ,' does not belong to
1127           resx(ip)  = 1000.           resx(ip)  = 1000.
1128           resy(ip)  = resyPAM           resy(ip)  = resyPAM
1129    
1130    cPP --- old ---
1131    c$$$         xm(ip) = -100.
1132    c$$$         ym(ip) = -100.
1133    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1134    c$$$         xm_A(ip) = xPAM_A
1135    c$$$         ym_A(ip) = yPAM_A
1136    c$$$         xm_B(ip) = xPAM_B
1137    c$$$         ym_B(ip) = yPAM_B
1138    cPP --- new ---
1139           xm(ip) = -100.           xm(ip) = -100.
1140           ym(ip) = -100.           ym(ip) = -100.
1141           zm(ip) = (zPAM_A+zPAM_B)/2.           zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1142           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1143           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1144             zm_A(ip) = zPAM_A
1145           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
1146           ym_B(ip) = yPAM_B           ym_B(ip) = yPAM_B
1147             zm_B(ip) = zPAM_B
1148    cPP -----------
1149    
1150  c         zv(ip) = (zPAM_A+zPAM_B)/2.  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1151                    
# Line 1242  c$$$     $        ,' does not belong to Line 1170  c$$$     $        ,' does not belong to
1170           resx(ip)  = resxPAM           resx(ip)  = resxPAM
1171           resy(ip)  = 1000.           resy(ip)  = 1000.
1172    
1173    cPP --- old ---
1174    c$$$         xm(ip) = -100.
1175    c$$$         ym(ip) = -100.
1176    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1177    c$$$         xm_A(ip) = xPAM_A
1178    c$$$         ym_A(ip) = yPAM_A
1179    c$$$         xm_B(ip) = xPAM_B
1180    c$$$         ym_B(ip) = yPAM_B
1181    cPP --- new ---
1182           xm(ip) = -100.           xm(ip) = -100.
1183           ym(ip) = -100.           ym(ip) = -100.
1184           zm(ip) = (zPAM_A+zPAM_B)/2.           zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1185           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1186           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1187             zm_A(ip) = zPAM_A
1188           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
1189           ym_B(ip) = yPAM_B           ym_B(ip) = yPAM_B
1190                     zm_B(ip) = zPAM_B
1191    cPP -----------
1192    
1193  c         zv(ip) = (zPAM_A+zPAM_B)/2.  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1194    
1195        else        else
# Line 1277  c         zv(ip) = z_mech_sensor(nplanes Line 1217  c         zv(ip) = z_mech_sensor(nplanes
1217    
1218        endif        endif
1219    
1220        if(DEBUG)then        if(DEBUG.EQ.1)then
1221  c         print*,'----------------------------- track coord'  c         print*,'----------------------------- track coord'
1222  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1223           write(*,22222)ip,zm(ip),xm(ip),ym(ip)           write(*,22222)ip,zm(ip),xm(ip),ym(ip)
# Line 1307  c$$$         print*,'------------------- Line 1247  c$$$         print*,'-------------------
1247  *      *    
1248  ********************************************************************************  ********************************************************************************
1249    
1250        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1251    
1252        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1253    
# Line 1316  c$$$         print*,'------------------- Line 1256  c$$$         print*,'-------------------
1256  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1257  *     -----------------------------------  *     -----------------------------------
1258    
1259          real rXPP,rYPP
1260          double precision XPP,YPP
1261        double precision distance,RE        double precision distance,RE
1262        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1263    
1264          XPP=DBLE(rXPP)
1265          YPP=DBLE(rYPP)
1266    
1267  *     ----------------------  *     ----------------------
1268        if (        if (
1269       +     xPAM.eq.0.and.       +     xPAM.eq.0.and.
# Line 1445  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1390  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1390        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1391        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1392        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1393        real*8 yvvv,xvvv        double precision yvvv,xvvv
1394        double precision xi,yi,zi        double precision xi,yi,zi
1395        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1396        real AA,BB        real AA,BB
1397        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1398    
1399  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1400        real ptoll        real ptoll
1401        data ptoll/150./          !um        data ptoll/150./          !um
1402    
1403        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1404    
1405        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1406        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1470  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1415  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1415  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1416  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1417  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1418                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)  c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1419       $              .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...
1420  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...
1421  c                  print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1422  c     $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1423                 endif  c               endif
1424                 xi = acoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1425                 yi = acoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1426                 zi = 0.                 zi = 0.D0
1427  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1428  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1429  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1728  c      include 'level1.f' Line 1673  c      include 'level1.f'
1673        integer iflag        integer iflag
1674    
1675        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1676        
1677          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1678    
1679  *     init variables  *     init variables
1680        ncp_tot=0        ncp_tot=0
# Line 1760  c      include 'level1.f' Line 1707  c      include 'level1.f'
1707           if( ncl_view(iv).gt. nclusterlimit)then           if( ncl_view(iv).gt. nclusterlimit)then
1708  c            mask_view(iv) = 1  c            mask_view(iv) = 1
1709              mask_view(iv) = mask_view(iv) + 2**0              mask_view(iv) = mask_view(iv) + 2**0
1710              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG.EQ.1)
1711       $           ,nclusterlimit,' on view ', iv,' --> masked!'       $        print*,' * WARNING * cl_to_couple: n.clusters > '
1712         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1713           endif           endif
1714        enddo        enddo
1715    
# Line 1899  c                  cut = chcut * sch(npl Line 1847  c                  cut = chcut * sch(npl
1847                 endif                 endif
1848    
1849                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1850                    if(verbose)print*,                    if(verbose.eq.1)print*,
1851       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1852       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1853       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
# Line 1937  c                  mask_view(nviewy(nply Line 1885  c                  mask_view(nviewy(nply
1885        enddo        enddo
1886                
1887                
1888        if(DEBUG)then        if(DEBUG.EQ.1)then
1889           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1890           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1891           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
1892           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1893        endif        endif
1894                
# Line 2000  c      double precision xm3,ym3,zm3 Line 1948  c      double precision xm3,ym3,zm3
1948        real xc,zc,radius        real xc,zc,radius
1949  *     -----------------------------  *     -----------------------------
1950    
1951          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1952    
1953  *     --------------------------------------------  *     --------------------------------------------
1954  *     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 1986  c     print*,'***',is1,xm1,ym1,zm1
1986                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1987                    if(  mask_view(nviewx(ip2)).ne.0 .or.                    if(  mask_view(nviewx(ip2)).ne.0 .or.
1988       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1989                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
                       
1990                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1991                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1992                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
# Line 2057  c     $                       (icx2,icy2 Line 2005  c     $                       (icx2,icy2
2005  *     (2 couples needed)  *     (2 couples needed)
2006  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2007                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2008                             if(verbose)print*,                             if(verbose.eq.1)print*,
2009       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2010       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2011       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
# Line 2121  c     $                               (i Line 2069  c     $                               (i
2069                                   ym3=yPAM                                   ym3=yPAM
2070                                   zm3=zPAM                                   zm3=zPAM
2071  *     find the circle passing through the three points  *     find the circle passing through the three points
2072    c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
2073    c$$$     $                                ,xc,zc,radius,iflag)
2074                                     iflag = DEBUG
2075                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2076       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag)
2077  c     print*,xc,zc,radius  c     print*,xc,zc,radius
# Line 2137  c     $                                 Line 2088  c     $                                
2088  *     (3 couples needed)  *     (3 couples needed)
2089  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2090                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2091                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2092       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2093       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2094       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
# Line 2199  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2150  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2150   10   continue   10   continue
2151        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2152                
2153        if(DEBUG)then        if(DEBUG.EQ.1)then
2154           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2155           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2156        endif        endif
# Line 2246  c      include 'momanhough_init.f' Line 2197  c      include 'momanhough_init.f'
2197        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2198        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2199    
2200          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2201    
2202  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2203  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2372  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2324  c         if(ncpused.lt.ncpyz_min)goto 2
2324  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2325    
2326           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2327              if(verbose)print*,              if(verbose.eq.1)print*,
2328       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2329       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2330       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
# Line 2400  c     ptcloud_yz_nt(nclouds_yz)=npt Line 2352  c     ptcloud_yz_nt(nclouds_yz)=npt
2352  c     print*,'>> ',ipt,db_all(ipt)  c     print*,'>> ',ipt,db_all(ipt)
2353           enddo             enddo  
2354           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2355           if(DEBUG)then           if(DEBUG.EQ.1)then
2356              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2357              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2358              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2359              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2360              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2361              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2362              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2363         $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2364                print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
2365  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
2366  c$$$     $           ,ptcloud_yz(nclouds_yz)  c$$$     $           ,ptcloud_yz(nclouds_yz)
2367  c$$$            print*,'nt-uple: db_cloud(...) = '  c$$$            print*,'nt-uple: db_cloud(...) = '
# Line 2425  c$$$     $           ,(db_cloud(iii),iii Line 2379  c$$$     $           ,(db_cloud(iii),iii
2379          goto 90                          goto 90                
2380        endif                            endif                    
2381                
2382        if(DEBUG)then        if(DEBUG.EQ.1)then
2383           print*,'---------------------- '           print*,'---------------------- '
2384           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
2385           print*,' '           print*,' '
# Line 2474  c      include 'momanhough_init.f' Line 2428  c      include 'momanhough_init.f'
2428        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2429        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2430    
2431          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2432    
2433  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2434  *     classification of TRIPLETS  *     classification of TRIPLETS
2435  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2531  c         tr_temp(1)=itr1 Line 2487  c         tr_temp(1)=itr1
2487       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2488                 distance = sqrt(distance)                 distance = sqrt(distance)
2489                                
2490                 if(distance.lt.cutdistxz)then  *     ------------------------------------------------------------------------
2491    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2492    *     ------------------------------------------------------------------------
2493    *     (added in august 2007)
2494                   istrimage=0
2495                   if(
2496         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2497         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2498         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2499         $              .true.)istrimage=1
2500    
2501                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2502  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
2503                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2504                    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 2562  c         if(ncpused.lt.ncpxz_min)goto 2
2562  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2563  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2564           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2565              if(verbose)print*,              if(verbose.eq.1)print*,
2566       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2567       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2568       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
# Line 2622  c               mask_view(iv) = 6 Line 2589  c               mask_view(iv) = 6
2589           enddo           enddo
2590           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2591                    
2592           if(DEBUG)then           if(DEBUG.EQ.1)then
2593              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2594              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2595              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2596              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2597              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2598              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2599              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2600                print*,'cpcloud_xz '
2601         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2602              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2603  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2604  c$$$     $           ,ptcloud_xz(nclouds_xz)  c$$$     $           ,ptcloud_xz(nclouds_xz)
# Line 2647  c$$$     $           ,(tr_cloud(iii),iii Line 2616  c$$$     $           ,(tr_cloud(iii),iii
2616           goto 91                           goto 91                
2617         endif                             endif                    
2618                
2619        if(DEBUG)then        if(DEBUG.EQ.1)then
2620           print*,'---------------------- '           print*,'---------------------- '
2621           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
2622           print*,' '           print*,' '
# Line 2699  c$$$     $           ,(tr_cloud(iii),iii Line 2668  c$$$     $           ,(tr_cloud(iii),iii
2668  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2669  *     variables for track fitting  *     variables for track fitting
2670        double precision AL_INI(5)        double precision AL_INI(5)
 c      double precision tath  
2671  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2672    
2673          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2674    
2675    
2676        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
# Line 2724  c      real fitz(nplanes)        !z coor Line 2692  c      real fitz(nplanes)        !z coor
2692                 enddo                 enddo
2693              enddo              enddo
2694              ncp_ok=0              ncp_ok=0
2695              do icp=1,ncp_tot    !loop on couples              do icp=1,ncp_tot    !loop over couples
2696  *     get info on  *     get info on
2697                 cpintersec(icp)=min(                 cpintersec(icp)=min(
2698       $              cpcloud_yz(iyz,icp),       $              cpcloud_yz(iyz,icp),
# Line 2733  c      real fitz(nplanes)        !z coor Line 2701  c      real fitz(nplanes)        !z coor
2701       $    (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.
2702       $    (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.
2703       $              .false.)cpintersec(icp)=0       $              .false.)cpintersec(icp)=0
2704    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2705                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
2706                    ncp_ok=ncp_ok+1                      ncp_ok=ncp_ok+1  
2707                                        
# Line 2765  c      real fitz(nplanes)        !z coor Line 2734  c      real fitz(nplanes)        !z coor
2734                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2735              enddo              enddo
2736                            
 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  
2737                            
2738              if(DEBUG)then              if(DEBUG.EQ.1)then
2739                 print*,'Combination ',iyz,ixz                 print*,'Combination ',iyz,ixz
2740       $              ,' db ',ptcloud_yz(iyz)       $              ,' db ',ptcloud_yz(iyz)
2741       $              ,' tr ',ptcloud_xz(ixz)       $              ,' tr ',ptcloud_xz(ixz)
2742       $              ,'  -----> # matching couples ',ncp_ok       $              ,'  -----> # matching couples ',ncp_ok
2743              endif              endif
2744    
2745    c            if(nplused.lt.nplxz_min)goto 888 !next combination
2746                if(nplused.lt.nplyz_min)goto 888 !next combination
2747                if(ncp_ok.lt.ncpok)goto 888 !next combination
2748    
2749  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2750  c$$$  print*,'Configurazione cluster XZ'  c$$$  print*,'Configurazione cluster XZ'
2751  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 2791  c$$$            AL_INI(5) = (1.e2*alfaxz
2791  c$$$              c$$$            
2792  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2793                                                    
2794              if(DEBUG)then              if(DEBUG.EQ.1)then
2795                   print*,'track candidate', ntracks+1
2796                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2797                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2798                 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 2825  c$$$            if(AL_INI(5).gt.defmax)g
2825                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2826                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2827                                                                
2828                                  *                             ---------------------------------------
2829    *                             check if this group of couples has been
2830    *                             already fitted    
2831    *                             ---------------------------------------
2832                                  do ica=1,ntracks
2833                                     isthesame=1
2834                                     do ip=1,NPLANES
2835                                        if(hit_plane(ip).ne.0)then
2836                                           if(  CP_STORE(nplanes-ip+1,ica)
2837         $                                      .ne.
2838         $                                      cp_match(ip,hit_plane(ip)) )
2839         $                                      isthesame=0
2840                                        else
2841                                           if(  CP_STORE(nplanes-ip+1,ica)
2842         $                                      .ne.
2843         $                                      0 )
2844         $                                      isthesame=0
2845                                        endif
2846                                     enddo
2847                                     if(isthesame.eq.1)then
2848                                        if(DEBUG.eq.1)
2849         $                                   print*,'(already fitted)'
2850                                        goto 666 !jump to next combination
2851                                     endif
2852                                  enddo
2853    
2854                                call track_init !init TRACK common                                call track_init !init TRACK common
2855    
2856                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2857                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2858                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2859                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2900  c$$$                              enddo Line 2897  c$$$                              enddo
2897                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2898                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2899                                iprint=0                                iprint=0
2900  c                              if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
2901                                if(DEBUG)iprint=2                                if(DEBUG.EQ.1)iprint=2
2902                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2903                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2904                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
2905                                      print *,                                      print *,
2906       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2907       $                              //'(clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
# Line 2929  c                                 chi2=- Line 2926  c                                 chi2=-
2926  *     --------------------------  *     --------------------------
2927                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2928                                                                    
2929                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
2930       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2931       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2932       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
# Line 2945  c                                    mas Line 2942  c                                    mas
2942                                                                
2943                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2944                                                                
2945                                do ip=1,nplanes                                do ip=1,nplanes !top to bottom
2946    
2947                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2948                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
# Line 2962  c                                    mas Line 2959  c                                    mas
2959                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
2960                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
2961                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
2962    *                                NB! hit_plane is defined from bottom to top
2963                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2964                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2965       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
# Line 2977  c                                    mas Line 2975  c                                    mas
2975                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0                                      SENSOR_STORE(nplanes-ip+1,ntracks)=0
2976                                      LADDER_STORE(nplanes-ip+1,ntracks)=0                                      LADDER_STORE(nplanes-ip+1,ntracks)=0
2977                                   endif                                   endif
2978                                   BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now                                   BX_STORE(ip,ntracks)=0!I dont need it now
2979                                   BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now                                   BY_STORE(ip,ntracks)=0!I dont need it now
2980                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(ip,ntracks)=0
2981                                   do i=1,5                                   do i=1,5
2982                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2983                                   enddo                                   enddo
# Line 3008  c                                    mas Line 3006  c                                    mas
3006           return           return
3007        endif        endif
3008                
3009  c$$$      if(DEBUG)then  c$$$      if(DEBUG.EQ.1)then
3010  c$$$         print*,'****** TRACK CANDIDATES ***********'  c$$$         print*,'****** TRACK CANDIDATES ***********'
3011  c$$$         print*,'#         R. chi2        RIG'  c$$$         print*,'#         R. chi2        RIG'
3012  c$$$         do i=1,ntracks  c$$$         do i=1,ntracks
# Line 3017  c$$$     $           ,1./abs(AL_STORE(5, Line 3015  c$$$     $           ,1./abs(AL_STORE(5,
3015  c$$$         enddo  c$$$         enddo
3016  c$$$         print*,'***********************************'  c$$$         print*,'***********************************'
3017  c$$$      endif  c$$$      endif
3018        if(DEBUG)then        if(DEBUG.EQ.1)then
3019          print*,'****** TRACK CANDIDATES *****************'          print*,'****** TRACK CANDIDATES *****************'
3020          print*,'#         R. chi2        RIG         ndof'          print*,'#         R. chi2        RIG         ndof'
3021          do i=1,ntracks          do i=1,ntracks
# Line 3069  c$$$      endif Line 3067  c$$$      endif
3067        real xyzp(3),bxyz(3)        real xyzp(3),bxyz(3)
3068        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3069    
3070          if(DEBUG.EQ.1)print*,'refine_track:'
3071  *     =================================================  *     =================================================
3072  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3073  *                          and  *                          and
# Line 3077  c$$$      endif Line 3076  c$$$      endif
3076        call track_init        call track_init
3077        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3078    
3079             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3080    
3081           xP=XV_STORE(nplanes-ip+1,ibest)           xP=XV_STORE(nplanes-ip+1,ibest)
3082           yP=YV_STORE(nplanes-ip+1,ibest)           yP=YV_STORE(nplanes-ip+1,ibest)
3083           zP=ZV_STORE(nplanes-ip+1,ibest)           zP=ZV_STORE(nplanes-ip+1,ibest)
# Line 3150  c     $           AYV_STORE(nplanes-ip+1 Line 3151  c     $           AYV_STORE(nplanes-ip+1
3151              LADDER_STORE(nplanes-ip+1,IBEST)=nldt              LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3152  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3153    
3154              if(DEBUG)then              if(DEBUG.EQ.1)then
3155                 print*,                 print*,
3156       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3157       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3161  c     $           AYV_STORE(nplanes-ip+1 Line 3162  c     $           AYV_STORE(nplanes-ip+1
3162  *     ===========================================  *     ===========================================
3163  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3164  *     ===========================================  *     ===========================================
3165  c            if(DEBUG)print*,'>>>> try to include a new couple'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
3166              distmin=1000000.              distmin=1000000.
3167              xmm = 0.              xmm = 0.
3168              ymm = 0.              ymm = 0.
# Line 3178  c            if(DEBUG)print*,'>>>> try t Line 3179  c            if(DEBUG)print*,'>>>> try t
3179                 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
3180  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
3181  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
3182       $              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
3183       $              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
3184       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3185  *            *          
3186                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3193  c     $              cl_used(icy).eq.1.o Line 3194  c     $              cl_used(icy).eq.1.o
3194                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3195  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3196                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3197                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)print*,'( couple ',id
3198       $              ,' ) distance ',distance       $              ,' ) distance ',distance
3199                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3200                    xmm = xPAM                    xmm = xPAM
# Line 3225  c            if(distmin.le.clinc)then   Line 3226  c            if(distmin.le.clinc)then  
3226                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3227  *              -----------------------------------  *              -----------------------------------
3228                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3229                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3230       $              ,' (dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3231                 goto 133         !next plane                 goto 133         !next plane
3232              endif              endif
# Line 3233  c            if(distmin.le.clinc)then   Line 3234  c            if(distmin.le.clinc)then  
3234  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3235  *                     either from a couple or single  *                     either from a couple or single
3236  *     ================================================  *     ================================================
3237  c            if(DEBUG)print*,'>>>> try to include a new cluster'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
3238              distmin=1000000.              distmin=1000000.
3239              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3240              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3270  c     $              AXV_STORE(nplanes-i Line 3271  c     $              AXV_STORE(nplanes-i
3271       $              )                     $              )              
3272                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3273  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3274                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)print*,'( cl-X ',icx
3275       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3276                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3277                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3303  c     $              0.,AYV_STORE(nplane Line 3304  c     $              0.,AYV_STORE(nplane
3304       $              )       $              )
3305                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3306  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3307                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy
3308       $              ,' in cp ',id,' ) distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3309                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3310                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 3325  c                 dedxmm = sgnl(icy)  !( Line 3326  c                 dedxmm = sgnl(icy)  !(
3326  *----- single clusters -----------------------------------------------    *----- single clusters -----------------------------------------------  
3327  c            print*,'## ncls(',ip,') ',ncls(ip)  c            print*,'## ncls(',ip,') ',ncls(ip)
3328              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3329    c               print*,'-',ic,'-'
3330                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3331  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
3332                 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 3350  c               if(cl_used(icl).eq.1.or.
3350    
3351                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3352  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3353                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)print*,'( cl-s ',icl
3354       $              ,' ) distance ',distance,'<',distmin,' ?'       $              ,' ) distance ',distance
3355                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3356                    if(DEBUG)print*,'YES'  c                  if(DEBUG.EQ.1)print*,'YES'
3357                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3358                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3359                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3397  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3399  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3399                    if(mod(VIEW(iclm),2).eq.0)then                    if(mod(VIEW(iclm),2).eq.0)then
3400                       XGOOD(nplanes-ip+1)=1.                       XGOOD(nplanes-ip+1)=1.
3401                       resx(nplanes-ip+1)=rxmm                       resx(nplanes-ip+1)=rxmm
3402                       if(DEBUG)print*,'%%%% included X-cl ',iclm                       if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
3403       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3404       $                    ,', dist.= ',distmin       $                    ,', dist.= ',distmin
3405       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clinc,' )'
3406                    else                    else
3407                       YGOOD(nplanes-ip+1)=1.                       YGOOD(nplanes-ip+1)=1.
3408                       resy(nplanes-ip+1)=rymm                       resy(nplanes-ip+1)=rymm
3409                       if(DEBUG)print*,'%%%% included Y-cl ',iclm                       if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
3410       $                    ,'( chi^2, ',RCHI2_STORE(ibest)       $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3411       $                    ,', dist.= ', distmin       $                    ,', dist.= ', distmin
3412       $                    ,', cut ',clinc,' )'       $                    ,', cut ',clinc,' )'
# Line 3413  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3415  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3415  *     ----------------------------  *     ----------------------------
3416                    xm_A(nplanes-ip+1) = xmm_A                    xm_A(nplanes-ip+1) = xmm_A
3417                    ym_A(nplanes-ip+1) = ymm_A                    ym_A(nplanes-ip+1) = ymm_A
3418                      zm_A(nplanes-ip+1) = zmm_A
3419                    xm_B(nplanes-ip+1) = xmm_B                    xm_B(nplanes-ip+1) = xmm_B
3420                    ym_B(nplanes-ip+1) = ymm_B                    ym_B(nplanes-ip+1) = ymm_B
3421                      zm_B(nplanes-ip+1) = zmm_B
3422                    zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.                    zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3423    c$$$                  zm(nplanes-ip+1) =
3424    c$$$     $                 z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
3425                    dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<                    dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3426                    dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<                    dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3427  *     ----------------------------  *     ----------------------------
# Line 3446  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3452  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3452        include 'common_momanhough.f'        include 'common_momanhough.f'
3453        include 'level2.f'              include 'level2.f'      
3454    
3455          if(DEBUG.EQ.1)print*,'clean_XYclouds:'
3456    
3457        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3458    
3459           id=CP_STORE(nplanes-ip+1,ibest)           id=CP_STORE(nplanes-ip+1,ibest)
# Line 3454  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3462  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3462              if(id.ne.0)then              if(id.ne.0)then
3463                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3464                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3465                 cl_used(iclx)=ntrk  !tag used clusters  c$$$               cl_used(iclx)=ntrk  !tag used clusters
3466                 cl_used(icly)=ntrk  !tag used clusters  c$$$               cl_used(icly)=ntrk  !tag used clusters
3467              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3468                 cl_used(icl)=ntrk   !tag used clusters  c$$$               cl_used(icl)=ntrk   !tag used clusters
3469              endif              endif
3470                            
3471  *     -----------------------------  *     -----------------------------
# Line 3476  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3484  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3484       $              cly(ip,icp).eq.icl       $              cly(ip,icp).eq.icl
3485       $              )then       $              )then
3486                    id=id_cp(ip,icp,1)                    id=id_cp(ip,icp,1)
3487                    if(DEBUG)then                    if(DEBUG.EQ.1)then
3488                       print*,ip,' <<< cp ',id                       print*,ip,' <<< cp ',id
3489       $                    ,' ( cl-x '       $                    ,' ( cl-x '
3490       $                    ,clx(ip,icp)       $                    ,clx(ip,icp)
# Line 3557  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3565  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3565              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3566              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3567              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3568                multmaxx(ip,it) = 0
3569                seedx(ip,it)    = 0
3570                xpu(ip,it)      = 0
3571                multmaxy(ip,it) = 0
3572                seedy(ip,it)    = 0
3573                ypu(ip,it)      = 0
3574           enddo           enddo
3575           do ipa=1,5           do ipa=1,5
3576              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3576  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3590  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3590          ys(1,ip)=0          ys(1,ip)=0
3591          ys(2,ip)=0          ys(2,ip)=0
3592          sgnlys(ip)=0          sgnlys(ip)=0
3593            sxbad(ip)=0
3594            sybad(ip)=0
3595            multmaxsx(ip)=0
3596            multmaxsy(ip)=0
3597        enddo        enddo
3598        end        end
3599    
# Line 3736  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3754  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3754           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3755           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3756           ayv_nt(ip,ntr)   = sngl(ayv(ip))             ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3757  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
3758           factor = sqrt(           factor = sqrt(
3759       $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3760       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3761       $        1. )       $        1. )
3762  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
3763           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3764           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3765        
# Line 3749  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Line 3767  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3767           ay   = ayv_nt(ip,ntr)           ay   = ayv_nt(ip,ntr)
3768           bfx  = BX_STORE(ip,IDCAND)           bfx  = BX_STORE(ip,IDCAND)
3769           bfy  = BY_STORE(ip,IDCAND)           bfy  = BY_STORE(ip,IDCAND)
3770           if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)  c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3771           if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)  c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3772           tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001  c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3773           angx     = 180.*atan(tgtemp)/acos(-1.)  c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3774           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001          c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3775           angy    = 180.*atan(tgtemp)/acos(-1.)  c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3776    
3777             angx = effectiveangle(ax,2*ip,bfy)
3778             angy = effectiveangle(ay,2*ip-1,bfx)
3779            
3780                    
3781  c         print*,'* ',ip,bfx,bfy,angx,angy  c         print*,'* ',ip,bfx,bfy,angx,angy
3782    
# Line 3771  c         print*,'* ',ip,bfx,bfy,angx,an Line 3793  c         print*,'* ',ip,bfx,bfy,angx,an
3793  c           >>> is a couple  c           >>> is a couple
3794              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3795              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3796                
3797  c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)              cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3798  c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)                          cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3799  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)))  
3800              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3801              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3802    
# Line 3785  c$$$            ybad(ip,ntr)= nbadstrips Line 3806  c$$$            ybad(ip,ntr)= nbadstrips
3806              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3807       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3808    
3809                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3810         $                         +10000*mult(cltrx(ip,ntr))
3811                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3812         $           /clsigma(indmax(cltrx(ip,ntr)))
3813                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3814                xpu(ip,ntr)      = corr
3815    
3816                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3817         $                         +10000*mult(cltry(ip,ntr))
3818                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3819         $           /clsigma(indmax(cltry(ip,ntr)))
3820                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3821                ypu(ip,ntr)      = corr
3822    
3823           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3824    
3825    
3826                cl_used(icl) = 1    !tag used clusters          
3827    
3828              if(mod(VIEW(icl),2).eq.0)then              if(mod(VIEW(icl),2).eq.0)then
3829                 cltrx(ip,ntr)=icl                 cltrx(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angx)  
 c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)  
3830                 xbad(ip,ntr) = nbadstrips(4,icl)                 xbad(ip,ntr) = nbadstrips(4,icl)
3831    
3832                 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)
3833    
3834                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3835         $                         +10000*mult(cltrx(ip,ntr))
3836                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3837         $           /clsigma(indmax(cltrx(ip,ntr)))
3838                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3839                   xpu(ip,ntr)      = corr
3840    
3841              elseif(mod(VIEW(icl),2).eq.1)then              elseif(mod(VIEW(icl),2).eq.1)then
3842                 cltry(ip,ntr)=icl                 cltry(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angy)  
 c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)  
3843                 ybad(ip,ntr) = nbadstrips(4,icl)                 ybad(ip,ntr) = nbadstrips(4,icl)
3844    
3845                 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)
3846    
3847                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3848         $                         +10000*mult(cltry(ip,ntr))
3849                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3850         $           /clsigma(indmax(cltry(ip,ntr)))
3851                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3852                   ypu(ip,ntr)      = corr
3853                  
3854              endif              endif
3855    
3856           endif                     endif          
3857    
3858        enddo        enddo
3859    
3860          if(DEBUG.eq.1)then
3861             print*,'> STORING TRACK ',ntr
3862             print*,'clusters: '
3863             do ip=1,6
3864                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3865             enddo
3866          endif
3867    
3868  c$$$      print*,(xgood(i),i=1,6)  c$$$      print*,(xgood(i),i=1,6)
3869  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 3898  c         if( mask_view(iv).ne.0 )good2(
3898           good2(iv) = good2(iv) + mask_view(iv)*2**8           good2(iv) = good2(iv) + mask_view(iv)*2**8
3899        enddo        enddo
3900    
3901          if(DEBUG.eq.1)then
3902             print*,'> STORING SINGLETS '
3903          endif
3904    
3905        do icl=1,nclstr1        do icl=1,nclstr1
3906    
3907             ip=nplanes-npl(VIEW(icl))+1            
3908            
3909           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
3910              ip=nplanes-npl(VIEW(icl))+1              
3911    cc            print*,' ic ',icl,' --- stored '
3912    
3913              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3914    
3915                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3916                 planex(nclsx) = ip                 planex(nclsx) = ip
3917                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3918                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3919                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3920                   sxbad(nclsx)  = nbadstrips(1,icl)
3921                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
3922                  
3923    cc               print*,icl,' >>>> ',sxbad(nclsx)
3924    
3925                 do is=1,2                 do is=1,2
3926  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3927  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 3939  c$$$               print*,'xs(2,nclsx)  
3939                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3940                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3941                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3942                   sybad(nclsy)  = nbadstrips(1,icl)
3943                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
3944    
3945    cc               print*,icl,' >>>> ',sybad(nclsy)
3946    
3947                 do is=1,2                 do is=1,2
3948  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3949  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 3960  c$$$               print*,'ys(2,nclsy)  
3960    
3961  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3962           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
3963    *     --------------------------------------------------
3964    *     per non perdere la testa...
3965    *     whichtrack(icl) e` una variabile del common level1
3966    *     che serve solo per sapere quali cluster sono stati
3967    *     associati ad una traccia, e permettere di salvare
3968    *     solo questi nell'albero di uscita
3969    *     --------------------------------------------------
3970                    
3971        enddo        enddo
3972        end        end
3973    
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.23