+PATCH,$CORR. +DECK,CORR1. Updates version: 2.00/00 to 2.00/04 +REP,*TITLE*,TITLE,1. LOCATER 2.00/04 06/11/95 15.32.16 +REP,$VERSION,V2_00,1-2. *CMZ : 2.00/04 06/11/95 15.32.29 by Christian Voelcker *CMZ : 2.00/03 26/10/95 17.49.17 by Christian Voelcker *CMZ : 2.00/02 26/10/95 17.15.05 by Christian Voelcker *CMZ : 2.00/01 05/09/95 14.44.37 by Christian Voelcker *CMZ : 2.00/00 24/08/95 17.44.50 by Christian Voelcker *-- Author : Christian Voelcker 24/08/95 * *::> VERSION 2.00/04 06/11/95 15.32.16 * (Chr. Voelcker) * remove one of the two SAVE IERRL statements in TCVRHX * * (M. Kunze) * rename ZBFORM to TCFORM, as there is a routine with * the same name in the cern library since this year. * * (C.A. Meyer) TSRPHI: Mark Lakata found a subtle bug in the old * TSGAUS routine, so I have replaced all calls with calls * to RNDNOR. I have then removed TSGAUS from the code. * * (J. Adomeit) * corrected TJDCGT to properly use of ILEN * variable to compute number of hits in bank. * * (H. Stoeck) * fix TCRAFU to use RADON2 on VMS systems * * (Ch. Voelcker) * add a 'select $OS' for VMS systems in the install file * * (Ch. Voelcker) * limit the number of hits per track to 50, * including the pwc or svtx hits. * *::> VERSION 2.00/03 26/10/95 17.49.17 * new multi-vertex fitter tcver3 completed * see cbnote 285 for details. * *::> VERSION 2.00/02 26/10/95 17.15.05 * /TC_RAWS/TSMUSH: * C.A.Meyer: Require that the drift time be no * smaller than 0.0001, rather than 0.00. This will * inhibit an error logged in TJTIME. * *::> VERSION 2.00/01 05/09/95 14.44.37 * remove printing in TCMATR and TCMATS +ADD,$KUMACS,INSTALL,0. *CMZ : 2.00/04 01/11/95 17.12.44 by Christian Voelcker +ADD,$KUMACS,INSTALL,25. * 01/11/95 Christian Voelcker * select $OS for VMS systems * +ADD,$KUMACS,INSTALL,131. Select $OS +ADD,TCCOMMON,TC_MACRO,0. *CMZU: 2.00/03 24/10/95 18.23.17 by M.T. Lakata +REP,TCCOMMON,TC_MACRO,870-882. * V3CONF - TCVER3: minimum conf. level to merge vertices * V3DIST - TCVER3: minimum distance to merge vertices * V3BACK - TCVER3: chi^2 penalty for interverted vertex momenta * V3MNHT - TCVER3: min number of hits on a track * VHMCH - TCVRHX: maximum rchi at any time during fit * VHGCH - TCVRHX: maximum rchi at end of fit * VHCFL - TCVRHX: minumum confidence level prob at end of fit * VHERC - TCVRHX: error scaling factor * VHCT - TCVRHX: rchisq threshold * * The parameters that end in D are the Defaults * INTEGER V3MNHD REAL V3COND, V3DISD, V3BACD REAL VHMCHD,VHGCHD,VHCFLD,VHERCD,VHCTD PARAMETER ( V3COND = 0.01, V3DISD=3.0, V3BACD=100.0 ) PARAMETER ( V3MNHD = 3 ) PARAMETER ( VHMCHD = 50.0, VHGCHD = 3.0, VHCFLD = 3.0 ) PARAMETER ( VHERCD = 1.0, VHCTD = 0.001 ) * INTEGER VERTVR,VFITVR,V3MNHT REAL V3CONF, V3DIST, V3BACK REAL VHMCH,VHGCH,VHCFL,VHERC, VHCT * COMMON /TCVTCM/ & VERTVR, VFITVR, & V3CONF, V3DIST, V3BACK, V3MNHT, & VHMCH, VHGCH, VHCFL, VHERC, VHCT +ADD,TC_MAIN,CJINIT,0. *CMZU: 2.00/03 24/10/95 17.57.34 by M.T. Lakata +REP,TC_MAIN,CJINIT,332-334. V3CONF = V3COND V3DIST = V3DISD V3BACK = V3BACD V3MNHT = V3MNHD VHCFL = VHCFLD VHERC = VHERCD VHGCH = VHGCHD VHMCH = VHMCHD VHCT = VHCTD +REP,TC_MAIN,CJINIT,497-499. * For user with TCVER3 CALL FFKEY('V3CL',V3CONF, 1,'R') ! Min Conf Level to combine 2 verts CALL FFKEY('V3DI',V3DIST, 1,'R') ! max distance bwt 2 vertices CALL FFKEY('V3BA',V3BACK, 1,'R') ! introversion penalty CALL FFKEY('V3MH',V3MNHT, 1,'I') ! minimum number of hits * For use with TCVRHX CALL FFKEY('VHCL',VHCFL, 1,'R') ! min confidence level CALL FFKEY('VHER',VHERC, 1,'R') ! Error scaling factor CALL FFKEY('VHGC',VHGCH, 1,'R') ! maximum chi during fit CALL FFKEY('VHMC',VHMCH, 1,'R') ! maximum chi after fit CALL FFKEY('VHCT',VHCT, 1,'R') ! chi threshold *STOPMTL +REP,TC_MAIN,ZBFORM,0. &DECK,ZBFORM,IF=NEVER. *CMZ : 2.00/04 01/11/95 15.24.14 by Christian Voelcker * replaced by SUBROUTINE TCFORM * this routine needed to be renamed to avoid ambiguities with * a new cernlib function (M. Kunze) +ADD,TC_MAIN,TCFORM,*. &DECK,TCFORM. *CMZ : 2.00/04 01/11/95 14.56.33 by Christian Voelcker *CMZU: 1.99/11 06/07/95 23.17.19 by Curtis A. Meyer *CMZU: 1.99/00 10/06/95 00.17.20 by Christian Voelcker *CMZU: 1.61/00 02/03/94 15.38.06 by Curtis A. Meyer *CMZ : 1.56/00 25/09/93 10.58.50 by Curtis A. Meyer *CMZU: 1.50/16 22/11/91 11.47.09 by Curtis A. Meyer *CMZ : 1.45/06 06/03/91 16.58.36 by Curtis A. Meyer *CMZ : 1.45/00 15/01/91 13.27.57 by Curtis A. Meyer *CMZ : 1.44/03 08/10/90 15.28.56 by Curtis A. Meyer *CMZ : 1.43/02 30/08/90 09.26.36 by Curtis A. Meyer *CMZ : 1.42/00 11/07/90 08.27.25 by Curtis A. Meyer *-- Author : SUBROUTINE TCFORM * * Author: Curtis A Meyer * * Creation date: 21 January 1988 * Modified: 22 November 1991 - Add TPCL banks. * Modified: 25 September 1993 - Add TCVH bank. * Modified: 01 March 1994 - Add new JDC calibration banks. * Modified: 10,June 1995 R.Ouared/Chr. Voelcker * NTVXT, NRSVX arrays added for svtx. * Modified: 01, Nov 1995 C. meyer, Chr. Voelcker * rename the routine to TCFORM, as there appeared * a routine ZBFORM in the cernlib! * * References: The ZEBRA user guide, V.3.53(January 1987) * * Description: This subroutine formats all the the data * structures used in the Crystal Barrel * analysis. The format arrays are stored * in the CBFORM common. It is called once at * the beginning of each run. * * **************************************************************** * &SEQ,IMPNONE. * INTEGER I INTEGER NRJDC(5),NTJDC(5),NTJDS(5),NTCHT(5),NTCLY(5) INTEGER NTCTK(5),NTCSG(5),NTCTR(5),NTCTX(5),NTCVH(5) INTEGER NTCHX(5),NTCHP(5),NTCVX(5),NTCVT(5),NTCVP(5) INTEGER NTPWC(5),NTPCH(5),NTPCL(5),NHTJD(5),NTJRF(5) INTEGER NTJCT(5),NTJXX(5),NTJCP(5),NTJST(5),NTJT0(5) INTEGER NTJCZ(5),NTJCE(5),NTJZ0(5),NTJZL(5),NTJWR(5) INTEGER NTJSL(5),NTJSR(5),NTJSZ(5),NTVXT(5),NRSVX(5) * &CDE,TCLIFT. &CDE,TRKPRM. * * Tracking banks: * DATA NRJDC /4HRJDC, 0, 0, 0,1/ DATA NTJDC /4HTJDC,30,30, 30,2/ DATA NTJDS /4HTJDS, 0, 0, 0,0/ DATA NTCHT /4HTCHT,23,23, 23,2/ DATA NTCLY /4HTCLY, 0, 0, 0,0/ DATA NTCTK /4HTCTK,20,20, 1,2/ DATA NTCSG /4HTCSG, 0, 0,LENTK,0/ DATA NTCTR /4HTCTR, 0, 0, 1,2/ DATA NTCTX /4HTCTX, 0, 0,LENTR,0/ DATA NTCVH /4HTCVH, 0, 0, 0,3/ DATA NTCHX /4HTCHX, 0, 0, 1,2/ DATA NTCHP /4HTCHP, 0, 0, 0,3/ DATA NTCVX /4HTCVX,10,10, 1,2/ DATA NTCVT /4HTCVT, 1, 1,LENVX,0/ DATA NTCVP /4HTCVP, 0, 0, 0,0/ DATA NTPWC /4HTPWC, 4, 4, 4,2/ DATA NTPCH /4HTPCH, 0, 0, 0,0/ DATA NTPCL /4HTPCL, 0, 0, 0,0/ DATA NHTJD /4HHTJD,10,10, 3,2/ DATA NTVXT /4HTVXT, 3, 3, 5,2/ DATA NRSVX /4HRSVX, 0, 0, 0,0/ * * Calibration and geometry banks. * DATA NTJRF /4HTJRF, 0, 0, 200,0/ DATA NTJCT /4HTJCT,23,23, 23,2/ DATA NTJXX /4HTJXX, 0, 0, 0,3/ DATA NTJCP /4HTJCP, 0, 0, 460,3/ DATA NTJST /4HTJST, 0, 0, 690,3/ DATA NTJSL /4HTJSL, 0, 0, 1380,3/ DATA NTJSR /4HTJSR, 0, 0, 690,3/ DATA NTJSZ /4HTJSZ, 0, 0, 690,3/ DATA NTJT0 /4HTJT0, 0, 0, 690,2/ DATA NTJCZ /4HTJCZ, 0, 0, 690,3/ DATA NTJCE /4HTJCE, 0, 0, 690,3/ DATA NTJZ0 /4HTJZ0, 0, 0, 690,3/ DATA NTJZL /4HTJZL, 0, 0, 690,3/ DATA NTJWR /4HTJWR, 0, 0, 23,1/ * * __________________________________________ * * Copy the N**** data into the TCLIFT common block. * DO 100 I = 1,5 MRJDC(I) = NRJDC(I) MTJDC(I) = NTJDC(I) MTJDS(I) = NTJDS(I) MTCHT(I) = NTCHT(I) MTCLY(I) = NTCLY(I) MTCTK(I) = NTCTK(I) MTCSG(I) = NTCSG(I) MTCTR(I) = NTCTR(I) MTCTX(I) = NTCTX(I) MTCVH(I) = NTCVH(I) MTCHX(I) = NTCHX(I) MTCHP(I) = NTCHP(I) MTCVX(I) = NTCVX(I) MTCVT(I) = NTCVT(I) MTCVP(I) = NTCVP(I) MTPWC(I) = NTPWC(I) MTPCH(I) = NTPCH(I) MTPCL(I) = NTPCL(I) MHTJD(I) = NHTJD(I) MTVXT(I) = NTVXT(I) MRSVX(I) = NRSVX(I) * MTJRF(I) = NTJRF(I) MTJCT(I) = NTJCT(I) MTJXX(I) = NTJXX(I) MTJCP(I) = NTJCP(I) MTJST(I) = NTJST(I) MTJSL(I) = NTJSL(I) MTJSR(I) = NTJSR(I) MTJSZ(I) = NTJSZ(I) MTJT0(I) = NTJT0(I) MTJCZ(I) = NTJCZ(I) MTJCE(I) = NTJCE(I) MTJZ0(I) = NTJZ0(I) MTJZL(I) = NTJZL(I) MTJWR(I) = NTJWR(I) * 100 CONTINUE * * Set up the format for all the data banks with unusual structures. * CALL MZFORM('TJDS','/4I 12F', MTJDS(5)) CALL MZFORM('TJRF','2I -F', MTJRF(5)) CALL MZFORM('TPCH','/2I 8F', MTPCH(5)) CALL MZFORM('TPCL','/3I 9F', MTPCL(5)) CALL MZFORM('TCLY','/8I 15F', MTCLY(5)) CALL MZFORM('TCSG','7I -F', MTCSG(5)) CALL MZFORM('TCTX','6I 1F 1I -F',MTCTX(5)) CALL MZFORM('TCVT','4I 10F', MTCVT(5)) CALL MZFORM('TCVP','/2I 11F', MTCVP(5)) CALL MZFORM('RSVX','/2I 8F', MRSVX(5)) * RETURN END +ADD,TC_RAWS,TJDCGT,0. *CMZ : 2.00/04 01/11/95 14.50.28 by Christian Voelcker +ADD,TC_RAWS,TJDCGT,133. * **> 01/11095 J.Adomeit/C. Meyer * corrected code to properly use of ILEN * variable to compute number of hits in bank. * +REP,TC_RAWS,TJDCGT,323-333. NHITS = (JBYT(IQ(LRJDC+5),1,16)-2 )/(MAX(2,ILEN)*2) &SELF,IF=ALT,SUN,VAX. ILEN = MAX(2,IBITS(IQ(LRJDC+2),16,16)) LTIM = MAX(1,IBITS(IQ(LRJDC+2), 0,16)) IF(IBITS(IQ(LRJDC+3),16,16).NE.0) LGT0 = .TRUE. NHITS = (IBITS(IQ(LRJDC+5),0,16)-2 )/(MAX(2,ILEN)*2) &SELF,IF=DECS,IBM. ILEN = MAX(2,JBYT(IQ(LRJDC+2),17,16)) LTIM = MAX(1,JBYT(IQ(LRJDC+2), 1,16)) IF(JBYT(IQ(LRJDC+3),17,16).NE.0) LGT0 = .TRUE. NHITS = (JBYT(IQ(LRJDC+5),1,16)-2 )/(MAX(2,ILEN)*2) +ADD,TC_RAWS,TSMUSH,0. *CMZU: 2.00/04 06/09/95 07.53.23 by Curtis A. Meyer +REP,TC_RAWS,TSMUSH,19. * Called By: TSMINI +ADD,TC_RAWS,TSMUSH,31. * 06.9.95: C.A.Meyer: Require that the drift time be no * smaller than 0.0001, rather than 0.00. This will * inhibit an error logged in TJTIME. +REP,TC_RAWS,TSMUSH,112-117. Q(KTJDC+5) = MAX(0.0001,Q(KTJDC+5)*(1.00+DIST)) END IF * * Add the global time_zero jitter to the hit. * Q(KTJDC+5) = MAX(0.0001,Q(KTJDC+5)+TIM0MC) +ADD,TC_RAWS,TSRPHI,0. *CMZU: 2.00/04 13/10/95 08.30.43 by Curtis A. Meyer +ADD,TC_RAWS,TSRPHI,17. * October 13,1995 C.A. Meyer * Mark Lakata found a subtle bug in the old TSGAUS routine, so * I have replaced all calls with calls to RNDNOR. I have then * removed TSGAUS from the code. * +ADD,TC_RAWS,TSRPHI,27. REAL RNDNOR EXTERNAL RNDNOR * +REP,TC_RAWS,TSRPHI,86-89. TSRPHI=C(2)+RNDNOR(0.0,C(3)) * use first gaussian ELSE TSRPHI=C(5)+RNDNOR(0.0,C(6)) +REP,TC_RAWS,TSRPHI,95-97. TSRPHI=CB(2)+RNDNOR(0.0,CB(3)) ELSE TSRPHI=CB(5)+RNDNOR(0.0,CB(6)) +REP,TC_RAWS,TSRPHI,112-150. TSRPHI=TSRPHI+C(2)+RNDNOR(0.0,C(3)) ELSE TSRPHI=TSRPHI+C(5)+RNDNOR(0.0,C(6)) ENDIF ENDIF RETURN END * +ADD,TC_FFUZ,TCMATR,0. *CMZ : 2.00/01 05/09/95 14.42.20 by Christian Voelcker +ADD,TC_FFUZ,TCMATR,27. C C Modifications: 05/09/95 M. Benayoun/ Chr. Voelcker C remove printing. +DEL,TC_FFUZ,TCMATR,54-59. +REP,TC_FFUZ,TCMATR,107-126. A12=(Y1-Y2)/(X1-X2) IF(DABS(A12).LT.1.0D-09) A12=1.0D-09 B12=Y1-A12*X1 C12=(A12+1.D0/A12)*(X1+X2)/2.D0+B12 A23=(Y2-Y3)/(X2-X3) IF(DABS(A23).LT.1.0D-09) A23=1.0D-09 B23=Y2-A23*X2 C23=(A23+1.D0/A23)*(X2+X3)/2.D0+B23 +DEL,TC_FFUZ,TCMATR,142. +DEL,TC_FFUZ,TCMATR,183-190. +DEL,TC_FFUZ,TCMATR,199-203. +ADD,TC_FFUZ,TCMATS,0. *CMZ : 2.00/01 05/09/95 14.44.28 by Christian Voelcker +ADD,TC_FFUZ,TCMATS,10. C Modifications : 05/09/95 M. Benayoun/ Chr. Voelcker C remove pronting C +DEL,TC_FFUZ,TCMATS,54-59. +REP,TC_FFUZ,TCMATS,107-126. A12=(Y1-Y2)/(X1-X2) IF(DABS(A12).LT.1.0D-09) A12=1.0D-09 B12=Y1-A12*X1 C12=(A12+1.D0/A12)*(X1+X2)/2.D0+B12 A23=(Y2-Y3)/(X2-X3) IF(DABS(A23).LT.1.0D-09) A23=1.0D-09 B23=Y2-A23*X2 C23=(A23+1.D0/A23)*(X2+X3)/2.D0+B23 +DEL,TC_FFUZ,TCMATS,142. +DEL,TC_FFUZ,TCMATS,184-191. +DEL,TC_FFUZ,TCMATS,200-204. +REP,TC_RADON2,00_PATCH,0-1. &PATCH,TC_RADON2,IF=CCODE. *CMZ : 2.00/04 16/09/95 22.17.27 by Holger Stoeck +ADD,TC_RADON2,TCRAFU,0. *CMZU: 2.00/04 19/09/95 01.22.17 by Holger Stoeck +ADD,TC_RADON2,TCRAFU,44. &SELF,IF=VMS. #include &SELF. +REP,TC_RADON2,TCRAFU,267-274. &SELF,IF=VMS. $DESCRIPTOR(output,""); output.dsc$a_pointer="SET PRINTOUT"; output.dsc$w_length=strlen("SET PRINTOUT"); mnexcm_(fcnfit_,&output,&arglis,&narg,&ierflg,&futil); &SELF,IF=-VMS. strcpy(output,"SET PRINTOUT"); mnexcm_(fcnfit_,output,&arglis,&narg,&ierflg,&futil,strlen(output)); &SELF. arglis = 0 ; narg = 0 ; futil = 0 ; &SELF,IF=VMS. $DESCRIPTOR(output,""); output.dsc$a_pointer="SET NOWARNINGS"; output.dsc$w_length=strlen("SET NOWARNINGS"); mnexcm_(fcnfit_,&output,&arglis,&narg,&ierflg,&futil); &SELF,IF=-VMS. strcpy(output,"SET NOWARNINGS"); mnexcm_(fcnfit_,output,&arglis,&narg,&ierflg,&futil,strlen(output)); &SELF. +REP,TC_RADON2,TCRAFU,282-283. &SELF,IF=VMS. $DESCRIPTOR(strategy,""); strategy.dsc$a_pointer="SET STRATEGY"; strategy.dsc$w_length=strlen("SET STRATEGY"); mnexcm_(fcnfit_,&strategy,&arglis,&narg,&ierflg,&futil); &SELF,IF=-VMS. strcpy(strategy,"SET STRATEGY"); mnexcm_(fcnfit_,strategy,&arglis,&narg,&ierflg,&futil,strlen(strategy)); &SELF. +REP,TC_RADON2,TCRAFU,294-301. &SELF,IF=VMS. $DESCRIPTOR(pnam,""); pnam.dsc$a_pointer="KAPPA"; pnam.dsc$w_length=strlen("KAPPA"); mnexcm_(&nprm,&pnam,&vstrt,&stp,&zero,&zero,&ierflg); &SELF,IF=-VMS. strcpy(pnam,"KAPPA") ; mnparm_(&nprm,pnam,&vstrt,&stp,&zero,&zero,&ierflg,strlen(pnam)); &SELF. nprm = 2 ; vstrt = param[1] ; stp = 1000.0 * steps[1]; &SELF,IF=VMS. $DESCRIPTOR(pnam,""); pnam.dsc$a_pointer="PHI"; pnam.dsc$w_length=strlen("PHI"); mnexcm_(&nprm,&pnam,&vstrt,&stp,&zero,&zero,&ierflg); &SELF,IF=-VMS. strcpy(pnam,"PHI") ; mnparm_(&nprm,pnam,&vstrt,&stp,&zero,&zero,&ierflg,strlen(pnam)); &SELF. +REP,TC_RADON2,TCRAFU,309-310. &SELF,IF=VMS. $DESCRIPTOR(chcom,""); chcom.dsc$a_pointer="MIGRAD"; chcom.dsc$w_length=strlen("MIGRAD"); mnexcm_(fcnfit_,&chcom,&arglis,&narg,&ierflg,&futil); &SELF,IF=-VMS. strcpy(chcom,"MIGRAD"); mnexcm_(fcnfit_,chcom,&arglis,&narg,&ierflg,&futil,strlen(chcom)); &SELF. +ADD,TC_CIRC,TXTGET,0. *CMZU: 2.00/03 24/10/95 18.29.56 by M.T. Lakata +REP,TC_CIRC,TXTGET,44. * Added PWC/SVX loading for vertex fitting +ADD,TC_HELX,TCHXLD,0. *CMZ : 2.00/04 01/11/95 17.29.39 by Christian Voelcker *CMZU: 2.00/03 24/10/95 17.59.01 by M.T. Lakata +ADD,TC_HELX,TCHXLD,129. * 01 November, 1995 Chr. Voelcker * limit the number of hits per track to 50, * including the pwc or svtx hits. * +ADD,TC_HELX,TCHXLD,181. * * ERRLOG variables * INTEGER IMAXHIT * +REP,TC_HELX,TCHXLD,228. * MTL: I commented this out C IF(NPTS .LE. 5) IERR = 10 +ADD,TC_HELX,TCHXLD,507. IF(NPTS.GT.50) THEN CALL ERRLOG(IMAXHIT, & ' maximum hits on track exceeded (50) ') NPTS = 50 ENDIF +ADD,TC_HELX,TCHXLD,530. IF(NPTS.GT.50) THEN CALL ERRLOG(IMAXHIT, & ' maximum hits on track exceeded (50) ') NPTS = 50 ENDIF +ADD,TC_VERT,TCVER3,0. *CMZU: 2.00/03 24/10/95 18.01.14 by M.T. Lakata +REP,TC_VERT,TCVER3,10-13. * reference CBNote 285 +REP,TC_VERT,TCVER3,27-32. * CALL FFKEY('V3CL',V3CONF, 1,'I') ! Min Conf Level to combine 2 verts * CALL FFKEY('V3DI',V3DIST, 1,'I') ! max distance bwt 2 vertices * CALL FFKEY('V3BA',V3BACK, 1,'I') ! introversion penalty * &SEQ,IMPNONE. &SEQ,CBLINK. &SEQ,TCLIFT. &SEQ,TCPRMS. &SEQ,TCSTAT. &SEQ,TCCUTS. &SEQ,CBUNIT. &SEQ,TRKPRM. &SEQ,CFCOMS. &SEQ,USKEYS. &SEQ,TCVTCM. +ADD,TC_VERT,TCVER3,44. * NHITS number of hits on a track +REP,TC_VERT,TCVER3,53-65. INTEGER NTRAK,IVRTX,ITRAK,IERR,I,J,K,NHITS INTEGER IV,IV1,IV2 INTEGER NV,CIRERR,HLXERR * * NTR number of tracks to try for fit * TRKS list of tracks to try * NDROP number of tracks dropped by fit * IDROP list of dropped tracks * VALID true = use track, false=ignore track INTEGER NTR,TRKS(MAXTRK) LOGICAL VALID(MAXTRK) +REP,TC_VERT,TCVER3,72. INTEGER ITCTR,ITCVT1,ITCVT2,ITCVP,ITCVH +REP,TC_VERT,TCVER3,96-98. INTEGER NPLUS, CHARGE,NMINUS, NVLIST,NJUNK EQUIVALENCE (NTYPE(1),NMINUS),(NTYPE(2),NJUNK),(NTYPE(3),NPLUS) +ADD,TC_VERT,TCVER3,114. * NDF number of degrees of freedom +REP,TC_VERT,TCVER3,124. INTEGER BESTCF,ITCVT,NUNUSE,UNUSED(MAXTRK),NDF +REP,TC_VERT,TCVER3,159-175. CHARACTER*80 CERR * * * for use with TCVRTX INTEGER NDROP,IDROP(MAXTRK) LOGICAL FIRST SAVE FIRST * DATA IERRL /10*0/ DATA FIRST /.TRUE./ +REP,TC_VERT,TCVER3,186-200. IF (V3CONF.NE.V3COND.OR. & V3DIST.NE.V3DISD.OR. & V3BACK.NE.V3BACD.OR. & V3MNHT.NE.V3MNHD) THEN WRITE(LLOG,295 ) V3CONF,V3COND,V3DIST,V3DISD, & V3BACK,V3BACD,V3MNHT,V3MNHD ELSE WRITE(LLOG,294) V3CONF, V3DIST, V3BACK, V3MNHT ENDIF 294 FORMAT( & ' *** TCVER3 Initialized',/, & ' V3CONF is default ',F,/, & ' V3DIST is default ',F,/, & ' V3BACK is default ',F,/, & ' V3MNHT is default ',I2) 295 FORMAT( & ' *** TCVER3 Initialized with changed parameters',/ & ' V3CONF is now ',F,' default was ',F,/, & ' V3DIST is now ',F,' default was ',F,/, & ' V3BACK is now ',F,' default was ',F,/, & ' V3MNHT is now ',I2,' default was ',I2) +REP,TC_VERT,TCVER3,231-245. NHITS = IQ(ITCTR+1) IERR = IQ(ITCTR+6) * * Make sure that we ignore bad tracks * maybe these should be totally ignored... * allowed circle errors: * 0 = ok * 1 = no iterations * 2 = chi^2 too big, but tcsplt says OK * 3 = chi^2 converged but above cutoff * allow helix errors * 00 = ok * 10 = only 3 points * * Also, a minimum number of hits is required, V3MNHT * * note that the circle fit could have a large chi^2 (error=3), but * if the helix fit is ok, then the track is accepted. * CIRERR = MOD(IERR,10) HLXERR = MOD(IERR/10,10) IF (CIRERR.LE.3.AND.HLXERR.LE.1.AND.NHITS.GE.V3MNHT) THEN CHARGE = NINT(Q(ITCTR+7)) NTYPE(CHARGE+2) = NTYPE(CHARGE+2) +1 TYPELI(NTYPE(CHARGE+2),CHARGE+2) = ITRAK VALID(ITRAK) = .TRUE. ELSE NJUNK = NJUNK + 1 TYPELI(NJUNK,2) = ITRAK VALID(ITRAK) = .FALSE. ENDIF 100 CONTINUE IF (NMINUS.EQ.0.OR.NPLUS.EQ.0) THEN BESTCF = 1 NPAIR(BESTCF) = 0 COMASK(BESTCF) = 0 GOTO 500 ENDIF +REP,TC_VERT,TCVER3,275. IF(IERR.GT.100) THEN +REP,TC_VERT,TCVER3,390-397. NDF = 0 DO 410 J=1,NPAIR(I) ITCVT1 = XXTCVT(COLIST(J,I)) CHI = CHI + Q(ITCVT1+14) NDF = NDF + IQ(ITCVT1+4) IF (V3BACK.GT.0) THEN DO 420 K=J+1,NPAIR(I) ITCVT2 = XXTCVT(COLIST(K,I)) * optional step: if momentum of vertex is pointing towards * another vertex, then add a penalty of V3BACK to chi +REP,TC_VERT,TCVER3,404-418. ITCVP = LQ(ITCVT1-1) PX = Q(ITCVP+4)+Q(ITCVP+4+LENVP) PY = Q(ITCVP+5)+Q(ITCVP+5+LENVP) PZ = Q(ITCVP+6)+Q(ITCVP+6+LENVP) IF (X*PX+Y*PY+Z*PZ.LT.0.0) CHI=CHI+V3BACK ITCVP = LQ(ITCVT2-1) PX = Q(ITCVP+4)+Q(ITCVP+4+LENVP) PY = Q(ITCVP+5)+Q(ITCVP+5+LENVP) PZ = Q(ITCVP+6)+Q(ITCVP+6+LENVP) IF (X*PX+Y*PY+Z*PZ.GT.0.0) CHI=CHI+V3BACK 420 CONTINUE ENDIF 410 CONTINUE IF (NDF.LE.0) NDF = 1 CHI = CHI/NDF +DEL,TC_VERT,TCVER3,426-427. +REP,TC_VERT,TCVER3,448. IF(IBITS(COMASK(BESTCF),I,1).EQ.0.AND. & VALID(I)) THEN +REP,TC_VERT,TCVER3,515. IF(R.GT.V3DIST) GOTO 850 +REP,TC_VERT,TCVER3,522. & 2.0*( +REP,TC_VERT,TCVER3,533. & 2.0*( +REP,TC_VERT,TCVER3,555-561. IF (CHI.GT.0.0) THEN P = PROB(CHI,3) ELSE P = 0.0 ENDIF * if vertices are not close together, skip this part IF (P.GE.V3CONF) THEN +REP,TC_VERT,TCVER3,574-576. XGUES(1) = 0.0 XGUES(2) = 0.0 XGUES(3) = 0.0 +REP,TC_VERT,TCVRHX,1-17. *CMZ : 2.00/04 27/10/95 11.45.17 by Christian Voelcker *CMZ : 27/10/95 11.42.49 by Christian Voelcker * remove one of the two SAVE IERRL statements *CMZU: 2.00/03 24/10/95 18.02.38 by M.T. Lakata *CMZU: 1.99/22 10/08/95 20.49.51 by Christian Voelcker *-- Author : Mark Lakata 10/08/95 SUBROUTINE TCVRHX(NTRKS,NTRAK,ITCVT1,IERR1,XGUES) * * Author: Mark Lakata * * Creation date: 05 October 1995. * * References: Siegmund Brandt, Data Analysis, Chapter 9 * editions: 1970, 1976 or 1992 * * Description: This routine will attemps to refit the tracks * listed in NTRAK, with the extra constraint that they * all pass through a common vertex. It does a similar * job as TCVRTX, returning the same zebra banks and * error codes(and 900). However, instead of simply fitting +ADD,TC_VERT,TCVRHX,24. * helix(alpha,tanl,phi) = 3 per track * * For full documentation consult the CB Note 285 by Mark Lakata * * * Called by: TCVER3 * * Subroutines Referenced: * DINV (cernlib F010). * TXTGET (locater) gets hits * MZLIFT (cernlib, zebra) * ERRLOG (cboff) error message logging * * Error codes returned in IERR1 and also stored in zebra banks: * * IERR1 = 0 No error encountered. Normal completion. * IERR1 = 100 Only one track to this vertex. * IERR1 = * IERR1 = 300 Converged to too large a value. * IERR1 = 400 Reached iteration limit without convergence. * IERR1 = 500 Began to diverge. * IERR1 = 600 There are fewer than 2 tracks, no vertex * fit was possible. * IERR1 = 700 Inversion of a singular matrix, routine * failed to extract a vertex. * IERR1 = 800 (only tcvrtx) * IERR1 = 900 Software Error * * All errors 300 and above are considered "bad" and should not be * used in regular data analysis. It is up to the user to check the * return value of ERR1. * * Zebra banks created: * * The format of the TCVT bank is: * * Word Type Description * ================================================ * LQ(-1) I*4 Pointer to the TCVP data bank * ------------------------------------------------ * +1 I*4 Number of tracks to this vertex. * +2 I*4 Number of neutrals. * +3 I*4 Error code of fit. * +4 I*4 Number of degrees of freedom. * +5 R*4 x-vertex (cm). * +6 R*4 y-vertex (cm). * +7 R*4 z-vertex (cm). * +8 R*4 covarianve matrix for vertex. * * +14 R*4 Chisquare of the fit. * * The format of the TCVP bank is: * * Word Type Description * ------------------------------------------------- * +1 I*4 Track number in TCTR * +2 I*4 Quality word. * +3 R*4 Charge of track. * +4 R*4 x-momentum (MeV/c) * +5 R*4 y-momentum (MeV/c) * +6 R*4 z-momentum (MeV/c) * +7 R*4 Energy (assuming pion mass) (Mev). * +8 R*4 Covariance matrix for the momentum * +9 R*4 " * +10 R*4 " | 8 | | s_xx | * +11 R*4 " | 9 10 | = | s_xy s_yy | * +12 R*4 " | 11 12 13 | | s_xz s_yz s_zz | * +13 R*4 " * === Repeat the last 13 for every track at this vertex === * * * Bank Name : TCVH "T"ra"C"king,"V"ertex fitted "H"elix * PREVIOUS bank link : ITCVP i.e. ITCVH = LQ(ITCVP) * * TCVH bank format: * * Word Type Description equivalent to TCTX word * ------------------------------------------- * +1 R*4 r0 +11 -- * +2 R*4 z0 +12 | * +3 R*4 alpha +13 | * +4 R*4 tan lambda +14 |- first track at vertex * +5 R*4 psi0 +15 | * +6 R*4 s +16 | * +7 R*4 Chi**2 +17 -- * The structure then repeats every 7=LENVH words for the rest of the * tracks at that vertex. i.e. Offset +8 is the r0 of track 2. * * ******************************************************************** * &SEQ,IMPNONE. &SEQ,PI2PI. &SEQ,CBUNIT. &SEQ,TCPRMS. &SEQ,TCCUTS. &SEQ,TCSTAT. &SEQ,CBLINK. &SEQ,TRKPRM. &SEQ,TCLIFT. &SEQ,CBSTOP. &SEQ,TCVTCM. * * CONSTANTS: -------------------------------------------------------- * * MAXTRK maximum number of tracks * MAXPT maximum number of hits on a track * MAX2PT MAXPT*2 * MAX3PT MAXPT*3 * INTEGER MAXTRK,MAXPT,MAX2PT,MAX3PT,MAXPAR PARAMETER (MAXTRK=10,MAXPT=30,MAX2PT=MAXPT*2,MAX3PT=MAXPT*3) PARAMETER (MAXPAR=3+3*MAXTRK) * ALPHAX offset of alpha parameter in X() * TANLAX offset of tanlambda parameter in X() * PHI0XX offset of phi0 parameter in X() * XA = x of 1st point * XB = x of 2nd point * YA = y of 1st point * YB = y of 2nd point * ZA = z of 1st point * ZB = z of 2nd point INTEGER ALPHAX,TANLAX,PHI0XX,VX,VY,VZ PARAMETER (ALPHAX=1,TANLAX=2,PHI0XX=3) PARAMETER (VX=1,VY=2,VZ=3) INTEGER XA,XB,YA,YB,ZA,ZB PARAMETER (XA = 1, XB = 4, YA = 2, YB = 5,ZA=3,ZB=6) * SQMPI Mass of the pion squared (MeV). REAL SQMPI PARAMETER (SQMPI=19479.03) * Passed variables:------------------------------------------ * * NTRKS number of tracks to fit to this vertex. * NTRAK array containing the track numbers of the NTRKS. * ITCVT vertex bank to fill (already exists) * IERR returned error code. * XGUES is a guess to the vertex psoition. * INTEGER NTRKS,NTRAK(*),ITCVT1,IERR1 REAL XGUES(3) * Internal variables:---------------------------------------- * IT number of the present track. * IPT number of the present hit * ITER iteration number. * IWORK working area for the DINV routine. * IOFF offset caused by tracks, 0,3,6 * IERR error code of various things * NPAR number of fitting parameters * NPTS number of points in the tracks (for "raw" hits) * NDF number of degrees of freedom * TOTPT total number of points of all tracks used in fitting INTEGER IT,IPT,ITER,IERR INTEGER IOFF,IWORK(MAXPAR) INTEGER NPAR,NDF,NPTS(MAXTRK),TOTPT * ITCHX pointer to helix bank * ITCVP pointer to track subbank of ITCVT * ITCVH pointer to helix subbank of ITCVT * ITCTR Reference link to the TCTR data bank for track IT. * ITCVT pointer to vertex bank INTEGER ITCHX,ITCVP, ITCVH, ITCTR, ITCVT * These are all used in the call to TXGET REAL R1(MAXPT),X1(MAXPT),Y1(MAXPT),Z1(MAXPT), & D1(MAXPT),D2(MAXPT),D3(MAXPT) REAL AMPL(MAXPT),AMPR(MAXPT),SR(MAXPT),SP(MAXPT) INTEGER ILAY(MAXPT),ISC(MAXPT),IRS(MAXPT),TIM(MAXPT) * I,J,K,II,JJ,KK loop indicies. * J1,J2,KX,KY,KZ,IX,IY,IZ offsets INTEGER I,J,K,J1,J2,KX,KY,KZ,IX,IY,IZ,JJ,KK * These are all temporary values (the names are obvious) * BETA beta at each (x,y,z) point. * SINB "sin(beta_i)" = * COSB "cos(beta_i)" = * RP temp value * XX,YY,ZZ position of a hit * SGN sign direction parameter. * INVALP 1/alpha * FACT0 temporary factor * FACT1 temporary factor DOUBLE PRECISION ALPHA,TANLA,PHI0,SGN1,COSPHI,SINPHI,C2S2B DOUBLE PRECISION R0, Z0, INVALP,FACT0,FACT1 DOUBLE PRECISION BETA,SINB,COSB,RP,XX,YY,ZZ,SGN(MAXTRK) * * CHISQ chisquare of the total event * OLDCH chisquare from the previous iteration (of vertex) * DELCH change in a chisquare * OLDBET The old value of beta * RCHISQ reduced chi squared (chi/ndf) DOUBLE PRECISION CHISQ ,OLDCH ,DELCH, OLDBET, RCHISQ * X Contains the values of the 6 helix parameters being fit. * XSI is a correction to X. * Y Contains the (x,y,z) triplets for each point on the track. * Y0 The initial values of Y * GYINV is the covariance matrix for the (x,y,z) points. * A matrix of derivatives of the constraints by fit values. * B matrix of derivatives of the constraints by measurements. * C The evaluations of the equations of constraint. * GB (B[transpose]*GYINV*B)[-1] * ATGBA is the A[trans]*GB*A matrix * ATGBC is A[trans}*GB*C * ATGB A^(transpose) X G_B * DC1..2 C-A*del_X * DY1..3 Corretion to the Y vector. * M1,M2 Used in computing the updated GYINV matrix. * GBxx The (2x2) subunit to the GB matrix * DET determinant of a matrix DOUBLE PRECISION X(MAXPAR),XSI(MAXPAR),X0(MAXPAR), & ATGBA(MAXPAR,MAXPAR),ATGBC(MAXPAR) DOUBLE PRECISION Y(MAX3PT,MAXTRK),Y0(MAX3PT,MAXTRK), & GYINV(MAX3PT,MAXTRK),C(MAX2PT,MAXTRK),GB(3,MAXPT,MAXTRK) DOUBLE PRECISION A(MAX2PT,6,MAXTRK),B(MAX2PT,3,MAXTRK) DOUBLE PRECISION DC1,DC2,DY1,DY2,DY3 DOUBLE PRECISION ATGB(6,2) DOUBLE PRECISION GB11,GB12,GB22,DET * PX,PY,PZ Momentum components at vertex of each track * PP Total momentum * SPP Momentum error matrix (upper triangle) DOUBLE PRECISION PX,PY,PZ,PP,SPP(6) * PSI - old parameterization of PSI * DB - beta difference between vertex and closest approach DOUBLE PRECISION PSI, DB * * AA - coefficients of line equation of first 2 hits * BB - constant of line equation of first 2 hits * AAINV - inverse of AA DOUBLE PRECISION AA(2,2), BB(2), AAINV(2,2) * Error message numbers needed for cboff routine ERRLOG INTEGER IERRL(10) CHARACTER*80 CERR * Confidence level function PROB REAL PROB EXTERNAL PROB * FIRST - true on startup LOGICAL FIRST * this may be uncommented later, if the updated hits are desired * DOUBLE PRECISION M1(3,2),M2(2,2) DATA VHMCH / VHMCHD / DATA VHGCH / VHGCHD / DATA VHCFL / VHCFLD / DATA VHERC / VHERCD / SAVE IERRL, FIRST DATA IERRL /10*0/ DATA FIRST /.TRUE./ * * ______________________________________________________________-- * IF (FIRST) THEN FIRST=.FALSE. * If any changes to defaults, print out warning, otherwise * WELCOME to TCVRHX ! IF ( VHMCH.NE.VHMCHD .OR. & VHGCH.NE.VHGCHD .OR. & VHCFL.NE.VHCFLD .OR. & VHERC.NE.VHERCD) THEN WRITE(LLOG,295 ) & VHMCH,VHMCHD, & VHGCH,VHGCHD, & VHCFL,VHCFLD, & VHERC,VHERCD ELSE WRITE(LLOG,294) VHMCH,VHGCH,VHCFL,VHERC 294 FORMAT(' *** TCVRHX Initialized',/, & ' VHMCH is default ',F,/, & ' VHGCH is default ',F,/, & ' VHCFL is default ',F,/, & ' VHERC is default ',F) 295 FORMAT(' *** TCVRHX Initialized with changed parameters',/ & ' VHMCH is now ',F,' default was ',F,/, & ' VHGCH is now ',F,' default was ',F,/, & ' VHCFL is now ',F,' default was ',F,/, & ' VHERC is now ',F,' default was ',F) ENDIF ENDIF * If ITCVT1 is negative, then it is the address of vertex * if it is positive, then it is an existing vertex number IF (ITCVT1.LT.0) THEN ITCVT = -ITCVT1 ELSE ITCVT = LQ(LTCVX-ITCVT1) ENDIF * ################################################## * # Case of no tracks, or more than MAXTRK tracks # * ################################################## 100 IF((NTRKS .LE. 0).OR.(NTRKS.GT.MAXTRK)) THEN WRITE(CERR,296) NTRKS,MAXTRK 296 FORMAT(' called with ',I3', tracks, out of range ', & '[0,',I2,']') CALL ERRLOG(IERRL(7),CERR) IERR1 = 900 GOTO 5000 ENDIF * ########################## * # Case of only one track # * ########################## * Now treat the special case where there is exactly one track. * This actually has two cases. If the first layer of the found * track is not 1 or 2, then we use the first point on the track * as the vertex. Otherwise the point closest to the point * (0,0,0) is chosesn IF(NTRKS .EQ. 1) THEN * Get the pointer to the track data in the TCHX data bank. ITCTR = LQ(LTCTR-NTRAK(1)) IF(LTCHX.LE.0) THEN ITCHX = 0 ELSE ITCHX = LQ(LTCHX-NTRAK(1)) ENDIF * no iterations ITER = 0 * perfect chisq ??? CHISQ = 0.0D0 * error code 100 = only 1 track at vertex IERR1 = 100 * NDF = 1 (to avoid divide by zero errors) NDF = 1 * zero diagonal elements ??? this can be improved ATGBA(1,2) = 0.0D0 ATGBA(1,3) = 0.0D0 ATGBA(2,3) = 0.0D0 ATGBA(3,2) = 0.0D0 ATGBA(2,1) = 0.0D0 ATGBA(3,1) = 0.0D0 R0 = DBLE(Q(ITCTR+11)) Z0 = DBLE(Q(ITCTR+12)) ALPHA = DBLE(Q(ITCTR+13)) TANLA = DBLE(Q(ITCTR+14)) PSI = DBLE(Q(ITCTR+15)) SGN1 = DBLE(Q(ITCTR+16)) INVALP = 1.0D0/ALPHA * * If the first layer of the track is layer 1,2 or 3, then * take the point of closest approach to (0,0,0) as the * vertex. (also, if there is a pwc hit in the track, take * the point of closest approach.) If the TCHX banks do * not exist, it is impossible to take the first point, * use the dca in this case. * IF((IQ(ITCTR+4).LE.3).OR.(IQ(ITCTR+8).GT.0) & .OR.(LTCHX.LE.0)) THEN * * Here we want to take the point closest to (0,0,0). This * point is given as : * * x = r_0 * sin(psi_0) * y =-r_0 * cos(psi_0) * z = z_0 * X(VX) = R0 * DSIN(PSI) X(VY) = -R0 * DCOS(PSI) X(VZ) = Z0 * * Now create the covariance matrix for this point. * * copy error matrix IF (ITCHX.GT.0) THEN ATGBA(1,1) = DBLE(Q(ITCHX+4)) ATGBA(2,2) = DBLE(Q(ITCHX+5)) ATGBA(3,3) = DBLE(Q(ITCHX+6)) ELSE ATGBA(1,1) = DBLE(Q(ITCTR+18)) ATGBA(2,2) = DBLE(Q(ITCTR+18)) ATGBA(3,3) = DBLE(Q(ITCTR+20)) ATGBA(3,1) = DBLE(Q(ITCTR+19)) ATGBA(3,2) = DBLE(Q(ITCTR+19)) ENDIF BETA = PSI + DPIHLF*SGN1 IF(BETA .LT. 0.D0) BETA = BETA + DPI2 IF(BETA .GE. DPI2) BETA = BETA - DPI2 COSB = DCOS(BETA) SINB = DSIN(BETA) * ELSE * * Here we want to take the coordinate of the first point * on the track. * X(VX) = DBLE(Q(ITCHX+1)) X(VY) = DBLE(Q(ITCHX+2)) X(VZ) = DBLE(Q(ITCHX+3)) * * For the covariance matrix, we want to just use the * measured errors in the points. This is wrong, but it * is a guess. * ATGBA(1,1) = DBLE(Q(ITCHX+4)) ATGBA(2,2) = DBLE(Q(ITCHX+5)) ATGBA(3,3) = DBLE(Q(ITCHX+6)) * * cos(beta_1) = x_1*alpha - (alpha*r_0+s)*sin(psi_0) * sin(beta_1) = y_1*alpha + (alpha*r_0+s)*cos(psi_0) * COSB = ALPHA * X(VX) - ((R0 * ALPHA + SGN1) * DSIN(PSI)) SINB = ALPHA * X(VY) + ((R0 * ALPHA + SGN1) * DCOS(PSI)) * BETA = DATAN2(SINB,COSB) IF(BETA.LT.0.D0) BETA = BETA + DPI2 SINB = DSIN(BETA) COSB = DCOS(BETA) * ENDIF * Fill vertex bank with result GOTO 3000 * ENDIF * * ################################################# * # CASE of between 2 and MAXTRK tracks, (normal) # * ################################################# * TOTPT = 0 IERR1 = 0 * DO 200 IT = 1,NTRKS * * TCTR bank hold helix parameters * IF (NTRAK(IT).EQ.0) THEN CALL ERRLOG(IERRL(2),' NTRAK parameter is corrupt') STOPCB = .TRUE. IERR1 = 900 GOTO 5000 ENDIF * extract circle fit hits(icode=12), use these as CONSTRAINTS CALL TXTGET(12,NTRAK(IT),NPTS(IT),R1,X1,Y1,Z1,D1,D2,D3, & ILAY,ISC,IRS,TIM,AMPL,AMPR,SR,SP) IX = 1 DO 505 IPT=1,NPTS(IT) IY = IX +1 IZ = IX +2 Y0 (IX,IT) = DBLE(X1(IPT)) Y (IX,IT) = Y0(IX,IT) GYINV(IX,IT) = VHERC*DBLE(D1(IPT)**2) Y0 (IY,IT) = DBLE(Y1(IPT)) Y (IY,IT) = Y0(IY,IT) GYINV(IY,IT) = VHERC*DBLE(D2(IPT)**2) Y0 (IZ,IT) = DBLE(Z1(IPT)) Y (IZ,IT) = Y0(IZ,IT) GYINV(IZ,IT) = VHERC*DBLE(D3(IPT)**2) IX = IX + 3 505 CONTINUE TOTPT = TOTPT + NPTS(IT) 200 CONTINUE * The Number of degrees of freedom in this fit is * N_DATA_POINTS - N_PARAMETERS * * However, I am not sure that the N_DATA_POINTS = 3* N_TOTAL_POINTS * (3 for (x,y,z)), since there are only two equations of constraint * per data point, and also, the Chi^2 seem too low for this case * In order to calculate the confidence level, we must have the * correct degrees of freedom. * * for now, I assume N_DATA_POINTS = 2*TOTPT NPAR = 3 + 3*NTRKS NDF = 2*TOTPT - NPAR * Now we need to have an initial guess for the common vertex. * Since TCVER3 will be handing us 2-prong vertices, this will be * a good start... IF (NTRKS.EQ.2) THEN * find intersection of lines formed by first 2 hits of each track * z dependence, just take average of initial hit and origin, * for each track) * need to solve these equations for (a,b) * The equation of a line that passes through (x,y)a and (x,y)b * is * (yb-ya) x + (xa-xb) y = (xa yb - ya xb) * * Thus we need to solve the system, which is the * intersection of two tracks * * (yb1-ya1) x + (xa1-xb1) y = (xa1 yb1 - ya1 xb1) * (yb2-ya2) x + (xa2-xb2) y = (xa2 yb2 - ya2 xb2) * * AA * (x,y) = BB * * (x,y) = AA^inverse * BB * DO 4564 IT =1,2 AA(IT,1) = Y(YB,IT) - Y(YA,IT) AA(IT,2) = Y(XA,IT) - Y(XB,IT) BB(IT) = Y(XA,IT)*Y(YB,IT)-Y(YA,IT)*Y(XB,IT) 4564 CONTINUE * invert AA DET = AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1) IF(DET .EQ. 0.D0) THEN * An attempt was made to invert a singular matrix. CALL ERRLOG(IERRL(4),' AA is singular') IERR1 = 700 GOTO 5000 ENDIF AAINV(1,1) = AA(2,2)/DET AAINV(2,2) = AA(1,1)/DET AAINV(1,2) = -AA(1,2)/DET AAINV(2,1) = -AA(2,1)/DET * X(VX) = DBLE(AAINV(1,1)*BB(1)+AAINV(1,2)*BB(2)) X(VY) = DBLE(AAINV(2,1)*BB(1)+AAINV(2,2)*BB(2)) * now that we know vx, we can extrapolate the Z position at vx * from each track, and then just average the two together. The z position * is not as important as the (x,y) position. IF (AA(1,2).NE.0D0.AND.AA(2,2).NE.0D0) THEN X(VZ) = ( & Y(ZA,1)-(Y(ZB,1)-Y(ZA,1))*(X(VX)-Y(XA,1))/AA(1,2)+ & Y(ZA,2)-(Y(ZB,2)-Y(ZA,2))*(X(VX)-Y(XA,2))/AA(2,2) & )*0.5D0 ELSEIF (AA(1,1).NE.0D0.AND.AA(2,1).NE.0D0) THEN X(VZ) = ( & Y(ZA,1)-(Y(ZB,1)-Y(ZA,1))*(X(VY)-Y(YA,1))/AA(1,1)+ & Y(ZA,2)-(Y(ZB,2)-Y(ZA,2))*(X(VY)-Y(YA,2))/AA(2,1) & )*0.5D0 ELSE X(VZ) = 0.D0 ENDIF IF (X(1)**2+X(2)**2.GT.2.D0.AND.DABS(X(3)).GT.2.D0) THEN DO 111 I=1,3 IF (DABS(X(I)).GT.30.D0.OR. & DABS(X(I)).GT.DABS(Y(I,1)).OR. & DABS(X(I)).GT.DABS(Y(I,2))) THEN * * The guess was bad, but don't give up... try 0,0,0 ! * X(1) = 0.D0 X(2) = 0.D0 X(3) = 0.D0 ENDIF 111 CONTINUE ENDIF ELSE X(VX) = DBLE(XGUES(1)) X(VY) = DBLE(XGUES(2)) X(VZ) = DBLE(XGUES(3)) ENDIF * Initialization. * set Chi-square to a very large value. * OLDCH = 1.0D10 DELCH = 1.0D10 ITER = 1 * Load the helix fit information from the TCTR bank for each * track. * IOFF = 0 DO 800 IT = 1,NTRKS ITCTR = LQ(LTCTR-NTRAK(IT)) IOFF = IOFF + 3 * * initial guesses SGN(IT) = DBLE(Q(ITCTR+16)) * alpha, parameter1 X(ALPHAX+IOFF) = DBLE(Q(ITCTR+13)) * tan (lambda) parameter2 X(TANLAX+IOFF) = DBLE(Q(ITCTR+14)) * phi0, parameter3 RP = (SGN(IT)*DBLE(Q(ITCTR+11))+1/X(ALPHAX+IOFF)) PHI0 = DATAN2( & X(2)-RP*DSIN(DBLE(Q(ITCTR+15))-SGN(IT)*DPIHLF), & X(1)-RP*DCOS(DBLE(Q(ITCTR+15))-SGN(IT)*DPIHLF)) IF (PHI0.LT.0.0) PHI0 = PHI0 + DPI2 X(PHI0XX+IOFF) = PHI0 * 800 CONTINUE **************************************************************************** **************************************************************************** **************************************************************************** **************************************************************************** * Save the initial values, and compare the updated ones- * make sure they don't get out of hand!! DO 905 I=1,NPAR X0(I) = X(I) 905 CONTINUE * beginning of main iteration loop. * control comes back here if the first pass is not good enough 1000 CONTINUE * zero some things CHISQ = 0.0D0 DO 1100 J = 1,NPAR XSI(J) = 0.D0 ATGBC(J) = 0.D0 DO 1120 K=1,NPAR ATGBA(J,K) = 0.D0 1120 CONTINUE 1100 CONTINUE IOFF = 0 DO 2800 IT = 1,NTRKS * Compute repeatedly used quantities in the loop: IOFF = IOFF + 3 ALPHA = X(ALPHAX+IOFF) TANLA = X(TANLAX+IOFF) PHI0 = X(PHI0XX+IOFF) SGN1 = SGN(IT) INVALP = 1.D0/ALPHA COSPHI = DCOS(PHI0) SINPHI = DSIN(PHI0) * Now loop through and load the A and B matricies, and * the C vector, and various products: * loop over all the point on this track, IT DO 2000 IPT = 1,NPTS(IT) * Compute offset pointers to the hit data in the Y vector, * we have x_i = Y(k1), y_i = Y(KY), and z_i = Y(KZ). The * J1,J2 pointers address data in the GB matrix: J2 = 2*IPT J1 = J2-1 KZ = 3*IPT KY = KZ-1 KX = KZ-2 XX = Y(KX,IT) YY = Y(KY,IT) ZZ = Y(KZ,IT) * Calculate sin and cos of beta_i at this point: COSB = ALPHA*(XX-X(VX)) + COSPHI SINB = ALPHA*(YY-X(VY)) + SINPHI * Force beta in the range -pi,pi. * note: by definition, beta = 0 at vertex BETA = DATAN2(SINB,COSB) - PHI0 IF (BETA.GT. DPI) BETA=BETA-DPI2 IF (BETA.LT.-DPI) BETA=BETA+DPI2 BETA = -BETA*SGN1 * This may occur in weird occasions, that the beta of a track * is negative, which means that something is corrupt in the * helix parameters. This seems to be very rare, so I * won't worry about the cause of it now. IF(BETA.LT.-0.1D0) THEN IF (IPT.GT.1) THEN IF (BETA-OLDBET.LT.0.D0) THEN WRITE(CERR,309) SNGL(BETA),SNGL(PHI0), & NINT(SNGL(SGN1)), & IT,IPT,SNGL(DATAN2(SINB,COSB)) 309 FORMAT(' Beta <0', & F10.3,F10.3,I,I2,I2,F10.3) CALL ERRLOG(IERRL(3),CERR) IERR1 = 700 GOTO 5000 ENDIF ENDIF ENDIF OLDBET = BETA * The C vector just contains the evaluations of the 2n * constraint equations. C2S2B = COSB**2 + SINB**2 C(J1,IT) = INVALP * ( C2S2B - 1.D0) C(J2,IT) = X(VZ) - ZZ + INVALP*TANLA*BETA * Compute the derivatives of F1 and F2 with respect to the * measurements, (x,y,z). * dF1/dx = 2 * cos(beta_i) * dF1/dy = 2 * sin(beta_i) * dF1/dz = 0 B(J1,1,IT) = 2.D0 * COSB B(J1,2,IT) = 2.D0 * SINB B(J1,3,IT) = 0.D0 * dF2/dx = sgn * tan(lambda) * sin(beta) * dF2/dy = -sgn * tan(lambda) * cos(beta) * dF2/dz = -1.0 B(J2,1,IT) = SGN1 * TANLA * SINB/C2S2B B(J2,2,IT) = -SGN1 * TANLA * COSB/C2S2B B(J2,3,IT) = -1.D0 * Now compute the GB matrix. GB = (B(transpose)Gy(-1)B)(-1) * This is a symetric matrix, so we only need to compute * three of the four elements. Start with GB(-1) * (remember that B3(J1)=dF1/dz = 0, B3(J2) = dF2/dz = -1). GB11 = GYINV(KX,IT)*B(J1,1,IT)*B(J1,1,IT) & + GYINV(KY,IT)*B(J1,2,IT)*B(J1,2,IT) GB12 = GYINV(KX,IT)*B(J1,1,IT)*B(J2,1,IT) & + GYINV(KY,IT)*B(J1,2,IT)*B(J2,2,IT) GB22 = GYINV(KX,IT)*B(J2,1,IT)*B(J2,1,IT) & + GYINV(KY,IT)*B(J2,2,IT)*B(J2,2,IT) & + GYINV(KZ,IT) * Make sure that the determinant is not zero, and then * invert the matrix. DET = GB11*GB22 - GB12*GB12 IF(DET .EQ. 0.D0) THEN * An attempt was made to invert a singular matrix. CALL ERRLOG(IERRL(4),' GB is singular') IERR1 = 700 GOTO 5000 ENDIF DET = 1.D0/DET GB(1,IPT,IT) = GB22*DET GB(2,IPT,IT) =-GB12*DET GB(3,IPT,IT) = GB11*DET GB11 = GB(1,IPT,IT) GB12 = GB(2,IPT,IT) GB22 = GB(3,IPT,IT) * Compute the elements of the (2n,6) A Matrix. The A matrix * is defined as A(i,j) = d(Fi)/d(Xj) * consult the CBnotes for pretty formatted versions of these * derivatives... * a(j,1) = dF_j / dV_x <-- these are easy, since they are the same as B, * a(j,2) = dF_j / dV_y < / with an overall minus * a(j,3) = dF_j / dV_z ATGBA is singular') IERR1 = 700 GOTO 5000 ENDIF * Now we start updating the X-helix parameters and Y-hits DO 2600 I = 1,NPAR DO 2500 J = 1,NPAR XSI(I) = XSI(I) - ATGBA(I,J) * ATGBC(J) 2500 CONTINUE X(I) = X(I) + XSI(I) 2600 CONTINUE * If the value of alpha has become negative, then change its * sign to positive, and change the value of SGN as well. IOFF = 0 DO 2650 IT=1,NTRKS IOFF= IOFF+3 IF(X(ALPHAX+IOFF) .LT. 0.D0) THEN I = alphax+ioff IF (DSQRT(ATGBA(i,i)).GT.DABS(X(I))) THEN * undo update !! X(I) = X(I) - XSI(I) ENDIF C X(ALPHAX+IOFF) = -X(ALPHAX+IOFF) C SGN(IT) = -SGN(IT) ISTKTC(48) = ISTKTC(48) + 1 ENDIF * The correction to the Y vector. DY = GY(-1) * B^{T} * G_{B} * (C + A XSI) DO 2700 IPT = 1,NPTS(IT) J2 = 2*IPT J1 = J2 - 1 KZ = 3*IPT KY = KZ - 1 KX = KY - 1 * Compute C + AdX DC1 = C(J1,IT) & + A(J1,1,IT)*XSI(1) + A(J1,2,IT)*XSI(2) & + A(J1,3,IT)*XSI(3) + A(J1,4,IT)*XSI(1+IOFF) & + A(J1,5,IT)*XSI(2+IOFF) + A(J1,6,IT)*XSI(3+IOFF) DC2 = C(J2,IT) & + A(J2,1,IT)*XSI(1) + A(J2,2,IT)*XSI(2) & + A(J2,3,IT)*XSI(3) + A(J2,4,IT)*XSI(1+IOFF) & + A(J2,5,IT)*XSI(2+IOFF) + A(J2,6,IT)*XSI(3+IOFF) * Now compute GB * (C+Adx) (2 elements) FACT0 = GB(1,IPT,IT)*DC1 + GB(2,IPT,IT)*DC2 FACT1 = GB(2,IPT,IT)*DC1 + GB(3,IPT,IT)*DC2 * Now compute delta_y: DY1 = GYINV(KX,IT) *(B(J1,1,IT) *FACT0 + B(J2,1,IT)*FACT1) DY2 = GYINV(KY,IT) *(B(J1,2,IT) *FACT0 + B(J2,2,IT)*FACT1) DY3 = GYINV(KZ,IT) *(0.D0 - 1.0D0 *FACT1) * Update the current values of y. Y(KX,IT) = Y(KX,IT) - DY1 Y(KY,IT) = Y(KY,IT) - DY2 Y(KZ,IT) = Y(KZ,IT) - DY3 * Compute the total change in Y wrt the original hit positions (y0) DY1 = Y(KX,IT) - Y0(KX,IT) DY2 = Y(KY,IT) - Y0(KY,IT) DY3 = Y(KZ,IT) - Y0(KZ,IT) * Compute the contribution to chisquare from this point. * chi_square = (B*DY)[transpose]*GB*(B*DY) * Calculate B*DY, (remember that dF1/dz = 0, dF2/dz = -1). FACT0 = B(J1,1,IT)*DY1 + B(J1,2,IT)*DY2 FACT1 = B(J2,1,IT)*DY1 + B(J2,2,IT)*DY2 - DY3 * * CHISQ = CHISQ & + GB(1,IPT,IT)*FACT0*FACT0 & + 2.D0*GB(2,IPT,IT)*FACT0*FACT1 & + GB(3,IPT,IT)*FACT1*FACT1 * end of loop over hits on a track 2700 CONTINUE * end of loop over tracks 2650 CONTINUE * done with all matrix stuff now check the results * Check for convergence. RCHISQ = CHISQ/NDF IF(RCHISQ .LT. VHCT) GOTO 3000 DELCH = DABS(RCHISQ - OLDCH) IF(DABS(DELCH) .LT. VHCT) GOTO 3000 * Exceeded iteration limit: IF(ITER .GT. NVTXTC) THEN IERR1 = 500 GOTO 5000 ENDIF * Check for divergence of the series. IF((ITER.GT.3).AND.(DELCH .LT. 0.D0).AND. & (RCHISQ .GT. OLDCH*1.10D0)) THEN IERR1 = 400 GOTO 5000 ENDIF * Is chisquare just too large? IF(RCHISQ .GT. VHMCH) THEN IERR1 = 400 GOTO 5000 ENDIF DO 3111 I=1,3 IF (DABS(X(I)).GT.30.D0) THEN IERR1 = 400 GOTO 5000 ENDIF 3111 CONTINUE ITER = ITER + 1 OLDCH = RCHISQ * The fit is not good enough. Go back for another * iteration. GOTO 1000 * Normal convergence of the loop. We now need to compute what the * new errors in Y are. First check to see if the convergence * value was too large. 3000 CONTINUE C IF( PROB(SNGL(CHISQ),NDF).LT.VHCFL ) THEN C IERR1 = 300 C GOTO 5000 C ENDIF * Normal convergence of the routine. First copy the fit vertex * information into the vertex bank. * Compute the total change in Y wrt the original hit positions (y0) IF(RCHISQ .GT. VHGCH) IERR1 = 300 * IQ(ITCVT+1) = NTRKS IQ(ITCVT+2) = 0 IQ(ITCVT+3) = IERR1 IQ(ITCVT+4) = NDF Q(ITCVT+5) = SNGL(X(1)) Q(ITCVT+6) = SNGL(X(2)) Q(ITCVT+7) = SNGL(X(3)) * copy error matrix Q(ITCVT+8) = SNGL(ATGBA(1,1)) Q(ITCVT+9) = SNGL(ATGBA(2,1)) Q(ITCVT+10) = SNGL(ATGBA(2,2)) Q(ITCVT+11) = SNGL(ATGBA(3,1)) Q(ITCVT+12) = SNGL(ATGBA(3,2)) Q(ITCVT+13) = SNGL(ATGBA(3,3)) Q(ITCVT+14) = SNGL(CHISQ) * If any diagonal elements of the covariance matrix are 0, this * is nonsense. IF((Q(ITCVT+8).LT.0.0).OR.(Q(ITCVT+10).LT.0.0).OR. & (Q(ITCVT+13).LT.0.0)) THEN CALL ERRLOG(IERRL(6), & ' Diagonal of cov matrix is all zeros') IERR1 = 700 GOTO 5000 ENDIF * ############################################## * # Store the fit data into the TCVP data bank # * ############################################## * Create the TCVP data bank under the TCVT data bank to hold the * fit momentum at the vertex. The bank has NTRKS rows, each of * which is LENVP columns wide. * * HACK/DEBUG: add this feature: check to see if the bank * exists already, if so then use it. also * make sure that the ntrks is the same. MTCVP(4) = LENVP * NTRKS * ITCVT = LQ(LTCVX-IVRTX) CALL MZLIFT(IXTCVX,ITCVP,ITCVT,-1,MTCVP,0) * Create the TCVH data bank next to the TCVP data bank to hold the * newly fitted helix parameters. The bank has NTRKS rows, each of * which is LENVH columns wide. MTCVH(4) = LENVH * NTRKS CALL MZLIFT(IXTCVX,ITCVH,ITCVP,0,MTCVH,0) IOFF = 0 DO 6000 IT = 1,NTRKS IOFF= IOFF+3 * Calculate the new hit errors * * This is almost never going to be used, so it is * just copied from helix fit... If it is ever used, it needs to be * debugged for sure! * ** _Datenanalyse or Data Analysis_, Siegmund Brandt ** p. 183 (1st ed.) ** p. 271? (2nd ed.) ** p. 318 (3rd ed.) ** ** GyINV += M1 * B * GyINV + M1 * M2 * M1^T ** ** M1 = GyINV B^T G_B ** M2 = A ATGBA A^T * DO 4000 I = 1,NPTS(IT) ** * KZ = 3*I * KY = KZ-1 * KX = KY-1 * J2 = 2*I * J1 = J2-1 ** ** The M1 matrix is GYINV * B(transpose) * G_{B}. ** * M1(1,1) = GYINV(KX)*(B(J1,1,IT)*GB(1,I)+B(J2,1,IT)*GB(2,I)) * M1(1,2) = GYINV(KX)*(B(J1,1,IT)*GB(2,I)+B(J2,1,IT)*GB(3,I)) * M1(2,1) = GYINV(KY)*(B(J1,2,IT)*GB(1,I)+B(J2,2,IT)*GB(2,I)) * M1(2,2) = GYINV(KY)*(B(J1,2,IT)*GB(2,I)+B(J2,2,IT)*GB(3,I)) * M1(3,1) = -GYINV(KZ)*GB(2,I) * M1(3,2) = -GYINV(KZ)*GB(3,I) ** ** The M2 matrix is the (2x2) along-the-diagonal element of ** A * Covr_x * A(transpose) ** * M2(1,1) = 0.D0 * M2(2,1) = 0.D0 * M2(2,2) = 0.D0 * DO 3200 J = 1,6 * DO 3100 K = 1,6 * M2(1,1) = M2(1,1)+A(J1,J,IT)*A(J1,K,IT)*ATGBA(K,J,IT) * M2(2,1) = M2(2,1)+A(J1,J,IT)*A(J2,K,IT)*ATGBA(K,J,IT) * M2(2,2) = M2(2,2)+A(J2,J,IT)*A(J2,K,IT)*ATGBA(K,J,IT) * 3100 CONTINUE * 3200 CONTINUE * M2(1,2) = M2(2,1) ** ** Now compute the diagonal elements of M1 * M2 * M1[transpose] ** * DY1 = M1(1,1)*M1(1,1)*M2(1,1) * & + 2.D0*M1(1,1)*M1(1,2)*M2(2,1) * & + M1(1,2)*M1(1,2)*M2(2,2) ** * DY2 = M1(2,1)*M1(2,1)*M2(1,1) * & + 2.D0*M1(2,1)*M1(2,2)*M2(2,1) * & + M1(2,2)*M1(2,2)*M2(2,2) ** * DY3 = M1(3,1)*M1(3,1)*M2(1,1) * & + 2.D0*M1(3,1)*M1(3,2)*M2(2,1) * & + M1(3,2)*M1(3,2)*M2(2,2) ** ** Now compute the diagonal elements of ** GY[-1]*B[trans]*GB*B*GY[-1] ** * DY1 = DY1 - GYINV(KX,IT)*(B(J1,1,IT)*M1(1,1) + B(J2,1,IT)*M1(1,2)) * DY2 = DY2 - GYINV(KY,IT)*(B(J1,2,IT)*M1(2,1) + B(J2,2,IT)*M1(2,2)) * DY3 = DY3 + GYINV(KZ,IT)*M1(3,2) ** ** Now compute the new diagonals along the covariance matrix. ** * GYINV(KX,IT) = GYINV(KX,IT) + DY1 * GYINV(KY,IT) = GYINV(KY,IT) + DY2 * GYINV(KZ,IT) = GYINV(KZ,IT) + DY3 ** ** * extract the momentum information at this point. This can be * done using the fact that the beta for each track has already * been obtained for these tracks. Note that in the case of exactly * one track, we do not improve the covariance matrix. * store the vertex number with the track data * (note that this is probably modified by TCVER3) ITCTR = LQ(LTCTR-NTRAK(IT)) IF (ITCVT1.GT.0) THEN IQ(ITCTR+3) = ITCVT1 ELSE IQ(ITCTR+3) = 0 ENDIF * store the track information with the vertex IQ(ITCVP+1) = NTRAK(IT) IQ(ITCVP+2) = IQ(ITCTR+1)+100*IQ(ITCTR+4)+10000*IQ(ITCTR+6) * Now we want to store the particle information. This is done * For all number of tracks, including exactly one. * Fill the new TCVH bank with the newly fitted helix parameters * convert the 3+3n parameters to 5n style IF (NTRKS.GT.1) THEN SGN1 = SGN(IT) INVALP = 1.D0/X(ALPHAX+IOFF) SINB = DSIN(X(PHI0XX+IOFF)) COSB = DCOS(X(PHI0XX+IOFF)) PSI = DATAN2(X(2)-INVALP*SINB,X(1)-INVALP*COSB) + & SGN1*DPIHLF IF(PSI.LT.0.0D0) PSI = PSI + DPI2 DB = PSI - SGN1*DPIHLF - X(PHI0XX+IOFF) + DPI IF (DB.GT.DPI) DB = DB-DPI2 DB = -SGN1*DB R0 = SGN1*(SQRT((X(1)-INVALP*COSB)**2 & +(X(2)-INVALP*SINB)**2)-INVALP) Z0 = X(3)+X(TANLAX+IOFF)*INVALP*DB ALPHA = X(ALPHAX+IOFF) TANLA = X(TANLAX+IOFF) ENDIF Q(ITCVH+1) = SNGL(R0) Q(ITCVH+2) = SNGL(Z0) Q(ITCVH+3) = SNGL(ALPHA) Q(ITCVH+4) = SNGL(TANLA) Q(ITCVH+5) = SNGL(PSI) Q(ITCVH+6) = SNGL(SGN1) Q(ITCVH+7) = SNGL(CHISQ) * Store the particle charge. Q(ITCVP+3) = SNGL(SGN1)*BSGNTC * Compute the three components of momentum along the track. PP = DBLE(BMGCTC) * INVALP PX = PP * SINB * SGN1 PY = -PP * COSB * SGN1 PZ = PP * TANLA * Store the momentum and energy of the particles. Q(ITCVP+4) = SNGL(PX) Q(ITCVP+5) = SNGL(PY) Q(ITCVP+6) = SNGL(PZ) Q(ITCVP+7) = SQRT(SNGL(PP*PP + PZ*PZ) + SQMPI) * Compute the three by three covariance matrix for this * track. It is stored as: * | 1 | * | 2 3 | * | 4 5 6 | SPP(1) = -PX * INVALP * GYINV( 6,IT) & + PY * GYINV(13,IT) * PX*INVALP & - PX * INVALP * GYINV(13,IT) & + PY * GYINV(15,IT) * PY SPP(2) = -PY * INVALP * GYINV( 6,IT) & - PX * GYINV(13,IT) * PX*INVALP & - PY * INVALP * GYINV(13,IT) & - PX * GYINV(15,IT) * PY SPP(3) = -PY * INVALP * GYINV( 6,IT) & - PX * GYINV(13,IT) * PY*INVALP & - PY * INVALP * GYINV(13,IT) & + PX * GYINV(15,IT) * PX SPP(4) = -PZ * INVALP * GYINV( 6,IT) & - PP * GYINV( 9,IT) * PX*INVALP & - PZ * INVALP * GYINV(13,IT) & - PP * GYINV(14,IT) * PY SPP(5) = -PZ * INVALP * GYINV( 6,IT) & - PP * GYINV( 9,IT) * PY*INVALP & - PZ * INVALP * GYINV(13,IT) & + PP * GYINV(14,IT) * PX SPP(6) = -PZ * INVALP * GYINV( 6,IT) & - PP * GYINV( 9,IT) * PZ*INVALP & - PZ * INVALP * GYINV( 9,IT) & + PP * GYINV(10,IT) * PP * Store the three by three covariance matrix for the momentum in * the TCVP bank. DO 8400 I = 1,6 Q(ITCVP+7+I) = SNGL(SPP(I)) 8400 CONTINUE * Move out to the next row in the bank. ITCVP = ITCVP + LENVP ITCVH = ITCVH + LENVH 6000 CONTINUE * Update the statistics on the error codes. C DO 6200 IT = 1,NTRKS C ITCTR = LQ(LTCTR-IT) C IERR = IQ(ITCTR+6)/10 C ISTKTC(31+IERR) = ISTKTC(31+IERR) + 1 C 6200 CONTINUE RETURN 5000 CONTINUE * FILL IN ERROR CODE * If the NTRKS = 0, then the vertex is a junk vertex, and none of * the value except for the IERR1 will be updated. IF (ITCVT.GT.0) THEN NTRKS = 0 IQ(ITCVT+1) = NTRKS IQ(ITCVT+2) = 0 IQ(ITCVT+3) = IERR1 ENDIF RETURN END +ADD,TC_VERT,TCVRHY,*. &DECK,TCVRHY. *CMZU: 2.00/03 24/10/95 18.18.35 by M.T. Lakata *-- Author : M.T. Lakata 24/10/95 SUBROUTINE TCVRHY(NVERT,NPRONG,IERR) * * ********************************************************************** * * SUBROUTINE TCVRHY * * Author : M.Lakata * Date : Oct 24, 1995 * * Returns number of vertices in NVERT, and the multiplicity * of n-pronged vertices in NPRONG(n). IERR .ne. 0 if there is * any error in any vertex. * * see CB Note 285 for documentation * * Input : - * Output : - * * Commons : /none/ * * ********************************************************************** &SEQ,IMPNONE. &SEQ,CBLINK. INTEGER NVERT,NPRONG(10),IERR INTEGER ITCVT,IV,NTRK NVERT = 0 IERR = 0 DO 100 IV=1,10 NPRONG(IV) = 0 100 CONTINUE IF (LTCVX.LE.0) RETURN DO 200 IV=1,IQ(LTCVX+1) ITCVT = LQ(LTCVX-IV) IF (ITCVT.GT.0) THEN NTRK = IQ(ITCVT+1) IF (NTRK.GT.0) THEN NVERT=NVERT+1 IF (IQ(ITCVT+3).GT.100) IERR = IERR + 1 IF (NTRK.LE.10) THEN NPRONG(NTRK)=NPRONG(NTRK)+1 ELSE IERR = IERR + 10 ENDIF ENDIF ENDIF 200 CONTINUE RETURN END +ADD,TC_VERT,TCVRHZ,*. &DECK,TCVRHZ. *CMZU: 2.00/03 24/10/95 18.18.44 by M.T. Lakata *-- Author : M.T. Lakata 24/10/95 SUBROUTINE TCVRHZ(NPRONG,IVERT,E,PX,PY,PZ,CHRG,MASS,X,T,CHI, & ITCVT,IERR) * * ********************************************************************** * * SUBROUTINE TCVRHZ * * Author : M.Lakata * Date: Oct 24, 1995 * * for a given value of NPRONG and IVERT (the i'th vertex with NPRONG * tracks), the routine returns several values. See CB Note 285. * * Commons : /none/ * * This routine unpacks the vertex banks. See the documentation for * TCVER3/TCVRHX * * input: NPRONG - how many prongs desired * IVERT - Ith vertex with that many prongs * output: E - array of Energies of outgoing tracks (1..NPRONG) * PX - array of momentum x * PY - array of momentum y * PZ - array of momentum z * CHRG - array of charge * MASS - mass of vertex * X - X(3) position of vertex * T - T(4) momentum of vertex (4=energy) * CHI - chi^2 of vertex * ITCVT - pointer to vertex zebra bank * IERR - return code * ********************************************************************** &SEQ,IMPNONE. &SEQ,CBLINK. &SEQ,TRKPRM. INTEGER MAXTRK PARAMETER (MAXTRK=10) INTEGER NPRONG,IVERT,IERR,ITCVP REAL E(MAXTRK),PX(MAXTRK),PY(MAXTRK),PZ(MAXTRK) REAL CHRG(MAXTRK),MASS,X(3),CHI,T(4) INTEGER ITCVT,IV,NTRK,I,J IERR = 1 I = 0 IF (NPRONG.LE.0.OR.NPRONG.GT.MAXTRK) RETURN IF (LTCVX.LE.0) RETURN DO 200 IV=1,IQ(LTCVX+1) ITCVT = LQ(LTCVX-IV) IF (ITCVT.GT.0) THEN NTRK = IQ(ITCVT+1) IF (NTRK.EQ.NPRONG) THEN I = I + 1 IF (I.EQ.IVERT) THEN X(1) = Q(ITCVT+5) X(2) = Q(ITCVT+6) X(3) = Q(ITCVT+7) IF (Q(ITCVT+4).GT.0) THEN CHI = Q(ITCVT+14)/Q(ITCVT+4) ELSE CHI = 0.0 ENDIF T(1) = 0.0 T(2) = 0.0 T(3) = 0.0 T(4) = 0.0 DO 100 J=1,NPRONG ITCVP = LQ(ITCVT-1)+LENVP*(J-1) CHRG(J) = Q(ITCVP+3) PX(J) = Q(ITCVP+4) PY(J) = Q(ITCVP+5) PZ(J) = Q(ITCVP+6) E (J) = Q(ITCVP+7) T(1) = T(1) + PX(J) T(2) = T(2) + PY(J) T(3) = T(3) + PZ(J) T(4) = T(4) + E (J) 100 CONTINUE MASS = T(4)**2 - T(1)**2 - T(2)**2 - T(3)**2 IF (MASS.GT.0.0) MASS = SQRT(MASS) IERR = 0 RETURN ENDIF ENDIF ENDIF 200 CONTINUE RETURN END