From 4c306063471084a080e5a869e63fed008ab03e9e Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Wed, 6 Mar 2024 02:09:14 +0000
Subject: [PATCH] BBA: ignore an unmeasured coordinate of Trd-1d
---
algo/ca/core/tracking/CaTrackFit.cxx | 10 ++--
algo/ca/core/tracking/CaTrackFit.h | 2 +-
reco/KF/CbmKfTrackFitter.cxx | 26 ++++++---
reco/KF/CbmKfTrackFitter.h | 6 ++-
reco/alignment/CbmBbaAlignmentTask.cxx | 75 ++++++++++++++++----------
reco/alignment/CbmBbaAlignmentTask.h | 11 ++--
6 files changed, 85 insertions(+), 45 deletions(-)
diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx
index f1a33a5c5b..e0118834ef 100644
--- a/algo/ca/core/tracking/CaTrackFit.cxx
+++ b/algo/ca/core/tracking/CaTrackFit.cxx
@@ -186,7 +186,7 @@ namespace cbm::algo::ca
}
- void TrackFit::FilterXY(const ca::MeasurementXy<fvec>& mxy)
+ void TrackFit::FilterXY(const ca::MeasurementXy<fvec>& mxy, bool skipUnmeasuredCoordinates)
{
{
ca::MeasurementU<fvec> mx;
@@ -203,8 +203,12 @@ namespace cbm::algo::ca
mu.SetDu2(mxy.Dy2() - mxy.Dxy() * mxy.Dxy() / mxy.Dx2());
mu.SetNdf(mxy.NdfY());
- Filter1d(mx);
- Filter1d(mu);
+ if (!(skipUnmeasuredCoordinates && mxy.NdfX()[0] == 0)) {
+ Filter1d(mx);
+ }
+ if (!(skipUnmeasuredCoordinates && mxy.NdfY()[0] == 0)) {
+ Filter1d(mu);
+ }
return;
}
diff --git a/algo/ca/core/tracking/CaTrackFit.h b/algo/ca/core/tracking/CaTrackFit.h
index b9c9b5f1e7..f0f1b12924 100644
--- a/algo/ca/core/tracking/CaTrackFit.h
+++ b/algo/ca/core/tracking/CaTrackFit.h
@@ -97,7 +97,7 @@ namespace cbm::algo::ca
void Filter1d(const ca::MeasurementU<fvec>& m);
/// filter the track with the XY measurement
- void FilterXY(const ca::MeasurementXy<fvec>& m);
+ void FilterXY(const ca::MeasurementXy<fvec>& m, bool skipUnmeasuredCoordinates = false);
/// filter the track with the time measurement
void FilterTime(fvec t, fvec dt2, fvec timeInfo);
diff --git a/reco/KF/CbmKfTrackFitter.cxx b/reco/KF/CbmKfTrackFitter.cxx
index 6c163cfc4b..9e61f26bb9 100644
--- a/reco/KF/CbmKfTrackFitter.cxx
+++ b/reco/KF/CbmKfTrackFitter.cxx
@@ -213,10 +213,10 @@ bool CbmKfTrackFitter::CreateGlobalTrack(CbmKfTrackFitter::Track& kfTrack, const
}
// lambda to set the node from the pixel hit
- auto setNode = [&](const CbmPixelHit& h, int stIdx, int hitIdx, ca::EDetectorID detId) {
+ auto setNode = [&](const CbmPixelHit& h, int stIdx, int hitIdx, ca::EDetectorID detId) -> CbmKfTrackFitter::FitNode* {
int ista = caPar.GetStationIndexActive(stIdx, detId);
if (ista < 0) {
- return;
+ return nullptr;
}
assert(ista < nStations);
CbmKfTrackFitter::FitNode& n = t.fNodes[ista];
@@ -238,6 +238,7 @@ bool CbmKfTrackFitter::CreateGlobalTrack(CbmKfTrackFitter::Track& kfTrack, const
n.fHitSystemId = cbm::L1Util::ConvertDetSystemId(detId);
n.fHitAddress = h.GetAddress();
n.fHitIndex = hitIdx;
+ return &n;
};
// Read MVD & STS hits
@@ -368,8 +369,12 @@ bool CbmKfTrackFitter::CreateGlobalTrack(CbmKfTrackFitter::Track& kfTrack, const
return false;
}
const auto& hit = *dynamic_cast<const CbmTrdHit*>(fInputTrdHits->At(hitIndex));
- setNode(hit, CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(&hit), hitIndex,
- ca::EDetectorID::kTrd);
+
+ auto* node = setNode(hit, CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(&hit), hitIndex,
+ ca::EDetectorID::kTrd);
+ if (node != nullptr && hit.GetClassType() == 1) {
+ node->fHitSystemId = ECbmModuleId::kTrd2d;
+ }
}
}
}
@@ -437,8 +442,13 @@ void CbmKfTrackFitter::FilterFirstMeasurement(const FitNode& n)
auto& tr = fFit.Tr();
tr.ResetErrors(mxy.Dx2(), mxy.Dy2(), 1., 1., 1., 1.e4, 1.e2);
tr.SetC10(mxy.Dxy());
- tr.SetX(mxy.X());
- tr.SetY(mxy.Y());
+
+ if (!(fSkipUnmeasuredCoordinates && mxy.NdfX()[0] == 0)) {
+ tr.SetX(mxy.X());
+ }
+ if (!(fSkipUnmeasuredCoordinates && mxy.NdfY()[0] == 0)) {
+ tr.SetY(mxy.Y());
+ }
tr.SetChiSq(0.);
tr.SetChiSqTime(0.);
tr.SetNdf(-5 + mxy.NdfX() + mxy.NdfY());
@@ -517,7 +527,7 @@ void CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
n.fTrack = fFit.Tr();
n.fIsFitted = false;
if (n.fIsXySet) {
- fFit.FilterXY(n.fMxy);
+ fFit.FilterXY(n.fMxy, fSkipUnmeasuredCoordinates);
}
if (n.fIsTimeSet) {
fFit.FilterTime(n.fMt);
@@ -541,7 +551,7 @@ void CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
fFit.Extrapolate(n.fZ, field);
if (n.fIsXySet) {
- fFit.FilterXY(n.fMxy);
+ fFit.FilterXY(n.fMxy, fSkipUnmeasuredCoordinates);
}
if (n.fIsTimeSet) {
fFit.FilterTime(n.fMt);
diff --git a/reco/KF/CbmKfTrackFitter.h b/reco/KF/CbmKfTrackFitter.h
index fa1ec0e32f..5f426fe81c 100644
--- a/reco/KF/CbmKfTrackFitter.h
+++ b/reco/KF/CbmKfTrackFitter.h
@@ -110,6 +110,10 @@ class CbmKfTrackFitter {
/// set electron flag (bremmstrallung will be applied)
void SetElectronFlag(bool isElectron);
+ /// skip unmeasured coordinates
+ void SetSkipUnmeasuredCoordinates(bool skip = true) { fSkipUnmeasuredCoordinates = skip; }
+
+ /// set the input data arrays
bool CreateMvdStsTrack(Track& kfTrack, int stsTrackIndex);
bool CreateGlobalTrack(Track& kfTrack, int globalTrackIndex);
bool CreateGlobalTrack(Track& kfTrack, const CbmGlobalTrack& globalTrack);
@@ -143,7 +147,7 @@ class CbmKfTrackFitter {
//
bool fIsInitialized = {false}; // is the fitter initialized
-
+ bool fSkipUnmeasuredCoordinates{false};
ca::TrackFit fFit; // track fit object
double fMass{ca::constants::phys::PionMass}; // mass hypothesis for the fit
diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index c8d4f6f4c4..4fe3516ad5 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/reco/alignment/CbmBbaAlignmentTask.cxx
@@ -51,17 +51,19 @@ namespace
void CbmBbaAlignmentTask::TrackContainer::MakeConsistent()
{
- fNmvdHits = 0;
- fNstsHits = 0;
- fNmuchHits = 0;
- fNtrdHits = 0;
- fNtofHits = 0;
+ fNmvdHits = 0;
+ fNstsHits = 0;
+ fNmuchHits = 0;
+ fNtrd1dHits = 0;
+ fNtrd2dHits = 0;
+ fNtofHits = 0;
for (const auto& n : fUnalignedTrack.fNodes) {
switch (n.fHitSystemId) {
case ECbmModuleId::kMvd: fNmvdHits++; break;
case ECbmModuleId::kSts: fNstsHits++; break;
case ECbmModuleId::kMuch: fNmuchHits++; break;
- case ECbmModuleId::kTrd: fNtrdHits++; break;
+ case ECbmModuleId::kTrd: fNtrd1dHits++; break;
+ case ECbmModuleId::kTrd2d: fNtrd2dHits++; break;
case ECbmModuleId::kTof: fNtofHits++; break;
default: break;
}
@@ -116,6 +118,7 @@ InitStatus CbmBbaAlignmentTask::Init()
//fNthreads = 1;
fFitter.Init();
+ fFitter.SetSkipUnmeasuredCoordinates(true);
//Get ROOT Manager
FairRootManager* ioman = FairRootManager::Instance();
@@ -197,10 +200,10 @@ InitStatus CbmBbaAlignmentTask::Init()
else if (fTrackingMode == kMcbm) {
if (i == -1) {
- sizeX = 20.;
- sizeY = 20.;
- sizePX = 20.;
- sizePY = 20.;
+ sizeX = 5.;
+ sizeY = 5.;
+ sizePX = 10.;
+ sizePY = 10.;
}
if (i == 0 || i == 1) {
@@ -217,8 +220,8 @@ InitStatus CbmBbaAlignmentTask::Init()
sizePX = 10.;
}
if (i == 4) {
- sizeX = 5.;
- sizeY = 10.;
+ sizeX = 20.;
+ sizeY = 5.;
sizePY = 10.;
}
if (i >= 5) {
@@ -319,12 +322,20 @@ 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]) {
+ n.fMxy.SetNdfY(0);
+ }
+ else {
+ n.fMxy.SetNdfX(0);
+ }
+ }
n.fTrack.SetQp(0.);
}
//if (t.fNstsHits < 1) continue;
//if (t.fNtrdHits < 2) continue;
- if (t.fNstsHits + t.fNmuchHits + t.fNtrdHits + t.fNtofHits < fNtrackingStations - 1) continue;
+ if (t.fNstsHits + t.fNmuchHits + t.fNtrd1dHits + t.fNtrd2dHits + t.fNtofHits < fNtrackingStations - 1) continue;
//if (t.fNtrdHits < 3) continue;
if (t.fNstsHits < 2) continue;
if (t.fNtofHits < 2) continue;
@@ -611,16 +622,21 @@ void CbmBbaAlignmentTask::Finish()
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];
- hResidualsBeforeAlignmentX[0]->Fill(x_residual);
- hResidualsBeforeAlignmentY[0]->Fill(y_residual);
+ 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.fMxy.NdfX()[0] > 0) {
+ hResidualsBeforeAlignmentX[0]->Fill(x_residual);
+ hPullsBeforeAlignmentX[0]->Fill(x_pull);
+ }
+ if (node.fMxy.NdfY()[0] > 0) {
+ hResidualsBeforeAlignmentY[0]->Fill(y_residual);
+ hPullsBeforeAlignmentY[0]->Fill(y_pull);
+ }
+
hResidualsBeforeAlignmentX[s.fTrackingStation + 1]->Fill(x_residual);
hResidualsBeforeAlignmentY[s.fTrackingStation + 1]->Fill(y_residual);
- 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]);
-
- hPullsBeforeAlignmentX[0]->Fill(x_pull);
- hPullsBeforeAlignmentY[0]->Fill(y_pull);
hPullsBeforeAlignmentX[s.fTrackingStation + 1]->Fill(x_pull);
hPullsBeforeAlignmentY[s.fTrackingStation + 1]->Fill(y_pull);
}
@@ -677,16 +693,21 @@ void CbmBbaAlignmentTask::Finish()
double x_residual = node.fMxy.X()[0] - node.fTrack.X()[0];
double y_residual = node.fMxy.Y()[0] - node.fTrack.Y()[0];
- hResidualsAfterAlignmentX[0]->Fill(x_residual);
- hResidualsAfterAlignmentY[0]->Fill(y_residual);
+ 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.fMxy.NdfX()[0] > 0) {
+ hResidualsAfterAlignmentX[0]->Fill(x_residual);
+ hPullsAfterAlignmentX[0]->Fill(x_pull);
+ }
+ if (node.fMxy.NdfY()[0] > 0) {
+ hResidualsAfterAlignmentY[0]->Fill(y_residual);
+ hPullsAfterAlignmentY[0]->Fill(y_pull);
+ }
+
hResidualsAfterAlignmentX[s.fTrackingStation + 1]->Fill(x_residual);
hResidualsAfterAlignmentY[s.fTrackingStation + 1]->Fill(y_residual);
- 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]);
-
- hPullsAfterAlignmentX[0]->Fill(x_pull);
- hPullsAfterAlignmentY[0]->Fill(y_pull);
hPullsAfterAlignmentX[s.fTrackingStation + 1]->Fill(x_pull);
hPullsAfterAlignmentY[s.fTrackingStation + 1]->Fill(y_pull);
}
diff --git a/reco/alignment/CbmBbaAlignmentTask.h b/reco/alignment/CbmBbaAlignmentTask.h
index af78145a7b..6f46befcbb 100644
--- a/reco/alignment/CbmBbaAlignmentTask.h
+++ b/reco/alignment/CbmBbaAlignmentTask.h
@@ -68,11 +68,12 @@ class CbmBbaAlignmentTask : public FairTask {
CbmKfTrackFitter::Track fUnalignedTrack; // track before alignment
CbmKfTrackFitter::Track fAlignedTrack; // track after alignment
- int fNmvdHits{0}; // number of MVD hits
- int fNstsHits{0}; // number of STS hits
- int fNmuchHits{0}; // number of MUCH hits
- int fNtrdHits{0}; // number of TRD hits
- int fNtofHits{0}; // number of TOF hits
+ int fNmvdHits{0}; // number of MVD hits
+ int fNstsHits{0}; // number of STS hits
+ int fNmuchHits{0}; // number of MUCH hits
+ int fNtrd1dHits{0}; // number of TRD hits
+ int fNtrd2dHits{0}; // number of TRD hits
+ int fNtofHits{0}; // number of TOF hits
bool fIsActive{1}; // is the track active
void MakeConsistent();
--
GitLab