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