+PATCH,$CORR. +DECK,CORR1. Updates version: 2.00/00 to 2.00/05 +REP,*TITLE*,TITLE,0-1. &TITLE. LOCATER 2.00/05 12/12/95 17.56.02 +REP,$VERSION,V2_00,1-4. *CMZU: 2.00/05 12/12/95 17.56.05 by Christian Voelcker *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/05 12/12/95 17.56.02 * CJINIT (C. A. Meyer) 12/12/95 * tiny bug initializing DLZZMC * * TCHXLD (F. H. Heinsius) 28/11/95 * reconstructing data without svtx and pwc still produces * crashes... (why??) * reorganize the handling of the bank length in TCHXLD. * * TCHXLD and TCVRHX: (M.T. Lakata) 21/11/95 * compiler strictness problems. * None of the changes should have any real effect. * * TXTGET: (M.T. Lakata) 21/11/95 * Added PWC/SVX loading for vertex fitting, * which got lost in correction cradle for * locater_20003. * * TCTRAK: (C.A.Meyer) 17/11/95 * In order to retrack the vertex information without * retracking the helix information, (something one * needs to do in repeating global tracking), it is * necessary to "reset" word 3 of all the helix * banks to zero. I have implemented this change here. * * TSZERO: (Chr. Voelcker) 15.Nov. 1995 * just put an error message in more suggestive words. * * TCFFUZ: (Chr. Voelcker) 15.Nov. 1995 * NART was only defined for the first call, when QONCE is .TRUE. * I SAVEed the variable. * *::> 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 * *::> VERSION 2.00/00 24/08/95 17.44.50 +REP,$KUMACS,INSTALL,1-11. *CMZU: 2.00/04 01/11/95 17.12.44 by Christian Voelcker *CMZ : 1.99/23 23/08/95 18.48.02 by Christian Voelcker *CMZU: 1.99/19 28/07/95 10.43.52 by Curtis A. Meyer *CMZU: 1.99/11 06/07/95 23.58.32 by Curtis A. Meyer *CMZU: 1.99/05 21/06/95 16.56.41 by Christian Voelcker *CMZU: 1.97/02 08/05/95 16.40.45 by Christian Voelcker *CMZU: 1.96/03 10/04/95 15.53.43 by Curtis A. Meyer *CMZU: 1.95/06 16/03/95 09.41.00 by Curtis A. Meyer *CMZU: 1.94/07 14/02/95 16.11.33 by Curtis A. Meyer *CMZ : 1.94/04 10/02/95 17.12.24 by Curtis A. Meyer *CMZ : 1.94/02 02/02/95 08.47.08 by Curtis A. Meyer *-- Author : Curtis A. Meyer 02.Nov.1994 +REP,$KUMACS,INSTALL,26-31. * 01/11/95 Christian Voelcker * select $OS for VMS systems * Macro INSTALL src=X * * [src] = src directory * [2-4] = options * * options: CALIB = Create Calibration Version +REP,$KUMACS,INSTALL,45-72. if [src]='X' Then Exec help Goto END Endif If [src]='?' Then Exec help Goto END Endif * * CODEDIR Source code directory * CODEDIR = $Lower([src]) message Putting source code files in [CODEDIR] * * set defaults * CCODEFLAG = TRUE DEBUGFLAG = FALSE TEXTFLAG = FALSE LIBFLAG = TRUE GNUFLAG = FALSE Select . * * check options * FAIL = FALSE Do I = 2,[#] +REP,$KUMACS,INSTALL,115-132. Enddo If [FAIL] = TRUE then exec help goto end Endif If [CCODEFLAG] = TRUE then Select CCODE Endif Case $OS In (UNIX) File $CBCMZ/cboff.cmz -R (VMS) File CB$CMZ:cboff.cmz -R Endcase Select $MACHINE Select $OS Case $Machine In +REP,$KUMACS,INSTALL,141-150. Endcase * Mess ' ' Mess Installation of LOCATER on $MACHINE $OS starting at $DATE $TIME . Select Mess ' ' * * Set up the FORTRAN compiling extensions * Case $MACHINE In +REP,$KUMACS,INSTALL,162-168. Endcase * Set FORTRAN -LAN * * Set up filenames for the extraction of the fortran code. * Case $OS In +REP,$KUMACS,INSTALL,175-179. Endcase * * Set up the debug options * If [DEBUGFLAG] = TRUE Then +REP,$KUMACS,INSTALL,194. Else +REP,$KUMACS,INSTALL,201-220. Endif * * choose fortran compiler options and other stuff * If [LIBFLAG] = TRUE Then Case $MACHINE In (HPUX) Set 'f77 -c'//[DEBUGOPT]//' +ppu $compfile' -C f77 Set 'cc -c -O -w '//[DEBUGCCC]//' $compfile' -C C (SUN) Set 'f77 -c'//[DEBUGOPT]//' -fnonstd -Nl100 -w $compfile' -C f77 Set '`which gcc` -c'//[DEBUGCCC]//' -w $compfile' -C C (IBMRT) Set 'xlf -c'//[DEBUGOPT]//' -qextname $compfile' -C f77 Set 'cc -c -w '//[DEBUGCCC]//' $compfile' -C C (SGI) Set 'f77 -c'//[DEBUGOPT]//' $compfile G 3' -C f77 Set 'cc -c -w '//[DEBUGCCC]//' $compfile' -C C (DECS) Set 'f77 -c -w -Nl99 '//[DEBUGOPT]//' $compfile' -C f77 Set 'cc -c -w '//[DEBUGCCC]//' $compfile' -C C (VAX) Set 'FOR/OBJ=$compfile.OBJ'//[DEBUGOPT]//' $compfile' -C f77 Set 'CC/OBJ=$compfile.obj'//[DEBUGCCC]//' $compfile' -C C * (ALPHA) Case $OS In +REP,$KUMACS,INSTALL,229-236. IF [GNUFLAG] = TRUE THEN Set '`which g77` -c'//[DEBUGOPT]//' -w $compfile' -C f77 Set '`which gcc` -c'//[DEBUGCCC]//' -w $compfile' -C C Endif * Sequence //cboff/commcb * * make the library... +REP,$KUMACS,INSTALL,267-271. Endif * END: * Return +REP,TCCOMMON,TC_MACRO,1-16. *CMZU: 2.00/03 24/10/95 18.23.17 by M.T. Lakata *CMZU: 1.99/22 08/08/95 19.37.14 by Christian Voelcker *CMZU: 1.99/21 02/08/95 15.57.59 by Christian Voelcker *CMZU: 1.99/08 22/06/95 18.05.16 by Christian Voelcker *CMZU: 1.99/00 09/06/95 21.37.00 by Christian Voelcker *CMZU: 1.98/00 08/06/95 18.52.58 by Christian Voelcker *CMZ : 1.97/02 08/05/95 15.47.03 by Christian Voelcker *CMZU: 1.96/06 25/04/95 10.46.23 by Curtis A. Meyer *CMZ : 1.96/03 10/04/95 15.20.02 by Curtis A. Meyer *CMZU: 1.96/00 06/04/95 17.05.10 by Christian Voelcker *CMZU: 1.95/04 15/03/95 17.47.50 by Curtis A. Meyer *-- Author : Curtis A. Meyer * *-- Modified: 24-May-1994, C.A.Meyer * Add OLDJDC flag to TCFLAG, and OLDJLG to LGHOLD. * *-- Modified: 29-AUG-1993, R.Bossingham +REP,TCCOMMON,TC_MACRO,562. * Here starts the statistic of the fuzzy Radon pattern recognition +REP,TCCOMMON,TC_MACRO,750. * PWC cluster information (added by M.Doser) +REP,TCCOMMON,TC_MACRO,867-882. * Common block TCVTCM, for vertex fit parameters * VERTVR - TCVER? version * VFITVR - 1=TCVRTX 2=TCVRHX * 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 +REP,TC_MAIN,CJINIT,1-33. *CMZU: 2.00/05 12/12/95 17.45.41 by Christian Voelcker *CMZU: 2.00/03 24/10/95 17.57.34 by M.T. Lakata *CMZU: 1.99/22 08/08/95 19.17.09 by Christian Voelcker *CMZU: 1.99/00 09/06/95 21.58.30 by Christian Voelcker *CMZU: 1.98/00 09/06/95 19.09.30 by Christian Voelcker *CMZ : 1.97/03 08/05/95 16.54.38 by Christian Voelcker *CMZU: 1.96/06 25/04/95 10.45.32 by Curtis A. Meyer *CMZ : 1.96/03 10/04/95 15.21.06 by Curtis A. Meyer *CMZU: 1.95/03 15/03/95 16.07.20 by Curtis A. Meyer *CMZ : 1.95/02 13/03/95 18.21.51 by Curtis A. Meyer *CMZU: 1.93/01 18/07/94 13.26.56 by Curtis A. Meyer *CMZU: 1.92/04 02/07/94 15.06.20 by Curtis A. Meyer *CMZ : 1.92/03 02/07/94 14.53.06 by Curtis A. Meyer *CMZU: 1.91/03 13/06/94 10.43.03 by Curtis A. Meyer *CMZ : 1.90/00 25/05/94 16.52.16 by Curtis A. Meyer *CMZU: 1.61/00 03/03/94 08.17.57 by Curtis A. Meyer *CMZ : 1.58/00 16/01/94 13.36.18 by Curtis A. Meyer *CMZ : 1.53/00 20/09/93 11.32.26 by Curtis A. Meyer *CMZU: 1.51/01 06/04/92 13.49.13 by Curtis A. Meyer *CMZ : 1.51/00 13/01/92 14.30.09 by Curtis A. Meyer *CMZ : 1.50/19 13/12/91 13.16.39 by Curtis A. Meyer *CMZ : 1.49/01 16/08/91 13.46.01 by Curtis A. Meyer *CMZ : 1.49/00 06/08/91 13.30.49 by Curtis A. Meyer *CMZ : 1.47/00 25/06/91 16.12.58 by Curtis A. Meyer *CMZU: 1.46/03 20/03/91 17.11.55 by Curtis A. Meyer *CMZ : 1.45/05 27/02/91 11.32.08 by Curtis A. Meyer *CMZ : 1.45/02 05/02/91 18.18.07 by Curtis A. Meyer *CMZ : 1.45/01 03/02/91 12.53.21 by Curtis A. Meyer *CMZ : 1.45/00 15/01/91 13.27.52 by Curtis A. Meyer *CMZ : 1.44/09 21/11/90 11.46.20 by Curtis A. Meyer *CMZ : 1.44/02 02/10/90 16.15.21 by CURTIS A. MEYER *CMZ : 1.43/02 29/08/90 18.51.52 by Curtis A. Meyer *CMZ : 1.43/01 14/08/90 18.08.10 by Curtis A. Meyer *CMZ : 1.42/00 11/07/90 08.27.22 by Curtis A. Meyer *-- Author : +ADD,TC_MAIN,CJINIT,87. * 12 December 1995 C. Meyer/ Ch. Voelcker * tiny bug initializing DLZZMC * +REP,TC_MAIN,CJINIT,128. * hollerith strings for PREC card +REP,TC_MAIN,CJINIT,244-245. *---Under UNIX, check to see if a file is assigned to the unit. If *---so, then get the file name. +REP,TC_MAIN,CJINIT,253-276. *---On a SUN, the environment is FORTii, on the ALLIANT it is *---FOR0ii. * &SELF,IF=ALT. 198 FORMAT('FOR00',I1) 199 FORMAT('FOR0',I2) &SELF,IF=DECS,HPUX,IBMRT,NXT,SGI,SUN,UNIX. 198 FORMAT('FORT0',I2) 199 FORMAT('FORT',I2) CLOSE(UNIT=LJCAL,ERR=200) &SELF,IF=ALT,DECS,HPUX,IBMRT,NXT,SGI,SUN,UNIX. * *---Get the environment variable for this unit. * 200 CALL GETENV(IENVIR,FILENM) * *---Loop over the returned name, and find the file. * ILENF = LENOCC(FILENM) IF(ILENF .LE. 0) RETURN * *---Open the file. * 220 OPEN(UNIT=LJCAL,FILE=FILENM(1:ILENF),STATUS='OLD', +REP,TC_MAIN,CJINIT,284-286. *---Check to see if the file exists. Under UNIX, this is redundent. * 250 INQUIRE(LJCAL,EXIST=LEXIST) +REP,TC_MAIN,CJINIT,299. 100 CONTINUE +REP,TC_MAIN,CJINIT,326-334. * set nominal values for TCVTCM common * * default for VERTVR is old, 1 vertex fit VERTVR = 1 * default of fitter is TCVRTX (1=tcvrtx, 2=tcvrhx) VFITVR = 1 V3CONF = V3COND V3DIST = V3DISD V3BACK = V3BACD V3MNHT = V3MNHD VHCFL = VHCFLD VHERC = VHERCD VHGCH = VHGCHD VHMCH = VHMCHD VHCT = VHCTD +REP,TC_MAIN,CJINIT,469. * user control of FAST FUZZY pattern recognition parameters +REP,TC_MAIN,CJINIT,480-484. * user control of which Pattern RECognition shall be used * CALL FFKEY('PREC',PRECTC,1,'MIXED') * * user control of svtx parameters +REP,TC_MAIN,CJINIT,493-499. * Vertex fitting parameters * CALL FFKEY('VERT',VERTVR, 1,'I') ! which TCVER? to use CALL FFKEY('VFIT',VFITVR, 1,'I') ! which TCVRTX/TCVRHX to use * 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,CJINIT,506-511. * check Vertex values * IF (VERTVR.LT.1.OR.VERTVR.GT.3) THEN WRITE(LLOG,294 ) VERTVR WRITE(LERR,294 ) VERTVR 294 FORMAT(' Locater Card VERT ',I,' is out range [1,3]') +REP,TC_MAIN,CJINIT,517-521. 295 FORMAT(' Locater Card VFIT ',I,' is out range [1,2]') STOPCB = .TRUE. ENDIF * * decode PRECTC +REP,TC_MAIN,CJINIT,542-547. 293 FORMAT(' WARNING: Smearing of Monte Carlo JDC', & /,' Data has been USER Disabled!!!!',/, & /,' WARNING: Results will be too good!!') * IF ( DLZZ(2) .NE. 0.0) THEN * DLZZMC = DLZZ(2) DLZZMC = DLZZ(1) +REP,TC_MAIN,CJINIT,596. 300 CONTINUE +REP,TC_MAIN,CJINIT,606. 315 CONTINUE +REP,TC_MAIN,CJINIT,612. 325 CONTINUE +REP,TC_MAIN,CJINIT,649. 500 CONTINUE +REP,TC_MAIN,CJINIT,655. 600 CONTINUE +REP,TC_MAIN,CJINIT,706-707. 9998 CLOSE(UNIT=LJCAL,ERR=9999) 9999 RETURN +REP,TC_MAIN,TCTRAK,1-28. *CMZU: 2.00/05 17/11/95 16.06.23 by Curtis A. Meyer *CMZU: 1.99/22 08/08/95 19.14.57 by Christian Voelcker *CMZU: 1.99/00 10/06/95 00.06.06 by Christian Voelcker *CMZU: 1.98/00 09/06/95 15.31.07 by Christian Voelcker *CMZU: 1.96/03 10/04/95 15.07.55 by Curtis A. Meyer *CMZU: 1.96/00 05/04/95 17.12.20 by Christian Voelcker *CMZU: 1.95/02 14/03/95 11.47.42 by Curtis A. Meyer *CMZU: 1.94/03 09/02/95 18.40.43 by Christian Voelcker *CMZ : 1.94/02 02/02/95 08.12.21 by Curtis A. Meyer *CMZU: 1.94/00 19/01/95 10.38.31 by Christian Voelcker *CMZU: 1.93/04 06/09/94 11.20.33 by Curtis A. Meyer *CMZU: 1.92/05 05/07/94 20.27.54 by Curtis A. Meyer *CMZ : 1.90/01 26/05/94 10.45.33 by Curtis A. Meyer *CMZ : 1.62/04 19/05/94 10.39.48 by Curtis A. Meyer *CMZU: 1.51/06 01/06/92 11.38.31 by Curtis A. Meyer *CMZU: 1.50/12 06/11/91 09.47.22 by Curtis A. Meyer *CMZ : 1.48/00 31/07/91 11.36.23 by Curtis A. Meyer *CMZ : 1.45/06 06/03/91 16.52.16 by Curtis A. Meyer *CMZ : 1.45/00 15/01/91 13.27.53 by Curtis A. Meyer *CMZ : 1.44/10 16/12/90 18.30.27 by Curtis A. Meyer *CMZ : 1.44/09 31/10/90 11.26.57 by Curtis A. Meyer *CMZ : 1.44/06 27/10/90 13.09.29 by Curtis A. Meyer *CMZ : 1.44/04 23/10/90 08.55.22 by Curtis A. Meyer *CMZ : 1.44/03 10/10/90 11.57.53 by Curtis A. Meyer *CMZ : 1.44/00 24/09/90 08.32.55 by Curtis A. Meyer *CMZ : 1.43/02 31/08/90 18.06.08 by Curtis A. Meyer *CMZ : 1.43/01 15/08/90 18.51.33 by Curtis A. Meyer *CMZ : 1.42/00 11/07/90 08.27.20 by Curtis A. Meyer *-- Author : Curtis A. Meyer +REP,TC_MAIN,TCTRAK,41-42. * Modified: 17 November, 1995 * * References: The locater offline manual, and the locater * USER's manual. +ADD,TC_MAIN,TCTRAK,146. * 17/11/95 C.A.Meyer * In order to retrack the vertex information without * retracking the helix information, (something one * needs to do in repeating global tracking), it is * necessary to "reset" word 3 of all the helix * banks to zero. I have implemented this change * here. * +REP,TC_MAIN,TCTRAK,207. * +SEQ,TCVERT. +REP,TC_MAIN,TCTRAK,225-230. 9191 FORMAT(' Called for event ',I9) &SELF. &SELF,IF=JDCALIBR. DO 100 I = 1,NESTCJ IESTCJ(I) = 0 100 CONTINUE +REP,TC_MAIN,TCTRAK,445. 200 CONTINUE +REP,TC_MAIN,TCTRAK,495. 300 CONTINUE +REP,TC_MAIN,TCTRAK,576. 310 CONTINUE +REP,TC_MAIN,TCTRAK,591. 350 CONTINUE +ADD,TC_MAIN,TCTRAK,755. * If we did NOT redo the helix fit, but want to redo the * vertex fit, we need to loop over the helix banks, and * set word 3 to zero for each track. * IF ( .NOT. HELXTC ) THEN * DO 375 I = 1 , IQ(LTCTR+1) IQ(LQ(LTCTR-I)+3) = 0 375 CONTINUE * ENDIF * +REP,TC_MAIN,TCTRAK,792-796. * Call the desired vertex fitter. * * TCVERT = default one vertex fit * TCVER2 = finds all possible vertices, for CBKFIT * TCVER3 = finds best configuration of multiple vertices +REP,TC_MAIN,TCTRAK,838. 400 CONTINUE +REP,TC_MAIN,ZBFORM,0-11. &DECK,ZBFORM,IF=NEVER. *CMZU: 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) *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 : +REP,TC_MAIN,ZBFORM,130. 100 CONTINUE +ADD,TC_MAIN,TCFORM,*. &DECK,TCFORM. *CMZU: 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 +REP,TC_RAWS,TJDCGT,1-28. *CMZU: 2.00/04 01/11/95 14.50.28 by Christian Voelcker *CMZU: 1.95/05 15/03/95 18.30.10 by Curtis A. Meyer *CMZ : 1.95/02 14/03/95 11.54.58 by Curtis A. Meyer *CMZU: 1.92/00 14/06/94 16.00.40 by Curtis A. Meyer *CMZ : 1.91/03 13/06/94 09.49.56 by Curtis A. Meyer *CMZU: 1.52/02 24/07/92 17.23.12 by Michael Doser *CMZ : 1.52/01 20/07/92 13.38.33 by Michael Doser *CMZU: 1.51/07 04/06/92 13.45.19 by Curtis A. Meyer *CMZ : 1.51/06 01/06/92 13.33.55 by Curtis A. Meyer *CMZU: 1.51/05 05/05/92 09.25.30 by Michael Doser *CMZU: 1.51/04 10/04/92 14.07.35 by Curtis A. Meyer *CMZ : 1.51/03 09/04/92 17.31.38 by Curtis A. Meyer *CMZ : 1.51/02 08/04/92 08.07.48 by Curtis A. Meyer *CMZU: 1.50/17 13/11/91 15.57.44 by Michael Doser *CMZ : 13/11/91 15.56.42 by Michael Doser * introduce an error check (actually for TCSGMT) to reject * events with more than 400 hits *CMZU: 1.50/02 20/08/91 12.10.05 by Gunter Folger *CMZ : 1.49/00 06/08/91 15.45.22 by Curtis A. Meyer *CMZ : 1.46/04 30/04/91 15.08.49 by Curtis A. Meyer *CMZU: 1.46/01 15/03/91 12.37.53 by Curtis A. Meyer *CMZ : 1.46/00 07/03/91 11.56.17 by Curtis A. Meyer *CMZ : 1.45/06 06/03/91 16.50.30 by Curtis A. Meyer *CMZ : 1.45/00 15/01/91 13.33.22 by Curtis A. Meyer *CMZ : 1.44/09 11/11/90 14.04.06 by Curtis A. Meyer *CMZ : 1.44/00 24/09/90 08.32.55 by Curtis A. Meyer *CMZ : 1.43/05 21/09/90 07.49.58 by Curtis A. Meyer *CMZ : 1.43/02 30/08/90 10.14.17 by Curtis A. Meyer *-- Author : Curtis A. Meyer +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,313. 100 CONTINUE +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) +REP,TC_RAWS,TJDCGT,357. 200 CONTINUE +REP,TC_RAWS,TJDCGT,448-450. 500 IRJDC = IRJDC + ILEN * 600 CONTINUE +REP,TC_RAWS,TJDCGT,545-550. 1000 CONTINUE * * Loop over all the hits, starting at the outermost ones. * I = LAYER(LYR,ISEC) 1100 I = I - 1 +REP,TC_RAWS,TJDCGT,605. 1250 CONTINUE +REP,TC_RAWS,TJDCGT,620. 1300 DO 1400 J = 1,LAYER(LYR,ISEC) +REP,TC_RAWS,TJDCGT,663-665. 1400 CONTINUE ENDIF 1500 CONTINUE +REP,TC_RAWS,TJDCGT,671. 2000 CONTINUE +REP,TC_RAWS,TSMUSH,1-3. *CMZU: 2.00/04 06/09/95 07.53.23 by Curtis A. Meyer *CMZU: 1.95/03 15/03/95 16.06.41 by Curtis A. Meyer *CMZ : 1.95/02 15/03/95 08.53.37 by Curtis A. Meyer *-- Author : Curtis A. Meyer 14/03/95 +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,76-84. *---Make sure that these are Monte Carlo Data. * IF(.NOT.LTYPMC) RETURN * *---Make sure that the TJDC bank exists. * IF(ITJDC .LE. 0) RETURN * *---Create a local copy of the pointer. +REP,TC_RAWS,TSMUSH,112-127. 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) * *---Smear z---In 76% of the cases, we want to smear with a *---gaussian of width DELZMC centered at -0.01. In 24% *---we will smear with a Gaussian of width 1.72 centered *---at 0.13. Tune the 1.72 to 1.42 to better agree with *---data. Change 76% to 74.5% * CALL GRNDM(P,1) IF (P.LE.0.745) THEN 850 DELZ=RNDNOR(-0.01,DELZMC) +REP,TC_RAWS,TSMUSH,133-134. *---It is now necessary to recompute the amplitudes based on *---this new z-position. +REP,TC_RAWS,TSRPHI,1-3. *CMZU: 2.00/04 13/10/95 08.30.43 by Curtis A. Meyer *CMZU: 1.99/11 07/07/95 12.38.25 by Curtis A. Meyer *CMZU: 1.95/03 15/03/95 13.37.00 by Curtis A. Meyer *CMZ : 1.95/02 14/03/95 11.42.22 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 * +REP,TC_RAWS,TSRPHI,167. * emulate Rannor (taken from GEANT 3.13 source) +REP,TC_RAWS,TSRPHI,209-215. 100 CONTINUE ISEED(1,ISEQ) = ISEED1 ISEED(2,ISEQ) = ISEED2 END CDECK ID>, GRNDMQ. *CMZ : 3.15/01 08/05/91 11.49.15 by Federico Carminati *-- Author : Rene Brun 23/02/89 +REP,TC_RAWS,TSRPHI,484. 999 RETURN +REP,TC_FFUZ,TCFFUZ,1-5. *CMZU: 2.00/05 15/11/95 06.59.55 by Christian Voelcker *CMZU: 1.99/18 27/07/95 14.52.00 by Christian Voelcker *CMZU: 1.96/06 25/04/95 10.38.09 by Curtis A. Meyer *CMZ : 1.96/05 18/04/95 10.55.50 by Curtis A. Meyer *CMZ : 1.96/03 10/04/95 16.11.14 by Curtis A. Meyer *-- Author : Bernd Kaemmle +ADD,TC_FFUZ,TCFFUZ,33. * 15.Nov. 1995 Chr. Voelcker * NART was only defined for the first call, when QONCE is .TRUE. * I SAVEed the variable. * +REP,TC_FFUZ,TCFFUZ,58. SAVE IDUMMY,NART +REP,TC_FFUZ,TCFFUZ,86. 121 FORMAT(' **************************************************** ',/, +REP,TC_FFUZ,TCFFUZ,166. 1000 CONTINUE +REP,TC_FFUZ,TCFFUZ,209. ** changed by Rouared CALL TCFSEL(QFIRST,NBLOCK,IBLOCK) +REP,TC_FFUZ,TCFFUZ,224. 1500 CONTINUE +REP,TC_FFUZ,TCFFUZ,246-252. *** IF(QFAST_FIT .AND. .NOT.QGOOD) THEN * * If the fast fuzzy Radon fit did not succeed try a long one. * ***** QFAST_FIT = .FALSE. ***** GO TO 1500 ***** END IF +REP,TC_FFUZ,TCFFUZ,258. 2000 CONTINUE +REP,TC_FFUZ,TCFFUZ,295. 999 RETURN +REP,TC_FFUZ,TCMATR,1-6. *CMZU: 2.00/01 05/09/95 14.42.20 by Christian Voelcker *CMZU: 1.99/21 01/08/95 15.16.00 by F.-H.Heinsius *CMZU: 1.96/05 18/04/95 10.49.46 by Curtis A. Meyer *CMZ : 1.96/03 10/04/95 15.32.57 by Curtis A. Meyer *-- Author : Maurice Benayoun (BNY) 12/02/95 SUBROUTINE TCMATR(INDEX,KHIT,XX,YY,XX0,YY0,RR0) C----------------------------------------------------------------- +REP,TC_FFUZ,TCMATR,28. C C Modifications: 05/09/95 M. Benayoun/ Chr. Voelcker C remove printing. C----------------------------------------------------------------- +REP,TC_FFUZ,TCMATR,43-63. cccccc IR(3) ! for matrix inversion (absolutely needed) REAL XX0,YY0,RR0 LOGICAL LECR C--- for checks -------------------------- c INTEGER ICOUNT1,ICOUNT2, ICOUNT3 c DATA ICOUNT1,ICOUNT2,ICOUNT3/0,0,0/ SAVE C----------------------------------------------- LECR=KY25TC.EQ.1 KEVT=IEHDCB(5) CCCC**** LECR=KEVT.EQ.199 ! get printout for evt 199 C--- for checks -------------------------- NP=KHIT CFHH CMN=0. CFHH CSG=0. +REP,TC_FFUZ,TCMATR,90-94. 1000 CONTINUE ! Entry point from IFLAG=-1 (see below) * !NP>=3 X3=XX(LL) Y3=YY(LL) 1200 CONTINUE ! entry point for NP=2 +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 +REP,TC_FFUZ,TCMATR,138-142. C-------------------------------------------------------------- C-- here is general case np>3 C-------------------------------------------------------------- C ICOUNT2=ICOUNT2+1 +REP,TC_FFUZ,TCMATR,169-170. C S0,SX,SY,SX2,SY2,SXY,SX3,SY3,SXY2,SX2Y C loading VEC and AMAT +REP,TC_FFUZ,TCMATR,183-203. C--> EQUATION IS AMAT*PARAM=VEC C WITH PARAM(1)=2*X0, PARAM(2)=2*Y0 AND PARAM(3)=R^2-X0^2-Y0^2 C------------------------------------------------------------------ C SOLVE THE EQUATION C C------------------------------------------------------------------ CALL DEQINV(3,AMAT,3,IR,IFAIL,1,VEC) C ICOUNT3=ICOUNT3+1 +REP,TC_FFUZ,TCMATR,230-236. C------------------------------------------------------------------ C--- CF F010 on return AMAT is replaced by AMAT**-1 C VEC=(2*x0,2*y0,c) C NORMAL RETURN CODE IS IFAIL=0 C IF IFAIL=-1 AMAT is singular then AMAT IS CHANGED C and VEC IS UNCHANGED C------------------------------------------------------------------ +REP,TC_FFUZ,TCMATS,1-9. *CMZU: 2.00/01 05/09/95 14.44.28 by Christian Voelcker *CMZU: 1.99/21 01/08/95 16.48.21 by F.-H.Heinsius *CMZU: 1.96/05 18/04/95 10.49.46 by Curtis A. Meyer *CMZ : 1.96/03 10/04/95 15.33.09 by Curtis A. Meyer *-- Author : Maurice Benayoun (BNY) 12/02/95 SUBROUTINE TCMATS(INDEX,KHIT,XX,YY,DELR2,XX0,YY0,RR0) C----------------------------------------------------------------- C Author: Maurice Benayoun C C Creation date : february 1995 C C Modifications : 05/09/95 M. Benayoun/ Chr. Voelcker C remove pronting +REP,TC_FFUZ,TCMATS,28. C----------------------------------------------------------------- +REP,TC_FFUZ,TCMATS,43-63. cccccc IR(3) ! for matrix inversion (absolutely needed) REAL XX0,YY0,RR0 LOGICAL LECR C--- for checks -------------------------- c INTEGER ICOUNT1,ICOUNT2, ICOUNT3 c DATA ICOUNT1,ICOUNT2,ICOUNT3/0,0,0/ SAVE C----------------------------------------------- LECR=KY25TC.EQ.1 KEVT=IEHDCB(5) CCCC**** LECR=KEVT.EQ.199 ! get printout for evt 199 C--- for checks -------------------------- NP=KHIT CFHH CMN=0. CFHH CSG=0. +REP,TC_FFUZ,TCMATS,90-94. 1000 CONTINUE ! Entry point from IFLAG=-1 (see below) * !NP>=3 X3=XX(LL) Y3=YY(LL) 1200 CONTINUE ! entry point for NP=2 +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 +REP,TC_FFUZ,TCMATS,138-142. C-------------------------------------------------------------- C-- here is general case np>3 C-------------------------------------------------------------- CFHH ICOUNT2=ICOUNT2+1 +REP,TC_FFUZ,TCMATS,170-171. C S0,SX,SY,SX2,SY2,SXY,SX3,SY3,SXY2,SX2Y C loading VEC and AMAT +REP,TC_FFUZ,TCMATS,184-204. C--> EQUATION IS AMAT*PARAM=VEC C WITH PARAM(1)=2*X0, PARAM(2)=2*Y0 AND PARAM(3)=R^2-X0^2-Y0^2 C------------------------------------------------------------------ C SOLVE THE EQUATION C C------------------------------------------------------------------ CALL DEQINV(3,AMAT,3,IR,IFAIL,1,VEC) C ICOUNT3=ICOUNT3+1 +REP,TC_FFUZ,TCMATS,231-237. C------------------------------------------------------------------ C--- CF F010 on return AMAT is replaced by AMAT**-1 C VEC=(2*x0,2*y0,c) C NORMAL RETURN CODE IS IFAIL=0 C IF IFAIL=-1 AMAT is singular then AMAT IS CHANGED C and VEC IS UNCHANGED C------------------------------------------------------------------ +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 +REP,TC_RADON2,TCRAFU,1-129. *CMZU: 2.00/04 19/09/95 01.22.17 by Holger Stoeck *CMZU: 1.95/01 07/03/95 21.59.24 by Holger Stoeck *CMZU: 1.94/05 21/02/95 10.40.10 by Holger Stoeck *CMZU: 1.94/03 10/02/95 16.27.06 by Christian Voelcker *CMZU: 1.94/00 08/01/95 19.14.03 by Holger Stoeck *-- Author : Holger Stoeck 27/04/94 /***************************************************************************** * * Author: Holger Stoeck * * Creation date: 27/4/1994 * Modified: * * References: * * Description: This file contains various functions for the * fuzzy radon patt. recog. RADON2. * * Subroutines referenced: * * Edit History: 25/1/1995 H. Stoeck * Fix the bug which produced the compiler warnings. * Changed the binning values for the start * parameters calculation. * Changed also the order of looping. * * 10/02/1995 Chr. Voelcker (on request from H.S.) * Set MINUIT output to zero (function minuit()). * * 20/2/1995 H. Stoeck * Changed the 3 dimensional fuzzy Radon transform * to a 2 dimensional one. * Now only the r and phi information is used. * * 7/3/1995 H. Stoeck * Fixed a bug in the calculation of the starting * values for the fit. * Built in some more securities. * *****************************************************************************/ /************************************************ * Include files * ************************************************/ #include &SELF,IF=VMS. #include &SELF. /************************************************ * Some global variables * ************************************************/ int max_hit; &SELF,IF=DECS. int radon_minuit; &SELF. double HIT_X[200], HIT_Y[200], SIN[200], COS[200]; double k_top, p_top, SIGMA[2]; /************************************************ * Function radon_hit_density * ************************************************/ double radon_hit_density(double *TRACK,double *ANGLE,double *X) { double eta_g; eta_g = sqrt((((TRACK[0] * X[0]) + ANGLE[1]) * ((TRACK[0] * X[0]) + ANGLE[1])) + (((TRACK[0] * X[1]) - ANGLE[2]) * ((TRACK[0] * X[1]) - ANGLE[2]))); return exp(-0.5 * (eta_g - 1) * (eta_g - 1) / (TRACK[0] * TRACK[0] * SIGMA[0] * SIGMA[0])); } /************************************************ * Function getStartPara * ************************************************/ void getStartPara(double START_PARA[],double charge) { int index; int k_bin, k_bin_min, k_bin_max, k_bin_step, k_bin_top; int p_bin, p_bin_min, p_bin_max, p_bin_step, p_bin_top; double k, k_min, k_max, k_step, k_bin_max_float; double p, p_min, p_max, p_step, p_bin_max_float; double track_density, track_density_top, ip; double TRACK[2], ANGLE[3], X[2]; /* * Set the looping and binning values. * CHARGE is the charge of the reconstructed track. * Consider only the solution -CHARGE. */ if (charge == 1.0) { k_min = -2.0; k_max = 0.0; } else { k_min = 0.0; k_max = 2.0; } k_step = 0.4; p_min = 0.0; p_max = 360.0; p_step = 10.0; /* * Transfer the float values to integer values. */ k_bin_min = 1; k_bin_max_float = (k_max - k_min) / k_step; k_bin_step = 1; if (modf(k_bin_max_float,&ip) > 0.5) {k_bin_max = ip + 1;} else {k_bin_max = ip;} p_bin_min = 1; p_bin_max_float = (p_max - p_min) / p_step; p_bin_step = 1; if (modf(p_bin_max_float,&ip) > 0.5) {p_bin_max = ip + 1;} else {p_bin_max = ip;} /* * Begin calculation of the start parameters */ track_density_top = 0.0; /* * First Loop over all PHI values. */ for (p_bin = (p_bin_min - 1); p_bin <= p_bin_max; p_bin++) { +REP,TC_RADON2,TCRAFU,136-138. /* * Loop over all KAPPA values. */ +REP,TC_RADON2,TCRAFU,153-155. /* * At last loop over all hits. */ +REP,TC_RADON2,TCRAFU,164-167. /* * If the calculated track density > 0.0 * and > old track density, store the parameter. */ +REP,TC_RADON2,TCRAFU,175-310. } /* * Get the real track parameters from the binning values * and check the found starting values. */ START_PARA[0] = (p_bin_top * p_step) + p_min; START_PARA[1] = (k_bin_top * k_step) + k_min; if (fabs(START_PARA[0]) < 0.1) { START_PARA[0] = 99; START_PARA[1] = 99; } } /************************************************ * Function fcnfit * ************************************************/ int fcnfit_(npar, gin, f, x, iflag) long int *npar, *iflag; double *gin, *f, *x; { /* * This is the function to be fitted by MINUIT. * The MINUIT manual suggests to fit the negative of the logarithm * of the sum of the hit densities. * This is to get a more parabel like form of the fit function. */ int i; double eta_g, hit_density, radon_density; radon_density = 0.0; &SELF,IF=DECS. radon_minuit = -1; &SELF. --x ; --gin; for (i = 0; i < max_hit; i++) { eta_g = sqrt((((x[1] * HIT_X[i]) + sin(x[2])) * ((x[1] * HIT_X[i]) + sin(x[2]))) + (((x[1] * HIT_Y[i]) - cos(x[2])) * ((x[1] * HIT_Y[i]) - cos(x[2])))); hit_density = exp(-0.5 * (eta_g - 1) * (eta_g - 1) / (x[1] * x[1] * SIGMA[0] * SIGMA[0])); radon_density += hit_density; } &SELF,IF=DECS. /* * This is to prevent a coredump by MINUIT on DECStations. */ if (radon_density < 0.001) { radon_minuit = 1; radon_density = 1; *f = -1; } else { *f = -log(radon_density); k_top = x[1] ; p_top = x[2] ; } &SELF,IF=UNIX. *f = -log(radon_density); &SELF. return 0; } /************************************************ * Function minuit * ************************************************/ void minuit(double *param,double *steps) { static char strategy[256],pnam[256],chcom[256],output[256]; int nprm,narg,ierflg,futil; double zero,vstrt,stp,arglis; /* * Set the output of MINUIT to the lowest level. */ arglis = -1 ; narg = 1 ; futil = 0 ; &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. /* * Set the strategy. */ arglis = 2 ; narg = 1 ; futil = 0 ; &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. /* * Set the parameters and steps. * At STRATEGY = 2 the steps are reduced by MINUIT by a factor of 1000. * To compensate that the steps are multiplied by a factor of 1000. */ zero = 0 ; nprm = 1 ; vstrt = param[0] ; stp = 1000.0 * steps[0]; &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. /* * Begin fitting using MIGRAD. */ arglis = 0 ; narg = 0 ; futil = 0 ; &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. +REP,TC_CIRC,TXTGET,1-10. *CMZU: 2.00/05 21/11/95 17.41.12 by M.T. Lakata *CMZU: 1.99/11 06/07/95 23.25.56 by Curtis A. Meyer *CMZU: 1.58/00 21/01/94 11.51.37 by F.-H.Heinsius *CMZ : 1.55/00 23/09/93 11.21.39 by R.BOSSINGHAM *-- Author : R.Bossingham *-- Derives in part from Mark Burchell's TCMPB routine. * SUBROUTINE TXTGET(ICODE,ITRK,NPTS,R1,X1,Y1,Z1,D1,D2,D3, & ILAY,ISC,IRS,TIM,AMPL,AMPR,SR,SP) * * Version 0.10 +REP,TC_CIRC,TXTGET,28-30. * 2 or 3 hits taken from TCHT banks. * 4 hits taken from TCTR banks. * +12 or 13 hits taken from TCHT and PWC/SVX banks +REP,TC_CIRC,TXTGET,43-45. * 17 Oct 1995, M. Lakata * Added PWC/SVX loading for vertex fitting * * ********************************************************************* * &SEQ,IMPNONE. &SEQ,CBLINK. &SEQ,CBUNIT. &SEQ,PI2PI. &SEQ,TCANGL. &SEQ,TCPRMS. &SEQ,TRKPRM. +REP,TC_CIRC,TXTGET,54-66. * R1 : R coordinates * X1 : x coordinates * Y1 : y coordinates * Z1 : z coordinates * D1 : x coordinate 1-sigma errors * D2 : y coordinate 1-sigma errors * D3 : z coordinate 1-sigma errors * ILAY : Hit layer * TIM : Drift time * ISC : Sector * IRS : Left/right resolution of hit * SR : Radial error * SP : Phi error +REP,TC_CIRC,TXTGET,74-97. * ITCHT : Location in the TCHT bank of the point being copied * ITCTK : TCTK bank pointer * ITCTR : Location in the TCTR bank of the track copied * ITCHX : Location of the hits bank for the helix fit * INTEGER ITCHT,ITCTK,ITCTR,ITCHX * * I : Loop index * LYR1 : Layer containing the present point * IHIT : Hit number of the present point * ISECT : Sector number * INTEGER I,LYR1,IHIT,ISECT * * TEMPORARY VARIABLES IL1,IS1,ITJ * INTEGER IL1,IS1,ITJ * * IPWC : 10 if pwc hits are to be loaded * IJDC : 2,3,4 * ITPWC : pointer to pwc bank * J : loop index * ITVXT : pointer to svx bank INTEGER IPWC,IJDC,ITPWC,J,ITVXT +REP,TC_CIRC,TXTGET,103-107. * Identify the operation code: * IPWC = ICODE/10 IJDC = MOD(ICODE,10) IF (IJDC.EQ.2 .OR. IJDC.EQ.3) THEN * * Code 2 gets data out of the TCHT data banks. +REP,TC_CIRC,TXTGET,121-127. 800 FORMAT(' TXTGET: ITRK, IQ(LTCTK+1) ',2I20) RETURN ENDIF ********************* * this code was stolen from tchxld * load any PWC hits IF(LTPWC.LE.0.OR.IPWC.EQ.0) GOTO 450 * * Loop over the two layers in the PWC. Down links 1 and 2 are * the hits, whereas down links 3 and 4 are the clusters. We want * to use the cluster information. * DO 400 I = 3,4 * * Make sure that this pwc has data: * IF((IQ(LTPWC+I).LE.0).OR.(LQ(LTPWC-I).LE.0)) GOTO 400 ITPWC = LQ(LTPWC-I) * * Loop over the hits in this layer until we find one * connected to track ITRK: * DO 300 J = 1,IQ(LTPWC+I) IF(IQ(ITPWC+1).EQ.ITRK) THEN * * Add this hit to the track: * NPTS = NPTS + 1 X1(NPTS) = Q(ITPWC+5) Y1(NPTS) = Q(ITPWC+6) Z1(NPTS) = Q(ITPWC+7) * * Store the errors squared in (x,y,z). Note that * the error is z is +/- 25cm 1 sigma! * D1(NPTS) = Q(ITPWC+10)**2 D2(NPTS) = Q(ITPWC+11)**2 D3(NPTS) = 25.0 * * As we can only have one PWC hit per chamber per * track, go to the next pwc. * GOTO 400 * ENDIF * ITPWC = ITPWC + LENCL * 300 CONTINUE 400 CONTINUE 450 CONTINUE * Add any SVX hits IF(LTVXT.LE.0.OR.IPWC.EQ.0) GOTO 4500 IF((IQ(LTVXT+2).LE.0).OR.(LQ(LTVXT-2).LE.0)) GOTO 4500 ITVXT = LQ(LTVXT-2) * * Loop over the clusters in this layer until we find one * connected to track ITRK: * DO 3000 J = 1,IQ(LTVXT+2) * IF(IQ(ITVXT+1).EQ.ITRK) THEN * * Add this hit to the track: * X1(NPTS) = Q(ITVXT+4) Y1(NPTS) = Q(ITVXT+5) Z1(NPTS) = Q(ITVXT+6) * * Store the errors squared in (x,y,z). Note that * the error is z is +/- 3.5cm 1 sigma! * D1(NPTS) = Q(ITVXT+10)**2 D2(NPTS) = Q(ITVXT+11)**2 D3(NPTS) = 3.5 ENDIF * ITVXT = ITVXT + LENCLS * 3000 CONTINUE * 4500 CONTINUE * Now add JDC hits ******************************** * ITCTK = LQ(LTCTK-ITRK) * * Extract the layer and hit number of first hit in track from TCTK bank +REP,TC_CIRC,TXTGET,133-135. * Now compute where this hit is stored in TCHT bank. * 200 CONTINUE +REP,TC_CIRC,TXTGET,141-143. * Now extract the chosen x and y coordinates. * The resolution is specified by the value in the data bank * at IQ(ITCHT+2): (0=Unresolved, 1=Left and 2=Right). +REP,TC_CIRC,TXTGET,173. * Use S.Spanier's error estimate (as used in helix-fit) +REP,TC_CIRC,TXTGET,183-191. * Get the layer and hit number of the next point in this track * from the forward pointer for this track. * 225 CONTINUE LYR1 = IQ(ITCHT+3) IHIT = IQ(ITCHT+4) * * Check to make sure that the pointer is pointing to a hit. * If not then we are done with this track. +REP,TC_CIRC,TXTGET,200-202. ELSEIF(IJDC.EQ.4) THEN * * Make sure that the needed banks exist, and contain data. +REP,TC_CIRC,TXTGET,230. 500 CONTINUE +REP,TC_HELX,TCHXLD,1-16. *CMZU: 2.00/05 28/11/95 16.25.27 by F.-H.Heinsius *CMZU: 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 *CMZ : 1.99/23 23/08/95 18.32.47 by Christian Voelcker *CMZU: 1.99/20 01/08/95 00.04.54 by Christian Voelcker *CMZU: 1.99/19 28/07/95 10.46.25 by Curtis A. Meyer *CMZU: 1.99/00 10/06/95 00.33.32 by Christian Voelcker *CMZU: 1.95/03 15/03/95 12.06.57 by Curtis A. Meyer *CMZU: 1.62/01 28/03/94 23.21.30 by Curtis A. Meyer *CMZ : 1.50/20 09/01/92 14.38.42 by Curtis A. Meyer *CMZU: 1.50/16 22/11/91 13.39.43 by Curtis A. Meyer *CMZ : 1.50/04 18/09/91 12.57.34 by Curtis A. Meyer *CMZ : 1.50/03 05/09/91 16.38.15 by Curtis A. Meyer *CMZ : 1.45/00 15/01/91 13.43.30 by Curtis A. Meyer *CMZ : 1.44/09 07/11/90 16.33.29 by Curtis A. Meyer *CMZ : 1.44/08 29/10/90 16.39.27 by Curtis A. Meyer *CMZ : 1.44/07 27/10/90 14.42.53 by Curtis A. Meyer *CMZ : 1.42/00 11/07/90 08.27.44 by Curtis A. Meyer *-- Author : Curtis A. Meyer +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. * * 27 Nov, 1995 F.-H.Heinsius * Fix bug in previous correction. * +REP,TC_HELX,TCHXLD,159. INTEGER ITCTK,ITCHT,ITCTR,ITCHX,ITPWC,ITVXT +ADD,TC_HELX,TCHXLD,181. * * ERRLOG variables * INTEGER IMAXHIT * +REP,TC_HELX,TCHXLD,228. * MTL: I commented this out IF(NPTS .LE. 5) IERR = 10 +REP,TC_HELX,TCHXLD,284-288. * Now store the charge of the particle, and the * dE/dx information in the new data bank. * Q(ITCTR+7) = Q(ITCTK+8) +REP,TC_HELX,TCHXLD,320. 100 CONTINUE +REP,TC_HELX,TCHXLD,346-351. IPNT = 0 NPWC = 0 NVTX = 0 IF(GPWCTC) THEN * Check to see if there is PWC data. If there is, then * connect it to the helix tracks: +REP,TC_HELX,TCHXLD,416-418. 300 CONTINUE 400 CONTINUE 450 CONTINUE +REP,TC_HELX,TCHXLD,428-435. IF(LTVXT.LE.0) GOTO 4500 * * Make sure that the svtx has data: * IF((IQ(LTVXT+2).LE.0).OR.(LQ(LTVXT-2).LE.0)) GOTO 4500 +ADD,TC_HELX,TCHXLD,446. * Do not add more than two vertex hits. IF (NVTX.GE.2) GOTO 4500 +REP,TC_HELX,TCHXLD,485-500. 3000 CONTINUE 4500 CONTINUE * ENDIF * * * If two PWC/SVX clusters were NOT connected to the track, then we * need to shrink the TCHX bank. * IRES = 0 IF(GPWCTC) THEN IF(LTPWC.GT.0) THEN IF(NPWC.LT.2) THEN IRES = IRES + LENHP*(NPWC-2) ENDIF ENDIF END IF IF(GVTXTC) THEN IF(LTVXT.GT.0) THEN IF(NVTX.LT.2) THEN IRES = IRES + LENHP*(NVTX-2) ENDIF ENDIF ENDIF * IF (IRES.GT.0) THEN ITCHX = LQ(LTCHX-ITRK) CALL MZPUSH(IXTCHX,ITCHX,0,IRES,'R') ITCHX = ITCHX + (NVTX+NPWC)*LENHP +REP,TC_HELX,TCHXLD,507-535. NPTS = NPTS + NVTX + NPWC IF(NPTS.GT.50) THEN CALL ERRLOG(IMAXHIT, & ' maximum hits on track exceeded (50) ') NPTS = 50 ENDIF IQ(ITCTR+8) = NVTX + NPWC IQ(ITCTR+1) = NPTS +REP,TC_HELX,TCHXLD,541-543. * this 500 is for down loop! * 500 IPNT = IPNT + 1 +REP,TC_HELX,TCHXLD,611. 2100 CONTINUE +REP,TC_HELX,TCHXLD,628-629. 2200 CONTINUE 2300 CONTINUE +REP,TC_HELX,TCHXLD,646. 2400 CONTINUE +REP,TC_HELX,TCHXLD,695-700. 2500 CONTINUE * * Store CHISQ in the TCTR bank. * Q(ITCTR+16) = SNGL(SGN) Q(ITCTR+17) = SNGL(CHISQ) +REP,TC_HELX,TCHXLD,743. 2900 CONTINUE +REP,TC_HELX,TCHXLD,765. 4600 RETURN +REP,TC_VERT,TCVER3,1-2. *CMZU: 2.00/03 24/10/95 18.01.14 by M.T. Lakata *CMZU: 1.99/22 11/08/95 14.06.06 by Christian Voelcker *-- Author : +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,156-175. * Error message numbers INTEGER IERRL(10) SAVE IERRL 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,183-200. * If any changes to defaults, print out warning, otherwise * WELCOME to TCVER3 ! 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,217. * first, seperate + from - tracks +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,255-260. * initial guess is average of first two hits ITCHX1 = LQ(LTCHX-TRKS(1)) ITCHX2 = LQ(LTCHX-TRKS(2)) DO 150 K=1,3 XGUES(K)=(Q(ITCHX1+K)+Q(ITCHX2+K))/2.0 150 CONTINUE +REP,TC_VERT,TCVER3,266. * call desired fitter +REP,TC_VERT,TCVER3,275-278. IF(IERR.GT.100) THEN * Check the result and drop the vertex if unsatisfactory * or VTCHSQ=Q(LTCVT+14) +REP,TC_VERT,TCVER3,291-312. 250 CONTINUE 200 CONTINUE * If no vertices, clean up and exit IF (IVRTX.LE.0) GOTO 4000 * now make valid configuration of pairs * this loops over all the vertices. * If it can add the vertex to the configuration without overlapping * with an already present vertex, the current configuration * is pushed onto the stack. If a sufficient number of vertices are in the * current configuration, the configuration is saved. * Once the vertices run out, the stack is * popped, and further configurations are tried. * When the stack is empty, all configurations have been tried * and the valid ones are given in COLIST() * this routine should work for any number of vertices or tracks, DESPAI = NTRAK/2 MINPAI = DESPAI 305 CONTINUE +REP,TC_VERT,TCVER3,321. 300 CONTINUE +REP,TC_VERT,TCVER3,332-338. 307 FORMAT('Too many vertex configs,Increase MAXCNF',I,I) CALL ERRLOG(IERRL(1),CERR) DO 315 I=1,IVRTX TRKSVT(I)=0 315 CONTINUE * Empty the stack 320 CALL TCPOP(IV) +REP,TC_VERT,TCVER3,349-353. 310 CONTINUE ENDIF ENDIF ENDIF 325 CONTINUE +REP,TC_VERT,TCVER3,365-367. 350 CONTINUE * humph, I'm glad that's over :) +REP,TC_VERT,TCVER3,383. 450 CONTINUE +REP,TC_VERT,TCVER3,390-426. 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 * x,y,z is vector difference between vertex positions * px,py,pz is vector difference between vertex momenta X = Q(ITCVT1+5)-Q(ITCVT2+5) Y = Q(ITCVT1+6)-Q(ITCVT2+6) Z = Q(ITCVT1+7)-Q(ITCVT2+7) 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 IF (BESTCH.GT.CHI) THEN BESTCF = I BESTCH = CHI ENDIF 400 CONTINUE +REP,TC_VERT,TCVER3,438-441. 460 CONTINUE ENDIF 500 CONTINUE +REP,TC_VERT,TCVER3,448-452. IF(IBITS(COMASK(BESTCF),I,1).EQ.0.AND. & VALID(I)) THEN NUNUSE = NUNUSE+1 UNUSED(NUNUSE) = I ENDIF 600 CONTINUE +REP,TC_VERT,TCVER3,460-464. * initial guess is average of first two hits ITCHX1 = LQ(LTCHX-TRKS(1)) DO 755 K=1,3 XGUES(K)=Q(ITCHX1+K) 755 CONTINUE +REP,TC_VERT,TCVER3,470. * call desired fitter +REP,TC_VERT,TCVER3,491-505. 700 CONTINUE ENDIF * now try to combine vertices. Try combine each vertex with everyother * vertex. DO 800 IV1=1,IVRTX-1 * don't bother with the vertex if it has been disabled (==0) IF (TRKSVT(IV1).EQ.0) GOTO 800 NVLIST=1 VLIST(1) = IV1 DO 850 IV2=IV1+1,IVRTX * don't bother with the vertex if it has been disabled (==0) IF(TRKSVT(IV2).EQ.0) GOTO 850 * don't bother with the vertices if both are single tracks +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,569-576. 850 CONTINUE IF (NVLIST.GT.1) THEN * aha , we found some close vertices NTR = 0 XGUES(1) = 0.0 XGUES(2) = 0.0 XGUES(3) = 0.0 +REP,TC_VERT,TCVER3,583-593. 920 CONTINUE * initial guess is average of all vertices DO 930 K=1,3 XGUES(K)=XGUES(K)+Q(ITCVT+4+K)/TRKSVT(IV) 930 CONTINUE 900 CONTINUE * Create a vertex and momentum bank for this vertex. IVRTX = IVRTX + 1 CALL MZLIFT(IXTCVX,XXTCVT(IVRTX),0,2,MTCVT,0) * call desired fitter +REP,TC_VERT,TCVER3,613-621. 950 CONTINUE TRKSVT(IVRTX) = NTR ENDIF ENDIF 800 CONTINUE * clear up, finish everything 4000 CONTINUE +REP,TC_VERT,TCVER3,633-637. * assign a vertex # to all the tracks DO 1250 J=1,IQ(ITCVT+1) ITCTR = LQ(LTCTR-IQ(LQ(ITCVT-1)+1+(J-1)*LENVP)) IQ(ITCTR+3) = IGOOD 1250 CONTINUE +REP,TC_VERT,TCVER3,652. 470 CONTINUE +REP,TC_VERT,TCVER3,658. 9999 CONTINUE +REP,TC_VERT,TCVRHX,1-17. *CMZU: 2.00/05 22/11/95 14.13.57 by M.T. Lakata *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(50),X1(50),Y1(50),Z1(50), & D1(50),D2(50),D3(50),TIM(50) REAL AMPL(50),AMPR(50),SR(50),SP(50) INTEGER ILAY(50),ISC(50),IRS(50) * 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./ DATA OLDBET /0.0D0/ * * ______________________________________________________________-- * 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) = DBLE(VHERC*D1(IPT)**2) Y0 (IY,IT) = DBLE(Y1(IPT)) Y (IY,IT) = Y0(IY,IT) GYINV(IY,IT) = DBLE(VHERC*D2(IPT)**2) Y0 (IZ,IT) = DBLE(Z1(IPT)) Y (IZ,IT) = Y0(IZ,IT) GYINV(IZ,IT) = DBLE(VHERC*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 +REP,TC_SVTX,TSZERO,1-8. *CMZU: 2.00/05 15/11/95 06.52.28 by Christian Voelcker *CMZU: 1.99/12 10/07/95 15.00.37 by RAFIK OUARED * new variables OVER(16),VALI(16) to store logical info for * backplane ADCs and new dimensions for TAMP; *CMZU: 1.99/08 22/06/95 17.55.39 by Christian Voelcker *CMZ : 1.99/06 22/06/95 11.01.01 by Christian Voelcker *CMZ : 1.99/01 14/06/95 15.47.05 by RAFIK OUARED *CMZ : 1.99/00 10/06/95 00.49.04 by Christian Voelcker *-- Author : +ADD,TC_SVTX,TSZERO,22. * Modifications: 15/11/95 CV * just put an error message in more suggestive words. * +REP,TC_SVTX,TSZERO,102. * fill up arrays IWIR (storing strip numbers(0-1920)) and IAMP (amplitudes +REP,TC_SVTX,TSZERO,124-126. 200 CONTINUE * * NOW WORK OUT CLUSTERS +REP,TC_SVTX,TSZERO,158. WRITE(LLOG,*) 'TSZERO: number of strips gt 128.. +REP,TC_SVTX,TSZERO,165. 202 CONTINUE +REP,TC_SVTX,TSZERO,185. * COMPUTE TOTAL NUMBER OF SHORT WORDS AND RSVX LENGTH +REP,TC_SVTX,TSZERO,199. * LIFT RSVX BANK +REP,TC_SVTX,TSZERO,210-212. * FILL UP RSVX BANK * * START WITH THE HEADER +REP,TC_SVTX,TSZERO,219. * THEN THE BACKPLANES +REP,TC_SVTX,TSZERO,233. * THEN CLUSTERS..START WITH HEADER AND CONTINUE WITH DATA +REP,TC_SVTX,TSZERO,260. 1000 RETURN