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