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