store digits on the QA structure to complete the history of

reconstruction

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();
   }
 

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
 };