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   delete [] fX;
  92   fX = new Float_t[fNPoints];
  93   delete [] fY;
  94   fY = new Float_t[fNPoints];
  95   delete [] fZ;
  96   fZ = new Float_t[fNPoints];
  97   delete [] fVolumeID;
  98   fVolumeID = new UShort_t[fNPoints];
  99   delete [] fCov;
 100   fCov = new Float_t[fSize];
 101   memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
 102   memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
 103   memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
 104   memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
 105   memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
 106
 107   return *this;
 108 }
 109
 110 //______________________________________________________________________________
 111 AliTrackPointArray::~AliTrackPointArray()
 112 {
 113   // Destructor
 114   //
 115   delete [] fX;
 116   delete [] fY;
 117   delete [] fZ;
 118   delete [] fVolumeID;
 119   delete [] fCov;
 120 }
 121
 122
 123 //______________________________________________________________________________
 124 Bool_t AliTrackPointArray::AddPoint(Int_t i, const AliTrackPoint *p)
 125 {
 126   // Add a point to the array at position i
 127   //
 128   if (i >= fNPoints) return kFALSE;
 129   fX[i] = p->GetX();
 130   fY[i] = p->GetY();
 131   fZ[i] = p->GetZ();
 132   fVolumeID[i] = p->GetVolumeID();
 133   memcpy(&fCov[6*i],p->GetCov(),6*sizeof(Float_t));
 134   return kTRUE;
 135 }
 136
 137
 138 //______________________________________________________________________________
 139 Bool_t AliTrackPointArray::GetPoint(AliTrackPoint &p, Int_t i) const
 140 {
 141   // Get the point at position i
 142   //
 143   if (i >= fNPoints) return kFALSE;
 144   p.SetXYZ(fX[i],fY[i],fZ[i],&fCov[6*i]);
 145   p.SetVolumeID(fVolumeID[i]);
 146   return kTRUE;
 147 }
 148
 149 //______________________________________________________________________________
 150 Bool_t AliTrackPointArray::HasVolumeID(UShort_t volid) const
 151 {
 152   // This method checks if the array
 153   // has at least one hit in the detector
 154   // volume defined by volid
 155   Bool_t check = kFALSE;
 156   for (Int_t ipoint = 0; ipoint < fNPoints; ipoint++)
 157     if (fVolumeID[ipoint] == volid) check = kTRUE;
 158
 159   return check;
 160 }
 161
 162 ClassImp(AliTrackPoint)
 163
 164 //______________________________________________________________________________
 165 AliTrackPoint::AliTrackPoint() :
 166   TObject(),
 167   fX(0),
 168   fY(0),
 169   fZ(0),
 170   fVolumeID(0)
 171 {
 172   // Default constructor
 173   //
 174   memset(fCov,0,6*sizeof(Float_t));
 175 }
 176
 177
 178 //______________________________________________________________________________
 179 AliTrackPoint::AliTrackPoint(Float_t x, Float_t y, Float_t z, const Float_t *cov, UShort_t volid) :
 180   TObject(),
 181   fX(0),
 182   fY(0),
 183   fZ(0),
 184   fVolumeID(0)
 185 {
 186   // Constructor
 187   //
 188   SetXYZ(x,y,z,cov);
 189   SetVolumeID(volid);
 190 }
 191
 192 //______________________________________________________________________________
 193 AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid) :
 194   TObject(),
 195   fX(0),
 196   fY(0),
 197   fZ(0),
 198   fVolumeID(0)
 199 {
 200   // Constructor
 201   //
 202   SetXYZ(xyz[0],xyz[1],xyz[2],cov);
 203   SetVolumeID(volid);
 204 }
 205
 206 //______________________________________________________________________________
 207 AliTrackPoint::AliTrackPoint(const AliTrackPoint &p):
 208   TObject(p),
 209   fX(0),
 210   fY(0),
 211   fZ(0),
 212   fVolumeID(0)
 213 {
 214   // Copy constructor
 215   //
 216   SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
 217   SetVolumeID(p.fVolumeID);
 218 }
 219
 220 //_____________________________________________________________________________
 221 AliTrackPoint &AliTrackPoint::operator =(const AliTrackPoint& p)
 222 {
 223   // assignment operator
 224   //
 225   if(this==&p) return *this;
 226   ((TObject *)this)->operator=(p);
 227
 228   SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
 229   SetVolumeID(p.fVolumeID);
 230
 231   return *this;
 232 }
 233
 234 //______________________________________________________________________________
 235 void AliTrackPoint::SetXYZ(Float_t x, Float_t y, Float_t z, const Float_t *cov)
 236 {
 237   // Set XYZ coordinates and their cov matrix
 238   //
 239   fX = x;
 240   fY = y;
 241   fZ = z;
 242   if (cov)
 243     memcpy(fCov,cov,6*sizeof(Float_t));
 244 }
 245
 246 //______________________________________________________________________________
 247 void AliTrackPoint::SetXYZ(const Float_t *xyz, const Float_t *cov)
 248 {
 249   // Set XYZ coordinates and their cov matrix
 250   //
 251   SetXYZ(xyz[0],xyz[1],xyz[2],cov);
 252 }
 253
 254 //______________________________________________________________________________
 255 void AliTrackPoint::GetXYZ(Float_t *xyz, Float_t *cov) const
 256 {
 257   xyz[0] = fX;
 258   xyz[1] = fY;
 259   xyz[2] = fZ;
 260   if (cov)
 261     memcpy(cov,fCov,6*sizeof(Float_t));
 262 }
 263
 264 //______________________________________________________________________________
 265 Float_t AliTrackPoint::GetResidual(const AliTrackPoint &p, Bool_t weighted) const
 266 {
 267   // This method calculates the track to space-point residuals. The track
 268   // interpolation is also stored as AliTrackPoint. Using the option
 269   // 'weighted' one can calculate the residual either with or without
 270   // taking into account the covariance matrix of the space-point and
 271   // track interpolation. The second case the residual becomes a pull.
 272
 273   Float_t res = 0;
 274
 275   if (!weighted) {
 276     Float_t xyz[3],xyzp[3];
 277     GetXYZ(xyz);
 278     p.GetXYZ(xyzp);
 279     res = (xyz[0]-xyzp[0])*(xyz[0]-xyzp[0])+
 280           (xyz[1]-xyzp[1])*(xyz[1]-xyzp[1])+
 281           (xyz[2]-xyzp[2])*(xyz[2]-xyzp[2]);
 282   }
 283   else {
 284     Float_t xyz[3],xyzp[3];
 285     Float_t cov[6],covp[6];
 286     GetXYZ(xyz,cov);
 287     TMatrixDSym mcov(3);
 288     mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
 289     mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
 290     mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
 291     p.GetXYZ(xyzp,covp);
 292     TMatrixDSym mcovp(3);
 293     mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
 294     mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
 295     mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
 296     TMatrixDSym msum = mcov + mcovp;
 297     msum.Invert();
 298     //    mcov.Print(); mcovp.Print(); msum.Print();
 299     if (msum.IsValid()) {
 300       for (Int_t i = 0; i < 3; i++)
 301         for (Int_t j = 0; j < 3; j++)
 302           res += (xyz[i]-xyzp[i])*(xyz[j]-xyzp[j])*msum(i,j);
 303     }
 304   }
 305
 306   return res;
 307 }
 308
 309 //_____________________________________________________________________________
 310 Bool_t AliTrackPoint::GetPCA(const AliTrackPoint &p, AliTrackPoint &out) const
 311 {
 312   //
 313   // Get the intersection point between this point and
 314   // the point "p" belongs to.
 315   // The result is stored as a point 'out'
 316   // return kFALSE in case of failure.
 317   out.SetXYZ(0,0,0);
 318
 319   TMatrixD t(3,1);
 320   t(0,0)=GetX();
 321   t(1,0)=GetY();
 322   t(2,0)=GetZ();
 323
 324   TMatrixDSym tC(3);
 325   {
 326   const Float_t *cv=GetCov();
 327   tC(0,0)=cv[0]; tC(0,1)=cv[1]; tC(0,2)=cv[2];
 328   tC(1,0)=cv[1]; tC(1,1)=cv[3]; tC(1,2)=cv[4];
 329   tC(2,0)=cv[2]; tC(2,1)=cv[4]; tC(2,2)=cv[5];
 330   }
 331
 332   TMatrixD m(3,1);
 333   m(0,0)=p.GetX();
 334   m(1,0)=p.GetY();
 335   m(2,0)=p.GetZ();
 336
 337   TMatrixDSym mC(3);
 338   {
 339   const Float_t *cv=p.GetCov();
 340   mC(0,0)=cv[0]; mC(0,1)=cv[1]; mC(0,2)=cv[2];
 341   mC(1,0)=cv[1]; mC(1,1)=cv[3]; mC(1,2)=cv[4];
 342   mC(2,0)=cv[2]; mC(2,1)=cv[4]; mC(2,2)=cv[5];
 343   }
 344
 345   TMatrixDSym tmW(tC);
 346   tmW+=mC;
 347   tmW.Invert();
 348   if (!tmW.IsValid()) return kFALSE;
 349
 350   TMatrixD mW(tC,TMatrixD::kMult,tmW);
 351   TMatrixD tW(mC,TMatrixD::kMult,tmW);
 352
 353   TMatrixD mi(mW,TMatrixD::kMult,m);
 354   TMatrixD ti(tW,TMatrixD::kMult,t);
 355   ti+=mi;
 356
 357   TMatrixD iC(tC,TMatrixD::kMult,tmW);
 358   iC*=mC;
 359
 360   out.SetXYZ(ti(0,0),ti(1,0),ti(2,0));
 361   UShort_t id=p.GetVolumeID();
 362   out.SetVolumeID(id);
 363
 364   return kTRUE;
 365 }
 366
 367 //______________________________________________________________________________
 368 Float_t AliTrackPoint::GetAngle() const
 369 {
 370   // The method uses the covariance matrix of
 371   // the space-point in order to extract the
 372   // orientation of the detector plane.
 373   // The rotation in XY plane only is calculated.
 374
 375   Float_t phi= TMath::ATan2(TMath::Sqrt(fCov[0]),TMath::Sqrt(fCov[3]));
 376   if (fCov[1] > 0) {
 377     phi = TMath::Pi() - phi;
 378     if ((fY-fX) < 0) phi += TMath::Pi();
 379   }
 380   else {
 381     if ((fX+fY) < 0) phi += TMath::Pi();
 382   }
 383
 384   return phi;
 385
 386 }
 387
 388 //_____________________________________________________________________________
 389 AliTrackPoint& AliTrackPoint::Rotate(Float_t alpha) const
 390 {
 391   // Transform the space-point coordinates
 392   // and covariance matrix from global to
 393   // local (detector plane) coordinate system
 394   // XY plane rotation only
 395
 396   static AliTrackPoint p;
 397   p = *this;
 398
 399   Float_t xyz[3],cov[6];
 400   GetXYZ(xyz,cov);
 401
 402   Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
 403
 404   Float_t newxyz[3],newcov[6];
 405   newxyz[0] = cos*xyz[0] + sin*xyz[1];
 406   newxyz[1] = cos*xyz[1] - sin*xyz[0];
 407   newxyz[2] = xyz[2];
 408
 409   newcov[0] = cov[0]*cos*cos+
 410             2*cov[1]*sin*cos+
 411               cov[3]*sin*sin;
 412   newcov[1] = cov[1]*(cos*cos-sin*sin)+
 413              (cov[3]-cov[0])*sin*cos;
 414   newcov[2] = cov[2]*cos+
 415               cov[4]*sin;
 416   newcov[3] = cov[0]*sin*sin-
 417             2*cov[1]*sin*cos+
 418               cov[3]*cos*cos;
 419   newcov[4] = cov[4]*cos-
 420               cov[2]*sin;
 421   newcov[5] = cov[5];
 422
 423   p.SetXYZ(newxyz,newcov);
 424   p.SetVolumeID(GetVolumeID());
 425
 426   return p;
 427 }
 428
 429 //_____________________________________________________________________________
 430 AliTrackPoint& AliTrackPoint::MasterToLocal() const
 431 {
 432   // Transform the space-point coordinates
 433   // and the covariance matrix from the
 434   // (master) to the local (tracking)
 435   // coordinate system
 436
 437   Float_t alpha = GetAngle();
 438   return Rotate(alpha);
 439 }
 440
 441 //_____________________________________________________________________________
 442 void AliTrackPoint::Print(Option_t *) const
 443 {
 444   // Print the space-point coordinates and
 445   // covariance matrix
 446
 447   printf("VolumeID=%d\n", GetVolumeID());
 448   printf("X = %12.6f    Tx = %12.6f%12.6f%12.6f\n", fX, fCov[0], fCov[1], fCov[2]);
 449   printf("Y = %12.6f    Ty = %12.6f%12.6f%12.6f\n", fY, fCov[1], fCov[3], fCov[4]);
 450   printf("Z = %12.6f    Tz = %12.6f%12.6f%12.6f\n", fZ, fCov[2], fCov[4], fCov[5]);
 451
 452 }
 453
 454
 455 //________________________________
 456 void AliTrackPoint::SetAlignCovMatrix(const TMatrixDSym alignparmtrx){
 457   // Add the uncertainty on the cluster position due to alignment
 458   // (using the 6x6 AliAlignObj Cov. Matrix alignparmtrx) to the already
 459   // present Cov. Matrix
 460
 461   TMatrixDSym cov(3);
 462   TMatrixD coval(3,3);
 463   TMatrixD jacob(3,6);
 464   Float_t newcov[6];
 465
 466   cov(0,0)=fCov[0];
 467   cov(1,0)=cov(0,1)=fCov[1];
 468   cov(2,0)=cov(0,2)=fCov[2];
 469   cov(1,1)=fCov[3];
 470   cov(2,1)=cov(1,2)=fCov[4];
 471   cov(2,2)=fCov[5];
 472
 473   jacob(0,0) = 1;      jacob(1,0) = 0;       jacob(2,0) = 0;
 474   jacob(0,1) = 0;      jacob(1,1) = 1;       jacob(2,1) = 0;
 475   jacob(0,2) = 0;      jacob(1,2) = 0;       jacob(2,2) = 1;
 476   jacob(0,3) = 0;      jacob(1,3) =-fZ;      jacob(2,3) = fY;
 477   jacob(0,4) = fZ;     jacob(1,4) = 0;       jacob(2,4) =-fX;
 478   jacob(0,5) = fY;     jacob(1,5) = fX;      jacob(2,5) = 0;
 479
 480   TMatrixD jacobT=jacob.T();jacob.T();
 481
 482   coval=jacob*alignparmtrx*jacobT+cov;
 483
 484
 485   newcov[0]=coval(0,0);
 486   newcov[1]=coval(1,0);
 487   newcov[2]=coval(2,0);
 488   newcov[3]=coval(1,1);
 489   newcov[4]=coval(2,1);
 490   newcov[5]=coval(2,2);
 491
 492   SetXYZ(fX,fY,fZ,newcov);
 493
 494 }