* * $Id$ * * $Log$ * Revision 1.1 1996/03/06 15:31:14 mclareni * Add geane321 history, CMZ and doc files * * *CMZ : 3.21/02 29/03/94 15.41.49 by S.Giani *-- Author : *-- Author : 1. Introduction 2. Definitions 3. Description of the User Routines and COMMONs 4. Examples of application 5. Interface with GEANT 6. Acknowledgements 1. Introduction ============ The present Package allows the user to calculate the average trajectories of particles and to calculate the transport matrix as well as the propagated error covariance matrix. It makes use of a set of routines worked out by the European Muon Collaboration [1] and it is integrated to the GEANT3 system [2] with expected applications in both simulation and reconstruction context. The package is available as a PAM-file. It contains two basic PATCHes, one with the original EMC routines, the other with new GEANT-type tracking routines and interface between them and the EMC routines. This second PATCH contains those routines the user should invoke to carry out the tracking (routines ERTRAK, EUFILL, EUFILP, EUFILV). In addition to this a series of utilities are available for the user (e.g. to transform the track representation from one system to another or to carry out 5 X 5 matrix multiplication in an optimum way). In Section 2 we give the definitions of the mathematical quantities to be dealt with. In Section 3 the description of the user routines are given. In Section 4 we illustrate the application of the program by several examples. Finally in Section 5 the GEANT part is described. Further development will concern the improvement of the error matrix by taking into account the Landau tail in the fluctuation of the energy loss, the bremsstrahlung and direct pair production of the muons. 2. Definitions =========== 2.1 Track variables, Representations -------------------------------- The particle trajectory is characterized by 5 independent variables as a function of one parameter (e.g. the pathlength). Among the 5 variables 1 is related to the curvature (to the absolute value of the momentum, p), 2 are related to the direction of the particle and the other 2 are related to the spatial location. The most usual representation of these 5 parameters are: I. 1/p, lambda, phi, y_perp, z_perp where lambda and phi are the dip and azimuthal angles related to the momentum components in the following way: p_x = p cos(lambda) cos(phi) p_y = p cos(lambda) sin(phi) p_z = p sin(lambda) y_perp and z_perp are the coordinates of the trajectory in a local orthonormal reference frame with the x_perp axis along the particle direction, the y_perp being parallel to the x-y plane. This representation is usually applied in the overall reference frame. (In the EMC code this reference frame is labelled by 'SC' since the overall system was identified with that of the Streamer Chamber.) II. 1/p, y', z', y, z where y'=dy/dx and z'=dz/dx. This representation is particularly useful in fixed target experiments, where the trajectory is evaluated on successive parallel planes (which are perpendicular to the x-axis). (In the EMC code this representation is labelled by 'SP' since a convenient mathematical description of a trajectory being approxima- tely parallel to the x-axis is a 'spline'.) III. 1/p, v', w', v, w where v'=dv/du and w'=dw/du in an orthonormal coordinate system with axis u, v and w. This representation is paricularly useful when the trajectory has to be evaluated on different detector planes in a colliding beam experiment, where the planes can take a great variety of directions.(In the EMC code this representation is labelled by 'SD' as System of Detection.) Of course, all the above representations of the trajectory are equivalent and one can go from one representation to the other by calculating the corresponding Jacobian. These Jacobians are provided by the following EMC routines: S/R TRSCSP from I to II S/R TRSPSC from II to I S/R TRSCSD from I to III S/R TRSDSC from III to I 2.2 Error Propagation ----------------- Let us denote in the following the 5 independent variables at a given value of parameter l_0 (e.g. pathlength) by x_i(l_0), (i=1,...,5). In many applications we are interested in the evolution of the average value of x_i for l>l_0: E(x_i). This is calculated by GEANT as will be outlined in Section 5. The knowledge on the avarage trajectory is characterized by the 5 X 5 covariance matrix of the variables: sigma_$ij(l_0) = E(x_i(l_0).x_j(l_0)) - E(x_i(l_0)).E(x_j(l_0)) We are also interested in the evolution of sigma_$ij for l>l_0, which we call error propagation. If the particle is propagating in a deterministic way, i.e. without any random process involved ( e.g. in vacuum) then the propagation of sigma is simply described by the so called transport matrix in the following way: sigma_$ij(l) = T_$jm(l,l_0).sigma_$mn(l_0)T_$in(l,l_0) where the transport matrix expresses the infinitesimal change in the parameters at l with respect to the change of parameters at l_0: T_$ij(l,l_0) = delta (x_i(l))/ delta (x_j(l_0)). In a realistic detector, however, the particles undergo random processes as well, like Multiple Coulomb scattering, energy loss due to delta ray production, etc. therefore the error propagation should contain an additional term: sigma_$ij(l) = T_$jm(l,l_0).sigma_$mn(l_0).T_$in(l,l_0) + + sigma_$ij:$random(l). The program calculates sigma_$ij(l) step by step using the above recursive formula, where T_$ij and sigma_$ij:$random refers to the actual step and sigma_$mn is the cumulative error for all previous steps. For the mathematical formulae to calculate T_$ij for a finite step the reader is referred to Ref [1]. By invoking the subroutine ERTRAK (see next Section) the user will have access to the average trajectory, to the full error matrix represented and in addition to this the program makes available also the transport matrix given by which can be useful in several applications (see Section 4.) 3. Description of the User Routines and COMMONs ============================================ To run the program the user should first initialize GEANT, (set-up the geometry and initialize the appropriate physics processes). This procedure will be described in Section 5. The tracking with error propagation is carried out by invoking subr. ERTRAK. However, before calling ERTRAK the user should provide informations to the program in two commons: /ERTRIO/ and the pair /EROPTS/ and /EROPTC/. For this purpose a series of user routines are forseen (routines EUFILL, EUFILP, EUFILV), which should be called by the user. The result of the tracking is partly returned in the arguments of the routine ERTRAK and partly can be accessed through the common /ERTRIO/. In the following we give a description of the user routines (subr. ERTRAK, EUFIL,L,P,V) and that of the commons /ERTRIO/, /EROPTS/ and /EROPTC/. 3.1 User Routines ------------- The output parameters are denoted by asterisk in the calling sequence. SUBROUTINE ERTRAK (X1, P1, X2*, P2*, IPA, CHOPT) ================================================ Performs the tracking of the track from point X1 to point X2 (Before calling this routine the user should also provide the input informations in /EROPTS/, /EROPTC/ and /ERTRIO/ using subr. EUFIL(L/P/V) X1 - Starting coordinates (Cartesian) P1 - Starting 3-momentum (Cartesian) X2 - Final coordinates (Cartesian) P2 - Final 3-momentum (Cartesian) IPA - Particle code (a la GEANT) of the track CHOPT 'B' 'Backward tracking' - i.e. energy loss added to the current energy 'E' 'Exact' calculation of errors assuming helix (i.e. pathlength not assumed as infinitesimal) 'L' Tracking upto prescribed Lengths reached 'M' 'Mixed' prediction (not yet coded) 'O' Tracking 'Only' without calculating errors 'P' Tracking upto prescribed Planes reached 'V' Tracking upto prescribed Volumes reached 'X' Tracking upto prescribed Point approached SUBROUTINE EUFILL (N, EIN, XLF) =============================== User routine to fill the input values of the commons : /EROPTS/, /EROPTC/ and /ERTRIO/ for CHOPT = 'L' N Number of predictions where to store results EIN Input error matrix XLF Defines the tracklengths which if passed the result should be stored SUBROUTINE EUFILP (N, EIN, PLI, PLF) ==================================== User routine to fill the input values of the commons : /EROPTS/, /EROPTC/ and /ERTRIO/ for CHOPT = 'P' N Number of predictions where to store results EIN Input error matrix (in the 'Plane' system ) PLI Defines the start plane PLI(3,1) - and PLI(3,2) - 2 unit vectors in the plane PLF Defines the end plane PLF(3,1,I) - and PLF(3,2,I) - 2 unit vectors in the plane PLF(3,3,I) - point on the plane at intermediate point I SUBROUTINE EUFILV (N, EIN, CNAMV, NUMV, IOVL) ============================================ User routine to fill the input values of the commons : /EROPTS/, /EROPTC/ and /ERTRIO/ for CHOPT = 'V' N Number of predictions where to store results EIN Input error matrix CNAMV Volume name of the prediction NUMV Volume number (if 0 = all volumes) IOVL = 1 prediction when entering in the volume = 2 prediction when leaving the volume 2.2 User COMMONs ------------ CHARACTER*8 CHOPTI PARAMETER (MXPRED = 10) DOUBLE PRECISION ERDTRP LOGICAL LEEXAC, LELENG, LEONLY, LEPLAN, LEPOIN, LEVOLU COMMON /EROPTS/ ERPLI(3,2), ERPLO(3,4,MXPRED), ERLENG(MXPRED), , NAMEER(MXPRED), NUMVER(MXPRED), IOVLER(MXPRED), , LEEXAC, LELENG, LEONLY, LEPLAN, LEPOIN, LEVOLU COMMON /EROPTC/CHOPTI COMMON /ERTRIO/ ERDTRP(5,5,MXPRED), ERRIN(15), ERROUT(15,MXPRED), , ERTRSP(5,5,MXPRED), ERXIN( 3), ERXOUT( 3,MXPRED), , ERPIN(3),ERPOUT(3,MXPRED),NEPRED,INLIST, ILPRED, , IEPRED(MXPRED) LEEXAC = .TRUE. if CHOPT = 'E' in ERTRAK LELENG = .TRUE. if CHOPT = 'L' in ERTRAK LEONLY = .TRUE. if CHOPT = 'O' in ERTRAK LEPLAN = .TRUE. if CHOPT = 'P' in ERTRAK LEPOIN = .TRUE. if CHOPT = 'X' in ERTRAK LEVOLU = .TRUE. if CHOPT = 'V' in ERTRAK IOPTER(I) = 1 if the Ith letter of the alphabet is occuring in CHOPT in ERTRAK (else = 0) NEPRED Number of predictions (c.f. N in EUFILL,P,V) ERPLI Initial plane descriptor (c.f. PLI in EUFILP) ERPLO(,,INLIST) Final plane descriptor - first 3 vectors are identic with PLF in EUFILP, the 4th vector is the cross-product of the first two vectors (plane normal) ERLENG(INLIST) Lengths to store results (c.f. XLF in EUFILL) NAMEER(INLIST) Volume names to store results (c.f. CNAMV in EUFILV) NUMVER(INLIST) Volume numbers to store results (c.f. NUMV in EUFILV) IOVLER(INLIST) (c.f. IOVL in EUFILV) ILPRED Current number of prediction IEPRED(ILPRED) = INLIST if the ILPREDth prediction reached (else = 0) ERDTRP(,,ILPRED) Double precision value of the Transport Matrix at the prediction ILPRED ERRIN Input Error Matrix in Triangular form ERROUT(,ILPRED) Output Error Matrix in Triangular form at the prediction ILPRED ERTRSP(,,ILPRED) Single precision value of the Transport Matrix at the prediction ILPRED ERXIN Starting coordinates (c.f. X1 in ERTRAK) ERXOUT(,ILPRED) Output coordinates at the prediction ILPRED ERPIN Starting momentum ERPOUT(,ILPRED) Output momentum at the prediction ILPRED Note that ERRIN, ERROUT, ERPIN, ERPOUT, ERTRSP and ERDTSP are given by the program in the representation which is requested by CHOPT in subr. ERTRAK. (E.g. if CHOPT='P', all the above quantites are given in the representation III.) 4. Examples of Application ======================= 4.1 The simplest case: Representing the trajectory at another point --------------------------------------------------------------- Usually the particle trajectory is not measured at the point of production where its physical parameters are of interest. Therefore the measurement has to be extrapolated back close to the origin. This can be achieved by the simple call: CALL ERTRAK(X1,P1,X2,P2,IT,CHOPT). Since this extrapolation is opposite to the particle direction, CHOPT should contain the letter 'B'. If the tracking should be stopped on a prescribed plane, CHOPT should also contain 'P', and before invoking ERTRAK the user should call subr. EUFILP. This extrapolation can be carried out simultaneousely onto several planes, in this case the 1st argument of EUFILP is greater than 1. The result can be retrieved from the common /ERTRIO/ as described in Section 3. 4.2 Joining track elements in different parts of the detector --------------------------------------------------------- It happens frequently that one measures a part of a trajectory in a downstream detector and would like to join this information to another one obtained in a detector close to the interaction point. Since there are usually several trajectories which could be a priori joined the first task is to find the one which matches the best. The next task is to improve the trajectory parameters. One chooses a plane near to the interaction point and extrapolates onto this plane all candidate trajectories as described in the preceeding section. For the i-th trajectory one obtains an avarage value x_i and a covariance: sigma_i. (In this discussion the indices will represent the trajectory numbers.) Next one extrapolates back the trajectory from the downstream detector to the same plane and obtains x_d and sigma_d. One can then construct a chi:2 for each track i: chi:2_i = (x_i-x_d)(sigma_i+sigma_d):$-1(x_i-x_d) The matching condition can be defined as: chi:2_i.leq.chi:2_0, where chi:2_0 is some prescribed value. Having chosen trajectory i for the matching the improved track parameters can be obtained by minimizing chi:2 = (x-x_d).sigma_d:$-1.(x-x_d)+(x-x_i).sigma_i:$-1.(x-x_i) w.r.t. x resulting in: x_$impr = (sigma_d:$-1+sigma_i:$-1):$-1. (sigma_d:$-1.x_d+sigma_i:$-1.x_i) The covariance of x_$impr sigma_$impr = (sigma_d:$-1+sigma_i:$-1):$-1 shows explicitely the improvement of the trajectory parameters. This procedure can be easily generalized to join more than 2 measurements on the particle trajectory. If e.g. between the two above planes there is another detection plane, one can first merge the informations of the downstream and intermediate plane and continue the backtracking from the intermediate plane to the plane close to the interaction point with the improved trajectory parameters. The procedure can be used in principle also if not all the five parameters are measured (e.g. if only the coordinate informations are available). In this case one starts the back-tracking with some initial values of the unmeasured parameters and assigns an error to these parameters which is much larger than the difference between the true and the initial value. The user is however has to ensure that the result is stable against the choice of the starting value of parameters and errors (e.g. by performing several iterations). These problems can be overcome by a fitting procedure which is described in Section 4.4. 4.3 Prediction of the trajectory ---------------------------- It is often needed to predict the particle trajectory in a detector plane at a certain confidence level in order to perform pattern recognition. An example is to find hits from a penetrating particle inside a segmented calorimeter when the particle trajectory is well measured at the two extrems of the calorimeter. In the case of 1 intermediate plane inside the calorimeter the solution can be obtained by combining the methods outlined in the previous two sections. One extrapolates the measured track parameters from the two endplanes of the calorimeter onto the plane set up inside the calorimeter (Section 3.1) and one joins the two informations on that plane (Section 3.2). This procedure of course can be carried out on any number of intermediate planes. However, if there is a large number of planes, it is advantageous to carry out the tracking in one direction and in one go, for which case a method is outlined below. Let's start the tracking from one end of the calorimeter and denote by x_i and by x_e the average track parameters on the intermediate plane i and on the other end-plane e, respectively. Let's denote the true track parameters on the same planes by x and by y, respectively. The corresponding chi:2, which we should minimize w.r.t. x is: chi:2 = (x-x_i).sigma_i:$-1.(x-x_i) + + (y(x)-x_m).sigma_$em:$-1.(y(x)-x_m) where x_m is the measured trajectory at the end-plane with covariance matrix sigma_m, sigma_i is the propagated error matrix from the starting plane to plane i and sigma_$em = sigma_m + sigma_$ei where sigma_$ei is the propagated error from plane i to plane e (n o t including the error on plane i itself). The minimization results in the following equation: sigma_i:$-1.(x-x_i) + dy/dx.sigma_$em:$-1.(y(x)-x_m) = 0, which we solve by linearization: x = x_i + Delta x y = y(x_i)+dy(x_i)/dx.Delta x = x_e+T(e,i).Delta x. The result is: Delta x = sigma_x.T:T(e,i).sigma_$em:$-1.(x_m-x_e) where: sigma_x = (sigma_i:$-1+T:T(e,i).sigma_$em:$-1.T(e,i)):$-1 is the covariance matrix of the trajectory prediction at the plane i (T:T means the transpose of T). The following glossary gives the correspondance between the mathematical quantities used in the above equations and the varibales in the user common /ERTRIO/: x_i ERXOUT(,I), ERPOUT(,I) (I standing for prediction i) x_ ERXOUT(,IE), ERPOUT(,IE) (IE standing for prediction e) sigma_i ERROUT(,I) T(e,i) T(e,1).T:$-1(i,1) T(e,1) ERTRSP(,,IE) T(i,1) ERTRSP(,,I) sigma_$ei sigma_e-T(e,i).sigma_iT(e,i) sigma_e ERROUT(,IE) 4.4 Fitting trajectory parameters ----------------------------- In the above examples all of the 5 variables of the trajectory have been known at least in one space point. However, in most of the cases direct mesurements yield only the coordinate informations, from which one should reconstruct the curvature and the direction. The following example shows how to use the program package for this pur- pose. This tool can be applied in the most general case: in inhomoge- neous magnetic field and even if the particle passes through a great amount of material. Suppose we would like to reconstruct the particle trajectory x_0 at plane 0 by measuring the coordinate informations x_i:m at N different detector planes (i=1,...,N). If in the formation of the trajectory random processes can be neglected, then the average trajectory can be obtained by minimizing the following chi:2 w.r.t. x_0: chi:2 = Sum_$i=1:N[(x_i(x_0)-x_i:m).sigma_i:m:$-1.(x_i(x_0)-x_i:m)] where x_i are the true track parameters at plane i, and sigma_i:m is the 2 X 2 covariance matrix of the measurement on plane i. This results in the following equation: Sum_$i=1:NT:T(i,0).sigma_i:m:$-1.(x_i(x_0)-x_i:m) = 0 where T(i,0) is the transport matrix between plane 0 and plane i. This equation is again solved by linearization. In the first approximation we calculate the trial trajectory x_i:t on plane i starting with a value x_0:t. The true value is then obtained by: x_0 = x_0:t + Delta x_0 with Delta x_0 = sigma_$x_0.Sum_$i=1:NT:T(i,0).sigma_i:m:$-1.(x_i:m-x_i:t), where the covariance matrix of x_0 is given by: sigma_$x_0 = Sum_$i=1:NT:T(i,0).sigma_i:m:$-1.T(i,0)]:$-1. If in the formation of the particle trajectory random processes, like multiple Coulomb scattering, cannot be neglected then obviousely there are correlations in the error matrix between different planes and therefore the above chi:2 should be written as chi:2 = (x(x_0)-x:m).sigma:$-1.(x(x_0)-x:m) where x is a vector of length 2 X N containing the coordinate values (xi,eta) of the average trajectory plane by plane: (xi_1,eta_1,xi_2,eta_2,...,xi_N,eta_N), x:m is the corresponding vector of the measured coordinates. sigma is a 2N X 2N matrix, whose 2 X 2 diagonal submatrices can be written as sigma_$ii = sigma_i:m + sigma(2)_i:r where sigma(2)_i:r is the 2 X 2 part of the covariance matrix sigma_i:r due to random processes. The off-diagonal 2 X 2 matrices give the the correlations between planes: sigma_$ij = T(j,i).sigma_i:r (i GUFLD | ----> TRSCSD | ----> TRSDSC | ----> ERBCER | ----> GEKBIN | ----> ERPINI ----> TRPROP | ----> ERTRGO ====== | | ----> GMEDIA | ----> GUFLD | ----> EVOLIO | ----> ERSTOR | ====== | | | | ----> ERBCER | | ----> ERBCTR | | ----> TRSCSD | | ----> DMMSS | | | ----> ERTRCH -| ----> GTNEXT | ====== | | ----> ERTRNT -| ----> GUSWIM ----> GUFLD | ====== | ----> GINVOL | | ----> ERLAND | | ----> GEKBIN | | ----> ERPROP | | ====== | ----> EUSTEP | | | ----> EVOLIO | | ----> GUFLD l | | ----> TRPROP l ----> GTMEDI | | ----> TRPRFN | | ----> SSMTST | | ----> ERMCSC | | ----> ERSTOR ====== 6. Acknowledgements ================ The authors of the present interface benefitted numerous critical remarks and useful suggestions from the authors of the GEANT3 Package, especially from F. Bruyant (CERN), which are greatly acknowledged here. Complaints and suggest must be sent to one of the authors : Innocent@cernvm Maire@cernvm Nagy@cernvm References: =========== [1] W.Wittek, EMC Internal Reports: EMC/80/15, EMCSW/80/39, EMCSW/81/13, EMCSW/81/18 A.Haas, The EMC Utility Package: UTIL42 [2] R.Brun, F.Bruyant, M.Maire, A.C.McPherson, P.Zanarini DD/EE/84-1, May 1986