From 69b3fd077ce1676698ffae5896af2baeffe49dc7 Mon Sep 17 00:00:00 2001
From: Alexandru Bercuci <abercuci@niham.nipne.ro>
Date: Thu, 5 May 2022 15:42:26 +0200
Subject: [PATCH] store digits on the QA structure to complete the history of
reconstruction
---
reco/detectors/trd/qa/data/CbmTrdHitMC.cxx | 77 +++++++++++++++++++---
reco/detectors/trd/qa/data/CbmTrdHitMC.h | 57 +++++++++++++---
2 files changed, 114 insertions(+), 20 deletions(-)
diff --git a/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx b/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
index f93565e2ad..db3f1ea430 100644
--- a/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
+++ b/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
@@ -4,12 +4,16 @@
#include "CbmTrdHitMC.h"
+#include "CbmTrdDigi.h"
+
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TVector3.h>
#include <sstream>
+// Error parametrization
+int CbmTrdHitMC::fSx[5][2] = {{190, 2321}, {147, 1920}, {50, 1461}, {22, 2094}, {21, 4297}};
//_____________________________________________________________________
CbmTrdHitMC::CbmTrdHitMC() : CbmTrdHit() {}
@@ -20,10 +24,7 @@ CbmTrdHitMC::CbmTrdHitMC(const CbmTrdHit& h) : CbmTrdHit(h) {}
CbmTrdHitMC::~CbmTrdHitMC() {}
//_____________________________________________________________________
-void CbmTrdHitMC::AddCluster(const CbmTrdCluster* c)
-{
- fCluster = *c;
-}
+void CbmTrdHitMC::AddCluster(const CbmTrdCluster* c) { fCluster = *c; }
//_____________________________________________________________________
size_t CbmTrdHitMC::AddPoint(const CbmTrdPoint* p, double t, int id)
@@ -33,10 +34,27 @@ size_t CbmTrdHitMC::AddPoint(const CbmTrdPoint* p, double t, int id)
}
//_____________________________________________________________________
-size_t CbmTrdHitMC::AddSignal(double s, int t)
+size_t CbmTrdHitMC::AddSignal(const CbmTrdDigi* digi, uint64_t t0)
+{
+ if (GetClassType() == 1) { // TRD2D type
+ int dt;
+ double t, r = digi->GetCharge(t, dt);
+ fTrdSignals.push_back(std::make_pair(t, digi->GetTimeDAQ() - t0));
+ fTrdSignals.push_back(std::make_pair(r, digi->GetTimeDAQ() + dt - t0));
+ }
+ else // TRD1D type
+ fTrdSignals.push_back(std::make_pair(digi->GetCharge(), digi->GetTime() - t0));
+
+ return fTrdSignals.size();
+}
+
+//_____________________________________________________________________
+size_t CbmTrdHitMC::PurgeSignals()
{
- printf("hit%dD : s(%f) t(%d)\n", (GetClassType() ? 2 : 1), s, t);
- fTrdSignals.push_back(std::make_pair(s, t));
+ if (!GetNSignals()) return 0;
+ if (fTrdSignals.front().first < 1.e-3) fTrdSignals.erase(fTrdSignals.begin());
+ if (fTrdSignals.back().first < 1.e-3) fTrdSignals.pop_back();
+
return fTrdSignals.size();
}
@@ -47,6 +65,28 @@ const CbmTrdPoint* CbmTrdHitMC::GetPoint(uint idx) const
return &std::get<0>(fTrdPoints[idx]);
}
+//_____________________________________________________________________
+double CbmTrdHitMC::GetSignal(uint idx) const
+{
+ if (idx >= fTrdSignals.size()) return 0;
+ return fTrdSignals[idx].first;
+}
+
+//_____________________________________________________________________
+CbmTrdHitMC::eCbmTrdHitMCshape CbmTrdHitMC::GetClShape() const
+{
+ if (fCluster.HasOpenStart()) {
+ if (fCluster.HasOpenStop()) return eCbmTrdHitMCshape::kRT;
+ else
+ return eCbmTrdHitMCshape::kRR;
+ }
+ else {
+ if (fCluster.HasOpenStop()) return eCbmTrdHitMCshape::kTT;
+ else
+ return eCbmTrdHitMCshape::kTR;
+ }
+}
+
//_____________________________________________________________________
double CbmTrdHitMC::GetDx() const
{
@@ -56,6 +96,20 @@ double CbmTrdHitMC::GetDx() const
return GetX() - x;
}
+//_____________________________________________________________________
+double CbmTrdHitMC::GetSx() const
+{
+ const CbmTrdPoint* p(nullptr);
+ if (!(p = GetPoint())) return 1;
+ double phi(p->GetPx() / p->GetPz());
+ int isz = GetNSignals() - 2;
+
+ if (isz < 0) isz = 0;
+ if (isz >= 5) isz = 4;
+
+ return 1.e-4 * (fSx[isz][0] + phi * phi * fSx[isz][1]); // error in cm
+}
+
//_____________________________________________________________________
double CbmTrdHitMC::GetDy() const
{
@@ -65,6 +119,9 @@ double CbmTrdHitMC::GetDy() const
return GetY() - y;
}
+//_____________________________________________________________________
+double CbmTrdHitMC::GetSy() const { return GetDy(); }
+
//_____________________________________________________________________
double CbmTrdHitMC::GetDt() const
{
@@ -77,7 +134,7 @@ double CbmTrdHitMC::GetDt() const
TParticlePDG* pmc = (TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg);
if (pdg < 9999999 && pmc) mass = pmc->Mass();
-
+
TVector3 mom3;
p->Momentum(mom3);
t += dz / (p->GetPz() * speedOfLight) * sqrt(mass * mass + mom3.Mag2());
@@ -86,10 +143,10 @@ double CbmTrdHitMC::GetDt() const
//_____________________________________________________________________
std::string CbmTrdHitMC::ToString() const
-{
+{
std::stringstream ss;
for (auto mcp : fTrdPoints) {
- ss << "Event time(ns)=" << std::get<1>(mcp) << " partId=" << std::get<2>(mcp) << "\n";
+ ss << "Event time(ns)=" << std::get<1>(mcp) << " partId=" << std::get<2>(mcp) << "\n";
ss << std::get<0>(mcp).ToString();
}
diff --git a/reco/detectors/trd/qa/data/CbmTrdHitMC.h b/reco/detectors/trd/qa/data/CbmTrdHitMC.h
index a2cff1f130..1be20cd029 100644
--- a/reco/detectors/trd/qa/data/CbmTrdHitMC.h
+++ b/reco/detectors/trd/qa/data/CbmTrdHitMC.h
@@ -5,12 +5,12 @@
#ifndef CBMTRDHITMC_H
#define CBMTRDHITMC_H
+#include "CbmTrdCluster.h"
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"
-#include "CbmTrdCluster.h"
-#include <vector> // for fTrdPoints
#include <string> // for ToString
+#include <vector> // for fTrdPoints
/** \class CbmTrdHitMC
* \brief TRD hit to MC point correlation class
@@ -24,8 +24,17 @@
* To describe main functionality ...
*/
+class CbmTrdDigi;
class CbmTrdHitMC : public CbmTrdHit {
public:
+ enum eCbmTrdHitMCshape
+ {
+ kRT = 0, // open left - open right shape
+ kRR, // open left - closed right shape
+ kTT, // closed left - open right shape
+ kTR // closed left - closed right shape
+ };
+
/** \brief Default constructor.*/
CbmTrdHitMC();
@@ -46,29 +55,55 @@ public:
*/
size_t AddPoint(const CbmTrdPoint* p, double t, int id);
/** \brief Add signal values in the increasing order of pad index
- * \param s signal from ch/pad
- * \param t relative time in the cluster
+ * \param d digi from ch/pad
+ * \param t0 relative time in the cluster
* \return the number of signals in the cluster
*/
- size_t AddSignal(double s, int t);
+ size_t AddSignal(const CbmTrdDigi* d, uint64_t t0);
+
+ /** \brief return MC pile-up size*/
+ std::string GetErrorMsg() const { return fErrMsg; }
/** \brief Register a MC point
* \param idx index of point being requested. by default the best fit is returned.
*/
const CbmTrdPoint* GetPoint(uint idx = 0) const;
+ /** \brief return signal at position
+ * \param idx index of signal counting from left to right looking upstream (from FEE).
+ */
+ double GetSignal(uint idx = 0) const;
+
+ /** \brief return MC pile-up size*/
+ size_t GetNPoints() const { return fTrdPoints.size(); }
+
+ /** \brief return cluster size*/
+ size_t GetNSignals() const { return fTrdSignals.size(); }
+
+ /** \brief return cluster shape according to the eCbmTrdHitMCshape definitions*/
+ eCbmTrdHitMCshape GetClShape() const;
+
/** \brief Calculate residuals in the bending plane.*/
double GetDx() const;
+ /** \brief Calculate error in the bending plane.*/
+ double GetSx() const;
+
/** \brief Calculate residuals for the azimuth direction.*/
double GetDy() const;
+ /** \brief Calculate error for the azimuth direction.*/
+ double GetSy() const;
+
/** \brief Calculate residuals for time.*/
double GetDt() const;
/** \brief Store error message.*/
void SetErrorMsg(std::string msg) { fErrMsg = msg; }
-
+
+ /** \brief Applies to TRD2D and remove 0 charges from the boundaries of the cluster.**/
+ size_t PurgeSignals();
+
/** \brief Verbosity functionality.**/
virtual std::string ToString() const;
@@ -78,10 +113,12 @@ private:
/** \brief Assignment operator.*/
CbmTrdHitMC& operator=(const CbmTrdHitMC&) = default;
- std::string fErrMsg = ""; //< error message from the QA task
- std::vector<std::pair<double, int>> fTrdSignals = {}; //< list of signal/time in cluster
- std::vector<std::tuple<CbmTrdPoint, double, int>> fTrdPoints = {}; //< list of MC points together with the event time and particle PDG code producing them
- CbmTrdCluster fCluster; //< data from the cluster
+ std::string fErrMsg = ""; //< error message from the QA task
+ std::vector<std::pair<double, int>> fTrdSignals = {}; //< list of signal/time in cluster
+ std::vector<std::tuple<CbmTrdPoint, double, int>> fTrdPoints =
+ {}; //< list of MC points together with the event time and particle PDG code producing them
+ CbmTrdCluster fCluster; //< data from the cluster
+ static int fSx[5][2]; //< x error parametrization as function of cluster size and incident direction
ClassDef(CbmTrdHitMC, 1) // Hit to MC point data correlation
};
--
GitLab