C=============================================================================== SUBROUTINE CTVMFT (PRINT,LEVPRT, IERR) C=============================================================================== C Authors: John Marriner & J.P. Berge C Original Creation: 22 Oct 4004 B.C.E. J.P. Marriner C revision: 1 Sep 2746 A.U.C. J.P. Berge Multiple vertices C revision: 10 Rabi` 1415 A.H. W.J. Ashmanskas Conversion fitting, C 1-track vertex C CDF multiple vertex and mass constrained fit package. C This package is documented in CDF/DOC/SEC_VTX/PUBLIC/1996 C (which needs revision) C Allows the specification of tracks at one or more vertices, and performs C (geometrical) fits such that at each vertex all tracks at that vertex are C constrained to pass through the same space point ("Vertex fits"). C There is a special case of the vertex fit that is designed to treat vertices C with two tracks where the two tracks originate from a photon conversion. C In this case, the two tracks have such a small space opening angle that C the usual vertex fit becomes unstable. For a vertex flagged as a conversion C candidate, extra constraints are imposed. The r-phi opening at the point C where they meet in space. The further requirement that the r-z opening angle C is zero may be optionally imposed. C These fits are not optional. They are always performed. C Allows constraining sub-groups of tracks to specified masses ("mass fit"). C Tracks in such a mass fit may be attached to different vertices (consider C the case of the decay Bd -> (psj)(k0). where the Bd mass is constra=ined). C These fits are optional. C Allows constraints such that the momentum sum of tracks at a vertex "points" C to another vertex ("pointing fits"). C "Target" vertices may be treated as "primary"; In these cases, vertex C coordinates and covariances are required input and are treated as C measurements in the fit. Such vertices have no attached tracks. C Alternatively, target vertices may be one of the vertices formed from C track intersection fitting in the current fit. Such vertices are considered C "secondary". A secondary vertex that is pointed at a secondary vertex is C a tertiary vertex. There is allowed a special case of a "one-track" vertex C (consider decays Xi -> (lambda)(pi) or B* -> (D0)(pi)). Here the fit C constrains the momentum vector from the tertiary vertex to intersect the C single track at the secondary vertex. C Note that "conversion vertices may point to parent vertices. Note that the C two tracks at such a vertex may appear in a multi-vertex mass-constrained C fit (consider the case Chi -> (psj)(gamma->ee)), but may not appear in any C mass constraint fit restricted only to this vertex. C These fits are optional. C NOTA BENE: Considerable care has been exercized in designing a formatted C dump of the CTVMFT fit results. This dump is activated through setting the C imput parameter PRINT to some non-zero value. If this effort has achieved C its purpose, many points which may remain unclear to the first time user C should be clarified through inspection of dumps of a few fits of interesting C fits. This formattewd output should be printed with 124 columns. C--Input parameters C PRINT If non-zero, dumps formatted fit results to unit PRINT C LEVPRT Controls the level of (optional) printing. Inoperative if PRINT=0. C Avoid this (other than 0,1) unless you know what you are doing. C--Input data communicated through the include file CTVMFT (COMMON /CTVMFC/) C Gross fit specifications (INTEGER); C NTRACK Number of tracks in this fit. (must be at least 2) C NVERTX Number of vertices in this fit (at least one is required) C NMASSC Number of mass constraints (may be zero) C Required topological and constraint information; C TrkVtx(MaxTrk,MaxVtx) LOGICAL C Specifies vertex to which each track is associated. C TrkVtx(Nt,Nv) is .TRUE. if track Nt is attached to vertex Nv. C VtxPnt(MaxVtx,2) INTEGER C Vertex association fit specification. C If Mv=VtxPnt(Nv,1) < 0, no pointing fit is done. C If Mv=VtxPnt(Nv,1) (Nv>Mv) is .GE.0, vertex Nv "points" to vertex Mv. C If Mv=0, the target vertex Mv is supplied as a measurement; the vertex C (x,y,z) together with its 3x3 covariance matrix must be furnished by C the user through the REAL arrays XYZPV0,EXYZPV. "Primary" vertices C are the most usual type of such vertices. C If Mv > 0, the target vertex Mv is a secondary vertex and its coord- C inates are parameters resulting from a track intersection fit. C The kind of "pointing" of vertex Nv to vertex Mv is specified by C setting VtxPnt(Nv,2) to 1, 2, or 3 as desired. C Note that more than one vertex can "point" to the same target vertex. C The constraint is that the sum of the momenta of all the tracks at C vertex Nv together with the summed momenta of all tracks at vertices C which point to Nv must be parallel OR anti-parallel to the vector C from the vertex Nv to vertex Mv. Allowed values of VtxPnt(Nv,2) are; C 1 = constrain pointing of Nv towards Mv in the (R,Phi) plane. C 2 = ALSO constrain pointing in (R,Z) plane (3 dimensional pointing). C 3 = point Nv toward single-track vertex Mv C Note that vertices connected via a charged track are not correctly C handled. C Cvtx(MaxVtx) INTEGER C Conversion constraint fit specification. (The conversion fit C constrains the two tracks at a vertex to be parallel at the vertex.) C 0 = no conversion fit is performed, the vertex is treated normally. C 1 = constrain the track directions in the r-phi plane to be the same C at the radius of the vertex (at the point where the tracks are C copunctual, they are clinear). C 2 = ALSO constrain opening angle to zero in r-z. C TrkMcn(MaxTrk,MaxMcn) LOGICAL C Specifies the tracks that participate in each mass constraint. C Track Nt participates in constraint Nm if TrkMcn(Nt,Nm) is .TRUE. C CMASS(MaxMcn) REAL C The Nmth cluster of tracks is constrained to have mass CMASS(Nm). C Track bank and data specifications; C TKBANK(MaxTrk) CHARACTER*4 C Bank name (exempli gratia, 'QTRK','TRKS','SVXS','CTCS') from which C to get each track's helix parameter and covariance data. C These bank names are used in any substantive way only in the data C access subroutine. In CDF, this routine is C$TRK:GETTRK. C LIST (MaxTrk) INTEGER C Bank numbers of the tracks to be used. C TMASS (MaxTrk) REAL C Masses (GeV/c**2) to be associated with the tracks, for purposes of C dE/dX and multiple Coulomb scattering calculations, and for any mass- C constrained fitting. C Output parameter C IERR = 0 (No error), or error code C Returned through the common /CTVMFC/ C PAR(5,NTRACK) = The track parameter fit results C XYZVRT(3,0:MaxVtx) The vertex x,y,z fit results. C VMAT(Maxdim,Maxdim) The covariance matrix on the fit parameters C FMCDIF(NMASSC) Residuals of the mass constraints (if used) C PCON(MaxVtx,1) Residual of the (R,Phi) pointing constraint (if used) C PCON(MaxVtx,2) Residual of the (R,Z) pointing constraint (if used) C SANG(MaxVtx,1) sin(ang) between vertex direction and momentum in x-y C SANG(MaxVtx,2) sin(ang) between vertex direction and momentum in rho-z C CHISQR(0:MaxItr) fit Chi Square result(s) C ... other goodies ... Go look at the include file. C Convenient access to VMAT entries for the fit results for vertex and track C parameters is provided by the offset pointer lists VOFF and TOFF. C VOFF(Nv+i) (i=1:3) gives row/column indices for vertex Nv (x,y,z). C TOFF(Nt+i) (i=1:3) gives row/column indices for track Nt (Crv,Phi,Ctg) C Subroutines; C CTVM00 Checks the input for completeness and consistency C CTVMFA Makes the vertex first approximations (uses GETTRK). C CTVM01 Collects the track data (uses GETTRK). C CTVMVF Calculates the derivative contributions for vertex fitting. C CTVMCF conversion fitting. C CTVMPF pointing fittting. C CTVMMF mass constraint fitting. C CTVMPR Prints the fit results, if desired, to output unit PRINT. C------------------------------------------------------------------------------- IMPLICIT NONE #include external ctvmft_bd c INTEGER PRINT, IERR LOGICAL PRIMARY DOUBLE PRECISION VXI(MAXDIM) DOUBLE PRECISION WXYZPV(3,3), SUM REAL C,PHI,COTH, D,Z, PT, RADV, TWOC, S,TS,TRAD, T1,T2,T3,T4 REAL SP,CP, XCYS,XSYC,XC2YC2, SINCS,COSCS,TWOCS, PMAG, XMAG REAL XSVT,YSVT, WORK(MAXDIM) REAL PARFAC,CUTFAC, DELVS(3,MaxVtx) REAL PCSAVE(MaxVtx,2), MFSAVE(MaxMcn), MFTEST(MaxMcn) REAL STPTK,STEP(5), PARDJF(5,MaxTrk) INTEGER Nv,Mv,Nt,Np,Nm, NvF,NtF,LpF,LmF,LcF INTEGER NITER, NMAXCT,ISCUT,NPC INTEGER I,J,K,L, II(2), I1,I2, IFAIL, IMAX, VERTEX,VRTX INTEGER KTRACK,NPOINT,NCONV, KPRINT C Local variables pertaining to single-track vertices INTEGER STVCOUNT,STVFLAG(MaxVtx) INTEGER BEAM LOGICAL FIRST INTEGER MRUN,MNEV, me CHARACTER*6 NAME CHARACTER*80 STRING C LEVPRT = print level (0=print disabled) INTEGER LEVPRT REAL ZERO,MINUS,SHOUT, TEST SAVE FIRST,BEAM SAVE NMAXCT, MRUN,MNEV SAVE ZERO,MINUS DATA ME /0/ DATA ZERO /0.0/, MINUS /-1.0/ C Maximum number of cut steps allowed in any step DATA NMAXCT / 7 / DATA FIRST /.TRUE./ DATA MRUN /-1/, MNEV /-1/ C******************************************************************************* DBGPRT = IABS(LEVPRT) C avoid terminating on error in CTVM00 EXCUSE = ME IF (LEVPRT.LT.0) THEN EXCUSE =-1 END IF IF (PRINT.EQ.0) DBGPRT=0 KPRINT = 0 IF (PRINT.GT.0 .AND. DBGPRT.GT.10) KPRINT=PRINT C First entrance initializations IF (FIRST) THEN FIRST =.FALSE. C do NOT use beam constrained parameters BEAM = 0 END IF IF (MRUN.NE.RUNNUM .OR. MNEV.NE.TRGNUM) THEN C New event, reset track "memory" DO I=1,MaxTrk TKERR(I) =-1 END DO MRUN = RUNNUM MNEV = TRGNUM END IF IERR = 0 DO I=1,3 IJKERR(I) = 0 END DO NTSCUT = 0 CALL VZERO(CHISQR,MaxItr+1) CHISQR(0) =-1. CALL VZERO(CHIV ,MaxVtx+1) CALL VZERO(CHIT ,MaxTrk) CALL VZERO(PARDIF,MaxTrk*5) CALL VZERO(PARDJF,MaxTrk*5) CALL VZERO(PARERR,MaxTrk*5) CALL VZERO(TrkP4, MaxTrk*6) CALL VZERO(DDA, MaxTrk*8) CALL VZERO(XYZVRT,MaxVtx*3+3) CALL VZERO(DELVS, MaxVtx*3) CALL VZERO(PCSAVE,MaxVtx*2) DO Nm=1,MaxMcn MFSAVE(Nm) = 0.0 MFTEST(Nm) = 0.0 END DO CALL VZERO(STVFLAG,MaxVtx) STVCOUNT = 0 C Check Input Conditions for illegal conditions C------------------------------------------------------------------------------- C Check input specifications CALL CTVM00 (IERR,WXYZPV) IF (IERR.NE.0) THEN RETURN END IF TRAD = ABS(XYZPV0(1)) + ABS(XYZPV0(2)) C flag the appearance of a primary vertex PRIMARY =.FALSE. C check for a primary vertex DO Nv=1,NVERTX IF (VtxPnt(Nv,1).EQ.0) THEN PRIMARY =.TRUE. END IF END DO C primary vertex coordinates IF (PRIMARY .OR. TRAD.GT.0.0) THEN DO I=1,3 XYZVRT(I,0) = XYZPV0(I) C primary vertex fit displacement DXYZPV(I) = 0.0 END DO END IF c IF (DBGPRT.GT.0) THEN c WRITE(PRINT,1000) RUNNUM,TRGNUM, NTRACK,NVERTX,NMASSC c IF (PRIMARY) THEN c WRITE(PRINT,1005) XYZPV0 c DO I=1,3 c WORK(I) = SQRT(EXYZPV(I,I)) c END DO c DO I=1,3 c DO J=1,3 c VMAT(I,J) = EXYZPV(I,J)/(WORK(I)*WORK(J)) c END DO c END DO c WRITE(PRINT,1006) (WORK(I),I=1,3) c WRITE(PRINT,1006) ((VMAT(I,J),I=1,3),J=1,3) c END IF c WRITE(PRINT,1001) (I, LIST(I),TkBank(I),Tmass(I),I=1,NTRACK) c1000 FORMAT('1Event ',I5,'.',I6,'; CTVMFT fit results,' c &,' (Tracks,Vertices,Mass Constraints) ',2(I2,','),I2, /) c1001 FORMAT(' Track',I3,';',I5,A5,F8.3 ) c1005 FORMAT( ' Primary vertex ',3F10.4) c1006 FORMAT(16X,3F10.4) c END IF C Vertex fit matrix dimensions; C 3 vertex parameters (x,y,z) per (possible) primary vertex C 3 vertex parameters per secondary vertex C + 0, 1, or 2 "pointing" parameters C + 3 parameters (curvature,phi0,coth(theta0)) per track C + N mass constraints C VOFF(Nv),TOFF(Nt) are offset pointers to ((Vertex,Pointing),Track) parameters C in (the 1st derivative vector, 2nd derivative matrix) VXI, VMAT C VOFF(Nv) points to the start of the 3, 4, or 5 variables for vertex Nv; C vertex Nv (x,y,z, p1,p2) C TOFF(Nt) points to the start of the 3 parameters; track Nt (Crv,Phi,Ctg) C POFF(Lp) points to the start of the Lagrangian multiplier values for the C pointing constraint Lp. C COFF(Lc) points to the start of the Lagrangian multiplier values for the C conversion constraint Lc. C MOFF points to the start of the Lagrangian multiplier values for the C mass constraints C initialize the parameter offset pointer NvF = 0 IF (PRIMARY) NvF = 3 C offsets to vertex parameters DO Nv=1,NVERTX VOFF(Nv) = NvF NvF = NvF + 3 END DO NtF = NvF C offsets to track parameters DO Nt=1,NTRACK TOFF(Nt) = NtF NtF = NtF + 3 END DO C count the number of "pointing" constraints NPOINT = 0 LpF = NtF C offsets to pointing Lagrangian multipliers DO Nv=1,NVERTX IF (VtxPnt(Nv,1).GE.0) THEN POFF(Nv) = LpF IF (VtxPnt(Nv,2).EQ.1.OR.VtxPnt(Nv,2).EQ.2) THEN NPOINT = NPOINT + VtxPnt(Nv,2) LpF = LpF + VtxPnt(Nv,2) C "one-track vertex" ELSE IF (VtxPnt(Nv,2).EQ.3) THEN IF (STVFLAG(VtxPnt(Nv,1)).EQ.0) THEN C Count number of unique single-track vertices STVFLAG(VtxPnt(Nv,1)) = 1 STVCOUNT = STVCOUNT+1 ENDIF NPOINT = NPOINT + 2 LpF = LpF + 2 END IF END IF END DO C count number of "conversion" constraints NCONV = 0 LcF = LpF C offsets to zero-opening-angle constraint DO Nv=1,NVERTX C Lagrangian multipliers IF (Cvtx(Nv).GT.0) THEN COFF(Nv) = LcF NCONV = NCONV + Cvtx(Nv) LcF = LcF + Cvtx(Nv) END IF END DO C offset to mass cnst Lagrangian multipliers MOFF = LcF C Fit dimension (number of degrees of freedom) NDOF = 2*NTRACK - 3*NVERTX + NPOINT+NMASSC+NCONV C Dimension of matrix equations MATDIM = NtF + NMASSC+NPOINT+NCONV C Data Collection and Checking, Fit Initialization, Vertex first approximation C------------------------------------------------------------------------------- C Loop (backward) over secondary vertices DO Nv=NVERTX,1,-1 CALL CTVMFA (Nv, PRINT,IERR) IF (IERR.NE.0) THEN IF (PRINT.GT.0) then call print_v('vertex first approximation failure %4i', $ ijkerr(1)) call print_v('%3i',ijkerr(2)) call print_v('%3i\n',ijkerr(3)) ENDIF GO TO 990 END IF END DO ITER = 0 NITER = MAXITR c----------------------------------------------------------------------- C fit iteration (re)entrance point c----------------------------------------------------------------------- 100 CONTINUE C (re)initialize VXI, VMAT DO I=1,MAXDIM VXI(I) = 0.0D0 DO J=1,MAXDIM VMAT(J,I) = 0.0D0 END DO END DO C stuff the primary Vtx covariance into VMAT IF (PRIMARY) THEN DO I=1,3 DO J=1,3 VMAT(I,J) = WXYZPV(I,J) C (measured - this fit) VXI(I) = VXI(I) - VMAT(I,J)*DXYZPV(J) ENDDO ENDDO ENDIF c----------------------------------------------------------------------- C Caclulate terms for the vertex fits via track inentersection. c----------------------------------------------------------------------- DO Nv=NVERTX,1,-1 CALL CTVMVF (KPRINT, Nv, VXI) ENDDO c----------------------------------------------------------------------- C accumulate vertex momentum sums, for (possible) pointing constraint c fitting 4-momentum sums for vertex pointing c----------------------------------------------------------------------- CALL VZERO(VtxP4,MaxVtx*4) DO Nv=NVERTX,1,-1 DO Nt=1,NTRACK IF (TrkVtx(Nt,Nv)) THEN DO I=1,4 VtxP4(I,Nv) = VtxP4(I,Nv) + TrkP4(Nt,I) ENDDO ENDIF ENDDO C does Nv point? Mv = VtxPnt(Nv,1) C tertiary; add P(Nv) to P(Mv) IF (Mv.GT.0) THEN DO I=1,4 VtxP4(I,Mv) = VtxP4(I,Mv) + VtxP4(I,Nv) ENDDO ENDIF ENDDO c----------------------------------------------------------------------- C calculate terms for "conversion" (zero opening angle) constraint c----------------------------------------------------------------------- DO Nv=1,NVERTX C conversion constraint, this vertex IF (Cvtx(Nv).GT.0) THEN CALL CTVMCF (KPRINT, Nv, VXI) ENDIF ENDDO c----------------------------------------------------------------------- C calculate terms for pointing from this vertex towards its parent c----------------------------------------------------------------------- DO Nv=NVERTX,1,-1 C pointing constraint, this vertex IF (VtxPnt(Nv,1).GE.0) THEN CALL CTVMPF (KPRINT, Nv, VXI) C save the last step quality LpF = POFF(Nv) PCSAVE(Nv,1) = ABS(VXI(LpF+1)) PCSAVE(Nv,2) = ABS(VXI(LpF+2)) ENDIF ENDDO c----------------------------------------------------------------------- C accumulate momentum sums, for (possible) mass constraint fitting C 4-momentum sums for mass constraints c----------------------------------------------------------------------- CALL VZERO(McnP4,MaxMcn*4) DO Nm=1,NMASSC DO Nt=1,NTRACK IF (TrkMcn(Nt,Nm)) THEN DO I=1,4 McnP4(I,Nm) = McnP4(I,Nm) + TrkP4(Nt,I) ENDDO ENDIF ENDDO ENDDO c----------------------------------------------------------------------- C loop over mass constraints c----------------------------------------------------------------------- DO Nm=1,NMASSC CALL CTVMMF (KPRINT, Nm, VXI) LmF = MOFF MFSAVE(Nm) = ABS(VXI(LmF+Nm)) ENDDO c----------------------------------------------------------------------- C End of the derivative accumulation phase; Now solve matrix c equations. Find solution without finding inverse c----------------------------------------------------------------------- IF (ITER.EQ.0) THEN CALL DEQN(MATDIM,VMAT,MAXDIM,WORK,IFAIL,1,VXI) C Find solution + covariance matrix ELSE CALL DEQINV(MATDIM,VMAT,MAXDIM,WORK,IFAIL,1,VXI) ENDIF IF (IFAIL .NE. 0) THEN IERR = 9 IJKERR(2) = 1 call print_s('CTVMFT error: singular solution matrix\n') RETURN ENDIF c----------------------------------------------------------------------- C Check the step just taken for acceptability and for possible c convergence. Initialize cut steps counter c----------------------------------------------------------------------- ISCUT = 0 C Step size scale factor - normally 1 CUTFAC = 1.0 C Fraction of the step already taken PARFAC = 0.0 C Check that track turning angles to the new vertex points are acceptable 110 CONTINUE C x,y for this vertex, this step cut factor DO 150 Nv=1,NVERTX NvF = VOFF(Nv) XSVT = XYZVRT(1,Nv) + CUTFAC*VXI(NvF+1) YSVT = XYZVRT(2,Nv) + CUTFAC*VXI(NvF+2) C track loop, checking turning angle DO 140 Nt=1,NTRACK IF (.NOT.TrkVtx(Nt,Nv)) GO TO 140 NtF = TOFF(Nt) C = PAR0(1,Nt) + PARDIF(1,Nt) + CUTFAC*VXI(NtF+1) PHI = PAR0(2,Nt) + PARDIF(2,Nt) + CUTFAC*VXI(NtF+2) TS = 2.0*C*(XSVT*COS(PHI)+YSVT*SIN(PHI)) C track turns too much, cut step size IF (ABS(TS).GT.TRNMAX) THEN ISCUT = ISCUT+1 NTSCUT = NTSCUT + 1 C maximum number of cut steps exceeded IF (ISCUT.GT.NMAXCT) THEN IERR = 9 IJKERR(2) = 2 IJKERR(3) = ITER call print_v('too many cut steps : %5i',ijkerr(1)) call print_v('%3i',ijkerr(2)) call print_v('%3i\n',ijkerr(3)) GO TO 990 END IF CUTFAC = 0.5*(CUTFAC-PARFAC) GO TO 110 END IF 140 CONTINUE 150 CONTINUE PARFAC = PARFAC+CUTFAC c----------------------------------------------------------------------- C Provisionally acceptable step; C Update all the fit parameters and calculate Chi Square, this step C Re-initialize, recompute momentum sums with the new step results C Update primary vertex coordinates c----------------------------------------------------------------------- IF (PRIMARY) THEN DO I=1,3 DXYZPV(I) = DXYZPV(I) + CUTFAC*VXI(I) XYZVRT(I,0) = XYZPV0(I) + DXYZPV(I) END DO END IF C Contribution from the "primary" vertex CHIV(0) = 0.0 IF (PRIMARY) THEN DO I=1,3 DO J=1,3 CHIV(0)= CHIV(0) + DXYZPV(I)*WXYZPV(I,J)*DXYZPV(J) END DO END DO END IF CHISQR(ITER) = CHIV(0) C Update secondary vertex coordinates DO Nv=1,NVERTX Nvf = VOFF(Nv) DO I=1,3 XYZVRT(I,Nv) = XYZVRT(I,Nv) + CUTFAC*VXI(NvF+I) Mv = VtxPnt(Nv,1) C vertex separation vector IF (Mv.GE.0) THEN DELVS(I,Nv) = XYZVRT(I,Mv) - XYZVRT(I,Nv) END IF END DO END DO DO 210 Nv=1,NVERTX CHIV(Nv) = 0.0 C Update the track parameters DO 210 Nt=1,NTRACK IF (.NOT.TrkVtx(Nt,Nv)) GO TO 210 NtF = TOFF(Nt) C half curvature step PARDIF(1,Nt) = PARDIF(1,Nt) + CUTFAC*VXI(NtF+1) C phi step PARDIF(2,Nt) = PARDIF(2,Nt) + CUTFAC*VXI(NtF+2) C cot(theta) step PARDIF(3,Nt) = PARDIF(3,Nt) + CUTFAC*VXI(NtF+3) C new half curvature C = PAR0(1,Nt) + PARDIF(1,Nt) C new phi PHI = PAR0(2,Nt) + PARDIF(2,Nt) C new coth(theta) COTH = PAR0(3,Nt) + PARDIF(3,Nt) SP = SIN(PHI) CP = COS(PHI) C Xsv*cos(phi) + Ysv*sin(phi), -Xsv*sin(phi) + Ysv*cos(phi), 2C XCYS = XYZVRT(1,Nv)*CP + XYZVRT(2,Nv)*SP XSYC =-XYZVRT(1,Nv)*SP + XYZVRT(2,Nv)*CP TWOC = 2.0*C S = ASIN(TWOC*XCYS)/TWOC TS = TWOC*XCYS TRAD = SQRT((1.-TS)*(1.+TS)) C Constrained d0 D = XSYC - SIN(C*S)**2/C C Constrained z0 Z = XYZVRT(3,Nv) - COTH*S C fit d0 - CTC fit d0 PARDIF(4,Nt) = D - PAR0(4,Nt) C fit z0 - CTC fit z0 PARDIF(5,Nt) = Z - PAR0(5,Nt) PT = PSCALE/ABS(C) C x component of momentum TrkP4(Nt,1) = PT*(CP*TRAD-SP*TS) C y component of momentum TrkP4(Nt,2) = PT*(SP*TRAD+CP*TS) C z component of momentum TrkP4(Nt,3) = PT*COTH C Transverse momentum TrkP4(Nt,5) = PT C total momentum TrkP4(Nt,6) = SQRT(PT**2+TrkP4(Nt,3)**2) C energy TrkP4(Nt,4) = SQRT(TrkP4(Nt,6)**2+TMASS(NT)**2) C Calculate contribution to chi-squared for this track CHIT(Nt) = 0.0 C Loop over rows of the weight matrix DO I=1,5 C Loop over columns of the weight matrix DO J=1,5 CHIT(Nt) = CHIT(Nt) + PARDIF(I,Nt)*G(J,I,Nt)*PARDIF(J,Nt) END DO END DO CHIV(Nv) = CHIV(Nv) + CHIT(Nt) CHISQR(ITER) = CHISQR(ITER) + CHIT(Nt) C End the loop on tracks 210 CONTINUE CALL VZERO(VtxP4,MaxVtx*4) CALL VZERO(McnP4,MaxMcn*4) DO Nt=1,NTRACK DO Nv=NVERTX,1,-1 IF (TrkVtx(Nt,Nv)) THEN C accumulate vertex sums DO I=1,4 VtxP4(I,Nv) = VtxP4(I,Nv) + TrkP4(Nt,I) ENDDO ENDIF ENDDO c----------------------------------------------------------------------- C accumulate mass constraint sums c----------------------------------------------------------------------- DO Nm=1,NMASSC IF (TrkMcn(Nt,Nm)) THEN DO I=1,4 McnP4(I,Nm) = McnP4(I,Nm) + TrkP4(Nt,I) ENDDO ENDIF ENDDO ENDDO c DO Nv=1,NVERTX C does Nv point? Mv = VtxPnt(Nv,1) C tertiary; add P(Nv) to P(Mv) IF (Mv.GT.0) THEN DO I=1,4 VtxP4(I,Mv) = VtxP4(I,Mv) + VtxP4(I,Nv) ENDDO ENDIF ENDDO c----------------------------------------------------------------------- C Check the quality of satisfaction of the pointing constraints, c this step c----------------------------------------------------------------------- DO 220 Nv=1,NVERTX Mv = VtxPnt(Nv,1) IF (Mv.LT.0) GO TO 220 C IMAX points to larger PSUM component (1=x,2=y) IMAX = 1 IF (ABS(VtxP4(1,Nv)).LT.ABS(VtxP4(2,Nv))) IMAX = 2 NPC = VtxPnt(Nv,2) IF (NPC.EQ.3) NPC = 2 DO Np=1,NPC IF (Np.EQ.1) THEN C Residual of first pointing constraint, this step SUM = VtxP4(1,Nv)*DELVS(2,Nv) - VtxP4(2,Nv)*DELVS(1,Nv) PMAG = SQRT(VtxP4(1,Nv)**2+VtxP4(2,Nv)**2) XMAG = SQRT(DELVS(1,Nv)**2+DELVS(2,Nv)**2) ELSE C Residual of second pointing constraint, this step SUM = DELVS(3,Nv)*VtxP4(IMAX,Nv) - DELVS(IMAX,Nv)*VtxP4(3,Nv) PMAG = SQRT(VtxP4(IMAX,Nv)**2 + VtxP4(3,Nv)**2) XMAG = SQRT(DELVS(IMAX,Nv)**2+DELVS(3,Nv)**2) END IF PCON(Nv,Np) = SUM C Check for very small vertex separation IF (XMAG .GT. .001) THEN C OK. Get sin(angle) SANG(Nv,Np) = PCON(Nv,Np)/(PMAG*XMAG) ELSE SANG(Nv,Np) = 2.0 ENDIF END DO C end vertex loop 220 CONTINUE C Check the quality of satisfaction pf the mass constraints, this step DO Nm=1,NMASSC C Residual of mass constraint. SUM = SQRT(McnP4(1,Nm)**2+McnP4(2,Nm)**2+McnP4(3,Nm)**2) SUM = (McnP4(4,Nm)+SUM) * (McnP4(4,Nm)-SUM) IF (SUM.LT.0.0) THEN SUM = 0.0 ELSE SUM = SQRT(SUM) END IF FMCDIF(Nm) = (SUM + CMASS(Nm)) * (SUM - CMASS(Nm)) MFTEST(Nm) = 0.5 * FMCDIF(Nm) C Convert to Delta m in MeV from Delta m**2 in Gev FMCDIF(Nm) = 500.*FMCDIF(Nm)/CMASS(Nm) END DO c----------------------------------------------------------------------- C Check the improvements(?) in pointing and mass constraints, this step C If step size is too large the constraints will not be satisfied. C If there is a problem, the step size will be halved C Check pointing constraints c----------------------------------------------------------------------- DO 230 Nv=1,NVERTX Mv = VtxPnt(Nv,1) IF (Mv.LT.0) GO TO 230 NPC = VtxPnt(Nv,2) IF (NPC.EQ.3) NPC = 2 DO I=1,NPC C = PCON(Nv,I) / (PCSAVE(Nv,I) + 0.001) IF ((ABS(C).GT.1.0) .AND. (ABS(SANG(Nv,I)).GT..01)) THEN GO TO 250 END IF END DO C End vertex loop 230 CONTINUE c----------------------------------------------------------------------- C Check mass constraints c----------------------------------------------------------------------- DO Nm=1,NMASSC C = ABS(MFTEST(Nm) / (MFSAVE(Nm) + 0.1)) IF (C.GT.1.0) THEN c$$$ IF (DBGPRT.GT.1) THEN c$$$ WRITE(PRINT,1052) ITER,ISCUT,Nm,C,MFTEST(Nm) c$$$1052 FORMAT(' Iter',I2,' Cut',I2,', Nm',I2,2F10.5) c$$$ END IF GO TO 250 END IF END DO GO TO 300 250 CONTINUE C Count a cut step ISCUT = ISCUT + 1 NTSCUT = NTSCUT + 1 C maximum number of cut steps exceeded IF (ISCUT .GT. NMAXCT) THEN IERR = 9 IJKERR(2) = 3 IJKERR(3) = ITER call print_v('too many cut steps : %5i',ijkerr(1)) call print_v('%3i',ijkerr(2)) call print_v('%3i\n',ijkerr(3)) GO TO 990 END IF C Prepare to subtract off 1/2 step size taken CUTFAC = -0.5*PARFAC GO TO 110 C Accept this step 300 CONTINUE CHISQR(0) = CHISQR(ITER) STPTK = 1.0 IF (ITER.GT.0) THEN STPTK = 0.0 DO NT=1,NTRACK C make the difference between this step and the previous DO I=1,5 PARDJF(I,NT) = PARDIF(I,NT) - PARDJF(I,NT) STEP(I) = PARDJF(I,NT)/PARERR(I,NT) STPTK = AMAX1(STPTK,ABS(STEP(I))) END DO END DO END IF D IF (DBGPRT.GT.1) THEN D C = AMIN1(CHISQR(0),9999.0) D L = 1 D IF (PRIMARY) L=0 D DO I=L,NVERTX D WORK(I) = AMIN1(CHIV(I),999.0) D END DO D IF (NVERTX.EQ.1 .AND. L.EQ.1) THEN D WRITE(PRINT,1070) stptk, C D1070 FORMAT(/,' StpTk ',F8.3,'; ChiSqr =',F7.2) D ELSE D WRITE(PRINT,1071) stptk, C, (I,WORK(I),I=L,NVERTX) D1071 FORMAT(/,' StpTk ',F8.3,'; ChiSqr =',F7.2, 4(3X,'ChiV',I2,F7.2) ) D END IF D IF (ITER.EQ.0) GO TO 500 D DO NT=1,NTRACK D WORK(Nt) = AMIN1(CHIT(Nt),999.0) D DO I=1,5 D STEP(I) = PARDJF(I,NT)/PARERR(I,NT) D END DO D WRITE(PRINT,1075) nt, WORK(Nt), STEP D1075 FORMAT(i5,F10.4,2x,5F10.6) D END DO D END IF C Check for step acceptability 500 CONTINUE IF (ITER.EQ.0) GO TO 550 IF (STPTK.LT.0.1) THEN C accept the convergence GO TO 900 ELSE IF (ITER.GE.MAXITR) THEN C "convergence", sort of. Accept the stinker? IF (STPTK.LT.0.2) THEN GO TO 900 ELSE STPTK = AMIN1(STPTK,9999.0) IERR = 9 IJKERR(2) = 4 IF (PRINT.GT.0) then call print_v('; Convergence failure, Itr = %2i',iter) call print_v('; StpTk %7.2f',stptk) call print_v(' ChiSqr(i) %9.2e\n',chisqr(0)) endif GO TO 990 END IF ELSE IF (ITER.GT.2) THEN C after the 2nd step, if the "ChiSquare" TEST = FLOAT(NDOF) * 225.0 C is >15 "standard deviations", and if IF (CHISQR(0).GT.TEST) THEN C the last step was did not do much good IF (STPTK.GT.0.5) THEN STPTK = AMIN1(STPTK,9999.0) C give up in disgust IERR = 9 IJKERR(2) = 5 if (PRINT.GT.0) then call print_v('; Bad convergence, Itr = %2i',iter) call print_v('; StpTk %7.2f',stptk) call print_v(' ChiSqr(i) %9.2e\n',chisqr(0)) endif GO TO 990 END IF END IF END IF END IF 550 CONTINUE C return for the next iteration step DO NT=1,NTRACK C save the last step on this track DO I=1,5 PARDJF(I,NT) = PARDIF(I,NT) END DO END DO ITER = ITER + 1 GO TO 100 C Successful Return C------------------------------------------------------------------------------- 900 CONTINUE C final track parameters DO Nt=1,NTRACK DO I=1,5 PAR(I,NT) = PAR0(I,Nt) + PARDIF(I,Nt) END DO END DO IF (PRINT.GT.0) THEN CALL CTVMPR (PRINT,DBGPRT, PARERR) END IF c----------------------------------------------------------------------- C Check for peculiar case of non-positive diagonal matrix elements, which C should never happen, but if it happens, better to report it than to C let it slip by unnoticed. All errors of this sort should be reported C immediately to the guilty parties. c----------------------------------------------------------------------- DO I = 1,TOFF(NTRACK)+3 IF (VMAT(I,I).LE.0) THEN STRING = 'SERIOUS problem with ill-formed covariance matrix' IERR = 9 IJKERR(1) = 9 IJKERR(2) = 9 CALL ctvmft_ERROR('CTVMFT',IERR,STRING) IF (PRINT.GT.0) then call print_v('CTVMFT error: %2i',mrun) call print_v(' %2i',mnev) call print_v(' %12.5e \n',vmat(i,j)) endif GOTO 999 ENDIF ENDDO RETURN c----------------------------------------------------------------------- C Error Return C----------------------------------------------------------------------- 990 CONTINUE IF (IERR.NE.0) THEN IJKERR(1) = IERR IF (PRINT.GT.0) then call print_v('CTVMFT error: IERR = %3i =>\n',ierr) endif c IF (PRINT.GT.0 .OR. IERR.GT.50) THEN c CALL ctvmft_ERROR (NAME,IERR,STRING) c ENDIF ENDIF 999 CONTINUE RETURN END