*file member=ft1_fit_event library=snoman language=fortran77
*file date=24:Nov:1994
subroutine ft1_fit_event(iretc)
:
:
:
HEADER
:
:
:
* Local Variable Declarations:-
* ===========================
integer I_FIT_MASK_BIT
parameter (I_FIT_MASK_BIT = KPF_FT_MASK + KFT_USER1)
integer num_pmt_avail , num_pmt_used ,
+ i_pmt_dr_max , i_pmt , i_pmt_id
real sum_x , sum_y , sum_z , sum_t , sum_dr2 ,
+ fit_xyz(3) , fit_r , dr , dr_max , rms ,
+ vdist
--
iretc = KSU_NO_DATA |
|
*** If no EV quit with KSU_NO_DATA. |
|
if (lev .le. 0) return |
|
*** Create FT1, also creating FT if necessary, but quit if not |
*** not enough data to start fit. |
|
lft = lq(lev - KEV_FT) 1
if (lft .eq. 0) |
+ call mzlift( IDIV_EVENT , lft , lev , -KEV_FT , mmft , 0 ) |
call mzlift( IDIV_EVENT , lft1 , lft , -KFT_FT1 , mmft1 , 0 ) |
iq(lft1 + KFTX_METHOD) = KFT_USER1 |
iq(lft1 + KFTX_RETC) = iretc |
|
lpf = lq(lev - KEV_PF) |
lpn = lq(lev - KEV_PN) |
lpt = lq(lev - KEV_PT) |
if (lpf .eq. 0 .or. lpn .eq. 0 .or. lpt .eq. 0) return |
--
num_pmt_avail = 0 --
sum_x = 0. |
sum_y = 0. 2
sum_z = 0. |
sum_t = 0. |
--
do i_pmt = 1 , iq(lpf-1) --
* The FT1 mask bit ought not to be set, but just in case... |
iq(lpf + i_pmt) = ibclr( iq(lpf + i_pmt) , I_FIT_MASK_BIT ) |
if (.not. btest( iq(lpf + i_pmt) , KPF_DIS ) ) then |
num_pmt_avail = num_pmt_avail + 1 |
i_pmt_id = iq(lpn + i_pmt) |
sum_x = sum_x + rge_pmt_xyz(1,i_pmt_id) |
sum_y = sum_y + rge_pmt_xyz(2,i_pmt_id) 3
sum_z = sum_z + rge_pmt_xyz(3,i_pmt_id) |
sum_t = sum_t + rq(lpt + i_pmt) |
endif |
enddo |
|
if (num_pmt_avail .lt. icons(ldtft1 + KTFT1_MIN_PMT) ) return |
--
--
*** Fit vertex and iterate dropping points until RMS acceptable |
*** or too many points dropped. |
|
num_pmt_used = num_pmt_avail |
iretc = KSU_NOOP |
|
do while (iretc .eq. KSU_NOOP) |
|
** Compute fit vertex as the C of G. of all PMT hits. |
|
fit_xyz(1) = sum_x / num_pmt_used |
fit_xyz(2) = sum_y / num_pmt_used |
fit_xyz(3) = sum_z / num_pmt_used 4
|
** Find mean radial distance from fit to each hit. |
|
fit_r = 0. |
do i_pmt = 1 , iq(lpf-1) |
if ( .not. btest( iq(lpf + i_pmt) , KPF_DIS ) |
+ .and. .not. btest( iq(lpf + i_pmt) , I_FIT_MASK_BIT ) |
+ ) then |
i_pmt_id = iq(lpn + i_pmt) |
fit_r= fit_r + vdist(fit_xyz,rge_pmt_xyz(1,i_pmt_id),3) |
endif |
enddo |
fit_r = fit_r / num_pmt_used |
|
** Find fit RMS and worst residual. |
|
sum_dr2 = 0. |
dr_max = 0. |
do i_pmt = 1 , iq(lpf-1) |
if ( .not. btest( iq(lpf + i_pmt) , KPF_DIS ) |
+ .and. .not. btest( iq(lpf + i_pmt) , I_FIT_MASK_BIT ) |
+ ) then |
i_pmt_id = iq(lpn + i_pmt) |
dr = vdist(fit_xyz,rge_pmt_xyz(1,i_pmt_id),3) - fit_r |
sum_dr2 = sum_dr2 + dr**2 |
if (abs(dr) .gt. dr_max) then |
dr_max = abs(dr) |
i_pmt_dr_max = i_pmt |
endif |
endif |
enddo |
rms = sqrt( sum_dr2 / num_pmt_used ) |
|
** If fit RMS O.K., accept fit. 4
|
if (RMS .le. rcons(ldtft1 + KTFT1_MAX_RMS) ) then |
iretc = KSU_OK |
|
** If enough hits left, drop worst and continue. |
|
elseif (num_pmt_used .gt. icons(ldtft1+KTFT1_MIN_PMT)) then |
i_pmt = i_pmt_dr_max |
i_pmt_id = iq(lpn + i_pmt) |
num_pmt_used = num_pmt_used - 1 |
sum_x = sum_x - rge_pmt_xyz(1,i_pmt_id) |
sum_y = sum_y - rge_pmt_xyz(2,i_pmt_id) |
sum_z = sum_z - rge_pmt_xyz(3,i_pmt_id) |
sum_t = sum_t - rq(lpt + i_pmt) |
iq(lpf + i_pmt) = ibset(iq(lpf + i_pmt), I_FIT_MASK_BIT) |
|
** Otherwise fail. |
|
else |
iretc = KSU_FAILED |
|
endif |
|
enddo --
--
iq(lft1 + KFTX_RETC) = iretc |
if (iretc .ne. KSU_OK) return |
|
*** Fit O.K.; complete FT1, store results in an FT1V and FT1T and |
*** point currently selected FTx to FT1. |
|
|
iq(lft1 + KFTX_PMT_AVAIL) = num_pmt_avail |
iq(lft1 + KFTX_PMT_USED) = num_pmt_used |
iq(lft1 + KFTX_ITER) = num_pmt_avail - num_pmt_used |
* At least this word is meaningful! |
rq(lft1 + KFTX_PROB) = 0.0 |
* Store data in fitter specific extension. 5
rq(lft1 + KFT1_RMS) = RMS |
|
call mzlift( IDIV_EVENT, lft1v, lft1, -KFTX_FTXV, mmft1v, 0 ) |
rq(lft1v + KFTXV_X) = sum_x / num_pmt_used |
rq(lft1v + KFTXV_Y) = sum_y / num_pmt_used |
rq(lft1v + KFTXV_Z) = sum_z / num_pmt_used |
rq(lft1v + KFTXV_T) = sum_t / num_pmt_used |
|
call mzlift(IDIV_EVENT, lft1t, lft1v, -KFTXV_FTXT, mmft1t, 0) |
rq(lft1t + KFTXT_U) = 0.57735 |
rq(lft1t + KFTXT_V) = 0.57735 |
rq(lft1t + KFTXT_W) = 0.57735 |
|
lq(lft - KFT_FTX_CUR) = lft1 |
|
return |
end |
--
*endfile member=ft1_fit_event
To facilitate explanation the code above has been split into 5 parts as indicated by the bracketing.
PART 1: Here a few checks are made, pointers set up and the FT1 bank is lifted. The user should not change anything.
PART 2: This is just some initialisation specific to this particular fitting method. It should be removed.
PART 3: This part comprises a loop over the hits (iq(lpf-1) holds
the number of hits in the trigger as opposed to the number of hit PMTs).
Users may find it useful to keep and modify this section, iq(lpn + N)
holds the PMT ID of the Nth hit and the array rge_pmt_xyz holds the XYZ
coordinates of all PMTs, with the time of the Nth hit being held in
rq(lpt+ N). The variables sum_x, sum_y etc. are specific to this
particular fitting method and can be discarded. The PF bank contains a set of
flags for each hit and it is possible to set one of these flags so that the hit
is not used by a particular fitter or to set another flag to discard it
entirely for all fits. The second line within the loop carries out a check to
see if a hit is to be discarded generally and the line above it clears all the
PF flags that indicate discarding for this particular fitter. Notice also the
check on the number of hits available at the end of this part using the titles
bank address mnemonic set up in FT1_COM.INC.
PART 4: Here is the real mathematics of the fit. To describe this particular fitting method Nick West's comments are appropriate:
``Calling this fitter Micky Mouse would be deeply insulting to Walt Disney for there is no justification for the method, which attempts to find a point equidistant between all hit PMTs (ignoring time differences!), throwing away points with a large RMS to the fitted radial distance. However the fitter does demonstrate all the features that a real fitter must have.''This part should be ripped out entirely and replaced by the fitting method of the user's choice with calls to subroutines if required and access to the PMT data having been already obtained in the loop of part 3. Two things should be noted, however. Firstly the setting of the bit flag towards the end of part 4 to discard a hit may be useful and secondly the call argument
iretc should be set appropriately before returning to the calling
routine FT1_FITTER. This is to enable run statistics to be printed by
FT1_FITTER_TRM. The available values for iretc, which are defined in
SU_MNEMONICS.INC, are:
KSU_OK |
Fit OK |
KSU_NOOP |
No operation (output already exists) |
KSU_UNSUPPORTED |
Unsupported operation |
KSU_NO_DATA |
Not enough data to fit |
KSU_TOO_MUCH_DATA |
Too much data to fit |
KSU_FAILED |
Fit failed |
PART 5: In this final part the output from the fitter is loaded into the data
structure. See the FTX, FTXV and FTXT documentation refered to in the
BANK_FT1.INC section for the predefined slots and mnemonics. It is here also
that the FT1 specific fit bank mnemonics that were defined in BANK_FT1.INC are
used. Remember also to store each data word in the appropriate array iq
or rq depending on whether it is integer or floating point. If a fit
vertex or fit track bank is required it should be lifted with an MZLIFT
call. As can be seen in the code above both are currently lifted and the call
arguments do not need to be changed. The final line assigns to currently
selected fitter to FT1 and should not be changed.