diff --git a/algo/ca/core/tracking/CaCloneMerger.cxx b/algo/ca/core/tracking/CaCloneMerger.cxx
index 561a0f5ea8dfab75ef482fe7434000ff09e22028..9332f72dc5b83bca274c9ef521922a69ab3bd0b9 100644
--- a/algo/ca/core/tracking/CaCloneMerger.cxx
+++ b/algo/ca/core/tracking/CaCloneMerger.cxx
@@ -80,12 +80,12 @@ void CloneMerger::Exec(const ca::InputData& input, WindowData& wData)
isDownstreamNeighbour[iTr] = false;
}
- ca::TrackFit fitB;
+ ca::TrackFit<fvec> fitB;
fitB.SetParticleMass(fDefaultMass);
fitB.SetMask(fmask::One());
fitB.SetQp0(fvec(0.));
- ca::TrackFit fitF;
+ ca::TrackFit<fvec> fitF;
fitF.SetParticleMass(fDefaultMass);
fitF.SetMask(fmask::One());
fitF.SetQp0(fvec(0.));
@@ -263,8 +263,8 @@ void CloneMerger::InvertCholesky(fvec a[15])
fvec smallval(1.e-12);
uud = iif(uud >= smallval, uud, smallval);
- d[i] = uud / abs(uud);
- u[i][i] = sqrt(abs(uud));
+ d[i] = uud / utils::fabs(uud);
+ u[i][i] = sqrt(utils::fabs(uud));
for (int j = i + 1; j < 5; j++) {
uud = 0.f;
diff --git a/algo/ca/core/tracking/CaTrackExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx
index a92a5b0193576043ea313ce7420c71b3cab16e4c..2b7e606765e5077fab15e44bdcb7758c48e0380c 100644
--- a/algo/ca/core/tracking/CaTrackExtender.cxx
+++ b/algo/ca/core/tracking/CaTrackExtender.cxx
@@ -42,7 +42,7 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
{
CBMCA_DEBUG_ASSERT(t.NHits >= 3);
- ca::TrackFit fit;
+ ca::TrackFit<fvec> fit;
fit.SetParticleMass(fDefaultMass);
fit.SetMask(fmask::One());
fit.SetTrack(Tout);
@@ -119,7 +119,7 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
const ca::Station<fvec>& sta = fParameters.GetStation(ista);
fit.Extrapolate(hit.Z(), fld);
- fit.FilterHit(sta, hit);
+ fit.FilterHit(hit, fmask(sta.timeInfo));
auto radThick = fParameters.GetMaterialThickness(ista, T.X(), T.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec(upstream ? 1. : -1.));
@@ -153,7 +153,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
Vector<ca::HitIndex_t> newHits{"ca::TrackExtender::newHits"};
newHits.reserve(fParameters.GetNstationsActive());
- ca::TrackFit fit;
+ ca::TrackFit<fvec> fit;
fit.SetParticleMass(fDefaultMass);
fit.SetMask(fmask::One());
fit.SetTrack(Tout);
@@ -215,8 +215,8 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
const auto& grid = frWData->Grid(ista);
ca::GridArea area(grid, tr.X()[0], tr.Y()[0],
- (sqrt(pickGather * tr.C00()) + grid.GetMaxRangeX() + maxDZ * abs(tr.Tx()))[0],
- (sqrt(pickGather * tr.C11()) + grid.GetMaxRangeY() + maxDZ * abs(tr.Ty()))[0]);
+ (sqrt(pickGather * tr.C00()) + grid.GetMaxRangeX() + maxDZ * utils::fabs(tr.Tx()))[0],
+ (sqrt(pickGather * tr.C11()) + grid.GetMaxRangeY() + maxDZ * utils::fabs(tr.Ty()))[0]);
if (fParameters.DevIsIgnoreHitSearchAreas()) {
area.DoLoopOverEntireGrid();
@@ -271,7 +271,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
newHits.push_back(hit.Id());
fit.Extrapolate(hit.Z(), fld);
- fit.FilterHit(sta, hit);
+ fit.FilterHit(hit, fmask(sta.timeInfo));
auto radThick = fParameters.GetMaterialThickness(ista, tr.X(), tr.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, upstream ? fvec(1.f) : fvec(-1.f));
diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx
index 3b7328c2b69adbaf801a5950c0f95c30661ad35f..3666eef74930a847ddf9a15c5109ef802cd4d449 100644
--- a/algo/ca/core/tracking/CaTrackFit.cxx
+++ b/algo/ca/core/tracking/CaTrackFit.cxx
@@ -12,13 +12,12 @@
namespace cbm::algo::ca
{
- typedef const fvec cnst;
-
- void TrackFit::Filter1d(const ca::MeasurementU<fvec>& m)
+ template<typename DataT>
+ void TrackFit<DataT>::Filter1d(const ca::MeasurementU<DataT>& m)
{
- fvec zeta, HCH;
- fvec F0, F1, F2, F3, F4, F5, F6;
- fvec K1, K2, K3, K4, K5, K6;
+ DataT zeta, HCH;
+ DataT F0, F1, F2, F3, F4, F5, F6;
+ DataT K1, K2, K3, K4, K5, K6;
zeta = m.CosPhi() * fTr.X() + m.SinPhi() * fTr.Y() - m.U();
@@ -34,17 +33,17 @@ namespace cbm::algo::ca
F5 = m.CosPhi() * fTr.C50() + m.SinPhi() * fTr.C51();
F6 = m.CosPhi() * fTr.C60() + m.SinPhi() * fTr.C61();
- constexpr bool doProtect = !std::is_same<fscal, double>::value;
+ constexpr bool doProtect = !std::is_same<DataTscal, double>::value;
- const fmask maskDoFilter = doProtect ? (HCH < m.Du2() * 16.f) : fmask::One();
+ const DataTmask maskDoFilter = doProtect ? (HCH < m.Du2() * 16.f) : DataTmask(true);
// correction to HCH is needed for the case when sigma2 is so small
// with respect to HCH that it disappears due to the roundoff error
//
- fvec wi = doProtect ?(fMaskF / (m.Du2() + fvec(1.0000001) * HCH)) : (fMaskF / (m.Du2() + HCH));
- fvec zetawi = fMaskF * zeta / (iif(maskDoFilter, m.Du2(), fvec::Zero()) + HCH);
+ DataT wi = doProtect ? (fMaskF / (m.Du2() + DataT(1.0000001) * HCH)) : (fMaskF / (m.Du2() + HCH));
+ DataT zetawi = fMaskF * zeta / (utils::iif(maskDoFilter, m.Du2(), DataT(0.)) + HCH);
- wi = iif(m.Du2() > fvec::Zero(), wi, fvec::Zero());
+ wi = utils::iif(m.Du2() > DataT(0.), wi, DataT(0.));
fTr.ChiSq() += m.Ndf() * zeta * zeta * wi;
fTr.Ndf() += m.Ndf() * fMaskF;
@@ -100,45 +99,44 @@ namespace cbm::algo::ca
fTr.C66() -= K6 * F6;
}
- void TrackFit::FilterTime(fvec t, fvec dt2, fvec timeInfo)
+ template<typename DataT>
+ void TrackFit<DataT>::FilterTime(DataT t, DataT dt2, const DataTmask& timeInfo)
{
// filter track with a time measurement
// F = CH'
- fvec F0 = fTr.C50();
- fvec F1 = fTr.C51();
- fvec F2 = fTr.C52();
- fvec F3 = fTr.C53();
- fvec F4 = fTr.C54();
- fvec F5 = fTr.C55();
- fvec F6 = fTr.C65();
-
- fvec HCH = fTr.C55();
+ DataT F0 = fTr.C50();
+ DataT F1 = fTr.C51();
+ DataT F2 = fTr.C52();
+ DataT F3 = fTr.C53();
+ DataT F4 = fTr.C54();
+ DataT F5 = fTr.C55();
+ DataT F6 = fTr.C65();
- fvec w = fMaskF;
+ DataT HCH = fTr.C55();
- w.setZero(timeInfo <= fvec::Zero());
+ DataT w = utils::iif(timeInfo, fMaskF, DataT(0.));
// when dt0 is much smaller than current time error,
// set track time exactly to the measurement value without filtering
// it helps to keep the initial time errors reasonably small
// the calculations in the covariance matrix are not affected
- const fmask maskDoFilter = (HCH < dt2 * 16.f);
+ const DataTmask maskDoFilter = (HCH < dt2 * 16.f);
- fvec wi = w / (dt2 + fvec(1.0000001) * HCH);
- fvec zeta = fTr.Time() - t;
- fvec zetawi = w * zeta / (iif(maskDoFilter, dt2, fvec::Zero()) + HCH);
+ DataT wi = w / (dt2 + DataT(1.0000001) * HCH);
+ DataT zeta = fTr.Time() - t;
+ DataT zetawi = w * zeta / (utils::iif(maskDoFilter, dt2, DataT(0.)) + HCH);
- fTr.ChiSqTime()(maskDoFilter) += zeta * zeta * wi;
+ fTr.ChiSqTime() += utils::iif(maskDoFilter, zeta * zeta * wi, DataT(0.));
fTr.NdfTime() += w;
- fvec K1 = F1 * wi;
- fvec K2 = F2 * wi;
- fvec K3 = F3 * wi;
- fvec K4 = F4 * wi;
- fvec K5 = F5 * wi;
- fvec K6 = F6 * wi;
+ DataT K1 = F1 * wi;
+ DataT K2 = F2 * wi;
+ DataT K3 = F3 * wi;
+ DataT K4 = F4 * wi;
+ DataT K5 = F5 * wi;
+ DataT K6 = F6 * wi;
fTr.X() -= F0 * zetawi;
fTr.Y() -= F1 * zetawi;
@@ -185,42 +183,48 @@ namespace cbm::algo::ca
}
- void TrackFit::FilterXY(const ca::MeasurementXy<fvec>& mxy, bool skipUnmeasuredCoordinates)
+ template<typename DataT>
+ void TrackFit<DataT>::FilterXY(const ca::MeasurementXy<DataT>& mxy, bool skipUnmeasuredCoordinates)
{
{
- ca::MeasurementU<fvec> mx;
- mx.SetCosPhi(fvec::One());
- mx.SetSinPhi(fvec::Zero());
+ ca::MeasurementU<DataT> mx;
+ mx.SetCosPhi(DataT(1.));
+ mx.SetSinPhi(DataT(0.));
mx.SetU(mxy.X());
mx.SetDu2(mxy.Dx2());
mx.SetNdf(mxy.NdfX());
- ca::MeasurementU<fvec> mu;
+ ca::MeasurementU<DataT> mu;
mu.SetCosPhi(-mxy.Dxy() / mxy.Dx2());
- mu.SetSinPhi(fvec::One());
+ mu.SetSinPhi(DataT(1.));
mu.SetU(mu.CosPhi() * mxy.X() + mxy.Y());
mu.SetDu2(mxy.Dy2() - mxy.Dxy() * mxy.Dxy() / mxy.Dx2());
mu.SetNdf(mxy.NdfY());
- if (!(skipUnmeasuredCoordinates && mxy.NdfX()[0] == 0)) {
- Filter1d(mx);
+ auto maskOld = fMask;
+ if (skipUnmeasuredCoordinates) {
+ fMask = maskOld & (mxy.NdfX() > DataT(0.));
}
- if (!(skipUnmeasuredCoordinates && mxy.NdfY()[0] == 0)) {
- Filter1d(mu);
+ Filter1d(mx);
+ if (skipUnmeasuredCoordinates) {
+ fMask = maskOld & (mxy.NdfY() > DataT(0.));
}
+ Filter1d(mu);
+ fMask = maskOld;
+
return;
}
//----------------------------------------------------------------------------------------------
// the other way: filter 2 dimensions at once
- cnst TWO(2.);
+ const DataT TWO(2.);
- fvec zeta0, zeta1, S00, S10, S11, si;
- fvec F00, F10, F20, F30, F40, F50, F60;
- fvec F01, F11, F21, F31, F41, F51, F61;
- fvec K00, K10, K20, K30, K40, K50, K60;
- fvec K01, K11, K21, K31, K41, K51, K61;
+ DataT zeta0, zeta1, S00, S10, S11, si;
+ DataT F00, F10, F20, F30, F40, F50, F60;
+ DataT F01, F11, F21, F31, F41, F51, F61;
+ DataT K00, K10, K20, K30, K40, K50, K60;
+ DataT K01, K11, K21, K31, K41, K51, K61;
zeta0 = fTr.X() - mxy.X();
zeta1 = fTr.Y() - mxy.Y();
@@ -246,11 +250,11 @@ namespace cbm::algo::ca
S10 = F10 + mxy.Dxy();
S11 = F11 + mxy.Dy2();
- si = 1.f / (S00 * S11 - S10 * S10);
- fvec S00tmp = S00;
- S00 = si * S11;
- S10 = -si * S10;
- S11 = si * S00tmp;
+ si = 1.f / (S00 * S11 - S10 * S10);
+ DataT S00tmp = S00;
+ S00 = si * S11;
+ S10 = -si * S10;
+ S11 = si * S00tmp;
fTr.ChiSq() += zeta0 * zeta0 * S00 + 2.f * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
fTr.Ndf() += TWO;
@@ -320,23 +324,25 @@ namespace cbm::algo::ca
fTr.C66() -= K60 * F60 + K61 * F61;
}
- void TrackFit::FilterHit(const ca::Station<fvec>& sta, const ca::Hit& hit)
+ template<typename DataT>
+ void TrackFit<DataT>::FilterHit(const ca::Hit& hit, const DataTmask& timeInfo)
{
- ca::MeasurementXy<fvec> m;
+ ca::MeasurementXy<DataT> m;
m.SetDx2(hit.dX2());
m.SetDy2(hit.dY2());
m.SetDxy(hit.dXY());
m.SetX(hit.X());
m.SetY(hit.Y());
- m.SetNdfX(fvec::One());
- m.SetNdfY(fvec::One());
+ m.SetNdfX(DataT(1.));
+ m.SetNdfY(DataT(1.));
FilterXY(m);
- FilterTime(hit.T(), hit.dT2(), sta.timeInfo);
+ FilterTime(hit.T(), hit.dT2(), timeInfo);
}
- void TrackFit::FilterExtrapolatedXY(const ca::MeasurementXy<fvec>& m, fvec extrX, fvec extrY,
- const std::array<fvec, ca::TrackParamV::kNtrackParam>& Jx,
- const std::array<fvec, ca::TrackParamV::kNtrackParam>& Jy)
+ template<typename DataT>
+ void TrackFit<DataT>::FilterExtrapolatedXY(const ca::MeasurementXy<DataT>& m, DataT extrX, DataT extrY,
+ const std::array<DataT, ca::TrackParamBase<DataT>::kNtrackParam>& Jx,
+ const std::array<DataT, ca::TrackParamBase<DataT>::kNtrackParam>& Jy)
{
// add a 2-D measurenent (x,y) at some z, that differs from fTr.GetZ()
// extrX, extrY are extrapolated track parameters at z, Jx, Jy are derivatives of the extrapolation
@@ -344,55 +350,55 @@ namespace cbm::algo::ca
// ! it is assumed that in the track covariance matrix all non-diagonal covariances are 0
// ! except of C10
- TrackParamV& T = fTr;
+ auto& T = fTr;
//zeta0 = T.x + Jx[2]*T.Tx() + Jx[3]*T.Ty() + Jx[4]*T.qp - x;
//zeta1 = T.y + Jy[2]*T.Tx() + Jy[3]*T.Ty() + Jy[4]*T.qp - y;
- fvec zeta0 = extrX - m.X();
- fvec zeta1 = extrY - m.Y();
+ DataT zeta0 = extrX - m.X();
+ DataT zeta1 = extrY - m.Y();
// H = 1 0 Jx[2] Jx[3] Jx[4] 0
// 0 1 Jy[2] Jy[3] Jy[4] 0
// F = CH'
- fvec F00 = T.C00();
- fvec F01 = T.C10();
- fvec F10 = T.C10();
- fvec F11 = T.C11();
- fvec F20 = Jx[2] * T.C22();
- fvec F21 = Jy[2] * T.C22();
- fvec F30 = Jx[3] * T.C33();
- fvec F31 = Jy[3] * T.C33();
- fvec F40 = Jx[4] * T.C44();
- fvec F41 = Jy[4] * T.C44();
+ DataT F00 = T.C00();
+ DataT F01 = T.C10();
+ DataT F10 = T.C10();
+ DataT F11 = T.C11();
+ DataT F20 = Jx[2] * T.C22();
+ DataT F21 = Jy[2] * T.C22();
+ DataT F30 = Jx[3] * T.C33();
+ DataT F31 = Jy[3] * T.C33();
+ DataT F40 = Jx[4] * T.C44();
+ DataT F41 = Jy[4] * T.C44();
// Jx[5,6] and Jy[5,6] are 0.
- fvec S00 = m.Dx2() + F00 + Jx[2] * F20 + Jx[3] * F30 + Jx[4] * F40;
- fvec S10 = m.Dxy() + F10 + Jy[2] * F20 + Jy[3] * F30 + Jy[4] * F40;
- fvec S11 = m.Dy2() + F11 + Jy[2] * F21 + Jy[3] * F31 + Jy[4] * F41;
+ DataT S00 = m.Dx2() + F00 + Jx[2] * F20 + Jx[3] * F30 + Jx[4] * F40;
+ DataT S10 = m.Dxy() + F10 + Jy[2] * F20 + Jy[3] * F30 + Jy[4] * F40;
+ DataT S11 = m.Dy2() + F11 + Jy[2] * F21 + Jy[3] * F31 + Jy[4] * F41;
- fvec si = fvec(1.) / (S00 * S11 - S10 * S10);
+ DataT si = DataT(1.) / (S00 * S11 - S10 * S10);
- fvec S00tmp = S00;
- S00 = si * S11;
- S10 = -si * S10;
- S11 = si * S00tmp;
+ DataT S00tmp = S00;
+ S00 = si * S11;
+ S10 = -si * S10;
+ S11 = si * S00tmp;
- T.ChiSq() += zeta0 * zeta0 * S00 + fvec(2.) * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
+ T.ChiSq() += zeta0 * zeta0 * S00 + DataT(2.) * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
T.Ndf() += m.NdfX() + m.NdfY();
- fvec K00 = F00 * S00 + F01 * S10;
- fvec K01 = F00 * S10 + F01 * S11;
- fvec K10 = F10 * S00 + F11 * S10;
- fvec K11 = F10 * S10 + F11 * S11;
- fvec K20 = F20 * S00 + F21 * S10;
- fvec K21 = F20 * S10 + F21 * S11;
- fvec K30 = F30 * S00 + F31 * S10;
- fvec K31 = F30 * S10 + F31 * S11;
- fvec K40 = F40 * S00 + F41 * S10;
- fvec K41 = F40 * S10 + F41 * S11;
+ DataT K00 = F00 * S00 + F01 * S10;
+ DataT K01 = F00 * S10 + F01 * S11;
+ DataT K10 = F10 * S00 + F11 * S10;
+ DataT K11 = F10 * S10 + F11 * S11;
+ DataT K20 = F20 * S00 + F21 * S10;
+ DataT K21 = F20 * S10 + F21 * S11;
+ DataT K30 = F30 * S00 + F31 * S10;
+ DataT K31 = F30 * S10 + F31 * S11;
+ DataT K40 = F40 * S00 + F41 * S10;
+ DataT K41 = F40 * S10 + F41 * S11;
T.X() -= K00 * zeta0 + K01 * zeta1;
T.Y() -= K10 * zeta0 + K11 * zeta1;
@@ -418,23 +424,24 @@ namespace cbm::algo::ca
}
- void TrackFit::MeasureVelocityWithQp()
+ template<typename DataT>
+ void TrackFit<DataT>::MeasureVelocityWithQp()
{
// measure velocity using measured qp
// assuming particle mass == fMass;
- cnst kClightNsInv = fvec(constants::phys::SpeedOfLightInv);
+ const DataT kClightNsInv = DataT(constants::phys::SpeedOfLightInv);
- fvec zeta, HCH;
- fvec F0, F1, F2, F3, F4, F5, F6;
- fvec K1, K2, K3, K4, K5, K6;
+ DataT zeta, HCH;
+ DataT F0, F1, F2, F3, F4, F5, F6;
+ DataT K1, K2, K3, K4, K5, K6;
- //FilterVi(sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv);
+ //FilterVi(sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv);
//return;
- fvec vi0 = sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
+ DataT vi0 = sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
- fvec h = fMass2 * fQp0 / sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
+ DataT h = fMass2 * fQp0 / sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
zeta = vi0 + h * (fTr.Qp() - fQp0) - fTr.Vi();
@@ -454,8 +461,8 @@ namespace cbm::algo::ca
HCH = F4 * h - F6;
- fvec wi = fMaskF / HCH;
- fvec zetawi = fMaskF * zeta / HCH;
+ DataT wi = fMaskF / HCH;
+ DataT zetawi = fMaskF * zeta / HCH;
fTr.ChiSqTime() += zeta * zeta * wi;
fTr.NdfTime() += fMaskF;
@@ -509,16 +516,17 @@ namespace cbm::algo::ca
fTr.C65() -= K6 * F5;
fTr.C66() -= K6 * F6;
- // fTr.Vi()( fTr.Vi() < fvec(TrackParamV::kClightNsInv) ) = fvec(TrackParamV::kClightNsInv);
+ // fTr.Vi()( fTr.Vi() < DataT(TrackParamV::kClightNsInv) ) = DataT(TrackParamV::kClightNsInv);
}
- void TrackFit::FilterVi(fvec vi)
+ template<typename DataT>
+ void TrackFit<DataT>::FilterVi(DataT vi)
{
// set inverse velocity to vi
- fvec zeta, HCH;
- fvec F0, F1, F2, F3, F4, F5, F6;
- fvec K1, K2, K3, K4, K5, K6;
+ DataT zeta, HCH;
+ DataT F0, F1, F2, F3, F4, F5, F6;
+ DataT K1, K2, K3, K4, K5; //, K6;
zeta = fTr.Vi() - vi;
@@ -536,8 +544,8 @@ namespace cbm::algo::ca
HCH = F6;
- fvec wi = fMaskF / HCH;
- fvec zetawi = fMaskF * zeta / HCH;
+ DataT wi = fMaskF / HCH;
+ DataT zetawi = fMaskF * zeta / HCH;
fTr.ChiSqTime() += zeta * zeta * wi;
fTr.NdfTime() += fMaskF;
@@ -546,7 +554,7 @@ namespace cbm::algo::ca
K3 = F3 * wi;
K4 = F4 * wi;
K5 = F5 * wi;
- K6 = F6 * wi;
+ // K6 = F6 * wi;
fTr.X() -= F0 * zetawi;
fTr.Y() -= F1 * zetawi;
@@ -591,18 +599,20 @@ namespace cbm::algo::ca
//fTr.C64() -= K6 * F4;
//fTr.C65() -= K6 * F5;
//fTr.C66() -= K6 * F6;
- fTr.C60() = fvec(0.);
- fTr.C61() = fvec(0.);
- fTr.C62() = fvec(0.);
- fTr.C63() = fvec(0.);
- fTr.C64() = fvec(0.);
- fTr.C65() = fvec(0.);
- fTr.C66() = fvec(1.e-8); // just for a case..
+ fTr.C60() = DataT(0.);
+ fTr.C61() = DataT(0.);
+ fTr.C62() = DataT(0.);
+ fTr.C63() = DataT(0.);
+ fTr.C64() = DataT(0.);
+ fTr.C65() = DataT(0.);
+ fTr.C66() = DataT(1.e-8); // just for a case..
}
- void TrackFit::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (fvec z_out, // extrapolate to this z position
- const ca::FieldRegion<fvec>& F)
+ template<typename DataT>
+ void
+ TrackFit<DataT>::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
+ (DataT z_out, // extrapolate to this z position
+ const ca::FieldRegion<DataT>& F)
{
// use Q/p linearisation at fQp0
@@ -610,18 +620,20 @@ namespace cbm::algo::ca
ExtrapolateLineNoField(z_out);
}
else {
- fvec sgn = iif(fTr.GetZ() < z_out, fvec(1.), fvec(-1.));
- while (!(fMaskF * abs(z_out - fTr.GetZ()) <= fvec(1.e-6)).isFull()) {
- fvec zNew = fTr.GetZ() + sgn * fMaxExtraplationStep; // max. 50 cm step
- zNew(sgn * (z_out - zNew) <= fvec(0.)) = z_out;
+ DataT sgn = utils::iif(fTr.GetZ() < z_out, DataT(1.), DataT(-1.));
+ while (!utils::isFull(fMaskF * utils::fabs(z_out - fTr.GetZ()) <= DataT(1.e-6))) {
+ DataT zNew = fTr.GetZ() + sgn * fMaxExtraplationStep; // max. 50 cm step
+ zNew = utils::iif(sgn * (z_out - zNew) <= DataT(0.), z_out, zNew);
ExtrapolateStep(zNew, F);
}
}
}
- void TrackFit::ExtrapolateStep // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (fvec zOut, // extrapolate to this z position
- const ca::FieldRegion<fvec>& Field)
+ template<typename DataT>
+ void TrackFit<
+ DataT>::ExtrapolateStep // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
+ (DataT zOut, // extrapolate to this z position
+ const ca::FieldRegion<DataT>& Field)
{
// use Q/p linearisation at fQp0
// implementation of the Runge-Kutta method without optimization
@@ -662,48 +674,48 @@ namespace cbm::algo::ca
//----------------------------------------------------------------
- cnst zMasked = iif(fMask, zOut, fTr.GetZ());
+ cnst zMasked = utils::iif(fMask, zOut, fTr.GetZ());
cnst h = (zMasked - fTr.GetZ());
- cnst stepDz[5] = {0., 0., h * fvec(0.5), h * fvec(0.5), h};
+ cnst stepDz[5] = {0., 0., h * DataT(0.5), h * DataT(0.5), h};
- fvec f[5][7] = {{0.}}; // ( d*/dz ) [step]
- fvec F[5][7][7] = {{{0.}}}; // ( d *new [step] / d *old )
+ DataT f[5][7] = {{DataT(0.)}}; // ( d*/dz ) [step]
+ DataT F[5][7][7] = {{{DataT(0.)}}}; // ( d *new [step] / d *old )
// Runge-Kutta steps
//
- fvec r0[7] = {fTr.X(), fTr.Y(), fTr.Tx(), fTr.Ty(), fQp0, fTr.Time(), fTr.Vi()};
- fvec R0[7][7] = {{0.}};
+ DataT r0[7] = {fTr.X(), fTr.Y(), fTr.Tx(), fTr.Ty(), fQp0, fTr.Time(), fTr.Vi()};
+ DataT R0[7][7] = {{DataT(0.)}};
for (int i = 0; i < 7; ++i) {
R0[i][i] = 1.;
}
for (int step = 1; step <= 4; ++step) {
- fvec rstep[7] = {{0.}};
+ DataT rstep[7] = {DataT(0.)};
for (int i = 0; i < 7; ++i) {
rstep[i] = r0[i] + stepDz[step] * f[step - 1][i];
}
- fvec z = fTr.GetZ() + stepDz[step];
+ DataT z = fTr.GetZ() + stepDz[step];
ca::FieldValue B = Field.Get(rstep[0], rstep[1], z);
- fvec Bx = B.x;
- fvec By = B.y;
- fvec Bz = B.z;
+ DataT Bx = B.x;
+ DataT By = B.y;
+ DataT Bz = B.z;
// NOTE: SZh 13.08.2024: The rstep[0] and rstep[1] make no effect on the B value, if Field is an approximation
- fvec tx = rstep[2];
- fvec ty = rstep[3];
- fvec tx2 = tx * tx;
- fvec ty2 = ty * ty;
- fvec txty = tx * ty;
- fvec L2 = fvec(1.) + tx2 + ty2;
- fvec L2i = fvec(1.) / L2;
- fvec L = sqrt(L2);
- fvec cL = c_light * L;
- fvec cLqp0 = cL * fQp0;
+ DataT tx = rstep[2];
+ DataT ty = rstep[3];
+ DataT tx2 = tx * tx;
+ DataT ty2 = ty * ty;
+ DataT txty = tx * ty;
+ DataT L2 = DataT(1.) + tx2 + ty2;
+ DataT L2i = DataT(1.) / L2;
+ DataT L = sqrt(L2);
+ DataT cL = c_light * L;
+ DataT cLqp0 = cL * fQp0;
f[step][0] = tx;
F[step][0][2] = 1.;
@@ -711,23 +723,23 @@ namespace cbm::algo::ca
f[step][1] = ty;
F[step][1][3] = 1.;
- fvec f2tmp = txty * Bx - (fvec(1.) + tx2) * By + ty * Bz;
- f[step][2] = cLqp0 * f2tmp;
+ DataT f2tmp = txty * Bx - (DataT(1.) + tx2) * By + ty * Bz;
+ f[step][2] = cLqp0 * f2tmp;
- F[step][2][2] = cLqp0 * (tx * f2tmp * L2i + ty * Bx - fvec(2.) * tx * By);
+ F[step][2][2] = cLqp0 * (tx * f2tmp * L2i + ty * Bx - DataT(2.) * tx * By);
F[step][2][3] = cLqp0 * (ty * f2tmp * L2i + tx * Bx + Bz);
F[step][2][4] = cL * f2tmp;
- fvec f3tmp = -txty * By - tx * Bz + (fvec(1.) + ty2) * Bx;
+ DataT f3tmp = -txty * By - tx * Bz + (DataT(1.) + ty2) * Bx;
f[step][3] = cLqp0 * f3tmp;
F[step][3][2] = cLqp0 * (tx * f3tmp * L2i - ty * By - Bz);
- F[step][3][3] = cLqp0 * (ty * f3tmp * L2i + fvec(2.) * ty * Bx - tx * By);
+ F[step][3][3] = cLqp0 * (ty * f3tmp * L2i + DataT(2.) * ty * Bx - tx * By);
F[step][3][4] = cL * f3tmp;
f[step][4] = 0.;
if (fDoFitVelocity) {
- fvec vi = rstep[6];
+ DataT vi = rstep[6];
f[step][5] = vi * L;
F[step][5][2] = vi * tx / L;
F[step][5][3] = vi * ty / L;
@@ -736,12 +748,12 @@ namespace cbm::algo::ca
F[step][5][6] = L;
}
else {
- fvec vi = sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * fvec(constants::phys::SpeedOfLightInv);
+ DataT vi = sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * DataT(constants::phys::SpeedOfLightInv);
f[step][5] = vi * L;
F[step][5][2] = vi * tx / L;
F[step][5][3] = vi * ty / L;
F[step][5][4] =
- fMass2 * fQp0 * L / sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * fvec(constants::phys::SpeedOfLightInv);
+ fMass2 * fQp0 * L / sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * DataT(constants::phys::SpeedOfLightInv);
F[step][5][5] = 0.;
F[step][5][6] = 0.;
}
@@ -750,12 +762,12 @@ namespace cbm::algo::ca
} // end of Runge-Kutta step
- fvec r[7] = {{0.}}; // extrapolated parameters
- fvec R[7][7] = {{0.}}; // Jacobian of the extrapolation
+ DataT r[7] = {DataT(0.)}; // extrapolated parameters
+ DataT R[7][7] = {{DataT(0.)}}; // Jacobian of the extrapolation
- cnst stepW[5] = {0., h / fvec(6.), h / fvec(3.), h / fvec(3.), h / fvec(6.)};
+ cnst stepW[5] = {0., h / DataT(6.), h / DataT(3.), h / DataT(3.), h / DataT(6.)};
- fvec k[5][7][7] = {{0.}};
+ DataT k[5][7][7] = {{{DataT(0.)}}};
for (int step = 1; step <= 4; ++step) {
for (int i = 0; i < 7; i++) {
for (int j = 0; j < 7; j++) {
@@ -776,7 +788,7 @@ namespace cbm::algo::ca
}
}
- fvec dqp = fTr.Qp() - fQp0;
+ DataT dqp = fTr.Qp() - fQp0;
for (int i = 0; i < 7; i++) {
r[i] = r0[i];
@@ -797,19 +809,19 @@ namespace cbm::algo::ca
fTr.Time() = r[5];
fTr.Vi() = r[6];
- //fTr.Vi()( fTr.Vi() < fvec(TrackParamV::kClightNsInv) ) = fvec(TrackParamV::kClightNsInv);
+ //fTr.Vi()( fTr.Vi() < DataT(TrackParamV::kClightNsInv) ) = DataT(TrackParamV::kClightNsInv);
fTr.Z() = zMasked;
// covariance matrix transport
- fvec C[7][7];
+ DataT C[7][7];
for (int i = 0; i < 7; i++) {
for (int j = 0; j < 7; j++) {
C[i][j] = fTr.C(i, j);
}
}
- fvec RC[7][7];
+ DataT RC[7][7];
for (int i = 0; i < 7; i++) {
for (int j = 0; j < 7; j++) {
RC[i][j] = 0.;
@@ -820,7 +832,7 @@ namespace cbm::algo::ca
}
for (int i = 0; i < 7; i++) {
for (int j = 0; j < 7; j++) {
- fvec Cij = 0.;
+ DataT Cij = 0.;
for (int m = 0; m < 7; m++) {
Cij += RC[i][m] * R[j][m];
}
@@ -830,18 +842,20 @@ namespace cbm::algo::ca
}
- void TrackFit::ExtrapolateLine(fvec z_out, const ca::FieldRegion<fvec>& F)
+ template<typename DataT>
+ void TrackFit<DataT>::ExtrapolateLine(DataT z_out, const ca::FieldRegion<DataT>& F)
{
// extrapolate the track assuming fQp0 == 0
// TODO: write special simplified procedure
//
auto qp0 = fQp0;
- fQp0 = fvec(0.);
+ fQp0 = DataT(0.);
Extrapolate(z_out, F);
fQp0 = qp0;
}
- void TrackFit::ExtrapolateLineNoField(fvec zOut)
+ template<typename DataT>
+ void TrackFit<DataT>::ExtrapolateLineNoField(DataT zOut)
{
// extrapolate the track assuming no field
@@ -850,26 +864,26 @@ namespace cbm::algo::ca
// t += dz * sqrt ( 1 + tx*tx + ty*ty ) * vi
if (0) { // debug: full Runge-Kutta extrapolation with zero field values
- FieldRegion<fvec> F;
+ FieldRegion<DataT> F;
ExtrapolateStep(zOut, F);
return;
}
auto& t = fTr; // use reference to shorten the text
- cnst zMasked = iif(fMask, zOut, t.GetZ());
+ cnst zMasked = utils::iif(fMask, zOut, t.GetZ());
cnst dz = (zMasked - t.GetZ());
- fvec tx = t.GetTx();
- fvec ty = t.GetTy();
- fvec vi = t.GetVi();
+ DataT tx = t.GetTx();
+ DataT ty = t.GetTy();
+ DataT vi = t.GetVi();
- fvec L = sqrt(fvec(1.) + tx * tx + ty * ty);
+ DataT L = sqrt(DataT(1.) + tx * tx + ty * ty);
- fvec j52 = dz * tx * vi / L;
- fvec j53 = dz * ty * vi / L;
- fvec j56 = dz * L;
+ DataT j52 = dz * tx * vi / L;
+ DataT j53 = dz * ty * vi / L;
+ DataT j56 = dz * L;
// transport parameters
@@ -891,34 +905,34 @@ namespace cbm::algo::ca
// JC = J * C
- fvec jc00 = t.C00() + dz * t.C20();
- //fvec jc01 = t.C01() + dz * t.C21();
- fvec jc02 = t.C02() + dz * t.C22();
- //fvec jc03 = t.C03() + dz * t.C23();
- //fvec jc04 = t.C04() + dz * t.C24();
- //fvec jc05 = t.C05() + dz * t.C25();
- //fvec jc06 = t.C06() + dz * t.C26();
+ DataT jc00 = t.C00() + dz * t.C20();
+ //DataT jc01 = t.C01() + dz * t.C21();
+ DataT jc02 = t.C02() + dz * t.C22();
+ //DataT jc03 = t.C03() + dz * t.C23();
+ //DataT jc04 = t.C04() + dz * t.C24();
+ //DataT jc05 = t.C05() + dz * t.C25();
+ //DataT jc06 = t.C06() + dz * t.C26();
- fvec jc10 = t.C10() + dz * t.C30();
- fvec jc11 = t.C11() + dz * t.C31();
- fvec jc12 = t.C12() + dz * t.C32();
- fvec jc13 = t.C13() + dz * t.C33();
- //fvec jc14 = t.C14() + dz * t.C34();
- //fvec jc15 = t.C15() + dz * t.C35();
- //fvec jc16 = t.C16() + dz * t.C36();
+ DataT jc10 = t.C10() + dz * t.C30();
+ DataT jc11 = t.C11() + dz * t.C31();
+ DataT jc12 = t.C12() + dz * t.C32();
+ DataT jc13 = t.C13() + dz * t.C33();
+ //DataT jc14 = t.C14() + dz * t.C34();
+ //DataT jc15 = t.C15() + dz * t.C35();
+ //DataT jc16 = t.C16() + dz * t.C36();
// jc2? = t.C2?
// jc3? = t.C3?
// jc4? = t.C4?
- fvec jc50 = t.C50() + j52 * t.C20() + j53 * t.C30() + j56 * t.C60();
- fvec jc51 = t.C51() + j52 * t.C21() + j53 * t.C31() + j56 * t.C61();
- fvec jc52 = t.C52() + j52 * t.C22() + j53 * t.C32() + j56 * t.C62();
- fvec jc53 = t.C53() + j52 * t.C23() + j53 * t.C33() + j56 * t.C63();
- fvec jc54 = t.C54() + j52 * t.C24() + j53 * t.C34() + j56 * t.C64();
- fvec jc55 = t.C55() + j52 * t.C25() + j53 * t.C35() + j56 * t.C65();
- fvec jc56 = t.C56() + j52 * t.C26() + j53 * t.C36() + j56 * t.C66();
+ DataT jc50 = t.C50() + j52 * t.C20() + j53 * t.C30() + j56 * t.C60();
+ DataT jc51 = t.C51() + j52 * t.C21() + j53 * t.C31() + j56 * t.C61();
+ DataT jc52 = t.C52() + j52 * t.C22() + j53 * t.C32() + j56 * t.C62();
+ DataT jc53 = t.C53() + j52 * t.C23() + j53 * t.C33() + j56 * t.C63();
+ DataT jc54 = t.C54() + j52 * t.C24() + j53 * t.C34() + j56 * t.C64();
+ DataT jc55 = t.C55() + j52 * t.C25() + j53 * t.C35() + j56 * t.C65();
+ DataT jc56 = t.C56() + j52 * t.C26() + j53 * t.C36() + j56 * t.C66();
// jc6? = t.C6?
@@ -972,90 +986,93 @@ namespace cbm::algo::ca
}
- void TrackFit::MultipleScattering(fvec radThick, fvec qp0)
+ template<typename DataT>
+ void TrackFit<DataT>::MultipleScattering(DataT radThick, DataT qp0)
{
- cnst ONE = 1.;
-
- fvec tx = fTr.Tx();
- fvec ty = fTr.Ty();
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
+ cnst kONE = 1.;
+
+ DataT tx = fTr.Tx();
+ DataT ty = fTr.Ty();
+ DataT txtx = tx * tx;
+ DataT tyty = ty * ty;
+ DataT txtx1 = txtx + kONE;
+ DataT h = txtx + tyty;
+ DataT t = sqrt(txtx1 + tyty);
+ DataT h2 = h * h;
+ DataT qp0t = qp0 * t;
cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
+ DataT s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //DataT a = ( (kONE+mass2*qp0*qp0t)*radThick*s0*s0 );
+ DataT a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
- fTr.C22()(fMask) += txtx1 * a;
- fTr.C32()(fMask) += tx * ty * a;
- fTr.C33()(fMask) += (ONE + tyty) * a;
+ fTr.C22() += fMaskF * txtx1 * a;
+ fTr.C32() += fMaskF * tx * ty * a;
+ fTr.C33() += fMaskF * (kONE + tyty);
}
- void TrackFit::MultipleScatteringInThickMaterial(fvec radThick, fvec thickness, bool fDownstream)
+ template<typename DataT>
+ void TrackFit<DataT>::MultipleScatteringInThickMaterial(DataT radThick, DataT thickness, bool fDownstream)
{
- cnst ONE = 1.;
-
- fvec tx = fTr.Tx();
- fvec ty = fTr.Ty();
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = fQp0 * t;
-
- cnst c1(0.0136), c2 = c1 * fvec(0.038), c3 = c2 * fvec(0.5), c4 = -c3 / fvec(2.0), c5 = c3 / fvec(3.0),
- c6 = -c3 / fvec(4.0);
-
- fvec s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + fMass2 * fQp0 * qp0t) * radThick * s0 * s0);
-
- fvec D = (fDownstream) ? fvec(1.) : fvec(-1.);
- fvec T23 = (thickness * thickness) / fvec(3.0);
- fvec T2 = thickness / fvec(2.0);
-
- fTr.C00()(fMask) += txtx1 * a * T23;
- fTr.C10()(fMask) += tx * ty * a * T23;
- fTr.C20()(fMask) += txtx1 * a * D * T2;
- fTr.C30()(fMask) += tx * ty * a * D * T2;
-
- fTr.C11()(fMask) += (ONE + tyty) * a * T23;
- fTr.C21()(fMask) += tx * ty * a * D * T2;
- fTr.C31()(fMask) += (ONE + tyty) * a * D * T2;
-
- fTr.C22()(fMask) += txtx1 * a;
- fTr.C32()(fMask) += tx * ty * a;
- fTr.C33()(fMask) += (ONE + tyty) * a;
+ cnst kONE = 1.;
+
+ DataT tx = fTr.Tx();
+ DataT ty = fTr.Ty();
+ DataT txtx = tx * tx;
+ DataT tyty = ty * ty;
+ DataT txtx1 = txtx + kONE;
+ DataT h = txtx + tyty;
+ DataT t = sqrt(txtx1 + tyty);
+ DataT h2 = h * h;
+ DataT qp0t = fQp0 * t;
+
+ cnst c1(0.0136), c2 = c1 * DataT(0.038), c3 = c2 * DataT(0.5), c4 = -c3 / DataT(2.0), c5 = c3 / DataT(3.0),
+ c6 = -c3 / DataT(4.0);
+
+ DataT s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //DataT a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
+ DataT a = ((t + fMass2 * fQp0 * qp0t) * radThick * s0 * s0);
+
+ DataT D = (fDownstream) ? DataT(1.) : DataT(-1.);
+ DataT T23 = (thickness * thickness) / DataT(3.0);
+ DataT T2 = thickness / DataT(2.0);
+
+ fTr.C00() += fMaskF * txtx1 * a * T23;
+ fTr.C10() += fMaskF * tx * ty * a * T23;
+ fTr.C20() += fMaskF * txtx1 * a * D * T2;
+ fTr.C30() += fMaskF * tx * ty * a * D * T2;
+
+ fTr.C11() += fMaskF * (kONE + tyty) * a * T23;
+ fTr.C21() += fMaskF * tx * ty * a * D * T2;
+ fTr.C31() += fMaskF * (kONE + tyty) * a * D * T2;
+
+ fTr.C22() += fMaskF * txtx1 * a;
+ fTr.C32() += fMaskF * tx * ty * a;
+ fTr.C33() += fMaskF * (kONE + tyty) * a;
}
- void TrackFit::EnergyLossCorrection(fvec radThick, fvec upstreamDirection)
+ template<typename DataT>
+ void TrackFit<DataT>::EnergyLossCorrection(DataT radThick, DataT upstreamDirection)
{
cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
- cnst qp02 = max(fQp0 * fQp0, qp2cut);
- cnst p2 = fvec(1.) / qp02;
+ cnst qp02 = utils::max(fQp0 * fQp0, qp2cut);
+ cnst p2 = DataT(1.) / qp02;
cnst E2 = fMass2 + p2;
cnst bethe = ApproximateBetheBloch(p2 / fMass2);
- fvec tr = sqrt(fvec(1.f) + fTr.Tx() * fTr.Tx() + fTr.Ty() * fTr.Ty());
+ DataT tr = sqrt(DataT(1.f) + fTr.Tx() * fTr.Tx() + fTr.Ty() * fTr.Ty());
cnst dE = bethe * radThick * tr * 2.33f * 9.34961f;
cnst ECorrected = sqrt(E2) + upstreamDirection * dE;
cnst E2Corrected = ECorrected * ECorrected;
- fvec corr = sqrt(p2 / (E2Corrected - fMass2));
- fmask ok = (corr == corr) && fMask;
- corr = iif(ok, corr, fvec::One());
+ DataT corr = sqrt(p2 / (E2Corrected - fMass2));
+ DataTmask ok = (corr == corr) && fMask;
+ corr = utils::iif(ok, corr, DataT(1.));
fQp0 *= corr;
fTr.Qp() *= corr;
@@ -1068,15 +1085,16 @@ namespace cbm::algo::ca
}
- void TrackFit::EnergyLossCorrection(int atomicZ, fscal atomicA, fscal rho, fscal radLen, fvec radThick,
- fvec upstreamDirection)
+ template<typename DataT>
+ void TrackFit<DataT>::EnergyLossCorrection(int atomicZ, DataTscal atomicA, DataTscal rho, DataTscal radLen,
+ DataT radThick, DataT upstreamDirection)
{
cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
- cnst qp02 = max(fQp0 * fQp0, qp2cut);
- cnst p2 = fvec(1.) / qp02;
+ cnst qp02 = utils::max(fQp0 * fQp0, qp2cut);
+ cnst p2 = DataT(1.) / qp02;
cnst E2 = fMass2 + p2;
- fvec i;
+ DataT i;
if (atomicZ < 13)
i = (12. * atomicZ + 7.) * 1.e-9;
else
@@ -1084,121 +1102,124 @@ namespace cbm::algo::ca
cnst bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
- fvec tr = sqrt(fvec(1.f) + fTr.Tx() * fTr.Tx() + fTr.Ty() * fTr.Ty());
+ DataT tr = sqrt(DataT(1.f) + fTr.Tx() * fTr.Tx() + fTr.Ty() * fTr.Ty());
- fvec dE = bethe * radThick * tr * radLen * rho;
+ DataT dE = bethe * radThick * tr * radLen * rho;
cnst E2Corrected = (sqrt(E2) + upstreamDirection * dE) * (sqrt(E2) + upstreamDirection * dE);
- fvec corr = sqrt(p2 / (E2Corrected - fMass2));
- fmask ok = (corr == corr) && fMask;
- corr = iif(ok, corr, fvec::One());
+ DataT corr = sqrt(p2 / (E2Corrected - fMass2));
+ DataTmask ok = (corr == corr) && fMask;
+ corr = utils::iif(ok, corr, DataT(1.));
fQp0 *= corr;
fTr.Qp() *= corr;
- fvec P(sqrt(p2)); // GeV
+ DataT P(sqrt(p2)); // GeV
- fvec Z(atomicZ);
- fvec A(atomicA);
- fvec RHO(rho);
+ DataT Z(atomicZ);
+ DataT A(atomicA);
+ DataT RHO(rho);
- fvec STEP = radThick * tr * radLen;
- fvec EMASS(0.511 * 1e-3); // GeV
+ DataT STEP = radThick * tr * radLen;
+ DataT EMASS(0.511 * 1e-3); // GeV
- fvec BETA = P / sqrt(E2Corrected);
- fvec GAMMA = sqrt(E2Corrected) / fMass;
+ DataT BETA = P / sqrt(E2Corrected);
+ DataT GAMMA = sqrt(E2Corrected) / fMass;
// Calculate xi factor (KeV).
- fvec XI = (fvec(153.5) * Z * STEP * RHO) / (A * BETA * BETA);
+ DataT XI = (DataT(153.5) * Z * STEP * RHO) / (A * BETA * BETA);
// Maximum energy transfer to atomic electron (KeV).
- fvec ETA = BETA * GAMMA;
- fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / fMass;
- fvec F1 = fvec(2.) * EMASS * ETASQ;
- fvec F2 = fvec(1.) + fvec(2.) * RATIO * GAMMA + RATIO * RATIO;
- fvec EMAX = fvec(1e6) * F1 / F2;
+ DataT ETA = BETA * GAMMA;
+ DataT ETASQ = ETA * ETA;
+ DataT RATIO = EMASS / fMass;
+ DataT F1 = DataT(2.) * EMASS * ETASQ;
+ DataT F2 = DataT(1.) + DataT(2.) * RATIO * GAMMA + RATIO * RATIO;
+ DataT EMAX = DataT(1e6) * F1 / F2;
- fvec DEDX2 = XI * EMAX * (fvec(1.) - (BETA * BETA / fvec(2.))) * fvec(1e-12);
+ DataT DEDX2 = XI * EMAX * (DataT(1.) - (BETA * BETA / DataT(2.))) * DataT(1e-12);
- fvec P2 = P * P;
- fvec SDEDX = (E2 * DEDX2) / (P2 * P2 * P2);
+ DataT P2 = P * P;
+ DataT SDEDX = (E2 * DEDX2) / (P2 * P2 * P2);
// T.fTr.C40() *= corr;
// T.fTr.C41() *= corr;
// T.fTr.C42() *= corr;
// T.fTr.C43() *= corr;
// T.fTr.C44() *= corr*corr;
- fTr.C44() += abs(SDEDX);
+ fTr.C44() += utils::fabs(SDEDX);
}
- void TrackFit::GetExtrapolatedXYline(fvec z, const ca::FieldRegion<fvec>& F, fvec& extrX, fvec& extrY,
- std::array<fvec, ca::TrackParamV::kNtrackParam>& Jx,
- std::array<fvec, ca::TrackParamV::kNtrackParam>& Jy) const
+ template<typename DataT>
+ void TrackFit<DataT>::GetExtrapolatedXYline(DataT z, const ca::FieldRegion<DataT>& F, DataT& extrX, DataT& extrY,
+ std::array<DataT, ca::TrackParamBase<DataT>::kNtrackParam>& Jx,
+ std::array<DataT, ca::TrackParamBase<DataT>::kNtrackParam>& Jy) const
{
// extrapolate track assuming it is straight (qp==0)
// return the extrapolated X, Y and the derivatives of the extrapolated X and Y
cnst c_light(0.000299792458), c1(1.), c2i(0.5), c6i(1. / 6.), c12i(1. / 12.);
- cnst tx = fTr.GetTx();
- cnst ty = fTr.GetTy();
- fvec dz = z - fTr.GetZ();
+ cnst tx = fTr.GetTx();
+ cnst ty = fTr.GetTy();
+ DataT dz = z - fTr.GetZ();
- fvec dz2 = dz * dz;
- fvec dzc6i = dz * c6i;
- fvec dz2c12i = dz2 * c12i;
+ DataT dz2 = dz * dz;
+ DataT dzc6i = dz * c6i;
+ DataT dz2c12i = dz2 * c12i;
- fvec xx = tx * tx;
- fvec yy = ty * ty;
- fvec xy = tx * ty;
+ DataT xx = tx * tx;
+ DataT yy = ty * ty;
+ DataT xy = tx * ty;
- fvec Ay = -xx - c1;
- fvec Bx = yy + c1;
+ DataT Ay = -xx - c1;
+ DataT Bx = yy + c1;
- fvec ctdz2 = c_light * sqrt(c1 + xx + yy) * dz2;
+ DataT ctdz2 = c_light * sqrt(c1 + xx + yy) * dz2;
- fvec Sx = F.cx0 * c2i + F.cx1 * dzc6i + F.cx2 * dz2c12i;
- fvec Sy = F.cy0 * c2i + F.cy1 * dzc6i + F.cy2 * dz2c12i;
- fvec Sz = F.cz0 * c2i + F.cz1 * dzc6i + F.cz2 * dz2c12i;
+ DataT Sx = F.cx0 * c2i + F.cx1 * dzc6i + F.cx2 * dz2c12i;
+ DataT Sy = F.cy0 * c2i + F.cy1 * dzc6i + F.cy2 * dz2c12i;
+ DataT Sz = F.cz0 * c2i + F.cz1 * dzc6i + F.cz2 * dz2c12i;
extrX = fTr.GetX() + tx * dz;
extrY = fTr.GetY() + ty * dz;
- Jx.fill(fvec(0.)); // for a case
- Jy.fill(fvec(0.));
+ Jx.fill(DataT(0.)); // for a case
+ Jy.fill(DataT(0.));
- Jx[0] = fvec(1.);
- Jx[1] = fvec(0.);
+ Jx[0] = DataT(1.);
+ Jx[1] = DataT(0.);
Jx[2] = dz;
- Jx[3] = fvec(0.);
+ Jx[3] = DataT(0.);
Jx[4] = ctdz2 * (Sx * xy + Sy * Ay + Sz * ty);
- Jx[5] = fvec(0.);
- Jx[6] = fvec(0.);
+ Jx[5] = DataT(0.);
+ Jx[6] = DataT(0.);
- Jy[0] = fvec(0.);
- Jy[1] = fvec(1.);
- Jy[2] = fvec(0.);
+ Jy[0] = DataT(0.);
+ Jy[1] = DataT(1.);
+ Jy[2] = DataT(0.);
Jy[3] = dz;
Jy[4] = ctdz2 * (Sx * Bx - Sy * xy - Sz * tx);
- Jy[5] = fvec(0.);
- Jx[6] = fvec(0.);
+ Jy[5] = DataT(0.);
+ Jx[6] = DataT(0.);
}
- void TrackFit::FilterWithTargetAtLine(fvec targZ, const ca::MeasurementXy<fvec>& targXY,
- const ca::FieldRegion<fvec>& F)
+ template<typename DataT>
+ void TrackFit<DataT>::FilterWithTargetAtLine(DataT targZ, const ca::MeasurementXy<DataT>& targXY,
+ const ca::FieldRegion<DataT>& F)
{
// Add the target constraint to a straight line track
- fvec eX, eY;
- std::array<fvec, ca::TrackParamV::kNtrackParam> Jx, Jy;
+ DataT eX, eY;
+ std::array<DataT, ca::TrackParamV::kNtrackParam> Jx, Jy;
GetExtrapolatedXYline(targZ, F, eX, eY, Jx, Jy);
FilterExtrapolatedXY(targXY, eX, eY, Jx, Jy);
}
- fvec TrackFit::ApproximateBetheBloch(fvec bg2)
+ template<typename DataT>
+ DataT TrackFit<DataT>::ApproximateBetheBloch(DataT bg2)
{
//
// This is the parameterization of the Bethe-Bloch formula inspired by Geant.
@@ -1220,31 +1241,32 @@ namespace cbm::algo::ca
cnst kp3 = 173e-9f;
cnst kp4 = 0.49848f;
- constexpr fscal mK = 0.307075e-3f; // [GeV*cm^2/g]
- constexpr fscal _2me = 1.022e-3f; // [GeV/c^2]
- cnst rho = kp0;
- cnst x0 = kp1 * 2.303f;
- cnst x1 = kp2 * 2.303f;
- cnst mI = kp3;
- cnst mZA = kp4;
- cnst maxT = _2me * bg2; // neglecting the electron mass
+ constexpr DataTscal mK = 0.307075e-3f; // [GeV*cm^2/g]
+ constexpr DataTscal _2me = 1.022e-3f; // [GeV/c^2]
+ cnst rho = kp0;
+ cnst x0 = kp1 * 2.303f;
+ cnst x1 = kp2 * 2.303f;
+ cnst mI = kp3;
+ cnst mZA = kp4;
+ cnst maxT = _2me * bg2; // neglecting the electron mass
//*** Density effect
- fvec d2(0.f);
+ DataT d2(0.f);
cnst x = 0.5f * log(bg2);
cnst lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
- fmask init = x > x1;
- d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
- cnst r = (x1 - x) / (x1 - x0);
- init = (x > x0) & (x1 > x);
- d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
+ DataTmask init = x > x1;
+ d2 = utils::iif(init, lhwI + x - 0.5f, DataT(0.));
+ cnst r = (x1 - x) / (x1 - x0);
+ init = (x > x0) & (x1 > x);
+ d2 = utils::iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
- return mK * mZA * (fvec(1.f) + bg2) / bg2
- * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
+ return mK * mZA * (DataT(1.f) + bg2) / bg2
+ * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (DataT(1.f) + bg2) - d2);
}
- fvec TrackFit::ApproximateBetheBloch(fvec bg2, fvec kp0, fvec kp1, fvec kp2, fvec kp3, fvec kp4)
+ template<typename DataT>
+ DataT TrackFit<DataT>::ApproximateBetheBloch(DataT bg2, DataT kp0, DataT kp1, DataT kp2, DataT kp3, DataT kp4)
{
//
// This is the parameterization of the Bethe-Bloch formula inspired by Geant.
@@ -1266,51 +1288,52 @@ namespace cbm::algo::ca
// cnst &kp3 = 173e-9f;
// cnst &kp4 = 0.49848f;
- constexpr fscal mK = 0.307075e-3f; // [GeV*cm^2/g]
- constexpr fscal _2me = 1.022e-3f; // [GeV/c^2]
- fvec rho = kp0;
- cnst x0 = kp1 * 2.303f;
- cnst x1 = kp2 * 2.303f;
- fvec mI = kp3;
- fvec mZA = kp4;
- cnst maxT = _2me * bg2; // neglecting the electron mass
+ constexpr DataTscal mK = 0.307075e-3f; // [GeV*cm^2/g]
+ constexpr DataTscal _2me = 1.022e-3f; // [GeV/c^2]
+ DataT rho = kp0;
+ cnst x0 = kp1 * 2.303f;
+ cnst x1 = kp2 * 2.303f;
+ DataT mI = kp3;
+ DataT mZA = kp4;
+ cnst maxT = _2me * bg2; // neglecting the electron mass
//*** Density effect
- fvec d2(0.f);
+ DataT d2(0.f);
cnst x = 0.5f * log(bg2);
cnst lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
- fmask init = x > x1;
- d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
- cnst r = (x1 - x) / (x1 - x0);
- init = (x > x0) & (x1 > x);
- d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
+ DataTmask init = x > x1;
+ d2 = utils::iif(init, lhwI + x - 0.5f, DataT(0.));
+ cnst r = (x1 - x) / (x1 - x0);
+ init = (x > x0) & (x1 > x);
+ d2 = utils::iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
- return mK * mZA * (fvec(1.f) + bg2) / bg2
- * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
+ return mK * mZA * (DataT(1.f) + bg2) / bg2
+ * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (DataT(1.f) + bg2) - d2);
}
- std::tuple<fvec, fvec> TrackFit::GetChi2XChi2U(ca::MeasurementXy<fvec> m, fvec x, fvec y, fvec C00, fvec C10,
- fvec C11)
+ template<typename DataT>
+ std::tuple<DataT, DataT> TrackFit<DataT>::GetChi2XChi2U(ca::MeasurementXy<DataT> m, DataT x, DataT y, DataT C00,
+ DataT C10, DataT C11)
{
- fvec chi2x{0.};
+ DataT chi2x{0.};
{ // filter X measurement
- fvec zeta = x - m.X();
+ DataT zeta = x - m.X();
// F = CH'
- fvec F0 = C00;
- fvec F1 = C10;
+ DataT F0 = C00;
+ DataT F1 = C10;
- fvec HCH = F0;
+ DataT HCH = F0;
- fvec wi = fvec(1.) / (m.Dx2() + HCH);
- fvec zetawi = zeta * wi;
- chi2x = m.NdfX() * zeta * zetawi;
+ DataT wi = DataT(1.) / (m.Dx2() + HCH);
+ DataT zetawi = zeta * wi;
+ chi2x = m.NdfX() * zeta * zetawi;
- fvec K1 = F1 * wi;
+ DataT K1 = F1 * wi;
x -= F0 * zetawi;
y -= F1 * zetawi;
@@ -1320,68 +1343,70 @@ namespace cbm::algo::ca
C11 -= K1 * F1;
}
- fvec chi2u{0.};
+ DataT chi2u{0.};
{ // filter U measurement, we need only chi2 here
- fvec cosPhi = -m.Dxy() / m.Dx2();
- fvec u = cosPhi * m.X() + m.Y();
- fvec du2 = m.Dy2() + cosPhi * m.Dxy();
+ DataT cosPhi = -m.Dxy() / m.Dx2();
+ DataT u = cosPhi * m.X() + m.Y();
+ DataT du2 = m.Dy2() + cosPhi * m.Dxy();
- fvec zeta = cosPhi * x + y - u;
+ DataT zeta = cosPhi * x + y - u;
// F = CH'
- fvec F0 = cosPhi * C00 + C10;
- fvec F1 = cosPhi * C10 + C11;
+ DataT F0 = cosPhi * C00 + C10;
+ DataT F1 = cosPhi * C10 + C11;
- fvec HCH = (F0 * cosPhi + F1);
+ DataT HCH = (F0 * cosPhi + F1);
chi2u += m.NdfY() * zeta * zeta / (du2 + HCH);
}
- return std::tuple<fvec, fvec>(chi2x, chi2u);
+ return std::tuple<DataT, DataT>(chi2x, chi2u);
}
- void TrackFit::GuessTrack(const fvec& trackZ, const fvec hitX[], const fvec hitY[], const fvec hitZ[],
- const fvec hitT[], const fvec By[], const fmask hitW[], const fmask hitWtime[], int NHits)
+ template<typename DataT>
+ void TrackFit<DataT>::GuessTrack(const DataT& trackZ, const DataT hitX[], const DataT hitY[], const DataT hitZ[],
+ const DataT hitT[], const DataT By[], const DataTmask hitW[],
+ const DataTmask hitWtime[], int NHits)
{
// gives nice initial approximation for x,y,tx,ty - almost same as KF fit. qp - is shifted by 4%, resid_ual - ~3.5% (KF fit resid_ual - 1%).
- const fvec c_light(0.000299792458), c_light_i(fvec(1.) / c_light);
+ const DataT c_light(0.000299792458), c_light_i(DataT(1.) / c_light);
- fvec A0 = 0., A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0.;
- fvec a0 = 0., a1 = 0., a2 = 0.;
- fvec b0 = 0., b1 = 0., b2 = 0.;
+ DataT A0 = 0., A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0.;
+ DataT a0 = 0., a1 = 0., a2 = 0.;
+ DataT b0 = 0., b1 = 0., b2 = 0.;
- fvec time = 0.;
- fmask isTimeSet = fmask::Zero();
+ DataT time = 0.;
+ DataTmask isTimeSet = DataTmask(false);
- fvec prevZ = 0.;
- fvec sy = 0., Sy = 0.; // field integrals
+ DataT prevZ = 0.;
+ DataT sy = 0., Sy = 0.; // field integrals
for (int i = 0; i < NHits; i++) {
- fvec w = iif(hitW[i], fvec::One(), fvec::Zero());
+ DataT w = utils::iif(hitW[i], DataT(1.), DataT(0.));
- fmask setTime = (!isTimeSet) && hitWtime[i] && hitW[i];
- time(setTime) = hitT[i];
- isTimeSet = isTimeSet || setTime;
+ DataTmask setTime = (!isTimeSet) && hitWtime[i] && hitW[i];
+ time = utils::iif(setTime, hitT[i], time);
+ isTimeSet = isTimeSet || setTime;
- fvec x = hitX[i];
- fvec y = hitY[i];
- fvec z = hitZ[i] - trackZ;
+ DataT x = hitX[i];
+ DataT y = hitY[i];
+ DataT z = hitZ[i] - trackZ;
{
- fvec dZ = z - prevZ;
- Sy(hitW[i]) += dZ * sy + fvec(0.5) * dZ * dZ * By[i];
- sy(hitW[i]) += dZ * By[i];
- prevZ(hitW[i]) = z;
+ DataT dZ = z - prevZ;
+ Sy += w * dZ * sy + DataT(0.5) * dZ * dZ * By[i];
+ sy += w * dZ * By[i];
+ prevZ = utils::iif(hitW[i], z, prevZ);
}
- fvec S = Sy;
+ DataT S = Sy;
- fvec wz = w * z;
- fvec wS = w * S;
+ DataT wz = w * z;
+ DataT wS = w * S;
A0 += w;
A1 += wz;
@@ -1399,16 +1424,16 @@ namespace cbm::algo::ca
b2 += wS * y;
}
- fvec m00 = A0;
- fvec m01 = A1;
- fvec m02 = A3;
+ DataT m00 = A0;
+ DataT m01 = A1;
+ DataT m02 = A3;
- fvec m11 = A2;
- fvec m12 = A4;
+ DataT m11 = A2;
+ DataT m12 = A4;
- fvec m22 = A5;
+ DataT m22 = A5;
- fvec m21 = m12;
+ DataT m21 = m12;
// { m00 m01 m02 } ( a0 )
// { m01 m11 m12 } = x * ( a1 )
@@ -1429,9 +1454,9 @@ namespace cbm::algo::ca
a2 = m11 * a2 - m21 * a1;
}
- fvec L = 0.;
+ DataT L = 0.;
{ // diagonalization row 2
- L(abs(m22) > fvec(1.e-4)) = a2 / m22;
+ L = utils::iif((utils::fabs(m22) > DataT(1.e-4)), a2 / m22, DataT(0.));
a1 -= L * m12;
a0 -= L * m02;
}
@@ -1445,9 +1470,9 @@ namespace cbm::algo::ca
fTr.X() = a0 / m00;
}
- fvec txtx1 = fvec(1.) + fTr.Tx() * fTr.Tx();
- L = L / txtx1;
- fvec L1 = L * fTr.Tx();
+ DataT txtx1 = DataT(1.) + fTr.Tx() * fTr.Tx();
+ L = L / txtx1;
+ DataT L1 = L * fTr.Tx();
A1 = A1 + A3 * L1;
A2 = A2 + (A4 + A4 + A5 * L1) * L1;
@@ -1469,4 +1494,8 @@ namespace cbm::algo::ca
fQp0 = fTr.Qp();
}
+ template class TrackFit<fvec>;
+ template class TrackFit<float>;
+ template class TrackFit<double>;
+
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaTrackFit.h b/algo/ca/core/tracking/CaTrackFit.h
index 619eabf470fedacc8a1e18069db47c214db2224d..a4f6b1f8cf7b8d8f5c240d02badba8977828a7fe 100644
--- a/algo/ca/core/tracking/CaTrackFit.h
+++ b/algo/ca/core/tracking/CaTrackFit.h
@@ -16,6 +16,9 @@
#include "CaMeasurementXy.h"
#include "CaSimd.h"
#include "KfTrackParam.h"
+#include "CaUtils.h"
+
+#include <type_traits>
namespace cbm::algo::ca
{
@@ -28,18 +31,19 @@ namespace cbm::algo::ca
} // namespace
class Hit;
- template<typename DataT>
- class Station;
-
/// Track fit utilities for the CA tracking based on the Kalman Filter
///
+ template<typename DataT>
class TrackFit {
+ using DataTscal = utils::scaltype<DataT>;
+ using DataTmask = utils::masktype<DataT>;
+
public:
TrackFit() = default;
- TrackFit(const TrackParamV& t) { SetTrack(t); }
+ TrackFit(const TrackParamBase<DataT>& t) { SetTrack(t); }
template<typename T>
TrackFit(const TrackParamBase<T>& t)
@@ -47,10 +51,15 @@ namespace cbm::algo::ca
SetTrack(t);
}
- void SetMask(const fmask& m)
+ void SetMask(const DataTmask& m)
{
- fMask = m;
- fMaskF = iif(m, fvec::One(), fvec::Zero());
+ fMask = m;
+ if constexpr (std::is_same<DataT, fvec>::value) {
+ fMaskF = iif(m, fvec::One(), fvec::Zero());
+ }
+ else {
+ fMaskF = (m ? 1. : 0.);
+ }
}
template<typename T>
@@ -60,11 +69,11 @@ namespace cbm::algo::ca
fQp0 = fTr.GetQp();
}
- void SetQp0(fvec qp0) { fQp0 = qp0; }
+ void SetQp0(DataT qp0) { fQp0 = qp0; }
- TrackParamV& Tr() { return fTr; }
+ TrackParamBase<DataT>& Tr() { return fTr; }
- fvec& Qp0() { return fQp0; }
+ DataT& Qp0() { return fQp0; }
void SetDoFitVelocity(bool v) { fDoFitVelocity = v; }
@@ -76,41 +85,41 @@ namespace cbm::algo::ca
/// Fit utilities
/// set particle mass for the fit
- void SetParticleMass(fvec mass)
+ void SetParticleMass(DataT mass)
{
fMass = mass;
fMass2 = mass * mass;
}
/// get the particle mass
- fvec GetParticleMass() const { return fMass; }
+ DataT GetParticleMass() const { return fMass; }
/// get the particle mass squared
- fvec GetParticleMass2() const { return fMass2; }
+ DataT GetParticleMass2() const { return fMass2; }
/// set max extrapolation step [cm]
- void SetMaxExtrapolationStep(fscal step) { fMaxExtraplationStep = fvec(step); }
+ void SetMaxExtrapolationStep(DataTscal step) { fMaxExtraplationStep = DataT(step); }
/// get the particle mass
- fvec GetMaxExtrapolationStep() const { return fMaxExtraplationStep; }
+ DataT GetMaxExtrapolationStep() const { return fMaxExtraplationStep; }
/// filter the track with the 1d measurement
- void Filter1d(const ca::MeasurementU<fvec>& m);
+ void Filter1d(const ca::MeasurementU<DataT>& m);
/// filter the track with the XY measurement
- void FilterXY(const ca::MeasurementXy<fvec>& m, bool skipUnmeasuredCoordinates = false);
+ void FilterXY(const ca::MeasurementXy<DataT>& m, bool skipUnmeasuredCoordinates = false);
/// filter the track with the time measurement
- void FilterTime(fvec t, fvec dt2, fvec timeInfo);
+ void FilterTime(DataT t, DataT dt2, const DataTmask& m);
/// filter the track with the time measurement
- void FilterTime(MeasurementTime<fvec> mt) { FilterTime(mt.T(), mt.Dt2(), mt.NdfT()); }
+ void FilterTime(MeasurementTime<DataT> mt) { FilterTime(mt.T(), mt.Dt2(), DataTmask(mt.NdfT() > DataT(0.))); }
/// filter the track with the hit
- void FilterHit(const ca::Station<fvec>& s, const ca::Hit& h);
+ void FilterHit(const ca::Hit& h, const DataTmask& timeInfo);
/// filter the inverse speed
- void FilterVi(fvec vi);
+ void FilterVi(DataT vi);
/// measure the track velocity with the track Qp and the mass
void MeasureVelocityWithQp();
@@ -119,22 +128,22 @@ namespace cbm::algo::ca
/// it can do several extrapolation steps if the Z is far away
/// \param z - Z coordinate to extrapolate to
/// \param F - field region
- void Extrapolate(fvec z, const ca::FieldRegion<fvec>& F);
+ void Extrapolate(DataT z, const ca::FieldRegion<DataT>& F);
/// extrapolate the track to the given Z using the field F
/// it does extrapolation in one step
- void ExtrapolateStep(fvec z_out, const ca::FieldRegion<fvec>& F);
+ void ExtrapolateStep(DataT z_out, const ca::FieldRegion<DataT>& F);
/// extrapolate the track to the given Z using linearization at the straight line
- void ExtrapolateLine(fvec z_out, const ca::FieldRegion<fvec>& F);
+ void ExtrapolateLine(DataT z_out, const ca::FieldRegion<DataT>& F);
/// extrapolate the track to the given Z assuming no magnetic field
- void ExtrapolateLineNoField(fvec z_out);
+ void ExtrapolateLineNoField(DataT z_out);
/// apply energy loss correction to the track
/// \param radThick - radiation length of the material
/// \param upstreamDirection - direction of the track (1 or 0)
- void EnergyLossCorrection(fvec radThick, fvec upstreamDirection);
+ void EnergyLossCorrection(DataT radThick, DataT upstreamDirection);
/// apply energy loss correction to the track
/// more accurate formula using material atomic numbers
@@ -144,26 +153,26 @@ namespace cbm::algo::ca
/// \param radLen - radiation length of the material
/// \param radThick - radiation length of the material
/// \param upstreamDirection - direction of the track (1 or 0)
- void EnergyLossCorrection(int atomicZ, fscal atomicA, fscal rho, fscal radLen, fvec radThick,
- fvec upstreamDirection);
+ void EnergyLossCorrection(int atomicZ, DataTscal atomicA, DataTscal rho, DataTscal radLen, DataT radThick,
+ DataT upstreamDirection);
/// apply multiple scattering correction to the track with the given Qp0
- void MultipleScattering(fvec radThick, fvec qp0);
+ void MultipleScattering(DataT radThick, DataT qp0);
/// apply multiple scattering correction to the track
- void MultipleScattering(fvec radThick) { MultipleScattering(radThick, fQp0); }
+ void MultipleScattering(DataT radThick) { MultipleScattering(radThick, fQp0); }
/// apply multiple scattering correction in thick material to the track
- void MultipleScatteringInThickMaterial(fvec radThick, fvec thickness, bool fDownstream);
+ void MultipleScatteringInThickMaterial(DataT radThick, DataT thickness, bool fDownstream);
///------------------------------------------------------------------
/// special utilities needed by the combinatorial track finder
/// extrapolate track as a line, return the extrapolated X, Y and the Jacobians
- void GetExtrapolatedXYline(fvec z, const ca::FieldRegion<fvec>& F, fvec& extrX, fvec& extrY,
- std::array<fvec, TrackParamV::kNtrackParam>& Jx,
- std::array<fvec, TrackParamV::kNtrackParam>& Jy) const;
+ void GetExtrapolatedXYline(DataT z, const kf::FieldRegion<DataT>& F, DataT& extrX, DataT& extrY,
+ std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jx,
+ std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jy) const;
/// filter the track with the XY measurement placed at different Z
/// \param m - measurement
@@ -171,32 +180,33 @@ namespace cbm::algo::ca
/// \param extrY - extrapolated Y of the track
/// \param Jx - Jacobian of the extrapolated X
/// \param Jy - Jacobian of the extrapolated Y
- void FilterExtrapolatedXY(const ca::MeasurementXy<fvec>& m, fvec extrX, fvec extrY,
- const std::array<fvec, TrackParamV::kNtrackParam>& Jx,
- const std::array<fvec, TrackParamV::kNtrackParam>& Jy);
+ void FilterExtrapolatedXY(const kf::MeasurementXy<DataT>& m, DataT extrX, DataT extrY,
+ const std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jx,
+ const std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jy);
/// extrapolate the track to the given Z using linearization at the straight line,
/// \param z_out - Z coordinate to extrapolate to
/// \return pair of the extrapolated X, and dX2 - the rms^2 of the extrapolated x
- std::pair<fvec, fvec> ExtrapolateLineXdX2(fvec z_out) const;
+ std::pair<DataT, DataT> ExtrapolateLineXdX2(DataT z_out) const;
/// extrapolate the track to the given Z using linearization at the straight line,
/// \param z_out - Z coordinate to extrapolate to
/// \return pair of the extrapolated Y, and dY2 - the rms^2 of the extrapolated y
- std::pair<fvec, fvec> ExtrapolateLineYdY2(fvec z_out) const;
+ std::pair<DataT, DataT> ExtrapolateLineYdY2(DataT z_out) const;
/// extrapolate the track to the given Z using linearization at the straight line,
/// \param z_out - Z coordinate to extrapolate to
/// \return extrapolated correlation cov<x,y>
- fvec ExtrapolateLineDxy(fvec z_out) const;
+ DataT ExtrapolateLineDxy(DataT z_out) const;
/// add target measuremet to the track using linearisation at a straight line
- void FilterWithTargetAtLine(fvec targZ, const ca::MeasurementXy<fvec>& targXYInfo, const ca::FieldRegion<fvec>& F);
+ void FilterWithTargetAtLine(DataT targZ, const ca::MeasurementXy<DataT>& targXYInfo,
+ const ca::FieldRegion<DataT>& F);
/// \brief Approximate mean energy loss with Bethe-Bloch formula
/// \param bg2 (beta*gamma)^2
/// \return mean energy loss
- static fvec ApproximateBetheBloch(fvec bg2);
+ static DataT ApproximateBetheBloch(DataT bg2);
/// \brief Approximate mean energy loss with Bethe-Bloch formula
/// \param bg2 (beta*gamma)^2
@@ -206,7 +216,7 @@ namespace cbm::algo::ca
/// \param kp3 mean excitation energy [GeV]
/// \param kp4 mean Z/A
/// \return mean energy loss
- static fvec ApproximateBetheBloch(fvec bg2, fvec kp0, fvec kp1, fvec kp2, fvec kp3, fvec kp4);
+ static DataT ApproximateBetheBloch(DataT bg2, DataT kp0, DataT kp1, DataT kp2, DataT kp3, DataT kp4);
/// \brief git two chi^2 components of the track fit to measurement
/// \param m - measurement
@@ -217,8 +227,8 @@ namespace cbm::algo::ca
/// \param C11 - track covariance C11
/// \return pair of (chi^2_x, chi^2_u) components of the chi^2.
/// chi^2_u is calculated after track is fit to the X measurement
- static std::tuple<fvec, fvec> GetChi2XChi2U(ca::MeasurementXy<fvec> m, fvec x, fvec y, fvec C00, fvec C10,
- fvec C11);
+ static std::tuple<DataT, DataT> GetChi2XChi2U(ca::MeasurementXy<DataT> m, DataT x, DataT y, DataT C00, DataT C10,
+ DataT C11);
/// \brief fast guess of track parameterts based on its hits
/// \param trackZ - Z coordinate of the track
@@ -230,23 +240,25 @@ namespace cbm::algo::ca
/// \param hitW - hit weight
/// \param hitWtime - hit weight for the time measurement
/// \param NHits - number of hits
- void GuessTrack(const fvec& trackZ, const fvec hitX[], const fvec hitY[], const fvec hitZ[], const fvec hitT[],
- const fvec By[], const fmask hitW[], const fmask hitWtime[], int NHits);
+ void GuessTrack(const DataT& trackZ, const DataT hitX[], const DataT hitY[], const DataT hitZ[], const DataT hitT[],
+ const DataT By[], const DataTmask hitW[], const DataTmask hitWtime[], int NHits);
private:
+ typedef const DataT cnst;
+
///--------------------------
/// Data members
- fmask fMask{fmask::One()}; ///< mask of active elements in simd vectors
- fvec fMaskF{fvec::One()}; ///< floating point representation of fMask
+ DataTmask fMask{true}; ///< mask of active elements in simd vectors
+ DataT fMaskF{1.}; ///< floating point representation of fMask
- TrackParamV fTr{}; ///< track parameters
- fvec fQp0{0.};
+ TrackParamBase<DataT> fTr{}; ///< track parameters
+ DataT fQp0{0.};
- fvec fMass{0.10565800}; ///< particle mass (muon mass by default)
- fvec fMass2{fMass * fMass}; ///< mass squared
+ DataT fMass{0.10565800}; ///< particle mass (muon mass by default)
+ DataT fMass2{fMass * fMass}; ///< mass squared
- fvec fMaxExtraplationStep{50.}; ///< max extrapolation step [cm]
+ DataT fMaxExtraplationStep{50.}; ///< max extrapolation step [cm]
bool fDoFitVelocity{0}; // should the track velocity be fitted as an independent parameter
@@ -254,30 +266,38 @@ namespace cbm::algo::ca
// =============================================================================================
- inline std::string TrackFit::ToString(int i) { return fTr.ToString(i); }
+ template<typename DataT>
+ inline std::string TrackFit<DataT>::ToString(int i)
+ {
+ return fTr.ToString(i);
+ }
- inline void TrackFit::SetOneEntry(const int i0, const TrackFit& T1, const int i1)
+ template<typename DataT>
+ inline void TrackFit<DataT>::SetOneEntry(const int i0, const TrackFit& T1, const int i1)
{
fTr.SetOneEntry(i0, T1.fTr, i1);
- fQp0[i0] = T1.fQp0[i1];
+ utils::VecCopy<DataT, DataT, false, false>::CopyEntries(fQp0, i0, T1.fQp0, i1);
}
- inline std::pair<fvec, fvec> TrackFit::ExtrapolateLineXdX2(fvec z_out) const
+ template<typename DataT>
+ inline std::pair<DataT, DataT> TrackFit<DataT>::ExtrapolateLineXdX2(DataT z_out) const
{
- fvec dz = (z_out - fTr.GetZ());
+ DataT dz = (z_out - fTr.GetZ());
return std::pair(fTr.GetX() + fTr.GetTx() * dz, fTr.C00() + dz * (2 * fTr.C20() + dz * fTr.C22()));
}
- inline std::pair<fvec, fvec> TrackFit::ExtrapolateLineYdY2(fvec z_out) const
+ template<typename DataT>
+ inline std::pair<DataT, DataT> TrackFit<DataT>::ExtrapolateLineYdY2(DataT z_out) const
{
- fvec dz = (z_out - fTr.GetZ());
+ DataT dz = (z_out - fTr.GetZ());
return std::pair(fTr.GetY() + fTr.GetTy() * dz, fTr.C11() + dz * (2 * fTr.C31() + dz * fTr.C33()));
}
- inline fvec TrackFit::ExtrapolateLineDxy(fvec z_out) const
+ template<typename DataT>
+ inline DataT TrackFit<DataT>::ExtrapolateLineDxy(DataT z_out) const
{
- fvec dz = (z_out - fTr.GetZ());
+ DataT dz = (z_out - fTr.GetZ());
return fTr.C10() + dz * (fTr.C21() + fTr.C30() + dz * fTr.C32());
}
diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
index 63f401fae12ccd1707007eb780aa3db0bda5266d..b9dc278886f39b00972c068d3e9bfeb4454c4ec1 100644
--- a/algo/ca/core/tracking/CaTrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -51,7 +51,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
const int nStations = fParameters.GetNstationsActive();
int nTracks_SIMD = fvec::size();
- ca::TrackFit fit; // fit parameters coresponding to the current track
+ ca::TrackFit<fvec> fit; // fit parameters coresponding to the current track
TrackParamV& tr = fit.Tr();
fit.SetParticleMass(fDefaultMass);
fit.SetDoFitVelocity(true);
@@ -269,7 +269,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
- fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
+ fit.FilterTime(time[ista], dt2[ista], fmask(sta[ista].timeInfo));
fldB2 = fldB1;
@@ -372,7 +372,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(-1.f));
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
- fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
+ fit.FilterTime(time[ista], dt2[ista], fmask(sta[ista].timeInfo));
fldB2 = fldB1;
fldZ2 = fldZ1;
diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
index 08a39ca921813511a84ec1e2f776a3188ba1a964..5a1227b5829b71e978e3dd4853e6c0d8002e8db1 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -235,7 +235,7 @@ void TripletConstructor::FitDoublets()
fFit.ExtrapolateLineNoField(z_2);
fFit.FilterXY(m_2);
- fFit.FilterTime(t_2, dt2_2, fStaM->timeInfo);
+ fFit.FilterTime(t_2, dt2_2, fmask(fStaM->timeInfo));
fFit.SetQp0(isMomentumFitted ? fFit.Tr().GetQp() : maxQp);
fFit.MultipleScattering(fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y()));
fFit.SetQp0(fFit.Tr().Qp());
@@ -440,7 +440,7 @@ void TripletConstructor::FitTriplets()
// "triplets", "ev:iter:i0:x0:y0:z0:i1:x1:y1:z1:i2:x2:y2:z2:mc:sta:p:vx:vy:vz:chi2:ndf:chi2time:ndfTime");
}
- ca::TrackFit fit;
+ ca::TrackFit<fvec> fit;
fit.SetMask(fmask::One());
if (frWData.CurrentIteration()->GetElectronFlag()) {
@@ -555,7 +555,7 @@ void TripletConstructor::FitTriplets()
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec(energyLossSign));
fit.FilterXY(mxy[ih]);
- fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
+ fit.FilterTime(t[ih], dt2[ih], fmask(sta[ih].timeInfo));
}
};
@@ -678,12 +678,13 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
const auto& grid = frWData.Grid(iSta);
const fvec maxDZ = frWData.CurrentIteration()->GetMaxDZ();
- ca::GridArea area(grid, T.X()[0], T.Y()[0], (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * abs(T.Tx()))[0],
- (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * abs(T.Ty()))[0]);
+ ca::GridArea area(grid, T.X()[0], T.Y()[0],
+ (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * utils::fabs(T.Tx()))[0],
+ (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * utils::fabs(T.Ty()))[0]);
if constexpr (fDebugCollectHits) {
LOG(info) << "search area: " << T.X()[0] << " " << T.Y()[0] << " "
- << (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * abs(T.Tx()))[0] << " "
- << (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * abs(T.Ty()))[0];
+ << (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * utils::fabs(T.Tx()))[0] << " "
+ << (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * utils::fabs(T.Ty()))[0];
}
if (fParameters.DevIsIgnoreHitSearchAreas()) {
area.DoLoopOverEntireGrid();
@@ -755,7 +756,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
ca::MeasurementXy<fvec> mxy(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
const fvec C10 = fFit.ExtrapolateLineDxy(z);
- const auto [chi2x, chi2u] = ca::TrackFit::GetChi2XChi2U(mxy, x, y, C00, C10, C11);
+ const auto [chi2x, chi2u] = ca::TrackFit<fvec>::GetChi2XChi2U(mxy, x, y, C00, C10, C11);
// TODO: adjust the cut, cut on chi2x & chi2u separately
if (!frWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
diff --git a/algo/ca/core/tracking/CaTripletConstructor.h b/algo/ca/core/tracking/CaTripletConstructor.h
index f5b51a14ba665f87b1ac4fc030bd59b2ef97a059..ba7e26cbd6d2bfa87ead7147299b05892f3e4770 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.h
+++ b/algo/ca/core/tracking/CaTripletConstructor.h
@@ -119,7 +119,7 @@ namespace cbm::algo::ca
Vector<ca::Triplet> fTriplets{"TripletConstructor::fTriplets"};
- ca::TrackFit fFit;
+ ca::TrackFit<fvec> fFit;
private:
static constexpr bool fDebugDublets = false; // print debug info for dublets
diff --git a/algo/kf/core/utils/KfSimdPseudo.h b/algo/kf/core/utils/KfSimdPseudo.h
index bfdb741408860539cf22cde6499571ca63db4e7d..9edf285e77fb36e662f43ff0a6acea34089ec654 100644
--- a/algo/kf/core/utils/KfSimdPseudo.h
+++ b/algo/kf/core/utils/KfSimdPseudo.h
@@ -237,7 +237,7 @@ namespace cbm::algo::kf
friend fvec max(const fvec& a, const fvec& b) { _f2(a, b, max) }
friend fvec asgnb(const fvec& a, const fvec& b) { _f2(a, b, asgnb) }
friend fvec sqrt(const fvec& a) { _f1(a, sqrt) }
- friend fvec abs(const fvec& a) { _f1(a, fabs) }
+ // friend fvec abs(const fvec& a) { _f1(a, fabs) } // disabled to avoid confusions, use fabs() instead
friend fvec sgn(const fvec& a) { _f1(a, sgn) }
friend fvec exp(const fvec& a) { _f1(a, exp) }
friend fvec log(const fvec& a) { _f1(a, log) }
diff --git a/algo/kf/core/utils/KfUtils.h b/algo/kf/core/utils/KfUtils.h
index 560bad4cc62129e26329ee16a30124ec2d38ea79..2f4c277f46804f90ed57b18b623ea0f16f17b41e 100644
--- a/algo/kf/core/utils/KfUtils.h
+++ b/algo/kf/core/utils/KfUtils.h
@@ -20,6 +20,12 @@
namespace cbm::algo::kf::utils
{
+ template<typename T>
+ using scaltype = typename std::conditional<std::is_same<T, fvec>::value, fscal, T>::type;
+
+ template<typename T>
+ using masktype = typename std::conditional<std::is_same<T, fvec>::value, fmask, bool>::type;
+
inline fvec iif(const fmask& m, const fvec& t, const fvec& f) { return Vc::iif(m, t, f); }
inline fvec fabs(const fvec& v) { return Vc::abs(v); }
@@ -35,6 +41,28 @@ namespace cbm::algo::kf::utils
return std::fabs(v);
}
+ template<typename T>
+ inline bool isFull(const T& m)
+ {
+ if constexpr (std::is_same_v<T, fmask>) {
+ return m.isFull();
+ }
+ else {
+ return m;
+ }
+ }
+
+ template<typename T>
+ inline T max(const T& a, const T& b)
+ {
+ if constexpr (std::is_same_v<T, fvec>) {
+ return Vc::max(a, b);
+ }
+ else {
+ return std::max(a, b);
+ }
+ }
+
/// tell the CPU that the bool condition is likely to be true
#if defined(__GNUC__) && __GNUC__ - 0 >= 3
__attribute__((always_inline)) inline bool IsLikely(bool b) { return __builtin_expect(!!(b), 1); }
diff --git a/reco/KF/CbmKfTrackFitter.cxx b/reco/KF/CbmKfTrackFitter.cxx
index 811412fcf685ae6e5eea6ad2432dca739cd6e87d..afe921dc9916082690bb7574f6baa5964a0d21cb 100644
--- a/reco/KF/CbmKfTrackFitter.cxx
+++ b/reco/KF/CbmKfTrackFitter.cxx
@@ -434,10 +434,10 @@ void CbmKfTrackFitter::FilterFirstMeasurement(const FitNode& n)
tr.ResetErrors(mxy.Dx2(), mxy.Dy2(), 100., 100., 1., 1.e4, 1.e2);
tr.SetC10(mxy.Dxy());
- if (!(fSkipUnmeasuredCoordinates && mxy.NdfX()[0] == 0)) {
+ if (!(fSkipUnmeasuredCoordinates && mxy.NdfX() == 0)) {
tr.SetX(mxy.X());
}
- if (!(fSkipUnmeasuredCoordinates && mxy.NdfY()[0] == 0)) {
+ if (!(fSkipUnmeasuredCoordinates && mxy.NdfY() == 0)) {
tr.SetY(mxy.Y());
}
tr.SetChiSq(0.);
@@ -468,14 +468,14 @@ void CbmKfTrackFitter::AddMaterialEffects(CbmKfTrackFitter::FitNode& n, CbmKfTra
// calculate the radiation thickness from the current track
if (!n.fIsRadThickFixed) {
n.fRadThick =
- CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThicknessScal(n.fMaterialLayer, tr.GetX()[0], tr.GetY()[0]);
+ CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThicknessScal(n.fMaterialLayer, tr.GetX(), tr.GetY());
}
- fvec msQp0 = fIsQpForMsFixed ? fDefaultQpForMs : fFit.Qp0();
+ double msQp0 = fIsQpForMsFixed ? fDefaultQpForMs : fFit.Qp0();
fFit.MultipleScattering(n.fRadThick, msQp0);
if (!fIsQpForMsFixed) {
- fFit.EnergyLossCorrection(n.fRadThick, (direction == kUpstream) ? fvec::One() : fvec::Zero());
+ fFit.EnergyLossCorrection(n.fRadThick, direction);
}
}
@@ -527,10 +527,9 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
LOG(warning) << "CbmKfTrackFitter: the nodes are not ordered in Z!";
}
- fFit.SetMask(fmask::One());
fFit.SetParticleMass(fMass);
- ca::FieldRegion<ca::fvec> field _fvecalignment;
+ ca::FieldRegion<double> field;
field.SetUseOriginalField();
{ // fit downstream
@@ -541,7 +540,7 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
if (n.fIsFitted) { // The initial parameters are taken from the previous fit
// Linearisation at the previously fitted momentum
- assert(n.fZ == n.fTrack.GetZ()[0]);
+ assert(n.fZ == n.fTrack.GetZ());
fFit.SetTrack(n.fTrack);
}
else { // The initial parameters are calculated from the first and the last measurements
@@ -562,9 +561,9 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
FilterFirstMeasurement(n);
if (fVerbosityLevel > 0) {
- std::cout << " fit downstream: node " << firstHitNode << " chi2 " << fFit.Tr().GetChiSq()[0] << " x "
- << fFit.Tr().GetX()[0] << " y " << fFit.Tr().GetY()[0] << " tx " << fFit.Tr().GetTx()[0] << " ty "
- << fFit.Tr().GetTy()[0] << std::endl;
+ std::cout << " fit downstream: node " << firstHitNode << " chi2 " << fFit.Tr().GetChiSq() << " x "
+ << fFit.Tr().GetX() << " y " << fFit.Tr().GetY() << " tx " << fFit.Tr().GetTx() << " ty "
+ << fFit.Tr().GetTy() << std::endl;
}
}
@@ -593,9 +592,9 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
}
if (fVerbosityLevel > 0) {
- std::cout << " fit downstream: node " << iNode << " chi2 " << fFit.Tr().GetChiSq()[0] << " x "
- << fFit.Tr().GetX()[0] << " y " << fFit.Tr().GetY()[0] << " tx " << fFit.Tr().GetTx()[0] << " ty "
- << fFit.Tr().GetTy()[0] << std::endl;
+ std::cout << " fit downstream: node " << iNode << " chi2 " << fFit.Tr().GetChiSq() << " x " << fFit.Tr().GetX()
+ << " y " << fFit.Tr().GetY() << " tx " << fFit.Tr().GetTx() << " ty " << fFit.Tr().GetTy()
+ << std::endl;
}
if (iNode >= lastHitNode) { // the track is fully fitted, store it
n.fTrack = fFit.Tr();
@@ -605,8 +604,8 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
}
}
- double dnChi2 = fFit.Tr().GetChiSq()[0];
- double dnNdf = fFit.Tr().GetNdf()[0];
+ double dnChi2 = fFit.Tr().GetChiSq();
+ double dnNdf = fFit.Tr().GetNdf();
{ // fit upstream
@@ -616,9 +615,9 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
fFit.SetTrack(n.fTrack);
FilterFirstMeasurement(n);
if (fVerbosityLevel > 0) {
- std::cout << "\n\n up node " << lastHitNode << " chi2 " << fFit.Tr().GetChiSq()[0] << " x "
- << fFit.Tr().GetX()[0] << " y " << fFit.Tr().GetY()[0] << " tx " << fFit.Tr().GetTx()[0] << " ty "
- << fFit.Tr().GetTy()[0] << std::endl;
+ std::cout << "\n\n up node " << lastHitNode << " chi2 " << fFit.Tr().GetChiSq() << " x " << fFit.Tr().GetX()
+ << " y " << fFit.Tr().GetY() << " tx " << fFit.Tr().GetTx() << " ty " << fFit.Tr().GetTy()
+ << std::endl;
}
}
@@ -636,9 +635,8 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
fFit.FilterTime(n.fMt);
}
if (fVerbosityLevel > 0) {
- std::cout << " up node " << iNode << " chi2 " << fFit.Tr().GetChiSq()[0] << " x " << fFit.Tr().GetX()[0]
- << " y " << fFit.Tr().GetY()[0] << " tx " << fFit.Tr().GetTx()[0] << " ty " << fFit.Tr().GetTy()[0]
- << std::endl;
+ std::cout << " up node " << iNode << " chi2 " << fFit.Tr().GetChiSq() << " x " << fFit.Tr().GetX() << " y "
+ << fFit.Tr().GetY() << " tx " << fFit.Tr().GetTx() << " ty " << fFit.Tr().GetTy() << std::endl;
}
// combine partially fitted downstream and upstream tracks
@@ -664,8 +662,8 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
}
}
- double upChi2 = fFit.Tr().GetChiSq()[0];
- double upNdf = fFit.Tr().GetNdf()[0];
+ double upChi2 = fFit.Tr().GetChiSq();
+ double upNdf = fFit.Tr().GetNdf();
if (!fDoSmooth) {
if (upNdf != dnNdf) {
@@ -696,19 +694,19 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
return ok;
}
-bool CbmKfTrackFitter::Smooth(kf::TrackParamV& t1, const kf::TrackParamV& t2)
+bool CbmKfTrackFitter::Smooth(kf::TrackParamBase<double>& t1, const kf::TrackParamBase<double>& t2)
{
// TODO: move to the CaTrackFit class
constexpr int nPar = kf::TrackParamV::kNtrackParam;
constexpr int nCov = kf::TrackParamV::kNcovParam;
- double r[nPar] = {t1.X()[0], t1.Y()[0], t1.Tx()[0], t1.Ty()[0], t1.Qp()[0], t1.Time()[0], t1.Vi()[0]};
- double m[nPar] = {t2.X()[0], t2.Y()[0], t2.Tx()[0], t2.Ty()[0], t2.Qp()[0], t2.Time()[0], t2.Vi()[0]};
+ double r[nPar] = {t1.X(), t1.Y(), t1.Tx(), t1.Ty(), t1.Qp(), t1.Time(), t1.Vi()};
+ double m[nPar] = {t2.X(), t2.Y(), t2.Tx(), t2.Ty(), t2.Qp(), t2.Time(), t2.Vi()};
double S[nCov], S1[nCov], Tmp[nCov];
for (int i = 0; i < nCov; i++) {
- S[i] = (t1.GetCovMatrix()[i] + t2.GetCovMatrix()[i])[0];
+ S[i] = (t1.GetCovMatrix()[i] + t2.GetCovMatrix()[i]);
}
int nullty = 0;
@@ -729,7 +727,7 @@ bool CbmKfTrackFitter::Smooth(kf::TrackParamV& t1, const kf::TrackParamV& t2)
double Si[nPar][nPar];
for (int i = 0, k = 0; i < nPar; i++) {
for (int j = 0; j <= i; j++, k++) {
- C[i][j] = t1.GetCovMatrix()[k][0];
+ C[i][j] = t1.GetCovMatrix()[k];
Si[i][j] = S1[k];
C[j][i] = C[i][j];
Si[j][i] = Si[i][j];
@@ -753,7 +751,7 @@ bool CbmKfTrackFitter::Smooth(kf::TrackParamV& t1, const kf::TrackParamV& t2)
for (int l = 0; l < nPar; l++) {
kc += K[i][l] * C[l][j];
}
- t1.CovMatrix()[k] = t1.CovMatrix()[k][0] - kc;
+ t1.CovMatrix()[k] = t1.CovMatrix()[k] - kc;
}
}
@@ -786,9 +784,9 @@ bool CbmKfTrackFitter::Smooth(kf::TrackParamV& t1, const kf::TrackParamV& t2)
chi2Time += dzeta[i] * SiDzeta;
}
}
- t1.SetChiSq(chi2 + t1.GetChiSq()[0] + t2.GetChiSq()[0]);
- t1.SetChiSqTime(chi2Time + t1.GetChiSqTime()[0] + t2.GetChiSqTime()[0]);
- t1.SetNdf(5 + t1.GetNdf()[0] + t2.GetNdf()[0]);
- t1.SetNdfTime(2 + t1.GetNdfTime()[0] + t2.GetNdfTime()[0]);
+ t1.SetChiSq(chi2 + t1.GetChiSq() + t2.GetChiSq());
+ t1.SetChiSqTime(chi2Time + t1.GetChiSqTime() + t2.GetChiSqTime());
+ t1.SetNdf(5 + t1.GetNdf() + t2.GetNdf());
+ t1.SetNdfTime(2 + t1.GetNdfTime() + t2.GetNdfTime());
return true;
}
diff --git a/reco/KF/CbmKfTrackFitter.h b/reco/KF/CbmKfTrackFitter.h
index b610d2a748ade45b351689152439560239772500..73786e32167c49b6858754abba7ca807b63f80ce 100644
--- a/reco/KF/CbmKfTrackFitter.h
+++ b/reco/KF/CbmKfTrackFitter.h
@@ -56,7 +56,7 @@ class CbmKfTrackFitter {
double fZ{0.}; ///< Z coordinate
- kf::TrackParamV fTrack{}; ///< fitted track
+ kf::TrackParamBase<double> fTrack{}; ///< fitted track
/// == Material information (if present)
@@ -69,9 +69,9 @@ class CbmKfTrackFitter {
/// == Hit information ( if present )
- ca::MeasurementXy<ca::fvec> fMxy{}; ///< XY-measurement at fZ
+ ca::MeasurementXy<double> fMxy{}; ///< XY-measurement at fZ
- ca::MeasurementTime<ca::fvec> fMt{}; ///< time measurement at fZ
+ ca::MeasurementTime<double> fMt{}; ///< time measurement at fZ
/// == Flags etc
bool fIsXySet{false}; ///< true if the XY measurement is set
@@ -145,7 +145,7 @@ class CbmKfTrackFitter {
void AddMaterialEffects(FitNode& n, Direction direction);
// combine two tracks
- bool Smooth(kf::TrackParamV& t1, const kf::TrackParamV& t2);
+ bool Smooth(kf::TrackParamBase<double>& t1, const kf::TrackParamBase<double>& t2);
private:
// input data arrays
@@ -165,7 +165,7 @@ class CbmKfTrackFitter {
bool fIsInitialized = {false}; // is the fitter initialized
bool fSkipUnmeasuredCoordinates{false};
- ca::TrackFit fFit; // track fit object
+ ca::TrackFit<double> fFit; // track fit object
/// externally defined inverse momentum for the Multiple Scattering calculation.
/// It is used for the tracks in field-free regions.
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index 99343d32b1d120837f3d8ae1c8789a8270ccd097..ad8f858906b9105e17c15198e733ad8bb7c3ba20 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -42,8 +42,6 @@
using namespace cbm::algo::ca;
using cbm::algo::kf::TrackParamD;
-//typedef ca::TrackFit1 ca::TrackFit;
-
using std::vector;
using namespace std;
@@ -152,7 +150,7 @@ inline int CbmL1PFFitter::GetStsStationIndex(const CbmStsHit* hit)
}
-void FilterFirst(ca::TrackFit& fit, ca::MeasurementXy<fvec>& mxy, fvec& t, fvec& dt2)
+void FilterFirst(ca::TrackFit<fvec>& fit, ca::MeasurementXy<fvec>& mxy, fvec& t, fvec& dt2)
{
TrackParamV& tr = fit.Tr();
tr.ResetErrors(mxy.Dx2(), mxy.Dy2(), 1., 1., 1., dt2, 1.e2);
@@ -178,7 +176,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
static int nHits = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive();
int nTracks_SIMD = fvec::size();
- ca::TrackFit fit;
+ ca::TrackFit<fvec> fit;
fit.SetParticleMass(CbmL1::Instance()->fpAlgo->GetDefaultParticleMass());
TrackParamV& T = fit.Tr(); // fitting parametr coresponding to current track
@@ -380,7 +378,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.SetMask(initialised && w[i]);
fit.FilterXY(mxy[i]);
- fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
+ fit.FilterTime(t[i], dt2[i], fmask(sta[i].timeInfo));
b2 = b1;
z2 = z1;
@@ -441,7 +439,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.SetMask(initialised && w[i]);
fit.FilterXY(mxy[i]);
- fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
+ fit.FilterTime(t[i], dt2[i], fmask(sta[i].timeInfo));
b2 = b1;
z2 = z1;
@@ -505,7 +503,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
int nTracks_SIMD = fvec::size();
- ca::TrackFit fit;
+ ca::TrackFit<fvec> fit;
TrackParamV& T = fit.Tr(); // fitting parametr coresponding to current track
CbmStsTrack* tr[fvec::size()]{nullptr};
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index b5766f6efd44b3294cdb1db21fb2b947f28080fb..ff8367828d5352a807346139451d9a8ca467b6c1 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1079,7 +1079,7 @@ void CbmL1::TrackFitPerformance()
static bool first_call = 1;
- ca::TrackFit fit;
+ ca::TrackFit<fvec> fit;
fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
fit.SetMask(fmask::One());
//fit.SetMaxExtrapolationStep(10.);
@@ -1506,7 +1506,7 @@ void CbmL1::TrackFitPerformance()
void CbmL1::FillFitHistos(TrackParamV& track, const cbm::ca::tools::MCPoint& mcP, bool isTimeFitted, TH1F* h[])
{
- ca::TrackFit fit;
+ ca::TrackFit<fvec> fit;
fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
fit.SetMask(fmask::One());
//fit.SetMaxExtrapolationStep(10.);
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index e4a646149649aeef6d93fe04f24a7913c88635dd..a56a5303aef0eed2f298c81ac77ac02fc95e51cc 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -445,7 +445,7 @@ void OutputQa::ExecQa()
// ** Fill distributions for reconstructed tracks (including ghost ones) **
// ************************************************************************
- for (size_t iTrkReco = 0; iTrkReco < static_cast<size_t>(nRecoTracks); ++iTrkReco) {
+ for (int iTrkReco = 0; iTrkReco < nRecoTracks; ++iTrkReco) {
const auto& recoTrk = fvRecoTracks[iTrkReco];
// Reject tracks, which do not contain hits
diff --git a/reco/L1/qa/CbmCaTrackFitQa.cxx b/reco/L1/qa/CbmCaTrackFitQa.cxx
index d50786a97e6580b469d95d034a1b89618b593f4c..357bae8080ea908c543c9b2b2026455871cc09c2 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.cxx
+++ b/reco/L1/qa/CbmCaTrackFitQa.cxx
@@ -160,7 +160,7 @@ void TrackFitQa::Init()
void TrackFitQa::Fill(const TrackParamV& trPar, const tools::MCPoint& mcPoint, bool bTimeMeasured, double /*weight*/)
{
// Probably, a bottleneck
- ca::TrackFit fitter;
+ ca::TrackFit<ca::fvec> fitter;
fitter.SetParticleMass(fMass);
fitter.SetMask(fmask::One());
fitter.SetDoFitVelocity(true);
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.h b/reco/L1/qa/CbmCaTrackTypeQa.h
index 8b4ab81f5f5e11405efc5267705f079f19056e57..393bbcae55ff9ed41178d98b900644f8260cd31c 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.h
+++ b/reco/L1/qa/CbmCaTrackTypeQa.h
@@ -278,7 +278,7 @@ namespace cbm::ca
TProfile* fph_stations_point = nullptr; ///< Average number of stations with MC point
TProfile* fph_stations_hit = nullptr; ///< Average number of stations with hit
- ca::TrackFit fTrackFit; ///< Track fitter
+ ca::TrackFit<ca::fvec> fTrackFit; ///< Track fitter
ca::FieldRegion<ca::fvec> fFieldRegion; ///< Magnetic field
int fCounterMC = 0; ///< Counter of MC tracks
@@ -293,9 +293,9 @@ namespace cbm::ca
TString fsTitle = ""; ///< Title of the track category
// TODO: SZh 20.03.2024: Maybe replace CbmL1Track with CaTrack?
- ca::Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to vector of reconstructed tracks
- ca::Vector<CbmL1HitDebugInfo>* fpvHits = nullptr; ///< Pointer to vector of reconstructed hits
- tools::MCData* fpMCData = nullptr; ///< Pointer to MC data object
+ ca::Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to vector of reconstructed tracks
+ ca::Vector<CbmL1HitDebugInfo>* fpvHits = nullptr; ///< Pointer to vector of reconstructed hits
+ tools::MCData* fpMCData = nullptr; ///< Pointer to MC data object
std::shared_ptr<ca::Parameters<ca::fvec>> fpParameters = nullptr; ///< Pointer to parameters object
diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index fb0a71ea5e2db9b65c7bcd459bb6d18848edfa11..536a380bbcb92c56337d8e4c70850ac79ecb8a02 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/reco/alignment/CbmBbaAlignmentTask.cxx
@@ -323,8 +323,8 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
const auto& par = t.fUnalignedTrack.fNodes.front().fTrack;
if (t.fNstsHits < l1Par.GetNstationsActive(ca::EDetectorID::kSts)) continue;
- if (!std::isfinite((ca::fscal) par.GetQp()[0])) continue;
- if (fabs(par.GetQp()[0]) > 1.) continue; // select tracks with min 1 GeV momentum
+ if (!std::isfinite(par.GetQp())) continue;
+ if (fabs(par.GetQp()) > 1.) continue; // select tracks with min 1 GeV momentum
t.fAlignedTrack = t.fUnalignedTrack;
fTracks.push_back(t);
@@ -355,7 +355,7 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
for (auto& n : t.fUnalignedTrack.fNodes) {
n.fIsTimeSet = false;
if (n.fHitSystemId == ECbmModuleId::kTrd) {
- if (n.fMxy.Dx2()[0] < n.fMxy.Dy2()[0]) {
+ if (n.fMxy.Dx2() < n.fMxy.Dy2()) {
n.fMxy.SetNdfY(0);
}
else {
@@ -413,12 +413,12 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
isGood = false;
break;
}
- double x_residual = node.fMxy.X()[0] - node.fTrack.X()[0]; // this should be the difference vector
- double y_residual = node.fMxy.Y()[0] - node.fTrack.Y()[0];
- double x_pull = x_residual / sqrt(node.fMxy.Dx2()[0] + node.fTrack.GetCovariance(0, 0)[0]);
- double y_pull = y_residual / sqrt(node.fMxy.Dy2()[0] + node.fTrack.GetCovariance(1, 1)[0]);
- if (!node.fIsFitted || (node.fMxy.NdfX()[0] > 0. && fabs(x_pull) > cutX)
- || (node.fMxy.NdfY()[0] > 0. && fabs(y_pull) > cutY)) {
+ double x_residual = node.fMxy.X() - node.fTrack.X(); // this should be the difference vector
+ double y_residual = node.fMxy.Y() - node.fTrack.Y();
+ double x_pull = x_residual / sqrt(node.fMxy.Dx2() + node.fTrack.GetCovariance(0, 0));
+ double y_pull = y_residual / sqrt(node.fMxy.Dy2() + node.fTrack.GetCovariance(1, 1));
+ if (!node.fIsFitted || (node.fMxy.NdfX() > 0. && fabs(x_pull) > cutX)
+ || (node.fMxy.NdfY() > 0. && fabs(y_pull) > cutY)) {
isGood = false;
break;
}
@@ -505,8 +505,8 @@ void CbmBbaAlignmentTask::ApplyAlignment(const std::vector<double>& par)
double dy = parConstrained[3 * s.fAlignmentBody + 1];
double dz = parConstrained[3 * s.fAlignmentBody + 2];
- nodeAligned.fMxy.SetX(node.fMxy.X() + ca::fvec(dx));
- nodeAligned.fMxy.SetY(node.fMxy.Y() + ca::fvec(dy));
+ nodeAligned.fMxy.SetX(node.fMxy.X() + dx);
+ nodeAligned.fMxy.SetY(node.fMxy.Y() + dy);
nodeAligned.fZ = node.fZ + dz;
}
}
@@ -546,8 +546,8 @@ double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
auto& t = fTracks[itr];
if (!t.fIsActive) continue;
bool ok = fit.FitTrack(t.fAlignedTrack);
- double chi2 = t.fAlignedTrack.fNodes.back().fTrack.GetChiSq()[0];
- double ndf = t.fAlignedTrack.fNodes.back().fTrack.GetNdf()[0];
+ double chi2 = t.fAlignedTrack.fNodes.back().fTrack.GetChiSq();
+ double ndf = t.fAlignedTrack.fNodes.back().fTrack.GetNdf();
if (!ok || !std::isfinite(chi2) || chi2 <= 0. || ndf <= 0.) {
LOG(fatal) << "BBA: fit failed! ";
assert(0);
@@ -822,8 +822,8 @@ void CbmBbaAlignmentTask::Finish()
fCostInitial = CostFunction(parMisaligned);
for (auto& t : fTracks) {
- double ndf = t.fAlignedTrack.fNodes.back().fTrack.GetNdf()[0];
- double chi2 = t.fAlignedTrack.fNodes.back().fTrack.GetChiSq()[0];
+ double ndf = t.fAlignedTrack.fNodes.back().fTrack.GetNdf();
+ double chi2 = t.fAlignedTrack.fNodes.back().fTrack.GetChiSq();
if (ndf < 0. || chi2 < 0. || !std::isfinite(ndf) || !std::isfinite(chi2)) {
t.fIsActive = false;
}
@@ -847,16 +847,16 @@ void CbmBbaAlignmentTask::Finish()
Sensor& s = fSensors[node.fReference1];
- double x_residual = node.fMxy.X()[0] - node.fTrack.X()[0]; // this should be the difference vector
- double y_residual = node.fMxy.Y()[0] - node.fTrack.Y()[0];
- double x_pull = x_residual / sqrt(node.fMxy.Dx2()[0] + node.fTrack.GetCovariance(0, 0)[0]);
- double y_pull = y_residual / sqrt(node.fMxy.Dy2()[0] + node.fTrack.GetCovariance(1, 1)[0]);
+ double x_residual = node.fMxy.X() - node.fTrack.X(); // this should be the difference vector
+ double y_residual = node.fMxy.Y() - node.fTrack.Y();
+ double x_pull = x_residual / sqrt(node.fMxy.Dx2() + node.fTrack.GetCovariance(0, 0));
+ double y_pull = y_residual / sqrt(node.fMxy.Dy2() + node.fTrack.GetCovariance(1, 1));
- if (node.fMxy.NdfX()[0] > 0) {
+ if (node.fMxy.NdfX() > 0) {
hResidualsBeforeAlignmentX[0]->Fill(x_residual);
hPullsBeforeAlignmentX[0]->Fill(x_pull);
}
- if (node.fMxy.NdfY()[0] > 0) {
+ if (node.fMxy.NdfY() > 0) {
hResidualsBeforeAlignmentY[0]->Fill(y_residual);
hPullsBeforeAlignmentY[0]->Fill(y_pull);
}
@@ -997,16 +997,16 @@ void CbmBbaAlignmentTask::Finish()
Sensor& s = fSensors[node.fReference1];
- double x_residual = node.fMxy.X()[0] - node.fTrack.X()[0];
- double y_residual = node.fMxy.Y()[0] - node.fTrack.Y()[0];
- double x_pull = x_residual / sqrt(node.fMxy.Dx2()[0] + node.fTrack.GetCovariance(0, 0)[0]);
- double y_pull = y_residual / sqrt(node.fMxy.Dy2()[0] + node.fTrack.GetCovariance(1, 1)[0]);
+ double x_residual = node.fMxy.X() - node.fTrack.X();
+ double y_residual = node.fMxy.Y() - node.fTrack.Y();
+ double x_pull = x_residual / sqrt(node.fMxy.Dx2() + node.fTrack.GetCovariance(0, 0));
+ double y_pull = y_residual / sqrt(node.fMxy.Dy2() + node.fTrack.GetCovariance(1, 1));
- if (node.fMxy.NdfX()[0] > 0) {
+ if (node.fMxy.NdfX() > 0) {
hResidualsAfterAlignmentX[0]->Fill(x_residual);
hPullsAfterAlignmentX[0]->Fill(x_pull);
}
- if (node.fMxy.NdfY()[0] > 0) {
+ if (node.fMxy.NdfY() > 0) {
hResidualsAfterAlignmentY[0]->Fill(y_residual);
hPullsAfterAlignmentY[0]->Fill(y_pull);
}
diff --git a/reco/qa/CbmRecoQaTask.cxx b/reco/qa/CbmRecoQaTask.cxx
index 655ab5de0fb0e9b16aaafea71e3711b83a7ad782..57ddfe9a6b4afccf381f9a32d012186c5e1adfb7 100644
--- a/reco/qa/CbmRecoQaTask.cxx
+++ b/reco/qa/CbmRecoQaTask.cxx
@@ -293,11 +293,10 @@ bool CbmRecoQaTask::Detector::View::Load(const CbmHit* h, const CbmEvent* ev)
bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::FitNode* n)
{
- double dx = n->fMxy.X()[0] - n->fTrack.X()[0], dy = n->fMxy.Y()[0] - n->fTrack.Y()[0],
- dt = n->fMt.T()[0] - n->fTrack.Time()[0],
- pullx = dx / sqrt(n->fMxy.Dx2()[0] + n->fTrack.GetCovariance(0, 0)[0]),
- pully = dy / sqrt(n->fMxy.Dy2()[0] + n->fTrack.GetCovariance(1, 1)[0]);
- // pullt = dt / sqrt(n->fMt.Dt2()[0] + n->fTrack.GetCovariance(5, 5)[0]);
+ double dx = n->fMxy.X() - n->fTrack.X(), dy = n->fMxy.Y() - n->fTrack.Y(), dt = n->fMt.T() - n->fTrack.Time(),
+ pullx = dx / sqrt(n->fMxy.Dx2() + n->fTrack.GetCovariance(0, 0)),
+ pully = dy / sqrt(n->fMxy.Dy2() + n->fTrack.GetCovariance(1, 1));
+ // pullt = dt / sqrt(n->fMt.Dt2() + n->fTrack.GetCovariance(5, 5));
for (auto& projection : fProjection) {
int scale = get<0>(projection.second);
@@ -305,13 +304,13 @@ bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::FitNode* n)
if (!hh) continue;
switch (projection.first) {
- case eViewProjection::kXYa: hh->Fill(n->fMxy.X()[0], n->fMxy.Y()[0]); break;
- case eViewProjection::kXYp: hh->Fill(n->fTrack.X()[0], n->fTrack.Y()[0]); break;
- case eViewProjection::kXdX: hh->Fill(n->fMxy.X()[0], scale * dx); break;
- case eViewProjection::kYdY: hh->Fill(n->fMxy.Y()[0], scale * dy); break;
- case eViewProjection::kWdT: hh->Fill(n->fMxy.X()[0], scale * dt); break;
- case eViewProjection::kXpX: hh->Fill(n->fMxy.X()[0], pullx); break;
- case eViewProjection::kYpY: hh->Fill(n->fMxy.Y()[0], pully); break;
+ case eViewProjection::kXYa: hh->Fill(n->fMxy.X(), n->fMxy.Y()); break;
+ case eViewProjection::kXYp: hh->Fill(n->fTrack.X(), n->fTrack.Y()); break;
+ case eViewProjection::kXdX: hh->Fill(n->fMxy.X(), scale * dx); break;
+ case eViewProjection::kYdY: hh->Fill(n->fMxy.Y(), scale * dy); break;
+ case eViewProjection::kWdT: hh->Fill(n->fMxy.X(), scale * dt); break;
+ case eViewProjection::kXpX: hh->Fill(n->fMxy.X(), pullx); break;
+ case eViewProjection::kYpY: hh->Fill(n->fMxy.Y(), pully); break;
default: break;
}
}
@@ -393,7 +392,7 @@ void CbmRecoQaTask::Exec(Option_t*)
auto det = fDetQa[n.fHitSystemId];
// det.Print();
- auto view = det.FindView(n.fMxy.X()[0], n.fMxy.Y()[0], n.fZ);
+ auto view = det.FindView(n.fMxy.X(), n.fMxy.Y(), n.fZ);
if (!view) {
LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
continue;
@@ -431,8 +430,8 @@ void CbmRecoQaTask::Exec(Option_t*)
// mod->Load(&n);
//
// if (mod->hdXx.second)
- // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X()[0], mod->hdXx.first * n.fTrack.Y()[0]);
- // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X()[0], n.fTrack.Time()[0]);
+ // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X(), mod->hdXx.first * n.fTrack.Y());
+ // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X(), n.fTrack.Time());
}
itrack++;
} // END TRACK LOOP
@@ -505,7 +504,7 @@ void CbmRecoQaTask::Exec(Option_t*)
auto det = fDetQa[n.fHitSystemId];
// det.Print();
- auto view = det.FindView(n.fMxy.X()[0], n.fMxy.Y()[0], n.fZ);
+ auto view = det.FindView(n.fMxy.X(), n.fMxy.Y(), n.fZ);
if (!view) {
LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
continue;
@@ -541,8 +540,8 @@ void CbmRecoQaTask::Exec(Option_t*)
// View* mod = &fProjs[n.fReference1 + nnodes][0];
// mod->Load(&n);
// if (mod->hdXx.second)
- // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X()[0], mod->hdXx.first * n.fTrack.Y()[0]);
- // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X()[0], n.fTrack.Time()[0]);
+ // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X(), mod->hdXx.first * n.fTrack.Y());
+ // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X(), n.fTrack.Time());
}
itrack++;
} // END TRACK LOOP