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