From 99d6505f4295c6d35a1da3eb5b52ce6c656f841b Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Thu, 2 Nov 2023 14:11:53 +0100
Subject: [PATCH] CA QA: Track fit QA in track vertex
---
reco/L1/catools/CaToolsMCTrack.cxx | 38 ++++++++++++++++-
reco/L1/catools/CaToolsMCTrack.h | 3 ++
reco/L1/qa/CbmCaTrackFitQa.cxx | 29 +++++++------
reco/L1/qa/CbmCaTrackFitQa.h | 9 ++--
reco/L1/qa/CbmCaTrackTypeQa.cxx | 68 ++++++++++++++++++++++++------
reco/L1/qa/CbmCaTrackTypeQa.h | 6 +++
6 files changed, 123 insertions(+), 30 deletions(-)
diff --git a/reco/L1/catools/CaToolsMCTrack.cxx b/reco/L1/catools/CaToolsMCTrack.cxx
index 36afa39172..abbf2aaf62 100644
--- a/reco/L1/catools/CaToolsMCTrack.cxx
+++ b/reco/L1/catools/CaToolsMCTrack.cxx
@@ -16,9 +16,9 @@
#include "CaToolsMCPoint.h"
+using cbm::ca::tools::MCPoint;
using cbm::ca::tools::MCTrack;
-
// ---------------------------------------------------------------------------------------------------------------------
//
void MCTrack::Clear()
@@ -30,6 +30,42 @@ void MCTrack::Clear()
ClearTouchTrackIndexes();
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+MCPoint MCTrack::GetVertexPoint() const
+{
+ MCPoint point;
+ point.SetCharge(this->GetCharge());
+ point.SetEventId(this->GetEventId());
+ point.SetFileId(this->GetFileId());
+ point.SetExternalId(-1); // not in the array
+ point.SetId(-1); // not in the array
+ point.SetStationId(-1);
+ point.SetMass(this->GetMass());
+ point.SetMotherId(this->GetMotherId());
+ point.SetPdgCode(this->GetPdgCode());
+ point.SetPx(this->GetPx());
+ point.SetPy(this->GetPy());
+ point.SetPz(this->GetPz());
+ point.SetPxIn(this->GetPx());
+ point.SetPyIn(this->GetPy());
+ point.SetPzIn(this->GetPz());
+ point.SetPxOut(this->GetPx());
+ point.SetPyOut(this->GetPy());
+ point.SetPzOut(this->GetPz());
+ point.SetX(this->GetStartX());
+ point.SetY(this->GetStartY());
+ point.SetZ(this->GetStartZ());
+ point.SetXIn(this->GetStartX());
+ point.SetYIn(this->GetStartY());
+ point.SetZIn(this->GetStartZ());
+ point.SetXOut(this->GetStartX());
+ point.SetYOut(this->GetStartY());
+ point.SetZOut(this->GetStartZ());
+ point.SetTime(this->GetStartT());
+ return point;
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void MCTrack::InitHitsInfo(const ca::Vector<CbmL1HitDebugInfo>& vHits)
diff --git a/reco/L1/catools/CaToolsMCTrack.h b/reco/L1/catools/CaToolsMCTrack.h
index 3ec771e864..aeda05ea67 100644
--- a/reco/L1/catools/CaToolsMCTrack.h
+++ b/reco/L1/catools/CaToolsMCTrack.h
@@ -211,6 +211,9 @@ namespace cbm::ca::tools
/// Gets track slope along y-axis
double GetTy() const { return fMom[1] / fMom[2]; }
+ /// \brief Creates an MC point from the track vertex
+ MCPoint GetVertexPoint() const;
+
/// Gets a reference to vector of reconstructed track indexes, not associated with this MC track but containing some
/// hits, produced by this MC track
const auto& GetTouchTrackIndexes() const { return fvTouchTrackIndexes; }
diff --git a/reco/L1/qa/CbmCaTrackFitQa.cxx b/reco/L1/qa/CbmCaTrackFitQa.cxx
index 97875b5c4b..fd8cbbe0ea 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.cxx
+++ b/reco/L1/qa/CbmCaTrackFitQa.cxx
@@ -149,17 +149,24 @@ void TrackFitQa::Fill(const TrackParamV& trPar, const tools::MCPoint& mcPoint, b
const TrackParamV& trParExtr = fitter.Tr(); // Track parameters extrapolated to given MC point
// ** Time-independent measurements **
- FillResAndPull(ETrackParType::kX, trParExtr.GetX()[0], trParExtr.GetXError()[0], mcPoint.GetXOut(), weight);
- FillResAndPull(ETrackParType::kY, trParExtr.GetY()[0], trParExtr.GetYError()[0], mcPoint.GetYOut(), weight);
- FillResAndPull(ETrackParType::kTX, trParExtr.GetTx()[0], trParExtr.GetTxError()[0], mcPoint.GetTxOut(), weight);
- FillResAndPull(ETrackParType::kTY, trParExtr.GetTy()[0], trParExtr.GetTyError()[0], mcPoint.GetTyOut(), weight);
- FillResAndPull(ETrackParType::kQP, trParExtr.GetQp()[0], trParExtr.GetQpError()[0], mcPoint.GetQpOut(), weight);
+ FillResAndPull(ETrackParType::kX, trPar.GetX()[0], trParExtr.GetXError()[0], mcPoint.GetXOut());
+ FillResAndPull(ETrackParType::kY, trPar.GetY()[0], trParExtr.GetYError()[0], mcPoint.GetYOut());
+ FillResAndPull(ETrackParType::kTX, trPar.GetTx()[0], trParExtr.GetTxError()[0], mcPoint.GetTxOut());
+ FillResAndPull(ETrackParType::kTY, trPar.GetTy()[0], trParExtr.GetTyError()[0], mcPoint.GetTyOut());
+ FillResAndPull(ETrackParType::kQP, trPar.GetQp()[0], trParExtr.GetQpError()[0], mcPoint.GetQpOut());
- // TODO: in which point do we calculate MC parameters (central, entrance, exit)??
+ // Momentum resolution
double recoP = std::fabs(mcPoint.GetCharge() / trParExtr.GetQp()[0]); // reco mom. (with MC charge)
double resP = recoP - mcPoint.GetPOut(); // residual of total momentum
- double resPhi = trParExtr.GetPhi()[0] - mcPoint.GetPhiOut(); // residual of azimuthal angle
- double resTheta = trParExtr.GetTheta()[0] - mcPoint.GetThetaOut(); // residual of polar angle
+
+ // Phi resolution
+ double mcPhi = mcPoint.GetPhiOut();
+ double recoTx = trPar.Tx()[0];
+ double recoTy = trPar.Ty()[0];
+ double resPhi = atan2(recoTx * cos(mcPhi) + recoTy * sin(mcPhi), -recoTx * cos(mcPhi) + recoTy * sin(mcPhi));
+
+ // Theta resolution
+ double resTheta = trPar.GetTheta()[0] - mcPoint.GetThetaOut(); // residual of polar angle
resPhi = std::atan2(std::sin(resPhi), std::cos(resPhi)); // overflow over (-pi, pi] protection
@@ -169,10 +176,8 @@ void TrackFitQa::Fill(const TrackParamV& trPar, const tools::MCPoint& mcPoint, b
// ** Time-dependent measurements **
if (bTimeMeasured) {
- FillResAndPull(ETrackParType::kTIME, trParExtr.GetTime()[0], trParExtr.GetTimeError()[0], mcPoint.GetTime(),
- weight);
- FillResAndPull(ETrackParType::kVI, trParExtr.GetVi()[0], trParExtr.GetViError()[0], mcPoint.GetInvSpeedOut(),
- weight);
+ FillResAndPull(ETrackParType::kTIME, trPar.GetTime()[0], trPar.GetTimeError()[0], mcPoint.GetTime());
+ FillResAndPull(ETrackParType::kVI, trPar.GetVi()[0], trPar.GetViError()[0], mcPoint.GetInvSpeedOut());
}
}
diff --git a/reco/L1/qa/CbmCaTrackFitQa.h b/reco/L1/qa/CbmCaTrackFitQa.h
index 1491ee5537..270bbf8171 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.h
+++ b/reco/L1/qa/CbmCaTrackFitQa.h
@@ -181,8 +181,7 @@ namespace cbm::ca
/// @param recoVal Reconstructed error of quantity
/// @param recoErr Error of quantity
/// @param trueVal True value of quantity
- /// @param weight Weight
- void FillResAndPull(ETrackParType type, double recoVal, double recoErr, double trueVal, double w);
+ void FillResAndPull(ETrackParType type, double recoVal, double recoErr, double trueVal);
TrackParArray_t<TH1F*> fvphResiduals = {{0}}; ///< Residuals for different track parameters
TrackParArray_t<TH1F*> fvphPulls = {{0}}; ///< Pulls for different track parameters
@@ -209,13 +208,13 @@ namespace cbm::ca
// ---------------------------------------------------------------------------------------------------------------------
// TODO: Test this function for performance penalties
- inline void TrackFitQa::FillResAndPull(ETrackParType type, double recoVal, double recoErr, double trueVal, double w)
+ inline void TrackFitQa::FillResAndPull(ETrackParType type, double recoVal, double recoErr, double trueVal)
{
if (fvbParIgnored[type]) { return; }
double res = recoVal - trueVal;
double pull = res / recoErr;
- fvphResiduals[type]->Fill(res, w);
- fvphPulls[type]->Fill(pull, w);
+ fvphResiduals[type]->Fill(res);
+ fvphPulls[type]->Fill(pull);
}
} // namespace cbm::ca
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.cxx b/reco/L1/qa/CbmCaTrackTypeQa.cxx
index de706d1e42..f5c9c09be1 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.cxx
+++ b/reco/L1/qa/CbmCaTrackTypeQa.cxx
@@ -14,6 +14,8 @@
#include "CaToolsMCData.h"
+using cbm::algo::ca::FieldRegion;
+using cbm::algo::ca::TrackFit;
using cbm::ca::TrackTypeQa;
using cbm::ca::tools::MCPoint;
using cbm::ca::tools::MCTrack;
@@ -172,6 +174,11 @@ void TrackTypeQa::Init()
fpFitQaVertex->SetTitle("Vertex");
fpFitQaVertex->Init();
}
+
+ // Track fitter
+ fTrackFit.SetDoFitVelocity(true); // TODO: set flag to the configuration
+ fTrackFit.SetMask(fmask::One());
+ fFieldRegion.SetUseOriginalField();
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -211,21 +218,58 @@ void TrackTypeQa::FillRecoTrack(int iTrkReco, double weight)
bool isTimeFitted = (nTimeMeasurements > 1);
// ** First hit **
- int iHfst = recoTrack.GetFirstHitIndex();
- int iPfst = (*fpvHits)[iHfst].GetBestMcPointId();
- if (iPfst > -1) {
- const auto& mcPoint = fpMCData->GetPoint(iPfst);
- TrackParamV trPar(recoTrack);
- fpFitQaFirstHit->Fill(trPar, mcPoint, isTimeFitted, weight);
+ fTrackFit.SetParticleMass(mcTrack.GetMass());
+ {
+ int iHfst = recoTrack.GetFirstHitIndex();
+ int iPfst = (*fpvHits)[iHfst].GetBestMcPointId();
+ if (iPfst > -1) {
+ const auto& mcPoint = fpMCData->GetPoint(iPfst);
+ TrackParamV trPar(recoTrack);
+ fTrackFit.SetTrack(trPar);
+ fTrackFit.Extrapolate(mcPoint.GetZOut(), fFieldRegion);
+ fpFitQaFirstHit->Fill(fTrackFit.Tr(), mcPoint, isTimeFitted);
+ }
}
// ** Last hit **
- int iHlst = recoTrack.GetLastHitIndex();
- int iPlst = (*fpvHits)[iHlst].GetBestMcPointId();
- if (iPlst > -1) {
- const auto& mcPoint = fpMCData->GetPoint(iPlst);
- TrackParamV trPar(recoTrack.TLast);
- fpFitQaLastHit->Fill(trPar, mcPoint, isTimeFitted, weight);
+ {
+ int iHlst = recoTrack.GetLastHitIndex();
+ int iPlst = (*fpvHits)[iHlst].GetBestMcPointId();
+ if (iPlst > -1) {
+ const auto& mcPoint = fpMCData->GetPoint(iPlst);
+ TrackParamV trPar(recoTrack.TLast);
+ fTrackFit.SetTrack(trPar);
+ fTrackFit.Extrapolate(mcPoint.GetZOut(), fFieldRegion);
+ fpFitQaLastHit->Fill(trPar, mcPoint, isTimeFitted);
+ }
+ }
+
+ // ** Vertex **
+ {
+ // Create an MC point for track vertex
+ MCPoint mcTrkVertex = mcTrack.GetVertexPoint();
+ TrackParamV trPar(recoTrack);
+ fTrackFit.SetTrack(trPar);
+ fTrackFit.Extrapolate(mcTrkVertex.GetZ(), fFieldRegion);
+ // Add material
+ int iStFst = (*fpvHits)[recoTrack.GetFirstHitIndex()].GetStationId();
+ double dZ = mcTrkVertex.GetZ() - fpParameters->GetStation(iStFst).fZ[0];
+ int direction = static_cast<int>(dZ / std::fabs(dZ));
+ for (int iSt = iStFst; (iSt >= 0) && (iSt < fpParameters->GetNstationsActive())
+ && (direction * (mcTrkVertex.GetZ() - fpParameters->GetStation(iSt).fZ[0]) > 0);
+ iSt += direction) {
+ fTrackFit.Extrapolate(fpParameters->GetStation(iSt).fZ, fFieldRegion);
+ auto radLength = fpParameters->GetMaterialThickness(iSt, fTrackFit.Tr().GetX(), fTrackFit.Tr().GetY());
+ fTrackFit.MultipleScattering(radLength);
+ fTrackFit.EnergyLossCorrection(radLength, fvec::One());
+ }
+ fTrackFit.Extrapolate(mcTrkVertex.GetZ(), fFieldRegion);
+ const TrackParamV& trParExtr = fTrackFit.Tr();
+ if (mcTrkVertex.GetZ() == trParExtr.GetZ()[0]) { fpFitQaVertex->Fill(trParExtr, mcTrkVertex, isTimeFitted); }
+ else {
+ LOG(warn) << "iTrkReco = " << iTrkReco << ", mcTrkVertex = " << mcTrkVertex.GetZ()
+ << ", par = " << trParExtr.GetZ()[0];
+ }
}
}
}
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.h b/reco/L1/qa/CbmCaTrackTypeQa.h
index 470cbb5a01..2e1a8de3bf 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.h
+++ b/reco/L1/qa/CbmCaTrackTypeQa.h
@@ -18,7 +18,9 @@
#include <map>
#include <string>
+#include "CaField.h"
#include "CaParameters.h"
+#include "CaTrackFit.h"
// Forward declarations
namespace cbm::ca::tools
@@ -275,6 +277,9 @@ 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::FieldRegion fFieldRegion; ///< Magnetic field
+
int fCounterMC = 0; ///< Counter of MC tracks
int fCounterClones = 0; ///< Counter of clone tracks
int fCounterRecoTotal = 0; ///< Counter of reco tracks (total = reco + ghost + clones)
@@ -291,6 +296,7 @@ namespace cbm::ca
tools::MCData* fpMCData = nullptr; ///< Pointer to MC data object
std::shared_ptr<ca::Parameters> fpParameters = nullptr; ///< Pointer to parameters object
+
// ** Cuts on tracks for a given track class **
/// Cut function on MC tracks
--
GitLab