From 2d5d47a0811a4e08dacd2fe6bcbcda2ad0083343 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Thu, 27 Apr 2023 15:23:24 +0200
Subject: [PATCH] CA: QA: Added resolution profiles for momentum, theta and phi
---
reco/L1/L1Algo/L1TrackPar.h | 63 +++++++++++++++++++++++++++++++-
reco/L1/catools/CaToolsMCPoint.h | 18 +++++++++
reco/L1/qa/CbmCaTrackFitQa.cxx | 21 +++++++++++
reco/L1/qa/CbmCaTrackFitQa.h | 5 +++
4 files changed, 105 insertions(+), 2 deletions(-)
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index a8a118b67a..af9eca6e0e 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -109,6 +109,18 @@ public:
/// @brief Gets inverse speed error [ns/cm]
fscal GetInvSpeedErr() const { return (std::isfinite(C66[0]) && C66[0] > 0) ? std::sqrt(C66[0]) : undef::kF; }
+ /// @brief Gets azimuthal angle [rad]
+ fscal GetPhi() const { return std::atan2(-GetTy(), -GetTx()); }
+
+ /// @brief Gets azimuthal angle error [rad]
+ fscal GetPhiErr() const;
+
+ /// @brief Gets polar angle [rad]
+ fscal GetTheta() const { return std::acos(1. / sqrt(GetTx() * GetTx() + GetTy() * GetTy() + 1.)); }
+
+ /// @brief Gets polar angle error [rad]
+ fscal GetThetaErr() const;
+
/// @brief Resets variances of track parameters
/// @param c00 Variance of x-position [cm2]
/// @param c11 Variance of y-position [cm2]
@@ -149,8 +161,46 @@ public:
} _fvecalignment;
-// =============================================================================================
+// ***************************************************
+// ** Inline and template function implementation **
+// ***************************************************
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+inline fscal L1TrackPar::GetPhiErr() const
+{
+ fscal phiDFactor = 1. / (GetTx() * GetTx() + GetTy() * GetTy());
+ fscal phiDTx = -phiDFactor * GetTy(); // partial derivative of phi over Tx
+ fscal phiDTy = +phiDFactor * GetTx(); // partial derivative of phi over Ty
+
+ fscal varTx = (std::isfinite(C22[0]) && C22[0] > 0) ? C22[0] : undef::kF; // variance of Tx
+ fscal varTy = (std::isfinite(C33[0]) && C33[0] > 0) ? C33[0] : undef::kF; // variance of Ty
+ fscal covTxTy = std::isfinite(C32[0]) ? C32[0] : undef::kF; // covariance of Tx and Ty
+
+ fscal varPhi = phiDTx * phiDTx * varTx + phiDTy * phiDTy * varTy + 2 * phiDTx * phiDTy * covTxTy;
+ return std::sqrt(varPhi);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+inline fscal L1TrackPar::GetThetaErr() const
+{
+ fscal sumSqSlopes = GetTx() * GetTx() + GetTy() * GetTy();
+ fscal thetaDFactor = 1. / ((sumSqSlopes + 1) * sqrt(sumSqSlopes));
+ fscal thetaDTx = thetaDFactor * GetTx();
+ fscal thetaDTy = thetaDFactor * GetTy();
+
+ fscal varTx = (std::isfinite(C22[0]) && C22[0] > 0) ? C22[0] : undef::kF; // variance of Tx
+ fscal varTy = (std::isfinite(C33[0]) && C33[0] > 0) ? C33[0] : undef::kF; // variance of Ty
+ fscal covTxTy = std::isfinite(C32[0]) ? C32[0] : undef::kF; // covariance of Tx and Ty
+ fscal varTheta = thetaDTx * thetaDTx * varTx + thetaDTy * thetaDTy * varTy + 2 * thetaDTx * thetaDTy * covTxTy;
+ return std::sqrt(varTheta);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
template<typename T>
inline void L1TrackPar::copyToArrays(const int iVec, T p[kNparTr], T c[kNparCov]) const
{
@@ -169,6 +219,8 @@ inline void L1TrackPar::copyToArrays(const int iVec, T p[kNparTr], T c[kNparCov]
}
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
template<typename T>
inline void L1TrackPar::copyFromArrays(const int iVec, const T p[kNparTr], const T c[kNparCov])
{
@@ -190,6 +242,8 @@ inline void L1TrackPar::copyFromArrays(const int iVec, const T p[kNparTr], const
NDF = 0.;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
template<typename T>
inline void L1TrackPar::copyFromArrays(const T p[kNparTr], const T c[kNparCov])
{
@@ -211,7 +265,8 @@ inline void L1TrackPar::copyFromArrays(const T p[kNparTr], const T c[kNparCov])
NDF = 0.;
}
-
+// ---------------------------------------------------------------------------------------------------------------------
+//
inline void L1TrackPar::SetOneEntry(const int i0, const L1TrackPar& T1, const int i1)
{
x[i0] = T1.x[i1];
@@ -256,6 +311,8 @@ inline void L1TrackPar::SetOneEntry(const int i0, const L1TrackPar& T1, const in
NDF[i0] = T1.NDF[i1];
} // SetOneEntry
+// ---------------------------------------------------------------------------------------------------------------------
+//
inline void L1TrackPar::ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec c44, fvec c55, fvec c66)
{
C10 = 0.;
@@ -278,6 +335,8 @@ inline void L1TrackPar::ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec
nTimeMeasurements = 0.;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
inline void L1TrackPar::InitVelocityRange(fscal minP)
{
// initialise the velocity range with respect to the minimal momentum minP {GeV/c}
diff --git a/reco/L1/catools/CaToolsMCPoint.h b/reco/L1/catools/CaToolsMCPoint.h
index d1476138d4..158f35e6e3 100644
--- a/reco/L1/catools/CaToolsMCPoint.h
+++ b/reco/L1/catools/CaToolsMCPoint.h
@@ -104,6 +104,15 @@ namespace ca::tools
/// @brief Gets PDG code of the particle
int GetPdgCode() const { return fPdgCode; }
+ /// @brief Gets track azimuthal angle at reference z of station [rad]
+ double GetPhi() const { return std::atan2(-fMom[1], fMom[0]); }
+
+ /// @brief Gets track azimuthal angle at entrance to station [rad]
+ double GetPhiIn() const { return std::atan2(-fMomIn[1], fMomIn[0]); }
+
+ /// @brief Gets track azimuthal angle at exit of station [rad]
+ double GetPhiOut() const { return std::atan2(-fMomOut[1], fMomOut[0]); }
+
/// @brief Gets track momentum absolute value at entrance to station [GeV/c]
double GetPIn() const { return std::sqrt(fMomIn[0] * fMomIn[0] + fMomIn[1] * fMomIn[1] + fMomIn[2] * fMomIn[2]); }
@@ -152,6 +161,15 @@ namespace ca::tools
/// @brief Gets global ID of the active tracking station
int GetStationId() const { return fStationId; }
+ /// @brief Gets polar angle at reference z of station [rad]
+ double GetTheta() const { return std::acos(GetPz() / GetP()); }
+
+ /// @brief Gets polar angle at entrance to station [rad]
+ double GetThetaIn() const { return std::acos(GetPzIn() / GetPIn()); }
+
+ /// @brief Gets polar angle at exit of station [rad]
+ double GetThetaOut() const { return std::acos(GetPzOut() / GetPOut()); }
+
/// @brief Gets time [ns]
double GetTime() const { return fTime; }
diff --git a/reco/L1/qa/CbmCaTrackFitQa.cxx b/reco/L1/qa/CbmCaTrackFitQa.cxx
index 7769050a36..21b23c0c4c 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.cxx
+++ b/reco/L1/qa/CbmCaTrackFitQa.cxx
@@ -14,6 +14,8 @@
#include "CaToolsMCData.h"
#include "L1Field.h"
#include "L1Fit.h"
+#include "TProfile.h"
+#include "TH1.h"
using cbm::ca::TrackFitQa;
@@ -37,6 +39,11 @@ void TrackFitQa::Init()
fph_res_qp = MakeHisto<TH1F>("res_qp", "", fBinsRESQP, fLoRESQP, fUpRESQP);
fph_res_vi = MakeHisto<TH1F>("res_vi", "", fBinsRESVI, fLoRESVI, fUpRESVI);
+ // FIXME: Replace hardcoded parameters with variables
+ fph_res_p_vs_pMC = MakeHisto<TProfile>("res_p_vs_pMC", "", 100, 0.0, 10.0, -2., 2.);
+ fph_res_phi_vs_phiMC = MakeHisto<TProfile>("res_phi_vs_phiMC", "", 100, -3.2, 3.2, -2., 2.);
+ fph_res_theta_vs_thetaMC = MakeHisto<TProfile>("res_theta_vs_phiMC", "", 100, 0., 3.2, -2., 2.);
+
fph_pull_x = MakeHisto<TH1F>("pull_x", "", fBinsPULLX, fLoPULLX, fUpPULLX);
fph_pull_y = MakeHisto<TH1F>("pull_y", "", fBinsPULLY, fLoPULLY, fUpPULLY);
fph_pull_t = MakeHisto<TH1F>("pull_t", "", fBinsPULLT, fLoPULLT, fUpPULLT);
@@ -55,6 +62,10 @@ void TrackFitQa::Init()
fph_res_qp->SetTitle(sPrefix + " residual for q/p;(q/p)_{reco} - (q/p)_{MC} [ec/GeV]");
fph_res_vi->SetTitle(sPrefix + " residual for inverse speed;v^{-1}_{reco} - v^{-1}_{MC} [c^{-1}]");
+ fph_res_p_vs_pMC->SetTitle(sPrefix + " resolution of momentum;p^{MC} [GeV/c];#delta p [GeV/c]");
+ fph_res_phi_vs_phiMC->SetTitle(sPrefix + " resolution of polar angle;#phi^{MC} [rad];#delta#phi [rad]");
+ fph_res_theta_vs_thetaMC->SetTitle(sPrefix + " resolution of polar angle;#theta^{MC} [rad];#delta#theta [rad]");
+
fph_pull_x->SetTitle(sPrefix + " pull for x-coordinate;(x_{reco} - x_{MC}) / #sigma_{x}");
fph_pull_y->SetTitle(sPrefix + " pull for y-coordinate;(y_{reco} - y_{MC}) / #sigma_{y}");
fph_pull_t->SetTitle(sPrefix + " pull for time;(t_{reco} - t_{MC}) / #sigma_{t}");
@@ -87,6 +98,12 @@ void TrackFitQa::Fill(const L1TrackPar& trPar, const ::ca::tools::MCPoint& mcPoi
double resTx = trParExtr.GetTx() - mcPoint.GetTxOut(); // residual of slope along x-axis [1]
double resTy = trParExtr.GetTy() - mcPoint.GetTyOut(); // residual of slope along y-axis [1]
double resQp = trParExtr.GetQp() - mcPoint.GetQpOut(); // residual of q/p [ec/GeV]
+
+ // TODO: in which point do we calculate MC parameters (center, in, out)??
+ double recoP = std::fabs(mcPoint.GetCharge() / trParExtr.GetQp()); // reco mom. (with MC charge)
+ double resP = recoP - mcPoint.GetPOut(); // residual of total momentum
+ double resPhi = trParExtr.GetPhi() - mcPoint.GetPhiOut(); // residual of azimuthal angle
+ double resTheta = trParExtr.GetTheta() - mcPoint.GetThetaOut(); // residual of polar angle
double pullX = resX / trParExtr.GetXErr(); // pull of x-position
double pullY = resY / trParExtr.GetYErr(); // pull of y-position
@@ -100,6 +117,10 @@ void TrackFitQa::Fill(const L1TrackPar& trPar, const ::ca::tools::MCPoint& mcPoi
fph_res_ty->Fill(resTy, weight);
fph_res_qp->Fill(resQp, weight);
+ fph_res_p_vs_pMC->Fill(mcPoint.GetPOut(), resP);
+ fph_res_phi_vs_phiMC->Fill(mcPoint.GetPhiOut(), resPhi);
+ fph_res_theta_vs_thetaMC->Fill(mcPoint.GetThetaOut(), resTheta);
+
fph_pull_x->Fill(pullX, weight);
fph_pull_y->Fill(pullY, weight);
fph_pull_tx->Fill(pullTx, weight);
diff --git a/reco/L1/qa/CbmCaTrackFitQa.h b/reco/L1/qa/CbmCaTrackFitQa.h
index 7211c325f8..837c7e3eee 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.h
+++ b/reco/L1/qa/CbmCaTrackFitQa.h
@@ -28,6 +28,7 @@ class CbmL1HitDebugInfo;
class TFolder;
class TH1F;
+class TProfile;
namespace cbm::ca
{
@@ -100,6 +101,10 @@ namespace cbm::ca
TH1F* fph_pull_t = nullptr; ///< Pull of time
TH1F* fph_pull_vi = nullptr; ///< Pull of inverse speed
+ // ** Resolution profiles **
+ TProfile* fph_res_p_vs_pMC = nullptr; ///< Resolution of momentum [GeV/c]
+ TProfile* fph_res_theta_vs_thetaMC = nullptr; ///< Resolution of polar angle [rad]
+ TProfile* fph_res_phi_vs_phiMC = nullptr; ///< Resolution of azimuthal angle [rad]
// **************************
// ** Histogram properties **
--
GitLab