1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////////
20 // Author: Artur Szostak
21 // Email: artur@alice.phy.uct.ac.za | artursz@iafrica.com
23 ////////////////////////////////////////////////////////////////////////////////
25 #include "AliHLTMUONCalculations.h"
26 #include "AliHLTMUONUtils.h"
27 #include "AliHLTMUONTriggerRecordsBlockStruct.h"
30 AliHLTFloat32_t AliHLTMUONCalculations::fgZf = -975.0; // cm
32 AliHLTFloat32_t AliHLTMUONCalculations::fgQBLScaled
33 = 3.0 * 2.99792458e8 / 1e9; // T.m.*c/1e9
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
40 AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaX2 = 1.; // cm^2
41 AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaY2 = 1.; // cm^2
43 AliHLTFloat32_t AliHLTMUONCalculations::fgMzx = 0;
44 AliHLTFloat32_t AliHLTMUONCalculations::fgMzy = 0;
45 AliHLTFloat32_t AliHLTMUONCalculations::fgCzx = 0;
46 AliHLTFloat32_t AliHLTMUONCalculations::fgCzy = 0;
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
56 bool AliHLTMUONCalculations::ComputeMomentum(
58 AliHLTFloat32_t y1, AliHLTFloat32_t y2,
59 AliHLTFloat32_t z1, AliHLTFloat32_t z2
62 /// Computes the momentum components based on the equations given in the
63 /// ALICE dimuon spectrometer Technical Design Report (TDR-5): trigger section.
66 /// "CERN/LHCC 2000-046
67 /// Addendum 1 to ALICE TDR 5
69 /// Section 3.1.2 pages 144 and 145.
71 /// Input can be in meters, cm or mm. Output is in GeV/c.
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.
82 AliHLTFloat64_t z2mz1 = z2 - z1;
83 if (z2mz1 == 0 or z1 == 0)
85 fgSign = kSignUnknown;
86 fgPx = fgPy = fgPz = 0;
89 AliHLTFloat64_t thetaTimesZf = (y1*z2 - y2*z1) / z2mz1;
90 AliHLTFloat64_t xf = x1 * fgZf / z1;
91 AliHLTFloat64_t yf = y2 - ((y2-y1) * (z2-fgZf)) / z2mz1;
93 if (thetaTimesZf == 0)
95 fgSign = kSignUnknown;
96 fgPx = fgPy = fgPz = 0;
99 AliHLTFloat64_t pDivZf = (fgQBLScaled / thetaTimesZf);
100 AliHLTFloat64_t p = pDivZf * fgZf;
101 pDivZf = fabs(pDivZf);
108 fgSign = kSignUnknown;
110 fgPx = AliHLTFloat32_t( pDivZf * xf );
111 fgPy = AliHLTFloat32_t( pDivZf * yf );
112 AliHLTFloat64_t k = p*p - fgPx*fgPx - fgPy*fgPy;
114 fgPz = AliHLTFloat32_t( sqrt(k) );
117 // fgPz must be the same sign as fgZf else it could not have been measured.
118 if (fgZf < 0) fgPz = -fgPz;
124 AliHLTFloat32_t AliHLTMUONCalculations::QBL()
126 // We have to convert back into Tesla metres.
127 return fgQBLScaled * 1e9 / 2.99792458e8;
131 void AliHLTMUONCalculations::QBL(AliHLTFloat32_t value)
133 // Note: 2.99792458e8/1e9 is the conversion factor for GeV.
134 // It is c/1e9, where c is the speed of light.
135 fgQBLScaled = value * 2.99792458e8 / 1e9;
139 AliHLTFloat32_t AliHLTMUONCalculations::ComputeMass(
140 AliHLTFloat32_t massA,
144 AliHLTFloat32_t massB,
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.
162 AliHLTFloat32_t massA2 = massA*massA;
163 AliHLTFloat32_t massB2 = massB*massB;
164 AliHLTFloat32_t energyA = sqrt(massA2 + pxA*pxA + pyA*pyA + pzA*pzA);
165 AliHLTFloat32_t energyB = sqrt(massB2 + pxB*pxB + pyB*pyB + pzB*pzB);
166 AliHLTFloat32_t mass2 = massA2 + massB2 + 2. * (energyA*energyB - pxA*pxB - pyA*pyB - pzA*pzB);
167 if (mass2 < 0.) return -1.;
172 bool AliHLTMUONCalculations::FitLineToTriggerRecord(
173 const AliHLTMUONTriggerRecordStruct& trigger
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.
191 AliHLTMUONParticleSign sign;
193 AliHLTMUONUtils::UnpackTriggerRecordFlags(trigger.fFlags, sign, hitset);
194 DebugTrace("hitset = {" << hitset[0] << ", " << hitset[1] << ", "
195 << hitset[2] << ", " << hitset[3] << "}"
198 return FitLineToTriggerRecord(trigger, hitset);
202 bool AliHLTMUONCalculations::FitLineToTriggerRecord(
203 const AliHLTMUONTriggerRecordStruct& trigger,
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.
213 AliHLTFloat32_t sumX = 0;
214 AliHLTFloat32_t sumY = 0;
215 AliHLTFloat32_t sumZ = 0;
217 for (int i = 0; i < 4; i++)
221 sumX += trigger.fHit[i].fX;
222 sumY += trigger.fHit[i].fY;
223 sumZ += trigger.fHit[i].fZ;
227 if (n < 2) return false;
228 AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
229 AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
230 AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
232 AliHLTFloat32_t vSSzz = 0;
233 AliHLTFloat32_t vSSzx = 0;
234 AliHLTFloat32_t vSSzy = 0;
235 for (int i = 0; i < 4; i++)
239 vSSzz += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fZ - meanZ);
240 vSSzx += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fX - meanX);
241 vSSzy += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fY - meanY);
245 // Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
246 if (vSSzz == 0) return false;
247 fgMzx = vSSzx / vSSzz;
248 fgMzy = vSSzy / vSSzz;
249 fgCzx = meanX - fgMzx * meanZ;
250 fgCzy = meanY - fgMzy * meanZ;
252 // Calculate ideal points on chambers 11 and 13:
253 fgIdealX1 = fgMzx * fgIdealZ1 + fgCzx;
254 fgIdealY1 = fgMzy * fgIdealZ1 + fgCzy;
255 fgIdealX2 = fgMzx * fgIdealZ2 + fgCzx;
256 fgIdealY2 = fgMzy * fgIdealZ2 + fgCzy;
262 bool AliHLTMUONCalculations::FitLineToData(
263 const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
264 const AliHLTFloat32_t* z, AliHLTUInt32_t n
267 /// Straight lines are fitted in the ZX and ZY planes using a least
268 /// squares fit for the (x, y, z) data points.
269 /// http://mathworld.wolfram.com/LeastSquaresFitting.html
270 /// If this method returns true, then the fitted parameters can fetched
271 /// using the method calls Mzx(), Mzy(), Czx() and Czy(). The lines are
272 /// then given by: x = Mzx() * z + Czx() and y = Mzy() * z + Czy()
273 /// \param x This must point to the array of x data values.
274 /// \param y This must point to the array of y data values.
275 /// \param z This must point to the array of z data values.
276 /// \param n Specifies the number of data points in the x, y and z arrays.
277 /// \return true if the line could be fitted or false otherwise.
278 /// The reason for failure could be either too few data points or the
279 /// slopes Mzx() or Mzy() would be infinite, implying a line that is
280 /// perpendicular to the z axis.
282 if (n < 2) return false;
284 AliHLTFloat32_t sumX = 0;
285 AliHLTFloat32_t sumY = 0;
286 AliHLTFloat32_t sumZ = 0;
287 for (AliHLTUInt32_t i = 0; i < n; i++)
293 AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
294 AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
295 AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
297 AliHLTFloat32_t vSSzz = 0;
298 AliHLTFloat32_t vSSzx = 0;
299 AliHLTFloat32_t vSSzy = 0;
300 for (AliHLTUInt32_t i = 0; i < n; i++)
302 vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
303 vSSzx += (z[i] - meanZ)*(x[i] - meanX);
304 vSSzy += (z[i] - meanZ)*(y[i] - meanY);
307 // Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
308 if (vSSzz == 0) return false;
309 fgMzx = vSSzx / vSSzz;
310 fgMzy = vSSzy / vSSzz;
311 fgCzx = meanX - fgMzx * meanZ;
312 fgCzy = meanY - fgMzy * meanZ;
318 bool AliHLTMUONCalculations::FitLineToData(
319 const AliHLTFloat32_t* x, const AliHLTFloat32_t* z, AliHLTUInt32_t n
322 /// A straight line is fitted in the X, Z data points using a least squares fit.
323 /// http://mathworld.wolfram.com/LeastSquaresFitting.html
324 /// If this method returns true, then the fitted parameters can fetched using the
325 /// method calls Mzx() and Czx(). The line is then given by: x = Mzx() * z + Czx()
326 /// \param x This must point to the array of x data values.
327 /// \param z This must point to the array of z data values.
328 /// \param n Specifies the number of data points in the x and z arrays.
329 /// \return true if the line could be fitted or false otherwise.
330 /// The reason for failure could be either too few data points or the slopes
331 /// Mzx() would be infinite, implying a line that is perpendicular to the z axis.
333 if (n < 2) return false;
335 AliHLTFloat32_t sumX = 0;
336 AliHLTFloat32_t sumZ = 0;
337 for (AliHLTUInt32_t i = 0; i < n; i++)
342 AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
343 AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
345 AliHLTFloat32_t vSSzz = 0;
346 AliHLTFloat32_t vSSzx = 0;
347 for (AliHLTUInt32_t i = 0; i < n; i++)
349 vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
350 vSSzx += (z[i] - meanZ)*(x[i] - meanX);
353 // Calculate params for line x = fgMzx * z + fgCzx.
354 if (vSSzz == 0) return false;
355 fgMzx = vSSzx / vSSzz;
356 fgCzx = meanX - fgMzx * meanZ;
362 AliHLTFloat32_t AliHLTMUONCalculations::AliHLTMUONCalculations::ComputeChi2(
363 const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
364 const AliHLTFloat32_t* z, AliHLTUInt32_t n
367 /// Calculates the chi squared value for the set of data points given
368 /// the fitted slope and coefficient parameters previously fitted by
369 /// LineFit(x, y, z, n);
370 /// The fgSigmaX2 and fgSigmaY2 are used as the variance for the X and
371 /// Y coordinates respectively. Note we assume that the covariance terms
373 /// \param x This must point to the array of x data values.
374 /// \param y This must point to the array of y data values.
375 /// \param z This must point to the array of z data values.
376 /// \param n Specifies the number of data points in the x, y and z arrays.
377 /// \return The chi squared value.
379 AliHLTFloat32_t chi2 = 0;
380 for (AliHLTUInt32_t i = 0; i < n; i++)
382 AliHLTFloat32_t residualX = fgMzx * z[i] + fgCzx - x[i];
383 AliHLTFloat32_t residualY = fgMzy * z[i] + fgCzy - y[i];
384 chi2 += residualX*residualX/fgSigmaX2 + residualY*residualY/fgSigmaY2;