Analyzer/HTRACKING/h_track_fit.f - annotate

Return to h_track_fit.f CVS log

Up to [HallC] / Analyzer / HTRACKING

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.11 * Revision 1.10 1995/08/30 16:11:39 cdaq 13 * (JRA) Don't fill single_residual arrray 14 *
15 cdaq 1.10 * Revision 1.9 1995/05/22 19:39:31 cdaq 16 * (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts" 17 *
18 cdaq 1.9 * Revision 1.8 1995/04/06 19:32:58 cdaq 19 * (JRA) Rename residuals variables 20 *
21 cdaq 1.8 * Revision 1.7 1995/01/27 20:26:50 cdaq 22 * (JRA) Remove Mack's personal focalplane diceamatic (z slicer) code 23 *
24 cdaq 1.7 * Revision 1.6 1994/12/06 15:45:27 cdaq 25 * (DJM) Take slices in Z to look for best focus 26 *
27 cdaq 1.6 * Revision 1.5 1994/10/12 18:52:06 cdaq 28 * (DJM) Initialize some variables 29 * (SAW) Prettify indentation 30 *
31 cdaq 1.5 * Revision 1.4 1994/09/01 12:29:07 cdaq 32 * (DJM) Make registered versions of residuals 33 *
34 cdaq 1.4 * Revision 1.3 1994/08/18 02:45:53 cdaq 35 * (DM) Add calculation of residuals 36 *
37 cdaq 1.3 * Revision 1.2 1994/02/22 05:25:00 cdaq 38 * (SAW) Remove dfloat calls with floating args 39 *
40 cdaq 1.2 * Revision 1.1 1994/02/19 06:20:53 cdaq 41 * Initial revision
42 cdaq 1.1 *
43 cdaq 1.2 *
44 cdaq 1.1 implicit none
45 cdaq 1.9 include "hms_data_structures.cmn"
46 cdaq 1.1 include "hms_tracking.cmn" 47 include "hms_geometry.cmn"
48 cdaq 1.11 external h_dpsifun
49 cdaq 1.1 * 50 * 51 * local variables 52 * 53 logical ABORT
54 cdaq 1.11 character*11 here
55 cdaq 1.1 parameter (here='H_TRACK_FIT') 56 character() err
57 cdaq 1.11 integer*4 itrk ! track loop index
58 cdaq 1.1 integer*4 ihit,ierr
59 cdaq 1.11 integer*4 hit,pln
60 cdaq 1.8 integer*4 i,j ! loop index
61 cdaq 1.7 * real*4 z_slice
62 cdaq 1.6
63 cdaq 1.11 real*8 h_dpsifun
64 cdaq 1.3 real8 pos 65 real8 ray1(4) 66 real*8 ray2(4)
67 cdaq 1.1 real*8 TT(hnum_fpray_param)
68 cdaq 1.3 real*8 AA(hnum_fpray_param,hnum_fpray_param)
69 cdaq 1.1 real8 dray(hnum_fpray_param) 70 real4 chi2,dummychi2 71 parameter (dummychi2 = 1.E4) 72 * array to remap hplane_coeff to param number 73 integer4 remap(hnum_fpray_param) 74 data remap/5,6,3,4/ 75 save remap 76 77 ABORT= .FALSE. 78 ierr=0
79 cdaq 1.3 * initailize residuals 80
81 cdaq 1.11 do pln=1,hdc_num_planes 82 do itrk=1,hntracks_fp 83 hdc_double_residual(itrk,pln)=1000 84 hdc_single_residual(itrk,pln)=1000
85 cdaq 1.3 enddo
86 cdaq 1.10 c fill the 1d arrays from the 2d arrays for good track (in h_physics)
87 cdaq 1.11 c hdc_sing_res(pln)=1000 88 hdc_dbl_res(pln)=1000
89 cdaq 1.3 enddo 90
91 cdaq 1.1 * test for no tracks
92 cdaq 1.11 if(hntracks_fp.ge.1) then 93 do itrk=1,hntracks_fp
94 cdaq 1.1 chi2= dummychi2
95 cdaq 1.11 htrack_fit_num=itrk
96 cdaq 1.3
97 cdaq 1.1 * are there enough degrees of freedom
98 cdaq 1.11 hnfree_fp(itrk)=hntrack_hits(itrk,1)-hnum_fpray_param 99 if(hnfree_fp(itrk).gt.0) then
100 cdaq 1.3
101 cdaq 1.1 * initialize parameters
102 cdaq 1.5 do i=1,hnum_fpray_param 103 TT(i)=0.
104 cdaq 1.11 do ihit=2,hntrack_hits(itrk,1)+1 105 hit=hntrack_hits(itrk,ihit) 106 pln=hdc_plane_num(hit) 107 TT(i)=TT(i)+((hdc_wire_coord(hit)* 108 & hplane_coeff(remap(i),pln)) 109 & /(hdc_sigma(pln)*hdc_sigma(pln)))
110 cdaq 1.5 enddo 111 enddo 112 do i=1,hnum_fpray_param 113 do j=1,hnum_fpray_param 114 AA(i,j)=0. 115 if(j.lt.i)then 116 AA(i,j)=AA(j,i) 117 else
118 cdaq 1.11 do ihit=2,hntrack_hits(itrk,1)+1 119 hit=hntrack_hits(itrk,ihit) 120 pln=hdc_plane_num(hit)
121 cdaq 1.5 AA(i,j)=AA(i,j) + (
122 cdaq 1.11 & hplane_coeff(remap(i),pln)hplane_coeff(remap(j) 123 $ ,pln)/(hdc_sigma(pln)hdc_sigma(pln)))
124 cdaq 1.5 enddo ! end loop on ihit 125 endif ! end test on j .lt. i 126 enddo ! end loop on j 127 enddo ! end loop on i
128 cdaq 1.1 * 129 * solve four by four equations
130 cdaq 1.5 call solve_four_by_four(TT,AA,dray,ierr)
131 cdaq 1.1 *
132 cdaq 1.5 if(ierr.ne.0) then 133 dray(1)=10000. 134 dray(2)=10000. 135 dray(3)=2. 136 dray(4)=2. 137 else 138 * calculate chi2 139 chi2=0. 140
141 cdaq 1.11 * calculate hit coord at each plane for chisquared and efficiency calculations. 142 do pln=1,hdc_num_planes 143 hdc_track_coord(itrk,pln)=hplane_coeff(remap(1),pln)dray(1) 144 & +hplane_coeff(remap(2),pln)dray(2) 145 & +hplane_coeff(remap(3),pln)dray(3) 146 & +hplane_coeff(remap(4),pln)dray(4) 147 enddo 148 149 do ihit=2,hntrack_hits(itrk,1)+1 150 hit=hntrack_hits(itrk,ihit) 151 pln=hdc_plane_num(hit)
152 cdaq 1.3
153 cdaq 1.1 * note chi2 is single precision
154 cdaq 1.11 155 hdc_single_residual(itrk,pln)= 156 & hdc_wire_coord(hit)-hdc_track_coord(itrk,pln) 157 chi2=chi2+ 158 & (hdc_single_residual(itrk,pln)/hdc_sigma(pln))**2
159 cdaq 1.5 enddo 160 endif 161
162 cdaq 1.11 hx_fp(itrk)=dray(1) 163 hy_fp(itrk)=dray(2) 164 hz_fp(itrk)=0. ! z=0 of tracking. 165 hxp_fp(itrk)=dray(3) 166 hyp_fp(itrk)=dray(4)
167 cdaq 1.5 endif ! end test on degrees of freedom
168 cdaq 1.11 hchi2_fp(itrk)=chi2
169 cdaq 1.5 enddo ! end loop over tracks
170 cdaq 1.6 endif 171
172 cdaq 1.3 * calculate residuals for each chamber if in single stub mode 173 * and there were 2 tracks found one in first chanber and one in the second 174 175 if (hsingle_stub.ne.0) then
176 cdaq 1.11 if (hntracks_fp.eq.2) then 177 itrk=1
178 cdaq 1.3 ihit=2
179 cdaq 1.11 hit=hntrack_hits(itrk,ihit) 180 pln=hdc_plane_num(hit) 181 if (pln.le.6) then 182 itrk=2 183 hit=hntrack_hits(itrk,ihit) 184 pln=hdc_plane_num(hit) 185 if (pln.ge.7) then
186 cdaq 1.3 187 * condition of above met calculating residuals 188 * assigning rays to tracks in each chamber 189 * ray1 is ray from first chamber fit 190 * ray2 is ray from second chamber fit 191
192 cdaq 1.11 ray1(1)=dble(hx_fp(1)) 193 ray1(2)=dble(hy_fp(1)) 194 ray1(3)=dble(hxp_fp(1)) 195 ray1(4)=dble(hyp_fp(1)) 196 ray2(1)=dble(hx_fp(2)) 197 ray2(2)=dble(hy_fp(2)) 198 ray2(3)=dble(hxp_fp(2)) 199 ray2(4)=dble(hyp_fp(2))
200 cdaq 1.3
201 cdaq 1.11 itrk=1
202 cdaq 1.3 * loop over hits in second chamber
203 cdaq 1.11 do ihit=1,hntrack_hits(itrk+1,1)
204 cdaq 1.3 205 * calculate residual in second chamber from first chamber track
206 cdaq 1.11 hit=hntrack_hits(itrk+1,ihit+1) 207 pln=hdc_plane_num(hit) 208 pos=h_dpsifun(ray1,pln) 209 hdc_double_residual(itrk,pln)=hdc_wire_coord(hit)-pos
210 cdaq 1.4 * djm 8/31/94 stuff this variable into 1d array we can register
211 cdaq 1.11 hdc_dbl_res(pln) = hdc_double_residual(1,pln)
212 cdaq 1.3
213 cdaq 1.5 enddo
214 cdaq 1.3
215 cdaq 1.11 itrk=2
216 cdaq 1.3 * loop over hits in first chamber
217 cdaq 1.11 do ihit=1,hntrack_hits(itrk-1,1)
218 cdaq 1.3
219 cdaq 1.4 * calculate residual in first chamber from second chamber track
220 cdaq 1.11 hit=hntrack_hits(itrk-1,ihit+1) 221 pln=hdc_plane_num(hit) 222 pos=h_dpsifun(ray2,pln) 223 hdc_double_residual(itrk,pln)=hdc_wire_coord(hit)-pos
224 cdaq 1.4 * djm 8/31/94 stuff this variable into 1d array we can register
225 cdaq 1.11 hdc_dbl_res(pln) = hdc_double_residual(2,pln)
226 cdaq 1.3
227 cdaq 1.5 enddo
228 cdaq 1.11 endif ! end pln ge 7 229 endif ! end pln le 6 230 endif ! end hntracks_fp eq 2
231 cdaq 1.4 endif ! end hsignle_stub .ne. 0
232 cdaq 1.3
233 cdaq 1.1 * test if we want to dump out trackfit results 234 if(hdebugtrackprint.ne.0) then
235 cdaq 1.5 call h_print_tracks 236 endif ! end test on zero tracks 237 1000 return
238 cdaq 1.1 end 239 240 *

Analyzer/Replay: Mark Jones, Documents: Stephen Wood