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