GPU QA: Add option to dump / read also O2 MC track information

Author: davidrohr

GPU/GPUTracking/Global/GPUChainTracking.cxx

@@ -405,6 +405,9 @@ int GPUChainTracking::PrepareEvent()
 
 int GPUChainTracking::ForceInitQA()
 {
+  if (!mQA) {
+    mQA.reset(new GPUQA(this));
+  }
   return mQA->InitQA();
 }
 

GPU/GPUTracking/Interface/GPUO2Interface.cxx

@@ -98,6 +98,13 @@ int GPUO2Interface::RunTracking(GPUTrackingInOutPointers* data, GPUInterfaceOutp
     mChain->DumpData(fname);
     if (nEvent == 0) {
       mRec->DumpSettings();
+#ifdef GPUCA_BUILD_QA
+      if (mConfig->configProcessing.runMC) {
+        mChain->ForceInitQA();
+        sprintf(fname, "mc.%d.dump", nEvent);
+        mChain->GetQA()->DumpO2MCData(fname);
+      }
+#endif
     }
     nEvent++;
     if (mConfig->configInterface.dumpEvents >= 2) {

GPU/GPUTracking/Standalone/Benchmark/standalone.cxx

@@ -492,6 +492,18 @@ int ReadEvent(int n)
   if (r) {
     return r;
   }
+#if defined(GPUCA_TPC_GEOMETRY_O2) && defined(GPUCA_BUILD_QA) && !defined(GPUCA_O2_LIB)
+  if (configStandalone.proc.runQA && !configStandalone.QA.noMC) {
+    chainTracking->ForceInitQA();
+    snprintf(filename, 256, "events/%s/mc.%d.dump", configStandalone.EventsDir, n);
+    if (chainTracking->GetQA()->ReadO2MCData(filename)) {
+      snprintf(filename, 256, "events/%s/mc.%d.dump", configStandalone.EventsDir, 0);
+      if (chainTracking->GetQA()->ReadO2MCData(filename)) {
+        throw std::runtime_error("Error reading O2 MC dump");
+      }
+    }
+  }
+#endif
   if (chainTracking->mIOPtrs.clustersNative && (configStandalone.TF.bunchSim || configStandalone.TF.nMerge || !configStandalone.runTransformation)) {
     if (configStandalone.proc.debugLevel >= 2) {
       printf("Converting Native to Legacy ClusterData for overlaying - WARNING: No raw clusters produced - Compression etc will not run!!!\n");

GPU/GPUTracking/qa/GPUQA.cxx

@@ -204,12 +204,12 @@ static constexpr Color_t defaultColorNUms[COLORCOUNT] = {kRed, kBlue, kGreen, kM
 
 #define TRACK_EXPECTED_REFERENCE_X_DEFAULT 81
 #ifdef GPUCA_TPC_GEOMETRY_O2
-inline unsigned int GPUQA::GetNMCCollissions()
+inline unsigned int GPUQA::GetNMCCollissions() const
 {
   return mNColTracks.size();
 }
-inline unsigned int GPUQA::GetNMCTracks(int iCol) { return mNColTracks[iCol]; }
-inline unsigned int GPUQA::GetNMCLabels() { return mClNative->clustersMCTruth ? mClNative->clustersMCTruth->getIndexedSize() : 0; }
+inline unsigned int GPUQA::GetNMCTracks(int iCol) const { return mNColTracks[iCol]; }
+inline unsigned int GPUQA::GetNMCLabels() const { return mClNative->clustersMCTruth ? mClNative->clustersMCTruth->getIndexedSize() : 0; }
 inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(unsigned int iTrk, unsigned int iCol) { return mMCInfos[iCol][iTrk]; }
 inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(const mcLabel_t& label) { return mMCInfos[label.getEventID()][label.getTrackID()]; }
 inline GPUQA::mcLabels_t GPUQA::GetMCLabel(unsigned int i) { return mClNative->clustersMCTruth->getLabels(i); }
@@ -231,9 +231,9 @@ inline GPUQA::mcLabelI_t::mcLabelI_t(const GPUQA::mcLabel_t& l) : track(l.fMCID)
 {
 }
 inline bool GPUQA::mcLabelI_t::operator==(const GPUQA::mcLabel_t& l) { return AbsLabelID(track) == l.fMCID; }
-inline unsigned int GPUQA::GetNMCCollissions() { return 1; }
-inline unsigned int GPUQA::GetNMCTracks(int iCol) { return mTracking->mIOPtrs.nMCInfosTPC; }
-inline unsigned int GPUQA::GetNMCLabels() { return mTracking->mIOPtrs.nMCLabelsTPC; }
+inline unsigned int GPUQA::GetNMCCollissions() const { return 1; }
+inline unsigned int GPUQA::GetNMCTracks(int iCol) const { return mTracking->mIOPtrs.nMCInfosTPC; }
+inline unsigned int GPUQA::GetNMCLabels() const { return mTracking->mIOPtrs.nMCLabelsTPC; }
 inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(unsigned int iTrk, unsigned int iCol) { return mTracking->mIOPtrs.mcInfosTPC[AbsLabelID(iTrk)]; }
 inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(const mcLabel_t& label) { return GetMCTrack(label.fMCID, 0); }
 inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(const mcLabelI_t& label) { return GetMCTrack(label.track, 0); }
@@ -555,6 +555,147 @@ int GPUQA::loadHistograms(std::vector<TH1F>& i1, std::vector<TH2F>& i2, std::vec
   return 0;
 }
 
+void GPUQA::DumpO2MCData(const char* filename) const
+{
+  FILE* fp = fopen(filename, "w+b");
+  if (fp == nullptr) {
+    return;
+  }
+  unsigned int n = GetNMCCollissions();
+  fwrite(&n, sizeof(n), 1, fp);
+  for (unsigned int i = 0; i < n; i++) {
+    unsigned int nn = GetNMCTracks(i);
+    fwrite(&nn, sizeof(nn), 1, fp);
+    fwrite(mMCInfos[i].data(), sizeof(mMCInfos[i][0]), nn, fp);
+  }
+  fclose(fp);
+}
+
+int GPUQA::ReadO2MCData(const char* filename)
+{
+  FILE* fp = fopen(filename, "rb");
+  if (fp == nullptr) {
+    return 1;
+  }
+  unsigned int n;
+  unsigned int x;
+  if ((x = fread(&n, sizeof(n), 1, fp)) != 1) {
+    fclose(fp);
+    return 1;
+  }
+  mNColTracks.resize(n);
+  mMCInfos.resize(GetNMCCollissions());
+  for (unsigned int i = 0; i < n; i++) {
+    unsigned int nn = GetNMCTracks(i);
+    if (fread(&nn, sizeof(nn), 1, fp) != 1) {
+      fclose(fp);
+      return 1;
+    }
+    mNColTracks[i] = nn;
+    mMCInfos[i].resize(nn);
+    if (fread(mMCInfos[i].data(), sizeof(mMCInfos[i][0]), nn, fp) != nn) {
+      fclose(fp);
+      return 1;
+    }
+  }
+  fclose(fp);
+  return 0;
+}
+
+void GPUQA::InitO2MCData()
+{
+#ifdef GPUCA_O2_LIB
+  HighResTimer timer;
+  if (mTracking && mTracking->GetProcessingSettings().debugLevel) {
+    GPUInfo("Start reading O2 Track MC information");
+    timer.Start();
+  }
+  static constexpr float PRIM_MAX_T = 0.01f;
+
+  o2::steer::MCKinematicsReader mcReader("collisioncontext.root");
+  int nSimEvents = mcReader.getNEvents(0);
+  mNColTracks.resize(nSimEvents);
+  mMCInfos.resize(GetNMCCollissions());
+  std::vector<int> refId;
+
+  auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root");
+  auto evrec = dc->getEventRecords();
+
+  for (int i = 0; i < nSimEvents; i++) {
+    auto ir = evrec[i];
+    auto ir0 = o2::raw::HBFUtils::Instance().getFirstIRofTF(ir);
+    float timebin = (float)ir.differenceInBC(ir0) / o2::tpc::constants::LHCBCPERTIMEBIN;
+
+    const std::vector<o2::MCTrack>& tracks = mcReader.getTracks(0, i);
+    const std::vector<o2::TrackReference>& trackRefs = mcReader.getTrackRefsByEvent(0, i);
+
+    refId.resize(tracks.size());
+    std::fill(refId.begin(), refId.end(), -1);
+    for (unsigned int j = 0; j < trackRefs.size(); j++) {
+      if (trackRefs[j].getDetectorId() == o2::detectors::DetID::TPC) {
+        int trkId = trackRefs[j].getTrackID();
+        if (refId[trkId] == -1) {
+          refId[trkId] = j;
+        }
+      }
+    }
+    mNColTracks[i] = tracks.size();
+    mMCInfos[i].resize(mNColTracks[i]);
+    for (unsigned int j = 0; j < tracks.size(); j++) {
+      auto& info = mMCInfos[i][j];
+      const auto& trk = tracks[j];
+      TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(trk.GetPdgCode());
+      Int_t pid = -1;
+      if (abs(trk.GetPdgCode()) == kElectron) {
+        pid = 0;
+      }
+      if (abs(trk.GetPdgCode()) == kMuonMinus) {
+        pid = 1;
+      }
+      if (abs(trk.GetPdgCode()) == kPiPlus) {
+        pid = 2;
+      }
+      if (abs(trk.GetPdgCode()) == kKPlus) {
+        pid = 3;
+      }
+      if (abs(trk.GetPdgCode()) == kProton) {
+        pid = 4;
+      }
+
+      info.charge = particle ? particle->Charge() : 0;
+      info.prim = trk.T() < PRIM_MAX_T;
+      info.primDaughters = 0;
+      if (trk.getFirstDaughterTrackId() != -1) {
+        for (int k = trk.getFirstDaughterTrackId(); k <= trk.getLastDaughterTrackId(); k++) {
+          if (tracks[k].T() < PRIM_MAX_T) {
+            info.primDaughters = 1;
+            break;
+          }
+        }
+      }
+      info.pid = pid;
+      info.t0 = timebin;
+      if (refId[j] >= 0) {
+        const auto& trkRef = trackRefs[refId[j]];
+        info.x = trkRef.X();
+        info.y = trkRef.Y();
+        info.z = trkRef.Z();
+        info.pX = trkRef.Px();
+        info.pY = trkRef.Py();
+        info.pZ = trkRef.Pz();
+        info.genRadius = std::sqrt(trk.GetStartVertexCoordinatesX() * trk.GetStartVertexCoordinatesX() + trk.GetStartVertexCoordinatesY() * trk.GetStartVertexCoordinatesY() + trk.GetStartVertexCoordinatesZ() * trk.GetStartVertexCoordinatesZ());
+      } else {
+        info.x = info.y = info.z = info.pX = info.pY = info.pZ = 0;
+        info.genRadius = 0;
+      }
+    }
+  }
+  if (mTracking && mTracking->GetProcessingSettings().debugLevel) {
+    GPUInfo("Finished reading O2 Track MC information (%f seconds)", timer.GetCurrentElapsedTime());
+  }
+#endif
+}
+
 int GPUQA::InitQA(int tasks)
 {
   if (mQAInitialized) {
@@ -583,98 +724,7 @@ int GPUQA::InitQA(int tasks)
 
 #ifdef GPUCA_O2_LIB
   if (!mConfig.noMC) {
-    HighResTimer timer;
-    if (mTracking && mTracking->GetProcessingSettings().debugLevel) {
-      GPUInfo("Start reading O2 Track MC information");
-      timer.Start();
-    }
-    static constexpr float PRIM_MAX_T = 0.01f;
-
-    o2::steer::MCKinematicsReader mcReader("collisioncontext.root");
-    int nSimEvents = mcReader.getNEvents(0);
-    mTrackMCLabelsReverse.resize(nSimEvents);
-    mRecTracks.resize(nSimEvents);
-    mFakeTracks.resize(nSimEvents);
-    mMCParam.resize(nSimEvents);
-    mMCInfos.resize(nSimEvents);
-    mNColTracks.resize(nSimEvents);
-    std::vector<int> refId;
-
-    auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root");
-    auto evrec = dc->getEventRecords();
-
-    for (int i = 0; i < nSimEvents; i++) {
-      auto ir = evrec[i];
-      auto ir0 = o2::raw::HBFUtils::Instance().getFirstIRofTF(ir);
-      float timebin = (float)ir.differenceInBC(ir0) / o2::tpc::constants::LHCBCPERTIMEBIN;
-
-      const std::vector<o2::MCTrack>& tracks = mcReader.getTracks(0, i);
-      const std::vector<o2::TrackReference>& trackRefs = mcReader.getTrackRefsByEvent(0, i);
-
-      refId.resize(tracks.size());
-      std::fill(refId.begin(), refId.end(), -1);
-      for (unsigned int j = 0; j < trackRefs.size(); j++) {
-        if (trackRefs[j].getDetectorId() == o2::detectors::DetID::TPC) {
-          int trkId = trackRefs[j].getTrackID();
-          if (refId[trkId] == -1) {
-            refId[trkId] = j;
-          }
-        }
-      }
-      mNColTracks[i] = tracks.size();
-      mMCInfos[i].resize(tracks.size());
-      for (unsigned int j = 0; j < tracks.size(); j++) {
-        auto& info = mMCInfos[i][j];
-        const auto& trk = tracks[j];
-        TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(trk.GetPdgCode());
-        Int_t pid = -1;
-        if (abs(trk.GetPdgCode()) == kElectron) {
-          pid = 0;
-        }
-        if (abs(trk.GetPdgCode()) == kMuonMinus) {
-          pid = 1;
-        }
-        if (abs(trk.GetPdgCode()) == kPiPlus) {
-          pid = 2;
-        }
-        if (abs(trk.GetPdgCode()) == kKPlus) {
-          pid = 3;
-        }
-        if (abs(trk.GetPdgCode()) == kProton) {
-          pid = 4;
-        }
-
-        info.charge = particle ? particle->Charge() : 0;
-        info.prim = trk.T() < PRIM_MAX_T;
-        info.primDaughters = 0;
-        if (trk.getFirstDaughterTrackId() != -1) {
-          for (int k = trk.getFirstDaughterTrackId(); k <= trk.getLastDaughterTrackId(); k++) {
-            if (tracks[k].T() < PRIM_MAX_T) {
-              info.primDaughters = 1;
-              break;
-            }
-          }
-        }
-        info.pid = pid;
-        info.t0 = timebin;
-        if (refId[j] >= 0) {
-          const auto& trkRef = trackRefs[refId[j]];
-          info.x = trkRef.X();
-          info.y = trkRef.Y();
-          info.z = trkRef.Z();
-          info.pX = trkRef.Px();
-          info.pY = trkRef.Py();
-          info.pZ = trkRef.Pz();
-          info.genRadius = std::sqrt(trk.GetStartVertexCoordinatesX() * trk.GetStartVertexCoordinatesX() + trk.GetStartVertexCoordinatesY() * trk.GetStartVertexCoordinatesY() + trk.GetStartVertexCoordinatesZ() * trk.GetStartVertexCoordinatesZ());
-        } else {
-          info.x = info.y = info.z = info.pX = info.pY = info.pZ = 0;
-          info.genRadius = 0;
-        }
-      }
-    }
-    if (mTracking && mTracking->GetProcessingSettings().debugLevel) {
-      GPUInfo("Finished reading O2 Track MC information (%f seconds)", timer.GetCurrentElapsedTime());
-    }
+    InitO2MCData();
   }
 #endif
 
@@ -729,6 +779,22 @@ void GPUQA::RunQA(bool matchOnly, const std::vector<o2::tpc::TrackTPC>* tracksEx
   }
   mClNative = clNative;
 
+#ifdef GPUCA_TPC_GEOMETRY_O2
+  unsigned int nSimEvents = GetNMCCollissions();
+  if (mTrackMCLabelsReverse.size() < nSimEvents) {
+    mTrackMCLabelsReverse.resize(nSimEvents);
+  }
+  if (mRecTracks.size() < nSimEvents) {
+    mRecTracks.resize(nSimEvents);
+  }
+  if (mFakeTracks.size() < nSimEvents) {
+    mFakeTracks.resize(nSimEvents);
+  }
+  if (mMCParam.size() < nSimEvents) {
+    mMCParam.resize(nSimEvents);
+  }
+#endif
+
   // Initialize Arrays
   unsigned int nReconstructedTracks = 0;
   if (tracksExternal) {

GPU/GPUTracking/qa/GPUQA.h

@@ -52,6 +52,8 @@ class GPUQA
   int HitAttachStatus(int iHit) const { return false; }
   int GetMCTrackLabel(unsigned int trackId) const { return -1; }
   bool clusterRemovable(int attach, bool prot) const { return false; }
+  void DumpO2MCData(const char* filename) const {}
+  int ReadO2MCData(const char* filename) { return 1; }
   static bool QAAvailable() { return false; }
   static bool IsInitialized() { return false; }
 };
@@ -104,6 +106,8 @@ class GPUQA
   int HitAttachStatus(int iHit) const;
   int GetMCTrackLabel(unsigned int trackId) const;
   bool clusterRemovable(int attach, bool prot) const;
+  void DumpO2MCData(const char* filename) const;
+  int ReadO2MCData(const char* filename);
   static bool QAAvailable() { return true; }
   bool IsInitialized() { return mQAInitialized; }
 
@@ -140,6 +144,7 @@ class GPUQA
   };
 
   int InitQACreateHistograms();
+  void InitO2MCData();
   int DoClusterCounts(unsigned long long int* attachClusterCounts, int mode = 0);
   void PrintClusterCount(int mode, int& num, const char* name, unsigned long long int n, unsigned long long int normalization);
 
@@ -188,9 +193,9 @@ class GPUQA
   template <class T>
   static auto& GetMCTrackObj(T& obj, const mcLabelI_t& l);
 
-  unsigned int GetNMCCollissions();
-  unsigned int GetNMCTracks(int iCol);
-  unsigned int GetNMCLabels();
+  unsigned int GetNMCCollissions() const;
+  unsigned int GetNMCTracks(int iCol) const;
+  unsigned int GetNMCLabels() const;
   const mcInfo_t& GetMCTrack(unsigned int iTrk, unsigned int iCol);
   const mcInfo_t& GetMCTrack(const mcLabel_t& label);
   int GetMCLabelNID(const mcLabels_t& label);
@@ -222,14 +227,14 @@ class GPUQA
   std::vector<std::vector<int>> mRecTracks;
   std::vector<std::vector<int>> mFakeTracks;
   std::vector<std::vector<additionalMCParameters>> mMCParam;
-  std::vector<std::vector<mcInfo_t>> mMCInfos;
-  std::vector<int> mNColTracks;
 #else
   std::vector<int> mTrackMCLabelsReverse[1];
   std::vector<int> mRecTracks[1];
   std::vector<int> mFakeTracks[1];
   std::vector<additionalMCParameters> mMCParam[1];
 #endif
+  std::vector<int> mNColTracks;
+  std::vector<std::vector<mcInfo_t>> mMCInfos;
   std::vector<additionalClusterParameters> mClusterParam;
   int mNTotalFakes = 0;