1 cdaq  1.1       subroutine H_TRACK_FIT(ABORT,err,ierr)
  2           *     primary track fitting routine for the HMS spectrometer
  3           *
  4           *     Called by H_TRACK
  5           *
  6           *     d.f. geesaman         17 January 1994
  7           *     modified 
  8           *                           17 Feb 1994        dfg
  9           *                              remove minuit. Make fit linear
 10           *                              still does not do errors properly

 11 cdaq  1.2 * $Log: h_track_fit.f,v $

 12 cdaq  1.4 * Revision 1.3  1994/08/18  02:45:53  cdaq
 13           * (DM) Add calculation of residuals
 14           *

 15 cdaq  1.3 * Revision 1.2  1994/02/22  05:25:00  cdaq
 16           * (SAW) Remove dfloat calls with floating args
 17           *

 18 cdaq  1.2 * Revision 1.1  1994/02/19  06:20:53  cdaq
 19           * Initial revision

 20 cdaq  1.1 *

 21 cdaq  1.2 *

 22 cdaq  1.1       implicit none
 23                 include "gen_data_structures.cmn"
 24                 include "hms_tracking.cmn"
 25                 include "hms_geometry.cmn"

 26 cdaq  1.3       external H_DPSIFUN

 27 cdaq  1.1 *
 28           *
 29           *     local variables
 30           *
 31                 logical ABORT
 32                 character*50 here
 33                 parameter (here='H_TRACK_FIT')
 34                 character*(*) err
 35                 integer*4 itrack                        ! track loop index
 36                 integer*4 ihit,ierr
 37                 integer*4 hit,plane
 38                 integer*4 i,j                             ! loop index

 39 cdaq  1.3       real*8   H_DPSIFUN
 40                 real*8   pos

 41 cdaq  1.1       real*4   ray(hnum_fpray_param)

 42 cdaq  1.3       real*8   ray1(4)
 43                 real*8   ray2(4)

 44 cdaq  1.1       real*8 TT(hnum_fpray_param)

 45 cdaq  1.3       real*8 AA(hnum_fpray_param,hnum_fpray_param)

 46 cdaq  1.1       real*8 dray(hnum_fpray_param)
 47                 real*8 error(hnum_fpray_param)
 48                 real*4 chi2,dummychi2
 49                 parameter (dummychi2 = 1.E4)
 50           *     array to remap hplane_coeff to param number
 51                 integer*4 remap(hnum_fpray_param)
 52                 data remap/5,6,3,4/
 53                 save remap
 54           *
 55                 ABORT= .FALSE.
 56                 ierr=0

 57 cdaq  1.3 *  initailize residuals
 58           
 59                 do itrack=1,HNTRACKS_MAX
 60                   do plane=1,HMAX_NUM_DC_PLANES
 61                     hdc_residual(itrack,plane)=1000
 62                     hdc_sing_res(itrack,plane)=1000
 63                   enddo
 64                 enddo
 65           
 66           

 67 cdaq  1.1 *     test for no tracks
 68                 if(HNTRACKS_FP.ge.1) then
 69                   do itrack=1,HNTRACKS_FP
 70                     chi2= dummychi2
 71                     htrack_fit_num=itrack

 72 cdaq  1.3 

 73 cdaq  1.1 *     are there enough degrees of freedom
 74                     HNFREE_FP(itrack)=HNTRACK_HITS(itrack,1)-hnum_fpray_param
 75                     if(HNFREE_FP(itrack).gt.0) then

 76 cdaq  1.3 

 77 cdaq  1.1 *     initialize parameters
 78                      do i=1,hnum_fpray_param
 79                       TT(i)=0.
 80                       do ihit=2,HNTRACK_HITS(itrack,1)+1
 81                         hit=HNTRACK_HITS(itrack,ihit)
 82                         plane=HDC_PLANE_NUM(hit)

 83 cdaq  1.2               TT(i)=TT(i)+((HDC_WIRE_COORD(hit)*

 84 cdaq  1.1      &          hplane_coeff(remap(i),plane))
 85                &          /(hdc_sigma(plane)*hdc_sigma(plane)))
 86                       enddo 
 87                      enddo
 88                     do i=1,hnum_fpray_param
 89                      do j=1,hnum_fpray_param
 90                      AA(i,j)=0.
 91                       do ihit=2,HNTRACK_HITS(itrack,1)+1
 92                         hit=HNTRACK_HITS(itrack,ihit)
 93                         plane=HDC_PLANE_NUM(hit)

 94 cdaq  1.2               AA(i,j)=AA(i,j) + (

 95 cdaq  1.1      &           hplane_coeff(remap(i),plane)*hplane_coeff(remap(j),plane)
 96                &          /(hdc_sigma(plane)*hdc_sigma(plane)))
 97                       enddo               ! end loop on ihit
 98                      enddo                ! end loop on j
 99                     enddo                 ! end loop on i
100           *
101           *     solve four by four equations
102                    call solve_four_by_four(TT,AA,dray,ierr)
103           *
104                 if(ierr.ne.0) then
105                    dray(1)=10000.
106                    dray(2)=10000.
107                    dray(3)=2.
108                    dray(4)=2.
109                 else
110           *      calculate chi2
111                   chi2=0.

112 cdaq  1.3 

113 cdaq  1.1         ray(1)=dray(1)
114                   ray(2)=dray(2)
115                   ray(3)=dray(3)
116                   ray(4)=dray(4)
117                    do ihit=2,HNTRACK_HITS(itrack,1)+1
118                         hit=HNTRACK_HITS(itrack,ihit)
119                         plane=HDC_PLANE_NUM(hit)

120 cdaq  1.3 

121 cdaq  1.1 * note chi2 is single precision

122 cdaq  1.3            hdc_sing_res(itrack,plane)=

123 cdaq  1.1      &         (HDC_WIRE_COORD(hit)
124                &         -hplane_coeff(remap(1),plane)*ray(1)
125                &         -hplane_coeff(remap(2),plane)*ray(2)
126                &         -hplane_coeff(remap(3),plane)*ray(3)

127 cdaq  1.3      &         -hplane_coeff(remap(4),plane)*ray(4))
128                      chi2=chi2+((hdc_sing_res(itrack,plane))**2
129                &         )/(hdc_sigma(plane)*hdc_sigma(plane))

130 cdaq  1.4 *djm
131                       if(itrack.eq.1 .and. plane.ge.1 .and. plane.le.6)
132                &        hdc1_sing_res(plane) = hdc_sing_res(itrack,plane)
133                       if(itrack.eq.2 .and. plane.ge.7 .and. plane.le.12)
134                &        hdc2_sing_res(plane) = hdc_sing_res(itrack,plane)

135 cdaq  1.1         enddo
136                 endif
137           
138                      HX_FP(itrack)=ray(1)
139                      HY_FP(itrack)=ray(2)
140                      HZ_FP(itrack)=0.                   ! z=0 of tracking.
141                      HXP_FP(itrack)=ray(3)
142                      HYP_FP(itrack)=ray(4)
143                     endif                               ! end test on degrees of freedom
144                     HCHI2_FP(itrack)=chi2
145                   enddo                                 ! end loop over tracks
146                 endif

147 cdaq  1.3 
148           * calculate residuals for each chamber if in single stub mode
149           * and there were 2 tracks found one in first chanber and one in the second
150           
151                 if (hsingle_stub.ne.0) then
152                   if (HNTRACKS_FP.eq.2) then
153                     itrack=1
154                     ihit=2
155                     hit=HNTRACK_HITS(itrack,ihit)
156                     plane=HDC_PLANE_NUM(hit)
157                     if (plane.le.6) then
158                       itrack=2
159                       hit=HNTRACK_HITS(itrack,ihit)
160                       plane=HDC_PLANE_NUM(hit)
161                       if (plane.ge.7) then
162           
163           * condition of above met calculating residuals  
164           * assigning rays to tracks in each chamber
165           * ray1 is ray from first chamber fit
166           * ray2 is ray from second chamber fit
167           
168 cdaq  1.3             ray1(1)=dble(HX_FP(1))
169                       ray1(2)=dble(HY_FP(1))
170                       ray1(3)=dble(HXP_FP(1))
171                       ray1(4)=dble(HYP_FP(1))
172                       ray2(1)=dble(HX_FP(2))
173                       ray2(2)=dble(HY_FP(2))
174                       ray2(3)=dble(HXP_FP(2))
175                       ray2(4)=dble(HYP_FP(2))
176           
177                       itrack=1
178           * loop over hits in second chamber
179                       do ihit=1,HNTRACK_HITS(itrack+1,1)
180           
181           * calculate residual in second chamber from first chamber track
182                         hit=HNTRACK_HITS(itrack+1,ihit+1)
183                         plane=HDC_PLANE_NUM(hit)
184                         pos=H_DPSIFUN(ray1,plane)
185                         hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos

186 cdaq  1.4 * djm 8/31/94 stuff this variable into 1d array we can register
187                         hdc2_dbl_res(plane) = hdc_residual(1,plane)

188 cdaq  1.3 
189                       enddo
190           
191                       itrack=2
192           * loop over hits in first chamber
193                       do ihit=1,HNTRACK_HITS(itrack-1,1)
194           

195 cdaq  1.4 * calculate residual in first chamber from second chamber track

196 cdaq  1.3               hit=HNTRACK_HITS(itrack-1,ihit+1)
197                         plane=HDC_PLANE_NUM(hit)
198                         pos=H_DPSIFUN(ray2,plane)
199                         hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos

200 cdaq  1.4 * djm 8/31/94 stuff this variable into 1d array we can register
201                         hdc1_dbl_res(plane) = hdc_residual(2,plane)

202 cdaq  1.3 
203                       enddo

204 cdaq  1.4             endif                       ! end plane ge 7
205                     endif                         ! end plane le 6
206                   endif                           ! end HNTRACKS_FP eq 2
207                 endif                             ! end hsignle_stub .ne. 0

208 cdaq  1.3 

209 cdaq  1.1 *     test if we want to dump out trackfit results
210                 if(hdebugtrackprint.ne.0) then
211                    call h_print_tracks
212                 endif                                   ! end test on zero tracks
213           1000  return
214                 end
215           
216           *