SUBROUTINE TRKTIM_RD( ITR, IPRNT ) C********************************************************************* C For a range stack track in standard RSTRAK format, find the C track time. This routine is intended to be as efficient as C possible. If there is ANY timing information available, it should C return a time. Given some timing info, it will try to use only C the best. C C Input variables: MOD ......... List of hit counters in track in C same format as ITRK of RSTRAK or C MOD_RS in RANGE_RS.CMN C NMOD ........ Number of counters in track C IPRNT ....... Print flag = 0, 1, 2 for no, C summary, and full printing C C Output variables: TIME ........ Track time, defaulted to -999. C SIGMA ....... Error in track time C NT .......... Number of times used C ICODE ....... Quality/error code C =0 no time found C >0 time found with inner layer hits C <0 time found with NO inner layer hits C For ICODE.NE.0, the absolute value is the C "quality", the score assigned to the C proto-cluster, see below. C C The routine uses the track as passed in MOD as a starting point. C It is assumed that MOD(1) is a T-counter. Only end-to-end coincidences C in layers A and outward are used in the time. These coincidences C are found using only zero-tolerence counters if possible. After C end-matching, proto-clusters (which may be as small as a single hit) C are found in layers A...11. These are ranked in quality (including C a coincidence with any single-end T-counter time as an important C virtue). The best proto-cluster is extended into the outer layers C to get the best average time. C PDM 7-OCT-1989 C********************************************************************* C 12-Dec-89: PDM -- Fixed bug: no-time branch was inside PRINT IF. C 16-Dec-89: PDM -- Changed TDCUNP window to +- 200 ns. C 6-Dec-90: TRI -- Got from TRIUMF with mods for new TD unpacking. C 17-Jan-91: PDM -- Added TRKTIM_MOD storage. C 5-Mar-91: PDM -- Modify to use 1990 TD tolerances and codes. C Append FUNCTION TT_SIG to source code. C Use tolerance-weighted statistics. C Keep 1989 code for NRUN<6600, but commented out C until TDTOLR.CMN provided. C 5-Oct-92: PDM -- Modify to get unified version for all years, C including 1989 and demux era. May have trouble C with 1992 data. C - NLAY_RS controls demux. This assumes that the C relevant arrays in RS_TRK will be handled this C way. HARDWIRED TO 15 FOR NOW. C - For 1989, combine standard tolerance codes with C Akerib bifurcation tolerances. The latter refer C to meantimes. C********************************************************************* C 5-Aug-94: PDM -- CHANGE TO TRKTIM_RD C - Abandon hope of having same version for pre- and post- C 1992. Info for demux RD is in different arrays, etc. C - NLAY_RS = 21 C - RSxxx --> RDxxx. C - Comment out all 1989 special code. C - Enlarge max cluster size from 20 to 30. C - Add constant LASTIN to be the last layer considered to C be the "inner region". This is currently set to 8. That C is, layers 2-8 (the old A and B), will be favored in C determining the time. In the multiplex version of TRKTIM, C this region was A-11. The new region does not go as deep, c but includes more counters. C 13-Jan-95 PDM: Copied from 1994 Pass1 directory. Change TRKTIM_MOD C variables from suffix TT to suffix TM to avoid confusion C with Target 2. C 95-Feb-03 MB: rename TT_SIG to TT_SIG_RD to avoid conflict with C function in TRKTIM C 95-May AK: call for tdunp is commented out to save CPU time, C assuming that it is already called by rd_trk. C 95-Jun-26 AK: add a check to avoid rare array out of bound problem C 95-Sep-11 JRS: New calling arguments: ITR = Track number (from RD_TRK) C for which to find time. C IPRNT = Same as before. C Results are now returned through the common block. C 96-Apr-24 GR: Tdunp is now called with a window corresponding to C to the FERA gate. C Added support for the array IHIT_TM in trktim_mod.cmn C********************************************************************* IMPLICIT NONE C #include #include #include #include #include C INTEGER MOD(50), NMOD, NT, ICODE, IPRNT, ITR REAL TIME, SIGMA C INTEGER I, IBEST, IC, IE, IL, IMOD(120), IS, IT, I1, I2, $ JMOD(30,120), LASTIN, NCLUST(120), NLAY_RS, $ NMAT, NMATMX, NMATIN, NMAX, $ NSTOP, NTOT, NTRY, YEAR, IHIT(2,120), IHITMOD(30,2,120) REAL ADDWIN, COWIN, EEC, ENDCON, ENDCOW, RN, $ S(120), SCLUST(120), SCORE(120), SIG(2), SMAX, SMOD(30,120), $ T(120), TCLUST(120), TMOD(30,120), TSOFAR, TT_SIG_RD, TWIN, $ WCLUST(120), WEIGHT, WMAX C DATA ENDCON, ENDCOW / 7., 20. / ! Narrow and wide half-windows for C end-to-end matching. DATA NMATMX / 120 / ! Max number of end-to-end matches. DATA COWIN / 3. / ! Half-window to match clusters. DATA TWIN / 10. / ! Half-window to match T-counter. DATA ADDWIN / 4. / ! Half-window to add outer layer DATA LASTIN / 8 / ! Last layer in "inner region" C C C Demultiplexed RS? C NLAY_RS = 21 C ICODE = 0 NT = 0 TIME = -999. NEV_TM = NEVT NRN_TM = NRUN NTRK_TM = ITR NHITS_TM = 0 TIM_TM = TIME SIG_TM = 0. c......Get list of hit counters from RD_TRK common block (JRS 9/11/95) NMOD = NHTRK(NTRK_TM) DO I=1,NMOD MOD(I)=ITRK(I,NTRK_TM) ENDDO CCC CALL TDCUNP( 'RS', 'CAL', -200., 200. ) C Call TDUNP with a window corresponding to the FERA gate. The gate C width is known to be 100ns, but there is some uncertainty in exactly C when the gate opens relative to prompt (it is ~20ns early at the C FERA input, but there are internal delays in the FERA, for example, C which complicate matters). We thus make our window a little bigger. TMIN_TM = -30. TMAX_TM = 90. c..k3 CALL TDUNP( 'RD', 'CAL', TMIN_TM, TMAX_TM ) CALL TDCUNP( 'RD', 'CAL', TMIN_TM, TMAX_TM ) C C Find end-to-end matches in layers beyond T-counter. Try two times C to find at least one match in inner layers (2-LASTIN) before continuing. C NTRY = 1 ............. Use only zero-tolerance counters. C 2 ............. Widen coincidence window. c 3 ............. include all counters, all layers. C NTRY = 0 100 NTRY = NTRY + 1 NMAT = 0 ! Total end-to-end matches NMATIN = 0 ! Matches in inner region NTOT = 0 DO 500 IC=1,NMOD IS = (MOD(IC)-1)/NLAY_RS + 1 ! Sector IL = MOD(IC) - NLAY_RS*(IS-1) ! Layer IF(IL .LT. 2) GO TO 500 ! Skip T for now. c.k3 NTOT = NTOT + RDMODHIT(1,IL,IS) + RDMODHIT(2,IL,IS) NTOT = NTOT + RDTHIT(1,IL,IS) + RDTHIT(2,IL,IS) C C On first two tries, require at least one match be found in A, B, C, or 11. C IF((NTRY.LT.3) .AND. (IL.GT.LASTIN) .AND. (NMAT.EQ.0))GO TO 100 C C Select window. C EEC = ENDCON IF(IL .GT. LASTIN)EEC = EEC + 2. ! Broader in z in outer layers IF(NTRY .EQ. 2)EEC = ENDCOW C CC IF(YEAR(NRUN) .EQ. 1989)THEN CC IF(TDTOLL(IL,IS)+TDTOLH(IL,IS) .GT. 0.1)THEN CC IF(NTRY .LT. 3)GO TO 490 ! Skip bad cntrs on first 2 tries CC EEC = EEC + 2.5*AMAX1(TDTOLL(IL,IS),TDTOLH(IL,IS)) CCC (Of the factor of 2.5, 2 is because the tolerances are for meantimes.) CC ENDIF CC ENDIF C C Loop through hits at both ends, accumulating matches. Allow multiple C use of hits and several matches per counter and sort them out later C with meantime coincidences because time differences are perhaps not as C well calibrated as meantimes. TT_SIG_RD returns the tolerance, coded to C indicate when, if ever, to use the hit. C c.k3 DO 300 I1=1,RDMODHIT(1,IL,IS) DO 300 I1=1,RDTHIT(1,IL,IS) c.k3 SIG(1) = TT_SIG_RD( RDFIDC(1,IL,IS,I1) ) c.k3 IF(SIG(1) .LT. -990.)GO TO 300 ! Bad hit, don't use c.k3 IF((SIG(1).LT.0.) .AND. (NTRY.LT.3))GO TO 300 ! Poor hit. c.k3 SIG(1) = ABS(SIG(1)) SIG(1) = 1.0 ! Temp.DV 0.5/SQRT(6.0) ! Assuming 0.5ns time step for TDC. DV. c.dv DO 200 I2=1,RDMODHIT(2,IL,IS) DO 200 I2=1,RDTHIT (2,IL,IS) c.k3 SIG(2) = TT_SIG_RD( RDFIDC(2,IL,IS,I2) ) c.k3 IF(SIG(2) .LT. -990.)GO TO 200 ! Bad hit, don't use c.k3 IF((SIG(2).LT.0.) .AND. (NTRY.LT.3))GO TO 200 ! Poor hit. c.k3 SIG(2) = ABS(SIG(2)) SIG(2) = 1.0 ! Temp.DV 0.5/SQRT(6.0) ! Assuming 0.5ns time step for TDC. DV. c.dv if(abs(rdtdtim(1,il,is,i1)-rdtdtim(2,il,is,i2)) ! eec includes if(abs(rdttdc (1,il,is,i1)-rdttdc (2,il,is,i2)) ! eec includes & .LT.EEC+SQRT(SIG(1)**2+SIG(2)**2+.000001))THEN ! 1989 toler NMAT = NMAT + 1 IF(IL .LE. LASTIN)NMATIN = NMATIN + 1 ! In inner region c.dv T(NMAT) = (RDTDTIM(1,IL,IS,I1)+RDTDTIM(2,IL,IS,I2))/2. T(NMAT) = (RDTTDC (1,IL,IS,I1)+RDTTDC (2,IL,IS,I2))/2. S(NMAT) = SQRT(SIG(1)**2+SIG(2)**2+.000001)/2. CC IF(YEAR(NRUN) .EQ. 1989)S(NMAT) = CC & S(NMAT) + AMAX1(TDTOLL(IL,IS),TDTOLH(IL,IS)) ! Crude IMOD(NMAT) = NLAY_RS*(IS-1) + IL IHIT(1,NMAT) = I1 IHIT(2,NMAT) = I2 IF(NMAT .GE. NMATMX)THEN WRITE(LLOG,510)NRUN, NEVT, IL GO TO 520 ENDIF ENDIF 200 CONTINUE 300 CONTINUE C C For tracks ending in inner region, check that all 3 tries have been C exhausted before giving up. C 490 IF((NTRY.LT.3) .AND. (IC.EQ.NMOD) .AND. (NMAT.EQ.0))GO TO 100 500 CONTINUE 510 FORMAT(' ***TRKTIM--RUN',I5,' EVENT',I6,'--MAX MATCHES BY LAY',I3) 520 CONTINUE IF(NMAT .EQ. 0)THEN IF(IPRNT.GE.1)WRITE(LLOG,530)NRUN, NEVT, $ NTOT, (MOD(1)-1)/NLAY_RS + 1, NMOD GO TO 9000 ENDIF 530 FORMAT(' ***TRKTIM--RUN',I5,' EVENT',I6, $'--NO END-TO-END MATCHES, TOTAL HITS,T-SECTOR,NMOD=',I4,2I3) IF(IPRNT .GE. 2)WRITE(LLOG,540)NRUN, NEVT, (MOD(1)-1)/NLAY_RS + 1, $ NMOD, NTRY, NMAT, NMATIN, NTOT 540 FORMAT(' TRKTIM--RUN',I5,' EVENT',I6,' T-SECTOR, NMOD=',2I3, $ ', TRIES, MATCHES, INNER MATCHES, TOTAL HITS=',3I3,I4) C C Use this list to find counter-to-counter coincidences. Seed a track C with hits from inner region only. Every such hit is a (short, squat) C cluster that will be expanded by finding coincidences. Thus the C number of clusters is NMATIN. Each hit can seed only one track, made C of hits further down the list. If no matches in inner layers, use C the same procedure on the outer layers. C NMAX = 0 ! Size of largest cluster. NSTOP = NMATIN ! Use inner layers. IF(NSTOP .EQ. 0)NSTOP = NMAT ! Use outer layers. DO 650 I1=1,NSTOP NCLUST(I1) = 1 TCLUST(I1) = T(I1)/S(I1)**2 SCLUST(I1) = T(I1)**2/S(I1)**2 WCLUST(I1) = 1./S(I1)**2 C C For later storage: C JMOD(1,I1) = IMOD(I1) TMOD(1,I1) = T(I1) SMOD(1,I1) = S(I1) IHITMOD(1,1,I1) = IHIT(1,I1) IHITMOD(1,2,I1) = IHIT(2,I1) c.dv DO 600 I2=I1+1,NSTOP DO 600 I2=1,NSTOP IF(I2.EQ.I1) GOTO 600 IF(ABS(T(I2)-T(I1)) .LE. COWIN)THEN NCLUST(I1) = NCLUST(I1) + 1 C C Keep track of the largest cluster. C IF(NCLUST(I1) .GT. NMAX)NMAX = NCLUST(I1) TCLUST(I1) = TCLUST(I1) + T(I2)/S(I2)**2 SCLUST(I1) = SCLUST(I1) + T(I2)**2/S(I2)**2 WCLUST(I1) = WCLUST(I1) + 1./S(I2)**2 JMOD(NCLUST(I1),I1) = IMOD(I2) TMOD(NCLUST(I1),I1) = T(I2) SMOD(NCLUST(I1),I1) = S(I2) IHITMOD(NCLUST(I1),1,I1) = IHIT(1,I2) IHITMOD(NCLUST(I1),2,I1) = IHIT(2,I2) ENDIF 600 CONTINUE 650 CONTINUE C C Pick best starting cluster by rating with this scoring system: C coincidence with one end of T-counter: +100 C N=hits in cluster (up to LASTIN): +10N C M=hits in cluster (over LASTIN+1): -10M C T=ABS(Time from prompt): -T/10 C S=sigma(T): -S C (Over LASTIN+1 real hits is not possible, sigma is the width of the hit C distribution used in the cluster.) C IT = (MOD(1)-1)/NLAY_RS + 1 ! Sector of T-counter hit SMAX = -999. DO 900 I=1,NSTOP SCORE(I) = 0. C C Strobe with T-counter. C TIME = TCLUST(I) / WCLUST(I) DO 800 IE=1,2 c.dv DO 700 I1=1,RDMODHIT(IE,1,IT) DO 700 I1=1,RDTHIT (IE,1,IT) c.dv IF(ABS(RDTDTIM(IE,1,IT,I1)-TIME) .LT. TWIN)THEN IF(ABS(RDTTDC (IE,1,IT,I1)-TIME) .LT. TWIN)THEN SCORE(I) = SCORE(I) + 100. GO TO 810 ENDIF 700 CONTINUE 800 CONTINUE 810 CONTINUE C C Number of hits. C SCORE(I) = SCORE(I) + 10.*MIN(NCLUST(I), LASTIN) IF(NCLUST(I) .GT. 6)SCORE(I) = $ SCORE(I) - 10.*(NCLUST(I)-(LASTIN+1)) C C Time from prompt. C SCORE(I) = SCORE(I) - ABS(TIME)/10. C C Width. C SIGMA = SQRT(AMAX1(SCLUST(I)/WCLUST(I) - TIME**2, 0.0001)) IF(NCLUST(I) .EQ. 1)SIGMA = SMOD(1,I) SCORE(I) = SCORE(I) - SIGMA IF(SCORE(I) .GT. SMAX)THEN IBEST = I SMAX = SCORE(I) ENDIF IF(IPRNT .GE. 2)WRITE(LLOG,910)I, TIME, SIGMA, NCLUST(I), $ SCORE(I) 900 CONTINUE 910 FORMAT(' TRKTIM--cluster',I3,' time,width,hits,score=', $ 2F8.2,I3,F8.2) C C Load this prototrack into final variables. C NT = NCLUST(IBEST) TIME = TCLUST(IBEST) SIGMA = SCLUST(IBEST) WEIGHT = WCLUST(IBEST) DO I=1,NT MOD_TM(I) = JMOD(I,IBEST) TMOD_TM(I) = TMOD(I,IBEST) SMOD_TM(I) = SMOD(I,IBEST) IHIT_TM(1,I) = IHITMOD(I,1,IBEST) IHIT_TM(2,I) = IHITMOD(I,2,IBEST) ENDDO TSOFAR = TIME / WEIGHT C C Avoid confusion with ICODE=0 by making minimum SMAX=1.1. C IF(SMAX .LT. 1.1)SMAX = 1.1 IF(NMATIN .EQ. 0)SMAX = -SMAX ! Tag case of no inner hits. C c now include any matches from outer layers that are in coincidence. C DO 1000 I=NSTOP+1,NMAT IF(ABS(T(I)-TSOFAR) .LT. ADDWIN+S(I))THEN IF(NT.GE.50) GO TO 1000 ! avoid array out of bound NT = NT + 1 TIME = TIME + T(I)/S(I)**2 SIGMA = SIGMA + T(I)**2/S(I)**2 WEIGHT = WEIGHT + 1./S(I)**2 MOD_TM(NT) = IMOD(I) TMOD_TM(NT) = T(I) SMOD_TM(NT) = S(I) IHIT_TM(1,NT) = IHIT(1,I) IHIT_TM(2,NT) = IHIT(2,I) IF(IPRNT.GE.2)WRITE(LLOG,1020) I,T(I),S(i),NT,TIME/WEIGHT 1020 FORMAT(' TRKTIM--addlaye',I3,' time,width,hits,ttime=', $ 2F8.2,I3,F8.2) ENDIF 1000 CONTINUE C TIME = TIME / WEIGHT C IF(NT .GT. 1)THEN C C For uncertainty in mean track time, count highest-weight hit as 1, C scale others. PDM LOG IV-188 C WMAX = 0. DO I=1,NT IF(1./SMOD_TM(I)**2 .GT. WMAX)WMAX = 1./SMOD_TM(I)**2 ENDDO RN = 0. DO I=1,NT RN = RN + 1./SMOD_TM(I)**2/WMAX ENDDO SIGMA = SQRT(AMAX1(SIGMA/WEIGHT - TIME**2, 0.0001))/SQRT(RN) ELSE ! If only one hit, use its tolerance. SIGMA = SMOD_TM(1) ENDIF NHITS_TM = NT TIM_TM = TIME SIG_TM = SIGMA IF(IPRNT .GE. 1)WRITE(LLOG,1010)NRUN, NEVT, TIME, SIGMA, NT, $ SMAX, TSOFAR 1010 FORMAT(' TRKTIM--RUN',I5,' EVENT',I6,' -- T=',F8.2,'+-',F6.2, $ I5,' HITS, INNER SCORE AND TIME=',2F8.2) C ICODE = IFIX(SMAX) 9000 CONTINUE C ICODE_TM = ICODE RETURN END REAL FUNCTION TT_SIG_RD( ICODE ) C************************************************************************* C Make decision on if, when, and how to use individual hits in TRKTIM, C based on TD code. This code contains information about the quality C of the fiducial pulse, the leading edge fit, and the flag pattern. C The routine FIDDEC unpacks these codes and assigns a timing tolerance C to the pulse. This routine returns a "sigma" to be used to weight C times in TRKTIM: C C TT_SIG_RD = SIGOOD if TOL = 0. C = TOL for usable non-zero TOL hits C = -TOL for hits not to be used until the C final (desperation) pass C = -999. for hits not to be used by TRKTIM. C C The decisions of what hits to use and when are value judgements made C by the author and are hard-wired here. C PDM 3/5/91 C************************************************************************* C Modified: C 5-Oct-92 - PDM: Skip fiducial tolerance for 1989, must be handled by C user, since Akerib tolerances apply to meantimes. C 13-Nov-92 - MMI, PDM: Change interpretation of EVTCOD. EVTCOD=4 is OK. 89 C data, no fiducials. C 95-Jun-20 AK: change the condition for HITCOD to accept double pulse C in fiducial. C 95-Jun-22 JRS: Add TD2 fiducial codes to EVTCOD list. C 95-Jul-31 JRS: Change the condition for FLGCOD to accept double pulses. C************************************************************************* IMPLICIT NONE C #include #include #include C INTEGER ICODE, EVTCOD, RUNCOD, HITCOD, FLGCOD, YEAR REAL SIGOOD, TOLER DATA SIGOOD / 0.36 / ! From Kpi2: PDM LOG IV-181,183 C TT_SIG_RD = -999. C C Unpack ICODE C CALL FIDDEC( ICODE, EVTCOD, RUNCOD, HITCOD, FLGCOD, TOLER ) C IF(TOLER .LT. 0.01)THEN ! good codes, zero tolerance TT_SIG_RD = SIGOOD GO TO 900 ENDIF TT_SIG_RD = TOLER IF(YEAR(NRUN) .EQ. 1989)GO TO 100 ! No fiducials C C To be safe, explicitly list the codes that SHOULD be used (on last pass). C 1990 value: IF(EVTCOD .GT. 5)THEN IF( (EVTCOD .EQ. IFID_ERR_AREA_TDDT) .OR. ! 13 $ (EVTCOD .EQ. IFID_ERR_AREA_TDPT) .OR. ! 14 IV-168 $ (EVTCOD .EQ. IFID_ERR_AREA_TDT0) .OR. ! 18 $ (EVTCOD .EQ. IFID_ERR_TD2_NO_ZEROS) .OR. ! 38 \ $ (EVTCOD .EQ. IFID_ERR_TD2_MAX_OVERFLOW) .OR. ! 41 > JRS $ (EVTCOD .EQ. IFID_ERR_TD2_SLOPE_BUMP) ) THEN ! 42 / 6/22/95 TT_SIG_RD = -TOLER ELSE TT_SIG_RD = -999. GO TO 900 ENDIF ENDIF IF(RUNCOD .GT. 0)THEN IF( (RUNCOD .EQ. IFID_ERR_NORMAL_BIFURCATION) .OR. ! 22 $ (RUNCOD .EQ. IFID_ERR_2NS_ZONE_SHIFT) .OR. ! 23 $ (RUNCOD .EQ. IFID_ERR_BAD_BIF_SIZE) .OR. ! 24 $ (RUNCOD .EQ. IFID_ERR_AGREE_WITH_HIGH_PEAK) .OR. ! 25 $ (RUNCOD .EQ. IFID_ERR_NO_BIF_DISAGREE_TDPT) .OR. ! 26 IV-171 $ (RUNCOD .EQ. IFID_ERR_NOPT_LAST_NOBIFNOSHIFT) .OR. ! 27 $ (RUNCOD .EQ. IFID_ERR_NOPT_LAST_NOBIFSHIFT) .OR. ! 28 $ (RUNCOD .EQ. IFID_ERR_NOPT_LAST_BIFNOSHIFT) .OR. ! 29 $ (RUNCOD .EQ. IFID_ERR_NOPT_LAST_BIFSHIFT) .OR. ! 30 $ (RUNCOD .EQ. IFID_ERR_NO_GOOD_PEAK_INFO) .OR. ! 32 $ (RUNCOD .EQ. IFID_ERR_BIF_NO_TDPT_RUN) )THEN ! 33 TT_SIG_RD = -TOLER ELSE TT_SIG_RD = -999. GO TO 900 ENDIF ENDIF 100 CONTINUE C C For HITCOD (leading-edge-finder), use good or no-rise codes on all passes. C Do not use the rest. IV-184 C Double pulse cases are accepted as well. June 20, 1995 AK C CC IF(HITCOD .GT. 1)THEN IF(HITCOD .GT. 3)THEN TT_SIG_RD = -999. GO TO 900 ENDIF C C For flag codes, use good or normal-sector-crossing codes on all passes C (FLGCOD=1 has a tolerance). Use the rest on pass 3. IV-183 C Multiple leading edge times in a given pulse also accepted. JRS 7/31/95 C CC IF(FLGCOD .GT. 1)THEN IF(FLGCOD .GT. 1 .AND. FLGCOD.NE.10)THEN TT_SIG_RD = -TOLER ENDIF C 900 CONTINUE RETURN END