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         AliHLTFloat32_t sumX = 0;
 214         AliHLTFloat32_t sumY = 0;
 215         AliHLTFloat32_t sumZ = 0;
 216         int n = 0;
 217         for (int i = 0; i < 4; i++)
 218         {
 219                 if (hitset[i])
 220                 {
 221                         sumX += trigger.fHit[i].fX;
 222                         sumY += trigger.fHit[i].fY;
 223                         sumZ += trigger.fHit[i].fZ;
 224                         n++;
 225                 }
 226         }
 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);
 231
 232         AliHLTFloat32_t vSSzz = 0;
 233         AliHLTFloat32_t vSSzx = 0;
 234         AliHLTFloat32_t vSSzy = 0;
 235         for (int i = 0; i < 4; i++)
 236         {
 237                 if (hitset[i])
 238                 {
 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);
 242                 }
 243         }
 244
 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;
 251
 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;
 257
 258         return true;
 259 }
 260
 261
 262 bool AliHLTMUONCalculations::FitLineToData(
 263                 const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
 264                 const AliHLTFloat32_t* z, AliHLTUInt32_t n
 265         )
 266 {
 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.
 281
 282         if (n < 2) return false;
 283
 284         AliHLTFloat32_t sumX = 0;
 285         AliHLTFloat32_t sumY = 0;
 286         AliHLTFloat32_t sumZ = 0;
 287         for (AliHLTUInt32_t i = 0; i < n; i++)
 288         {
 289                 sumX += x[i];
 290                 sumY += y[i];
 291                 sumZ += z[i];
 292         }
 293         AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
 294         AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
 295         AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
 296
 297         AliHLTFloat32_t vSSzz = 0;
 298         AliHLTFloat32_t vSSzx = 0;
 299         AliHLTFloat32_t vSSzy = 0;
 300         for (AliHLTUInt32_t i = 0; i < n; i++)
 301         {
 302                 vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
 303                 vSSzx += (z[i] - meanZ)*(x[i] - meanX);
 304                 vSSzy += (z[i] - meanZ)*(y[i] - meanY);
 305         }
 306
 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;
 313
 314         return true;
 315 }
 316
 317
 318 AliHLTFloat32_t AliHLTMUONCalculations::AliHLTMUONCalculations::ComputeChi2(
 319                 const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
 320                 const AliHLTFloat32_t* z, AliHLTUInt32_t n
 321         )
 322 {
 323         /// Calculates the chi squared value for the set of data points given
 324         /// the fitted slope and coefficient parameters previously fitted by
 325         /// LineFit(x, y, z, n);
 326         /// The fgSigmaX2 and fgSigmaY2 are used as the variance for the X and
 327         /// Y coordinates respectively. Note we assume that the covariance terms
 328         /// are zero.
 329         /// \param x  This must point to the array of x data values.
 330         /// \param y  This must point to the array of y data values.
 331         /// \param z  This must point to the array of z data values.
 332         /// \param n  Specifies the number of data points in the x, y and z arrays.
 333         /// \return  The chi squared value.
 334
 335         AliHLTFloat32_t chi2 = 0;
 336         for (AliHLTUInt32_t i = 0; i < n; i++)
 337         {
 338                 AliHLTFloat32_t residualX = fgMzx * z[i] + fgCzx - x[i];
 339                 AliHLTFloat32_t residualY = fgMzy * z[i] + fgCzy - y[i];
 340                 chi2 += residualX*residualX/fgSigmaX2 + residualY*residualY/fgSigmaY2;
 341         }
 342         return chi2;
 343 }