/[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.14 by pam-fi, Tue Nov 21 14:00:40 2006 UTC revision 1.27 by pam-fi, Fri Aug 17 14:36:05 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 233  c$$$            endif Line 266  c$$$            endif
266  c$$$         enddo  c$$$         enddo
267  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
268    
269           rchi2best=1000000000.  *     -------------------------------------------------------
270           ndofbest=0             !(1)  *     order track-candidates according to:
271    *     1st) decreasing n.points
272    *     2nd) increasing chi**2
273    *     -------------------------------------------------------
274             rchi2best=1000000000.
275             ndofbest=0            
276           do i=1,ntracks           do i=1,ntracks
277             if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0              
278       $          RCHI2_STORE(i).gt.0)then             do ii=1,nplanes    
279               ndof=0             !(1)               ndof=ndof        
280               do ii=1,nplanes    !(1)       $            +int(xgood_store(ii,i))
281                 ndof=ndof        !(1)       $            +int(ygood_store(ii,i))
282       $              +int(xgood_store(ii,i)) !(1)             enddo              
283       $              +int(ygood_store(ii,i)) !(1)             if(ndof.gt.ndofbest)then
284               enddo              !(1)               ibest=i
285               if(ndof.ge.ndofbest)then !(1)               rchi2best=RCHI2_STORE(i)
286                 ndofbest=ndof    
287               elseif(ndof.eq.ndofbest)then
288                 if(RCHI2_STORE(i).lt.rchi2best.and.
289         $            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    
297    c$$$         rchi2best=1000000000.
298    c$$$         ndofbest=0             !(1)
299    c$$$         do i=1,ntracks
300    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
301    c$$$     $          RCHI2_STORE(i).gt.0)then
302    c$$$             ndof=0             !(1)
303    c$$$             do ii=1,nplanes    !(1)
304    c$$$               ndof=ndof        !(1)
305    c$$$     $              +int(xgood_store(ii,i)) !(1)
306    c$$$     $              +int(ygood_store(ii,i)) !(1)
307    c$$$             enddo              !(1)
308    c$$$             if(ndof.ge.ndofbest)then !(1)
309    c$$$               ibest=i
310    c$$$               rchi2best=RCHI2_STORE(i)
311    c$$$               ndofbest=ndof    !(1)
312    c$$$             endif              !(1)
313    c$$$           endif
314    c$$$         enddo
315    
316           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
317  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
318  *     The best track candidate (ibest) is selected and a new fitting is performed.  *     The best track candidate (ibest) is selected and a new fitting is performed.
# Line 269  c$$$         if(ibest.eq.0)goto 880 !>> Line 331  c$$$         if(ibest.eq.0)goto 880 !>>
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 285  c$$$         if(ibest.eq.0)goto 880 !>> Line 349  c$$$         if(ibest.eq.0)goto 880 !>>
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 297  c$$$         if(ibest.eq.0)goto 880 !>> Line 362  c$$$         if(ibest.eq.0)goto 880 !>>
362           iprint=0           iprint=0
363  c         if(DEBUG)iprint=1  c         if(DEBUG)iprint=1
364           if(VERBOSE)iprint=1           if(VERBOSE)iprint=1
365             if(DEBUG)iprint=2
366           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
367           if(ifail.ne.0) then           if(ifail.ne.0) then
368              if(.true.)then              if(VERBOSE)then
369                 print *,                 print *,
370       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
371       $              ,iev       $              ,iev
# Line 446  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 484  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 537  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          if(sensor.lt.1.or.sensor.gt.2)then
600             print*,'xyz_PAM   ***ERROR*** wrong input '
601             print*,'sensor ',sensor
602             icx=0
603             icy=0
604          endif
605    
606  *     -----------------  *     -----------------
607  *     CLUSTER X  *     CLUSTER X
608  *     -----------------  *     -----------------      
   
609        if(icx.ne.0)then        if(icx.ne.0)then
610           viewx = VIEW(icx)  
611           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
612           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
613           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
614                     resxPAM = RESXAV
615           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
616           if(PFAx.eq.'COG1')then  !(1)  
617              stripx = stripx      !(1)           if(
618              resxPAM = resxPAM    !(1)       $        viewx.lt.1.or.        
619         $        viewx.gt.12.or.
620         $        nldx.lt.1.or.
621         $        nldx.gt.3.or.
622         $        stripx.lt.1.or.
623         $        stripx.gt.3072.or.
624         $        .false.)then
625                print*,'xyz_PAM   ***ERROR*** wrong input '
626                print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
627                icx = 0
628                goto 10
629             endif
630    
631    *        --------------------------
632    *        magnetic-field corrections
633    *        --------------------------
634             angtemp  = ax
635             bfytemp  = bfy
636    *        /////////////////////////////////
637    *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!
638    *        *grvzkkjsdgjhhhgngbn###>:(
639    *        /////////////////////////////////
640    c         if(nplx.eq.6) angtemp = -1. * ax
641    c         if(nplx.eq.6) bfytemp = -1. * bfy
642             if(viewx.eq.12) angtemp = -1. * ax
643             if(viewx.eq.12) bfytemp = -1. * bfy
644             tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
645             angx     = 180.*atan(tgtemp)/acos(-1.)
646             stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
647    c$$$         print*,nplx,ax,bfy/10.
648    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
649    c$$$         print*,'========================'
650    c$$$         if(bfy.ne.0.)print*,viewx,'-x- '
651    c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ
652    *        --------------------------
653    
654    c$$$         print*,'--- x-cl ---'
655    c$$$         istart = INDSTART(ICX)
656    c$$$         istop  = TOTCLLENGTH
657    c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
658    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
659    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
660    c$$$         print*,INDMAX(icx)
661    c$$$         print*,cog(4,icx)
662    c$$$         print*,fbad_cog(4,icx)
663            
664    
665             if(PFAx.eq.'COG1')then
666    
667                stripx  = stripx
668                resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM
669    
670           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
671              stripx = stripx + cog(2,icx)              
672                stripx  = stripx + cog(2,icx)            
673                resxPAM = risx_cog(abs(angx))!TEMPORANEO              
674              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
675    
676             elseif(PFAx.eq.'COG3')then
677    
678                stripx  = stripx + cog(3,icx)            
679                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
680                resxPAM = resxPAM*fbad_cog(3,icx)
681    
682             elseif(PFAx.eq.'COG4')then
683    
684                stripx  = stripx + cog(4,icx)            
685                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
686                resxPAM = resxPAM*fbad_cog(4,icx)
687    
688           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
689  c            cog2 = cog(2,icx)  
690  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                          stripx  = stripx + pfaeta2(icx,angx)          
691  c            stripx = stripx + etacorr              resxPAM = risxeta2(abs(angx))
692              stripx = stripx + pfaeta2(icx,angx)            !(3)              resxPAM = resxPAM*fbad_cog(2,icx)
             resxPAM = risx_eta2(angx)                       !   (4)  
693              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
694       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
695              resxPAM = resxPAM*fbad_cog(2,icx)  
696           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                        
697              stripx = stripx + pfaeta3(icx,angx)            !(3)  
698              resxPAM = risx_eta3(angx)                       !   (4)              stripx  = stripx + pfaeta3(icx,angx)          
699              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              resxPAM = risxeta3(abs(angx))                      
700       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
701              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)  c            if(DEBUG.and.fbad_cog(3,icx).ne.1)            
702           elseif(PFAx.eq.'ETA4')then                         !(3)  c     $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
703              stripx = stripx + pfaeta4(icx,angx)            !(3)  
704              resxPAM = risx_eta4(angx)                       !   (4)           elseif(PFAx.eq.'ETA4')then                        
705              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
706       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)              stripx  = stripx + pfaeta4(icx,angx)            
707              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = risxeta4(abs(angx))                      
708           elseif(PFAx.eq.'ETA')then                          !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
709              stripx = stripx + pfaeta(icx,angx)             !(3)  c            if(DEBUG.and.fbad_cog(4,icx).ne.1)              
710              resxPAM = ris_eta(icx,angx)                     !   (4)  c     $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
711              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
712       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)           elseif(PFAx.eq.'ETA')then  
713  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
714              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)              stripx  = stripx + pfaeta(icx,angx)            
715           elseif(PFAx.eq.'COG')then           !(2)  c            resxPAM = riseta(icx,angx)                    
716              stripx = stripx + cog(0,icx)     !(2)                      resxPAM = riseta(viewx,angx)                    
717              resxPAM = risx_cog(angx)                        !   (4)              resxPAM = resxPAM*fbad_eta(icx,angx)            
718              resxPAM = resxPAM*fbad_cog(0,icx)!(2)  c            if(DEBUG.and.fbad_cog(2,icx).ne.1)              
719    c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
720    
721             elseif(PFAx.eq.'ETAL')then  
722    
723                stripx  = stripx + pfaetal(icx,angx)            
724                resxPAM = riseta(viewx,angx)                    
725                resxPAM = resxPAM*fbad_eta(icx,angx)            
726    c            if(DEBUG.and.fbad_cog(2,icx).ne.1)              
727    c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
728    
729             elseif(PFAx.eq.'COG')then          
730    
731                stripx  = stripx + cog(0,icx)            
732                resxPAM = risx_cog(abs(angx))                    
733                resxPAM = resxPAM*fbad_cog(0,icx)
734    
735           else           else
736              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
737           endif           endif
738    
739        endif  
740  c      if(icy.eq.0.and.icx.ne.0)  *     ======================================
741  c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  *     temporary patch for saturated clusters
742          *     ======================================
743             if( nsatstrips(icx).gt.0 )then
744                stripx  = stripx + cog(4,icx)            
745                resxPAM = pitchX*1e-4/sqrt(12.)
746                cc=cog(4,icx)
747    c$$$            print*,icx,' *** ',cc
748    c$$$            print*,icx,' *** ',resxPAM
749             endif
750    
751     10   endif
752    
753        
754  *     -----------------  *     -----------------
755  *     CLUSTER Y  *     CLUSTER Y
756  *     -----------------  *     -----------------
757    
758        if(icy.ne.0)then        if(icy.ne.0)then
759    
760           viewy = VIEW(icy)           viewy = VIEW(icy)
761           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
762           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
763           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
764             stripy = float(MAXS(icy))
765    
766             if(
767         $        viewy.lt.1.or.        
768         $        viewy.gt.12.or.
769         $        nldy.lt.1.or.
770         $        nldy.gt.3.or.
771         $        stripy.lt.1.or.
772         $        stripy.gt.3072.or.
773         $        .false.)then
774                print*,'xyz_PAM   ***ERROR*** wrong input '
775                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
776                icy = 0
777                goto 20
778             endif
779    
780           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
781              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG) then
782       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
783         $              ,icx,icy
784                endif
785              goto 100              goto 100
786           endif           endif
787            *        --------------------------
788           stripy = float(MAXS(icy))  *        magnetic-field corrections
789           if(PFAy.eq.'COG1')then !(1)  *        --------------------------
790              stripy = stripy     !(1)           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
791              resyPAM = resyPAM   !(1)           angy    = 180.*atan(tgtemp)/acos(-1.)
792             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
793    c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
794    c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
795    *        --------------------------
796            
797    c$$$         print*,'--- y-cl ---'
798    c$$$         istart = INDSTART(ICY)
799    c$$$         istop  = TOTCLLENGTH
800    c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
801    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
802    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
803    c$$$         print*,INDMAX(icy)
804    c$$$         print*,cog(4,icy)
805    c$$$         print*,fbad_cog(4,icy)
806    
807             if(PFAy.eq.'COG1')then
808    
809                stripy  = stripy    
810                resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  
811    
812           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
813              stripy = stripy + cog(2,icy)  
814                stripy  = stripy + cog(2,icy)
815                resyPAM = risy_cog(abs(angy))!TEMPORANEO
816              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
817    
818             elseif(PFAy.eq.'COG3')then
819    
820                stripy  = stripy + cog(3,icy)
821                resyPAM = risy_cog(abs(angy))!TEMPORANEO
822                resyPAM = resyPAM*fbad_cog(3,icy)
823    
824             elseif(PFAy.eq.'COG4')then
825    
826                stripy  = stripy + cog(4,icy)
827                resyPAM = risy_cog(abs(angy))!TEMPORANEO
828                resyPAM = resyPAM*fbad_cog(4,icy)
829    
830           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
831  c            cog2 = cog(2,icy)  
832  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)              stripy  = stripy + pfaeta2(icy,angy)          
833  c            stripy = stripy + etacorr              resyPAM = risyeta2(abs(angy))              
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
834              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
835              if(DEBUG.and.fbad_cog(2,icy).ne.1)  c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
836       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
837           elseif(PFAy.eq.'ETA3')then                         !(3)  
838              stripy = stripy + pfaeta3(icy,angy)            !(3)           elseif(PFAy.eq.'ETA3')then                      
839              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
840              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              stripy  = stripy + pfaeta3(icy,angy)
841       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)              resyPAM = resyPAM*fbad_cog(3,icy)  
842           elseif(PFAy.eq.'ETA4')then                         !(3)  c            if(DEBUG.and.fbad_cog(3,icy).ne.1)
843              stripy = stripy + pfaeta4(icy,angy)            !(3)  c     $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
844              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
845              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA4')then  
846       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
847           elseif(PFAy.eq.'ETA')then                          !(3)              stripy  = stripy + pfaeta4(icy,angy)
848              stripy = stripy + pfaeta(icy,angy)             !(3)              resyPAM = resyPAM*fbad_cog(4,icy)
849              resyPAM = ris_eta(icy,angy)                     !   (4)  c            if(DEBUG.and.fbad_cog(4,icy).ne.1)
850  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c     $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
851              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
852              if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA')then
853       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
854                stripy  = stripy + pfaeta(icy,angy)
855    c            resyPAM = riseta(icy,angy)  
856                resyPAM = riseta(viewy,angy)  
857                resyPAM = resyPAM*fbad_eta(icy,angy)
858    c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
859    c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
860    
861             elseif(PFAy.eq.'ETAL')then
862    
863                stripy  = stripy + pfaetal(icy,angy)
864                resyPAM = riseta(viewy,angy)  
865                resyPAM = resyPAM*fbad_eta(icy,angy)
866    c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
867    c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
868    
869           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
870              stripy = stripy + cog(0,icy)              
871              resyPAM = risy_cog(angy)                        !   (4)              stripy  = stripy + cog(0,icy)            
872  c            resyPAM = ris_eta(icy,angy)                    !   (4)              resyPAM = risy_cog(abs(angy))
873              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
874    
875           else           else
876              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
877           endif           endif
878    
       endif  
879    
880          *     ======================================
881    *     temporary patch for saturated clusters
882    *     ======================================
883             if( nsatstrips(icy).gt.0 )then
884                stripy  = stripy + cog(4,icy)            
885                resyPAM = pitchY*1e-4/sqrt(12.)
886                cc=cog(4,icy)
887    c$$$            print*,icy,' *** ',cc
888    c$$$            print*,icy,' *** ',resyPAM
889             endif
890    
891    
892     20   endif
893    
894    c$$$      print*,'## stripx,stripy ',stripx,stripy
895    
896  c===========================================================  c===========================================================
897  C     COUPLE  C     COUPLE
898  C===========================================================  C===========================================================
# Line 667  c     (xi,yi,zi) = mechanical coordinate Line 903  c     (xi,yi,zi) = mechanical coordinate
903  c------------------------------------------------------------------------  c------------------------------------------------------------------------
904           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
905       $        .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...
906              print*,'xyz_PAM (couple):',              if(DEBUG) then
907       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
908         $              ' WARNING: false X strip: strip ',stripx
909                endif
910           endif           endif
911           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
912           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
# Line 760  c            print*,'X-singlet ',icx,npl Line 998  c            print*,'X-singlet ',icx,npl
998  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...
999              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1000       $           .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...
1001                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG) then
1002       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
1003         $                 ' WARNING: false X strip: strip ',stripx
1004                   endif
1005              endif              endif
1006              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
1007    
# Line 783  c            print*,'X-cl ',icx,stripx,' Line 1023  c            print*,'X-cl ',icx,stripx,'
1023  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
1024    
1025           else           else
1026                if(DEBUG) then
1027              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
1028              print *,'icx = ',icx                 print *,'icx = ',icx
1029              print *,'icy = ',icy                 print *,'icy = ',icy
1030                endif
1031              goto 100              goto 100
1032                            
1033           endif           endif
# Line 851  c--------------------------------------- Line 1092  c---------------------------------------
1092  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1093    
1094        else        else
1095                       if(DEBUG) then
1096           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1097           print *,'icx = ',icx              print *,'icx = ',icx
1098           print *,'icy = ',icy              print *,'icy = ',icy
1099                         endif
1100        endif        endif
1101                    
1102    
1103    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1104    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1105    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1106    
1107   100  continue   100  continue
1108        end        end
1109    
1110    ************************************************************************
1111    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1112    *     (done to be called from c/c++)
1113    ************************************************************************
1114    
1115          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1116    
1117          include 'commontracker.f'
1118          include 'level1.f'
1119          include 'common_mini_2.f'
1120          include 'common_xyzPAM.f'
1121          include 'common_mech.f'
1122          include 'calib.f'
1123          
1124    *     flag to chose PFA
1125    c$$$      character*10 PFA
1126    c$$$      common/FINALPFA/PFA
1127    
1128          integer icx,icy           !X-Y cluster ID
1129          integer sensor
1130          character*4 PFAx,PFAy     !PFA to be used
1131          real ax,ay                !X-Y geometric angle
1132          real bfx,bfy              !X-Y b-field components
1133    
1134          ipx=0
1135          ipy=0      
1136          
1137    c$$$      PFAx = 'COG4'!PFA
1138    c$$$      PFAy = 'COG4'!PFA
1139    
1140          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1141                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1142         $           ,' does not exists (nclstr1=',nclstr1,')'
1143                icx = -1*icx
1144                icy = -1*icy
1145                return
1146            
1147          endif
1148          
1149          call idtoc(pfaid,PFAx)
1150          call idtoc(pfaid,PFAy)
1151    
1152    c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1153    
1154    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1155          
1156          if(icx.ne.0.and.icy.ne.0)then
1157    
1158             ipx=npl(VIEW(icx))
1159             ipy=npl(VIEW(icy))
1160    c$$$         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1161    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1162    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1163            
1164             if( (nplanes-ipx+1).ne.ip )then
1165                print*,'xyzpam: ***WARNING*** cluster ',icx
1166         $           ,' does not belong to plane: ',ip
1167                icx = -1*icx
1168                return
1169             endif
1170             if( (nplanes-ipy+1).ne.ip )then
1171                print*,'xyzpam: ***WARNING*** cluster ',icy
1172         $           ,' does not belong to plane: ',ip
1173                icy = -1*icy
1174                return
1175             endif
1176    
1177             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1178    
1179             xgood(ip) = 1.
1180             ygood(ip) = 1.
1181             resx(ip)  = resxPAM
1182             resy(ip)  = resyPAM
1183    
1184             xm(ip) = xPAM
1185             ym(ip) = yPAM
1186             zm(ip) = zPAM
1187             xm_A(ip) = 0.
1188             ym_A(ip) = 0.
1189             xm_B(ip) = 0.
1190             ym_B(ip) = 0.
1191    
1192    c         zv(ip) = zPAM
1193    
1194          elseif(icx.eq.0.and.icy.ne.0)then
1195    
1196             ipy=npl(VIEW(icy))
1197    c$$$         if((nplanes-ipy+1).ne.ip)
1198    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1199    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1200             if( (nplanes-ipy+1).ne.ip )then
1201                print*,'xyzpam: ***WARNING*** cluster ',icy
1202         $           ,' does not belong to plane: ',ip
1203                icy = -1*icy
1204                return
1205             endif
1206    
1207             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1208            
1209             xgood(ip) = 0.
1210             ygood(ip) = 1.
1211             resx(ip)  = 1000.
1212             resy(ip)  = resyPAM
1213    
1214             xm(ip) = -100.
1215             ym(ip) = -100.
1216             zm(ip) = (zPAM_A+zPAM_B)/2.
1217             xm_A(ip) = xPAM_A
1218             ym_A(ip) = yPAM_A
1219             xm_B(ip) = xPAM_B
1220             ym_B(ip) = yPAM_B
1221    
1222    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1223            
1224          elseif(icx.ne.0.and.icy.eq.0)then
1225    
1226             ipx=npl(VIEW(icx))
1227    c$$$         if((nplanes-ipx+1).ne.ip)
1228    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1229    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1230    
1231             if( (nplanes-ipx+1).ne.ip )then
1232                print*,'xyzpam: ***WARNING*** cluster ',icx
1233         $           ,' does not belong to plane: ',ip
1234                icx = -1*icx
1235                return
1236             endif
1237    
1238             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1239          
1240             xgood(ip) = 1.
1241             ygood(ip) = 0.
1242             resx(ip)  = resxPAM
1243             resy(ip)  = 1000.
1244    
1245             xm(ip) = -100.
1246             ym(ip) = -100.
1247             zm(ip) = (zPAM_A+zPAM_B)/2.
1248             xm_A(ip) = xPAM_A
1249             ym_A(ip) = yPAM_A
1250             xm_B(ip) = xPAM_B
1251             ym_B(ip) = yPAM_B
1252            
1253    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1254    
1255          else
1256    
1257             il = 2
1258             if(lad.ne.0)il=lad
1259             is = 1
1260             if(sensor.ne.0)is=sensor
1261    c         print*,nplanes-ip+1,il,is
1262    
1263             xgood(ip) = 0.
1264             ygood(ip) = 0.
1265             resx(ip)  = 1000.
1266             resy(ip)  = 1000.
1267    
1268             xm(ip) = -100.
1269             ym(ip) = -100.          
1270             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1271             xm_A(ip) = 0.
1272             ym_A(ip) = 0.
1273             xm_B(ip) = 0.
1274             ym_B(ip) = 0.
1275    
1276    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1277    
1278          endif
1279    
1280          if(DEBUG)then
1281    c         print*,'----------------------------- track coord'
1282    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1283             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1284         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1285         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1286    c$$$         print*,'-----------------------------'
1287          endif
1288          end
1289    
1290  ********************************************************************************  ********************************************************************************
1291  ********************************************************************************  ********************************************************************************
# Line 936  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1361  c         print*,'A-(',xPAM_A,yPAM_A,')
1361           endif                   endif        
1362    
1363           distance=           distance=
1364       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1365    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1366           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1367    
1368  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 961  c$$$         print*,' resolution ',re Line 1387  c$$$         print*,' resolution ',re
1387  *     ----------------------  *     ----------------------
1388                    
1389           distance=           distance=
1390       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1391       $        +       $       +
1392       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1393    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1394    c$$$     $        +
1395    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1396           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1397    
1398  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 972  c$$$         print*,' resolution ',resxP Line 1401  c$$$         print*,' resolution ',resxP
1401                    
1402        else        else
1403                    
1404           print*  c         print*
1405       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1406           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1407           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 '
1408       $        ,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
1409        endif          endif  
1410    
1411        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1044  c--------------------------------------- Line 1473  c---------------------------------------
1473                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1474       $              .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...
1475  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...
1476                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1477       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1478                 endif                 endif
1479                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1480                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1200  c      include 'common_analysis.f' Line 1629  c      include 'common_analysis.f'
1629        is_cp=0        is_cp=0
1630        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1631        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1632        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1633    
1634        return        return
1635        end        end
# Line 1329  c      include 'level1.f' Line 1758  c      include 'level1.f'
1758  *     mask views with too many clusters  *     mask views with too many clusters
1759        do iv=1,nviews        do iv=1,nviews
1760           if( ncl_view(iv).gt. nclusterlimit)then           if( ncl_view(iv).gt. nclusterlimit)then
1761              mask_view(iv) = 1  c            mask_view(iv) = 1
1762                mask_view(iv) = mask_view(iv) + 2**0
1763              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1764       $           ,nclusterlimit,' on view ', iv,' --> masked!'       $           ,nclusterlimit,' on view ', iv,' --> masked!'
1765           endif           endif
# Line 1348  c      include 'level1.f' Line 1778  c      include 'level1.f'
1778  *     ----------------------------------------------------  *     ----------------------------------------------------
1779  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1780  *     ----------------------------------------------------  *     ----------------------------------------------------
1781           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1782              cl_single(icx)=0              cl_single(icx)=0
1783              goto 10              goto 10
1784           endif           endif
# Line 1398  c     endif Line 1828  c     endif
1828  *     ----------------------------------------------------  *     ----------------------------------------------------
1829  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1830  *     ----------------------------------------------------  *     ----------------------------------------------------
1831              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1832                 cl_single(icy)=0                 cl_single(icy)=0
1833                 goto 20                 goto 20
1834              endif              endif
# Line 1444  c     endif Line 1874  c     endif
1874  *     charge correlation  *     charge correlation
1875  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1876    
1877                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1878       $              .and.       $              .and.
1879       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1880       $              .and.       $              .and.
1881       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1882       $              .and.       $              .and.
1883       $              .true.)then       $              .true.)then
1884    
1885                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1886       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1887                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1888    
1889  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1890    
1891                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1892       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1893                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1894                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1467  c                  cut = chcut * sch(npl Line 1897  c                  cut = chcut * sch(npl
1897                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1898                    endif                    endif
1899                 endif                 endif
                 
 *     ------------------> COUPLE <------------------  
                ncp_plane(nplx) = ncp_plane(nplx) + 1  
                clx(nplx,ncp_plane(nplx))=icx  
                cly(nply,ncp_plane(nplx))=icy  
                cl_single(icx)=0  
                cl_single(icy)=0  
1900    
1901                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1902                    if(verbose)print*,                    if(verbose)print*,
1903       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1904       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1905       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1906       $                 ,'( ',ncouplemax,' ) --> masked!'       $                 ,'( ',ncouplemax,' ) --> masked!'
1907                    mask_view(nviewx(nplx)) = 2  c                  mask_view(nviewx(nplx)) = 2
1908                    mask_view(nviewy(nply)) = 2  c                  mask_view(nviewy(nply)) = 2
1909                      mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1910                      mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1911                      goto 10
1912                 endif                 endif
1913                  
1914    *     ------------------> COUPLE <------------------
1915                   ncp_plane(nplx) = ncp_plane(nplx) + 1
1916                   clx(nplx,ncp_plane(nplx))=icx
1917                   cly(nply,ncp_plane(nplx))=icy
1918                   cl_single(icx)=0
1919                   cl_single(icy)=0
1920  *     ----------------------------------------------  *     ----------------------------------------------
1921    
1922              endif                                            endif                              
# Line 1575  c      double precision xm3,ym3,zm3 Line 2008  c      double precision xm3,ym3,zm3
2008  *     --------------------------------------------  *     --------------------------------------------
2009        do ip=1,nplanes        do ip=1,nplanes
2010           if(ncp_plane(ip).gt.ncouplelimit)then           if(ncp_plane(ip).gt.ncouplelimit)then
2011              mask_view(nviewx(ip)) = 8  c            mask_view(nviewx(ip)) = 8
2012              mask_view(nviewy(ip)) = 8  c            mask_view(nviewy(ip)) = 8
2013                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
2014                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
2015           endif           endif
2016        enddo        enddo
2017    
# Line 1592  c      double precision xm3,ym3,zm3 Line 2027  c      double precision xm3,ym3,zm3
2027                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
2028                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
2029  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
2030                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
2031                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
2032                 xm1=xPAM                 xm1=xPAM
2033                 ym1=yPAM                 ym1=yPAM
2034                 zm1=zPAM                                   zm1=zPAM                  
# Line 1608  c     print*,'***',is1,xm1,ym1,zm1 Line 2044  c     print*,'***',is1,xm1,ym1,zm1
2044                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
2045  c                        call xyz_PAM  c                        call xyz_PAM
2046  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
2047    c                        call xyz_PAM
2048    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2049                          call xyz_PAM                          call xyz_PAM
2050       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2051                          xm2=xPAM                          xm2=xPAM
2052                          ym2=yPAM                          ym2=yPAM
2053                          zm2=zPAM                          zm2=zPAM
# Line 1626  c     $                       (icx2,icy2 Line 2064  c     $                       (icx2,icy2
2064  c                           good2=.false.  c                           good2=.false.
2065  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
2066                             do iv=1,12                             do iv=1,12
2067                                mask_view(iv) = 3  c                              mask_view(iv) = 3
2068                                  mask_view(iv) = mask_view(iv)+ 2**2
2069                             enddo                             enddo
2070                             iflag=1                             iflag=1
2071                             return                             return
# Line 1673  c$$$ Line 2112  c$$$
2112                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2113  c                                 call xyz_PAM  c                                 call xyz_PAM
2114  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2115    c                                 call xyz_PAM
2116    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2117                                   call xyz_PAM                                   call xyz_PAM
2118       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2119         $                                ,0.,0.,0.,0.)
2120                                   xm3=xPAM                                   xm3=xPAM
2121                                   ym3=yPAM                                   ym3=yPAM
2122                                   zm3=zPAM                                   zm3=zPAM
# Line 1702  c     $                                 Line 2144  c     $                                
2144  c                                    good2=.false.  c                                    good2=.false.
2145  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2146                                      do iv=1,nviews                                      do iv=1,nviews
2147                                         mask_view(iv) = 4  c                                       mask_view(iv) = 4
2148                                           mask_view(iv)=mask_view(iv)+ 2**3
2149                                      enddo                                      enddo
2150                                      iflag=1                                      iflag=1
2151                                      return                                      return
# Line 1936  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2379  c         if(ncpused.lt.ncpyz_min)goto 2
2379  c               good2=.false.  c               good2=.false.
2380  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2381              do iv=1,nviews              do iv=1,nviews
2382                 mask_view(iv) = 5  c               mask_view(iv) = 5
2383                   mask_view(iv) = mask_view(iv) + 2**4
2384              enddo              enddo
2385              iflag=1              iflag=1
2386              return              return
# Line 2146  c     print*,'check cp_used' Line 2590  c     print*,'check cp_used'
2590           enddo           enddo
2591  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2592           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2593           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2594                    
2595  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2596  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
# Line 2158  c         if(ncpused.lt.ncpxz_min)goto 2 Line 2602  c         if(ncpused.lt.ncpxz_min)goto 2
2602  c     good2=.false.  c     good2=.false.
2603  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2604              do iv=1,nviews              do iv=1,nviews
2605                 mask_view(iv) = 6  c               mask_view(iv) = 6
2606                   mask_view(iv) =  mask_view(iv) + 2**5
2607              enddo              enddo
2608              iflag=1              iflag=1
2609              return              return
# Line 2223  c$$$     $           ,(tr_cloud(iii),iii Line 2668  c$$$     $           ,(tr_cloud(iii),iii
2668  **************************************************  **************************************************
2669    
2670        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2671    
2672        include 'commontracker.f'        include 'commontracker.f'
2673        include 'level1.f'        include 'level1.f'
# Line 2233  c*************************************** Line 2675  c***************************************
2675        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2676        include 'common_mini_2.f'        include 'common_mini_2.f'
2677        include 'common_mech.f'        include 'common_mech.f'
2678  c      include 'momanhough_init.f'  
2679    
2680    
2681  *     output flag  *     output flag
# Line 2323  c      real fitz(nplanes)        !z coor Line 2765  c      real fitz(nplanes)        !z coor
2765                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2766              enddo              enddo
2767                            
2768              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2769                if(nplused.lt.nplyz_min)goto 888 !next doublet
2770              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2771                            
2772              if(DEBUG)then              if(DEBUG)then
# Line 2427  c$$$            if(AL_INI(5).gt.defmax)g Line 2870  c$$$            if(AL_INI(5).gt.defmax)g
2870  *                                   *************************  *                                   *************************
2871  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2872  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2873    c                                    call xyz_PAM(icx,icy,is, !(1)
2874    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2875                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2876       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2877  *                                   *************************  *                                   *************************
2878  *                                   -----------------------------  *                                   -----------------------------
2879                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 2456  c$$$                              enddo Line 2901  c$$$                              enddo
2901                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2902                                iprint=0                                iprint=0
2903  c                              if(DEBUG)iprint=1  c                              if(DEBUG)iprint=1
2904                                if(DEBUG)iprint=1                                if(DEBUG)iprint=2
2905                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2906                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2907                                   if(DEBUG)then                                   if(DEBUG)then
# Line 2491  c                                 chi2=- Line 2936  c                                 chi2=-
2936  c                                 good2=.false.  c                                 good2=.false.
2937  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2938                                   do iv=1,nviews                                   do iv=1,nviews
2939                                      mask_view(iv) = 7  c                                    mask_view(iv) = 7
2940                                        mask_view(iv) = mask_view(iv) + 2**6
2941                                   enddo                                   enddo
2942                                   iflag=1                                   iflag=1
2943                                   return                                   return
# Line 2499  c                                 goto 8 Line 2945  c                                 goto 8
2945                                                                
2946                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2947                                                                
 c$$$                              ndof=0                                  
2948                                do ip=1,nplanes                                do ip=1,nplanes
2949  c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2950                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2951                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2952                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 2522  c$$$     $                               Line 2965  c$$$     $                              
2965                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2966                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2967       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2968                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2969         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2970                                        LADDER_STORE(nplanes-ip+1,ntracks)
2971         $                                   = LADDER(
2972         $                                   clx(ip,icp_cp(
2973         $                                   cp_match(ip,hit_plane(ip)
2974         $                                   ))));
2975                                   else                                   else
2976                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2977                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2978                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2979                                   endif                                   endif
2980                                     BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2981                                     BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2982                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(nplanes-ip+1,ntracks)=0
2983                                   do i=1,5                                   do i=1,5
2984                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2985                                   enddo                                   enddo
2986                                enddo                                enddo
2987                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2988                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2989                                                                
2990  *     --------------------------------  *     --------------------------------
# Line 2558  c$$$  rchi2=chi2/dble(ndof) Line 3008  c$$$  rchi2=chi2/dble(ndof)
3008           return           return
3009        endif        endif
3010                
3011    c$$$      if(DEBUG)then
3012    c$$$         print*,'****** TRACK CANDIDATES ***********'
3013    c$$$         print*,'#         R. chi2        RIG'
3014    c$$$         do i=1,ntracks
3015    c$$$            print*,i,' --- ',rchi2_store(i),' --- '
3016    c$$$     $           ,1./abs(AL_STORE(5,i))
3017    c$$$         enddo
3018    c$$$         print*,'***********************************'
3019    c$$$      endif
3020        if(DEBUG)then        if(DEBUG)then
3021           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
3022           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
3023           do i=1,ntracks          do i=1,ntracks
3024              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
3025       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
3026           enddo              ndof=ndof           !(1)
3027           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
3028         $           +int(ygood_store(ii,i)) !(1)
3029              enddo                 !(1)
3030              print*,i,' --- ',rchi2_store(i),' --- '
3031         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3032            enddo
3033            print*,'*****************************************'
3034        endif        endif
3035                
3036                
# Line 2584  c$$$  rchi2=chi2/dble(ndof) Line 3049  c$$$  rchi2=chi2/dble(ndof)
3049    
3050        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3051    
 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******************************************************  
3052    
3053        include 'commontracker.f'        include 'commontracker.f'
3054        include 'level1.f'        include 'level1.f'
# Line 2596  c*************************************** Line 3056  c***************************************
3056        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3057        include 'common_mini_2.f'        include 'common_mini_2.f'
3058        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
3059        include 'calib.f'        include 'calib.f'
3060    
   
3061  *     flag to chose PFA  *     flag to chose PFA
3062        character*10 PFA        character*10 PFA
3063        common/FINALPFA/PFA        common/FINALPFA/PFA
3064    
3065          real k(6)
3066          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3067    
3068          real xp,yp,zp
3069          real xyzp(3),bxyz(3)
3070          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3071    
3072  *     =================================================  *     =================================================
3073  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3074  *                          and  *                          and
# Line 2613  c      include 'level1.f' Line 3077  c      include 'level1.f'
3077        call track_init        call track_init
3078        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3079    
3080             xP=XV_STORE(nplanes-ip+1,ibest)
3081             yP=YV_STORE(nplanes-ip+1,ibest)
3082             zP=ZV_STORE(nplanes-ip+1,ibest)
3083             call gufld(xyzp,bxyz)
3084             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3085             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3086    c$$$  bxyz(1)=0
3087    c$$$         bxyz(2)=0
3088    c$$$         bxyz(3)=0
3089  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3090  *     -------------------------------------------------  *     -------------------------------------------------
3091  *     If the plane has been already included, it just  *     If the plane has been already included, it just
# Line 2633  c      include 'level1.f' Line 3106  c      include 'level1.f'
3106       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3107              icx=clx(ip,icp)              icx=clx(ip,icp)
3108              icy=cly(ip,icp)              icy=cly(ip,icp)
3109    c            call xyz_PAM(icx,icy,is,
3110    c     $           PFA,PFA,
3111    c     $           AXV_STORE(nplanes-ip+1,ibest),
3112    c     $           AYV_STORE(nplanes-ip+1,ibest))
3113              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3114       $           PFA,PFA,       $           PFA,PFA,
3115       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3116       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3117  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3118  c$$$  $              'COG2','COG2',       $           bxyz(2)
3119  c$$$  $              0.,       $           )
3120  c$$$  $              0.)  
3121              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3122              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3123              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 2650  c$$$  $              0.) Line 3126  c$$$  $              0.)
3126              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3127              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3128    
3129  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3130              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)  
3131                            
3132  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3133  *     -------------------------------------------------  *     -------------------------------------------------
# Line 2667  c            dedxtrk(nplanes-ip+1) = (de Line 3142  c            dedxtrk(nplanes-ip+1) = (de
3142                                
3143  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3144  *     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)  
3145              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3146  *     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
3147              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3148    
3149                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3150                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3151  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3152    
3153              if(DEBUG)then              if(DEBUG)then
# Line 2709  c     $              cl_used(icy).eq.1.o Line 3184  c     $              cl_used(icy).eq.1.o
3184  *            *          
3185                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3186       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3187       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3188       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3189         $              bxyz(1),
3190         $              bxyz(2)
3191         $              )
3192                                
3193                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3194                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3195                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3196                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
3197       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3198                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3199                    xmm = xPAM                    xmm = xPAM
3200                    ymm = yPAM                    ymm = yPAM
# Line 2726  c     $              'ETA2','ETA2', Line 3203  c     $              'ETA2','ETA2',
3203                    rymm = resyPAM                    rymm = resyPAM
3204                    distmin = distance                    distmin = distance
3205                    idm = id                                      idm = id                  
3206  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3207                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3208                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3209                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3210         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3211                 endif                 endif
3212   1188          continue   1188          continue
3213              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3214              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3215                if(distmin.le.clincnewc)then     !QUIQUI              
3216  *              -----------------------------------  *              -----------------------------------
3217                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3218                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3219                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3220                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3221                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3222                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3223                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3224  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3225                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3226  *              -----------------------------------  *              -----------------------------------
3227                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3228                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG)print*,'%%%% included couple ',idm
3229       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3230                 goto 133         !next plane                 goto 133         !next plane
3231              endif              endif
3232  *     ================================================  *     ================================================
# Line 2780  c            if(DEBUG)print*,'>>>> try t Line 3259  c            if(DEBUG)print*,'>>>> try t
3259  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
3260                 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)
3261  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3262    c               call xyz_PAM(icx,0,ist,
3263    c     $              PFA,PFA,
3264    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3265                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3266       $              PFA,PFA,       $              PFA,PFA,
3267       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3268         $              bxyz(1),
3269         $              bxyz(2)
3270         $              )              
3271                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3272                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3273                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3274       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3275                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3276                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3277                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2799  c     $              'ETA2','ETA2', Line 3283  c     $              'ETA2','ETA2',
3283                    rymm = resyPAM                    rymm = resyPAM
3284                    distmin = distance                    distmin = distance
3285                    iclm = icx                    iclm = icx
3286  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3287                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3288                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3289                 endif                                   endif                  
3290  11881          continue  11881          continue
# Line 2808  c                  dedxmm = dedx(icx) !( Line 3292  c                  dedxmm = dedx(icx) !(
3292  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
3293                 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)
3294  *                                              !jump to the next couple  *                                              !jump to the next couple
3295    c               call xyz_PAM(0,icy,ist,
3296    c     $              PFA,PFA,
3297    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3298                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3299       $              PFA,PFA,       $              PFA,PFA,
3300       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3301         $              bxyz(1),
3302         $              bxyz(2)
3303         $              )
3304                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3305                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3306                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3307       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3308                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3309                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3310                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2827  c     $              'ETA2','ETA2', Line 3316  c     $              'ETA2','ETA2',
3316                    rymm = resyPAM                    rymm = resyPAM
3317                    distmin = distance                    distmin = distance
3318                    iclm = icy                    iclm = icy
3319  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3320                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3321                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3322                 endif                                   endif                  
3323  11882          continue  11882          continue
3324              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3325  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3326    c            print*,'## ncls(',ip,') ',ncls(ip)
3327              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3328                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3329  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 2842  c               if(cl_used(icl).eq.1.or. Line 3332  c               if(cl_used(icl).eq.1.or.
3332       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3333                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3334                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3335       $                 PFA,PFA,       $                 PFA,PFA,
3336       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3337         $                 bxyz(1),
3338         $                 bxyz(2)
3339         $                 )
3340                 else                         !<---- Y view                 else                         !<---- Y view
3341                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3342       $                 PFA,PFA,       $                 PFA,PFA,
3343       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3344         $                 bxyz(1),
3345         $                 bxyz(2)
3346         $                 )
3347                 endif                 endif
3348    
3349                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3350                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3351                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3352       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance,'<',distmin,' ?'
3353                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3354                      if(DEBUG)print*,'YES'
3355                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3356                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3357                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 2867  c     $                 'ETA2','ETA2', Line 3362  c     $                 'ETA2','ETA2',
3362                    rymm = resyPAM                    rymm = resyPAM
3363                    distmin = distance                      distmin = distance  
3364                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3365                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3366                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3367                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3368                    else          !<---- Y view                    else          !<---- Y view
3369                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3370                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3371                    endif                    endif
3372                 endif                                   endif                  
3373  18882          continue  18882          continue
3374              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3375    c            print*,'## distmin ', distmin,' clinc ',clinc
3376    
3377              if(distmin.le.clinc)then                    c     QUIQUI------------
3378                  c     anche qui: non ci vuole clinc???
3379                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3380                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3381                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3382                    resx(nplanes-ip+1)=rxmm       $                 20*
3383                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3384  c                  if(.true.)print*,'%%%% included X-cl ',iclm                 else if(mod(VIEW(iclm),2).ne.0)then
3385       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    clincnew=
3386       $                 ,', norm.dist.= ',distmin       $                 10*
3387       $                 ,', cut ',clinc,' )'       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3388                 else                 endif
3389                    YGOOD(nplanes-ip+1)=1.  c     QUIQUI------------
3390                    resy(nplanes-ip+1)=rymm                
3391                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 if(distmin.le.clincnew)then   !QUIQUI
3392  c                  if(.true.)print*,'%%%% included Y-cl ',iclm  c     if(distmin.le.clinc)then          !QUIQUI          
3393       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    
3394       $                 ,', norm.dist.= ', distmin                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3395       $                 ,', cut ',clinc,' )'  *     ----------------------------
3396    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3397                      if(mod(VIEW(iclm),2).eq.0)then
3398                         XGOOD(nplanes-ip+1)=1.
3399                         resx(nplanes-ip+1)=rxmm
3400                         if(DEBUG)print*,'%%%% included X-cl ',iclm
3401         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3402         $                    ,', dist.= ',distmin
3403         $                    ,', cut ',clinc,' )'
3404                      else
3405                         YGOOD(nplanes-ip+1)=1.
3406                         resy(nplanes-ip+1)=rymm
3407                         if(DEBUG)print*,'%%%% included Y-cl ',iclm
3408         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3409         $                    ,', dist.= ', distmin
3410         $                    ,', cut ',clinc,' )'
3411                      endif
3412    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3413    *     ----------------------------
3414                      xm_A(nplanes-ip+1) = xmm_A
3415                      ym_A(nplanes-ip+1) = ymm_A
3416                      xm_B(nplanes-ip+1) = xmm_B
3417                      ym_B(nplanes-ip+1) = ymm_B
3418                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3419                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3420                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3421    *     ----------------------------
3422                 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)  
 *              ----------------------------  
3423              endif              endif
3424           endif           endif
3425   133     continue   133     continue
# Line 2930  c              dedxtrk(nplanes-ip+1) = d Line 3438  c              dedxtrk(nplanes-ip+1) = d
3438  *                                                 *  *                                                 *
3439  *                                                 *  *                                                 *
3440  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3441  *  *
3442        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3443    
3444        include 'commontracker.f'        include 'commontracker.f'
3445        include 'level1.f'        include 'level1.f'
3446        include 'common_momanhough.f'        include 'common_momanhough.f'
3447  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
   
3448    
3449        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3450    
# Line 2952  c      include 'level1.f' Line 3454  c      include 'level1.f'
3454              if(id.ne.0)then              if(id.ne.0)then
3455                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3456                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3457  c               cl_used(iclx)=1  !tag used clusters                 cl_used(iclx)=ntrk  !tag used clusters
3458  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)  
3459              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3460  c               cl_used(icl)=1   !tag used clusters                 cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3461              endif              endif
3462                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3463  *     -----------------------------  *     -----------------------------
3464  *     remove the couple from clouds  *     remove the couple from clouds
3465  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3024  c               endif Line 3520  c               endif
3520        include 'level1.f'        include 'level1.f'
3521        include 'common_momanhough.f'        include 'common_momanhough.f'
3522        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3523    
3524    *     ---------------------------------
3525    *     variables initialized from level1
3526    *     ---------------------------------
3527        do i=1,nviews        do i=1,nviews
3528           good2(i)=good1(i)           good2(i)=good1(i)
3529             do j=1,nva1_view
3530                vkflag(i,j)=1
3531                if(cnnev(i,j).le.0)then
3532                   vkflag(i,j)=cnnev(i,j)
3533                endif
3534             enddo
3535        enddo        enddo
3536    *     ----------------
3537    *     level2 variables
3538    *     ----------------
3539        NTRK = 0        NTRK = 0
3540        do it=1,NTRKMAX        do it=1,NTRKMAX
3541           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3041  c      include 'level1.f' Line 3546  c      include 'level1.f'
3546              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3547              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3548              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3549                TAILX_nt(IP,IT) = 0
3550                TAILY_nt(IP,IT) = 0
3551                XBAD(IP,IT) = 0
3552                YBAD(IP,IT) = 0
3553              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3554              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3555                LS(IP,IT) = 0
3556              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3557              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3558              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
# Line 3175  c      include 'level1.f' Line 3685  c      include 'level1.f'
3685    
3686            
3687        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3688        include 'level1.f'        include 'level1.f'
3689        include 'common_momanhough.f'        include 'common_momanhough.f'
3690        include 'level2.f'        include 'level2.f'
3691        include 'common_mini_2.f'        include 'common_mini_2.f'
3692        real sinth,phi,pig              include 'calib.f'
3693    
3694          character*10 PFA
3695          common/FINALPFA/PFA
3696    
3697          real sinth,phi,pig
3698          integer ssensor,sladder
3699        pig=acos(-1.)        pig=acos(-1.)
3700    
3701    *     -------------------------------------
3702        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3703        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3704    *     -------------------------------------
3705        phi   = al(4)                  phi   = al(4)          
3706        sinth = al(3)                    sinth = al(3)            
3707        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3198  c      include 'level1.f' Line 3714  c      include 'level1.f'
3714       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3715        al(4) = phi                      al(4) = phi              
3716        al(3) = sinth                    al(3) = sinth            
   
3717        do i=1,5        do i=1,5
3718           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3719           do j=1,5           do j=1,5
3720              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3721           enddo           enddo
3722        enddo        enddo
3723          *     -------------------------------------      
3724        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3725           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3726           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3214  c      include 'level1.f' Line 3729  c      include 'level1.f'
3729           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3730           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3731           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3732             TAILX_nt(IP,ntr) = 0.
3733             TAILY_nt(IP,ntr) = 0.
3734           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3735           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3736           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3737           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3738           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3739           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3740           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3741         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3742         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3743         $        1. )
3744    c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3745             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3746             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3747        
3748           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3749             ay   = ayv_nt(ip,ntr)
3750             bfx  = BX_STORE(ip,IDCAND)
3751             bfy  = BY_STORE(ip,IDCAND)
3752             if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3753             if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3754             tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3755             angx     = 180.*atan(tgtemp)/acos(-1.)
3756             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3757             angy    = 180.*atan(tgtemp)/acos(-1.)
3758            
3759    c         print*,'* ',ip,bfx,bfy,angx,angy
3760    
3761             id  = CP_STORE(ip,IDCAND) ! couple id
3762           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3763             ssensor = -1
3764             sladder = -1
3765             ssensor = SENSOR_STORE(ip,IDCAND)
3766             sladder = LADDER_STORE(ip,IDCAND)
3767             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3768             LS(IP,ntr)      = ssensor+10*sladder
3769    
3770           if(id.ne.0)then           if(id.ne.0)then
3771    c           >>> is a couple
3772              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3773              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3774  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)              
3775    c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3776    c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3777    c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
3778    c$$$            ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
3779                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3780                ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3781    
3782    
3783                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3784         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3785                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3786         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3787    
3788           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3789              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3790              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl              if(mod(VIEW(icl),2).eq.0)then
3791  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)                 cltrx(ip,ntr)=icl
3792    c$$$               nnnnn = npfastrips(icl,PFA,angx)
3793    c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)
3794                   xbad(ip,ntr) = nbadstrips(4,icl)
3795    
3796                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3797                elseif(mod(VIEW(icl),2).eq.1)then
3798                   cltry(ip,ntr)=icl
3799    c$$$               nnnnn = npfastrips(icl,PFA,angy)
3800    c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)
3801                   ybad(ip,ntr) = nbadstrips(4,icl)
3802                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3803                endif
3804    
3805           endif                     endif          
3806    
3807        enddo        enddo
3808    
3809    
3810    c$$$      print*,(xgood(i),i=1,6)
3811    c$$$      print*,(ygood(i),i=1,6)
3812    c$$$      print*,(ls(i,ntr),i=1,6)
3813    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3814    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3815    c$$$      print*,'-----------------------'
3816    
3817        end        end
3818    
3819        subroutine fill_level2_siglets        subroutine fill_level2_siglets
# Line 3248  c            print*,ip,' ',cltrx(ip,ntr) Line 3825  c            print*,ip,' ',cltrx(ip,ntr)
3825  *     -------------------------------------------------------  *     -------------------------------------------------------
3826    
3827        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3828        include 'calib.f'        include 'calib.f'
3829        include 'level1.f'        include 'level1.f'
3830        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 3256  c      include 'level1.f' Line 3832  c      include 'level1.f'
3832        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3833    
3834  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3835        nclsx = 0        nclsx = 0
3836        nclsy = 0        nclsy = 0
3837    
3838        do iv = 1,nviews        do iv = 1,nviews
3839           if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)  c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3840             good2(iv) = good2(iv) + mask_view(iv)*2**8
3841        enddo        enddo
3842    
3843        do icl=1,nclstr1        do icl=1,nclstr1
# Line 3270  c      good2=1!.true. Line 3846  c      good2=1!.true.
3846              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3847                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3848                 planex(nclsx) = ip                 planex(nclsx) = ip
3849                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3850                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3851                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3852                 do is=1,2                 do is=1,2
3853  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3854                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3855                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3856                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3857                 enddo                 enddo
3858  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3285  c$$$               print*,'xs(2,nclsx)   Line 3863  c$$$               print*,'xs(2,nclsx)  
3863              else                          !=== Y views              else                          !=== Y views
3864                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3865                 planey(nclsy) = ip                 planey(nclsy) = ip
3866                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3867                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3868                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3869                 do is=1,2                 do is=1,2
3870  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3871                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3872                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3873                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3874                 enddo                 enddo
3875  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 3299  c$$$               print*,'ys(1,nclsy)   Line 3879  c$$$               print*,'ys(1,nclsy)  
3879  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3880              endif              endif
3881           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3882    
3883  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3884           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.23