/[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.29 by pam-fi, Tue Aug 21 15:21:52 2007 UTC revision 1.34 by pam-fi, Wed Mar 5 17:00:20 2008 UTC
# Line 688  c     $        rchi2best.lt.15..and. Line 688  c     $        rchi2best.lt.15..and.
688    
689        real stripx,stripy        real stripx,stripy
690    
691          double precision xi,yi,zi
692          double precision xi_A,yi_A,zi_A
693          double precision xi_B,yi_B,zi_B
694        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
695        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
696        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
# Line 701  c$$$      print*,icx,icy,sensor,PFAx,PFA Line 704  c$$$      print*,icx,icy,sensor,PFAx,PFA
704        resxPAM = 0        resxPAM = 0
705        resyPAM = 0        resyPAM = 0
706    
707        xPAM = 0.        xPAM = 0.D0
708        yPAM = 0.        yPAM = 0.D0
709        zPAM = 0.        zPAM = 0.D0
710        xPAM_A = 0.        xPAM_A = 0.D0
711        yPAM_A = 0.        yPAM_A = 0.D0
712        zPAM_A = 0.        zPAM_A = 0.D0
713        xPAM_B = 0.        xPAM_B = 0.D0
714        yPAM_B = 0.        yPAM_B = 0.D0
715        zPAM_B = 0.        zPAM_B = 0.D0
716  c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy  cc      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
717    
718        if(sensor.lt.1.or.sensor.gt.2)then        if(sensor.lt.1.or.sensor.gt.2)then
719           print*,'xyz_PAM   ***ERROR*** wrong input '           print*,'xyz_PAM   ***ERROR*** wrong input '
# Line 727  c      print*,'## xyz_PAM: ',icx,icy,sen Line 730  c      print*,'## xyz_PAM: ',icx,icy,sen
730           viewx   = VIEW(icx)           viewx   = VIEW(icx)
731           nldx    = nld(MAXS(icx),VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
732           nplx    = npl(VIEW(icx))           nplx    = npl(VIEW(icx))
733           resxPAM = RESXAV  c         resxPAM = RESXAV
734           stripx  = float(MAXS(icx))           stripx  = float(MAXS(icx))
735    
736           if(           if(
# Line 747  c      print*,'## xyz_PAM: ',icx,icy,sen Line 750  c      print*,'## xyz_PAM: ',icx,icy,sen
750  *        --------------------------  *        --------------------------
751  *        magnetic-field corrections  *        magnetic-field corrections
752  *        --------------------------  *        --------------------------
753           angtemp  = ax           stripx  = stripx +  fieldcorr(viewx,bfy)      
754           bfytemp  = bfy           angx    = effectiveangle(ax,viewx,bfy)
 *        /////////////////////////////////  
 *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!  
 *        *grvzkkjsdgjhhhgngbn###>:(  
 *        /////////////////////////////////  
 c         if(nplx.eq.6) angtemp = -1. * ax  
 c         if(nplx.eq.6) bfytemp = -1. * bfy  
          if(viewx.eq.12) angtemp = -1. * ax  
          if(viewx.eq.12) bfytemp = -1. * bfy  
          tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001  
          angx     = 180.*atan(tgtemp)/acos(-1.)  
          stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX  
 c$$$         print*,nplx,ax,bfy/10.  
 c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX  
 c$$$         print*,'========================'  
 c$$$         if(bfy.ne.0.)print*,viewx,'-x- '  
 c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ  
 *        --------------------------  
   
 c$$$         print*,'--- x-cl ---'  
 c$$$         istart = INDSTART(ICX)  
 c$$$         istop  = TOTCLLENGTH  
 c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1  
 c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)  
 c$$$         print*,(CLSIGNAL(i),i=istart,istop)  
 c$$$         print*,INDMAX(icx)  
 c$$$         print*,cog(4,icx)  
 c$$$         print*,fbad_cog(4,icx)  
755                    
756             call applypfa(PFAx,icx,angx,corr,res)
757           if(PFAx.eq.'COG1')then           stripx  = stripx + corr
758             resxPAM = res
             stripx  = stripx  
             resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM  
   
          elseif(PFAx.eq.'COG2')then  
   
             stripx  = stripx + cog(2,icx)              
             resxPAM = risx_cog(abs(angx))!TEMPORANEO                
             resxPAM = resxPAM*fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'COG3')then  
   
             stripx  = stripx + cog(3,icx)              
             resxPAM = risx_cog(abs(angx))!TEMPORANEO                        
             resxPAM = resxPAM*fbad_cog(3,icx)  
   
          elseif(PFAx.eq.'COG4')then  
   
             stripx  = stripx + cog(4,icx)              
             resxPAM = risx_cog(abs(angx))!TEMPORANEO                        
             resxPAM = resxPAM*fbad_cog(4,icx)  
   
          elseif(PFAx.eq.'ETA2')then  
   
             stripx  = stripx + pfaeta2(icx,angx)            
             resxPAM = risxeta2(abs(angx))  
             resxPAM = resxPAM*fbad_cog(2,icx)  
 c$$$            if(DEBUG.EQ.1.and.fbad_cog(2,icx).ne.1)  
 c$$$     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'ETA3')then                          
   
             stripx  = stripx + pfaeta3(icx,angx)            
             resxPAM = risxeta3(abs(angx))                        
             resxPAM = resxPAM*fbad_cog(3,icx)                
 c            if(DEBUG.and.fbad_cog(3,icx).ne.1)              
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)  
   
          elseif(PFAx.eq.'ETA4')then                          
   
             stripx  = stripx + pfaeta4(icx,angx)              
             resxPAM = risxeta4(abs(angx))                        
             resxPAM = resxPAM*fbad_cog(4,icx)                
 c            if(DEBUG.and.fbad_cog(4,icx).ne.1)                
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)  
   
          elseif(PFAx.eq.'ETA')then    
   
             stripx  = stripx + pfaeta(icx,angx)              
 c            resxPAM = riseta(icx,angx)                      
             resxPAM = riseta(viewx,angx)                      
             resxPAM = resxPAM*fbad_eta(icx,angx)              
 c            if(DEBUG.and.fbad_cog(2,icx).ne.1)                
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'ETAL')then    
   
             stripx  = stripx + pfaetal(icx,angx)              
             resxPAM = riseta(viewx,angx)                      
             resxPAM = resxPAM*fbad_eta(icx,angx)              
 c            if(DEBUG.and.fbad_cog(2,icx).ne.1)                
 c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
   
          elseif(PFAx.eq.'COG')then            
   
             stripx  = stripx + cog(0,icx)              
             resxPAM = risx_cog(abs(angx))                      
             resxPAM = resxPAM*fbad_cog(0,icx)  
   
          else  
             if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx  
          endif  
   
   
 *     ======================================  
 *     temporary patch for saturated clusters  
 *     ======================================  
          if( nsatstrips(icx).gt.0 )then  
             stripx  = stripx + cog(4,icx)              
             resxPAM = pitchX*1e-4/sqrt(12.)  
 cc            cc=cog(4,icx)  
 c$$$            print*,icx,' *** ',cc  
 c$$$            print*,icx,' *** ',resxPAM  
          endif  
759    
760   10   endif   10   endif
   
761            
762  *     -----------------  *     -----------------
763  *     CLUSTER Y  *     CLUSTER Y
# Line 876  c$$$            print*,icx,' *** ',resxP Line 768  c$$$            print*,icx,' *** ',resxP
768           viewy = VIEW(icy)           viewy = VIEW(icy)
769           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
770           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
771           resyPAM = RESYAV  c         resyPAM = RESYAV
772           stripy = float(MAXS(icy))           stripy = float(MAXS(icy))
773    
774           if(           if(
# Line 900  c$$$            print*,icx,' *** ',resxP Line 792  c$$$            print*,icx,' *** ',resxP
792              endif              endif
793              goto 100              goto 100
794           endif           endif
795    
796  *        --------------------------  *        --------------------------
797  *        magnetic-field corrections  *        magnetic-field corrections
798  *        --------------------------  *        --------------------------
799           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001                   stripy  = stripy +  fieldcorr(viewy,bfx)      
800           angy    = 180.*atan(tgtemp)/acos(-1.)           angy    = effectiveangle(ay,viewy,bfx)
          stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY  
 c$$$         if(bfx.ne.0.)print*,viewy,'-y- '  
 c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ  
 *        --------------------------  
801                    
802  c$$$         print*,'--- y-cl ---'           call applypfa(PFAy,icy,angy,corr,res)
803  c$$$         istart = INDSTART(ICY)           stripy  = stripy + corr
804  c$$$         istop  = TOTCLLENGTH           resyPAM = res
 c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1  
 c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)  
 c$$$         print*,(CLSIGNAL(i),i=istart,istop)  
 c$$$         print*,INDMAX(icy)  
 c$$$         print*,cog(4,icy)  
 c$$$         print*,fbad_cog(4,icy)  
   
          if(PFAy.eq.'COG1')then  
   
             stripy  = stripy      
             resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM    
   
          elseif(PFAy.eq.'COG2')then  
   
             stripy  = stripy + cog(2,icy)  
             resyPAM = risy_cog(abs(angy))!TEMPORANEO  
             resyPAM = resyPAM*fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'COG3')then  
   
             stripy  = stripy + cog(3,icy)  
             resyPAM = risy_cog(abs(angy))!TEMPORANEO  
             resyPAM = resyPAM*fbad_cog(3,icy)  
   
          elseif(PFAy.eq.'COG4')then  
   
             stripy  = stripy + cog(4,icy)  
             resyPAM = risy_cog(abs(angy))!TEMPORANEO  
             resyPAM = resyPAM*fbad_cog(4,icy)  
   
          elseif(PFAy.eq.'ETA2')then  
   
             stripy  = stripy + pfaeta2(icy,angy)            
             resyPAM = risyeta2(abs(angy))                
             resyPAM = resyPAM*fbad_cog(2,icy)  
 c            if(DEBUG.and.fbad_cog(2,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'ETA3')then                        
   
             stripy  = stripy + pfaeta3(icy,angy)  
             resyPAM = resyPAM*fbad_cog(3,icy)    
 c            if(DEBUG.and.fbad_cog(3,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)  
   
          elseif(PFAy.eq.'ETA4')then    
   
             stripy  = stripy + pfaeta4(icy,angy)  
             resyPAM = resyPAM*fbad_cog(4,icy)  
 c            if(DEBUG.and.fbad_cog(4,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)  
   
          elseif(PFAy.eq.'ETA')then  
   
             stripy  = stripy + pfaeta(icy,angy)  
 c            resyPAM = riseta(icy,angy)    
             resyPAM = riseta(viewy,angy)    
             resyPAM = resyPAM*fbad_eta(icy,angy)  
 c            if(DEBUG.and.fbad_cog(2,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'ETAL')then  
   
             stripy  = stripy + pfaetal(icy,angy)  
             resyPAM = riseta(viewy,angy)    
             resyPAM = resyPAM*fbad_eta(icy,angy)  
 c            if(DEBUG.and.fbad_cog(2,icy).ne.1)  
 c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
   
          elseif(PFAy.eq.'COG')then  
   
             stripy  = stripy + cog(0,icy)              
             resyPAM = risy_cog(abs(angy))  
             resyPAM = resyPAM*fbad_cog(0,icy)  
   
          else  
             if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx  
          endif  
   
   
 *     ======================================  
 *     temporary patch for saturated clusters  
 *     ======================================  
          if( nsatstrips(icy).gt.0 )then  
             stripy  = stripy + cog(4,icy)              
             resyPAM = pitchY*1e-4/sqrt(12.)  
 cc            cc=cog(4,icy)  
 c$$$            print*,icy,' *** ',cc  
 c$$$            print*,icy,' *** ',resyPAM  
          endif  
   
805    
806   20   endif   20   endif
807    
808  c$$$      print*,'## stripx,stripy ',stripx,stripy  cc      print*,'## stripx,stripy ',stripx,stripy
809    
810  c===========================================================  c===========================================================
811  C     COUPLE  C     COUPLE
# Line 1024  c--------------------------------------- Line 822  c---------------------------------------
822       $              ' WARNING: false X strip: strip ',stripx       $              ' WARNING: false X strip: strip ',stripx
823              endif              endif
824           endif           endif
825           xi = acoordsi(stripx,viewx)           xi = dcoordsi(stripx,viewx)
826           yi = acoordsi(stripy,viewy)           yi = dcoordsi(stripy,viewy)
827           zi = 0.           zi = 0.D0
828                    
   
829  c------------------------------------------------------------------------  c------------------------------------------------------------------------
830  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
831  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1062  c--------------------------------------- Line 859  c---------------------------------------
859           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
860           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
861    
862           xPAM_A = 0.           xPAM_A = 0.D0
863           yPAM_A = 0.           yPAM_A = 0.D0
864           zPAM_A = 0.           zPAM_A = 0.D0
865    
866           xPAM_B = 0.           xPAM_B = 0.D0
867           yPAM_B = 0.           yPAM_B = 0.D0
868           zPAM_B = 0.           zPAM_B = 0.D0
869    
870        elseif(        elseif(
871       $        (icx.ne.0.and.icy.eq.0).or.       $        (icx.ne.0.and.icy.eq.0).or.
# Line 1088  C======================================= Line 885  C=======================================
885              nldx = nldy              nldx = nldy
886              viewx = viewy + 1              viewx = viewy + 1
887    
888              yi   = acoordsi(stripy,viewy)              yi   = dcoordsi(stripy,viewy)
889    
890              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
891              yi_A = yi              yi_A = yi
# Line 1119  c            if((stripx.le.3).or.(stripx Line 916  c            if((stripx.le.3).or.(stripx
916       $                 ' WARNING: false X strip: strip ',stripx       $                 ' WARNING: false X strip: strip ',stripx
917                 endif                 endif
918              endif              endif
919              xi   = acoordsi(stripx,viewx)              xi   = dcoordsi(stripx,viewx)
920    
921              xi_A = xi              xi_A = xi
922              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 1192  c     (xPAM,yPAM,zPAM) = measured coordi Line 989  c     (xPAM,yPAM,zPAM) = measured coordi
989  c                        in PAMELA reference system  c                        in PAMELA reference system
990  c------------------------------------------------------------------------  c------------------------------------------------------------------------
991    
992           xPAM = 0.           xPAM = 0.D0
993           yPAM = 0.           yPAM = 0.D0
994           zPAM = 0.           zPAM = 0.D0
995    
996           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
997           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 1256  c$$$      PFAy = 'COG4'!PFA Line 1053  c$$$      PFAy = 'COG4'!PFA
1053    
1054        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then        if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1055              print*,'xyzpam: ***WARNING*** clusters ',icx,icy              print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1056       $           ,' does not exists (nclstr1=',nclstr1,')'       $           ,' do not exists (n.clusters=',nclstr1,')'
1057              icx = -1*icx              icx = -1*icx
1058              icy = -1*icy              icy = -1*icy
1059              return              return
# Line 1303  c$$$     $        ,' does not belong to Line 1100  c$$$     $        ,' does not belong to
1100           xm(ip) = xPAM           xm(ip) = xPAM
1101           ym(ip) = yPAM           ym(ip) = yPAM
1102           zm(ip) = zPAM           zm(ip) = zPAM
1103           xm_A(ip) = 0.           xm_A(ip) = 0.D0
1104           ym_A(ip) = 0.           ym_A(ip) = 0.D0
1105           xm_B(ip) = 0.           xm_B(ip) = 0.D0
1106           ym_B(ip) = 0.           ym_B(ip) = 0.D0
1107    
1108  c         zv(ip) = zPAM  c         zv(ip) = zPAM
1109    
# Line 1426  c$$$         print*,'------------------- Line 1223  c$$$         print*,'-------------------
1223  *      *    
1224  ********************************************************************************  ********************************************************************************
1225    
1226        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1227    
1228        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1229    
# Line 1435  c$$$         print*,'------------------- Line 1232  c$$$         print*,'-------------------
1232  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1233  *     -----------------------------------  *     -----------------------------------
1234    
1235          real rXPP,rYPP
1236          double precision XPP,YPP
1237        double precision distance,RE        double precision distance,RE
1238        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1239    
1240          XPP=DBLE(rXPP)
1241          YPP=DBLE(rYPP)
1242    
1243  *     ----------------------  *     ----------------------
1244        if (        if (
1245       +     xPAM.eq.0.and.       +     xPAM.eq.0.and.
# Line 1564  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1366  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1366        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1367        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1368        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1369        real*8 yvvv,xvvv        double precision yvvv,xvvv
1370        double precision xi,yi,zi        double precision xi,yi,zi
1371        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1372        real AA,BB        real AA,BB
1373        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1374    
1375  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1376        real ptoll        real ptoll
1377        data ptoll/150./          !um        data ptoll/150./          !um
1378    
1379        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1380    
1381        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1382        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1589  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1391  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1391  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1392  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1393  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1394                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)  c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1395       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...  c     $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1396  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...
1397  c                  print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1398  c     $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1399                 endif  c               endif
1400                 xi = acoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1401                 yi = acoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1402                 zi = 0.                 zi = 0.D0
1403  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1404  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1405  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 2243  c     $                               (i Line 2045  c     $                               (i
2045                                   ym3=yPAM                                   ym3=yPAM
2046                                   zm3=zPAM                                   zm3=zPAM
2047  *     find the circle passing through the three points  *     find the circle passing through the three points
2048    c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
2049    c$$$     $                                ,xc,zc,radius,iflag)
2050                                     iflag = DEBUG
2051                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2052       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag)
2053  c     print*,xc,zc,radius  c     print*,xc,zc,radius
# Line 3247  c$$$      endif Line 3052  c$$$      endif
3052        call track_init        call track_init
3053        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3054    
3055             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3056    
3057           xP=XV_STORE(nplanes-ip+1,ibest)           xP=XV_STORE(nplanes-ip+1,ibest)
3058           yP=YV_STORE(nplanes-ip+1,ibest)           yP=YV_STORE(nplanes-ip+1,ibest)
3059           zP=ZV_STORE(nplanes-ip+1,ibest)           zP=ZV_STORE(nplanes-ip+1,ibest)
# Line 3348  c            if(DEBUG.EQ.1)print*,'>>>> Line 3155  c            if(DEBUG.EQ.1)print*,'>>>>
3155                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3156  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
3157  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
3158       $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)       $              cl_used(icx).ne.0.or. !or the X cluster is already used
3159       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)       $              cl_used(icy).ne.0.or. !or the Y cluster is already used
3160       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3161  *            *          
3162                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3495  c                 dedxmm = sgnl(icy)  !( Line 3302  c                 dedxmm = sgnl(icy)  !(
3302  *----- single clusters -----------------------------------------------    *----- single clusters -----------------------------------------------  
3303  c            print*,'## ncls(',ip,') ',ncls(ip)  c            print*,'## ncls(',ip,') ',ncls(ip)
3304              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3305    c               print*,'-',ic,'-'
3306                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3307  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
3308                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
# Line 3729  c$$$               cl_used(icl)=ntrk   ! Line 3537  c$$$               cl_used(icl)=ntrk   !
3537              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3538              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3539              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3540                multmaxx(ip,it) = 0
3541                seedx(ip,it)    = 0
3542                xpu(ip,it)      = 0
3543                multmaxy(ip,it) = 0
3544                seedy(ip,it)    = 0
3545                ypu(ip,it)      = 0
3546           enddo           enddo
3547           do ipa=1,5           do ipa=1,5
3548              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3748  c$$$               cl_used(icl)=ntrk   ! Line 3562  c$$$               cl_used(icl)=ntrk   !
3562          ys(1,ip)=0          ys(1,ip)=0
3563          ys(2,ip)=0          ys(2,ip)=0
3564          sgnlys(ip)=0          sgnlys(ip)=0
3565            sxbad(ip)=0
3566            sybad(ip)=0
3567            multmaxsx(ip)=0
3568            multmaxsy(ip)=0
3569        enddo        enddo
3570        end        end
3571    
# Line 3921  c$$$               cl_used(icl)=ntrk   ! Line 3739  c$$$               cl_used(icl)=ntrk   !
3739           ay   = ayv_nt(ip,ntr)           ay   = ayv_nt(ip,ntr)
3740           bfx  = BX_STORE(ip,IDCAND)           bfx  = BX_STORE(ip,IDCAND)
3741           bfy  = BY_STORE(ip,IDCAND)           bfy  = BY_STORE(ip,IDCAND)
3742           if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)  c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3743           if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)  c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3744           tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001  c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3745           angx     = 180.*atan(tgtemp)/acos(-1.)  c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3746           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001          c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3747           angy    = 180.*atan(tgtemp)/acos(-1.)  c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3748    
3749             angx = effectiveangle(ax,2*ip,bfy)
3750             angy = effectiveangle(ay,2*ip-1,bfx)
3751            
3752                    
3753  c         print*,'* ',ip,bfx,bfy,angx,angy  c         print*,'* ',ip,bfx,bfy,angx,angy
3754    
# Line 3947  c           >>> is a couple Line 3769  c           >>> is a couple
3769              cl_used(cltrx(ip,ntr)) = 1     !tag used clusters                        cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3770              cl_used(cltry(ip,ntr)) = 1     !tag used clusters                        cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3771    
 c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)  
 c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)              
 c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))  
 c$$$            ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))  
3772              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3773              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3774    
# Line 3960  c$$$            ybad(ip,ntr)= nbadstrips Line 3778  c$$$            ybad(ip,ntr)= nbadstrips
3778              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3779       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3780    
3781                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3782         $                         +10000*mult(cltrx(ip,ntr))
3783                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3784         $           /clsigma(indmax(cltrx(ip,ntr)))
3785                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3786                xpu(ip,ntr)      = corr
3787    
3788                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3789         $                         +10000*mult(cltry(ip,ntr))
3790                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3791         $           /clsigma(indmax(cltry(ip,ntr)))
3792                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3793                ypu(ip,ntr)      = corr
3794    
3795           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3796    
3797    
3798              cl_used(icl) = 1    !tag used clusters                        cl_used(icl) = 1    !tag used clusters          
3799    
3800              if(mod(VIEW(icl),2).eq.0)then              if(mod(VIEW(icl),2).eq.0)then
3801                 cltrx(ip,ntr)=icl                 cltrx(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angx)  
 c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)  
3802                 xbad(ip,ntr) = nbadstrips(4,icl)                 xbad(ip,ntr) = nbadstrips(4,icl)
3803    
3804                 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)                 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3805    
3806                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3807         $                         +10000*mult(cltrx(ip,ntr))
3808                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3809         $           /clsigma(indmax(cltrx(ip,ntr)))
3810                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3811                   xpu(ip,ntr)      = corr
3812    
3813              elseif(mod(VIEW(icl),2).eq.1)then              elseif(mod(VIEW(icl),2).eq.1)then
3814                 cltry(ip,ntr)=icl                 cltry(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angy)  
 c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)  
3815                 ybad(ip,ntr) = nbadstrips(4,icl)                 ybad(ip,ntr) = nbadstrips(4,icl)
3816    
3817                 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)                 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3818    
3819                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3820         $                         +10000*mult(cltry(ip,ntr))
3821                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3822         $           /clsigma(indmax(cltry(ip,ntr)))
3823                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3824                   ypu(ip,ntr)      = corr
3825                  
3826              endif              endif
3827    
3828           endif                     endif          
# Line 4033  c         if( mask_view(iv).ne.0 )good2( Line 3879  c         if( mask_view(iv).ne.0 )good2(
3879           ip=nplanes-npl(VIEW(icl))+1                       ip=nplanes-npl(VIEW(icl))+1            
3880                    
3881           if(cl_used(icl).eq.0)then !cluster not included in any track           if(cl_used(icl).eq.0)then !cluster not included in any track
3882    
3883              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3884    
3885                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3886                 planex(nclsx) = ip                 planex(nclsx) = ip
3887                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3888                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3889                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3890                   sxbad(nclsx)  = nbadstrips(1,icl)
3891                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
3892                  
3893    cc               print*,icl,' >>>> ',sxbad(nclsx)
3894    
3895                 do is=1,2                 do is=1,2
3896  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3897  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
# Line 4056  c$$$               print*,'xs(2,nclsx)   Line 3909  c$$$               print*,'xs(2,nclsx)  
3909                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3910                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3911                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3912                   sybad(nclsy)  = nbadstrips(1,icl)
3913                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
3914    
3915    cc               print*,icl,' >>>> ',sybad(nclsy)
3916    
3917                 do is=1,2                 do is=1,2
3918  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3919  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
# Line 4079  c$$$               print*,'ys(2,nclsy)   Line 3937  c$$$               print*,'ys(2,nclsy)  
3937  *     associati ad una traccia, e permettere di salvare  *     associati ad una traccia, e permettere di salvare
3938  *     solo questi nell'albero di uscita  *     solo questi nell'albero di uscita
3939  *     --------------------------------------------------  *     --------------------------------------------------
3940                            
   
 c$$$         print*,' cl ',icl,' --> ',cl_used(icl)  
 c$$$  
 c$$$         if( cl_used(icl).ne.0 )then  
 c$$$            if(  
 c$$$     $           mod(VIEW(icl),2).eq.0.and.  
 c$$$     $           cltrx(ip,whichtrack(icl)).ne.icl )  
 c$$$     $           print*,'**WARNING** cltrx(',ip,',',whichtrack(icl)  
 c$$$     $           ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl  
 c$$$            if(  
 c$$$     $           mod(VIEW(icl),2).eq.1.and.  
 c$$$     $           cltry(ip,whichtrack(icl)).ne.icl )  
 c$$$     $           print*,'**WARNING** cltry(',ip,',',whichtrack(icl)  
 c$$$     $           ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl  
 c$$$         endif  
           
   
3941        enddo        enddo
3942        end        end
3943    
Legend:
Removed from v.1.29  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.23