// ----------------------------------------------------------------------------
// $Id: LCDVToplTRK.cxx,v 1.2 2001/06/19 03:47:27 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDVToplTRK.cxx,v $
// Revision 1.2 2001/06/19 03:47:27 toshi
// Changes to allow direct LCDVToplVRT pointers in LCDVToplTRK.
//
// Revision 1.1 2001/05/10 19:32:12 toshi
// C++ file name convention has been changed from .C to .cxx .
//
//
#include "LCDVToplTRK.h"
////////////////////////////////////////////////////////////
//
// LCDVToplTRK
//
////////////////////////////////////////////////////////////
ClassImp(LCDVToplTRK)
//__________________________________________________________________________
//Constructor&deconstructor
LCDVToplTRK::LCDVToplTRK() {
Init();
}
LCDVToplTRK::LCDVToplTRK(Int_t tkid, Double_t* tkpar, Double_t* e_tkpar,
Double_t chrg,
Double_t mfld, Int_t mc_ind,const TVector3& ipv) {
Init(tkid,tkpar,e_tkpar,chrg,mfld,mc_ind,ipv);
}
LCDVToplTRK::LCDVToplTRK(TObjArray* track_list,
Int_t trk_id,
const TVector3& ipv) {
LCDTrack* trk=(LCDTrack*)(track_list->At(trk_id));
Double_t* trk_parm=trk->GetTrackParameters();
Double_t emtrk[15];
Int_t i=0,j=0,k=0;
for(i=0; i < 5 ; i++) {
for(j=0; j <= i ; j++ ) {
emtrk[k++] = trk->GetErrorMatrixElement(i,j);
}
}
Init(trk_id,trk_parm,emtrk,trk->GetCharge(),trk->GetMagneticField(),
trk->GetParticle(),ipv);
}
//Initialize.
void LCDVToplTRK::Init(){
m_nvrt=0;
m_track_id=-1;
Int_t i=0;
for(i=0; i < 5 ; i++) {
m_trackPar[i] = 0.0;
}
Int_t k=0;
for(k=0; k < 15 ; k++) {
m_trackErrorMatrix[k++] = 0;
}
m_charge=0;
m_magneticfield=0;
m_index = -1;
m_trkval = 0;
m_rmax=2.0;
Int_t j;
for (j=0 ; j < 3 ; j++) {
m_err_tsi[j]=1E10;
m_err_eta[j]=1E10;
}
m_v0flag=0;
}
void LCDVToplTRK::Init(Int_t tkid,Double_t* tkpar, Double_t* e_tkpar,
Double_t chrg,
Double_t mfld,Int_t mc_ind,
const TVector3& ipv){
Init();
SetUp(tkid,tkpar,e_tkpar,chrg,mfld,mc_ind,ipv);
}
void LCDVToplTRK::SetUp(Int_t tkid, Double_t* tkpar, Double_t* e_tkpar,
Double_t chrg,
Double_t mfld, Int_t mc_ind, const TVector3& ipv) {
m_track_id=tkid;
Int_t i=0;
for(i=0; i < 5 ; i++) {
m_trackPar[i] = tkpar[i];
}
Int_t k=0;
for(k=0; k < 15 ; k++) {
m_trackErrorMatrix[k] = e_tkpar[k];
}
m_charge =chrg;
m_magneticfield=mfld;
m_index =mc_ind;
Double_t s1i=TMath::Sqrt(GetErrorTsi(0.0 ,ipv));
Double_t s2i=TMath::Sqrt(GetErrorEta(0.0 ,ipv));
Double_t s1m=TMath::Sqrt(GetErrorTsi(m_rmax/2.0,ipv));
Double_t s2m=TMath::Sqrt(GetErrorEta(m_rmax/2.0,ipv));
Double_t s1f=TMath::Sqrt(GetErrorTsi(m_rmax ,ipv));
Double_t s2f=TMath::Sqrt(GetErrorEta(m_rmax ,ipv));
m_err_tsi[0] = s1i;
m_err_tsi[2] = (s1f + s1i)/2.0 - s1m;
m_err_tsi[1] = (s1m - s1i) - m_err_tsi[2];
m_err_eta[0] = s2i;
m_err_eta[2] = (s2f + s2i)/2.0 - s2m;
m_err_eta[1] = (s2m - s2i) - m_err_eta[2];
}
Double_t LCDVToplTRK::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 LCDVToplTRK::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 LCDVToplTRK::GetPositionZ(Double_t l) {
Double_t z0 = m_trackPar[3];
Double_t tanl = m_trackPar[4];
return z0 + l*tanl;
}
TVector3 LCDVToplTRK::GetPositionVector(Double_t l) {
return TVector3(GetPositionX(l),GetPositionY(l),GetPositionZ(l));
}
Double_t LCDVToplTRK::GetMomentumPx(Double_t l) {
Double_t phi0 = m_trackPar[1];
Double_t omega= m_trackPar[2];
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 LCDVToplTRK::GetMomentumPy(Double_t l) {
Double_t phi0 = m_trackPar[1];
Double_t omega= m_trackPar[2];
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 LCDVToplTRK::GetMomentumPz(Double_t l) {
Double_t omega= m_trackPar[2];
Double_t tanl = m_trackPar[4];
Double_t pt = m_magneticfield/333.567/TMath::Abs(omega);
return pt*tanl;
}
TVector3 LCDVToplTRK::GetMomentumVector(Double_t l) {
return TVector3(GetMomentumPx(l),GetMomentumPy(l),GetMomentumPz(l));
}
Double_t LCDVToplTRK::GetTsi(Double_t l, const TVector3& ipv) {
Double_t bz = GetMagneticField();
Double_t phi0 = GetTrackParameter(1);
Double_t omega = GetTrackParameter(2);
Double_t phi = phi0 + omega*l;
TVector3 n(TMath::Sin(phi),-TMath::Cos(phi),0.0);
if (bz < 0) { n*=-1; }
TVector3 dx=GetPositionVector(l) - ipv;
return n*dx;
}
Double_t LCDVToplTRK::GetEta(Double_t l, const TVector3& ipv) {
Double_t bz = GetMagneticField();
Double_t phi0 = GetTrackParameter(1);
Double_t omega = GetTrackParameter(2);
Double_t tanl = GetTrackParameter(4);
Double_t cosl = 1.0/TMath::Sqrt(1.0 + tanl*tanl);
Double_t phi = phi0 + omega*l;
TVector3 q(-TMath::Cos(phi)*tanl*cosl,-TMath::Sin(phi)*tanl*cosl,cosl);
if (bz < 0) { q *= -1; }
TVector3 dx= GetPositionVector(l) - ipv;
return q*dx;
}
Double_t LCDVToplTRK::GetTrackChiSqContrib(const TVector3& xvrt,
const TVector3& ipv,
Int_t f_debug) {
Double_t l=0;
TVector2 x0(GetPositionX(l)-ipv[0],GetPositionY(l)-ipv[1]);
Double_t omega = GetTrackParameter(2);
Double_t phi = GetTrackParameter(1) - TMath::Pi()/2.0;
TVector2 x0p = x0.Rotate(-phi);
TVector2 x0pr = (xvrt-ipv).XYvector().Rotate(-phi);
Double_t y0pr = x0pr.Y() - x0p.Y();
Double_t lp = y0pr;
Double_t z0p = GetPositionZ(l);
Double_t tanl = GetTrackParameter(4);
Double_t s1s = GetErrorTsiApprox(lp);
Double_t s2s = GetErrorEtaApprox(lp)*(1.0 + tanl*tanl);
Double_t dx0pr = x0pr.X() - (x0p.X() - omega/2.0*y0pr*y0pr);
Double_t dz0pr = xvrt[2] - (z0p + tanl *lp );
if (f_debug) {
printf(" x0pvrt =%f z0pvrt=%fn",x0pr.X(),xvrt[2]);
printf(" x0ptk =%f z0ptk =%fn",
x0p.X() - omega/2.0*y0pr*y0pr,z0p+tanl*lp);
printf(" x0p.X()=%f y0pr=%f omega=%f z0p=%f tanl=%fn",
x0p.X(),y0pr,omega,z0p,tanl);
}
return dx0pr/s1s*dx0pr + dz0pr/s2s*dz0pr;
}
Double_t LCDVToplTRK::GetTrackProbabirity(const TVector3& xvrt,
const TVector3& ipv) {
Double_t u = GetTrackChiSqContrib(xvrt,ipv);
Double_t prfu = 0.0;
// protect against floating overflow
if (u < 100.0) {
prfu = TMath::Exp(-0.5*u);
}
return prfu;
}
Double_t LCDVToplTRK::GetDistance(Double_t l, const TVector3& xvrt) {
return (GetPositionVector(l) - xvrt).Mag();
}
Double_t LCDVToplTRK::GetDistance(Double_t l, Double_t x, Double_t y, Double_t z) {
TVector3 x0(x,y,z);
return GetDistance(l,x0);
}
Double_t LCDVToplTRK::GetDistance2D(Double_t l, Double_t x, Double_t y) {
Double_t x0=GetPositionX(l);
Double_t y0=GetPositionY(l);
return TMath::Sqrt((x0 - x)*(x0 - x) + (y0 - y)*(y0 - y));
}
Double_t LCDVToplTRK::GetTrdi(const TVector3& ipv, const TVector3& xvrt) {
Double_t l0=GetLPOCA2(xvrt);
Double_t l=GetLminTrdi(l0,ipv,xvrt);
return GetTrdi(l,ipv,xvrt);
}
Double_t LCDVToplTRK::GetTrdi(Double_t l, const TVector3& ipv, const TVector3& xvrt) {
TVector3 ipv2xvrt=xvrt - ipv;
TVector3 ipv2xtrk=GetPositionVector(l) - ipv;
return ipv2xtrk.Perp(ipv2xvrt);
}
Double_t LCDVToplTRK::GetLodi(const TVector3& ipv, const TVector3& xvrt) {
Double_t l0=GetLPOCA2(xvrt);
Double_t l=GetLminTrdi(l0,ipv,xvrt);
return GetLodi(l,ipv,xvrt);
}
Double_t LCDVToplTRK::GetLodi(Double_t l, const TVector3& ipv, const TVector3& xvrt) {
TVector3 ipv2xvrt=xvrt - ipv;
TVector3 ipv2xtrk=GetPositionVector(l) - ipv;
Double_t lodi=0;
Double_t aipv2xvrt=ipv2xvrt.Mag();
if (aipv2xvrt > 0) lodi=(ipv2xvrt*ipv2xtrk)/aipv2xvrt;
return lodi;
}
void LCDVToplTRK::GetTrdiLodi(Double_t l, const TVector3& ipv, const TVector3& xvt,
Double_t* trdi, Double_t* lodi) {
*lodi=0;
*trdi=0;
TVector3 xtk=GetPositionVector(l);
TVector3 ptk=GetMomentumVector(l);
TVector3 vx=xvt - ipv;
//Double_t anta=vx.Angle(ptk);
TVector3 it2ip=xtk-ipv;
Double_t aa=it2ip*vx;
Double_t dd=it2ip*ptk;
Double_t bb=ptk*vx;
Double_t cc=-vx.Mag2();
Double_t ee=ptk.Mag2();
Double_t lambda=0;
Double_t mu =0;
if (cc+bb*bb/ee == 0.0) { // no solution, track, axis parallel.
*lodi =-1000.0;
lambda= 0.0;
mu = -dd/ee;
} else {
lambda = (dd*bb/ee - aa)/(cc + bb*bb/ee);
mu = (bb*lambda - dd)/ee;
}
// xu = point of normal hit on vertex axis
// xt = point of normal hit on track vector
TVector3 xu=ipv+lambda*vx ;
TVector3 xt=xtk+mu *ptk;
TVector3 xu2ip=xu - ipv;
TVector3 xu2xt=xu - xt;
*lodi = xu2ip.Mag();
//if we're behing IP...
if (xu2ip*vx < 0.0) { *lodi=-*lodi; }
*trdi=xu2xt.Mag();
if (cc+bb*bb/ee == 0.0) {//no solution, track, axis parallel: reset lodi
*lodi=-1000.0;
}
}
void LCDVToplTRK::GetTrdiLodiAnta(Double_t l, const TVector3& ipv, const TVector3& xvt,
Double_t* trdi,
Double_t* lodi,
Double_t* anta) {
*lodi=0;
*trdi=0;
*anta=0;
TVector3 xtk=GetPositionVector(l);
TVector3 ptk=GetMomentumVector(l);
TVector3 vx=xvt - ipv;
*anta=vx.Angle(ptk);
TVector3 it2ip=xtk-ipv;
Double_t aa=it2ip*vx;
Double_t dd=it2ip*ptk;
Double_t bb=ptk*vx;
Double_t cc=-vx.Mag2();
Double_t ee=ptk.Mag2();
Double_t lambda=0;
Double_t mu =0;
if (cc+bb*bb/ee == 0.0) { // no solution, track, axis parallel.
*lodi =-1000.0;
lambda= 0.0;
mu = -dd/ee;
} else {
lambda = (dd*bb/ee - aa)/(cc + bb*bb/ee);
mu = (bb*lambda - dd)/ee;
}
// xu = point of normal hit on vertex axis
// xt = point of normal hit on track vector
TVector3 xu=ipv+lambda*vx ;
TVector3 xt=xtk+mu *ptk;
TVector3 xu2ip=xu - ipv;
TVector3 xu2xt=xu - xt;
*lodi = xu2ip.Mag();
//if we're behing IP...
if (xu2ip*vx < 0.0) { *lodi=-*lodi; }
*trdi=xu2xt.Mag();
if (cc+bb*bb/ee == 0.0) {//no solution, track, axis parallel: reset lodi
*lodi=-1000.0;
}
}
Double_t LCDVToplTRK::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 LCDVToplTRK::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 LCDVToplTRK::GetErrorMatrixPosition(Double_t l,TMatrixD& var2) {
TMatrixD var0(5,5); GetErrorMatrix(var0);
TMatrixD dxdaT(5,3); GetdxdaT(l,dxdaT);
TMatrixD var1(var0,TMatrixD::kMult,dxdaT);
TMatrixD dxda(TMatrixD::kTransposed,dxdaT);
var2.Mult(dxda,var1);
}
void LCDVToplTRK::GetErrorMatrixMomentum(Double_t l,TMatrixD& var2) {
TMatrixD var0(5,5); GetErrorMatrix(var0);
TMatrixD dpdaT(5,3); GetdpdaT(l,dpdaT);
TMatrixD var1(var0,TMatrixD::kMult,dpdaT);
TMatrixD dpda(TMatrixD::kTransposed,dpdaT);
var2.Mult(dpda,var1);
}
Double_t LCDVToplTRK::GetErrorTsi(Double_t l, const TVector3& ipv) {
Double_t bz = GetMagneticField();
if (bz != 0) bz /= TMath::Abs(bz);
Double_t d0 = GetTrackParameter(0);
Double_t phi0 = GetTrackParameter(1);
Double_t omega = GetTrackParameter(2);
Double_t omegal = omega*l;
Double_t phi = phi0 + omegal;
Double_t csomegal = TMath::Cos(omegal);
Double_t snomegal = TMath::Sin(omegal);
Double_t csphi = TMath::Cos(phi);
Double_t snphi = TMath::Sin(phi);
TMatrixD dtsida(5,1);
dtsida(0,0)=-bz*csomegal;
dtsida(1,0)=-bz*(ipv[0]*csphi + ipv[1]*snphi);
dtsida(2,0)=-bz*((csomegal - 1.0)/omega/omega + (1.0/omega + d0)*l*snomegal
- ipv[0]*l*csphi - ipv[1]*l*snphi);
dtsida(3,0)= 0.0;
dtsida(4,0)= 0.0;
TMatrixD trkparmerror(5,5); GetErrorMatrix(trkparmerror);
TMatrixD tmp0(trkparmerror,TMatrixD::kMult,dtsida);
TMatrixD tmp1(dtsida,TMatrixD::kTransposeMult,tmp0);
return tmp1(0,0);
}
Double_t LCDVToplTRK::GetErrorEta(Double_t l, const TVector3& ipv) {
Double_t bz = GetMagneticField();
if (bz != 0) bz /= TMath::Abs(bz);
Double_t d0 = GetTrackParameter(0);
Double_t phi0 = GetTrackParameter(1);
Double_t omega = GetTrackParameter(2);
Double_t z0 = GetTrackParameter(3);
Double_t tanl = GetTrackParameter(4);
Double_t cosl = 1.0/TMath::Sqrt(1.0 + tanl*tanl);
Double_t sinl = tanl*cosl;
Double_t omegal = omega*l;
Double_t phi = phi0 + omegal;
Double_t csomegal = TMath::Cos(omegal);
Double_t snomegal = TMath::Sin(omegal);
Double_t csphi = TMath::Cos(phi);
Double_t snphi = TMath::Sin(phi);
TMatrixD detadaT(5,1);
detadaT(0,0)=-bz*sinl*snomegal;
detadaT(1,0)= bz*sinl*(-ipv[0]*snphi + ipv[1]*csphi);
detadaT(2,0)= bz*sinl*(snomegal/omega/omega - (1.0/omega + d0)*l*csomegal
- ipv[0]*l*snphi - ipv[1]*l*csphi);
detadaT(3,0)= bz*cosl;
detadaT(4,0)= bz/cosl*(-(1.0/omega + d0)*snomegal + l + ipv[0]*csphi
+ipv[1]*snphi - (z0 - ipv[2])*tanl);
TMatrixD trkparmerror(5,5); GetErrorMatrix(trkparmerror);
TMatrixD tmp0(trkparmerror,TMatrixD::kMult,detadaT);
TMatrixD tmp1(detadaT,TMatrixD::kTransposeMult,tmp0);
return tmp1(0,0);
}
Double_t LCDVToplTRK::GetErrorTsiApprox(Double_t l) {
Double_t ll=TMath::Min(l,m_rmax);
ll=TMath::Max(ll,0.0);
return TMath::Power(m_err_tsi[0] + m_err_tsi[1]*ll + m_err_tsi[2]*ll*ll,2);
}
Double_t LCDVToplTRK::GetErrorEtaApprox(Double_t l) {
Double_t ll=TMath::Min(l,m_rmax);
ll=TMath::Max(ll,0.0);
return TMath::Power(m_err_eta[0] + m_err_eta[1]*ll + m_err_eta[2]*ll*ll,2);
}
Double_t LCDVToplTRK::GetErrorDistance(Double_t l, const TVector3& xvrt) {
return GetErrorDistance(l,xvrt[0],xvrt[1],xvrt[2]);
}
Double_t LCDVToplTRK::GetErrorDistance(Double_t l,
Double_t x, Double_t y, Double_t z) {
TMatrixD dddaT(5,1); GetdddaT(l,x,y,z,dddaT);
TMatrixD em(5,5); GetErrorMatrix(em);
TMatrixD tmp0(em,TMatrixD::kMult,dddaT);
TMatrixD tmp1(dddaT,TMatrixD::kTransposeMult,tmp0);
return tmp1(0,0);
}
Double_t LCDVToplTRK::GetErrorDistance(Double_t l, const TVector3& x_from,
TMatrixD& ex_from) {
Double_t svrt=GetErrorDistance(l,x_from);
TVector3 xvrt=GetPositionVector(l);
TVector3 dx =xvrt - x_from;
Double_t d =dx.Mag();
TMatrixD dddxipT(3,1);
if (d > 0) {
dddxipT(0,0)=-dx[0]/d;
dddxipT(1,0)=-dx[1]/d;
dddxipT(2,0)=-dx[2]/d;
} else {
dddxipT(0,0)=0.0;
dddxipT(1,0)=0.0;
dddxipT(2,0)=0.0;
}
TMatrixD tmp0(ex_from,TMatrixD::kMult,dddxipT);
TMatrixD tmp1(dddxipT,TMatrixD::kTransposeMult,tmp0);
return tmp1(0,0)+svrt;
}
Double_t LCDVToplTRK::GetErrorDistance(Double_t l,
Double_t x, Double_t y, Double_t z,
TMatrixD& ex_from) {
TVector3 xv(x,y,z);
return GetErrorDistance(l,xv,ex_from);
}
Double_t LCDVToplTRK::GetErrorDistance2D(Double_t l, Double_t x, Double_t y) {
TMatrixD dddaT(5,1); GetdddaT2D(l,x,y,dddaT);
TMatrixD em(5,5); GetErrorMatrix(em);
TMatrixD tmp0(em,TMatrixD::kMult,dddaT);
TMatrixD tmp1(dddaT,TMatrixD::kTransposeMult,tmp0);
return tmp1(0,0);
}
void LCDVToplTRK::GetdxdaT(Double_t l, TMatrixD& dxdaT) {
Double_t d0 = m_trackPar[0];
Double_t phi0 = m_trackPar[1];
Double_t omega = m_trackPar[2];
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);
dxdaT(0,0) =-snphi0;
dxdaT(1,0) = 1.0/omega*csphi - (1.0/omega + d0)*csphi0;
dxdaT(2,0) = (1.0/omega*(-snphi + snphi0) + csphi*l)/omega;
dxdaT(3,0) = 0.0;
dxdaT(4,0) = 0.0;
dxdaT(0,1) = csphi0;
dxdaT(1,1) = 1.0/omega*snphi - (1.0/omega + d0)*snphi0;
dxdaT(2,1) = (1.0/omega*( csphi - csphi0) + snphi*l)/omega;
dxdaT(3,1) = 0.0;
dxdaT(4,1) = 0.0;
dxdaT(0,2) = 0.0;
dxdaT(1,2) = 0.0;
dxdaT(2,2) = 0.0;
dxdaT(3,2) = 1.0;
dxdaT(4,2) = l;
}
void LCDVToplTRK::GetdpdaT(Double_t l, TMatrixD& dpdaT) {
Double_t phi0 = m_trackPar[1];
Double_t omega = m_trackPar[2];
Double_t tanl = m_trackPar[4];
Double_t phi = phi0 + omega*l;
Double_t csphi = TMath::Cos(phi);
Double_t snphi = TMath::Sin(phi);
Double_t r2p = TMath::Abs(m_magneticfield)/333.567;
dpdaT(0,0) = 0.0;
dpdaT(1,0) =-r2p/omega*snphi;
dpdaT(2,0) =-r2p/omega/omega*csphi - r2p/omega*snphi*l;
dpdaT(3,0) = 0.0;
dpdaT(4,0) = 0.0;
dpdaT(0,1) = 0.0;
dpdaT(1,1) = r2p/omega*csphi;
dpdaT(2,1) =-r2p/omega/omega*snphi + r2p/omega*csphi*l;
dpdaT(3,1) = 0.0;
dpdaT(4,1) = 0.0;
dpdaT(0,2) = 0.0;
dpdaT(1,2) = 0.0;
dpdaT(2,2) =-r2p/omega/omega*tanl;
dpdaT(3,2) = 0.0;
dpdaT(4,2) = r2p/omega;
}
void LCDVToplTRK::GetdddaT(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(3,1);
dx(0,0) = -(x - xp)/d;
dx(1,0) = -(y - yp)/d;
dx(2,0) = -(z - zp)/d;
TMatrixD dxdaT(5,3); GetdxdaT(l,dxdaT);
ddda.Mult(dxdaT,dx);
}
void LCDVToplTRK::GetdddaT2D(Double_t l,
Double_t x, Double_t y,
TMatrixD& ddda) {
Double_t xp=GetPositionX(l);
Double_t yp=GetPositionY(l);
Double_t d =GetDistance2D(l,x,y);
TMatrixD dx(3,1);
dx(0,0) = -(x - xp)/d;
dx(1,0) = -(y - yp)/d;
dx(2,0) = 0.0;
TMatrixD dxdaT(5,3); GetdxdaT(l,dxdaT);
ddda.Mult(dxdaT,dx);
}
void LCDVToplTRK::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 LCDVToplTRK::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;
}
TVector3 LCDVToplTRK::GetdxdlVector(Double_t l) {
Double_t phi0=GetTrackParameter(1);
Double_t omega=GetTrackParameter(2);
Double_t tanl=GetTrackParameter(4);
Double_t phi=phi0+omega*l;
return TVector3(TMath::Cos(phi),TMath::Sin(phi),tanl);
}
Double_t LCDVToplTRK::GetLPOCA2(const TVector3& x_given) {
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[2]=0.0;
TVector3 nxxc=xxc.Unit();
TVector3 nxc=-xc.Unit();
TVector3 nxcXnxxc=nxc.Cross(nxxc);
Double_t nxcXnxxcZ=nxcXnxxc[2];
if (nxcXnxxcZ < -1.0) { nxcXnxxcZ=-1.0; }
else if (nxcXnxxcZ > +1.0) { nxcXnxxcZ=+1.0; }
return TMath::ASin(nxcXnxxcZ)/omega;
}
// Get Closest Approach Point to beam axis...
Double_t LCDVToplTRK::CalcPOCA3(Double_t l0,const TVector3& x_given) {
const Int_t MAXTRY=50;
const Double_t DLEN2TH=1E-10;
Double_t l=l0;
Double_t l_best=0.0,ddlen2dl_best=0.0,d2dlen2d2l_best=0.0;
Double_t dlen2_c;
Double_t dlen2_p=1e10;
Double_t a;
Double_t f=1.0;
for (Int_t i=0; i < MAXTRY ; i++) {
dlen2_c=dlen2(l,x_given);
if ( TMath::Abs(dlen2_c - dlen2_p) < DLEN2TH && i != 0) {
break;
} else if ( i >= MAXTRY-1 ) {
l = l_best;
break;
} else if ( dlen2_c > dlen2_p && i != 0 ) {
f /= 2.0;
l = l_best;
l -= ddlen2dl_best/d2dlen2d2l_best*f;
} 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 /= 2.0;
l = l_best;
}
l -= ddlen2dl_best/d2dlen2d2l_best*f;
}
}
return l;
}
Double_t LCDVToplTRK::CalcPOCA3(const TVector3& x_given) {
Double_t l=GetLPOCA2(x_given);
TVector3 a=x_given-GetPositionVector(l);
TVector3 b=GetMomentumVector(l);
Double_t tanl=GetTrackParameter(4);
l += a*b/b.Mag()*TMath::Sqrt(1.0/(1.0 + tanl*tanl));
return CalcPOCA3(l,x_given);
}
// void LCDVToplTRK::AddVertex(TObject* vtx) {
void LCDVToplTRK::AddVertex(LCDVToplVRT* vtx) {
if (m_nvrt < 99) {
vertex_list[m_nvrt++]=vtx;
} else {
fprintf(stderr,"LCDVToplTRK::AddVertex # of vertices is over 100!n");
}
}
// void LCDVToplTRK::AddVertexAt(TObject* vtx,Int_t i) {
void LCDVToplTRK::AddVertexAt(LCDVToplVRT* vtx,Int_t i) {
if (i < m_nvrt+1) {
if (vertex_list[i] == 0) m_nvrt++;
vertex_list[i]=vtx;
} else {
fprintf(stderr,
"LCDVToplTRK::AddVertex vertex index is over 100! idx=%dn nv=%d",
i,m_nvrt);
}
}
// void LCDVToplTRK::RemoveVertex(TObject* vrt){
void LCDVToplTRK::RemoveVertex(LCDVToplVRT* vrt){
for (Int_t i=0 ; i < m_nvrt ; i++) {
if(vertex_list[i] == vrt) {
vertex_list[i] = 0;
}
}
CompressVerticesList();
}
void LCDVToplTRK::RemoveVertexAt(Int_t ivrt){
if (ivrt < m_nvrt) {
vertex_list[ivrt]=0;
CompressVerticesList();
} else {
fprintf(stderr,
"LCDVToplTRK::RemoveVertexAt ver. idx is over nv! idx=%d nv=%dn",
ivrt,m_nvrt);
}
}
void LCDVToplTRK::CompressVerticesList() {
Int_t nzero=0;
for (Int_t i=0 ; i < m_nvrt ; i++) {
if(vertex_list[i] == 0) {
nzero++;
continue;
}
if (nzero == 0) continue;
vertex_list[i-nzero]=vertex_list[i];
}
m_nvrt -= nzero;
}
Double_t LCDVToplTRK::dlen2(Double_t l, const TVector3& x_given) {
return (x_given - GetPositionVector(l)).Mag2();
}
Double_t LCDVToplTRK::ddlen2dl(Double_t l,const 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[0];
Double_t y0p =x_given[1];
Double_t z0p =x_given[2];
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 LCDVToplTRK::d2dlen2d2l(Double_t l,const TVector3& x_given)
{
Double_t d0 =m_trackPar[0];
Double_t phi0 =m_trackPar[1];
Double_t omega=m_trackPar[2];
Double_t tanl =m_trackPar[4];
Double_t x0p =x_given[0];
Double_t y0p =x_given[1];
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);
}
// Get Closest Approach Point to given axis...
Double_t LCDVToplTRK::GetLminTrdiApprox(const TVector3& xfrom,const TVector3& xto) {
Double_t l0=GetLPOCA2(xto);
return GetLminTrdiApprox(l0,xfrom,xto);
}
Double_t LCDVToplTRK::GetLminTrdiApprox(Double_t l,
const TVector3& ipv,
const TVector3& xvrt) {
TVector3 ipv2xvrt=xvrt - ipv;
TVector3 xtk=GetPositionVector(l);
TVector3 ptk=GetMomentumVector(l);
TVector3 ipv2xtk=xtk - ipv;
Double_t a=ipv2xtk*ipv2xvrt;
Double_t d=ipv2xtk*ptk;
Double_t b=ptk*ipv2xvrt;
Double_t c=-ipv2xvrt.Mag2();
Double_t e= ptk.Mag2();
Double_t lambda=0;
Double_t mu=0;
if (c+b*b/e == 0.0) { // no solution, track, axis parallel.
lambda= 0.0;
mu = -d/e;
} else {
lambda = (d*b/e - a)/(c + b*b/e);
mu = (b*lambda - d)/e;
}
return l+mu*ptk.Pt();
}
Double_t LCDVToplTRK::GetLminTrdi(const TVector3& xfrom,const TVector3& xto) {
Double_t l0=GetLPOCA2(xto);
return GetLminTrdi(l0,xfrom,xto);
}
Double_t LCDVToplTRK::GetLminTrdi(Double_t l0,const TVector3& xfrom,const TVector3& xto) {
const Int_t MAXTRY=50;
const Double_t TRDITH=1E-10;
Double_t l=l0;
Double_t l_best=0.0;
Double_t dtrdidl_best=0.0;
Double_t d2trdid2l_best=0.0;
Double_t trdi_c=1e10;
Double_t trdi_p=1e10;
Double_t a=0.0;
Double_t f=1.0;
for (Int_t i=0; i < MAXTRY ; i++) {
trdi_c=TMath::Power(GetTrdi(l,xfrom,xto),2);
if ( TMath::Abs(trdi_c - trdi_p) < TRDITH && i != 0) {
break;
} else if ( i >= MAXTRY-1 ) {
l = l_best;
break;
} else if ( trdi_c > trdi_p && i != 0 ) {
f /= 2.0;
l = l_best;
l -= dtrdidl_best/d2trdid2l_best*f;
} else {
a= Getd2Trdi2d2l(l,xfrom,xto);
if (a != 0) {
trdi_p = trdi_c;
f = 1.0;
l_best = l;
dtrdidl_best = GetdTrdi2dl(l,xfrom,xto);
d2trdid2l_best = a;
} else {
f /= 2.0;
l = l_best;
}
l -= dtrdidl_best/d2trdid2l_best*f;
}
}
return l;
}
Double_t LCDVToplTRK::GetdTrdi2dl(Double_t l, const TVector3& xfrom, const TVector3& xto) {
TVector3 n = xto - xfrom; n.Unit();
TVector3 dxdl = GetdxdlVector(l);
TVector3 x = GetPositionVector(l);
TVector3 xfrom2x = x - xfrom;
return 2*(xfrom2x*dxdl) - 2*(x*n)*(n*dxdl);
}
Double_t LCDVToplTRK::Getd2Trdi2d2l(Double_t l, const TVector3& xfrom, const TVector3& xto) {
Double_t d0 = GetTrackParameter(0);
Double_t phi0 = GetTrackParameter(1);
Double_t omega = GetTrackParameter(2);
Double_t tanl = GetTrackParameter(4);
Double_t omegal = omega*l;
Double_t phi = phi0 + omegal;
Double_t csphi = TMath::Cos(phi);
Double_t snphi = TMath::Sin(phi);
Double_t csomegal= TMath::Cos(omegal);
Double_t snomegal= TMath::Sin(omegal);
TVector3 n = xto - xfrom; n.Unit();
TVector3 x = GetPositionVector(l);
TVector3 xfrom2x = x - xfrom;
return 2*(d0 + 1.0/omega)*omega*csomegal + 2*xfrom[0]*omega*snphi
- 2*xfrom[1]*omega*csphi + 2*tanl*tanl
- 2*n[0]*n[0]*csphi*csphi + 2*n[0]*n[0]*omega*xfrom2x[0]*snphi
- 2*n[1]*n[1]*snphi*snphi - 2*n[1]*n[1]*omega*xfrom2x[1]*csphi
- 2*n[2]*n[2]*tanl*tanl
+ 2*n[0]*n[1]*(d0 + 1.0/omega)*omega*snomegal
+ 2*n[0]*n[1]*(-xfrom[1]*snphi + xfrom[0]*snphi)
- 2*n[0]*n[2]*tanl*csphi + 2*n[0]*n[2]*xfrom2x[2]*omega*snphi
- 2*n[0]*n[2]*snphi*tanl
- 2*n[1]*n[2]*tanl*snphi - 2*n[1]*n[2]*xfrom2x[2]*omega*csphi
- 2*n[1]*n[2]*snphi*tanl;
}
Double_t LCDVToplTRK::GetSignedDistance(const TVector3& xfrom,
const TVector3& pjet) {
Double_t l=CalcPOCA3(xfrom);
return GetSignedDistance(l,xfrom,pjet);
}
Double_t LCDVToplTRK::GetSignedDistance(Double_t l,
const TVector3& xfrom,
const TVector3& pjet) {
TVector3 dx(GetPositionVector(l) - xfrom);
Double_t dlen=dx.Mag();
return (dx*pjet > 0) ? dlen : -dlen;
}
Double_t LCDVToplTRK::GetDistanceNormByError(const TVector3& xfrom) {
Double_t l=CalcPOCA3(xfrom);
return GetDistanceNormByError(l,xfrom);
}
Double_t LCDVToplTRK::GetDistanceNormByError(Double_t l, const TVector3& xfrom) {
Double_t dist=GetDistance(l,xfrom);
Double_t sdst=TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom)));
return dist/sdst;
}
Double_t LCDVToplTRK::GetDistanceNormByError(const TVector3& xfrom,
TMatrixD& ex_from) {
Double_t l=CalcPOCA3(xfrom);
return GetDistanceNormByError(l,xfrom,ex_from);
}
Double_t LCDVToplTRK::GetDistanceNormByError(Double_t l, const TVector3& xfrom,
TMatrixD& ex_from) {
Double_t dist=GetDistance(l,xfrom);
Double_t sdst=
TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom,ex_from)));
return dist/sdst;
}
Double_t LCDVToplTRK::GetSignedDistanceNormByError(const TVector3& xfrom,
const TVector3& pjet) {
Double_t l=CalcPOCA3(xfrom);
return GetSignedDistanceNormByError(l,xfrom,pjet);
}
Double_t LCDVToplTRK::GetSignedDistanceNormByError(Double_t l,
const TVector3& xfrom,
const TVector3& pjet) {
Double_t dist=GetSignedDistance(l,xfrom,pjet);
Double_t sdst=TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom)));
return dist/sdst;
}
Double_t LCDVToplTRK::GetSignedDistanceNormByError(const TVector3& xfrom,
const TVector3& pjet,
TMatrixD& ex_from) {
Double_t l=CalcPOCA3(xfrom);
return GetSignedDistanceNormByError(l,xfrom,pjet,ex_from);
}
Double_t LCDVToplTRK::GetSignedDistanceNormByError(Double_t l, const TVector3& xfrom,
const TVector3& pjet,
TMatrixD& ex_from) {
Double_t dist=GetSignedDistance(l,xfrom,pjet);
Double_t sdst=
TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom,ex_from)));
return dist/sdst;
}
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.