/[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.36 by pam-fi, Tue Nov 25 14:41:38 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 1330  c$$$     $        ,' does not belong to Line 1127  c$$$     $        ,' does not belong to
1127           resx(ip)  = 1000.           resx(ip)  = 1000.
1128           resy(ip)  = resyPAM           resy(ip)  = resyPAM
1129    
1130    cPP --- old ---
1131    c$$$         xm(ip) = -100.
1132    c$$$         ym(ip) = -100.
1133    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1134    c$$$         xm_A(ip) = xPAM_A
1135    c$$$         ym_A(ip) = yPAM_A
1136    c$$$         xm_B(ip) = xPAM_B
1137    c$$$         ym_B(ip) = yPAM_B
1138    cPP --- new ---
1139           xm(ip) = -100.           xm(ip) = -100.
1140           ym(ip) = -100.           ym(ip) = -100.
1141           zm(ip) = (zPAM_A+zPAM_B)/2.           zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1142           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1143           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1144             zm_A(ip) = zPAM_A
1145           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
1146           ym_B(ip) = yPAM_B           ym_B(ip) = yPAM_B
1147             zm_B(ip) = zPAM_B
1148    cPP -----------
1149    
1150  c         zv(ip) = (zPAM_A+zPAM_B)/2.  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1151                    
# Line 1361  c$$$     $        ,' does not belong to Line 1170  c$$$     $        ,' does not belong to
1170           resx(ip)  = resxPAM           resx(ip)  = resxPAM
1171           resy(ip)  = 1000.           resy(ip)  = 1000.
1172    
1173    cPP --- old ---
1174    c$$$         xm(ip) = -100.
1175    c$$$         ym(ip) = -100.
1176    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1177    c$$$         xm_A(ip) = xPAM_A
1178    c$$$         ym_A(ip) = yPAM_A
1179    c$$$         xm_B(ip) = xPAM_B
1180    c$$$         ym_B(ip) = yPAM_B
1181    cPP --- new ---
1182           xm(ip) = -100.           xm(ip) = -100.
1183           ym(ip) = -100.           ym(ip) = -100.
1184           zm(ip) = (zPAM_A+zPAM_B)/2.           zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1185           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1186           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1187             zm_A(ip) = zPAM_A
1188           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
1189           ym_B(ip) = yPAM_B           ym_B(ip) = yPAM_B
1190                     zm_B(ip) = zPAM_B
1191    cPP -----------
1192    
1193  c         zv(ip) = (zPAM_A+zPAM_B)/2.  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1194    
1195        else        else
# Line 1426  c$$$         print*,'------------------- Line 1247  c$$$         print*,'-------------------
1247  *      *    
1248  ********************************************************************************  ********************************************************************************
1249    
1250        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1251    
1252        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1253    
# Line 1435  c$$$         print*,'------------------- Line 1256  c$$$         print*,'-------------------
1256  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1257  *     -----------------------------------  *     -----------------------------------
1258    
1259          real rXPP,rYPP
1260          double precision XPP,YPP
1261        double precision distance,RE        double precision distance,RE
1262        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1263    
1264          XPP=DBLE(rXPP)
1265          YPP=DBLE(rYPP)
1266    
1267  *     ----------------------  *     ----------------------
1268        if (        if (
1269       +     xPAM.eq.0.and.       +     xPAM.eq.0.and.
# Line 1564  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1390  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1390        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1391        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1392        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1393        real*8 yvvv,xvvv        double precision yvvv,xvvv
1394        double precision xi,yi,zi        double precision xi,yi,zi
1395        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1396        real AA,BB        real AA,BB
1397        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1398    
1399  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1400        real ptoll        real ptoll
1401        data ptoll/150./          !um        data ptoll/150./          !um
1402    
1403        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1404    
1405        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1406        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1589  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA Line 1415  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPA
1415  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1416  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1417  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1418                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)  c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1419       $              .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...
1420  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...
1421  c                  print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1422  c     $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1423                 endif  c               endif
1424                 xi = acoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1425                 yi = acoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1426                 zi = 0.                 zi = 0.D0
1427  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1428  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1429  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 2243  c     $                               (i Line 2069  c     $                               (i
2069                                   ym3=yPAM                                   ym3=yPAM
2070                                   zm3=zPAM                                   zm3=zPAM
2071  *     find the circle passing through the three points  *     find the circle passing through the three points
2072    c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
2073    c$$$     $                                ,xc,zc,radius,iflag)
2074                                     iflag = DEBUG
2075                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2076       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag)
2077  c     print*,xc,zc,radius  c     print*,xc,zc,radius
# Line 3247  c$$$      endif Line 3076  c$$$      endif
3076        call track_init        call track_init
3077        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3078    
3079             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3080    
3081           xP=XV_STORE(nplanes-ip+1,ibest)           xP=XV_STORE(nplanes-ip+1,ibest)
3082           yP=YV_STORE(nplanes-ip+1,ibest)           yP=YV_STORE(nplanes-ip+1,ibest)
3083           zP=ZV_STORE(nplanes-ip+1,ibest)           zP=ZV_STORE(nplanes-ip+1,ibest)
# Line 3348  c            if(DEBUG.EQ.1)print*,'>>>> Line 3179  c            if(DEBUG.EQ.1)print*,'>>>>
3179                 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
3180  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
3181  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
3182       $              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
3183       $              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
3184       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3185  *            *          
3186                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3495  c                 dedxmm = sgnl(icy)  !( Line 3326  c                 dedxmm = sgnl(icy)  !(
3326  *----- single clusters -----------------------------------------------    *----- single clusters -----------------------------------------------  
3327  c            print*,'## ncls(',ip,') ',ncls(ip)  c            print*,'## ncls(',ip,') ',ncls(ip)
3328              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3329    c               print*,'-',ic,'-'
3330                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3331  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
3332                 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 3583  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~ Line 3415  c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~
3415  *     ----------------------------  *     ----------------------------
3416                    xm_A(nplanes-ip+1) = xmm_A                    xm_A(nplanes-ip+1) = xmm_A
3417                    ym_A(nplanes-ip+1) = ymm_A                    ym_A(nplanes-ip+1) = ymm_A
3418                      zm_A(nplanes-ip+1) = zmm_A
3419                    xm_B(nplanes-ip+1) = xmm_B                    xm_B(nplanes-ip+1) = xmm_B
3420                    ym_B(nplanes-ip+1) = ymm_B                    ym_B(nplanes-ip+1) = ymm_B
3421                      zm_B(nplanes-ip+1) = zmm_B
3422                    zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.                    zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3423    c$$$                  zm(nplanes-ip+1) =
3424    c$$$     $                 z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
3425                    dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<                    dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3426                    dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<                    dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3427  *     ----------------------------  *     ----------------------------
# Line 3729  c$$$               cl_used(icl)=ntrk   ! Line 3565  c$$$               cl_used(icl)=ntrk   !
3565              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3566              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3567              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3568                multmaxx(ip,it) = 0
3569                seedx(ip,it)    = 0
3570                xpu(ip,it)      = 0
3571                multmaxy(ip,it) = 0
3572                seedy(ip,it)    = 0
3573                ypu(ip,it)      = 0
3574           enddo           enddo
3575           do ipa=1,5           do ipa=1,5
3576              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3748  c$$$               cl_used(icl)=ntrk   ! Line 3590  c$$$               cl_used(icl)=ntrk   !
3590          ys(1,ip)=0          ys(1,ip)=0
3591          ys(2,ip)=0          ys(2,ip)=0
3592          sgnlys(ip)=0          sgnlys(ip)=0
3593            sxbad(ip)=0
3594            sybad(ip)=0
3595            multmaxsx(ip)=0
3596            multmaxsy(ip)=0
3597        enddo        enddo
3598        end        end
3599    
# Line 3921  c$$$               cl_used(icl)=ntrk   ! Line 3767  c$$$               cl_used(icl)=ntrk   !
3767           ay   = ayv_nt(ip,ntr)           ay   = ayv_nt(ip,ntr)
3768           bfx  = BX_STORE(ip,IDCAND)           bfx  = BX_STORE(ip,IDCAND)
3769           bfy  = BY_STORE(ip,IDCAND)           bfy  = BY_STORE(ip,IDCAND)
3770           if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)  c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3771           if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)  c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3772           tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001  c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3773           angx     = 180.*atan(tgtemp)/acos(-1.)  c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3774           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001          c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3775           angy    = 180.*atan(tgtemp)/acos(-1.)  c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3776    
3777             angx = effectiveangle(ax,2*ip,bfy)
3778             angy = effectiveangle(ay,2*ip-1,bfx)
3779            
3780                    
3781  c         print*,'* ',ip,bfx,bfy,angx,angy  c         print*,'* ',ip,bfx,bfy,angx,angy
3782    
# Line 3947  c           >>> is a couple Line 3797  c           >>> is a couple
3797              cl_used(cltrx(ip,ntr)) = 1     !tag used clusters                        cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3798              cl_used(cltry(ip,ntr)) = 1     !tag used clusters                        cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3799    
 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)))  
3800              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))              xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3801              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))              ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3802    
# Line 3960  c$$$            ybad(ip,ntr)= nbadstrips Line 3806  c$$$            ybad(ip,ntr)= nbadstrips
3806              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)              if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3807       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)       $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3808    
3809                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3810         $                         +10000*mult(cltrx(ip,ntr))
3811                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3812         $           /clsigma(indmax(cltrx(ip,ntr)))
3813                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3814                xpu(ip,ntr)      = corr
3815    
3816                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3817         $                         +10000*mult(cltry(ip,ntr))
3818                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3819         $           /clsigma(indmax(cltry(ip,ntr)))
3820                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3821                ypu(ip,ntr)      = corr
3822    
3823           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3824    
3825    
3826              cl_used(icl) = 1    !tag used clusters                        cl_used(icl) = 1    !tag used clusters          
3827    
3828              if(mod(VIEW(icl),2).eq.0)then              if(mod(VIEW(icl),2).eq.0)then
3829                 cltrx(ip,ntr)=icl                 cltrx(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angx)  
 c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)  
3830                 xbad(ip,ntr) = nbadstrips(4,icl)                 xbad(ip,ntr) = nbadstrips(4,icl)
3831    
3832                 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)
3833    
3834                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3835         $                         +10000*mult(cltrx(ip,ntr))
3836                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3837         $           /clsigma(indmax(cltrx(ip,ntr)))
3838                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3839                   xpu(ip,ntr)      = corr
3840    
3841              elseif(mod(VIEW(icl),2).eq.1)then              elseif(mod(VIEW(icl),2).eq.1)then
3842                 cltry(ip,ntr)=icl                 cltry(ip,ntr)=icl
 c$$$               nnnnn = npfastrips(icl,PFA,angy)  
 c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)  
3843                 ybad(ip,ntr) = nbadstrips(4,icl)                 ybad(ip,ntr) = nbadstrips(4,icl)
3844    
3845                 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)
3846    
3847                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3848         $                         +10000*mult(cltry(ip,ntr))
3849                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3850         $           /clsigma(indmax(cltry(ip,ntr)))
3851                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3852                   ypu(ip,ntr)      = corr
3853                  
3854              endif              endif
3855    
3856           endif                     endif          
# Line 4033  c         if( mask_view(iv).ne.0 )good2( Line 3907  c         if( mask_view(iv).ne.0 )good2(
3907           ip=nplanes-npl(VIEW(icl))+1                       ip=nplanes-npl(VIEW(icl))+1            
3908                    
3909           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
3910    
3911    cc            print*,' ic ',icl,' --- stored '
3912    
3913              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3914    
3915                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3916                 planex(nclsx) = ip                 planex(nclsx) = ip
3917                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3918                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)                 if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3919                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3920                   sxbad(nclsx)  = nbadstrips(1,icl)
3921                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
3922                  
3923    cc               print*,icl,' >>>> ',sxbad(nclsx)
3924    
3925                 do is=1,2                 do is=1,2
3926  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3927  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 3939  c$$$               print*,'xs(2,nclsx)  
3939                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3940                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)                 if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3941                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3942                   sybad(nclsy)  = nbadstrips(1,icl)
3943                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
3944    
3945    cc               print*,icl,' >>>> ',sybad(nclsy)
3946    
3947                 do is=1,2                 do is=1,2
3948  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3949  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 3967  c$$$               print*,'ys(2,nclsy)  
3967  *     associati ad una traccia, e permettere di salvare  *     associati ad una traccia, e permettere di salvare
3968  *     solo questi nell'albero di uscita  *     solo questi nell'albero di uscita
3969  *     --------------------------------------------------  *     --------------------------------------------------
3970                            
   
 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  
           
   
3971        enddo        enddo
3972        end        end
3973    
Legend:
Removed from v.1.29  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.23