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.6 * Revision 1.5 1994/10/12 18:52:06 cdaq 13 * (DJM) Initialize some variables 14 * (SAW) Prettify indentation 15 *
16 cdaq 1.5 * Revision 1.4 1994/09/01 12:29:07 cdaq 17 * (DJM) Make registered versions of residuals 18 *
19 cdaq 1.4 * Revision 1.3 1994/08/18 02:45:53 cdaq 20 * (DM) Add calculation of residuals 21 *
22 cdaq 1.3 * Revision 1.2 1994/02/22 05:25:00 cdaq 23 * (SAW) Remove dfloat calls with floating args 24 *
25 cdaq 1.2 * Revision 1.1 1994/02/19 06:20:53 cdaq 26 * Initial revision
27 cdaq 1.1 *
28 cdaq 1.2 *
29 cdaq 1.1 implicit none 30 include "gen_data_structures.cmn" 31 include "hms_tracking.cmn" 32 include "hms_geometry.cmn"
33 cdaq 1.3 external H_DPSIFUN
34 cdaq 1.1 * 35 * 36 * local variables 37 * 38 logical ABORT 39 character50 here 40 parameter (here='H_TRACK_FIT') 41 character() err 42 integer4 itrack ! track loop index 43 integer4 ihit,ierr 44 integer4 hit,plane
45 cdaq 1.6 integer4 i,j,k ! loop index 46 real4 z_slice 47
48 cdaq 1.3 real8 H_DPSIFUN 49 real8 pos
50 cdaq 1.1 real*4 ray(hnum_fpray_param)
51 cdaq 1.3 real8 ray1(4) 52 real8 ray2(4)
53 cdaq 1.1 real*8 TT(hnum_fpray_param)
54 cdaq 1.3 real*8 AA(hnum_fpray_param,hnum_fpray_param)
55 cdaq 1.1 real8 dray(hnum_fpray_param) 56 real8 error(hnum_fpray_param) 57 real4 chi2,dummychi2 58 parameter (dummychi2 = 1.E4) 59 array to remap hplane_coeff to param number 60 integer4 remap(hnum_fpray_param) 61 data remap/5,6,3,4/ 62 save remap 63 64 ABORT= .FALSE. 65 ierr=0
66 cdaq 1.3 * initailize residuals 67 68 do itrack=1,HNTRACKS_MAX 69 do plane=1,HMAX_NUM_DC_PLANES 70 hdc_residual(itrack,plane)=1000 71 hdc_sing_res(itrack,plane)=1000
72 cdaq 1.5 hdc1_sing_res(plane)=1000 73 hdc2_sing_res(plane)=1000 74 hdc1_dbl_res(plane)=1000 75 hdc2_dbl_res(plane)=1000 76
77 cdaq 1.3 enddo 78 enddo 79 80
81 cdaq 1.1 * test for no tracks 82 if(HNTRACKS_FP.ge.1) then 83 do itrack=1,HNTRACKS_FP 84 chi2= dummychi2 85 htrack_fit_num=itrack
86 cdaq 1.3
87 cdaq 1.1 * are there enough degrees of freedom 88 HNFREE_FP(itrack)=HNTRACK_HITS(itrack,1)-hnum_fpray_param 89 if(HNFREE_FP(itrack).gt.0) then
90 cdaq 1.3
91 cdaq 1.1 * initialize parameters
92 cdaq 1.5 do i=1,hnum_fpray_param 93 TT(i)=0. 94 do ihit=2,HNTRACK_HITS(itrack,1)+1 95 hit=HNTRACK_HITS(itrack,ihit) 96 plane=HDC_PLANE_NUM(hit) 97 TT(i)=TT(i)+((HDC_WIRE_COORD(hit)* 98 & hplane_coeff(remap(i),plane)) 99 & /(hdc_sigma(plane)hdc_sigma(plane))) 100 enddo 101 enddo 102 do i=1,hnum_fpray_param 103 do j=1,hnum_fpray_param 104 AA(i,j)=0. 105 if(j.lt.i)then 106 AA(i,j)=AA(j,i) 107 else 108 do ihit=2,HNTRACK_HITS(itrack,1)+1 109 hit=HNTRACK_HITS(itrack,ihit) 110 plane=HDC_PLANE_NUM(hit) 111 AA(i,j)=AA(i,j) + ( 112 & hplane_coeff(remap(i),plane)hplane_coeff(remap(j) 113 cdaq 1.5 $ ,plane)/(hdc_sigma(plane)*hdc_sigma(plane))) 114 enddo ! end loop on ihit 115 endif ! end test on j .lt. i 116 enddo ! end loop on j 117 enddo ! end loop on i
118 cdaq 1.1 * 119 * solve four by four equations
120 cdaq 1.5 call solve_four_by_four(TT,AA,dray,ierr)
121 cdaq 1.1 *
122 cdaq 1.5 if(ierr.ne.0) then 123 dray(1)=10000. 124 dray(2)=10000. 125 dray(3)=2. 126 dray(4)=2. 127 else 128 * calculate chi2 129 chi2=0. 130 131 ray(1)=dray(1) 132 ray(2)=dray(2) 133 ray(3)=dray(3) 134 ray(4)=dray(4) 135 do ihit=2,HNTRACK_HITS(itrack,1)+1 136 hit=HNTRACK_HITS(itrack,ihit) 137 plane=HDC_PLANE_NUM(hit)
138 cdaq 1.3
139 cdaq 1.1 * note chi2 is single precision
140 cdaq 1.5 hdc_sing_res(itrack,plane)= 141 & (HDC_WIRE_COORD(hit) 142 & -hplane_coeff(remap(1),plane)ray(1) 143 & -hplane_coeff(remap(2),plane)ray(2) 144 & -hplane_coeff(remap(3),plane)ray(3) 145 & -hplane_coeff(remap(4),plane)ray(4)) 146 chi2=chi2+((hdc_sing_res(itrack,plane))*2 147 & )/(hdc_sigma(plane)hdc_sigma(plane))
148 cdaq 1.4 *djm
149 cdaq 1.5 if(itrack.eq.1 .and. plane.ge.1 .and. plane.le.6) 150 & hdc1_sing_res(plane) = hdc_sing_res(itrack,plane) 151 if(itrack.eq.2 .and. plane.ge.7 .and. plane.le.12) 152 & hdc2_sing_res(plane) = hdc_sing_res(itrack,plane) 153 enddo 154 endif 155 156 HX_FP(itrack)=ray(1) 157 HY_FP(itrack)=ray(2) 158 HZ_FP(itrack)=0. ! z=0 of tracking. 159 HXP_FP(itrack)=ray(3) 160 HYP_FP(itrack)=ray(4) 161 endif ! end test on degrees of freedom
162 cdaq 1.1 HCHI2_FP(itrack)=chi2
163 cdaq 1.5 enddo ! end loop over tracks
164 cdaq 1.6 endif 165 166 * A reasonable selection of slices is presently -80,-60,-40,-20,0,20,40 167 * ,60,80. Zero is the nominal midplane between the chambers, -80 168 * corresponds closely to the exit flange position. This slice pattern 169 * is created with hz_wild=-80., hdelta_z_wild=20., and hnum_zslice=9 170 if(hz_slice_enable.ne.0)then 171 do k=1,hnum_zslice 172 z_slice = hz_wild + (k-1)hdelta_z_wild 173 hx_fp_wild(k) = hx_fp(1) + hxp_fp(1)z_slice 174 hy_fp_wild(k) = hy_fp(1) + hyp_fp(1)*z_slice 175 enddo
176 cdaq 1.1 endif
177 cdaq 1.3 178 * calculate residuals for each chamber if in single stub mode 179 * and there were 2 tracks found one in first chanber and one in the second 180 181 if (hsingle_stub.ne.0) then 182 if (HNTRACKS_FP.eq.2) then 183 itrack=1 184 ihit=2 185 hit=HNTRACK_HITS(itrack,ihit) 186 plane=HDC_PLANE_NUM(hit) 187 if (plane.le.6) then 188 itrack=2 189 hit=HNTRACK_HITS(itrack,ihit) 190 plane=HDC_PLANE_NUM(hit) 191 if (plane.ge.7) then 192 193 * condition of above met calculating residuals 194 * assigning rays to tracks in each chamber 195 * ray1 is ray from first chamber fit 196 * ray2 is ray from second chamber fit 197
198 cdaq 1.5 ray1(1)=dble(HX_FP(1)) 199 ray1(2)=dble(HY_FP(1)) 200 ray1(3)=dble(HXP_FP(1)) 201 ray1(4)=dble(HYP_FP(1)) 202 ray2(1)=dble(HX_FP(2)) 203 ray2(2)=dble(HY_FP(2)) 204 ray2(3)=dble(HXP_FP(2)) 205 ray2(4)=dble(HYP_FP(2))
206 cdaq 1.3
207 cdaq 1.5 itrack=1
208 cdaq 1.3 * loop over hits in second chamber
209 cdaq 1.5 do ihit=1,HNTRACK_HITS(itrack+1,1)
210 cdaq 1.3 211 * calculate residual in second chamber from first chamber track
212 cdaq 1.5 hit=HNTRACK_HITS(itrack+1,ihit+1) 213 plane=HDC_PLANE_NUM(hit) 214 pos=H_DPSIFUN(ray1,plane) 215 hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos
216 cdaq 1.4 * djm 8/31/94 stuff this variable into 1d array we can register
217 cdaq 1.5 hdc2_dbl_res(plane) = hdc_residual(1,plane)
218 cdaq 1.3
219 cdaq 1.5 enddo
220 cdaq 1.3
221 cdaq 1.5 itrack=2
222 cdaq 1.3 * loop over hits in first chamber
223 cdaq 1.5 do ihit=1,HNTRACK_HITS(itrack-1,1)
224 cdaq 1.3
225 cdaq 1.4 * calculate residual in first chamber from second chamber track
226 cdaq 1.5 hit=HNTRACK_HITS(itrack-1,ihit+1) 227 plane=HDC_PLANE_NUM(hit) 228 pos=H_DPSIFUN(ray2,plane) 229 hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos
230 cdaq 1.4 * djm 8/31/94 stuff this variable into 1d array we can register
231 cdaq 1.5 hdc1_dbl_res(plane) = hdc_residual(2,plane)
232 cdaq 1.3
233 cdaq 1.5 enddo
234 cdaq 1.4 endif ! end plane ge 7 235 endif ! end plane le 6 236 endif ! end HNTRACKS_FP eq 2 237 endif ! end hsignle_stub .ne. 0
238 cdaq 1.3
239 cdaq 1.1 * test if we want to dump out trackfit results 240 if(hdebugtrackprint.ne.0) then
241 cdaq 1.5 call h_print_tracks 242 endif ! end test on zero tracks 243 1000 return
244 cdaq 1.1 end 245 246 *

Analyzer/Replay: Mark Jones, Documents: Stephen Wood