Return to h_track_fit.f CVS log

Up to [HallC] / Analyzer / HTRACKING

(JRA) Remove slices code, misc fixes, reindent.

      subroutine H_TRACK_FIT(ABORT,err,ierr)
*     primary track fitting routine for the HMS spectrometer
*
*     Called by H_TRACK
*
*     d.f. geesaman         17 January 1994
*     modified 
*                           17 Feb 1994        dfg
*                              remove minuit. Make fit linear
*                              still does not do errors properly
* $Log: h_track_fit.f,v $
* 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
* (JRA) Don't fill single_residual arrray
*
* Revision 1.9  1995/05/22  19:39:31  cdaq
* (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts"
*
* Revision 1.8  1995/04/06  19:32:58  cdaq
* (JRA) Rename residuals variables
*
* Revision 1.7  1995/01/27  20:26:50  cdaq
* (JRA) Remove Mack's personal focalplane diceamatic (z slicer) code
*
* Revision 1.6  1994/12/06  15:45:27  cdaq
* (DJM) Take slices in Z to look for best focus
*
* Revision 1.5  1994/10/12  18:52:06  cdaq
* (DJM) Initialize some variables
* (SAW) Prettify indentation
*
* Revision 1.4  1994/09/01  12:29:07  cdaq
* (DJM) Make registered versions of residuals
*
* Revision 1.3  1994/08/18  02:45:53  cdaq
* (DM) Add calculation of residuals
*
* Revision 1.2  1994/02/22  05:25:00  cdaq
* (SAW) Remove dfloat calls with floating args
*
* Revision 1.1  1994/02/19  06:20:53  cdaq
* Initial revision
*
*
      implicit none
      include "hms_data_structures.cmn"
      include "hms_tracking.cmn"
      include "hms_geometry.cmn"
      external h_dpsifun
*
*
*     local variables
*
      logical ABORT
      character*11 here
      parameter (here='H_TRACK_FIT')
      character*(*) err
      integer*4 itrk                        ! track loop index
      integer*4 ihit,ierr
      integer*4 hit,pln
      integer*4 i,j                             ! loop index
*      real*4 z_slice

      real*8   h_dpsifun
      real*8   pos
      real*8   ray1(4)
      real*8   ray2(4)
      real*8 TT(hnum_fpray_param)
      real*8 AA(hnum_fpray_param,hnum_fpray_param)
      real*8 dray(hnum_fpray_param)
      real*4 chi2,dummychi2
      parameter (dummychi2 = 1.E4)
*     array to remap hplane_coeff to param number
      integer*4 remap(hnum_fpray_param)
      data remap/5,6,3,4/
      save remap
*
      ABORT= .FALSE.
      ierr=0
*  initailize residuals

      do pln=1,hdc_num_planes
        do itrk=1,hntracks_fp
          hdc_double_residual(itrk,pln)=1000
          hdc_single_residual(itrk,pln)=1000
        enddo
c fill the 1d arrays from the 2d arrays for good track (in h_physics)
c        hdc_sing_res(pln)=1000
        hdc_dbl_res(pln)=1000
      enddo

*     test for no tracks
      if(hntracks_fp.ge.1) then
        do itrk=1,hntracks_fp
          chi2= dummychi2
          htrack_fit_num=itrk

*     are there enough degrees of freedom
          hnfree_fp(itrk)=hntrack_hits(itrk,1)-hnum_fpray_param
          if(hnfree_fp(itrk).gt.0) then

*     initialize parameters
            do i=1,hnum_fpray_param
              TT(i)=0.
              do ihit=2,hntrack_hits(itrk,1)+1
                hit=hntrack_hits(itrk,ihit)
                pln=hdc_plane_num(hit)
                TT(i)=TT(i)+((hdc_wire_coord(hit)*
     &               hplane_coeff(remap(i),pln))
     &               /(hdc_sigma(pln)*hdc_sigma(pln)))
              enddo 
            enddo
            do i=1,hnum_fpray_param
              do j=1,hnum_fpray_param
                AA(i,j)=0.
                if(j.lt.i)then
                  AA(i,j)=AA(j,i)
                else 
                  do ihit=2,hntrack_hits(itrk,1)+1
                    hit=hntrack_hits(itrk,ihit)
                    pln=hdc_plane_num(hit)
                    AA(i,j)=AA(i,j) + (
     &                   hplane_coeff(remap(i),pln)*hplane_coeff(remap(j)
     $                   ,pln)/(hdc_sigma(pln)*hdc_sigma(pln)))
                  enddo                 ! end loop on ihit
                endif                   ! end test on j .lt. i
              enddo                     ! end loop on j
            enddo                       ! end loop on i
*
*     solve four by four equations
            call solve_four_by_four(TT,AA,dray,ierr)
*
            if(ierr.ne.0) then
              dray(1)=10000.
              dray(2)=10000.
              dray(3)=2.
              dray(4)=2.
            else
*     calculate chi2
              chi2=0.
              
* calculate hit coord at each plane for chisquared and efficiency calculations.
              do pln=1,hdc_num_planes
                hdc_track_coord(itrk,pln)=hplane_coeff(remap(1),pln)*dray(1)
     &               +hplane_coeff(remap(2),pln)*dray(2)
     &               +hplane_coeff(remap(3),pln)*dray(3)
     &               +hplane_coeff(remap(4),pln)*dray(4)
              enddo

              do ihit=2,hntrack_hits(itrk,1)+1
                hit=hntrack_hits(itrk,ihit)
                pln=hdc_plane_num(hit)

* note chi2 is single precision

                hdc_single_residual(itrk,pln)=
     &              hdc_wire_coord(hit)-hdc_track_coord(itrk,pln)
                chi2=chi2+
     &              (hdc_single_residual(itrk,pln)/hdc_sigma(pln))**2
              enddo
            endif

            hx_fp(itrk)=dray(1)
            hy_fp(itrk)=dray(2)
            hz_fp(itrk)=0.            ! z=0 of tracking.
            hxp_fp(itrk)=dray(3)
            hyp_fp(itrk)=dray(4)
          endif                         ! end test on degrees of freedom
          hchi2_fp(itrk)=chi2
        enddo                           ! end loop over tracks
      endif

* 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

      if (hsingle_stub.ne.0) then
        if (hntracks_fp.eq.2) then
          itrk=1
          ihit=2
          hit=hntrack_hits(itrk,ihit)
          pln=hdc_plane_num(hit)
          if (pln.le.6) then
            itrk=2
            hit=hntrack_hits(itrk,ihit)
            pln=hdc_plane_num(hit)
            if (pln.ge.7) then

* condition of above met calculating residuals  
* assigning rays to tracks in each chamber
* ray1 is ray from first chamber fit
* ray2 is ray from second chamber fit

              ray1(1)=dble(hx_fp(1))
              ray1(2)=dble(hy_fp(1))
              ray1(3)=dble(hxp_fp(1))
              ray1(4)=dble(hyp_fp(1))
              ray2(1)=dble(hx_fp(2))
              ray2(2)=dble(hy_fp(2))
              ray2(3)=dble(hxp_fp(2))
              ray2(4)=dble(hyp_fp(2))

              itrk=1
* loop over hits in second chamber
              do ihit=1,hntrack_hits(itrk+1,1)

* calculate residual in second chamber from first chamber track
                hit=hntrack_hits(itrk+1,ihit+1)
                pln=hdc_plane_num(hit)
                pos=h_dpsifun(ray1,pln)
                hdc_double_residual(itrk,pln)=hdc_wire_coord(hit)-pos
* djm 8/31/94 stuff this variable into 1d array we can register
                hdc_dbl_res(pln) = hdc_double_residual(1,pln)

              enddo

              itrk=2
* loop over hits in first chamber
              do ihit=1,hntrack_hits(itrk-1,1)

* calculate residual in first chamber from second chamber track
                hit=hntrack_hits(itrk-1,ihit+1)
                pln=hdc_plane_num(hit)
                pos=h_dpsifun(ray2,pln)
                hdc_double_residual(itrk,pln)=hdc_wire_coord(hit)-pos
* djm 8/31/94 stuff this variable into 1d array we can register
                hdc_dbl_res(pln) = hdc_double_residual(2,pln)

              enddo
            endif                       ! end pln ge 7
          endif                         ! end pln le 6
        endif                           ! end hntracks_fp eq 2
      endif                             ! end hsignle_stub .ne. 0

*     test if we want to dump out trackfit results
      if(hdebugtrackprint.ne.0) then
        call h_print_tracks
      endif                             ! end test on zero tracks
 1000 return
      end

*