(file) Return to h_track_fit.f CVS log (file) (dir) Up to [HallC] / Analyzer / HTRACKING

Diff for /Analyzer/HTRACKING/h_track_fit.f between version 1.10 and 1.11

version 1.10, 1995/08/30 16:11:39 version 1.11, 1996/01/16 21:42:18
Line 9 
Line 9 
 *                              remove minuit. Make fit linear *                              remove minuit. Make fit linear
 *                              still does not do errors properly *                              still does not do errors properly
 * $Log$ * $Log$
   * Revision 1.11  1996/01/16 21:42:18  cdaq
   * (JRA) Remove slices code, misc fixes, reindent.
   *
 * Revision 1.10  1995/08/30 16:11:39  cdaq * Revision 1.10  1995/08/30 16:11:39  cdaq
 * (JRA) Don't fill single_residual arrray * (JRA) Don't fill single_residual arrray
 * *
Line 45 
Line 48 
       include "hms_data_structures.cmn"       include "hms_data_structures.cmn"
       include "hms_tracking.cmn"       include "hms_tracking.cmn"
       include "hms_geometry.cmn"       include "hms_geometry.cmn"
       external H_DPSIFUN        external h_dpsifun
 * *
 * *
 *     local variables *     local variables
 * *
       logical ABORT       logical ABORT
       character*50 here        character*11 here
       parameter (here='H_TRACK_FIT')       parameter (here='H_TRACK_FIT')
       character*(*) err       character*(*) err
       integer*4 itrack                        ! track loop index        integer*4 itrk                        ! track loop index
       integer*4 ihit,ierr       integer*4 ihit,ierr
       integer*4 hit,plane        integer*4 hit,pln
       integer*4 i,j                             ! loop index       integer*4 i,j                             ! loop index
 *      real*4 z_slice *      real*4 z_slice
  
       real*8   H_DPSIFUN        real*8   h_dpsifun
       real*8   pos       real*8   pos
       real*4   ray(hnum_fpray_param)  
       real*8   ray1(4)       real*8   ray1(4)
       real*8   ray2(4)       real*8   ray2(4)
       real*8 TT(hnum_fpray_param)       real*8 TT(hnum_fpray_param)
       real*8 AA(hnum_fpray_param,hnum_fpray_param)       real*8 AA(hnum_fpray_param,hnum_fpray_param)
       real*8 dray(hnum_fpray_param)       real*8 dray(hnum_fpray_param)
 *      real*8 error(hnum_fpray_param)  
       real*4 chi2,dummychi2       real*4 chi2,dummychi2
       parameter (dummychi2 = 1.E4)       parameter (dummychi2 = 1.E4)
 *     array to remap hplane_coeff to param number *     array to remap hplane_coeff to param number
Line 80 
Line 81 
       ierr=0       ierr=0
 *  initailize residuals *  initailize residuals
  
       do plane=1,HDC_NUM_PLANES        do pln=1,hdc_num_planes
         do itrack=1,HNTRACKS_MAX          do itrk=1,hntracks_fp
           hdc_double_residual(itrack,plane)=1000            hdc_double_residual(itrk,pln)=1000
           hdc_single_residual(itrack,plane)=1000            hdc_single_residual(itrk,pln)=1000
         enddo         enddo
 c fill the 1d arrays from the 2d arrays for good track (in h_physics) c fill the 1d arrays from the 2d arrays for good track (in h_physics)
 c        hdc_sing_res(plane)=1000  c        hdc_sing_res(pln)=1000
         hdc_dbl_res(plane)=1000          hdc_dbl_res(pln)=1000
       enddo       enddo
  
 *     test for no tracks *     test for no tracks
       if(HNTRACKS_FP.ge.1) then        if(hntracks_fp.ge.1) then
         do itrack=1,HNTRACKS_FP          do itrk=1,hntracks_fp
           chi2= dummychi2           chi2= dummychi2
           htrack_fit_num=itrack            htrack_fit_num=itrk
  
 *     are there enough degrees of freedom *     are there enough degrees of freedom
           HNFREE_FP(itrack)=HNTRACK_HITS(itrack,1)-hnum_fpray_param            hnfree_fp(itrk)=hntrack_hits(itrk,1)-hnum_fpray_param
           if(HNFREE_FP(itrack).gt.0) then            if(hnfree_fp(itrk).gt.0) then
  
 *     initialize parameters *     initialize parameters
             do i=1,hnum_fpray_param             do i=1,hnum_fpray_param
               TT(i)=0.               TT(i)=0.
               do ihit=2,HNTRACK_HITS(itrack,1)+1                do ihit=2,hntrack_hits(itrk,1)+1
                 hit=HNTRACK_HITS(itrack,ihit)                  hit=hntrack_hits(itrk,ihit)
                 plane=HDC_PLANE_NUM(hit)                  pln=hdc_plane_num(hit)
                 TT(i)=TT(i)+((HDC_WIRE_COORD(hit)*                  TT(i)=TT(i)+((hdc_wire_coord(hit)*
      &               hplane_coeff(remap(i),plane))       &               hplane_coeff(remap(i),pln))
      &               /(hdc_sigma(plane)*hdc_sigma(plane)))       &               /(hdc_sigma(pln)*hdc_sigma(pln)))
               enddo               enddo
             enddo             enddo
             do i=1,hnum_fpray_param             do i=1,hnum_fpray_param
Line 117 
Line 118 
                 if(j.lt.i)then                 if(j.lt.i)then
                   AA(i,j)=AA(j,i)                   AA(i,j)=AA(j,i)
                 else                 else
                   do ihit=2,HNTRACK_HITS(itrack,1)+1                    do ihit=2,hntrack_hits(itrk,1)+1
                     hit=HNTRACK_HITS(itrack,ihit)                      hit=hntrack_hits(itrk,ihit)
                     plane=HDC_PLANE_NUM(hit)                      pln=hdc_plane_num(hit)
                     AA(i,j)=AA(i,j) + (                     AA(i,j)=AA(i,j) + (
      &                   hplane_coeff(remap(i),plane)*hplane_coeff(remap(j)       &                   hplane_coeff(remap(i),pln)*hplane_coeff(remap(j)
      $                   ,plane)/(hdc_sigma(plane)*hdc_sigma(plane)))       $                   ,pln)/(hdc_sigma(pln)*hdc_sigma(pln)))
                   enddo                 ! end loop on ihit                   enddo                 ! end loop on ihit
                 endif                   ! end test on j .lt. i                 endif                   ! end test on j .lt. i
               enddo                     ! end loop on j               enddo                     ! end loop on j
Line 140 
Line 141 
 *     calculate chi2 *     calculate chi2
               chi2=0.               chi2=0.
  
               ray(1)=dray(1)  * calculate hit coord at each plane for chisquared and efficiency calculations.
               ray(2)=dray(2)                do pln=1,hdc_num_planes
               ray(3)=dray(3)                  hdc_track_coord(itrk,pln)=hplane_coeff(remap(1),pln)*dray(1)
               ray(4)=dray(4)       &               +hplane_coeff(remap(2),pln)*dray(2)
               do ihit=2,HNTRACK_HITS(itrack,1)+1       &               +hplane_coeff(remap(3),pln)*dray(3)
                 hit=HNTRACK_HITS(itrack,ihit)       &               +hplane_coeff(remap(4),pln)*dray(4)
                 plane=HDC_PLANE_NUM(hit)                enddo
   
                 do ihit=2,hntrack_hits(itrk,1)+1
                   hit=hntrack_hits(itrk,ihit)
                   pln=hdc_plane_num(hit)
  
 * note chi2 is single precision * note chi2 is single precision
                 hdc_single_residual(itrack,plane)=  
      &               (HDC_WIRE_COORD(hit)                  hdc_single_residual(itrk,pln)=
      &               -hplane_coeff(remap(1),plane)*ray(1)       &              hdc_wire_coord(hit)-hdc_track_coord(itrk,pln)
      &               -hplane_coeff(remap(2),plane)*ray(2)                  chi2=chi2+
      &               -hplane_coeff(remap(3),plane)*ray(3)       &              (hdc_single_residual(itrk,pln)/hdc_sigma(pln))**2
      &               -hplane_coeff(remap(4),plane)*ray(4))  
                 chi2=chi2+((hdc_single_residual(itrack,plane))**2  
      &               )/(hdc_sigma(plane)*hdc_sigma(plane))  
   
 * Fill single_residual array.  Note that due to ihit loop, the plane  
 * will always contain a wire on the track being examined.  
 c                hdc_sing_res(plane) = hdc_single_residual(itrack,plane)  
               enddo               enddo
             endif             endif
  
             HX_FP(itrack)=ray(1)              hx_fp(itrk)=dray(1)
             HY_FP(itrack)=ray(2)              hy_fp(itrk)=dray(2)
             HZ_FP(itrack)=0.            ! z=0 of tracking.              hz_fp(itrk)=0.            ! z=0 of tracking.
             HXP_FP(itrack)=ray(3)              hxp_fp(itrk)=dray(3)
             HYP_FP(itrack)=ray(4)              hyp_fp(itrk)=dray(4)
           endif                         ! end test on degrees of freedom           endif                         ! end test on degrees of freedom
           HCHI2_FP(itrack)=chi2            hchi2_fp(itrk)=chi2
         enddo                           ! end loop over tracks         enddo                           ! end loop over tracks
       endif       endif
  
 **     A reasonable selection of slices is presently -80,-60,-40,-20,0,20,40  
 **     ,60,80.   Zero is the nominal midplane between the chambers, -80  
 **     corresponds closely to the exit flange position.  This slice pattern  
 **     is created with hz_wild=-80., hdelta_z_wild=20., and hnum_zslice=9  
 *      if(hz_slice_enable.ne.0)then  
 *        do k=1,hnum_zslice  
 *          z_slice = hz_wild + (k-1)*hdelta_z_wild  
 *          hx_fp_wild(k) = hx_fp(1) + hxp_fp(1)*z_slice  
 *          hy_fp_wild(k) = hy_fp(1) + hyp_fp(1)*z_slice  
 *        enddo  
 *      endif  
   
 * calculate residuals for each chamber if in single stub mode * calculate residuals for each chamber if in single stub mode
 * and there were 2 tracks found one in first chanber and one in the second * and there were 2 tracks found one in first chanber and one in the second
  
       if (hsingle_stub.ne.0) then       if (hsingle_stub.ne.0) then
         if (HNTRACKS_FP.eq.2) then          if (hntracks_fp.eq.2) then
           itrack=1            itrk=1
           ihit=2           ihit=2
           hit=HNTRACK_HITS(itrack,ihit)            hit=hntrack_hits(itrk,ihit)
           plane=HDC_PLANE_NUM(hit)            pln=hdc_plane_num(hit)
           if (plane.le.6) then            if (pln.le.6) then
             itrack=2              itrk=2
             hit=HNTRACK_HITS(itrack,ihit)              hit=hntrack_hits(itrk,ihit)
             plane=HDC_PLANE_NUM(hit)              pln=hdc_plane_num(hit)
             if (plane.ge.7) then              if (pln.ge.7) then
  
 * condition of above met calculating residuals * condition of above met calculating residuals
 * assigning rays to tracks in each chamber * assigning rays to tracks in each chamber
 * ray1 is ray from first chamber fit * ray1 is ray from first chamber fit
 * ray2 is ray from second chamber fit * ray2 is ray from second chamber fit
  
               ray1(1)=dble(HX_FP(1))                ray1(1)=dble(hx_fp(1))
               ray1(2)=dble(HY_FP(1))                ray1(2)=dble(hy_fp(1))
               ray1(3)=dble(HXP_FP(1))                ray1(3)=dble(hxp_fp(1))
               ray1(4)=dble(HYP_FP(1))                ray1(4)=dble(hyp_fp(1))
               ray2(1)=dble(HX_FP(2))                ray2(1)=dble(hx_fp(2))
               ray2(2)=dble(HY_FP(2))                ray2(2)=dble(hy_fp(2))
               ray2(3)=dble(HXP_FP(2))                ray2(3)=dble(hxp_fp(2))
               ray2(4)=dble(HYP_FP(2))                ray2(4)=dble(hyp_fp(2))
  
               itrack=1                itrk=1
 * loop over hits in second chamber * loop over hits in second chamber
               do ihit=1,HNTRACK_HITS(itrack+1,1)                do ihit=1,hntrack_hits(itrk+1,1)
  
 * calculate residual in second chamber from first chamber track * calculate residual in second chamber from first chamber track
                 hit=HNTRACK_HITS(itrack+1,ihit+1)                  hit=hntrack_hits(itrk+1,ihit+1)
                 plane=HDC_PLANE_NUM(hit)                  pln=hdc_plane_num(hit)
                 pos=H_DPSIFUN(ray1,plane)                  pos=h_dpsifun(ray1,pln)
                 hdc_double_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos                  hdc_double_residual(itrk,pln)=hdc_wire_coord(hit)-pos
 * djm 8/31/94 stuff this variable into 1d array we can register * djm 8/31/94 stuff this variable into 1d array we can register
                 hdc_dbl_res(plane) = hdc_double_residual(1,plane)                  hdc_dbl_res(pln) = hdc_double_residual(1,pln)
  
               enddo               enddo
  
               itrack=2                itrk=2
 * loop over hits in first chamber * loop over hits in first chamber
               do ihit=1,HNTRACK_HITS(itrack-1,1)                do ihit=1,hntrack_hits(itrk-1,1)
  
 * calculate residual in first chamber from second chamber track * calculate residual in first chamber from second chamber track
                 hit=HNTRACK_HITS(itrack-1,ihit+1)                  hit=hntrack_hits(itrk-1,ihit+1)
                 plane=HDC_PLANE_NUM(hit)                  pln=hdc_plane_num(hit)
                 pos=H_DPSIFUN(ray2,plane)                  pos=h_dpsifun(ray2,pln)
                 hdc_double_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos                  hdc_double_residual(itrk,pln)=hdc_wire_coord(hit)-pos
 * djm 8/31/94 stuff this variable into 1d array we can register * djm 8/31/94 stuff this variable into 1d array we can register
                 hdc_dbl_res(plane) = hdc_double_residual(2,plane)                  hdc_dbl_res(pln) = hdc_double_residual(2,pln)
  
               enddo               enddo
             endif                       ! end plane ge 7              endif                       ! end pln ge 7
           endif                         ! end plane le 6            endif                         ! end pln le 6
         endif                           ! end HNTRACKS_FP eq 2          endif                           ! end hntracks_fp eq 2
       endif                             ! end hsignle_stub .ne. 0       endif                             ! end hsignle_stub .ne. 0
  
 *     test if we want to dump out trackfit results *     test if we want to dump out trackfit results

Legend:
Removed from v.1.10  
changed lines
  Added in v.1.11
Analyzer/Replay: Mark Jones, Documents: Stephen Wood
 Powered by
ViewCVS 0.9.2-cvsgraph-1.4.0