diff --git a/algo/kf/core/KfTrackKalmanFilter.cxx b/algo/kf/core/KfTrackKalmanFilter.cxx
index f9bdb65f257a85338e1a043eb0c2b88ec3e73a62..7230a22f136d33fea3e1a0945b7d2024cb3a54a2 100644
--- a/algo/kf/core/KfTrackKalmanFilter.cxx
+++ b/algo/kf/core/KfTrackKalmanFilter.cxx
@@ -4,6 +4,7 @@
#include "KfTrackKalmanFilter.h"
+#include "AlgoFairloggerCompat.h"
#include "KfMeasurementU.h"
#include "KfMeasurementXy.h"
@@ -46,7 +47,7 @@ namespace cbm::algo::kf
wi = kf::utils::iif(m.Du2() > DataT(0.), wi, DataT(0.));
- fTr.ChiSq() += m.Ndf() * zeta * zeta * wi;
+ fTr.ChiSq() += zeta * zeta * wi;
fTr.Ndf() += kf::utils::iif(fMask, m.Ndf(), DataT(0.));
K1 = F1 * wi;
@@ -977,29 +978,37 @@ namespace cbm::algo::kf
template<typename DataT>
- void TrackKalmanFilter<DataT>::MultipleScattering(DataT radThick, DataT qp0)
+ void TrackKalmanFilter<DataT>::MultipleScattering(DataT radThick, DataT tx, DataT ty, DataT qp)
{
- cnst kONE = 1.;
+ // 1974 - Highland (PDG) correction for multiple scattering
+ // In the formula there is a replacement:
+ // log(X/X0) == log( RadThick * sqrt(1+h) ) -> log(RadThick) + 0.5 * log(1+h);
+ // then we use an approximation:
+ // log(1+h) -> h - h^2/2 + h^3/3 - h^4/4
- DataT tx = fTr.Tx();
- DataT ty = fTr.Ty();
DataT txtx = tx * tx;
DataT tyty = ty * ty;
- DataT txtx1 = txtx + kONE;
+ DataT txtx1 = DataT(1.) + txtx;
+ DataT tyty1 = DataT(1.) + tyty;
DataT h = txtx + tyty;
DataT t = sqrt(txtx1 + tyty);
DataT h2 = h * h;
- DataT qp0t = qp0 * t;
+ DataT qpt = qp * 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;
+ DataT lg = DataT(.0136) * (DataT(1.) + DataT(0.038) * log(radThick * t));
+ lg = kf::utils::iif(lg > DataT(0.), lg, DataT(0.));
- DataT s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ DataT s0 = lg * qp * t;
+ DataT a = (DataT(1.) + fMass2 * qp * qp) * s0 * s0 * t * radThick;
+
+ //cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
+ //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);
+ //DataT a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
fTr.C22() += kf::utils::iif(fMask, txtx1 * a, DataT(0.));
fTr.C32() += kf::utils::iif(fMask, tx * ty * a, DataT(0.));
- fTr.C33() += kf::utils::iif(fMask, (kONE + tyty) * a, DataT(0.));
+ fTr.C33() += kf::utils::iif(fMask, tyty1 * a, DataT(0.));
}
template<typename DataT>
@@ -1087,10 +1096,12 @@ namespace cbm::algo::kf
cnst E2 = fMass2 + p2;
DataT i;
- if (atomicZ < 13)
+ if (atomicZ < 13) {
i = (12. * atomicZ + 7.) * 1.e-9;
- else
+ }
+ else {
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
+ }
cnst bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
diff --git a/algo/kf/core/KfTrackKalmanFilter.h b/algo/kf/core/KfTrackKalmanFilter.h
index 525ef228e98caa00ea649b40c8fd2f810f001c95..ad18cab1061dbbd0d492481c97d5e7286c4e697f 100644
--- a/algo/kf/core/KfTrackKalmanFilter.h
+++ b/algo/kf/core/KfTrackKalmanFilter.h
@@ -151,10 +151,10 @@ namespace cbm::algo::kf
/// apply multiple scattering correction to the track with the given Qp0
- void MultipleScattering(DataT radThick, DataT qp0);
+ void MultipleScattering(DataT radThick, DataT tx0, DataT ty0, DataT qp0);
/// apply multiple scattering correction to the track
- void MultipleScattering(DataT radThick) { MultipleScattering(radThick, fQp0); }
+ void MultipleScattering(DataT radThick) { MultipleScattering(radThick, fTr.GetTx(), fTr.GetTy(), fQp0); }
/// apply multiple scattering correction in thick material to the track
void MultipleScatteringInThickMaterial(DataT radThick, DataT thickness, bool fDownstream);
diff --git a/reco/KF/CbmKfTrackFitter.cxx b/reco/KF/CbmKfTrackFitter.cxx
index 28c2940f4f84443292f70424cf24781d93d0b7fe..fe974c9d069d80d20d803ab5b60fa4f747eb6814 100644
--- a/reco/KF/CbmKfTrackFitter.cxx
+++ b/reco/KF/CbmKfTrackFitter.cxx
@@ -428,7 +428,7 @@ void CbmKfTrackFitter::FilterFirstMeasurement(const FitNode& n)
auto& tr = fFit.Tr();
- tr.ResetErrors(mxy.Dx2(), mxy.Dy2(), 100., 100., 1., 1.e4, 1.e2);
+ tr.ResetErrors(mxy.Dx2(), mxy.Dy2(), 100., 100., 10., 1.e4, 1.e2);
tr.SetC10(mxy.Dxy());
if (!(fSkipUnmeasuredCoordinates && mxy.NdfX() == 0)) {
@@ -437,6 +437,8 @@ void CbmKfTrackFitter::FilterFirstMeasurement(const FitNode& n)
if (!(fSkipUnmeasuredCoordinates && mxy.NdfY() == 0)) {
tr.SetY(mxy.Y());
}
+ tr.SetZ(n.fZ);
+
tr.SetChiSq(0.);
tr.SetChiSqTime(0.);
tr.SetNdf(-5 + mxy.NdfX() + mxy.NdfY());
@@ -462,15 +464,18 @@ void CbmKfTrackFitter::AddMaterialEffects(CbmKfTrackFitter::FitNode& n, kf::FitD
auto& tr = n.fIsFitted ? n.fTrack : fFit.Tr();
+ double msQp0 = n.fIsFitted ? n.fTrack.GetQp() : fFit.Qp0();
+ if (fIsQpForMsFixed) {
+ msQp0 = fDefaultQpForMs;
+ }
+
// calculate the radiation thickness from the current track
if (!n.fIsRadThickFixed) {
n.fRadThick =
CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThicknessScal(n.fMaterialLayer, tr.GetX(), tr.GetY());
}
- double msQp0 = fIsQpForMsFixed ? fDefaultQpForMs : fFit.Qp0();
-
- fFit.MultipleScattering(n.fRadThick, msQp0);
+ fFit.MultipleScattering(n.fRadThick, tr.GetTx(), tr.GetTy(), msQp0);
if (!fIsQpForMsFixed) {
fFit.EnergyLossCorrection(n.fRadThick, direction);
}
@@ -479,7 +484,7 @@ void CbmKfTrackFitter::AddMaterialEffects(CbmKfTrackFitter::FitNode& n, kf::FitD
bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
{
- // fit the track
+ // (re)fit the track
if (fVerbosityLevel > 0) {
std::cout << "FitTrack ... " << std::endl;
}
@@ -521,11 +526,13 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
}
if (!isOrdered) {
- LOG(warning) << "CbmKfTrackFitter: the nodes are not ordered in Z!";
+ LOG(warning) << "CbmKfTrackFitter: track nodes are not ordered in Z!";
}
fFit.SetParticleMass(fMass);
+ // kf::GlobalField::fgOriginalFieldType = kf::EFieldType::Null;
+
kf::FieldRegion<double> field(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
{ // fit downstream
@@ -536,10 +543,16 @@ 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());
+ //if (n.fZ != n.fTrack.GetZ()) {
+ //LOG(ERROR) << "Z " << n.fZ << " track Z " << n.fTrack.GetZ() << " dz " << n.fZ - n.fTrack.GetZ();
+ //}
+ //assert(n.fZ == n.fTrack.GetZ());
fFit.SetTrack(n.fTrack);
}
- else { // The initial parameters are calculated from the first and the last measurements
+ else { // The initial parameters are calculated from the first and the last measurements
+ if (!fDoSmooth) { //SG!!!
+ LOG(fatal) << "CbmKfTrackFitter: Debug: downstream: the first measurement is not fitted";
+ }
auto& tr = fFit.Tr();
auto& n1 = t.fNodes[lastHitNode];
tr.X() = n.fMxy.X();
@@ -558,8 +571,8 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
if (fVerbosityLevel > 0) {
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;
+ << fFit.Tr().GetX() << " y " << fFit.Tr().GetY() << " z " << fFit.Tr().GetZ() << " tx "
+ << fFit.Tr().GetTx() << " ty " << fFit.Tr().GetTy() << std::endl;
}
}
@@ -567,20 +580,25 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
FitNode& n = t.fNodes[iNode];
+ if (!fDoSmooth && !n.fIsFitted) { //SG!!!
+ LOG(fatal) << "CbmKfTrackFitter: Debug: downstream: node " << iNode << " is not fitted";
+ }
+
fFit.Extrapolate(n.fZ, field);
+ if (n.fIsFitted && (iNode <= lastHitNode)) { // linearisation is taken from the previous fit
+ fFit.SetQp0(n.fTrack.GetQp());
+ }
+
if (fDoSmooth) { // store the partially fitted track
n.fTrack = fFit.Tr();
n.fIsFitted = false; // the stored track is not fully fitted
}
- if (n.fIsFitted && (iNode <= lastHitNode)) { // linearisation is taken from the previous fit
- fFit.SetQp0(n.fTrack.GetQp());
- }
-
AddMaterialEffects(n, kf::FitDirection::kDownstream);
if (n.fIsXySet) {
+ assert(iNode <= lastHitNode);
fFit.FilterXY(n.fMxy, fSkipUnmeasuredCoordinates);
}
if (n.fIsTimeSet) {
@@ -589,10 +607,10 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
if (fVerbosityLevel > 0) {
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;
+ << " y " << fFit.Tr().GetY() << " z " << fFit.Tr().GetZ() << " tx " << fFit.Tr().GetTx() << " ty "
+ << fFit.Tr().GetTy() << std::endl;
}
- if (iNode >= lastHitNode) { // the track is fully fitted, store it
+ if (fDoSmooth && iNode >= lastHitNode) { // the track is fully fitted, store it
n.fTrack = fFit.Tr();
n.fIsFitted = true;
fFit.SetQp0(fFit.Tr().GetQp()); // set the linearisation of q/p to the fitted value
@@ -607,22 +625,28 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
{
FitNode& n = t.fNodes[lastHitNode];
+ if (!fDoSmooth && !n.fIsFitted) {
+ LOG(fatal) << "CbmKfTrackFitter: Debug: upstream: the last measurement is not fitted";
+ }
+
assert(n.fIsFitted);
fFit.SetTrack(n.fTrack);
FilterFirstMeasurement(n);
if (fVerbosityLevel > 0) {
- 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;
+ std::cout << "\n\n fit upstream: node " << lastHitNode << " chi2 " << fFit.Tr().GetChiSq() << " x "
+ << fFit.Tr().GetX() << " y " << fFit.Tr().GetY() << " z " << fFit.Tr().GetZ() << " tx "
+ << fFit.Tr().GetTx() << " ty " << fFit.Tr().GetTy() << std::endl;
}
}
for (int iNode = lastHitNode - 1; iNode >= 0; iNode--) {
FitNode& n = t.fNodes[iNode];
- fFit.Extrapolate(n.fZ, field);
+ if (!fDoSmooth && !n.fIsFitted) {
+ LOG(fatal) << "CbmKfTrackFitter: Debug: upstream: the node " << iNode << "is not fitted";
+ }
- AddMaterialEffects(n, kf::FitDirection::kUpstream);
+ fFit.Extrapolate(n.fZ, field);
if (n.fIsXySet) {
fFit.FilterXY(n.fMxy, fSkipUnmeasuredCoordinates);
@@ -630,31 +654,35 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
if (n.fIsTimeSet) {
fFit.FilterTime(n.fMt);
}
+
+ AddMaterialEffects(n, kf::FitDirection::kUpstream);
+
if (fVerbosityLevel > 0) {
- 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;
+ std::cout << " fit upstream: node " << iNode << " chi2 " << fFit.Tr().GetChiSq() << " x " << fFit.Tr().GetX()
+ << " y " << fFit.Tr().GetY() << " z " << fFit.Tr().GetZ() << " tx " << fFit.Tr().GetTx() << " ty "
+ << fFit.Tr().GetTy() << std::endl;
}
// combine partially fitted downstream and upstream tracks
- n.fIsFitted = true;
if (iNode > firstHitNode) {
- n.fIsFitted = false;
if (fDoSmooth) {
n.fIsFitted = Smooth(n.fTrack, fFit.Tr());
}
}
else {
- n.fIsFitted = true;
- n.fTrack = fFit.Tr();
+ if (fDoSmooth) {
+ n.fIsFitted = true;
+ n.fTrack = fFit.Tr();
+ }
}
+
if (fDoSmooth && !n.fIsFitted) {
ok = false;
}
if (n.fIsFitted) {
fFit.SetQp0(n.fTrack.GetQp());
}
- AddMaterialEffects(n, kf::FitDirection::kUpstream);
}
}
@@ -665,9 +693,13 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
if (upNdf != dnNdf) {
LOG(fatal) << "CbmKfTrackFitter: ndf mismatch: dn " << dnNdf << " != up " << upNdf;
}
- if (fabs(upChi2 - dnChi2) > 1.e-2) {
- LOG(fatal) << "CbmKfTrackFitter: chi2 mismatch: dn " << dnChi2 << " != up " << upChi2 << " first node "
- << firstHitNode << " last node " << lastHitNode << " of " << nNodes;
+ if (fabs(upChi2 - dnChi2) > 1.e-1) {
+ //if (upChi2 > dnChi2 + 1.e-2) {
+ LOG(error) << "CbmKfTrackFitter: " << fDebugInfo << ": chi2 mismatch: dn " << dnChi2 << " != up " << upChi2
+ << " first node " << firstHitNode << " last node " << lastHitNode << " of " << nNodes;
+ //if (fVerbosityLevel > 0) {
+ exit(1);
+ //}
}
}
// distribute the final chi2, ndf to all nodes
diff --git a/reco/KF/CbmKfTrackFitter.h b/reco/KF/CbmKfTrackFitter.h
index 00a5f11e97c317c97afecd2748b0b86c1fd4e04a..50f6182c3cefd5123da40fb10b04af60f9c00169 100644
--- a/reco/KF/CbmKfTrackFitter.h
+++ b/reco/KF/CbmKfTrackFitter.h
@@ -51,7 +51,10 @@ class CbmKfTrackFitter {
/// When fIsFitted flag is set (e.g. by the fitter), the node track parameters are used
/// for the trajectory linearisation during the fit.
///
- // TODO: proper interface
+ /// When fitting w/o smoother we assume that smoothed parameters are always better than partially fitted parameters
+ ///
+ /// TODO: proper interface and description
+ ///
struct FitNode {
double fZ{0.}; ///< Z coordinate
@@ -115,7 +118,13 @@ class CbmKfTrackFitter {
void SetSkipUnmeasuredCoordinates(bool skip = true) { fSkipUnmeasuredCoordinates = skip; }
/// set the default inverse momentum for the Multiple Scattering calculation
- void SetDefaultMomentumForMs(double p) { fDefaultQpForMs = 1. / p; }
+ void SetDefaultMomentumForMs(double p) { fDefaultQpForMs = fabs(1. / p); }
+
+ /// set the default inverse momentum for the Multiple Scattering calculation
+ void SetDefaultInverseMomentumForMs(double invP) { fDefaultQpForMs = fabs(invP); }
+
+ /// set the default inverse momentum for the Multiple Scattering calculation
+ void SetNoMultipleScattering() { fDefaultQpForMs = 0.; }
/// fix the inverse momentum for the Multiple Scattering calculation
void FixMomentumForMs(bool fix = true) { fIsQpForMsFixed = fix; }
@@ -134,6 +143,9 @@ class CbmKfTrackFitter {
/// set verbosity level
void SetVerbosityLevel(int level) { fVerbosityLevel = level; }
+ /// set information about the track for debug output
+ void SetDebugInfo(const std::string &info) { fDebugInfo = info; }
+
private:
void FilterFirstMeasurement(const FitNode& n);
@@ -173,4 +185,5 @@ class CbmKfTrackFitter {
bool fIsElectron{false}; // fit track as an electron (with the bermsstrallung effect)
bool fDoSmooth{true}; // do the KF-smoothing to define track pars at all the nodes
int fVerbosityLevel{0}; // verbosity level
+ std::string fDebugInfo{}; // debug info
};