/[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.44 by pamelats, Wed Dec 30 14:22:23 2009 UTC revision 1.47 by pam-ts, Wed Jun 4 07:57:04 2014 UTC
# Line 1046  c$$$      PFAy = 'COG4'!PFA Line 1046  c$$$      PFAy = 'COG4'!PFA
1046           ipy=npl(VIEW(icy))           ipy=npl(VIEW(icy))
1047                    
1048           if( (nplanes-ipx+1).ne.ip )then           if( (nplanes-ipx+1).ne.ip )then
1049              print*,'xyzpam: ***WARNING*** cluster ',icx              print*,'xyzpam: ***WARNING*** cluster icx=',icx
1050       $           ,' does not belong to plane: ',ip       $           ,' belongs to plane ',(nplanes-ipx+1)            
1051         $           ,' and not ',ip
1052              icx = -1*icx              icx = -1*icx
1053              return              return
1054           endif           endif
1055           if( (nplanes-ipy+1).ne.ip )then           if( (nplanes-ipy+1).ne.ip )then
1056              print*,'xyzpam: ***WARNING*** cluster ',icy              print*,'xyzpam: ***WARNING*** cluster icy=',icy
1057       $           ,' does not belong to plane: ',ip       $           ,' belongs to plane ',(nplanes-ipy+1)            
1058              icy = -1*icy       $           ,' and not ',ip
1059                 icy = -1*icy
1060              return              return
1061           endif           endif
1062    
# Line 1070  c$$$      PFAy = 'COG4'!PFA Line 1072  c$$$      PFAy = 'COG4'!PFA
1072           zm(ip) = zPAM           zm(ip) = zPAM
1073           xm_A(ip) = 0.D0           xm_A(ip) = 0.D0
1074           ym_A(ip) = 0.D0           ym_A(ip) = 0.D0
1075             zm_A(ip) = 0.D0
1076           xm_B(ip) = 0.D0           xm_B(ip) = 0.D0
1077           ym_B(ip) = 0.D0           ym_B(ip) = 0.D0
1078             zm_B(ip) = 0.D0
1079    
1080  c         zv(ip) = zPAM  c         zv(ip) = zPAM
1081    
# Line 1079  c         zv(ip) = zPAM Line 1083  c         zv(ip) = zPAM
1083    
1084           ipy=npl(VIEW(icy))           ipy=npl(VIEW(icy))
1085           if( (nplanes-ipy+1).ne.ip )then           if( (nplanes-ipy+1).ne.ip )then
1086              print*,'xyzpam: ***WARNING*** cluster ',icy              print*,'xyzpam: ***WARNING*** cluster icy=',icy
1087       $           ,' does not belong to plane: ',ip       $           ,' belongs to plane ',(nplanes-ipy+1)            
1088         $           ,' and not ',ip
1089              icy = -1*icy              icy = -1*icy
1090              return              return
1091           endif           endif
# Line 1100  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2. Line 1105  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1105           zm(ip) = zPAM           zm(ip) = zPAM
1106           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1107           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1108             zm_A(ip) = zPAM_A
1109           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
1110           ym_B(ip) = yPAM_B           ym_B(ip) = yPAM_B
1111             zm_B(ip) = zPAM_B
1112    
1113  c         zv(ip) = (zPAM_A+zPAM_B)/2.  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1114                    
# Line 1110  c         zv(ip) = (zPAM_A+zPAM_B)/2. Line 1117  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1117           ipx=npl(VIEW(icx))           ipx=npl(VIEW(icx))
1118    
1119           if( (nplanes-ipx+1).ne.ip )then           if( (nplanes-ipx+1).ne.ip )then
1120              print*,'xyzpam: ***WARNING*** cluster ',icx              print*,'xyzpam: ***WARNING*** cluster icx=',icx
1121       $           ,' does not belong to plane: ',ip       $           ,' belongs to plane ',(nplanes-ipx+1)            
1122         $           ,' and not ',ip
1123              icx = -1*icx              icx = -1*icx
1124              return              return
1125           endif           endif
# Line 1131  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2. Line 1139  c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1139           zm(ip) = zPAM           zm(ip) = zPAM
1140           xm_A(ip) = xPAM_A           xm_A(ip) = xPAM_A
1141           ym_A(ip) = yPAM_A           ym_A(ip) = yPAM_A
1142             zm_A(ip) = zPAM_A
1143           xm_B(ip) = xPAM_B           xm_B(ip) = xPAM_B
1144           ym_B(ip) = yPAM_B           ym_B(ip) = yPAM_B
1145             zm_B(ip) = zPAM_B
1146                    
1147  c         zv(ip) = (zPAM_A+zPAM_B)/2.  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1148    
# Line 1153  c         zv(ip) = (zPAM_A+zPAM_B)/2. Line 1163  c         zv(ip) = (zPAM_A+zPAM_B)/2.
1163           zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4           zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1164           xm_A(ip) = 0.           xm_A(ip) = 0.
1165           ym_A(ip) = 0.           ym_A(ip) = 0.
1166             zm_A(ip) = 0.
1167           xm_B(ip) = 0.           xm_B(ip) = 0.
1168           ym_B(ip) = 0.           ym_B(ip) = 0.
1169             zm_B(ip) = 0.
1170    
1171  c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4  c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1172    
# Line 1379  c--------------------------------------- Line 1391  c---------------------------------------
1391                 iv1=iside                 iv1=iside
1392                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1393  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1394                 AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))                 AA = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7
1395                 BB = yvv(iv1) - AA*xvv(iv1)                 BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7
1396  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1397                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
1398                 yoo = AA*xoo + BB                 yoo = AA*xoo + BB
# Line 1392  c--------------------------------------- Line 1404  c---------------------------------------
1404                 iv1=iside                 iv1=iside
1405                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1406  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1407                 AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))                 AA = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7
1408                 BB = xvv(iv1) - AA*yvv(iv1)                 BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7
1409  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1410                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
1411                 xoo = AA*yoo + BB                 xoo = AA*yoo + BB
# Line 1976  c$$$               print*,'(1) ip ',ip1, Line 1988  c$$$               print*,'(1) ip ',ip1,
1988  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1989  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1990                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1991                 xm1=xPAM                 xm1=REAL(xPAM) !EM GCC4.7
1992                 ym1=yPAM                 ym1=REAL(yPAM) !EM GCC4.7
1993                 zm1=zPAM                                   zm1=REAL(zPAM) !EM GCC4.7
1994    
1995                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1996  c$$$                  print*,'(2) ip ',ip2  c$$$                  print*,'(2) ip ',ip2
# Line 1999  c                        call xyz_PAM Line 2011  c                        call xyz_PAM
2011  c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)  c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2012                          call xyz_PAM                          call xyz_PAM
2013       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2014                          xm2=xPAM                          xm2=REAL(xPAM) !EM GCC4.7
2015                          ym2=yPAM                          ym2=REAL(yPAM) !EM GCC4.7
2016                          zm2=zPAM                          zm2=REAL(zPAM) !EM GCC4.7
2017    
2018  *                       ---------------------------------------------------  *                       ---------------------------------------------------
2019  *                       both couples must have a y-cluster  *                       both couples must have a y-cluster
# Line 2042  ccc                        print*,'<doub Line 2054  ccc                        print*,'<doub
2054  *     tg(th_yz)  *     tg(th_yz)
2055                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2056  *     y0 (cm)  *     y0 (cm)
2057                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1                      alfayz1(ndblt)=alfayz2(ndblt)*(REAL(zini)-zm1)+ym1! EM GCC4.7 zm1, ym1 and alfayz1/2 are REAL
2058                                                    
2059  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2060  ****  reject non phisical couples                    ****  ****  reject non phisical couples                    ****
# Line 2105  c     $                               (i Line 2117  c     $                               (i
2117                                   call xyz_PAM                                   call xyz_PAM
2118       $                                (icx3,icy3,is3,PFAdef,PFAdef       $                                (icx3,icy3,is3,PFAdef,PFAdef
2119       $                                ,0.,0.,0.,0.)       $                                ,0.,0.,0.,0.)
2120                                   xm3=xPAM                                   xm3=REAL(xPAM) !EM GCC4.7
2121                                   ym3=yPAM                                   ym3=REAL(yPAM) !EM GCC4.7
2122                                   zm3=zPAM                                   zm3=REAL(zPAM) !EM GCC4.7
2123    
2124    
2125  *     find the circle passing through the three points  *     find the circle passing through the three points
# Line 2144  cc                                 if(if Line 2156  cc                                 if(if
2156                                      DET=SZZ*S1-SZ*SZ                                      DET=SZZ*S1-SZ*SZ
2157                                      AX=(SZX*S1-SZ*SSX)/DET                                      AX=(SZX*S1-SZ*SSX)/DET
2158                                      BX=(SZZ*SSX-SZX*SZ)/DET                                      BX=(SZZ*SSX-SZX*SZ)/DET
2159                                      X0  = AX*ZINI+BX                                      X0  = AX*REAL(ZINI)+BX ! EM GCC4.7
2160                                                                            
2161                                   endif                                   endif
2162    
# Line 2183  ccc                                 prin Line 2195  ccc                                 prin
2195                                                                    
2196                                   if(radius.ne.0.and.xc.lt.0)then                                   if(radius.ne.0.and.xc.lt.0)then
2197  *************POSITIVE DEFLECTION  *************POSITIVE DEFLECTION
2198               alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc+sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2199               alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = (REAL(ZINI)-zc)/
2200               alfaxz3(ntrpt) = 1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2201               alfaxz3(ntrpt) = 1/radius
2202                                  else if(radius.ne.0.and.xc.ge.0)then                                  else if(radius.ne.0.and.xc.ge.0)then
2203  *************NEGATIVE DEFLECTION  *************NEGATIVE DEFLECTION
2204               alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc-sqrt(radius**2-(REAL(ZINI)-zc)**2)
2205               alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/
2206               alfaxz3(ntrpt) = -1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2207               alfaxz3(ntrpt) = -1/radius
2208                                  else if(radius.eq.0)then                                  else if(radius.eq.0)then
2209  *************straight fit  *************straight fit
2210               alfaxz1(ntrpt) = X0               alfaxz1(ntrpt) = X0
# Line 3104  c                                    mas Line 3118  c                                    mas
3118                                   enddo                                   enddo
3119                                enddo                                enddo
3120                                                                
3121                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=REAL(chi2)
3122                                                                
3123  *     --------------------------------  *     --------------------------------
3124  *     STORE candidate TRACK INFO - end  *     STORE candidate TRACK INFO - end
# Line 3172  cc         return Line 3186  cc         return
3186        character*10 PFA        character*10 PFA
3187        common/FINALPFA/PFA        common/FINALPFA/PFA
3188    
3189          double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7
3190          double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7
3191          double precision zmm,zmm_A,zmm_B !EM GCC4.7
3192          double precision clincnewc !EM GCC4.7
3193          double precision clincnew !EM GCC4.7
3194    
3195        real k(6)        real k(6)
3196        DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/        DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3197    
# Line 3390  c               distance = distance / RC Line 3410  c               distance = distance / RC
3410                    idm = id                                      idm = id                  
3411                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3412                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3413                    clincnewc=10*sqrt(rymm**2+rxmm**2                    clincnewc=10.*dsqrt(rymm**2+rxmm**2
3414       $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))       $            +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))))! EM GCC4.7
3415                 endif                 endif
3416   1188          continue   1188          continue
3417              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
# Line 3559  c               distance = distance / RC Line 3579  c               distance = distance / RC
3579              if(iclm.ne.0)then              if(iclm.ne.0)then
3580                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3581                    clincnew=                    clincnew=
3582       $                 20*       $            20.*     !EM GCC4.7
3583       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))       $            dsqrt(rxmm**2 +
3584         $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1))
3585                 else if(mod(VIEW(iclm),2).ne.0)then                 else if(mod(VIEW(iclm),2).ne.0)then
3586                    clincnew=                    clincnew=
3587       $                 10*       $            10.* !EM GCC4.7
3588       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))       $            dsqrt(rymm**2 +
3589         $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2))
3590                 endif                 endif
3591    
3592                 if(distmin.le.clincnew)then                   if(distmin.le.clincnew)then  
# Line 3813  c               distance = distance / RC Line 3835  c               distance = distance / RC
3835        character*10 PFA        character*10 PFA
3836        common/FINALPFA/PFA        common/FINALPFA/PFA
3837    
3838        real sinth,phi,pig        real sinth,phi,pig, npig ! EM GCC4.7
3839        integer ssensor,sladder        integer ssensor,sladder
3840        pig=acos(-1.)        pig=acos(-1.)
3841    
# Line 3823  c               distance = distance / RC Line 3845  c               distance = distance / RC
3845        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3846        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3847  *     -------------------------------------  *     -------------------------------------
3848        phi   = al(4)                  phi   = REAL(al(4))
3849        sinth = al(3)                    sinth = REAL(al(3))
3850        if(sinth.lt.0)then              if(sinth.lt.0)then      
3851           sinth = -sinth                   sinth = -sinth        
3852           phi = phi + pig                 phi = phi + pig      
# Line 4047  c         if( mask_view(iv).ne.0 )good2( Line 4069  c         if( mask_view(iv).ne.0 )good2(
4069                 do is=1,2                 do is=1,2
4070  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4071  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4072                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)                    call xyz_PAM(icl,0,is,PFAdef,'    ',0.,0.,0.,0.)
4073                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7
4074                 enddo                 enddo
4075              else                          !=== Y views              else                          !=== Y views
4076                 nclsy = nclsy + 1                 nclsy = nclsy + 1
# Line 4063  c                  call xyz_PAM(icl,0,is Line 4085  c                  call xyz_PAM(icl,0,is
4085                 do is=1,2                 do is=1,2
4086  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4087  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4088                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)                    call xyz_PAM(0,icl,is,'    ',PFAdef,0.,0.,0.,0.)
4089                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7
4090                 enddo                 enddo
4091              endif              endif
4092           endif           endif
Legend:
Removed from v.1.44  
changed lines
  Added in v.1.47
  ViewVC Help
Powered by ViewVC 1.1.23