*-- Author : R. Henderson
SUBROUTINE FPDG32(IP)
**: FPDG32 40000 RH. New steering parameter.
*------------------------------------------------------------------
**: FPDG32 30907 RH. Bug fix in cluster finding.
C------------------------------------------------------------------
C
C--- This routine finds clusters from 3 digitizings
C--- wires 1 2 3 only
C--- This routine should always be called from within FPDG4
C--- after FPDG31
C
C---
*KEEP,FPPRAM.
C
C--- MAXSEG is maximum number of segments per supermodule
C--- MAXCON is maximum number of amibiguous segments associatable with
C--- one segment
C--- LIMSTO is maximum number of 2 cluster planes intersections to be
C--- stored per supermodule
C--- MSEGLM is maximum number of clusters that can be found before
C--- connectivity considered
C--- MAXCLU is maximum number of clusters that can be found after
C--- forming non-connected set MUST BE 50 IF RUN WITH OLD RCW
C--- (cluster = 3/4 digits found in a straight line in one
C--- 4-wire orientation)
C
PARAMETER (MAXSEG = 200)
PARAMETER (MAXCON = 100)
PARAMETER (LIMSTO = 5000)
PARAMETER (MSEGLM = 150)
PARAMETER (MAXCLU = 50)
C---
*KEND.
C---
*KEEP,FPLGEO.
C---
COMMON /FPLGEO/ ZPLAN(36) , TP(9) , YP(26) , PLANE(3,9),
1 RMAX , RMIN , YSTART , YSPACE ,
2 X0 , Y0 , PZSTRU (8), STAGER ,
3 RESOL , ACUT , CTP(9) , STP(9)
C---
*KEND.
C---
*KEEP,FPCLUS.
COMMON /FPCLUS/ TC(3,9,MAXCLU) , NTC(9) , TOC(3,9,MAXCLU) ,
2 TCYUV(4,9,MAXCLU), TCYUVW(4,9,MAXCLU)
C---
*KEND.
C---
*KEEP,FPH1WRK.
PARAMETER (MAXHTS=200)
PARAMETER (NUMWPL=36)
PARAMETER (MAXTRK=200)
PARAMETER (MXTTRK=900)
PARAMETER (MAXTR3=200)
PARAMETER (MAXHPW=2)
PARAMETER (MAXDIG=2000)
PARAMETER (NUMRWR=1727)
PARAMETER (NUMPWR=1151)
LOGICAL DRMASK
COMMON /H1WORK/
C-- *KEEP,FPCSEG.
C---
3 TPNORM(3,9,MAXCLU), PCONST(9,MAXCLU) ,
4 SMLS(4,2,LIMSTO,3) ,
C---
C-- *KEEP,FPDIGI.
5 DRSTO(MSEGLM,4),NDRSTO(4),
6 IDIGST(4,MSEGLM),
7 SEGTAB(MSEGLM,MSEGLM),DRMASK(MSEGLM,4),
8 IDCELL(MSEGLM,4),
9 NSGTAB(MSEGLM),ASGTAB(MSEGLM),
A RESSTO(MSEGLM,4) ,
C---
C-- *KEEP,FPDGI.
B IDGIST(MSEGLM,4),IDGISM(4,9,MAXCLU)
C ,RCHI(MAXSEG,3) ,
C---
C-- *KEEP,FPSTID.
D IDRSTO(MSEGLM,4),IDYUV(4,9,MSEGLM),
E IDYUVS(12,MAXSEG,3),FREQ(MSEGLM+1) ,
C---
C-- *interface to real data
F NDPW(NUMWPL),DW(MAXHTS,NUMWPL),DWG(MAXHTS,NUMWPL),
G DRIWP(MAXHTS,NUMWPL),DRIWM(MAXHTS,NUMWPL),
G DRIW(MAXHTS,NUMWPL),IPHOLE(MAXHTS,NUMWPL),
H IPFRPE(MAXHTS,36), IPPIOS(MAXDIG,2)
C---.
*KEND.
C---
DIMENSION COVSLZ(2,2)
C---
LOGICAL LGKS,LSPLIT
REAL LINEX(2),LINEY(2),LOOP(50),W(4),PROBL(50)
DOUBLE PRECISION Y(4),ONE
C-------------------------------------------------------------------
C DIGITS ANALYSIS
C
C The following arrays are used in the remainder of this routine
C
C drsto(drift,wire_in_orientation):
C This has both drift and reflection per orientation
C
C ndrsto(wire_in_orientation):
C Number of drifts + reflections on wire (Naughty)
C
C drmask(drift,wire_in_orientation):
C 1.0 If digit or reflection used
C 0.0 If digit or reflection unused
C
C idigst(wire_in_orientation,found_segment_candidate):
C index of digitising associated with a particular
C candidate segment for a particular wire
C
C nsgtab(candidate_segment):
C Number of candidate segments with which current segment
C shares digitizings
C
C segtab(candidate_segment,associated_ambiguous_segments):
C segment index ambiguous with current segment
C
C-------------------------------------------------------------------
C
C--- tolerance for finding digitizings
C
TOLER = RESOL*3.0
C
C--- define maximum slope to be found for segment
C
SLMAX = 3.0
SLYMIN = 0.005
SLYMAX = 0.015
C
C--- Section to find 2d line segments
C
C--- loop over first and last wire in orientation
C--- but allow one of the intermediate points to be missing
C
C
C--- entry point for more relaxed segment search
C--- ifirst = 1 on first pass through
C
IFIRST = 1
1100 CONTINUE
NSEG = 0
DO 200 IO1 = 1 , NDRSTO(1)
IF(DRMASK(IO1,1)) GO TO 200
C
C--- determine if dealing with a split cell
C
IF( IDCELL(IO1,1) .EQ. 1 .OR.
1 IDCELL(IO1,1) .EQ. -1 ) THEN
LSPLIT = .TRUE.
ELSE
LSPLIT = .FALSE.
ENDIF
DO 201 IO3 = 1 , NDRSTO(3)
IF(DRMASK(IO3,3)) GO TO 201
IF( LSPLIT .AND. IDCELL(IO1,1) .NE. IDCELL(IO3,3)
1 .AND. IDCELL(IO3,3) .NE. 0 )
1 GO TO 201
LINEX(1) = 9.0
LINEY(1) = DRSTO(IO1,1)
LINEX(2) = ( ( ZPLAN(IP+2) - ZPLAN(IP) ) / 10.0 )
1 + 9.0
LINEY(2) = DRSTO(IO3,3)
GRAD = (LINEY(2) - LINEY(1)) / (LINEX(2) - LINEX(1))
C
C--- On first pass filter out large slopes
C
IF ( ABS(GRAD) .GT. SLMAX ) GO TO 201
C---
IF ( (ABS(LINEY(1)) .GT. 100.0 .AND.
1 GRAD*LINEY(1) .LT. 0.0 ).OR.
2 ABS(GRAD) .GT. ABS(LINEY(1))*SLYMAX .OR.
3 ABS(GRAD) .LT. ABS(LINEY(1))*SLYMIN ) GO TO 201
C
C--- use tolerance to find digitizings to form segments
C
C
C--- plane 2
C
DO 202 IO2 = 1 , NDRSTO(2)
IF(DRMASK(IO2,2)) GO TO 202
IF( LSPLIT .AND. IDCELL(IO1,1) .NE. IDCELL(IO2,2)
1 .AND. IDCELL(IO2,2) .NE. 0 )
1 GO TO 202
PRED2 = LINEY(1) + GRAD*
1 ( ( ZPLAN(IP+1) - ZPLAN(IP) ) / 10.0 )
IF( ABS( PRED2 - DRSTO(IO2,2)) .GT. TOLER)GO TO 202
C
C--- store the digitizing per segment found
C
NSEG = NSEG + 1
IF (NSEG .GT. MSEGLM) THEN
CALL ERRLOG(207,'W:FPDG32: NSEG > MSEGLM')
NSEG = NSEG - 1
GO TO 205
ENDIF
IDIGST(1,NSEG) = IO1
IDIGST(2,NSEG) = IO2
IDIGST(3,NSEG) = IO3
IDIGST(4,NSEG) = -1
202 CONTINUE
201 CONTINUE
200 CONTINUE
205 CONTINUE
C
C--- Now sort out which initial segments to keep
C
C
C--- Create segtab showing interconnectivity
C
CALL VZERO(NSGTAB,MSEGLM)
DO 300 ISEG = 1,NSEG
IF(NSGTAB(ISEG) .EQ. MSEGLM)GO TO 300
C
C--- loop over remaining segments
C
DO 302 KSEG = ISEG+1,NSEG
IF(NSGTAB(KSEG) .EQ. MSEGLM)GO TO 302
C
C--- comparison loop over each wire in turn
C
DO 301 ID = 1,4
ID1 = IDIGST(ID,ISEG)
IF(ID1 .EQ. -1) GO TO 301
IF ( MOD(ID1,2) .EQ. 0 )THEN
ID2 = ID1 - 1
ELSE
ID2 = ID1 + 1
ENDIF
C
C--- has the same wire the same digit
C
IF( IDIGST(ID,KSEG) .NE. ID1 .AND.
1 IDIGST(ID,KSEG) .NE. ID2) GO TO 301
NSGTAB(ISEG) = NSGTAB(ISEG) + 1
SEGTAB(ISEG,NSGTAB(ISEG)) = KSEG
NSGTAB(KSEG) = NSGTAB(KSEG) + 1
SEGTAB(KSEG,NSGTAB(KSEG)) = ISEG
GO TO 302
301 CONTINUE
302 CONTINUE
300 CONTINUE
C
C--- remove nodes greater than or equal to 3
C
400 CONTINUE
IF( NSEG .LT. 1)GO TO 500
CALL VFLOAT(NSGTAB,ASGTAB,NSEG)
MXSEG = LVMAX(ASGTAB,NSEG)
VMXSEG = ASGTAB(MXSEG)
IF ( VMXSEG .LE. 2.0 ) GO TO 500
C
C--- greater than 2.0 so remove by setting nsgtab = -1.0
C
NSGTAB(MXSEG) = -1.0
C
C--- Now remove any reference to this node in the remaining nodes
C
DO 401 ISEG = 1, NSEG
IF ( NSGTAB(ISEG) .EQ. -1) GO TO 401
DO 402 ID = 1,NSGTAB(ISEG)
IF( SEGTAB(ISEG,ID) .NE. MXSEG ) GO TO 402
SEGTAB(ISEG,ID) = SEGTAB(ISEG,NSGTAB(ISEG))
NSGTAB(ISEG) = NSGTAB(ISEG) - 1
402 CONTINUE
401 CONTINUE
GO TO 400
500 CONTINUE
C
C--- Now try to find loops and angles and eliminate
C
DO 700 ISEG = 1,NSEG
IF( NSGTAB(ISEG) .EQ. -1 ) GO TO 700
VALSEG = NSGTAB(ISEG)
C
C--- Find first candidate with 2 links
C
IF ( VALSEG .LT. 2.0 ) GO TO 700
C
C--- Now trace its path
C
ILOOP = 1
LOOP(ILOOP) = ISEG
LSTSEG = ISEG
NXTSEG = SEGTAB(ISEG,1)
C
C--- Entry point for step along chain
C
703 CONTINUE
C
C--- Test if path at end
C
IF( NSGTAB(NXTSEG) .LT. 2) GO TO 701
C
C--- Skip link if pointing back
C
NEWSEG = SEGTAB(NXTSEG,1)
IF ( NEWSEG .EQ. LSTSEG )THEN
NEWSEG = SEGTAB(NXTSEG,2)
ENDIF
C
C--- Store next element of chain
C
ILOOP = ILOOP + 1
LOOP(ILOOP) = NXTSEG
LSTSEG = NXTSEG
NXTSEG = NEWSEG
C
C--- Test if loop complete
C
IF (NXTSEG .EQ. ISEG)GO TO 702
C
C--- Points to next element in chain
C
GO TO 703
C
C--- End of branch one
C
701 CONTINUE
C
C--- This cannot be a loop so kill off node
C
NSGTAB(ISEG) = -1
C
C--- Now remove any reference to this node in the remaining nodes
C
DO 801 LSEG = 1, NSEG
IF ( NSGTAB(LSEG) .EQ. -1) GO TO 801
DO 802 ID = 1,NSGTAB(LSEG)
IF( SEGTAB(LSEG,ID) .NE. ISEG ) GO TO 802
SEGTAB(LSEG,ID) = SEGTAB(LSEG,NSGTAB(LSEG))
NSGTAB(LSEG) = NSGTAB(LSEG) - 1
802 CONTINUE
801 CONTINUE
C
C--- Now continue on node loop
C
GO TO 700
C
C--- Loop complete
C
702 CONTINUE
C
C--- Perform fits and eliminate all but one node in loop
C
DO 860 KLOOP = 1,ILOOP
KSEG = LOOP(KLOOP)
DO 861 IWIRE = 1,4
IF( IDIGST(IWIRE,KSEG) .EQ. -1)THEN
Y(IWIRE) = 0.0
W(IWIRE) = 0.0
ELSE
Y(IWIRE) = DRSTO( IDIGST(IWIRE,KSEG) , IWIRE )
W(IWIRE) = 1.0/RESSTO( IDIGST(IWIRE,KSEG) , IWIRE)**2
ENDIF
861 CONTINUE
CALL FPCFIT(IP,Y,W,SLOPE,ZERO,COVSLZ,CHISQ,PBCHI)
PROBL(KLOOP) = PBCHI
860 CONTINUE
C
C--- Keep only highest probability iff prob gt 10**-5
C
LPSAVE = LVMAX(PROBL,ILOOP)
VALPRB = PROBL(LPSAVE)
IF( VALPRB .GT. 0.00001)THEN
KSSAVE = LOOP(LPSAVE)
ELSE
KSSAVE = -1
ENDIF
C
C--- now loop over segments and remove all but the saved
C
DO 870 KLOOP = 1,ILOOP
KSEG = LOOP(KLOOP)
IF( KSEG .EQ. KSSAVE ) GO TO 870
NSGTAB(KSEG) = -1
C
C--- Now remove any reference to this node in the remaining nodes
C
DO 871 LSEG = 1, NSEG
IF ( NSGTAB(LSEG) .EQ. -1) GO TO 871
DO 872 ID = 1,NSGTAB(LSEG)
IF( SEGTAB(LSEG,ID) .NE. KSEG ) GO TO 872
SEGTAB(LSEG,ID) = SEGTAB(LSEG,NSGTAB(LSEG))
NSGTAB(LSEG) = NSGTAB(LSEG) - 1
872 CONTINUE
871 CONTINUE
870 CONTINUE
875 CONTINUE
C
700 CONTINUE
C
C--- now remove any pairs by fitting
C
DO 900 ISEG = 1,NSEG
IF ( NSGTAB(ISEG) .NE. 1)GO TO 900
C
C--- Found a pair so find partner
C
ISEGP = SEGTAB(ISEG,1)
C
C--- fit the first possiblity
C
DO 901 IWIRE = 1,4
IF( IDIGST(IWIRE,ISEG) .EQ. -1)THEN
Y(IWIRE) = 0.0
W(IWIRE) = 0.0
ELSE
Y(IWIRE) = DRSTO( IDIGST(IWIRE,ISEG) , IWIRE )
W(IWIRE) = 1.0/RESSTO( IDIGST(IWIRE,ISEG) , IWIRE)**2
ENDIF
901 CONTINUE
CALL FPCFIT(IP,Y,W,SLOPE,ZERO,COVSLZ,CHISQ1,PBCHI1)
C
C--- Fit the second possiblity
C
DO 902 IWIRE = 1,4
IF( IDIGST(IWIRE,ISEGP) .EQ. -1)THEN
Y(IWIRE) = 0.0
W(IWIRE) = 0.0
ELSE
Y(IWIRE) = DRSTO( IDIGST(IWIRE,ISEGP) , IWIRE )
W(IWIRE) = 1.0/RESSTO( IDIGST(IWIRE,ISEGP) , IWIRE)**2
ENDIF
902 CONTINUE
CALL FPCFIT(IP,Y,W,SLOPE,ZERO,COVSLZ,CHISQ2,PBCHI2)
C
C--- Now remove the smaller probability segment
C
IF ( PBCHI1 .GT. PBCHI2) THEN
KSEG = ISEGP
ELSE
KSEG = ISEG
ENDIF
NSGTAB(KSEG) = -1
C
C--- Now remove any reference to this node in the remaining nodes
C
DO 911 LSEG = 1, NSEG
IF ( NSGTAB(LSEG) .EQ. -1) GO TO 911
DO 912 ID = 1,NSGTAB(LSEG)
IF( SEGTAB(LSEG,ID) .NE. KSEG ) GO TO 912
SEGTAB(LSEG,ID) = SEGTAB(LSEG,NSGTAB(LSEG))
NSGTAB(LSEG) = NSGTAB(LSEG) - 1
912 CONTINUE
911 CONTINUE
900 CONTINUE
C
C--- Now draw remaining segments
C
DO 600 ISEG = 1,NSEG
IF( NSGTAB(ISEG) .EQ. -1 ) GO TO 600
C
C--- fit remaining segments
C
DO 670 IWIRE = 1,4
IF( IDIGST(IWIRE,ISEG) .EQ. -1)THEN
Y(IWIRE) = 0.0
W(IWIRE) = 0.0
ELSE
Y(IWIRE) = DRSTO( IDIGST(IWIRE,ISEG) , IWIRE )
W(IWIRE) = 1.0/RESSTO( IDIGST(IWIRE,ISEG) , IWIRE)**2
ENDIF
670 CONTINUE
CALL FPCFIT(IP,Y,W,SLOPE,ZERO,COVSLZ,CHISQ,PBCHI)
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C--- find the absolute coordinate normal current wire orientation
C
COWIRE = 0.0
C
C--- find which orientation plane (1-9) is current
C
IPLANE = IP/4 + 1
C
C--- find fitted 'y' position at begining and end of 4 wire set
C--- fdrsta and fdrend respectively
C
FDRSTA = COWIRE + ZERO
FDREND = FDRSTA + SLOPE*(
1 ZPLAN( (IPLANE)*4 ) - ZPLAN( (IPLANE-1)*4 + 1) )
C
C--- Transform from orientation drift coordinates to global coordinates
C--- assuming that drift x coordinate is zero
C
XREAL1= -STP(IPLANE) * FDRSTA
YREAL1= CTP(IPLANE) * FDRSTA
XREAL2= -STP(IPLANE) * FDREND
YREAL2= CTP(IPLANE) * FDREND
C
C--- Fill track cluster banks and banks counter
C
C
C--- Increment cluster counter per plane
C
IF( NTC(IPLANE) .GE. MAXCLU)THEN
CALL ERRLOG(208,'W:FPDG32: NTC(IPLANE) > MAXCLU')
ELSE
NTC(IPLANE) = NTC(IPLANE) + 1
ENDIF
TC(1,IPLANE,NTC(IPLANE)) = XREAL2 - XREAL1
TC(2,IPLANE,NTC(IPLANE)) = YREAL2 - YREAL1
TC(3,IPLANE,NTC(IPLANE)) =
1 ZPLAN((IPLANE)*4 ) - ZPLAN((IPLANE-1)*4 + 1)
C
C--- store toc
C
TOC(1,IPLANE,NTC(IPLANE))=XREAL1
TOC(2,IPLANE,NTC(IPLANE))=YREAL1
TOC(3,IPLANE,NTC(IPLANE))=ZPLAN( (IPLANE-1)*4 + 1)
C
C--- store the digitisings associated with plane/track for final
C fit
C
DO 695 IWW = 1,4
IF(IDIGST(IWW,ISEG) .NE. -1)THEN
IDGISM(IWW,IPLANE,NTC(IPLANE)) =
1 IDGIST( IDIGST(IWW,ISEG) , IWW )
ELSE
IDGISM(IWW,IPLANE,NTC(IPLANE)) = 0
ENDIF
TCYUV(IWW,IPLANE,NTC(IPLANE)) = COWIRE + Y(IWW)
TCYUVW(IWW,IPLANE,NTC(IPLANE)) = W(IWW)
695 CONTINUE
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
600 CONTINUE
C
C--- Now plot unused digitisizings
C
C
C--- Set drmask to 1.0 for used digits
C
DO 750 ISEG = 1,NSEG
IF(NSGTAB(ISEG) .EQ. -1) GO TO 750
DO 751 ID = 1,4
ID1 = IDIGST(ID,ISEG)
IF (ID1 .EQ. -1) GO TO 751
IF ( MOD(ID1,2) .EQ. 0 )THEN
ID2 = ID1 - 1
ELSE
ID2 = ID1 + 1
ENDIF
DRMASK(ID1,ID) = .TRUE.
DRMASK(ID2,ID) = .TRUE.
751 CONTINUE
750 CONTINUE
IDUNUS = 0
DO 760 IWIRE = 1,4
DO 761 ID = 1,NDRSTO(IWIRE)
IF( DRMASK(ID,IWIRE) ) GO TO 761
DRIFT = DRSTO(ID,IWIRE)
IDUNUS = IDUNUS + 1
761 CONTINUE
760 CONTINUE
C
C--- If first pass now loop back with 2*tolerance and 3*maxang
C
IF (IDUNUS .NE. 0 .AND. IFIRST .EQ. 1)THEN
IFIRST = 0
TOLER = RESOL*8.0
SLYMAX = 1.0
SLYMIN = 0.0
SLMAX = 5.0
GO TO 1100
ENDIF
END
*