]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1 1996/03/06 15:31:14 mclareni | |
6 | * Add geane321 history, CMZ and doc files | |
7 | * | |
8 | * | |
9 | *CMZ : 3.21/02 29/03/94 15.41.49 by S.Giani | |
10 | *-- Author : | |
11 | ||
12 | *-- Author : | |
13 | 1. Introduction | |
14 | ||
15 | 2. Definitions | |
16 | ||
17 | 3. Description of the User Routines and COMMONs | |
18 | ||
19 | 4. Examples of application | |
20 | ||
21 | 5. Interface with GEANT | |
22 | ||
23 | 6. Acknowledgements | |
24 | ||
25 | ||
26 | ||
27 | 1. Introduction | |
28 | ============ | |
29 | ||
30 | The present Package allows the user to calculate the | |
31 | average trajectories of particles and to calculate the | |
32 | transport matrix as well as the propagated error covariance | |
33 | matrix. It makes use of a set of routines worked out | |
34 | by the European Muon Collaboration [1] and it is integrated | |
35 | to the GEANT3 system [2] with expected applications in both | |
36 | simulation and reconstruction context. | |
37 | ||
38 | The package is available as a PAM-file. It contains | |
39 | two basic PATCHes, one with the original EMC routines, | |
40 | the other with new GEANT-type tracking routines and interface | |
41 | between them and the EMC routines. This second PATCH | |
42 | contains those routines the user should invoke | |
43 | to carry out the tracking (routines ERTRAK, EUFILL, EUFILP, EUFILV). | |
44 | In addition to this a series of utilities are available for | |
45 | the user (e.g. to transform the track representation from one | |
46 | system to another or to carry out 5 X 5 matrix multiplication | |
47 | in an optimum way). | |
48 | ||
49 | In Section 2 we give the definitions of the mathematical | |
50 | quantities to be dealt with. In Section 3 the description of | |
51 | the user routines are given. In Section 4 we illustrate the | |
52 | application of the program by several examples. Finally in | |
53 | Section 5 the GEANT part is described. | |
54 | ||
55 | Further development will concern the improvement of the | |
56 | error matrix by taking into account the Landau tail in the | |
57 | fluctuation of the energy loss, the bremsstrahlung and direct | |
58 | pair production of the muons. | |
59 | ||
60 | ||
61 | 2. Definitions | |
62 | =========== | |
63 | ||
64 | 2.1 Track variables, Representations | |
65 | -------------------------------- | |
66 | ||
67 | The particle trajectory is characterized by 5 independent | |
68 | variables as a function of one parameter (e.g. the pathlength). | |
69 | Among the 5 variables 1 is related to the curvature (to the absolute | |
70 | value of the momentum, p), 2 are related to the direction of the | |
71 | particle and the other 2 are related to the spatial location. | |
72 | The most usual representation of these 5 parameters are: | |
73 | ||
74 | I. | |
75 | ||
76 | 1/p, lambda, phi, y_perp, z_perp | |
77 | ||
78 | ||
79 | where lambda and phi are the dip and azimuthal angles related | |
80 | to the momentum components in the following way: | |
81 | ||
82 | p_x = p cos(lambda) cos(phi) | |
83 | p_y = p cos(lambda) sin(phi) | |
84 | p_z = p sin(lambda) | |
85 | ||
86 | ||
87 | y_perp and z_perp are the coordinates of the trajectory in a | |
88 | local orthonormal reference frame with the x_perp axis along the | |
89 | particle direction, the y_perp being parallel to the x-y plane. | |
90 | This representation is usually applied in the overall reference | |
91 | frame. (In the EMC code this reference frame is labelled by 'SC' | |
92 | since the overall system was identified with that of the Streamer | |
93 | Chamber.) | |
94 | ||
95 | II. | |
96 | ||
97 | 1/p, y', z', y, z | |
98 | ||
99 | ||
100 | where y'=dy/dx and z'=dz/dx. This representation is particularly | |
101 | useful in fixed target experiments, where the trajectory is evaluated | |
102 | on successive parallel planes (which are perpendicular to the x-axis). | |
103 | (In the EMC code this representation is labelled by 'SP' since a | |
104 | convenient mathematical description of a trajectory being approxima- | |
105 | tely parallel to the x-axis is a 'spline'.) | |
106 | ||
107 | III. | |
108 | ||
109 | 1/p, v', w', v, w | |
110 | ||
111 | ||
112 | where v'=dv/du and w'=dw/du in an orthonormal coordinate system with | |
113 | axis u, v and w. This representation is paricularly useful when the | |
114 | trajectory has to be evaluated on different detector planes | |
115 | in a colliding beam experiment, where the planes can take a great | |
116 | variety of directions.(In the EMC code this representation is | |
117 | labelled by 'SD' as System of Detection.) | |
118 | ||
119 | Of course, all the above representations of the trajectory | |
120 | are equivalent and one can go from one representation to the | |
121 | other by calculating the corresponding Jacobian. These Jacobians | |
122 | are provided by the following EMC routines: | |
123 | ||
124 | S/R TRSCSP from I to II | |
125 | S/R TRSPSC from II to I | |
126 | S/R TRSCSD from I to III | |
127 | S/R TRSDSC from III to I | |
128 | ||
129 | ||
130 | ||
131 | 2.2 Error Propagation | |
132 | ----------------- | |
133 | ||
134 | ||
135 | Let us denote in the following the 5 independent variables at | |
136 | a given value of parameter l_0 (e.g. pathlength) by x_i(l_0), | |
137 | (i=1,...,5). In many applications we are interested in the evolution | |
138 | of the average value of x_i for l>l_0: E(x_i). This is calculated | |
139 | by GEANT as will be outlined in Section 5. | |
140 | ||
141 | The knowledge on the avarage trajectory is characterized by the | |
142 | 5 X 5 covariance matrix of the variables: | |
143 | ||
144 | sigma_$ij\a(l_0) = E(x_i(l_0).x_j(l_0)) - E(x_i(l_0)).E(x_j(l_0)) | |
145 | ||
146 | We are also interested in the evolution of sigma_$ij\a for l>l_0, | |
147 | which we call error propagation. If the particle is propagating in a | |
148 | deterministic way, i.e. without any random process involved ( | |
149 | e.g. in vacuum) then the propagation of sigma is simply described by | |
150 | the so called transport matrix in the following way: | |
151 | ||
152 | sigma_$ij\a(l) = T_$jm\a(l,l_0).sigma_$mn\a(l_0)T_$in\a(l,l_0) | |
153 | ||
154 | where the transport matrix expresses the infinitesimal change | |
155 | in the parameters at l with respect to the change of parameters | |
156 | at l_0: | |
157 | ||
158 | T_$ij\a(l,l_0) = delta (x_i(l))/ delta (x_j(l_0)). | |
159 | ||
160 | In a realistic detector, however, the particles undergo random | |
161 | processes as well, like Multiple Coulomb scattering, energy loss | |
162 | due to delta ray production, etc. therefore the error propagation | |
163 | should contain an additional term: | |
164 | ||
165 | sigma_$ij\a(l) = T_$jm\a(l,l_0).sigma_$mn\a(l_0).T_$in\a(l,l_0) + | |
166 | + sigma_$ij\a:$random\a(l). | |
167 | ||
168 | The program calculates sigma_$ij\a(l) step by step using the above | |
169 | recursive formula, where T_$ij\a and sigma_$ij\a:$random\a refers to the | |
170 | actual step and sigma_$mn\a is the cumulative error for all previous | |
171 | steps. For the mathematical formulae to calculate T_$ij\a for a finite | |
172 | step the reader is referred to Ref [1]. | |
173 | ||
174 | By invoking the subroutine ERTRAK (see next Section) the user | |
175 | will have access to the average trajectory, to the full error matrix | |
176 | represented and in addition to this the program makes available also | |
177 | the transport matrix given by which can be useful in several | |
178 | applications (see Section 4.) | |
179 | ||
180 | ||
181 | ||
182 | 3. Description of the User Routines and COMMONs | |
183 | ============================================ | |
184 | ||
185 | To run the program the user should first initialize GEANT, | |
186 | (set-up the geometry and initialize the appropriate physics | |
187 | processes). This procedure will be described in Section 5. | |
188 | The tracking with error propagation is carried out by invoking | |
189 | subr. ERTRAK. However, before calling ERTRAK the user should | |
190 | provide informations to the program in two commons: | |
191 | /ERTRIO/ and the pair /EROPTS/ and /EROPTC/. For this purpose | |
192 | a series of user routines | |
193 | are forseen (routines EUFILL, EUFILP, EUFILV), which should be | |
194 | called by the user. The result of the tracking is partly | |
195 | returned in the arguments of the routine ERTRAK and partly can be | |
196 | accessed through the common /ERTRIO/. In the following we give | |
197 | a description of the user routines (subr. ERTRAK, EUFIL,L,P,V) | |
198 | and that of the commons /ERTRIO/, /EROPTS/ and /EROPTC/. | |
199 | ||
200 | ||
201 | 3.1 User Routines | |
202 | ------------- | |
203 | ||
204 | The output parameters are denoted by asterisk in the calling | |
205 | sequence. | |
206 | ||
207 | SUBROUTINE ERTRAK (X1, P1, X2*, P2*, IPA, CHOPT) | |
208 | ================================================ | |
209 | ||
210 | ||
211 | Performs the tracking of the track from point X1 to | |
212 | point X2 | |
213 | (Before calling this routine the user should also provide | |
214 | the input informations in /EROPTS/, /EROPTC/ and /ERTRIO/ | |
215 | using subr. EUFIL(L/P/V) | |
216 | ||
217 | X1 - Starting coordinates (Cartesian) | |
218 | P1 - Starting 3-momentum (Cartesian) | |
219 | X2 - Final coordinates (Cartesian) | |
220 | P2 - Final 3-momentum (Cartesian) | |
221 | IPA - Particle code (a la GEANT) of the track | |
222 | ||
223 | CHOPT | |
224 | 'B' 'Backward tracking' - i.e. energy loss | |
225 | added to the current energy | |
226 | 'E' 'Exact' calculation of errors assuming | |
227 | helix (i.e. pathlength not | |
228 | assumed as infinitesimal) | |
229 | 'L' Tracking upto prescribed Lengths reached | |
230 | 'M' 'Mixed' prediction (not yet coded) | |
231 | 'O' Tracking 'Only' without calculating errors | |
232 | 'P' Tracking upto prescribed Planes reached | |
233 | 'V' Tracking upto prescribed Volumes reached | |
234 | 'X' Tracking upto prescribed Point approached | |
235 | ||
236 | ||
237 | SUBROUTINE EUFILL (N, EIN, XLF) | |
238 | =============================== | |
239 | ||
240 | ||
241 | User routine to fill the input values of the commons : | |
242 | /EROPTS/, /EROPTC/ and /ERTRIO/ for CHOPT = 'L' | |
243 | ||
244 | N Number of predictions where to store results | |
245 | EIN Input error matrix | |
246 | XLF Defines the tracklengths which if passed the | |
247 | result should be stored | |
248 | ||
249 | SUBROUTINE EUFILP (N, EIN, PLI, PLF) | |
250 | ==================================== | |
251 | ||
252 | ||
253 | User routine to fill the input values of the commons : | |
254 | /EROPTS/, /EROPTC/ and /ERTRIO/ for CHOPT = 'P' | |
255 | N Number of predictions where to store results | |
256 | EIN Input error matrix (in the 'Plane' system ) | |
257 | PLI Defines the start plane | |
258 | PLI(3,1) - and | |
259 | PLI(3,2) - 2 unit vectors in the plane | |
260 | PLF Defines the end plane | |
261 | PLF(3,1,I) - and | |
262 | PLF(3,2,I) - 2 unit vectors in the plane | |
263 | PLF(3,3,I) - point on the plane | |
264 | at intermediate point I | |
265 | ||
266 | ||
267 | SUBROUTINE EUFILV (N, EIN, CNAMV, NUMV, IOVL) | |
268 | ============================================ | |
269 | ||
270 | ||
271 | User routine to fill the input values of the commons : | |
272 | /EROPTS/, /EROPTC/ and /ERTRIO/ for CHOPT = 'V' | |
273 | N Number of predictions where to store results | |
274 | EIN Input error matrix | |
275 | CNAMV Volume name of the prediction | |
276 | NUMV Volume number (if 0 = all volumes) | |
277 | IOVL = 1 prediction when entering in the volume | |
278 | = 2 prediction when leaving the volume | |
279 | ||
280 | ||
281 | 2.2 User COMMONs | |
282 | ------------ | |
283 | ||
284 | CHARACTER*8 CHOPTI | |
285 | PARAMETER (MXPRED = 10) | |
286 | DOUBLE PRECISION ERDTRP | |
287 | LOGICAL LEEXAC, LELENG, LEONLY, LEPLAN, LEPOIN, LEVOLU | |
288 | ||
289 | COMMON /EROPTS/ ERPLI(3,2), ERPLO(3,4,MXPRED), ERLENG(MXPRED), | |
290 | , NAMEER(MXPRED), NUMVER(MXPRED), IOVLER(MXPRED), | |
291 | , LEEXAC, LELENG, LEONLY, LEPLAN, LEPOIN, LEVOLU | |
292 | COMMON /EROPTC/CHOPTI | |
293 | ||
294 | COMMON /ERTRIO/ ERDTRP(5,5,MXPRED), ERRIN(15), ERROUT(15,MXPRED), | |
295 | , ERTRSP(5,5,MXPRED), ERXIN( 3), ERXOUT( 3,MXPRED), | |
296 | , ERPIN(3),ERPOUT(3,MXPRED),NEPRED,INLIST, ILPRED, | |
297 | , IEPRED(MXPRED) | |
298 | ||
299 | ||
300 | LEEXAC = .TRUE. if CHOPT = 'E' in ERTRAK | |
301 | LELENG = .TRUE. if CHOPT = 'L' in ERTRAK | |
302 | LEONLY = .TRUE. if CHOPT = 'O' in ERTRAK | |
303 | LEPLAN = .TRUE. if CHOPT = 'P' in ERTRAK | |
304 | LEPOIN = .TRUE. if CHOPT = 'X' in ERTRAK | |
305 | LEVOLU = .TRUE. if CHOPT = 'V' in ERTRAK | |
306 | IOPTER(I) = 1 if the Ith letter of the alphabet is | |
307 | occuring in CHOPT in ERTRAK (else = 0) | |
308 | ||
309 | NEPRED Number of predictions (c.f. N in EUFILL,P,V) | |
310 | ERPLI Initial plane descriptor (c.f. PLI in EUFILP) | |
311 | ERPLO(,,INLIST) Final plane descriptor - first 3 vectors are identic | |
312 | with PLF in EUFILP, the 4th vector is the cross-product | |
313 | of the first two vectors (plane normal) | |
314 | ERLENG(INLIST) Lengths to store results (c.f. XLF in EUFILL) | |
315 | NAMEER(INLIST) Volume names to store results (c.f. CNAMV in EUFILV) | |
316 | NUMVER(INLIST) Volume numbers to store results (c.f. NUMV in EUFILV) | |
317 | IOVLER(INLIST) (c.f. IOVL in EUFILV) | |
318 | ||
319 | ILPRED Current number of prediction | |
320 | IEPRED(ILPRED) = INLIST if the ILPREDth prediction reached (else = 0) | |
321 | ERDTRP(,,ILPRED) Double precision value of the Transport Matrix | |
322 | at the prediction ILPRED | |
323 | ERRIN Input Error Matrix in Triangular form | |
324 | ERROUT(,ILPRED) Output Error Matrix in Triangular form | |
325 | at the prediction ILPRED | |
326 | ERTRSP(,,ILPRED) Single precision value of the Transport Matrix | |
327 | at the prediction ILPRED | |
328 | ERXIN Starting coordinates (c.f. X1 in ERTRAK) | |
329 | ERXOUT(,ILPRED) Output coordinates at the prediction ILPRED | |
330 | ERPIN Starting momentum | |
331 | ERPOUT(,ILPRED) Output momentum at the prediction ILPRED | |
332 | ||
333 | ||
334 | ||
335 | Note that ERRIN, ERROUT, ERPIN, ERPOUT, ERTRSP and ERDTSP are | |
336 | given by the program in the representation which is requested by | |
337 | CHOPT in subr. ERTRAK. (E.g. if CHOPT='P', all the above quantites | |
338 | are given in the representation III.) | |
339 | ||
340 | ||
341 | ||
342 | 4. Examples of Application | |
343 | ======================= | |
344 | ||
345 | ||
346 | 4.1 The simplest case: Representing the trajectory at another point | |
347 | --------------------------------------------------------------- | |
348 | ||
349 | Usually the particle trajectory is not measured at the point of | |
350 | production where its physical parameters are of interest. Therefore | |
351 | the measurement has to be extrapolated back close to the origin. This | |
352 | can be achieved by the simple call: | |
353 | ||
354 | CALL ERTRAK(X1,P1,X2,P2,IT,CHOPT). | |
355 | ||
356 | Since this extrapolation is opposite to the particle direction, | |
357 | CHOPT should contain the letter 'B'. If the tracking should be | |
358 | stopped on a prescribed plane, CHOPT should also contain 'P', | |
359 | and before invoking ERTRAK the user should call subr. EUFILP. | |
360 | This extrapolation can be carried out simultaneousely onto several | |
361 | planes, in this case the 1st argument of EUFILP is greater than 1. | |
362 | The result can be retrieved from the common /ERTRIO/ as described | |
363 | in Section 3. | |
364 | ||
365 | 4.2 Joining track elements in different parts of the detector | |
366 | --------------------------------------------------------- | |
367 | ||
368 | It happens frequently that one measures a part of a trajectory | |
369 | in a downstream detector and would like to join this information | |
370 | to another one obtained in a detector close to the interaction | |
371 | point. Since there are usually several trajectories which could be | |
372 | a priori joined the first task is to find the one which matches | |
373 | the best. The next task is to improve the trajectory parameters. | |
374 | ||
375 | One chooses a plane near to the interaction point and extrapolates | |
376 | onto this plane all candidate trajectories as described in the | |
377 | preceeding section. For the i-th trajectory one obtains an | |
378 | avarage value x_i and a covariance: sigma_i. (In this discussion | |
379 | the indices will represent the trajectory numbers.) Next one | |
380 | extrapolates back the trajectory from the downstream detector | |
381 | to the same plane and obtains x_d and sigma_d. One can then construct | |
382 | a chi:2 for each track i: | |
383 | ||
384 | chi:2_i = (x_i-x_d)(sigma_i+sigma_d):$-1\a(x_i-x_d) | |
385 | ||
386 | The matching condition can be defined as: | |
387 | ||
388 | chi:2_i.leq.chi:2_0, | |
389 | ||
390 | where chi:2_0 is some prescribed value. | |
391 | ||
392 | Having chosen trajectory i for the matching the improved track | |
393 | parameters can be obtained by minimizing | |
394 | ||
395 | chi:2 = (x-x_d).sigma_d:$-1\a.(x-x_d)+(x-x_i).sigma_i:$-1\a.(x-x_i) | |
396 | ||
397 | w.r.t. x resulting in: | |
398 | ||
399 | x_$impr\a = (sigma_d:$-1\a+sigma_i:$-1\a):$-1\a. | |
400 | (sigma_d:$-1\a.x_d+sigma_i:$-1\a.x_i) | |
401 | ||
402 | The covariance of x_$impr\a | |
403 | ||
404 | sigma_$impr\a = (sigma_d:$-1\a+sigma_i:$-1\a):$-1\a | |
405 | ||
406 | shows explicitely the improvement of the trajectory parameters. | |
407 | ||
408 | This procedure can be easily generalized to join more than 2 | |
409 | measurements on the particle trajectory. If e.g. between the two | |
410 | above planes there is another detection plane, one can first | |
411 | merge the informations of the downstream and intermediate | |
412 | plane and continue the backtracking from the intermediate | |
413 | plane to the plane close to the interaction point with the | |
414 | improved trajectory parameters. | |
415 | ||
416 | The procedure can be used in principle also if not all | |
417 | the five parameters are measured (e.g. if only the | |
418 | coordinate informations are available). In this case one | |
419 | starts the back-tracking with some initial values of the | |
420 | unmeasured parameters and assigns an error to these | |
421 | parameters which is much larger than the difference between | |
422 | the true and the initial value. The user is however has to | |
423 | ensure that the result is stable against the choice of | |
424 | the starting value of parameters and errors (e.g. by performing | |
425 | several iterations). These problems can be overcome | |
426 | by a fitting procedure which is described in Section 4.4. | |
427 | ||
428 | 4.3 Prediction of the trajectory | |
429 | ---------------------------- | |
430 | ||
431 | It is often needed to predict the particle trajectory in a | |
432 | detector plane at a certain confidence level in order to | |
433 | perform pattern recognition. An example is to find hits from | |
434 | a penetrating particle inside a segmented calorimeter when the | |
435 | particle trajectory is well measured at the two extrems of the | |
436 | calorimeter. | |
437 | ||
438 | In the case of 1 intermediate plane inside the calorimeter | |
439 | the solution can be obtained by combining the methods outlined | |
440 | in the previous two sections. One extrapolates the measured track | |
441 | parameters from the two endplanes of the calorimeter onto the plane | |
442 | set up inside the calorimeter (Section 3.1) and one joins the two | |
443 | informations on that plane (Section 3.2). This procedure of course | |
444 | can be carried out on any number of intermediate planes. However, | |
445 | if there is a large number of planes, it is advantageous to carry | |
446 | out the tracking in one direction and in one go, for which case a | |
447 | method is outlined below. | |
448 | ||
449 | Let's start the tracking from one end of the calorimeter | |
450 | and denote by x_i and by x_e the average track parameters on the | |
451 | intermediate plane i and on the other end-plane e, respectively. | |
452 | Let's denote the true track parameters on the same planes by x | |
453 | and by y, respectively. The corresponding chi:2, which we should | |
454 | minimize w.r.t. x is: | |
455 | ||
456 | chi:2 = (x-x_i).sigma_i:$-1\a.(x-x_i) + | |
457 | + (y(x)-x_m).sigma_$em\a:$-1\a.(y(x)-x_m) | |
458 | ||
459 | where x_m is the measured trajectory at the end-plane with | |
460 | covariance matrix sigma_m, sigma_i is the propagated error | |
461 | matrix from the starting plane to plane i and | |
462 | ||
463 | sigma_$em\a = sigma_m + sigma_$ei\a | |
464 | ||
465 | where sigma_$ei\a is the propagated error from plane i to plane e | |
466 | (n o t including the error on plane i itself). | |
467 | ||
468 | The minimization results in the following equation: | |
469 | ||
470 | sigma_i:$-1\a.(x-x_i) + dy/dx.sigma_$em\a:$-1\a.(y(x)-x_m) = 0, | |
471 | ||
472 | which we solve by linearization: | |
473 | ||
474 | x = x_i + Delta x | |
475 | ||
476 | y = y(x_i)+dy(x_i)/dx.Delta x = x_e+T(e,i).Delta x. | |
477 | ||
478 | The result is: | |
479 | ||
480 | Delta x = sigma_x.T:T(e,i).sigma_$em\a:$-1\a.(x_m-x_e) | |
481 | ||
482 | where: | |
483 | ||
484 | sigma_x = (sigma_i:$-1\a+T:T(e,i).sigma_$em\a:$-1.\aT(e,i)):$-1\a | |
485 | ||
486 | is the covariance matrix of the trajectory prediction at the | |
487 | plane i (T:T means the transpose of T). | |
488 | ||
489 | The following glossary gives the correspondance between the | |
490 | mathematical quantities used in the above equations and the | |
491 | varibales in the user common /ERTRIO/: | |
492 | ||
493 | x_i ERXOUT(,I), ERPOUT(,I) (I standing for prediction i) | |
494 | x_ ERXOUT(,IE), ERPOUT(,IE) (IE standing for prediction e) | |
495 | sigma_i ERROUT(,I) | |
496 | T(e,i) T(e,1).T:$-1\a(i,1) | |
497 | T(e,1) ERTRSP(,,IE) | |
498 | T(i,1) ERTRSP(,,I) | |
499 | sigma_$ei\a sigma_e-T(e,i).sigma_iT(e,i) | |
500 | sigma_e ERROUT(,IE) | |
501 | ||
502 | ||
503 | ||
504 | ||
505 | 4.4 Fitting trajectory parameters | |
506 | ----------------------------- | |
507 | ||
508 | ||
509 | In the above examples all of the 5 variables of the trajectory | |
510 | have been known at least in one space point. However, in most of the | |
511 | cases direct mesurements yield only the coordinate informations, from | |
512 | which one should reconstruct the curvature and the direction. The | |
513 | following example shows how to use the program package for this pur- | |
514 | pose. This tool can be applied in the most general case: in inhomoge- | |
515 | neous magnetic field and even if the particle passes through a great | |
516 | amount of material. | |
517 | ||
518 | Suppose we would like to reconstruct the particle trajectory | |
519 | x_0 at plane 0 by measuring the coordinate informations x_i:m | |
520 | at N different detector planes (i=1,...,N). If in the formation | |
521 | of the trajectory random processes can be neglected, then | |
522 | the average trajectory can be obtained by | |
523 | minimizing the following chi:2 w.r.t. x_0: | |
524 | ||
525 | chi:2 = Sum_$i=1\a:N[(x_i(x_0)-x_i:m).sigma_i:m:$-1\a.(x_i(x_0)-x_i:m)] | |
526 | ||
527 | where x_i are the true track parameters at plane i, and sigma_i:m is | |
528 | the 2 X 2 covariance matrix of the measurement on plane i. This results | |
529 | in the following equation: | |
530 | ||
531 | Sum_$i=1\a:NT:T(i,0).sigma_i:m:$-1\a.(x_i(x_0)-x_i:m) = 0 | |
532 | ||
533 | where T(i,0) is the transport matrix between plane 0 and plane i. | |
534 | This equation is again solved by linearization. In the first | |
535 | approximation we calculate the trial trajectory x_i:t on plane i | |
536 | starting with a value x_0:t. The true value is then obtained by: | |
537 | ||
538 | x_0 = x_0:t + Delta x_0 | |
539 | with | |
540 | ||
541 | Delta x_0 = sigma_$x_0\a.Sum_$i=1\a:NT:T(i,0).sigma_i:m:$-1\a.(x_i:m-x_i:t), | |
542 | ||
543 | where the covariance matrix of x_0 is given by: | |
544 | ||
545 | sigma_$x_0\a = Sum_$i=1\a:NT:T(i,0).sigma_i:m:$-1\a.T(i,0)]:$-1\a. | |
546 | ||
547 | If in the formation of the particle trajectory random processes, | |
548 | like multiple Coulomb scattering, cannot be neglected then obviousely | |
549 | there are correlations in the error matrix between different planes | |
550 | and therefore the above chi:2 should be written as | |
551 | ||
552 | chi:2 = (x(x_0)-x:m).sigma:$-1\a.(x(x_0)-x:m) | |
553 | ||
554 | where x is a vector of length 2 X N containing the coordinate values | |
555 | (xi,eta) of the average trajectory plane by plane: | |
556 | ||
557 | (xi_1,eta_1,xi_2,eta_2,...,xi_N,eta_N), | |
558 | ||
559 | x:m is the corresponding vector of the measured coordinates. | |
560 | sigma is a 2N X 2N matrix, whose 2 X 2 diagonal submatrices can be | |
561 | written as | |
562 | ||
563 | sigma_$ii\a = sigma_i:m + sigma(2)_i:r | |
564 | ||
565 | where sigma(2)_i:r is the 2 X 2 part of the covariance matrix sigma_i:r | |
566 | due to random processes. The off-diagonal 2 X 2 matrices give the the | |
567 | correlations between planes: | |
568 | ||
569 | sigma_$ij\a = T(j,i).sigma_i:r (i<j) | |
570 | ||
571 | where T(j,i) is the 5 X 5 transport matrix between plane i and j and | |
572 | only the 2 X 2 part of sigma_$ij\a is considered. | |
573 | ||
574 | The minimization procedure is formally the same as before: | |
575 | ||
576 | x_0 = x_0:t + Delta x_0 | |
577 | ||
578 | with | |
579 | ||
580 | Delta x_0 = sigma_$x_0\a.tau:T.sigma:$-1\a.(x:m-x:t), | |
581 | ||
582 | where the 5 X 5 covariance matrix of x_0 is given by: | |
583 | ||
584 | sigma_$x_0\a = (tau:T.sigma:$-1\a.tau):$-1\a. | |
585 | ||
586 | Here tau is the joint transport matrix of dimension 5 X 2N | |
587 | containing the transport matrices T(i,0) (5 variables for | |
588 | plane 0 and 2 coordinate variables for plane i). | |
589 | ||
590 | Again the correspondance between the above symbols and variables | |
591 | calculated by the program package is given by the following glossary: | |
592 | ||
593 | In case of tracking from plane 0 to plane i in one go: | |
594 | ||
595 | ||
596 | x_i:t ERXOUT(,I), ERPOUT(,I) (I standing for prediction i) | |
597 | sigma_i:r ERROUT(,I) (starting with 0 error) | |
598 | T(j,i) T(j,0).T:$-1\a(i,0) | |
599 | T(i,0) ERTRSP(,,I) | |
600 | ||
601 | ||
602 | ||
603 | ||
604 | 5. Interface with GEANT | |
605 | ==================== | |
606 | ||
607 | ||
608 | The following diagram shows the program flow-chart of ERTRAK. | |
609 | Two running examples show how to build up a complete application. | |
610 | (see patch EREXAM1 and EREXAM2) | |
611 | ||
612 | As can be seen from the flow-chart the particle is propagated | |
613 | through the experimental setup by the routine ERTRGO which is a | |
614 | simplified version of the GEANT routines GTRACK and GTVOL. In fact, | |
615 | the calculation of the average trajectory is independent of the basic | |
616 | GEANT Tracking Package. | |
617 | ||
618 | On the other hand, the GEANT Geometry Package is invoked | |
619 | by calling GMEDIA,GTMEDI,GINVOL,GTNEXT (see the flowchart). | |
620 | Therefore the GEANT Data Structures: JVOLU, JMATE, JTMED, JPART | |
621 | must be set up at initialization time (see patch erexample below). | |
622 | ||
623 | In fact, in the JMATE data structure, only the energy loss | |
624 | tables (JLOSS and JRANG) are required. The energy loss calculation must | |
625 | include the contributions coming from the undetectable random | |
626 | processes: delta rays, bremsstrahlung and direct pair production | |
627 | of muons. The user should check that the values of the parameters | |
628 | DCUTE/M, BCUTE/M, PPCUTM have been set coherently with his | |
629 | application. (In most of the cases these parameters should be set | |
630 | to their maximum value (say 10 TeV).) | |
631 | ||
632 | No other GEANT facilities are required to run ERTRAK. | |
633 | No GEANT data structurs are overwritten. Only the commons /GCKINE/ | |
634 | and /GCTRAK/ are used to keep the current values of the parameters | |
635 | of the average trajectory. | |
636 | ||
637 | For debugging purposes the control is given to the user at each | |
638 | tracking step via the routine EUSTEP. This is the equivalent of | |
639 | the GEANT GUSTEP routine. The routine ERXYZC is a copy of GPCXYZ. | |
640 | ||
641 | After the initialization a call to ERTRAK can be done at any | |
642 | place of the user's code both in a context of Reconstruction or | |
643 | Simulation, without interfering with other eventual GEANT calls | |
644 | in the same event. | |
645 | ||
646 | ||
647 | ||
648 | Flow-Chart for subroutine ERTRAK | |
649 | ================================ | |
650 | ||
651 | ||
652 | ERTRAK | |
653 | ====== | |
654 | | | |
655 | | ----> GUFLD | |
656 | | ----> TRSCSD | |
657 | | ----> TRSDSC | |
658 | | ----> ERBCER | |
659 | | ----> GEKBIN | |
660 | | ----> ERPINI ----> TRPROP | |
661 | | ----> ERTRGO | |
662 | ====== | |
663 | | | |
664 | | ----> GMEDIA | |
665 | | ----> GUFLD | |
666 | | ----> EVOLIO | |
667 | | ----> ERSTOR | |
668 | | ====== | |
669 | | | | |
670 | | | ----> ERBCER | |
671 | | | ----> ERBCTR | |
672 | | | ----> TRSCSD | |
673 | | | ----> DMMSS | |
674 | | | |
675 | | | |
676 | | ----> ERTRCH -| ----> GTNEXT | |
677 | | ====== | | |
678 | | ----> ERTRNT -| ----> GUSWIM ----> GUFLD | |
679 | | ====== | ----> GINVOL | |
680 | | | ----> ERLAND | |
681 | | | ----> GEKBIN | |
682 | | | ----> ERPROP | |
683 | | | ====== | |
684 | | ----> EUSTEP | | | |
685 | | ----> EVOLIO | | ----> GUFLD | |
686 | l | | ----> TRPROP | |
687 | l ----> GTMEDI | | ----> TRPRFN | |
688 | | | ----> SSMTST | |
689 | | | ----> ERMCSC | |
690 | | | |
691 | | ----> ERSTOR | |
692 | ====== | |
693 | ||
694 | ||
695 | ||
696 | 6. Acknowledgements | |
697 | ================ | |
698 | ||
699 | The authors of the present interface benefitted numerous critical | |
700 | remarks and useful suggestions from the authors of the GEANT3 Package, | |
701 | especially from F. Bruyant (CERN), which are greatly acknowledged here. | |
702 | ||
703 | Complaints and suggest must be sent to one of the authors : | |
704 | Innocent@cernvm Maire@cernvm Nagy@cernvm | |
705 | ||
706 | ||
707 | ||
708 | References: | |
709 | =========== | |
710 | ||
711 | ||
712 | [1] W.Wittek, EMC Internal Reports: EMC/80/15, EMCSW/80/39, | |
713 | EMCSW/81/13, EMCSW/81/18 | |
714 | A.Haas, The EMC Utility Package: UTIL42 | |
715 | ||
716 | [2] R.Brun, F.Bruyant, M.Maire, A.C.McPherson, P.Zanarini | |
717 | DD/EE/84-1, May 1986 | |
718 | ||
719 | ||
720 | ||
721 |