GPU QA: Several improvements

GPU/GPUTracking/Base/GPUReconstruction.cxx

@@ -100,9 +100,6 @@ GPUReconstruction::GPUReconstruction(const GPUSettingsDeviceBackend& cfg) : mHos
     processors()->tpcNNClusterer[i].mISector = i;
 #endif
   }
-#ifndef GPUCA_NO_ROOT
-  mROOTDump = GPUROOTDumpCore::getAndCreate();
-#endif
 }
 
 GPUReconstruction::~GPUReconstruction()
@@ -135,6 +132,11 @@ int32_t GPUReconstruction::Init()
   if (mMaster) {
     throw std::runtime_error("Must not call init on slave!");
   }
+#ifndef GPUCA_NO_ROOT
+  if (!mROOTDump) {
+    mROOTDump = GPUROOTDumpCore::getAndCreate(GetProcessingSettings().ROOTDumpFile.c_str());
+  }
+#endif
   int32_t retVal = InitPhaseBeforeDevice();
   if (retVal) {
     return retVal;

GPU/GPUTracking/Debug/GPUROOTDump.h

@@ -50,8 +50,11 @@ struct internal_Branch<TTree> {
 };
 } // namespace
 
+template <typename... Args>
+class GPUROOTDump;
+
 template <class T, typename... Args>
-class GPUROOTDump : public GPUROOTDump<Args...>
+class GPUROOTDump<T, Args...> : public GPUROOTDump<Args...>
 {
  public:
   template <typename... Names>
@@ -65,59 +68,53 @@ class GPUROOTDump : public GPUROOTDump<Args...>
   {
     return GPUROOTDump<T, Args...>(name1, names...);
   }
-  void Fill(const T& o, Args... args)
+
+  void Fill(const T& o, const Args&... args)
+  {
+    stdspinlock spinlock(GPUROOTDumpBase::mMutex);
+    FillInternal(o, args...);
+  }
+
+ protected:
+  void FillInternal(const T& o, const Args&... args)
   {
     mObj = o;
     GPUROOTDump<Args...>::Fill(args...);
   }
 
- protected:
   using GPUROOTDump<Args...>::mTree;
   template <typename... Names>
   GPUROOTDump(const char* name1, Names... names) : GPUROOTDump<Args...>(names...)
   {
+    stdspinlock spinlock(GPUROOTDumpBase::mMutex);
     mTree->Branch(name1, &mObj);
   }
 
  private:
   T mObj;
 };
 
-template <class T>
-class GPUROOTDump<T> : public GPUROOTDumpBase
+template <>
+class GPUROOTDump<> : public GPUROOTDumpBase
 {
  public:
-  static GPUROOTDump<T>& get(const char* name) // return always the same instance, identified by template
-  {
-    static GPUROOTDump<T> instance(name);
-    return instance;
-  }
-  static GPUROOTDump<T> getNew(const char* name) // return new individual instance
-  {
-    return GPUROOTDump<T>(name);
-  }
-
   void write() override { mTree->Write(); }
 
-  void Fill(const T& o)
+ protected:
+  void Fill()
   {
-    mObj = o;
     mTree->Fill();
   }
 
- protected:
   GPUROOTDump(const char* name1, const char* nameTree = nullptr)
   {
     if (nameTree == nullptr) {
       nameTree = name1;
     }
+    stdspinlock spinlock(GPUROOTDumpBase::mMutex);
     mTree = new TTree(nameTree, nameTree);
-    mTree->Branch(name1, &mObj);
   }
   TTree* mTree = nullptr;
-
- private:
-  T mObj;
 };
 
 template <>
@@ -137,14 +134,16 @@ class GPUROOTDump<TNtuple> : public GPUROOTDumpBase
   void write() override { mNTuple->Write(); }
 
   template <typename... Args>
-  void Fill(Args... args)
+  void Fill(const Args&... args)
   {
+    stdspinlock spinlock(GPUROOTDumpBase::mMutex);
     mNTuple->Fill(args...);
   }
 
  private:
   GPUROOTDump(const char* name, const char* options)
   {
+    stdspinlock spinlock(GPUROOTDumpBase::mMutex);
     mNTuple = new TNtuple(name, name, options);
   }
   TNtuple* mNTuple;
@@ -155,18 +154,18 @@ class GPUROOTDump
 {
  public:
   template <typename... Names>
-  GPUd() void Fill(Args... args) const
+  GPUd() static void Fill(Args... args)
   {
   }
   template <typename... Names>
   GPUd() static GPUROOTDump<Args...>& get(Args... args)
   {
-    return *(GPUROOTDump<Args...>*)(size_t)(1024); // Will never be used, return just some reference, which must not be nullptr by specification
+    return GPUROOTDump<Args...>();
   }
   template <typename... Names>
   GPUd() static GPUROOTDump<Args...>& getNew(Args... args)
   {
-    return *(GPUROOTDump<Args...>*)(size_t)(1024); // Will never be used, return just some reference, which must not be nullptr by specification
+    return GPUROOTDump<Args...>();
   }
 };
 #endif

GPU/GPUTracking/Debug/GPUROOTDumpCore.cxx

@@ -22,6 +22,7 @@
 using namespace o2::gpu;
 
 std::weak_ptr<GPUROOTDumpCore> GPUROOTDumpCore::sInstance;
+std::atomic_flag GPUROOTDumpBase::mMutex = ATOMIC_FLAG_INIT;
 
 GPUROOTDumpCore::GPUROOTDumpCore(GPUROOTDumpCore::GPUROOTDumpCorePrivate)
 {
@@ -38,7 +39,7 @@ GPUROOTDumpCore::~GPUROOTDumpCore()
   }
 }
 
-std::shared_ptr<GPUROOTDumpCore> GPUROOTDumpCore::getAndCreate()
+std::shared_ptr<GPUROOTDumpCore> GPUROOTDumpCore::getAndCreate(const char* filename)
 {
   static std::atomic_flag lock = ATOMIC_FLAG_INIT;
   while (lock.test_and_set(std::memory_order_acquire)) {
@@ -48,6 +49,10 @@ std::shared_ptr<GPUROOTDumpCore> GPUROOTDumpCore::getAndCreate()
     retVal = std::make_shared<GPUROOTDumpCore>(GPUROOTDumpCorePrivate());
     sInstance = retVal;
   }
+  if (*filename != 0 && retVal->mFileName != "" && retVal->mFileName != filename) {
+    throw std::runtime_error("GPUROOTDump reinitialized with different file name");
+  }
+  retVal->mFileName = filename;
   lock.clear(std::memory_order_release);
   return retVal;
 }
@@ -60,8 +65,11 @@ GPUROOTDumpBase::GPUROOTDumpBase()
   }
   p->mBranches.emplace_back(this);
   if (!p->mFile) {
-    std::remove("gpudebug.root");
-    p->mFile.reset(new TFile("gpudebug.root", "recreate"));
+    if (p->mFileName == "") {
+      throw std::runtime_error("GPUROOTDump output file name not set");
+    }
+    std::remove(p->mFileName.c_str());
+    p->mFile.reset(new TFile(p->mFileName.c_str(), "recreate"));
   }
   p->mFile->cd();
 }

GPU/GPUTracking/Debug/GPUROOTDumpCore.h

@@ -16,8 +16,10 @@
 #define GPUROOTDUMPCORE_H
 
 #include "GPUCommonDef.h"
+#include "utils/stdspinlock.h"
 #include <memory>
 #include <vector>
+#include <string>
 
 class TFile;
 
@@ -32,6 +34,7 @@ class GPUROOTDumpBase
 
  protected:
   GPUROOTDumpBase();
+  static std::atomic_flag mMutex;
   std::weak_ptr<GPUROOTDumpCore> mCore;
 };
 
@@ -52,9 +55,10 @@ class GPUROOTDumpCore
   ~GPUROOTDumpCore();
 
  private:
-  static std::shared_ptr<GPUROOTDumpCore> getAndCreate();
+  static std::shared_ptr<GPUROOTDumpCore> getAndCreate(const char* filename);
   static std::weak_ptr<GPUROOTDumpCore> get() { return sInstance; }
   static std::weak_ptr<GPUROOTDumpCore> sInstance;
+  std::string mFileName;
   std::unique_ptr<TFile> mFile;
   std::vector<GPUROOTDumpBase*> mBranches;
 #endif

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -414,6 +414,7 @@ AddOption(debugOnFailureMaxN, uint32_t, 1, "", 0, "Max number of times to run th
 AddOption(debugOnFailureMaxFiles, uint32_t, 0, "", 0, "Max number of files to have in the target folder")
 AddOption(debugOnFailureMaxSize, uint32_t, 0, "", 0, "Max size of existing dumps in the target folder in GB")
 AddOption(debugOnFailureDirectory, std::string, ".", "", 0, "Target folder for debug / dump")
+AddOption(ROOTDumpFile, std::string, "gpudebug.root", "", 0, "File name for ROOT dump (default gpudebug.root)")
 AddOption(memoryStat, bool, false, "", 0, "Print memory statistics")
 AddVariable(eventDisplay, o2::gpu::GPUDisplayFrontendInterface*, nullptr)
 AddSubConfig(GPUSettingsProcessingRTC, rtc)
@@ -546,7 +547,7 @@ AddOption(filterCharge, int32_t, 0, "", 0, "Filter for positive (+1) or negative
 AddOption(filterPID, int32_t, -1, "", 0, "Filter for Particle Type (0 Electron, 1 Muon, 2 Pion, 3 Kaon, 4 Proton)")
 AddOption(nativeFitResolutions, bool, false, "", 0, "Create resolution histograms in the native fit units (sin(phi), tan(lambda), Q/Pt)")
 AddOption(enableLocalOutput, bool, true, "", 0, "Enable normal output to local PDF files / console")
-AddOption(dumpToROOT, int32_t, 0, "", 0, "Dump all clusters and tracks to a ROOT file, 1 = combined TNTUple dump, 2 = also individual cluster / track branch dump")
+AddOption(dumpToROOTLevel, int32_t, -1, "", 0, "Dump all clusters and tracks to a ROOT file, 1 = combined TNTUple dump, 2 = also individual cluster / track branch dump", def(1))
 AddOption(writeFileExt, std::string, "", "", 0, "Write extra output file with given extension (default ROOT Canvas)", def("root"))
 AddOption(writeMCLabels, bool, false, "", 0, "Store mc labels to file for later matching")
 AddOptionVec(matchMCLabels, std::string, "", 0, "Read labels from files and match them, only process tracks where labels differ")

GPU/GPUTracking/Standalone/Benchmark/standalone.cxx

@@ -212,7 +212,7 @@ int32_t ReadConfiguration(int argc, char** argv)
     configStandalone.rundEdx = false;
     configStandalone.noEvents = true;
   }
-  if (configStandalone.QA.dumpToROOT) {
+  if (configStandalone.QA.dumpToROOTLevel >= 1) {
     configStandalone.proc.outputSharedClusterMap = true;
   }
   if (configStandalone.eventDisplay) {

GPU/GPUTracking/qa/GPUQA.cxx

@@ -539,6 +539,10 @@ int32_t GPUQA::InitQACreateHistograms()
     }
     std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], PT_MIN_CLUST, PT_MAX)};
     createHist(mTrackPt, "tracks_pt", "tracks_pt", AXIS_BINS[4], binsPt.get());
+    for (int32_t i = 0; i < 2; i++) {
+      snprintf(name, 2048, i ? "tracks_dedx_max" : "tracks_dedx_tot");
+      createHist(mTrackdEdx[i], name, name, 200, -3, 1, 1000, 0, i ? 1500 : 5000);
+    }
     const uint32_t maxTime = (mTracking && mTracking->GetParam().continuousMaxTimeBin > 0) ? mTracking->GetParam().continuousMaxTimeBin : constants::TPC_MAX_TIME_BIN_TRIGGERED;
     createHist(mT0[0], "tracks_t0", "tracks_t0", (maxTime + 1) / 10, 0, maxTime);
     createHist(mT0[1], "tracks_t0_res", "tracks_t0_res", 1000, -100, 100);
@@ -1281,7 +1285,7 @@ void GPUQA::RunQA(bool matchOnly, const std::vector<o2::tpc::TrackTPC>* tracksEx
           float s = std::sin(alpha);
           float localY = -info.x * s + info.y * c;
 
-          if (mConfig.dumpToROOT) {
+          if (mConfig.dumpToROOTLevel >= 1) {
             static auto effdump = GPUROOTDump<TNtuple>::getNew("eff", "alpha:x:y:z:mcphi:mceta:mcpt:rec:fake:findable:prim:ncls");
             float localX = info.x * c + info.y * s;
             effdump.Fill(alpha, localX, localY, info.z, mcphi, mceta, mcpt, mRecTracks[iCol][i], mFakeTracks[iCol][i], findable, info.prim, mc2.nWeightCls);
@@ -1745,6 +1749,12 @@ void GPUQA::RunQA(bool matchOnly, const std::vector<o2::tpc::TrackTPC>* tracksEx
         continue;
       }
       mTrackPt->Fill(1.f / fabsf(track.GetParam().GetQPt()));
+      if (mParam->par.dodEdx && mParam->dodEdxEnabled && track.NClusters() >= 60) {
+        const GPUdEdxInfo& trackdEdx = mTracking->GetProcessors()->tpcMerger.MergedTracksdEdx()[i];
+        const float logp = logf(1.f / fabsf(track.GetParam().GetQPt()) * sqrtf(1.f + track.GetParam().GetDzDs() * track.GetParam().GetDzDs()));
+        mTrackdEdx[0]->Fill(logp, trackdEdx.dEdxTotTPC);
+        mTrackdEdx[1]->Fill(logp, trackdEdx.dEdxMaxTPC);
+      }
       mNCl[0]->Fill(track.NClustersFitted());
       uint32_t nClCorrected = 0;
       const auto& trackClusters = mTracking->mIOPtrs.mergedTrackHits;
@@ -1919,7 +1929,7 @@ void GPUQA::RunQA(bool matchOnly, const std::vector<o2::tpc::TrackTPC>* tracksEx
     GPUInfo("QA Time: Cluster Counts:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
   }
 
-  if (mConfig.dumpToROOT && !tracksExternal) {
+  if (mConfig.dumpToROOTLevel >= 1 && !tracksExternal) {
     if (!clNative || !mTracking || !mTracking->mIOPtrs.mergedTrackHitAttachment || !mTracking->mIOPtrs.mergedTracks) {
       throw std::runtime_error("Cannot dump non o2::tpc::clusterNative clusters, need also hit attachmend and GPU tracks");
     }
@@ -1938,7 +1948,7 @@ void GPUQA::RunQA(bool matchOnly, const std::vector<o2::tpc::TrackTPC>* tracksEx
           }
           uint32_t extState = mTracking->mIOPtrs.mergedTrackHitStates ? mTracking->mIOPtrs.mergedTrackHitStates[clid] : 0;
 
-          if (mConfig.dumpToROOT >= 2) {
+          if (mConfig.dumpToROOTLevel >= 2) {
             GPUTPCGMMergedTrack trk;
             GPUTPCGMMergedTrackHit trkHit;
             memset((void*)&trk, 0, sizeof(trk));
@@ -2255,12 +2265,15 @@ int32_t GPUQA::DrawQAHistograms(TObjArray* qcout)
 
     // Create Canvas for track statistic histos
     if (mQATasks & taskTrackStatistics) {
-      mCTrackPt = createGarbageCollected<TCanvas>("ctrackspt", "ctrackspt", 0, 0, 700, 700. * 2. / 3.);
-      mCTrackPt->cd();
-      mPTrackPt = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
-      mPTrackPt->Draw();
-      mLTrackPt = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.5, 0.93 - (0.93 - 0.86) / 2. * (float)ConfigNumInputs, 0.98, 0.949);
-      SetLegend(mLTrackPt, true);
+      for (int32_t i = 0; i < 3; i++) {
+        snprintf(name, 2048, "ctracks%s", i ? (i == 2 ? "dedxmax" : "dedxtot") : "pt");
+        mCTracks[i] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
+        mCTracks[i]->cd();
+        mPTracks[i] = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
+        mPTracks[i]->Draw();
+        mLTracks[i] = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.5, 0.93 - (0.93 - 0.86) / 2. * (float)ConfigNumInputs, 0.98, 0.949);
+        SetLegend(mLTracks[i], true);
+      }
 
       for (int32_t i = 0; i < 2; i++) {
         snprintf(name, 2048, "ctrackst0%d", i);
@@ -2916,6 +2929,7 @@ int32_t GPUQA::DrawQAHistograms(TObjArray* qcout)
 
   if (mQATasks & taskTrackStatistics) {
     // Process track statistic histograms
+
     float tmpMax = 0.;
     for (int32_t k = 0; k < ConfigNumInputs; k++) { // TODO: Simplify this drawing, avoid copy&paste
       TH1F* e = mTrackPt;
@@ -2927,8 +2941,8 @@ int32_t GPUQA::DrawQAHistograms(TObjArray* qcout)
         tmpMax = e->GetMaximum();
       }
     }
-    mPTrackPt->cd();
-    mPTrackPt->SetLogx();
+    mPTracks[0]->cd();
+    mPTracks[0]->SetLogx();
     for (int32_t k = 0; k < ConfigNumInputs; k++) {
       TH1F* e = mTrackPt;
       if (GetHist(e, tin, k, nNewInput) == nullptr) {
@@ -2952,17 +2966,40 @@ int32_t GPUQA::DrawQAHistograms(TObjArray* qcout)
       e->SetLineColor(colorNums[k % COLORCOUNT]);
       e->Draw(k == 0 ? "" : "same");
       GetName(fname, k, mConfig.inputHistogramsOnly);
-      mLTrackPt->AddEntry(e, Form(mConfig.inputHistogramsOnly ? "%s" : "%sTrack #it{p}_{T}", fname), "l");
+      mLTracks[0]->AddEntry(e, Form(mConfig.inputHistogramsOnly ? "%s" : "%sTrack #it{p}_{T}", fname), "l");
     }
-    mLTrackPt->Draw();
+    mLTracks[0]->Draw();
     doPerfFigure(0.63, 0.7, 0.030);
-    mCTrackPt->cd();
-    mCTrackPt->Print(Form("%s/tracks.pdf", mConfig.plotsDir.c_str()));
+    mCTracks[0]->cd();
+    mCTracks[0]->Print(Form("%s/tracks%s.pdf", mConfig.plotsDir.c_str(), "pt"));
     if (mConfig.writeFileExt != "") {
-      mCTrackPt->Print(Form("%s/tracks.%s", mConfig.plotsDir.c_str(), mConfig.writeFileExt.c_str()));
+      mCTracks[0]->Print(Form("%s/tracks%s.%s", mConfig.plotsDir.c_str(), "pt", mConfig.writeFileExt.c_str()));
     }
 
     for (int32_t i = 0; i < 2; i++) {
+      mPTracks[1 + i]->cd();
+      {
+        TH2F* e = mTrackdEdx[i];
+        if (tout && !mConfig.inputHistogramsOnly) {
+          e->Write();
+        }
+        // e->SetStats(kFALSE);
+        e->SetTitle(mConfig.plotsNoTitle ? "" : (i ? "Track dE/dx (Max)" : "Track dE/dx (Tot)"));
+        e->GetYaxis()->SetTitle(i ? "dE/dx (max)" : "dE/dx (tot)");
+        e->GetXaxis()->SetTitle("log(#it{p})");
+        e->GetXaxis()->SetTitleOffset(1.2);
+        if (qcout) {
+          qcout->Add(e);
+        }
+        e->SetOption("colz");
+        e->Draw();
+      }
+      mCTracks[1 + i]->cd();
+      mCTracks[1 + i]->Print(Form("%s/tracks%s.pdf", mConfig.plotsDir.c_str(), i ? "dedx_max" : "dedx_tot"));
+      if (mConfig.writeFileExt != "") {
+        mCTracks[1 + i]->Print(Form("%s/tracks%s.%s", mConfig.plotsDir.c_str(), i ? "dedx_max" : "dedx_tot", mConfig.writeFileExt.c_str()));
+      }
+
       tmpMax = 0.;
       for (int32_t k = 0; k < ConfigNumInputs; k++) {
         TH1F* e = mT0[i];