subroutine prmain(ndig)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Main routine to call pattern recognition
C-
C- Controls:
C-
C-
C-
C-
C- Created 20-DEC-1994 Alexander Andryakov
C-
C----------------------------------------------------------------------
INTEGER NDIG
c
c Prepare data for the p.r. package
c
call ptrnprep(ndig)
c
c 1-st stereo view ( + )
c
call dtptrn(1)
call dtsort(1)
c
c 2-nd stereo view ( - )
c
call dtptrn(2)
call dtsort(2)
call dtdist
c
c create YBOS banks for each t.c. (temporary put here, will go inside p.r)
c
CALL DPRARG
c
c merge results of 1-view p.r. to get space cand.
c
call v_merge
c
c 0-approximation fit to get z0 and cot(theta)
c
call z0_ctgth
c call ch_restor
end
SUBROUTINE PTRNPREP(n_hits)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Created 21-APR-1993 Alexander Andryakov
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDTHIT.INC'
$$INCLUDE 'K$ITRK:CTRACK.INC'
$$INCLUDE 'K$ITRK:MCKINE.INC'
INTEGER N_HITS
INTEGER IHIT_N,IHITS_CH,NLAYER_C
INTEGER I_ST_SIGN,I_LAY,NLAY_C
INTEGER ibuf,ibufidt,IBUFHIT_N
REAL bufp,bufd,bufe
DIMENSION ibuf(numhit),bufp(numhit),bufd(numhit),bufe(numhit),
& ibufidt(numhit),ibufhit_n(numhit)
integer old_hit_num
REAL btidra,btiddp,btidst,btidsc
common/chorod/
& btidra(nlayero),btiddp(nlayero),btidst(nlayero),btidsc(nlayero)
C
C---- Clean CWORK working area again ----------------
C
call clean_cwork
C
C----------------------------------------------------
c
c Prepare key-list for sorting
c
ihit_n=0
do ihits_ch=1,n_hits
if(LDT(ihits_ch).lt.1) then
call INERR(210,1)
ibuf(ihits_ch)=3*10**5
go to 23
endif
c LDT(ihits_ch)=NLAYER-LDT(ihits_ch)+1
nlayer_c=LDT(ihits_ch)
if(TIDST(nlayer_c).lt.0) then
i_st_sign=1
else
i_st_sign=2
endif
ibuf(ihits_ch)=i_st_sign*10**5+LDT(ihits_ch)*1000+n_wir(ihits_ch)
ihit_n=ihit_n+1
ibufidt(ihits_ch)=IDT(ihits_ch)
bufd(ihits_ch)=DTD(ihits_ch)
bufe(ihits_ch)=DTE(ihits_ch)
bufp(ihits_ch)=DTP(ihits_ch)
ibufhit_n(ihits_ch)=MDT(ihits_ch)
23 continue
enddo
call sortzv(ibuf,list_hit,n_hits,-1,0,0)
c
c Fill or refill array with new order
c
do ihits_ch=1,ihit_n
old_hit_num=list_hit(ihits_ch)
idtr(ihits_ch)=ibufidt(old_hit_num)
c IDT(old_hit_num)=-999
IDT(ihits_ch)=-999
dtdsr(ihits_ch)=bufd(old_hit_num)
LDT(ihits_ch)=mod(ibuf(old_hit_num),100000)/1000
DTD(ihits_ch)=abs(dtdsr(ihits_ch))
DTE(ihits_ch)=bufe(old_hit_num)
DTP(ihits_ch)=bufp(old_hit_num)
dtds(ihits_ch)=-5000.
MDT(ihits_ch)=ibufhit_n(old_hit_num)
enddo
c
c Fill NLDTB-array
c
do i_lay=1,nlayer
NLDTB(1,i_lay)=1
NLDTB(2,i_lay)=0
enddo
nlay_c=1
do ihits_ch=1,ihit_n
if(LDT(ihits_ch).ne.nlay_c) then
NLDTB(1,LDT(ihits_ch))=ihits_ch
NLDTB(2,LDT(ihits_ch))=ihits_ch
nlay_c=LDT(ihits_ch)
else
NLDTB(2,LDT(ihits_ch))=ihits_ch
endif
enddo
c
c Check if NLDTB has been filled correctly
c
* write(6,'(''1 layer first wire last wire'')')
* do ipr=1,nlayer
* write(6,'(3x,i2,7x,i3,8x,i3))') ipr,NLDTB(1,ipr),NLDTB(2,ipr)
* enddo
c if(n_hits.gt.ihit_n) then
CALL VZERO(IDT(ihit_n+1),numhit-ihit_n)
CALL VZERO(LDT(ihit_n+1),numhit-ihit_n)
CALL VZERO(DTP(ihit_n+1),numhit-ihit_n)
CALL VZERO(DTD(ihit_n+1),numhit-ihit_n)
n_hits=ihit_n
c endif
c
c Erase all arrays which have to be filled in DTPTRN
c
LWMDTC = 0 ! current position in MB
LWMDTL = 0 ! last usable position in MB
LWKDTC = 0 ! last bank in KB
LWIDTC = 0 ! last bank in IB
LWIDT1 = 0 ! pointer to 1st track in IB
LWIDTL = 0 ! pointer to last track in IB
LWJDTC = 0 ! last bank in JB
call vzero(KB(1,1),4*LWKDT)
call vzero(IB(1,1),8*LWIDT)
call vzero(JB(1,1),10*LWJDT)
END
subroutine clean_cwork
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : clean CWORK area
C-
C- Created 11-JUL-1994 Alexander Andryakov
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CDTWORK.INC'
call vzero(MB,18000)
end
SUBROUTINE CH_RESTOR
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Created 26-APR-1993 Alexander Andryakov
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CTIDRI.INC'
INTEGER INOR
REAL BTIDRA,BTIDDP,BTIDST,BTIDSC
common/chorod/
& btidra(nlayero),btiddp(nlayero),btidst(nlayero),btidsc(nlayero)
c
c Restore original order of the chamber constant arrays
c
do inor=1,NLAYER
tidra(inor)=btidra(inor)
tiddp(inor)=btiddp(inor)
tidst(inor)=btidst(inor)
tidsc(inor)=btidsc(inor)
enddo
END
C 27/11/89 MEMBER NAME DTPTRN (ES) FVS
C 13/03/81 MEMBER NAME DTPTRN (ES) FORTRAN
C
SUBROUTINE DTPTRN(i_view)
$$IMPLICIT
C
C PATTERN RECOGNITION
C
C --- HARTWIG ALBRECHT --- 80-04 --- VERSION 82-10-11
C
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CDTHIT.INC'
$$INCLUDE 'K$ITRK:CTRACK.INC'
$$INCLUDE 'K$ITRK:CDTWORK.INC'
$$INCLUDE 'K$ITRK:CDTCUT.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
INTEGER I_VIEW
INTEGER MOPT,MN1,MN0,LABELT,JJ,JM,NB3,NY3,NT3,NT2,NUSED
INTEGER NT1,LABEL8,NTX,NTX2,NTD,NTA,LABELX
REAL CH
REAL PB,RT2,RT3,B2,B3,GB1,GB2,DT1,DT2,DT3,V1,V2
REAL CUR1,CUR2
REAL CUTV0,CUTV1,CUTV2,CUTC
INTEGER LABELM,MAZ,MR,L,NY,NZ,M,LXTR,LXT2,LXT3
C ------------------ my changes ---- A.Andryakov ----------------
C Continuous fraction presentation of sqrt is used in original
C code. I've changed it to standard call to intrinsic fortran
C function sqrt. But of course I keep the old presentation as well.
C Logical variable cont_fract control the usage of these two poss.
C cont_fract = .true. - continuous fraction presentation
C = .false. - standard fortran
C
logical cont_fract
data cont_fract /.false./
C
C ------------------ my changes ---- A.Andryakov ----------------
C
REAL G1_cut,G2_cut
REAL PD2,PD3,D1,F1,G2,F2,R4,R2,RD,B1,SDP
REAL SDPT,CDP,CUTV,S1,G1,CHIM,C1,S1T
REAL PD,R,RSQ,CUTP0,CUTP3,CUTPS,CUTP,CUTPHI,PIADD
INTEGER N1,NX,ND,K2,N2,K3,N3,N4
INTEGER NB,J,J2,LPREV,MAA,MAC,MZC,MC,LABELA,NY2,NB2
INTEGER MMAIN,MTRACK,MDISCO,MPREV,MCURR,INEXT,IDTFL,LABELY
LOGICAL FLSOL,FLSOL1,FLAGT
DIMENSION GB1(2),GB2(2),MAZ(4,5)
DIMENSION JJ(64),PB(64),CUTP(5)
DATA CUTV0,CUTV1,CUTV2 /0.0,-1.0,-4.0/,
+ CUTC /-0.02/
logical entry_flag
C ------------------------------------------------------------------
C ------------------------------------------------- INITIALIZATION
C ------------------------------------------------------------------
entry_flag=.true.
LABELM = 1
MAZ(1,1) = 0
MAZ(1,2) = 0
MAZ(1,3) = 0
MAZ(1,4) = 0
MAZ(1,5) = 0
MAZ(2,1) = 0
MAZ(2,2) = 0
MAZ(2,3) = 0
MAZ(2,4) = 0
MAZ(2,5) = 0
MAZ(3,1) = 0
MAZ(3,2) = 0
MAZ(3,3) = 0
MAZ(3,4) = 0
MAZ(3,5) = 0
LWMDTL = LWMDT
MR = 0
L = i_view-2
NY = 0
NZ = 0
M = 1
c CUTDNB = DTCDNB ! 0.4(7) * cell_size
c CUTDNB = 0.477777 * cell_size
c CUTD1 = cell_size - 0.01
c CUTD1 = cell_size - 0.01
LXTR = LDTXTR
LXT2 = LDTXT2
LXT3 = LDTXT3
CUR1 = DTXUR1
CUR2 = DTXUR2
C ------------------------------------------------------------------
C -------------------------------------------- LOOP OVER ALL WIRES
C ------------------------------------------------------------------
100 NB = NY + 1
IF (NB.LE.NZ) GO TO 130
C ------------------------------------------------------------------
C ------------------------------------------------------ NEW LAYER
C ------------------------------------------------------------------
110 IF (L.EQ.NLAYER) RETURN ! stereo chamber
L = L + 2
IF (L.gt.NLAYER) RETURN
MAZ(2,5) = MAZ(2,4)
MAZ(2,4) = MAZ(2,3)
MAZ(2,3) = MAZ(2,2)
MAZ(2,2) = MAZ(2,1)
MAZ(3,5) = MAZ(3,4)
MAZ(3,4) = MAZ(3,3)
MAZ(3,3) = MAZ(3,2)
MAZ(3,2) = MAZ(3,1)
MAZ(4,5) = MAZ(4,4)
MAZ(4,4) = MAZ(4,3)
MAZ(4,3) = MAZ(4,2)
MAZ(4,2) = MAZ(4,1)
NZ = NLDTB(2,L)
IF (NZ.GE.NLDTB(1,L)) GO TO 120
MAZ(2,1) = 0
GO TO 110
C-----
120 ny=nldtb(1,l)-1
c if(entry_flag) then
c entry_flag=.false.
nb=ny+1
c endif
c write(6,'('' l,nb,ny: '',3i4)') l,nb,ny
MAZ(2,1) = M
MAZ(3,1) = M
MAZ(2,2) = MAZ(3,2)
MAZ(2,3) = MAZ(3,3)
MAZ(2,4) = MAZ(3,4)
MAZ(2,5) = MAZ(3,5)
PD = TIDDP(L)
R = TIDRA(L)
RSQ = R * R
C ----------------------------------------------------------- CUTS
CUTP0 = PD * DTXPN0
CUTP3 = PD * DTXPN3
CUTPS = PD * 1.1
CUTP(1) = PD * DTXPN1
CUTP(2) = PD * DTXPN2
CUTP(3) = CUTP(2)
CUTP(4) = CUTP(3)
CUTP(5) = CUTP(4)
C ------------------------------------------------------------------
C ------------------------------------------- GET CLUSTER OF WIRES
C ------------------------------------------------------------------
130 NY = NY + 1
IF (DTP(NY+1)-DTP(NY).GE.CUTP0
+ . OR . NY.EQ.NZ . OR.NY.GE.NB+15) GO TO 150
IF (DTP(NY+1)-DTP(NY).GE.CUTPS) GO TO 130
IF (DTD(NY+1)+DTD(NY).GE.
+ max(cell_size(ldt(ny)),cell_size(ldt(ny+1))) - 0.01) GO TO 135
IF (DTD(NY).LT.0.477777 * cell_size(ldt(ny)) . OR .
+ DTD(NY+1).LT.0.477777 * cell_size(ldt(ny+1))) GO TO 130
135 IF (DTD(NY).GT.DTD(NY+1)) GO TO 140
IDT(NY+1) = 0
GO TO 130
C-----
140 IDT(NY) = 0
GO TO 130
C-----
150 J = 0
J2 = 0
LPREV = 1
C ------------------------------------------- SET POINTERS IN MBUF
c write(6,'('' l,nb,ny: '',3i4)') l,nb,ny
MAZ(4,1) = M
MB(1,M) = NB
MB(2,M) = NY
MB(1,M+2) = 0
MB(2,M+2) = 0
LWMDTC = M + 3
IF (LWMDTC.LT.LWMDTL) GO TO 210
GO TO 910
C
C-----------------------------------------------------------------------
C------------------------------------------- LOOP OVER PREVIOUS LAYERS
C-----------------------------------------------------------------------
200 MAZ(2,LPREV+1) = MAA
IF (LPREV.EQ.LXTR) GO TO 600
LPREV = LPREV + 1
210 CUTPHI = CUTP(LPREV)
MAA = MAZ(2,LPREV+1)
C EMPTY LAYER ?
IF (MAA.EQ.0) GO TO 200
MAC = MAZ(3,LPREV+1)
MZC = MAZ(4,LPREV+1)
IF (MAC.NE.MAA) GO TO 260
C ------------ NEIGHBOURS IN THE PREVIOUS LAYERS AT PHI = 0 / 2*PI
MC = MZC
PIADD = -PI2
ASSIGN 580 TO LABELA
220 NY2 = MB(2,MC)
IF (DTP(NB) - DTP(NY2) + PI2 . GT . CUTPHI) GO TO 260
NB2 = MB(1,MC)
GO TO 300
C-----
230 IF (MAC.EQ.MZC) GO TO 200
MC = MAC
240 NB2 = MB(1,MC)
IF (DTP(NB2) - DTP(NY) + PI2 . GT . CUTPHI) GO TO 200
NY2 = MB(2,MC)
PIADD = PI2
ASSIGN 590 TO LABELA
GO TO 300
C
C--------------------------------------------------------- NO ACTION :
C DELTAPHI = 0
260 PIADD = 0.
ASSIGN 570 TO LABELA
C 1ST NEIGHBOUR :
MC = MAA
270 NB2 = MB(1,MC)
NY2 = MB(2,MC)
IF (DTP(NB)-DTP(NY2).LE.CUTPHI) GO TO 280
IF (MC.EQ.MZC) GO TO 230
MC = MB(1,MC+1)
MAA = MC
GO TO 270
C-----
280 IF (DTP(NB2)-DTP(NY).GT.CUTPHI) GO TO 200
C ------------------------------------------------------------------
C --------------------------------------------------------- TRACKS
C ------------------------------------------------------------------
300 MMAIN = MC + 2
MTRACK = MB(1,MMAIN)
IF (MTRACK.EQ.0) GO TO 550
310 MDISCO = MMAIN
MPREV = MTRACK
MCURR = MTRACK
INEXT = -1
IDTFL = -2
FLSOL = .FALSE.
FLAGT = .TRUE.
C ------------------------------------------------------------------
C --------------------------------- LOOP OVER WIRES OF THE CLUSTER
C ------------------------------------------------------------------
320 ASSIGN 420 TO LABELY
IF (BM(2,MPREV+2).LT.0.) GO TO 321
N1 = NY
NX = NB
ND = -1
GO TO 330
C-----
321 N1 = NB
NX = NY
ND = 1
330 IF (IDT(N1).GE.0) GO TO LABELY, (420,430,440)
K2 = MB(2,MPREV)
N2 = KB(2,K2)
340 PD = DTP(N2) - DTP(N1) + PIADD
IF (ABS(PD).GT.CUTPHI) GO TO 410
PD2 = BM(2,MPREV+1)
PD3 = PD + PD2
D1 = DTD(N1)
K3 = KB(1,K2)
F1 = D1 - BK(3,K2)
G2 = BK(4,K2) + BK(4,K3)
N3 = KB(2,K3)
K3 = KB(1,K3)
F2 = BK(3,K2) - BK(3,K3)
N4 = KB(2,K3)
R4 = TIDRA(LDT(N4))
R2 = TIDRA(LDT(N2))
RD = R2 - R
IF (RD.EQ.0.) GO TO 350
C -------------------------------- W1 AND W2 NOT IN THE SAME LAYER
B1 = RD * RD + R * R2 * (PD * PD)
SDP = B1 * BM(1,MPREV+2)
sdpt=sqrt(SDP)
CDP = ((R - R4)**2 + R * R4 * (PD3 * PD3) -
+ B1 - BM(1,MPREV+2)) * 0.5 / SDP
IF (CDP.LE.0.) GO TO 410
SDP = (R * R2 * PD * (1. - (PD * PD ) * .177778) +
+ R2 * R4 * PD2 * (1. - (PD2 * PD2) * .177778) -
+ R * R4 * PD3 * (1. - (PD3 * PD3) * .177778)) / SDP
CUTV = CUTV1
IF (RD.GT.4.) CUTV = CUTV2
S1 = B1 - F1 * F1
G1 = 0.5 * (RD + S1 / RD)
GO TO 360
C
C----------------------------------------- W1 AND W2 IN THE SAME LAYER
350 CONTINUE
B1 = (R * PD) * (R * PD)
SDP = B1 * BM(1,MPREV+2)
sdpt=sqrt(SDP)
if (abs(sdp).lt.1.e-10) then
write(6,1000)
& sdp, b1, bm(1,mprev+2), pd,
& n1, ldt(n1), dtp(n1),
& n2, ldt(n2), dtp(n2)
1000 FORMAT(' sdp, b1, bm(1,mprev+2), pd:',4g13.5/
& ' n1, ldt(n1), dtp(n1): ',2i5,g13.5/
& ' n2, ldt(n2), dtp(n2): ',2i5,g13.5)
go to 410
endif
CDP = ((R - R4)**2 + R * R4 * (PD+PD2)**2 -
+ B1 - BM(1,MPREV+2)) * 0.5 / SDP
IF (CDP.LE.CUTC) GO TO 410
SDP = R * (R4 * (PD2 * (1. - PD2 * PD2 * .177778) -
+ PD3 * (1. - PD3 * PD3 * .177778)) +
+ R * PD * (1. - PD * PD * .177778)) / SDP
CUTV = CUTV0
S1 = B1 - F1 * F1
G1 = 0.5 * (1.3 + S1)
360 CHIM = 1.E38
IF (LWKDTC.GT.LWKDT) GO TO 960
IF (LWMDTC.GE.LWMDTL) GO TO 911
C ------------------------------------------------------------------
C ----------------------------------------- TRY RIGHT AND LEFT AMB
C ------------------------------------------------------------------
370 F1 = D1 - BK(3,K2)
BM(2,LWMDTC+2) = (PD * R - F1)
C1 = BM(2,LWMDTC+2) * BM(2,MPREV+2)
IF (C1.LT.CUTV) GO TO 400
if(cont_fract) then
G1 = 0.5 * (G1 + S1 / G1)
G1 = 0.5 * (G1 + S1 / G1)
G1 = 0.5 * (G1 + S1 / G1)
else
if(s1.lt.0.) then
call INERR(211,1)
go to 400
else
G1 = sqrt(s1) ! I switched off sqrt calc. via cont.fract.
endif
endif
S1 = G1 + G2
s1t=s1
C ------------------------------------------------------ CURVATURE
C1 = ((G1*F2 - F1*G2) * CDP -
+ (F1*F2 + G1*G2) * SDP) / S1
IF (ABS(C1).GT.CUR1) GO TO 400
C ------------------------------------------- 1 / ERROR**2 ON CURV
S1 = (S1*S1) * (G1*G1) * (G2*G2) /
+ (DTE(N1)*(G2*G2) + DTE(N4)*(G1*G1) + DTE(N2)*(S1*S1))
C ------------------------------------------------ SUM 1/ERR(N)**2
BM(2,LWMDTC+3) = S1 + BM(2,MPREV+3)
C ------------------------------------------ CHISQ OF THE NEW WIRE
CH = S1 * BM(2,MPREV+3) *
+ (C1 - BM(1,MPREV+3))**2 / BM(2,LWMDTC+3)
IF (CH.GT.chi2_pr_cut) GO TO 400
IF (CH.GT.CHIM) GO TO 380
c call hf1(1,c1*s1t-sdp*sdpt,1.)
c call hf1(2,c1*s1t,1.)
c call hf1(3,sdp*sdpt,1.)
CHIM = CH
MOPT = LWMDTC
C --------------------------------------------------------- <CURV>
380 BM(1,LWMDTC+3) = (C1 * S1 +
+ BM(1,MPREV+3) * BM(2,MPREV+3)) / BM(2,LWMDTC+3)
C ------------------------------------------------ FILL K - BUFFER
LWKDTC = LWKDTC + 1
KB(1,LWKDTC) = K2
KB(2,LWKDTC) = N1
BK(3,LWKDTC) = D1
BK(4,LWKDTC) = G1
dtds(n1)=d1
C ------------------------------------------------ FILL TEMP BUFFER
MB(2,LWMDTC) = LWKDTC
MB(1,LWMDTC+1) = MPREV
BM(2,LWMDTC+1) = BM(2,MPREV+5) + PD
BM(1,LWMDTC+2) = (R - TIDRA(LDT(N3)))**2 +
+ R * TIDRA(LDT(N3)) * BM(2,LWMDTC+1)**2
MB(1,LWMDTC+4) = 0
MB(2,LWMDTC+4) = LDT(N2)
BM(1,LWMDTC+5) = BM(1,MPREV+5) + CH - 1.
BM(2,LWMDTC+5) = PD
C ----------------------------------------------- POINTERS IN MB :
IF (.NOT.FLSOL) GO TO 384
MN1 = MB(1,M+2)
IF (BM(1,MN1+5).GE.BM(1,LWMDTC+5)) GO TO 382
MN0 = MN1
381 MN1 = MB(1,MN0+1)
IF (BM(1,MN1+5).GE.BM(1,LWMDTC+5)) GO TO 383
MN0 = MN1
GO TO 381
C-----
382 MB(1,LWMDTC) = MB(1,MN1)
MB(1,MN1) = 0
MB(1,MN1+4) = 0
MB(1,M+2) = LWMDTC
MB(1,LWMDTC+1) = MN1
GO TO 390
C-----
383 MB(1,LWMDTC) = 0
MB(1,MN0+1) = LWMDTC
MB(1,LWMDTC+1) = MN1
GO TO 390
C-----
384 FLSOL = .TRUE.
IF (MDISCO.EQ.0) GO TO 386
MB(1,LWMDTC+1) = MB(1,MDISCO)
MB(1,MDISCO) = MB(1,MB(1,MDISCO))
MDISCO = 0
GO TO 387
386 MB(1,LWMDTC+1) = MPREV
387 MB(1,LWMDTC) = MB(1,M+2)
MB(1,MPREV) = 0
MB(1,M+2) = LWMDTC
390 ASSIGN 430 TO LABELY
INEXT = MAX0 (INEXT,MB(1,MPREV+4))
MB(1,MPREV+4) = 0
LWMDTC = LWMDTC + 6
400 IF (D1.LE.0.) GO TO 410
D1 = -D1
F1 = D1 - BK(3,K2)
S1 = B1 - F1 * F1
G1 = 0.5 * (G1 + S1 / G1)
GO TO 370
C
C -----------------------------------------------------------------
C ----------------------------------------- END OF LOOP OVER AMB'S
C----------------------------------------------------------------------
410 CONTINUE
C
GO TO LABELY, (420,430,440)
C --------------------------------- WIRE DOES NOT FIT TO THE TRACK
420 IF (N1.EQ.NX) GO TO 490
N1 = N1 + ND
GO TO 330
C
C----------------------------------------------------------- WIRE FITS
430 ASSIGN 440 TO LABELY
IF (.NOT.FLAGT . OR .
+ BM(1,MPREV+5).GT.-1.5) GO TO 431
IDT(N1) = -1
GO TO 432
C-----
431 IDT(N1) = MAX0 (IDT(N1),IDTFL)
IDT(N2) = MAX0 (IDT(N2),IDTFL)
IDT(N3) = MAX0 (IDT(N3),IDTFL)
IDT(N4) = MAX0 (IDT(N4),IDTFL)
432 IF (N1.EQ.NX) GO TO 480
N1 = N1 + ND
MTRACK = MPREV
MPREV = MOPT
GO TO 330
C
C------------------------------------- WIRE DOES NOT FIT TO THE TRACK,
C BUT ANOTHER ONE DID.
440 IF (.NOT.FLAGT) GO TO 460
IF (BM(1,MCURR+5).LT.0.) GO TO 450
IF (MPREV.EQ.MCURR) GO TO 491
MPREV = MCURR
GO TO 330
C-----
450 IF (MPREV.EQ.MCURR) GO TO 530
MPREV = MCURR
GO TO 330
C-----
460 IF (MPREV.EQ.MTRACK) GO TO 470
MPREV = MTRACK
GO TO 330
C-----
470 IF (N1.EQ.NX) GO TO 510
MPREV = MOPT
N1 = N1 + ND
GO TO 330
C
C -----------------------------------------------------------------
C ----------------------------------------- END OF LOOP OVER WIRES
C----------------------------------------------------------------------
480 IF (BM(1,MCURR+5).GE.0.) GO TO 490
IF (FLAGT) GO TO 530
GO TO 510
C
C -----------------------------------------------------------------
C ------------------------------------------------------- GO BACKW
C----------------------------------------------------------------------
490 IF (.NOT.FLAGT) GO TO 510
491 MPREV = MB(1,MCURR+1)
IF (MPREV.EQ.0) GO TO 530
IF (L-MB(2,MCURR+4).GT.LXTR) GO TO 530
IF (BM(1,MCURR+5).LT.0.) MDISCO = MCURR + 1
MCURR = MPREV
GO TO 320
C
C -----------------------------------------------------------------
C ------------------------------------------------ EXTEND TRIPLETS
C----------------------------------------------------------------------
500 MCURR = M + 1
INEXT = -1
C NEXT TRIPLET
510 MPREV = MB(1,MCURR+1)
IF (MPREV.EQ.0) GO TO 520
MCURR = MPREV
ASSIGN 420 TO LABELY
K2 = MB(2,MPREV)
N2 = KB(2,K2)
C ------------------------------------------------------------------
C --------------------------------- LOOP OVER WIRES OF THE CLUSTER
C ------------------------------------------------------------------
IF (BM(2,MPREV+2).LT.0.) GO TO 511
IF (N2.EQ.NB) GO TO 510
NX = NB
ND = -1
N1 = N2 + ND
IF (IDT(N1).LT.0) GO TO 340
GO TO 420
C-----
511 IF (N2.EQ.NY) GO TO 510
NX = NY
ND = 1
N1 = N2 + ND
IF (IDT(N1).LT.0) GO TO 340
GO TO 420
C
C----------------------------------------------------- FILL I - BUFFER
520 IF (INEXT.LT.0) GO TO LABELT, (691,758,818,867)
IF (LWIDTC.EQ.LWIDT) GO TO 970
LWIDTC = LWIDTC + 1
MOPT = MB(1,M+2)
MB(1,MOPT+4) = LWIDTC
IB(2,LWIDTC) = 0
IB(3,LWIDTC) = MB(2,MOPT)
IB(7,LWIDTC) = MB(1,MOPT+3)
IB(8,LWIDTC) = MB(1,MOPT+5)
GO TO LABELT, (691,758,818,867)
C-----
530 IF (INEXT) 540,531,532
531 IF (LWIDTC.EQ.LWIDT) GO TO 970
LWIDTC = LWIDTC + 1
INEXT = LWIDTC
532 MOPT = MB(1,M+2)
MB(1,MOPT+4) = INEXT
IB(2,INEXT) = 0
IB(3,INEXT) = MB(2,MOPT)
IB(7,INEXT) = MB(1,MOPT+3)
IB(8,INEXT) = MB(1,MOPT+5)
C ----------------------------------------------------- NEXT TRACK
540 IF (MB(1,MMAIN).EQ.MTRACK) MMAIN = MTRACK
MTRACK = MB(1,MMAIN)
IF (MTRACK.NE.0) GO TO 310
550 IF (LPREV.GT.LXT3) GO TO 561
J = J + 1
if(j.ge.64) write (6,'('' @@@ DTPTRN @@@ j ='',i4)') j
IF (LPREV.LE.LXT2) J2 = J
JJ(J) = MC
PB(J) = PIADD
C -----------------------------------------------------------------
C ---------------------------------------------------------- PAIRS
C -----------------------------------------------------------------
IF (LPREV.GT.LXT2) GO TO 561
IF (LWMDTC.GE.LWMDTL) GO TO 912
560 MB(1,LWMDTC) = MB(2,M+2)
MB(2,M+2) = LWMDTC
BM(2,LWMDTC) = PIADD
MB(1,LWMDTC+1) = MB(1,MC)
MB(2,LWMDTC+1) = MB(2,MC)
LWMDTC = LWMDTC + 2
561 GO TO LABELA, (580,570,590)
C -----------------------------------------------------------------
C ------------------------------------------------- NEXT CLUSTER :
C -----------------------------------------------------------------
570 IF (MC.EQ.MZC) GO TO 230
MC = MB(1,MC+1)
GO TO 270
C-----
580 IF (MC . EQ . MAC) GO TO 200
MC = MB(2,MC+1)
GO TO 220
C-----
590 MC = MB(1,MC+1)
GO TO 240
C
C----------------------------------------------------------------------
C------------------------------------------------------------ TRIPLETS
C----------------------------------------------------------------------
600 PIADD = 0.
FLAGT = .FALSE.
IF (J.EQ.0) GO TO 850
JM = 1
630 MC = JJ(JM)
NB2 = MB(1,MC)
NY2 = MB(2,MC)
RT2 = TIDRA(LDT(NB2))
MC = MB(2,MC+2)
IF (MC.EQ.0) GO TO 740
C -----------------------------------------------------------------
C --------------------------- A : WIRES IN THREE DIFFERENT LAYERS
C -----------------------------------------------------------------
IDTFL = -200
ASSIGN 691 TO LABELT
CUTPHI = CUTP(1)
CUTV = CUTV1
640 NB3 = MB(1,MC+1)
NY3 = MB(2,MC+1)
RT3 = TIDRA(LDT(NB3))
NT3 = NB3
650 NT2 = NB2
660 NUSED = -2
IF (IDT(NT2).GT.IDTFL) NUSED = -1
IF (IDT(NT3).GT.IDTFL) NUSED = NUSED + 1
IF (NUSED.EQ.0) GO TO 692
NT1 = NB
690 IF (IDT(NT1).GE.IDTFL) NUSED = NUSED + 1
IF (NUSED.LT.0) GO TO 695
691 IF (NT1.EQ.NY) GO TO 692
NT1 = NT1 + 1
GO TO 690
C-----
692 NT2 = NT2 + 1
IF (NT2.LE.NY2) GO TO 660
693 NT3 = NT3 + 1
IF (NT3.LE.NY3) GO TO 650
MC = MB(1,MC)
IF (MC.NE.0) GO TO 640
GO TO 740
C-----
695 PD = DTP(NT2) - DTP(NT1) + PB(JM)
PD2 = DTP(NT3) - DTP(NT2) + BM(2,MC)
PD3 = PD + PD2
IF (ABS(PD3).GT.CUTP(2)) GO TO 691
B1 = (RT2 - R)**2 + R * RT2 * (PD*PD)
B2 = (RT3 - RT2)**2 + RT2 * RT3 * (PD2*PD2)
B3 = (R - RT3)**2 + R * RT3 * PD3 * PD3
CDP = (B3 - B1 - B2) * 0.5 / (B1 * B2)
IF (CDP.LE.0.) GO TO 691
SDP = (R * RT2 * PD * (1. - (PD * PD ) * .177778) +
+ RT2 * RT3 * PD2 * (1. - (PD2 * PD2) * .177778) -
+ R * RT3 * PD3 * (1. - (PD3 * PD3) * .177778)) / (B1 * B2)
GB1(1) = RT2 - R
GB1(2) = GB1(1)
GB2(1) = RT3 - RT2
GB2(2) = GB2(1)
700 IF (LWMDTC.GE.LWMDTL) GO TO 913
701 FLSOL = .FALSE.
FLSOL1 = .FALSE.
G1_cut = G1_frac * max(cell_size(ldt(nt1)),cell_size(ldt(nt2)))
G2_cut = G2_frac * max(cell_size(ldt(nt2)),cell_size(ldt(nt3)))
DT1 = DTD(NT1)
DT2 = DTD(NT2)
DT3 = DTD(NT3)
F1 = DT1 - DT2
V1 = PD * R - F1
S1 = B1 - F1 * F1
if(cont_fract) then
GB1(1) = 0.5 * (GB1(1) + S1 / GB1(1))
GB1(1) = 0.5 * (GB1(1) + S1 / GB1(1))
GB1(1) = 0.5 * (GB1(1) + S1 / GB1(1))
else
if(s1.lt.0.) then
cw write(6,'('' !!! 701 s1,gb1(1) '',2g12.5)') s1,gb1(1)
call INERR(212,1)
s1=abs(s1)
go to 730
else
GB1(1)=sqrt(s1) ! I switched off sqrt calc. via cont.fract.
endif
endif
G1 = GB1(1)
F2 = DT2 - DT3
V2 = PD2 * RT2 - F2
S1 = B2 - F2 * F2
if(cont_fract) then
GB2(1) = 0.5 * (GB2(1) + S1 / GB2(1))
GB2(1) = 0.5 * (GB2(1) + S1 / GB2(1))
GB2(1) = 0.5 * (GB2(1) + S1 / GB2(1))
else
if(s1.lt.0.) then
cw write(6,'('' !!! 701 s1,gb2(1) '',2g12.5)') s1,gb2(1)
call INERR(213,1)
s1=abs(s1)
go to 730
else
GB2(1)=sqrt(s1) ! I switched off sqrt calc. via cont.fract.
endif
endif
G2 = GB2(1)
C -----------------------------------------------------------------
C -------------------------------------------------- 8 AMBIGUITIES
C -----------------------------------------------------------------
ASSIGN 702 TO LABEL 8
GO TO 710
C-----
702 DT1 = - DT1
F1 = DT1 - DT2
V1 = PD * R - F1
S1 = B1 - F1 * F1
if(cont_fract) then
GB1(2) = 0.5 * (GB1(2) + S1 / GB1(2))
GB1(2) = 0.5 * (GB1(2) + S1 / GB1(2))
GB1(2) = 0.5 * (GB1(2) + S1 / GB1(2))
else
if(s1.lt.0.) then
cw write(6,'('' !!! 702 s1,gb1(2) '',2g12.5)') s1,gb1(2)
call INERR(214,1)
s1=abs(s1)
go to 730
else
GB1(2)=sqrt(s1) ! I switched off sqrt calc. via cont.fract.
endif
endif
G1 = GB1(2)
ASSIGN 703 TO LABEL 8
GO TO 710
C-----
703 FLSOL1 = .FALSE.
DT3 = - DT3
F2 = DT2 - DT3
V2 = PD2 * RT2 - F2
S1 = B2 - F2 * F2
if(cont_fract) then
GB2(2) = 0.5 * (GB2(2) + S1 / GB2(2))
GB2(2) = 0.5 * (GB2(2) + S1 / GB2(2))
GB2(2) = 0.5 * (GB2(2) + S1 / GB2(2))
else
if(s1.lt.0.) then
cw write(6,'('' 703 s1,gb2(2):'',2g12.5)') s1,gb2(2)
call INERR(215,1)
s1=abs(s1)
go to 730
endif
GB2(2)=sqrt(s1) ! I switched off sqrt calc. via cont.fract.
endif
G2 = GB2(2)
ASSIGN 704 TO LABEL 8
GO TO 710
C-----
704 DT1 = - DT1
F1 = DT1 - DT2
V1 = PD * R - F1
G1 = GB1(1)
ASSIGN 705 TO LABEL 8
GO TO 710
C-----
705 FLSOL1 = .FALSE.
DT2 = - DT2
F1 = DT1 - DT2
V1 = PD * R - F1
G1 = GB1(2)
F2 = DT2 - DT3
V2 = PD2 * RT2 - F2
G2 = GB2(1)
ASSIGN 706 TO LABEL 8
GO TO 710
C-----
706 DT1 = - DT1
F1 = DT1 - DT2
V1 = PD * R - F1
G1 = GB1(1)
ASSIGN 707 TO LABEL 8
GO TO 710
C-----
707 FLSOL1 = .FALSE.
DT3 = - DT3
F2 = DT2 - DT3
V2 = PD2 * RT2 - F2
G2 = GB2(2)
ASSIGN 708 TO LABEL 8
GO TO 710
C-----
708 DT1 = - DT1
F1 = DT1 - DT2
V1 = PD * R - F1
G1 = GB1(1)
ASSIGN 730 TO LABEL 8
710 IF (V1 * V2 .LT.CUTV) GO TO 720
S1 = G1 + G2
C1 = ((G1*F2 - F1*G2) * CDP -
+ (F1*F2 + G1*G2) * SDP) / S1
IF (ABS(C1).GT.CUR1) GO TO 720
IF (G1.LE.G1_cut) GO TO 711
IF (G2.GT.G2_cut) GO TO 712
711 IF (ABS(C1).GT.CUR2) GO TO 720
712 S1 = (S1*S1) * (G1*G1) * (G2*G2) /
+ (DTE(NT1)*(G2*G2) + DTE(NT3)*(G1*G1) + DTE(NT2)*(S1*S1))
C ------------------------------------------------ FILL K - BUFFER
IF (.NOT.FLSOL1) GO TO 715
LWKDTC = LWKDTC + 1
KB(1,LWKDTC) = LWKDTC - 2
GO TO 716
C-----
715 LWKDTC = LWKDTC + 3
IF (LWKDTC.GT.LWKDT) GO TO 960
KB(1,LWKDTC-2) = 0
KB(2,LWKDTC-2) = NT3
BK(3,LWKDTC-2) = DT3
dtds(nt3)=dt3
KB(4,LWKDTC-2) = 0
KB(1,LWKDTC-1) = LWKDTC - 2
KB(2,LWKDTC-1) = NT2
BK(3,LWKDTC-1) = DT2
dtds(nt2)=dt2
BK(4,LWKDTC-1) = G2
KB(1,LWKDTC) = LWKDTC - 1
716 KB(2,LWKDTC) = NT1
BK(3,LWKDTC) = DT1
dtds(nt1)=dt1
BK(4,LWKDTC) = G1
C ---------------------------------------------------------- FLAGS
IDT(NT3) = MAX0 (IDT(NT3),IDTFL)
IDT(NT2) = MAX0 (IDT(NT2),IDTFL)
IDT(NT1) = MAX0 (IDT(NT1),IDTFL)
C ------------------------------------------------ FILL TEMP BUFFER
MB(2,LWMDTC) = LWKDTC
BM(2,LWMDTC+1) = PD3
BM(1,LWMDTC+2) = B3
BM(2,LWMDTC+2) = V1
BM(1,LWMDTC+3) = C1
BM(2,LWMDTC+3) = S1
MB(1,LWMDTC+4) = 0
MB(2,LWMDTC+4) = LDT(NT2)
BM(1,LWMDTC+5) = 0.
BM(2,LWMDTC+5) = PD
C --------------------------------------------- SET POINTERS IN MB
IF (FLSOL) GO TO 718
MB(1,LWMDTC) = MB(1,M+2)
MB(1,LWMDTC+1) = 0
MB(1,M+2) = LWMDTC
GO TO 719
C-----
718 MB(1,LWMDTC) = MB(1,MB(1,M+2))
MB(1,MB(1,M+2)) = 0
MB(1,LWMDTC+1) = MB(1,M+2)
MB(1,M+2) = LWMDTC
719 FLSOL = .TRUE.
FLSOL1 = .TRUE.
LWMDTC = LWMDTC + 6
720 GO TO LABEL 8, (702,703,704,705,706,707,708,730)
730 CONTINUE
IF (FLSOL . AND . NT1.NE.NTX) GO TO 500
GO TO LABELT, (691,758,818,867)
C
C -----------------------------------------------------------------
C --------------------- TRIPLETS B : W2 AND W3 IN THE SAME LAYER
C----------------------------------------------------------------------
740 CUTPHI = CUTP(2)
CUTV = CUTV0
IF (JM.GT.J2) GO TO 840
IF (NB2.EQ.NY2) GO TO 800
ASSIGN 758 TO LABELT
IDTFL = JM - 400
NT2 = NB2
NTX2 = NY2
NTD = -1
NTX = NB
NTA = NY
750 NT3 = NT2 - NTD
IF (IDT(NT2).LE.IDTFL) GO TO 751
IF (IDT(NT3).GT.IDTFL) GO TO 759
751 IF (IDT(NT2).GE.-1) GO TO 759
IF (IDT(NT3).GE.-1) GO TO 759
PD2 = DTP(NT3) - DTP(NT2)
NT1 = NTA
752 IF (IDT(NT1).GT.IDTFL) GO TO 758
c absence of this line seems to be a bug (A.Andryakov)
PD2 = DTP(NT3) - DTP(NT2)
c
PD = DTP(NT2) - DTP(NT1) + PB(JM)
IF (ABS(PD).GT.CUTP3) GO TO 758
PD3 = PD + PD2
B1 = (RT2 - R)**2 + RT2 * (R*PD) * PD
B2 = (RT2 * PD2)**2
B3 = (R-RT2)**2 + R * RT2 * PD3 * PD3
C
C THE FOLLOWING SHOULD NEVER HAPPEN.....WHAT DOES IT MEAN??
IF (B1.EQ.0.) THEN
B1 = 1.
CALL ERLOGR('DTPTRNS',ERWARN,0,0,
+ 'B1=0. AT LINE 1739')
ENDIF
IF (B2.EQ.0.) THEN
B2 = 1.
CALL ERLOGR('DTPTRNS',ERWARN,0,0,
+ 'B2=0. AT LINE 1739')
ENDIF
CDP = (B3 - B1 - B2) * 0.5 / (B1 * B2)
IF (CDP.LE.CUTC) GO TO 758
C
SDP = RT2 * (R * (PD * (1. - PD * PD * .177778) -
+ PD3 * (1. - PD3 * PD3 * .177778)) +
+ RT2 * PD2 * (1. - PD2 * PD2 * .177778)) / (B1 * B2)
GB1(1) = RT2 - R
GB1(2) = GB1(1)
GB2(1) = max(cell_size(ldt(nt1)),cell_size(ldt(nt2)))
GB2(2) = 0.5*GB2(1)/3.
GO TO 700
C-----
758 IF (NT1.EQ.NTX) GO TO 759
NT1 = NT1 + NTD
GO TO 752
C-----
759 IF (NT3.EQ.NTX2) GO TO 760
NT2 = NT3
GO TO 750
C-----
760 IF (NTD.GT.0) GO TO 800
NT2 = NY2
NTX2 = NB2
NTD = 1
NTX = NY
NTA = NB
GO TO 750
C
C -----------------------------------------------------------------
C -------------------- TRIPLETS C : W1 AND W2 IN THE SAME LAYER
C----------------------------------------------------------------------
800 IF (NB.EQ.NY) GO TO 840
ASSIGN 818 TO LABELT
IDTFL = -600
NT3 = NB2
NTX2 = NY2
NTD = -1
NTX = NB
NTA = NY
810 IF (IDT(NT3).GT.-300) GO TO 819
NT2 = NTA
811 NT1 = NT2 + NTD
IF (IDT(NT1).GE.-1) GO TO 818
IF (IDT(NT2).GE.-1) GO TO 818
IF (IDT(NT1).GE.JM-400) GO TO 812
IF (IDT(NT2).LT.JM-400) GO TO 813
812 IF (IDT(NT1).LE.IDTFL) GO TO 813
IF (IDT(NT2).GT.IDTFL) GO TO 818
813 PD2 = DTP(NT3) - DTP(NT2) + PB(JM)
IF (ABS(PD2).GT.CUTP3) GO TO 818
PD = DTP(NT2) - DTP(NT1)
PD3 = PD + PD2
B1 = RSQ * (PD*PD)
B2 = (RT2 - R)**2 + R * RT2 * PD2 * PD2
B3 = (RT2-R)**2 + R * (RT2*PD3) * PD3
CDP = (B3 - B1 - B2) * 0.5 / (B1 * B2)
IF (CDP.LE.CUTC) GO TO 818
SDP = R * (RT2 * (PD2 * (1. - PD2 * PD2 * .177778) -
+ PD3 * (1. - PD3 * PD3 * .177778)) +
+ R * PD * (1. - PD * PD * .177778)) / (B1 * B2)
GB1(1) = max(cell_size(ldt(nt1)),cell_size(ldt(nt3)))
GB1(2) = 0.5*GB1(1)/3.
GB2(1) = RT2 - R
GB2(2) = GB1(1)
GO TO 700
C-----
818 IF (NT1.EQ.NTX) GO TO 819
NT2 = NT1
GO TO 811
C-----
819 IF (NT3.EQ.NTX2) GO TO 820
NT3 = NT3 - NTD
GO TO 810
C-----
820 IF (NTD.GT.0) GO TO 840
NT3 = NY2
NTX2 = NB2
NTD = 1
NTX = NY
NTA = NB
GO TO 810
C-----
840 JM = JM + 1
IF (JM.LE.J) GO TO 630
C -----------------------------------------------------------------
C ------------------------------- D : ALL WIRES IN THE SAME LAYER
C -----------------------------------------------------------------
850 IF (NY-NB.LE.1) GO TO 900
ASSIGN 867 TO LABELT
CUTPHI = CUTP(1)
CUTV = CUTV0
IDTFL = -900
NT3 = NB
NTX2 = NY - 2
NTD = 1
NTX = NY
NTA = 2
RT2 = R
860 IF (IDT(NT3).GE.-1) GO TO 870
NT1 = NT3 + NTA
IF (IDT(NT1).GE.-1) GO TO 869
861 NT2 = NT3 + NTD
IF (IDT(NT2).GE.-1) GO TO 868
862 NUSED = 0
IF (IDT(NT3).GE.IDTFL) NUSED = NUSED + 1
IF (IDT(NT2).GE.IDTFL) NUSED = NUSED + 1
IF (IDT(NT1).GE.IDTFL) NUSED = NUSED + 1
IF (NUSED.GT.1) GO TO 868
PD = DTP(NT2) - DTP(NT1)
PD2 = DTP(NT3) - DTP(NT2)
PD3 = PD + PD2
IF (ABS(PD3).GT.CUTP3) GO TO 868
GB1(1) = cell_size(ldt(nt1))
GB2(1) = GB1(1)
GB1(2) = 0.5*GB1(1)/3.
GB2(2) = GB1(2)
B1 = (RSQ*PD) * PD
B2 = RSQ * PD2 * PD2
B3 = RSQ * PD3 * PD3
CDP = (B3 - B1 - B2) * 0.5 / (B1 * B2)
SDP = RSQ * (PD * (1. - PD * PD * .177778) +
+ PD2 * (1. - PD2 * PD2 * .177778) -
+ PD3 * (1. - PD3 * PD3 * .177778)) / (B1 * B2)
GO TO 700
C-----
867 IDTFL = IDTFL - 1
868 NT2 = NT2 + NTD
IF (NT2.NE.NT1) GO TO 862
869 IF (NT1.EQ.NTX) GO TO 870
NT1 = NT1 + NTD
GO TO 861
C-----
870 IF (NT3.EQ.NTX2) GO TO 880
NT3 = NT3 + NTD
GO TO 860
C-----
880 IF (NTD.LT.0) GO TO 900
IDTFL = -800
NT3 = NY
NTX2 = NB + 2
NTD = -1
NTX = NB
NTA = -2
GO TO 860
C
C-----------------------------------------------------------------------
C----------------------------------------------------------- NEXT WIRE
C-----------------------------------------------------------------------
900 MB(1,M+1) = LWMDTC
MB(2,M+1) = MR
MR = M
M = LWMDTC
GO TO 100
C
C ------------------------------------------------------------------
C ------------------------------------------- OVERWRITE M - BUFFER
C-----------------------------------------------------------------------
910 LABELM = LABELM + 1
911 LABELM = LABELM + 1
912 LABELM = LABELM + 1
913 LABELX = LABELM
LABELM = 1
IF (LWMDTL.NE.LWMDT) GO TO 930
LWMDTC = 1
920 LWMDTL = MAZ(3,LXTR+1) - 16
GO TO (701,560,370,210), LABELX
C-----
930 IF (LWMDTL-MAZ(3,LXTR+1) + 16) 920,950,940
940 LWMDTL = LWMDT
GO TO (701,560,370,210), LABELX
C
C ------------------------------------------------------------------
C ------------------------------------------------ BUFFER OVERFLOW
C-----------------------------------------------------------------------
950 CALL INERR (203,1)
RETURN
C-----
960 CALL INERR (204,1)
RETURN
C-----
970 CALL INERR (205,1)
RETURN
9410 FORMAT (5I5,L5,9I5,L5)
END
C 27/11/89 MEMBER NAME DTSORT (ES) FVS
SUBROUTINE DTSORT(i_view)
$$IMPLICIT
C
C SELECT THE BEST TRACK CANDIDATES FOUND BY DTPTRN
C
C --- HARTWIG ALBRECHT --- 80-05
C --- HARTWIG ALBRECHT --- 82-01
C --- HARTWIG ALBRECHT --- 83-10
C
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:CDTHIT.INC'
$$INCLUDE 'K$ITRK:CDTCUT.INC'
$$INCLUDE 'K$ITRK:CTRACK.INC'
$$INCLUDE 'K$ITRK:CUWORK.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
INTEGER I_VIEW
INTEGER NWORK
COMMON /CUWORK/ NWORK(512)
INTEGER I_VIEW_S,K,I,I1,I2,I3,I4,I0,KC,LF,N,LL
C
IF (LWIDTC.EQ.0) RETURN
if(tidst(i_view).lt.0) then
i_view_s = 6
else
i_view_s = 4
endif
C------------------------------------------------- SORT WORKING SPACE
K = 0
DO 100 I=1,LWIDTC
IF (IB(2,I).NE.0 . OR . IB(3,I).EQ.0) GO TO 100
K = K + 1
NWORK(K) = I
100 CONTINUE
IF (K.EQ.0) RETURN
C--------------------------------------------------------------- SORT
I1 = K
110 I1 = I1 / 2
IF (I1.EQ.0) GO TO 140
I2 = K - I1
DO 130 I3=1,I2
I4 = I3
120 IF (I4.LT.1) GO TO 130
IF (BI(8,NWORK(I4+I1)).GE.BI(8,NWORK(I4))) GO TO 130
I0 = NWORK(I4)
NWORK(I4) = NWORK(I4+I1)
NWORK(I4+I1) = I0
I4 = I4 - I1
GO TO 120
130 CONTINUE
GO TO 110
C-------------------------------------------------------- SET POINTERS
140 IF (LWIDTL.NE.0) GO TO 150
LWIDT1 = NWORK(1)
GO TO 160
150 IB(1,LWIDTL) = NWORK(1)
I1 = LWIDTL
160 IF (K.EQ.1) GO TO 180
DO 170 I=2,K
170 IB(1,NWORK(I-1)) = NWORK(I)
180 IB(1,NWORK(K)) = 0
C-------------------------------------------------- SELECT GOOD TRACKS
I = NWORK(1)
200 I2 = IB(1,I)
IF (IB(2,I).NE.0) GO TO 330
IB(5,I) = i_view
IB(i_view_s,I) = 0
KC = 0
K = IB(3,I)
LF = 0
210 N = KB(2,K)
IF (IDT(N).GT.0) GO TO 250
IF (LF.EQ.0) LF = LDT(N)
IB(i_view_s,I) = IB(i_view_s,I) + 1
KC = K
K = KB(1,K)
IF (K.NE.0) GO TO 210
GO TO 300
C----------------------------------------------- WIRE ALREADY ASSIGNED
250 BI(8,I) = BI(8,I) + 1.
IF (KC.EQ.0) GO TO 260
BK(4,KC) = BK(4,KC) + BK(4,K)
K = KB(1,K)
KB(1,KC) = K
IF (K.NE.0) GO TO 210
GO TO 300
260 IB(3,I) = KB(1,K)
K = KB(1,K)
IF (K.NE.0) GO TO 210
300 LL = LDT(N)
c IF (LL.EQ.LF) GO TO 400
IF (IB(i_view_s,I).LT.4) GO TO 400
IF (I2.EQ.0) GO TO 310
IF (BI(8,I).GT.BI(8,I2)) GO TO 340
C--------------------------------------------- MARK THE ASSIGNED WIRES
310 LWJDTC = LWJDTC + 1
LWIDTL = I
IB(2,I) = LWJDTC
JB(1,LWJDTC) = 0
JB(9,LWJDTC) = 1
c IB(4,I) = 0
c IB(6,I) = 0
K = IB(3,I)
320 IDT(KB(2,K)) = I
K = KB(1,K)
IF (K.NE.0) GO TO 320
330 I1 = I
GO TO 430
C------------------------------------------------------------- REORDER
340 I3 = IB(1,I)
350 I4 = I3
I3 = IB(1,I4)
IF (I3.EQ.0) GO TO 360
IF (BI(8,I).GT.BI(8,I3)) GO TO 350
360 IB(1,I4) = I
IB(1,I) = I3
GO TO 410
400 IB(3,I) = 0
410 IF (I1.NE.0) GO TO 420
LWIDT1 = I2
GO TO 430
420 IB(1,I1) = I2
430 IF (LWJDTC.GE.LWJDT) GO TO 600
I = I2
IF (I.NE.0) GO TO 200
500 IF (LWIDT1.LE.0) LWIDTC = 0
RETURN
600 CALL INERR (206,1)
IB(1,I) = 0
LWIDTL = 0
GO TO 500
END
C 27/11/89 MEMBER NAME DTDIST (ES) FVS
SUBROUTINE DTDIST
$$IMPLICIT
C
C H.SCHROEDER VERSION 20-03-81
C H.ALBRECHT VERSION 26-01-82
C H.SCHROEDER VERSION 05-03-83
C H.SCHROEDER VERSION 28-08-83
C-----CALCULATES THE TRACK PARAMETERS :
C----- D : DISTANCE OF CLOSEST APPROACH
C----- X0: X VALUE AT D
C----- Y0: Y VALUE AT D
C----- XR: X VALUE AT FIRST LAYER WITH RADIUS R
C----- YR: Y VALUE AT FIRST LAYER WITH RADIUS R
C-----THESE VALUES ARE STORED IN THE BJ-BANK
C----------------------------------------------------------------------
C- Updated 3-DEC-1993 Alexander Andryakov
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CTRACK.INC'
$$INCLUDE 'K$ITRK:CDTHIT.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
C
c COMMON /CUWORK/ DTXH(256),DTYH(256)
INTEGER I,ILAST,J,NHIT,K1,K2,NHIM,IN,IC,I2
REAL R2,P2,R,P,PD,SP,CP,X1,Y1,CPD,X2,Y2
REAL CEN,RHONEW,RHOSQ,QD,X0,Y0,RADSQ,RT
REAL XQ,YQ,T0X,T0Y,TQX,TQY,TTQ0,Q,D,PSI,SPP,CPP
REAL D1,D2,FAKY,FAKX,DTXH,DTYH,XP,YP,XC,YC
dimension DTXH(256),DTYH(256)
DOUBLE PRECISION XA,YB,DET,SXX,SXXP,SXY,SXPXP,SXPY,SIGMA,FY
+,DX1,DY1
C
C-----LOOP OVER ALL TRACKS
IF (LWIDTC.EQ.0) RETURN
I = LWIDT1
ILAST = 0
10 CONTINUE
J = IB(2,I)
if(j.ge.LWJDT) write (6,'('' @@@ DTDIST @@@ j ='',i4)') j
IF (JB(9,J).EQ.0) GO TO 90
C
C-----LOOP OVER THE WIRES OF THE TRACK
NHIT = 0
K1 = IB(3,I)
20 K2 = KB(1,K1)
IF (K2.EQ.0) GO TO 30
R2 = TIDRA(LDT(KB(2,K2)))
P2 = DTP(KB(2,K2))
R = TIDRA(LDT(KB(2,K1)))
P = DTP(KB(2,K1))
PD = P2 - P
IF (PD.GT.PI) PD = PD - PI2
IF (PD.LT.-PI) PD = PD + PI2
SP = SIN(P)
CP = COS(P)
X1 = R * CP
Y1 = R * SP
CPD = 1. - 0.5*PD*PD
X2 = R2 * (CPD*CP-SP*PD)
Y2 = R2 * (CPD*SP+CP*PD)
D1 = BK(3,K1)
D2 = BK(3,K2)
FAKY = (R - R2)**2 + R * R2 * PD * PD
FAKX = ((D2-D1) * (X2-X1) + BK(4,K1) * (Y2-Y1)) / FAKY
FAKY = (BK(4,K1) * (X2-X1) - (D2-D1) * (Y2-Y1)) / FAKY
NHIT = NHIT + 2
if(nhit.gt.256) then
nhit = 256
call INERR(207,1)
endif
NHIM = NHIT-1
DTXH(NHIM) = X1 - D1 * FAKX
DTYH(NHIM) = Y1 + D1 * FAKY
DTXH(NHIT) = X2 - D2 * FAKX
DTYH(NHIT) = Y2 + D2 * FAKY
c
c Save x,y of the hit in /CDTHIT/ for 3-d pt.r.
c
DTXH_K(LIST_HIT(KB(2,K1))) = DTXH(NHIM)
DTYH_K(LIST_HIT(KB(2,K1))) = DTYH(NHIM)
DTXH_K(LIST_HIT(KB(2,K2))) = DTXH(NHIT)
DTYH_K(LIST_HIT(KB(2,K2))) = DTYH(NHIT)
c
c
c
25 CONTINUE
K1 = KB(1,K2)
IF (K1.NE.0) GO TO 20
30 CONTINUE
C
C-----DETERMINE FIXED POINT ON THE TRACK
IN = NHIT
XP = DTXH(IN)
YP = DTYH(IN)
SXX = 0.
SXXP = 0.
SXY = 0.
SXPXP = 0.
SXPY = 0.
c write (6,'('' dtdist track can :'',i4/'' x y'')') i
DO 50 IC=1,NHIT
c write (6,'(2g11.4)') DTXH(IC),DTYH(IC)
DX1 = DTXH(IC) - XP
DY1 = DTYH(IC) - YP
IF ( IC.EQ.IN ) GOTO 50
FY = DX1**2 + DY1**2
IF (FY.EQ.0.) GOTO 50
SIGMA = 0.0001*FY
SXX = SXX + DX1*DX1*SIGMA
SXXP = SXXP + DX1*DY1*SIGMA
SXY = SXY + DX1*FY*SIGMA
SXPXP = SXPXP + DY1*DY1*SIGMA
SXPY = SXPY + DY1*FY*SIGMA
50 CONTINUE
DET = SXX*SXPXP - SXXP*SXXP
XA = -0.5 * ( SXY*SXPXP - SXXP*SXPY ) / DET
YB = -0.5 * ( SXX*SXPY -SXY*SXXP ) / DET
XC = XP - XA
YC = YP - YB
CEN = SQRT( XC*XC + YC*YC )
RHONEW = DSQRT( XA*XA + YB*YB )
RHOSQ = RHONEW*RHONEW
QD = RHONEW - CEN
X0 = - QD*XC/CEN
Y0 = - QD*YC/CEN
I2 = NHIT/2
RADSQ = DTXH(I2)**2 + DTYH(I2)**2
RT = RHOSQ - RADSQ - CEN*CEN
DX1 = 0.5 * DTYH(I2) * ( RT + 2.*XC*DTXH(I2) +2.*YC*DTYH(I2) ) /
+ ( YC*DTXH(I2) - XC*DTYH(I2) )
DY1 = - DTXH(I2) * DX1 / DTYH(I2)
XQ = DTXH(I2) + DX1
YQ = DTYH(I2) + DY1
T0X = X0 - XC
T0Y = Y0 - YC
TQX = XQ - XC
TQY = YQ - YC
TTQ0 = T0Y*TQX - T0X*TQY
Q = SIGN(1.,TTQ0)
D = Q*QD
C
C-----FILL J-BUFFER
C
BJ(2,J) = D
RADSQ = DTXH(1)**2 + DTYH(1)**2
RT = RHOSQ - RADSQ - CEN*CEN
DX1 = 0.5 * DTYH(1) * ( RT + 2.*XC*DTXH(1) +2.*YC*DTYH(1) ) /
+ ( YC*DTXH(1) - XC*DTYH(1) )
IF (DABS(DX1).GT.0.3) DX1=0.
DY1 = - DTXH(1) * DX1 / DTYH(1)
BJ(4,J) = DTXH(1) + DX1
BJ(5,J) = DTYH(1) + DY1
BJ(3,J) = 0.
BJ(6,J) = 0.
BJ(8,J) = 0.
BJ(10,J) = (BJ(4,J) -X0 )**2 + (BJ(5,J) - Y0)**2
PSI = .25 * BJ(10,J) / RHOSQ
PSI = RHONEW * PSI * (0.5 + 0.125 * PSI)
BJ(10,J) = SQRT( BJ(10,J) + (5.33333 * PSI * PSI) )
PSI = BJ(10,J) / RHONEW
SPP = Q*(XC-BJ(4,J))/CEN
CPP = -Q*(YC-BJ(5,J))/CEN
BJ(7,J) = ATAN2 (SPP,CPP)
JB(9,J) = 1
ILAST = I
C
C------------ for debugging only -------------
C
c BI(8,I) = BI(7,I)
C
C------------ for debugging only -------------
C
BI(7,I) = 0.5 * Q / RHONEW
90 I = IB(1,I)
IF (I.NE.0) GO TO 10
RETURN
END
[KLOE]
[Offline Doc]
[TRK Files]
Generated with Light on Thu Apr 8 13:00:16 MET DST 1999
.
Mail comments and suggestions.
