STEER/AliTrackPointArray.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //////////////////////////////////////////////////////////////////////////////
  17 //                          Class AliTrackPointArray                        //
  18 //   This class contains the ESD track space-points which are used during   //
  19 //   the alignment procedures. Each space-point consist of 3 coordinates    //
  20 //   (and their errors) and the index of the sub-detector which contains    //
  21 //   the space-point.                                                       //
  22 //   cvetan.cheshkov@cern.ch 3/11/2005                                      //
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25 #include <TMath.h>
  26 #include <TMatrixD.h>
  27 #include <TMatrixDSym.h>
  28
  29 #include "AliTrackPointArray.h"
  30
  31 ClassImp(AliTrackPointArray)
  32
  33 //______________________________________________________________________________
  34 AliTrackPointArray::AliTrackPointArray() :
  35   TObject(),
  36   fNPoints(0),
  37   fX(0),
  38   fY(0),
  39   fZ(0),
  40   fSize(0),
  41   fCov(0),
  42   fVolumeID(0)
  43 {
  44 }
  45
  46 //______________________________________________________________________________
  47 AliTrackPointArray::AliTrackPointArray(Int_t npoints):
  48   TObject(),
  49   fNPoints(npoints),
  50   fX(new Float_t[npoints]),
  51   fY(new Float_t[npoints]),
  52   fZ(new Float_t[npoints]),
  53   fSize(6*npoints),
  54   fCov(new Float_t[fSize]),
  55   fVolumeID(new UShort_t[npoints])
  56 {
  57   // Constructor
  58   //
  59 }
  60
  61 //______________________________________________________________________________
  62 AliTrackPointArray::AliTrackPointArray(const AliTrackPointArray &array):
  63   TObject(array),
  64   fNPoints(array.fNPoints),
  65   fX(new Float_t[fNPoints]),
  66   fY(new Float_t[fNPoints]),
  67   fZ(new Float_t[fNPoints]),
  68   fSize(array.fSize),
  69   fCov(new Float_t[fSize]),
  70   fVolumeID(new UShort_t[fNPoints])
  71 {
  72   // Copy constructor
  73   //
  74   memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
  75   memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
  76   memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
  77   memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
  78   memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
  79 }
  80
  81 //_____________________________________________________________________________
  82 AliTrackPointArray &AliTrackPointArray::operator =(const AliTrackPointArray& array)
  83 {
  84   // assignment operator
  85   //
  86   if(this==&array) return *this;
  87   ((TObject *)this)->operator=(array);
  88
  89   fNPoints = array.fNPoints;
  90   fSize = array.fSize;
  91   fX = new Float_t[fNPoints];
  92   fY = new Float_t[fNPoints];
  93   fZ = new Float_t[fNPoints];
  94   fVolumeID = new UShort_t[fNPoints];
  95   fCov = new Float_t[fSize];
  96   memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
  97   memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
  98   memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
  99   memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
 100   memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
 101
 102   return *this;
 103 }
 104
 105 //______________________________________________________________________________
 106 AliTrackPointArray::~AliTrackPointArray()
 107 {
 108   // Destructor
 109   //
 110   delete [] fX;
 111   delete [] fY;
 112   delete [] fZ;
 113   delete [] fVolumeID;
 114   delete [] fCov;
 115 }
 116
 117
 118 //______________________________________________________________________________
 119 Bool_t AliTrackPointArray::AddPoint(Int_t i, const AliTrackPoint *p)
 120 {
 121   // Add a point to the array at position i
 122   //
 123   if (i >= fNPoints) return kFALSE;
 124   fX[i] = p->GetX();
 125   fY[i] = p->GetY();
 126   fZ[i] = p->GetZ();
 127   fVolumeID[i] = p->GetVolumeID();
 128   memcpy(&fCov[6*i],p->GetCov(),6*sizeof(Float_t));
 129   return kTRUE;
 130 }
 131
 132 //______________________________________________________________________________
 133 Bool_t AliTrackPointArray::GetPoint(AliTrackPoint &p, Int_t i) const
 134 {
 135   // Get the point at position i
 136   //
 137   if (i >= fNPoints) return kFALSE;
 138   p.SetXYZ(fX[i],fY[i],fZ[i],&fCov[6*i]);
 139   p.SetVolumeID(fVolumeID[i]);
 140   return kTRUE;
 141 }
 142
 143 //______________________________________________________________________________
 144 Bool_t AliTrackPointArray::HasVolumeID(UShort_t volid) const
 145 {
 146   // This method checks if the array
 147   // has at least one hit in the detector
 148   // volume defined by volid
 149   Bool_t check = kFALSE;
 150   for (Int_t ipoint = 0; ipoint < fNPoints; ipoint++)
 151     if (fVolumeID[ipoint] == volid) check = kTRUE;
 152
 153   return check;
 154 }
 155
 156 ClassImp(AliTrackPoint)
 157
 158 //______________________________________________________________________________
 159 AliTrackPoint::AliTrackPoint() :
 160   TObject(),
 161   fX(0),
 162   fY(0),
 163   fZ(0),
 164   fVolumeID(0)
 165 {
 166   // Default constructor
 167   //
 168   memset(fCov,0,6*sizeof(Float_t));
 169 }
 170
 171
 172 //______________________________________________________________________________
 173 AliTrackPoint::AliTrackPoint(Float_t x, Float_t y, Float_t z, const Float_t *cov, UShort_t volid) :
 174   TObject(),
 175   fX(0),
 176   fY(0),
 177   fZ(0),
 178   fVolumeID(0)
 179 {
 180   // Constructor
 181   //
 182   SetXYZ(x,y,z,cov);
 183   SetVolumeID(volid);
 184 }
 185
 186 //______________________________________________________________________________
 187 AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid) :
 188   TObject(),
 189   fX(0),
 190   fY(0),
 191   fZ(0),
 192   fVolumeID(0)
 193 {
 194   // Constructor
 195   //
 196   SetXYZ(xyz[0],xyz[1],xyz[2],cov);
 197   SetVolumeID(volid);
 198 }
 199
 200 //______________________________________________________________________________
 201 AliTrackPoint::AliTrackPoint(const AliTrackPoint &p):
 202   TObject(p),
 203   fX(0),
 204   fY(0),
 205   fZ(0),
 206   fVolumeID(0)
 207 {
 208   // Copy constructor
 209   //
 210   SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
 211   SetVolumeID(p.fVolumeID);
 212 }
 213
 214 //_____________________________________________________________________________
 215 AliTrackPoint &AliTrackPoint::operator =(const AliTrackPoint& p)
 216 {
 217   // assignment operator
 218   //
 219   if(this==&p) return *this;
 220   ((TObject *)this)->operator=(p);
 221
 222   SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
 223   SetVolumeID(p.fVolumeID);
 224
 225   return *this;
 226 }
 227
 228 //______________________________________________________________________________
 229 void AliTrackPoint::SetXYZ(Float_t x, Float_t y, Float_t z, const Float_t *cov)
 230 {
 231   // Set XYZ coordinates and their cov matrix
 232   //
 233   fX = x;
 234   fY = y;
 235   fZ = z;
 236   if (cov)
 237     memcpy(fCov,cov,6*sizeof(Float_t));
 238 }
 239
 240 //______________________________________________________________________________
 241 void AliTrackPoint::SetXYZ(const Float_t *xyz, const Float_t *cov)
 242 {
 243   // Set XYZ coordinates and their cov matrix
 244   //
 245   SetXYZ(xyz[0],xyz[1],xyz[2],cov);
 246 }
 247
 248 //______________________________________________________________________________
 249 void AliTrackPoint::GetXYZ(Float_t *xyz, Float_t *cov) const
 250 {
 251   xyz[0] = fX;
 252   xyz[1] = fY;
 253   xyz[2] = fZ;
 254   if (cov)
 255     memcpy(cov,fCov,6*sizeof(Float_t));
 256 }
 257
 258 //______________________________________________________________________________
 259 Float_t AliTrackPoint::GetResidual(const AliTrackPoint &p, Bool_t weighted) const
 260 {
 261   // This method calculates the track to space-point residuals. The track
 262   // interpolation is also stored as AliTrackPoint. Using the option
 263   // 'weighted' one can calculate the residual either with or without
 264   // taking into account the covariance matrix of the space-point and
 265   // track interpolation. The second case the residual becomes a pull.
 266
 267   Float_t res = 0;
 268
 269   if (!weighted) {
 270     Float_t xyz[3],xyzp[3];
 271     GetXYZ(xyz);
 272     p.GetXYZ(xyzp);
 273     res = (xyz[0]-xyzp[0])*(xyz[0]-xyzp[0])+
 274           (xyz[1]-xyzp[1])*(xyz[1]-xyzp[1])+
 275           (xyz[2]-xyzp[2])*(xyz[2]-xyzp[2]);
 276   }
 277   else {
 278     Float_t xyz[3],xyzp[3];
 279     Float_t cov[6],covp[6];
 280     GetXYZ(xyz,cov);
 281     TMatrixDSym mcov(3);
 282     mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
 283     mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
 284     mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
 285     p.GetXYZ(xyzp,covp);
 286     TMatrixDSym mcovp(3);
 287     mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
 288     mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
 289     mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
 290     TMatrixDSym msum = mcov + mcovp;
 291     msum.Invert();
 292     //    mcov.Print(); mcovp.Print(); msum.Print();
 293     if (msum.IsValid()) {
 294       for (Int_t i = 0; i < 3; i++)
 295         for (Int_t j = 0; j < 3; j++)
 296           res += (xyz[i]-xyzp[i])*(xyz[j]-xyzp[j])*msum(i,j);
 297     }
 298   }
 299
 300   return res;
 301 }
 302
 303 //_____________________________________________________________________________
 304 Bool_t AliTrackPoint::GetPCA(const AliTrackPoint &p, AliTrackPoint &out) const
 305 {
 306   //
 307   // Get the intersection point between this point and
 308   // the point "p" belongs to.
 309   // The result is stored as a point 'out'
 310   // return kFALSE in case of failure.
 311   out.SetXYZ(0,0,0);
 312
 313   TMatrixD t(3,1);
 314   t(0,0)=GetX();
 315   t(1,0)=GetY();
 316   t(2,0)=GetZ();
 317
 318   TMatrixDSym tC(3);
 319   {
 320   const Float_t *cv=GetCov();
 321   tC(0,0)=cv[0]; tC(0,1)=cv[1]; tC(0,2)=cv[2];
 322   tC(1,0)=cv[1]; tC(1,1)=cv[3]; tC(1,2)=cv[4];
 323   tC(2,0)=cv[2]; tC(2,1)=cv[4]; tC(2,2)=cv[5];
 324   }
 325
 326   TMatrixD m(3,1);
 327   m(0,0)=p.GetX();
 328   m(1,0)=p.GetY();
 329   m(2,0)=p.GetZ();
 330
 331   TMatrixDSym mC(3);
 332   {
 333   const Float_t *cv=p.GetCov();
 334   mC(0,0)=cv[0]; mC(0,1)=cv[1]; mC(0,2)=cv[2];
 335   mC(1,0)=cv[1]; mC(1,1)=cv[3]; mC(1,2)=cv[4];
 336   mC(2,0)=cv[2]; mC(2,1)=cv[4]; mC(2,2)=cv[5];
 337   }
 338
 339   TMatrixDSym tmW(tC);
 340   tmW+=mC;
 341   tmW.Invert();
 342   if (!tmW.IsValid()) return kFALSE;
 343
 344   TMatrixD mW(tC,TMatrixD::kMult,tmW);
 345   TMatrixD tW(mC,TMatrixD::kMult,tmW);
 346
 347   TMatrixD mi(mW,TMatrixD::kMult,m);
 348   TMatrixD ti(tW,TMatrixD::kMult,t);
 349   ti+=mi;
 350
 351   TMatrixD iC(tC,TMatrixD::kMult,tmW);
 352   iC*=mC;
 353
 354   out.SetXYZ(ti(0,0),ti(1,0),ti(2,0));
 355   UShort_t id=p.GetVolumeID();
 356   out.SetVolumeID(id);
 357
 358   return kTRUE;
 359 }
 360
 361 //______________________________________________________________________________
 362 Float_t AliTrackPoint::GetAngle() const
 363 {
 364   // The method uses the covariance matrix of
 365   // the space-point in order to extract the
 366   // orientation of the detector plane.
 367   // The rotation in XY plane only is calculated.
 368
 369   Float_t phi= TMath::ATan2(TMath::Sqrt(fCov[0]),TMath::Sqrt(fCov[3]));
 370   if (fCov[1] > 0) {
 371     phi = TMath::Pi() - phi;
 372     if ((fY-fX) < 0) phi += TMath::Pi();
 373   }
 374   else {
 375     if ((fX+fY) < 0) phi += TMath::Pi();
 376   }
 377
 378   return phi;
 379
 380 }
 381
 382 //_____________________________________________________________________________
 383 AliTrackPoint& AliTrackPoint::Rotate(Float_t alpha) const
 384 {
 385   // Transform the space-point coordinates
 386   // and covariance matrix from global to
 387   // local (detector plane) coordinate system
 388   // XY plane rotation only
 389
 390   static AliTrackPoint p;
 391   p = *this;
 392
 393   Float_t xyz[3],cov[6];
 394   GetXYZ(xyz,cov);
 395
 396   Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
 397
 398   Float_t newxyz[3],newcov[6];
 399   newxyz[0] = cos*xyz[0] + sin*xyz[1];
 400   newxyz[1] = cos*xyz[1] - sin*xyz[0];
 401   newxyz[2] = xyz[2];
 402
 403   newcov[0] = cov[0]*cos*cos+
 404             2*cov[1]*sin*cos+
 405               cov[3]*sin*sin;
 406   newcov[1] = cov[1]*(cos*cos-sin*sin)+
 407              (cov[3]-cov[0])*sin*cos;
 408   newcov[2] = cov[2]*cos+
 409               cov[4]*sin;
 410   newcov[3] = cov[0]*sin*sin-
 411             2*cov[1]*sin*cos+
 412               cov[3]*cos*cos;
 413   newcov[4] = cov[4]*cos-
 414               cov[2]*sin;
 415   newcov[5] = cov[5];
 416
 417   p.SetXYZ(newxyz,newcov);
 418   p.SetVolumeID(GetVolumeID());
 419
 420   return p;
 421 }
 422
 423 //_____________________________________________________________________________
 424 AliTrackPoint& AliTrackPoint::MasterToLocal() const
 425 {
 426   // Transform the space-point coordinates
 427   // and the covariance matrix from the
 428   // (master) to the local (tracking)
 429   // coordinate system
 430
 431   Float_t alpha = GetAngle();
 432   return Rotate(alpha);
 433 }
 434
 435 //_____________________________________________________________________________
 436 void AliTrackPoint::Print(Option_t *) const
 437 {
 438   // Print the space-point coordinates and
 439   // covariance matrix
 440
 441   printf("VolumeID=%d\n", GetVolumeID());
 442   printf("X = %12.6f    Tx = %12.6f%12.6f%12.6f\n", fX, fCov[0], fCov[1], fCov[2]);
 443   printf("Y = %12.6f    Ty = %12.6f%12.6f%12.6f\n", fY, fCov[1], fCov[3], fCov[4]);
 444   printf("Z = %12.6f    Tz = %12.6f%12.6f%12.6f\n", fZ, fCov[2], fCov[4], fCov[5]);
 445
 446 }