// ----------------------------------------------------------------------------
// $Id: LCDTrack.cxx,v 1.1 2001/05/10 17:25:09 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDTrack.cxx,v $
// Revision 1.1 2001/05/10 17:25:09 toshi
// Rename *.C to *.cxx.
//
// Revision 1.2 2001/04/28 23:46:37 toshi
// Start to use CVS.
//
//
///////////////////////////////////////////////////////////////////////////
//
// LCDTrack
//
// The Track class is used to store the data from the Track System of the
// detector. This includes the momentum, position, charge as well as a
// pointer to the corresponding McPart object.
//
///////////////////////////////////////////////////////////////////////////
#include "LCDTrack.h"
#include "LCDMcPart.h"
ClassImp(LCDTrack)
///________________________________________________________________________
LCDTrack::LCDTrack(Int_t index,
Int_t index0,
Double_t bfld_z,
Int_t ndf,
Double_t chi2,
Double_t* par,
Double_t* errMat ) {
SetLCDTrack(index,index0,bfld_z,ndf,chi2,par,errMat);
}
LCDTrack::LCDTrack(Int_t index,
Int_t index0,
Double_t bfld_z,
Int_t ndf,
Double_t chi2,
Double_t qMC,
TVector3 pMC,
TVector3 xMC,
Double_t* errMat ) {
Double_t par[5];
CalcTrackParameters(qMC,bfld_z,pMC,xMC,par);
SetLCDTrack(index,index0,bfld_z,ndf,chi2,par,errMat);
}
LCDTrack::LCDTrack(TClonesArray* mclist,
Int_t index,
Double_t bfld_z,
Int_t ndf,
Double_t chi2,
Double_t* errMat) {
LCDMcPart* mc =(LCDMcPart*)mclist->At(index );
LCDMcPart* mcp=(LCDMcPart*)mclist->At(mc->GetParentIdx());
Int_t index0=0;
TVector3 xMC(0.0,0.0,0.0);
if (mcp != 0 ) {
index0=mclist->IndexOf((TObject*)mcp);
xMC = mcp->GetTermPosition();
}
Double_t qMC=mc->GetCharge();
TVector3 pMC=(mc->Get4Momentum()).Vect();
Double_t par[5];
CalcTrackParameters(qMC,bfld_z,pMC,xMC,par);
SetLCDTrack(index,index0,bfld_z,ndf,chi2,par,errMat);
}
LCDTrack::LCDTrack(LCDTrack* trk) {
SetLCDTrack(trk->GetParticle(),trk->GetParent(),trk->GetMagneticField(),
trk->GetNdf(),trk->GetChi2(),trk->GetTrackParameters(),
trk->GetErrorMatrix());
}
LCDTrack::~LCDTrack() {
// delete trackErrorMatrix;
}
void LCDTrack::SetLCDTrack(Int_t index,
Int_t index0,
Double_t bfld_z,
Int_t ndf,
Double_t chi2,
Double_t* par,
Double_t* errMat ) {
// Create a record of a track in the track system
m_magneticfield=bfld_z;
m_charge = ( par[2] < 0. ) ? 1. : -1.;
m_index = index;
// Temp
m_index0 = index0;
m_ndf = ndf;
m_chisquared = chi2;
for (Int_t j=0 ; j < 5 ; j++) {
m_trackPar[j] = par[j];
}
for (Int_t k=0 ; k < 15 ; k++) {
m_trackErrorMatrix[k] = errMat[k];
}
}
void LCDTrack::CalcPOCA(Double_t qMC,
Double_t bfld_z,
TVector3 p_orig,
TVector3 x_orig,
TVector3& p_poca,
TVector3& x_poca) {
Double_t pT = p_orig.Pt();
Double_t tanL = p_orig.Z()/pT;
Double_t arho = 333.567/bfld_z*pT;
TVector3 nzv(0.0,0.0,qMC/TMath::Abs(qMC));
TVector3 xxc = p_orig.Cross(nzv);
xxc.SetZ(0.0);
TVector3 nxxc = xxc.Unit();
TVector3 xc = x_orig + arho*nxxc;
xc.SetZ(0.0);
TVector3 nxc = xc.Unit();
TVector3 catMC=nzv.Cross(nxc);
catMC.SetZ(0.0);
TVector3 ncatMC=catMC.Unit();
p_poca = pT*ncatMC;
p_poca.SetZ(p_orig.Z());
Double_t dphi = TMath::ASin(nxxc.Cross(nxc).Z());
Double_t dl =-dphi*arho*qMC;
x_poca = xc - arho*nxc;
x_poca.SetZ(x_orig.Z() + dl*tanL);
}
// get track parameter
void LCDTrack::CalcTrackParameters(Double_t qMC,Double_t bfld_z,
TVector3 pMC,TVector3 xMC, Double_t* tkpar){
TVector3 capMC;
TVector3 caxMC;
CalcPOCA(qMC,bfld_z,pMC,xMC,capMC,caxMC);
Double_t pT= capMC.Pt();
Double_t d0= caxMC.Perp();
if (caxMC.Cross(capMC).Z() > 0) {
d0 = -d0;
}
Double_t phi0 = capMC.Phi();
Double_t z0 = caxMC.Z();
Double_t tanL = capMC.Z()/pT;
tkpar[0]= d0;
tkpar[1]= phi0;
tkpar[2]=-qMC*bfld_z/333.567/pT;
tkpar[3]= z0;
tkpar[4]= tanL;
}
void LCDTrack::WriteTrack(FILE* ofile)const {
fprintf(ofile,
"%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %i %e %e %i %in",
m_trackPar[0],
m_trackPar[1],
m_trackPar[2],
m_trackPar[3],
m_trackPar[4],
m_trackErrorMatrix[ 0],
m_trackErrorMatrix[ 1],
m_trackErrorMatrix[ 2],
m_trackErrorMatrix[ 3],
m_trackErrorMatrix[ 4],
m_trackErrorMatrix[ 5],
m_trackErrorMatrix[ 6],
m_trackErrorMatrix[ 7],
m_trackErrorMatrix[ 8],
m_trackErrorMatrix[ 9],
m_trackErrorMatrix[10],
m_trackErrorMatrix[11],
m_trackErrorMatrix[12],
m_trackErrorMatrix[13],
m_trackErrorMatrix[14],
m_ndf,
m_chisquared,
m_magneticfield,
(int)m_index,
(int)m_index0
);
}
Double_t LCDTrack::GetTrackParameter(Int_t i) {
return m_trackPar[i];
}
Double_t* LCDTrack::GetTrackParameters() {
return m_trackPar;
}
Double_t LCDTrack::GetPositionX(Double_t l) {
Double_t d0 = m_trackPar[0];
Double_t phi0 = m_trackPar[1];
Double_t omega= m_trackPar[2];
Double_t phi = phi0+omega*l;
return 1.0/omega*TMath::Sin(phi) - (1.0/omega + d0)*TMath::Sin(phi0);
}
Double_t LCDTrack::GetPositionY(Double_t l) {
Double_t d0 = m_trackPar[0];
Double_t phi0 = m_trackPar[1];
Double_t omega= m_trackPar[2];
Double_t phi = phi0+omega*l;
return -1.0/omega*TMath::Cos(phi) + (1.0/omega + d0)*TMath::Cos(phi0);
}
Double_t LCDTrack::GetPositionZ(Double_t l) {
Double_t z0 = m_trackPar[3];
Double_t tanl = m_trackPar[4];
return z0 + l*tanl;
}
Double_t* LCDTrack::GetPosition() {
m_position[0]=GetPositionX(0.0);
m_position[1]=GetPositionY(0.0);
m_position[2]=GetPositionZ(0.0);
return m_position;
}
Double_t* LCDTrack::GetPosition(Double_t l) {
m_position[0]=GetPositionX(l);
m_position[1]=GetPositionY(l);
m_position[2]=GetPositionZ(l);
return m_position;
}
TVector3 LCDTrack::GetPositionVector(Double_t l) {
return TVector3(GetPositionX(l),GetPositionY(l),GetPositionZ(l));
}
Double_t LCDTrack::GetMomentumPx(Double_t l) {
//Double_t d0 = m_trackPar[0];
Double_t phi0 = m_trackPar[1];
Double_t omega= m_trackPar[2];
//Double_t z0 = m_trackPar[3];
// Double_t tanl = m_trackPar[4];
Double_t phi = phi0+omega*l;
Double_t pt = m_magneticfield/333.567/TMath::Abs(m_trackPar[2]);
return pt*TMath::Cos(phi);
}
Double_t LCDTrack::GetMomentumPy(Double_t l) {
//Double_t d0 = m_trackPar[0];
Double_t phi0 = m_trackPar[1];
Double_t omega= m_trackPar[2];
//Double_t z0 = m_trackPar[3];
// Double_t tanl = m_trackPar[4];
Double_t phi = phi0+omega*l;
Double_t pt = m_magneticfield/333.567/TMath::Abs(m_trackPar[2]);
return pt*TMath::Sin(phi);
}
Double_t LCDTrack::GetMomentumPz(Double_t l) {
//Double_t d0 = m_trackPar[0];
//Double_t phi0 = m_trackPar[1];
Double_t omega= m_trackPar[2];
//Double_t z0 = m_trackPar[3];
Double_t tanl = m_trackPar[4];
//Double_t phi = phi0+omega*l;
Double_t pt = m_magneticfield/333.567/TMath::Abs(omega);
return pt*tanl;
}
TVector3 LCDTrack::GetMomentumVector(Double_t l) {
return TVector3(GetMomentumPx(l),GetMomentumPy(l),GetMomentumPz(l));
}
Double_t* LCDTrack::GetMomentum() {
m_momentum[0]=GetMomentumPx(0.0);
m_momentum[1]=GetMomentumPy(0.0);
m_momentum[2]=GetMomentumPz(0.0);
return m_momentum;
}
Double_t* LCDTrack::GetMomentum(Double_t l) {
m_momentum[0]=GetMomentumPx(l);
m_momentum[1]=GetMomentumPy(l);
m_momentum[2]=GetMomentumPz(l);
return m_momentum;
}
Double_t LCDTrack::GetDistance(Double_t l, Double_t x, Double_t y, Double_t z) {
Double_t x0=GetPositionX(l);
Double_t y0=GetPositionY(l);
Double_t z0=GetPositionZ(l);
return TMath::Sqrt((x0 - x)*(x0 - x) + (y0 - y)*(y0 - y) + (z0 - z)*(z0 - z));
}
Double_t LCDTrack::GetDistance(TVector3 xpt) {
Double_t l =CalcPOCA3(xpt);
TVector3 dx=xpt-GetPositionVector(l);
return dx.Mag();
}
Double_t LCDTrack::GetDistance(Double_t x, Double_t y, Double_t z) {
TVector3 xpt(x,y,z);
return GetDistance(xpt);
}
Double_t LCDTrack::GetErrorMatrixElement(Int_t i, Int_t j) {
Int_t ij[5][5]= { { 0, 1, 3, 6,10},
{ 1, 2, 4, 7,11},
{ 3, 4, 5, 8,12},
{ 6, 7, 8, 9,13},
{10,11,12,13,14}};
return m_trackErrorMatrix[ij[i][j]];
}
void LCDTrack::GetErrorMatrix(TMatrixD& trackErrorMatrix){
Int_t ij[5][5]= { { 0, 1, 3, 6,10},
{ 1, 2, 4, 7,11},
{ 3, 4, 5, 8,12},
{ 6, 7, 8, 9,13},
{10,11,12,13,14}};
Int_t ip=0,jp=0;
for (ip=0 ; ip < 5 ; ip++) {
for (jp=0 ; jp < 5 ; jp++) {
trackErrorMatrix(ip,jp) = m_trackErrorMatrix[ij[ip][jp]];
}
}
}
void LCDTrack::GetErrorMatrixPosition(Double_t l,TMatrixD& var2) {
TMatrixD var0(5,5); GetErrorMatrix(var0);
TMatrixD dxda(3,5); Getdxda(l,dxda);
TMatrixD dxdaT(TMatrixD::kTransposed,dxda);
TMatrixD var1(var0,TMatrixD::kMult,dxdaT);
var2.Mult(dxda,var1);
}
void LCDTrack::GetErrorMatrixMomentum(Double_t l,TMatrixD& var2) {
TMatrixD var0(5,5); GetErrorMatrix(var0);
TMatrixD dpda(3,5); Getdpda(l,dpda);
TMatrixD dpdaT(TMatrixD::kTransposed,dpda);
TMatrixD var1(var0,TMatrixD::kMult,dpdaT);
var2.Mult(dpda,var1);
}
Double_t LCDTrack::GetErrorDistance(Double_t l, Double_t x, Double_t y, Double_t z) {
TMatrixD ddda(1,5); Getddda(l,x,y,z,ddda);
TMatrixD dddaT(TMatrixD::kTransposed,ddda);
TMatrixD em(5,5); GetErrorMatrix(em);
TMatrixD tmp0(em,TMatrixD::kMult,dddaT);
TMatrixD tmp1(ddda,TMatrixD::kMult,tmp0);
return tmp1(0,0);
}
Double_t LCDTrack::GetErrorDistance(TVector3 xpt) {
Double_t l=CalcPOCA3(xpt);
return GetErrorDistance(l,xpt[0],xpt[1],xpt[2]);
}
Double_t LCDTrack::GetErrorDistance(Double_t x, Double_t y, Double_t z) {
TVector3 xpt(x,y,z);
return GetErrorDistance(xpt);
}
void LCDTrack::Getdxda(Double_t l, TMatrixD& dxda) {
Double_t d0 = m_trackPar[0];
Double_t phi0 = m_trackPar[1];
Double_t omega = m_trackPar[2];
//Double_t z0 = m_trackPar[3];
//Double_t tanl = m_trackPar[4];
Double_t phi = phi0 + omega*l;
Double_t csphi0= TMath::Cos(phi0);
Double_t snphi0= TMath::Sin(phi0);
Double_t csphi = TMath::Cos(phi);
Double_t snphi = TMath::Sin(phi);
dxda(0,0) =-snphi0;
dxda(0,1) = 1.0/omega*csphi - (1.0/omega + d0)*csphi0;
dxda(0,2) = (1.0/omega*(-snphi + snphi0) + csphi*l)/omega;
dxda(0,3) = 0.0;
dxda(0,4) = 0.0;
dxda(1,0) = csphi0;
dxda(1,1) = 1.0/omega*snphi - (1.0/omega + d0)*snphi0;
dxda(1,2) = (1.0/omega*( csphi - csphi0) + snphi*l)/omega;
dxda(1,3) = 0.0;
dxda(1,4) = 0.0;
dxda(2,0) = 0.0;
dxda(2,1) = 0.0;
dxda(2,2) = 0.0;
dxda(2,3) = 1.0;
dxda(2,4) = l;
}
void LCDTrack::Getdpda(Double_t l, TMatrixD& dxda) {
//Double_t d0 = m_trackPar[0];
Double_t phi0 = m_trackPar[1];
Double_t omega = m_trackPar[2];
//Double_t z0 = m_trackPar[3];
Double_t tanl = m_trackPar[4];
Double_t phi = phi0 + omega*l;
// Double_t csphi0= TMath::Cos(phi0);
// Double_t snphi0= TMath::Sin(phi0);
Double_t csphi = TMath::Cos(phi);
Double_t snphi = TMath::Sin(phi);
Double_t r2p = TMath::Abs(m_magneticfield)/333.567;
TMatrixD dpda(5,3);
dpda(0,0) = 0.0;
dpda(1,0) =-r2p/omega*snphi;
dpda(2,0) =-r2p/omega/omega*csphi - r2p/omega*snphi*l;
dpda(3,0) = 0.0;
dpda(4,0) = 0.0;
dpda(0,1) = 0.0;
dpda(1,1) = r2p/omega*csphi;
dpda(2,1) =-r2p/omega/omega*snphi + r2p/omega*csphi*l;
dpda(3,1) = 0.0;
dpda(4,1) = 0.0;
dpda(0,2) = 0.0;
dpda(1,2) = 0.0;
dpda(2,2) =-r2p/omega/omega*tanl;
dpda(3,2) = 0.0;
dpda(4,2) = r2p/omega;
}
void LCDTrack::Getddda(Double_t l,
Double_t x, Double_t y, Double_t z,
TMatrixD& ddda) {
Double_t xp=GetPositionX(l);
Double_t yp=GetPositionY(l);
Double_t zp=GetPositionZ(l);
Double_t d =GetDistance(l,x,y,z);
TMatrixD dx(1,3);
dx(0,0) = -(x - xp)/d;
dx(0,1) = -(y - yp)/d;
dx(0,2) = -(z - zp)/d;
TMatrixD dxda(3,5); Getdxda(l,dxda);
ddda.Mult(dx,dxda);
}
void LCDTrack::Getd2xd2l(Double_t l, TMatrixD& d2xd2l) {
Double_t phi0 = m_trackPar[1];
Double_t omega = m_trackPar[2];
Double_t phi = phi0 + omega*l;
d2xd2l(0,0)=-omega*TMath::Sin(phi);
d2xd2l(1,0)= omega*TMath::Cos(phi);
d2xd2l(2,0)= 0.0;
}
void LCDTrack::Getdxdl(Double_t l,TMatrixD& dxdl) {
Double_t phi0=GetTrackParameter(1);
Double_t omega=GetTrackParameter(2);
Double_t tanl=GetTrackParameter(4);
Double_t phi=phi0+omega*l;
dxdl(0,0)=TMath::Cos(phi);
dxdl(1,0)=TMath::Sin(phi);
dxdl(2,0)=tanl;
}
// Get Closest Approach Point to beam axis...
Double_t LCDTrack::CalcPOCA3(TVector3 x_given) {
const Int_t MAXTRY=50;
const Double_t DLEN2TH=1E-10;
Double_t d0 =m_trackPar[0];
Double_t phi0 =m_trackPar[1];
Double_t omega=m_trackPar[2];
TVector3 xc=(-(d0+1.0/omega)*TMath::Sin(phi0),
+(d0+1.0/omega)*TMath::Cos(phi0),
+0.0);
TVector3 xxc=x_given-xc;
xxc.SetZ(0.0);
TVector3 nxxc=xxc; nxxc.Unit();
TVector3 nxc=xc ; nxc.Unit();
nxc *=-1.0;
TVector3 nxcXnxxc=nxc; nxcXnxxc.Cross(nxxc);
Double_t u=nxcXnxxc.Z();
if (u <-1.0) { u=-1.0; }
else if (u > 1.0) { u= 1.0; }
Double_t l=1.0/omega*TMath::ASin(u);
Double_t l_best=0.0,ddlen2dl_best=0.0,d2dlen2d2l_best=0.0;
Double_t dlen2_c=0;
Double_t dlen2_p=1e10;
Double_t a=0;
Double_t f=1.0;
Int_t i=0;
for (i=0; i < MAXTRY ; i++) {
dlen2_c=dlen2(l,x_given);
if ( i != 0 && TMath::Abs(dlen2_c - dlen2_p) < DLEN2TH) {
break;
} else if ( i >= MAXTRY-1 ) {
l = l_best;
break;
} else if ( i != 0 && dlen2_c > dlen2_p ) {
f *= 3.0;
l = l_best;
} else {
a = d2dlen2d2l(l,x_given);
if (a != 0) {
dlen2_p = dlen2_c;
f = 1.0;
l_best = l;
ddlen2dl_best = ddlen2dl(l,x_given);
d2dlen2d2l_best = a;
} else {
f *= 3.0;
l = l_best;
}
}
l -= ddlen2dl_best/d2dlen2d2l_best/(1.0+f);
}
return l;
}
Double_t LCDTrack::dlen2(Double_t l, TVector3 x_given) {
TVector3 dx=x_given - GetPositionVector(l);
return dx.Mag2();
}
Double_t LCDTrack::ddlen2dl(Double_t l,TVector3 x_given)
{
Double_t d0 =m_trackPar[0];
Double_t phi0 =m_trackPar[1];
Double_t omega=m_trackPar[2];
Double_t z0 =m_trackPar[3];
Double_t tanl =m_trackPar[4];
Double_t x0p =x_given.X();
Double_t y0p =x_given.Y();
Double_t z0p =x_given.Z();
Double_t tmp=2*(d0 + 1.0/omega)*TMath::Sin(omega*l);
tmp += 2.0*z0*tanl + 2.0*l*tanl*tanl;
tmp += -2*x0p*TMath::Cos(phi0+omega*l);
tmp += - 2*y0p*TMath::Sin(phi0+omega*l);
tmp += -2*z0p*tanl;
return tmp;
// return 2*(d0 + 1.0/omega)*TMath::Sin(omega*l) + 2.0*z0*tanl + 2.0*l*tanl*tanl
// -2*x0p*TMath::Cos(phi0+omega*l) - 2*y0p*TMath::Sin(phi0+omega*l)-2*z0p*tanl;
}
Double_t LCDTrack::d2dlen2d2l(Double_t l,TVector3 x_given)
{
Double_t d0 =m_trackPar[0];
Double_t phi0 =m_trackPar[1];
Double_t omega=m_trackPar[2];
// Double_t z0 =m_trackPar[3];
Double_t tanl =m_trackPar[4];
Double_t x0p =x_given.X();
Double_t y0p =x_given.Y();
// Double_t z0p =x_given.Z();
Double_t tmp= 2*omega*(d0 + 1.0/omega)*TMath::Cos(omega*l);
tmp += 2.0*tanl*tanl;
tmp += 2*omega*x0p*TMath::Sin(phi0+omega*l);
tmp -= 2*omega*y0p*TMath::Cos(phi0+omega*l);
return tmp;
// return 2*omega*(d0 + 1.0/omega)*TMath::Cos(omega*l) + 2.0*tanl*tanl
// + 2*omega*x0p*TMath::Sin(phi0+omega*l)-2*omega*y0p*TMath::Cos(phi0+omega*l);
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.