/[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.20 by pam-fi, Fri Apr 27 10:39:58 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 285  c$$$         if(ibest.eq.0)goto 880 !>> Line 347  c$$$         if(ibest.eq.0)goto 880 !>>
347           do i=1,5           do i=1,5
348              AL_GUESS(i)=AL(i)              AL_GUESS(i)=AL(i)
349           enddo           enddo
350    c         print*,'## guess: ',al
351    
352           do i=1,5           do i=1,5
353              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
# Line 297  c$$$         if(ibest.eq.0)goto 880 !>> Line 360  c$$$         if(ibest.eq.0)goto 880 !>>
360           iprint=0           iprint=0
361  c         if(DEBUG)iprint=1  c         if(DEBUG)iprint=1
362           if(VERBOSE)iprint=1           if(VERBOSE)iprint=1
363             if(DEBUG)iprint=2
364           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
365           if(ifail.ne.0) then           if(ifail.ne.0) then
366              if(.true.)then              if(VERBOSE)then
367                 print *,                 print *,
368       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
369       $              ,iev       $              ,iev
# Line 446  c     $        rchi2best.lt.15..and. Line 510  c     $        rchi2best.lt.15..and.
510  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
511  *     angx   - Projected angle in x  *     angx   - Projected angle in x
512  *     angy   - Projected angle in y  *     angy   - Projected angle in y
513    *     bfx    - x-component of magnetci field
514    *     bfy    - y-component of magnetic field
515  *  *
516  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
517  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 484  c     $        rchi2best.lt.15..and. Line 550  c     $        rchi2best.lt.15..and.
550  *  *
551  *  *
552    
553        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
554    
 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*****************************************************  
555                
556        include 'commontracker.f'        include 'commontracker.f'
557        include 'level1.f'        include 'level1.f'
558        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
559        include 'common_align.f'        include 'common_align.f'
560        include 'common_mech.f'        include 'common_mech.f'
561        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
 c      include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
562    
563        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
564        integer sensor        integer sensor
565        integer viewx,viewy        integer viewx,viewy
566        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
567        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
568          real angx,angy            !X-Y effective angle
569          real bfx,bfy              !X-Y b-field components
570    
571        real stripx,stripy        real stripx,stripy
572    
573        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
574        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
575        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  
576                
577    
578        parameter (ndivx=30)        parameter (ndivx=30)
# Line 537  c      double precision xi_B,yi_B,zi_B Line 589  c      double precision xi_B,yi_B,zi_B
589        xPAM_B = 0.        xPAM_B = 0.
590        yPAM_B = 0.        yPAM_B = 0.
591        zPAM_B = 0.        zPAM_B = 0.
592    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
593    
594  *     -----------------  *     -----------------
595  *     CLUSTER X  *     CLUSTER X
596  *     -----------------  *     -----------------
597    
598        if(icx.ne.0)then        if(icx.ne.0)then
599           viewx = VIEW(icx)  
600           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
601           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
602           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
603                     resxPAM = RESXAV
604           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
605           if(PFAx.eq.'COG1')then  !(1)  *        --------------------------
606              stripx = stripx      !(1)  *        magnetic-field corrections
607              resxPAM = resxPAM    !(1)  *        --------------------------
608             angtemp  = ax
609             bfytemp  = bfy
610             if(nplx.eq.6) angtemp = -1. * ax
611             if(nplx.eq.6) bfytemp = -1. * bfy
612             tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
613             angx     = 180.*atan(tgtemp)/acos(-1.)
614             stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
615    c$$$         print*,nplx,ax,bfy/10.
616    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
617    c$$$         print*,'========================'
618    *        --------------------------
619    
620    c$$$         print*,'--- x-cl ---'
621    c$$$         istart = INDSTART(ICX)
622    c$$$         istop  = TOTCLLENGTH
623    c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
624    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
625    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
626    c$$$         print*,INDMAX(icx)
627    c$$$         print*,cog(4,icx)
628    c$$$         print*,fbad_cog(4,icx)
629            
630    
631             if(PFAx.eq.'COG1')then
632    
633                stripx  = stripx
634                resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM
635    
636           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
637              stripx = stripx + cog(2,icx)              
638                stripx  = stripx + cog(2,icx)            
639                resxPAM = risx_cog(abs(angx))!TEMPORANEO              
640              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
641    
642             elseif(PFAx.eq.'COG3')then
643    
644                stripx  = stripx + cog(3,icx)            
645                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
646                resxPAM = resxPAM*fbad_cog(3,icx)
647    
648             elseif(PFAx.eq.'COG4')then
649    
650                stripx  = stripx + cog(4,icx)            
651                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
652                resxPAM = resxPAM*fbad_cog(4,icx)
653    
654           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
655  c            cog2 = cog(2,icx)  
656  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                          stripx  = stripx + pfaeta2(icx,angx)          
657  c            stripx = stripx + etacorr              resxPAM = risx_eta2(abs(angx))
658              stripx = stripx + pfaeta2(icx,angx)            !(3)              resxPAM = resxPAM*fbad_cog(2,icx)
             resxPAM = risx_eta2(angx)                       !   (4)  
659              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
660       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
661              resxPAM = resxPAM*fbad_cog(2,icx)  
662           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                        
663              stripx = stripx + pfaeta3(icx,angx)            !(3)  
664              resxPAM = risx_eta3(angx)                       !   (4)              stripx  = stripx + pfaeta3(icx,angx)          
665              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              resxPAM = risx_eta3(abs(angx))                      
666       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
667              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)            
668           elseif(PFAx.eq.'ETA4')then                         !(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
669              stripx = stripx + pfaeta4(icx,angx)            !(3)  
670              resxPAM = risx_eta4(angx)                       !   (4)           elseif(PFAx.eq.'ETA4')then                        
671              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
672       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)              stripx  = stripx + pfaeta4(icx,angx)            
673              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = risx_eta4(abs(angx))                      
674           elseif(PFAx.eq.'ETA')then                          !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
675              stripx = stripx + pfaeta(icx,angx)             !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              
676              resxPAM = ris_eta(icx,angx)                     !   (4)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
677              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
678       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)           elseif(PFAx.eq.'ETA')then  
679  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
680              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)              stripx  = stripx + pfaeta(icx,angx)            
681           elseif(PFAx.eq.'COG')then           !(2)              resxPAM = ris_eta(icx,angx)                    
682              stripx = stripx + cog(0,icx)     !(2)                      resxPAM = resxPAM*fbad_eta(icx,angx)            
683              resxPAM = risx_cog(angx)                        !   (4)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              
684              resxPAM = resxPAM*fbad_cog(0,icx)!(2)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
685    
686             elseif(PFAx.eq.'COG')then          
687    
688                stripx  = stripx + cog(0,icx)            
689                resxPAM = risx_cog(abs(angx))                    
690                resxPAM = resxPAM*fbad_cog(0,icx)
691    
692           else           else
693              print*,'*** Non valid p.f.a. (x) --> ',PFAx              print*,'*** Non valid p.f.a. (x) --> ',PFAx
694           endif           endif
695    
696    
697    *     ======================================
698    *     temporary patch for saturated clusters
699    *     ======================================
700             if( nsatstrips(icx).gt.0 )then
701                stripx  = stripx + cog(4,icx)            
702                resxPAM = pitchX*1e-4/sqrt(12.)
703                cc=cog(4,icx)
704    c$$$            print*,icx,' *** ',cc
705    c$$$            print*,icx,' *** ',resxPAM
706             endif
707    
708        endif        endif
 c      if(icy.eq.0.and.icx.ne.0)  
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
709                
710  *     -----------------  *     -----------------
711  *     CLUSTER Y  *     CLUSTER Y
712  *     -----------------  *     -----------------
713    
714        if(icy.ne.0)then        if(icy.ne.0)then
715    
716           viewy = VIEW(icy)           viewy = VIEW(icy)
717           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
718           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
719           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
720             stripy = float(MAXS(icy))
721    
722           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
723              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
724       $           ,icx,icy       $           ,icx,icy
725              goto 100              goto 100
726           endif           endif
727            *        --------------------------
728           stripy = float(MAXS(icy))  *        magnetic-field corrections
729           if(PFAy.eq.'COG1')then !(1)  *        --------------------------
730              stripy = stripy     !(1)           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
731              resyPAM = resyPAM   !(1)           angy    = 180.*atan(tgtemp)/acos(-1.)
732             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
733    *        --------------------------
734            
735    c$$$         print*,'--- y-cl ---'
736    c$$$         istart = INDSTART(ICY)
737    c$$$         istop  = TOTCLLENGTH
738    c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
739    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
740    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
741    c$$$         print*,INDMAX(icy)
742    c$$$         print*,cog(4,icy)
743    c$$$         print*,fbad_cog(4,icy)
744    
745             if(PFAy.eq.'COG1')then
746    
747                stripy  = stripy    
748                resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  
749    
750           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
751              stripy = stripy + cog(2,icy)  
752                stripy  = stripy + cog(2,icy)
753                resyPAM = risy_cog(abs(angy))!TEMPORANEO
754              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
755    
756             elseif(PFAy.eq.'COG3')then
757    
758                stripy  = stripy + cog(3,icy)
759                resyPAM = risy_cog(abs(angy))!TEMPORANEO
760                resyPAM = resyPAM*fbad_cog(3,icy)
761    
762             elseif(PFAy.eq.'COG4')then
763    
764                stripy  = stripy + cog(4,icy)
765                resyPAM = risy_cog(abs(angy))!TEMPORANEO
766                resyPAM = resyPAM*fbad_cog(4,icy)
767    
768           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
769  c            cog2 = cog(2,icy)  
770  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)              stripy  = stripy + pfaeta2(icy,angy)          
771  c            stripy = stripy + etacorr              resyPAM = risy_eta2(abs(angy))              
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
772              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
773              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
774       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
775           elseif(PFAy.eq.'ETA3')then                         !(3)  
776              stripy = stripy + pfaeta3(icy,angy)            !(3)           elseif(PFAy.eq.'ETA3')then                      
777              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
778              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              stripy  = stripy + pfaeta3(icy,angy)
779       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)              resyPAM = resyPAM*fbad_cog(3,icy)  
780           elseif(PFAy.eq.'ETA4')then                         !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)
781              stripy = stripy + pfaeta4(icy,angy)            !(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
782              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
783              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA4')then  
784       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
785           elseif(PFAy.eq.'ETA')then                          !(3)              stripy  = stripy + pfaeta4(icy,angy)
786              stripy = stripy + pfaeta(icy,angy)             !(3)              resyPAM = resyPAM*fbad_cog(4,icy)
787              resyPAM = ris_eta(icy,angy)                     !   (4)              if(DEBUG.and.fbad_cog(4,icy).ne.1)
788  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
789              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
790              if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA')then
791       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
792                stripy  = stripy + pfaeta(icy,angy)
793                resyPAM = ris_eta(icy,angy)  
794                resyPAM = resyPAM*fbad_eta(icy,angy)
795                if(DEBUG.and.fbad_cog(2,icy).ne.1)
796         $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
797    
798           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
799              stripy = stripy + cog(0,icy)              
800              resyPAM = risy_cog(angy)                        !   (4)              stripy  = stripy + cog(0,icy)            
801  c            resyPAM = ris_eta(icy,angy)                    !   (4)              resyPAM = risy_cog(abs(angy))
802              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
803    
804           else           else
805              print*,'*** Non valid p.f.a. (x) --> ',PFAx              print*,'*** Non valid p.f.a. (x) --> ',PFAx
806           endif           endif
807    
808    
809    *     ======================================
810    *     temporary patch for saturated clusters
811    *     ======================================
812             if( nsatstrips(icy).gt.0 )then
813                stripy  = stripy + cog(4,icy)            
814                resyPAM = pitchY*1e-4/sqrt(12.)
815                cc=cog(4,icy)
816    c$$$            print*,icy,' *** ',cc
817    c$$$            print*,icy,' *** ',resyPAM
818             endif
819    
820    
821        endif        endif
822    
823          c      print*,'## stripx,stripy ',stripx,stripy
824    
825  c===========================================================  c===========================================================
826  C     COUPLE  C     COUPLE
827  C===========================================================  C===========================================================
# Line 858  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1023  c         print*,'A-(',xPAM_A,yPAM_A,')
1023                            
1024        endif        endif
1025                    
1026    
1027    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1028    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1029    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1030    
1031   100  continue   100  continue
1032        end        end
1033    
# Line 936  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1106  c         print*,'A-(',xPAM_A,yPAM_A,')
1106           endif                   endif        
1107    
1108           distance=           distance=
1109       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1110    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1111           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1112    
1113  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 1132  c$$$         print*,' resolution ',re
1132  *     ----------------------  *     ----------------------
1133                    
1134           distance=           distance=
1135       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1136       $        +       $       +
1137       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1138    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1139    c$$$     $        +
1140    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1141           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1142    
1143  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1348  c      include 'level1.f' Line 1522  c      include 'level1.f'
1522  *     ----------------------------------------------------  *     ----------------------------------------------------
1523  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1524  *     ----------------------------------------------------  *     ----------------------------------------------------
1525           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1526              cl_single(icx)=0              cl_single(icx)=0
1527              goto 10              goto 10
1528           endif           endif
# Line 1398  c     endif Line 1572  c     endif
1572  *     ----------------------------------------------------  *     ----------------------------------------------------
1573  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1574  *     ----------------------------------------------------  *     ----------------------------------------------------
1575              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1576                 cl_single(icy)=0                 cl_single(icy)=0
1577                 goto 20                 goto 20
1578              endif              endif
# Line 1444  c     endif Line 1618  c     endif
1618  *     charge correlation  *     charge correlation
1619  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1620    
1621                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1622       $              .and.       $              .and.
1623       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1624       $              .and.       $              .and.
1625       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1626       $              .and.       $              .and.
1627       $              .true.)then       $              .true.)then
1628    
1629                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1630       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1631                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1632    
1633  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1634    
1635                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1636       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1637                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1638                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1467  c                  cut = chcut * sch(npl Line 1641  c                  cut = chcut * sch(npl
1641                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1642                    endif                    endif
1643                 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  
1644    
1645                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1646                    if(verbose)print*,                    if(verbose)print*,
1647       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1648       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
# Line 1483  c                  cut = chcut * sch(npl Line 1650  c                  cut = chcut * sch(npl
1650       $                 ,'( ',ncouplemax,' ) --> masked!'       $                 ,'( ',ncouplemax,' ) --> masked!'
1651                    mask_view(nviewx(nplx)) = 2                    mask_view(nviewx(nplx)) = 2
1652                    mask_view(nviewy(nply)) = 2                    mask_view(nviewy(nply)) = 2
1653                      goto 10
1654                 endif                 endif
1655                  
1656    *     ------------------> COUPLE <------------------
1657                   ncp_plane(nplx) = ncp_plane(nplx) + 1
1658                   clx(nplx,ncp_plane(nplx))=icx
1659                   cly(nply,ncp_plane(nplx))=icy
1660                   cl_single(icx)=0
1661                   cl_single(icy)=0
1662  *     ----------------------------------------------  *     ----------------------------------------------
1663    
1664              endif                                            endif                              
# Line 1592  c      double precision xm3,ym3,zm3 Line 1767  c      double precision xm3,ym3,zm3
1767                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1768                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1769  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1770                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1771                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1772                 xm1=xPAM                 xm1=xPAM
1773                 ym1=yPAM                 ym1=yPAM
1774                 zm1=zPAM                                   zm1=zPAM                  
# Line 1608  c     print*,'***',is1,xm1,ym1,zm1 Line 1784  c     print*,'***',is1,xm1,ym1,zm1
1784                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1785  c                        call xyz_PAM  c                        call xyz_PAM
1786  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1787    c                        call xyz_PAM
1788    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1789                          call xyz_PAM                          call xyz_PAM
1790       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1791                          xm2=xPAM                          xm2=xPAM
1792                          ym2=yPAM                          ym2=yPAM
1793                          zm2=zPAM                          zm2=zPAM
# Line 1673  c$$$ Line 1851  c$$$
1851                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
1852  c                                 call xyz_PAM  c                                 call xyz_PAM
1853  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
1854    c                                 call xyz_PAM
1855    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
1856                                   call xyz_PAM                                   call xyz_PAM
1857       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
1858         $                                ,0.,0.,0.,0.)
1859                                   xm3=xPAM                                   xm3=xPAM
1860                                   ym3=yPAM                                   ym3=yPAM
1861                                   zm3=zPAM                                   zm3=zPAM
# Line 2146  c     print*,'check cp_used' Line 2327  c     print*,'check cp_used'
2327           enddo           enddo
2328  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2329           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2330           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2331                    
2332  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2333  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
# Line 2223  c$$$     $           ,(tr_cloud(iii),iii Line 2404  c$$$     $           ,(tr_cloud(iii),iii
2404  **************************************************  **************************************************
2405    
2406        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2407    
2408        include 'commontracker.f'        include 'commontracker.f'
2409        include 'level1.f'        include 'level1.f'
# Line 2233  c*************************************** Line 2411  c***************************************
2411        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2412        include 'common_mini_2.f'        include 'common_mini_2.f'
2413        include 'common_mech.f'        include 'common_mech.f'
2414  c      include 'momanhough_init.f'  
2415    
2416    
2417  *     output flag  *     output flag
# Line 2323  c      real fitz(nplanes)        !z coor Line 2501  c      real fitz(nplanes)        !z coor
2501                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2502              enddo              enddo
2503                            
2504              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2505                if(nplused.lt.nplyz_min)goto 888 !next doublet
2506              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2507                            
2508              if(DEBUG)then              if(DEBUG)then
# Line 2427  c$$$            if(AL_INI(5).gt.defmax)g Line 2606  c$$$            if(AL_INI(5).gt.defmax)g
2606  *                                   *************************  *                                   *************************
2607  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2608  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2609    c                                    call xyz_PAM(icx,icy,is, !(1)
2610    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2611                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2612       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2613  *                                   *************************  *                                   *************************
2614  *                                   -----------------------------  *                                   -----------------------------
2615                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 2456  c$$$                              enddo Line 2637  c$$$                              enddo
2637                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2638                                iprint=0                                iprint=0
2639  c                              if(DEBUG)iprint=1  c                              if(DEBUG)iprint=1
2640                                if(DEBUG)iprint=1                                if(DEBUG)iprint=2
2641                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2642                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2643                                   if(DEBUG)then                                   if(DEBUG)then
# Line 2499  c                                 goto 8 Line 2680  c                                 goto 8
2680                                                                
2681                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2682                                                                
 c$$$                              ndof=0                                  
2683                                do ip=1,nplanes                                do ip=1,nplanes
2684  c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2685                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2686                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2687                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 2522  c$$$     $                               Line 2700  c$$$     $                              
2700                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2701                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2702       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2703                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2704         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2705                                        LADDER_STORE(nplanes-ip+1,ntracks)
2706         $                                   = LADDER(
2707         $                                   clx(ip,icp_cp(
2708         $                                   cp_match(ip,hit_plane(ip)
2709         $                                   ))));
2710                                   else                                   else
2711                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2712                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2713                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2714                                   endif                                   endif
2715                                     BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2716                                     BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2717                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(nplanes-ip+1,ntracks)=0
2718                                   do i=1,5                                   do i=1,5
2719                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2720                                   enddo                                   enddo
2721                                enddo                                enddo
2722                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2723                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2724                                                                
2725  *     --------------------------------  *     --------------------------------
# Line 2558  c$$$  rchi2=chi2/dble(ndof) Line 2743  c$$$  rchi2=chi2/dble(ndof)
2743           return           return
2744        endif        endif
2745                
2746    c$$$      if(DEBUG)then
2747    c$$$         print*,'****** TRACK CANDIDATES ***********'
2748    c$$$         print*,'#         R. chi2        RIG'
2749    c$$$         do i=1,ntracks
2750    c$$$            print*,i,' --- ',rchi2_store(i),' --- '
2751    c$$$     $           ,1./abs(AL_STORE(5,i))
2752    c$$$         enddo
2753    c$$$         print*,'***********************************'
2754    c$$$      endif
2755        if(DEBUG)then        if(DEBUG)then
2756           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
2757           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
2758           do i=1,ntracks          do i=1,ntracks
2759              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
2760       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
2761           enddo              ndof=ndof           !(1)
2762           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
2763         $           +int(ygood_store(ii,i)) !(1)
2764              enddo                 !(1)
2765              print*,i,' --- ',rchi2_store(i),' --- '
2766         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
2767            enddo
2768            print*,'*****************************************'
2769        endif        endif
2770                
2771                
# Line 2584  c$$$  rchi2=chi2/dble(ndof) Line 2784  c$$$  rchi2=chi2/dble(ndof)
2784    
2785        subroutine refine_track(ibest)        subroutine refine_track(ibest)
2786    
 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******************************************************  
2787    
2788        include 'commontracker.f'        include 'commontracker.f'
2789        include 'level1.f'        include 'level1.f'
# Line 2596  c*************************************** Line 2791  c***************************************
2791        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2792        include 'common_mini_2.f'        include 'common_mini_2.f'
2793        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
2794        include 'calib.f'        include 'calib.f'
2795    
   
2796  *     flag to chose PFA  *     flag to chose PFA
2797        character*10 PFA        character*10 PFA
2798        common/FINALPFA/PFA        common/FINALPFA/PFA
2799    
2800          real k(6)
2801          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
2802    
2803          real xp,yp,zp
2804          real xyzp(3),bxyz(3)
2805          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
2806    
2807  *     =================================================  *     =================================================
2808  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
2809  *                          and  *                          and
# Line 2613  c      include 'level1.f' Line 2812  c      include 'level1.f'
2812        call track_init        call track_init
2813        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2814    
2815             xP=XV_STORE(nplanes-ip+1,ibest)
2816             yP=YV_STORE(nplanes-ip+1,ibest)
2817             zP=ZV_STORE(nplanes-ip+1,ibest)
2818             call gufld(xyzp,bxyz)
2819             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
2820             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
2821    c$$$  bxyz(1)=0
2822    c$$$         bxyz(2)=0
2823    c$$$         bxyz(3)=0
2824  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
2825  *     -------------------------------------------------  *     -------------------------------------------------
2826  *     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 2841  c      include 'level1.f'
2841       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
2842              icx=clx(ip,icp)              icx=clx(ip,icp)
2843              icy=cly(ip,icp)              icy=cly(ip,icp)
2844    c            call xyz_PAM(icx,icy,is,
2845    c     $           PFA,PFA,
2846    c     $           AXV_STORE(nplanes-ip+1,ibest),
2847    c     $           AYV_STORE(nplanes-ip+1,ibest))
2848              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
2849       $           PFA,PFA,       $           PFA,PFA,
2850       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
2851       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
2852  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
2853  c$$$  $              'COG2','COG2',       $           bxyz(2)
2854  c$$$  $              0.,       $           )
2855  c$$$  $              0.)  
2856              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
2857              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
2858              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 2650  c$$$  $              0.) Line 2861  c$$$  $              0.)
2861              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
2862              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
2863    
2864  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
2865              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)  
2866                            
2867  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
2868  *     -------------------------------------------------  *     -------------------------------------------------
# Line 2667  c            dedxtrk(nplanes-ip+1) = (de Line 2877  c            dedxtrk(nplanes-ip+1) = (de
2877                                
2878  *     --------------------------------------------------------------  *     --------------------------------------------------------------
2879  *     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)  
2880              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
2881  *     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
2882              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
2883    
2884                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
2885                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
2886  *     --------------------------------------------------------------  *     --------------------------------------------------------------
2887    
2888              if(DEBUG)then              if(DEBUG)then
# Line 2709  c     $              cl_used(icy).eq.1.o Line 2919  c     $              cl_used(icy).eq.1.o
2919  *            *          
2920                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
2921       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
2922       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
2923       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
2924         $              bxyz(1),
2925         $              bxyz(2)
2926         $              )
2927                                
2928                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2929                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
2930                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
2931                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
2932       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
2933                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
2934                    xmm = xPAM                    xmm = xPAM
2935                    ymm = yPAM                    ymm = yPAM
# Line 2726  c     $              'ETA2','ETA2', Line 2938  c     $              'ETA2','ETA2',
2938                    rymm = resyPAM                    rymm = resyPAM
2939                    distmin = distance                    distmin = distance
2940                    idm = id                                      idm = id                  
2941  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2942                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2943                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
2944                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
2945         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
2946                 endif                 endif
2947   1188          continue   1188          continue
2948              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
2949              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
2950                if(distmin.le.clincnewc)then     !QUIQUI              
2951  *              -----------------------------------  *              -----------------------------------
2952                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
2953                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
2954                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
2955                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
2956                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
2957                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
2958                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
2959  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
2960                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
2961  *              -----------------------------------  *              -----------------------------------
2962                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
2963                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG)print*,'%%%% included couple ',idm
2964       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
2965                 goto 133         !next plane                 goto 133         !next plane
2966              endif              endif
2967  *     ================================================  *     ================================================
# Line 2780  c            if(DEBUG)print*,'>>>> try t Line 2994  c            if(DEBUG)print*,'>>>> try t
2994  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
2995                 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)
2996  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
2997    c               call xyz_PAM(icx,0,ist,
2998    c     $              PFA,PFA,
2999    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3000                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3001       $              PFA,PFA,       $              PFA,PFA,
3002       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3003         $              bxyz(1),
3004         $              bxyz(2)
3005         $              )              
3006                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3007                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3008                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3009       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3010                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3011                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3012                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2799  c     $              'ETA2','ETA2', Line 3018  c     $              'ETA2','ETA2',
3018                    rymm = resyPAM                    rymm = resyPAM
3019                    distmin = distance                    distmin = distance
3020                    iclm = icx                    iclm = icx
3021  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3022                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3023                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3024                 endif                                   endif                  
3025  11881          continue  11881          continue
# Line 2808  c                  dedxmm = dedx(icx) !( Line 3027  c                  dedxmm = dedx(icx) !(
3027  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
3028                 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)
3029  *                                              !jump to the next couple  *                                              !jump to the next couple
3030    c               call xyz_PAM(0,icy,ist,
3031    c     $              PFA,PFA,
3032    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3033                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3034       $              PFA,PFA,       $              PFA,PFA,
3035       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3036         $              bxyz(1),
3037         $              bxyz(2)
3038         $              )
3039                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3040                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3041                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3042       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3043                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3044                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3045                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2827  c     $              'ETA2','ETA2', Line 3051  c     $              'ETA2','ETA2',
3051                    rymm = resyPAM                    rymm = resyPAM
3052                    distmin = distance                    distmin = distance
3053                    iclm = icy                    iclm = icy
3054  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3055                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3056                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3057                 endif                                   endif                  
3058  11882          continue  11882          continue
3059              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3060  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3061    c            print*,'## ncls(',ip,') ',ncls(ip)
3062              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3063                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3064  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 3067  c               if(cl_used(icl).eq.1.or.
3067       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3068                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3069                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3070       $                 PFA,PFA,       $                 PFA,PFA,
3071       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3072         $                 bxyz(1),
3073         $                 bxyz(2)
3074         $                 )
3075                 else                         !<---- Y view                 else                         !<---- Y view
3076                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3077       $                 PFA,PFA,       $                 PFA,PFA,
3078       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3079         $                 bxyz(1),
3080         $                 bxyz(2)
3081         $                 )
3082                 endif                 endif
3083    
3084                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3085                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3086                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3087       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance,'<',distmin,' ?'
3088                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3089                      if(DEBUG)print*,'YES'
3090                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3091                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3092                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 2867  c     $                 'ETA2','ETA2', Line 3097  c     $                 'ETA2','ETA2',
3097                    rymm = resyPAM                    rymm = resyPAM
3098                    distmin = distance                      distmin = distance  
3099                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3100                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3101                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3102                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3103                    else          !<---- Y view                    else          !<---- Y view
3104                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3105                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3106                    endif                    endif
3107                 endif                                   endif                  
3108  18882          continue  18882          continue
3109              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3110    c            print*,'## distmin ', distmin,' clinc ',clinc
3111    
3112              if(distmin.le.clinc)then                    c     QUIQUI------------
3113                  c     anche qui: non ci vuole clinc???
3114                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3115                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3116                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3117                    resx(nplanes-ip+1)=rxmm       $                 20*
3118                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3119  c                  if(.true.)print*,'%%%% included X-cl ',iclm                 else if(mod(VIEW(iclm),2).ne.0)then
3120       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    clincnew=
3121       $                 ,', norm.dist.= ',distmin       $                 10*
3122       $                 ,', cut ',clinc,' )'       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3123                 else                 endif
3124                    YGOOD(nplanes-ip+1)=1.  c     QUIQUI------------
3125                    resy(nplanes-ip+1)=rymm                
3126                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 if(distmin.le.clincnew)then   !QUIQUI
3127  c                  if(.true.)print*,'%%%% included Y-cl ',iclm  c     if(distmin.le.clinc)then          !QUIQUI          
3128       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    
3129       $                 ,', norm.dist.= ', distmin                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3130       $                 ,', cut ',clinc,' )'  *     ----------------------------
3131    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3132                      if(mod(VIEW(iclm),2).eq.0)then
3133                         XGOOD(nplanes-ip+1)=1.
3134                         resx(nplanes-ip+1)=rxmm
3135                         if(DEBUG)print*,'%%%% included X-cl ',iclm
3136         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3137         $                    ,', dist.= ',distmin
3138         $                    ,', cut ',clinc,' )'
3139                      else
3140                         YGOOD(nplanes-ip+1)=1.
3141                         resy(nplanes-ip+1)=rymm
3142                         if(DEBUG)print*,'%%%% included Y-cl ',iclm
3143         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3144         $                    ,', dist.= ', distmin
3145         $                    ,', cut ',clinc,' )'
3146                      endif
3147    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3148    *     ----------------------------
3149                      xm_A(nplanes-ip+1) = xmm_A
3150                      ym_A(nplanes-ip+1) = ymm_A
3151                      xm_B(nplanes-ip+1) = xmm_B
3152                      ym_B(nplanes-ip+1) = ymm_B
3153                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3154                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3155                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3156    *     ----------------------------
3157                 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)  
 *              ----------------------------  
3158              endif              endif
3159           endif           endif
3160   133     continue   133     continue
# Line 2930  c              dedxtrk(nplanes-ip+1) = d Line 3173  c              dedxtrk(nplanes-ip+1) = d
3173  *                                                 *  *                                                 *
3174  *                                                 *  *                                                 *
3175  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3176  *  *
3177        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3178    
3179        include 'commontracker.f'        include 'commontracker.f'
3180        include 'level1.f'        include 'level1.f'
3181        include 'common_momanhough.f'        include 'common_momanhough.f'
3182  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
   
3183    
3184        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3185    
# Line 2952  c      include 'level1.f' Line 3189  c      include 'level1.f'
3189              if(id.ne.0)then              if(id.ne.0)then
3190                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3191                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3192  c               cl_used(iclx)=1  !tag used clusters                 cl_used(iclx)=ntrk  !tag used clusters
3193  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)  
3194              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3195  c               cl_used(icl)=1   !tag used clusters                 cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3196              endif              endif
3197                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3198  *     -----------------------------  *     -----------------------------
3199  *     remove the couple from clouds  *     remove the couple from clouds
3200  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3024  c               endif Line 3255  c               endif
3255        include 'level1.f'        include 'level1.f'
3256        include 'common_momanhough.f'        include 'common_momanhough.f'
3257        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3258    
3259    *     ---------------------------------
3260    *     variables initialized from level1
3261    *     ---------------------------------
3262        do i=1,nviews        do i=1,nviews
3263           good2(i)=good1(i)           good2(i)=good1(i)
3264             do j=1,nva1_view
3265                vkflag(i,j)=1
3266                if(cnnev(i,j).le.0)then
3267                   vkflag(i,j)=cnnev(i,j)
3268                endif
3269             enddo
3270        enddo        enddo
3271    *     ----------------
3272    *     level2 variables
3273    *     ----------------
3274        NTRK = 0        NTRK = 0
3275        do it=1,NTRKMAX        do it=1,NTRKMAX
3276           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3041  c      include 'level1.f' Line 3281  c      include 'level1.f'
3281              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3282              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3283              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3284                TAILX_nt(IP,IT) = 0
3285                TAILY_nt(IP,IT) = 0
3286                XBAD(IP,IT) = 0
3287                YBAD(IP,IT) = 0
3288              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3289              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3290                LS(IP,IT) = 0
3291              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3292              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3293              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
# Line 3175  c      include 'level1.f' Line 3420  c      include 'level1.f'
3420    
3421            
3422        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3423        include 'level1.f'        include 'level1.f'
3424        include 'common_momanhough.f'        include 'common_momanhough.f'
3425        include 'level2.f'        include 'level2.f'
3426        include 'common_mini_2.f'        include 'common_mini_2.f'
3427        real sinth,phi,pig              include 'calib.f'
3428    
3429          character*10 PFA
3430          common/FINALPFA/PFA
3431    
3432          real sinth,phi,pig
3433          integer ssensor,sladder
3434        pig=acos(-1.)        pig=acos(-1.)
3435    
3436    *     -------------------------------------
3437        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3438        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3439    *     -------------------------------------
3440        phi   = al(4)                  phi   = al(4)          
3441        sinth = al(3)                    sinth = al(3)            
3442        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3198  c      include 'level1.f' Line 3449  c      include 'level1.f'
3449       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3450        al(4) = phi                      al(4) = phi              
3451        al(3) = sinth                    al(3) = sinth            
   
3452        do i=1,5        do i=1,5
3453           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3454           do j=1,5           do j=1,5
3455              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3456           enddo           enddo
3457        enddo        enddo
3458          *     -------------------------------------      
3459        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3460           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3461           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3214  c      include 'level1.f' Line 3464  c      include 'level1.f'
3464           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3465           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3466           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3467             TAILX_nt(IP,ntr) = 0.
3468             TAILY_nt(IP,ntr) = 0.
3469           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3470           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3471           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3472           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3473           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3474           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3475           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3476         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3477         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3478         $        1. )
3479    c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3480             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3481             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3482        
3483           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3484             ay   = ayv_nt(ip,ntr)
3485             bfx  = BX_STORE(ip,IDCAND)
3486             bfy  = BY_STORE(ip,IDCAND)
3487             if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3488             if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3489             tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3490             angx     = 180.*atan(tgtemp)/acos(-1.)
3491             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3492             angy    = 180.*atan(tgtemp)/acos(-1.)
3493            
3494    c         print*,'* ',ip,bfx,bfy,angx,angy
3495    
3496             id  = CP_STORE(ip,IDCAND) ! couple id
3497           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3498             ssensor = -1
3499             sladder = -1
3500             ssensor = SENSOR_STORE(ip,IDCAND)
3501             sladder = LADDER_STORE(ip,IDCAND)
3502             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3503             LS(IP,ntr)      = ssensor+10*sladder
3504    
3505           if(id.ne.0)then           if(id.ne.0)then
3506    c           >>> is a couple
3507              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3508              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3509  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3510    c$$$            if(is_cp(id).ne.ssensor)
3511    c$$$     $           print*,'ERROR is sensor assignment (couple)'
3512    c$$$     $           ,is_cp(id),ssensor
3513    c$$$            if(LADDER(clx(nplanes-ip+1,icp_cp(id))).ne.sladder)
3514    c$$$     $           print*,'ERROR is ladder assignment (couple)'
3515    c$$$     $           ,LADDER(clx(nplanes-ip+1,icp_cp(id))),sladder
3516                
3517                nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3518                nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3519                xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
3520                ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
3521    
3522                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3523         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3524                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3525         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3526    
3527           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3528              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  c           >>> is a singlet
3529              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl  c$$$            if(LADDER(icl).ne.sladder)
3530  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  c$$$     $           print*,'ERROR is ladder assignment (single)'
3531    c$$$     $           ,LADDER(icl),sladder
3532                if(mod(VIEW(icl),2).eq.0)then
3533                   cltrx(ip,ntr)=icl
3534                   nnnnn = npfastrips(icl,PFA,angx)
3535                   xbad(ip,ntr) = nbadstrips(nnnnn,icl)
3536                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3537                elseif(mod(VIEW(icl),2).eq.1)then
3538                   cltry(ip,ntr)=icl
3539                   nnnnn = npfastrips(icl,PFA,angy)
3540                   ybad(ip,ntr) = nbadstrips(nnnnn,icl)
3541                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3542                endif
3543           endif                     endif          
3544    
3545        enddo        enddo
3546    
3547    
3548    c$$$      print*,(xgood(i),i=1,6)
3549    c$$$      print*,(ygood(i),i=1,6)
3550    c$$$      print*,(ls(i,ntr),i=1,6)
3551    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3552    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3553    c$$$      print*,'-----------------------'
3554    
3555        end        end
3556    
3557        subroutine fill_level2_siglets        subroutine fill_level2_siglets
# Line 3248  c            print*,ip,' ',cltrx(ip,ntr) Line 3563  c            print*,ip,' ',cltrx(ip,ntr)
3563  *     -------------------------------------------------------  *     -------------------------------------------------------
3564    
3565        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3566        include 'calib.f'        include 'calib.f'
3567        include 'level1.f'        include 'level1.f'
3568        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 3256  c      include 'level1.f' Line 3570  c      include 'level1.f'
3570        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3571    
3572  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3573        nclsx = 0        nclsx = 0
3574        nclsy = 0        nclsy = 0
3575    
# Line 3270  c      good2=1!.true. Line 3583  c      good2=1!.true.
3583              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3584                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3585                 planex(nclsx) = ip                 planex(nclsx) = ip
3586                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3587                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3588                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3589                 do is=1,2                 do is=1,2
3590  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3591                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3592                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3593                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3594                 enddo                 enddo
3595  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3285  c$$$               print*,'xs(2,nclsx)   Line 3600  c$$$               print*,'xs(2,nclsx)  
3600              else                          !=== Y views              else                          !=== Y views
3601                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3602                 planey(nclsy) = ip                 planey(nclsy) = ip
3603                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3604                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3605                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3606                 do is=1,2                 do is=1,2
3607  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3608                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3609                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3610                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3611                 enddo                 enddo
3612  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 3299  c$$$               print*,'ys(1,nclsy)   Line 3616  c$$$               print*,'ys(1,nclsy)  
3616  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3617              endif              endif
3618           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3619    
3620  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3621           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23