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