*-- Author : I.O. Skillicorn
SUBROUTINE FTFIT
c ftfit4.f on graphics
c ftfit6.f on h1rec
c see later version ftfit7.f for isolation flags
C fit LINE SEGMENTS in phi-z,r-z
C use mean r in phi calculation
c modify to compare and join segments if
c consistent in phi/r at centre of radial chamber.
c objective: join segments that cross cell boundary.
c
c add section to check if a better line segment can be obtained by
c reversing drift sign and rejecting the worst point.
c
C AUTHOR I.O.SKILLICORN
*KEEP,FRDIMS.
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)
*KEEP,FH1WORK.
COMMON/FGMIOS/
* Planar geometry
+ ZPP(140),C(140),S(140),ITYPE(140),WZERO(140),WSPACE,
*
* Radial geometry
+ ZP(36),PHW(36),WS(36)
*
COMMON/H1WORK/
* Radial data...
+ WW(MAXHTS,36),DRI(MAXHTS,36),RM(MAXHTS,36),
+ NDP(36), NW(MAXHTS,36), DWS(MAXHTS,36),
*
* Planar Data
+ NDPW(NUMWPL),DW(MAXHTS,NUMWPL),
+ DRIW(MAXHTS,NUMWPL),NDW(MAXHTS,NUMWPL),
+ WWP(MAXHTS,NUMWPL),
+ IPHOLE(MAXHTS,NUMWPL),
*
* Pointers into DIGI bank for IOS labelled hits
+ IPFRRE(MAXHTS,36),IPFRPE(MAXHTS,36),NFRRE,NFRPE,
+ IRPIOS(MAXDIG,2), IPPIOS(MAXDIG,2),
*
* Track segment data
+ NTRAKS(3),IRPT(12,MAXTRK,3),SDRFT(12,MAXTRK,3),
*
* Fit data
+ PCOSL(MAXTRK,3),PSINL(MAXTRK,3),PHZL(MAXTRK,3),
+ DPCOSL(MAXTRK,3),DPSINL(MAXTRK,3),
+ DPHZL(MAXTRK,3),CHSQ(MAXTRK,3),RZI(MAXTRK,3),
+ RPCOSG(MAXTRK),RPSING(MAXTRK),
+ PHZG(MAXTRK),CC(3,MAXTRK),ZIG(MAXTRK),
+ IRADG(36,MAXTRK),PHIG(36,MAXTRK),
+ IG,SDRADG(36,MAXTRK),
+ R1,Z1,RFIT(MAXTRK,3),
+ CHG(MAXTRK),
+ PPA(MAXTRK,3), ZZA(MAXTRK,3),
+ GPA(MAXTRK,3),GZA(MAXTRK,3)
*
*
*KEEP,FPTVTX.
COMMON/VERTVV/ZV ,XVV,YVV
**the common/VERTEX/ becomes /VERTVV/ (in analogy to /VERTFF/) on the
** 17/6/91, since it is in conflict with the VERTEX module (g.bernardi)
** (note that all these common names should start by F in this deck...)
*KEND.
*
DIMENSION IPA(36),SD(36)
dimension xx(12),yy(12),zz(12),xxx(12),zzz(12)
data istart/0/
* segments to be joined agree in r within rcut
* agree in phi within pcut / r ie within 0.5 cm
* in drift. segments have no more than one z-plane in common
if(istart.eq.0)then
istart=1
call stext(3040,4,'ftfit rms original')
call bhs(3040,0,50,0.,0.1)
call stext(3041,4,'ftfit rms point rejected')
call bhs(3041,0,50,0.,0.1)
call stext(3042,4,'ftfit rms reversed sign')
call bhs(3042,0,50,0.,0.2)
call stext(3043,4,'ftfit rms reversed sign pt rej')
call bhs(3043,0,50,0.,0.2)
call stext(3044,4,'ftfit rms rev. sign pt rej sel')
call bhs(3044,0,50,0.,0.1)
call stext(3045,4,'ftfit res max original ')
call bhs(3045,0,50,0.,0.2)
call stext(3046,4,'ftfit res max reversed sel')
call bhs(3046,0,50,0.,0.2)
call stext(3047,4,'ftfit delta phi*r - proj')
call bhs(3047,0,50,-5.00,5.00)
call stext(3048,4,'ftfit delta phi - proj')
call bhs(3048,0,50,-0.25,0.25)
endif
rcut=20.
pcut=1.0
c write(*,*)' ftfit1 entered '
*
c check for wrong drift-sign line segments
do 100 ism=1,3
do 110 j=1,ntraks(ism)
ll=0
nww=0
do 120 k=1,12
i=irpt(k,j,ism)
if(i.eq.0)goto120
ip=k+(ism-1)*12
if(nww.eq.0)nww=nw(i,ip)
if(nww.ne.nw(i,ip))goto110
c same wire
ll=ll+1
xx(ll)=zp(ip)
yy(ll)= sdrft(k,j,ism)*dri(i,ip)+dws(i,ip)
zz(ll)=-sdrft(k,j,ism)*dri(i,ip)+dws(i,ip)
120 continue
if(ll.le.5)goto110
call ftlft(xx,yy,ll,0,at,bt,ee)
ee=sqrt(abs(ee))
call shs(3040,0,ee)
c if rms is good do not check
* if(ee.lt.0.04)goto110
c
c find worst point
c remove and refit
c this makes posible a direct comparison with
c the reversed sign remove and refit
dmax=0.
do 130 l=1,ll
diff=abs(yy(l)-(at*xx(l)+bt))
if(diff.gt.dmax)then
dmax=diff
lmax=l
endif
130 continue
call shs(3045,0,dmax)
lll=0
do 135 l=1,ll
if(l.eq.lmax)goto135
lll=lll+1
xxx(lll)=xx(l)
zzz(lll)=yy(l)
135 continue
call ftlft(xxx,zzz,lll,0,at1,bt1,ee1)
ee1=sqrt(abs(ee1))
call shs(3041,0,ee1)
c if rms is good do not check
if(ee1.lt.0.03)goto110
c
c fit oposite sign
call ftlft(xx,zz,ll,0,at2,bt2,ee2)
ee2=sqrt(abs(ee2))
call shs(3042,0,ee2)
c find worst point
dmax=0.
do 140 l=1,ll
diff=abs(zz(l)-(at2*xx(l)+bt2))
if(diff.gt.dmax)then
dmax=diff
lmax=l
endif
140 continue
lll=0
do 150 l=1,ll
if(l.eq.lmax)goto150
lll=lll+1
xxx(lll)=xx(l)
zzz(lll)=zz(l)
150 continue
call ftlft(xxx,zzz,lll,0,at3,bt3,ee3)
c drift sign may not be determined after point rejection
c so add 0.0001 to counter rounding errors
ee3=sqrt(abs(ee3)) + 0.0001
call shs(3043,0,ee3)
if(ee3.lt.ee1)then
call shs(3044,0,ee3)
call shs(3046,0,dmax)
c write(*,*)' ftfit4: k j ism res',lmax,j,ism,ee,ee1,ee2
c 1 ,ee3
c rejection of one point gives
c improved residuals with reversed drift sign.
c reset drift sign, remove this point
ll=0
do 160 k=1,12
i=irpt(k,j,ism)
if(i.eq.0)goto160
ll=ll+1
if(ll.eq.lmax)then
irpt(k,j,ism)=0
endif
sdrft(k,j,ism)=-sdrft(k,j,ism)
160 continue
endif
110 continue
100 continue
c
c
c end of section for correcting drift sign
c
c
DO 1001 I=1,3
FN=NTRAKS(I)
DO 1002 J=1,NTRAKS(I)
DO 1005 II=1,36
SD(II)=0.
1005 IPA(II)=0
IF (I.EQ.1) THEN
IFP=1
ILP=12
ELSE IF (I.EQ.2) THEN
IFP=13
ILP=24
ELSE
IFP=25
ILP=36
ENDIF
IC=1
DO 1003 K=IFP,ILP
SD(K)=SDRFT(IC,J,I)
IPA(K)=IRPT(IC,J,I)
IC=IC+1
1003 CONTINUE
*
*------------------------------------------
C PRINT5000,I,J,IPA
5000 FORMAT(2I5,3X,36I2 )
CDEB PRINT 5001,SD
5001 FORMAT(36F3.0)
*------------------------------------------
*
ZVV=ZV
NPLA=36
C1=0.
C2=0.
C3=0.
CALL FTFHQQ(NPLA,IPA,SD,ZVV,C1,C2,C3,
1 PCOS,PSIN,PHZ,DPCOS,DPSIN,DPHZ,RI,CH,IT)
C WRITE(*,*)PCOS,PSIN,PHZ,DPCOS,RI,R1
PCOSL(J,I)=PCOS
PSINL(J,I)=PSIN
PHZL(J,I)=PHZ
DPCOSL(J,I)=DPCOS
DPSINL(J,I)=DPSIN
DPHZL(J,I)=DPHZ
RZI(J,I)=RI
C STORE MEAN R OF LINE SEGMENT
RFIT(J,I)=R1
CHSQ(J,I)=CH
1002 CONTINUE
if(ntraks(i).le.1)goto1001
c compare line segments and join if compatible
if(i.eq.1)im=6
if(i.eq.2)im=18
if(i.eq.3)im=30
c write(*,*)' module ,# segs',i,ntraks(i)
do 10 j=1,ntraks(i)-1
if(chsq(j,i).gt.1000.)goto10
phi1=zp(im)*pcosl(j,i)+phzl(j,i)
if(phi1.gt.6.28318)phi1=phi1-6.28318
do 20 k=j+1,ntraks(i)
if(chsq(k,i).gt.1000.)goto20
c write(*,*)' rfit',rfit(j,i),rfit(k,i),k,j
if(abs(rfit(j,i)-rfit(k,i)).gt.rcut)goto20
phi2=zp(im)*pcosl(k,i)+phzl(k,i)
if(phi2.gt.6.28318)phi2=phi2-6.28318
rr=0.5*(rfit(j,i)+rfit(k,i))
dphi12=amod((phi1-phi2),6.28318)
c write(*,*)' phi ',phi1,phi2,dphi12,pcut/rr
if(abs(dphi12).gt.5.*pcut/rr)goto20
c similar r and phi
ifw1=0
ifw2=0
do 30 l=1,12
if(irpt(l,j,i).ne.0.and.ifw1.eq.0)ifw1=l
if(irpt(l,k,i).ne.0.and.ifw2.eq.0)ifw2=l
if(irpt(l,j,i).ne.0)ilw1=l
if(irpt(l,k,i).ne.0)ilw2=l
30 continue
c write(*,*)' wires f,l 1,2',ifw1,ilw1,ifw2,ilw2
jj=irpt(ifw1,j,i)
kk=irpt(ifw2,k,i)
jp=ifw1+(i-1)*12
kp=ifw2+(i-1)*12
c check adjacent wedges
c write(*,*)' nw ',nw(jj,jp),nw(kk,kp)
if(nw(jj,jp).eq.1.and.nw(kk,kp).eq.48)goto60
if(nw(kk,kp).eq.1.and.nw(jj,jp).eq.48)goto60
if(iabs(nw(jj,jp)-nw(kk,kp)).ne.1)goto20
60 continue
if(ilw1.le.ifw2.or.ilw2.le.ifw1)then
c overlap by at most one hit
c histogram here
call shs(3047,0,dphi12*rr)
call shs(3048,0,dphi12)
if(abs(dphi12).gt.pcut/rr)goto20
do40 ii=1,36
sd(ii)=0.
40 ipa(ii)=0
c write(*,*)' JOIN '
c print 3001,i,j,(irpt(ll,j,i),ll=1,12),chsq(j,i)
c print 3001,i,k,(irpt(ll,k,i),ll=1,12),chsq(k,i)
3001 format(' trk, irpt ',2i3,4x,12i2,f5.0)
ic=1
do 41 ii=12*(i-1)+1,12*i
if(irpt(ic,j,i).ne.0)sd(ii)=sdrft(ic,j,i)
if(irpt(ic,j,i).ne.0)ipa(ii)=irpt(ic,j,i)
if(irpt(ic,k,i).ne.0)sd(ii)=sdrft(ic,k,i)
if(irpt(ic,k,i).ne.0)ipa(ii)=irpt(ic,k,i)
ic=ic+1
41 continue
c write(*,*)' module ',i
c print 3000,ipa
3000 format(' ipa ',36i2)
ZVV=ZV
NPLA=36
C1=0.
C2=0.
C3=0.
CALL FTFHQQ(NPLA,IPA,SD,ZVV,C1,C2,C3,
1 PCOS,PSIN,PHZ,DPCOS,DPSIN,DPHZ,RI,CH,IT)
C WRITE(*,*)PCOS,PSIN,PHZ,DPCOS,RI,R1
PCOSL(J,I)=PCOS
PSINL(J,I)=PSIN
PHZL(J,I)=PHZ
DPCOSL(J,I)=DPCOS
DPSINL(J,I)=DPSIN
DPHZL(J,I)=DPHZ
RZI(J,I)=RI
C STORE MEAN R OF LINE SEGMENT
RFIT(J,I)=R1
CHSQ(J,I)=CH
c replace points
do 42 l=1,12
irpt(l,j,i) =ipa(l+12*(i-1))
sdrft(l,j,i)=sd(l+12*(i-1))
irpt(l,k,i)=0
42 continue
c remove joined line segment
chsq(k,i)=5000.
c print 3001,i,j,(irpt(ll,j,i),ll=1,12),chsq(j,i)
c print 3001,i,k,(irpt(ll,k,i),ll=1,12),chsq(k,i)
endif
20 continue
10 continue
1001 CONTINUE
RETURN
END
*
*
*
*