/[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.24 by pam-fi, Tue May 15 16:22:18 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 678  c$$$         print*,fbad_cog(4,icx) Line 691  c$$$         print*,fbad_cog(4,icx)
691           elseif(PFAx.eq.'ETA')then             elseif(PFAx.eq.'ETA')then  
692    
693              stripx  = stripx + pfaeta(icx,angx)                          stripx  = stripx + pfaeta(icx,angx)            
694              resxPAM = ris_eta(icx,angx)                      c            resxPAM = riseta(icx,angx)                    
695                resxPAM = riseta(viewx,angx)                    
696              resxPAM = resxPAM*fbad_eta(icx,angx)                          resxPAM = resxPAM*fbad_eta(icx,angx)            
697              if(DEBUG.and.fbad_cog(2,icx).ne.1)                            if(DEBUG.and.fbad_cog(2,icx).ne.1)              
698       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
# Line 690  c$$$         print*,fbad_cog(4,icx) Line 704  c$$$         print*,fbad_cog(4,icx)
704              resxPAM = resxPAM*fbad_cog(0,icx)              resxPAM = resxPAM*fbad_cog(0,icx)
705    
706           else           else
707              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
708           endif           endif
709    
710    
# Line 720  c$$$            print*,icx,' *** ',resxP Line 734  c$$$            print*,icx,' *** ',resxP
734           stripy = float(MAXS(icy))           stripy = float(MAXS(icy))
735    
736           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
737              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG) then
738       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
739         $              ,icx,icy
740                endif
741              goto 100              goto 100
742           endif           endif
743  *        --------------------------  *        --------------------------
# Line 730  c$$$            print*,icx,' *** ',resxP Line 746  c$$$            print*,icx,' *** ',resxP
746           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001                   tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
747           angy    = 180.*atan(tgtemp)/acos(-1.)           angy    = 180.*atan(tgtemp)/acos(-1.)
748           stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY           stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
749    c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
750    c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
751  *        --------------------------  *        --------------------------
752                    
753  c$$$         print*,'--- y-cl ---'  c$$$         print*,'--- y-cl ---'
# Line 790  c$$$         print*,fbad_cog(4,icy) Line 808  c$$$         print*,fbad_cog(4,icy)
808           elseif(PFAy.eq.'ETA')then           elseif(PFAy.eq.'ETA')then
809    
810              stripy  = stripy + pfaeta(icy,angy)              stripy  = stripy + pfaeta(icy,angy)
811              resyPAM = ris_eta(icy,angy)    c            resyPAM = riseta(icy,angy)  
812                resyPAM = riseta(viewy,angy)  
813              resyPAM = resyPAM*fbad_eta(icy,angy)              resyPAM = resyPAM*fbad_eta(icy,angy)
814              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
815       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
# Line 802  c$$$         print*,fbad_cog(4,icy) Line 821  c$$$         print*,fbad_cog(4,icy)
821              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
822    
823           else           else
824              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
825           endif           endif
826    
827    
# Line 820  c$$$            print*,icy,' *** ',resyP Line 839  c$$$            print*,icy,' *** ',resyP
839    
840        endif        endif
841    
842  c      print*,'## stripx,stripy ',stripx,stripy  c$$$      print*,'## stripx,stripy ',stripx,stripy
843    
844  c===========================================================  c===========================================================
845  C     COUPLE  C     COUPLE
# Line 832  c     (xi,yi,zi) = mechanical coordinate Line 851  c     (xi,yi,zi) = mechanical coordinate
851  c------------------------------------------------------------------------  c------------------------------------------------------------------------
852           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
853       $        .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...
854              print*,'xyz_PAM (couple):',              if(DEBUG) then
855       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
856         $              ' WARNING: false X strip: strip ',stripx
857                endif
858           endif           endif
859           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
860           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
# Line 925  c            print*,'X-singlet ',icx,npl Line 946  c            print*,'X-singlet ',icx,npl
946  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...
947              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
948       $           .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...
949                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG) then
950       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
951         $                 ' WARNING: false X strip: strip ',stripx
952                   endif
953              endif              endif
954              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
955    
# Line 948  c            print*,'X-cl ',icx,stripx,' Line 971  c            print*,'X-cl ',icx,stripx,'
971  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
972    
973           else           else
974                if(DEBUG) then
975              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
976              print *,'icx = ',icx                 print *,'icx = ',icx
977              print *,'icy = ',icy                 print *,'icy = ',icy
978                endif
979              goto 100              goto 100
980                            
981           endif           endif
# Line 1016  c--------------------------------------- Line 1040  c---------------------------------------
1040  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1041    
1042        else        else
1043                       if(DEBUG) then
1044           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1045           print *,'icx = ',icx              print *,'icx = ',icx
1046           print *,'icy = ',icy              print *,'icy = ',icy
1047                         endif
1048        endif        endif
1049                    
1050    
# Line 1031  c      print*,'## xPAM_B,yPAM_B,zPAM_B ' Line 1055  c      print*,'## xPAM_B,yPAM_B,zPAM_B '
1055   100  continue   100  continue
1056        end        end
1057    
1058    ************************************************************************
1059    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1060    *     (done to be called from c/c++)
1061    ************************************************************************
1062    
1063          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1064    
1065          include 'commontracker.f'
1066          include 'level1.f'
1067          include 'common_mini_2.f'
1068          include 'common_xyzPAM.f'
1069          include 'common_mech.f'
1070          include 'calib.f'
1071          
1072    *     flag to chose PFA
1073    c$$$      character*10 PFA
1074    c$$$      common/FINALPFA/PFA
1075    
1076          integer icx,icy           !X-Y cluster ID
1077          integer sensor
1078          character*4 PFAx,PFAy     !PFA to be used
1079          real ax,ay                !X-Y geometric angle
1080          real bfx,bfy              !X-Y b-field components
1081    
1082          ipx=0
1083          ipy=0      
1084          
1085    c$$$      PFAx = 'COG4'!PFA
1086    c$$$      PFAy = 'COG4'!PFA
1087          
1088          call idtoc(pfaid,PFAx)
1089          call idtoc(pfaid,PFAy)
1090    
1091          call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1092    
1093    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1094          
1095          if(icx.ne.0.and.icy.ne.0)then
1096    
1097             ipx=npl(VIEW(icx))
1098             ipy=npl(VIEW(icy))
1099             if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1100         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1101         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1102    
1103             xgood(ip) = 1.
1104             ygood(ip) = 1.
1105             resx(ip)  = resxPAM
1106             resy(ip)  = resyPAM
1107    
1108             xm(ip) = xPAM
1109             ym(ip) = yPAM
1110             zm(ip) = zPAM
1111             xm_A(ip) = 0.
1112             ym_A(ip) = 0.
1113             xm_B(ip) = 0.
1114             ym_B(ip) = 0.
1115    
1116    c         zv(ip) = zPAM
1117    
1118          elseif(icx.eq.0.and.icy.ne.0)then
1119    
1120             ipy=npl(VIEW(icy))
1121             if((nplanes-ipy+1).ne.ip)
1122         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1123         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1124    
1125             xgood(ip) = 0.
1126             ygood(ip) = 1.
1127             resx(ip)  = 1000.
1128             resy(ip)  = resyPAM
1129    
1130             xm(ip) = -100.
1131             ym(ip) = -100.
1132             zm(ip) = (zPAM_A+zPAM_B)/2.
1133             xm_A(ip) = xPAM_A
1134             ym_A(ip) = yPAM_A
1135             xm_B(ip) = xPAM_B
1136             ym_B(ip) = yPAM_B
1137    
1138    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1139            
1140          elseif(icx.ne.0.and.icy.eq.0)then
1141    
1142             ipx=npl(VIEW(icx))
1143             if((nplanes-ipx+1).ne.ip)
1144         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1145         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1146    
1147             xgood(ip) = 1.
1148             ygood(ip) = 0.
1149             resx(ip)  = resxPAM
1150             resy(ip)  = 1000.
1151    
1152             xm(ip) = -100.
1153             ym(ip) = -100.
1154             zm(ip) = (zPAM_A+zPAM_B)/2.
1155             xm_A(ip) = xPAM_A
1156             ym_A(ip) = yPAM_A
1157             xm_B(ip) = xPAM_B
1158             ym_B(ip) = yPAM_B
1159            
1160    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1161    
1162          else
1163    
1164             il = 2
1165             if(lad.ne.0)il=lad
1166             is = 1
1167             if(sensor.ne.0)is=sensor
1168    c         print*,nplanes-ip+1,il,is
1169    
1170             xgood(ip) = 0.
1171             ygood(ip) = 0.
1172             resx(ip)  = 1000.
1173             resy(ip)  = 1000.
1174    
1175             xm(ip) = -100.
1176             ym(ip) = -100.          
1177             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1178             xm_A(ip) = 0.
1179             ym_A(ip) = 0.
1180             xm_B(ip) = 0.
1181             ym_B(ip) = 0.
1182    
1183    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1184    
1185          endif
1186    
1187          if(DEBUG)then
1188    c         print*,'----------------------------- track coord'
1189    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1190             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1191         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1192         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1193    c$$$         print*,'-----------------------------'
1194          endif
1195          end
1196    
1197  ********************************************************************************  ********************************************************************************
1198  ********************************************************************************  ********************************************************************************
# Line 1146  c$$$         print*,' resolution ',resxP Line 1308  c$$$         print*,' resolution ',resxP
1308                    
1309        else        else
1310                    
1311           print*  c         print*
1312       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1313           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1314           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 '
1315       $        ,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
1316        endif          endif  
1317    
1318        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1218  c--------------------------------------- Line 1380  c---------------------------------------
1380                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1381       $              .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...
1382  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...
1383                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1384       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1385                 endif                 endif
1386                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1387                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1374  c      include 'common_analysis.f' Line 1536  c      include 'common_analysis.f'
1536        is_cp=0        is_cp=0
1537        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1538        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1539        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1540    
1541        return        return
1542        end        end
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.23