/[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.20 by pam-fi, Fri Apr 27 10:39:58 2007 UTC revision 1.23 by pam-fi, Mon May 14 11:03:06 2007 UTC
# Line 331  c$$$         enddo Line 331  c$$$         enddo
331              iimage=0              iimage=0
332           endif           endif
333           if(icand.eq.0)then           if(icand.eq.0)then
334              print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand              if(VERBOSE)then
335       $           ,ibest,iimage                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
336         $              ,ibest,iimage
337                endif
338              return              return
339           endif           endif
340    
# Line 576  c     $        rchi2best.lt.15..and. Line 578  c     $        rchi2best.lt.15..and.
578                
579    
580        parameter (ndivx=30)        parameter (ndivx=30)
581    
582    
583    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
584                
585        resxPAM = 0        resxPAM = 0
586        resyPAM = 0        resyPAM = 0
# Line 607  c      print*,'## xyz_PAM: ',icx,icy,sen Line 612  c      print*,'## xyz_PAM: ',icx,icy,sen
612  *        --------------------------  *        --------------------------
613           angtemp  = ax           angtemp  = ax
614           bfytemp  = bfy           bfytemp  = bfy
615           if(nplx.eq.6) angtemp = -1. * ax  *        /////////////////////////////////
616           if(nplx.eq.6) bfytemp = -1. * bfy  *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!
617    *        *grvzkkjsdgjhhhgngbn###>:(
618    *        /////////////////////////////////
619    c         if(nplx.eq.6) angtemp = -1. * ax
620    c         if(nplx.eq.6) bfytemp = -1. * bfy
621             if(viewx.eq.12) angtemp = -1. * ax
622             if(viewx.eq.12) bfytemp = -1. * bfy
623           tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001           tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
624           angx     = 180.*atan(tgtemp)/acos(-1.)           angx     = 180.*atan(tgtemp)/acos(-1.)
625           stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX           stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
626  c$$$         print*,nplx,ax,bfy/10.  c$$$         print*,nplx,ax,bfy/10.
627  c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX  c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
628  c$$$         print*,'========================'  c$$$         print*,'========================'
629    c$$$         if(bfy.ne.0.)print*,viewx,'-x- '
630    c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ
631  *        --------------------------  *        --------------------------
632    
633  c$$$         print*,'--- x-cl ---'  c$$$         print*,'--- x-cl ---'
# Line 690  c$$$         print*,fbad_cog(4,icx) Line 703  c$$$         print*,fbad_cog(4,icx)
703              resxPAM = resxPAM*fbad_cog(0,icx)              resxPAM = resxPAM*fbad_cog(0,icx)
704    
705           else           else
706              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
707           endif           endif
708    
709    
# Line 720  c$$$            print*,icx,' *** ',resxP Line 733  c$$$            print*,icx,' *** ',resxP
733           stripy = float(MAXS(icy))           stripy = float(MAXS(icy))
734    
735           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
736              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG) then
737       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
738         $              ,icx,icy
739                endif
740              goto 100              goto 100
741           endif           endif
742  *        --------------------------  *        --------------------------
# Line 730  c$$$            print*,icx,' *** ',resxP Line 745  c$$$            print*,icx,' *** ',resxP
745           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001                   tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
746           angy    = 180.*atan(tgtemp)/acos(-1.)           angy    = 180.*atan(tgtemp)/acos(-1.)
747           stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY           stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
748    c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
749    c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
750  *        --------------------------  *        --------------------------
751                    
752  c$$$         print*,'--- y-cl ---'  c$$$         print*,'--- y-cl ---'
# Line 802  c$$$         print*,fbad_cog(4,icy) Line 819  c$$$         print*,fbad_cog(4,icy)
819              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
820    
821           else           else
822              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
823           endif           endif
824    
825    
# Line 820  c$$$            print*,icy,' *** ',resyP Line 837  c$$$            print*,icy,' *** ',resyP
837    
838        endif        endif
839    
840  c      print*,'## stripx,stripy ',stripx,stripy  c$$$      print*,'## stripx,stripy ',stripx,stripy
841    
842  c===========================================================  c===========================================================
843  C     COUPLE  C     COUPLE
# Line 832  c     (xi,yi,zi) = mechanical coordinate Line 849  c     (xi,yi,zi) = mechanical coordinate
849  c------------------------------------------------------------------------  c------------------------------------------------------------------------
850           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
851       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
852              print*,'xyz_PAM (couple):',              if(DEBUG) then
853       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
854         $              ' WARNING: false X strip: strip ',stripx
855                endif
856           endif           endif
857           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
858           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
# Line 925  c            print*,'X-singlet ',icx,npl Line 944  c            print*,'X-singlet ',icx,npl
944  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...
945              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
946       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
947                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG) then
948       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
949         $                 ' WARNING: false X strip: strip ',stripx
950                   endif
951              endif              endif
952              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
953    
# Line 948  c            print*,'X-cl ',icx,stripx,' Line 969  c            print*,'X-cl ',icx,stripx,'
969  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
970    
971           else           else
972                if(DEBUG) then
973              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
974              print *,'icx = ',icx                 print *,'icx = ',icx
975              print *,'icy = ',icy                 print *,'icy = ',icy
976                endif
977              goto 100              goto 100
978                            
979           endif           endif
# Line 1016  c--------------------------------------- Line 1038  c---------------------------------------
1038  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1039    
1040        else        else
1041                       if(DEBUG) then
1042           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1043           print *,'icx = ',icx              print *,'icx = ',icx
1044           print *,'icy = ',icy              print *,'icy = ',icy
1045                         endif
1046        endif        endif
1047                    
1048    
# Line 1031  c      print*,'## xPAM_B,yPAM_B,zPAM_B ' Line 1053  c      print*,'## xPAM_B,yPAM_B,zPAM_B '
1053   100  continue   100  continue
1054        end        end
1055    
1056    ************************************************************************
1057    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1058    *     (done to be called from c/c++)
1059    ************************************************************************
1060    
1061          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1062    
1063          include 'commontracker.f'
1064          include 'level1.f'
1065          include 'common_mini_2.f'
1066          include 'common_xyzPAM.f'
1067          include 'common_mech.f'
1068          include 'calib.f'
1069          
1070    *     flag to chose PFA
1071    c$$$      character*10 PFA
1072    c$$$      common/FINALPFA/PFA
1073    
1074          integer icx,icy           !X-Y cluster ID
1075          integer sensor
1076          character*4 PFAx,PFAy     !PFA to be used
1077          real ax,ay                !X-Y geometric angle
1078          real bfx,bfy              !X-Y b-field components
1079    
1080          ipx=0
1081          ipy=0      
1082          
1083    c$$$      PFAx = 'COG4'!PFA
1084    c$$$      PFAy = 'COG4'!PFA
1085          
1086          call idtoc(pfaid,PFAx)
1087          call idtoc(pfaid,PFAy)
1088    
1089          call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1090    
1091    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1092          
1093          if(icx.ne.0.and.icy.ne.0)then
1094    
1095             ipx=npl(VIEW(icx))
1096             ipy=npl(VIEW(icy))
1097             if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1098         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1099         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1100    
1101             xgood(ip) = 1.
1102             ygood(ip) = 1.
1103             resx(ip)  = resxPAM
1104             resy(ip)  = resyPAM
1105    
1106             xm(ip) = xPAM
1107             ym(ip) = yPAM
1108             zm(ip) = zPAM
1109             xm_A(ip) = 0.
1110             ym_A(ip) = 0.
1111             xm_B(ip) = 0.
1112             ym_B(ip) = 0.
1113    
1114    c         zv(ip) = zPAM
1115    
1116          elseif(icx.eq.0.and.icy.ne.0)then
1117    
1118             ipy=npl(VIEW(icy))
1119             if((nplanes-ipy+1).ne.ip)
1120         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1121         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1122    
1123             xgood(ip) = 0.
1124             ygood(ip) = 1.
1125             resx(ip)  = 1000.
1126             resy(ip)  = resyPAM
1127    
1128             xm(ip) = -100.
1129             ym(ip) = -100.
1130             zm(ip) = (zPAM_A+zPAM_B)/2.
1131             xm_A(ip) = xPAM_A
1132             ym_A(ip) = yPAM_A
1133             xm_B(ip) = xPAM_B
1134             ym_B(ip) = yPAM_B
1135    
1136    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1137            
1138          elseif(icx.ne.0.and.icy.eq.0)then
1139    
1140             ipx=npl(VIEW(icx))
1141             if((nplanes-ipx+1).ne.ip)
1142         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1143         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1144    
1145             xgood(ip) = 1.
1146             ygood(ip) = 0.
1147             resx(ip)  = resxPAM
1148             resy(ip)  = 1000.
1149    
1150             xm(ip) = -100.
1151             ym(ip) = -100.
1152             zm(ip) = (zPAM_A+zPAM_B)/2.
1153             xm_A(ip) = xPAM_A
1154             ym_A(ip) = yPAM_A
1155             xm_B(ip) = xPAM_B
1156             ym_B(ip) = yPAM_B
1157            
1158    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1159    
1160          else
1161    
1162             il = 2
1163             if(lad.ne.0)il=lad
1164             is = 1
1165             if(sensor.ne.0)is=sensor
1166    c         print*,nplanes-ip+1,il,is
1167    
1168             xgood(ip) = 0.
1169             ygood(ip) = 0.
1170             resx(ip)  = 1000.
1171             resy(ip)  = 1000.
1172    
1173             xm(ip) = -100.
1174             ym(ip) = -100.          
1175             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1176             xm_A(ip) = 0.
1177             ym_A(ip) = 0.
1178             xm_B(ip) = 0.
1179             ym_B(ip) = 0.
1180    
1181    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1182    
1183          endif
1184    
1185          if(DEBUG)then
1186    c         print*,'----------------------------- track coord'
1187    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1188             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1189         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1190         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1191    c$$$         print*,'-----------------------------'
1192          endif
1193          end
1194    
1195  ********************************************************************************  ********************************************************************************
1196  ********************************************************************************  ********************************************************************************
# Line 1146  c$$$         print*,' resolution ',resxP Line 1306  c$$$         print*,' resolution ',resxP
1306                    
1307        else        else
1308                    
1309           print*  c         print*
1310       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1311           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1312           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 '
1313       $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
1314        endif          endif  
1315    
1316        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1218  c--------------------------------------- Line 1378  c---------------------------------------
1378                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1379       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1380  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...
1381                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1382       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1383                 endif                 endif
1384                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1385                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1374  c      include 'common_analysis.f' Line 1534  c      include 'common_analysis.f'
1534        is_cp=0        is_cp=0
1535        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1536        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1537        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1538    
1539        return        return
1540        end        end
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23