C*HE 03/22/95 17:52:39
C*DK DFMAIN
C**********************************************************************
C The track fit package of ARGUS collaboration
C Written primarily by Hartwig Albrecht in 1982 and then updated by
C him and R.Kutschke, H.Kapitza, Ubi, Roland Waldi
C The description of the algorithm can be found in NIM A275, p.1-48.
C and in unpublished and unfinished version of H.Kapitza's note
C (at LNF vax cluster d19:[alex.argus.write-up]).
C Formulas from NIM paper and Kapitza note are referenced in the code.
C The original ARGUS code has been obtained from H.Kapitza
C-----------------------------------------------------------------------
C KLOE history
C ------------
C Adopted for KLOE by A.Andryakov in 1992
C Updated - 24.11.94 by A.Andryakov
C Updated - 16.02.95 by G.Capon (style changes)
C Updated - 17.03.95 by A.Andryakov (cleaned include file)
C effects all parts p.r. + t.f. + v.f
C Updated - 3.04.95 by G.Capon (ARGUS dynamical track candidates
C structure is replaced with YBOS bank structure)
C by A.Antonelli & A.Calcaterra (new DFSTOR to switch
C to YBOS banks)
C***********************************************************************
C
C=======================================================================
SUBROUTINE DFMAIN
C=======================================================================
$$IMPLICIT
C
C --- MAIN ROUTINE FOR TRACK FITTING
C-----------------------------------------------------------------------
C
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:OFERCO.INC'
C
$$INCLUDE 'C$INC:JOBSTA.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DTFS.INC'
$$INCLUDE 'K$ITRK:DCHD.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
$$INCLUDE 'K$ITRK:STUMPCOM.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
INTEGER SUPERSTATUS,DPVER,IHIT
INTEGER BLOCAT,BCOPYG,BNEXT,BDATA,BNUMB
INTEGER STATUS,STATUS_DPRS,STATUS_DTFS
INTEGER NPASS,I
INTEGER IND,INDAT,IND_DPRS,INDAT_DPRS,IND_DTFS,INDAT_DTFS
INTEGER INDR,INDDATR,DHRENCO,NDCHITS
INTEGER PTRCAN,LOOP
REAL ALONE,PRHITS,FTHITS,CHINOR,PCHI2
NPASS = 1
C Before starting the fit
C create the second (running) version of banks DTCE and DHRE
STATUS = BCOPYG(IW,'DTCE',1,'DTCE',2,IND,INDAT)
STATUS = BCOPYG(IW,'DHRE',1,'DHRE',2,IND,INDAT)
C Now locate the drift chamber header bank
STATUS = BLOCAT(IW,'DCHD',1,IND,INDAT)
IF (STATUS.EQ.YESUCC) THEN
IND = INDAT + IW(INDAT+DCHHDS)
IW(IND+DCHTRK) = 0
IW(IND+DCHTCN) = 0
ELSE
CALL ERLOGR('DFMAIN',ERWARN,0,STATUS,
+ 'DC Header Bank not found!')
RETURN
ENDIF
C WRITE(6,*)'Event nb :',NEV
C Loop over p.r. track candidates
20 continue
STATUS_DPRS = BLOCAT(IW,'DPRS',-1,IND_DPRS,INDAT_DPRS)
DO WHILE (IND_DPRS.GT.0)
IF (ldfjoi.GE.0) THEN
DPVER=IW(INDAT_DPRS+DPRVRN)
IF(NPASS.eq.1) THEN
if((DPVER.lt.10).and.DPVER.ne.NPASS) goto 100
ELSE
if(DPVER.ne.NPASS) goto 100
ENDIF
ENDIF
JDPRS = INDAT_DPRS + IW(INDAT_DPRS+DPRHDS)
*******************************************
IF (IW(JDPRS+DPRPOK) .EQ. 3) THEN
*******************************************
JDPRS1 = JDPRS + DPRCEN
STATUS = BNUMB(IW,IND_DPRS,ITR)
PRHITS = IW(JDPRS+DPRPOS)+IW(JDPRS+DPRNEG)
CALL HFILL(1,PRHITS,0.,1.)
CALL VZERO(NDF,MXHIT)
CALL DFFIT (NPASS)
IF (ldfjoi.GE.0) THEN
do IHIT=1,MXHIT
KDFS(ITR,IHIT)=KDF(IHIT)
DFSS(ITR,IHIT)=DFS(IHIT)
enddo
if(NPASS.eq.2) then
call DTFSKILL(ITR)
endif
ENDIF
C If the track fit succeeded, then the following bank must exist
STATUS_DTFS = BLOCAT(IW,'DTFS',ITR,IND_DTFS,INDAT_DTFS)
IND_DTFS = INDAT_DTFS + IW(INDAT_DTFS+DTFHDS)
IF (STATUS_DTFS.EQ.YESUCC) THEN
IW(IND+DCHTRK) = IW(IND+DCHTRK)+1
IF (IW(IND_DTFS+DTFCEN).EQ.1) THEN
IW(IND+DCHTCN) = IW(IND+DCHTCN)+1
ENDIF
FTHITS = NFITPOS+NFITNEG
CALL HFILL(4,FTHITS,0.,1.)
CHINOR=RW(IND_DTFS+DTFCHI)/(FTHITS-5)
CALL HFILL(9,CHINOR,0.,1.)
PCHI2 = RW(IND_DTFS+DTFPCH)
CALL HFILL(10,PCHI2,0.,1.)
ELSE
C If not, then the track fit failed. Histogram the number of hits for
C the candidate, and the last track fit error code
CALL HFILL(11,PRHITS,0.,1.)
CALL HFILL(12,FLOAT(NDFERR),0.,1.)
ENDIF
******************************************
ENDIF
******************************************
100 continue
STATUS_DPRS=BNEXT(IW,IND_DPRS,IND_DPRS)
STATUS_DTFS=BDATA(IW,IND_DPRS,INDAT_DPRS)
ENDDO
C At the end of the first step of track-fitting do some control histograms
CALL HFILL(5,FLOAT(IW(IND+DCHTRK)),0.,1.)
CALL HFILL(6,FLOAT(IW(IND+DCHTCN)),0.,1.)
STATUS = BLOCAT(IW,'DHRE',2,INDR,INDDATR)
IF (STATUS .NE. YESUCC) RETURN
INDR = INDDATR+IW(INDDATR+DHRHDS)
DHRENCO = IW(INDDATR+DHRNCO)
NDCHITS = IW(INDDATR+DHRNRO)
ALONE = 0.
PTRCAN = INDR+DHRTRN
DO LOOP=1,NDCHITS
IF (IW(PTRCAN).LE.0) ALONE = ALONE+1.
PTRCAN = PTRCAN+DHRENCO
ENDDO
CALL HFILL(7,ALONE,0.,1.)
CALL HFILL(8,ALONE/FLOAT(NDCHITS),0.,1.)
C
C Decides if a track is to be splitted in two
IF(KINK_RECO.and.(NPASS.eq.1)) THEN
CALL DFSPLIT
ENDIF
C
C
C Identify and join track stumps
IF (ldfjoi.GE.0) THEN
if(NPASS.eq.2) then
goto 900 !<COM6-11 for the moment only one join trk at time
endif
CALL DFJOIN(SUPERSTATUS)
if(SUPERSTATUS.eq.OFSUCC) then
NPASS = NPASS+1
goto 20
endif
ENDIF
900 RETURN
END
C
C*DK DFADD
C
C=======================================================================
C
SUBROUTINE DFADD (NADD,NERRF1)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Add a single hit to a partially fitted track
C
C Input parameters
C ================
INTEGER NADD ! index of the hit in the DHRE bank
INTEGER NERRF1 ! error multiplier from DFFIT via DFNEWW to DFITER
C
C Output parameters
C =================
C None
C
C Author:
C =======
C-
C- Originally from ARGUS collaboration
C- Modified June 1996 A. Calcaterra
C-
C----------------------------------------------------------------------
C
C Global Declarations:-
C =====================
C
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
C External Functions:-
C ====================
INTEGER BLOCAT
C Local Declarations:-
C ====================
REAL MAXTRP
PARAMETER (MAXTRP=10.)
LOGICAL CLOSE_E
INTEGER STATUS,INDR,INDDATR,DHRENCO
INTEGER N,NLST
INTEGER NC,I,INDWIRE,ITRSAV,RSAV
REAL D2D,STEP,MIN_STEP,EXTRAP,PHIZ0
REAL A1,A2,CHIMLT
REAL C,C1,C2,CMIN,AX,AY,AZ
C------- SPACE IN /CDFWORK/ EXHAUSTED. ADD NOTHING.
IF (NDFC .GE. MXHIT-1) THEN
CALL ERLOGR('DFADD',ERWARN,0,STATUS,
+ 'Working space exhausted.')
RETURN
ENDIF
C------- LOCATE RUNNING VERSION OF BANK DHRE
STATUS = BLOCAT(IW,'DHRE',2,INDR,INDDATR)
IF (STATUS.NE.YESUCC) THEN
CALL ERLOGR('DFADD',ERWARN,0,STATUS,
+ 'Bank DHRE not found.')
RETURN
ENDIF
INDR = INDDATR+IW(INDDATR+DHRHDS)
DHRENCO = IW(INDDATR+DHRNCO)
INDWIRE = INDR + (NADD-1)*DHRENCO
C------- FILL DF..(NC) ARRAYS. FIND FIRST EMPTY SLOT IN CDFWORK
C------- THERE *MUST* BE ONE BECAUSE OF THE TEST MADE ABOVE ON NDFC/MXHIT-1
NC = 1
DO WHILE (NDF(NC).GT.0 .OR. NC.EQ.NDFL)
NC = NC+1
ENDDO
KDF(NC) = NADD
DFD(NC) = ABS(RW(INDWIRE+DHRRAD)) !? BUT ...
MDF(NC) = 0
LDFL(NC) = IW(INDWIRE+DHRSLR)
DFE(NC) = RW(INDWIRE+DHRERR)**2
DFA(5,NC) = 0.
IDF(NC) = IW(INDWIRE+DHRWNR)
PHIZ0 = TIDDP(LDFL(NC)) * (IDF(NC)-1)
DFVX(NC) = -TIDST(LDFL(NC)) * SIN(PHIZ0)
DFVY(NC) = TIDST(LDFL(NC)) * COS(PHIZ0)
DFWX(NC) = TIDRA(LDFL(NC)) * COS(PHIZ0)
DFWY(NC) = TIDRA(LDFL(NC)) * SIN(PHIZ0)
C------- TRY TO FIND THE RIGHT PLACE IN THE PRESENT CHAIN OF HITS.
C------- "NLST" IS THE INDEX OF THE LAST HIT IN THE CHAIN BEFORE "STEP"
C------- BECOMES NEGATIVE. IF "NLST" TURNS OUT TO BE 0, THEN THE HIT TO
C------- BE ADDED IS THE FIRST ONE.
C------- BUT BEFORE TRUSTING DFPNT ONE MUST FIRST GET "CLOSE ENOUGH"
MIN_STEP = 1.E+14
STEP = 1.
NLST = 0
N = NDF1
CLOSE_E = .FALSE.
DO WHILE (N.GT.0 .AND. STEP.GT.0.)
D2D = SQRT((DFX(N)-DFWX(NC))**2+(DFY(N)-DFWY(NC))**2)
IF (D2D.LT.50.) THEN
CLOSE_E = .TRUE.
CALL DFPNT(N,NC,STEP)
IF (STEP.GT.0. .AND. STEP.LT.MIN_STEP) THEN
MIN_STEP = STEP
NLST = N
ENDIF
ENDIF
N = NDF(N)
ENDDO
IF (.NOT. CLOSE_E) THEN
CALL INERR (453)
RETURN
ENDIF
IF (NLST.GT.0) THEN
CALL DFNEXT(NLST,NC,MIN_STEP)
EXTRAP = SQRT((DFX(NC)-DFX(NLST))**2+(DFY(NC)-DFY(NLST))**2)
CHIMLT = 1.
ELSE
CALL DFNEXT(NDF1,NC,STEP)
EXTRAP = SQRT((DFX(NC)-DFX(NDF1))**2+(DFY(NC)-DFY(NDF1))**2)
CHIMLT = 5.
ENDIF
C------- REJECT EXTRAPOLATIONS TOO BIG
IF (EXTRAP .GT. MAXTRP) THEN
RETURN
ENDIF
C------- CHISQ CONTRIBUTION OF THE NEW WIRE; FIND SIGN OF DRIFT DISTANCE
A1 = TIDSC(LDFL(NC))
AX = DFVY(NC) * DFPZ(NC) - A1 * DFPY(NC)
AY = A1 * DFPX(NC) - DFVX(NC) * DFPZ(NC)
AZ = DFVX(NC) * DFPY(NC) - DFVY(NC) * DFPX(NC)
A2 = AX**2 + AY**2 + AZ**2
A1 = SQRT(A2)
C = ((DFWX(NC) - DFX(NC)) * AX +
+ (DFWY(NC) - DFY(NC)) * AY - DFZ(NC) * AZ ) / A1
C1 = (C + DFD(NC))**2 / DFE(NC)
C2 = (C - DFD(NC))**2 / DFE(NC)
C------- CUT ON CHI2: INCREASE IT BY FACTOR 5 IF EXTRAPOLATING BACKWARD
CMIN = AMIN1(C1,C2)
IF (CMIN .GT. CHIMLT*CDFCAD) RETURN
IF (CMIN.EQ.C2) THEN
DFD(NC) = -DFD(NC)
ELSE
ENDIF
C------- WE HAVE A CANDIDATE: IT MUST BE MARKED, WITH THE RIGHT SIGN, BOTH
C------- IN DHRE (BECAUSE DFITER--> DFBCOR WILL READ FROM THERE) AND IN DFD.
C------- IF THE FIT IS ATTEMPTED, AND FAILS, THE ORIGINAL VALUES WILL
C------- BE RECOVERED FROM ITRSAV AND RSAV.
ITRSAV = IW(INDWIRE+DHRTRN)
IW(INDWIRE+DHRTRN) = ITR
RSAV = RW(INDWIRE+DHRRAD)
RW(INDWIRE+DHRRAD) = DFD(NC)
NDFC = NDFC+1
IF (TIDST(LDFL(NC)) .LT. 0.) NFITNEG = NFITNEG+1
IF (TIDST(LDFL(NC)) .GT. 0.) NFITPOS = NFITPOS+1
IF (NLST .EQ. 0) THEN
C------- THIS HIT MUST BE ADDED BEFORE THE FIRST ONE
NDF(NC) = NDF1
DFS(NC) = -STEP
DFL(NC) = 0.
NDF1 = NC
N = NDF1
DO WHILE (NDF(N).GT.0)
DFL(NDF(N)) = DFL(N) + DFS(N)
N = NDF(N)
ENDDO
NDFMS(1) = NC
MDF(NC) = 0
DFCTGT = DFPZ(NC) / SQRT (DFPX(NC)**2 + DFPY(NC)**2)
DFSTHE = 1. / SQRT (1. + DFCTGT**2)
DFPHI = ATAN2 (DFPY(NC) , DFPX(NC))
ELSEIF (NLST .EQ. NDFL) THEN
C------- THIS HIT MUST BE ADDED AFTER THE LAST ONE
NDF(NDFL) = NC
NDF(NC) = 0
DFS(NDFL) = MIN_STEP
DFL(NC) = DFL(NDFL) + DFS(NDFL)
NDFL = NC
ELSE
C------- THIS HIT MUST BE ADDED IN THE MIDDLE OF THE CHAIN
NDF(NC) = NDF(NLST)
NDF(NLST) = NC
DFS(NC) = DFS(NLST) - MIN_STEP
DFS(NLST) = MIN_STEP
DFL(NC) = DFL(NLST) + DFS(NLST)
ENDIF
C------- CLEAR DERIVATIVES
DO I=1,NDFDIM
DFA(I,NC) = 0.
ENDDO
C------- ATTEMPT THE FIT NOW (IF THE POINT IS NOT REALLY TOO CLOSE: CDFDC1)
C------- IF IT WORKS, SAVE NEW HIT INCREASING NDFC
C------- IF IT DOES NOT, CALL DFKILL TO CLEAN UP VECTORS
C------- IF THE HIT *IS* REALLY TOO CLOSE DO NOT WASTE TIME: LEAVE ALL SUCH
C------- HITS IN THEIR VECTORS, FOR THE LAST ITERATION IN DFFIT.
IF (CMIN .GT. CDFDC1) THEN
C------- Let's save current situation
NDFERR = 0
IF(NLST.eq.0) NLST = NDFL
CALL DFSAVE(NLST,.FALSE.)
CALL DFITER (NERRF1*CDFDC1,1.)
IF (NDFERR.EQ.0) THEN
ELSE
C------- RECOVER ORIGINAL VALUES IN DHRE
IW(INDWIRE+DHRTRN) = ITRSAV
RW(INDWIRE+DHRRAD) = RSAV
C------- let's restore all vital track info !
CALL DFRESTORE(NLST,.FALSE.)
CALL DFKILL(NC)
IF (NDFERR.NE.0) RETURN
ENDIF
ELSE
ENDIF
RETURN
END
C
C*DK DFBCOR
C
C=======================================================================
C
SUBROUTINE DFBCOR
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : To correct the determination of drift distances
C applying fine space-time relations, depending
C on the phi impact angle of the track and on
C the particular local shape of the cells.
C Author: A. Calcaterra
C =======
C-
C- Created 17-1-96
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:BPOYBS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
INTEGER STATUS,BLOCAT,FINEST
INTEGER INDC,INDDATC,DTCENCO,INDR,INDDATR,DHRENCO
INTEGER N,LGEANFI,CELLDIM
REAL*4 BETA,PHITILDE
REAL*4 TZCOR,PREV_R,DRIFTR,ERRD
STATUS = BLOCAT(IW,'DTCE',2,INDC,INDDATC)
IF (STATUS.NE.YESUCC) RETURN
DTCENCO = IW(INDDATC+DTCNCO)
INDC = INDDATC+IW(INDDATC+DTCHDS)
STATUS = BLOCAT(IW,'DHRE',2,INDR,INDDATR)
IF (STATUS.NE.YESUCC) RETURN
DHRENCO = IW(INDDATR+DHRNCO)
INDR = INDDATR+IW(INDDATR+DHRHDS)
N = NDF1
DO WHILE (N.GT.0)
IF (LDFL(N) .GT. 0) THEN
C Translate layer n. from the ARGUS to the GEANFI convention
LGEANFI = 1+2*(59-LDFL(N))
CELLDIM = 3
IF (LDFL(N).GT.46) CELLDIM = 2
C Now calculate beta and phitilde for this particular hit,cell and track
CALL DFCBET(IDF(N),LGEANFI,
. DFX(N),DFY(N),DFZ(N),
. BETA,PHITILDE)
IF (BETA.LT.65.) BETA = 65.
IF (BETA.GT.130.) BETA = 130.
DFB(N) = BETA
DFF(N) = PHITILDE
C Finally recalculate drift distance from the Z-corrected drift time
TZCOR = ABS(RW(INDC+(KDF(N)-1)*DTCENCO+DTCTIM))
STATUS = FINEST(CELLDIM,TZCOR,
. BETA,PHITILDE*180./PI,
. DRIFTR,ERRD)
PREV_R = RW(INDR+(KDF(N)-1)*DHRENCO+DHRRAD)
RW(INDR+(KDF(N)-1)*DHRENCO+DHRRAD) = SIGN(DRIFTR,PREV_R)
RW(INDR+(KDF(N)-1)*DHRENCO+DHRERR) = ERRD
C Assume that big changes should not change signs.
C Small ones don't matter.
DFD(N) = SIGN(DRIFTR,PREV_R)
DFE(N) = ERRD**2
ENDIF
N = NDF(N)
ENDDO
RETURN
END
C
C*DK DFTCOR
C
C=======================================================================
C
SUBROUTINE DFTCOR
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : To correct the determination of drift times
C taking into account time propagation of the
C signal along the wire, which depends on Z,
C and the particle's time of flight.
C Author: A. Calcaterra
C =======
C-
C- Created 17-1-96
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:BPOYBS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
INTEGER BLOCAT
INTEGER STATUS,IND1,INDDAT1,IND2,INDDAT2,NHITC,DTCENCO
INTEGER N,LGEANFI
REAL*4 TMEAS,TPROP,TZCOR
LOGICAL ODDLAY
REAL*4 VINV,ALENGBY2
DATA VINV/0.035/,ALENGBY2/155./
STATUS = BLOCAT(IW,'DTCE',1,IND1,INDDAT1)
IF (STATUS.NE.YESUCC) RETURN
STATUS = BLOCAT(IW,'DTCE',2,IND2,INDDAT2)
IF (STATUS.NE.YESUCC) RETURN
NHITC = IW(INDDAT1+DTCNRO)
DTCENCO = IW(INDDAT1+DTCNCO)
IND1 = INDDAT1+IW(INDDAT1+DTCHDS)
IND2 = INDDAT2+IW(INDDAT2+DTCHDS)
N = NDF1
DO WHILE (N.GT.0)
IF (LDFL(N) .GT. 0) THEN
TMEAS = RW(IND1+(KDF(N)-1)*DTCENCO+DTCTIM)
LGEANFI = 59-LDFL(N)
ODDLAY = BTEST(LGEANFI,0)
IF (ODDLAY) THEN
TPROP = (ALENGBY2 - DFZ(N))*VINV
ELSE
TPROP = (ALENGBY2 + DFZ(N))*VINV
ENDIF
C Times from MC are signed like drift distances: so take the ABS()
TZCOR = ABS(TMEAS) - TPROP
TZCOR = TZCOR-(DFL(N)+RW(JDPRS+DPRCEN))/(DFBETA*30.)
C Assume that big changes should not change signs.
C Small ones don't matter.
RW(IND2+(KDF(N)-1)*DTCENCO+DTCTIM) = SIGN(TZCOR,TMEAS)
ENDIF
N = NDF(N)
ENDDO
RETURN
END
C
C*DK DFCBET
C=======================================================================
C
SUBROUTINE DFCBET(ISENSE,NPLANE,XHIT,YHIT,ZHIT,BETA,PHITILDE)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Finds the cell shape
C
C Input parameters
C ================
C INTEGER ISENSE ! index of the hit sense wire
C INTEGER NPLANE ! index of the hit sense wire plane
C REAL*4 XHIT ! X coordinate of the point of closest approach
C REAL*4 YHIT ! Y coordinate of the point of closest approach
C REAL*4 ZHIT ! Z coordinate of the point of closest approach
C
C Output parameters
C =================
C REAL*4 BETA ! angle that identifies the shape of the cell
C REAL*4 PHITILDE ! angle of the point of closest approach
C ! in the local reference system of the hit cell
C
C Author: P. de Simone
C =======
C-
C- Created 8-1-96
C-
C----------------------------------------------------------------------
INTEGER ISENSE, NPLANE
REAL*4 XHIT, YHIT, ZHIT
REAL*4 BETA, PHITILDE
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'RTDB$LIBRARY:GFGEOM.INC'
C
C Local declaration
C
INTEGER NF, IDFIELD, ISFIELD
REAL*4 PHI1, PHI2, PHI3, PHI4, PHI5
REAL*4 R1Z, R3Z, R4Z
REAL*4 X1, Y1, X2, Y2, X3, Y3, X4, Y4, X5, Y5
REAL*4 XP, YP, XF, YF
REAL*4 PHIZ, PHIF, PHID, PHIS
REAL*4 XD, YD, XS, YS
REAL*4 DD, DS
REAL*4 DX24, DY24, DX34, DY34, D24, D34
REAL*4 COSB
REAL*4 DX, DY, D, DX51, DY51, D51
REAL*4 COSP, RAGG
C
C--Definition of BETA and PHITILDE
C
C 4 3
C o-----o o
C /_/ beta
C /
C 2 o 1 + o
C
C
C o o o
C
C Calculate the PHITILDE angle, i.e. the "angle" of closest approach,
C which is an approximation of the direction of approach of the first
C cluster. The cells are divided in 36 PHITILDE wedges of +-5 degrees
C so that the appropriate STR can be chosen. For each of the 36 PHITILDE
C wedges there are 6 (all?)different STRs computed for 6 different BETAs
C
C o o o
C
C x (hit)
C \
C \--\ <-- phitilde
C o X__\____o --> wire 5
C
C
C
C
C o o o
C
C ----------------------------------
C find the phi angle of the sense wire on the plane z = 0
PHI1 = (ISENSE-1)*(2.*DTETWI(NPLANE))
C find the coordinates of the sense wire on the plane z = hit(3)
R1Z = RP0(NPLANE)*RP0(NPLANE) +
> ZHIT*ZHIT*TSTE(NPLANE)*TSTE(NPLANE)
R1Z = SQRT(R1Z)
X1 = RP0(NPLANE)*COS(PHI1) - ZHIT*SIN(PHI1)*TSTE(NPLANE)
Y1 = RP0(NPLANE)*SIN(PHI1) + ZHIT*COS(PHI1)*TSTE(NPLANE)
C find the phi angle of the sense wire on the plane z = hit(3)
PHIZ = ATAN2(Y1,X1)
IF ( PHIZ.LT.0. ) PHIZ = PHIZ + PI2
C Calculate Phitilde. Find coordinate of field wire 5 (see above)
C on the plane Z=Zhit
PHI5 = PHIZ - DTETWI(NPLANE)
IF (PHI5.LT.0.) PHI5 = PHI5 + PI2
X5 = R1Z*COS(PHI5)
Y5 = R1Z*SIN(PHI5)
C vector from the hit to the sense wire
DX = XHIT - X1
DY = YHIT - Y1
D = DX*DX + DY*DY
C vector from the field wire 5 (see drawing) to the sense wire
DX51 = X5 - X1
DY51 = Y5 - Y1
D51 = DX51*DX51 + DY51*DY51
IF ((D.NE.0.0).AND.(D51.NE.0.0)) THEN
COSP = ( DX*DX51 + DY*DY51 )/SQRT(D*D51)
ELSE
COSP = 0.0
ENDIF
IF (ABS(COSP).GT.0.99999) COSP = SIGN(0.99999,COSP)
RAGG = XHIT*XHIT + YHIT*YHIT
RAGG = SQRT(RAGG)
IF( R1Z.LT.RAGG ) THEN
PHITILDE = ACOS(COSP)
ELSE
PHITILDE = PI2 - ACOS(COSP)
ENDIF
C find the coordinates of field wire 2 on the plane z = hit(3)
PHI2 = PHIZ + DTETWI(NPLANE)
IF (PHI2.GT.PI2) PHI2 = PHI2 - PI2
X2 = R1Z*COS(PHI2)
Y2 = R1Z*SIN(PHI2)
C find the coordinates of field wire 3 on the plane z = hit(3)
NF = NPLANE + 1
R3Z = RP0(NF)*RP0(NF) + ZHIT*ZHIT*TSTE(NF)*TSTE(NF)
R3Z = SQRT(R3Z)
C 1. project the coord X1 and Y1 on the hyperboloid nf at z = hit(3)
XP = X1*(R3Z/R1Z)
YP = Y1*(R3Z/R1Z)
C 2. calculate the phi angle PHIF at the plane z = hit(3)
C of the first wire at R = RP0(NF) that has phi = 0 at z = 0
XF = RP0(NF)
YF = ZHIT*TSTE(NF)
PHIF = ATAN2(YF,XF)
C 3. index and phi angle of the field wires on the right
C and on the left of (XP,YP) on the plane z = hit(3)
IF( PHIF.LT.0 ) THEN
PHIF = - PHIF
IDFIELD = INT( (PHIZ+PHIF)/DTETWI(NF) )
PHID = DTETWI(NF)*IDFIELD - PHIF
ISFIELD = IDFIELD + 1
PHIS = DTETWI(NF)*ISFIELD - PHIF
ELSE
IDFIELD = INT( (PHIZ-PHIF)/DTETWI(NF) )
PHID = DTETWI(NF)*IDFIELD + PHIF
ISFIELD = IDFIELD + 1
PHIS = DTETWI(NF)*ISFIELD + PHIF
ENDIF
C 4. find the coordinate of the two found field wires and take the
C closest to the projection of the sense wire (XP,YP)
XD = R3Z*COS(PHID)
YD = R3Z*SIN(PHID)
XS = R3Z*COS(PHIS)
YS = R3Z*SIN(PHIS)
DD = SQRT((XP-XD)**2 + (YP-YD)**2)
DS = SQRT((XP-XS)**2 + (YP-YS)**2)
IF ( DD.GT.DS ) THEN
X3 = XS
Y3 = YS
PHI3 = PHIS
ELSE
X3 = XD
Y3 = YD
PHI3 = PHID
ENDIF
C find the coordinate of the field wire 4 on the plane z = hit(3)
PHI4 = PHI3 + DTETWI(NF)
R4Z = R3Z
X4 = R4Z*COS(PHI4)
Y4 = R4Z*SIN(PHI4)
C define the auxiliary vectors
10 DX24 = X2 - X4
DY24 = Y2 - Y4
DX34 = X3 - X4
DY34 = Y3 - Y4
D24 = DX24*DX24 + DY24*DY24
D34 = DX34*DX34 + DY34*DY34
C .. and calculate the beta angle
COSB = (DX24*DX34 + DY24*DY34)/SQRT(D24*D34)
BETA = ACOS(COSB)*180./PI
IF ( BETA.GE.65. ) RETURN
PHI4 = PHI3
PHI3 = PHI3 - DTETWI(NF)
X3 = R3Z*COS(PHI3)
Y3 = R3Z*SIN(PHI3)
X4 = R4Z*COS(PHI4)
Y4 = R4Z*SIN(PHI4)
GO TO 10
END
C
C*DK DFDRV
C
C=======================================================================
C
SUBROUTINE DFDRV
$$IMPLICIT
C------- CALCULATE DERIVATIVES, CHISQ, AND FILL MATRIX
$$INCLUDE 'C$INC:JOBSTA.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$include 'K$ITRK:CDFWORK.INC'
REAL X1,Y1,Z1
REAL PX1,PY1,PZ1
REAL PXN1,PYN1,XN,YN
REAL A1,A2,AX,AY,AZ,H1
REAL BTHETA,N1A,B1A,PHIMS
REAL THETAMS,TMS
INTEGER I,J,ND,N,ND1,M
INTEGER NPRHITS
X1 = DFX(NDF1)
Y1 = DFY(NDF1)
Z1 = DFZ(NDF1)
PX1 = DFPX(NDF1)
PY1 = DFPY(NDF1)
PZ1 = DFPZ(NDF1)
PXN1 = -PY1 / DFSTHE
PYN1 = PX1 / DFSTHE
C------- CLEAR XX AND YY
DFCHI2 = 0.0
DO I=1,NDFDIM+2
YY(I) = 0.D0
DO J=1,NDFDIM+2
XX(I,J) = 0.D0
ENDDO
ENDDO
ND = 6
N = NDF1
DO WHILE (N.NE.0)
IF (LDFL(N).LE.0) GO TO 400
XN = DFX(N)
YN = DFY(N)
C------- FORM VECTOR A = WIRE CROSS TANGENT (CF NIM(42))
A1 = TIDSC(LDFL(N))
AX = DFVY(N) * DFPZ(N) - A1 * DFPY(N)
AY = A1 * DFPX(N) - DFVX(N) * DFPZ(N)
AZ = DFVX(N) * DFPY(N) - DFVY(N) * DFPX(N)
A2 = AX**2 + AY**2 + AZ**2
A1 = SQRT(A2)
AX = AX / A1
AY = AY / A1
C------- CALCULATE DERIVATIVES (STEREO WIRES)
C------- CF NIM(50-54)
C --- DFA(1,N) = D DIST / D D
C --- DFA(2,N) = D DIST / D Z0
C --- DFA(3,N) = D DIST / D CURV
C --- DFA(4,N) = D DIST / D CTG(THETA)
C --- DFA(5,N) = D DIST / D (DP/DX) dummy !!!!
C --- DFA(6,N) = D DIST / D PHI
DFA(1,N) = PXN1 * AX + PYN1 * AY
DFA(2,N) = AZ / A1
H1 = (X1 - XN)*AX + (Y1 - YN)*AY + (Z1 - DFZ(N))*DFA(2,N)
DFA(3,N) = H1 / DFKURV
DFA(4,N) = DFL(N) * DFA(2,N) * DFSTHE
DFA(5,N) = 0.
DFA(6,N) = (XN-X1)*AY - (YN-Y1)*AX
ND1 = MAX0 (6,ND)
DO M=7,ND1,2
DFA(M,N) = (XN-DFX(NDFMS(M-5)))*AY
+ - (YN-DFY(NDFMS(M-5)))*AX
PHIMS=ATAN2(DFPY(NDFMS(M-4)),DFPX(NDFMS(M-4)))
N1A=SIN(PHIMS)*AX-COS(PHIMS)*AY
B1A=COS(PHIMS)*AX+SIN(PHIMS)*AY
THETAMS=ASIN(DFSTHE)
BTHETA=(DFL(N)-DFL(NDFMS(M-4)))
+ *DFH(N)*DFKURV*SIN(THETAMS)
DFA(M+1,N)=(DFL(N)-DFL(NDFMS(M-4)))
+ *COS(THETAMS)*( SIN(BTHETA)*N1A+COS(BTHETA)*B1A
+ - TAN(THETAMS)*DFA(2,N) )
ENDDO
C------- CHISQ CONTRIBUTION - CF NIM(21,43)
DFCH(N) = DFD(N) + (DFWX(N) - XN) * AX +
+ (DFWY(N) - YN) * AY - DFZ(N) * DFA(2,N)
DFCHI2 = DFCHI2 + DFCH(N)**2 / DFE(N)
C------- FILL MATRIX AND VECTOR - CF NIM(28)
DO I=1,ND
YY(I) = YY(I) + DFA(I,N) * DFCH(N) / DFE(N)
DO J=1,I
XX(I,J) = XX(I,J) + DFA(I,N) * DFA(J,N) / DFE(N)
ENDDO
ENDDO
400 IF (NDFMS(ND-4).EQ.N) ND = ND + 2
N = NDF(N)
ENDDO
C------- MULTIPLE SCATTERING CONSTRAINT
C------- DFEMS=1/SIGMA_SQUARE_MS - CF NIM(56)
XX(5,5) = 1.
IF (NDFDIM .GT. 5) THEN
DO I=7,NDFDIM
XX(I,I) = XX(I,I) + DFEMS(I-5)
YY(I) = YY(I) - (DFPMS(I-5)) * DFEMS(I-5)
ENDDO
ENDIF
RETURN
END
C
C*DK DFFIT
C
C=======================================================================
C
SUBROUTINE DFFIT (NPASS)
$$IMPLICIT
INTEGER NPASS
C======= FIT ONE TRACK
$$INCLUDE 'C$INC:JOBSTA.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$include 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
integer npass2
common /NP/npass2 !<COM9-11
C External Functions:-
C ====================
C Local Declarations:-
C ====================
INTEGER KINKPOS,KINK_FLAG
INTEGER NPRHITS
INTEGER NHITS1,NHITS2,NHITS3,NHITS4,
+ SHITS1,SHITS2
INTEGER NERRF1,NERRF2,IFLREV,N,I,M,I1
REAL RESID
REAL CH2A,CH2B
REAL ANG1,ANG2
COMMON /ANGLES/ ANG1(64),ANG2(64)
COMMON /NBHITS/ NPRHITS
COMMON /CH2AB/ CH2A(MXHIT),CH2B(MXHIT)
C======= SET FLAGS
NERRF1=1
NERRF2=1
KINK_FLAG = 0
c TOP OF LOOP OF NEXT TRIAL WITH LOOSER CHI2 CUTS
IFLREV = -1
NDFERR = 0
c NPRHITS : initial number of hits
NPRHITS = IW(JDPRS+DPRPOS) + IW(JDPRS+DPRNEG)
C LOOP BACK TO THIS POINT IN CASE OF WIRE-REVERSING
10 CONTINUE
C======= 1ST ITERATION
npass2= npass
DO M=1,64
ANG1(M) = 0.
ANG2(M) = 0.
ENDDO
CALL DFINIT (1,NPRHITS)
C Check that the track left a sufficient number of hits, with a sufficient
C number of stereo layers etc. etc.
IF (NFITPOS+NFITNEG .LE. NCUTSTEREO) THEN
NDFERR = 8
CALL INERR(416)
GOTO 888
ELSE
IF (NFITPOS.LT.3 .OR. NFITNEG.LT.3) THEN
NDFERR = 8
CALL INERR(416)
GOTO 888
ENDIF
ENDIF
C Check that the track crossed a sufficient number of different layers.
IF (NL_CR .LT. 2) THEN
NDFERR = 9
CALL INERR(407)
GOTO 888
ENDIF
CALL DFTRAC(.FALSE.)
FDFAST = .FALSE.
NDFINV = 0
NDFDIM = 6
100 CALL DFITER (NERRF1*CDFDC1,CDFM1L)
IF (NDFERR.NE.0) THEN
NERRF1=NERRF1+1
CALL INERR(408)
C Some error conditions: (Chi2 increasing) and (too many iterations)
C invite to try again the fit relaxing the first cut. The same could
C apply in principle for (pred chi2 wors), NDFERR=3
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) THEN
GOTO 888
ENDIF
IF (NERRF1.LT.15) GOTO 100
ENDIF
FDFAST = .TRUE.
C======= TRY FLIPPING SMALL DRIFT DISTANCES
CALL DFLIP(NERRF1)
C======= TRY TO FIND A KINK ALONG THE TRACK
C
C
IF(KINK_RECO) THEN
CALL DFKINK(ITR,KINKPOS,KINK_FLAG)
ENDIF
C
C dE/dx part (particles are assumed to be pions)
IF (ldfdx.GE.0) THEN
C+++++++++++++++++++++++++++++++++
C if(nprhits.gt.80) then ! for n<80 no dedx correction !<COM9-11
C+++++++++++++++++++++++++++++++++
150 CALL DFDEDX(1)
CALL DFITER(NERRF2*CDFDC2,CDFM1L)
IF (NDFERR.NE.0) THEN
NERRF2=NERRF2+1
CALL INERR(417)
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) then
C write(6,*) 'dedx - rej track, ndferr=',NDFERR
goto 888
ENDIF
IF (NERRF2.LT.15) THEN
NDFERR = 0
GOTO 150
ENDIF
ENDIF
C+++++++++++++++++++++++++++++++++
C endif !<COM9-11
C+++++++++++++++++++++++++++++++++
ENDIF
C======= MULTIPLE SCATTERING
IF (npass.eq.2.and.LDFMUS.GE.0) THEN
CALL DFMUSC
ENDIF
200 IF (NDFDIM.GT.6) THEN
CALL DFITER (NERRF1*CDFDC1,CDFM1L)
IF (NDFERR.NE.0) THEN
NERRF1=NERRF1+1
CALL INERR(409)
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) THEN
GOTO 888
ENDIF
IF (NERRF1.LT.15) GOTO 200
ENDIF
ENDIF
C======= REJECT WIRE(S) OR SWITCH SIGNS OF DRIFT TIME
C NHITS1 : number of hits after drift distance sign flipping
NHITS1 = NFITNEG+NFITPOS
IF (LDFREJ.GE.0) CALL DFREJ (NERRF1,60.,4.) !sv
IF (NDFERR.NE.0) THEN
CALL INERR(420)
GOTO 888
ENDIF
C======= FIND NEW WIRES
c NHITS2 : number of hits after hit rejection (DFREJ)
c SHITS1 : number of rejected hits after 1st call to DFREJ
c NHITS3 : number of hits after hit addition
NHITS2 = NFITNEG+NFITPOS
SHITS1 = NHITS1-NHITS2
IF (LDFADD.GE.0) THEN
CALL DFNEWW(NERRF1)
NHITS3 = NFITNEG+NFITPOS
CALL HFILL(2,FLOAT(NHITS3-NHITS2),0.,1.)
ENDIF
C======= REJECT WIRES, AGAIN
IF (LDFREJ.GE.0) CALL DFREJ (NERRF1,60.,4.) !sv
IF (NDFERR.NE.0) THEN
CALL INERR(421)
GOTO 888
ENDIF
c NHITS4 : number of hits after hit rejection (DFREJ)
c SHITS2 : number of rejected hits after 2nd call to DFREJ
NHITS4 = NFITNEG+NFITPOS
SHITS2 = NHITS3-NHITS4
C======= MULTIPLE SCATTERING
IF (npass.eq.1.and.LDFMUS.GE.0) CALL DFMUSC
210 IF (NDFDIM.GT.6) THEN
CALL DFITER (NERRF1*CDFDC1,CDFM1L)
IF (NDFERR.NE.0) THEN
NERRF1=NERRF1+1
CALL INERR(409)
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) THEN
GOTO 888
ENDIF
IF (NERRF1.LT.15) GOTO 210
ENDIF
ENDIF
C======= LAST ITERATION ; GET ERROR MATRIX
NDFINV = 6
300 CALL DFITER (NERRF2*CDFDC2,CDFM1L)
IF (NDFERR.NE.0) THEN
NERRF2=NERRF2+1
CALL INERR(410)
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) THEN
GOTO 888
ENDIF
IF (NERRF2.LT.15) GOTO 300
ENDIF
C======= IF EVERYTHING WORKED, PLOT RESIDUALS AND STORE RESULTS
C nhi=0 !<COMM6-11
IF (NDFERR .EQ. 0) THEN
IF (ABS(DFCTGT) .LT. 1000.) THEN
N = NDF1
DO WHILE (N.GT.0)
IF (LDFL(N).GT.0) THEN
RESID = - DFCH(N)
DO I=1,NDFDIM
RESID = RESID + XX(I,NDFDIM+2) * DFA(I,N)
ENDDO
CALL HFILL(450,RESID/SQRT(DFE(N)),0.,1.)
ENDIF
N = NDF(N)
ENDDO
CALL DFSTOR(NPASS)
ELSE ! THE FIT CONVERGED ALL RIGHT, BUT TO WHAT, GOD KNOWS.
CALL INERR(492)
ENDIF
ENDIF
IF (NDFERR .EQ. 0) THEN
CALL HFILL(3,FLOAT(SHITS1+SHITS2),0.,1.)
ENDIF
888 Continue
C
C===== PREPARES NEW DPRS BANKS IF A KINK HAS BEEN FOUND
C
IF(KINK_RECO) THEN
IF(KINK_FLAG.EQ.-1) CALL DFTWOPR(KINKPOS,ITR)
ENDIF
C
C
RETURN
END
C
C*DK DFINIT
C
C=======================================================================
C
SUBROUTINE DFINIT (NFIRST,NLAST)
$$IMPLICIT
C------ INITIALIZE TRACK FIT FOR ONE TRACK
C Call Parameters:-
C ================
INTEGER NFIRST,NLAST
C Global Declarations:-
C =====================
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
C External Functions:-
C ====================
INTEGER BLOCAT
C Local Declarations:-
C ====================
INTEGER STATUS,INDR,INDDATR,DHRENCO,IND
INTEGER NDIFF(MXHIT),TEST
LOGICAL DIFFER
REAL X(3),H(3)
REAL PHIZ0, PGEV, DIVSTEP, DIST2D
INTEGER NC,IWIRE
INTEGER NADD,NADDST,L
LOGICAL ITERATE
CALL VZERO(DFX,MXHIT)
CALL VZERO(DFY,MXHIT)
CALL VZERO(DFZ,MXHIT)
CALL VZERO(DFWX,MXHIT)
CALL VZERO(DFWY,MXHIT)
CALL VZERO(DFPX,MXHIT)
CALL VZERO(DFPY,MXHIT)
CALL VZERO(DFPZ,MXHIT)
CALL VZERO(DFS,MXHIT)
CALL VZERO(DFL,MXHIT)
CALL VZERO(LDFL,MXHIT)
* dhre bank (INDR,INDDATR)
STATUS = BLOCAT(IW,'DHRE',2,INDR,INDDATR)
INDR = INDDATR+IW(INDDATR+DHRHDS)
DHRENCO = IW(INDDATR+DHRNCO)
C------ GET HELIX STARTING PARAMETERS
C------ DFKURV = 1/(RHO*HFIELD) - CF NIM(37)
NDF1 = 1
DFKURV = RW(JDPRS+DPRCUR)
DFCTGT = RW(JDPRS+DPRCTT)
DFSTHE = 1. / SQRT (1. + DFCTGT**2)
PGEV = 1./ABS(DFKURV*DFSTHE)
C------ BETA IN THE PION MASS HYPOTHESIS
DFBETA = PGEV/SQRT(PGEV**2+PION_MASS**2)
DFPHI = RW(JDPRS+DPRPHI)
NDFMS(1) = NDF1
NDFMS(2) = 0
DFDPDX = 0.
DFCMSO = 0.
DFCMS2 = 0.
DFX(1) = RW(JDPRS+DPRXCO)
DFY(1) = RW(JDPRS+DPRYCO)
DFZ(1) = RW(JDPRS+DPRZCO)
DFPY(1) = SIN (DFPHI) * DFSTHE
DFPX(1) = COS (DFPHI) * DFSTHE
DFPZ(1) = DFCTGT * DFSTHE
DFSTLC = SQRT (CDFACU/ABS(DFKURV)) / DFSTHE
C DFSTCL = SQRT (2 * RHO * ACCURACY)/SIN(THETA)
C-------- IS MAX STEP LENGTH FOR WHICH STRAIGHT LINE EXTRAPOLATION
C IS ACCURATE WITHIN CDFACU = 200 MICRON (JAN 95)
C SEE KAPITZA (65)
DFH(1) = TIFLAV
X(1) = DFX(1)
X(2) = DFY(1)
X(3) = DFZ(1)
CALL MGFLD(X,H)
CALL VSCALE(H,CONVERS,H,3)
DFHX(1) = H(1) / DFH(1)
DFHY(1) = H(2) / DFH(1)
DFHZ(1) = H(3) / DFH(1)
NC = 0
NFITPOS = 0
NFITNEG = 0
NL_CR = 0
L = NFIRST
ITERATE = .TRUE.
DO WHILE (L.LE.NLAST .AND. ITERATE)
NC = NC + 1
C-------- DFINIT MAY BE CALLED AGAIN, AFTER A HIT-REJECTION STEP.
C-------- FOR NOW, WE DECIDE TO RE-ACCEPT THE HIT ...
IWIRE = IABS(IW(JDPRS1+DPRLWI*(L-1)+DPRIWI))
NDF(NC) = NC + 1
DIST2D = RW(JDPRS1+DPRLWI*(L-1)+DPRWTR)
C--------------------- 2D distances greater than 20 cm are not plausible
C--------------------- and they may make DFTRAC to loop forever
DFS(NC) = AMIN1(DIST2D,20.) / DFSTHE
DFL(NC) = 0.
C--------------------- FILL DF..(NC) ARRAYS
IND = INDR + (IWIRE-1)*DHRENCO
KDF(NC) = IWIRE
IDF(NC) = IW(IND+DHRWNR)
MDF(NC) = 0
DFD(NC) = RW(IND+DHRRAD)
LDFL(NC) = IW(IND+DHRSLR)
DFE(NC) = RW(IND+DHRERR)**2
DFA(5,NC) = 0.
PHIZ0 = TIDDP(LDFL(NC)) * (IW(IND+DHRWNR)-1)
DFVX(NC) = -TIDST(LDFL(NC)) * SIN(PHIZ0)
DFVY(NC) = TIDST(LDFL(NC)) * COS(PHIZ0)
DFWX(NC) = TIDRA(LDFL(NC)) * COS(PHIZ0)
DFWY(NC) = TIDRA(LDFL(NC)) * SIN(PHIZ0)
IF (TIDST(LDFL(NC)) .LT. 0.) NFITNEG = NFITNEG + 1
IF (TIDST(LDFL(NC)) .GT. 0.) NFITPOS = NFITPOS + 1
DIFFER = .TRUE.
DO TEST=1,NL_CR
IF (LDFL(NC) .EQ. NDIFF(TEST)) DIFFER = .FALSE.
ENDDO
IF (DIFFER) THEN
NL_CR = NL_CR + 1
NDIFF(NL_CR) = LDFL(NC)
ENDIF
C-------- INSERT ADDITIONAL STEPS IF STEP LENGTH TOO LARGE
IF (DFS(NC) .GT. CDFSTL) THEN
NADDST = DFS(NC) / CDFSTL
DIVSTEP = DFS(NC) / (NADDST+1)
DFS(NC) = DIVSTEP
DO NADD=NC+1,MIN0(NC+NADDST,MXHIT)
DFS(NADD) = DIVSTEP
NDF(NADD) = NADD + 1
KDF(NADD) = 0
MDF(NADD) = 0
DFCH(NADD) = 0.
LDFL(NADD) = 0
ENDDO
NC = NC+NADDST
ENDIF
C Here the condition for iterating the loop is calculated
ITERATE = NC+1 .LE. MXHIT-1
L = L + 1
ENDDO
C Now flag premature termination of the loop, just to know
IF (L.LT.NLAST) THEN
CALL INERR (405)
ENDIF
DFS(NC) = 0.
NDF(NC) = 0
NDFL = NC
NDFC = NC
RETURN
END
C
C*DK DFITER
C
C=======================================================================
C
SUBROUTINE DFITER (CUTCHI,CMAX1L)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : This routine minimizes the track chi-square
C- doing a maximum of "LDFITE" iterations.
C- The outputs are : a chi-2 value in DFCHIN,
C- the fitted set of 5 helyx parameters, and
C- the full series of points of closest approach
C- to the wires (indirectly, via the cumulative
C- flight distances DFL for example)
C- (All of this meaningful only if NDFERR = 0 ....)
C-
C- Created
C- Updated
C-
C----------------------------------------------------------------------
$$INCLUDE 'C$INC:JOBSTA.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
INTEGER NDFITE,M,N,I,NPRHITS
REAL CUTCHI,CMAX1L
REAL X,CH2,CH2A,CH2B,DELTA
LOGICAL TOLERATE,FAST_SWIM
COMMON /NBHITS/ NPRHITS
COMMON /CH2AB/ CH2A(MXHIT),CH2B(MXHIT)
NDFERR = 0
DFCHIO = 1.E38
DFNDGF = NFITNEG + NFITPOS - 5.
TOLERATE = .TRUE.
DO NDFITE=1,LDFITE
C------- TRACE PARTICLE: FIND POINTS OF CLOSEST APPROACH TO THE WIRES
FAST_SWIM = NDFITE .GT. 1
CALL DFTRAC(FAST_SWIM)
IF (NDFERR.NE.0) THEN
C------- THE ONLY ERROR IN DFTRAC IS SIGNALLED IN DFTRAC ITSELF
RETURN
ENDIF
C------- USE POINTS OF C.A. TO IMPROVE CALCULATION OF DRIFT DISTANCES.
C------- DO 1) Z-DEPENDENT CORRECTION OF TIME FOR PROPAGATION ALONG WIRE
C------- 2) CHOOSE PHI AND BETA-DEPENDENT SPACE-TIME RELATION
IF (NDFITE.EQ.1 .AND. TS_RAWREL) CALL DFTCOR
IF (NDFITE.EQ.1 .AND. TS_FINEREL)CALL DFBCOR
C------- CALCULATE DERIVATIVES AND DFCHI2, DIAGONAL CHI-SQUARE
C------- WITH CONTRIBUTIONS ONLY FROM MEASUREMENTS.
C------- DFCMS2 IS CALCULATED BELOW, HERE IN DFITER.
CALL DFDRV
X = DFCHI2 + DFCMS2 - DFCHIO - DFCMSO
IF (X.GT.0.1 .AND. NDFITE.GT.1) THEN
C-------- CHISQ INCREASED SINCE LAST ITERATION.
C-------- DURING A GIVEN FIT THIS IS ACTUALLY ALLOWED, BUT ONLY ONCE.
IF (TOLERATE) THEN
TOLERATE = .FALSE.
ELSE
CALL INERR (412)
NDFERR = 2
RETURN
ENDIF
ENDIF
DFCHIO = DFCHI2
DFCMSO = DFCMS2
C-------- INVERT MATRIX; IF SINGULAR MATRIX, SIGNAL ERROR AND RETURNS
CALL LEQU64 (XX,XI,YY,NDFDIM,NDFINV,NDFERR)
IF (NDFERR.NE.0) THEN
CALL INERR (411)
RETURN
ENDIF
C-------- CALCULATE THE NEW TRACK PARAMETERS
DFX(NDF1) = DFX(NDF1) - DFPY(NDF1) * XX(1,NDFDIM+2) / DFSTHE
DFY(NDF1) = DFY(NDF1) + DFPX(NDF1) * XX(1,NDFDIM+2) / DFSTHE
DFZ(NDF1) = DFZ(NDF1) + XX(2,NDFDIM+2)
DFKURV = DFKURV + XX(3,NDFDIM+2)
IF (ABS(DFKURV).LE.1.E-10) THEN
NDFERR=10
RETURN
ENDIF
DFCTGT = DFCTGT + XX(4,NDFDIM+2)
IF (ABS(DFCTGT).GT.1.E+08) THEN
NDFERR=10
RETURN
ENDIF
DFSTHE = 1. / SQRT (1. + DFCTGT**2)
DFPMS(1) = DFPMS(1) + XX(6,NDFDIM+2)
DO M=7,NDFDIM,2
DFPMS(M-5) = DFPMS(M-5) + XX(M,NDFDIM+2)
DFPMS(M-4) = DFPMS(M-4) + XX(M+1,NDFDIM+2)
ENDDO
DFPY(NDF1) = SIN (DFPHI) * DFSTHE
DFPX(NDF1) = COS (DFPHI) * DFSTHE
DFPZ(NDF1) = DFCTGT * DFSTHE
IF (ABS (DFZ(NDF1)).GT.CDFZMA . OR .
+ ABS (DFZ(NDFL)).GT.CDFZMA) THEN
C-------- TRACK OUTSIDE VISIBLE VOLUME
NDFERR = 5
CALL INERR (431)
RETURN
ENDIF
C-------- NEW CHISQ - KAPITZA(24)
C-------- NDFMCH IS THE WIRE WITH THE BIGGEST CHISQ (=DFMCH)
DFCHIN = 0.
DFMCH = 0.
CALL VZERO(CH2A,MXHIT)
CALL VZERO(CH2B,MXHIT)
N = NDF1
C IT_LOOP=0 !<COM 6-11
DO WHILE (N.GT.0)
IF (LDFL(N).GT.0) THEN
CH2A(N) = - DFCH(N)
DO I=1,NDFDIM
CH2B(N) = CH2B(N) + XX(I,NDFDIM+2) * DFA(I,N)
ENDDO
CH2 = CH2A(N) + CH2B(N)
CH2 = CH2**2 / DFE(N)
DFCHIN = DFCHIN + CH2
IF (N.EQ.NDF1 . OR . N.EQ.NDFL) CH2 = CH2 * CMAX1L
C if(CH2.GE.25) IT_LOOP = IT_LOOP + 1 !<COM 6-11
IF (CH2.GE.DFMCH) THEN
DFMCH = CH2
NDFMCH = N
ENDIF
ENDIF
N = NDF(N)
ENDDO
C------- MULT. SCATT. CONTRIBUTION TO CHISQ - CF NIM(55)
DFCMS2 = 0.
DO M=7,NDFDIM
DFCMS2 = DFCMS2 + DFPMS(M-5)**2 * DFEMS(M-5)
ENDDO
DELTA=0.05*NPRHITS
IF (DFCHIN+DFCMS2 .GT. DFCHI2+DFCMSO+DELTA) THEN
C------- PREDICTED CHISQ WORSE THAN CHISQ OF THIS ITERATION
CALL INERR (413)
NDFERR = 3
RETURN
ENDIF
C CHECK IF THE ITERATION PROCEDURE HAS CONVERGED
IF (DFCHI2+DFCMSO-DFCHIN-DFCMS2 .LE. CUTCHI) THEN
RETURN
ENDIF
ENDDO
C------- TOO MANY ITERATIONS
CALL INERR (414)
NDFERR = 6
RETURN
END
C
C*DK DFJOIN !< COM6-11 (METTERE DFJOIN)
C=======================================================================
C
SUBROUTINE DFJOIN(SUPERSTATUS)
C
C----------------------------------------------------------------------
C-
C- Purpose and Methods : This routine scans all track pairs to identify
C- cases in which a track has been split in 2
C- (or more) stumps by the old ARGUS pattern
C- routines.
C-
C- Author:
C =======
C G.Venanzoni 29-9-97
C
C----------------------------------------------------------------------
C
C------ Global Declarations -------------------------------------------
C
$$implicit
C
C------ Include Files -------------------------------------------------
C
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:JOBSTA.INC'
$$include 'K$INC:OFERCO.INC'
$$INCLUDE 'K$ITRK:DTFS.INC'
$$INCLUDE 'K$ITRK:STUMPCOM.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
C
C------ External Declarations ------------------------------------------
C
C------ Local Declarations ---------------------------------------------
integer SUPERSTATUS,STATUS
integer COMSTATUS,IS,NDUPL,NSTRNG,ICOP
integer TRK1,TRK2,VCOP(5),FCOP,NCOPCON
integer V_FCOP(5),V_NCOP(5),I
integer A_NCOP(5,5),ISTR,JSTR
integer NSLD,NSLD_ICH,S_NCOP,S_VCOP(5)
logical SFL(5),DUPFL,LFLG
real pi
DATA pi /3.1415/
C------ Start of code -------------------------------------------------
CStart with success
SUPERSTATUS=OFSUCC
call DFJSTUMP(COMSTATUS) !Fill STUMPCOM Common Block
if (COMSTATUS.ne.OFSUCC) NCOP =0
if(NCOP.eq.0) goto 900
do IS=1,5
SFL(IS)=.true.
enddo
NDUPL = 0
NSTRNG = 0
do 100 ICOP =1, NCOP
TRK1 = FIRSTUMP(ICOP)
TRK2 = SECSTUMP(ICOP)
if (FHIT(ICOP).eq.1) then !icop is broken on FIRST HIT
DUPFL = FHIT(ICOP).eq.1.and.LASTHIT(TRK1).eq.0
& .and.LASTHIT(TRK2).eq.0 !flag on duplet
if(DUPFL) then ! icop is a DUBLET
NDUPL = NDUPL + 1
call DFJMRG(TRK1,TRK2,STATUS)
goto 100
else ! icop is a point of a TRIPLET or more
NSTRNG = NSTRNG + 1
LFLG = .true.
call DFJTRP(LFLG,TRK1,TRK2,VCOP,FCOP,NCOPCON)
V_FCOP(NSTRNG)=FCOP
V_NCOP(NSTRNG)=NCOPCON
do I=1,NCOPCON
A_NCOP(NSTRNG,I) =VCOP(I)
enddo
goto 100
endif
else ! icop is broken on LAST HIT
DUPFL = LHIT(ICOP).eq.1.and.FIRSTHIT(TRK1).eq.0
& .and.FIRSTHIT(TRK2).eq.0 !flag on duplet
if(DUPFL) then
NDUPL = NDUPL + 1
call DFJMRG(TRK1,TRK2,STATUS)
goto 100
else ! icop is a point of a TRIPLET or more
NSTRNG = NSTRNG + 1
LFLG = .false.
call DFJTRP(LFLG,TRK1,TRK2,VCOP,FCOP,NCOPCON)
V_FCOP(NSTRNG)=FCOP
V_NCOP(NSTRNG)=NCOPCON
do I=1,NCOPCON
A_NCOP(NSTRNG,I) =VCOP(I)
enddo
goto 100
endif
endif
100 continue
CLook how many different configurations:
if(NSTRNG.le.1) goto 200
do 170 ISTR = 1,NSTRNG-1
if(.not.SFL(ISTR)) goto 170
do JSTR = ISTR+1,NSTRNG
if(V_FCOP(ISTR).eq.V_FCOP(JSTR)) SFL(JSTR) = .false.
enddo
170 continue
C ncop = nsing + ntrip..
200 continue
NSLD = 0
do ISTR =1,NSTRNG
if (SFL(ISTR)) then
NSLD_ICH = NCOPCON
NSLD = NSLD + NSLD_ICH
endif
enddo
if (NCOP.ne.NDUPL+NSLD) then
print*, 'ERROR in calculation NCOP !!!!'
endif
do ISTR =1,NSTRNG
if (SFL(ISTR)) then
S_NCOP=V_NCOP(ISTR)
do I=1,S_NCOP
S_VCOP(I)=A_NCOP(ISTR,I)
enddo
call DFJSTR(ISTR,S_NCOP,S_VCOP)
endif
enddo
goto 990
900 SUPERSTATUS=OFFAIL
990 continue
return
end
C=======================================================================
C
SUBROUTINE DFJTRP(LFLG,TRK1,TRK2,VCOP,FCOP,NCOPCON)
C
C=======================================================================
C
C Description:
C ============
C
C
C Author:
C =======
C G.Venanzoni 22-9-97
C------------------------------------------------------------------------
C----------------------------------------------------------------------
C
C------ Global Declarations -------------------------------------------
C
$$IMPLICIT
C
C
C------ Include Files -------------------------------------------------
C
$$INCLUDE 'K$INC:ERLEVL.INC'
$$include 'k$inc:oferco.inc'
$$INCLUDE 'K$ITRK:STUMPCOM.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
C
C------ External Declaration --------------------------------------------
C
integer GETCOP
C
C------ Local Declarations -------------------------------------------
C
logical LFLG,LFLG_FIRST
integer TRK,TRK1,TRK2,NTRK,IT,TRKNUM(6)
integer COPARRAY(5),HIT,TRK_OLD,COPMIN
integer ICOP,NCOPCON,FCOP,VCOP(5)
C
C----- Start of Code --------------------------------------------------
C
LFLG_FIRST = LFLG
TRK = TRK1
NTRK= 1
IT = 1
TRKNUM(1) = TRK1
COPARRAY(1)= GETCOP(TRK1,TRK2)
if (COPARRAY(1).eq.-1) then
call ERLOGR('DFJTRP',ERWARN,0,0,'Cannot find COP')
goto 900
endif
20 continue
if (LFLG) then
HIT = LASTHIT(TRK)
else
HIT = FIRSTHIT(TRK)
endif
do while (HIT.eq.1)
IT = IT + 1 ! # of trks connected with FIRSTUMP (included the first)
if(LFLG) then
TRK_OLD = TRK
TRK = COMP_L(TRK)
TRKNUM(IT) = TRK
COPARRAY(IT) = GETCOP(TRK_OLD,TRK)
if(COPARRAY(IT).eq.-1) then
call ERLOGR('DFJTRP',ERWARN,0,0,'Cannot find COP')
goto 900
endif
HIT = FIRSTHIT(TRK)
LFLG = .false.
else
TRK_OLD = TRK
TRK = COMP_F(TRK)
TRKNUM(IT) = TRK
COPARRAY(IT) = GETCOP(TRK_OLD,TRK)
if(COPARRAY(IT).eq.-1) then
call ERLOGR('DFJTRP',ERWARN,0,0,'Cannot find COP')
goto 900
endif
HIT = LASTHIT(TRK)
LFLG = .true.
endif
enddo
if(NTRK.eq.2) then
goto 100
else
NTRK=2
TRK = TRK2
LFLG = LFLG_FIRST
goto 20
endif
100 continue
TRKNUM(IT+1) = TRK2
COPMIN=10
do ICOP=1,IT
if(COPARRAY(ICOP).lt.COPMIN) then
COPMIN=COPARRAY(ICOP)
endif
enddo
FCOP=COPMIN !NUMBER OF FIRST COP
NCOPCON=IT !NUMBER OF COPS CONNECTED
do ICOP=1,NCOPCON
VCOP(ICOP)=COPARRAY(ICOP)
enddo
call VZERO(COPARRAY,NCOPCON)
900 return
end
C================================================================
C
Integer Function GETCOP(TRK1,TRK2)
C
C================================================================
C
C Description:
C ============
CIT retrives the Number of COP connecting TRK1 AND TRK2
C
C Author:
C =======
C G.Venanzoni 22-9
C
C---------------------------------------------------------------------
C------ Global Declarations -------------------------------------------
C
$$IMPLICIT
C
C
C------ Include Files -------------------------------------------------
C
$$INCLUDE 'K$INC:ERLEVL.INC'
$$include 'k$inc:oferco.inc'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:STUMPCOM.INC'
C
C------ External Declaration --------------------------------------------
c
C------ Local Declarations -------------------------------------------
C
integer IC,TRK1,TRK2
C----- Start of Code --------------------------------------------------
C
GETCOP=-1
do IC=1,NCOP
if((FIRSTUMP(IC).eq.TRK1.and.SECSTUMP(IC).eq.TRK2).or.
& (FIRSTUMP(IC).eq.TRK2.and.SECSTUMP(IC).eq.TRK1)) then
GETCOP=IC
goto 900
endif
enddo
900 continue
return
end
C===================================================================
C
SUBROUTINE DFJSTR(ISTRING,S_NCOP,S_VCOP)
C
C===================================================================
C
C Description:
C ============
CIT breaks String #ISTRING according to the Quality Factor of each COP
C
C Author:
C =======
C G.Venanzoni 23-9
C
C------ Global Declarations -------------------------------------------
C
$$IMPLICIT
C
C
C------ Include Files -------------------------------------------------
C
$$INCLUDE 'K$INC:ERLEVL.INC'
$$include 'k$inc:oferco.inc'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:STUMPCOM.INC'
C
C------ External Declaration --------------------------------------------
c
C------ Local Declarations -------------------------------------------
C
integer ICOP,S_NCOP,S_VCOP(5),SNCOP,ISTRING
integer BESTICOP,TRK1,TRK2,TRK1_OLD,TRK2_OLD
real QF(5),QFGI,QF_OLD
integer IDCOP_OLD,STATUS
C----- Start of Code --------------------------------------------------
C
do ICOP=1,S_NCOP
QF(ICOP)=QF0(S_VCOP(ICOP))+QF1(S_VCOP(ICOP))+QF2(S_VCOP(ICOP))+
+ QF3(S_VCOP(ICOP))+QF4(S_VCOP(ICOP)) !FOR THE MOMENT THE SIMPLEST COMB.
enddo
Cget the best QF and start to connect the corrispective ICOP
SNCOP = S_NCOP
do while(SNCOP.gt.0)
QFGI =0.
do 20 ICOP =1,SNCOP
if(QF(ICOP).gt.QFGI) then
BESTICOP = ICOP
QFGI = QF(ICOP)
endif
20 continue
TRK1=FIRSTUMP(S_VCOP(BESTICOP))
TRK2=SECSTUMP(S_VCOP(BESTICOP))
TRK1_OLD = 0
TRK2_OLD = 0
if (TRK1.eq.TRK1_OLD.or.TRK1.eq.TRK2_OLD.or.
& TRK2.eq.TRK1_OLD.or.TRK1.eq.TRK2_OLD) goto 100
call DFJMRG(TRK1,TRK2,STATUS)
TRK1_OLD = TRK1
TRK2_OLD = TRK2
100 continue
QF_OLD = QF(NCOP) !Store Q.F.
IDCOP_OLD = S_VCOP(NCOP) !Store Identity of COP
QF(NCOP) = QF(BESTICOP)
S_VCOP(NCOP) = S_VCOP(BESTICOP)
QF(BESTICOP)=QF_OLD
S_VCOP(BESTICOP)= IDCOP_OLD
SNCOP = SNCOP-1
enddo
900 continue
RETURN
END
C=======================================================================
C
SUBROUTINE DFJMRG(TRK1,TRK2,STATUS)
C=======================================================================
C
C DESCRIPTION:
C ===========
CIt makes the merging of the two DPRS banks, creating the # BIGGGEST+1
Cnew DPRS Bank
C
C Author:
C =======
C G.Venanzoni 23-9-97
C
C------ Global declaration ---------------------------------------------
c
$$implicit
c
c------ Include Files --------------------------------------------------
c
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$include 'K$INC:OFERCO.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$include 'K$ITRK:DHRE.INC'
$$include 'K$ITRK:DPRS.INC'
$$include 'K$ITRK:DTFS.INC'
$$include 'K$ITRK:STUMPCOM.INC'
c
c------ External Declaration -------------------------------------------
c
integer BIGEST,BTYMKG,BLOCAT,GETCOP
real GETLENG
c------ Local Declaration ----------------------------------------------
c
integer JOIN_NUM,ISTAT,STATUS
integer DPRSBNK,BNKSIZ,IND_DP,INDDAT_DP
integer ISTAT_DTF1,ISTAT_DTF2,IND_DTF1,IND_DTF2
integer TRK1,TRK2,INDDAT_DTF1,INDDAT_DTF2
integer POINT_TR1,POINT_TR2,POINT_FIRSTRK,POINT_SECTRK
integer FTRK,STRK,NHITPOS_FTRK,NHITNEG_FTRK,NHIT_FTRK
integer NHITPOS_STRK,NHITNEG_STRK,NHIT_STRK,NHITD
integer I,IWIRE,W_PREV,ICOP
integer IND_DHR,INDDAT_DHR,DHRECO
real XCOOR_HIT1,YCOOR_HIT1,XCOOR_HIT2,YCOOR_HIT2
real DISTXY1,DISTXY2,F1,F2,COTFIRST,SINFIRST
real COSFIRST,COTSEC,SINSEC,COSSEC,COSTHETA
real ZCOP,SCOP,PI,DRIFT_P
data pi /3.1415/
logical DFLG,FH
character*80 BNKTYP
data BNKTYP/'2I4,9R4,2I4,R4,2I4,2R4,I4,XXX(I4,R4)'/
real ZCOOR_HIT1,ZCOOR_HIT2,DISTR1,DISTR2 !for joined trk close to IR
c
C------ Start of code --------------------------------------------------
C
CSTART with success
STATUS = OFSUCC
c
C----- Take the BIGGEST number of DPRS ---------------------------------
C
ISTAT = BIGEST (IW,'DPRS',JOIN_NUM)
if (ISTAT.ne.YESUCC) then
call ERLOGR ('DFJOIN',ERWARN,0,0,'Cannot find DPRS')
goto 990
endif
c------ Locate The two bank --------------------------------------------
c
ISTAT = BLOCAT (IW,'DTFS',TRK1,IND_DTF1,INDDAT_DTF1)
if (ISTAT.ne.YESUCC) then
call ERLOGR ('DPRMERGE',ERWARN,0,0,'Cannot find DTFS')
goto 990
endif
ISTAT = BLOCAT (IW,'DTFS',TRK2,IND_DTF2,INDDAT_DTF2)
if (ISTAT.ne.YESUCC) then
call ERLOGR ('DPRMERGE',ERWARN,0,0,'Cannot find DTFS')
goto 990
endif
C----- Locate also DHRE bank for reversing the measured drift time ------
c
ISTAT = BLOCAT(IW,'DHRE',2,IND_DHR,INDDAT_DHR)
if (ISTAT.ne.YESUCC) then
call ERLOGR ('DPRMERGE',ERWARN,0,0,'Cannot locate DHRE')
goto 990
endif
DHRECO=IW(INDDAT_DHR+DHRNCO)
INDDAT_DHR=INDDAT_DHR+IW(INDDAT_DHR+DHRHDS)
C Get COP Number:
ICOP= GETCOP(TRK1,TRK2)
POINT_TR1=INDDAT_DTF1+IW(INDDAT_DTF1+DTFHDS)
POINT_TR2=INDDAT_DTF2+IW(INDDAT_DTF2+DTFHDS)
if(FHIT(ICOP).eq.1) then !ICOP is broken in the First hit
FH=.true.
XCOOR_HIT1= RW(POINT_TR1+DTFXLP)
YCOOR_HIT1= RW(POINT_TR1+DTFYLP)
XCOOR_HIT2= RW(POINT_TR2+DTFXLP)
YCOOR_HIT2= RW(POINT_TR2+DTFYLP)
else !ICOP is broken in the Last hit
FH=.false.
XCOOR_HIT1= RW(POINT_TR1+DTFXCO)
YCOOR_HIT1= RW(POINT_TR1+DTFYCO)
ZCOOR_HIT1= RW(POINT_TR1+DTFZCO)
XCOOR_HIT2= RW(POINT_TR2+DTFXCO)
YCOOR_HIT2= RW(POINT_TR2+DTFYCO)
ZCOOR_HIT2= RW(POINT_TR2+DTFZCO)
endif
DISTXY1=sqrt(XCOOR_HIT1**2+YCOOR_HIT1**2)
DISTXY2=sqrt(XCOOR_HIT2**2+YCOOR_HIT2**2)
C! tracks with endpoints close to the IR- compare tridim r:
if(DISTXY1.lt.50.and.DISTXY2.lt.50) then
DISTR1 = sqrt(DISTXY1**2+ZCOOR_HIT1**2)
DISTR2 = sqrt(DISTXY2**2+ZCOOR_HIT2**2)
DFLG =DISTR1 .lt. DISTR2
else
DFLG = DISTXY1.lt.DISTXY2
endif
if (DFLG) then !lastpoint of trk1 is the first hit !!!
POINT_FIRSTRK= POINT_TR1
POINT_SECTRK = POINT_TR2
FTRK=TRK1
STRK=TRK2
else !lastpoint of trk2 is the first hit !!!
POINT_FIRSTRK= POINT_TR2
POINT_SECTRK = POINT_TR1
FTRK=TRK2
STRK=TRK1
endif
NHITPOS_FTRK=IW(POINT_FIRSTRK+DTFNHF)/10000
NHITNEG_FTRK=mod(IW(POINT_FIRSTRK+DTFNHF),10000)
NHIT_FTRK=NHITPOS_FTRK+NHITNEG_FTRK
NHITPOS_STRK=IW(POINT_SECTRK+DTFNHF)/10000
NHITNEG_STRK=mod(IW(POINT_SECTRK+DTFNHF),10000)
NHIT_STRK=NHITPOS_STRK+NHITNEG_STRK
NHITD=NHIT_FTRK+NHIT_STRK
F1=float(NHIT_FTRK)/NHITD
F2=float(NHIT_STRK)/NHITD
C----- Create the 2nd version of DPRS -----------------------------------
c
C Now make the corresponding DPRS bank, and fill the header
WRITE (BNKTYP(27:29),'(I3)') NHITD
STATUS = BTYMKG(IW,'DPRS',JOIN_NUM+1,
& BNKTYP,BNKSIZ,IND_DP,INDDAT_DP)
if (STATUS.NE.YESUCC) then
CALL ERLOGR('DFJOIN',ERWARN,0,0,
+ 'Error creating bank DPRS')
goto 990
endif
BNKS(ICOP)=JOIN_NUM+1
if(FH) then
C Fill bank DPRS first (for the first HIT)
C ====================
IW(INDDAT_DP+DPRHDS) = 2
IW(INDDAT_DP+DPRVRN) = 2
INDDAT_DP = INDDAT_DP + IW(INDDAT_DP+DPRHDS)
RW(INDDAT_DP+DPRXCO) = RW(POINT_FIRSTRK+DTFXLP) !
RW(INDDAT_DP+DPRYCO) = RW(POINT_FIRSTRK+DTFYLP) !
RW(INDDAT_DP+DPRZCO) = RW(POINT_FIRSTRK+DTFZLP) !
RW(INDDAT_DP+DPRXLS) = RW(POINT_SECTRK+DTFXLP)
RW(INDDAT_DP+DPRYLS) = RW(POINT_SECTRK+DTFYLP)
RW(INDDAT_DP+DPRZLS) = RW(POINT_SECTRK+DTFZLP)
RW(INDDAT_DP+DPRCUR) = -RW(POINT_FIRSTRK+DTFCUR)*F1+
& RW(POINT_SECTRK+DTFCUR)*F2 ! invert first track sign
RW(INDDAT_DP+DPRPHI) = RW(POINT_FIRSTRK+DTFPLP) - PI !F(-PI;PI) for FHIT
RW(INDDAT_DP+DPRCTT) = -RW(POINT_FIRSTRK+DTFTLP)*F1 +
& RW(POINT_SECTRK+DTFCTT)*F2
IW(INDDAT_DP+DPRPOS) = NHITPOS_FTRK+NHITPOS_STRK !
IW(INDDAT_DP+DPRNEG) = NHITNEG_FTRK+NHITNEG_STRK !
RW(INDDAT_DP+DPRQUA) = 0.
IW(INDDAT_DP+DPRPOK) = 3
IW(INDDAT_DP+DPRLRG) = 0
RW(INDDAT_DP+DPRPHF) = 0.
RW(INDDAT_DP+DPRPHL) = 0.
IW(INDDAT_DP+DPRCEN) = 0
INDDAT_DP=INDDAT_DP+16 !point to the hits block
C Now First TRK
do I=1,NHIT_FTRK-1 !I in the DPRS sense
IWIRE=IW(POINT_FIRSTRK+DTFFHA+NHIT_FTRK-I)
IW(INDDAT_DP+DPRLWI*(I-1)+DPRIWI)=IWIRE
W_PREV=IW(POINT_FIRSTRK+DTFFHA+NHIT_FTRK-I-1) !take prec. one
RW(INDDAT_DP+DPRLWI*(I-1)+DPRWTR)=GETLENG(FTRK,W_PREV,STATUS)
if(STATUS.ne.OFSUCC) then
CALL ERLOGR('DFJOIN',ERWARN,0,0,
+ 'Error filling DPRS')
goto 990
endif
C reverse drift distance sign:
DRIFT_P=RW(INDDAT_DHR+DHRECO*(IWIRE-1)+DHRRAD)
RW(INDDAT_DHR+DHRECO*(IWIRE-1)+DHRRAD)=-DRIFT_P
enddo
C Now the COP (wire position:NHIT_FTRK):
C For the moment the simple relation:s=dz/cos(theta)
IW(INDDAT_DP+DPRLWI*(NHIT_FTRK-1)+DPRIWI)=
& IW(POINT_FIRSTRK+DTFFHA)
DRIFT_P=RW(INDDAT_DHR+DHRECO*(IW(POINT_FIRSTRK+DTFFHA)-1)
& +DHRRAD)
RW(INDDAT_DHR+DHRECO*(IW(POINT_FIRSTRK+DTFFHA)-1)+
& +DHRRAD)=-DRIFT_P
ZCOP=RW(POINT_FIRSTRK+DTFZCO)-RW(POINT_SECTRK+
& DTFZCO)
COTFIRST = RW(POINT_FIRSTRK+DTFCTT)
SINFIRST = SQRT(1./(1.+COTFIRST**2))
COSFIRST = SINFIRST*COTFIRST
COTSEC = RW(POINT_SECTRK+DTFCTT)
SINSEC = SQRT(1./(1.+COTSEC**2))
COSSEC = SINSEC*COTSEC
COSTHETA=-COSFIRST*F1+COSSEC*F2
if(COSTHETA.gt.1.) COSTHETA=1.
if(COSTHETA.lt.-1.) COSTHETA=-1.
SCOP=ABS(ZCOP/COSTHETA*SINFIRST) !2D, AS IN DFINIT
RW(INDDAT_DP+DPRLWI*(NHIT_FTRK-1)+DPRWTR)=SCOP
C Now Second TRK
do I=NHIT_FTRK+1,NHITD
IWIRE=IW(POINT_SECTRK+DTFFHA+(I-NHIT_FTRK)-1)
IW(INDDAT_DP+DPRLWI*(I-1)+DPRIWI)=IWIRE
RW(INDDAT_DP+DPRLWI*(I-1)+DPRWTR)=GETLENG(STRK,IWIRE,STATUS)
if(STATUS.ne.OFSUCC) then
CALL ERLOGR('DFJOIN',ERWARN,0,0,
+ 'Error filling DPRS')
goto 990
endif
enddo
else
C Fill bank DPRS first (for the LAST HIT)
C ====================
IW(INDDAT_DP+DPRHDS) = 2
IW(INDDAT_DP+DPRVRN) = 2
INDDAT_DP = INDDAT_DP + IW(INDDAT_DP+DPRHDS)
RW(INDDAT_DP+DPRXCO) = RW(POINT_FIRSTRK+DTFXCO) !
RW(INDDAT_DP+DPRYCO) = RW(POINT_FIRSTRK+DTFYCO) !
RW(INDDAT_DP+DPRZCO) = RW(POINT_FIRSTRK+DTFZCO) !
RW(INDDAT_DP+DPRXLS) = RW(POINT_SECTRK+DTFXCO)
RW(INDDAT_DP+DPRYLS) = RW(POINT_SECTRK+DTFYCO)
RW(INDDAT_DP+DPRZLS) = RW(POINT_SECTRK+DTFZCO)
RW(INDDAT_DP+DPRCUR) = +RW(POINT_FIRSTRK+DTFCUR)*F1-
& RW(POINT_SECTRK+DTFCUR)*F2 ! invert the sign of the first trk
RW(INDDAT_DP+DPRPHI) = RW(POINT_FIRSTRK+DTFPHI)
RW(INDDAT_DP+DPRCTT) = RW(POINT_FIRSTRK+DTFCTT)*F1 -
& RW(POINT_SECTRK+DTFTLP)*F2
IW(INDDAT_DP+DPRPOS) = NHITPOS_FTRK+NHITPOS_STRK
IW(INDDAT_DP+DPRNEG) = NHITNEG_FTRK+NHITNEG_STRK
RW(INDDAT_DP+DPRQUA) = 0.
IW(INDDAT_DP+DPRPOK) = 3
IW(INDDAT_DP+DPRLRG) = 0
RW(INDDAT_DP+DPRPHF) = 0.
RW(INDDAT_DP+DPRPHL) = 0.
IW(INDDAT_DP+DPRCEN) = 0
INDDAT_DP=INDDAT_DP+16 !point to the hits block
C Now First TRK !I is in the DPRS direction
do I=1,NHIT_FTRK-1
IWIRE=IW(POINT_FIRSTRK+DTFFHA+I-1)
IW(INDDAT_DP+DPRLWI*(I-1)+DPRIWI)=IWIRE
RW(INDDAT_DP+DPRLWI*(I-1)+DPRWTR)=GETLENG(FTRK,IWIRE,STATUS)
if(STATUS.ne.OFSUCC) then
CALL ERLOGR('DFJOIN',ERWARN,0,0,
+ 'Error filling DPRS')
goto 990
endif
enddo
C Now the COP (wire position:NHIT_FTRK)
C For the moment the simple relation:s=dz/cos(theta)
ZCOP=RW(POINT_FIRSTRK+DTFZLP)-RW(POINT_SECTRK+DTFZLP)
IW(INDDAT_DP+DPRLWI*(NHIT_FTRK-1)+DPRIWI)=
& IW(POINT_FIRSTRK+DTFFHA+NHIT_FTRK-1)
COTFIRST = RW(POINT_FIRSTRK+DTFTLP)
SINFIRST = SQRT(1./(1.+COTFIRST**2))
COSFIRST = SINFIRST*COTFIRST
COTSEC = RW(POINT_SECTRK+DTFTLP)
SINSEC = SQRT(1./(1.+COTSEC**2))
COSSEC = SINSEC*COTSEC
COSTHETA=COSFIRST*F1-COSSEC*F2
if(COSTHETA.gt.1.) COSTHETA=1.
if(COSTHETA.lt.-1.) COSTHETA=-1.
SCOP=ABS(ZCOP/COSTHETA*SINFIRST) !2D, as in DFINIT
RW(INDDAT_DP+DPRLWI*(NHIT_FTRK-1)+DPRWTR)=SCOP
C Now second TRK
do I=NHIT_FTRK+1,NHITD-1
IWIRE=IW(POINT_SECTRK+DTFFHA+NHITD-I)
IW(INDDAT_DP+DPRLWI*(I-1)+DPRIWI)=IWIRE
W_PREV=IW(POINT_SECTRK+DTFFHA+NHITD-I-1) !take prec. one
RW(INDDAT_DP+DPRLWI*(I-1)+DPRWTR)=GETLENG(STRK,W_PREV,STATUS)
if(STATUS.ne.OFSUCC) then
CALL ERLOGR('DFJOIN',ERWARN,0,0,
+ 'Error filling DPRS')
goto 990
endif
C reverse drift distance sign:
DRIFT_P=RW(INDDAT_DHR+DHRECO*(IWIRE-1)+DHRRAD)
RW(INDDAT_DHR+DHRECO*(IWIRE-1)+DHRRAD)=-DRIFT_P
enddo
C Now Last HIT
C I=NHIT
IW(INDDAT_DP+DPRLWI*(NHITD-1)+DPRIWI)=
& IW(POINT_SECTRK+DTFFHA)
RW(INDDAT_DP+DPRLWI*(NHITD-1)+DPRWTR)=0.
DRIFT_P=RW(INDDAT_DHR+DHRECO*(IW(POINT_SECTRK+DTFFHA)-1)
& +DHRRAD)
RW(INDDAT_DHR+DHRECO*(IW(POINT_SECTRK+DTFFHA)-1)+
& DHRRAD)=-DRIFT_P
endif
900 continue
goto 999
990 continue
STATUS=OFFAIL
999 return
end
C=================================================================
C
REAL FUNCTION GETLENG(TRK,IWIRE,STATUS)
C
C=================================================================
C
C Description:
C ============
CIt retrives step lengh of the corresponding iwire
C
C Author:
C =======
C G.Venanzoni 26-9-97
C
C------ Global Declarations -------------------------------------------
C
$$IMPLICIT
C
C
C------ Include Files -------------------------------------------------
C
$$INCLUDE 'K$INC:ERLEVL.INC'
$$include 'k$inc:oferco.inc'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:STUMPCOM.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
C
C------ External Declaration --------------------------------------------
c
C------ Local Declarations -------------------------------------------
C
integer STATUS,KK,TRK,IWIRE
C----- Start of Code --------------------------------------------------
C
status=ofsucc
do 100 kk=1,mxhit
if(KDFS(TRK,KK).eq.IWIRE) then
GETLENG=DFSS(TRK,KK)
goto 999
endif
100 continue
900 continue
STATUS=OFFAIL
999 return
end
C=======================================================================
C
SUBROUTINE DFJSTUMP(COMSTATUS)
C----------------------------------------------------------------------
C IT scan all the track pairs to identify cases in which a track has been
C split in 2 (or more) stumps by the old ARGUS pattern routines.
C It fills the COPCOM & TRSCOM COmmon Block (in $TRK_INC:STUMPCOM.INC)
C
C Author:
C =======
C G.Venanzoni 19-9-97
C----------------------------------------------------------------------
C
C------ Global Declarations -------------------------------------------
C
$$IMPLICIT
C
C
C------ Include Files -------------------------------------------------
C
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$include 'k$inc:oferco.inc'
$$INCLUDE 'K$ITRK:DTFS.INC'
$$INCLUDE 'K$ITRK:STUMPCOM.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC' <ADD 11-11kink
C
C------ External Declaration --------------------------------------------
c
INTEGER BLOCAT,BNEXT,BDATA,BNUMB
REAL VDOTN
C------ Local Declarations -------------------------------------------
C
logical LFHIT,CUTFLG,CUTLLG,MASSFLG
integer STATUS,STATUS_DTFS,IND,JND,INDTFS,JNDTFS
integer COMSTATUS,NHITP
real VXFTRKI,VYFTRKI,VZFTRKI,CURVI,CHARGI
real VXFTRKJ,VYFTRKJ,VZFTRKJ,CURVJ,CHARGJ
real VXLTRKI,VYLTRKI,VZLTRKI
real VXLTRKJ,VYLTRKJ,VZLTRKJ
real DFTIJXY,DLTIJXY,RCOOI(3),COTI,PHII,SINTHI,TNGI(3)
real RCOOJ(3),COTJ,PHIJ,SINTHJ,TNGJ(3),DTIJZ
real COSX,XANG,RMED(2),ORTI,ORTJ,ORT
real PTTI,PXTI,PYTI,PZTI,ENI
real PTTJ,PXTJ,PYTJ,PZTJ,ENJ
real PIJ2,KLMASS,PI,MPI
DATA PI/3.1415/,MPI/139.570/
integer STATUS_PR,IND_PR,INDAT_PR,DPVERI !<ADD11-11 kink
integer JND_PR,JNDAT_PR,DPVERJ !<ADD11-11 kink
C----- Start of Code --------------------------------------------------
C
C Start with Success
NCOP=0
C NTI=0
C NTJ=0
COMSTATUS = OFSUCC
call VZERO(FIRSTUMP,MAXNUMNCOP)
call VZERO(SECSTUMP,MAXNUMNCOP)
call VZERO(FHIT,MAXNUMNCOP)
call VZERO(LHIT,MAXNUMNCOP)
call VZERO(BNKS,MAXNUMNCOP) !New entry 10-10
call VZERO(NHIT,MAXNUMTRS)
call VZERO(FIRSTHIT,MAXNUMTRS)
call VZERO(LASTHIT,MAXNUMTRS)
call VZERO(COMP_F,MAXNUMTRS)
call VZERO(COMP_L,MAXNUMTRS)
C Try to identify ARGUS-splits from Pattern Recognition
STATUS_DTFS = BLOCAT(IW,'DTFS',-1,IND,INDTFS)
if (STATUS_DTFS .NE. YESUCC) THEN
call ERLOGR('DFJSTUMP',ERWARN,0,0,'Cannot find DTFS')
goto 900
endif
DO WHILE (STATUS_DTFS.eq.YESUCC.and.IND.GT.0)
STATUS=BDATA(IW,IND,INDTFS)
IF (STATUS.NE.YESUCC) GOTO 900
INDTFS = INDTFS + IW(INDTFS+DTFHDS)
STATUS = BNUMB(IW,IND,NTI)
C Store DPRS #vers for kinks !
STATUS_PR = BLOCAT(IW,'DPRS',NTI,IND_PR,INDAT_PR) !<ADD11-11 kink
if(STATUS_PR.NE.YESUCC) then !<ADD11-11 kink
write(6,*)'Error with join-kink' !<ADD11-11 kink
GOTO 900 !<ADD11-11 kink
endif !<ADD11-11 kink
DPVERI=IW(INDAT_PR+DPRVRN) !<ADD11-11 kink
VXFTRKI = RW(INDTFS + DTFXCO)
VYFTRKI = RW(INDTFS + DTFYCO)
VXLTRKI = RW(INDTFS + DTFXLP)
VYLTRKI = RW(INDTFS + DTFYLP)
CURVI = RW(INDTFS + DTFCUR)
CHARGI = sign(1., CURVI)
STATUS_DTFS = BNEXT(IW,IND,JND)
IF (STATUS_DTFS.ne.YESUCC) GOTO 900
DO WHILE (STATUS_DTFS.eq.YESUCC.and.JND.GT.0)
STATUS = BDATA(IW,JND,JNDTFS)
IF (STATUS.NE.YESUCC) GOTO 900
JNDTFS = JNDTFS + IW(JNDTFS+DTFHDS)
STATUS = BNUMB(IW,JND,NTJ)
C Store DPRS #vers for kinks !
STATUS_PR = BLOCAT(IW,'DPRS',NTJ,JND_PR,JNDAT_PR) !<ADD11-11 kink
if(STATUS_PR.NE.YESUCC) then !<ADD11-11 kink
write(6,*)'Error with join-kink' !<ADD11-11 kink
GOTO 900 !<ADD11-11 kink
endif !<ADD11-11 kink
DPVERJ=IW(JNDAT_PR+DPRVRN) !<ADD11-11 kink
C Skip kinks !
if(dpverj.gt.10 .and. !<ADD11-11 kink
+ (dpverj .eq. dpveri)) goto 100 !<ADD11-11 kink
VXFTRKJ = RW(JNDTFS +DTFXCO)
VYFTRKJ = RW(JNDTFS + DTFYCO)
VXLTRKJ = RW(JNDTFS + DTFXLP)
VYLTRKJ = RW(JNDTFS + DTFYLP)
CURVJ = RW(JNDTFS + DTFCUR)
CHARGJ = sign(1., CURVJ)
C**** 0) First the selection on total charge.and.on curvature
if(ABS(CURVI+CURVJ).gt.0.5.or.(CHARGI+CHARGJ).ne.0) goto 100
***** 1) (x,y)distance
dftijxy=sqrt((vxftrkJ-vxftrkI)**2 +
+ (vyftrkJ-vyftrkI)**2)
dltijxy=sqrt((vxltrkJ-vxltrkI)**2 +
+ (vyltrkJ-vyltrkI)**2)
LFHIT = (DFTIJXY.lt.DLTIJXY)
if(LFHIT) then ! split point is in the first hit
RCOOI(1) = RW(INDTFS + DTFXCO)
RCOOI(2) = RW(INDTFS + DTFYCO)
RCOOI(3) = RW(INDTFS + DTFZCO)
COTI = RW(INDTFS + DTFCTT)
PHII = RW(INDTFS + DTFPHI)
SINTHI = SQRT(1./(1+COTI**2))
TNGI(1) = SINTHI*COS(PHII)
TNGI(2) = SINTHI*SIN(PHII)
TNGI(3) = SINTHI*COTI
PTTI = 1000./ABS(RW(INDTFS+DTFCUR))
RCOOJ(1) = RW(JNDTFS + DTFXCO)
RCOOJ(2) = RW(JNDTFS + DTFYCO)
RCOOJ(3) = RW(JNDTFS + DTFZCO)
COTJ = RW(JNDTFS + DTFCTT)
PHIJ = RW(JNDTFS + DTFPHI)
SINTHJ = SQRT(1./(1+COTJ**2))
TNGJ(1) = SINTHJ*COS(PHIJ)
TNGJ(2) = SINTHJ*SIN(PHIJ)
TNGJ(3) = SINTHJ*COTJ
PTTJ = 1000./ABS(RW(JNDTFS+DTFCUR))
else !LAST HIT
RCOOI(1) = RW(INDTFS + DTFXLP)
RCOOI(2) = RW(INDTFS + DTFYLP)
RCOOI(3) = RW(INDTFS + DTFZLP)
COTI = RW(INDTFS + DTFTLP)
PHII = RW(INDTFS + DTFPLP)
SINTHI = SQRT(1./(1+COTI**2))
TNGI(1) = SINTHI*COS(PHII)
TNGI(2) = SINTHI*SIN(PHII)
TNGI(3) = SINTHI*COTI
PTTI = 1000./ABS(RW(INDTFS+DTFCLP))
RCOOJ(1) = RW(JNDTFS + DTFXLP)
RCOOJ(2) = RW(JNDTFS + DTFYLP)
RCOOJ(3) = RW(JNDTFS + DTFZLP)
COTJ = RW(JNDTFS + DTFTLP)
PHIJ = RW(JNDTFS + DTFPLP)
SINTHJ = SQRT(1./(1+COTJ**2))
TNGJ(1) = SINTHJ*COS(PHIJ)
TNGJ(2) = SINTHJ*SIN(PHIJ)
TNGJ(3) = SINTHJ*COTJ
PTTJ = 1000./ABS(RW(JNDTFS+DTFCLP))
endif
C***** 2) Z Distance:
DTIJZ = RCOOJ(3)-RCOOI(3)
C***** 3)angle between <ti,-tj>
cosx = vdotn(tngi,tngj,3)
if(cosx.gt.1.) cosx=1.
if(cosx.lt.-1.) cosx=-1.
xang=(pi-acos(cosx))*180./pi
***** 4) ortogonality w.r.t (x,y)
RMED(1) =0.5*(RCOOI(1)+RCOOJ(1))
RMED(2) =0.5*(RCOOI(2)+RCOOJ(2))
ORTI = VDOTN(TNGI,RMED,2)
ORTJ = VDOTN(TNGJ,RMED,2)
ORT = ABS(ORTI+ORTJ)
***** 5) invariant mass
PXTI = PTTI * COS(PHII) ! Px PAC
PYTI = PTTI * SIN(PHII) ! Py PAC
PZTI = PTTI * COTI ! Pz PAC
PXTJ = PTTJ * COS(PHIJ) ! Px PAC
PYTJ = PTTJ * SIN(PHIJ) ! Py PAC
PZTJ = PTTJ * COTJ ! Pz PAC
ENI = SQRT(PXTI*PXTI+PYTI*PYTI+PZTI*PZTI+
& MPI*MPI)
ENJ = SQRT(PXTJ*PXTJ+PYTJ*PYTJ+PZTJ*PZTJ+
& MPI*MPI)
PIJ2 = (PXTI+PXTJ)*(PXTI+PXTJ)+
& (PYTI+PYTJ)*(PYTI+PYTJ)+(PZTI+PZTJ)*(PZTI+PZTJ)
KLMASS = SQRT((ENI+ENJ)*(ENI+ENJ)-PIJ2)
C cut on the i.m. (to be sure)
MASSFLG =KLMASS.lt.480.or.KLMASS.gt.550.
CUTFLG = LFHIT.and.(ABS(DTIJZ).lt.20.).and.(DFTIJXY.lt.20.).and.
+ (xang.lt.15.).and.(ort.lt.0.4).and.MASSFLG
CUTLLG = (.not.LFHIT).and.(ABS(DTIJZ).lt.20.).and.(DLTIJXY.lt.20.).and.
+ (xang.lt.20.).and.(ort.lt.0.3).and.MASSFLG
if (CUTFLG.or.CUTLLG) then ! these tracks are candidate stumps
NCOP = NCOP +1 !count the stumps
QF0(NCOP) = ABS(CURVI+CURVJ)/.5
NT = NTI
FIRSTUMP(NCOP) = NT
NHITP = IW(INDTFS+DTFNHF)
NHIT(NT) = NHITP/10000+MOD(NHITP,10000)
NT = NTJ
SECSTUMP(NCOP) = NT
NHITP = IW(JNDTFS+DTFNHF)
NHIT(NT) = NHITP/10000+MOD(NHITP,10000)
C FIRSTUMP(NCOP) = NTI
C NHIT(NTI) = IW(INDTFS+DTFNHF)
C SECSTUMP(NCOP) = NTJ
C NHIT(NTJ) = IW(JNDTFS+DTFNHF)
if (LFHIT) then ! FIRST HIT
FHIT(NCOP) = 1 ! ncop is in the first hit
QF1(NCOP) = ABS(DTIJZ)/15.
QF2(NCOP) = DFTIJXY/30.
QF3(NCOP) = XANG/15.
QF4(NCOP) = ORT/.4
NT = NTI
C FIRSTHIT(NTI) = 1
C FIRSTHIT(NTJ) = 1
C COMP_F(NTI) = NTJ
C COMP_F(NTJ) = NTI
FIRSTHIT(NT) = 1
COMP_F(NT) = NTJ
NT=NTJ
FIRSTHIT(NT) = 1
COMP_F(NT) = NTI
else
LHIT(NCOP) = 1 ! ncop is in the last hit
QF1(NCOP) = ABS(DTIJZ)/30.
QF2(NCOP) = DLTIJXY/30.
QF3(NCOP) = XANG/30.
QF4(NCOP) = ORT/.4
NT = NTI
LASTHIT(NT) = 1
COMP_L(NT) = NTJ
NT = NTJ
LASTHIT(NT) = 1
COMP_L(NT) = NTI
C LASTHIT(NTI) = 1
C LASTHIT(NTJ) = 1
C COMP_L(NTI) = NTJ
C COMP_L(NTJ) = NTI
endif
endif
100 continue
STATUS_DTFS = BNEXT(IW,JND,JND)
ENDDO
STATUS_DTFS = BNEXT(IW,IND,IND)
ENDDO
goto 990
900 COMSTATUS = OFFAIL
990 continue
return
end
C
C
C*DK DFKILL
C
C=======================================================================
C
SUBROUTINE DFKILL (NW)
C
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Mark wire "NW" as rejected during track fit, both
C- in the bank DPRS (by sign reversing) and in the
C- working vectors in /CDFWORK/ (by removing it from
C- the chain NDF)
C- NOTE : NW is NOT an index in the overall hit banks
C- DHRE,DTHA (which are to be left untouched)
C- but in the (smaller) one-track working banks
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:cdfcut.inc'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'C$INC:JOBSTA.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:cdfwork.inc'
integer npass2 !< AGG6-11
common /NP/ npass2 !< AGG6-11
INTEGER NW,N1,N0,NDIFF(MXHIT),TEST,STATUS
LOGICAL DIFFER
NDFERR = 0
C------- ENOUGH WIRES LEFT ?
C------- THIS TEST MAY FAIL ONLY WHEN CALLED FROM DFREJ; WHEN CALLED
C------- FROM DFADD WE REMOVE JUST A HIT PREVIOUSLY ADDED TO THE TRACK
IF (TIDST(LDFL(NW)) .LT. 0.) NFITNEG = NFITNEG-1
IF (TIDST(LDFL(NW)) .GT. 0.) NFITPOS = NFITPOS-1
IF (NFITNEG+NFITPOS .LE. MIN0(NCUTSTEREO,6)) THEN
CALL ERLOGR('DFKILL',ERWARN,0,STATUS,
+ 'Too many hits rejected : only 5 hits remain.')
NDFERR = 8
RETURN
ENDIF
C------- ENOUGH DIFFERENT LAYERS ?
N1 = NDF1
NL_CR = 0
DO WHILE (N1.GT.0 .AND. NL_CR.LT.2)
IF (N1.NE.NW) THEN
DIFFER = .TRUE.
DO TEST=1,NL_CR
IF (LDFL(N1) .EQ. NDIFF(TEST)) DIFFER = .FALSE.
ENDDO
IF (DIFFER) THEN
NL_CR = NL_CR + 1
NDIFF(NL_CR) = LDFL(N1)
ENDIF
ENDIF
N1 = NDF(N1)
ENDDO
IF (NL_CR .LT. 2) THEN
CALL ERLOGR('DFKILL',ERWARN,0,STATUS,
+ 'number of different crossed layer <2')
NDFERR = 9
RETURN
ENDIF
NDFC = NDFC-1
C------- MARK THE WIRE IN THE BANK DPRS, BUT CAREFUL: ONLY A POINT
C------- WHICH WAS ORIGINALLY THERE CAN BE REMOVED. POINTS ADDED BY
C------- DFADD ARE NOT IN DPRS BY DEFINITION.
IF (NW .LE. IW(JDPRS+DPRPOS)+IW(JDPRS+DPRNEG))
.IW(JDPRS1+DPRLWI*(NW-1)+DPRIWI) = - IW(JDPRS1+DPRLWI*(NW-1)+DPRIWI)
C------- MARK THE WIRE IN THE WORKING AREA
IF (NW .EQ. NDF1) THEN
C------- 1ST WIRE REJECTED : SET NDF1 TO THE FIRST POINT NOT YET REJECTED
C------- AND DECREASE THE DFL STEPS.
N0=NDF1
NDF1 = NDF(NW)
NDF(NW) = 0
N1 = NDF1
DO WHILE (N1.GT.0)
DFL(N1) = DFL(N1)-DFS(NW)
N1 = NDF(N1)
ENDDO
DFCTGT = DFPZ(NDF1) / SQRT (DFPX(NDF1)**2 + DFPY(NDF1)**2)
DFSTHE = 1. / SQRT (1. + DFCTGT**2)
DFPHI = ATAN2 (DFPY(NDF1) , DFPX(NDF1))
DO WHILE (NDF(NDF1).GT.0 .AND. LDFL(NDF1).EQ.0)
DFL(NDF(NDF1)) = DFL(NDF1) + DFS(NDF1)
CALL DFNEXT (NDF1,NDF(NDF1),DFS(NDF1))
NDF1 = NDF(NDF1)
ENDDO
NDFMS(1) = NDF1
ELSEIF (NW .EQ. NDFL) THEN
C------- LAST WIRE REJECTED ; DECREASE NUMBER OF STEPS (NDFL)
N1 = NDF1
DO WHILE (NDF(N1).NE.NW .AND. NDF(N1).GT.0)
N1 = NDF(N1)
ENDDO
NDFL = N1
NDF(NDFL) = 0
NDF(NW) = 0
DFS(NDFL) = 0.
N1 = NDFL
DO WHILE (LDFL(N1).EQ.0 .AND. N1.GE.1)
IF (MDF(N1).NE.0) NDFDIM = NDFDIM - 2
MDF(N1) = 0
N1 = N1-1
ENDDO
ELSE
C------- REJECTED WIRE IS IN THE MIDDLE; SHORT-CIRCUIT POINTERS "NDF"
N1 = NDF1
DO WHILE (NDF(N1).NE.NW .AND. NDF(N1).GT.0)
N1 = NDF(N1)
ENDDO
NDF(N1) = NDF(NW)
NDF(NW) = 0
DFS(N1) = DFS(N1)+DFS(NW)
ENDIF
IF (npass2.eq.2 .and. LDFMUS.GE.0) THEN
IF (MDF(NW).NE.0.OR.NW.EQ.N0) THEN
CALL DFMUSC
ENDIF
ENDIF
RETURN
END
C
C*DK DFMUSC
C
C=======================================================================
C
SUBROUTINE DFMUSC
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : To prepare the trck fit in DFITER for a treatment
C- of Multiple Coulomb Scattering. The lateral
C- average displacement after a STEP is given by
C-
C- STEP 13.6
C- -------- * ----------- * SQRT(STEP/X0)
C- SQRT(3.) P[MeV]*BETA
C-
C- When the cumulative track length since the last
C- kink (or since track beginning) given in DFL
C- is such that the above quantity is comparable
C- to the drift chamber accuracy, CDFACU, the kink
C- will be prepared.
C- Inputs : Attention to DFBETA: it is set in DFINIT in the
C- pion mass hypothesis.
C- Created
C- Updated
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
C
INTEGER N,M,NPRHITS
REAL PMEV,S,FACTOR,FACTM2
COMMON /NBHITS/ NPRHITS
PMEV = 1000./ABS(DFKURV*DFSTHE)
FACTOR = 13.6/PMEV/DFBETA
FACTM2 = 1./FACTOR/FACTOR
DFMLEN = (3.*CDFACU*CDFACU*FACTM2*TIDRL(1))**0.333333
DFMLEN = AMAX1(DFMLEN,DFL(NDFL)/AMIN1(CDFNMS,DFNDGF))
IF (DFMLEN .GT. .25*DFL(NDFL) .AND. NDFDIM.EQ.6) THEN
RETURN
ENDIF
NDFDIM = 6
CALL VZERO(MDF,MXHIT)
N = NDF(NDF1)
S = 0.
DO WHILE (N.GT.0.AND.N.NE.NDFL)
IF (DFL(N) - S .GE. DFMLEN) THEN
IF (NDFDIM.GT.60) THEN
CALL INERR(442)
GOTO 30
ENDIF
NDFDIM = NDFDIM + 2
MDF(N) = NDFDIM
NDFMS(NDFDIM-6) = N
NDFMS(NDFDIM-5) = N
DFPMS(NDFDIM-6) = 0.
DFPMS(NDFDIM-5) = 0.
C--- DFEMS = 1/(THETA M.S.)**2 - CF NIM(56)
DFEMS(NDFDIM-6) = FACTM2 * TIDRL(1)/(DFL(N)-S)
DFEMS(NDFDIM-5) = FACTM2 * TIDRL(1)/(DFL(N)-S)
S = DFL(N)
ENDIF
N = NDF(N)
ENDDO
MDF(NDFL) = 0
30 NDFMS(NDFDIM-4) = 0
DO N=1,MXHIT
DO M=6,NDFDIM
DFA(M,N) = 0.
ENDDO
ENDDO
RETURN
END
C
C*DK DFNEWW
C
C=======================================================================
C
SUBROUTINE DFNEWW (NERRF1)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : To add new hits to an already found track,
C- looking first in the direction of the particle,
C- and in a second iteration backward.
C
C Input parameters
C ================
INTEGER NERRF1 ! error multiplier from DFFIT to DFADD to DFITER
C
C Output parameters
C =================
C None
C
C Author:
C =======
C-
C- Originally from ARGUS collaboration
C- Modified June 1996 A. Calcaterra
C-
C----------------------------------------------------------------------
C
C Global Declarations:-
C =====================
C
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DHPT.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
C
C External Functions:-
C ====================
INTEGER BLOCAT
REAL PROXIM
C Local Declarations:-
C ====================
INTEGER MAXJMP
PARAMETER (MAXJMP=2)
INTEGER STATUS,INDR,INDDATR,DHRENCO,INDP,INDDATP
INTEGER NHITS
INTEGER MNEXT,WTEST
INTEGER ISIGN,N,LL
INTEGER MP,MM,M,LD,L1,L2,L
REAL P,PHIW,PHIHIT,PHIDEL,PHITST,CROSS,DOT
INTEGER MINC,SCAN,NDFCSV
INTEGER MSTRT(2),MSTOP(2),MSTEP(2)
INTEGER LFRST,LLAST
INTEGER ITRN
INTEGER NEWLR
LOGICAL LOG1,LOG2
C Locate bank DHRE, the running version of it
STATUS = BLOCAT(IW,'DHRE',2,INDR,INDDATR)
IF (STATUS.NE.YESUCC) THEN
CALL ERLOGR('DFNEWW',ERWARN,0,STATUS,
+ 'Bank DHRE not found.')
RETURN
ENDIF
INDR = INDDATR+IW(INDDATR+DHRHDS)
NHITS = IW(INDDATR+DHRNRO)
DHRENCO = IW(INDDATR+DHRNCO)
C Locate bank of DHPT pointers
STATUS = BLOCAT(IW,'DHPT',1,INDP,INDDATP)
IF (STATUS.NE.YESUCC) THEN
CALL ERLOGR('DFNEWW',ERWARN,0,STATUS,
+ 'Bank DHPT not found.')
RETURN
ENDIF
INDP = INDDATP+IW(INDDATP+DHPHDS)
C ===== INITIALISE LOOP OVER HIT WIRES
ISIGN = 1
C ===== AT ITERATION 1 THIS LOOP WORKS IN THE TRACK DIRECTION,
C ===== AND IT EXECUTES ONLY ONCE.
C ===== AT ITERATION 2 IT TRIES TO EXTEND THE TRACK AT THE INNER END
C AND IT REPEATS UNTIL NO NEW POINTS HAVE BEEN ADDED.
100 CONTINUE
NDFCSV = NDFC
N = NDF1
DO WHILE (N.GT.0)
LL = LDFL(N)
IF (LL.GT.0) THEN
C ===== MP = SIGN OF D AZIMUTH / D S (S = TRACK LENGTH)
MP = ISIGN
CROSS = (DFPY(N)*DFX(N)-DFPX(N)*DFY(N)) /
. SQRT(DFX(N)*DFX(N)+DFY(N)*DFY(N))
IF (CROSS .LT. 0.) MP = - MP
MM = KDF(N)
PHIW = TIDDP(LL)*(IDF(N)-1)
MINC = IW(INDP+DHPAPT+NHITS+MM-1)
LFRST = IW(INDP+DHPPT1+LL-1)
LLAST = IW(INDP+DHPPTL+LL-1)
IF (MP.EQ.1) THEN
MSTRT(1) = MINC+1
MSTOP(1) = LLAST
MSTEP(1) = 1
MSTRT(2) = MINC-1
MSTOP(2) = LFRST
MSTEP(2) = -1
ELSE
MSTRT(1) = MINC-1
MSTOP(1) = LFRST
MSTEP(1) = -1
MSTRT(2) = MINC+1
MSTOP(2) = LLAST
MSTEP(2) = 1
ENDIF
DO SCAN=1,2
DO M = MSTRT(SCAN),MSTOP(SCAN),MSTEP(SCAN)
MNEXT = IW(INDP+DHPAPT+M-1)
ITRN = IW(INDR+(MNEXT-1)*DHRENCO+DHRTRN)
IF (ITRN.LE.0) THEN
NEWLR = IW(INDR+(MNEXT-1)*DHRENCO+DHRSLR)
WTEST = IW(INDR+(MNEXT-1)*DHRENCO+DHRWNR)
PHITST = TIDDP(LL)*(WTEST-1)
PHIDEL = PROXIM(PHITST-PHIW,0.)
IF (ABS(PHIDEL) .LT. CDFPHI*TIDDP(LL)) THEN
IF (LL.NE.NEWLR) THEN
CALL ERLOGR('DFNEWW',ERWARN,0,STATUS,
+ 'Pointer crossing: a BAD error.')
ENDIF
CALL DFADD (MNEXT,NERRF1)
ENDIF
ENDIF
ENDDO
ENDDO
C ===== SEARCH IN THE NEXT LAYERS
DOT = (DFPX(N)*DFX(N)+DFPY(N)*DFY(N)) /
. SQRT(DFX(N)*DFX(N)+DFY(N)*DFY(N))
IF (DOT .GT. 0.) THEN ! TRACK IS PROPAGATING OUTWARD...
IF (ISIGN.EQ.1) THEN ! ...AND WE ARE LOOKING OUTWARD,
L1 = MAX0(LL-1,1) ! SO DECREASE LAYER INDEX!
L2 = MAX0(LL-MAXJMP,1) ! (ARGUS CONVENTION!)
LD = -1
ELSE ! ...BUT WE ARE LOOKING INWARD!
L1 = MIN0(LL+1,58) ! SO INCREASE LAYER INDEX INSTEAD
L2 = MIN0(LL+MAXJMP,58)
LD = +1
ENDIF
ELSE ! TRACK IS PROPAGATING INWARD
IF (ISIGN.EQ.1) THEN ! (AND SO WE LOOK FIRST)
L1 = MIN0(LL+1,58) ! SO DECREASE LAYER INDEX!
L2 = MIN0(LL+MAXJMP,58) ! (ARGUS CONVENTION!)
LD = +1
ELSE ! LATER, WILL LOOK OUTWARD
L1 = MAX0(LL-1,1) ! SO DECREASE LAYER INDEX. EASY. ;-)
L2 = MAX0(LL-MAXJMP,1)
LD = -1
ENDIF
ENDIF
PHIHIT = ATAN2(DFY(N),DFX(N)) - DFZ(N)*TIDST(LL)/TIDRA(LL)
IF (PHIHIT.LT.0.) P = P + PI2
C ===== NEXT LAYER = CHECK INNERMOST AND OUTERMOST LAYER
DO L=L1,L2,LD
IF (L.EQ.LL) GOTO 160
LFRST = IW(INDP+DHPPT1+L-1)
LLAST = IW(INDP+DHPPTL+L-1)
IF (MP.EQ.1) THEN
MSTRT(1) = LFRST
MSTOP(1) = LLAST
MSTEP(1) = 1
MSTRT(2) = LLAST
MSTOP(2) = LFRST
MSTEP(2) = -1
ELSE
MSTRT(1) = LLAST
MSTOP(1) = LFRST
MSTEP(1) = -1
MSTRT(2) = LFRST
MSTOP(2) = LLAST
MSTEP(2) = 1
ENDIF
DO M = MSTRT(1),MSTOP(1),MSTEP(1)
MNEXT = IW(INDP+DHPAPT+M-1)
ITRN = IW(INDR+(MNEXT-1)*DHRENCO+DHRTRN)
IF (ITRN.LE.0) THEN
WTEST = IW(INDR+(MNEXT-1)*DHRENCO+DHRWNR)
PHITST = TIDDP(L)*(WTEST-1)
PHIDEL = PROXIM(PHITST-PHIHIT,0.)
IF (MP.EQ.1) THEN
LOG1 = PHIDEL .GT. 0.
LOG2 = PHIDEL .LT. CDFPHI*TIDDP(L)/(1.2*ABS(DOT))
ELSE
LOG1 = PHIDEL .LT. 0.
LOG2 = PHIDEL .GT. -CDFPHI*TIDDP(L)/(1.2*ABS(DOT))
ENDIF
IF (LOG1 .AND. LOG2) THEN
CALL DFADD (MNEXT,NERRF1)
ENDIF
ENDIF
ENDDO
DO M = MSTRT(2),MSTOP(2),MSTEP(2)
MNEXT = IW(INDP+DHPAPT+M-1)
ITRN = IW(INDR+(MNEXT-1)*DHRENCO+DHRTRN)
IF (ITRN.LE.0) THEN
WTEST = IW(INDR+(MNEXT-1)*DHRENCO+DHRWNR)
PHITST = TIDDP(L)*(WTEST-1)
PHIDEL = PROXIM(PHITST-PHIHIT,0.)
IF (MP.EQ.1) THEN
LOG1 = PHIDEL .LT. 0.
LOG2 = PHIDEL .GT. -CDFPHI*TIDDP(L)/(1.2*ABS(DOT))
ELSE
LOG1 = PHIDEL .GT. 0.
LOG2 = PHIDEL .LT. CDFPHI*TIDDP(L)/(1.2*ABS(DOT))
ENDIF
IF (LOG1 .AND. LOG2) THEN
CALL DFADD (MNEXT,NERRF1)
ENDIF
ENDIF
ENDDO
160 ENDDO
ENDIF
170 N = NDF(N)
ENDDO
IF (ISIGN.EQ.-1) THEN
IF (NDFC .EQ. NDFCSV) RETURN
GO TO 100
ELSE
ISIGN = - 1
GO TO 100
ENDIF
END
C
C=======================================================================
SUBROUTINE DFDEDX(ISIGN)
C=======================================================================
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Calculation of the energy loss, put in DFDPDX,
C- used in the subroutine DFPNT
C- Created
C- Updated -G.Venanzoni 8/98
C- E. De Lucia
C- P. Pages
C-
C- INPUT PARAMETER : ISIGN (+ or - 1)
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
INTEGER ISIGN
REAL MOM,MASS,BETA
REAL ACF,ZCF,DAL,ALI
REAL DMOM,DENS
MOM=1./ABS(DFKURV*DFSTHE)
MASS=PION_MASS
BETA=MOM/SQRT(MASS**2+MOM**2)
DENS=0.427E-3
ACF=7.823
ZCF=3.965
DAL=0.000307*ZCF/ACF
ALI=62500./ZCF**0.9
DFDPDX=DAL/BETA**2
+ *(LOG(ALI*(MOM/MASS)**2)-BETA**2)*DENS
DFDPDX=DFDPDX/BETA
DMOM=MOM-DFL(NDFL)*DFDPDX
DFDPDX=SIGN(DFKURV**2*DFSTHE*DFDPDX,DFKURV)
if(isign.eq.-1) DFDPDX =-DFDPDX
RETURN
END
C
C*DK DFNEXT
SUBROUTINE DFNEXT (N1,N2,S)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Given the particle position and momentum
C- at a point "N1", a helyx segment of 3D-length S
C- is employed to compute position and momentum
C- at a point "N2". This routine implicitly takes
C- input from vectors DFX/Y/Z and DFPX/Y/Z at N1
C- and modifies the same vectors at N2.
C- IMPORTANT: "N1" and "N2" can also be the same.
C-
C- Created
C- Updated
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
INTEGER N1,N2
REAL S,SLAM,CLAM,C,PSI,SPSI,CPSI,A1
REAL VN(3)
REAL ANG1,ANG2,DFCTHE,PHIMS,THETMS
REAL P1,P2,NSLAM,NCLAM
COMMON /ANGLES/ ANG1(64),ANG2(64)
DFCTHE=SQRT(1-DFSTHE**2)
C-------- GET VECTOR VN = TANGENT CROSS H_FIELD ; CF KAPITZA(54)
VN(1) = DFPY(N1) * DFHZ(N1) - DFPZ(N1) * DFHY(N1)
VN(2) = DFPZ(N1) * DFHX(N1) - DFPX(N1) * DFHZ(N1)
VN(3) = DFPX(N1) * DFHY(N1) - DFPY(N1) * DFHX(N1)
C-------- LAM = LOCAL DIP ANGLE ; CF NIM(37)
SLAM = DFHX(N1)*DFPX(N1) + DFHY(N1)*DFPY(N1) + DFHZ(N1)*DFPZ(N1)
IF (ABS(SLAM).GT.0.999) THEN
CLAM = 0.
ELSE
CLAM = SQRT (1. - SLAM**2)
ENDIF
C-------- PSI = TURNING ANGLE OF ORBIT ; CF NIM(32,34)
C = DFH(N1) * (DFKURV + DFL(N2) * DFDPDX) * DFSTHE
PSI = C * S
SPSI = SIN (PSI)
CPSI = 1. - COS (PSI)
C-------- NEW XYZ ; CF KAPITZA(55)
A1 = (PSI - SPSI) * SLAM
DFX(N2) = DFX(N1) +
+ (SPSI*DFPX(N1) + CPSI*VN(1) + A1*DFHX(N1)) / C
DFY(N2) = DFY(N1) +
+ (SPSI*DFPY(N1) + CPSI*VN(2) + A1*DFHY(N1)) / C
DFZ(N2) = DFZ(N1) +
+ (SPSI*DFPZ(N1) + CPSI*VN(3) + A1*DFHZ(N1)) / C
C-------- NEW DIRECTION ; CF KAPITZA(56)
A1 = CPSI * SLAM
CPSI = 1. - CPSI
DFPX(N2) = CPSI * DFPX(N1) + SPSI * VN(1) + A1 * DFHX(N1)
DFPY(N2) = CPSI * DFPY(N1) + SPSI * VN(2) + A1 * DFHY(N1)
DFPZ(N2) = CPSI * DFPZ(N1) + SPSI * VN(3) + A1 * DFHZ(N1)
C-------- MULTIPLE SCATTERING
IF (MDF(N2).NE.0 .AND. N1.NE.N2) THEN
NCLAM=CLAM*COS(ANG2(MDF(N2)-6))
+ +DFPZ(N2)*SIN(ANG2(MDF(N2)-6))
IF (ABS(NCLAM).LT.0.99999)THEN
NSLAM=SQRT(1.-NCLAM**2)
ELSE
NSLAM=0.
ENDIF
CPSI=COS(PSI)
A1 = (1-CPSI) * NSLAM
DFPX(N2) = CPSI * DFPX(N1) + SPSI * VN(1)
+ + A1 * DFHX(N1)
DFPY(N2) = CPSI * DFPY(N1) + SPSI * VN(2)
+ + A1 * DFHY(N1)
DFPZ(N2) = DFPZ(N2)*COS(ANG2(MDF(N2)-6))
+ - CLAM*SIN(ANG2(MDF(N2)-6))
A1 = DFPX(N2)
DFPX(N2) = A1 * COS (ANG1(MDF(N2)-6)) -
+ DFPY(N2) * SIN (ANG1(MDF(N2)-6))
DFPY(N2) = A1 * SIN (ANG1(MDF(N2)-6)) +
+ DFPY(N2) * COS (ANG1(MDF(N2)-6))
A1 = DFPX(N2)**2 + DFPY(N2)**2 + DFPZ(N2)**2
A1 = 2. / (1. + A1)
DFPX(N2) = DFPX(N2) * A1
DFPY(N2) = DFPY(N2) * A1
DFPZ(N2) = DFPZ(N2) * A1
ENDIF
RETURN
END
C
C*DK DFPNT
SUBROUTINE DFPNT (N1,N2,S)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Given the particle position and unit vector
C- at the measured point "N1", the track is LINEARLY
C- extrapolated to the point of CA to wire N2.
C- Also, the magnetic field must be
C- recalculated, for now it is assumed uniform.
C-
C Input parameters
C ================
INTEGER N1,N2 ! index of the 2 hits in the 1-track common /CDFWORK/
C
C Output parameters
C =================
REAL S ! distance from point N1 to point of CA to N2.
C
C- Created
C- Updated
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
REAL WT,X(3),H(3)
REAL EPSILON /1.E-6/
DFH(N1) = TIFLAV
X(1) = DFX(N1)
X(2) = DFY(N1)
X(3) = DFZ(N1)
CALL MGFLD(X,H)
CALL VSCALE(H,CONVERS,H,3)
DFHX(N1) = H(1) / DFH(N1)
DFHY(N1) = H(2) / DFH(N1)
DFHZ(N1) = H(3) / DFH(N1)
IF (LDFL(N2).GT.0) THEN
C THIS IS A POINT WHERE A REAL DISTANCE MEASUREMENT EXISTS
WT = DFVX(N2)*DFPX(N1) + DFVY(N2)*DFPY(N1)
+ + TIDSC(LDFL(N2))*DFPZ(N1)
IF (ABS(1.-WT*WT) .LT. EPSILON) THEN
IF (WT .LE. -1.) THEN
WT = -1. - EPSILON
ELSEIF (WT .LE. 0.) THEN
WT = -1. + EPSILON
ELSEIF (WT .LE. 1.) THEN
WT = 1. - EPSILON
ELSE
WT = 1. + EPSILON
ENDIF
ENDIF
S = ((DFWX(N2) - DFX(N1)) * (DFPX(N1) - WT * DFVX(N2)) +
+ (DFWY(N2) - DFY(N1)) * (DFPY(N1) - WT * DFVY(N2)) -
+ DFZ(N1) * (DFPZ(N1) - WT * TIDSC(LDFL(N2)))) /
+ (1. - WT * WT)
ELSE
C THIS IS EITHER A POINT PREVIOUSLY REJECTED IN THE FIT ITERATIONS,
C OR A FICTITIOUS POINT, INTRODUCED TO TREAT MULTIPLE SCATTERING.
S = 0.
ENDIF
RETURN
END
C
C*DK DFLIP
C
C=======================================================================
C
SUBROUTINE DFLIP (NERRF1)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : In case where the drift distance is small
C (typically 3-4 mm) it is possible that
C the sign of the drift distance came out wrong
C from Pattern Recognition.
C This routine tries to flip all (small) drift
C distances, starting from the wire with the
C worst residual.
C
C Author: A. Calcaterra
C =======
C-
C- Created 10-14-96
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
INTEGER NERRF1
INTEGER STATUS,BLOCAT
INTEGER INDR,INDDATR,DHRENCO
INTEGER INDICES(MXHIT),FLIPPED(MXHIT)
LOGICAL IN_RANGE, BETTER, NOTFLP
LOGICAL SAV_HITS,ITERATE
INTEGER N,NSORT,NFLIP
REAL OTHER,PREV_R,CHISAV
REAL DMIN,DMAX
DATA DMIN/0.000/,DMAX/3./
REAL RSORT(MXHIT),DFOTH(MXHIT)
STATUS = BLOCAT(IW,'DHRE',2,INDR,INDDATR)
IF (STATUS.NE.YESUCC) RETURN
DHRENCO = IW(INDDATR+DHRNCO)
INDR = INDDATR+IW(INDDATR+DHRHDS)
CALL VZERO(FLIPPED,MXHIT)
C Select all wires with small drift distance in order of decreasing
C (absolute) normalized residual. Flip wires only once please !!
ITERATE = .TRUE.
*************************************************************
DO WHILE (ITERATE)
*************************************************************
NSORT = 0
N = NDF1
DO WHILE (N.GT.0)
OTHER = DFCH(N) - 2.*DFD(N)
IN_RANGE = ABS(DFD(N)) .LT. DMAX .AND.
+ ABS(DFD(N)) .GT. DMIN
BETTER = ABS(OTHER) .LT. ABS(DFCH(N))
NOTFLP = FLIPPED(N).EQ.0
IF (IN_RANGE .AND. BETTER .AND. NOTFLP) THEN
NSORT = NSORT+1
INDICES(NSORT) = N
RSORT(N) = ABS(DFCH(N))-ABS(OTHER)
DFOTH(N) = OTHER
ENDIF
N = NDF(N)
ENDDO
C Are there (more) candidates for sign flipping?
IF (NSORT.EQ.0) THEN
ITERATE = .FALSE.
GOTO 999
ENDIF
C Yes, do the sorting, flip the worst hit, and try to refit the track
CALL SORTZV(RSORT,INDICES,NSORT,1,-1,-1)
NFLIP = INDICES(1)
FLIPPED(NFLIP) = 1
C Let's save all info about track before calling DFITER....
SAV_HITS = .FALSE.
NDFERR = 0
CALL DFSAVE(NFLIP,SAV_HITS)
CHISAV = DFCHIN
PREV_R = RW(INDR+(KDF(NFLIP)-1)*DHRENCO+DHRRAD)
RW(INDR+(KDF(NFLIP)-1)*DHRENCO+DHRRAD) = -PREV_R
DFD(NFLIP) = -DFD(NFLIP)
100 CALL DFITER (NERRF1*CDFDC1,CDFM1L)
C Did something BAD happen ?
IF (NDFERR.NE.0) THEN
NERRF1=NERRF1+1
CALL INERR(408)
C Some error conditions: (Chi2 increasing) and (too many iterations)
C invite to try again the fit relaxing the first cut. The same could
C apply in principle for (pred chi2 wors), NDFERR=3
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) GO TO 200
IF (DFCHIN.LT.CHISAV) THEN
IF (NERRF1.LT.15) THEN
GOTO 100
ENDIF
ENDIF
200 CONTINUE
RW(INDR+(KDF(NFLIP)-1)*DHRENCO+DHRRAD) = PREV_R
C It's better to restore previous situation !
CALL DFRESTORE(NFLIP,SAV_HITS)
CALL HFILL(408,0.,0.,1.)
ELSE
CALL HFILL(409,(DFCHIN-CHISAV)/DFNDGF,0.,1.)
IF (DFCHIN.GT.CHISAV) THEN
RW(INDR+(KDF(NFLIP)-1)*DHRENCO+DHRRAD) = PREV_R
C It's better to restore previous situation!
CALL DFRESTORE(NFLIP,SAV_HITS)
CALL HFILL(408,1.,0.,1.)
ELSE
CALL HFILL(408,2.,0.,1.)
ENDIF
ENDIF
999 CONTINUE
*************************************************************
ENDDO
*************************************************************
RETURN
END
C
C*DK DFREJ
C
C=======================================================================
C
SUBROUTINE DFREJ (NERRF1,CHIMAX,CHIFAC)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : This deceptively short and simple routine does
C- one of the most complicate jobs I have ever seen:
C- it keeps trying to obtain better track fits by picking the point giving
C- the biggest contribution to the chi^2 and doing the following:
C- 1) test if the zeroth-order approximation to the new fit will give
C- a smaller overall chi^2.
C- 2) If yes, re-fit the track with DFITER
C- 3) If the worst point is now another one, do the same to the other point.
C- 4) Repeat until EITHER this worst-contribution is smaller than "CHIMAX"
C- OR it is smaller than "CHIFAC"*the average chi^2 / point.
C- 5) At that point, it is judged that the fit has already improved as much
C- as possible: ALL drift signs are switched to the most favorable case,
C- and a final fit will be remade in DFFIT.
C- 6) If, during the iteration, one of the reversals is anticipated as
C- UNSUCCESSFUL (NO at point 2) above), and the worst contribution STILL
C- is bigger than "CHIMAX", the point is killed and the fit is repeated
C- without that point.
C- 7) If, during the iteration, even after a successful fit the worst point
C- does not change, the test at point 4) is made. If it fails (meaning
C- still more improvement is requested) the point is killed and a new
C- "worst point" is calculated (via DFITER) and examined.
C- 8) In the original routine there was the possibility of an infinite loop
C- when the worst possible chi square contribution oscillates between
C- two numbers. This actually happened! So the ATTEMPTS are limited to
C- the actual number of points.
C-
C- Created
C- Updated
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'C$INC:JOBSTA.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' ! 9/6/98 MXHIT=250
INTEGER NERRF1
REAL CHIMAX,CHIFAC
INTEGER STATUS,BLOCAT
INTEGER INDR,INDDATR,DHRENCO
INTEGER N,NPRHITS
REAL CHIPART
REAL ATTEMPTS
REAL OTHER,PREV_R,CHISAV
LOGICAL SAV_HITS
COMMON /NBHITS/ NPRHITS
CALL VZERO(CHIPART,MXHIT)
STATUS = BLOCAT(IW,'DHRE',2,INDR,INDDATR)
IF (STATUS.NE.YESUCC) RETURN
DHRENCO = IW(INDDATR+DHRNCO)
INDR = INDDATR+IW(INDDATR+DHRHDS)
ATTEMPTS = 0.
10 IF (DFMCH.LT.CHIMAX .OR.
+ DFMCH.LT.DFCHIN*CHIFAC/(DFNDGF+5.) .OR.
+ ATTEMPTS.GT.DFNDGF+5.) GOTO 200
C------- TRY THE OTHER SIGN OF THE DRIFT DISTANCE
N = NDFMCH
CHISAV = DFCHIN
OTHER = DFCH(N) - 2.*DFD(N)
IF (ABS(DFCH(N)) .GE. ABS(OTHER)) THEN
ATTEMPTS = ATTEMPTS+1.
PREV_R = RW(INDR+(KDF(N)-1)*DHRENCO+DHRRAD)
RW(INDR+(KDF(N)-1)*DHRENCO+DHRRAD) = -PREV_R
C Let's save all info about track before calling DFITER....
SAV_HITS = .FALSE.
NDFERR = 0
CALL DFSAVE(N,SAV_HITS)
DFD(N) = -DFD(N)
DFCH(N) = OTHER
100 CALL DFITER (NERRF1*CDFDC1,CDFM1L)
IF (NDFERR.NE.0) THEN
NERRF1=NERRF1+1
CALL INERR(408)
C Some error conditions: (Chi2 increasing) and (too many iterations)
C invite to try again the fit relaxing the first cut. The same could
C apply in principle for (pred chi2 wors), NDFERR=3
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) GO TO 400
IF(DFCHIN.LT.CHISAV)THEN
IF (NERRF1.LT.15) THEN
GOTO 100
ENDIF
ENDIF
400 CONTINUE
RW(INDR+(KDF(N)-1)*DHRENCO+DHRRAD) = PREV_R
C It's better to restore previous situation !
CALL DFRESTORE(N,SAV_HITS)
ELSE
IF (DFCHIN.GT.CHISAV) THEN
RW(INDR+(KDF(N)-1)*DHRENCO+DHRRAD) = PREV_R
C It's better to restore previous situation !
CALL DFRESTORE(N,SAV_HITS)
ENDIF
ENDIF
IF (N.NE.NDFMCH) GOTO 10
ENDIF
C------- REJECT WIRE
IF (DFMCH.LT.CHIMAX .OR.
+ DFMCH.LT.DFCHIN*CHIFAC/(DFNDGF+5.)) GOTO 200
IF (DFMCH.LT.CHIMAX) GOTO 200
CHISAV = DFCHIN
CALL DFKILL(N)
IF (NDFERR.NE.0) THEN
RETURN
ENDIF
300 CALL DFITER (NERRF1*CDFDC1,CDFM1L)
IF (NDFERR.NE.0) THEN
NERRF1=NERRF1+1
CALL INERR(408)
C Some error conditions: (Chi2 increasing) and (too many iterations)
C invite to try again the fit relaxing the first cut. The same could
C apply in principle for (pred chi2 wors), NDFERR=3
IF (NDFERR.NE.2 .AND. NDFERR.NE.6) RETURN
IF (DFCHIN.LT.CHISAV)THEN
IF (NERRF1.LT.15) GOTO 300
ENDIF
ENDIF
GOTO 10
C------- TRY ALL OTHER SIGNS OF THE DRIFT DISTANCE
200 CONTINUE
RETURN
END
C
C*DK DFSTOR
C
C=======================================================================
C
SUBROUTINE DFSTOR(NPASS)
C
$$IMPLICIT
C FILL FIT RESULTS INTO BANK "DTFS"
C author : A. Antonelli + S. Calcaterra - march 1995
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC' !<
C
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DTFS.INC'
$$INCLUDE 'K$ITRK:DHSP.INC'
$$INCLUDE 'C$INC:JOBSTA.INC'
C External Functions:-
C ====================
INTEGER BTYMKG
REAL PROB
C Local Declarations:-
C ====================
INTEGER KINENUM,PIDNUM
INTEGER STATUS,IND,INDDAT,INDP,INDDATP,DHSPNCO
INTEGER BNKSIZ,II,N
INTEGER NHIPT,NHITF
INTEGER NPASS
REAL CH2N
LOGICAL REVERSE
REAL*4 XST(3),RTOBEAM,ZTOBEAM,S
CHARACTER*32 DTFSBNK
DATA DTFSBNK/'13I4,30R4,I4,R4,2I4,6R4,XXXI4'/
CHARACTER*15 DHSPBNK
DATA DHSPBNK/'4I4,XXX(I4,7R4)'/
C Executable Code:
C =================
NHIPT = IW(JDPRS+DPRPOS)+IW(JDPRS+DPRNEG)
NHITF = NFITNEG + NFITPOS
IF (NHITF.GT.999) THEN
CALL INERR(401)
RETURN
ENDIF
C Create the fitted track bank
C =============================
WRITE(DTFSBNK(25:27),'(I3)') NHITF
STATUS = BTYMKG(IW,'DTFS',ITR,DTFSBNK,BNKSIZ,IND,INDDAT)
IF (STATUS.NE.YESUCC) THEN
CALL INERR(490)
RETURN
ENDIF
C Create the space points bank
C =============================
WRITE(DHSPBNK(5:7),'(I3)') NHITF
STATUS = BTYMKG(IW,'DHSP',ITR,DHSPBNK,BNKSIZ,INDP,INDDATP)
IF (STATUS.NE.YESUCC) THEN
CALL INERR(491)
RETURN
ENDIF
C Fill bank DTFS first
C ====================
IW(INDDAT+DTFHDS) = 2
IW(INDDAT+DTFVRN) = NPASS
INDDAT = INDDAT + IW(INDDAT+DTFHDS)
IW(INDDAT+DTFINH) = 10000*IW(JDPRS+DPRPOS)+IW(JDPRS+DPRNEG)
IW(INDDAT+DTFNHF) = 10000*NFITPOS+NFITNEG
RW(INDDAT+DTFLNG) = DFL(NDFL)
RW(INDDAT+DTFXCO) = DFX(NDF1)
RW(INDDAT+DTFYCO) = DFY(NDF1)
RW(INDDAT+DTFZCO) = DFZ(NDF1)
RW(INDDAT+DTFCUR) = DFKURV
RW(INDDAT+DTFCTT) = DFCTGT
RW(INDDAT+DTFPHI) = DFPHI
RW(INDDAT+DTFCVM) = XI(1,1)
RW(INDDAT+DTFCVM+1) = XI(2,1)
RW(INDDAT+DTFCVM+2) = XI(3,1)
RW(INDDAT+DTFCVM+3) = XI(4,1)
RW(INDDAT+DTFCVM+4) = XI(6,1)
RW(INDDAT+DTFCVM+5) = XI(2,2)
RW(INDDAT+DTFCVM+6) = XI(3,2)
RW(INDDAT+DTFCVM+7) = XI(4,2)
RW(INDDAT+DTFCVM+8) = XI(6,2)
RW(INDDAT+DTFCVM+9) = XI(3,3)
RW(INDDAT+DTFCVM+10)= XI(4,3)
RW(INDDAT+DTFCVM+11)= XI(6,3)
RW(INDDAT+DTFCVM+12)= XI(4,4)
RW(INDDAT+DTFCVM+13)= XI(6,4)
RW(INDDAT+DTFCVM+14)= XI(6,6)
RW(INDDAT+DTFCHI) = DFCHIN
RW(INDDAT+DTFPCH) = PROB(DFCHIN,NHITF-5)
RW(INDDAT+DTFXLP) = DFX(NDFL)
RW(INDDAT+DTFYLP) = DFY(NDFL)
RW(INDDAT+DTFZLP) = DFZ(NDFL)
RW(INDDAT+DTFCLP) = ABS (DFKURV + DFL(NDFL) * DFDPDX)
IF (SQRT(DFPX(NDFL)**2 + DFPY(NDFL)**2).GT.1.E-20) THEN
RW(INDDAT+DTFTLP)= DFPZ(NDFL)/SQRT(DFPX(NDFL)**2+DFPY(NDFL)**2)
RW(INDDAT+DTFPLP)= ATAN2 (DFPY(NDFL),DFPX(NDFL))
ENDIF
IF(RW(INDDAT+DTFPLP).LT.0.)RW(INDDAT+DTFPLP)=RW(INDDAT+DTFPLP)+PI2
C------ MULTIPLE SCATTERING
IW(INDDAT+DTFMSK)= NDFDIM - 6
RW(INDDAT+DTFCHM)= 0.
IF (NDFDIM.GT.6) RW(INDDAT+DTFCHM)= DFCMS2
IW(INDDAT+DTFPVS)= 0
N = NDF1
II = 0
DO WHILE (N.GT.0)
IF (LDFL(N) .GT. 0) THEN
IW(INDDAT+DTFFHA+II)= KDF(N)
II = II+1
ENDIF
N = NDF(N)
ENDDO
C======= PROJECT THIS TRACK (IF IT HAS POINTS IN AN INNER LAYER)
C======= TOWARD THE BEAM AXIS. FLAG THE TRACK AS COMING FROM THE IP
ZTOBEAM = -1000.
RTOBEAM = 1000.
IF (LDFL(NDF1).EQ.57 .OR. LDFL(NDF1).EQ.58) THEN
REVERSE = .TRUE.
XST(1) = DFX(NDF1)
XST(2) = DFY(NDF1)
XST(3) = DFZ(NDF1)
CALL DFTOIP(XST,DFKURV,DFCTGT,DFPHI,REVERSE,ZTOBEAM,RTOBEAM,S)
IF (RTOBEAM.LE.3. .AND. ABS(ZTOBEAM).LE.10.) THEN
IW(INDDAT+DTFCEN) = 1
ENDIF
ENDIF
C And now bank DHSP
C =================
IW(INDDATP+DHSHDS) = 4
IW(INDDATP+DHSVRN) = 1
IW(INDDATP+DHSNRO) = NHITF
IW(INDDATP+DHSNCO) = 8
DHSPNCO = IW(INDDATP+DHSNCO)
INDP = INDDATP + IW(INDDATP+DHSHDS)
N = NDF1
DO WHILE (N.GT.0)
IF (LDFL(N) .GT. 0) THEN
IW(INDP+DHSADR) = IW(JDPRS1+DPRLWI*(N-1)+DPRIWI)
RW(INDP+DHSXCO) = DFX(N)
RW(INDP+DHSYCO) = DFY(N)
RW(INDP+DHSZCO) = DFZ(N)
CH2N = - DFCH(N)
DO II=1,NDFDIM
CH2N = CH2N + XX(II,NDFDIM+2)*DFA(II,N)
ENDDO
RW(INDP+DHSCHN) = CH2N
RW(INDP+DHSBET) = DFB(N)
RW(INDP+DHSPHT) = DFF(N)
RW(INDP+DHSDRD) = DFD(N)-DFCH(N)
INDP = INDP + DHSPNCO
ENDIF
N = NDF(N)
ENDDO
RETURN
END
C
C*DK DFTOIP
C
C=======================================================================
C
SUBROUTINE DFTOIP(XSTART,CURV,CTGT,PHI0,REVERSE,ZMIN,RMIN,S)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Given the particle position and unit vector
C- DFTOIP finds the point of closest approach of
C- the track to the beam axis.
C-
C Author: A. Calcaterra
C =======
C-
C- Created Mar 1st 1996
C-
C----------------------------------------------------------------------
REAL*4 XSTART(3),CURV,CTGT,PHI0,S
LOGICAL REVERSE
REAL*4 ZMIN,RMIN
$$INCLUDE 'K$ITRK:CCONST.INC'
REAL RHO,Q,XCEN,YCEN,STEP,DSTEP
REAL OLD_Z,OLD_R,XSTOP(3),NEW_R
LOGICAL HOMING
RHO = 1./ABS(CURV*TIFLAV)
Q = SIGN(1.,CURV)
XCEN = XSTART(1)+Q*RHO*SIN(PHI0)
YCEN = XSTART(2)-Q*RHO*COS(PHI0)
HOMING = .TRUE.
OLD_R = SQRT(XSTART(1)**2+XSTART(2)**2)
STEP = OLD_R/2.
IF (REVERSE) STEP = -STEP
DO WHILE (HOMING .AND. ABS(STEP).GT.0.1)
XSTOP(1) = XCEN-RHO*Q*SIN(PHI0-Q*STEP/RHO)
XSTOP(2) = YCEN+RHO*Q*COS(PHI0-Q*STEP/RHO)
XSTOP(3) = XSTART(3)+STEP*CTGT
NEW_R = SQRT(XSTOP(1)**2+XSTOP(2)**2)
DSTEP = NEW_R/2.
IF (REVERSE) DSTEP = -DSTEP
STEP = STEP+DSTEP
HOMING = NEW_R.LT.OLD_R
OLD_R = NEW_R
OLD_Z = XSTOP(3)
ENDDO
S = STEP
RMIN = OLD_R
ZMIN = OLD_Z
RETURN
END
C
C*DK DFTRAC
C
C=======================================================================
C
SUBROUTINE DFTRAC(FAST_SWIM)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Given the particle position and unit vector
C- at the measured point "NDF1", a succession
C- of helical segments is found, joining points
C- of closest approach to successive hit wires.
C- IMPORTANT : the length output from DFPNT has the intrinsic
C- meaning of a CORRECTION, and it should be "small".
C- The length input to DFNEXT is a CORRECTION, and
C- must come from DFPNT, when the first 2 arguments
C- are the same, is a STEP when they are different.
C-
C- Created
C- Updated
C-
C----------------------------------------------------------------------
LOGICAL FAST_SWIM
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
INTEGER N1,N2,M
REAL STEP, CORRECTION
REAL ANG1,ANG2
COMMON /ANGLES/ ANG1(64),ANG2(64)
IF (NDFDIM.GT.6) THEN
DO M=7,NDFDIM,2
ANG1(M-5) = DFPMS(M-5)
ANG2(M-5) = DFPMS(M-4)
ENDDO
ENDIF
NDFERR = 0
N1 = NDF1
N2 = NDF(N1)
IF (FDFAST .AND. FAST_SWIM) THEN
C-------- FAST TRACING
DO WHILE (N2.NE.0)
DFL(N2) = DFL(N1) + DFS(N1)
CALL DFNEXT(N1,N2,DFS(N1))
N1 = NDF(N1)
N2 = NDF(N2)
ENDDO
ELSE
C-------- SLOW TRACING : START CORRECTING THE POSITION OF THE 1ST POINT
CALL DFPNT (N1,N1,CORRECTION)
IF (ABS(CORRECTION) .GT. AMIN1 (0.2*DFSTLC,0.2)) THEN
DFL(N1) = DFL(N1) + CORRECTION
CALL DFNEXT (N1,N1,CORRECTION)
CALL DFPNT (N1,N1,CORRECTION)
DFCTGT = DFPZ(N1) / SQRT (DFPX(N1)**2 + DFPY(N1)**2)
DFSTHE = 1. / SQRT (1. + DFCTGT**2)
DFPHI = ATAN2 (DFPY(N1) , DFPX(N1))
ENDIF
DFL(N1) = 0.
STEP = DFS(N1)
DO WHILE (N2.NE.0)
DFL(N2) = DFL(N1) + STEP
CALL DFNEXT (N1,N2,STEP)
CALL DFPNT (N2,N2,CORRECTION)
IF (CORRECTION .LT. -CDFSTL) THEN
CALL ERLOGR('DFTRAC',ERWARN,0,0,
+ 'Large negative step length: wrong order in D3SORT.')
NDFERR = 7
RETURN
ENDIF
DFS(N1) = DFS(N1) + CORRECTION
IF (ABS(CORRECTION) .LE. DFSTLC) THEN
STEP = DFS(N2) + CORRECTION
ELSE
DFL(N2) = DFL(N2) + CORRECTION
CALL DFNEXT (N2,N2,CORRECTION)
CALL DFPNT (N2,N2,CORRECTION)
STEP = DFS(N2)
ENDIF
N1 = NDF(N1)
N2 = NDF(N2)
ENDDO
ENDIF
RETURN
END
C
C*DK INERR
C
C=======================================================================
C
SUBROUTINE INERR(IENUMB)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Definition of error codes
C-
C- Created 1-JUL-1992 Alexander Andryakov
C- Updated 11-OCT-1994 Alexander Andryakov
C-
C----------------------------------------------------------------------
INTEGER IENUMB,IHIST
REAL CODE
IHIST = 100*(IENUMB/100)
CODE = FLOAT(MOD(IENUMB,100))
CALL HFILL(IHIST,CODE,0.,1.)
RETURN
END
C
C*DK MGFLD
C
C=======================================================================
C
SUBROUTINE MGFLD(X,H)
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : Calculate magnetic field inside drift chamber
C- Magnetic field is supposed to be constant
C- H(0,0,Hz)
C-
C- Created 1-JUL-1992 Alexander Andryakov
C- Updated 7-jun-1995 A.Andryakov
C- Magnetic field in the track fit is in units of
C- kGauss*0.00029979 (this was the major change
C- of the present code update)
C- B U T the vertex fit uses magnetic field in units of
C- kGauss only. Current version of the package (p.r.
C- +t.f.+v.f.) is self-consistent in this sence.
C- !!!!! Be very careful when creating the real m.f. map
C- !!!!! routine and keep in mind the above difference.
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:CCONST.INC'
REAL X(3), H(3)
REAL BRADL/0.000/
H(1) = 0.
H(2) = TIFLAV*BRADL/CONVERS
H(3) = TIFLAV*(1.-0.5*BRADL*BRADL)/CONVERS
RETURN
END
C
C*DK DFSAVE
C
C=======================================================================
C
Subroutine DFSAVE(NH,HITS)
C==============================================================
C-
C- Purpose and methods:
C-
C- Stores vital information about partially fitted tracks.
C- It should be called before calling DFITER, so in case
C- something BAD happens a subsequent call to routine DFRESTORE
C- may be done to restore the previous situation.
C- By default it saves the TRACK PARAMETERS,the drift dist.
C- of hit NH and its error (useful when hit NH is going to
C- be flipped), the STEP LENGTH AT FIRST HIT, the error flag NDFERR
C- the guessed chisq DFCHIN, the pointer to the hit contributing
C- most to chisq NDFMCH, the maximum individual chisq
C- DFMCH and multiple scattering parameters DFPMS.
C-
C- Input parameters:
C-
C- Integer NH The current hit, whose Drift Distance is to be
C- saved together with its error
C-
C-
C- Logical HITS If .TRUE. saves also the hits chain NDF and
C- all related info: DFD,DFE,DFS,DFL vectors,
C- NFITPOS and NFITNEG and all multiple scattering
C- related infos: NDFDIM,MDF,NDFMS,DFEMS
C-
C- N.B.
C- DFSAVE DOES NOT DEAL WITH DRIFT DISTANCE IN DHRE2 !!!
C-
C- Created: Feb 1998 F. Ambrosino
C-
C===============================================================
$$IMPLICIT
$$INCLUDE 'K$ITRK:CDFWORK.INC'
$$INCLUDE 'K$ITRK:DFSAV.INC'
INTEGER N,M,NH
LOGICAL HITS
IF(HITS) THEN
NDF1SAVE = NDF1
NDFLSAVE = NDFL
NDFDIMSAVE = NDFDIM
N = NDF1
M = 0
DO WHILE(N.GT.0)
M = M + 1
NDFSAVE(N) = NDF(N)
DFDSAVE(N) = DFD(N)
DFERRSAVE(N) = DFE(N)
DFSSAVE(N) = DFS(N)
DFLSAVE(N) = DFL(N)
MDFSAVE(N) = MDF(N)
N = NDF(N)
ENDDO
NHITSAVE = M
DO M = 2,(NDFDIM-5)
NDFMSSAVE(M) = NDFMS(M)
DFEMSSAVE(M) = DFEMS(M)
ENDDO
NPOSSAVE = NFITPOS
NNEGSAVE = NFITNEG
SAVED_HITS = .TRUE.
ENDIF
NDFERRSAVE = NDFERR
NMCHSAVE = NDFMCH
DFMCHSAVE = DFMCH
DFSAV( 1) = DFKURV
DFSAV( 2) = DFCTGT
DFSAV( 3) = DFSTHE
DFSAV( 4) = DFPHI
DFSAV( 5) = DFX(NDF1)
DFSAV( 6) = DFY(NDF1)
DFSAV( 7) = DFZ(NDF1)
DFSAV( 8) = DFPX(NDF1)
DFSAV( 9) = DFPY(NDF1)
DFSAV(10) = DFPZ(NDF1)
DFSAV(11) = DFCHIN
DFSAV(12) = DFD(NH)
DFSAV(13) = DFE(NH)
DFSAV(14) = DFS(NDF1)
C------- Multiple scattering
DO M = 2,(NDFDIM-5)
DFPMSSAVE(M) = DFPMS(M)
ENDDO
SAVED_PARA = .TRUE.
Return
End
C
C*DK DFRESTORE
C
C=======================================================================
C
Subroutine DFRESTORE(NH,HITS)
C==============================================================
C-
C- Purpose and methods:
C-
C- Restores vital information about partially fitted tracks
C- previously stored by routine DFSAVE.
C- It should be called after a call to DFITER, (and in case
C- something BAD happened) if one wants to restore the
C- previous situation.
C- By default it restores the TRACK PARAMETERS, the STEP LENGTH at 1st hit
C- the drift dist. of hit NH and its error, the guessed chisq DFCHIN,
C- the pointer to the hit contributing most to chisq NDFMCH, and the
C- maximum individual chisq DFMCH, the error flag NDFERR and mult. scatt.
C- parameters DFPMS
C-
C-
C- Input parameters
C-
C- Integer NH The current hit, whose D.D. and err. is to be restored
C-
C- Logical HITS If .TRUE. restore the hits chain NDF and
C- all related info (i.e. vectors DFD,DFE,DFL,DFS)
C- NFITPOS and NFITNEG and all mult. scattering info:
C- NDFDIM,MDF,NDFMS,DFEMS.
C-
C- N.B.
C- DFRESTORE DOES NOT DEAL WITH DRIFT DISTANCE IN DHRE2 !!!
C-
C-
C- Created: Feb 1998 F. Ambrosino
C-
C===============================================================
$$IMPLICIT
$$INCLUDE 'K$ITRK:CDFWORK.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:DFSAV.INC'
INTEGER NH,M,K,K_OLD
LOGICAL HITS
IF(HITS) THEN
IF(.NOT. SAVED_HITS) THEN
CALL ERLOGR('DFRESTORE',ERWARN,0,0,
& 'Tried to restore hits, but they were not saved !!')
GO TO 100
ENDIF
IF(NHITSAVE.NE.(NPOSSAVE+NNEGSAVE)) THEN
CALL ERLOGR('DFRESTORE',ERWARN,0,0,
& 'Inconsistent number of hits from DFSAVE!?!')
GO TO 100
ENDIF
NDF1 = NDF1SAVE
NDFL = NDFLSAVE
NDFDIM = NDFDIMSAVE
K = NDF1
DO WHILE(K.GT.0)
NDF(K) = NDFSAVE(K)
DFD(K) = DFDSAVE(K)
DFE(K) = DFERRSAVE(K)
DFS(K) = DFSSAVE(K)
DFL(K) = DFLSAVE(K)
MDF(K) = MDFSAVE(K)
K = NDFSAVE(K)
ENDDO
NFITPOS = NPOSSAVE
NFITNEG = NNEGSAVE
DO M = 2,(NDFDIM-5)
NDFMS(M) = NDFMSSAVE(M)
DFEMS(M) = DFEMSSAVE(M)
ENDDO
ENDIF
100 Continue
IF(.NOT. SAVED_PARA) THEN
CALL ERLOGR('DFRESTORE',ERWARN,0,0,
& 'Tried to restore params, but were not saved !!')
GO TO 200
ENDIF
NDFERR = NDFERRSAVE
NDFMCH = NMCHSAVE
DFMCH = DFMCHSAVE
DFKURV = DFSAV( 1)
DFCTGT = DFSAV( 2)
DFSTHE = DFSAV( 3)
DFPHI = DFSAV( 4)
DFX(NDF1) = DFSAV( 5)
DFY(NDF1) = DFSAV( 6)
DFZ(NDF1) = DFSAV( 7)
DFPX(NDF1) = DFSAV( 8)
DFPY(NDF1) = DFSAV( 9)
DFPZ(NDF1) = DFSAV(10)
DFCHIN = DFSAV(11)
DFD(NH) = DFSAV(12)
DFE(NH) = DFSAV(13)
DFS(NDF1) = DFSAV(14)
C------ Multiple scattering
DO M = 2,(NDFDIM-5)
DFPMS(M) = DFPMSSAVE(M)
ENDDO
200 Continue
Return
End
C
C*DK DTFSKILL
C
C====================================================================
SUBROUTINE DTFSKILL(TRKMERGE)
C====================================================================
C
C DESCRIPTION:
C ===========
C BY comparing the chi2 of the new trk with the mean of old ones it kill
c the worse(s). WARNING:not always the trk is fitted so DTFSKILL
c can return without any action
C
$$implicit
c------ Include Files ----
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$include 'K$INC:OFERCO.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$include 'K$ITRK:DTFS.INC'
$$include 'K$ITRK:STUMPCOM.INC'
c------ External Declaration -----
integer BLOCAT,MDROP
c------ Local Declaration --------
integer STATUS,ISTAT,TRKMERGE
integer IND_DTM,INDDAT_DTM,IC,BESTICOP
integer IND_DTF1,IND_DTF2
integer TRK1,TRK2,INDDAT_DTF1,INDDAT_DTF2
integer POINT_TRM,POINT_TR1,POINT_TR2
integer NPOS_TRM,NNEG_TRM,NHIT_TRM
integer NPOS_TR1,NNEG_TR1,NHIT_TR1
integer NPOS_TR2,NNEG_TR2,NHIT_TR2
real CHI2_TRM,CHI2_TR1,CHI2_TR2,CHI2_MEAN
logical TRM_KILL,TRM_KILL0,HIT_CUT
real PTTRI,COTTRI,PHITRI,PTTRJ,COTTRJ,PHITRJ
real PXTRI,PYTRI,PZTRI,PXTRJ,PYTRJ,PZTRJ
real ENI,ENJ,PIJ2,KLMASS,MPI
data MPI/139.570/
C------ Start of code ------
STATUS = OFSUCC
c------ Firstly check if the DTFS TRK existed
c------ (or existed only the DPRS)
ISTAT = BLOCAT (IW,'DTFS',TRKMERGE,IND_DTM,INDDAT_DTM)
if (ISTAT.ne.YESUCC) then !this could happen
goto 900
endif
do IC=1,NCOP
if(TRKMERGE.eq.BNKS(IC)) then
BESTICOP =IC
endif
enddo
TRK1=FIRSTUMP(BESTICOP)
TRK2=SECSTUMP(BESTICOP)
ISTAT = BLOCAT (IW,'DTFS',TRK1,IND_DTF1,INDDAT_DTF1)
if (ISTAT.ne.YESUCC) then
call ERLOGR ('DTFSKILL',ERWARN,0,0,'Cannot find DTFS')
goto 990
endif
ISTAT = BLOCAT (IW,'DTFS',TRK2,IND_DTF2,INDDAT_DTF2)
if (ISTAT.ne.YESUCC) then
call ERLOGR ('DTFSKILL',ERWARN,0,0,'Cannot find DTFS')
goto 990
endif
POINT_TRM=INDDAT_DTM+IW(INDDAT_DTM+DTFHDS)
POINT_TR1=INDDAT_DTF1+IW(INDDAT_DTF1+DTFHDS)
POINT_TR2=INDDAT_DTF2+IW(INDDAT_DTF2+DTFHDS)
C---- Compare the CHI2-------------------------------
C First Merged TRK
NPOS_TRM = IW(POINT_TRM+DTFNHF)/10000
NNEG_TRM = MOD(IW(POINT_TRM+DTFNHF),10000)
NHIT_TRM = NPOS_TRM + NNEG_TRM
CHI2_TRM = RW(POINT_TRM+DTFCHI)/(NHIT_TRM-5)
C Then First TRK
NPOS_TR1 = IW(POINT_TR1+DTFNHF)/10000
NNEG_TR1 = MOD(IW(POINT_TR1+DTFNHF),10000)
NHIT_TR1 = NPOS_TR1 + NNEG_TR1
CHI2_TR1 = RW(POINT_TR1+DTFCHI)/(NHIT_TR1-5)
C Then Second TRK
NPOS_TR2 = IW(POINT_TR2+DTFNHF)/10000
NNEG_TR2 = MOD(IW(POINT_TR2+DTFNHF),10000)
NHIT_TR2 = NPOS_TR2 + NNEG_TR2
CHI2_TR2 = RW(POINT_TR2+DTFCHI)/(NHIT_TR2-5)
CHI2_MEAN=0.5*(CHI2_TR1+CHI2_TR2)
if(chi2_mean.lt.0.000001) chi2_mean = 0.000001
TRM_KILL0=abs((CHI2_TRM-CHI2_MEAN)/CHI2_MEAN).gt.0.6
HIT_CUT = (NHIT_TR1+ NHIT_TR2 - NHIT_TRM) .gt. 20
TRM_KILL = TRM_KILL0 .or. HIT_CUT
if(TRM_KILL) then ! kill new merged Trk
STATUS = MDROP(IW,'DTFS',TRKMERGE)
if (IW(2).ne.YESUCC) then
call ERLOGR ('DTFSKILL',ERWARN,0,0,'Cannot drop DTFS')
goto 990
endif
else
C Store also the #dtfs of joined tracks
IW(POINT_TRM+DTFVMT) = TRK1
IW(POINT_TRM+DTFVDT) = TRK2
RW(POINT_TRM+DTFCHM)=CHI2_MEAN
C ALSO for Monte Carlo store the number of hits of the two tracks
C
IW(POINT_TRM+DTFCEN) = NHIT_TR1+NHIT_TR2
STATUS = MDROP(IW,'DTFS',TRK1)
if (IW(2).ne.YESUCC) then
call ERLOGR ('DTFSKILL',ERWARN,0,0,'Cannot drop DTFS')
goto 990
endif
STATUS = MDROP(IW,'DTFS',TRK2)
if (IW(2).ne.YESUCC) then
call ERLOGR ('DTFSKILL',ERWARN,0,0,'Cannot drop DTFS')
goto 990
endif
endif
900 continue
goto 999
990 continue
STATUS=OFFAIL
999 return
end
C
C*DK HITCODE
C
C=======================================================================
C
SUBROUTINE HITCODE()
$$IMPLICIT
C Global Declarations:-
C =====================
$$INCLUDE 'C$INC:JOBSTA.INC'
$$INCLUDE 'C$INC:BPOYBS.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:BNKTYP.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
INTEGER BLOCAT
INTEGER STATUS,INDR,INDDATR,DHRENCO,IND,
+ IDHI,IDDATHI,DHIENCO,DHIENRO,
+ IDKI,INDKATH,
+ HITG,IPOINT,IPTYPE,IPCODE(50),
+ PACODE,JCODE,IWIRE,N,JJ
INTEGER IFLAGCHI
INTEGER MUP,MUM,PIONP,PIONM
COMMON/CHIINIT/IPTYPE(MXHIT),PACODE
COMMON/CHIINIT2/MUP,MUM,PIONP,PIONM
COMMON/IFCHI/IFLAGCHI
IFLAGCHI=0
CALL VZERO(IPTYPE,MXHIT)
CALL VZERO(IPCODE,50)
* dhit bank (IDHI,IDDATHI)
STATUS=BLOCAT(IW,'DHIT',1,IDHI,IDDATHI)
IF (STATUS.NE.YESUCC) RETURN
IDHI=IDDATHI+IW(IDDATHI+DHIHDS)
DHIENCO=IW(IDDATHI+DHINCO)
DHIENRO=IW(IDDATHI+DHINRO)
N = NDF1
DO WHILE (N.GT.0)
IF (LDFL(N).GT.0) THEN
IWIRE=IABS(IW(JDPRS1+DPRLWI*(N-1)+DPRIWI))
IF (IWIRE.GT.10000) THEN
IFLAGCHI=1
WRITE(6,*)'strange value of iwire'
WRITE(6,*)'iwire=',iwire
RETURN
ENDIF
IPTYPE(N)=IW(IDHI+(IWIRE-1)*DHIENCO+DHIPTY)
IPCODE(IPTYPE(N))=IPCODE(IPTYPE(N))+1
ENDIF
N = NDF(N)
ENDDO
JCODE=0
PACODE=0
DO JJ=1,50
IF (IPCODE(JJ).GE.JCODE) THEN
JCODE=IPCODE(JJ)
PACODE=JJ
ENDIF
ENDDO
MUP=0
MUM=0
PIONP=0
PIONM=0
DO JJ=1,50
IF (IPCODE(JJ).GT.0) THEN
IF (JJ.EQ.5) THEN
MUP=IPCODE(JJ)
ENDIF
IF (JJ.EQ.6) THEN
MUM=IPCODE(JJ)
ENDIF
IF (JJ.EQ.8) THEN
PIONP=IPCODE(JJ)
ENDIF
IF (JJ.EQ.9) THEN
PIONM=IPCODE(JJ)
ENDIF
ENDIF
ENDDO
END
C==========================================================================
SUBROUTINE DFTWOPR(KINK_POS,NTRACK)
C========================================================================
C
C INPUT PARAMETERS:
C
C Int KINK_POS ! Guessed position of kink in the track
C
C
C Int NTRACK ! Number of DPRS bank of the track
C
C
C OUTPUT PARAMETERS: None
C
C
C======================================================================
C=
C= Purpose and Methods : This routine splits one track, with kink,
C= (as signaled by DFKINK) in two separate tracks,
C= then writes two new DPRS bank with version number
C= BNKVRS = 10 + NTRACK (so to keep track of their
C= origin)
C=
C= Created 15-4-1998 F.Ambrosino C. Di Donato
C=========================================================================
C
C
$$IMPLICIT
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:BPOYBS.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$INC:OFERCO.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:MCKINE.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
C
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
C
$$INCLUDE 'C$INC:jobsta.INC'
C
C
C
C
C---------------- COMMON SPLITTED TRACK PARAMETERS ------------------------
C
INTEGER MAX_PAR
PARAMETER (MAX_PAR=9)
INTEGER MAX_HITS
PARAMETER (MAX_HITS=250)
C
REAL PAR1
INTEGER NPOS_1, NNEG_1, NHITS1
REAL PAR2
INTEGER NPOS_2, NNEG_2, NHITS2
REAL STEPS
C
COMMON /PAR/ PAR1(MAX_PAR), NPOS_1, NNEG_1, NHITS1,
+ PAR2(MAX_PAR), NPOS_2, NNEG_2, NHITS2,
+ STEPS(MAX_HITS)
C-----------------------------------------------------------------------
C
C
C------------------------------------------External functions
INTEGER BTYMKG, BLOCAT, BIGEST
C
C------------------------------------------Local Declarations
INTEGER STATUS
INTEGER KINK_POS, NTRACK
INTEGER IND_DHRE, INDAT_DHRE
INTEGER IND_DPRS, INDAT_DPRS
INTEGER NUM_ROW, NUM_DATA_ROW
INTEGER MAXTRACK, NHITSP
INTEGER NEW1, NEW2
INTEGER IND1, INDAT1, IND2, INDAT2
INTEGER B1NKSIZ, B2NKSIZ
INTEGER DHRE_ADDR, LINK
INTEGER IOFFSET, IOFFSET2
C
INTEGER I, II
C
CHARACTER*36 B1NKTYP
DATA B1NKTYP/'2I4,9R4,2I4,R4,2I4,2R4,I4,XYX(I4,R4)'/
C
CHARACTER*36 B2NKTYP
DATA B2NKTYP/'2I4,9R4,2I4,R4,2I4,2R4,I4,XXY(I4,R4)'/
C
C
C
C
C--------------------------------Executable code ---------------------
C
C Inizialization
C
I = 0
II = 0
NUM_ROW = 0
NUM_DATA_ROW = 0
NEW1 = 0
NEW2 = 0
C
C ====== Find the max bank number of DPRS chain
C
STATUS = BIGEST(IW,'DPRS',MAXTRACK)
C
NEW1 = MAXTRACK + 1
NEW2 = MAXTRACK + 2
C
C
C ======= Create the new DPRS bank for the daughter tracks
C
C
WRITE (B1NKTYP(27:29),'(I3)') NHITS1
WRITE (B2NKTYP(27:29),'(I3)') NHITS2
C
STATUS = BTYMKG(IW,'DPRS',NEW1,
+ B1NKTYP,B1NKSIZ,IND1,INDAT1)
IF (STATUS.NE.YESUCC) THEN
CALL ERLOGR('DSPLIT',ERWARN,0,STATUS,
+ 'Error creating new bank DPRS #1')
GOTO 999
ENDIF
C
STATUS = BTYMKG(IW,'DPRS',NEW2,
+ B2NKTYP,B2NKSIZ,IND2,INDAT2)
IF (STATUS.NE.YESUCC) THEN
CALL ERLOGR('DSPLIT',ERWARN,0,STATUS,
+ 'Error creating new bank DPRS #2')
GOTO 999
ENDIF
C
C
C
C===== Locate DHRE 2 bank
C
STATUS = BLOCAT(IW,'DHRE',2,IND_DHRE
& ,INDAT_DHRE)
IF(STATUS.NE.YESUCC) go to 999
NUM_ROW = IW(INDAT_DHRE +DHRNRO)
NUM_DATA_ROW = IW(INDAT_DHRE +DHRNCO)
INDAT_DHRE = INDAT_DHRE + IW(INDAT_DHRE+DHRHDS)
C
C======== Locate the DPRS bank to split
C
STATUS = BLOCAT(IW,'DPRS',NTRACK,IND_DPRS
& ,INDAT_DPRS)
IF(STATUS.NE.YESUCC) go to 999
INDAT_DPRS = INDAT_DPRS + IW(INDAT_DPRS + DPRHDS)
IND_DPRS = INDAT_DPRS + DPRCEN
NHITSP = IW(INDAT_DPRS+DPRPOS) + IW(INDAT_DPRS+DPRNEG)
C
C=================================================================
C Fill new DPRS banks => from DFKINK parameters vectors
C and related info
C
C NEW DPRS #1----------------------------------------------------
C
IW(INDAT1+DPRHDS) = 2
IW(INDAT1+DPRVRN) = 10+NTRACK
INDAT1 = INDAT1 + IW(INDAT1+DPRHDS)
RW(INDAT1+DPRXCO) = PAR1(1)
RW(INDAT1+DPRYCO) = PAR1(2)
RW(INDAT1+DPRZCO) = PAR1(3)
RW(INDAT1+DPRXLS) = PAR1(4)
RW(INDAT1+DPRYLS) = PAR1(5)
RW(INDAT1+DPRZLS) = PAR1(6)
RW(INDAT1+DPRCUR) = PAR1(7)
RW(INDAT1+DPRPHI) = PAR1(8)
RW(INDAT1+DPRCTT) = PAR1(9)
RW(INDAT1+DPRQUA) = 0.
IW(INDAT1+DPRPOK) = 3
IW(INDAT1+DPRPOS) = NPOS_1
IW(INDAT1+DPRNEG) = NNEG_1
IW(INDAT1+DPRLRG) = 0
RW(INDAT1+DPRPHF) = 0.
RW(INDAT1+DPRPHL) = 0.
C
C
IW(INDAT1+DPRCEN) = IW(INDAT_DPRS+DPRCEN)
IND1 = INDAT1 + DPRCEN
C
C Fills Variable length part of the Bank
C
DO I = 1,NHITS1
IOFFSET = (I-1)*DPRLWI
DHRE_ADDR = IABS(IW(IND_DPRS+IOFFSET+DPRIWI))
IW(IND1+IOFFSET+DPRIWI) = DHRE_ADDR
RW(IND1+IOFFSET+DPRWTR) = STEPS(I)
C
C Fills DHRE2 link database with the new DPRS bank number
C DHRE_ADDR is number of the hit in DHRE2 bank
C
IOFFSET = (DHRE_ADDR-1)*NUM_DATA_ROW
IW(INDAT_DHRE+IOFFSET+DHRTRN) = NEW1
ENDDO
C
C
C DPRS BANK #2--------------------------------------------------------
C
IW(INDAT2+DPRHDS) = 2
IW(INDAT2+DPRVRN) = 10+NTRACK
INDAT2 = INDAT2 + IW(INDAT2+DPRHDS)
RW(INDAT2+DPRXCO) = PAR2(1)
RW(INDAT2+DPRYCO) = PAR2(2)
RW(INDAT2+DPRZCO) = PAR2(3)
RW(INDAT2+DPRXLS) = PAR2(4)
RW(INDAT2+DPRYLS) = PAR2(5)
RW(INDAT2+DPRZLS) = PAR2(6)
RW(INDAT2+DPRCUR) = PAR2(7)
RW(INDAT2+DPRPHI) = PAR2(8)
RW(INDAT2+DPRCTT) = PAR2(9)
RW(INDAT2+DPRQUA) = 0.
IW(INDAT2+DPRPOK) = 3
IW(INDAT2+DPRPOS) = NPOS_2
IW(INDAT2+DPRNEG) = NNEG_2
IW(INDAT2+DPRLRG) = 0
RW(INDAT2+DPRPHF) = 0.
RW(INDAT2+DPRPHL) = 0.
C
C
IW(INDAT2+DPRCEN) = IW(INDAT_DPRS+DPRCEN)
IND2 = INDAT2 + DPRCEN
C
C Fills Variable length part of the Bank
C
DO I = 1,NHITS2
IOFFSET = (I+NHITS1-1)*DPRLWI
IOFFSET2 = (I-1)*DPRLWI
DHRE_ADDR = IABS(IW(IND_DPRS+IOFFSET+DPRIWI))
IW(IND2+IOFFSET2+DPRIWI) = DHRE_ADDR
RW(IND2+IOFFSET2+DPRWTR) = STEPS(I+NHITS1)
C
C Fills DHRE2 link database with the new DPRS bank number
C DHRE_ADDR is number of the hit in DHRE2 bank
C
IOFFSET = (DHRE_ADDR-1)*NUM_DATA_ROW
IW(INDAT_DHRE+IOFFSET+DHRTRN) = NEW2
ENDDO
C
C==================================================================
C
C Purge track candidate number-link (DHRTRN)
C for the ''old'' track (useful to set free the hits added to
C the old track and give the new tracks a chance to get them)
C
DO I = 1, NUM_ROW
IOFFSET = (I-1)*NUM_DATA_ROW
LINK = IW(INDAT_DHRE + IOFFSET + DHRTRN)
IF(LINK.EQ.NTRACK) THEN
IW(INDAT_DHRE + IOFFSET + DHRTRN) = 0
ENDIF
ENDDO
C
C-------------------------------------------------------------------
C
999 CONTINUE
RETURN
END
C
C========================================================================
SUBROUTINE DFSPLIT
C========================================================================
$$IMPLICIT
C------------------------------------------------------------------------
C-
C- Purpose and Methods : This routine is intended to perform the final
C- check wether a track is to be splitted or not
C- The real decision, however is taken in routine
C- DFCOMPARE, and DFSPLIT makes all needed operations
C- to clean DPRS and DTFS bank chains
C-
C- Author: F. Ambrosino C. Di Donato
C =======
C-
C- created 17 4 98
C- updated 16 9 98
C----------------------------------------------------------------------
C
C Global Declarations:-
C =====================
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:DTFS.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
C
C
C External Functions
C ==================
C
INTEGER BLOCAT, BNUMB, BDATA, BNEXT, MDROP
C
C Local Declarations:-
C ====================
C
C
INTEGER STATUS
INTEGER IND_DHRE, INDAT_DHRE, DHRE_DATA
INTEGER STATUS_DTFS, IND_DTFS, INDAT_DTFS
INTEGER STATUS_DPRS, IND_DPRS, INDAT_DPRS
INTEGER VERSION, IMOTHER, IDAUGHT1, IDAUGHT2
C
INTEGER IND_MOT, INDAT_MOT, NHITS
INTEGER IOFFSET, DHRE_ADDR, DHRE_HITS, LINK
C
LOGICAL KEEPMOTHER
LOGICAL KILL_D1,KILL_D2,CLEAN
C
C
INTEGER I
C
C
IF(.NOT. KINK_RECO) goto 1000
C
C Look for kink-track and splitted-track
C If bank Version Number in DPRS bank is >10,
C this is a stump of a splitted track; the mother
C bank number is (bank Version number - 10)
C
C
IMOTHER = 0
IDAUGHT1 = 0
IDAUGHT2 = 0
Keepmother = .false.
KILL_D1 = .FALSE.
KILL_D2 = .FALSE.
CLEAN = .FALSE.
C
C
C
STATUS = BLOCAT(IW,'DHRE',2,IND_DHRE,INDAT_DHRE)
IF(STATUS.NE.YESUCC) THEN
Call ERLOGR('DFSPLIT',ERWARN,0,0,
& 'DHRE bank NOT FOUND!!!')
GO TO 1000
ENDIF
C
DHRE_DATA = IW(INDAT_DHRE + DHRNCO)
DHRE_HITS = IW(INDAT_DHRE + DHRNRO)
IND_DHRE = INDAT_DHRE + IW(INDAT_DHRE + DHRHDS)
C
STATUS_DPRS = BLOCAT(IW,'DPRS',-1,IND_DPRS,INDAT_DPRS)
DO WHILE (IND_DPRS.GT.0)
VERSION = IW(INDAT_DPRS + DPRVRN)
IF (VERSION.lt.10) then
STATUS = BNEXT(IW,IND_DPRS,IND_DPRS)
go to 999
ENDIF
IMOTHER = VERSION - 10
C
C--------------------------------This bank SHOULD exist, but in case the
C mother is itself daughter of another track
C it may have been erased : in such cases
C we just drop all the 'discendence'
C
c
STATUS = BLOCAT(IW,'DPRS',IMOTHER,IND_MOT,INDAT_MOT)
IF(STATUS.NE.YESUCC) THEN
KEEPMOTHER = .TRUE.
STATUS = BNUMB(IW,IND_DPRS,IDAUGHT1)
STATUS_DPRS=BNEXT(IW,IND_DPRS,IND_DPRS)
STATUS = BNUMB(IW,IND_DPRS,IDAUGHT2)
GO TO 777
ENDIF
C
C--------------------------IDAUGHT1 is the first daughter bank number
C
STATUS = BNUMB(IW,IND_DPRS,IDAUGHT1)
C
C-------------------------This is the second daughter, this bank MUST exist
C
STATUS_DPRS=BNEXT(IW,IND_DPRS,IND_DPRS)
IF(IND_DPRS.EQ.0) THEN
Call Erlogr('DFSPLIT',erwarn,0,0,
& 'Second daughter not found !!!')
keepmother = .true.
Go to 777
ENDIF
STATUS_DPRS = BDATA(IW,IND_DPRS,INDAT_DPRS)
IF(IW(INDAT_DPRS+DPRVRN).NE.VERSION) THEN
Call Erlogr('DFSPLIT',erwarn,0,0,
& 'Second daughter wrong link !!!')
keepmother = .true.
Go to 777
ENDIF
C--------------------------IDAUGHT2 is the second daughter bank number
C
STATUS = BNUMB(IW,IND_DPRS,IDAUGHT2)
C-------------------------------------------------------
C
C Now we have links to mother track and to both daughters and we call
C DFCOMPARE to decide wether the mother bank must be kept (and the
c daughters killed) or vice-versa
C--------------------------------------------------------
Call DFCOMPARE(IMOTHER,IDAUGHT1,IDAUGHT2,KEEPMOTHER)
C
777 CONTINUE
CC
C This MUST be put here to be safe, to locate next bank the previous
C must exist, and maybe we are going to kill last DPRS bank located....
C
STATUS = BNEXT(IW,IND_DPRS,IND_DPRS)
C
C
IF(KEEPMOTHER) THEN
C
C-----------------------------------Kill daughters DPRS banks...
C
STATUS = MDROP(IW,'DPRS',IDAUGHT1)
STATUS = MDROP(IW,'DPRS',IDAUGHT2)
KILL_D1 = .TRUE.
KILL_D2 = .TRUE.
C
C------------------------------------..and then corresponding DTFS banks
C if they exist....
STATUS_DTFS = BLOCAT(IW,'DTFS',IDAUGHT1
& ,IND_DTFS,INDAT_DTFS)
IF(STATUS_DTFS.EQ.YESUCC) THEN
STATUS = MDROP(IW,'DTFS',IDAUGHT1)
ENDIF
C
STATUS_DTFS = BLOCAT(IW,'DTFS',IDAUGHT2
& ,IND_DTFS,INDAT_DTFS)
IF(STATUS_DTFS.EQ.YESUCC) THEN
STATUS = MDROP(IW,'DTFS',IDAUGHT2)
ENDIF
C------------------------------------...finally restore DHRE2 original links
C
STATUS = BLOCAT(IW,'DPRS',IMOTHER,IND_MOT,INDAT_MOT)
IF(STATUS.ne.YESUCC) THEN
Call erlogr('DFSPLIT',erwarn,0,0,
& 'DPRS mother bank not found ?!?!')
goto 999
ENDIF
IND_MOT = INDAT_MOT + IW(INDAT_MOT + DPRHDS)
NHITS = IW(IND_MOT+DPRPOS)+IW(IND_MOT+DPRNEG)
IND_MOT = IND_MOT + DPRCEN
DO I = 1,NHITS
IOFFSET = (I-1)*DPRLWI
DHRE_ADDR = IABS(IW(IND_MOT+IOFFSET+DPRIWI))
IOFFSET = (DHRE_ADDR-1)*DHRE_DATA
IW(IND_DHRE+IOFFSET+DHRTRN) = IMOTHER
ENDDO
ELSE
C------------------------------Kill mother banks
C
STATUS = MDROP(IW,'DPRS',IMOTHER)
STATUS_DTFS = BLOCAT(IW,'DTFS',IMOTHER
& ,IND_DTFS,INDAT_DTFS)
IF(STATUS_DTFS.EQ.YESUCC) THEN
STATUS = MDROP(IW,'DTFS',IMOTHER)
ENDIF
C
C--------------------------..and also the daughters if have not
C been fitted !!!
C
C
STATUS_DTFS = BLOCAT(IW,'DTFS',IDAUGHT1
& ,IND_DTFS,INDAT_DTFS)
IF(STATUS_DTFS.NE.YESUCC) THEN
STATUS = MDROP(IW,'DPRS',IDAUGHT1)
KILL_D1 = .true.
ENDIF
C
STATUS_DTFS = BLOCAT(IW,'DTFS',IDAUGHT2
& ,IND_DTFS,INDAT_DTFS)
IF(STATUS_DTFS.NE.YESUCC) THEN
STATUS = MDROP(IW,'DPRS',IDAUGHT2)
KILL_D2 = .TRUE.
C IF(KILL_D1) then
C Call erlogr('DFSPLIT',erwarn,0,0,
C & 'Both daughters killed')
C ENDIF
ENDIF
ENDIF
C
C-------------------------FINALLY CLEAN DHRE2 LINKS TO KILLED DAUGHTERS
C
IF(KILL_D1.OR.KILL_D2) THEN
DO I = 1,DHRE_HITS
IOFFSET = (I-1)*DHRE_DATA
LINK = IW(IND_DHRE + IOFFSET + DHRTRN)
CLEAN = (KILL_D1.AND.(LINK.EQ.IDAUGHT1)).OR.
& (KILL_D2.AND.(LINK.EQ.IDAUGHT2))
IF(CLEAN) THEN
IW(IND_DHRE + IOFFSET + DHRTRN) = 0
ENDIF
ENDDO
ENDIF
C
C
999 Continue
STATUS_DPRS=BDATA(IW,IND_DPRS,INDAT_DPRS)
C
ENDDO
C
C
C
1000 CONTINUE
RETURN
END
C
C===================================================================
SUBROUTINE DFCOMPARE(IMOTHER,IDAUGHT1,IDAUGHT2,KEEPMOTHER)
C===================================================================
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Purpose and Methods : This routine is intended to take the final decision
C- wether a track which has been splitted in two
C- has to be kept (and its daughters killed) or
C- vice-versa.
C-
C- Input parameters :
C-
C- Int IMOTHER # of Mother track's DPRS bank
C- Int Idaught1 # of 1st daughter track's DPRS bank
C- Int Idaught2 # of 2nd daughter track's DPRS bank
C-
C- Output parameters :
C-
C- Logical keepmother If .true. keep mother track
C-
C-
C- Author: F. Ambrosino C. Di Donato
C =======
C-
C- created 17 4 98
C- updated 16 9 98
C- updated 4 12 98 (New cuts definition) F. A.
C----------------------------------------------------------------------
C
C Global Declarations:-
C =====================
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$ITRK:DTFS.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
C
C External Functions:-
C ===================
C
INTEGER BLOCAT
C
C Local Declarations:-
C ====================
C
INTEGER IMOTHER
& ,idaught1
& ,idaught2
C
LOGICAL KEEPMOTHER
C
INTEGER STATUS
INTEGER IND_MOTH, INDAT_MOTH
INTEGER IND_DTFS1, INDAT_DTFS1
INTEGER IND_DTFS2, INDAT_DTFS2
C
REAL CHIMOTHER, CHIDAUGHT1, CHIDAUGHT2,BESTCHI
REAL CHIM1, CHIM2
REAL TF_CURV, TF_CTGT, TF_PHI
REAL PX1,PY1,PZ1,PX2,PY2,PZ2,P1,P2,DOT,PT1,PT2
REAL COSTHETA,DELTAPT,DELTAPZ
REAL PXMOT,PYMOT,PZMOT
REAL FRAC,FRAC2
INTEGER MAXHITS,MINHITS
INTEGER CODEH,CODEH1,CODEH2
INTEGER NHITPOS, NHITNEG
INTEGER NHITPOS1, NHITNEG1
INTEGER NHITPOS2, NHITNEG2
INTEGER NHITMOT,NHIT1,NHIT2
C
LOGICAL FITD1,FITD2,FITMOT
LOGICAL ONLY_1_FITD,ONLY_2_FITD,BOTH_FITD,NO_FITD
LOGICAL CUT1,CUT2,CUT3,CUT4,CUT5,CUT6
C
C----------------------------------------CUTS
INTEGER DIFFHIT
REAL MINCHI1,MINCHIDIFF1
PARAMETER(DIFFHIT=5)
PARAMETER(MINCHI1=9.)
PARAMETER(MINCHIDIFF1=6.)
REAL MINCOSTHETA,MAXDPT,MAXDPZ,MAXFRAC,MINCHIFRAC
PARAMETER(MINCOSTHETA=0.9975)
PARAMETER(MAXDPT=5.)
PARAMETER(MAXDPZ=5.)
PARAMETER(MAXFRAC=-0.26)
PARAMETER(MINCHIFRAC=0.16)
C
C--------------------------------------INITIALIZATIONS
FITD1 = .FALSE.
FITD2 = .FALSE.
FITMOT = .FALSE.
ONLY_1_FITD = .FALSE.
ONLY_2_FITD = .FALSE.
BOTH_FITD = .FALSE.
NO_FITD = .FALSE.
C
C
C
C---------------------------------LOCATE BANKS (IF THEY EXIST!)
C
STATUS = BLOCAT(IW,'DTFS',IMOTHER,IND_MOTH,INDAT_MOTH)
IF(STATUS.EQ.YESUCC) THEN
FITMOT = .TRUE.
INDAT_MOTH = INDAT_MOTH + IW(INDAT_MOTH+DTFHDS)
CHIMOTHER = RW(INDAT_MOTH+DTFCHI)
TF_CURV = RW(INDAT_MOTH+DTFCUR)
TF_PHI = RW(INDAT_MOTH+DTFPHI)
TF_CTGT = RW(INDAT_MOTH+DTFCTT)
PXMOT = (1000./ABS(TF_CURV))
& *COS(TF_PHI)
PYMOT = (1000./ABS(TF_CURV))
& *SIN(TF_PHI)
PZMOT = (1000./ABS(TF_CURV))
& *TF_CTGT
CODEH = IW(INDAT_MOTH+DTFNHF)
NHITPOS = (CODEH/10000)
NHITNEG = (CODEH-NHITPOS*10000)
NHITMOT = NHITPOS + NHITNEG
IF(NHITMOT.GT.5) THEN
CHIMOTHER = CHIMOTHER/FLOAT(NHITMOT - 5)
ENDIF
ENDIF
C
STATUS = BLOCAT(IW,'DTFS',IDAUGHT1,IND_DTFS1,INDAT_DTFS1)
IF(STATUS.EQ.YESUCC) THEN
FITD1 = .TRUE.
INDAT_DTFS1 = INDAT_DTFS1 + IW(INDAT_DTFS1+DTFHDS)
CHIDAUGHT1 = RW(INDAT_DTFS1+DTFCHI)
TF_CURV = RW(INDAT_DTFS1+DTFCUR)
TF_PHI = RW(INDAT_DTFS1+DTFPHI)
TF_CTGT = RW(INDAT_DTFS1+DTFCTT)
PX1 = (1000./ABS(TF_CURV))
& *COS(TF_PHI)
PY1 = (1000./ABS(TF_CURV))
& *SIN(TF_PHI)
PZ1 = (1000./ABS(TF_CURV))
& *TF_CTGT
CODEH1 = IW(INDAT_DTFS1+DTFNHF)
NHITPOS1 = (CODEH1/10000)
NHITNEG1 = (CODEH1-NHITPOS1*10000)
NHIT1 = NHITPOS1 + NHITNEG1
IF(NHIT1.GT.5) THEN
CHIDAUGHT1 = CHIDAUGHT1/FLOAT(NHIT1 - 5)
ENDIF
ENDIF
C
STATUS = BLOCAT(IW,'DTFS',IDAUGHT2,IND_DTFS2,INDAT_DTFS2)
IF(STATUS.EQ.YESUCC) THEN
FITD2 = .TRUE.
INDAT_DTFS2 = INDAT_DTFS2 + IW(INDAT_DTFS2+DTFHDS)
CHIDAUGHT2 = RW(INDAT_DTFS2+DTFCHI)
TF_CURV = RW(INDAT_DTFS2+DTFCUR)
TF_PHI = RW(INDAT_DTFS2+DTFPHI)
TF_CTGT = RW(INDAT_DTFS2+DTFCTT)
PX2 = (1000./ABS(TF_CURV))
& *COS(TF_PHI)
PY2 = (1000./ABS(TF_CURV))
& *SIN(TF_PHI)
PZ2 = (1000./ABS(TF_CURV))
& *TF_CTGT
CODEH2 = IW(INDAT_DTFS2+DTFNHF)
NHITPOS2 = (CODEH2/10000)
NHITNEG2 = (CODEH2-NHITPOS2*10000)
NHIT2 = NHITPOS2 + NHITNEG2
IF(NHIT2.GT.5) THEN
CHIDAUGHT2 = CHIDAUGHT2/FLOAT(NHIT2 - 5)
ENDIF
ENDIF
C
ONLY_1_FITD = (FITD1.AND..NOT.FITD2)
ONLY_2_FITD = (FITD2.AND..NOT.FITD1)
C
BOTH_FITD = FITD1.AND.FITD2
NO_FITD = .NOT.FITD1.AND..NOT.FITD2
C
C================================================NOW COMPARING TRACKS
C
IF(.NOT.FITMOT) THEN
KEEPMOTHER = .FALSE.
GO TO 999
ENDIF
C
IF(NO_FITD) THEN
KEEPMOTHER = .TRUE.
GO TO 999
ENDIF
C--------------------------------HERE THE CHOICE WAS VERY SIMPLE...
C
C
C--------------------------------....BUT NOW WE HAVE REALLY TO CHOOSE!
IF(ONLY_1_FITD) THEN
CUT1 = (NHITMOT - NHIT1).GT.DIFFHIT
CUT2 = CHIDAUGHT1.GT.MINCHI1
CUT3 = (CHIDAUGHT1 - CHIMOTHER).GT.MINCHIDIFF1
KEEPMOTHER = CUT1.OR.(CUT2.AND.CUT3)
ELSEIF(ONLY_2_FITD) THEN
CUT1 = (NHITMOT - NHIT2).GT.DIFFHIT
CUT2 = CHIDAUGHT2.GT.MINCHI1
CUT3 = (CHIDAUGHT2 - CHIMOTHER).GT.MINCHIDIFF1
KEEPMOTHER = CUT1.OR.(CUT2.AND.CUT3)
ELSEIF(BOTH_FITD) THEN
C
p1 = sqrt((px1**2)+(py1**2)+(pz1**2))
p2 = sqrt((px2**2)+(py2**2)+(pz2**2))
dot = (PX1*PX2)+(PY1*PY2)+(PZ1*PZ2)
COSTHETA = dot/(p1*p2)
pt1 = sqrt(px1**2+py1**2)
pt2 = sqrt(px2**2+py2**2)
DELTAPT = pt1 - pt2
DELTAPZ = pz1 - pz2
MAXHITS = MAX(NHIT1,NHIT2)
MINHITS = MIN(NHIT1,NHIT2)
BESTCHI = MIN(CHIDAUGHT1,CHIDAUGHT2)
FRAC = float(MAXHITS - NHITMOT)/float(NHITMOT)
FRAC2 = (BESTCHI - CHIMOTHER)/CHIMOTHER
C
CUT1 = FRAC2.GT.MINCHIFRAC
CUT2 = FRAC.LT.MAXFRAC
CUT3 = COSTHETA.GT.MINCOSTHETA
CUT4 = ABS(DELTAPT).LT.MAXDPT
CUT5 = ABS(DELTAPZ).LT.MAXDPZ
CUT6 = (CUT4.AND.CUT5).OR.CUT3
C
KEEPMOTHER = (CUT1.OR.CUT2).OR.CUT6
C
ENDIF
C
999 CONTINUE
RETURN
END
C
C
C
C*DFKINK
C==========================================================================
SUBROUTINE DFKINK(NTRACK,KINKPOS,KINK_FLAG)
C========================================================================
C
C INPUT PARAMETERS:
C
C Int ntrack ! number of track (= # of DPRS bank)
C
C OUTPUT PARAMETERS:
C
C Int kinkpos ! kink position from algoritm
C ! in PR candidate
C
C Int kink_flag ! 0 if there is no kink
C ! -1 if there is kink
C======================================================================
C=
C= Purpose and Methods : This routine scans the partially fitted
C= track 'ntrack' looking for a kink.
C= If a 'kink candidate' is found, the vectors
C= of parameters needed to write the new
C= DPRS banks in DFTWOPR are stored in
C= common /PAR/.
C= In any case the track 'ntrack' is left
C= unaffected; decision wether to keep it
C= or its daughters is taken later, in routine
C= DFSPLIT.
C=
C= Created 7-4-1998 F.Ambrosino C. Di Donato
C= Updated 8-5-1998 F.Ambrosino C. Di Donato
C= Updated 10-1998 F.Ambrosino
C- Updated 4-12-1998 (New cuts definition) F. A.
C=
C=========================================================================
C
C
C Global Declarations:-
C ====================
C
$$IMPLICIT
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'K$INC:BPOYBS.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$INC:OFERCO.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:CDFWORK.INC'
$$INCLUDE 'K$ITRK:CFLAGS.INC'
$$INCLUDE 'K$ITRK:CDFCUT.INC'
$$INCLUDE 'C$INC:jobsta.INC'
$$INCLUDE 'K$ITRK:ctidri.INC'
C
C
C
C
C External Functions:-
C ==================
C
INTEGER BLOCAT
INTEGER LVMAXA
REAL VSUM
C
C Local Declarations:-
C ===================
C
C---------------- COMMON FOR NEW TRACKS PARAMETERS ------------------------
INTEGER MAX_PAR
PARAMETER (MAX_PAR=9)
INTEGER MAX_HITS
PARAMETER (MAX_HITS=250)
C
REAL PAR1
INTEGER NPOS_1, NNEG_1, NHITS1
REAL PAR2
INTEGER NPOS_2, NNEG_2, NHITS2
REAL STEPS
C
COMMON /PAR/ PAR1(MAX_PAR), NPOS_1, NNEG_1, NHITS1,
+ PAR2(MAX_PAR), NPOS_2, NNEG_2, NHITS2,
+ STEPS(MAX_HITS)
C-----------------------------------------------------------------------
INTEGER SIGN_SRES(MAX_HITS)
& ,ZERO_CROSS(MAX_HITS)
& ,KPOS(MAX_HITS)
& ,CAND(MAX_HITS)
C
C
REAL RES(MAX_HITS)
& ,SMOOTHRES(MAX_HITS)
& ,SPOT_INI(MAX_HITS)
& ,SPOT_SIZE(MAX_HITS)
INTEGER NTRACK,KINK_FLAG,KINKPOS
INTEGER STATUS, IND_DPRS, INDAT_DPRS, I
INTEGER K_HIT, N , K
C
REAL CHISQ, SMEAN, RSIGMA, RDEV, AVRES
C
INTEGER KINK_FND, N_K, NZERO_CROSS
REAL sign,fraczero,diffres
C
INTEGER nsp, maxspotsize, last_cross
REAL AVSPOTSIZE, AV_GAP
REAL zero
INTEGER KKK,NN
real thisres
real pvres1, pvres2, nxres1, nxres2
real presm1, nresm1
real diff1, diff2
REAL MAXRES
C
INTEGER STUMP1, STUMP2, MINSTUMP
C
C
LOGICAL CONS
LOGICAL LOG1
LOGICAL LOG2
C
C---------------------------------------------Cuts definition
real sigma_min, chi_min,cut_cross
real cut_spsize, cut_zero, cut_gap
real sharp
PARAMETER(SIGMA_MIN = 0.055)
PARAMETER(CHI_MIN = 7.0)
PARAMETER(CUT_CROSS = 0.4)
PARAMETER(CUT_GAP = 3)
PARAMETER(CUT_SPSIZE = 0.17)
PARAMETER(CUT_ZERO = 20.)
PARAMETER(SHARP = 0.25)
C-----------------------------------------EXECUTABLE CODE
C
C-----------------------Please let us not split twice!!!-------
C
STATUS = BLOCAT(IW,'DPRS',NTRACK,IND_DPRS,INDAT_DPRS)
IF(STATUS.NE.YESUCC) THEN
CALL ERLOGR('DFKINK',ERWARN,0 ,0
$ ,'DPRS BANK NOT FOUND ??')
KINK_FLAG = 0
RETURN
ENDIF
C
IF(IW(INDAT_DPRS+DPRVRN).GT.1) THEN
KINK_FLAG = 0
RETURN
ENDIF
C
C===================================================Inititalization
KINK_FLAG = 0
KINKPOS = 0
K_HIT = 0
C
C
C--------------------Fill working vectors
N = NDF1
DO WHILE (N.GT.0)
k_hit = k_hit + 1
RES(k_hit) = DFCH(N)
N = NDF(N)
C
ENDDO
C
CHISQ = DFCHIN
C===============================================================IS THERE A KINK IN THIS TRACK?
C
C----------------------------------Number of hits
KINK_FND = -200
IF(K_HIT.LT.16) GO TO 999
KINK_FND = -300
C----------------------------------Chisq
C
C IF(CHISQ.LT.CHI_MIN) GO TO 999
IF(CHISQ/float(k_hit-5).LT.CHI_MIN) GO TO 999
KINK_FND = -400
C
C---------------------------------- Standard dev. (sigma) of residuals
rsigma = 0.
Avres = VSUM(res,k_hit)/float(k_hit)
C
do i = 1, k_hit
rdev = (Avres - res(i))**2
rsigma = rsigma + rdev
enddo
rsigma = sqrt(rsigma/float(k_hit - 1))
IF(RSIGMA.LT.SIGMA_MIN) GOTO 999
C
C
C
C------------------------------------ Smoothing the residuals distribution
C
SMEAN = 0.
do i = 1, k_hit
if(i.gt.1) then
SMEAN = (res(i) + res(i-1))/2.
else
SMEAN = res(i)
endif
SMOOTHRES(i) = SMEAN
enddo
C
C------------------------------------ Counting number of zero crossings
C
Call Vzero(zero_cross,k_hit)
nzero_cross = 0
do i=1,(k_hit-1)
sign = smoothres(i)*smoothres(i+1)
if(sign.lt.0) then
nzero_cross = nzero_cross + 1
zero_cross(i) = 1
endif
enddo
C
kink_fnd = -500
fraczero = float(nzero_cross)/float(k_hit)
If(fraczero.gt.cut_cross) go to 999
C
C----------------------------------------- Maximum spot size
C
Call VZERO(SPOT_SIZE,K_HIT)
CONS = .FALSE.
nsp = 0
maxspotsize = 0
last_cross = -99
C
Do i = 1,k_hit-1
If(zero_cross(i).gt.0.) then
cons = (i - last_cross).le.2
if(cons) then
spot_size(nsp) = spot_size(nsp) + 1
if(spot_size(nsp).gt.maxspotsize) then
maxspotsize = spot_size(nsp)
endif
last_cross = i
else
nsp = nsp + 1
spot_size(nsp) = spot_size(nsp) + 1
if(spot_size(nsp).gt.maxspotsize) then
maxspotsize = spot_size(nsp)
endif
last_cross = i
Spot_ini(nsp) = float(i)
endif
endif
Enddo
C
Avspotsize = float(maxspotsize)/float(k_hit)
C
kink_fnd = -600
If(Avspotsize.gt.cut_spsize) go to 999
C
C--------------------------------------------Average gap
Av_gap = 0.
If (nsp.lt.2) then
go to 999
endif
C
Av_gap = spot_ini(1) - 1.
do i = 2,nsp
Av_gap = Av_gap +
& (spot_ini(i)-spot_ini(i-1)) - spot_size(i-1)
enddo
Av_gap = Av_gap + (k_hit - spot_ini(nsp) -
& spot_size(nsp))
C
Av_gap = Av_gap/float(nsp+1)
C
kink_fnd = -700
If(Av_gap.lt.cut_gap) go to 999
C
C=============================================================YES, BUT WHERE IS IT ???
C
n_k = 0
C---------------------------Study of the derivative
C
zero = sqrt(2.)*rsigma/cut_zero
do i=1,(k_hit-1)
diffres = (smoothres(i+1)-smoothres(i))
if(diffres.lt.-zero) then
sign_sres(i) = -1
else if(diffres.lt.zero) then
sign_sres(i) = 0
else
sign_sres(i) = 1
endif
enddo
C
C
do i=1,(k_hit-1)
sign = (sign_sres(i+1)*sign_sres(i))
if(sign.lt.0.) then
n_k = n_k + 1
kpos(n_k) = i+1
endif
enddo
C
C
C--------------------------------At most one kink ?
IF(n_k.eq.0) THEN
kink_fnd = -800
go to 999
ELSEIF(n_k.eq.1) THEN
if(kpos(1).lt.5) then
kink_fnd = -830
go to 999
elseif(kpos(1).gt.(k_hit-5)) then
kink_fnd = -860
go to 999
else
kink_fnd = kpos(n_k)
go to 777
endif
ENDIF
C--------------------------We have more than one kink cand.
C
C-------------Filter out external kinks
C
do k = 1, n_k
if(kpos(k).lt.5) then
kpos(k) = -kpos(k)
go to 666
elseif(kpos(k)
& .gt.(k_hit-5)) then
kpos(k) = -kpos(k)
go to 666
endif
666 continue
enddo
C
nn = 0
do k = 1, n_k
if(kpos(k).gt.0) then
nn = nn + 1
cand(nn) = k
endif
enddo
C
if(nn.eq.0) then
kink_fnd = -900
go to 999
elseif(nn.eq.1) then
kink_fnd = kpos(cand(nn))
go to 777
endif
C
C------------Study sharpness of local maximum or minimum
C
C
maxres = -10000.
DO i = 1,nn
KKK = kpos(cand(i))
if(KKK.le.0) go to 111
thisres = res(kkk)
pvres1 = res(kkk - 1)
pvres2 = res(kkk - 2)
nxres1 = res(kkk + 1)
nxres2 = res(kkk + 2)
presm1 = (pvres1+pvres2)/2.
nresm1 = (nxres1+nxres2)/2.
diff1 = ((presm1)-(thisres))
diff2 = ((nresm1)-(thisres))
log1 = (abs(diff1).ge.(sharp*rsigma))
log2 = (abs(diff2).ge.(sharp*rsigma))
if(log1.or.log2) then
if (abs(res(KKK)).ge.maxres) then
maxres = abs(res(KKK))
kink_fnd = KKK
endif
endif
111 continue
ENDDO
C
C-------------If everything else fails, find out the absolute
C maximum (minimum) of residuals
C
IF(MAXRES.lt.0.) then
kink_fnd = lvmaxa(res(kkk),k_hit)
endif
C
C==============================================================WE FOUND IT ! IT`S KINK_FND
C
777 continue
C
STUMP1 = k_hit - kink_fnd
STUMP2 = kink_fnd
MINSTUMP = MIN(STUMP1,STUMP2)
C
C
if (MINSTUMP.lt.6) then
KINK_FLAG = 0
else
KINK_FLAG = -1
endif
C
IF(KINK_FLAG.ge.0) go to 999
C
KINKPOS = KINK_FND
C
C ==============================================================
C Before hit addition, hit rejection, and so on,
C store in a vector the information to be put
C in the new dprs banks
C
C
C
NPOS_1 = 0
NNEG_1 = 0
NPOS_2 = 0
NNEG_2 = 0
N = NDF1
DO WHILE(N.GT.0)
IF(N.LT.KINKPOS) THEN
STEPS(N) = DFS(N)
IF(TIDST(LDFL(N)).GT.0.) THEN
NPOS_1 = NPOS_1 + 1
ELSE
NNEG_1 = NNEG_1 + 1
ENDIF
ELSE
STEPS(N) = DFS(N)
IF(TIDST(LDFL(N)).GT.0.) THEN
NPOS_2 = NPOS_2 + 1
ELSE
NNEG_2 = NNEG_2 + 1
ENDIF
ENDIF
N=NDF(N)
ENDDO
C
NHITS1 = NPOS_1 + NNEG_1
NHITS2 = NPOS_2 + NNEG_2
C
C-----------------------------Last step of first stump must be zero!
STEPS(KINKPOS-1) = 0.
C
C
PAR1(1) = DFX(NDF1)
PAR1(2) = DFY(NDF1)
PAR1(3) = DFZ(NDF1)
PAR1(4) = DFX(KINKPOS-1) ! This is right only if this routine
PAR1(5) = DFY(KINKPOS-1) ! is called BEFORE any addition/rejection
PAR1(6) = DFZ(KINKPOS-1) ! of hits, so NDF is still SEQUENCIAL
PAR1(7) = DFKURV
PAR1(8) = DFPHI
PAR1(9) = DFCTGT
C
C
PAR2(1) = DFX(KINKPOS) ! This is right only if this routine
PAR2(2) = DFY(KINKPOS) ! is called BEFORE any addition/rejection
PAR2(3) = DFZ(KINKPOS) ! of hits, so NDF is still SEQUENCIAL
PAR2(4) = DFX(NDFL)
PAR2(5) = DFY(NDFL)
PAR2(6) = DFZ(NDFL)
PAR2(7) = DFKURV
PAR2(8) = ATAN2(DFPY(KINKPOS), DFPX(KINKPOS))
PAR2(9) = DFCTGT
C
C
999 CONTINUE
C
IF(KINK_FND.lt.0) then
Call HFILL(71,-FLOAT(KINK_FND),0.,1.)
ELSE
CALL HFILL(72,FLOAT(KINK_FND),0.,1.)
ENDIF
C
C
C
RETURN
END
[KLOE]
[Offline Doc]
[TRK Files]
Generated with Light on Thu Apr 8 13:00:16 MET DST 1999
.
Mail comments and suggestions.
