GPU: Improve debug-dump of TPC cluster to track residuals

Author: davidrohr

GPU/GPUTracking/Debug/GPUROOTDump.h

@@ -14,13 +14,15 @@
 #ifndef GPUROOTDUMP_H
 #define GPUROOTDUMP_H
 
-#include "GPUROOTDumpCore.h"
-
+#include "GPUCommonDef.h"
 #if !defined(GPUCA_NO_ROOT) && !defined(GPUCA_GPUCODE)
+#include "GPUROOTDumpCore.h"
 #include <TTree.h>
 #include <TNtuple.h>
 #include <memory>
 #include <stdexcept>
+#else
+class TNtuple;
 #endif
 
 namespace GPUCA_NAMESPACE
@@ -55,11 +57,6 @@ class GPUROOTDump : public GPUROOTDumpBase
   {
     static GPUROOTDump<T> instance(name);
   }
-  GPUROOTDump(const char* name)
-  {
-    mTree = new TTree(name, name);
-    mTree->Branch(name, &mObj);
-  }
 
   void write() override { mTree->Write(); }
 
@@ -70,6 +67,11 @@ class GPUROOTDump : public GPUROOTDumpBase
   }
 
  private:
+  GPUROOTDump(const char* name)
+  {
+    mTree = new TTree(name, name);
+    mTree->Branch(name, &mObj);
+  }
   TTree* mTree = nullptr;
   T mObj;
 };
@@ -83,10 +85,6 @@ class GPUROOTDump<TNtuple> : public GPUROOTDumpBase
     static GPUROOTDump<TNtuple> instance(name, options);
     return instance;
   }
-  GPUROOTDump(const char* name, const char* options)
-  {
-    mNTuple = new TNtuple(name, name, options);
-  }
 
   void write() override { mNTuple->Write(); }
 
@@ -97,14 +95,26 @@ class GPUROOTDump<TNtuple> : public GPUROOTDumpBase
   }
 
  private:
+  GPUROOTDump(const char* name, const char* options)
+  {
+    mNTuple = new TNtuple(name, name, options);
+  }
   TNtuple* mNTuple;
 };
 #else
 template <class T>
 class GPUROOTDump
 {
  public:
-  GPUROOTDump() = delete;
+  template <typename... Args>
+  GPUd() void Fill(Args... args) const
+  {
+  }
+  template <typename... Args>
+  GPUd() static GPUROOTDump<T>& get(Args... args)
+  {
+    return *(GPUROOTDump<T>*)(size_t)(1024); // Will never be used, return just some reference
+  }
 };
 #endif
 } // namespace gpu

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx

@@ -13,8 +13,12 @@
 
 #define GPUCA_CADEBUG 0
 #define DEBUG_SINGLE_TRACK -1
+#define EXTRACT_RESIDUALS 0
+
+#if EXTRACT_RESIDUALS == 1
+#include "GPUROOTDump.h"
+#endif
 
-#include "GPUTPCClusterErrorStat.h"
 #include "GPUTPCDef.h"
 #include "GPUTPCGMTrackParam.h"
 #include "GPUTPCGMPhysicalTrackModel.h"
@@ -58,8 +62,6 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int iT
 {
   const GPUParam& GPUrestrict() param = merger->Param();
 
-  GPUTPCClusterErrorStat errorStat(N);
-
   GPUdEdx dEdx;
   GPUTPCGMPropagator prop;
   gputpcgmmergertypes::InterpolationErrors interpolation;
@@ -266,14 +268,32 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int iT
         continue;
       }
       CADEBUG(printf("\n"));
-      errorStat.Fill(xx, yy, zz, prop.GetAlpha(), mX, mP, mC, ihit, iWay);
 
       int retVal;
       float threshold = 3.f + (lastUpdateX >= 0 ? (CAMath::Abs(mX - lastUpdateX) / 2) : 0.f);
       if (mNDF > 5 && (CAMath::Abs(yy - mP[0]) > threshold || CAMath::Abs(zz - mP[1]) > threshold)) {
         retVal = 2;
       } else {
         char rejectChi2 = attempt ? 0 : ((param.rec.mergerInterpolateErrors && CAMath::Abs(ihit - ihitMergeFirst) <= 1) ? (refit ? (2 + ((nWays - iWay) & 1)) : 0) : (allowModification && goodRows > 5));
+#if EXTRACT_RESIDUALS == 1
+        if (iWay == nWays - 1 && interpolation.hit[ihit].errorY > (GPUCA_MERGER_INTERPOLATION_ERROR_TYPE)0) {
+          const float Iz0 = interpolation.hit[ihit].posY - mP[0];
+          const float Iz1 = interpolation.hit[ihit].posZ - mP[1];
+          float Iw0 = mC[2] + (float)interpolation.hit[ihit].errorZ;
+          float Iw2 = mC[0] + (float)interpolation.hit[ihit].errorY;
+          float Idet1 = 1.f / CAMath::Max(1e-10f, Iw0 * Iw2 - mC[1] * mC[1]);
+          const float Ik00 = (mC[0] * Iw0 + mC[1] * mC[1]) * Idet1;
+          const float Ik01 = (mC[0] * mC[1] + mC[1] * Iw2) * Idet1;
+          const float Ik10 = (mC[1] * Iw0 + mC[2] * mC[1]) * Idet1;
+          const float Ik11 = (mC[1] * mC[1] + mC[2] * Iw2) * Idet1;
+          const float ImP0 = mP[0] + Ik00 * Iz0 + Ik01 * Iz1;
+          const float ImP1 = mP[1] + Ik10 * Iz0 + Ik11 * Iz1;
+          const float ImC0 = mC[0] - Ik00 * mC[0] + Ik01 * mC[1];
+          const float ImC2 = mC[2] - Ik10 * mC[1] + Ik11 * mC[2];
+          auto& tup = GPUROOTDump<TNtuple>::get("clusterres", "row:clX:clY:clZ:angle:trkX:trkY:trkZ:trkSinPhi:trkDzDs:trkQPt:trkSigmaY2:trkSigmaZ2trkSigmaQPt2");
+          tup.Fill((float)clusters[ihit].row, xx, yy, zz, clAlpha, mX, ImP0, ImP1, mP[2], mP[3], mP[4], ImC0, ImC2, mC[14]);
+        }
+#endif
         retVal = prop.Update(yy, zz, clusters[ihit].row, param, clusterState, rejectChi2, &interpolation.hit[ihit], refit);
       }
       // clang-format off