HLT/MUON/OnlineAnalysis/AliHLTMUONCalculations.cxx

   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         AliHLTFloat64_t z2mz1 = z2 - z1;
  83         if (z2mz1 == 0 or z1 == 0)
  84         {
  85                 fgSign = kSignUnknown;
  86                 fgPx = fgPy = fgPz = 0;
  87                 return false;
  88         }
  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;
  92
  93         if (thetaTimesZf == 0)
  94         {
  95                 fgSign = kSignUnknown;
  96                 fgPx = fgPy = fgPz = 0;
  97                 return false;
  98         }
  99         AliHLTFloat64_t pDivZf = (fgQBLScaled / thetaTimesZf);
 100         AliHLTFloat64_t p = pDivZf * fgZf;
 101         pDivZf = fabs(pDivZf);
 102
 103         if (p < 0)
 104                 fgSign = kSignMinus;
 105         else if (p > 0)
 106                 fgSign = kSignPlus;
 107         else
 108                 fgSign = kSignUnknown;
 109
 110         fgPx = AliHLTFloat32_t( pDivZf * xf );
 111         fgPy = AliHLTFloat32_t( pDivZf * yf );
 112         AliHLTFloat64_t k = p*p - fgPx*fgPx - fgPy*fgPy;
 113         if (k > 0)
 114                 fgPz = AliHLTFloat32_t( sqrt(k) );
 115         else
 116                 fgPz = 0;
 117         // fgPz must be the same sign as fgZf else it could not have been measured.
 118         if (fgZf < 0) fgPz = -fgPz;
 119
 120         return true;
 121 }
 122
 123
 124 AliHLTFloat32_t AliHLTMUONCalculations::QBL()
 125 {
 126         // We have to convert back into Tesla metres.
 127         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         fgQBLScaled = value * 2.99792458e8 / 1e9;
 136 }
 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         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.;
 168         return sqrt(mass2);
 169 }
 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         AliHLTMUONParticleSign sign;
 192         bool hitset[4];
 193         AliHLTMUONUtils::UnpackTriggerRecordFlags(trigger.fFlags, sign, hitset);
 194         DebugTrace("hitset = {" << hitset[0] << ", " << hitset[1] << ", "
 195                 << hitset[2] << ", " << hitset[3] << "}"
 196         );
 197
 198         return FitLineToTriggerRecord(trigger, hitset);
 199 }
 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         bool lineOk = FitLine(trigger, hitset);
 214         if (lineOk)
 215         {
 216                 // Calculate ideal points on chambers 11 and 13:
 217                 fgIdealX1 = fgMzx * fgIdealZ1 + fgCzx;
 218                 fgIdealY1 = fgMzy * fgIdealZ1 + fgCzy;
 219                 fgIdealX2 = fgMzx * fgIdealZ2 + fgCzx;
 220                 fgIdealY2 = fgMzy * fgIdealZ2 + fgCzy;
 221         }
 222         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         for (int i = 0; i < 4; i++)
 242         {
 243                 if (hitset[i])
 244                 {
 245                         sumX += trigger.fHit[i].fX;
 246                         sumY += trigger.fHit[i].fY;
 247                         sumZ += trigger.fHit[i].fZ;
 248                         n++;
 249                 }
 250         }
 251         if (n < 2) return false;
 252         AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
 253         AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
 254         AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
 255
 256         AliHLTFloat32_t vSSzz = 0;
 257         AliHLTFloat32_t vSSzx = 0;
 258         AliHLTFloat32_t vSSzy = 0;
 259         for (int i = 0; i < 4; i++)
 260         {
 261                 if (hitset[i])
 262                 {
 263                         vSSzz += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fZ - meanZ);
 264                         vSSzx += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fX - meanX);
 265                         vSSzy += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fY - meanY);
 266                 }
 267         }
 268
 269         // Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
 270         if (vSSzz == 0) return false;
 271         fgMzx = vSSzx / vSSzz;
 272         fgMzy = vSSzy / vSSzz;
 273         fgCzx = meanX - fgMzx * meanZ;
 274         fgCzy = meanY - fgMzy * meanZ;
 275
 276         return true;
 277 }
 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         if (n < 2) return false;
 301
 302         AliHLTFloat32_t sumX = 0;
 303         AliHLTFloat32_t sumY = 0;
 304         AliHLTFloat32_t sumZ = 0;
 305         for (AliHLTUInt32_t i = 0; i < n; i++)
 306         {
 307                 sumX += x[i];
 308                 sumY += y[i];
 309                 sumZ += z[i];
 310         }
 311         AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
 312         AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
 313         AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
 314
 315         AliHLTFloat32_t vSSzz = 0;
 316         AliHLTFloat32_t vSSzx = 0;
 317         AliHLTFloat32_t vSSzy = 0;
 318         for (AliHLTUInt32_t i = 0; i < n; i++)
 319         {
 320                 vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
 321                 vSSzx += (z[i] - meanZ)*(x[i] - meanX);
 322                 vSSzy += (z[i] - meanZ)*(y[i] - meanY);
 323         }
 324
 325         // Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
 326         if (vSSzz == 0) return false;
 327         fgMzx = vSSzx / vSSzz;
 328         fgMzy = vSSzy / vSSzz;
 329         fgCzx = meanX - fgMzx * meanZ;
 330         fgCzy = meanY - fgMzy * meanZ;
 331
 332         return true;
 333 }
 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         if (n < 2) return false;
 352
 353         AliHLTFloat32_t sumX = 0;
 354         AliHLTFloat32_t sumZ = 0;
 355         for (AliHLTUInt32_t i = 0; i < n; i++)
 356         {
 357                 sumX += x[i];
 358                 sumZ += z[i];
 359         }
 360         AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
 361         AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
 362
 363         AliHLTFloat32_t vSSzz = 0;
 364         AliHLTFloat32_t vSSzx = 0;
 365         for (AliHLTUInt32_t i = 0; i < n; i++)
 366         {
 367                 vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
 368                 vSSzx += (z[i] - meanZ)*(x[i] - meanX);
 369         }
 370
 371         // Calculate params for line x = fgMzx * z + fgCzx.
 372         if (vSSzz == 0) return false;
 373         fgMzx = vSSzx / vSSzz;
 374         fgCzx = meanX - fgMzx * meanZ;
 375
 376         return true;
 377 }
 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         for (AliHLTUInt32_t i = 0; i < n; i++)
 399         {
 400                 AliHLTFloat32_t residualX = fgMzx * z[i] + fgCzx - x[i];
 401                 AliHLTFloat32_t residualY = fgMzy * z[i] + fgCzy - y[i];
 402                 chi2 += residualX*residualX/fgSigmaX2 + residualY*residualY/fgSigmaY2;
 403         }
 404         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         if (not FitLine(trigger, hitset)) return -1;
 420         AliHLTFloat32_t chi2 = 0;
 421         for (AliHLTUInt32_t i = 0; i < 4; i++)
 422         {
 423                 if (hitset[i])
 424                 {
 425                         AliHLTFloat32_t residualX = fgMzx * trigger.fHit[i].fZ + fgCzx - trigger.fHit[i].fX;
 426                         AliHLTFloat32_t residualY = fgMzy * trigger.fHit[i].fZ + fgCzy - trigger.fHit[i].fY;
 427                         chi2 += residualX*residualX/fgSigmaX2 + residualY*residualY/fgSigmaY2;
 428                 }
 429         }
 430         return chi2;
 431 }