c=======================================================================
SUBROUTINE VTXFIT
c=======================================================================
c
c This is the driving routine of the vertex fitting package ;
c it looks for vertex candidates and calls the minimization routine.
c If the fit is succsefull the candidate vertex is temporary stored
c in routine VTXCAN for later use (in VTXMERG).
c
c Author: Marco Incagli
c =======
c
c Creation date: 22-9-1997
c ==============
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c External declarations
c
c
c Local declarations
c
INTEGER Iextrvc, Iextrbp, iii, Iextr, Ipoinvc(20), Ipoinbp(20)
INTEGER ifi, ise, Ilist(2), IERRCR, Ifit, Nflag
INTEGER itk, ktk, Nreg, Ntrfit, i
REAL Qfi, Qse, Xint(3), Xcn(3), Chisq, Xstor(3)
REAL X1F, Y1F, Z1F, X1L, Y1L, Z1L, X2F, Y2F ,Z2F, X2L, Y2L, Z2L
REAL D1F2F, Z1F2F, D1F2L, Z1F2L, D1L2F, Z1L2F, D1L2L, Z1L2L
LOGICAL Tangent, Constr
LOGICAL FlagIP(50) ! common VTXFIT-VTXEXTR
COMMON /VTXCLOSIP/ FlagIP
c
c=======================================================================
c
c Initialize variables
c
Ntrfit = 2 ! Assume a "2 tracks fit" as default
Nvxfit = 0 ! Number of fitted vertices
c
c ===================================================================
c
c Organize tracks according to their extrapolated inner point:
c IEXTRVC are the tracks whose innermost extrapolated point falls
c inside the Vertex Chamber region; although KLOE has not a Vertex
c Chamber (as now) this is the zone included between the inner DC and
c the BP walls
c IEXTRBP are the tracks whose innermost extrapolated point falls
c inside the Beam Pipe
c
iextrvc = 0
iextrbp = 0
do iii = 1, TOTTRK
IEXTR = IW(INDVWRK + VWRNCO*(iii-1) + VWREG)
if (iextr.eq.NVC) then
iextrvc = iextrvc+1
ipoinvc(iextrvc) = iii
elseif(iextr.gt.NVC) then
iextrbp = iextrbp+1
ipoinbp(iextrbp) = iii
endif
enddo
cmi
C
C === LOOK FOR VERTICES NOW ======================================
C
c =================================
c (1) start from vertices inside BP
c =================================
c
c loop on first track
do 200 ifi = 1, iextrbp-1
QFI = SIGN(1.,
+ RW(INDVWRK+VWRNCO*(IPOINBP(ifi)-1)+VWCUBI)) ! track charge
c check if track has already been used (unless it passes close to the
c Interaction Point; these tracks can be used for more than one vertex)
if ( .not.FlagIP(IPOINBP(ifi)) .and.
+ IW(INDVWRK+VWRNCO*(IPOINBP(ifi)-1)+VWFLF).eq.FITGOLD) goto 200
c loop on second track
do 180 ise = iextrbp, ifi+1, -1
QSE = SIGN(1.,
+ RW(INDVWRK+VWRNCO*(IPOINBP(ise)-1)+VWCUBI)) ! track charge
c Check total charge (must be zero).
if ( abs(QSE+QFI) .gt. 0.5 ) goto 180
c check if track has already been used (unless it passes close to the
c Interaction Point; these tracks can be used for more than one vertex)
if ( .not.FlagIP(IPOINBP(ise)) .and.
+ IW(INDVWRK+VWRNCO*(IPOINBP(ise)-1)+VWFLF).eq.FITGOLD) goto 180
c
ILIST(1) = IPOINBP(ifi)*ISWMBACK
ILIST(2) = IPOINBP(ise)*ISWMBACK
Nreg = 3
c
call VTXCRSS(Nreg,ILIST,XINT,IERRCR,TANGENT)
c don't look for vertex if VTXCRSS doesn't find any intersection nor a
c "tangency within cuts" (see VTXCRSS for this definition)
if ((IERRCR.gt.1) .or.
+ (IERRCR.eq.1 .and. .not. TANGENT)) goto 180
c call fitting routine with two tracks
c
c if the tracks are parallel then:
c 1 - assume that they come from the origin (K+K- or Bhabha ...)
c
if ( TANGENT ) then
CONSTR = TANGENT
do i = 1, 3
xstor(i) = XINT(i)
XCN(i) = 0.
XINT(i) = 0.
enddo
endif
c
c call minimization procedure
c
NFLAG = 1 ! convergence occurred in loop (1)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
c
c if (1) fails then:
c 2 - fix the constraint at the crossing point
c
if (IFIT.ne.FITGOLD .and. TANGENT) then
XCN(1) = xstor(1)
XCN(2) = xstor(2)
XCN(3) = xstor(3)
XINT(1) = xstor(1)
XINT(2) = xstor(2)
XINT(3) = xstor(3)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
endif
CONSTR = .false.
c
180 continue
200 continue
C
C ====================================================================
C
c =================================
c (2) now vertices in VC region
c =================================
c
c loop on first track; it has to be of type IEXTRVC
do 201 ifi = 1, iextrvc
QFI = SIGN(1.,
+ RW(INDVWRK+VWRNCO*(IPOINVC(ifi)-1)+VWCUDC)) ! track charge
c check if the track has aready been used
if ( .not.FlagIP(IPOINVC(ifi)) .and.
+ IW(INDVWRK+VWRNCO*(IPOINVC(ifi)-1)+VWFLF).eq.FITGOLD) goto 201
c loop on second track: it can be either of type IEXTRVC or IEXTRBP
c first pass: IEXTRVC track
do 181 ise = 1, ifi-1
QSE = SIGN(1.,
+ RW(INDVWRK+VWRNCO*(IPOINVC(ise)-1)+VWCUDC)) ! track charge
if ( abs(QSE+QFI) .gt. 1. ) goto 181
if ( .not.FlagIP(IPOINVC(ise)) .and.
+ IW(INDVWRK+VWRNCO*(IPOINVC(ise)-1)+VWFLF).eq.FITGOLD) goto 181
c
ILIST(1) = IPOINVC(ifi)*ISWMBACK
ILIST(2) = IPOINVC(ise)*ISWMBACK
Nreg = 2
c
call VTXCRSS(Nreg,ILIST,XINT,IERRCR,TANGENT)
c don't look for vertex if VTXCRSS doesn't find any intersection nor a
c "tangency within cuts" (see VTXCRSS for this definition)
if ((IERRCR.gt.1) .or.
+ (IERRCR.eq.1 .and. .not. TANGENT)) goto 181
c call fitting routine with two tracks
cmi
cmi Apply a constrained fit if two tangents are parallel
cmi ===============
cmi
if ( TANGENT ) then
CONSTR = TANGENT
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
endif
c
c call minimization procedure
c
NFLAG = 2 ! convergence occurred in loop (2)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
c
181 continue
c second pass: IEXTRBP tracks
do 182 ise = 1, iextrbp
QSE = SIGN(1.,
+ RW(INDVWRK+VWRNCO*(IPOINBP(ise)-1)+VWCUBI)) ! track charge
if ( abs(QSE+QFI) .gt. 1. ) goto 182
c check if track has already been used (unless it passes close to the
c Interaction Point; these tracks can be used for more than one vertex)
if ( .not.FlagIP(IPOINBP(ise)) .and.
+ IW(INDVWRK+VWRNCO*(IPOINBP(ise)-1)+VWFLF).eq.FITGOLD) goto 182
c
ILIST(1) = IPOINVC(ifi)*ISWMBACK
ILIST(2) = IPOINBP(ise)*ISWMBACK
Nreg = 2 ! Vertex in VC region
c
call VTXCRSS(Nreg,ILIST,XINT,IERRCR,TANGENT)
c don't look for vertex if VTXCRSS doesn't find any intersection nor a
c "tangency within cuts" (see VTXCRSS for this definition)
if ((IERRCR.gt.1) .or.
+ (IERRCR.eq.1 .and. .not. TANGENT)) goto 182
c call fitting routine with two tracks
cmi
cmi Apply a constrained fit if two tangents are parallel - MI sep 97
cmi ===============
cmi
if ( TANGENT ) then
CONSTR = TANGENT
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
endif
c
c call minimization procedure
c
NFLAG = 2 ! convergence occurred in loop (2)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
c
182 continue
201 continue
C
C ======================================================================
C
c =================================
c (3) now vertices in DC region
c =================================
C
c In this section ALL tracks are combined to look for vertices. The
c distance, in r-phi and in z, of the First points, Last points and
c First-to-Last points are evaluated; a fit is attempted if DeltaR<RCUTM
c and DeltaZ<ZCUTM . A succesfull fit in previous sections is treated
c by putting Ztrack to a fake value, unless the track passes close
c to the IP (FlagIP.eq..true.).
c
do 85 itk = 1, TOTTRK-1
QFI = SIGN(1.,RW(INDVWRK+VWRNCO*(itk-1)+VWCUFH)) ! track charge
Z1F = 999.
Z1L = 999.
if ( IW(INDVWRK+VWRNCO*(itk-1)+VWFLF).ne.FITGOLD .or.
+ (FlagIP(itk) .and.
+ IW(INDVWRK+VWRNCO*(itk-1)+VWFLF).eq.FITGOLD) ) then
X1F = RW(INDVWRK+VWRNCO*(itk-1)+VWXFH)
Y1F = RW(INDVWRK+VWRNCO*(itk-1)+VWYFH)
Z1F = RW(INDVWRK+VWRNCO*(itk-1)+VWZFH)
endif
if ( IW(INDVWRK+VWRNCO*(itk-1)+VWFLL).ne.FITGOLD ) then
X1L = RW(INDVWRK+VWRNCO*(itk-1)+VWXLH)
Y1L = RW(INDVWRK+VWRNCO*(itk-1)+VWYLH)
Z1L = RW(INDVWRK+VWRNCO*(itk-1)+VWZLH)
endif
c
do 82 ktk = ITK+1, TOTTRK
QSE = SIGN(1.,RW(INDVWRK+VWRNCO*(ktk-1)+VWCUFH)) ! track charge
Z2F = -999.
Z2L = -999.
if ( IW(INDVWRK+VWRNCO*(ktk-1)+VWFLF).ne.FITGOLD .or.
+ (FlagIP(ktk) .and.
+ IW(INDVWRK+VWRNCO*(ktk-1)+VWFLF).eq.FITGOLD) ) then
X2F = RW(INDVWRK+VWRNCO*(ktk-1)+VWXFH)
Y2F = RW(INDVWRK+VWRNCO*(ktk-1)+VWYFH)
Z2F = RW(INDVWRK+VWRNCO*(ktk-1)+VWZFH)
endif
if ( IW(INDVWRK+VWRNCO*(ktk-1)+VWFLL).ne.FITGOLD ) then
X2L = RW(INDVWRK+VWRNCO*(ktk-1)+VWXLH)
Y2L = RW(INDVWRK+VWRNCO*(ktk-1)+VWYLH)
Z2L = RW(INDVWRK+VWRNCO*(ktk-1)+VWZLH)
endif
c distances first-first first-last last-last
D1F2F = sqrt( (X1F-X2F)**2+(Y1F-Y2F)**2 )
Z1F2F = abs ( Z1F - Z2F )
D1F2L = sqrt( (X1F-X2L)**2+(Y1F-Y2L)**2 )
Z1F2L = abs ( Z1F - Z2L )
D1L2F = sqrt( (X1L-X2F)**2+(Y1L-Y2F)**2 )
Z1L2F = abs ( Z1L - Z2F )
D1L2L = sqrt( (X1L-X2L)**2+(Y1L-Y2L)**2 )
Z1L2L = abs ( Z1L - Z2L )
c check if tracks from same vertex ; if yes then try to fit and store
c the check is done with all 4 combinations First-First, First-Last, ...
c note that in these combinations I check also the charge: for
c First-First and Last-Last the net charge has to be 0 (opposite sign)
c for First-Last combination the net charge has to be 2 (same sign)
c (1) 1 FIRST - 2 FIRST
if ( D1F2F.lt.RCUTMI .and. Z1F2F.lt.ZCUTMI .and.
+ abs(QFI+QSE).lt.1. ) then
c
c initial conditions for fitting routine (including swim direction)
c
ILIST(1) = itk*ISWMBACK
ILIST(2) = ktk*ISWMBACK
Nreg = 1
call VTXCRSS(Nreg,ILIST,XINT,IERRCR,TANGENT)
c don't look for vertex if VTXCRSS doesn't find any intersection nor a
c "tangency within cuts" (see VTXCRSS for this definition)
if ((IERRCR.gt.1) .or.
+ (IERRCR.eq.1 .and. .not. TANGENT)) goto 82
c
c Apply a constrained fit if two tangents are parallel - MI sep 97
c ===============
c
if ( TANGENT ) then
CONSTR = TANGENT
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
endif
c
c call minimization procedure
c
NFLAG = 3 ! convergence occurred in loop (3)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
c If minimization procedure converged then put fake Z value (=exclude
c this track end for next vertex searches); if fit is SILVER skip the
c rest of this section, but don't exclude other combinations (VTXMERG
c will eventually select the best one)
if (IFIT.eq.FITGOLD) then
Z1F = 999.
Z2F = -999.
goto 82
elseif (IFIT.eq.FITSILVER) then
goto 82
endif
c
c (2) 1 FIRST - 2 LAST
c
elseif ( D1F2L.lt.RCUTMI .and. Z1F2L.lt.ZCUTMI .and.
+ abs(QFI+QSE).gt.1. ) then
c
c initial conditions for fitting routine (including swim direction)
c
ILIST(1) = itk*ISWMBACK
ILIST(2) = ktk*ISWMFORW
Nreg = 1
call VTXCRSS(Nreg,ILIST,XINT,IERRCR,TANGENT)
c don't look for vertex if VTXCRSS doesn't find any intersection nor a
c "tangency within cuts" (see VTXCRSS for this definition)
if ((IERRCR.gt.1) .or.
+ (IERRCR.eq.1 .and. .not. TANGENT)) goto 82
c
c Apply a constrained fit if two tangents are parallel - MI sep 97
c ===============
c
if ( TANGENT ) then
CONSTR = TANGENT
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
endif
c
c call minimization procedure
c
NFLAG = 3 ! convergence occurred in loop (3)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
c If minimization procedure converged then put fake Z value (=exclude
c this track end for next vertex searches); if fit is SILVER skip the
c rest of this section, but don't exclude other combinations (VTXMERG
c will eventually select the best one)
if (IFIT.eq.FITGOLD) then
Z1F = 999.
Z2L = -999.
goto 82
elseif (IFIT.eq.FITSILVER) then
goto 82
endif
c
c (3) 1 LAST - 2 FIRST
elseif ( D1L2F.lt.RCUTMI .and. Z1L2F.lt.ZCUTMI .and.
+ abs(QFI+QSE).gt.1. ) then
c
c initial conditions for fitting routine (including swim direction)
c
ILIST(1) = itk*ISWMFORW
ILIST(2) = ktk*ISWMBACK
Nreg = 1
call VTXCRSS(Nreg,ILIST,XINT,IERRCR,TANGENT)
c don't look for vertex if VTXCRSS doesn't find any intersection nor a
c "tangency within cuts" (see VTXCRSS for this definition)
if ((IERRCR.gt.1) .or.
+ (IERRCR.eq.1 .and. .not. TANGENT)) goto 82
c
c Apply a constrained fit if two tangents are parallel - MI sep 97
c ===============
c
if ( TANGENT ) then
CONSTR = TANGENT
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
endif
c
c call minimization procedure
c
NFLAG = 3 ! convergence occurred in loop (3)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
c If minimization procedure converged then put fake Z value (=exclude
c this track end for next vertex searches); if fit is SILVER skip the
c rest of this section, but don't exclude other combinations (VTXMERG
c will eventually select the best one)
if (IFIT.eq.FITGOLD) then
Z1L = 999.
Z2F = -999.
goto 82
elseif (IFIT.eq.FITSILVER) then
goto 82
endif
c
c (4) 1 LAST - 2 LAST
elseif ( D1L2L.lt.RCUTMI .and. Z1L2L.lt.ZCUTMI .and.
+ abs(QFI+QSE).lt.1. ) then
c
c initial conditions for fitting routine (including swim direction)
c
ILIST(1) = itk*ISWMFORW
ILIST(2) = ktk*ISWMFORW
Nreg = 1
call VTXCRSS(Nreg,ILIST,XINT,IERRCR,TANGENT)
c don't look for vertex if VTXCRSS doesn't find any intersection nor a
c "tangency within cuts" (see VTXCRSS for this definition)
if ((IERRCR.gt.1) .or.
+ (IERRCR.eq.1 .and. .not. TANGENT)) goto 82
c
c Apply a constrained fit if two tangents are parallel - MI sep 97
c ===============
c
if ( TANGENT ) then
CONSTR = TANGENT
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
endif
c
c call minimization procedure
c
NFLAG = 3 ! convergence occurred in loop (3)
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
c If minimization procedure converged then put fake Z value (=exclude
c this track end for next vertex searches); if fit is SILVER skip the
c rest of this section, but don't exclude other combinations (VTXMERG
c will eventually select the best one)
if (IFIT.eq.FITGOLD) then
Z1L = 999.
Z2L = -999.
goto 82
elseif (IFIT.eq.FITSILVER) then
goto 82
endif
c
endif
82 continue
85 continue
c
return
end
c
c=======================================================================
INTEGER FUNCTION VTXINIT(iii,TRKDTFS)
c=======================================================================
c
c This function fills the work bank for iii-th DTFS track (having track
c number TRKDTFS) by reading DTFS and by calling the extrapolation
c routine VTXEXTR
c
c Author: Marco Incagli
c =======
c
c Creation date: 18-9-1997
c ==============
c
c Input:
c ======
c iii Track scalar number
c TRKDTFS DTFS track number
c
c Output:
c =======
c VWRK bank is filled for this track
c
c=======================================================================
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:DTFS.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c External declarations
c
integer BOPTN, WBANK, BLOCAT, IWdata, Wsetd
c
c Local declarations
c
integer III, DIMVWRK, ISTAT, TRKDTFS, IND, INDTFS, I
real PHI, COT, SINTHE, COSTHE, Xfh, Yfh, Rfh
character*80 message
REAL Rdc
EQUIVALENCE (TIRGVX(1),Rdc)
c
c=======================================================================
c
c Set initial value to failure
c
VTXINIT = -1
IWdata = INDVWRK + VWRNCO*(iii-1)
c Start filling Work bank (the extrap region is the innermost region
c reached by the particle)
IW(IWdata + VWDTFS) = TRKDTFS ! DTFS trk number
IW(IWdata + VWREG) = NFH ! Extrapolation region
c
c Read DTFS and continue filling work bank
c
ISTAT = BLOCAT(IW,'DTFS',TRKDTFS,IND,INDTFS)
INDTFS = INDTFS + IW(INDTFS+DTFHDS)
c
c First hit position, tangent, Q/Pt and track length
c
RW(IWdata + VWXFH) = RW(INDTFS + DTFXCO)
RW(IWdata + VWYFH) = RW(INDTFS + DTFYCO)
RW(IWdata + VWZFH) = RW(INDTFS + DTFZCO)
Rfh = sqrt( RW(IWdata + VWXFH)**2 + RW(IWdata + VWYFH)**2 )
if (Rfh.lt.RDC .or. Rfh.gt.Rcal) return
phi = RW(INDTFS + DTFPHI)
cot = RW(INDTFS + DTFCTT)
sinthe = 1./sqrt(1.+cot**2)
costhe = SIGN(sqrt(1.-sinthe**2),cot)
RW(IWdata + VWTXFH) = sinthe * cos(phi)
RW(IWdata + VWTYFH) = sinthe * sin(phi)
RW(IWdata + VWTZFH) = costhe
RW(IWdata + VWCUFH) = RW(INDTFS + DTFCUR)
RW(IWdata + VWTLFH) = RW(INDTFS + DTFLNG)
c
c Last hit position, tangent, Q/Pt and track length
c
RW(IWdata + VWXLH) = RW(INDTFS + DTFXLP)
RW(IWdata + VWYLH) = RW(INDTFS + DTFYLP)
RW(IWdata + VWZLH) = RW(INDTFS + DTFZLP)
phi = RW(INDTFS + DTFPLP)
cot = RW(INDTFS + DTFTLP)
sinthe = 1./sqrt(1.+cot**2)
costhe = SIGN(sqrt(1.-sinthe**2),cot)
RW(IWdata + VWTXLH) = sinthe * cos(phi)
RW(IWdata + VWTYLH) = sinthe * sin(phi)
RW(IWdata + VWTZLH) = costhe
RW(IWdata + VWCULH) = SIGN(RW(INDTFS+DTFCLP),RW(INDTFS+DTFCUR))
RW(IWdata + VWTLLH) = RW(INDTFS + DTFLNG)
c
c Copy error matrix at FH
c
do i = 1, 15
RW(IWdata+VWERRF+i-1) = RW(INDTFS+DTFCVM+i-1)
enddo
c
c Now extrapolate track IF its First Hit is "close enough" to z axis
c
If ( Rfh .lt. RPCAMAX ) call VTXEXTR(iii)
c
c Set function to success
c
VTXINIT = YESUCC
c
999 return
end
c=======================================================================
INTEGER FUNCTION VTXEXTR(iii)
c=======================================================================
c
c This function extrapolates track iii toward the interaction point,
c calling the appropriate routine to evaluate the dE/dx when it crosses
c boundaries. The energy loss is evaluated in the pion mass hypothesis.
c
c Author: Marco Incagli
c =======
c
c Creation date: 18-9-1997
c ==============
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c External declarations
c
INTEGER IWDATA
c
c Local declarations
c
integer FLGDIR, Nsurf, iii
real COSTHE, SINTHE, COTTHE, COSPHI, SINPHI, RHO, Q
real Xc, Yc, D2, S2, Delta, X1, X2, Dist1, Dist2, Dist
real Xdc, Ydc, Txdc, Tydc, Alpha, Arc, Zdc, Tzdc
real Pin, Pout, Rmass, Xbp, Ybp, Zbp, Txbp, Tybp, Tzbp
REAL Rdc, Rbp, Bfield, Y1, Y2, XR(3)
EQUIVALENCE (TIRGVX(1),Rdc),(TIRGVX(2),Rbp)
LOGICAL FlagIP(50) ! common VTXFIT-VTXEXTR
COMMON /VTXCLOSIP/ FlagIP
c
c=======================================================================
c
c -------------------------
c 0 - Initialize parameters
c -------------------------
c
Bfield = 6. ! Magnetic field
Flgdir = -1 ! Backward direction
VTXEXTR = YESUCC
IWDATA = INDVWRK + VWRNCO*(iii-1) ! Data index for track iii in VWRK
FlagIP(iii) = .false. ! track is "distant" from IP (see below)
c
c ---------------------------------------
c 1 - Intersection with DC inner boundary
c ---------------------------------------
c
c Find track circle parameters (center and radius):
c
costhe = RW(IWDATA + VWTZFH)
sinthe = sqrt(1-costhe**2)
cotthe = costhe / sinthe
cosphi = RW(IWdata+VWTXFH) / sinthe
sinphi = RW(IWdata+VWTYFH) / sinthe
rho = CONVERS * Bfield * RW(IWdata+VWCUFH)
Q = SIGN(1.,rho)
rho = 1./rho
Xc = RW(IWdata+VWXFH) + rho*sinphi
Yc = RW(IWdata+VWYFH) - rho*cosphi
c intersect
D2 = Xc**2 + Yc**2
S2 = D2 + Rdc**2 - rho**2
Delta = D2 * Rdc**2 - S2**2 / 4.
If (Delta .lt. 0.) then ! No intersections with DC inner wall
IW(IWDATA + VWREG) = NFH
return
Endif
c
X1 = (Xc*S2 + 2.*Yc*sqrt(Delta)) / (2.*D2)
X2 = (Xc*S2 - 2.*Yc*sqrt(Delta)) / (2.*D2)
if (abs(Yc).lt.1.E-6) then
Y1 = sqrt( Rdc**2 - X1**2 )
Y2 = Y1
else
Y1 = S2/(2.*Yc) - X1*Xc/Yc
Y2 = S2/(2.*Yc) - X2*Xc/Yc
endif
c Select the appropriate intersection as the closest to First Hit
Dist1 = (X1 - RW(IWdata+VWXFH))**2 +
+ (Y1 - RW(IWdata+VWYFH))**2
Dist2 = (X2 - RW(IWdata+VWXFH))**2 +
+ (Y2 - RW(IWdata+VWYFH))**2
If (dist1.lt.Dist2) then
Dist = sqrt(Dist1)
Xdc = X1
Ydc = Y1
else
Dist = sqrt(Dist2)
Xdc = X2
Ydc = Y2
endif
c Now find rotation angle to get "z" and new tangent.
c Note that, since the track is swimmed back, this angle is POSITIVE
c (counterclockwise) for positive tracks and viceversa.
c Since rho is signed, this means that ALPHA HAS THE SAME SIGN AS RHO
c The "-" in the z evaluation has a similar motivation.
Alpha = 2 * ASIN(0.5*Dist/Rho)
c
Arc = rho*Alpha
Zdc = RW(IWdata+VWZFH) - Arc*Cotthe
c
Txdc = ( cosphi*cos(Alpha) - sinphi*sin(alpha) ) * sinthe
Tydc = ( sinphi*cos(Alpha) + cosphi*sin(alpha) ) * sinthe
Tzdc = costhe
c Write quantities in Work bank
IW(IWdata+VWREG) = NVC
RW(IWdata+VWXDC) = Xdc
RW(IWdata+VWYDC) = Ydc
RW(IWdata+VWZDC) = Zdc
RW(IWdata+VWTXDC) = Txdc
RW(IWdata+VWTYDC) = Tydc
RW(IWdata+VWTZDC) = Tzdc
RW(IWdata+VWTLDC) = arc
c And now dE/dx correction
Nsurf = NVC ! DC inner wall
Rmass = PION_MASS
call VTXDEDX(iii,Nsurf,Flgdir,Rmass)
c the following flag will be used in VTXFIT. It flags wether the track
c passes or not "close enough" to the Interaction Point (in particular
c it checks the z-coordinate). These tracks are of interest for the
c interference analysis, so they will be treated with particular care.
c So start finding closest point to origin
xr(1) = 0.
xr(2) = 0.
xr(3) = 0.
call VTXSWM(iii,NVC,xr)
c
c now check its (z) position
c
If ( abs(RW(IWdata+VWZVX)) .lt. ZclosIP ) FlagIP(iii) = .true.
c
c ---------------------------------------
c 2 - Inward intersection with BP
c ---------------------------------------
c
c The method is the same as before, however, due to the curvature variation
c caused by the energy loss, all quantities on transverse plane have
c to be revaluated.
c
c Find track circle parameters (center and radius):
c
cosphi = RW(IWdata+VWTXDC) / sinthe
sinphi = RW(IWdata+VWTYDC) / sinthe
rho = CONVERS * Bfield * RW(IWdata+VWCUDC)
rho = 1./rho
Xc = RW(IWdata+VWXDC) + rho*sinphi
Yc = RW(IWdata+VWYDC) - rho*cosphi
c intersect
D2 = Xc**2 + Yc**2
S2 = D2 + Rbp**2 - rho**2
Delta = D2 * Rbp**2 - S2**2 / 4.
If (Delta .lt. 0.) then ! No intersections with BP
IW(IWDATA + VWREG) = NVC
return
Endif
c
X1 = (Xc*S2 + 2.*Yc*sqrt(Delta)) / (2.*D2)
X2 = (Xc*S2 - 2.*Yc*sqrt(Delta)) / (2.*D2)
if (abs(Yc).lt.1.E-6) then
Y1 = sqrt( Rbp**2 - X1**2 )
Y2 = Y1
else
Y1 = S2/(2.*Yc) - X1*Xc/Yc
Y2 = S2/(2.*Yc) - X2*Xc/Yc
endif
c Select the appropriate intersection as the closest to First Hit
Dist1 = (X1 - RW(IWdata+VWXDC))**2 +
+ (Y1 - RW(IWdata+VWYDC))**2
Dist2 = (X2 - RW(IWdata+VWXDC))**2 +
+ (Y2 - RW(IWdata+VWYDC))**2
If (dist1.lt.Dist2) then
Dist = sqrt(Dist1)
Xbp = X1
Ybp = Y1
else
Dist = sqrt(Dist2)
Xbp = X2
Ybp = Y2
endif
c Now find rotation angle to get "z" and new tangent.
c Note that, since the track is swimmed back, this angle is POSITIVE
c (counterclockwise) for positive tracks and viceversa.
c Since rho is signed, this means that ALPHA HAS THE SAME SIGN AS RHO
c The "-" in the z evaluation has a similar motivation.
Alpha = 2 * ASIN(0.5*Dist/Rho)
c
Arc = rho*Alpha
Zbp = RW(IWdata+VWZDC) - Arc*Cotthe
c
Txbp = ( cosphi*cos(Alpha) - sinphi*sin(alpha) ) * sinthe
Tybp = ( sinphi*cos(Alpha) + cosphi*sin(alpha) ) * sinthe
Tzbp = costhe
c Write quantities in Work bank
IW(IWdata+VWREG) = NBI
RW(IWdata+VWXBI) = Xbp
RW(IWdata+VWYBI) = Ybp
RW(IWdata+VWZBI) = Zbp
RW(IWdata+VWTXBI) = Txbp
RW(IWdata+VWTYBI) = Tybp
RW(IWdata+VWTZBI) = Tzbp
RW(IWdata+VWTLBI) = arc
c And now dE/dx correction
Nsurf = NBI ! BP wall entrance
Rmass = PION_MASS
call VTXDEDX(iii,Nsurf,Flgdir,Rmass)
c
c ---------------------------------------
c 3 - Outward intersection with BP
c ---------------------------------------
c
c It may happen that the track crosses twice the Beam Pipe, so it is
c nevessary to evaluate also the second intersection of the track with
c the Beam Pipe (if one exist, than also a second one does). Since we
c are tracing back, this intersection is called the "outward intersection"
c even if, when the vertex is in the VC region and the track crosses
c twice the BP, the "outward intersection" corresponds to the charged
c particle "entering" the BP.
c As before, transverse quantities need to be re-evaluated.
c
c Find track circle parameters (center and radius):
c
cosphi = RW(IWdata+VWTXBI) / sinthe
sinphi = RW(IWdata+VWTYBI) / sinthe
rho = CONVERS * Bfield * RW(IWdata+VWCUBI)
rho = 1./rho
Xc = RW(IWdata+VWXBI) + rho*sinphi
Yc = RW(IWdata+VWYBI) - rho*cosphi
c intersect
D2 = Xc**2 + Yc**2
S2 = D2 + Rbp**2 - rho**2
Delta = D2 * Rbp**2 - S2**2 / 4.
If (Delta .lt. 0.) then !This shouldn't happen; left as protection
IW(IWDATA + VWREG) = NBI
return
Endif
c
X1 = (Xc*S2 + 2.*Yc*sqrt(Delta)) / (2.*D2)
X2 = (Xc*S2 - 2.*Yc*sqrt(Delta)) / (2.*D2)
if (abs(Yc).lt.1.E-6) then
Y1 = sqrt( Rbp**2 - X1**2 )
Y2 = Y1
else
Y1 = S2/(2.*Yc) - X1*Xc/Yc
Y2 = S2/(2.*Yc) - X2*Xc/Yc
endif
c Select the appropriate intersection as the FAREST to First Hit
Dist1 = (X1 - RW(IWdata+VWXBI))**2 +
+ (Y1 - RW(IWdata+VWYBI))**2
Dist2 = (X2 - RW(IWdata+VWXBI))**2 +
+ (Y2 - RW(IWdata+VWYBI))**2
If (Dist1.gt.Dist2) then
Dist = sqrt(Dist1)
Xbp = X1
Ybp = Y1
else
Dist = sqrt(Dist2)
Xbp = X2
Ybp = Y2
endif
c Now find rotation angle to get "z" and new tangent.
c Note that, since the track is swimmed back, this angle is POSITIVE
c (counterclockwise) for positive tracks and viceversa.
c Since rho is signed, this means that ALPHA HAS THE SAME SIGN AS RHO
c The "-" in the z evaluation has a similar motivation.
Alpha = 2 * ASIN(0.5*Dist/Rho)
c
Arc = rho*Alpha
Zbp = RW(IWdata+VWZBI) - Arc*Cotthe
c
Txbp = ( cosphi*cos(Alpha) - sinphi*sin(alpha) ) * sinthe
Tybp = ( sinphi*cos(Alpha) + cosphi*sin(alpha) ) * sinthe
Tzbp = costhe
c Write quantities in Work bank
IW(IWdata+VWREG) = NBU
RW(IWdata+VWXBU) = Xbp
RW(IWdata+VWYBU) = Ybp
RW(IWdata+VWZBU) = Zbp
RW(IWdata+VWTXBU) = Txbp
RW(IWdata+VWTYBU) = Tybp
RW(IWdata+VWTZBU) = Tzbp
RW(IWdata+VWTLBU) = arc
c And now dE/dx correction
Nsurf = NBU ! BP wall at exit
Rmass = PION_MASS
call VTXDEDX(iii,Nsurf,Flgdir,Rmass)
c
c ------------------------------------------------
c End of extrapolation procedure
c ------------------------------------------------
c
return
end
c=======================================================================
SUBROUTINE VTXDEDX(iii,Nsurf,Flgdir,Rmass)
c=======================================================================
c
c This routine evaluates dE/dx loss due to surface crossing
c
c Author: Marco Incagli
c =======
c
c Creation date: 18-9-1997
c ==============
c
c Input data:
c ===========
c iii Track scalar number (not DTFS)
c NSurf Surface number (2=DC; 3=BP in; 4=BP out)
c Flgdir -1 for backward swimming ; otherwise +1
c RMass Particle mass (normally pion)
c
c Output data:
c ============
c RW(...+VWCUxx) the new curvature (defined as Q/Pt) of track
c iii at surface Nsurf is evaluated
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
INTEGER IWdata
INTEGER i, iii, nsurf, flgdir
REAL Rdc, Rbp, Rmass, Pin, Ptin, Pfin, Ptfin, dgdx, Q, sinthe
REAL C1Epxy, C2Epxy, C1Al, C2Al, C1Fe, C2Fe, C1Cu, C2Cu
REAL C1Pb, C2Pb, C1Be, C2Be, C1Ag, C2Ag
REAL XEpxydc, XAldc, XBebp
REAL C1,C2, GAM, GAMOUT, Delx, rx, ry, tx, ty, cosgamma
REAL dgdx1,dgdx2,dgdx3,dgdx4
CHARACTER*80 Message
c
c Bethe-Block constants C1, C2 for seven materials are stored
c in CTIDRI include file (C-epoxy,Al,Fe,Cu,Ag,Pb,Be).
c Not all of them are of interest for KLOE, however.
c C1 = 4 * pi * N_av * rad_el**2 * m_el * c**2 * Z / A (Gev*cm**2/gr)
c C2 = 2 * m_el * c**2 / I (I=mean excitation energy - adimensional)
c
EQUIVALENCE
+(TIDEVX(1,1),C1Epxy),(TIDEVX(2,1),C2Epxy),
+(TIDEVX(1,2),C1Al), (TIDEVX(2,2),C2Al),
+(TIDEVX(1,3),C1Fe), (TIDEVX(2,3),C2Fe),
+(TIDEVX(1,4),C1Cu), (TIDEVX(2,4),C2Cu),
+(TIDEVX(1,5),C1Ag), (TIDEVX(2,5),C2Ag),
+(TIDEVX(1,6),C1Pb), (TIDEVX(2,6),C2Pb),
+(TIDEVX(1,7),C1Be), (TIDEVX(2,7),C2Be)
c
EQUIVALENCE (TIRGVX( 1),Rdc),(TIRGVX( 2),Rbp)
c
c Thickness, in gr/cm**2, of the materials that make up the walls
c
EQUIVALENCE (TIRGVX( 3),XEpxydc),(TIRGVX( 6),XAldc)
EQUIVALENCE (TIRGVX(11),XBebp)
c
c Temporary solution to be corrected as soon as possible!!! (19-feb-1998)
c
real xxxaldc
PARAMETER (XXXAldc=0.11025)
C
c dGAMMA/dx function to be integrated (statement function)
c
dgdx(GAM) = -C1*(GAM**2/(GAM**2-1)*log(C2*(GAM**2-1))-1)
c
c=======================================================================
c
c 1 - DC wall
c ----------
c
If (Nsurf.eq.NVC) then
c get pointer to data and thickness crossed by particle
IWdata = Indvwrk + Vwrnco*(iii-1)
Ptin = 1. / RW(IWdata+VWCUFH) ! Starting Pt in GeV (signed)
Q = SIGN(1.,Ptin)
Pin = Ptin / sqrt(1-RW(IWdata+VWTZFH)**2) ! Starting P in GeV
c get spatial angle momentum-surface
rx = RW(IWdata+VWXDC)
ry = RW(IWdata+VWYDC)
tx = RW(IWdata+VWTXDC)
ty = RW(IWdata+VWTYDC)
cosgamma = abs(rx*tx+ry*ty)/sqrt(rx*rx+ry*ry)
if (cosgamma.lt.0.0001) cosgamma=0.0001
c prepare variables for dE/dx integration (flgdir determines the
c integration sign)
C1 = float(Flgdir) * C1Al/Rmass
C2 = C2Al
GAM = sqrt(Pin**2+Rmass**2)/Rmass
DELX = XAldc/cosgamma
sinthe = sqrt(1-RW(IWdata+VWTZDC)**2)
c
c 2 - BP wall (enter)
c ------------------
c
Elseif (Nsurf.eq.NBI) then
c get pointer to data and thickness crossed by particle
IWdata = Indvwrk + Vwrnco*(iii-1)
Ptin = 1. / RW(IWdata+VWCUDC) ! Starting Pt in GeV (signed)
Q = SIGN(1.,Ptin)
Pin = Ptin / sqrt(1-RW(IWdata+VWTZDC)**2) ! Starting P in GeV
c get spatial angle momentum-surface
rx = RW(IWdata+VWXBI)
ry = RW(IWdata+VWYBI)
tx = RW(IWdata+VWTXBI)
ty = RW(IWdata+VWTYBI)
cosgamma = abs(rx*tx+ry*ty)/sqrt(rx*rx+ry*ry)
if (cosgamma.lt.0.0001) cosgamma=0.0001
c prepare variables for dE/dx integration (flgdir determines the
c integration sign)
C1 = float(Flgdir) * C1Be/Rmass
C2 = C2Be
GAM = sqrt(Pin**2+Rmass**2)/Rmass
DELX = XBebp/cosgamma
sinthe = sqrt(1-RW(IWdata+VWTZBI)**2)
c
c 3 - BP wall (exit)
c ------------------
c
Elseif (Nsurf.eq.NBU) then
c get pointer to data and thickness crossed by particle
IWdata = Indvwrk + Vwrnco*(iii-1)
Ptin = 1. / RW(IWdata+VWCUBI) ! Starting Pt in GeV (signed)
Q = SIGN(1.,Ptin)
Pin = Ptin / sqrt(1-RW(IWdata+VWTZBI)**2) ! Starting P in GeV
c get spatial angle momentum-surface
rx = RW(IWdata+VWXBU)
ry = RW(IWdata+VWYBU)
tx = RW(IWdata+VWTXBU)
ty = RW(IWdata+VWTYBU)
cosgamma = abs(rx*tx+ry*ty)/sqrt(rx*rx+ry*ry)
if (cosgamma.lt.0.0001) cosgamma=0.0001
c prepare variables for dE/dx integration (flgdir determines the
c integration sign)
C1 = float(Flgdir) * C1Be/Rmass
C2 = C2Be
GAM = sqrt(Pin**2+Rmass**2)/Rmass
DELX = XBebp/cosgamma
sinthe = sqrt(1-RW(IWdata+VWTZBU)**2)
c
Else
message = 'Non exixting surface'
call ERLOGR ('VTXDEDX',ERWARN,0,0,message)
Endif
c
c =========================================================================
c Now perform integration of dGdx function using the 4th order Runge-Kutta
c method (see "Numerical Recipes" by Press,Flannery,Teukolsky,Vetterling,
c 1986, pag.550 - Cambridge University Press)
c ==========================================================================
c
dgdx1 = DELX * DGDX(GAM)
dgdx2 = DELX * DGDX(GAM+dgdx1/2.)
dgdx3 = DELX * DGDX(GAM+dgdx2/2.)
dgdx4 = DELX * DGDX(GAM+dgdx3)
c
GAMOUT = GAM + dgdx1/6. + dgdx2/3. + dgdx3/3. + dgdx4/6.
if (GAMOUT.lt.1.) then
message = 'Gamma .lt. 1 !'
call ERLOGR ('VTXDEDX',ERWARN,0,0,message)
GAMOUT = GAM
endif
c
Pfin = Rmass * sqrt(GAMOUT**2-1.) ! Final momentum
Ptfin = Pfin * sinthe ! Final Pt (unsigned)
c
c Write out informations
c
if (Nsurf.eq.2) then
RW(IWdata+VWCUDC) = Q/Ptfin
elseif (Nsurf.eq.3) then
RW(IWdata+VWCUBI) = Q/Ptfin
elseif (Nsurf.eq.4) then
RW(IWdata+VWCUBU) = Q/Ptfin
endif
c
return
end
c
c=======================================================================
subroutine vtxmerg
c=======================================================================
c
c final part of vertex fitting alorythm: the candidate vertices are
c compared in order to:
c 1 - select the best in chi2 if 2 of them are within the bp and use
c the same tracks
c 2 - merge 2 tracks vertices and build up a 3 or 4 tracks vertex if
c two vertices have a track in common .or. if they are in region
c NDC and are close (i.e. four tracks vertex split in 2 two tracks
c vertices)
c 3 - Order the DVFS banks (=remove holes)
c 4 - try to assign the left over particles - to the DVFS banks obtained
c (no new vertex at this stage)
c
c author: marco incagli
c =======
c
c creation date: 3-10-1997
c ==============
c
c=======================================================================
c
c include files
c
$$implicit none
$$include 'k$inc:bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$include 'k$itrk:ctidri.inc'
$$include 'k$itrk:cconst.inc'
$$include 'K$ITRK:vwrk.inc'
$$include 'K$ITRK:vtxfin.inc'
$$include 'K$ITRK:dvcs.inc'
c
integer blocat, ind, inddat, idummy, ii1, ii2, ii3, ii4
integer status, canvtx, i1, i2, vtxctov, ntrfit, ifit
integer fltr(2,50), flvtx(50), idv(50), ilist(10), thisdtfs
integer dhds, dtrp, dtrb, indvcs(50), nflag, vtxflg(50)
integer ndtfs1(50), ndtfs2(50), flend1(50), flend2(50)
integer IWDATA
real c2dvcs(50), vxcan(50), vycan(50), vzcan(50)
real xcn(3), xint(3), chisq, qpt, px1, py1, pz1, ee1
real px2, py2, pz2, ee2, rmass, dist
character*80 message
logical constr
c
c=======================================================================
c
c number of candidate vertices
c
canvtx = nvxfit
c
c Get informations from each vertex
c
do i1 = 1, canvtx
IDV(i1) = i1 ! Pointer to vertex candidate
STATUS = BLOCAT(IW,'DVCS',i1,IND,INDDAT)
DHDS = IW(INDDAT+DVCHDS) ! header size
DTRP = IW(INDDAT+DVCTRK) ! pointer to track data
DTRB = IW(INDDAT+DVCBLK) ! track data block size
C2DVCS(i1) = RW(INDDAT+DHDS+DVCCSQ)! track Chi squared
c decrease track CHI2 if two particles invariant mass is close to Kl mass
Qpt = RW(INDDAT+DHDS+DTRP+DVCQPT)
Pz1 = RW(INDDAT+DHDS+DTRP+DVCCOT)/abs(Qpt)
Px1 = cos(RW(INDDAT+DHDS+DTRP+DVCPHI))/abs(Qpt)
Py1 = sin(RW(INDDAT+DHDS+DTRP+DVCPHI))/abs(Qpt)
Ee1 = sqrt( PION_MASS**2 + Px1**2 + Py1**2 + Pz1**2 )
Qpt = RW(INDDAT+DHDS+DTRP+DTRB+DVCQPT)
Pz2 = RW(INDDAT+DHDS+DTRP+DTRB+DVCCOT)/abs(Qpt)
Px2 = cos(RW(INDDAT+DHDS+DTRP+DTRB+DVCPHI))/abs(Qpt)
Py2 = sin(RW(INDDAT+DHDS+DTRP+DTRB+DVCPHI))/abs(Qpt)
Ee2 = sqrt( PION_MASS**2 + Px2**2 + Py2**2 + Pz2**2 )
Rmass = (Ee1+Ee2)**2-(Px1+Px2)**2-(Py1+Py2)**2-(Pz1+Pz2)**2
Rmass = sqrt(Rmass)
If (abs(Rmass-KAON_MASS).lt.0.05) C2DVCS(i1) = C2DVCS(i1)/2.
c
INDVCS(i1) = INDDAT + DHDS! index to data
VTXFLG(i1) = IW(INDDAT + DHDS + DVCVFL) ! Vtx flag (VTXFIN structure)
VXCAN(i1) = RW(INDDAT+DHDS+DVCVX) ! Vertex x coordinate
VYCAN(i1) = RW(INDDAT+DHDS+DVCVY) ! Vertex y coordinate
VZCAN(i1) = RW(INDDAT+DHDS+DVCVZ) ! Vertex z coordinate
NDTFS1(i1) = IW(INDDAT+DHDS+DTRP+DVCTRN) ! track number
FLEND1(i1) = IW(INDDAT+DHDS+DTRP+DVCFFL) ! "end" flag (FH vs LH)
NDTFS2(i1) = IW(INDDAT+DHDS+DTRP+DTRB+DVCTRN) ! as before for
FLEND2(i1) = IW(INDDAT+DHDS+DTRP+DTRB+DVCFFL) ! second track
c Set initial flags for vertex and tracks.
c For vertices FLVTX=1 means "good", while FLVTX=0 means "bad"
c for tracks FLTR( Track_end , DTFS_num ) = 1 means that this end of
c the track (First or Last hit) is used in more than one vertex candidate.
c Start by assuming vertices good and tracks used in just one vertex.
FLVTX(i1) = 1
ii1 = (FLEND1(i1)+3)/2 ! 1 if FH; 2 if LH
ii2 = (FLEND2(i1)+3)/2 ! 1 if FH; 2 if LH
FLTR(ii1,NDTFS1(i1)) = 0 ! flag for end ii1 of track 1
FLTR(ii2,NDTFS2(i1)) = 0 ! flag for end ii2 of track 2
enddo
c
c Sort vertices according to CHI2
c
do 32 i1 = 1, canvtx-1
do 33 i2 = i1+1, canvtx
if ( C2DVCS(IDV(i1)) .gt. C2DVCS(IDV(i2)) ) then !switch tracks
IDUMMY = IDV(i1)
IDV(i1) = IDV(i2)
IDV(i2) = IDUMMY
endif
33 enddo
32 enddo
c
c ---------- Start loop on vertices -------------
c
do 35 i1 = 1, canvtx-1
c First vertex MUST have both tracks NOT FLAGGED
II1 = (FLEND1(IDV(i1))+3)/2 ! this is 1 for FH and 2 for LH
II2 = NDTFS1(IDV(i1))
II3 = (FLEND2(IDV(i1))+3)/2
II4 = NDTFS2(IDV(i1))
if ( FLTR(ii1,ii2)+FLTR(ii3,ii4) .gt. 0) goto 35
do 36 i2 = i1+1 , canvtx
c check if this is the same vertex found in two different parts of the
c algorithm if this is the case kill second vertex
if ( ( (ndtfs1(IDV(i1))+ndtfs1(IDV(i2))).eq.
& (ndtfs2(IDV(i1))+ndtfs2(IDV(i2))) ).and.
& ( (ndtfs1(IDV(i1))-ndtfs1(IDV(i2))).eq.
& (ndtfs2(IDV(i2))-ndtfs2(IDV(i1))) ) ) then
c flag bad vertex and release tracks 1 and 2 of vertex 2
FLVTX(IDV(i2)) = 0 ! flag for bad vertex
II1 = (FLEND1(IDV(i2))+3)/2
II2 = (FLEND2(IDV(i2))+3)/2
FLTR(II1,NDTFS1(IDV(i2))) = 0 ! flag for "free" track
FLTR(II2,NDTFS2(IDV(i2))) = 0 ! flag for "free" track
c flag good vertex and "used" tracks (tracks 1 and 2 of vertex 1)
II1 = (FLEND1(IDV(i1))+3)/2
II2 = (FLEND2(IDV(i1))+3)/2
FLVTX(IDV(i1)) = 1 ! flag for good vertex
FLTR(ii1,NDTFS1(IDV(i1))) = 1 ! flag for tracks in common
FLTR(ii2,NDTFS2(IDV(i1))) = 1 ! with more VTCS vertices
goto 36
endif
c Second vertex MUST have both tracks NOT FLAGGED
II1 = (FLEND1(IDV(i2))+3)/2 ! this is 1 for FH and 2 for LH
II2 = NDTFS1(IDV(i2))
II3 = (FLEND2(IDV(i2))+3)/2
II4 = NDTFS2(IDV(i2))
if ( FLTR(ii1,ii2)+FLTR(ii3,ii4) .gt. 0) then
FLVTX(IDV(i2)) = 0
goto 36
endif
dist = sqrt( (vxcan(IDV(i1))-vxcan(IDV(i2)))**2 +
+ (vycan(IDV(i1))-vycan(IDV(i2)))**2 +
+ (vzcan(IDV(i1))-vzcan(IDV(i2)))**2 )
c --- Step 1 --- Check if there is a track in common
c (FLEND=-1 for FH and +1 for LH) between two vertices.
c Start with track 2 of vtx 2.
if ( NDTFS1(IDV(i1))*FLEND1(IDV(i1)) .eq.
+ NDTFS2(IDV(i2))*FLEND2(IDV(i2)) .or.
+ NDTFS2(IDV(i1))*FLEND2(IDV(i1)) .eq.
+ NDTFS2(IDV(i2))*FLEND2(IDV(i2)) ) then
c call fitting routine with three tracks; do a constrained fit
c note that the minimization routine needs the "track number" and
c not the "DTFS number"; this is the reason of the following loop
ilist(1) = 0
ilist(2) = 0
ilist(3) = 0
do ii1=1, tottrk
thisdtfs = IW(INDVWRK+VWRNCO*(ii1-1)+VWDTFS)
if ( thisdtfs.eq.NDTFS1(IDV(i1)) ) then
Ilist(1) = ii1*FLEND1(IDV(i1))
elseif ( thisdtfs.eq.NDTFS2(IDV(i1)) ) then
Ilist(2) = ii1*FLEND2(IDV(i1))
elseif ( thisdtfs.eq.NDTFS1(IDV(i2)) ) then
Ilist(3) = ii1*FLEND1(IDV(i2))
endif
enddo
c check above assignement
if ( Ilist(1)*Ilist(2)*Ilist(3) .eq. 0 ) then
message =
+ 'Error in DTFS number assignement for 3 tracks vertex'
call ERLOGR ('VTXMERG',ERWARN,0,0,message)
CHISQ = 1000.
goto 998
endif
c
ntrfit = 3
CONSTR = .true.
XINT(1) = VXCAN(IDV(i1))
XINT(2) = VYCAN(IDV(i1))
XINT(3) = VZCAN(IDV(i1))
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
c
c call minimization procedure
c
NFLAG = 4 ! convergence occurred in merging procedure
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
998 continue
c
c check if three is better than two by looking at REDUCED CHI2
c
if ( IFIT.eq.FITGOLD .and.
+ CHISQ/3..lt.C2DVCS(IDV(i1))/2. ) then
c write DVFS bank by merging DVCS banks i1 and i2 using the two tracks
c of bank i1 and track number 1 of bank i2
ii1 = IDV(i1)
ii2 = IDV(i2)
call VTXMRGCAN(ii1,ii2,1)
c flag the two original vertices as "bad"
FLVTX(IDV(i1)) = 0 ! flag for bad vertex
FLVTX(IDV(i2)) = 0 ! flag for bad vertex
else
c Here there is one track in common, but second vertex is bad (large
c CHI2); flag bad vertex and release track 1 of vertex 2
FLVTX(IDV(i2)) = 0 ! flag for bad vertex
II1 = (FLEND1(IDV(i2))+3)/2
FLTR(II1,NDTFS1(IDV(i2))) = 0 ! flag for "free" track
c
c change also flag inside work bank for this end of the track
c
IWDATA = INDVWRK+VWRNCO*(abs(Ilist(3))-1)+II1+1
IW(IWDATA) = 0
c flag good vertex and "used" tracks (tracks 1 and 2 of vertex 1)
II1 = (FLEND1(IDV(i1))+3)/2
II2 = (FLEND2(IDV(i1))+3)/2
FLVTX(IDV(i1)) = 1 ! flag for good vertex
FLTR(ii1,NDTFS1(IDV(i1))) = 1 ! flag for tracks in common
FLTR(ii2,NDTFS2(IDV(i1))) = 1 ! with more VTCS vertices
endif
c --- Step 2 --- Check if there is a track in common
c (FLEND=-1 for FH and +1 for LH) between two vertices.
c Now check track 1 of vtx 2.
elseif ( NDTFS1(IDV(i1))*FLEND1(IDV(i1)) .eq.
+ NDTFS1(IDV(i2))*FLEND1(IDV(i2)) .or.
+ NDTFS2(IDV(i1))*FLEND2(IDV(i1)) .eq.
+ NDTFS1(IDV(i2))*FLEND1(IDV(i2)) ) then
c call fitting routine with three tracks; do a constrained fit ;
c note that the minimization routine needs the "track number" and
c not the "DTFS number"; this is the reason for the following loop
ilist(1) = 0
ilist(2) = 0
ilist(3) = 0
do ii1=1, tottrk
thisdtfs = IW(INDVWRK+VWRNCO*(ii1-1)+VWDTFS)
if ( thisdtfs.eq.NDTFS1(IDV(i1)) ) then
Ilist(1) = ii1*FLEND1(IDV(i1))
elseif ( thisdtfs.eq.NDTFS2(IDV(i1)) ) then
Ilist(2) = ii1*FLEND2(IDV(i1))
elseif ( thisdtfs.eq.NDTFS2(IDV(i2)) ) then
Ilist(3) = ii1*FLEND2(IDV(i2))
endif
enddo
c check above assignement
if ( Ilist(1)*Ilist(2)*Ilist(3) .eq. 0 ) then
message =
+ 'Error in DTFS number assignement for 3 tracks vertex'
call ERLOGR ('VTXMERG',ERWARN,0,0,message)
CHISQ = 1000.
goto 997
endif
c
ntrfit = 3
CONSTR = .true.
XINT(1) = VXCAN(IDV(i1))
XINT(2) = VYCAN(IDV(i1))
XINT(3) = VZCAN(IDV(i1))
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
c
c call minimization procedure
c
NFLAG = 4 ! convergence occurred in merging procedure
CALL VTXMIN(Nflag,Ntrfit,ILIST,XINT,CHISQ,IFIT,XCN,CONSTR)
CONSTR = .false.
997 continue
c
c check if three is better than two by looking at REDUCED CHI2
c
if ( IFIT.eq.FITGOLD .and.
+ CHISQ/3..lt.C2DVCS(IDV(i1))/2. ) then
c write DVFS bank by merging DVCS banks i1 and i2 using the two tracks
c of bank i1 and track number 2 of bank i2
ii1 = IDV(i1)
ii2 = IDV(i2)
call VTXMRGCAN(ii1,ii2,2)
c flag the two original vertices as "bad";
FLVTX(IDV(i1)) = 0 ! flag for bad vertex
FLVTX(IDV(i2)) = 0 ! flag for bad vertex
else
c Here there is one track in common, but second vertex is bad (large
c CHI2); flag bad vertex and release track 2 of vertex 2
FLVTX(IDV(i2)) = 0 ! flag for bad vertex
II1 = (FLEND2(IDV(i2))+3)/2
FLTR(ii1,NDTFS2(IDV(i2))) = 0 ! flag for "free" track
c
c change also flag inside work bank for this end of the track
c
IWDATA = INDVWRK+VWRNCO*(abs(Ilist(3))-1)+II1+1
IW(IWDATA) = 0
c flag good vertex and "used" tracks (tracks 1 and 2 of vertex 1)
II1 = (FLEND1(IDV(i1))+3)/2
II2 = (FLEND2(IDV(i1))+3)/2
FLVTX(IDV(i1)) = 1 ! flag for good vertex
FLTR(ii1,NDTFS1(IDV(i1))) = 1 ! flag for tracks used in more
FLTR(ii2,NDTFS2(IDV(i1))) = 1 ! than one candidate vertex
endif
c --- Step 3 --- Check if the two vertices are inside DC .and.
c if they are close. When this happens try to build up a 4 tracks vertex.
c
elseif ( ( vtxflg(IDV(i1)).eq.3.and.vtxflg(IDV(i2)).eq.3 ) .and.
+ dist .lt. distvtx ) then
c Call fitting routine with four tracks and do a constrained fit.
c Note that the minimization routine needs the "track number" and
c not the "DTFS number"; this is the reason for the following loop
ilist(1) = 0
ilist(2) = 0
ilist(3) = 0
ilist(4) = 0
do ii1=1, tottrk
thisdtfs = IW(INDVWRK+VWRNCO*(ii1-1)+VWDTFS)
if ( thisdtfs.eq.NDTFS1(IDV(i1)) ) then
Ilist(1) = ii1*FLEND1(IDV(i1))
elseif ( thisdtfs.eq.NDTFS2(IDV(i1)) ) then
Ilist(2) = ii1*FLEND2(IDV(i1))
elseif ( thisdtfs.eq.NDTFS1(IDV(i2)) ) then
Ilist(3) = ii1*FLEND1(IDV(i2))
elseif ( thisdtfs.eq.NDTFS2(IDV(i2)) ) then
Ilist(4) = ii1*FLEND2(IDV(i2))
endif
enddo
c check above assignement
if ( Ilist(1)*Ilist(2)*Ilist(3)*Ilist(4) .eq. 0 ) then
message =
+ 'Error in DTFS number assignement for 4 tracks vertex'
call ERLOGR ('VTXMERG',ERWARN,0,0,message)
goto 779
endif
c Perform a constrained 4 tracks fit and compare CHI2s
ntrfit = 4
CONSTR = .true.
XINT(1) = ( VXCAN(IDV(i1)) + VXCAN(IDV(i2)) ) / 2.
XINT(2) = ( VYCAN(IDV(i1)) + VYCAN(IDV(i2)) ) / 2.
XINT(3) = ( VZCAN(IDV(i1)) + VZCAN(IDV(i2)) ) / 2.
XCN(1) = XINT(1)
XCN(2) = XINT(2)
XCN(3) = XINT(3)
c
c call minimization procedure
c
NFLAG = 4 ! convergence occurred in merging procedure
CALL VTXMIN(Nflag, Ntrfit, ILIST,
+ XINT, CHISQ, IFIT, XCN, CONSTR)
CONSTR = .false.
c
c check if FOUR is better than TWO
c
if ( IFIT.eq.FITGOLD .and.
+ CHISQ .lt. (C2DVCS(IDV(i1))+C2DVCS(IDV(i2))) ) then
c write DVFS bank by merging DVCS banks i1 and i2 using the two tracks
c of bank i1 and track number 2 of bank i2
ii1 = IDV(i1)
ii2 = IDV(i2)
call VTXMRGCAN(ii1,ii2,2)
c flag the two original vertices as "bad";
FLVTX(IDV(i1)) = 0 ! flag for bad vertex
FLVTX(IDV(i2)) = 0 ! flag for bad vertex
endif
c
779 continue
endif
c
36 enddo
35 enddo
c
c 2 - loop on DVCS bank and write out DVFS bank for "good vertices"
c
do i1 = 1, Canvtx
if (FLVTX(IDV(i1)) .eq. 1) then
STATUS = VTXCTOV(IDV(i1))
IF (STATUS.NE.YESUCC) THEN
message = 'Error in DVCS_to_DVFS routine'
call ERLOGR ('VTXMERG',ERWARN,0,0,message)
return
ENDIF
endif
enddo
c
c 3 - order the DVFS banks (=remove holes)
c
i1 = 0
STATUS = BLOCAT(IW,'DVFS',-1,IND,INDDAT)
if (STATUS .ne. YESUCC) return ! no DVFS bank present
c
10 i1 = i1 + 1
c change bank number
IW(IND-2) = i1
c jump to next bank
IND = IW(IND-1)
c stop if IND = 0 (last named bank)
if (IND.eq.0) then
NVXFIT = i1
goto 11
endif
c
goto 10
11 continue
c
c 4 - add tracks to vertices
c
c 4.1) loop on vertices
c
c
return
end
c
c=======================================================================
INTEGER FUNCTION VTXMRGCAN(i1,i2,itrk)
c=======================================================================
c
c Routine that merges DVCS banks number i1 and i2 into a DVFS bank.
c If itrk = 1 .or. itrk = 2 then the DVFS bank has three tracks:
c the 2 tracks of bank i1 and track "itrk" from bank i2.
c Else if itrk = 3, than bank DVFS is composed by four tracks, the two
c from bank i1 and the 2 from bank i2.
c
c Author: Marco Incagli
c =======
c
c Creation date: 25-11-1997
c ==============
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
$$INCLUDE 'K$ITRK:DVCS.INC'
$$INCLUDE 'K$ITRK:DVFS.INC'
c
INTEGER STATUS, iii, itrk, IND, INDVCS, INDVFS
INTEGER i1, i2, i3
INTEGER BNKSIZ, BLOCAT, BTYMKG, Idvfs, Idvcs, TRKPOI, BLKSIZ
REAL tx, ty
CHARACTER*80 message
c
c=======================================================================
c
c set status to failure
c
VTXMRGCAN = -1
c
c Open the Vertex Candidate bank DVCS i1
c
STATUS = BLOCAT(IW,'DVCS',i1,IND,INDVCS)
IF (STATUS.NE.YESUCC) THEN
message = 'Error in locating bank DVCS'
call ERLOGR ('VTXCTOV',ERWARN,0,0,message)
return
ENDIF
TRKPOI = IW(INDVCS+DVCTRK)
BLKSIZ = IW(INDVCS+DVCBLK)
INDVCS = INDVCS + IW( INDVCS+DVCHDS )
c
c create the corresponding Vertex bank DVFS
c at the moment only 3 or 4 tracks banks can be created
c
if (itrk.lt.3) then
STATUS = BTYMKG(IW,'DVFS',i1,'3I4,10R4,3I4,3(I4,16R4,2I4)',
+ BNKSIZ, IND, INDVFS)
else
STATUS = BTYMKG(IW,'DVFS',i1,'3I4,10R4,3I4,4(I4,16R4,2I4)',
+ BNKSIZ, IND, INDVFS)
endif
c
IF (STATUS.NE.YESUCC) THEN
message = 'Error in creating bank DVFS'
call ERLOGR ('VTXCTOV',ERWARN,0,0,message)
return
ENDIF
c fill header
IW(INDVFS+DVFHDS) = 2
IW(INDVFS+DVFVRN) = 1
INDVFS = INDVFS + IW(INDVFS+DVFHDS)
c fill body 1 (vertex)
IW(indvfs+dvfqfl) = IW(indvcs+dvcqfl)
RW(indvfs+dvfxco) = RW(indvcs+dvcvx)
RW(indvfs+dvfyco) = RW(indvcs+dvcvy)
RW(indvfs+dvfzco) = RW(indvcs+dvcvz)
RW(indvfs+dvfchi) = RW(indvcs+dvccsq)
do i1=0,5
RW(indvfs+dvfcvv+i1) = RW(indvcs+dvccov+i1)
enddo
IW(indvfs+dvfidf) = IW(indvcs+dvcvfl)
c fill body 2 (connected tracks). Start from tracks from bank i1.
do i1 = 1, 2
idvfs = indvfs + dvfntr + (i1-1)*19 ! pointer to body2 data
idvcs = indvcs + trkpoi + (i1-1)*blksiz ! pointer to body2 data
IW(idvfs+dvftra) = IW(idvcs+dvctrn)
RW(idvfs+dvfcrv) = RW(idvcs+dvcqpt)
RW(idvfs+dvfctv) = RW(idvcs+dvccot)
RW(idvfs+dvfphv) = RW(idvcs+dvcphi)
RW(idvfs+dvfchv) = RW(idvcs+dvctch)
RW(idvfs+dvflnv) = RW(idvcs+dvclen)
c
c now copy covariance matrix after swim; in DVFS this matrix is
c broken into two pieces (the r-z covariance - 3 elements - and
c the Kur-cot-phi covariance - 6 elements)
c
RW(idvfs+dvfc1v+0) = RW(idvcs+dvctcv+0) ! <r,r>
RW(idvfs+dvfc1v+1) = RW(idvcs+dvctcv+1) ! <r,z>
RW(idvfs+dvfc1v+2) = RW(idvcs+dvctcv+6) ! <z,z>
c
do i3 = 0, 5 ! covariance matrix after swim (momentum)
RW(idvfs+dvfcov+i3) = RW(idvcs+dvctcv+9+i3)
enddo
c
c now evaluate D0 and Z0 after fit and store them
c
tx = SIN(RW(idvcs+dvcphi))
ty = COS(RW(idvcs+dvcphi))
RW(idvfs+dvfD0v) = ( RW(idvcs+dvctx)-RW(indvcs+dvcvx) )*ty
& + ( RW(idvcs+dvcty)-RW(indvcs+dvcvy) )*tx
RW(idvfs+dvfdzv) = RW(idvcs+dvctz)-RW(indvcs+dvcvz)
c
c finally store the flag First Hit vs Last Hit in DVFPID word (currently
c not used)
c
IW(idvfs+dvfpid) = IW(idvcs+dvcffl)
c
enddo
c
c Now add track (or tracks) bank i2. Open this bank first.
c
STATUS = BLOCAT(IW,'DVCS',i2,IND,INDVCS)
IF (STATUS.NE.YESUCC) THEN
message = 'Error in locating bank DVCS'
call ERLOGR ('VTXCTOV',ERWARN,0,0,message)
return
ENDIF
TRKPOI = IW(INDVCS+DVCTRK)
BLKSIZ = IW(INDVCS+DVCBLK)
INDVCS = INDVCS + IW( INDVCS+DVCHDS )
c
c if itrk<3 than I have to add just one track from bank DVCS, the one
c corresponding to track # itrk; if itrk>=3, than both DVCS tracks have
c to be added to DVFS bank
c
If (itrk .lt. 3) then ! build a 3 tracks bank
c ===================== =====================
c
IW(indvfs+dvfntr) = 3 ! 3 tracks bank
c
i1 = 3 ! add third track to DVFS
i3 = itrk ! using track "itrk" from DVCS
idvfs = indvfs + dvfntr + (i1-1)*19 ! pointer to body2 data
idvcs = indvcs + trkpoi + (i3-1)*blksiz ! pointer to body2 data
IW(idvfs+dvftra) = IW(idvcs+dvctrn)
RW(idvfs+dvfcrv) = RW(idvcs+dvcqpt)
RW(idvfs+dvfctv) = RW(idvcs+dvccot)
RW(idvfs+dvfphv) = RW(idvcs+dvcphi)
RW(idvfs+dvfchv) = RW(idvcs+dvctch)
RW(idvfs+dvflnv) = RW(idvcs+dvclen)
c
c now copy covariance matrix after swim; in DVFS this matrix is
c broken into two pieces (the r-z covariance - 3 elements - and
c the Kur-cot-phi covariance - 6 elements)
c
RW(idvfs+dvfc1v+0) = RW(idvcs+dvctcv+0) ! <r,r>
RW(idvfs+dvfc1v+1) = RW(idvcs+dvctcv+1) ! <r,z>
RW(idvfs+dvfc1v+2) = RW(idvcs+dvctcv+6) ! <z,z>
c
do i3 = 0, 5 ! covariance matrix after swim (momentum)
RW(idvfs+dvfcov+i3) = RW(idvcs+dvctcv+9+i3)
enddo
c
c now evaluate D0 and Z0 after fit and store them
c
tx = SIN(RW(idvcs+dvcphi))
ty = COS(RW(idvcs+dvcphi))
RW(idvfs+dvfD0v) = ( RW(idvcs+dvctx)-RW(indvcs+dvcvx) )*ty
& + ( RW(idvcs+dvcty)-RW(indvcs+dvcvy) )*tx
RW(idvfs+dvfdzv) = RW(idvcs+dvctz)-RW(indvcs+dvcvz)
c
c finally store the flag First Hit vs Last Hit in DVFPID word (currently
c not used)
c
IW(idvfs+dvfpid) = IW(idvcs+dvcffl)
c
c
else ! build a 4 tracks bank
c ==== =====================
c
IW(indvfs+dvfntr) = 4 ! number of tracks
c
do i1 = 3, 4 ! add tracks 3 and 4
i3 = i1 - 2 ! using tracks 1 and 2 from DVCS
idvfs = indvfs + dvfntr + (i1-1)*19 ! pointer to body2 data
idvcs = indvcs + trkpoi + (i3-1)*blksiz ! pointer to body2 data
IW(idvfs+dvftra) = IW(idvcs+dvctrn)
RW(idvfs+dvfcrv) = RW(idvcs+dvcqpt)
RW(idvfs+dvfctv) = RW(idvcs+dvccot)
RW(idvfs+dvfphv) = RW(idvcs+dvcphi)
RW(idvfs+dvfchv) = RW(idvcs+dvctch)
RW(idvfs+dvflnv) = RW(idvcs+dvclen)
c
c now copy covariance matrix after swim; in DVFS this matrix is
c broken into two pieces (the r-z covariance - 3 elements - and
c the Kur-cot-phi covariance - 6 elements)
c
RW(idvfs+dvfc1v+0) = RW(idvcs+dvctcv+0) ! <r,r>
RW(idvfs+dvfc1v+1) = RW(idvcs+dvctcv+1) ! <r,z>
RW(idvfs+dvfc1v+2) = RW(idvcs+dvctcv+6) ! <z,z>
c
do i3 = 0, 5 ! covariance matrix after swim (momentum)
RW(idvfs+dvfcov+i3) = RW(idvcs+dvctcv+9+i3)
enddo
c
c now evaluate D0 and Z0 after fit and store them
c
tx = SIN(RW(idvcs+dvcphi))
ty = COS(RW(idvcs+dvcphi))
RW(idvfs+dvfD0v) = ( RW(idvcs+dvctx)-RW(indvcs+dvcvx) )*ty
& + ( RW(idvcs+dvcty)-RW(indvcs+dvcvy) )*tx
RW(idvfs+dvfdzv) = RW(idvcs+dvctz)-RW(indvcs+dvcvz)
c
c finally store the flag First Hit vs Last Hit in DVFPID word (currently
c not used)
c
IW(idvfs+dvfpid) = IW(idvcs+dvcffl)
c
enddo
c
endif
c
VTXMRGCAN = YESUCC
c
return
end
c=======================================================================
INTEGER FUNCTION VTXCTOV(iii)
c=======================================================================
c
c routine that translates DVCS into DVFS bank
c
c Author: Marco Incagli
c =======
c
c Creation date: 6-10-1997
c ==============
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
$$INCLUDE 'K$ITRK:DVCS.INC'
$$INCLUDE 'K$ITRK:DVFS.INC'
c
INTEGER STATUS, iii, i1, i2, IND, INDVCS, INDVFS, BLKSIZ
INTEGER BNKSIZ, BLOCAT, BTYMKG, Idvfs, Idvcs, TRKPOI
REAL tx, ty
CHARACTER*80 message
c
c=======================================================================
c
c set status to failure
c
VTXCTOV = -1
c
c Open the Vertex Candidate bank DVCS
c
STATUS = BLOCAT(IW,'DVCS',iii,IND,INDVCS)
IF (STATUS.NE.YESUCC) THEN
message = 'Error in locating bank DVCS'
call ERLOGR ('VTXCTOV',ERWARN,0,0,message)
return
ENDIF
TRKPOI = IW(INDVCS+DVCTRK)
BLKSIZ = IW(INDVCS+DVCBLK)
INDVCS = INDVCS + IW( INDVCS+DVCHDS )
c
c create the corresponding Vertex bank DVFS
c
STATUS = BTYMKG(IW,'DVFS',iii,'3I4,10R4,3I4,2(I4,16R4,2I4)',
+ BNKSIZ, IND, INDVFS)
IF (STATUS.NE.YESUCC) THEN
message = 'Error in creating bank DVFS'
call ERLOGR ('VTXCTOV',ERWARN,0,0,message)
return
ENDIF
c fill header
IW(INDVFS+DVFHDS) = 2
IW(INDVFS+DVFVRN) = 1
INDVFS = INDVFS + IW(INDVFS+DVFHDS)
c fill body 1 (vertex)
IW(indvfs+dvfqfl) = IW(indvcs+dvcqfl)
RW(indvfs+dvfxco) = RW(indvcs+dvcvx)
RW(indvfs+dvfyco) = RW(indvcs+dvcvy)
RW(indvfs+dvfzco) = RW(indvcs+dvcvz)
RW(indvfs+dvfchi) = RW(indvcs+dvccsq)
do i1=0,5
RW(indvfs+dvfcvv+i1) = RW(indvcs+dvccov+i1)
enddo
IW(indvfs+dvfidf) = IW(indvcs+dvcvfl)
IW(indvfs+dvfntr) = IW(indvcs+dvcntr)
c fill body 2 (connected tracks)
do i1 = 1, IW(indvcs+dvcntr)
idvfs = indvfs + dvfntr + (i1-1)*19 ! pointer to body2 data
idvcs = indvcs + trkpoi + (i1-1)*blksiz ! pointer to body2 data
IW(idvfs+dvftra) = IW(idvcs+dvctrn)
RW(idvfs+dvfcrv) = RW(idvcs+dvcqpt)
RW(idvfs+dvfctv) = RW(idvcs+dvccot)
RW(idvfs+dvfphv) = RW(idvcs+dvcphi)
RW(idvfs+dvfchv) = RW(idvcs+dvctch)
RW(idvfs+dvflnv) = RW(idvcs+dvclen)
c
c now copy covariance matrix after swim; in DVFS this matrix is
c broken into two pieces (the r-z covariance - 3 elements - and
c the Kur-cot-phi covariance - 6 elements)
c
RW(idvfs+dvfc1v+0) = RW(idvcs+dvctcv+0) ! <r,r>
RW(idvfs+dvfc1v+1) = RW(idvcs+dvctcv+1) ! <r,z>
RW(idvfs+dvfc1v+2) = RW(idvcs+dvctcv+6) ! <z,z>
c
do i2 = 0, 5 ! covariance matrix after swim (momentum)
RW(idvfs+dvfcov+i2) = RW(idvcs+dvctcv+9+i2)
enddo
c
c now evaluate D0 and Z0 after fit and store them
c
tx = SIN(RW(idvcs+dvcphi))
ty = COS(RW(idvcs+dvcphi))
RW(idvfs+dvfD0v) = ( RW(idvcs+dvctx)-RW(indvcs+dvcvx) )*ty
& + ( RW(idvcs+dvcty)-RW(indvcs+dvcvy) )*tx
RW(idvfs+dvfdzv) = RW(idvcs+dvctz)-RW(indvcs+dvcvz)
c
c finally store the flag First Hit vs Last Hit in DVFPID word (currently
c not used)
c
IW(idvfs+dvfpid) = IW(idvcs+dvcffl)
c
enddo
c
VTXCTOV = YESUCC
c
return
end
c=======================================================================
SUBROUTINE VTXCAN
c=======================================================================
c
c This routine stores candidate vertices in the candidate vertex bank
c DVCS. This bank has the same structure as the VERSION 2 of the
c DVFS bank, except that it always assumes a two tracks vertex.
c (In routine VTXMERG, that actually builds DVFS banks, this assumption
c is released).
c
c Author: Marco Incagli
c =======
c
c Creation date: 18-9-1997
c ==============
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
$$INCLUDE 'K$ITRK:DVCS.INC'
c
c External variables
c
INTEGER BTYMKG, IND, INDDAT, STATUS
c
c Local variables
c
INTEGER i, ii, Bnksiz, IWdata, ICdata
REAL sinphi, cosphi, sinthe, costhe
CHARACTER*80 message
c
c common passed to VTXCAN from VTXMIN for creating candidate bank DVCS
c
INTEGER Jfit, Nflg, Jlist(2), Flagfl(2)
REAL Vx, Vy, Vz, mixga(3,3), Vchi2, VTrkcsq(2)
COMMON /DVCS_COM/ Jfit, Nflg, Jlist, Flagfl,
+ Vx, Vy, Vz, mixga, Vchi2, VTrkcsq
c
c=======================================================================
c
c create the Vertex Candidate bank DVCS,
c
STATUS = BTYMKG(IW,'DVCS',Nvxfit,'5I4,10R4,2I4,2(2I4,23R4)',
+ BNKSIZ, IND, INDDAT)
IF (STATUS.NE.YESUCC) THEN
message = 'Error in creating bank DVCS'
call ERLOGR ('VTXCAN',ERWARN,0,0,message)
return
ENDIF
c
c fill bank header (parameters defined in DVCS.INC)
c
IW( INDDAT + DVCHDS ) = DVCHDSIZ ! Header size
IW( INDDAT + DVCVRN ) = DVCVERNUM ! Version number
IW( INDDAT + DVCTRK ) = DVCTRKPOI ! Pointer to connected track data
IW( INDDAT + DVCBLK ) = DVCBLKSIZ ! Size of track informations block
c
c fill bank Body 1 (infos on vertex)
c
INDDAT = INDDAT + IW(INDDAT+DVCHDS)
IW( INDDAT + DVCQFL ) = JFIT
RW( INDDAT + DVCVX ) = Vx
RW( INDDAT + DVCVY ) = Vy
RW( INDDAT + DVCVZ ) = Vz
RW( INDDAT + DVCCSQ ) = Vchi2
RW( INDDAT + DVCCOV + 0 ) = mixga(1,1)
RW( INDDAT + DVCCOV + 1 ) = mixga(1,2)
RW( INDDAT + DVCCOV + 2 ) = mixga(1,3)
RW( INDDAT + DVCCOV + 3 ) = mixga(2,2)
RW( INDDAT + DVCCOV + 4 ) = mixga(2,3)
RW( INDDAT + DVCCOV + 5 ) = mixga(3,3)
IW( INDDAT + DVCVFL ) = Nflg
IW( INDDAT + DVCNTR ) = 2 ! Assume two tracks at this level
c
c Fill bank Body 2 (informations on connected tracks)
c Note that the "track length" is defined as the length of the arc
c delimited by the Vertex and the track end used for extrapolation
c (First Hit or Last Hit)
c
do i = 1, IW(INDDAT+DVCNTR)
IWdata = INDVWRK+VWRNCO*(Jlist(i)-1) ! pointer into VWRK
ICdata = INDDAT + DVCTRKPOI + (i-1)*DVCBLKSIZ ! pointer into DVCS
IW( ICdata + DVCTRN ) = IW( IWdata + VWDTFS )
IW( ICdata + DVCFFL ) = Flagfl(i)
RW( ICdata + DVCTCH ) = VTrkcsq(i)
If ( Flagfl(i) .lt. 0 ) then ! FH of track connected to vertex
RW( ICdata + DVCTX ) = RW( IWdata + VWXVX )
RW( ICdata + DVCTY ) = RW( IWdata + VWYVX )
RW( ICdata + DVCTZ ) = RW( IWdata + VWZVX )
RW( ICdata + DVCQPT ) = RW( IWdata + VWCUVX )
costhe = RW( IWdata + VWTZVX )
sinthe = sqrt( 1 - costhe**2 )
cosphi = RW( IWdata + VWTXVX ) / sinthe
sinphi = RW( IWdata + VWTYVX ) / sinthe
RW( ICdata + DVCPHI ) = ATAN2(sinphi,cosphi)
RW( ICdata + DVCCOT ) = costhe/sinthe
RW( ICdata + DVCLEN ) = RW( IWdata + VWTLDC ) +
+ RW( IWdata + VWTLBU ) + RW( IWdata + VWTLBI ) +
+ RW( IWdata + VWTLVX )
do ii = 0, 14 ! covariance matrix after swim
RW( ICdata + DVCTCV + ii ) = RW( IWdata + VWERRV + ii)
enddo
else ! LH of track connected to vertex
RW( ICdata + DVCTX ) = RW( IWdata + VWXVL )
RW( ICdata + DVCTY ) = RW( IWdata + VWYVL )
RW( ICdata + DVCTZ ) = RW( IWdata + VWZVL )
RW( ICdata + DVCQPT ) = RW( IWdata + VWCUVL )
costhe = RW( IWdata + VWTZVL )
sinthe = sqrt( 1 - costhe**2 )
cosphi = RW( IWdata + VWTXVL ) / sinthe
sinphi = RW( IWdata + VWTYVL ) / sinthe
RW( ICdata + DVCPHI ) = ATAN2(sinphi,cosphi)
RW( ICdata + DVCCOT ) = costhe/sinthe
RW( ICdata + DVCLEN ) = RW( IWdata + VWTLVL )
do ii = 0, 14 ! covariance matrix after swim
RW( ICdata + DVCTCV + ii ) = RW( IWdata + VWERRV + ii)
enddo
endif
c
enddo
c
return
end
c
c=======================================================================
SUBROUTINE VTXCRSS(Nreg,ilist,XINT,IERR,TANGENT)
c=======================================================================
c
c This routine looks for an intersection between two helices and, if
c found (within cuts), it perform a geometrical iterative procedure
c to find the best CANDITATE VERTEX to be fed to the minimization routine
c (this iterative procedure is described below, in this routine).
c This routine has been obtained by "cannibalizing" some sections of
c a previous routine by D.Coppage, from ARGUS, modified by A.Calcaterra,
c that was operating just on the transverse plane.
c
c Author: Marco Incagli
c =======
c
c Creation date: 23-9-1997
c ==============
c
c Input data:
c ===========
c Nreg region of extrapolation (1 is DC ; >1 is VC.or.BP)
c ilist track numbers (with sign! "-" means "backward swim")
c
c Output data:
c ============
c XINT(3) coordinates of intersection
c IERR error code (see below)
c TANGENT flag set when the two circles on x-y plane are "tangent
c within cuts" or, for a positive crossing, when the
c two extrapolated tracks are "parallel within cuts"
c (this last check is used in "constrained fits").
c
c Error codes:
c ============
c IERR = 0 found crossing (check TANGENT for "parallelism")
c IERR = 1 circles on x-y don't intersecate (check TANGENT for
c "tangent within cuts")
c IERR = 2 OK intersection in x-y, but z separation above cuts
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c
c=======================================================================
c
C
REAL PROXIM ! routine from CERNLIB (MI)
C
LOGICAL EXTCRC,INTCRC,tangent
INTEGER LSTRG1,LSTRG2,INRMST,REG,REGION,NDX
INTEGER IL,IG
REAL RRMIN,RR,TRMIV
REAL DCC,RLCOS,ARGMNT,A_HALF
REAL DLTZ,SGNIL,SGNIG,DLTCPL,DLTCMN,ANGNW
REAL VS(3,2),TS(3,2),VH(3),RHO(2), XINT(3)
REAL ANGSTR(2),XCTR(2,2),CTGTHE(2),U(2),ANG(2),ZPL(2),ZMN(2)
REAL SGNQB(2),ANGCTR(2),TURNPL(2),TURNMN(2)
C
CMI
INTEGER LITER, IERR, LTK, LERRS, ilist(2), jlist(2)
INTEGER nreg, Offset
REAL ZSEP, XINTO(3), VA(3), VB(3), XINTDIS, AZIMUT(2)
REAL ZSEPPL, ZSEPMN, PARAL, AVCOTHE
CMI
C========================================================================
C
c Initialization of parameters (note that extrapolation sign, coded into
c ILIST, is not relevant for this routine). VWRBLK is the block dimension
c inside VWRK bank
c
IERR = 0
Offset = VWRBLK * (Nreg - 1) ! pointer to region of extrapolation
Jlist(1) = abs(Ilist(1))
Jlist(2) = abs(Ilist(2))
c
DO NDX=1,2
if ( ilist(ndx)*iswmback .gt. 0 ) then
VS(1,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWXFH)
VS(2,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWYFH)
VS(3,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWZFH)
TS(1,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWTXFH)
TS(2,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWTYFH)
TS(3,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWTZFH)
else
VS(1,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWXLH)
VS(2,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWYLH)
VS(3,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWZLH)
TS(1,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWTXLH)
TS(2,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWTYLH)
TS(3,NDX) = RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWTZLH)
endif
c
CALL MGFLD(VS(1,NDX),VH)
IF (ABS(VH(3)) .LE. 1.0E-4) THEN
IERR = 7
RETURN
ELSE
C
C KLOE convention: RHO is >0 for a positive particle, turning CLOCKWISE.
RHO(NDX) = 1./
+ (CONVERS*VH(3)*RW(Indvwrk+vwrnco*(Jlist(ndx)-1)+Offset+VWCUFH))
SGNQB(NDX) = SIGN(1.0,RHO(NDX))
C
C This is the angle from the center of curvature to the starting point
ANGSTR(NDX)=
+ ATAN2(SGNQB(NDX)*TS(1,NDX), -SGNQB(NDX)*TS(2,NDX))
C
C These are the coordinates of the center of curvature
TRMIV = 1./SQRT(TS(1,NDX)**2+TS(2,NDX)**2)
XCTR(1,NDX) = VS(1,NDX) + RHO(NDX)*TS(2,NDX)*TRMIV
XCTR(2,NDX) = VS(2,NDX) - RHO(NDX)*TS(1,NDX)*TRMIV
C
C This is the cotangent of theta (or the tangent of lambda, the dip angle)
CTGTHE(NDX) = TS(3,NDX)/SQRT(1.0-TS(3,NDX)*TS(3,NDX))
RHO(NDX) = ABS(RHO(NDX))
ENDIF
ENDDO
C ----------------->> ASSURE THAT RHO(IG) .GT. RHO(IL)
IF (RHO(1) .LT. RHO(2)) THEN
IL=1
IG=2
ELSE
IL=2
IG=1
ENDIF
C
C U is the vector joining the centers of the two circles in the tr. plane
U(1) = XCTR(1,IG)-XCTR(1,IL)
U(2) = XCTR(2,IG)-XCTR(2,IL)
DCC = SQRT(U(1)**2+U(2)**2)
C
C ----------------->> TEST FOR ALMOST TANGENT CIRCLES
C
EXTCRC = DCC .GT. (RHO(1)+RHO(2))
INTCRC = DCC .LT. (RHO(IG)-RHO(IL))
C
IF (EXTCRC .OR. INTCRC) THEN
C
C The two circles have *NO* intersections. Use the point of closest
C approach to both helices.
C Finally, if the circles are "tangent within cuts", flag TANGENT = .true.
CMI In this case set error flag
IERR = 1
CMI
ANGCTR(IL) = ATAN2(U(2),U(1))
ANGCTR(IG) = ATAN2(-U(2),-U(1))
IF (INTCRC) ANGCTR(IL) = ANGCTR(IG)
C
TURNPL(IL) = PROXIM(ANGCTR(IL)-ANGSTR(IL),0.)
ZPL(IL) = VS(3,IL) - SGNQB(IL)*TURNPL(IL)*RHO(IL)*CTGTHE(IL)
C
TURNPL(IG) = PROXIM(ANGCTR(IG)-ANGSTR(IG),0.)
ZPL(IG) = VS(3,IG) - SGNQB(IG)*TURNPL(IG)*RHO(IG)*CTGTHE(IG)
C
DLTZ = ABS(ZPL(IG)-ZPL(IL))
IF (INTCRC) THEN
A_HALF = (RHO(IG)-RHO(IL)-DCC)/2.
ELSE
A_HALF = (DCC-RHO(IG)-RHO(IL))/2.
ENDIF
XINT(1) = XCTR(1,IL) + (RHO(IL)+A_HALF)*COS(ANGCTR(IL))
XINT(2) = XCTR(2,IL) + (RHO(IL)+A_HALF)*SIN(ANGCTR(IL))
XINT(3) = 0.5*(ZPL(IG)+ZPL(IL))
ZSEP = abs(ZPL(IG)-ZPL(IL))
C
IF (INTCRC) THEN
TANGENT = RHO(IG)-DCC-RHO(IL) .LT. AMAX1(ABSTOLER,TOLER*DCC)
ELSE
TANGENT = DCC-RHO(1)-RHO(2) .LT.
+ AMAX1(ABSTOLER,TOLER*(RHO(1)+RHO(2)))
ENDIF
ELSE
C
C The two circles have exactly two intersections
TANGENT = .FALSE.
C
C This is the distance from the center of circle IL to the chord
C between the two intersections.
RLCOS = ((DCC**2-RHO(IG)**2)+RHO(IL)**2)/(2.*DCC)
ARGMNT = RLCOS/RHO(IL)
IF (ABS(ARGMNT) .GT. 1.) THEN
ARGMNT = SIGN(1.,ARGMNT)
ENDIF
C
C This is a half of the chord angle, defined positive
ANG(IL) = ACOS(ARGMNT)
C
C Same for circle IG....
ARGMNT = (DCC-RLCOS)/RHO(IG)
IF (ABS(ARGMNT) .GT. 1.) THEN
ARGMNT = SIGN(1.,ARGMNT)
ENDIF
ANG(IG) = ACOS(ARGMNT)
C
C This is the angle of center IG seen from IL, correctly in (-PI,PI)
ANGCTR(IL) = ATAN2(U(2),U(1))
SGNIL = SIGN(1.0,ANGCTR(IL))
C
C And here is the angle of center IL seen from IG, also in (-PI,PI)
ANGCTR(IG) = ATAN2(-U(2),-U(1))
SGNIG = SIGN(1.0,ANGCTR(IG))
C
C These 4 TURNs are the ang.distances from the 2 starting points up to the
C intersections (2 for each track). Signs are geometric: clockwise --> < 0
TURNPL(IL) = PROXIM(ANGCTR(IL)+SGNIL*ANG(IL)-ANGSTR(IL),0.)
C
C The 4 - signs below take into account the KLOE convention above
ZPL(IL) = VS(3,IL) - SGNQB(IL)*TURNPL(IL)*RHO(IL)*CTGTHE(IL)
C
TURNPL(IG) = PROXIM(ANGCTR(IG)+SGNIG*ANG(IG)-ANGSTR(IG),0.)
ZPL(IG) = VS(3,IG) - SGNQB(IG)*TURNPL(IG)*RHO(IG)*CTGTHE(IG)
C
TURNMN(IL) = PROXIM(ANGCTR(IL)-SGNIL*ANG(IL)-ANGSTR(IL),0.)
ZMN(IL) = VS(3,IL) - SGNQB(IL)*TURNMN(IL)*RHO(IL)*CTGTHE(IL)
C
TURNMN(IG) = PROXIM(ANGCTR(IG)-SGNIG*ANG(IG)-ANGSTR(IG),0.)
ZMN(IG) = VS(3,IG) - SGNQB(IG)*TURNMN(IG)*RHO(IG)*CTGTHE(IG)
C
C The choice of the good intersection is based on the Z-separation
c between the two helices. If both separations are small, check angular
c distance
c
ZSEPMN = abs(ZMN(IG)-ZMN(IL))
DLTCMN = TURNMN(IL)**2+TURNMN(IG)**2
ZSEPPL = abs(ZPL(IG)-ZPL(IL))
DLTCPL = TURNPL(IL)**2+TURNPL(IG)**2
if ( AMAX1(ZSEPPL,ZSEPMN) .lt. 2. ) then
if ( DLTCPL .lt. DLTCMN ) goto 112
else
if ( ZSEPPL .lt. ZSEPMN ) goto 112
endif
ZSEP = ZSEPMN
ANGNW = ANGCTR(IL) - SGNIL*ANG(IL)
XINT(3) = 0.5*(ZMN(IG)+ZMN(IL))
c Azimuthal angle at intersection point.
AZIMUT(IG) = ATAN2(TS(2,IG),TS(1,IG)) + TURNMN(IG)
AZIMUT(IL) = ATAN2(TS(2,IL),TS(1,IL)) + TURNMN(IL)
goto 113
112 continue
ZSEP = ZSEPPL
ANGNW = ANGCTR(IL) + SGNIL*ANG(IL)
XINT(3) = 0.5*(ZPL(IG)+ZPL(IL))
AZIMUT(IG) = ATAN2(TS(2,IG),TS(1,IG)) + TURNPL(IG)
AZIMUT(IL) = ATAN2(TS(2,IL),TS(1,IL)) + TURNPL(IL)
113 continue
cmi Here check for parallelism (used later in minimization procedure)
PARAL = abs( PI - abs(AZIMUT(IG)-AZIMUT(IL)) )
if ( PARAL .lt. PARALCUT ) TANGENT = .true.
cmi
XINT(1) = XCTR(1,IL) + RHO(IL)*COS(ANGNW)
XINT(2) = XCTR(2,IL) + RHO(IL)*SIN(ANGNW)
ENDIF
C
CMI Check wether z_separation is small enough. If not flag it out.
CMI
AVCOTHE = amax1(abs(CTGTHE(IL)),abs(CTGTHE(IG)))
If (ZSEP.gt.ZSEPMAX*(1+AVCOTHE)) then
IERR = 2
goto 999
Endif
CMI
cmi
cmi now an iterative procedure starts that tends to improve the
cmi initial conditions. The vectors VA and VB, joining XINT with the
cmi points closest to XINT, belonging to the two tracks, are built.
cmi The vector sum (VA+VB)/2 is evaluated and XINT is moved by this
cmi amount. The iteration ends either after a certain number of attempts
cmi (NITER) or when the new position is far from the old one less than
cmi CUTINTER cm.
cmi
cmi initialize counter
c LITER = 0
c
c do while (LITER .le. NITERCRSS)
c LITER = LITER + 1
c
cmi keep track of old interaction point
c XINTO(1) = XINT(1)
c XINTO(2) = XINT(2)
c XINTO(3) = XINT(3)
c
cmi get vector VA
c LTK = ILIST(1)
c call VXSWIM(XINTO, LTK, TRKBNK(1,4,LTK), LERRS)
c if (LERRS .ne. 0) goto 999 ! bad track swim
c VA(1) = XINTO(1) - TRKBNK(2,4,LTK)
c VA(2) = XINTO(2) - TRKBNK(3,4,LTK)
c VA(3) = XINTO(3) - TRKBNK(4,4,LTK)
c
cmi get vector VB
c LTK = ILIST(2)
c call VXSWIM(XINTO, LTK, TRKBNK(1,4,LTK), LERRS)
c if (LERRS .ne. 0) goto 999 ! bad track swim
c VB(1) = XINTO(1) - TRKBNK(2,4,LTK)
c VB(2) = XINTO(2) - TRKBNK(3,4,LTK)
c VB(3) = XINTO(3) - TRKBNK(4,4,LTK)
c
ccmi get new interaction point
c XINT(1) = XINTO(1) + (VA(1)+VB(1))/2.
c XINT(2) = XINTO(2) + (VA(2)+VB(2))/2.
c XINT(3) = XINTO(3) + (VA(3)+VB(3))/2.
c
ccmi check distance with old interaction point
c XINTDIS = ( XINT(1) - XINTO(1) )**2 +
c + ( XINT(2) - XINTO(2) )**2 +
c + ( XINT(3) - XINTO(3) )**2
c XINTDIS = sqrt(XINTDIS)
c
c if (XINTDIS .lt. CUTINTER) goto 999
c
c enddo
ccmi
c
999 continue
RETURN
END
c=============================================================
integer function VTXSWM(iii,Ireg,XR)
c=============================================================
c
c This routine finds the coordinates of the point XP, belonging to
c track iii, CLOSEST IN SPACE to a given point XR (for example a
c candidate vertex).
c The point is found via an iterative procedure.
c - If Ireg=1-4 the track is extrapolated moving backward starting
c from the corresponding "First Hit" (1=FH,2=DC,3=BI,4=BU - see
c VTXFIN.INC). The difference between using one of these "First
c hits" is the energy loss due to DC and BP boundaries.
c - If Ireg=6 the track is extrapolated from the last point, in the DC
c chain of hits, moving outward.
c - If Ireg=0 then the routine starts from the FH and it tries to find
c the direction of extrapolation by rotating using Phi and Phi+-2Pi,
c chosing the "best one" [FOR THE MOMENT THIS IS NOT WORKING].
c
c Author: Marco Incagli
c =======
c
c Creation date: 25-9-1997
c ==============
c
c Input:
c ======
c iii Track number
c Ireg Region of extrapolation (see VTXFIN.INC)
c XR(3) Coordinates of external point
c
c Output:
c =======
c the following variables are filled inside bank VWRK:
c
c RW(INDVWRK+...+VWXVi), VWYVi, VWZVi, VWTXVi, VWTYVi, VWTZVi,
c VWCUVi, VWTLVi
c
c (position, tangent, Q/Pt and Track Length from vertex to FH or LH)
c
c where i=L , if Ireg=NLH (Vertex at Last point) , or i=X otherwise.
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c Global declarations
c
INTEGER NITER ! number of iterations IN VTXSWM
PARAMETER (NITER=10)
REAL EPSILON ! End of iteration procedure
PARAMETER (EPSILON=0.1)
REAL ALPHATOL
PARAMETER (ALPHATOL=0.3) ! Tolerance on "wrong extrapolation"
c
c Local declarations
c
INTEGER Ireg, iii, I1, i2, IWDATA
REAL xfh(3), vh(3)
REAL XP(3), XR(3), X0, Y0, Z0, TX0, TY0, TZ0, Bfield, Alphatot
REAL rho, Q, costhe, sinthe, cotthe, sinphi, cosphi, Xc, Yc
REAL Rnorm, sinalpha, cosalpha, alpha, PRdist, Rcur, prdistold
REAL PX(NITER),PY(NITER),PZ(NITER),TPX(NITER),TPY(NITER),TPZ(NITER)
CHARACTER*80 message
c
c================================================================
c --- Get initial conditions according to IREG parameter
c
XP(1) = 0.
XP(2) = 0.
XP(3) = 0.
VTXSWM = -1 ! set initial status to failure
IWDATA = INDVWRK+VWRNCO*(iii-1) ! data pointer for track iii in VWRK
c check IREG with VWREG parameter in VWRK bank
if ( IREG.ne.NLH .and. IREG.gt.IW(IWDATA+VWREG) ) then
message = 'Track cannot be extrapolated in this region'
call ERLOGR ('VTXSWM',ERWARN,0,0,message)
return
endif
If (IREG.eq.NFH .OR. IREG.EQ.0) then
X0 = RW(IWdata+VWXFH)
Y0 = RW(IWdata+VWYFH)
Z0 = RW(IWdata+VWZFH)
TX0 = RW(IWdata+VWTXFH)
TY0 = RW(IWdata+VWTYFH)
TZ0 = RW(IWdata+VWTZFH)
c get B field at First Hit
xfh(1) = X0
xfh(2) = Y0
xfh(3) = Z0
call MGFLD(xfh,Vh)
Bfield = Vh(3)
rho = CONVERS*Bfield*RW(IWdata+VWCUFH)
if ( abs(rho).lt.0.0001 ) then
message = 'Error in swimming track'
call ERLOGR ('VTXSWM',ERWARN,0,0,message)
return
endif
Q = SIGN(1.,rho)
rho = 1./rho
elseif (IREG.eq.NVC) then
X0 = RW(IWdata+VWXDC)
Y0 = RW(IWdata+VWYDC)
Z0 = RW(IWdata+VWZDC)
TX0 = RW(IWdata+VWTXDC)
TY0 = RW(IWdata+VWTYDC)
TZ0 = RW(IWdata+VWTZDC)
c get B field at First Hit
xfh(1) = X0
xfh(2) = Y0
xfh(3) = Z0
call MGFLD(xfh,Vh)
Bfield = Vh(3)
rho = CONVERS*Bfield*RW(IWdata+VWCUDC)
if ( abs(rho).lt.0.0001 ) then
message = 'Error in swimming track'
call ERLOGR ('VTXSWM',ERWARN,0,0,message)
return
endif
Q = SIGN(1.,rho)
rho = 1./rho
elseif (IREG.eq.NBI) then
X0 = RW(IWdata+VWXBI)
Y0 = RW(IWdata+VWYBI)
Z0 = RW(IWdata+VWZBI)
TX0 = RW(IWdata+VWTXBI)
TY0 = RW(IWdata+VWTYBI)
TZ0 = RW(IWdata+VWTZBI)
c get B field at First Hit
xfh(1) = X0
xfh(2) = Y0
xfh(3) = Z0
call MGFLD(xfh,Vh)
Bfield = Vh(3)
rho = CONVERS*Bfield*RW(IWdata+VWCUBI)
if ( abs(rho).lt.0.0001 ) then
message = 'Error in swimming track'
call ERLOGR ('VTXSWM',ERWARN,0,0,message)
return
endif
Q = SIGN(1.,rho)
rho = 1./rho
elseif (IREG.eq.NBU) then
X0 = RW(IWdata+VWXBU)
Y0 = RW(IWdata+VWYBU)
Z0 = RW(IWdata+VWZBU)
TX0 = RW(IWdata+VWTXBU)
TY0 = RW(IWdata+VWTYBU)
TZ0 = RW(IWdata+VWTZBU)
c get B field at First Hit
xfh(1) = X0
xfh(2) = Y0
xfh(3) = Z0
call MGFLD(xfh,Vh)
Bfield = Vh(3)
rho = CONVERS*Bfield*RW(IWdata+VWCUBU)
if ( abs(rho).lt.0.0001 ) then
message = 'Error in swimming track'
call ERLOGR ('VTXSWM',ERWARN,0,0,message)
return
endif
Q = SIGN(1.,rho)
rho = 1./rho
elseif (IREG.eq.NLH) then
X0 = RW(IWdata+VWXLH)
Y0 = RW(IWdata+VWYLH)
Z0 = RW(IWdata+VWZLH)
TX0 = RW(IWdata+VWTXLH)
TY0 = RW(IWdata+VWTYLH)
TZ0 = RW(IWdata+VWTZLH)
c get B field at First Hit
xfh(1) = X0
xfh(2) = Y0
xfh(3) = Z0
call MGFLD(xfh,Vh)
Bfield = Vh(3)
rho = CONVERS*Bfield*RW(IWdata+VWCULH)
if ( abs(rho).lt.0.0001 ) then
message = 'Error in swimming track'
call ERLOGR ('VTXSWM',ERWARN,0,0,message)
return
endif
Q = SIGN(1.,rho)
rho = 1./rho
else
message = 'Unknown detector region'
call ERLOGR ('VTXSWM',ERWARN,0,0,message)
return
endif
c --- Get helix parameters (center, radius and cot(theta))
costhe = TZ0
sinthe = sqrt(1-costhe**2)
cotthe = costhe/sinthe
cosphi = TX0/sinthe
sinphi = TY0/sinthe
Xc = X0 + rho*sinphi
Yc = Y0 - rho*cosphi
c --- Get starting point of iterative procedure
c Starting point P1 of iterative procedure is on the line C-R. Find
c rotation angle "alpha" from 0 to P1.
c While doing a backward extrapolation, the sign of Alpha is positive
c (counterclockwise) for positive tracks and negative for negative
c tracks. The situation is the opposite for forward extrapolations (from
c the last point moving outward).
c HOWEVER this approach must have a CONTROLLED TOLERANCE. In fact, the
c candidate vertex may be put in the "wrong position with respect to X0";
c in this situation "alpha" has the wrong sign.
c 1 - get SIN of rotation angle using vector product:
rnorm = ( (X0-Xc)**2 + (Y0-Yc)**2 ) *
+ ( (Xr(1)-Xc)**2 + (Xr(2)-Yc)**2 )
rnorm = sqrt(rnorm)
sinalpha = (X0-Xc)*(Xr(2)-Yc) - (Y0-Yc)*(Xr(1)-Xc)
sinalpha = sinalpha / rnorm
If (abs(sinalpha) .ge. 1.0) sinalpha = SIGN(0.9999,sinalpha)
c 2 - check sign according to direction of extrapolation and track charge
c (the passage from sinalpha to cosalpha is made in order to correctly
c define the angle alpha in the interval [-pi;pi])
If ( Ireg.eq.NLH ) then ! sign of alpha opposite to charge
if ( sinalpha*Q .lt. 0 ) then ! normal situation
cosalpha = sqrt(1. - sinalpha**2)
alpha = -Q * ACOS(cosalpha) ! defined in [-pi:pi]
else ! wrong situation; check tolerance
cosalpha = sqrt(1. - sinalpha**2)
alpha = Q * ACOS(cosalpha)
if (abs(alpha).gt.alphatol) return
endif
else ! alpha same sign as charge
if ( sinalpha*Q .gt. 0 ) then ! normal situation
cosalpha = sqrt(1 - sinalpha**2)
alpha = Q * ACOS(cosalpha)
else ! wrong situation; check tolerance
cosalpha = sqrt(1 - sinalpha**2)
alpha = -Q * ACOS(cosalpha)
if (abs(alpha).gt.alphatol) return
endif
endif
c 3 - get tangent and coordinates at P1, start of iteration, through a
c rotation of "alpha"
TPX(1) = ( cosphi*cos(Alpha) - sinphi*sin(alpha) ) * sinthe
TPY(1) = ( sinphi*cos(Alpha) + cosphi*sin(alpha) ) * sinthe
TPZ(1) = costhe
PX(1) = Xc - rho*TPY(1)/sinthe
PY(1) = Yc + rho*TPX(1)/sinthe
PZ(1) = Z0 - rho*alpha*cotthe
c keep track of rotation angle
alphatot = alpha
c----------------------------------------------------------------------------
c
c Now start iterative procedure
c
PRdistold = 9999.
Do i1 = 2, NITER
i2 = i1
c get distance P(i-1) - R', R' being the projection of R on the 3-d tangent
c to the helix; this distance is signed.
PRdist = (PX(i1-1)-Xr(1))*TPX(i1-1)+(PY(i1-1)-Xr(2))*TPY(i1-1)
+ + (PZ(i1-1)-Xr(3))*TPZ(i1-1)
c check if distance increased in this iteration; if yes exit
if ( abs(PRdist) .gt. abs(PRdistold) ) then
i2 = i1-1
goto 999
endif
c check if distance is below cut
if (abs(PRdist) .lt. epsilon ) goto 999
c If not find new point and iterate. New point is found moving by an
c arc of length PRdist along the curve. Note that PRdist is a signed
c quantity (and RHO, too).
alpha = PRdist / rho
TPX(i1) = TPX(i1-1)*cos(Alpha) - TPY(i1-1)*sin(alpha)
TPY(i1) = TPY(i1-1)*cos(Alpha) + TPX(i1-1)*sin(alpha)
TPZ(i1) = TPZ(i1-1)
PX(i1) = Xc - rho*TPY(i1)/sinthe
PY(i1) = Yc + rho*TPX(i1)/sinthe
PZ(i1) = PZ(i1-1) - rho*alpha*cotthe
PRdistold = PRdist
c keep track of rotation angle
alphatot = alphatot + alpha
enddo
c------------------------------------------------------------------------
999 continue
c output results
Rcur = 1./(convers*Bfield*Rho)
If (Ireg .ne. NLH) then
RW(IWdata+VWXVX) = PX(I2-1)
RW(IWdata+VWYVX) = PY(I2-1)
RW(IWdata+VWZVX) = PZ(I2-1)
RW(IWdata+VWTXVX) = TPX(I2-1)
RW(IWdata+VWTYVX) = TPY(I2-1)
RW(IWdata+VWTZVX) = TPZ(I2-1)
RW(IWdata+VWCUVX) = Rcur
RW(IWdata+VWTLVX) = abs(alphatot*rho)
else
RW(IWdata+VWXVL) = PX(I2-1)
RW(IWdata+VWYVL) = PY(I2-1)
RW(IWdata+VWZVL) = PZ(I2-1)
RW(IWdata+VWTXVL) = TPX(I2-1)
RW(IWdata+VWTYVL) = TPY(I2-1)
RW(IWdata+VWTZVL) = TPZ(I2-1)
RW(IWdata+VWCUVL) = Rcur
RW(IWdata+VWTLVL) = abs(alphatot*rho)
ENDIF
VTXSWM = YESUCC
return
end
c=======================================================================
SUBROUTINE VTXMIN(Nflag,Ntrfit,ILIST,XVX,CSQ,IFIT,XCN,CONSTR)
c=======================================================================
c
c Routine perfoming the minimization of the distance tracks - candidate
c vertex.
c
c Author: Marco Incagli
c =======
c
c Creation date: 26-9-1997
c ==============
c
c Input:
c ======
c Nflag Section of program in which cand. vertex was found
c Ntrfit Number of tracks to be fitted
c Ilist(Ntrfit) Track numbers
c XVX(3) Coordinates of candidate vertex
c Constr Logical value; if true then do a constrained fit
c XCN(3) Coordinates of the constraint
c
c Output:
c =======
c XVX(3) Coordinates of fitted vertex
c CSQ CHI2 value of this fit
c Ifit flag (values in VTXFIN)
c = FITGOLD -> "successfull" fit
c = FITSILVER -> "acceptable" fit
c = FITBAD -> fit converged but CHI2 too high
c = FITNOCONV -> fit didn't converge
c = FITNEGCHI2 -> negative CHI2
c = FITMAXITER -> too many iterations
c = FITTRKEXTR -> error in track extrapolation (VTXSWM)
c = FITOUTDC -> cand vertex outside DC boundaries
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c External declarations
c
INTEGER Istat, VTXSWM, IWDATA
c
c Local declarations
c
INTEGER i1, i2, iii, Niter, ilist(10), ibadtk, Ntrfit, IWfitflg(5)
INTEGER Ireg, Iireg, Idum(3), nerror, Tolerate, Ifit, Nflag
REAL Xcn(3), miverr(3), midv(3), vershift, vershprev
REAL mixgb(3), vh(3), rho, theta, Parallel
REAL rpca(3), trpca(3), cotheta
REAL midx, midy, midt, midz, Xvx(3), txynorm, cosa, sina
REAL sig2t, sig2z, sigtz, denom, vart, varz, covar, csq
REAL Trkcsq(5), den, mitx, mity
REAL mixgb1, mixgb2, mixgb3, mixga11, mixga12, mixga13
REAL mixga22, mixga23, mixga33, xga(3,3)
CHARACTER*80 message
LOGICAL Constr, converg, accept, verterr
REAL Rdc, Rbp, rvertex, con1, con2, con3
EQUIVALENCE (TIRGVX(1),Rdc),(TIRGVX(2),Rbp)
c
c common passed to VTXCAN from VTXMIN for creating candidate bank DVCS
c
INTEGER Jfit, Nflg, Jlist(2), Flagfl(2)
REAL Vx, Vy, Vz, mixga(3,3), Vchi2, VTrkcsq(2)
COMMON /DVCS_COM/ Jfit, Nflg, Jlist, Flagfl,
+ Vx, Vy, Vz, mixga, Vchi2, VTrkcsq
c
c
c=======================================================================
c
c initialize variables
c
if ( Ntrfit .gt. 10 ) then
message = 'More than 10 tracks connected to a vertex'
call ERLOGR ('VTXMIN',ERWARN,0,0,message)
return
endif
c
verterr = .false. ! set to true if vertex outside detector boundaries
NITER = 0 ! iteration nummber
vershprev = 999999. ! vertex shift (decrease at each iteration)
Tolerate = 0 ! we tolerate one "increase" in vertex shift
Parallel = 1. ! if "parallesim" flag is up (CONSTR=.true.)
If (constr) Parallel=2. ! then loose convergence criteria.
c
c-----------------------------------------------------------------------------
c Find region of candidate vertex.
c Note that this region needs not to be the same for all tracks.
c For example, if the vertex is in the DC region, than one
c track may have the First Hit connected to it (Ireg=NFH),
c while the other may have the Last Hit (Ireg=NLH).
c
rvertex = sqrt( xvx(1)**2 + xvx(2)**2 )
if ( rvertex .lt. Rbp ) then
Ireg = NBI
elseif ( rvertex .lt. Rdc ) then
Ireg = NVC
elseif ( rvertex .lt. Rcal*1.05 ) then ! 50f tolerance
Ireg = NFH ! Assume First Hit (eventually change later)
else
message = 'Candidate vertex outside calorimeter boundaries'
call ERLOGR ('VTXMIN',ERWARN,0,0,message)
verterr = .true. ! raise a flag that will be checked later
endif
c Find pointers used later to store fit results, but first check if
c this is a BACKWORD or FORWARD extrapolation; in this last case the
c vertex is associated to the track Last Hit => change IWFITFLG and FLAGFL.
do 414 i1 = 1, Ntrfit
iii = abs( ilist(i1) )
IWdata = INDVWRK+VWRNCO*(iii-1) ! data pointer for track iii in VWRK
IWfitflg(i1) = IWdata + VWFLF ! pointer to fit flag (assume 1st point)
Flagfl(i1) = Flagfh ! flag First Hit
If ( Ilist(i1)*ISWMFORW .gt. 0. .or. verterr ) then
IWfitflg(i1) = IWdata + VWFLL ! change to last point
Flagfl(i1) = Flaglh ! flag First Hit
endif
414 continue
c check VERTERR flag (IWFITFLG pointers are needed, so this check is
c performed after loop 414)
If ( verterr ) goto 305
c
c Start of iterative process
c
888 continue
NITER = NITER + 1
c check number of iterations
if ( NITER .gt. Nitfitmax ) goto 303
c
c ===========================================================
c loop on tracks and build derivative matrix and error vector
c ===========================================================
c
c
c zero CHI2, derivative matrix and error vector
c
CSQ = 0.
do i1=1,3
mixgb(i1) = 0.
do i2=1,3
mixga(i1,i2) = 0.
enddo
enddo
c get point of closest approach to candidate vertex (VPCA)
do 415 i1 = 1, Ntrfit
Iireg = Ireg ! use Iireg that can be different from Ireg
If ( Ilist(i1)*ISWMFORW .gt. 0 ) Iireg = NLH
iii = abs( ilist(i1) )
IWdata = INDVWRK+VWRNCO*(iii-1) ! data pointer for track iii in VWRK
istat = VTXSWM(iii,Iireg,XVX)
if ( ISTAT .ne. YESUCC ) then
c
c If swimming routine returns error check wether the candidate vertex is
c close to the boarder between two regions. If yes try again swimming
c with a different Iireg value
c
if ( iireg.eq.nbi .and. rvertex.lt.0.8*rbp ) then
goto 304 ! not at boarder
elseif ( iireg.eq.nvc .and. rvertex.lt..8*rdc) then
goto 304 ! not at boarder
elseif ( iireg.eq.nfh .or. iireg.eq.nlh ) then
goto304
else
iireg = iireg-1 ! try swimming with external region
istat = VTXSWM(iii,Iireg,XVX)
if ( ISTAT .ne. YESUCC ) goto 304 !No hope; flag error and ret
endif
endif
c
If ( Ilist(i1)*ISWMBACK .gt. 0. ) then
rpca(1) = RW(IWdata+VWXVX)
rpca(2) = RW(IWdata+VWYVX)
rpca(3) = RW(IWdata+VWZVX)
trpca(1) = RW(IWdata+VWTXVX)
trpca(2) = RW(IWdata+VWTYVX)
trpca(3) = RW(IWdata+VWTZVX)
call mgfld(RW(IWdata+VWXVX),vh)
rho = 1./(convers * vh(3) *
+ abs(RW(IWdata+VWCUVX)) )
else
rpca(1) = RW(IWdata+VWXVL)
rpca(2) = RW(IWdata+VWYVL)
rpca(3) = RW(IWdata+VWZVL)
trpca(1) = RW(IWdata+VWTXVL)
trpca(2) = RW(IWdata+VWTYVL)
trpca(3) = RW(IWdata+VWTZVL)
call mgfld(RW(IWdata+VWXVL),vh)
rho = 1./(convers * vh(3) *
+ abs(RW(IWdata+VWCUVL)) )
endif
c if at first iteration, than evaluate error matrix at VPCA
IF ( NITER .EQ. 1 ) call vtxerr(iii, Iireg, rpca, trpca)
sig2t = RW(IWdata+VWERRV+0)
sigtz = RW(IWdata+VWERRV+1)
sig2z = RW(IWdata+VWERRV+5)
c this is the vector from PCA to VX
midx = XVX(1) - rpca(1)
midy = XVX(2) - rpca(2)
midz = XVX(3) - rpca(3)
c this is the distance PCA-VX on transverse plane
txynorm = 1/sqrt(trpca(1)**2+trpca(2)**2)
midt = txynorm * (midx*trpca(2) - midy*trpca(1))
c alpha is the "rotation angle" along the helix at PCA; cosa=cos(alpha); sina...
mitx = trpca(1)
mity = trpca(2)
cosa = mity * txynorm ! cos(alpha) = sin(phi)
sina = - mitx * txynorm ! sin(alpha) = -cos(phi)
cotheta = trpca(3) / sqrt( trpca(1)**2 + trpca(2)**2 )
c correct cot(theta) by vertex distance from helix (normally theta ~
c cot(theta), since midt<<rho)
theta = cotheta * rho/(rho-midt)
c before evaluating the matrix, calculate the total and the diagonal CHI2
DENOM = 1 - sigtz**2/(sig2t*sig2z)
VART = midt**2/sig2t
VARZ = midz**2/sig2z
COVAR = - 2 * abs(midt) * abs(midz) * sigtz / sig2t / sig2z
c CSQ = global diagonal CHI2 (the sqrt will be taken after the loop)
c Trkcsq = CHI2 of single track (including correlations)
CSQ = CSQ + VART + VARZ
Trkcsq(i1) = ( VART + COVAR + VARZ ) / DENOM
IF (Trkcsq(i1) .LT. 0.0) Trkcsq(i1) = 100000000.
Trkcsq(i1) = sqrt(Trkcsq(i1))
c now evaluate matrix
DEN = sig2t*sig2z - sigtz**2
c
mixgb1 = - cosa*sig2z*midt + theta*sina*sig2t*midz +
& ( cosa*midz - theta*sina*midt ) * sigtz
mixgb2 = - sina*sig2z*midt + theta*cosa*sig2t*midz +
& ( sina*midz - theta*cosa*midt ) * sigtz
mixgb3 = - sig2t*midz + sigtz*midt
mixgb(1) = mixgb(1) + mixgb1/DEN
mixgb(2) = mixgb(2) + mixgb2/DEN
mixgb(3) = mixgb(3) + mixgb3/DEN
*
mixga11 = cosa**2*sig2z + theta**2*sina**2*sig2t +
& 2*theta*sina*cosa*sigtz
mixga12 = sina*cosa*(theta**2*sig2t+sig2z) +
& theta*sigtz
mixga13 = - theta*sina*sig2t - cosa*sigtz
mixga22 = sina**2*sig2z + theta**2*cosa**2*sig2t +
& 2*theta*sina*cosa*sigtz
mixga23 = - theta*cosa*sig2t - sina*sigtz
mixga33 = sig2t
mixga(1,1) = mixga(1,1) + mixga11/DEN
mixga(1,2) = mixga(1,2) + mixga12/DEN
mixga(1,3) = mixga(1,3) + mixga13/DEN
mixga(2,2) = mixga(2,2) + mixga22/DEN
mixga(2,3) = mixga(2,3) + mixga23/DEN
mixga(3,3) = mixga(3,3) + mixga33/DEN
415 enddo
c Fix global CHI2 (sum in quadrature of all diagonal CHI2)
CSQ = sqrt(CSQ)
c Check CHI2 cut on single track; flag worst track (if above limit)
do i1= 1, ntrfit
Ibadtk = 0 ! no BAD tracks
if ( Trkcsq(i1) .gt. Cutrkcsq ) Ibadtk = iii
enddo
c now apply the constraint (if asked) => fake track with coordinates
c (con1,con2,con3), wher con1,2,3 are defined below
if ( .not.CONSTR ) goto 416
con1 = XVX(1) - XCN(1)
con2 = XVX(2) - XCN(2)
con3 = XVX(3) - XCN(3)
midt = sqrt(con1**2+con2**2)
midz = abs(con3)
if ( midt.lt.0.00001 ) then
cosa = 1.
sina = 0.
theta = 0.
else
cosa = COS( ATAN2(con2,con1) )
sina = SIN( ATAN2(con2,con1) )
theta = con3 / sqrt(midt**2+midz**2)
if ( (abs(theta)-1).lt..0001 ) theta = SIGN(0.999,theta)
theta = ACOS( theta )
theta = 1. / TAN(theta)
endif
sig2t = sig2tcon
sig2z = sig2zcon
sigtz = 0.
c now update matrix with constraint conditions
DEN = sig2t*sig2z - sigtz**2
c
mixgb1 = - cosa*sig2z*midt + theta*sina*sig2t*midz +
& ( cosa*midz - theta*sina*midt ) * sigtz
mixgb2 = - sina*sig2z*midt + theta*cosa*sig2t*midz +
& ( sina*midz - theta*cosa*midt ) * sigtz
mixgb3 = - sig2t*midz + sigtz*midt
mixgb(1) = mixgb(1) + mixgb1/DEN
mixgb(2) = mixgb(2) + mixgb2/DEN
mixgb(3) = mixgb(3) + mixgb3/DEN
*
mixga11 = cosa**2*sig2z + theta**2*sina**2*sig2t +
& 2*theta*sina*cosa*sigtz
mixga12 = sina*cosa*(theta**2*sig2t+sig2z) +
& theta*sigtz
mixga13 = - theta*sina*sig2t - cosa*sigtz
mixga22 = sina**2*sig2z + theta**2*cosa**2*sig2t +
& 2*theta*sina*cosa*sigtz
mixga23 = - theta*cosa*sig2t - sina*sigtz
mixga33 = sig2t
c
c before modifying the matrix, store it. The final covariance matrix
c is the inverse of this matrix BEFORE being modified by the constraint.
c
do i1=1,3
do i2=i1,3
xga(i1,i2) = mixga(i1,i2)
enddo
enddo
mixga(1,1) = mixga(1,1) + mixga11/DEN
mixga(1,2) = mixga(1,2) + mixga12/DEN
mixga(1,3) = mixga(1,3) + mixga13/DEN
mixga(2,2) = mixga(2,2) + mixga22/DEN
mixga(2,3) = mixga(2,3) + mixga23/DEN
mixga(3,3) = mixga(3,3) + mixga33/DEN
416 continue
miXGA(2,1)=miXGA(1,2)
miXGA(3,1)=miXGA(1,3)
miXGA(3,2)=miXGA(2,3)
c
c----------------------------------------------------------------------------
c
c Use derivatives to evaluate new vertex position
c
c 1 - keep track of error vector
c
do i1 = 1,3
miverr(i1)=mixgb(i1)
enddo
c
c 2 - solve linearized system
c
call reqinv(3,mixga,3,idum,nerror,1,mixgb)
c
c 3 - get displacement
c
do i1=1,3
midv(i1)=mixgb(i1)
enddo
c
c 4 - get "vertex shift" (= how large was this step, in the iterative
c process, normalized with the inverse error vector)
vershift = 0.
do i1=1,3
vershift = vershift + miverr(i1)*midv(i1)
enddo
c
c 5 - flag error if normalized vertex dshift is negative
if ( vershift .lt. 0. ) goto 302
c
c 6 - Allows for an increment in vertex shift from previous step to this one,
c (TOLERATE), but only once.
c
if ( vershift .gt. vershprev ) then
Tolerate = Tolerate + 1
if ( Tolerate .gt. 1 ) goto 301
endif
c
vershprev = vershift
c
c 7 - update vertex position
c
XVX(1) = XVX(1) + midv(1)
XVX(2) = XVX(2) + midv(2)
XVX(3) = XVX(3) + midv(3)
c
c 8 - Check for fit convergence (normalized vector shift smaller than a
c given cut). The "Parallel" factor looses the fit convergence
c requirements for quasi-parallel tracks.
c
if ( vershift .lt. minshift*parallel ) then
converg = .true.
c
c If constraint flag was raised, restore the old XGA matrix, invert it
c and pass it to the VTXCAN routine IF POSITIVELY DEFINED.
c
if (constr) then
xga(2,1) = xga(1,2)
xga(3,1) = xga(1,3)
xga(3,2) = xga(2,3)
call rinv(3,xga,3,idum,nerror)
if (xga(1,1).gt.0.and.xga(2,2).gt.0
& .and.xga(3,3).gt.0) then
do i1=1,3
do i2=1,3
mixga(i1,i2) = xga(i1,i2)
enddo
enddo
endif
endif
c
goto 300
endif
c
c 9 - Iterate
c
goto 888
c==================================================================
c
c end of iterative procedure; fill FIT code (see comment above for
c meaning of "parallel" parameter).
c
300 if ( Ibadtk .eq. 0 ) then
IFIT = FITGOLD
elseif ( CSQ .lt. CHI2limit*parallel ) then
IFIT = FITSILVER
else
IFIT = FITBAD
endif
goto 999
c
301 IFIT = FITNOCONV
goto 999
c
302 IFIT = FITNEGCHI2
goto 999
c
303 IFIT = FITMAXITER
goto 999
c
304 IFIT = FITTRKEXTR
goto 999
c
305 IFIT = FITOUTDC
goto 999
c
999 continue
c fill track codes
do i1 = 1, ntrfit
IW( IWfitflg(i1) ) = IFIT
enddo
c If fit was succesfull than prepare the common block for VTXCAN and
c call the routine to store this candidate vertex.
c Note that DVCS bank is created only for 2 tracks vertices. Vertices
c with more tracks are built in VTXMERG and written directly into
c DVFS bank.
Accept = Ntrfit.eq.2 .and.
+ ( IFIT.eq.FITGOLD .or. IFIT.eq.FITSILVER )
If ( Accept ) then
c fill common block with infos to be put into the candidate bank
Nflg = Nflag ! internal flag
Jfit = Ifit ! quality flag
cmi
cmi just for debugging
cmi
if (constr) jfit = -Jfit
cmi
Vchi2 = CSQ ! (sqrt of) CHI squared
Nvxfit = Nvxfit + 1 ! vertex index
Vx = XVX(1) !
Vy = XVX(2) ! Vertex position
Vz = XVX(3) !
do i1 = 1, 2
Jlist(i1) = abs(Ilist(i1)) ! track number
Vtrkcsq(i1) = Trkcsq(i1) ! trk contribution to CHI2
enddo
c call the routine that creates DVCS bank (Candidate bank)
call VTXCAN
endif
c
return
end
c=======================================================================
SUBROUTINE VTXERR(iii, Ireg, rpca, trpca)
c=======================================================================
c
c Routine perfoming the extrapolation to point RPCA(3), in region Ireg,
c of error matrix.
c
c Author: Marco Incagli
c =======
c
c Creation date: 30-9-1997
c ==============
c
c Input:
c ======
c iii track number
c Ireg region of extrapolation
c rpca(3) coordinates of extrapolation point
c trpca(3) tangent at extrapolation point
c
c Output:
c =======
c RW(...+VWERRV) vector (15 elements) filled up
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$Include 'K$inc:Bcs.inc'
$$include 'ybos$library:errcod.inc'
$$include 'k$inc:jobsta.inc'
$$include 'k$inc:erlevl.inc' ! error logger
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:VWRK.INC'
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c Local declarations
c
INTEGER iii, Ireg, Iregcon, IWDATA, IWDATAN
REAL rpca(3), trpca(3), sinthe, Vh(3), tBP, tDC, sDC, sBP
REAL cosOMEGA, tEFF, sEFF, SIG2phi, SIG2cot, SIG2z
REAL Ptot, Etot, Bfield, SIG2kur
CHARACTER*80 message
c Parameters in common block (VTXERR - VTXTRANS)
REAL rxi, ryi, rzi, txi, tyi, exi, eyi,
+ rxf, ryf, rzf, txf, tyf, exf, eyf,
+ Rho, Kur,
+ V11i, V12i, V13i, V14i, V15i,
+ V22i, V23i, V24i, V25i, V33i,
+ V34i, V35i, V44i, V45i, V55i,
+ V11f, V12f, V13f, V14f, V15f,
+ V22f, V23f, V24f, V25f, V33f,
+ V34f, V35f, V44f, V45f, V55f
COMMON / TRANSCOM / rxi, ryi, rzi, txi, tyi, exi, eyi,
+ rxf, ryf, rzf, txf, tyf, exf, eyf,
+ Rho, Kur,
+ V11i, V12i, V13i, V14i, V15i,
+ V22i, V23i, V24i, V25i, V33i,
+ V34i, V35i, V44i, V45i, V55i,
+ V11f, V12f, V13f, V14f, V15f,
+ V22f, V23f, V24f, V25f, V33f,
+ V34f, V35f, V44f, V45f, V55f
c
c Thickness, in gr/cm**2, of the materials that make up the walls
c
REAL XAldc, XBebp
EQUIVALENCE (TIRGVX( 6),XAldc),(TIRGVX(11),XBebp)
c
c Radiation length, in gr/cm**2, of the materials that make up the walls
c
REAL RlAl, RlBe
EQUIVALENCE (TIDEVX(3,2), RlAl), (TIDEVX(3,7), RlBe)
c
c Density, in gr/cm**3, of the materials
c
Real RhoCf, RhoBe
PARAMETER (RhoCf=1.57,RhoBe=1.848)
c
c=======================================================================
c
c Initialize parameters
c
tDC = XAldc/RlAl ! thickness of DC wall in Rad. Len. units
sDC = XAldc/RhoCf ! thickness of DC wall in cm
tBP = XBebp/RlBe ! thickness of BP wall in Rad. Len. units
sBP = XBebp/RhoBe ! thickness of BP wall in cm
IREGCON = 1 ! control word
IWDATA = INDVWRK+VWRNCO*(iii-1) ! data pointer for track iii in VWRK
c
c Covariant matrix at First point
c
V11i = RW( IWDATA + VWERRF + 0 )
V12i = RW( IWDATA + VWERRF + 1 )
V13i = RW( IWDATA + VWERRF + 2 )
V14i = RW( IWDATA + VWERRF + 3 )
V15i = RW( IWDATA + VWERRF + 4 )
V22i = RW( IWDATA + VWERRF + 5 )
V23i = RW( IWDATA + VWERRF + 6 )
V24i = RW( IWDATA + VWERRF + 7 )
V25i = RW( IWDATA + VWERRF + 8 )
V33i = RW( IWDATA + VWERRF + 9 )
V34i = RW( IWDATA + VWERRF + 10 )
V35i = RW( IWDATA + VWERRF + 11 )
V44i = RW( IWDATA + VWERRF + 12 )
V45i = RW( IWDATA + VWERRF + 13 )
V55i = RW( IWDATA + VWERRF + 14 )
c get B field at First Hit
call MGFLD(RW(IWDATA+VWXFH),Vh)
Bfield = Vh(3)
c === Now start extrapolation of covariant matrix ===
c Two comments on conventions used
c 1 - In all of the following :
c Txi, Tyi = tangent unit vector on trasverse plane at initial point
c Exi, Eyi = versor defined as E = T X z (z being the B-field direction)
c Rxi, Ryi, Rzi = Position of initial point of extrapolation
c Txf, Tyf, Exf, Eyf, Rxf, ... = as before for final point
c Rho, Kur = rafius of curvature and curvature (defined as Q/Pt)
c
c 2 - To get from VWRK bank the informations relative to First Point, DC,
c or BP, I use the fact that each subregion is described by VWRBLK entries
c in VWRK bank. IREGCON is the control word
c that keep tracks of what subregion we are in.
c
c ====> 1 - Forward extrapolation
c
If (Ireg .eq. NLH) then
c values of versors and position at initial point
IWDATAN = IWDATA + (IREGCON-1) * VWRBLK !elements per block in VWRK
sinthe = sqrt(1.-RW(IWDATAN+VWTZFH)**2)
Txi = RW( IWDATAN + VWTXFH ) / sinthe
Tyi = RW( IWDATAN + VWTYFH ) / sinthe
Exi = -Tyi
Eyi = Txi
Kur = RW( IWDATAN + VWCUFH )
Rho = 1/(Convers * Bfield * Kur)
rxi = RW( IWDATAN + VWXFH )
ryi = RW( IWDATAN + VWYFH )
rzi = RW( IWDATAN + VWZFH )
c values of versors and position at final point
Txf = trpca(1) / sinthe
Tyf = trpca(2) / sinthe
Exf = -Tyf
Eyf = Txf
rxf = rpca(1)
ryf = rpca(2)
rzf = rpca(3)
c
c call transforming routine
c
call VTXTRANS(ISWMFORW)
c
c ====> 2 - If here than backward extrapolation.
c Loop on IREGCON control word (stop loop if IREGCON>IREG)
c
else
do while (IREGCON-IREG .le. 0)
c values of versors and position at initial point
IWDATAN = IWDATA + (IREGCON-1) * VWRBLK !elements per block in VWRK
sinthe = sqrt(1.-RW(IWDATAN+VWTZFH)**2)
Txi = RW( IWDATAN + VWTXFH ) / sinthe
Tyi = RW( IWDATAN + VWTYFH ) / sinthe
Exi = -Tyi
Eyi = Txi
Kur = RW( IWDATAN + VWCUFH )
Rho = 1/(Convers * Bfield * Kur)
rxi = RW( IWDATAN + VWXFH )
ryi = RW( IWDATAN + VWYFH )
rzi = RW( IWDATAN + VWZFH )
c values of versors and position at final point
IREGCON = IREGCON + 1
if (IREGCON .gt. IREG ) then ! end of extrapolation
Txf = trpca(1) / sinthe
Tyf = trpca(2) / sinthe
Exf = -Tyf
Eyf = Txf
rxf = rpca(1)
ryf = rpca(2)
rzf = rpca(3)
call VTXTRANS(Iswmback)
goto 999
else
IWDATAN = IWDATA + (IREGCON-1) * VWRBLK !elements per block in VWRK
sinthe = sqrt(1.-RW(IWDATAN+VWTZFH)**2)
Txf = RW( IWDATAN + VWTXFH ) / sinthe
Tyf = RW( IWDATAN + VWTYFH ) / sinthe
Exf = -Tyf
Eyf = Txf
rxf = RW( IWDATAN + VWXFH )
ryf = RW( IWDATAN + VWYFH )
rzf = RW( IWDATAN + VWZFH )
call VTXTRANS(Iswmback)
c --> Here we are crossing a surface!!!
c So add Multiple Scattering effect in pion mass hypothesis.
c get total momentum and energy
Ptot = 1./abs(Kur)/sinthe
Etot = sqrt(PION_MASS**2 + Ptot**2)
c get angle OMEGA, between momentum versor and surface, and eff. thickness
cosOMEGA = sinthe*(rxf*txf+ryf*tyf)/sqrt(rxf**2+ryf**2)
if (IREGCON .eq. NVC) then
tEFF = tDC / abs(cosOMEGA) ! eff.thickness in rad.len.
sEFF = sDC / abs(cosOMEGA) ! eff.thickness in cm.
elseif (IREGCON.eq.NBI .or. IREGCON.eq.NBU ) then
tEFF = tBP / abs(cosOMEGA)
sEFF = sBP / abs(cosOMEGA)
else
message = 'Unknown extrapolation region'
call ERLOGR ('VTXERR',ERWARN,0,0,message)
return
endif
c get SIG2phi as [ 0.0136GeV/(Beta*p) * sqrt(tEFF) ]**2
SIG2phi = (0.0136 * Etot)/Ptot**2 * sqrt(tEFF)
SIG2phi = SIG2phi*SIG2phi
c derive SIG2kur, SIG2z and SIG2cottheta
SIG2z = SIG2phi * sEFF**2 / 3.
SIG2cot = SIG2phi / sinthe**4
SIG2kur = SIG2phi * Kur**2 * (1-sinthe**2) / sinthe**2
c update diagonal elements of covariant matrix
V22f = V22f + SIG2z
V33f = V33f + SIG2kur
V44f = V44f + SIG2cot
V55f = V55f + SIG2phi
c --> End of MULTIPLE SCATTERING
c update error matrix for next iteration
V11i = V11f
V12i = V12f
V13i = V13f
V14i = V14f
V15i = V15f
V22i = V22f
V23i = V23f
V24i = V24f
V25i = V25f
V33i = V33f
V34i = V34f
V35i = V35f
V44i = V44f
V45i = V45f
V55i = V55f
endif
enddo ! end of DO WHILE loop
c
c if here than error (from final extrapolation should jump to 999)
c
message = 'Error in error matrix extrapolation'
call ERLOGR ('VTXERR',ERWARN,0,0,message)
return
c
endif ! end of extrapolation
999 continue
c
c Covariant Matrix at Candidate Vertex
c
RW( IWDATA + VWERRV + 0 ) = V11f
RW( IWDATA + VWERRV + 1 ) = V12f
RW( IWDATA + VWERRV + 2 ) = V13f
RW( IWDATA + VWERRV + 3 ) = V14f
RW( IWDATA + VWERRV + 4 ) = V15f
RW( IWDATA + VWERRV + 5 ) = V22f
RW( IWDATA + VWERRV + 6 ) = V23f
RW( IWDATA + VWERRV + 7 ) = V24f
RW( IWDATA + VWERRV + 8 ) = V25f
RW( IWDATA + VWERRV + 9 ) = V33f
RW( IWDATA + VWERRV + 10 ) = V34f
RW( IWDATA + VWERRV + 11 ) = V35f
RW( IWDATA + VWERRV + 12 ) = V44f
RW( IWDATA + VWERRV + 13 ) = V45f
RW( IWDATA + VWERRV + 14 ) = V55f
return
end
c=======================================================================
SUBROUTINE VTXTRANS(Iswim)
c=======================================================================
c
c Routine that perform the matrix transformation
c
c Author: Marco Incagli
c =======
c
c Creation date: 11-11-1997
c ==============
c
c Input:
c ======
c Iswim Swim direction (see VTXFIN include file)
c V11i,..,V55i Initial covariant matrix (15 elements)
c TRANSCOM Common file with initial and final point
c
c Output:
c =======
c V11f,..,V55f Final covariant matrix (15 elements)
c
c=======================================================================
c
c include files
c
$$IMPLICIT NONE
$$INCLUDE 'K$ITRK:VTXFIN.INC'
c
c Local declarations
c
INTEGER Iswim
REAL V21i, V31i, V32i, V41i, V42i, V43i
REAL V51i, V52i, V53i, V54i, rif, Alf
REAL D11, D13, D15, D22, D23, D24, D33, D44, D53, D55
c Parameters in common block (VTXERR - VTXTRANS)
REAL rxi, ryi, rzi, txi, tyi, exi, eyi,
+ rxf, ryf, rzf, txf, tyf, exf, eyf,
+ Rho, Kur,
+ V11i, V12i, V13i, V14i, V15i,
+ V22i, V23i, V24i, V25i, V33i,
+ V34i, V35i, V44i, V45i, V55i,
+ V11f, V12f, V13f, V14f, V15f,
+ V22f, V23f, V24f, V25f, V33f,
+ V34f, V35f, V44f, V45f, V55f
COMMON / TRANSCOM / rxi, ryi, rzi, txi, tyi, exi, eyi,
+ rxf, ryf, rzf, txf, tyf, exf, eyf,
+ Rho, Kur,
+ V11i, V12i, V13i, V14i, V15i,
+ V22i, V23i, V24i, V25i, V33i,
+ V34i, V35i, V44i, V45i, V55i,
+ V11f, V12f, V13f, V14f, V15f,
+ V22f, V23f, V24f, V25f, V33f,
+ V34f, V35f, V44f, V45f, V55f
c
c=======================================================================
c
c simmetrization of input matrix
c
V21i = V12i
V31i = V13i
V32i = V23i
V41i = V14i
V42i = V24i
V43i = V34i
V51i = V15i
V52i = V25i
V53i = V35i
V54i = V45i
c
c get rotation angle (same sign as Rho for FORWARD extrapolation,
c opposite sign for BACKWARD extrapolation)
c rif = distance initial-final point
c
rif = sqrt( (rxi-rxf)**2 + (ryi-ryf)**2 )
Alf = 0.5*rif/rho
If ( abs(Alf).ge.1. ) alf = SIGN(1.,alf)
Alf = Iswim * 2 * ASIN( Alf )
c
c evaluate derivatives matrix (all other terms are 0)
c
D11 = Txi*txf + Tyi*Tyf
D13 = ((rxi-rxf)*Exf + (ryi-ryf)*Eyf) / Kur
D15 = (rxf-rxi)*Txf + (ryf-ryi)*Tyf
D22 = 1.
D23 = (rzi-rzf) / Kur
D24 = Rho*Alf
D33 = 1.
D44 = 1.
D53 = - Alf / Kur
D55 = 1.
c
c Final covariant matrix
c
V11f = (D11*V11i+D13*V31i+D15*V51i)*D11 +
+ (D11*V13i+D13*V33i+D15*V53i)*D13 +
+ (D11*V15i+D13*V35i+D15*V55i)*D15
V12f = (D11*V12i+D13*V32i+D15*V52i)*D22 +
+ (D11*V13i+D13*V33i+D15*V53i)*D23 +
+ (D11*V14i+D13*V34i+D15*V54i)*D24
V13f = (D11*V13i+D13*V33i+D15*V53i)*D33
V14f = (D11*V14i+D13*V34i+D15*V54i)*D44
V15f = (D11*V13i+D13*V33i+D15*V53i)*D53 +
+ (D11*V15i+D13*V35i+D15*V55i)*D55
V22f = (D22*V22i+D23*V32i+D24*V42i)*D22 +
+ (D22*V23i+D23*V33i+D24*V43i)*D23 +
+ (D22*V24i+D23*V34i+D24*V44i)*D24
V23f = (D22*V23i+D23*V33i+D24*V43i)*D33
V24f = (D22*V24i+D23*V34i+D24*V44i)*D44
V25f = (D22*V23i+D23*V33i+D24*V43i)*D53 +
+ (D22*V25i+D23*V35i+D24*V45i)*D55
V33f = D33*V33i*D33
V34f = D33*V34i*D44
V35f = D33*V33i*D53 + D33*V35i*D55
V44f = D44*V44i*D44
V45f = D44*V43i*D53 + D44*V45i*D55
V55f = (D53*V33i+D55*V53i)*D53 +
+ (D53*V35i+D55*V55i)*D55
return
end
[KLOE]
[Offline Doc]
[TRK Files]
Generated with Light on Thu Apr 8 13:00:16 MET DST 1999
.
Mail comments and suggestions.
