Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : // $Id$
17 :
18 : ////////////////////////////////////////////////////////////////////////////////
19 : //
20 : // Author: Artur Szostak
21 : // Email: artur@alice.phy.uct.ac.za | artursz@iafrica.com
22 : //
23 : ////////////////////////////////////////////////////////////////////////////////
24 :
25 : #include "AliHLTMUONCalculations.h"
26 : #include "AliHLTMUONUtils.h"
27 : #include "AliHLTMUONTriggerRecordsBlockStruct.h"
28 : #include <cmath>
29 :
30 : AliHLTFloat32_t AliHLTMUONCalculations::fgZf = -975.0; // cm
31 :
32 : AliHLTFloat32_t AliHLTMUONCalculations::fgQBLScaled
33 : = 3.0 * 2.99792458e8 / 1e9; // T.m.*c/1e9
34 :
35 : AliHLTMUONParticleSign AliHLTMUONCalculations::fgSign = kSignUnknown;
36 : AliHLTFloat32_t AliHLTMUONCalculations::fgPx = 0; // GeV/c
37 : AliHLTFloat32_t AliHLTMUONCalculations::fgPy = 0; // GeV/c
38 : AliHLTFloat32_t AliHLTMUONCalculations::fgPz = 0; // GeV/c
39 :
40 : AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaX2 = 1.; // cm^2
41 : AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaY2 = 1.; // cm^2
42 :
43 : AliHLTFloat32_t AliHLTMUONCalculations::fgMzx = 0;
44 : AliHLTFloat32_t AliHLTMUONCalculations::fgMzy = 0;
45 : AliHLTFloat32_t AliHLTMUONCalculations::fgCzx = 0;
46 : AliHLTFloat32_t AliHLTMUONCalculations::fgCzy = 0;
47 :
48 : AliHLTFloat32_t AliHLTMUONCalculations::fgIdealX1 = 0; // cm
49 : AliHLTFloat32_t AliHLTMUONCalculations::fgIdealY1 = 0; // cm
50 : AliHLTFloat32_t AliHLTMUONCalculations::fgIdealZ1 = -1603.5f; // cm
51 : AliHLTFloat32_t AliHLTMUONCalculations::fgIdealX2 = 0; // cm
52 : AliHLTFloat32_t AliHLTMUONCalculations::fgIdealY2 = 0; // cm
53 : AliHLTFloat32_t AliHLTMUONCalculations::fgIdealZ2 = -1703.5f; // cm
54 :
55 :
56 : bool AliHLTMUONCalculations::ComputeMomentum(
57 : AliHLTFloat32_t x1,
58 : AliHLTFloat32_t y1, AliHLTFloat32_t y2,
59 : AliHLTFloat32_t z1, AliHLTFloat32_t z2
60 : )
61 : {
62 : /// Computes the momentum components based on the equations given in the
63 : /// ALICE dimuon spectrometer Technical Design Report (TDR-5): trigger section.
64 : ///
65 : /// Reference:
66 : /// "CERN/LHCC 2000-046
67 : /// Addendum 1 to ALICE TDR 5
68 : /// 15 Dec 2000"
69 : /// Section 3.1.2 pages 144 and 145.
70 : ///
71 : /// Input can be in meters, cm or mm. Output is in GeV/c.
72 : ///
73 : /// \param x1 X coordinate of hit point 1 on the track.
74 : /// \param y1 Y coordinate of hit point 1 on the track.
75 : /// \param z1 Z coordinate of hit point 1 on the track.
76 : /// \param y2 Y coordinate of hit point 2 on the track.
77 : /// \param z2 Z coordinate of hit point 2 on the track.
78 : /// \return true if the momentum could be calculated and false otherwise.
79 : /// If true is returned then the estimated momentum can be fetched by the
80 : /// method calls: Px(), Py() and Pz() for the individual components.
81 :
82 0 : AliHLTFloat64_t z2mz1 = z2 - z1;
83 0 : if (z2mz1 == 0 or z1 == 0)
84 : {
85 0 : fgSign = kSignUnknown;
86 0 : fgPx = fgPy = fgPz = 0;
87 0 : return false;
88 : }
89 0 : AliHLTFloat64_t thetaTimesZf = (y1*z2 - y2*z1) / z2mz1;
90 0 : AliHLTFloat64_t xf = x1 * fgZf / z1;
91 0 : AliHLTFloat64_t yf = y2 - ((y2-y1) * (z2-fgZf)) / z2mz1;
92 :
93 0 : if (thetaTimesZf == 0)
94 : {
95 0 : fgSign = kSignUnknown;
96 0 : fgPx = fgPy = fgPz = 0;
97 0 : return false;
98 : }
99 0 : AliHLTFloat64_t pDivZf = (fgQBLScaled / thetaTimesZf);
100 0 : AliHLTFloat64_t p = pDivZf * fgZf;
101 0 : pDivZf = fabs(pDivZf);
102 :
103 0 : if (p < 0)
104 0 : fgSign = kSignMinus;
105 0 : else if (p > 0)
106 0 : fgSign = kSignPlus;
107 : else
108 0 : fgSign = kSignUnknown;
109 :
110 0 : fgPx = AliHLTFloat32_t( pDivZf * xf );
111 0 : fgPy = AliHLTFloat32_t( pDivZf * yf );
112 0 : AliHLTFloat64_t k = p*p - fgPx*fgPx - fgPy*fgPy;
113 0 : if (k > 0)
114 0 : fgPz = AliHLTFloat32_t( sqrt(k) );
115 : else
116 0 : fgPz = 0;
117 : // fgPz must be the same sign as fgZf else it could not have been measured.
118 0 : if (fgZf < 0) fgPz = -fgPz;
119 :
120 : return true;
121 0 : }
122 :
123 :
124 : AliHLTFloat32_t AliHLTMUONCalculations::QBL()
125 : {
126 : // We have to convert back into Tesla metres.
127 0 : return fgQBLScaled * 1e9 / 2.99792458e8;
128 : }
129 :
130 :
131 : void AliHLTMUONCalculations::QBL(AliHLTFloat32_t value)
132 : {
133 : // Note: 2.99792458e8/1e9 is the conversion factor for GeV.
134 : // It is c/1e9, where c is the speed of light.
135 0 : fgQBLScaled = value * 2.99792458e8 / 1e9;
136 0 : }
137 :
138 :
139 : AliHLTFloat32_t AliHLTMUONCalculations::ComputeMass(
140 : AliHLTFloat32_t massA,
141 : AliHLTFloat32_t pxA,
142 : AliHLTFloat32_t pyA,
143 : AliHLTFloat32_t pzA,
144 : AliHLTFloat32_t massB,
145 : AliHLTFloat32_t pxB,
146 : AliHLTFloat32_t pyB,
147 : AliHLTFloat32_t pzB
148 : )
149 : {
150 : /// Calculates the invariant mass for a pair of particles.
151 : /// \param massA Mmass in GeV/c of particle A.
152 : /// \param pxA X component of the momentum in GeV/c for particle A.
153 : /// \param pyA Y component of the momentum in GeV/c for particle A.
154 : /// \param pzA Z component of the momentum in GeV/c for particle A.
155 : /// \param massB Mass in GeV/c of particle B.
156 : /// \param pxB X component of the momentum in GeV/c for particle B.
157 : /// \param pyB Y component of the momentum in GeV/c for particle B.
158 : /// \param pzB Z component of the momentum in GeV/c for particle B.
159 : /// \return The invariant mass in GeV/c^2 or -1 if there was a problem
160 : /// in the calculation due to bad input parameters.
161 :
162 0 : AliHLTFloat32_t massA2 = massA*massA;
163 0 : AliHLTFloat32_t massB2 = massB*massB;
164 0 : AliHLTFloat32_t energyA = sqrt(massA2 + pxA*pxA + pyA*pyA + pzA*pzA);
165 0 : AliHLTFloat32_t energyB = sqrt(massB2 + pxB*pxB + pyB*pyB + pzB*pzB);
166 0 : AliHLTFloat32_t mass2 = massA2 + massB2 + 2. * (energyA*energyB - pxA*pxB - pyA*pyB - pzA*pzB);
167 0 : if (mass2 < 0.) return -1.;
168 0 : return sqrt(mass2);
169 0 : }
170 :
171 :
172 : bool AliHLTMUONCalculations::FitLineToTriggerRecord(
173 : const AliHLTMUONTriggerRecordStruct& trigger
174 : )
175 : {
176 : /// Straight lines are fitted in the ZX and ZY planes using a least
177 : /// squares fit to the coordinates in the trigger record.
178 : /// http://mathworld.wolfram.com/LeastSquaresFitting.html
179 : /// If this method returns true, then the fitted parameters can fetched
180 : /// using the method calls Mzx(), Mzy(), Czx() and Czy(). The lines are
181 : /// then given by: x = Mzx() * z + Czx() and y = Mzy() * z + Czy()
182 : /// The ideal coordinates are also calculated and can be fetched with
183 : /// the method calls: IdealX1(), IdealY1() and IdealZ1() for point on MT1,
184 : /// and IdealX2(), IdealY2() and IdealZ2() for point on MT2.
185 : /// \param trigger The trigger record structure to which we fit a line.
186 : /// \return true if the line could be fitted or false otherwise.
187 : /// The reason for failure could be either too few hits or the slopes
188 : /// Mzx() or Mzy() would be infinite, implying a line that is
189 : /// perpendicular to the z axis.
190 :
191 0 : AliHLTMUONParticleSign sign;
192 0 : bool hitset[4];
193 0 : AliHLTMUONUtils::UnpackTriggerRecordFlags(trigger.fFlags, sign, hitset);
194 : DebugTrace("hitset = {" << hitset[0] << ", " << hitset[1] << ", "
195 : << hitset[2] << ", " << hitset[3] << "}"
196 : );
197 :
198 0 : return FitLineToTriggerRecord(trigger, hitset);
199 0 : }
200 :
201 :
202 : bool AliHLTMUONCalculations::FitLineToTriggerRecord(
203 : const AliHLTMUONTriggerRecordStruct& trigger,
204 : const bool hitset[4]
205 : )
206 : {
207 : /// Performs a straight line fit like FitLineToTriggerRecord(trigger)
208 : /// but requires pree-decoded flags indicating which hits were set.
209 : /// \param trigger The trigger record structure to which we fit a line.
210 : /// \param hitset Flags indicating which hits were set in the trigger record.
211 : /// \return true if the line could be fitted or false otherwise.
212 :
213 0 : bool lineOk = FitLine(trigger, hitset);
214 0 : if (lineOk)
215 : {
216 : // Calculate ideal points on chambers 11 and 13:
217 0 : fgIdealX1 = fgMzx * fgIdealZ1 + fgCzx;
218 0 : fgIdealY1 = fgMzy * fgIdealZ1 + fgCzy;
219 0 : fgIdealX2 = fgMzx * fgIdealZ2 + fgCzx;
220 0 : fgIdealY2 = fgMzy * fgIdealZ2 + fgCzy;
221 0 : }
222 0 : return lineOk;
223 : }
224 :
225 :
226 : bool AliHLTMUONCalculations::FitLine(
227 : const AliHLTMUONTriggerRecordStruct& trigger,
228 : const bool hitset[4]
229 : )
230 : {
231 : /// Performs a straight line fit to the trigger record hits which are indicated
232 : /// by the hitset flags array.
233 : /// \param trigger The trigger record structure to which we fit a line.
234 : /// \param hitset Flags indicating which hits to use and were set in the trigger record.
235 : /// \return true if the line could be fitted or false otherwise.
236 :
237 : AliHLTFloat32_t sumX = 0;
238 : AliHLTFloat32_t sumY = 0;
239 : AliHLTFloat32_t sumZ = 0;
240 : int n = 0;
241 0 : for (int i = 0; i < 4; i++)
242 : {
243 0 : if (hitset[i])
244 : {
245 0 : sumX += trigger.fHit[i].fX;
246 0 : sumY += trigger.fHit[i].fY;
247 0 : sumZ += trigger.fHit[i].fZ;
248 0 : n++;
249 0 : }
250 : }
251 0 : if (n < 2) return false;
252 0 : AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
253 0 : AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
254 0 : AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
255 :
256 : AliHLTFloat32_t vSSzz = 0;
257 : AliHLTFloat32_t vSSzx = 0;
258 : AliHLTFloat32_t vSSzy = 0;
259 0 : for (int i = 0; i < 4; i++)
260 : {
261 0 : if (hitset[i])
262 : {
263 0 : vSSzz += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fZ - meanZ);
264 0 : vSSzx += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fX - meanX);
265 0 : vSSzy += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fY - meanY);
266 0 : }
267 : }
268 :
269 : // Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
270 0 : if (vSSzz == 0) return false;
271 0 : fgMzx = vSSzx / vSSzz;
272 0 : fgMzy = vSSzy / vSSzz;
273 0 : fgCzx = meanX - fgMzx * meanZ;
274 0 : fgCzy = meanY - fgMzy * meanZ;
275 :
276 0 : return true;
277 0 : }
278 :
279 :
280 : bool AliHLTMUONCalculations::FitLineToData(
281 : const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
282 : const AliHLTFloat32_t* z, AliHLTUInt32_t n
283 : )
284 : {
285 : /// Straight lines are fitted in the ZX and ZY planes using a least
286 : /// squares fit for the (x, y, z) data points.
287 : /// http://mathworld.wolfram.com/LeastSquaresFitting.html
288 : /// If this method returns true, then the fitted parameters can fetched
289 : /// using the method calls Mzx(), Mzy(), Czx() and Czy(). The lines are
290 : /// then given by: x = Mzx() * z + Czx() and y = Mzy() * z + Czy()
291 : /// \param x This must point to the array of x data values.
292 : /// \param y This must point to the array of y data values.
293 : /// \param z This must point to the array of z data values.
294 : /// \param n Specifies the number of data points in the x, y and z arrays.
295 : /// \return true if the line could be fitted or false otherwise.
296 : /// The reason for failure could be either too few data points or the
297 : /// slopes Mzx() or Mzy() would be infinite, implying a line that is
298 : /// perpendicular to the z axis.
299 :
300 0 : if (n < 2) return false;
301 :
302 : AliHLTFloat32_t sumX = 0;
303 : AliHLTFloat32_t sumY = 0;
304 : AliHLTFloat32_t sumZ = 0;
305 0 : for (AliHLTUInt32_t i = 0; i < n; i++)
306 : {
307 0 : sumX += x[i];
308 0 : sumY += y[i];
309 0 : sumZ += z[i];
310 : }
311 0 : AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
312 0 : AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
313 0 : AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
314 :
315 : AliHLTFloat32_t vSSzz = 0;
316 : AliHLTFloat32_t vSSzx = 0;
317 : AliHLTFloat32_t vSSzy = 0;
318 0 : for (AliHLTUInt32_t i = 0; i < n; i++)
319 : {
320 0 : vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
321 0 : vSSzx += (z[i] - meanZ)*(x[i] - meanX);
322 0 : vSSzy += (z[i] - meanZ)*(y[i] - meanY);
323 : }
324 :
325 : // Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
326 0 : if (vSSzz == 0) return false;
327 0 : fgMzx = vSSzx / vSSzz;
328 0 : fgMzy = vSSzy / vSSzz;
329 0 : fgCzx = meanX - fgMzx * meanZ;
330 0 : fgCzy = meanY - fgMzy * meanZ;
331 :
332 0 : return true;
333 0 : }
334 :
335 :
336 : bool AliHLTMUONCalculations::FitLineToData(
337 : const AliHLTFloat32_t* x, const AliHLTFloat32_t* z, AliHLTUInt32_t n
338 : )
339 : {
340 : /// A straight line is fitted in the X, Z data points using a least squares fit.
341 : /// http://mathworld.wolfram.com/LeastSquaresFitting.html
342 : /// If this method returns true, then the fitted parameters can fetched using the
343 : /// method calls Mzx() and Czx(). The line is then given by: x = Mzx() * z + Czx()
344 : /// \param x This must point to the array of x data values.
345 : /// \param z This must point to the array of z data values.
346 : /// \param n Specifies the number of data points in the x and z arrays.
347 : /// \return true if the line could be fitted or false otherwise.
348 : /// The reason for failure could be either too few data points or the slopes
349 : /// Mzx() would be infinite, implying a line that is perpendicular to the z axis.
350 :
351 0 : if (n < 2) return false;
352 :
353 : AliHLTFloat32_t sumX = 0;
354 : AliHLTFloat32_t sumZ = 0;
355 0 : for (AliHLTUInt32_t i = 0; i < n; i++)
356 : {
357 0 : sumX += x[i];
358 0 : sumZ += z[i];
359 : }
360 0 : AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
361 0 : AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
362 :
363 : AliHLTFloat32_t vSSzz = 0;
364 : AliHLTFloat32_t vSSzx = 0;
365 0 : for (AliHLTUInt32_t i = 0; i < n; i++)
366 : {
367 0 : vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
368 0 : vSSzx += (z[i] - meanZ)*(x[i] - meanX);
369 : }
370 :
371 : // Calculate params for line x = fgMzx * z + fgCzx.
372 0 : if (vSSzz == 0) return false;
373 0 : fgMzx = vSSzx / vSSzz;
374 0 : fgCzx = meanX - fgMzx * meanZ;
375 :
376 0 : return true;
377 0 : }
378 :
379 :
380 : AliHLTFloat32_t AliHLTMUONCalculations::AliHLTMUONCalculations::ComputeChi2(
381 : const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
382 : const AliHLTFloat32_t* z, AliHLTUInt32_t n
383 : )
384 : {
385 : /// Calculates the chi squared value for the set of data points given
386 : /// the fitted slope and coefficient parameters previously fitted by
387 : /// one of FitLine(x, y, z, n) or FitLineToTriggerRecord
388 : /// The fgSigmaX2 and fgSigmaY2 are used as the variance for the X and
389 : /// Y coordinates respectively. Note we assume that the covariance terms
390 : /// are zero.
391 : /// \param x This must point to the array of x data values.
392 : /// \param y This must point to the array of y data values.
393 : /// \param z This must point to the array of z data values.
394 : /// \param n Specifies the number of data points in the x, y and z arrays.
395 : /// \return The chi squared value.
396 :
397 : AliHLTFloat32_t chi2 = 0;
398 0 : for (AliHLTUInt32_t i = 0; i < n; i++)
399 : {
400 0 : AliHLTFloat32_t residualX = fgMzx * z[i] + fgCzx - x[i];
401 0 : AliHLTFloat32_t residualY = fgMzy * z[i] + fgCzy - y[i];
402 0 : chi2 += residualX*residualX/fgSigmaX2 + residualY*residualY/fgSigmaY2;
403 : }
404 0 : return chi2;
405 : }
406 :
407 :
408 : AliHLTFloat32_t AliHLTMUONCalculations::AliHLTMUONCalculations::ComputeChi2(
409 : const AliHLTMUONTriggerRecordStruct& trigger,
410 : const bool hitset[4]
411 : )
412 : {
413 : /// Calculates the chi squared value for trigger record using the hits
414 : /// indicated by the hitset array.
415 : /// \param trigger The trigger record structure for which we compute the chi squared value.
416 : /// \param hitset Flags indicating which hits to use and were set in the trigger record.
417 : /// \return The chi squared value or -1 if it could not be calculated.
418 :
419 0 : if (not FitLine(trigger, hitset)) return -1;
420 : AliHLTFloat32_t chi2 = 0;
421 0 : for (AliHLTUInt32_t i = 0; i < 4; i++)
422 : {
423 0 : if (hitset[i])
424 : {
425 0 : AliHLTFloat32_t residualX = fgMzx * trigger.fHit[i].fZ + fgCzx - trigger.fHit[i].fX;
426 0 : AliHLTFloat32_t residualY = fgMzy * trigger.fHit[i].fZ + fgCzy - trigger.fHit[i].fY;
427 0 : chi2 += residualX*residualX/fgSigmaX2 + residualY*residualY/fgSigmaY2;
428 0 : }
429 : }
430 : return chi2;
431 0 : }
|
