From 5758cc4c3972924f14ed0b405580e3d61cb12bd9 Mon Sep 17 00:00:00 2001
From: Alexandru Bercuci <abercuci@niham.nipne.ro>
Date: Mon, 2 May 2022 21:52:01 +0200
Subject: [PATCH] add TrdHit QA data class which packs the whole history of a
reconstructed hit from MC points to the actual data. To be used in fine
tuning the error parametrization, pile-ups etc.
---
reco/L1/CMakeLists.txt | 2 +
reco/L1/qa/CbmTrackerInputQaTrd.cxx | 28 +++++++--
reco/L1/qa/CbmTrackerInputQaTrd.h | 9 +--
reco/detectors/trd/CMakeLists.txt | 3 +
reco/detectors/trd/CbmTrdRecoLinkDef.h | 1 +
reco/detectors/trd/qa/data/CbmTrdHitMC.cxx | 66 ++++++++++++++++++++++
reco/detectors/trd/qa/data/CbmTrdHitMC.h | 60 ++++++++++++++++++++
7 files changed, 159 insertions(+), 10 deletions(-)
create mode 100644 reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
create mode 100644 reco/detectors/trd/qa/data/CbmTrdHitMC.h
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index d2602fd95d..5ac6e870c6 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -21,6 +21,7 @@ ${CBMROOT_SOURCE_DIR}/reco/base
${CBMROOT_SOURCE_DIR}/reco/detectors/sts
${CBMROOT_SOURCE_DIR}/reco/detectors/rich
${CBMROOT_SOURCE_DIR}/reco/detectors/rich/fitter
+${CBMROOT_SOURCE_DIR}/reco/detectors/trd/qa/data
${CBMBASE_DIR}
${CBMDATA_DIR}
@@ -240,6 +241,7 @@ Set(DEPENDENCIES
# CbmGeoSetup
CbmMuchBase
CbmTrdBase
+ CbmTrdReco
CbmTofBase
CbmStsBase
CbmRecoBase
diff --git a/reco/L1/qa/CbmTrackerInputQaTrd.cxx b/reco/L1/qa/CbmTrackerInputQaTrd.cxx
index 5ebe1fae51..8a38cc6143 100644
--- a/reco/L1/qa/CbmTrackerInputQaTrd.cxx
+++ b/reco/L1/qa/CbmTrackerInputQaTrd.cxx
@@ -22,6 +22,7 @@
#include "CbmTrdCluster.h"
#include "CbmTrdDigi.h"
#include "CbmTrdHit.h"
+#include "CbmTrdHitMC.h"
#include "CbmTrdParModDigi.h" // for CbmTrdModule
#include "CbmTrdParModGeo.h"
#include "CbmTrdParSetDigi.h" // for CbmTrdParSetDigi
@@ -118,6 +119,12 @@ void CbmTrackerInputQaTrd::DeInit()
fhPointsPerHit.clear();
fhHitsPerPoint.clear();
+
+ if (fHitsMc) {
+ fHitsMc->Clear("C");
+ fHitsMc->Delete();
+ delete fHitsMc;
+ }
}
// -------------------------------------------------------------------------
@@ -270,6 +277,10 @@ InitStatus CbmTrackerInputQaTrd::ReInit()
return kERROR;
}
+ // init output tree
+ fHitsMc = new TClonesArray("CbmTrdHitMC", 100);
+ manager->Register("TrdHitMC", "TRD", fHitsMc, IsOutputBranchPersistent("TrdHitMC"));
+
// initialise histograms
fOutFolder.SetOwner(false);
fHistFolder = fOutFolder.AddFolder("rawHist", "Raw Histograms");
@@ -425,9 +436,10 @@ void CbmTrackerInputQaTrd::ResolutionQa()
Int_t nDigis = fDigiManager->GetNofDigis(ECbmModuleId::kTrd);
int nMcEvents = fMcEventList->GetNofEvents();
+ int imc(0); // index of hit->MC QA objects
+ fHitsMc->Delete();
// Vector of integers parallel to mc points
-
std::vector<std::vector<int>> pointNhits;
pointNhits.resize(nMcEvents);
for (int iE = 0; iE < nMcEvents; iE++) {
@@ -439,7 +451,7 @@ void CbmTrackerInputQaTrd::ResolutionQa()
for (Int_t iHit = 0; iHit < nHits; iHit++) {
- CbmTrdHit* hit = dynamic_cast<CbmTrdHit*>(fHits->At(iHit));
+ const CbmTrdHit* hit = dynamic_cast<const CbmTrdHit*>(fHits->At(iHit));
if (!hit) {
LOG(error) << "TRD hit N " << iHit << " doesn't exist";
return;
@@ -553,6 +565,9 @@ void CbmTrackerInputQaTrd::ResolutionQa()
// skip hits from the noise digis
if (bestLink.GetIndex() < 0) { continue; }
+ // construct QA object
+ CbmTrdHitMC* hmc = new ((*fHitsMc)[imc++]) CbmTrdHitMC(*hit);
+
CbmTrdPoint* p = dynamic_cast<CbmTrdPoint*>(fMcPoints->Get(bestLink));
if (p == nullptr) {
LOG(error) << "link points to a non-existing MC point";
@@ -631,15 +646,16 @@ void CbmTrackerInputQaTrd::ResolutionQa()
mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
}
}
+ hmc->AddPoint(p, t0, mass);
constexpr double speedOfLight = 29.979246; // cm/ns
TVector3 mom3;
p->Momentum(mom3);
t += dz / (pz * speedOfLight) * sqrt(mass * mass + mom3.Mag2());
- double du = hit->GetX() - x;
- double dv = hit->GetY() - y;
- double dt = hit->GetTime() - t;
+ double du = hmc->GetDx(); //hit->GetX() - x;
+ double dv = hmc->GetDy(); //hit->GetY() - y;
+ double dt = hmc->GetDt(); //hit->GetTime() - t;
double su = hit->GetDx();
double sv = hit->GetDy();
double st = hit->GetTimeError();
@@ -673,7 +689,7 @@ void CbmTrackerInputQaTrd::ResolutionQa()
fh1DpullT.Fill(dt / st);
}
else {
- fh2DpullX.Fill(du / su);
+ fh2DpullX.Fill(6.2 * du / su);
fh2DpullY.Fill(dv / sv);
fh2DpullT.Fill(dt / st);
}
diff --git a/reco/L1/qa/CbmTrackerInputQaTrd.h b/reco/L1/qa/CbmTrackerInputQaTrd.h
index c7c6f5aad8..3404ed5b9e 100644
--- a/reco/L1/qa/CbmTrackerInputQaTrd.h
+++ b/reco/L1/qa/CbmTrackerInputQaTrd.h
@@ -106,6 +106,7 @@ private:
TClonesArray* fHitMatches {nullptr};
/// Output
+ TClonesArray* fHitsMc {nullptr};
TFolder fOutFolder {"TrackerInputQaTrd", "TrackerInputQaTrd"}; /// output folder with histos and canvases
@@ -118,8 +119,8 @@ private:
CbmQaHist<TH1D> fh1DresidualV {"h1DresidualV", "Trd1D: Residual V;(V_{reco} - V_{MC})(cm)", 100, -10, 10};
CbmQaHist<TH1D> fh1DresidualT {"h1DresidualT", "Trd1D: Residual T;(T_{reco} - T_{MC})(ns)", 100, -50, 50};
- CbmQaHist<TH1D> fh2DresidualX {"h2DresidualX", "Trd2D: Residual X;(X_{reco} - X_{MC})(cm)", 100, -5, 5};
- CbmQaHist<TH1D> fh2DresidualY {"h2DresidualY", "Trd2D: Residual Y;(Y_{reco} - Y_{MC})(cm)", 100, -5, 5};
+ CbmQaHist<TH1D> fh2DresidualX {"h2DresidualX", "Trd2D: Residual X;(X_{reco} - X_{MC})(cm)", 1000, -1, 1};
+ CbmQaHist<TH1D> fh2DresidualY {"h2DresidualY", "Trd2D: Residual Y;(Y_{reco} - Y_{MC})(cm)", 1000, -2, 2};
CbmQaHist<TH1D> fh2DresidualT {"h2DresidualT", "Trd2D: Residual T;(T_{reco} - T_{MC})(ns)", 100, -1000, 1000};
/// Histogram for PULL Distribution
@@ -127,8 +128,8 @@ private:
CbmQaHist<TH1D> fh1DpullV {"h1DpullV", "Trd1D: Pull V;(V_{reco} - V_{MC}) / #sigmaV_{reco}", 100, -5, 5};
CbmQaHist<TH1D> fh1DpullT {"h1DpullT", "Trd1D: Pull T;(T_{reco} - T_{MC}) / #sigmaT_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fh2DpullX {"h2DpullX", "Trd2D: Pull X;(X_{reco} - X_{MC}) / #sigmaX_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fh2DpullY {"h2DpullY", "Trd2D: Pull Y;(Y_{reco} - Y_{MC}) / #sigmaY_{reco}", 100, -5, 5};
+ CbmQaHist<TH1D> fh2DpullX {"h2DpullX", "Trd2D: Pull X;(X_{reco} - X_{MC}) / #sigmaX_{reco}", 1000, -5, 5};
+ CbmQaHist<TH1D> fh2DpullY {"h2DpullY", "Trd2D: Pull Y;(Y_{reco} - Y_{MC}) / #sigmaY_{reco}", 1000, -5, 5};
CbmQaHist<TH1D> fh2DpullT {"h2DpullT", "Trd2D: Pull T;(T_{reco} - T_{MC}) / #sigmaT_{reco}", 100, -5, 5};
/// List of the above histogramms
diff --git a/reco/detectors/trd/CMakeLists.txt b/reco/detectors/trd/CMakeLists.txt
index 7eed7e333f..87f1dd0b5c 100644
--- a/reco/detectors/trd/CMakeLists.txt
+++ b/reco/detectors/trd/CMakeLists.txt
@@ -5,6 +5,7 @@
Set(INCLUDE_DIRECTORIES
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/qa
+ ${CMAKE_CURRENT_SOURCE_DIR}/qa/data
${CMAKE_CURRENT_SOURCE_DIR}/pid
${CMAKE_CURRENT_SOURCE_DIR}/rawToDigiMethods
${CMAKE_CURRENT_SOURCE_DIR}/unpack
@@ -74,6 +75,8 @@ qa/CbmTrdQa.cxx
qa/CbmTrdRecoQa.cxx
qa/CbmTrdTracksPidQa.cxx
+qa/data/CbmTrdHitMC.cxx
+
pid/CbmTrdElectronsTrainAnn.cxx
pid/CbmTrdSetTracksPidWkn.cxx
pid/CbmTrdSetTracksPidModWkn.cxx
diff --git a/reco/detectors/trd/CbmTrdRecoLinkDef.h b/reco/detectors/trd/CbmTrdRecoLinkDef.h
index 122fb9cdcc..3c0121a4ad 100644
--- a/reco/detectors/trd/CbmTrdRecoLinkDef.h
+++ b/reco/detectors/trd/CbmTrdRecoLinkDef.h
@@ -27,6 +27,7 @@
#pragma link C++ class CbmTrdQa + ;
#pragma link C++ class CbmTrdRecoQa + ;
#pragma link C++ class CbmTrdTracksPidQa + ;
+#pragma link C++ class CbmTrdHitMC + ;
#pragma link C++ class std::vector<std::pair<unsigned long, unsigned long>> + ;
#pragma link C++ class CbmTrdUnpackAlgoBaseR + ;
diff --git a/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx b/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
new file mode 100644
index 0000000000..48406d78e6
--- /dev/null
+++ b/reco/detectors/trd/qa/data/CbmTrdHitMC.cxx
@@ -0,0 +1,66 @@
+/* Copyright (C) 2018-2022 Horia Hulubei National Institute of Physics and Nuclear Engineering, Bucharest
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Alexandru Bercuci [committer] */
+
+#include "CbmTrdHitMC.h"
+
+#include <TVector3.h>
+
+//_____________________________________________________________________
+CbmTrdHitMC::CbmTrdHitMC() : CbmTrdHit() {}
+
+//_____________________________________________________________________
+CbmTrdHitMC::CbmTrdHitMC(const CbmTrdHit& h) : CbmTrdHit(h) {}
+
+//_____________________________________________________________________
+CbmTrdHitMC::~CbmTrdHitMC() {}
+
+//_____________________________________________________________________
+size_t CbmTrdHitMC::AddPoint(const CbmTrdPoint* p, double t, double m)
+{
+ fTrdPoints.push_back(std::make_tuple(CbmTrdPoint(*p), t, m));
+ return fTrdPoints.size();
+}
+
+//_____________________________________________________________________
+const CbmTrdPoint* CbmTrdHitMC::GetPoint(uint idx) const
+{
+ if (idx >= fTrdPoints.size()) return nullptr;
+ return &std::get<0>(fTrdPoints[idx]);
+}
+
+//_____________________________________________________________________
+double CbmTrdHitMC::GetDx() const
+{
+ const CbmTrdPoint* p(nullptr);
+ if (!(p = GetPoint())) return -999;
+ double dz(GetZ() - p->GetZ()), x(p->GetX() + dz * p->GetPx() / p->GetPz());
+ return GetX() - x;
+}
+
+//_____________________________________________________________________
+double CbmTrdHitMC::GetDy() const
+{
+ const CbmTrdPoint* p(nullptr);
+ if (!(p = GetPoint())) return -999;
+ double dz(GetZ() - p->GetZ()), y(p->GetY() + dz * p->GetPy() / p->GetPz());
+ return GetY() - y;
+}
+
+//_____________________________________________________________________
+double CbmTrdHitMC::GetDt() const
+{
+ constexpr double speedOfLight = 29.979246; // cm/ns
+ 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());
+
+ TVector3 mom3;
+ p->Momentum(mom3);
+ t += dz / (p->GetPz() * speedOfLight) * sqrt(mass * mass + mom3.Mag2());
+ return GetTime() - t;
+}
+
+ClassImp(CbmTrdHitMC)
diff --git a/reco/detectors/trd/qa/data/CbmTrdHitMC.h b/reco/detectors/trd/qa/data/CbmTrdHitMC.h
new file mode 100644
index 0000000000..015b795d2b
--- /dev/null
+++ b/reco/detectors/trd/qa/data/CbmTrdHitMC.h
@@ -0,0 +1,60 @@
+/* Copyright (C) 2018-2022 Horia Hulubei National Institute of Physics and Nuclear Engineering, Bucharest
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Alexandru Bercuci [committer] */
+
+#ifndef CBMTRDHITMC_H
+#define CBMTRDHITMC_H
+
+#include "CbmTrdHit.h"
+#include "CbmTrdPoint.h"
+
+#include <vector> // for vector
+
+/** \class CbmTrdHitMC
+ * \brief TRD hit to MC point correlation class
+ * \author Alexandru Bercuci <abercuci@niham.nipne.ro>
+ * \date 02.05.2022
+ *
+ * The class
+ */
+
+class CbmTrdHitMC : public CbmTrdHit {
+public:
+ /** \brief Default constructor.*/
+ CbmTrdHitMC();
+ /** \brief Standard constructor.*/
+ CbmTrdHitMC(const CbmTrdHit& hit);
+ /** \brief Destructor. */
+ virtual ~CbmTrdHitMC();
+
+ /** \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);
+
+ /** \brief Retrieve 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;
+
+private:
+ /** \brief Copy Constructor.*/
+ CbmTrdHitMC(const CbmTrdHitMC&) = default;
+ /** \brief Assignment operator.*/
+ CbmTrdHitMC& operator=(const CbmTrdHitMC&) = default;
+
+ std::vector<std::tuple<CbmTrdPoint, double, double>> fTrdPoints = {};
+ ClassDef(CbmTrdHitMC, 1) // Hit to MC point data correlation
+};
+
+#endif
--
GitLab