From 7a09d69df5ef1387a4d9db95e52fd71ab153b5e2 Mon Sep 17 00:00:00 2001
From: Dominik Smith <d.smith@gsi.de>
Date: Mon, 7 Feb 2022 15:59:38 +0100
Subject: [PATCH] CbmSeedFinderQa: Completed implementation of time offset
 measurement.

---
 reco/eventbuilder/digis/CbmSeedFinderQa.cxx   | 58 ++++++++++---------
 reco/eventbuilder/digis/CbmSeedFinderQa.h     | 10 +++-
 .../digis/CbmSeedFinderSlidingWindow.cxx      |  5 ++
 .../digis/CbmSeedFinderSlidingWindow.h        |  3 +
 .../digis/CbmTaskBuildRawEvents.cxx           |  3 +
 5 files changed, 52 insertions(+), 27 deletions(-)

diff --git a/reco/eventbuilder/digis/CbmSeedFinderQa.cxx b/reco/eventbuilder/digis/CbmSeedFinderQa.cxx
index 62db95c0e1..9e8b0bf17d 100644
--- a/reco/eventbuilder/digis/CbmSeedFinderQa.cxx
+++ b/reco/eventbuilder/digis/CbmSeedFinderQa.cxx
@@ -29,7 +29,8 @@ CbmSeedFinderQa::CbmSeedFinderQa() : fOutFolder("SeedFinderQA", "Seed finder QA"
   fhCorrectVsFoundNoNoise =
     new TH2I("fhCorrectVsFoundNoNoise", "Correct digis  [pct] vs. Found digis [pct], no noise; Correct; Found ", 110,
              -5., 105., 110, -5., 105.);
-  fhTimeOffset = new TH1F("fhTimeOffset", "tSeed - tMCEvent [ns]", 40, -10, 10);
+  fhTimeOffset = new TH1F("fhTimeOffset", "tSeed - tMCEvent [ns]", 20, -5, 5);
+  fhTimeOffset->SetCanExtend(TH1::kAllAxes);
 
   histFolder->Add(fhCorrectDigiRatio);
   histFolder->Add(fhCorrectDigiRatioNoNoise);
@@ -39,8 +40,8 @@ CbmSeedFinderQa::CbmSeedFinderQa() : fOutFolder("SeedFinderQA", "Seed finder QA"
   histFolder->Add(fhCorrectVsFoundNoNoise);
   histFolder->Add(fhTimeOffset);
 
-  fCanv = new CbmQaCanvas("cSummary", "", 3 * 400, 2 * 400);
-  fCanv->Divide2D(6);
+  fCanv = new CbmQaCanvas("cSummary", "", 4 * 400, 2 * 400);
+  fCanv->Divide2D(7);
   fOutFolder.Add(fCanv);
 }
 
@@ -56,6 +57,15 @@ CbmSeedFinderQa::~CbmSeedFinderQa()
   delete fCanv;
 }
 
+void CbmSeedFinderQa::Init()
+{
+  if (!FairRootManager::Instance() || !FairRootManager::Instance()->GetObject("MCEventList.")) {
+    LOG(error) << "No MC event list found";
+    return;
+  }
+  fEventList = (CbmMCEventList*) FairRootManager::Instance()->GetObject("MCEventList.");
+}
+
 void CbmSeedFinderQa::FillQaInfo(const int32_t WinStart, const int32_t WinEnd, const std::vector<CbmMatch>* vDigiMatch,
                                  const double seedTime)
 {
@@ -98,8 +108,11 @@ void CbmSeedFinderQa::FillQaInfo(const int32_t WinStart, const int32_t WinEnd, c
   {
     seedMatch.AddLink(1.0, 0, -1, 0);
     fvEventMatches.push_back(seedMatch);
-    fvCorrectDigiRatio.push_back(-1.0);
-    fvFoundDigiRatio.push_back(-1.0);
+    //should never show up in histograms
+    fvCorrectDigiRatio.push_back(std::numeric_limits<double>::quiet_NaN());
+    fvCorrectDigiRatioNoNoise.push_back(std::numeric_limits<double>::quiet_NaN());
+    fvFoundDigiRatio.push_back(std::numeric_limits<double>::quiet_NaN());
+    fvTimeOffset.push_back(std::numeric_limits<double>::quiet_NaN());
     return;
   }
   else {
@@ -131,18 +144,22 @@ void CbmSeedFinderQa::FillQaInfo(const int32_t WinStart, const int32_t WinEnd, c
   }
   const double foundDigiRatio = (double) correctDigiCount / matchedEventDigiCount;
   fvFoundDigiRatio.push_back(foundDigiRatio);
+
+  //seed time offset to MC
+  const CbmLink& eventLink = seedMatch.GetMatchedLink();
+  const double timeDiff    = seedTime - fEventList->GetEventTime(eventLink.GetEntry(), eventLink.GetFile());
+  fvTimeOffset.push_back(timeDiff);
 }
 
 void CbmSeedFinderQa::FillHistos()
 {
-  if (!FairRootManager::Instance() || !FairRootManager::Instance()->GetObject("MCEventList.")) {
-    LOG(error) << "No MC event list found";
-    return;
-  }
-  CbmMCEventList* eventList = (CbmMCEventList*) FairRootManager::Instance()->GetObject("MCEventList.");
-  //eventList->Print();
-
   for (uint32_t iEvent = 0; iEvent < fvEventMatches.size(); iEvent++) {
+
+    const CbmMatch* match = &(fvEventMatches.at(iEvent));
+    const CbmLink& link   = match->GetMatchedLink();
+    if (link.GetEntry() == -1) { continue; }
+
+    fhTimeOffset->Fill(fvTimeOffset.at(iEvent));
     const int32_t digiCount                 = fvFullDigiCount.at(iEvent);
     const int32_t noiseCount                = fvNoiseDigiCount.at(iEvent);
     const double noiseDigisPercent          = 100. * (double) noiseCount / digiCount;
@@ -155,14 +172,6 @@ void CbmSeedFinderQa::FillHistos()
     fhFoundDigiRatio->Fill(foundDigisPercent);
     fhCorrectVsFound->Fill(correctDigisPercent, foundDigisPercent);
     fhCorrectVsFoundNoNoise->Fill(correctDigisPercentNoNoise, foundDigisPercent);
-
-    const CbmMatch* match = &(fvEventMatches.at(iEvent));
-    const CbmLink& link   = match->GetMatchedLink();
-
-    if (link.GetEntry() != -1) {
-      const double timeDiff = fvSeedTimesFull.at(iEvent) - eventList->GetEventTime(link.GetEntry(), link.GetFile());
-      fhTimeOffset->Fill(timeDiff);
-    }
   }
 }
 
@@ -194,12 +203,6 @@ void CbmSeedFinderQa::OutputQa()
 
 void CbmSeedFinderQa::WriteHistos()
 {
-  //output histograms
-  if (!FairRootManager::Instance() || !FairRootManager::Instance()->GetSink()) {
-    LOG(error) << "No sink found";
-    return;
-  }
-
   fCanv->cd(1);
   fhCorrectDigiRatio->DrawCopy("colz", "");
 
@@ -218,6 +221,9 @@ void CbmSeedFinderQa::WriteHistos()
   fCanv->cd(6);
   fhCorrectVsFoundNoNoise->DrawCopy("colz", "");
 
+  fCanv->cd(7);
+  fhTimeOffset->DrawCopy("colz", "");
+
   FairSink* sink = FairRootManager::Instance()->GetSink();
   sink->WriteObject(&fOutFolder, nullptr);
 }
diff --git a/reco/eventbuilder/digis/CbmSeedFinderQa.h b/reco/eventbuilder/digis/CbmSeedFinderQa.h
index b2d1729a5b..02c5a223f0 100644
--- a/reco/eventbuilder/digis/CbmSeedFinderQa.h
+++ b/reco/eventbuilder/digis/CbmSeedFinderQa.h
@@ -27,6 +27,7 @@
 
 class TH1F;
 class TH2I;
+class CbmMCEventList;
 class CbmQaCanvas;
 
 class CbmSeedFinderQa {
@@ -55,6 +56,8 @@ public:
   std::vector<double> fvCorrectDigiRatioNoNoise;
   /** @brief Ratio of digis of matched events that were included in event seed. */
   std::vector<double> fvFoundDigiRatio;
+  /** @brief Difference between true event time and seed time. */
+  std::vector<double> fvTimeOffset;
 
   /** @brief Gather QA Information. 
   * @params WinStart Starting position of seed window.
@@ -71,6 +74,9 @@ public:
   /** @brief Finalize histograms and canvases and write to file. */
   void WriteHistos();
 
+  /** @brief Initialize communication with FairRootManager (needed for MC events). */
+  void Init();
+
 private:
   TFolder* histFolder = nullptr;  /// subfolder for histograms
   TFolder fOutFolder;             /// output folder with histos and canvases
@@ -82,10 +88,12 @@ private:
   TH1F* fhFoundDigiRatio          = nullptr;  /// digis found per event
   TH2I* fhCorrectVsFound          = nullptr;  /// correct digis per event vs found digis per event
   TH2I* fhCorrectVsFoundNoNoise   = nullptr;  /// correct digis per event vs found digis per event, disregarding noise
-  TH1F* fhTimeOffset              = nullptr;  /// difference between true event time and seed
+  TH1F* fhTimeOffset              = nullptr;  /// difference between true event time and seed time
 
   CbmQaCanvas* fCanv;  ///summary canvas
 
+  CbmMCEventList* fEventList = nullptr;  // to access MC truth
+
   /** @brief Fill output histograms with data from current timeslice. */
   void FillHistos();
 };
diff --git a/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.cxx b/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.cxx
index 8dd3e7dbf7..0e4f2f72fa 100644
--- a/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.cxx
+++ b/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.cxx
@@ -126,6 +126,11 @@ void CbmSeedFinderSlidingWindow::SetQa(bool doQA)
   }
 }
 
+void CbmSeedFinderSlidingWindow::Init()
+{
+  if (fQa) { fQa->Init(); }
+}
+
 void CbmSeedFinderSlidingWindow::OutputQa()
 {
   if (fQa) { fQa->OutputQa(); }
diff --git a/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.h b/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.h
index 06e55390e0..d8dce3b971 100644
--- a/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.h
+++ b/reco/eventbuilder/digis/CbmSeedFinderSlidingWindow.h
@@ -65,6 +65,9 @@ public:
   /** @brief Returns number of seed times currently stored in buffer. */
   size_t GetNofSeeds() { return fvSeedTimes->size(); }
 
+  /** @brief Initializes QA object if set. */
+  void Init();
+
 private:
   /** @brief Processes QA info. */
   CbmSeedFinderQa* fQa = nullptr;
diff --git a/reco/eventbuilder/digis/CbmTaskBuildRawEvents.cxx b/reco/eventbuilder/digis/CbmTaskBuildRawEvents.cxx
index 647f28e3d3..7f106fee54 100644
--- a/reco/eventbuilder/digis/CbmTaskBuildRawEvents.cxx
+++ b/reco/eventbuilder/digis/CbmTaskBuildRawEvents.cxx
@@ -170,6 +170,9 @@ InitStatus CbmTaskBuildRawEvents::Init()
     LOG(info) << "CbmTaskBuildRawEvents::Init(): Real time " << rtime << " s, CPU time " << ctime << " s";
   }
 
+  // Init seed finder
+  if (fSeedFinderSlidingWindow) { fSeedFinderSlidingWindow->Init(); }
+
   /// Call Algo Init method
   if (kTRUE == fpAlgo->InitAlgo()) return kSUCCESS;
   else
-- 
GitLab