version 1.1, 1994/02/19 06:20:53
|
version 1.5, 1994/10/12 18:52:06
|
|
|
* remove minuit. Make fit linear | * remove minuit. Make fit linear |
* still does not do errors properly | * still does not do errors properly |
* $Log$ | * $Log$ |
|
* 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 | * Revision 1.1 1994/02/19 06:20:53 cdaq |
* Initial revision | * Initial revision |
* | * |
|
|
include "gen_data_structures.cmn" | include "gen_data_structures.cmn" |
include "hms_tracking.cmn" | include "hms_tracking.cmn" |
include "hms_geometry.cmn" | include "hms_geometry.cmn" |
|
external H_DPSIFUN |
* | * |
* | * |
* local variables | * local variables |
|
|
integer*4 ihit,ierr | integer*4 ihit,ierr |
integer*4 hit,plane | integer*4 hit,plane |
integer*4 i,j ! loop index | integer*4 i,j ! loop index |
|
real*8 H_DPSIFUN |
|
real*8 pos |
real*4 ray(hnum_fpray_param) | real*4 ray(hnum_fpray_param) |
|
real*8 ray1(4) |
|
real*8 ray2(4) |
real*8 TT(hnum_fpray_param) | real*8 TT(hnum_fpray_param) |
real*8 AA(hnum_fpray_param,hnum_fpray_param) | real*8 AA(hnum_fpray_param,hnum_fpray_param) |
real*8 dray(hnum_fpray_param) | real*8 dray(hnum_fpray_param) |
|
|
* | * |
ABORT= .FALSE. | ABORT= .FALSE. |
ierr=0 | ierr=0 |
|
* initailize residuals |
|
|
|
do itrack=1,HNTRACKS_MAX |
|
do plane=1,HMAX_NUM_DC_PLANES |
|
hdc_residual(itrack,plane)=1000 |
|
hdc_sing_res(itrack,plane)=1000 |
|
hdc1_sing_res(plane)=1000 |
|
hdc2_sing_res(plane)=1000 |
|
hdc1_dbl_res(plane)=1000 |
|
hdc2_dbl_res(plane)=1000 |
|
|
|
enddo |
|
enddo |
|
|
|
|
* test for no tracks | * test for no tracks |
if(HNTRACKS_FP.ge.1) then | if(HNTRACKS_FP.ge.1) then |
do itrack=1,HNTRACKS_FP | do itrack=1,HNTRACKS_FP |
chi2= dummychi2 | chi2= dummychi2 |
htrack_fit_num=itrack | htrack_fit_num=itrack |
|
|
* are there enough degrees of freedom | * are there enough degrees of freedom |
HNFREE_FP(itrack)=HNTRACK_HITS(itrack,1)-hnum_fpray_param | HNFREE_FP(itrack)=HNTRACK_HITS(itrack,1)-hnum_fpray_param |
if(HNFREE_FP(itrack).gt.0) then | if(HNFREE_FP(itrack).gt.0) then |
|
|
* initialize parameters | * initialize parameters |
do i=1,hnum_fpray_param | do i=1,hnum_fpray_param |
TT(i)=0. | TT(i)=0. |
do ihit=2,HNTRACK_HITS(itrack,1)+1 | do ihit=2,HNTRACK_HITS(itrack,1)+1 |
hit=HNTRACK_HITS(itrack,ihit) | hit=HNTRACK_HITS(itrack,ihit) |
plane=HDC_PLANE_NUM(hit) | plane=HDC_PLANE_NUM(hit) |
TT(i)=TT(i)+DFLOAT((HDC_WIRE_COORD(hit)* |
TT(i)=TT(i)+((HDC_WIRE_COORD(hit)* |
& hplane_coeff(remap(i),plane)) | & hplane_coeff(remap(i),plane)) |
& /(hdc_sigma(plane)*hdc_sigma(plane))) | & /(hdc_sigma(plane)*hdc_sigma(plane))) |
enddo | enddo |
|
|
do i=1,hnum_fpray_param | do i=1,hnum_fpray_param |
do j=1,hnum_fpray_param | do j=1,hnum_fpray_param |
AA(i,j)=0. | AA(i,j)=0. |
|
if(j.lt.i)then |
|
AA(i,j)=AA(j,i) |
|
else |
do ihit=2,HNTRACK_HITS(itrack,1)+1 | do ihit=2,HNTRACK_HITS(itrack,1)+1 |
hit=HNTRACK_HITS(itrack,ihit) | hit=HNTRACK_HITS(itrack,ihit) |
plane=HDC_PLANE_NUM(hit) | plane=HDC_PLANE_NUM(hit) |
AA(i,j)=AA(i,j) + DFLOAT( |
AA(i,j)=AA(i,j) + ( |
& hplane_coeff(remap(i),plane)*hplane_coeff(remap(j),plane) |
& hplane_coeff(remap(i),plane)*hplane_coeff(remap(j) |
& /(hdc_sigma(plane)*hdc_sigma(plane))) |
$ ,plane)/(hdc_sigma(plane)*hdc_sigma(plane))) |
enddo ! end loop on ihit | enddo ! end loop on ihit |
|
endif ! end test on j .lt. i |
enddo ! end loop on j | enddo ! end loop on j |
enddo ! end loop on i | enddo ! end loop on i |
* | * |
|
|
else | else |
* calculate chi2 | * calculate chi2 |
chi2=0. | chi2=0. |
|
|
ray(1)=dray(1) | ray(1)=dray(1) |
ray(2)=dray(2) | ray(2)=dray(2) |
ray(3)=dray(3) | ray(3)=dray(3) |
|
|
do ihit=2,HNTRACK_HITS(itrack,1)+1 | do ihit=2,HNTRACK_HITS(itrack,1)+1 |
hit=HNTRACK_HITS(itrack,ihit) | hit=HNTRACK_HITS(itrack,ihit) |
plane=HDC_PLANE_NUM(hit) | plane=HDC_PLANE_NUM(hit) |
|
|
* note chi2 is single precision | * note chi2 is single precision |
chi2=chi2+ |
hdc_sing_res(itrack,plane)= |
& (HDC_WIRE_COORD(hit) | & (HDC_WIRE_COORD(hit) |
& -hplane_coeff(remap(1),plane)*ray(1) | & -hplane_coeff(remap(1),plane)*ray(1) |
& -hplane_coeff(remap(2),plane)*ray(2) | & -hplane_coeff(remap(2),plane)*ray(2) |
& -hplane_coeff(remap(3),plane)*ray(3) | & -hplane_coeff(remap(3),plane)*ray(3) |
& -hplane_coeff(remap(4),plane)*ray(4))**2 |
& -hplane_coeff(remap(4),plane)*ray(4)) |
& /(hdc_sigma(plane)*hdc_sigma(plane)) |
chi2=chi2+((hdc_sing_res(itrack,plane))**2 |
|
& )/(hdc_sigma(plane)*hdc_sigma(plane)) |
|
*djm |
|
if(itrack.eq.1 .and. plane.ge.1 .and. plane.le.6) |
|
& hdc1_sing_res(plane) = hdc_sing_res(itrack,plane) |
|
if(itrack.eq.2 .and. plane.ge.7 .and. plane.le.12) |
|
& hdc2_sing_res(plane) = hdc_sing_res(itrack,plane) |
enddo | enddo |
endif | endif |
| |
|
|
HCHI2_FP(itrack)=chi2 | HCHI2_FP(itrack)=chi2 |
enddo ! end loop over tracks | enddo ! end loop over tracks |
endif | 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 |
|
itrack=1 |
|
ihit=2 |
|
hit=HNTRACK_HITS(itrack,ihit) |
|
plane=HDC_PLANE_NUM(hit) |
|
if (plane.le.6) then |
|
itrack=2 |
|
hit=HNTRACK_HITS(itrack,ihit) |
|
plane=HDC_PLANE_NUM(hit) |
|
if (plane.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)) |
|
|
|
itrack=1 |
|
* loop over hits in second chamber |
|
do ihit=1,HNTRACK_HITS(itrack+1,1) |
|
|
|
* calculate residual in second chamber from first chamber track |
|
hit=HNTRACK_HITS(itrack+1,ihit+1) |
|
plane=HDC_PLANE_NUM(hit) |
|
pos=H_DPSIFUN(ray1,plane) |
|
hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos |
|
* djm 8/31/94 stuff this variable into 1d array we can register |
|
hdc2_dbl_res(plane) = hdc_residual(1,plane) |
|
|
|
enddo |
|
|
|
itrack=2 |
|
* loop over hits in first chamber |
|
do ihit=1,HNTRACK_HITS(itrack-1,1) |
|
|
|
* calculate residual in first chamber from second chamber track |
|
hit=HNTRACK_HITS(itrack-1,ihit+1) |
|
plane=HDC_PLANE_NUM(hit) |
|
pos=H_DPSIFUN(ray2,plane) |
|
hdc_residual(itrack,plane)=HDC_WIRE_COORD(hit)-pos |
|
* djm 8/31/94 stuff this variable into 1d array we can register |
|
hdc1_dbl_res(plane) = hdc_residual(2,plane) |
|
|
|
enddo |
|
endif ! end plane ge 7 |
|
endif ! end plane le 6 |
|
endif ! end HNTRACKS_FP eq 2 |
|
endif ! end hsignle_stub .ne. 0 |
|
|
* test if we want to dump out trackfit results | * test if we want to dump out trackfit results |
if(hdebugtrackprint.ne.0) then | if(hdebugtrackprint.ne.0) then |
call h_print_tracks | call h_print_tracks |
|
|
end | end |
| |
* | * |
|
|