/[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.16 by pam-fi, Thu Nov 30 17:04:27 2006 UTC revision 1.24 by pam-fi, Tue May 15 16:22:18 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 314  c$$$         enddo Line 349  c$$$         enddo
349           do i=1,5           do i=1,5
350              AL_GUESS(i)=AL(i)              AL_GUESS(i)=AL(i)
351           enddo           enddo
352    c         print*,'## guess: ',al
353    
354           do i=1,5           do i=1,5
355              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
# Line 476  c     $        rchi2best.lt.15..and. Line 512  c     $        rchi2best.lt.15..and.
512  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
513  *     angx   - Projected angle in x  *     angx   - Projected angle in x
514  *     angy   - Projected angle in y  *     angy   - Projected angle in y
515    *     bfx    - x-component of magnetci field
516    *     bfy    - y-component of magnetic field
517  *  *
518  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
519  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 514  c     $        rchi2best.lt.15..and. Line 552  c     $        rchi2best.lt.15..and.
552  *  *
553  *  *
554    
555        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
556    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c     01/02/2006 modified by elena vannuccini --> (2)  
 c     02/02/2006 modified by Elena Vannuccini --> (3)  
 c                (implemented new p.f.a.)  
 c     03/02/2006 modified by Elena Vannuccini --> (4)  
 c                (implemented variable resolution)  
 c*****************************************************  
557                
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
567        integer viewx,viewy        integer viewx,viewy
568        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
569        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
570          real angx,angy            !X-Y effective angle
571          real bfx,bfy              !X-Y b-field components
572    
573        real stripx,stripy        real stripx,stripy
574    
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 567  c      double precision xi_B,yi_B,zi_B Line 594  c      double precision xi_B,yi_B,zi_B
594        xPAM_B = 0.        xPAM_B = 0.
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
598    
599  *     -----------------  *     -----------------
600  *     CLUSTER X  *     CLUSTER X
601  *     -----------------  *     -----------------
602    
603        if(icx.ne.0)then        if(icx.ne.0)then
604           viewx = VIEW(icx)  
605           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
606           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
607           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
608                     resxPAM = RESXAV
609           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
610           if(PFAx.eq.'COG1')then  !(1)  *        --------------------------
611              stripx = stripx      !(1)  *        magnetic-field corrections
612              resxPAM = resxPAM    !(1)  *        --------------------------
613             angtemp  = ax
614             bfytemp  = bfy
615    *        /////////////////////////////////
616    *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!
617    *        *grvzkkjsdgjhhhgngbn###>:(
618    *        /////////////////////////////////
619    c         if(nplx.eq.6) angtemp = -1. * ax
620    c         if(nplx.eq.6) bfytemp = -1. * bfy
621             if(viewx.eq.12) angtemp = -1. * ax
622             if(viewx.eq.12) bfytemp = -1. * bfy
623             tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
624             angx     = 180.*atan(tgtemp)/acos(-1.)
625             stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
626    c$$$         print*,nplx,ax,bfy/10.
627    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
628    c$$$         print*,'========================'
629    c$$$         if(bfy.ne.0.)print*,viewx,'-x- '
630    c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ
631    *        --------------------------
632    
633    c$$$         print*,'--- x-cl ---'
634    c$$$         istart = INDSTART(ICX)
635    c$$$         istop  = TOTCLLENGTH
636    c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
637    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
638    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
639    c$$$         print*,INDMAX(icx)
640    c$$$         print*,cog(4,icx)
641    c$$$         print*,fbad_cog(4,icx)
642            
643    
644             if(PFAx.eq.'COG1')then
645    
646                stripx  = stripx
647                resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM
648    
649           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
650              stripx = stripx + cog(2,icx)              
651                stripx  = stripx + cog(2,icx)            
652                resxPAM = risx_cog(abs(angx))!TEMPORANEO              
653              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
654    
655             elseif(PFAx.eq.'COG3')then
656    
657                stripx  = stripx + cog(3,icx)            
658                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
659                resxPAM = resxPAM*fbad_cog(3,icx)
660    
661             elseif(PFAx.eq.'COG4')then
662    
663                stripx  = stripx + cog(4,icx)            
664                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
665                resxPAM = resxPAM*fbad_cog(4,icx)
666    
667           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
668  c            cog2 = cog(2,icx)  
669  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                          stripx  = stripx + pfaeta2(icx,angx)          
670  c            stripx = stripx + etacorr              resxPAM = risx_eta2(abs(angx))
671              stripx = stripx + pfaeta2(icx,angx)            !(3)              resxPAM = resxPAM*fbad_cog(2,icx)
             resxPAM = risx_eta2(angx)                       !   (4)  
672              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
673       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
674              resxPAM = resxPAM*fbad_cog(2,icx)  
675           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                        
676              stripx = stripx + pfaeta3(icx,angx)            !(3)  
677              resxPAM = risx_eta3(angx)                       !   (4)              stripx  = stripx + pfaeta3(icx,angx)          
678              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              resxPAM = risx_eta3(abs(angx))                      
679       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
680              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)            
681           elseif(PFAx.eq.'ETA4')then                         !(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
682              stripx = stripx + pfaeta4(icx,angx)            !(3)  
683              resxPAM = risx_eta4(angx)                       !   (4)           elseif(PFAx.eq.'ETA4')then                        
684              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
685       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)              stripx  = stripx + pfaeta4(icx,angx)            
686              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = risx_eta4(abs(angx))                      
687           elseif(PFAx.eq.'ETA')then                          !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
688              stripx = stripx + pfaeta(icx,angx)             !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              
689              resxPAM = ris_eta(icx,angx)                     !   (4)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
690              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
691       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)           elseif(PFAx.eq.'ETA')then  
692  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
693              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)              stripx  = stripx + pfaeta(icx,angx)            
694           elseif(PFAx.eq.'COG')then           !(2)  c            resxPAM = riseta(icx,angx)                    
695              stripx = stripx + cog(0,icx)     !(2)                      resxPAM = riseta(viewx,angx)                    
696              resxPAM = risx_cog(angx)                        !   (4)              resxPAM = resxPAM*fbad_eta(icx,angx)            
697              resxPAM = resxPAM*fbad_cog(0,icx)!(2)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              
698         $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
699    
700             elseif(PFAx.eq.'COG')then          
701    
702                stripx  = stripx + cog(0,icx)            
703                resxPAM = risx_cog(abs(angx))                    
704                resxPAM = resxPAM*fbad_cog(0,icx)
705    
706           else           else
707              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
708             endif
709    
710    
711    *     ======================================
712    *     temporary patch for saturated clusters
713    *     ======================================
714             if( nsatstrips(icx).gt.0 )then
715                stripx  = stripx + cog(4,icx)            
716                resxPAM = pitchX*1e-4/sqrt(12.)
717                cc=cog(4,icx)
718    c$$$            print*,icx,' *** ',cc
719    c$$$            print*,icx,' *** ',resxPAM
720           endif           endif
721    
722        endif        endif
 c      if(icy.eq.0.and.icx.ne.0)  
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
723                
724  *     -----------------  *     -----------------
725  *     CLUSTER Y  *     CLUSTER Y
726  *     -----------------  *     -----------------
727    
728        if(icy.ne.0)then        if(icy.ne.0)then
729    
730           viewy = VIEW(icy)           viewy = VIEW(icy)
731           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
732           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
733           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
734             stripy = float(MAXS(icy))
735    
736           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
737              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG) then
738       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
739         $              ,icx,icy
740                endif
741              goto 100              goto 100
742           endif           endif
743            *        --------------------------
744           stripy = float(MAXS(icy))  *        magnetic-field corrections
745           if(PFAy.eq.'COG1')then !(1)  *        --------------------------
746              stripy = stripy     !(1)           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
747              resyPAM = resyPAM   !(1)           angy    = 180.*atan(tgtemp)/acos(-1.)
748             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
749    c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
750    c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
751    *        --------------------------
752            
753    c$$$         print*,'--- y-cl ---'
754    c$$$         istart = INDSTART(ICY)
755    c$$$         istop  = TOTCLLENGTH
756    c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
757    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
758    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
759    c$$$         print*,INDMAX(icy)
760    c$$$         print*,cog(4,icy)
761    c$$$         print*,fbad_cog(4,icy)
762    
763             if(PFAy.eq.'COG1')then
764    
765                stripy  = stripy    
766                resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  
767    
768           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
769              stripy = stripy + cog(2,icy)  
770                stripy  = stripy + cog(2,icy)
771                resyPAM = risy_cog(abs(angy))!TEMPORANEO
772              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
773    
774             elseif(PFAy.eq.'COG3')then
775    
776                stripy  = stripy + cog(3,icy)
777                resyPAM = risy_cog(abs(angy))!TEMPORANEO
778                resyPAM = resyPAM*fbad_cog(3,icy)
779    
780             elseif(PFAy.eq.'COG4')then
781    
782                stripy  = stripy + cog(4,icy)
783                resyPAM = risy_cog(abs(angy))!TEMPORANEO
784                resyPAM = resyPAM*fbad_cog(4,icy)
785    
786           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
787  c            cog2 = cog(2,icy)  
788  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)              stripy  = stripy + pfaeta2(icy,angy)          
789  c            stripy = stripy + etacorr              resyPAM = risy_eta2(abs(angy))              
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
790              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
791              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
792       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
793           elseif(PFAy.eq.'ETA3')then                         !(3)  
794              stripy = stripy + pfaeta3(icy,angy)            !(3)           elseif(PFAy.eq.'ETA3')then                      
795              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
796              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              stripy  = stripy + pfaeta3(icy,angy)
797       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)              resyPAM = resyPAM*fbad_cog(3,icy)  
798           elseif(PFAy.eq.'ETA4')then                         !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)
799              stripy = stripy + pfaeta4(icy,angy)            !(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
800              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
801              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA4')then  
802       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
803           elseif(PFAy.eq.'ETA')then                          !(3)              stripy  = stripy + pfaeta4(icy,angy)
804              stripy = stripy + pfaeta(icy,angy)             !(3)              resyPAM = resyPAM*fbad_cog(4,icy)
805              resyPAM = ris_eta(icy,angy)                     !   (4)              if(DEBUG.and.fbad_cog(4,icy).ne.1)
806  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
807              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
808              if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA')then
809       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
810                stripy  = stripy + pfaeta(icy,angy)
811    c            resyPAM = riseta(icy,angy)  
812                resyPAM = riseta(viewy,angy)  
813                resyPAM = resyPAM*fbad_eta(icy,angy)
814                if(DEBUG.and.fbad_cog(2,icy).ne.1)
815         $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
816    
817           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
818              stripy = stripy + cog(0,icy)              
819              resyPAM = risy_cog(angy)                        !   (4)              stripy  = stripy + cog(0,icy)            
820  c            resyPAM = ris_eta(icy,angy)                    !   (4)              resyPAM = risy_cog(abs(angy))
821              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
822    
823           else           else
824              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
825             endif
826    
827    
828    *     ======================================
829    *     temporary patch for saturated clusters
830    *     ======================================
831             if( nsatstrips(icy).gt.0 )then
832                stripy  = stripy + cog(4,icy)            
833                resyPAM = pitchY*1e-4/sqrt(12.)
834                cc=cog(4,icy)
835    c$$$            print*,icy,' *** ',cc
836    c$$$            print*,icy,' *** ',resyPAM
837           endif           endif
838    
839    
840        endif        endif
841    
842          c$$$      print*,'## stripx,stripy ',stripx,stripy
843    
844  c===========================================================  c===========================================================
845  C     COUPLE  C     COUPLE
846  C===========================================================  C===========================================================
# Line 697  c     (xi,yi,zi) = mechanical coordinate Line 851  c     (xi,yi,zi) = mechanical coordinate
851  c------------------------------------------------------------------------  c------------------------------------------------------------------------
852           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
853       $        .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...
854              print*,'xyz_PAM (couple):',              if(DEBUG) then
855       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
856         $              ' WARNING: false X strip: strip ',stripx
857                endif
858           endif           endif
859           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
860           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
# Line 790  c            print*,'X-singlet ',icx,npl Line 946  c            print*,'X-singlet ',icx,npl
946  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...
947              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
948       $           .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...
949                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG) then
950       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
951         $                 ' WARNING: false X strip: strip ',stripx
952                   endif
953              endif              endif
954              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
955    
# Line 813  c            print*,'X-cl ',icx,stripx,' Line 971  c            print*,'X-cl ',icx,stripx,'
971  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
972    
973           else           else
974                if(DEBUG) then
975              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
976              print *,'icx = ',icx                 print *,'icx = ',icx
977              print *,'icy = ',icy                 print *,'icy = ',icy
978                endif
979              goto 100              goto 100
980                            
981           endif           endif
# Line 881  c--------------------------------------- Line 1040  c---------------------------------------
1040  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1041    
1042        else        else
1043                       if(DEBUG) then
1044           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1045           print *,'icx = ',icx              print *,'icx = ',icx
1046           print *,'icy = ',icy              print *,'icy = ',icy
1047                         endif
1048        endif        endif
1049                    
1050    
1051    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1052    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1053    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1054    
1055   100  continue   100  continue
1056        end        end
1057    
1058    ************************************************************************
1059    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1060    *     (done to be called from c/c++)
1061    ************************************************************************
1062    
1063          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1064    
1065          include 'commontracker.f'
1066          include 'level1.f'
1067          include 'common_mini_2.f'
1068          include 'common_xyzPAM.f'
1069          include 'common_mech.f'
1070          include 'calib.f'
1071          
1072    *     flag to chose PFA
1073    c$$$      character*10 PFA
1074    c$$$      common/FINALPFA/PFA
1075    
1076          integer icx,icy           !X-Y cluster ID
1077          integer sensor
1078          character*4 PFAx,PFAy     !PFA to be used
1079          real ax,ay                !X-Y geometric angle
1080          real bfx,bfy              !X-Y b-field components
1081    
1082          ipx=0
1083          ipy=0      
1084          
1085    c$$$      PFAx = 'COG4'!PFA
1086    c$$$      PFAy = 'COG4'!PFA
1087          
1088          call idtoc(pfaid,PFAx)
1089          call idtoc(pfaid,PFAy)
1090    
1091          call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1092    
1093    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1094          
1095          if(icx.ne.0.and.icy.ne.0)then
1096    
1097             ipx=npl(VIEW(icx))
1098             ipy=npl(VIEW(icy))
1099             if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1100         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1101         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1102    
1103             xgood(ip) = 1.
1104             ygood(ip) = 1.
1105             resx(ip)  = resxPAM
1106             resy(ip)  = resyPAM
1107    
1108             xm(ip) = xPAM
1109             ym(ip) = yPAM
1110             zm(ip) = zPAM
1111             xm_A(ip) = 0.
1112             ym_A(ip) = 0.
1113             xm_B(ip) = 0.
1114             ym_B(ip) = 0.
1115    
1116    c         zv(ip) = zPAM
1117    
1118          elseif(icx.eq.0.and.icy.ne.0)then
1119    
1120             ipy=npl(VIEW(icy))
1121             if((nplanes-ipy+1).ne.ip)
1122         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1123         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1124    
1125             xgood(ip) = 0.
1126             ygood(ip) = 1.
1127             resx(ip)  = 1000.
1128             resy(ip)  = resyPAM
1129    
1130             xm(ip) = -100.
1131             ym(ip) = -100.
1132             zm(ip) = (zPAM_A+zPAM_B)/2.
1133             xm_A(ip) = xPAM_A
1134             ym_A(ip) = yPAM_A
1135             xm_B(ip) = xPAM_B
1136             ym_B(ip) = yPAM_B
1137    
1138    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1139            
1140          elseif(icx.ne.0.and.icy.eq.0)then
1141    
1142             ipx=npl(VIEW(icx))
1143             if((nplanes-ipx+1).ne.ip)
1144         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1145         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1146    
1147             xgood(ip) = 1.
1148             ygood(ip) = 0.
1149             resx(ip)  = resxPAM
1150             resy(ip)  = 1000.
1151    
1152             xm(ip) = -100.
1153             ym(ip) = -100.
1154             zm(ip) = (zPAM_A+zPAM_B)/2.
1155             xm_A(ip) = xPAM_A
1156             ym_A(ip) = yPAM_A
1157             xm_B(ip) = xPAM_B
1158             ym_B(ip) = yPAM_B
1159            
1160    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1161    
1162          else
1163    
1164             il = 2
1165             if(lad.ne.0)il=lad
1166             is = 1
1167             if(sensor.ne.0)is=sensor
1168    c         print*,nplanes-ip+1,il,is
1169    
1170             xgood(ip) = 0.
1171             ygood(ip) = 0.
1172             resx(ip)  = 1000.
1173             resy(ip)  = 1000.
1174    
1175             xm(ip) = -100.
1176             ym(ip) = -100.          
1177             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1178             xm_A(ip) = 0.
1179             ym_A(ip) = 0.
1180             xm_B(ip) = 0.
1181             ym_B(ip) = 0.
1182    
1183    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1184    
1185          endif
1186    
1187          if(DEBUG)then
1188    c         print*,'----------------------------- track coord'
1189    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1190             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1191         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1192         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1193    c$$$         print*,'-----------------------------'
1194          endif
1195          end
1196    
1197  ********************************************************************************  ********************************************************************************
1198  ********************************************************************************  ********************************************************************************
# Line 966  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1268  c         print*,'A-(',xPAM_A,yPAM_A,')
1268           endif                   endif        
1269    
1270           distance=           distance=
1271       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1272    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1273           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1274    
1275  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 991  c$$$         print*,' resolution ',re Line 1294  c$$$         print*,' resolution ',re
1294  *     ----------------------  *     ----------------------
1295                    
1296           distance=           distance=
1297       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1298       $        +       $       +
1299       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1300    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1301    c$$$     $        +
1302    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1303           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1304    
1305  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1002  c$$$         print*,' resolution ',resxP Line 1308  c$$$         print*,' resolution ',resxP
1308                    
1309        else        else
1310                    
1311           print*  c         print*
1312       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1313           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1314           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 '
1315       $        ,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
1316        endif          endif  
1317    
1318        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1074  c--------------------------------------- Line 1380  c---------------------------------------
1380                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1381       $              .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...
1382  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...
1383                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1384       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1385                 endif                 endif
1386                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1387                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1230  c      include 'common_analysis.f' Line 1536  c      include 'common_analysis.f'
1536        is_cp=0        is_cp=0
1537        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1538        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1539        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1540    
1541        return        return
1542        end        end
# Line 1378  c      include 'level1.f' Line 1684  c      include 'level1.f'
1684  *     ----------------------------------------------------  *     ----------------------------------------------------
1685  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1686  *     ----------------------------------------------------  *     ----------------------------------------------------
1687           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1688              cl_single(icx)=0              cl_single(icx)=0
1689              goto 10              goto 10
1690           endif           endif
# Line 1428  c     endif Line 1734  c     endif
1734  *     ----------------------------------------------------  *     ----------------------------------------------------
1735  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1736  *     ----------------------------------------------------  *     ----------------------------------------------------
1737              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1738                 cl_single(icy)=0                 cl_single(icy)=0
1739                 goto 20                 goto 20
1740              endif              endif
# Line 1474  c     endif Line 1780  c     endif
1780  *     charge correlation  *     charge correlation
1781  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1782    
1783                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1784       $              .and.       $              .and.
1785       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1786       $              .and.       $              .and.
1787       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1788       $              .and.       $              .and.
1789       $              .true.)then       $              .true.)then
1790    
1791                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1792       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1793                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1794    
1795  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1796    
1797                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1798       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1799                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1800                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1623  c      double precision xm3,ym3,zm3 Line 1929  c      double precision xm3,ym3,zm3
1929                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1930                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1931  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1932                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1933                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1934                 xm1=xPAM                 xm1=xPAM
1935                 ym1=yPAM                 ym1=yPAM
1936                 zm1=zPAM                                   zm1=zPAM                  
# Line 1639  c     print*,'***',is1,xm1,ym1,zm1 Line 1946  c     print*,'***',is1,xm1,ym1,zm1
1946                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1947  c                        call xyz_PAM  c                        call xyz_PAM
1948  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1949    c                        call xyz_PAM
1950    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1951                          call xyz_PAM                          call xyz_PAM
1952       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1953                          xm2=xPAM                          xm2=xPAM
1954                          ym2=yPAM                          ym2=yPAM
1955                          zm2=zPAM                          zm2=zPAM
# Line 1704  c$$$ Line 2013  c$$$
2013                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2014  c                                 call xyz_PAM  c                                 call xyz_PAM
2015  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2016    c                                 call xyz_PAM
2017    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2018                                   call xyz_PAM                                   call xyz_PAM
2019       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2020         $                                ,0.,0.,0.,0.)
2021                                   xm3=xPAM                                   xm3=xPAM
2022                                   ym3=yPAM                                   ym3=yPAM
2023                                   zm3=zPAM                                   zm3=zPAM
# Line 2177  c     print*,'check cp_used' Line 2489  c     print*,'check cp_used'
2489           enddo           enddo
2490  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2491           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2492           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2493                    
2494  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2495  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
# Line 2254  c$$$     $           ,(tr_cloud(iii),iii Line 2566  c$$$     $           ,(tr_cloud(iii),iii
2566  **************************************************  **************************************************
2567    
2568        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2569    
2570        include 'commontracker.f'        include 'commontracker.f'
2571        include 'level1.f'        include 'level1.f'
# Line 2264  c*************************************** Line 2573  c***************************************
2573        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2574        include 'common_mini_2.f'        include 'common_mini_2.f'
2575        include 'common_mech.f'        include 'common_mech.f'
2576  c      include 'momanhough_init.f'  
2577    
2578    
2579  *     output flag  *     output flag
# Line 2459  c$$$            if(AL_INI(5).gt.defmax)g Line 2768  c$$$            if(AL_INI(5).gt.defmax)g
2768  *                                   *************************  *                                   *************************
2769  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2770  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2771    c                                    call xyz_PAM(icx,icy,is, !(1)
2772    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2773                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2774       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2775  *                                   *************************  *                                   *************************
2776  *                                   -----------------------------  *                                   -----------------------------
2777                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 2531  c                                 goto 8 Line 2842  c                                 goto 8
2842                                                                
2843                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2844                                                                
 c$$$                              ndof=0                                  
2845                                do ip=1,nplanes                                do ip=1,nplanes
2846  c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2847                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2848                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2849                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 2554  c$$$     $                               Line 2862  c$$$     $                              
2862                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2863                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2864       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2865                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2866         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2867                                        LADDER_STORE(nplanes-ip+1,ntracks)
2868         $                                   = LADDER(
2869         $                                   clx(ip,icp_cp(
2870         $                                   cp_match(ip,hit_plane(ip)
2871         $                                   ))));
2872                                   else                                   else
2873                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2874                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2875                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2876                                   endif                                   endif
2877                                     BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2878                                     BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2879                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(nplanes-ip+1,ntracks)=0
2880                                   do i=1,5                                   do i=1,5
2881                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2882                                   enddo                                   enddo
2883                                enddo                                enddo
2884                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2885                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2886                                                                
2887  *     --------------------------------  *     --------------------------------
# Line 2590  c$$$  rchi2=chi2/dble(ndof) Line 2905  c$$$  rchi2=chi2/dble(ndof)
2905           return           return
2906        endif        endif
2907                
2908    c$$$      if(DEBUG)then
2909    c$$$         print*,'****** TRACK CANDIDATES ***********'
2910    c$$$         print*,'#         R. chi2        RIG'
2911    c$$$         do i=1,ntracks
2912    c$$$            print*,i,' --- ',rchi2_store(i),' --- '
2913    c$$$     $           ,1./abs(AL_STORE(5,i))
2914    c$$$         enddo
2915    c$$$         print*,'***********************************'
2916    c$$$      endif
2917        if(DEBUG)then        if(DEBUG)then
2918           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
2919           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
2920           do i=1,ntracks          do i=1,ntracks
2921              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
2922       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
2923           enddo              ndof=ndof           !(1)
2924           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
2925         $           +int(ygood_store(ii,i)) !(1)
2926              enddo                 !(1)
2927              print*,i,' --- ',rchi2_store(i),' --- '
2928         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
2929            enddo
2930            print*,'*****************************************'
2931        endif        endif
2932                
2933                
# Line 2616  c$$$  rchi2=chi2/dble(ndof) Line 2946  c$$$  rchi2=chi2/dble(ndof)
2946    
2947        subroutine refine_track(ibest)        subroutine refine_track(ibest)
2948    
 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******************************************************  
2949    
2950        include 'commontracker.f'        include 'commontracker.f'
2951        include 'level1.f'        include 'level1.f'
# Line 2628  c*************************************** Line 2953  c***************************************
2953        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2954        include 'common_mini_2.f'        include 'common_mini_2.f'
2955        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
2956        include 'calib.f'        include 'calib.f'
2957    
   
2958  *     flag to chose PFA  *     flag to chose PFA
2959        character*10 PFA        character*10 PFA
2960        common/FINALPFA/PFA        common/FINALPFA/PFA
2961    
2962          real k(6)
2963          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
2964    
2965          real xp,yp,zp
2966          real xyzp(3),bxyz(3)
2967          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
2968    
2969  *     =================================================  *     =================================================
2970  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
2971  *                          and  *                          and
# Line 2645  c      include 'level1.f' Line 2974  c      include 'level1.f'
2974        call track_init        call track_init
2975        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2976    
2977             xP=XV_STORE(nplanes-ip+1,ibest)
2978             yP=YV_STORE(nplanes-ip+1,ibest)
2979             zP=ZV_STORE(nplanes-ip+1,ibest)
2980             call gufld(xyzp,bxyz)
2981             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
2982             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
2983    c$$$  bxyz(1)=0
2984    c$$$         bxyz(2)=0
2985    c$$$         bxyz(3)=0
2986  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
2987  *     -------------------------------------------------  *     -------------------------------------------------
2988  *     If the plane has been already included, it just  *     If the plane has been already included, it just
# Line 2665  c      include 'level1.f' Line 3003  c      include 'level1.f'
3003       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3004              icx=clx(ip,icp)              icx=clx(ip,icp)
3005              icy=cly(ip,icp)              icy=cly(ip,icp)
3006    c            call xyz_PAM(icx,icy,is,
3007    c     $           PFA,PFA,
3008    c     $           AXV_STORE(nplanes-ip+1,ibest),
3009    c     $           AYV_STORE(nplanes-ip+1,ibest))
3010              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3011       $           PFA,PFA,       $           PFA,PFA,
3012       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3013       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3014  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3015  c$$$  $              'COG2','COG2',       $           bxyz(2)
3016  c$$$  $              0.,       $           )
3017  c$$$  $              0.)  
3018              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3019              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3020              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 2682  c$$$  $              0.) Line 3023  c$$$  $              0.)
3023              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3024              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3025    
3026  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3027              dedxtrk_x(nplanes-ip+1)=dedx(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)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  
3028                            
3029  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3030  *     -------------------------------------------------  *     -------------------------------------------------
# Line 2699  c            dedxtrk(nplanes-ip+1) = (de Line 3039  c            dedxtrk(nplanes-ip+1) = (de
3039                                
3040  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3041  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
             xP=XV_STORE(nplanes-ip+1,ibest)  
             yP=YV_STORE(nplanes-ip+1,ibest)  
             zP=ZV_STORE(nplanes-ip+1,ibest)  
3042              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3043  *     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
3044              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3045    
3046                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3047                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3048  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3049    
3050              if(DEBUG)then              if(DEBUG)then
# Line 2741  c     $              cl_used(icy).eq.1.o Line 3081  c     $              cl_used(icy).eq.1.o
3081  *            *          
3082                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3083       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3084       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3085       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3086         $              bxyz(1),
3087         $              bxyz(2)
3088         $              )
3089                                
3090                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3091                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3092                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3093                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
3094       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3095                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3096                    xmm = xPAM                    xmm = xPAM
3097                    ymm = yPAM                    ymm = yPAM
# Line 2758  c     $              'ETA2','ETA2', Line 3100  c     $              'ETA2','ETA2',
3100                    rymm = resyPAM                    rymm = resyPAM
3101                    distmin = distance                    distmin = distance
3102                    idm = id                                      idm = id                  
3103  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3104                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3105                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3106                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3107         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3108                 endif                 endif
3109   1188          continue   1188          continue
3110              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3111              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3112                if(distmin.le.clincnewc)then     !QUIQUI              
3113  *              -----------------------------------  *              -----------------------------------
3114                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3115                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3116                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3117                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3118                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3119                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3120                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3121  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3122                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3123  *              -----------------------------------  *              -----------------------------------
3124                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3125                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG)print*,'%%%% included couple ',idm
3126       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3127                 goto 133         !next plane                 goto 133         !next plane
3128              endif              endif
3129  *     ================================================  *     ================================================
# Line 2812  c            if(DEBUG)print*,'>>>> try t Line 3156  c            if(DEBUG)print*,'>>>> try t
3156  c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used  c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
3157                 if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)                 if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
3158  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3159    c               call xyz_PAM(icx,0,ist,
3160    c     $              PFA,PFA,
3161    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3162                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3163       $              PFA,PFA,       $              PFA,PFA,
3164       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3165         $              bxyz(1),
3166         $              bxyz(2)
3167         $              )              
3168                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3169                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3170                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3171       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3172                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3173                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3174                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2831  c     $              'ETA2','ETA2', Line 3180  c     $              'ETA2','ETA2',
3180                    rymm = resyPAM                    rymm = resyPAM
3181                    distmin = distance                    distmin = distance
3182                    iclm = icx                    iclm = icx
3183  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3184                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3185                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3186                 endif                                   endif                  
3187  11881          continue  11881          continue
# Line 2840  c                  dedxmm = dedx(icx) !( Line 3189  c                  dedxmm = dedx(icx) !(
3189  c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used  c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
3190                 if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)                 if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
3191  *                                              !jump to the next couple  *                                              !jump to the next couple
3192    c               call xyz_PAM(0,icy,ist,
3193    c     $              PFA,PFA,
3194    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3195                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3196       $              PFA,PFA,       $              PFA,PFA,
3197       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3198         $              bxyz(1),
3199         $              bxyz(2)
3200         $              )
3201                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3202                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3203                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3204       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3205                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3206                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3207                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2859  c     $              'ETA2','ETA2', Line 3213  c     $              'ETA2','ETA2',
3213                    rymm = resyPAM                    rymm = resyPAM
3214                    distmin = distance                    distmin = distance
3215                    iclm = icy                    iclm = icy
3216  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3217                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3218                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3219                 endif                                   endif                  
3220  11882          continue  11882          continue
3221              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3222  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3223    c            print*,'## ncls(',ip,') ',ncls(ip)
3224              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3225                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3226  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
# Line 2874  c               if(cl_used(icl).eq.1.or. Line 3229  c               if(cl_used(icl).eq.1.or.
3229       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3230                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3231                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3232       $                 PFA,PFA,       $                 PFA,PFA,
3233       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3234         $                 bxyz(1),
3235         $                 bxyz(2)
3236         $                 )
3237                 else                         !<---- Y view                 else                         !<---- Y view
3238                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3239       $                 PFA,PFA,       $                 PFA,PFA,
3240       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3241         $                 bxyz(1),
3242         $                 bxyz(2)
3243         $                 )
3244                 endif                 endif
3245    
3246                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3247                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3248                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3249       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance,'<',distmin,' ?'
3250                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3251                      if(DEBUG)print*,'YES'
3252                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3253                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3254                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 2899  c     $                 'ETA2','ETA2', Line 3259  c     $                 'ETA2','ETA2',
3259                    rymm = resyPAM                    rymm = resyPAM
3260                    distmin = distance                      distmin = distance  
3261                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3262                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3263                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3264                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3265                    else          !<---- Y view                    else          !<---- Y view
3266                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3267                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3268                    endif                    endif
3269                 endif                                   endif                  
3270  18882          continue  18882          continue
3271              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3272    c            print*,'## distmin ', distmin,' clinc ',clinc
3273    
3274              if(distmin.le.clinc)then                    c     QUIQUI------------
3275                  c     anche qui: non ci vuole clinc???
3276                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3277                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3278                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3279                    resx(nplanes-ip+1)=rxmm       $                 20*
3280                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3281  c                  if(.true.)print*,'%%%% included X-cl ',iclm                 else if(mod(VIEW(iclm),2).ne.0)then
3282       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    clincnew=
3283       $                 ,', norm.dist.= ',distmin       $                 10*
3284       $                 ,', cut ',clinc,' )'       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3285                 else                 endif
3286                    YGOOD(nplanes-ip+1)=1.  c     QUIQUI------------
3287                    resy(nplanes-ip+1)=rymm                
3288                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 if(distmin.le.clincnew)then   !QUIQUI
3289  c                  if(.true.)print*,'%%%% included Y-cl ',iclm  c     if(distmin.le.clinc)then          !QUIQUI          
3290       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    
3291       $                 ,', norm.dist.= ', distmin                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3292       $                 ,', cut ',clinc,' )'  *     ----------------------------
3293    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3294                      if(mod(VIEW(iclm),2).eq.0)then
3295                         XGOOD(nplanes-ip+1)=1.
3296                         resx(nplanes-ip+1)=rxmm
3297                         if(DEBUG)print*,'%%%% included X-cl ',iclm
3298         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3299         $                    ,', dist.= ',distmin
3300         $                    ,', cut ',clinc,' )'
3301                      else
3302                         YGOOD(nplanes-ip+1)=1.
3303                         resy(nplanes-ip+1)=rymm
3304                         if(DEBUG)print*,'%%%% included Y-cl ',iclm
3305         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3306         $                    ,', dist.= ', distmin
3307         $                    ,', cut ',clinc,' )'
3308                      endif
3309    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3310    *     ----------------------------
3311                      xm_A(nplanes-ip+1) = xmm_A
3312                      ym_A(nplanes-ip+1) = ymm_A
3313                      xm_B(nplanes-ip+1) = xmm_B
3314                      ym_B(nplanes-ip+1) = ymm_B
3315                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3316                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3317                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3318    *     ----------------------------
3319                 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)  
 *              ----------------------------  
3320              endif              endif
3321           endif           endif
3322   133     continue   133     continue
# Line 2962  c              dedxtrk(nplanes-ip+1) = d Line 3335  c              dedxtrk(nplanes-ip+1) = d
3335  *                                                 *  *                                                 *
3336  *                                                 *  *                                                 *
3337  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3338  *  *
3339        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3340    
3341        include 'commontracker.f'        include 'commontracker.f'
3342        include 'level1.f'        include 'level1.f'
3343        include 'common_momanhough.f'        include 'common_momanhough.f'
3344  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
   
3345    
3346        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3347    
# Line 2984  c      include 'level1.f' Line 3351  c      include 'level1.f'
3351              if(id.ne.0)then              if(id.ne.0)then
3352                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3353                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3354  c               cl_used(iclx)=1  !tag used clusters                 cl_used(iclx)=ntrk  !tag used clusters
3355  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)  
3356              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3357  c               cl_used(icl)=1   !tag used clusters                 cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3358              endif              endif
3359                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3360  *     -----------------------------  *     -----------------------------
3361  *     remove the couple from clouds  *     remove the couple from clouds
3362  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3056  c               endif Line 3417  c               endif
3417        include 'level1.f'        include 'level1.f'
3418        include 'common_momanhough.f'        include 'common_momanhough.f'
3419        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3420    
3421    *     ---------------------------------
3422    *     variables initialized from level1
3423    *     ---------------------------------
3424        do i=1,nviews        do i=1,nviews
3425           good2(i)=good1(i)           good2(i)=good1(i)
3426             do j=1,nva1_view
3427                vkflag(i,j)=1
3428                if(cnnev(i,j).le.0)then
3429                   vkflag(i,j)=cnnev(i,j)
3430                endif
3431             enddo
3432        enddo        enddo
3433    *     ----------------
3434    *     level2 variables
3435    *     ----------------
3436        NTRK = 0        NTRK = 0
3437        do it=1,NTRKMAX        do it=1,NTRKMAX
3438           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3073  c      include 'level1.f' Line 3443  c      include 'level1.f'
3443              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3444              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3445              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3446                TAILX_nt(IP,IT) = 0
3447                TAILY_nt(IP,IT) = 0
3448                XBAD(IP,IT) = 0
3449                YBAD(IP,IT) = 0
3450              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3451              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3452                LS(IP,IT) = 0
3453              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3454              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3455              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
# Line 3207  c      include 'level1.f' Line 3582  c      include 'level1.f'
3582    
3583            
3584        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3585        include 'level1.f'        include 'level1.f'
3586        include 'common_momanhough.f'        include 'common_momanhough.f'
3587        include 'level2.f'        include 'level2.f'
3588        include 'common_mini_2.f'        include 'common_mini_2.f'
3589        real sinth,phi,pig              include 'calib.f'
3590    
3591          character*10 PFA
3592          common/FINALPFA/PFA
3593    
3594          real sinth,phi,pig
3595          integer ssensor,sladder
3596        pig=acos(-1.)        pig=acos(-1.)
3597    
3598    *     -------------------------------------
3599        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3600        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3601    *     -------------------------------------
3602        phi   = al(4)                  phi   = al(4)          
3603        sinth = al(3)                    sinth = al(3)            
3604        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3230  c      include 'level1.f' Line 3611  c      include 'level1.f'
3611       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3612        al(4) = phi                      al(4) = phi              
3613        al(3) = sinth                    al(3) = sinth            
   
3614        do i=1,5        do i=1,5
3615           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3616           do j=1,5           do j=1,5
3617              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3618           enddo           enddo
3619        enddo        enddo
3620          *     -------------------------------------      
3621        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3622           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3623           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3246  c      include 'level1.f' Line 3626  c      include 'level1.f'
3626           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3627           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3628           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3629             TAILX_nt(IP,ntr) = 0.
3630             TAILY_nt(IP,ntr) = 0.
3631           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3632           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3633           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3634           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3635           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3636           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3637           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3638         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3639         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3640         $        1. )
3641    c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3642             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3643             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3644        
3645           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3646             ay   = ayv_nt(ip,ntr)
3647             bfx  = BX_STORE(ip,IDCAND)
3648             bfy  = BY_STORE(ip,IDCAND)
3649             if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3650             if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3651             tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3652             angx     = 180.*atan(tgtemp)/acos(-1.)
3653             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3654             angy    = 180.*atan(tgtemp)/acos(-1.)
3655            
3656    c         print*,'* ',ip,bfx,bfy,angx,angy
3657    
3658             id  = CP_STORE(ip,IDCAND) ! couple id
3659           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3660             ssensor = -1
3661             sladder = -1
3662             ssensor = SENSOR_STORE(ip,IDCAND)
3663             sladder = LADDER_STORE(ip,IDCAND)
3664             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3665             LS(IP,ntr)      = ssensor+10*sladder
3666    
3667           if(id.ne.0)then           if(id.ne.0)then
3668    c           >>> is a couple
3669              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3670              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3671  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3672    c$$$            if(is_cp(id).ne.ssensor)
3673    c$$$     $           print*,'ERROR is sensor assignment (couple)'
3674    c$$$     $           ,is_cp(id),ssensor
3675    c$$$            if(LADDER(clx(nplanes-ip+1,icp_cp(id))).ne.sladder)
3676    c$$$     $           print*,'ERROR is ladder assignment (couple)'
3677    c$$$     $           ,LADDER(clx(nplanes-ip+1,icp_cp(id))),sladder
3678                
3679                nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3680                nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3681                xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
3682                ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
3683    
3684                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3685         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3686                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3687         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3688    
3689           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3690              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  c           >>> is a singlet
3691              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl  c$$$            if(LADDER(icl).ne.sladder)
3692  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  c$$$     $           print*,'ERROR is ladder assignment (single)'
3693    c$$$     $           ,LADDER(icl),sladder
3694                if(mod(VIEW(icl),2).eq.0)then
3695                   cltrx(ip,ntr)=icl
3696                   nnnnn = npfastrips(icl,PFA,angx)
3697                   xbad(ip,ntr) = nbadstrips(nnnnn,icl)
3698                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3699                elseif(mod(VIEW(icl),2).eq.1)then
3700                   cltry(ip,ntr)=icl
3701                   nnnnn = npfastrips(icl,PFA,angy)
3702                   ybad(ip,ntr) = nbadstrips(nnnnn,icl)
3703                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3704                endif
3705           endif                     endif          
3706    
3707        enddo        enddo
3708    
3709    
3710    c$$$      print*,(xgood(i),i=1,6)
3711    c$$$      print*,(ygood(i),i=1,6)
3712    c$$$      print*,(ls(i,ntr),i=1,6)
3713    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3714    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3715    c$$$      print*,'-----------------------'
3716    
3717        end        end
3718    
3719        subroutine fill_level2_siglets        subroutine fill_level2_siglets
# Line 3280  c            print*,ip,' ',cltrx(ip,ntr) Line 3725  c            print*,ip,' ',cltrx(ip,ntr)
3725  *     -------------------------------------------------------  *     -------------------------------------------------------
3726    
3727        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3728        include 'calib.f'        include 'calib.f'
3729        include 'level1.f'        include 'level1.f'
3730        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 3288  c      include 'level1.f' Line 3732  c      include 'level1.f'
3732        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3733    
3734  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3735        nclsx = 0        nclsx = 0
3736        nclsy = 0        nclsy = 0
3737    
# Line 3302  c      good2=1!.true. Line 3745  c      good2=1!.true.
3745              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3746                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3747                 planex(nclsx) = ip                 planex(nclsx) = ip
3748                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3749                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3750                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3751                 do is=1,2                 do is=1,2
3752  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3753                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3754                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3755                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3756                 enddo                 enddo
3757  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3317  c$$$               print*,'xs(2,nclsx)   Line 3762  c$$$               print*,'xs(2,nclsx)  
3762              else                          !=== Y views              else                          !=== Y views
3763                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3764                 planey(nclsy) = ip                 planey(nclsy) = ip
3765                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3766                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3767                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3768                 do is=1,2                 do is=1,2
3769  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3770                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3771                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3772                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3773                 enddo                 enddo
3774  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 3331  c$$$               print*,'ys(1,nclsy)   Line 3778  c$$$               print*,'ys(1,nclsy)  
3778  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3779              endif              endif
3780           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3781    
3782  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3783           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.23