From 4c853da934b82531c90005da15c4bac9536eadcb Mon Sep 17 00:00:00 2001
From: Alexandru Bercuci <abercuci@niham.nipne.ro>
Date: Tue, 3 May 2022 15:36:10 +0200
Subject: [PATCH] extend the TrdHit Qa data class for digi signals and clusters
---
reco/L1/qa/CbmTrackerInputQaTrd.cxx | 123 ++++++++++++++++-----
reco/detectors/trd/qa/data/CbmTrdHitMC.cxx | 52 ++++++++-
reco/detectors/trd/qa/data/CbmTrdHitMC.h | 38 ++++++-
3 files changed, 177 insertions(+), 36 deletions(-)
diff --git a/reco/L1/qa/CbmTrackerInputQaTrd.cxx b/reco/L1/qa/CbmTrackerInputQaTrd.cxx
index 8a38cc6143..ecc587141e 100644
--- a/reco/L1/qa/CbmTrackerInputQaTrd.cxx
+++ b/reco/L1/qa/CbmTrackerInputQaTrd.cxx
@@ -497,13 +497,17 @@ void CbmTrackerInputQaTrd::ResolutionQa()
return;
}
+ // construct QA object
+ CbmTrdHitMC* hmc = new ((*fHitsMc)[imc++]) CbmTrdHitMC(*hit);
+ hmc->AddCluster(cluster);
+ uint64_t tdigi = 0;
+
// custom finder of the digi matches
-
CbmMatch myHitMatch;
for (Int_t iDigi = 0; iDigi < nClusterDigis; iDigi++) {
Int_t digiIdx = cluster->GetDigi(iDigi);
if (digiIdx < 0 || digiIdx >= nDigis) {
- LOG(error) << "TRD cluster: digi index " << digiIdx << " out of range ";
+ LOG(fatal) << "TRD cluster: digi index " << digiIdx << " out of range ";
return;
}
const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(digiIdx);
@@ -513,9 +517,31 @@ void CbmTrackerInputQaTrd::ResolutionQa()
}
if (digi->GetAddressModule() != address) {
- LOG(error) << "TRD hit address " << address << " differs from its digi address " << digi->GetAddressModule();
+ std::stringstream ss;
+ ss << "TRD hit address " << address << " differs from its digi address " << digi->GetAddressModule();
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
+ switch (digi->GetType()) {
+ case CbmTrdDigi::eCbmTrdAsicType::kFASP:
+ {
+ int dt;
+ double t, r = digi->GetCharge(t, dt);
+ if (!tdigi) tdigi = digi->GetTimeDAQ();
+ if (t > 0) hmc->AddSignal(t, digi->GetTimeDAQ() - tdigi);
+ if (r > 0) hmc->AddSignal(r, digi->GetTimeDAQ() - tdigi + dt);
+ break;
+ }
+ case CbmTrdDigi::eCbmTrdAsicType::kSPADIC:
+ hmc->AddSignal(digi->GetCharge(), digi->GetTime());
+ break;
+ default:
+ LOG(fatal) << "TRD digi type neither SPADIC or FASP";
+ return;
+ }
+
+
const CbmMatch* match = dynamic_cast<const CbmMatch*>(fDigiManager->GetMatch(ECbmModuleId::kTrd, digiIdx));
if (!match) {
LOG(fatal) << "TRD digi match " << digiIdx << " not found";
@@ -529,11 +555,19 @@ void CbmTrackerInputQaTrd::ResolutionQa()
{ // check if the hit match is correct
CbmMatch* hitMatch = dynamic_cast<CbmMatch*>(fHitMatches->At(iHit));
- if (hitMatch == nullptr) { LOG(error) << "hit match for TRD hit " << iHit << " doesn't exist"; }
+ if (hitMatch == nullptr) {
+ std::stringstream ss;
+ ss << "hit match for TRD hit " << iHit << " doesn't exist";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
+ }
else {
const CbmLink& link = hitMatch->GetMatchedLink();
if ((link != bestLink) && (link.GetWeight() != bestLink.GetWeight())) {
- LOG(error) << "hit match for TRD hit " << iHit << " doesn't correspond to digi matches";
+ std::stringstream ss;
+ ss << "hit match for TRD hit " << iHit << " doesn't correspond to digi matches";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
}
}
}
@@ -547,11 +581,17 @@ void CbmTrackerInputQaTrd::ResolutionQa()
nHitPoints++;
int iE = fMcEventList->GetEventIndex(link);
if (iE < 0 || iE >= nMcEvents) {
- LOG(error) << "link points to a non-existing MC event";
+ std::stringstream ss;
+ ss << "link points to a non-existing MC event";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
if (link.GetIndex() >= (int) pointNhits[iE].size()) {
- LOG(error) << "link points to a non-existing MC point";
+ std::stringstream ss;
+ ss << "link points to a non-existing MC index";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
pointNhits[iE][link.GetIndex()]++;
@@ -563,19 +603,27 @@ void CbmTrackerInputQaTrd::ResolutionQa()
// take corresponding MC point
// skip hits from the noise digis
- if (bestLink.GetIndex() < 0) { continue; }
-
- // construct QA object
- CbmTrdHitMC* hmc = new ((*fHitsMc)[imc++]) CbmTrdHitMC(*hit);
-
+ if (bestLink.GetIndex() < 0) {
+ std::stringstream ss;
+ ss << "hit from noise [INFO]";
+ hmc->SetErrorMsg(ss.str());
+ continue;
+ }
+
CbmTrdPoint* p = dynamic_cast<CbmTrdPoint*>(fMcPoints->Get(bestLink));
if (p == nullptr) {
- LOG(error) << "link points to a non-existing MC point";
+ std::stringstream ss;
+ ss << "link points to a non-existing MC point";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
if (p->GetModuleAddress() != (int) CbmTrdAddress::GetModuleAddress(address)) {
- LOG(error) << "mc point module address differs from the hit module address";
+ std::stringstream ss;
+ ss << "mc point module address differs from the hit module address";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
@@ -593,25 +641,39 @@ void CbmTrackerInputQaTrd::ResolutionQa()
double staZ = pGeo->GetZ(); // module->GetZ(); //+ 410;
if ((staZ < p->GetZIn() - 1.) || (staZ > p->GetZOut() + 1.)) {
- LOG(error) << "TRD station " << StationIndex << ": active material Z[" << p->GetZIn() << " cm," << p->GetZOut()
+ std::stringstream ss;
+ ss << "TRD station " << StationIndex << ": active material Z[" << p->GetZIn() << " cm," << p->GetZOut()
<< " cm] is too far from the nominal station Z " << staZ << " cm";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
// the cut of 1 cm is choosen arbitrary and can be changed
if (fabs(hit->GetZ() - staZ) > 1.) {
- LOG(error) << "TRD station " << StationIndex << ": hit Z " << hit->GetZ()
+ std::stringstream ss;
+ ss << "TRD station " << StationIndex << ": hit Z " << hit->GetZ()
<< " cm, is too far from the nominal station Z " << staZ << " cm";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
}
// residual and pull
- if (nHitPoints != 1) continue; // only check residual for non-mixed hits
+ if (nHitPoints != 1) {
+ std::stringstream ss;
+ ss << "hit from mixed hit [INFO] nPoints=" << nHitPoints;
+ hmc->SetErrorMsg(ss.str());
+ continue; // only check residual for non-mixed hits
+ }
double t0 = fMcEventList->GetEventTime(bestLink);
if (t0 < 0) {
- LOG(error) << "MC event time doesn't exist for a TRD link";
+ std::stringstream ss;
+ ss << "MC event time doesn't exist for a TRD link";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
@@ -632,30 +694,35 @@ void CbmTrackerInputQaTrd::ResolutionQa()
y += dz * py / pz;
// get particle mass
+ Int_t pdg(0);
double mass = 0;
{
CbmLink mcTrackLink = bestLink;
mcTrackLink.SetIndex(p->GetTrackID());
CbmMCTrack* mcTrack = dynamic_cast<CbmMCTrack*>(fMcTracks->Get(mcTrackLink));
if (!mcTrack) {
- LOG(error) << "MC track " << p->GetTrackID() << " doesn't exist";
+ std::stringstream ss;
+ ss << "MC track " << p->GetTrackID() << " doesn't exist";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
- Int_t pdg = mcTrack->GetPdgCode();
+ pdg = mcTrack->GetPdgCode();
if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
}
}
- hmc->AddPoint(p, t0, mass);
-
+ hmc->AddPoint(p, t0, pdg);
+ //std::cout << hmc->ToString() << "\n";
+
constexpr double speedOfLight = 29.979246; // cm/ns
TVector3 mom3;
p->Momentum(mom3);
t += dz / (pz * speedOfLight) * sqrt(mass * mass + mom3.Mag2());
- double du = hmc->GetDx(); //hit->GetX() - x;
- double dv = hmc->GetDy(); //hit->GetY() - y;
- double dt = hmc->GetDt(); //hit->GetTime() - t;
+ double du = hit->GetX() - x; //hmc->GetDx(); //hit->GetX() - x;
+ double dv = hit->GetY() - y; //hmc->GetDy(); //hit->GetY() - y;
+ double dt = hit->GetTime() - t; // hmc->GetDt(); //hit->GetTime() - t;
double su = hit->GetDx();
double sv = hit->GetDy();
double st = hit->GetTimeError();
@@ -678,8 +745,10 @@ void CbmTrackerInputQaTrd::ResolutionQa()
// pulls
if ((su < 1.e-5) || (sv < 1.e-5) || (st < 1.e-5)) {
- LOG(error) << "TRD hit errors are not properly set: errX " << hit->GetDx() << " errY " << hit->GetDy() << " errT "
- << st;
+ std::stringstream ss;
+ ss << "TRD hit errors are not properly set: errX " << hit->GetDx() << " errY " << hit->GetDy() << " errT ";
+ hmc->SetErrorMsg(ss.str());
+ LOG(error) << ss.str();
return;
}
diff --git a/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx b/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
index 48406d78e6..f93565e2ad 100644
--- a/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
+++ b/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
@@ -4,8 +4,12 @@
#include "CbmTrdHitMC.h"
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
#include <TVector3.h>
+#include <sstream>
+
//_____________________________________________________________________
CbmTrdHitMC::CbmTrdHitMC() : CbmTrdHit() {}
@@ -16,12 +20,26 @@ CbmTrdHitMC::CbmTrdHitMC(const CbmTrdHit& h) : CbmTrdHit(h) {}
CbmTrdHitMC::~CbmTrdHitMC() {}
//_____________________________________________________________________
-size_t CbmTrdHitMC::AddPoint(const CbmTrdPoint* p, double t, double m)
+void CbmTrdHitMC::AddCluster(const CbmTrdCluster* c)
+{
+ fCluster = *c;
+}
+
+//_____________________________________________________________________
+size_t CbmTrdHitMC::AddPoint(const CbmTrdPoint* p, double t, int id)
{
- fTrdPoints.push_back(std::make_tuple(CbmTrdPoint(*p), t, m));
+ fTrdPoints.push_back(std::make_tuple(CbmTrdPoint(*p), t, id));
return fTrdPoints.size();
}
+//_____________________________________________________________________
+size_t CbmTrdHitMC::AddSignal(double s, int t)
+{
+ printf("hit%dD : s(%f) t(%d)\n", (GetClassType() ? 2 : 1), s, t);
+ fTrdSignals.push_back(std::make_pair(s, t));
+ return fTrdSignals.size();
+}
+
//_____________________________________________________________________
const CbmTrdPoint* CbmTrdHitMC::GetPoint(uint idx) const
{
@@ -54,13 +72,39 @@ double CbmTrdHitMC::GetDt() const
const CbmTrdPoint* p(nullptr);
if (!(p = GetPoint())) return -999;
- double t0(std::get<1>(fTrdPoints[0])), mass(std::get<2>(fTrdPoints[0]));
- double dz(GetZ() - p->GetZ()), t(t0 + p->GetTime());
+ int pdg = std::get<2>(fTrdPoints[0]);
+ double t0(std::get<1>(fTrdPoints[0])), dz(GetZ() - p->GetZ()), t(t0 + p->GetTime()), mass(0);
+ 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());
return GetTime() - t;
}
+//_____________________________________________________________________
+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 << std::get<0>(mcp).ToString();
+ }
+
+ ss << "CbmTrdDigi: [" << fTrdSignals.size() << "] Signal / Relative Time\n ";
+ for (auto sgn : fTrdSignals) {
+ ss << sgn.first << "/" << sgn.second << " ";
+ }
+ ss << "\n";
+
+ ss << fCluster.ToString();
+
+ ss << CbmTrdHit::ToString();
+
+ if (fErrMsg != "") ss << "Error : " << fErrMsg << "\n";
+ return ss.str();
+}
+
ClassImp(CbmTrdHitMC)
diff --git a/reco/detectors/trd/qa/data/CbmTrdHitMC.h b/reco/detectors/trd/qa/data/CbmTrdHitMC.h
index 015b795d2b..a2cff1f130 100644
--- a/reco/detectors/trd/qa/data/CbmTrdHitMC.h
+++ b/reco/detectors/trd/qa/data/CbmTrdHitMC.h
@@ -7,53 +7,81 @@
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"
+#include "CbmTrdCluster.h"
-#include <vector> // for vector
+#include <vector> // for fTrdPoints
+#include <string> // for ToString
/** \class CbmTrdHitMC
* \brief TRD hit to MC point correlation class
* \author Alexandru Bercuci <abercuci@niham.nipne.ro>
* \date 02.05.2022
*
- * The class
+ * The class packs the whole history of a TrdHit from the set of MC points generating the signals,
+ * to information on the digits structure in the cluster and the reconstructed observables.
+ * The class should be used to make error parametrization, pile-up estimations, etc.
+ *
+ * To describe main functionality ...
*/
class CbmTrdHitMC : public CbmTrdHit {
public:
/** \brief Default constructor.*/
CbmTrdHitMC();
+
/** \brief Standard constructor.*/
CbmTrdHitMC(const CbmTrdHit& hit);
+
/** \brief Destructor. */
virtual ~CbmTrdHitMC();
+ /** \brief Copy cluster details.*/
+ void AddCluster(const CbmTrdCluster* c);
+
/** \brief Add MC points to the hit. The first time this function is called is for the best matched MC point
* \param p: pointer to the point being added
* \param t: time of the event to which the point belongs
* \param m: mass of the particle producing the point
* \return the number of points in the list
*/
- size_t AddPoint(const CbmTrdPoint* p, double t, double m);
+ 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
+ * \return the number of signals in the cluster
+ */
+ size_t AddSignal(double s, int t);
- /** \brief Retrieve a MC point
+ /** \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 Calculate residuals in the bending plane.*/
double GetDx() const;
+
/** \brief Calculate residuals for the azimuth direction.*/
double GetDy() const;
+
/** \brief Calculate residuals for time.*/
double GetDt() const;
+ /** \brief Store error message.*/
+ void SetErrorMsg(std::string msg) { fErrMsg = msg; }
+
+ /** \brief Verbosity functionality.**/
+ virtual std::string ToString() const;
+
private:
/** \brief Copy Constructor.*/
CbmTrdHitMC(const CbmTrdHitMC&) = default;
/** \brief Assignment operator.*/
CbmTrdHitMC& operator=(const CbmTrdHitMC&) = default;
- std::vector<std::tuple<CbmTrdPoint, double, double>> fTrdPoints = {};
+ 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
ClassDef(CbmTrdHitMC, 1) // Hit to MC point data correlation
};
--
GitLab