// cvetan.cheshkov@cern.ch 3/11/2005 //
//////////////////////////////////////////////////////////////////////////////
+#include <TMath.h>
+#include <TMatrixD.h>
+#include <TMatrixDSym.h>
+#include <TGeoMatrix.h>
+#include <TMatrixDSymEigen.h>
+
#include "AliTrackPointArray.h"
ClassImp(AliTrackPointArray)
//______________________________________________________________________________
-AliTrackPointArray::AliTrackPointArray()
+AliTrackPointArray::AliTrackPointArray() :
+ TObject(),
+ fSorted(kFALSE),
+ fNPoints(0),
+ fX(0),
+ fY(0),
+ fZ(0),
+ fCharge(0),
+ fDriftTime(0),
+ fChargeRatio(0),
+ fIsExtra(0),
+ fSize(0),
+ fCov(0),
+ fVolumeID(0)
{
- fNPoints = fSize = 0;
- fX = fY = fZ = 0;
- fVolumeID = 0;
- fCov = 0;
}
//______________________________________________________________________________
AliTrackPointArray::AliTrackPointArray(Int_t npoints):
- fNPoints(npoints)
+ TObject(),
+ fSorted(kFALSE),
+ fNPoints(npoints),
+ fX(new Float_t[npoints]),
+ fY(new Float_t[npoints]),
+ fZ(new Float_t[npoints]),
+ fCharge(new Float_t[npoints]),
+ fDriftTime(new Float_t[npoints]),
+ fChargeRatio(new Float_t[npoints]),
+ fIsExtra(new Bool_t[npoints]),
+ fSize(6*npoints),
+ fCov(new Float_t[fSize]),
+ fVolumeID(new UShort_t[npoints])
{
// Constructor
//
- fSize = 6*npoints;
- fX = new Float_t[npoints];
- fY = new Float_t[npoints];
- fZ = new Float_t[npoints];
- fVolumeID = new UShort_t[npoints];
- fCov = new Float_t[fSize];
+ for (Int_t ip=0; ip<npoints;ip++){
+ fX[ip]=0;
+ fY[ip]=0;
+ fZ[ip]=0;
+ fCharge[ip]=0;
+ fDriftTime[ip]=0;
+ fChargeRatio[ip]=0;
+ fIsExtra[ip]=kFALSE;
+ fVolumeID[ip]=0;
+ for (Int_t icov=0;icov<6; icov++)
+ fCov[6*ip+icov]=0;
+ }
}
//______________________________________________________________________________
AliTrackPointArray::AliTrackPointArray(const AliTrackPointArray &array):
- TObject(array)
+ TObject(array),
+ fSorted(array.fSorted),
+ fNPoints(array.fNPoints),
+ fX(new Float_t[fNPoints]),
+ fY(new Float_t[fNPoints]),
+ fZ(new Float_t[fNPoints]),
+ fCharge(new Float_t[fNPoints]),
+ fDriftTime(new Float_t[fNPoints]),
+ fChargeRatio(new Float_t[fNPoints]),
+ fIsExtra(new Bool_t[fNPoints]),
+ fSize(array.fSize),
+ fCov(new Float_t[fSize]),
+ fVolumeID(new UShort_t[fNPoints])
{
// Copy constructor
//
- fNPoints = array.fNPoints;
- fSize = array.fSize;
- fX = new Float_t[fNPoints];
- fY = new Float_t[fNPoints];
- fZ = new Float_t[fNPoints];
- fVolumeID = new UShort_t[fNPoints];
- fCov = new Float_t[fSize];
memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
+ if (array.fCharge) {
+ memcpy(fCharge,array.fCharge,fNPoints*sizeof(Float_t));
+ } else {
+ memset(fCharge, 0, fNPoints*sizeof(Float_t));
+ }
+ if (array.fDriftTime) {
+ memcpy(fDriftTime,array.fDriftTime,fNPoints*sizeof(Float_t));
+ } else {
+ memset(fDriftTime, 0, fNPoints*sizeof(Float_t));
+ }
+ if (array.fChargeRatio) {
+ memcpy(fChargeRatio,array.fChargeRatio,fNPoints*sizeof(Float_t));
+ } else {
+ memset(fChargeRatio, 0, fNPoints*sizeof(Float_t));
+ }
+ if (array.fIsExtra) {
+ memcpy(fIsExtra,array.fIsExtra,fNPoints*sizeof(Bool_t));
+ } else {
+ memset(fIsExtra, 0, fNPoints*sizeof(Bool_t));
+ }
memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
}
if(this==&array) return *this;
((TObject *)this)->operator=(array);
+ fSorted = array.fSorted;
fNPoints = array.fNPoints;
fSize = array.fSize;
+ delete [] fX;
fX = new Float_t[fNPoints];
+ delete [] fY;
fY = new Float_t[fNPoints];
+ delete [] fZ;
fZ = new Float_t[fNPoints];
+ delete [] fCharge;
+ fCharge = new Float_t[fNPoints];
+ delete [] fDriftTime;
+ fDriftTime = new Float_t[fNPoints];
+ delete [] fChargeRatio;
+ fChargeRatio = new Float_t[fNPoints];
+ delete [] fIsExtra;
+ fIsExtra = new Bool_t[fNPoints];
+ delete [] fVolumeID;
fVolumeID = new UShort_t[fNPoints];
+ delete [] fCov;
fCov = new Float_t[fSize];
memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
+ memcpy(fCharge,array.fCharge,fNPoints*sizeof(Float_t));
+ memcpy(fDriftTime,array.fDriftTime,fNPoints*sizeof(Float_t));
+ memcpy(fChargeRatio,array.fChargeRatio,fNPoints*sizeof(Float_t));
+ memcpy(fIsExtra,array.fIsExtra,fNPoints*sizeof(Bool_t));
memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
delete [] fX;
delete [] fY;
delete [] fZ;
+ delete [] fCharge;
+ delete [] fDriftTime;
+ delete [] fChargeRatio;
+ delete [] fIsExtra;
delete [] fVolumeID;
delete [] fCov;
}
fX[i] = p->GetX();
fY[i] = p->GetY();
fZ[i] = p->GetZ();
+ fCharge[i] = p->GetCharge();
+ fDriftTime[i] = p->GetDriftTime();
+ fChargeRatio[i]=p->GetChargeRatio();
+ fIsExtra[i] = p->IsExtra();
fVolumeID[i] = p->GetVolumeID();
memcpy(&fCov[6*i],p->GetCov(),6*sizeof(Float_t));
return kTRUE;
}
+
//______________________________________________________________________________
Bool_t AliTrackPointArray::GetPoint(AliTrackPoint &p, Int_t i) const
{
if (i >= fNPoints) return kFALSE;
p.SetXYZ(fX[i],fY[i],fZ[i],&fCov[6*i]);
p.SetVolumeID(fVolumeID[i]);
+ p.SetCharge(fCharge ? fCharge[i] : 0);
+ p.SetDriftTime(fDriftTime ? fDriftTime[i] : 0);
+ p.SetChargeRatio(fChargeRatio ? fChargeRatio[i] : 0);
+ p.SetExtra(fIsExtra ? fIsExtra[i] : kFALSE);
return kTRUE;
}
return check;
}
+//______________________________________________________________________________
+void AliTrackPointArray::Sort(Bool_t down)
+{
+ // Sort the array by the values of Y-coordinate of the track points.
+ // The order is given by "down".
+ // Optimized more for maintenance rather than for speed.
+
+ if (fSorted) return;
+
+ Int_t *index=new Int_t[fNPoints];
+ AliTrackPointArray a(*this);
+ TMath::Sort(fNPoints,a.GetY(),index,down);
+
+ AliTrackPoint p;
+ for (Int_t i = 0; i < fNPoints; i++) {
+ a.GetPoint(p,index[i]);
+ AddPoint(i,&p);
+ }
+
+ delete[] index;
+ fSorted=kTRUE;
+}
+
ClassImp(AliTrackPoint)
//______________________________________________________________________________
-AliTrackPoint::AliTrackPoint()
+AliTrackPoint::AliTrackPoint() :
+ TObject(),
+ fX(0),
+ fY(0),
+ fZ(0),
+ fCharge(0),
+ fDriftTime(0),
+ fChargeRatio(0),
+ fIsExtra(kFALSE),
+ fVolumeID(0)
{
// Default constructor
//
- fX = fY = fZ = 0;
- fVolumeID = 0;
memset(fCov,0,6*sizeof(Float_t));
}
//______________________________________________________________________________
-AliTrackPoint::AliTrackPoint(Float_t x, Float_t y, Float_t z, const Float_t *cov, UShort_t volid)
+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) :
+ TObject(),
+ fX(0),
+ fY(0),
+ fZ(0),
+ fCharge(0),
+ fDriftTime(0),
+ fChargeRatio(0),
+ fIsExtra(kFALSE),
+ fVolumeID(0)
{
// Constructor
//
SetXYZ(x,y,z,cov);
+ SetCharge(charge);
+ SetDriftTime(drifttime);
+ SetChargeRatio(chargeratio);
SetVolumeID(volid);
}
//______________________________________________________________________________
-AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid)
+AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid, Float_t charge, Float_t drifttime,Float_t chargeratio) :
+ TObject(),
+ fX(0),
+ fY(0),
+ fZ(0),
+ fCharge(0),
+ fDriftTime(0),
+ fChargeRatio(0),
+ fIsExtra(kFALSE),
+ fVolumeID(0)
{
// Constructor
//
SetXYZ(xyz[0],xyz[1],xyz[2],cov);
+ SetCharge(charge);
+ SetDriftTime(drifttime);
+ SetChargeRatio(chargeratio);
SetVolumeID(volid);
}
//______________________________________________________________________________
AliTrackPoint::AliTrackPoint(const AliTrackPoint &p):
- TObject(p)
+ TObject(p),
+ fX(0),
+ fY(0),
+ fZ(0),
+ fCharge(0),
+ fDriftTime(0),
+ fChargeRatio(0),
+ fIsExtra(kFALSE),
+ fVolumeID(0)
{
// Copy constructor
//
SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
+ SetCharge(p.fCharge);
+ SetDriftTime(p.fDriftTime);
+ SetChargeRatio(p.fChargeRatio);
+ SetExtra(p.fIsExtra);
SetVolumeID(p.fVolumeID);
}
((TObject *)this)->operator=(p);
SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
+ SetCharge(p.fCharge);
+ SetDriftTime(p.fDriftTime);
+ SetChargeRatio(p.fChargeRatio);
+ SetExtra(p.fIsExtra);
SetVolumeID(p.fVolumeID);
return *this;
if (cov)
memcpy(cov,fCov,6*sizeof(Float_t));
}
+
+//______________________________________________________________________________
+Float_t AliTrackPoint::GetResidual(const AliTrackPoint &p, Bool_t weighted) const
+{
+ // This method calculates the track to space-point residuals. The track
+ // interpolation is also stored as AliTrackPoint. Using the option
+ // 'weighted' one can calculate the residual either with or without
+ // taking into account the covariance matrix of the space-point and
+ // track interpolation. The second case the residual becomes a pull.
+
+ Float_t res = 0;
+
+ if (!weighted) {
+ Float_t xyz[3],xyzp[3];
+ GetXYZ(xyz);
+ p.GetXYZ(xyzp);
+ res = (xyz[0]-xyzp[0])*(xyz[0]-xyzp[0])+
+ (xyz[1]-xyzp[1])*(xyz[1]-xyzp[1])+
+ (xyz[2]-xyzp[2])*(xyz[2]-xyzp[2]);
+ }
+ else {
+ Float_t xyz[3],xyzp[3];
+ Float_t cov[6],covp[6];
+ GetXYZ(xyz,cov);
+ TMatrixDSym mcov(3);
+ mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
+ mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
+ mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
+ p.GetXYZ(xyzp,covp);
+ TMatrixDSym mcovp(3);
+ mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
+ mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
+ mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
+ TMatrixDSym msum = mcov + mcovp;
+ msum.Invert();
+ // mcov.Print(); mcovp.Print(); msum.Print();
+ if (msum.IsValid()) {
+ for (Int_t i = 0; i < 3; i++)
+ for (Int_t j = 0; j < 3; j++)
+ res += (xyz[i]-xyzp[i])*(xyz[j]-xyzp[j])*msum(i,j);
+ }
+ }
+
+ return res;
+}
+
+//_____________________________________________________________________________
+Bool_t AliTrackPoint::GetPCA(const AliTrackPoint &p, AliTrackPoint &out) const
+{
+ //
+ // Get the intersection point between this point and
+ // the point "p" belongs to.
+ // The result is stored as a point 'out'
+ // return kFALSE in case of failure.
+ out.SetXYZ(0,0,0);
+
+ TMatrixD t(3,1);
+ t(0,0)=GetX();
+ t(1,0)=GetY();
+ t(2,0)=GetZ();
+
+ TMatrixDSym tC(3);
+ {
+ const Float_t *cv=GetCov();
+ tC(0,0)=cv[0]; tC(0,1)=cv[1]; tC(0,2)=cv[2];
+ tC(1,0)=cv[1]; tC(1,1)=cv[3]; tC(1,2)=cv[4];
+ tC(2,0)=cv[2]; tC(2,1)=cv[4]; tC(2,2)=cv[5];
+ }
+
+ TMatrixD m(3,1);
+ m(0,0)=p.GetX();
+ m(1,0)=p.GetY();
+ m(2,0)=p.GetZ();
+
+ TMatrixDSym mC(3);
+ {
+ const Float_t *cv=p.GetCov();
+ mC(0,0)=cv[0]; mC(0,1)=cv[1]; mC(0,2)=cv[2];
+ mC(1,0)=cv[1]; mC(1,1)=cv[3]; mC(1,2)=cv[4];
+ mC(2,0)=cv[2]; mC(2,1)=cv[4]; mC(2,2)=cv[5];
+ }
+
+ TMatrixDSym tmW(tC);
+ tmW+=mC;
+ tmW.Invert();
+ if (!tmW.IsValid()) return kFALSE;
+
+ TMatrixD mW(tC,TMatrixD::kMult,tmW);
+ TMatrixD tW(mC,TMatrixD::kMult,tmW);
+
+ TMatrixD mi(mW,TMatrixD::kMult,m);
+ TMatrixD ti(tW,TMatrixD::kMult,t);
+ ti+=mi;
+
+ TMatrixD iC(tC,TMatrixD::kMult,tmW);
+ iC*=mC;
+
+ out.SetXYZ(ti(0,0),ti(1,0),ti(2,0));
+ UShort_t id=p.GetVolumeID();
+ out.SetVolumeID(id);
+
+ return kTRUE;
+}
+
+//______________________________________________________________________________
+Float_t AliTrackPoint::GetAngle() const
+{
+ // The method uses the covariance matrix of
+ // the space-point in order to extract the
+ // orientation of the detector plane.
+ // The rotation in XY plane only is calculated.
+
+ Float_t phi= TMath::ATan2(TMath::Sqrt(fCov[0]),TMath::Sqrt(fCov[3]));
+ if (fCov[1] > 0) {
+ phi = TMath::Pi() - phi;
+ if ((fY-fX) < 0) phi += TMath::Pi();
+ }
+ else {
+ if ((fX+fY) < 0) phi += TMath::Pi();
+ }
+
+ return phi;
+
+}
+
+//______________________________________________________________________________
+Bool_t AliTrackPoint::GetRotMatrix(TGeoRotation& rot) const
+{
+ // Returns the orientation of the
+ // sensitive layer (using cluster
+ // covariance matrix).
+ // Assumes that cluster has errors only in the layer's plane.
+ // Return value is kTRUE in case of success.
+
+ TMatrixDSym mcov(3);
+ {
+ const Float_t *cov=GetCov();
+ mcov(0,0)=cov[0]; mcov(0,1)=cov[1]; mcov(0,2)=cov[2];
+ mcov(1,0)=cov[1]; mcov(1,1)=cov[3]; mcov(1,2)=cov[4];
+ mcov(2,0)=cov[2]; mcov(2,1)=cov[4]; mcov(2,2)=cov[5];
+ }
+
+ TMatrixDSymEigen eigen(mcov);
+ TMatrixD eigenMatrix = eigen.GetEigenVectors();
+
+ rot.SetMatrix(eigenMatrix.GetMatrixArray());
+
+ return kTRUE;
+}
+
+
+//_____________________________________________________________________________
+AliTrackPoint& AliTrackPoint::Rotate(Float_t alpha) const
+{
+ // Transform the space-point coordinates
+ // and covariance matrix from global to
+ // local (detector plane) coordinate system
+ // XY plane rotation only
+
+ static AliTrackPoint p;
+ p = *this;
+
+ Float_t xyz[3],cov[6];
+ GetXYZ(xyz,cov);
+
+ Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
+
+ Float_t newxyz[3],newcov[6];
+ newxyz[0] = cos*xyz[0] + sin*xyz[1];
+ newxyz[1] = cos*xyz[1] - sin*xyz[0];
+ newxyz[2] = xyz[2];
+
+ newcov[0] = cov[0]*cos*cos+
+ 2*cov[1]*sin*cos+
+ cov[3]*sin*sin;
+ newcov[1] = cov[1]*(cos*cos-sin*sin)+
+ (cov[3]-cov[0])*sin*cos;
+ newcov[2] = cov[2]*cos+
+ cov[4]*sin;
+ newcov[3] = cov[0]*sin*sin-
+ 2*cov[1]*sin*cos+
+ cov[3]*cos*cos;
+ newcov[4] = cov[4]*cos-
+ cov[2]*sin;
+ newcov[5] = cov[5];
+
+ p.SetXYZ(newxyz,newcov);
+ p.SetVolumeID(GetVolumeID());
+
+ return p;
+}
+
+//_____________________________________________________________________________
+AliTrackPoint& AliTrackPoint::MasterToLocal() const
+{
+ // Transform the space-point coordinates
+ // and the covariance matrix from the
+ // (master) to the local (tracking)
+ // coordinate system
+
+ Float_t alpha = GetAngle();
+ return Rotate(alpha);
+}
+
+//_____________________________________________________________________________
+void AliTrackPoint::Print(Option_t *) const
+{
+ // Print the space-point coordinates and
+ // covariance matrix
+
+ printf("VolumeID=%d\n", GetVolumeID());
+ printf("X = %12.6f Tx = %12.6f%12.6f%12.6f\n", fX, fCov[0], fCov[1], fCov[2]);
+ printf("Y = %12.6f Ty = %12.6f%12.6f%12.6f\n", fY, fCov[1], fCov[3], fCov[4]);
+ printf("Z = %12.6f Tz = %12.6f%12.6f%12.6f\n", fZ, fCov[2], fCov[4], fCov[5]);
+ printf("Charge = %f\n", fCharge);
+ printf("Drift Time = %f\n", fDriftTime);
+ printf("Charge Ratio = %f\n", fChargeRatio);
+ if(fIsExtra) printf("This is an extra point\n");
+
+}
+
+
+//________________________________
+void AliTrackPoint::SetAlignCovMatrix(const TMatrixDSym alignparmtrx){
+ // Add the uncertainty on the cluster position due to alignment
+ // (using the 6x6 AliAlignObj Cov. Matrix alignparmtrx) to the already
+ // present Cov. Matrix
+
+ TMatrixDSym cov(3);
+ TMatrixD coval(3,3);
+ TMatrixD jacob(3,6);
+ Float_t newcov[6];
+
+ cov(0,0)=fCov[0];
+ cov(1,0)=cov(0,1)=fCov[1];
+ cov(2,0)=cov(0,2)=fCov[2];
+ cov(1,1)=fCov[3];
+ cov(2,1)=cov(1,2)=fCov[4];
+ cov(2,2)=fCov[5];
+
+ jacob(0,0) = 1; jacob(1,0) = 0; jacob(2,0) = 0;
+ jacob(0,1) = 0; jacob(1,1) = 1; jacob(2,1) = 0;
+ jacob(0,2) = 0; jacob(1,2) = 0; jacob(2,2) = 1;
+ jacob(0,3) = 0; jacob(1,3) =-fZ; jacob(2,3) = fY;
+ jacob(0,4) = fZ; jacob(1,4) = 0; jacob(2,4) =-fX;
+ jacob(0,5) = -fY; jacob(1,5) = fX; jacob(2,5) = 0;
+
+ TMatrixD jacobT=jacob.T();jacob.T();
+
+ coval=jacob*alignparmtrx*jacobT+cov;
+
+
+ newcov[0]=coval(0,0);
+ newcov[1]=coval(1,0);
+ newcov[2]=coval(2,0);
+ newcov[3]=coval(1,1);
+ newcov[4]=coval(2,1);
+ newcov[5]=coval(2,2);
+
+ SetXYZ(fX,fY,fZ,newcov);
+
+}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff