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