[u/mrichter/AliRoot.git] / STEER / AliTrackPointArray.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//                          Class AliTrackPointArray                        //
//   This class contains the ESD track space-points which are used during   //
//   the alignment procedures. Each space-point consist of 3 coordinates    //
//   (and their errors) and the index of the sub-detector which contains    //
//   the space-point.                                                       //
//   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() :
  TObject(),
  fSorted(kFALSE),
  fNPoints(0),
  fX(0),
  fY(0),
  fZ(0),
  fCharge(0),
  fDriftTime(0),
  fSize(0),
  fCov(0),
  fVolumeID(0)
{
}

//______________________________________________________________________________
AliTrackPointArray::AliTrackPointArray(Int_t 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]),
  fSize(6*npoints),
  fCov(new Float_t[fSize]),
  fVolumeID(new UShort_t[npoints])
{
  // Constructor
  //
  for (Int_t ip=0; ip<npoints;ip++){
    fX[ip]=0;
    fY[ip]=0;
    fZ[ip]=0;
    fCharge[ip]=0;
    fDriftTime[ip]=0;
    fVolumeID[ip]=0;
    for (Int_t icov=0;icov<6; icov++)
      fCov[6*ip+icov]=0;
  }
}

//______________________________________________________________________________
AliTrackPointArray::AliTrackPointArray(const AliTrackPointArray &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]),
  fSize(array.fSize),
  fCov(new Float_t[fSize]),
  fVolumeID(new UShort_t[fNPoints])
{
  // Copy constructor
  //
  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));
  }
  memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
  memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
}

//_____________________________________________________________________________
AliTrackPointArray &AliTrackPointArray::operator =(const AliTrackPointArray& array)
{
  // assignment operator
  //
  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 [] 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(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
  memcpy(fCov,array.fCov,fSize*sizeof(Float_t));

  return *this;
}

//______________________________________________________________________________
AliTrackPointArray::~AliTrackPointArray()
{
  // Destructor
  //
  delete [] fX;
  delete [] fY;
  delete [] fZ;
  delete [] fCharge;
  delete [] fDriftTime;
  delete [] fVolumeID;
  delete [] fCov;
}


//______________________________________________________________________________
Bool_t AliTrackPointArray::AddPoint(Int_t i, const AliTrackPoint *p)
{
  // Add a point to the array at position i
  //
  if (i >= fNPoints) return kFALSE;
  fX[i] = p->GetX();
  fY[i] = p->GetY();
  fZ[i] = p->GetZ();
  fCharge[i] = p->GetCharge();
  fDriftTime[i] = p->GetDriftTime();
  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
{
  // Get the point at position i
  //
  if (i >= fNPoints) return kFALSE;
  p.SetXYZ(fX[i],fY[i],fZ[i],&fCov[6*i]);
  p.SetVolumeID(fVolumeID[i]);
  p.SetCharge(fCharge[i]);
  p.SetDriftTime(fDriftTime[i]);
  return kTRUE;
}

//______________________________________________________________________________
Bool_t AliTrackPointArray::HasVolumeID(UShort_t volid) const
{
  // This method checks if the array
  // has at least one hit in the detector
  // volume defined by volid
  Bool_t check = kFALSE;
  for (Int_t ipoint = 0; ipoint < fNPoints; ipoint++)
    if (fVolumeID[ipoint] == volid) check = 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() :
  TObject(),
  fX(0),
  fY(0),
  fZ(0),
  fCharge(0),
  fDriftTime(0),
  fVolumeID(0)
{
  // Default constructor
  //
  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, Float_t charge, Float_t drifttime) :
  TObject(),
  fX(0),
  fY(0),
  fZ(0),
  fCharge(0),
  fDriftTime(0),
  fVolumeID(0)
{
  // Constructor
  //
  SetXYZ(x,y,z,cov);
  SetCharge(charge);
  SetDriftTime(drifttime);
  SetVolumeID(volid);
}

//______________________________________________________________________________
AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid, Float_t charge, Float_t drifttime) :
  TObject(),
  fX(0),
  fY(0),
  fZ(0),
  fCharge(0),
  fDriftTime(0),
  fVolumeID(0)
{
  // Constructor
  //
  SetXYZ(xyz[0],xyz[1],xyz[2],cov);
  SetCharge(charge);
  SetDriftTime(drifttime);
  SetVolumeID(volid);
}

//______________________________________________________________________________
AliTrackPoint::AliTrackPoint(const AliTrackPoint &p):
  TObject(p),
  fX(0),
  fY(0),
  fZ(0),
  fCharge(0),
  fDriftTime(0),
  fVolumeID(0)
{
  // Copy constructor
  //
  SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
  SetCharge(p.fCharge);
  SetDriftTime(p.fDriftTime);
  SetVolumeID(p.fVolumeID);
}

//_____________________________________________________________________________
AliTrackPoint &AliTrackPoint::operator =(const AliTrackPoint& p)
{
  // assignment operator
  //
  if(this==&p) return *this;
  ((TObject *)this)->operator=(p);

  SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
  SetCharge(p.fCharge);
  SetDriftTime(p.fDriftTime);
  SetVolumeID(p.fVolumeID);

  return *this;
}

//______________________________________________________________________________
void AliTrackPoint::SetXYZ(Float_t x, Float_t y, Float_t z, const Float_t *cov)
{
  // Set XYZ coordinates and their cov matrix
  //
  fX = x;
  fY = y;
  fZ = z;
  if (cov)
    memcpy(fCov,cov,6*sizeof(Float_t));
}

//______________________________________________________________________________
void AliTrackPoint::SetXYZ(const Float_t *xyz, const Float_t *cov)
{
  // Set XYZ coordinates and their cov matrix
  //
  SetXYZ(xyz[0],xyz[1],xyz[2],cov);
}

//______________________________________________________________________________
void AliTrackPoint::GetXYZ(Float_t *xyz, Float_t *cov) const
{
  xyz[0] = fX;
  xyz[1] = fY;
  xyz[2] = fZ;
  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);

}


//________________________________
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);

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