/[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.18 by pam-fi, Fri Feb 16 14:56:02 2007 UTC revision 1.26 by pam-fi, Thu May 24 16:45:48 2007 UTC
# Line 20  Line 20 
20        include 'calib.f'        include 'calib.f'
21        include 'level2.f'        include 'level2.f'
22    
23  c      include 'momanhough_init.f'  
24    
25    c      print*,'======================================================'
26    c$$$      do ic=1,NCLSTR1
27    c$$$         if(.false.
28    c$$$     $        .or.nsatstrips(ic).gt.0
29    c$$$c     $        .or.nbadstrips(0,ic).gt.0
30    c$$$c     $        .or.nbadstrips(4,ic).gt.0
31    c$$$c     $        .or.nbadstrips(3,ic).gt.0
32    c$$$     $        .or..false.)then
33    c$$$            print*,'--- cl-',ic,' ------------------------'
34    c$$$            istart = INDSTART(IC)
35    c$$$            istop  = TOTCLLENGTH
36    c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
37    c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
38    c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
39    c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
40    c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
41    c$$$            print*,'view ',VIEW(ic)
42    c$$$            print*,'maxs ',MAXS(ic)
43    c$$$            print*,'COG4 ',cog(4,ic)
44    c$$$            ff = fbad_cog(4,ic)
45    c$$$            print*,'fbad ',ff
46    c$$$            print*,(CLBAD(i),i=istart,istop)
47    c$$$            bb=nbadstrips(0,ic)
48    c$$$            print*,'#BAD (tot)',bb
49    c$$$            bb=nbadstrips(4,ic)
50    c$$$            print*,'#BAD (4)',bb
51    c$$$            bb=nbadstrips(3,ic)
52    c$$$            print*,'#BAD (3)',bb
53    c$$$            ss=nsatstrips(ic)
54    c$$$            print*,'#saturated ',ss
55    c$$$         endif
56    c$$$      enddo
57                
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
59  *     STEP 1  *     STEP 1
# Line 41  c      include 'momanhough_init.f' Line 74  c      include 'momanhough_init.f'
74  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
75  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
76    
77  c      iflag=0  
78        call cl_to_couples(iflag)        call cl_to_couples(iflag)
79        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
80           goto 880               !go to next event           goto 880               !go to next event
# Line 75  c      iflag=0 Line 108  c      iflag=0
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    
111  c      iflag=0  
112        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
113        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
114           goto 880               !go to next event                       goto 880               !go to next event            
# Line 239  c$$$         if(ibest.eq.0)goto 880 !>> Line 272  c$$$         if(ibest.eq.0)goto 880 !>>
272  *     2nd) increasing chi**2  *     2nd) increasing chi**2
273  *     -------------------------------------------------------  *     -------------------------------------------------------
274           rchi2best=1000000000.           rchi2best=1000000000.
275           ndofbest=0             !(1)           ndofbest=0            
276           do i=1,ntracks           do i=1,ntracks
277             ndof=0               !(1)             ndof=0              
278             do ii=1,nplanes      !(1)             do ii=1,nplanes    
279               ndof=ndof          !(1)               ndof=ndof        
280       $            +int(xgood_store(ii,i)) !(1)       $            +int(xgood_store(ii,i))
281       $            +int(ygood_store(ii,i)) !(1)       $            +int(ygood_store(ii,i))
282             enddo                !(1)             enddo              
283             if(ndof.gt.ndofbest)then !(1)             if(ndof.gt.ndofbest)then
284               ibest=i               ibest=i
285               rchi2best=RCHI2_STORE(i)               rchi2best=RCHI2_STORE(i)
286               ndofbest=ndof      !(1)               ndofbest=ndof    
287             elseif(ndof.eq.ndofbest)then !(1)             elseif(ndof.eq.ndofbest)then
288               if(RCHI2_STORE(i).lt.rchi2best.and.               if(RCHI2_STORE(i).lt.rchi2best.and.
289       $            RCHI2_STORE(i).gt.0)then       $            RCHI2_STORE(i).gt.0)then
290                 ibest=i                 ibest=i
291                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
292                 ndofbest=ndof    !(1)                 ndofbest=ndof  
293               endif              !(1)               endif            
294             endif             endif
295           enddo           enddo
296    
# Line 298  c$$$         enddo Line 331  c$$$         enddo
331              iimage=0              iimage=0
332           endif           endif
333           if(icand.eq.0)then           if(icand.eq.0)then
334              print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand              if(VERBOSE)then
335       $           ,ibest,iimage                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
336         $              ,ibest,iimage
337                endif
338              return              return
339           endif           endif
340    
# Line 523  c     $        rchi2best.lt.15..and. Line 558  c     $        rchi2best.lt.15..and.
558        include 'commontracker.f'        include 'commontracker.f'
559        include 'level1.f'        include 'level1.f'
560        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
561        include 'common_align.f'        include 'common_align.f'
562        include 'common_mech.f'        include 'common_mech.f'
563        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
 c      include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
564    
565        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
566        integer sensor        integer sensor
# Line 545  c      common/dbg/DEBUG Line 575  c      common/dbg/DEBUG
575        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
576        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
577        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
 c      double precision xi,yi,zi  
 c      double precision xi_A,yi_A,zi_A  
 c      double precision xi_B,yi_B,zi_B  
578                
579    
580        parameter (ndivx=30)        parameter (ndivx=30)
581    
582    
583    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
584                
585        resxPAM = 0        resxPAM = 0
586        resyPAM = 0        resyPAM = 0
# Line 565  c      double precision xi_B,yi_B,zi_B Line 595  c      double precision xi_B,yi_B,zi_B
595        yPAM_B = 0.        yPAM_B = 0.
596        zPAM_B = 0.        zPAM_B = 0.
597  c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy  c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
598    
599          if(sensor.lt.1.or.sensor.gt.2)then
600             print*,'xyz_PAM   ***ERROR*** wrong input '
601             print*,'sensor ',sensor
602             icx=0
603             icy=0
604          endif
605    
606  *     -----------------  *     -----------------
607  *     CLUSTER X  *     CLUSTER X
608  *     -----------------  *     -----------------      
   
609        if(icx.ne.0)then        if(icx.ne.0)then
610    
611           viewx = VIEW(icx)           viewx   = VIEW(icx)
612           nldx = nld(MAXS(icx),VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
613           nplx = npl(VIEW(icx))           nplx    = npl(VIEW(icx))
614           resxPAM = RESXAV           resxPAM = RESXAV
615           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
616    
617             if(
618         $        viewx.lt.1.or.        
619         $        viewx.gt.12.or.
620         $        nldx.lt.1.or.
621         $        nldx.gt.3.or.
622         $        stripx.lt.1.or.
623         $        stripx.gt.3072.or.
624         $        .false.)then
625                print*,'xyz_PAM   ***ERROR*** wrong input '
626                print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
627                icx = 0
628                goto 10
629             endif
630    
631  *        --------------------------  *        --------------------------
632  *        magnetic-field corrections  *        magnetic-field corrections
633  *        --------------------------  *        --------------------------
 c$$$         print*,nplx,ax,bfy/10.  
634           angtemp  = ax           angtemp  = ax
635           bfytemp  = bfy           bfytemp  = bfy
636           if(nplx.eq.6) angtemp = -1. * ax  *        /////////////////////////////////
637           if(nplx.eq.6) bfytemp = -1. * bfy  *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!
638           tgtemp = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001  *        *grvzkkjsdgjhhhgngbn###>:(
639           angx = 180.*atan(tgtemp)/acos(-1.)  *        /////////////////////////////////
640           stripx = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX  c         if(nplx.eq.6) angtemp = -1. * ax
641    c         if(nplx.eq.6) bfytemp = -1. * bfy
642             if(viewx.eq.12) angtemp = -1. * ax
643             if(viewx.eq.12) bfytemp = -1. * bfy
644             tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
645             angx     = 180.*atan(tgtemp)/acos(-1.)
646             stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
647    c$$$         print*,nplx,ax,bfy/10.
648  c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX  c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
649  c$$$         print*,'========================'  c$$$         print*,'========================'
650    c$$$         if(bfy.ne.0.)print*,viewx,'-x- '
651    c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ
652  *        --------------------------  *        --------------------------
653    
654    c$$$         print*,'--- x-cl ---'
655    c$$$         istart = INDSTART(ICX)
656    c$$$         istop  = TOTCLLENGTH
657    c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
658    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
659    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
660    c$$$         print*,INDMAX(icx)
661    c$$$         print*,cog(4,icx)
662    c$$$         print*,fbad_cog(4,icx)
663            
664    
665           if(PFAx.eq.'COG1')then           if(PFAx.eq.'COG1')then
666              stripx = stripx      
667              resxPAM = resxPAM                stripx  = stripx
668                resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM
669    
670           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
671              stripx = stripx + cog(2,icx)              
672                stripx  = stripx + cog(2,icx)            
673                resxPAM = risx_cog(abs(angx))!TEMPORANEO              
674              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
675    
676           elseif(PFAx.eq.'COG3')then           elseif(PFAx.eq.'COG3')then
677              stripx = stripx + cog(3,icx)              
678                stripx  = stripx + cog(3,icx)            
679                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
680              resxPAM = resxPAM*fbad_cog(3,icx)              resxPAM = resxPAM*fbad_cog(3,icx)
681    
682           elseif(PFAx.eq.'COG4')then           elseif(PFAx.eq.'COG4')then
683  c            print*,'COG4'  
684              stripx = stripx + cog(4,icx)                          stripx  = stripx + cog(4,icx)            
685                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
686              resxPAM = resxPAM*fbad_cog(4,icx)              resxPAM = resxPAM*fbad_cog(4,icx)
687    
688           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
689              stripx = stripx + pfaeta2(icx,angx)            
690              resxPAM = risx_eta2(abs(angx))                                  stripx  = stripx + pfaeta2(icx,angx)          
691                resxPAM = risxeta2(abs(angx))
692                resxPAM = resxPAM*fbad_cog(2,icx)
693              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
694       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
695              resxPAM = resxPAM*fbad_cog(2,icx)  
696           elseif(PFAx.eq.'ETA3')then                                   elseif(PFAx.eq.'ETA3')then                        
697              stripx = stripx + pfaeta3(icx,angx)            
698              resxPAM = risx_eta3(abs(angx))                                    stripx  = stripx + pfaeta3(icx,angx)          
699                resxPAM = risxeta3(abs(angx))                      
700                resxPAM = resxPAM*fbad_cog(3,icx)              
701              if(DEBUG.and.fbad_cog(3,icx).ne.1)                          if(DEBUG.and.fbad_cog(3,icx).ne.1)            
702       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
703              resxPAM = resxPAM*fbad_cog(3,icx)                
704           elseif(PFAx.eq.'ETA4')then                                   elseif(PFAx.eq.'ETA4')then                        
705              stripx = stripx + pfaeta4(icx,angx)              
706              resxPAM = risx_eta4(abs(angx))                                    stripx  = stripx + pfaeta4(icx,angx)            
707                resxPAM = risxeta4(abs(angx))                      
708                resxPAM = resxPAM*fbad_cog(4,icx)              
709              if(DEBUG.and.fbad_cog(4,icx).ne.1)                            if(DEBUG.and.fbad_cog(4,icx).ne.1)              
710       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
711              resxPAM = resxPAM*fbad_cog(4,icx)                
712           elseif(PFAx.eq.'ETA')then             elseif(PFAx.eq.'ETA')then  
713  c            print*,'ETA'  
714              stripx = stripx + pfaeta(icx,angx)                          stripx  = stripx + pfaeta(icx,angx)            
715              resxPAM = ris_eta(icx,angx)                      c            resxPAM = riseta(icx,angx)                    
716                resxPAM = riseta(viewx,angx)                    
717                resxPAM = resxPAM*fbad_eta(icx,angx)            
718              if(DEBUG.and.fbad_cog(2,icx).ne.1)                            if(DEBUG.and.fbad_cog(2,icx).ne.1)              
719       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
720              resxPAM = resxPAM*fbad_eta(icx,angx)              
721           elseif(PFAx.eq.'COG')then                     elseif(PFAx.eq.'COG')then          
722              stripx = stripx + cog(0,icx)              
723                stripx  = stripx + cog(0,icx)            
724              resxPAM = risx_cog(abs(angx))                                  resxPAM = risx_cog(abs(angx))                    
725              resxPAM = resxPAM*fbad_cog(0,icx)              resxPAM = resxPAM*fbad_cog(0,icx)
726    
727           else           else
728              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
729           endif           endif
730    
 c         print*,'%%%%%%%%%%%%'  
731    
732        endif  *     ======================================
733          *     temporary patch for saturated clusters
734    *     ======================================
735             if( nsatstrips(icx).gt.0 )then
736                stripx  = stripx + cog(4,icx)            
737                resxPAM = pitchX*1e-4/sqrt(12.)
738                cc=cog(4,icx)
739    c$$$            print*,icx,' *** ',cc
740    c$$$            print*,icx,' *** ',resxPAM
741             endif
742    
743     10   endif
744    
745        
746  *     -----------------  *     -----------------
747  *     CLUSTER Y  *     CLUSTER Y
748  *     -----------------  *     -----------------
# Line 653  c         print*,'%%%%%%%%%%%%' Line 755  c         print*,'%%%%%%%%%%%%'
755           resyPAM = RESYAV           resyPAM = RESYAV
756           stripy = float(MAXS(icy))           stripy = float(MAXS(icy))
757    
758             if(
759         $        viewy.lt.1.or.        
760         $        viewy.gt.12.or.
761         $        nldy.lt.1.or.
762         $        nldy.gt.3.or.
763         $        stripy.lt.1.or.
764         $        stripy.gt.3072.or.
765         $        .false.)then
766                print*,'xyz_PAM   ***ERROR*** wrong input '
767                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
768                icy = 0
769                goto 20
770             endif
771    
772           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
773              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG) then
774       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
775         $              ,icx,icy
776                endif
777              goto 100              goto 100
778           endif           endif
779  *        --------------------------  *        --------------------------
# Line 664  c         print*,'%%%%%%%%%%%%' Line 782  c         print*,'%%%%%%%%%%%%'
782           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001                   tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
783           angy    = 180.*atan(tgtemp)/acos(-1.)           angy    = 180.*atan(tgtemp)/acos(-1.)
784           stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY           stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
785    c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
786    c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
787  *        --------------------------  *        --------------------------
788                    
789           if(PFAy.eq.'COG1')then !(1)  c$$$         print*,'--- y-cl ---'
790              stripy = stripy     !(1)  c$$$         istart = INDSTART(ICY)
791              resyPAM = resyPAM   !(1)  c$$$         istop  = TOTCLLENGTH
792    c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
793    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
794    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
795    c$$$         print*,INDMAX(icy)
796    c$$$         print*,cog(4,icy)
797    c$$$         print*,fbad_cog(4,icy)
798    
799             if(PFAy.eq.'COG1')then
800    
801                stripy  = stripy    
802                resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  
803    
804           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
805              stripy = stripy + cog(2,icy)  
806                stripy  = stripy + cog(2,icy)
807                resyPAM = risy_cog(abs(angy))!TEMPORANEO
808              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
809    
810           elseif(PFAy.eq.'COG3')then           elseif(PFAy.eq.'COG3')then
811              stripy = stripy + cog(3,icy)  
812                stripy  = stripy + cog(3,icy)
813                resyPAM = risy_cog(abs(angy))!TEMPORANEO
814              resyPAM = resyPAM*fbad_cog(3,icy)              resyPAM = resyPAM*fbad_cog(3,icy)
815    
816           elseif(PFAy.eq.'COG4')then           elseif(PFAy.eq.'COG4')then
817              stripy = stripy + cog(4,icy)  
818                stripy  = stripy + cog(4,icy)
819                resyPAM = risy_cog(abs(angy))!TEMPORANEO
820              resyPAM = resyPAM*fbad_cog(4,icy)              resyPAM = resyPAM*fbad_cog(4,icy)
821    
822           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
823  c            cog2 = cog(2,icy)  
824  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)              stripy  = stripy + pfaeta2(icy,angy)          
825  c            stripy = stripy + etacorr              resyPAM = risyeta2(abs(angy))              
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(abs(angy))                       !   (4)  
826              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
827              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
828       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
829           elseif(PFAy.eq.'ETA3')then                         !(3)  
830              stripy = stripy + pfaeta3(icy,angy)            !(3)           elseif(PFAy.eq.'ETA3')then                      
831              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
832              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              stripy  = stripy + pfaeta3(icy,angy)
833       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)              resyPAM = resyPAM*fbad_cog(3,icy)  
834           elseif(PFAy.eq.'ETA4')then                         !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)
835              stripy = stripy + pfaeta4(icy,angy)            !(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
836              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
837              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA4')then  
838       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
839           elseif(PFAy.eq.'ETA')then                          !(3)              stripy  = stripy + pfaeta4(icy,angy)
840              stripy = stripy + pfaeta(icy,angy)             !(3)              resyPAM = resyPAM*fbad_cog(4,icy)
841              resyPAM = ris_eta(icy,angy)                     !   (4)              if(DEBUG.and.fbad_cog(4,icy).ne.1)
842  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
843              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
844              if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA')then
845       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
846                stripy  = stripy + pfaeta(icy,angy)
847    c            resyPAM = riseta(icy,angy)  
848                resyPAM = riseta(viewy,angy)  
849                resyPAM = resyPAM*fbad_eta(icy,angy)
850                if(DEBUG.and.fbad_cog(2,icy).ne.1)
851         $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
852    
853           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
854              stripy = stripy + cog(0,icy)              
855              resyPAM = risy_cog(abs(angy))                        !   (4)              stripy  = stripy + cog(0,icy)            
856  c            resyPAM = ris_eta(icy,angy)                    !   (4)              resyPAM = risy_cog(abs(angy))
857              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
858    
859           else           else
860              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
861             endif
862    
863    
864    *     ======================================
865    *     temporary patch for saturated clusters
866    *     ======================================
867             if( nsatstrips(icy).gt.0 )then
868                stripy  = stripy + cog(4,icy)            
869                resyPAM = pitchY*1e-4/sqrt(12.)
870                cc=cog(4,icy)
871    c$$$            print*,icy,' *** ',cc
872    c$$$            print*,icy,' *** ',resyPAM
873           endif           endif
874    
       endif  
875    
876  c      print*,'## stripx,stripy ',stripx,stripy   20   endif
877    
878    c$$$      print*,'## stripx,stripy ',stripx,stripy
879    
880  c===========================================================  c===========================================================
881  C     COUPLE  C     COUPLE
# Line 727  c     (xi,yi,zi) = mechanical coordinate Line 887  c     (xi,yi,zi) = mechanical coordinate
887  c------------------------------------------------------------------------  c------------------------------------------------------------------------
888           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
889       $        .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...
890              print*,'xyz_PAM (couple):',              if(DEBUG) then
891       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
892         $              ' WARNING: false X strip: strip ',stripx
893                endif
894           endif           endif
895           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
896           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
# Line 820  c            print*,'X-singlet ',icx,npl Line 982  c            print*,'X-singlet ',icx,npl
982  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...
983              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
984       $           .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...
985                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG) then
986       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
987         $                 ' WARNING: false X strip: strip ',stripx
988                   endif
989              endif              endif
990              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
991    
# Line 843  c            print*,'X-cl ',icx,stripx,' Line 1007  c            print*,'X-cl ',icx,stripx,'
1007  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
1008    
1009           else           else
1010                if(DEBUG) then
1011              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
1012              print *,'icx = ',icx                 print *,'icx = ',icx
1013              print *,'icy = ',icy                 print *,'icy = ',icy
1014                endif
1015              goto 100              goto 100
1016                            
1017           endif           endif
# Line 911  c--------------------------------------- Line 1076  c---------------------------------------
1076  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1077    
1078        else        else
1079                       if(DEBUG) then
1080           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1081           print *,'icx = ',icx              print *,'icx = ',icx
1082           print *,'icy = ',icy              print *,'icy = ',icy
1083                         endif
1084        endif        endif
1085                    
1086    
# Line 926  c      print*,'## xPAM_B,yPAM_B,zPAM_B ' Line 1091  c      print*,'## xPAM_B,yPAM_B,zPAM_B '
1091   100  continue   100  continue
1092        end        end
1093    
1094    ************************************************************************
1095    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1096    *     (done to be called from c/c++)
1097    ************************************************************************
1098    
1099          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1100    
1101          include 'commontracker.f'
1102          include 'level1.f'
1103          include 'common_mini_2.f'
1104          include 'common_xyzPAM.f'
1105          include 'common_mech.f'
1106          include 'calib.f'
1107          
1108    *     flag to chose PFA
1109    c$$$      character*10 PFA
1110    c$$$      common/FINALPFA/PFA
1111    
1112          integer icx,icy           !X-Y cluster ID
1113          integer sensor
1114          character*4 PFAx,PFAy     !PFA to be used
1115          real ax,ay                !X-Y geometric angle
1116          real bfx,bfy              !X-Y b-field components
1117    
1118          ipx=0
1119          ipy=0      
1120          
1121    c$$$      PFAx = 'COG4'!PFA
1122    c$$$      PFAy = 'COG4'!PFA
1123    
1124          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1125                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1126         $           ,' does not exists (nclstr1=',nclstr1,')'
1127                icx = -1*icx
1128                icy = -1*icy
1129                return
1130            
1131          endif
1132          
1133          call idtoc(pfaid,PFAx)
1134          call idtoc(pfaid,PFAy)
1135    
1136    c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1137    
1138    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1139          
1140          if(icx.ne.0.and.icy.ne.0)then
1141    
1142             ipx=npl(VIEW(icx))
1143             ipy=npl(VIEW(icy))
1144    c$$$         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1145    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1146    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1147            
1148             if( (nplanes-ipx+1).ne.ip )then
1149                print*,'xyzpam: ***WARNING*** cluster ',icx
1150         $           ,' does not belong to plane: ',ip
1151                icx = -1*icx
1152                return
1153             endif
1154             if( (nplanes-ipy+1).ne.ip )then
1155                print*,'xyzpam: ***WARNING*** cluster ',icy
1156         $           ,' does not belong to plane: ',ip
1157                icy = -1*icy
1158                return
1159             endif
1160    
1161             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1162    
1163             xgood(ip) = 1.
1164             ygood(ip) = 1.
1165             resx(ip)  = resxPAM
1166             resy(ip)  = resyPAM
1167    
1168             xm(ip) = xPAM
1169             ym(ip) = yPAM
1170             zm(ip) = zPAM
1171             xm_A(ip) = 0.
1172             ym_A(ip) = 0.
1173             xm_B(ip) = 0.
1174             ym_B(ip) = 0.
1175    
1176    c         zv(ip) = zPAM
1177    
1178          elseif(icx.eq.0.and.icy.ne.0)then
1179    
1180             ipy=npl(VIEW(icy))
1181    c$$$         if((nplanes-ipy+1).ne.ip)
1182    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1183    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1184             if( (nplanes-ipy+1).ne.ip )then
1185                print*,'xyzpam: ***WARNING*** cluster ',icy
1186         $           ,' does not belong to plane: ',ip
1187                icy = -1*icy
1188                return
1189             endif
1190    
1191             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1192            
1193             xgood(ip) = 0.
1194             ygood(ip) = 1.
1195             resx(ip)  = 1000.
1196             resy(ip)  = resyPAM
1197    
1198             xm(ip) = -100.
1199             ym(ip) = -100.
1200             zm(ip) = (zPAM_A+zPAM_B)/2.
1201             xm_A(ip) = xPAM_A
1202             ym_A(ip) = yPAM_A
1203             xm_B(ip) = xPAM_B
1204             ym_B(ip) = yPAM_B
1205    
1206    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1207            
1208          elseif(icx.ne.0.and.icy.eq.0)then
1209    
1210             ipx=npl(VIEW(icx))
1211    c$$$         if((nplanes-ipx+1).ne.ip)
1212    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1213    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1214    
1215             if( (nplanes-ipx+1).ne.ip )then
1216                print*,'xyzpam: ***WARNING*** cluster ',icx
1217         $           ,' does not belong to plane: ',ip
1218                icx = -1*icx
1219                return
1220             endif
1221    
1222             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1223          
1224             xgood(ip) = 1.
1225             ygood(ip) = 0.
1226             resx(ip)  = resxPAM
1227             resy(ip)  = 1000.
1228    
1229             xm(ip) = -100.
1230             ym(ip) = -100.
1231             zm(ip) = (zPAM_A+zPAM_B)/2.
1232             xm_A(ip) = xPAM_A
1233             ym_A(ip) = yPAM_A
1234             xm_B(ip) = xPAM_B
1235             ym_B(ip) = yPAM_B
1236            
1237    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1238    
1239          else
1240    
1241             il = 2
1242             if(lad.ne.0)il=lad
1243             is = 1
1244             if(sensor.ne.0)is=sensor
1245    c         print*,nplanes-ip+1,il,is
1246    
1247             xgood(ip) = 0.
1248             ygood(ip) = 0.
1249             resx(ip)  = 1000.
1250             resy(ip)  = 1000.
1251    
1252             xm(ip) = -100.
1253             ym(ip) = -100.          
1254             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1255             xm_A(ip) = 0.
1256             ym_A(ip) = 0.
1257             xm_B(ip) = 0.
1258             ym_B(ip) = 0.
1259    
1260    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1261    
1262          endif
1263    
1264          if(DEBUG)then
1265    c         print*,'----------------------------- track coord'
1266    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1267             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1268         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1269         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1270    c$$$         print*,'-----------------------------'
1271          endif
1272          end
1273    
1274  ********************************************************************************  ********************************************************************************
1275  ********************************************************************************  ********************************************************************************
# Line 1001  c      print*,'## xPAM_B,yPAM_B,zPAM_B ' Line 1345  c      print*,'## xPAM_B,yPAM_B,zPAM_B '
1345           endif                   endif        
1346    
1347           distance=           distance=
1348       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1349    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1350           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1351    
1352  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
# Line 1026  c$$$         print*,' resolution ',re Line 1371  c$$$         print*,' resolution ',re
1371  *     ----------------------  *     ----------------------
1372                    
1373           distance=           distance=
1374       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1375       $        +       $       +
1376       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1377    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1378    c$$$     $        +
1379    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1380           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1381    
1382  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1037  c$$$         print*,' resolution ',resxP Line 1385  c$$$         print*,' resolution ',resxP
1385                    
1386        else        else
1387                    
1388           print*  c         print*
1389       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1390           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1391           print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '  c         print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
1392       $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
1393        endif          endif  
1394    
1395        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1109  c--------------------------------------- Line 1457  c---------------------------------------
1457                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1458       $              .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...
1459  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...
1460                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1461       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1462                 endif                 endif
1463                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1464                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1265  c      include 'common_analysis.f' Line 1613  c      include 'common_analysis.f'
1613        is_cp=0        is_cp=0
1614        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1615        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1616        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1617    
1618        return        return
1619        end        end
# Line 1394  c      include 'level1.f' Line 1742  c      include 'level1.f'
1742  *     mask views with too many clusters  *     mask views with too many clusters
1743        do iv=1,nviews        do iv=1,nviews
1744           if( ncl_view(iv).gt. nclusterlimit)then           if( ncl_view(iv).gt. nclusterlimit)then
1745              mask_view(iv) = 1  c            mask_view(iv) = 1
1746                mask_view(iv) = mask_view(iv) + 2**0
1747              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1748       $           ,nclusterlimit,' on view ', iv,' --> masked!'       $           ,nclusterlimit,' on view ', iv,' --> masked!'
1749           endif           endif
# Line 1539  c                  cut = chcut * sch(npl Line 1888  c                  cut = chcut * sch(npl
1888       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1889       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1890       $                 ,'( ',ncouplemax,' ) --> masked!'       $                 ,'( ',ncouplemax,' ) --> masked!'
1891                    mask_view(nviewx(nplx)) = 2  c                  mask_view(nviewx(nplx)) = 2
1892                    mask_view(nviewy(nply)) = 2  c                  mask_view(nviewy(nply)) = 2
1893                      mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1894                      mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1895                    goto 10                    goto 10
1896                 endif                 endif
1897                                
# Line 1641  c      double precision xm3,ym3,zm3 Line 1992  c      double precision xm3,ym3,zm3
1992  *     --------------------------------------------  *     --------------------------------------------
1993        do ip=1,nplanes        do ip=1,nplanes
1994           if(ncp_plane(ip).gt.ncouplelimit)then           if(ncp_plane(ip).gt.ncouplelimit)then
1995              mask_view(nviewx(ip)) = 8  c            mask_view(nviewx(ip)) = 8
1996              mask_view(nviewy(ip)) = 8  c            mask_view(nviewy(ip)) = 8
1997                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1998                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1999           endif           endif
2000        enddo        enddo
2001    
# Line 1695  c     $                       (icx2,icy2 Line 2048  c     $                       (icx2,icy2
2048  c                           good2=.false.  c                           good2=.false.
2049  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
2050                             do iv=1,12                             do iv=1,12
2051                                mask_view(iv) = 3  c                              mask_view(iv) = 3
2052                                  mask_view(iv) = mask_view(iv)+ 2**2
2053                             enddo                             enddo
2054                             iflag=1                             iflag=1
2055                             return                             return
# Line 1774  c     $                                 Line 2128  c     $                                
2128  c                                    good2=.false.  c                                    good2=.false.
2129  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2130                                      do iv=1,nviews                                      do iv=1,nviews
2131                                         mask_view(iv) = 4  c                                       mask_view(iv) = 4
2132                                           mask_view(iv)=mask_view(iv)+ 2**3
2133                                      enddo                                      enddo
2134                                      iflag=1                                      iflag=1
2135                                      return                                      return
# Line 2008  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2363  c         if(ncpused.lt.ncpyz_min)goto 2
2363  c               good2=.false.  c               good2=.false.
2364  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2365              do iv=1,nviews              do iv=1,nviews
2366                 mask_view(iv) = 5  c               mask_view(iv) = 5
2367                   mask_view(iv) = mask_view(iv) + 2**4
2368              enddo              enddo
2369              iflag=1              iflag=1
2370              return              return
# Line 2218  c     print*,'check cp_used' Line 2574  c     print*,'check cp_used'
2574           enddo           enddo
2575  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2576           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2577           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2578                    
2579  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2580  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
# Line 2230  c         if(ncpused.lt.ncpxz_min)goto 2 Line 2586  c         if(ncpused.lt.ncpxz_min)goto 2
2586  c     good2=.false.  c     good2=.false.
2587  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2588              do iv=1,nviews              do iv=1,nviews
2589                 mask_view(iv) = 6  c               mask_view(iv) = 6
2590                   mask_view(iv) =  mask_view(iv) + 2**5
2591              enddo              enddo
2592              iflag=1              iflag=1
2593              return              return
# Line 2295  c$$$     $           ,(tr_cloud(iii),iii Line 2652  c$$$     $           ,(tr_cloud(iii),iii
2652  **************************************************  **************************************************
2653    
2654        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2655    
2656        include 'commontracker.f'        include 'commontracker.f'
2657        include 'level1.f'        include 'level1.f'
# Line 2305  c*************************************** Line 2659  c***************************************
2659        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2660        include 'common_mini_2.f'        include 'common_mini_2.f'
2661        include 'common_mech.f'        include 'common_mech.f'
2662  c      include 'momanhough_init.f'  
2663    
2664    
2665  *     output flag  *     output flag
# Line 2566  c                                 chi2=- Line 2920  c                                 chi2=-
2920  c                                 good2=.false.  c                                 good2=.false.
2921  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2922                                   do iv=1,nviews                                   do iv=1,nviews
2923                                      mask_view(iv) = 7  c                                    mask_view(iv) = 7
2924                                        mask_view(iv) = mask_view(iv) + 2**6
2925                                   enddo                                   enddo
2926                                   iflag=1                                   iflag=1
2927                                   return                                   return
# Line 2574  c                                 goto 8 Line 2929  c                                 goto 8
2929                                                                
2930                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2931                                                                
 c$$$                              ndof=0                                  
2932                                do ip=1,nplanes                                do ip=1,nplanes
2933  c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2934                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2935                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2936                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 2597  c$$$     $                               Line 2949  c$$$     $                              
2949                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2950                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2951       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2952                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2953         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2954                                        LADDER_STORE(nplanes-ip+1,ntracks)
2955         $                                   = LADDER(
2956         $                                   clx(ip,icp_cp(
2957         $                                   cp_match(ip,hit_plane(ip)
2958         $                                   ))));
2959                                   else                                   else
2960                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2961                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2962                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2963                                   endif                                   endif
2964                                     BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2965                                     BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2966                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(nplanes-ip+1,ntracks)=0
2967                                   do i=1,5                                   do i=1,5
2968                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2969                                   enddo                                   enddo
2970                                enddo                                enddo
2971                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2972                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2973                                                                
2974  *     --------------------------------  *     --------------------------------
# Line 2633  c$$$  rchi2=chi2/dble(ndof) Line 2992  c$$$  rchi2=chi2/dble(ndof)
2992           return           return
2993        endif        endif
2994                
2995    c$$$      if(DEBUG)then
2996    c$$$         print*,'****** TRACK CANDIDATES ***********'
2997    c$$$         print*,'#         R. chi2        RIG'
2998    c$$$         do i=1,ntracks
2999    c$$$            print*,i,' --- ',rchi2_store(i),' --- '
3000    c$$$     $           ,1./abs(AL_STORE(5,i))
3001    c$$$         enddo
3002    c$$$         print*,'***********************************'
3003    c$$$      endif
3004        if(DEBUG)then        if(DEBUG)then
3005           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
3006           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
3007           do i=1,ntracks          do i=1,ntracks
3008              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
3009       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
3010           enddo              ndof=ndof           !(1)
3011           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
3012         $           +int(ygood_store(ii,i)) !(1)
3013              enddo                 !(1)
3014              print*,i,' --- ',rchi2_store(i),' --- '
3015         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3016            enddo
3017            print*,'*****************************************'
3018        endif        endif
3019                
3020                
# Line 2659  c$$$  rchi2=chi2/dble(ndof) Line 3033  c$$$  rchi2=chi2/dble(ndof)
3033    
3034        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3035    
 c******************************************************  
 cccccc 06/10/2005 modified by elena vannuccini ---> (1)  
 cccccc 31/01/2006 modified by elena vannuccini ---> (2)  
 cccccc 12/08/2006 modified by elena vannucicni ---> (3)  
 c******************************************************  
3036    
3037        include 'commontracker.f'        include 'commontracker.f'
3038        include 'level1.f'        include 'level1.f'
# Line 2671  c*************************************** Line 3040  c***************************************
3040        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3041        include 'common_mini_2.f'        include 'common_mini_2.f'
3042        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
3043        include 'calib.f'        include 'calib.f'
3044    
3045  *     flag to chose PFA  *     flag to chose PFA
3046        character*10 PFA        character*10 PFA
3047        common/FINALPFA/PFA        common/FINALPFA/PFA
3048    
3049          real k(6)
3050          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3051    
3052        real xp,yp,zp        real xp,yp,zp
3053        real xyzp(3),bxyz(3)        real xyzp(3),bxyz(3)
3054        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))        equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
# Line 2695  c      include 'level1.f' Line 3065  c      include 'level1.f'
3065           yP=YV_STORE(nplanes-ip+1,ibest)           yP=YV_STORE(nplanes-ip+1,ibest)
3066           zP=ZV_STORE(nplanes-ip+1,ibest)           zP=ZV_STORE(nplanes-ip+1,ibest)
3067           call gufld(xyzp,bxyz)           call gufld(xyzp,bxyz)
3068  c$$$         bxyz(1)=0           BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3069             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3070    c$$$  bxyz(1)=0
3071  c$$$         bxyz(2)=0  c$$$         bxyz(2)=0
3072  c$$$         bxyz(3)=0  c$$$         bxyz(3)=0
3073  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
# Line 2729  c     $           AYV_STORE(nplanes-ip+1 Line 3101  c     $           AYV_STORE(nplanes-ip+1
3101       $           bxyz(1),       $           bxyz(1),
3102       $           bxyz(2)       $           bxyz(2)
3103       $           )       $           )
3104  c$$$  call xyz_PAM(icx,icy,is,  
 c$$$  $              'COG2','COG2',  
 c$$$  $              0.,  
 c$$$  $              0.)  
3105              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3106              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3107              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 2741  c$$$  $              0.) Line 3110  c$$$  $              0.)
3110              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3111              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3112    
3113  c            dedxtrk(nplanes-ip+1) = (sgnl(icx)+sgnl(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3114              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
             dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  
3115                            
3116  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3117  *     -------------------------------------------------  *     -------------------------------------------------
# Line 2758  c            dedxtrk(nplanes-ip+1) = (sg Line 3126  c            dedxtrk(nplanes-ip+1) = (sg
3126                                
3127  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3128  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
 c$$$            xP=XV_STORE(nplanes-ip+1,ibest)  
 c$$$            yP=YV_STORE(nplanes-ip+1,ibest)  
 c$$$            zP=ZV_STORE(nplanes-ip+1,ibest)  
3129              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3130  *     if the track hit the plane in a dead area, go to the next plane  *     if the track hit the plane in a dead area, go to the next plane
3131              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3132    
3133                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3134                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3135  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3136    
3137              if(DEBUG)then              if(DEBUG)then
# Line 2798  c     $              cl_used(icy).eq.1.o Line 3166  c     $              cl_used(icy).eq.1.o
3166       $              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 !(3)
3167       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3168  *            *          
 c               call xyz_PAM(icx,icy,ist,  
 c     $              PFA,PFA,  
 c     $              AXV_STORE(nplanes-ip+1,ibest),  
 c     $              AYV_STORE(nplanes-ip+1,ibest))  
3169                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3170       $              PFA,PFA,       $              PFA,PFA,
3171       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
# Line 2811  c     $              AYV_STORE(nplanes-i Line 3175  c     $              AYV_STORE(nplanes-i
3175       $              )       $              )
3176                                
3177                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3178                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3179                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3180                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
3181       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3182                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3183                    xmm = xPAM                    xmm = xPAM
3184                    ymm = yPAM                    ymm = yPAM
# Line 2823  c     $              AYV_STORE(nplanes-i Line 3187  c     $              AYV_STORE(nplanes-i
3187                    rymm = resyPAM                    rymm = resyPAM
3188                    distmin = distance                    distmin = distance
3189                    idm = id                                      idm = id                  
 c                 dedxmm = (sgnl(icx)+sgnl(icy))/2. !(1)  
3190                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3191                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3192    c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3193                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3194         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3195                 endif                 endif
3196   1188          continue   1188          continue
3197              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3198              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3199                if(distmin.le.clincnewc)then     !QUIQUI              
3200  *              -----------------------------------  *              -----------------------------------
3201                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3202                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3203                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3204                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3205                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3206                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3207                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3208  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3209                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3210  *              -----------------------------------  *              -----------------------------------
3211                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3212                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG)print*,'%%%% included couple ',idm
3213       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3214                 goto 133         !next plane                 goto 133         !next plane
3215              endif              endif
3216  *     ================================================  *     ================================================
# Line 2887  c     $              AXV_STORE(nplanes-i Line 3253  c     $              AXV_STORE(nplanes-i
3253       $              bxyz(2)       $              bxyz(2)
3254       $              )                     $              )              
3255                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3256                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3257                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3258       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3259                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3260                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3261                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2920  c     $              0.,AYV_STORE(nplane Line 3286  c     $              0.,AYV_STORE(nplane
3286       $              bxyz(2)       $              bxyz(2)
3287       $              )       $              )
3288                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3289                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3290                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3291       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3292                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3293                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3294                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2944  c                 dedxmm = sgnl(icy)  !( Line 3310  c                 dedxmm = sgnl(icy)  !(
3310  c            print*,'## ncls(',ip,') ',ncls(ip)  c            print*,'## ncls(',ip,') ',ncls(ip)
3311              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3312                 icl=cls(ip,ic)                 icl=cls(ip,ic)
 c              print*,'## ic ',ic,' ist ',ist  
3313  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
3314                 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)
3315       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3316       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3317                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
 c                  call xyz_PAM(icl,0,ist,  
 c     $                 PFA,PFA,  
 c     $                 AXV_STORE(nplanes-ip+1,ibest),0.)  
3318                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
3319       $                 PFA,PFA,       $                 PFA,PFA,
3320       $                 AXV_STORE(nplanes-ip+1,ibest),0.,       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
# Line 2960  c     $                 AXV_STORE(nplane Line 3322  c     $                 AXV_STORE(nplane
3322       $                 bxyz(2)       $                 bxyz(2)
3323       $                 )       $                 )
3324                 else                         !<---- Y view                 else                         !<---- Y view
 c                  call xyz_PAM(0,icl,ist,  
 c     $                 PFA,PFA,  
 c     $                 0.,AYV_STORE(nplanes-ip+1,ibest))  
3325                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
3326       $                 PFA,PFA,       $                 PFA,PFA,
3327       $                 0.,AYV_STORE(nplanes-ip+1,ibest),       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
# Line 2972  c     $                 0.,AYV_STORE(npl Line 3331  c     $                 0.,AYV_STORE(npl
3331                 endif                 endif
3332    
3333                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3334                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3335                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3336       $              ,' ) normalized distance ',distance,'<',distmin,' ?'       $              ,' ) distance ',distance,'<',distmin,' ?'
3337                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3338                    if(DEBUG)print*,'YES'                    if(DEBUG)print*,'YES'
3339                    xmm_A = xPAM_A                    xmm_A = xPAM_A
# Line 2987  c     $                 0.,AYV_STORE(npl Line 3346  c     $                 0.,AYV_STORE(npl
3346                    rymm = resyPAM                    rymm = resyPAM
3347                    distmin = distance                      distmin = distance  
3348                    iclm = icl                    iclm = icl
 c                  dedxmm = sgnl(icl)                   !(1)  
3349                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3350                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3351                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3352                    else          !<---- Y view                    else          !<---- Y view
3353                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3354                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3355                    endif                    endif
3356                 endif                                   endif                  
3357  18882          continue  18882          continue
3358              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3359  c            print*,'## distmin ', distmin,' clinc ',clinc  c            print*,'## distmin ', distmin,' clinc ',clinc
3360              if(distmin.le.clinc)then                    
3361                  c     QUIQUI------------
3362                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<      c     anche qui: non ci vuole clinc???
3363  *              ----------------------------              if(iclm.ne.0)then
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3364                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3365                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3366                    resx(nplanes-ip+1)=rxmm       $                 20*
3367                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3368  c                  if(.true.)print*,'%%%% included X-cl ',iclm                 else if(mod(VIEW(iclm),2).ne.0)then
3369       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    clincnew=
3370       $                 ,', norm.dist.= ',distmin       $                 10*
3371       $                 ,', cut ',clinc,' )'       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3372                 else                 endif
3373                    YGOOD(nplanes-ip+1)=1.  c     QUIQUI------------
3374                    resy(nplanes-ip+1)=rymm                
3375                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 if(distmin.le.clincnew)then   !QUIQUI
3376  c                  if(.true.)print*,'%%%% included Y-cl ',iclm  c     if(distmin.le.clinc)then          !QUIQUI          
3377       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    
3378       $                 ,', norm.dist.= ', distmin                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3379       $                 ,', cut ',clinc,' )'  *     ----------------------------
3380    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3381                      if(mod(VIEW(iclm),2).eq.0)then
3382                         XGOOD(nplanes-ip+1)=1.
3383                         resx(nplanes-ip+1)=rxmm
3384                         if(DEBUG)print*,'%%%% included X-cl ',iclm
3385         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3386         $                    ,', dist.= ',distmin
3387         $                    ,', cut ',clinc,' )'
3388                      else
3389                         YGOOD(nplanes-ip+1)=1.
3390                         resy(nplanes-ip+1)=rymm
3391                         if(DEBUG)print*,'%%%% included Y-cl ',iclm
3392         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3393         $                    ,', dist.= ', distmin
3394         $                    ,', cut ',clinc,' )'
3395                      endif
3396    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3397    *     ----------------------------
3398                      xm_A(nplanes-ip+1) = xmm_A
3399                      ym_A(nplanes-ip+1) = ymm_A
3400                      xm_B(nplanes-ip+1) = xmm_B
3401                      ym_B(nplanes-ip+1) = ymm_B
3402                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3403                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3404                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3405    *     ----------------------------
3406                 endif                 endif
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
 *              ----------------------------  
                xm_A(nplanes-ip+1) = xmm_A  
                ym_A(nplanes-ip+1) = ymm_A  
                xm_B(nplanes-ip+1) = xmm_B  
                ym_B(nplanes-ip+1) = ymm_B  
                zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.  
 c              dedxtrk(nplanes-ip+1) = dedxmm !<<<    !(1)  
                dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1)  
                dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1)  
 *              ----------------------------  
3407              endif              endif
3408           endif           endif
3409   133     continue   133     continue
# Line 3050  c              dedxtrk(nplanes-ip+1) = d Line 3422  c              dedxtrk(nplanes-ip+1) = d
3422  *                                                 *  *                                                 *
3423  *                                                 *  *                                                 *
3424  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3425  *  *
3426        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3427    
3428        include 'commontracker.f'        include 'commontracker.f'
3429        include 'level1.f'        include 'level1.f'
3430        include 'common_momanhough.f'        include 'common_momanhough.f'
3431  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
   
3432    
3433        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3434    
# Line 3072  c      include 'level1.f' Line 3438  c      include 'level1.f'
3438              if(id.ne.0)then              if(id.ne.0)then
3439                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3440                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3441  c               cl_used(iclx)=1  !tag used clusters                 cl_used(iclx)=ntrk  !tag used clusters
3442  c               cl_used(icly)=1  !tag used clusters                 cl_used(icly)=ntrk  !tag used clusters
                cl_used(iclx)=ntrk  !tag used clusters !(1)  
                cl_used(icly)=ntrk  !tag used clusters !(1)  
3443              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3444  c               cl_used(icl)=1   !tag used clusters                 cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3445              endif              endif
3446                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3447  *     -----------------------------  *     -----------------------------
3448  *     remove the couple from clouds  *     remove the couple from clouds
3449  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3144  c               endif Line 3504  c               endif
3504        include 'level1.f'        include 'level1.f'
3505        include 'common_momanhough.f'        include 'common_momanhough.f'
3506        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3507    
3508    *     ---------------------------------
3509    *     variables initialized from level1
3510    *     ---------------------------------
3511        do i=1,nviews        do i=1,nviews
3512           good2(i)=good1(i)           good2(i)=good1(i)
3513             do j=1,nva1_view
3514                vkflag(i,j)=1
3515                if(cnnev(i,j).le.0)then
3516                   vkflag(i,j)=cnnev(i,j)
3517                endif
3518             enddo
3519        enddo        enddo
3520    *     ----------------
3521    *     level2 variables
3522    *     ----------------
3523        NTRK = 0        NTRK = 0
3524        do it=1,NTRKMAX        do it=1,NTRKMAX
3525           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3161  c      include 'level1.f' Line 3530  c      include 'level1.f'
3530              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3531              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3532              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3533                TAILX_nt(IP,IT) = 0
3534                TAILY_nt(IP,IT) = 0
3535                XBAD(IP,IT) = 0
3536                YBAD(IP,IT) = 0
3537              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3538              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3539                LS(IP,IT) = 0
3540              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3541              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3542              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
# Line 3295  c      include 'level1.f' Line 3669  c      include 'level1.f'
3669    
3670            
3671        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3672        include 'level1.f'        include 'level1.f'
3673        include 'common_momanhough.f'        include 'common_momanhough.f'
3674        include 'level2.f'        include 'level2.f'
3675        include 'common_mini_2.f'        include 'common_mini_2.f'
3676        real sinth,phi,pig              include 'calib.f'
3677    
3678          character*10 PFA
3679          common/FINALPFA/PFA
3680    
3681          real sinth,phi,pig
3682          integer ssensor,sladder
3683        pig=acos(-1.)        pig=acos(-1.)
3684    
3685    *     -------------------------------------
3686        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3687        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3688    *     -------------------------------------
3689        phi   = al(4)                  phi   = al(4)          
3690        sinth = al(3)                    sinth = al(3)            
3691        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3318  c      include 'level1.f' Line 3698  c      include 'level1.f'
3698       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3699        al(4) = phi                      al(4) = phi              
3700        al(3) = sinth                    al(3) = sinth            
   
3701        do i=1,5        do i=1,5
3702           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3703           do j=1,5           do j=1,5
3704              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3705           enddo           enddo
3706        enddo        enddo
3707          *     -------------------------------------      
3708        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3709           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3710           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3334  c      include 'level1.f' Line 3713  c      include 'level1.f'
3713           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3714           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3715           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3716             TAILX_nt(IP,ntr) = 0.
3717             TAILY_nt(IP,ntr) = 0.
3718           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3719           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3720           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
# Line 3341  c      include 'level1.f' Line 3722  c      include 'level1.f'
3722           ayv_nt(ip,ntr)   = sngl(ayv(ip))             ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3723  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3724           factor = sqrt(           factor = sqrt(
3725       $        sin( acos(-1.) * sngl(axv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3726       $        sin( acos(-1.) * sngl(ayv(ip)) /180. )**2 +       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3727       $        1. )       $        1. )
3728  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3729           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3730           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3731        
3732           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3733             ay   = ayv_nt(ip,ntr)
3734             bfx  = BX_STORE(ip,IDCAND)
3735             bfy  = BY_STORE(ip,IDCAND)
3736             if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3737             if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3738             tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3739             angx     = 180.*atan(tgtemp)/acos(-1.)
3740             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3741             angy    = 180.*atan(tgtemp)/acos(-1.)
3742            
3743    c         print*,'* ',ip,bfx,bfy,angx,angy
3744    
3745             id  = CP_STORE(ip,IDCAND) ! couple id
3746           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3747             ssensor = -1
3748             sladder = -1
3749             ssensor = SENSOR_STORE(ip,IDCAND)
3750             sladder = LADDER_STORE(ip,IDCAND)
3751             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3752             LS(IP,ntr)      = ssensor+10*sladder
3753    
3754           if(id.ne.0)then           if(id.ne.0)then
3755    c           >>> is a couple
3756              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3757              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3758  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)              
3759    c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3760    c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3761    c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
3762    c$$$            ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
3763                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3764                ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3765    
3766    
3767                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3768         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3769                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3770         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3771    
3772           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3773              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3774              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl              if(mod(VIEW(icl),2).eq.0)then
3775  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)                 cltrx(ip,ntr)=icl
3776    c$$$               nnnnn = npfastrips(icl,PFA,angx)
3777    c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)
3778                   xbad(ip,ntr) = nbadstrips(4,icl)
3779    
3780                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3781                elseif(mod(VIEW(icl),2).eq.1)then
3782                   cltry(ip,ntr)=icl
3783    c$$$               nnnnn = npfastrips(icl,PFA,angy)
3784    c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)
3785                   ybad(ip,ntr) = nbadstrips(4,icl)
3786                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3787                endif
3788    
3789           endif                     endif          
3790    
3791        enddo        enddo
3792    
3793    
3794    c$$$      print*,(xgood(i),i=1,6)
3795    c$$$      print*,(ygood(i),i=1,6)
3796    c$$$      print*,(ls(i,ntr),i=1,6)
3797    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3798    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3799    c$$$      print*,'-----------------------'
3800    
3801        end        end
3802    
3803        subroutine fill_level2_siglets        subroutine fill_level2_siglets
# Line 3374  c            print*,ip,' ',cltrx(ip,ntr) Line 3809  c            print*,ip,' ',cltrx(ip,ntr)
3809  *     -------------------------------------------------------  *     -------------------------------------------------------
3810    
3811        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3812        include 'calib.f'        include 'calib.f'
3813        include 'level1.f'        include 'level1.f'
3814        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 3382  c      include 'level1.f' Line 3816  c      include 'level1.f'
3816        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3817    
3818  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3819        nclsx = 0        nclsx = 0
3820        nclsy = 0        nclsy = 0
3821    
3822        do iv = 1,nviews        do iv = 1,nviews
3823           if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)  c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3824             good2(iv) = good2(iv) + mask_view(iv)*2**8
3825        enddo        enddo
3826    
3827        do icl=1,nclstr1        do icl=1,nclstr1
# Line 3396  c      good2=1!.true. Line 3830  c      good2=1!.true.
3830              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3831                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3832                 planex(nclsx) = ip                 planex(nclsx) = ip
3833                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3834                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3835                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3836                 do is=1,2                 do is=1,2
3837  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
# Line 3412  c$$$               print*,'xs(2,nclsx)   Line 3847  c$$$               print*,'xs(2,nclsx)  
3847              else                          !=== Y views              else                          !=== Y views
3848                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3849                 planey(nclsy) = ip                 planey(nclsy) = ip
3850                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3851                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3852                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3853                 do is=1,2                 do is=1,2
3854  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
# Line 3427  c$$$               print*,'ys(1,nclsy)   Line 3863  c$$$               print*,'ys(1,nclsy)  
3863  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3864              endif              endif
3865           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3866    
3867  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3868           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.23