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.3 * Revision 1.2 1994/02/22 05:25:00 cdaq 13 * (SAW) Remove dfloat calls with floating args 14 *
15 cdaq 1.2 * Revision 1.1 1994/02/19 06:20:53 cdaq 16 * Initial revision
17 cdaq 1.1 *
18 cdaq 1.2 *
19 cdaq 1.1 implicit none 20 include "gen_data_structures.cmn" 21 include "hms_tracking.cmn" 22 include "hms_geometry.cmn"
23 cdaq 1.3 external H_DPSIFUN
24 cdaq 1.1 * 25 * 26 * local variables 27 * 28 logical ABORT 29 character50 here 30 parameter (here='H_TRACK_FIT') 31 character() err 32 integer4 itrack ! track loop index 33 integer4 ihit,ierr 34 integer4 hit,plane 35 integer*4 i,j ! loop index
36 cdaq 1.3 real8 H_DPSIFUN 37 real8 pos
38 cdaq 1.1 real*4 ray(hnum_fpray_param)
39 cdaq 1.3 real8 ray1(4) 40 real8 ray2(4)
41 cdaq 1.1 real*8 TT(hnum_fpray_param)
42 cdaq 1.3 real*8 AA(hnum_fpray_param,hnum_fpray_param)
43 cdaq 1.1 real8 dray(hnum_fpray_param) 44 real8 error(hnum_fpray_param) 45 real4 chi2,dummychi2 46 parameter (dummychi2 = 1.E4) 47 array to remap hplane_coeff to param number 48 integer4 remap(hnum_fpray_param) 49 data remap/5,6,3,4/ 50 save remap 51 52 ABORT= .FALSE. 53 ierr=0
54 cdaq 1.3 * initailize residuals 55 56 do itrack=1,HNTRACKS_MAX 57 do plane=1,HMAX_NUM_DC_PLANES 58 hdc_residual(itrack,plane)=1000 59 hdc_sing_res(itrack,plane)=1000 60 enddo 61 enddo 62 63
64 cdaq 1.1 * test for no tracks 65 if(HNTRACKS_FP.ge.1) then 66 do itrack=1,HNTRACKS_FP 67 chi2= dummychi2 68 htrack_fit_num=itrack
69 cdaq 1.3
70 cdaq 1.1 * are there enough degrees of freedom 71 HNFREE_FP(itrack)=HNTRACK_HITS(itrack,1)-hnum_fpray_param 72 if(HNFREE_FP(itrack).gt.0) then
73 cdaq 1.3
74 cdaq 1.1 * initialize parameters 75 do i=1,hnum_fpray_param 76 TT(i)=0. 77 do ihit=2,HNTRACK_HITS(itrack,1)+1 78 hit=HNTRACK_HITS(itrack,ihit) 79 plane=HDC_PLANE_NUM(hit)
80 cdaq 1.2 TT(i)=TT(i)+((HDC_WIRE_COORD(hit)*
81 cdaq 1.1 & hplane_coeff(remap(i),plane)) 82 & /(hdc_sigma(plane)*hdc_sigma(plane))) 83 enddo 84 enddo 85 do i=1,hnum_fpray_param 86 do j=1,hnum_fpray_param 87 AA(i,j)=0. 88 do ihit=2,HNTRACK_HITS(itrack,1)+1 89 hit=HNTRACK_HITS(itrack,ihit) 90 plane=HDC_PLANE_NUM(hit)
91 cdaq 1.2 AA(i,j)=AA(i,j) + (
92 cdaq 1.1 & hplane_coeff(remap(i),plane)hplane_coeff(remap(j),plane) 93 & /(hdc_sigma(plane)hdc_sigma(plane))) 94 enddo ! end loop on ihit 95 enddo ! end loop on j 96 enddo ! end loop on i 97 * 98 * solve four by four equations 99 call solve_four_by_four(TT,AA,dray,ierr) 100 * 101 if(ierr.ne.0) then 102 dray(1)=10000. 103 dray(2)=10000. 104 dray(3)=2. 105 dray(4)=2. 106 else 107 * calculate chi2 108 chi2=0.
109 cdaq 1.3
110 cdaq 1.1 ray(1)=dray(1) 111 ray(2)=dray(2) 112 ray(3)=dray(3) 113 ray(4)=dray(4) 114 do ihit=2,HNTRACK_HITS(itrack,1)+1 115 hit=HNTRACK_HITS(itrack,ihit) 116 plane=HDC_PLANE_NUM(hit)
117 cdaq 1.3
118 cdaq 1.1 * note chi2 is single precision
119 cdaq 1.3 hdc_sing_res(itrack,plane)=
120 cdaq 1.1 & (HDC_WIRE_COORD(hit) 121 & -hplane_coeff(remap(1),plane)ray(1) 122 & -hplane_coeff(remap(2),plane)ray(2) 123 & -hplane_coeff(remap(3),plane)*ray(3)
124 cdaq 1.3 & -hplane_coeff(remap(4),plane)ray(4)) 125 chi2=chi2+((hdc_sing_res(itrack,plane))2 126 & )/(hdc_sigma(plane)hdc_sigma(plane)) 127
128 cdaq 1.1 enddo 129 endif 130 131 HX_FP(itrack)=ray(1) 132 HY_FP(itrack)=ray(2) 133 HZ_FP(itrack)=0. ! z=0 of tracking. 134 HXP_FP(itrack)=ray(3) 135 HYP_FP(itrack)=ray(4) 136 endif ! end test on degrees of freedom 137 HCHI2_FP(itrack)=chi2 138 enddo ! end loop over tracks 139 endif
140 cdaq 1.3 141 * calculate residuals for each chamber if in single stub mode 142 * and there were 2 tracks found one in first chanber and one in the second 143 144 if (hsingle_stub.ne.0) then 145 if (HNTRACKS_FP.eq.2) then 146 itrack=1 147 ihit=2 148 hit=HNTRACK_HITS(itrack,ihit) 149 plane=HDC_PLANE_NUM(hit) 150 if (plane.le.6) then 151 itrack=2 152 hit=HNTRACK_HITS(itrack,ihit) 153 plane=HDC_PLANE_NUM(hit) 154 if (plane.ge.7) then 155 156 * condition of above met calculating residuals 157 * assigning rays to tracks in each chamber 158 * ray1 is ray from first chamber fit 159 * ray2 is ray from second chamber fit 160 161 cdaq 1.3 ray1(1)=dble(HX_FP(1)) 162 ray1(2)=dble(HY_FP(1)) 163 ray1(3)=dble(HXP_FP(1)) 164 ray1(4)=dble(HYP_FP(1)) 165 ray2(1)=dble(HX_FP(2)) 166 ray2(2)=dble(HY_FP(2)) 167 ray2(3)=dble(HXP_FP(2)) 168 ray2(4)=dble(HYP_FP(2)) 169 170 itrack=1 171 * loop over hits in second chamber 172 do ihit=1,HNTRACK_HITS(itrack+1,1) 173 174 * calculate residual in second chamber from first chamber track 175 hit=HNTRACK_HITS(itrack+1,ihit+1) 176 plane=HDC_PLANE_NUM(hit) 177 pos=H_DPSIFUN(ray1,plane) 178 hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos 179 180 enddo 181 182 cdaq 1.3 itrack=2 183 * loop over hits in first chamber 184 do ihit=1,HNTRACK_HITS(itrack-1,1) 185 186 * calculate residual in first chamber form second chamber track 187 hit=HNTRACK_HITS(itrack-1,ihit+1) 188 plane=HDC_PLANE_NUM(hit) 189 pos=H_DPSIFUN(ray2,plane) 190 hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos 191 192 enddo 193 endif 194 endif 195 endif 196 endif 197
198 cdaq 1.1 * test if we want to dump out trackfit results 199 if(hdebugtrackprint.ne.0) then 200 call h_print_tracks 201 endif ! end test on zero tracks 202 1000 return 203 end 204 205 *

Analyzer/Replay: Mark Jones, Documents: Stephen Wood