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