Add mTrd timeshift correction container

During mCbm 2020 timeshifts of the Trd digi time were observed.
This commit adds a parameter container that allows a first order
correction of these shifts.

fles/mcbm2018/CMakeLists.txt

@@ -68,6 +68,7 @@ Set(SRCS
    parameter/CbmMcbm2018HodoPar.cxx
    parameter/CbmMcbm2018ContFact.cxx
    parameter/CbmMcbm2018PsdPar.cxx
+   parameter/CbmMcbm2020TrdTshiftPar.cxx
 
    unpacker/CbmMcbm2018RawConverterSdpb.cxx
    unpacker/CbmMcbm2018RawConverterGdpb.cxx

fles/mcbm2018/CbmFlibMcbm2018LinkDef.h

@@ -15,6 +15,7 @@
 #pragma link C++ class CbmMcbm2018HodoPar;
 #pragma link C++ class CbmMcbm2018ContFact;
 #pragma link C++ class CbmMcbm2018PsdPar;
+#pragma link C++ class CbmMcbm2020TrdTshiftPar + ;
 
 #pragma link C++ class CbmMcbm2018RawConverterSdpb;
 #pragma link C++ class CbmMcbm2018RawConverterGdpb;

fles/mcbm2018/parameter/CbmMcbm2018ContFact.cxx

@@ -14,6 +14,7 @@
 #include "CbmMcbm2018RichPar.h"
 #include "CbmMcbm2018StsPar.h"
 #include "CbmMcbm2018TofPar.h"
+#include "CbmMcbm2020TrdTshiftPar.h"
 
 #include "FairRuntimeDb.h"
 
@@ -72,6 +73,13 @@ void CbmMcbm2018ContFact::setAllContainers() {
                                           "TestDefaultContext");
   pPsd->addContext("TestNonDefaultContext");
   containers->Add(pPsd);
+
+  FairContainer* pTrd =
+    new FairContainer("CbmMcbm2020TrdTshiftPar",
+                      "TRD timeshift unpacker parameters mCbm 2020",
+                      "Trd mCbm 2020 timeshift correction");
+  pTrd->addContext("TestDefaultContext");
+  containers->Add(pTrd);
 }
 
 FairParSet* CbmMcbm2018ContFact::createContainer(FairContainer* c) {
@@ -105,6 +113,10 @@ FairParSet* CbmMcbm2018ContFact::createContainer(FairContainer* c) {
     p = new CbmMcbm2018PsdPar(
       c->getConcatName().Data(), c->GetTitle(), c->getContext());
   }
+  if (strcmp(name, "CbmMcbm2020TrdTshiftPar") == 0) {
+    p = new CbmMcbm2020TrdTshiftPar(
+      c->getConcatName().Data(), c->GetTitle(), c->getContext());
+  }
 
   return p;
 }

fles/mcbm2018/parameter/CbmMcbm2018ContFact.h

@@ -16,7 +16,7 @@ public:
   CbmMcbm2018ContFact();
   ~CbmMcbm2018ContFact() {}
   FairParSet* createContainer(FairContainer*);
-  ClassDef(CbmMcbm2018ContFact, 0)  // Factory for all TRD parameter containers
+  ClassDef(CbmMcbm2018ContFact, 0)  // Factory for all mcbm parameter containers
 };
 
 #endif /* !CBMMCBM2018CONTFACT_H */

fles/mcbm2018/parameter/CbmMcbm2020TrdTshiftPar.cxx

+/*
+ * File: /CbmMcbm2020TrdTshiftPar.cxx
+ * Created Date: Thursday December 17th 2020
+ * Author: Pascal Raisig -- praisig@ikf.uni-frankfurt.de
+ * -----
+ * Last Modified: Tuesday December 22nd 2020 16:35:01
+ * Modified By: Pascal Raisig
+ * -----
+ * Purpose: This container is only to be used for Trd mCbm2020 data. It contains the
+ * corrections for the observed timeshifts within the 2020 data. 
+ * -----
+ */
+
+
+#include "CbmMcbm2020TrdTshiftPar.h"
+
+// FairRoot
+#include "FairLogger.h"
+#include "FairParamList.h"
+
+// ROOT
+#include "THnSparse.h"
+
+// C/C++
+#include <iomanip>
+#include <iostream>
+
+// #define CTAVM CbmTrdAnalysisVarManager // See comments in GetNevents() et al.
+// #define CTAH CbmTrdAnalysisHisto // See comments in GetNevents() et al.
+
+CbmMcbm2020TrdTshiftPar::CbmMcbm2020TrdTshiftPar(const char* name,
+                                                 const char* title,
+                                                 const char* context)
+  : FairParGenericSet(name, title, context), fNtimeslices(1) {
+  detName = "TRD";
+}
+
+CbmMcbm2020TrdTshiftPar::~CbmMcbm2020TrdTshiftPar() { clear(); }
+
+// ----- Public method clear ----
+void CbmMcbm2020TrdTshiftPar::clear()
+
+{
+  status = kFALSE;
+  resetInputVersions();
+}
+
+// ----- Public method printparams ----
+void CbmMcbm2020TrdTshiftPar::printparams()
+
+{
+  std::cout << "CbmMcbm2020TrdTshiftPar::printparams() " << &fTimeshifts
+            << std::endl;
+  auto ntimeshifts = (fTimeshifts.GetSize() / (fgNchannels + 1));
+  std::cout << "ParSet has " << ntimeshifts << " timeshift changes stored"
+            << std::endl;
+
+  for (size_t ishiftblock = 0; ishiftblock < ntimeshifts; ishiftblock++) {
+    std::cout << "Shiftblock of event "
+              << fTimeshifts[ishiftblock * (fgNchannels + 1)] << std::endl;
+    for (size_t ichannel = 0; ichannel < fgNchannels; ichannel++) {
+      if (ichannel % 6 == 0) std::cout << std::endl;
+      std::cout << " Ch " << ichannel << " shift "
+                << fTimeshifts[(ishiftblock * (fgNchannels + 1)) + ichannel];
+    }
+  }
+}
+
+// ---- putParams ----
+void CbmMcbm2020TrdTshiftPar::putParams(FairParamList* l) {
+  if (!l) return;
+  l->add(pararraynames[0].data(), fNtimeslices);
+  l->add(pararraynames[1].data(), fTimeshifts);
+}
+
+// ---- getParams ----
+Bool_t CbmMcbm2020TrdTshiftPar::getParams(FairParamList* l) {
+  if (!l) return kFALSE;
+
+  if (!l->fill(pararraynames[0].data(), &fNtimeslices)) return kFALSE;
+  if (!l->fill(pararraynames[1].data(), &fTimeshifts)) return kFALSE;
+
+  // Create a map from the list imported from the par file
+  size_t nthShift = 0;
+
+  size_t nentries   = fTimeshifts.GetSize();
+  size_t ichannel   = 0;
+  size_t itimeslice = 0;
+  for (size_t ientry = 0; ientry < nentries; ientry++) {
+    // We have a chain of one eventId value + fgNchannels shift values
+    nthShift = ientry / (fgNchannels + 1);
+
+    // Look for the first value in the chain, which is the eventId
+    if ((ientry - nthShift * (fgNchannels + 1)) == 0) {
+      // Before starting the extraction of the next chain, emplace the previous chain to the map
+      if (ientry != 0) {
+        auto tspair = std::pair<size_t, std::vector<Int_t>>(
+          itimeslice, fvecCurrentTimeshifts);
+        fmapTimeshifts.emplace(tspair);
+      }
+      itimeslice = fTimeshifts[ientry];
+      fvecCurrentTimeshifts.clear();
+      fvecCurrentTimeshifts.resize(fgNchannels);
+
+    } else {
+      ichannel = ientry - 1 - nthShift * (fgNchannels + 1);
+      fvecCurrentTimeshifts[ichannel] = fTimeshifts[ientry];
+    }
+  }
+  // Now we have to fill the blank spots in the map, since, we do not now if the timeslices are accessed in a completely ordered manor
+  for (itimeslice = 0; itimeslice < std::abs(fNtimeslices[0]); itimeslice++) {
+    auto itspair = fmapTimeshifts.find(itimeslice);
+
+    if (itspair != fmapTimeshifts.end()) {
+      continue;
+    } else {
+      // Get previous timeshift vector to add it also to the current pair
+      itspair--;
+
+      auto newtspair =
+        std::pair<size_t, std::vector<Int_t>>(itimeslice, itspair->second);
+      fmapTimeshifts.emplace(newtspair);
+    }
+  }
+  return kTRUE;
+}
+
+// ---- GetTimeshifts ----
+bool CbmMcbm2020TrdTshiftPar::GetTimeshifts(
+  std::shared_ptr<TH2> timeshiftHisto) {
+
+  ///< Extract the timeshift values from a histo containing the information and write them to fTimeshifts.
+
+  fNtimeslices[0] = timeshiftHisto->GetNbinsX();
+  size_t nshifts  = 1;
+  fvecCurrentTimeshifts.clear();
+  fvecCurrentTimeshifts.resize(0, 0);
+  fvecCurrentTimeshifts.resize(fgNchannels, 0);
+  Int_t tshift   = 0;
+  size_t tsidx   = 0;
+  size_t ientry  = 0;
+  size_t eventId = 0;
+
+  // We write at least one shift value to fTimeshifts. Thus, we resize it to one valueset
+  size_t nentries = nshifts * (fgNchannels + 1);
+  fTimeshifts.Set(nentries);
+
+  bool didChange = true;
+  for (size_t ievent = 1; ievent < std::abs(fNtimeslices[0]); ievent++) {
+    tsidx = timeshiftHisto->GetXaxis()->GetBinLowEdge(ievent);
+    for (size_t ichannel = 0; ichannel < fgNchannels; ichannel++) {
+      tshift = (Int_t) timeshiftHisto->GetBinContent(ievent, ichannel);
+      if (tshift != fvecCurrentTimeshifts[ichannel]) {
+        didChange                       = true;
+        fvecCurrentTimeshifts[ichannel] = tshift;
+      }
+    }
+    if (didChange) {
+      // First resize the fTimeshifts array
+      nshifts++;
+      nentries = nshifts * (fgNchannels + 1);
+      fTimeshifts.Set(nentries);
+
+      // Then write the eventId correlated to the tsIdx to the array
+      // For mCbm2020 a timeslice contains 10 µslices, the tsidx jumps accordingly by 10, e.g. tsidx0 = 7 -> tsidx = 17
+      eventId             = tsidx / 10;
+      fTimeshifts[ientry] = eventId;
+      ientry++;
+      // Followed by the channelwise shift values
+      for (size_t ichannel = 0; ichannel < fgNchannels; ichannel++) {
+        tshift = fvecCurrentTimeshifts[ichannel];
+        // The above represents the timeshift in multiples of a µSlice length (for mCbm 2020 1280 ns) since here we have time to do the calculation of the actual timeshift, we do it here and not in the unpacker algo.
+        tshift *= fgMicroSliceLength;
+        fTimeshifts.AddAt(tshift, ientry);
+        ientry++;
+      }
+    }
+    // reset didChange value
+    didChange = false;
+  }
+  return fTimeshifts.GetSize() > 0;
+}
+
+// ---- GetTimeshiftsVec ----
+std::vector<Int_t> CbmMcbm2020TrdTshiftPar::GetTimeshiftsVec(size_t tsidx) {
+  auto pair = fmapTimeshifts.find(tsidx);
+  return pair->second;
+}
+
+// ---- GetNEvents ----
+double CbmMcbm2020TrdTshiftPar::GetNEvents(std::shared_ptr<TFile> mcbmanafile) {
+  ///< Extract the number of events from a mcbmana task output file, which has to contain the required histogram.
+
+  TH1* histo = nullptr;
+
+  // We have to do some hardcoding here, since the TrdAnalysisFramework is not yet in the common master. The handling with TAF in place will be given commented out
+
+  //First get the required histogram from the mcbmana file
+  // std::vector<CTAVM::eVars> varvec = {CTAVM::eVars::kRunId, CTAVM::eVars::kEnd, CTAVM::eVars::kEnd};
+  // auto weight                      = CTAVM::eVars::kNEvents;
+  // auto fillstation = CTAH::eFillStation::kRunInfo;
+  // auto htype       = (int) CTAH::eHistoType::kGlobal;
+  // histo = CTAH::GetHistoFromFile(varvec, fillstation, htype, mcbmanafile, weight);
+
+  std::string hpath = "FillStation-RunInfo/FullTrd/RunId_wNEvents-RunInfo";
+  histo             = (TH1*) mcbmanafile->Get(hpath.data());
+  if (!histo) {
+    LOG(fatal) << " CbmMcbm2020TrdTshiftPar::GetNEvents " << hpath.data()
+               << " not found in the file" << std::endl;
+  }
+  double nevents = histo->GetBinContent(histo->GetMaximumBin());
+  return nevents;
+}
+
+// ---- GetCalibHisto ----
+std::shared_ptr<TH3>
+CbmMcbm2020TrdTshiftPar::GetCalibHisto(std::shared_ptr<TFile> mcbmanafile) {
+  ///< Extract the required base histogram from a mcbmana task output file.
+  THnSparse* hsparse = nullptr;
+
+  // std::vector<CTAVM::eVars> varvec = {CTAVM::eVars::kTsSourceTsIndex, CTAVM::eVars::kDigiTrdChannel, CTAVM::eVars::kDigiDtCorrSlice};
+  // auto fillstation = CTAH::eFillStation::kTrdT0Digi;
+  // auto htype       = (int) CTAH::eHistoType::kGlobal;
+  // histo =
+  //   (THnSparseD*) CTAH::GetHistoFromFile(varvec, fillstation, htype, mcbmanafile);
+
+  std::string hpath =
+    "FillStation-TrdT0Digi/FullTrd/"
+    "TsSourceTsIndex_DigiTrdChannel_DigiDtCorrSlice-TrdT0Digi";
+  hsparse = (THnSparse*) mcbmanafile->Get(hpath.data());
+
+  if (!hsparse) {
+    LOG(fatal) << " CbmMcbm2020TrdTshiftPar::GetCalibHisto " << hpath.data()
+               << " not found in the file" << std::endl;
+  }
+
+  auto nevents = GetNEvents(mcbmanafile);
+
+  // auto tsaxis = CTAH::GetVarAxis(hsparse, CTAVM::eVars::kTsSourceTsIndex);
+  auto tsaxis =
+    hsparse->GetAxis(0);  // For now we know that the TsIndex is on the X-Axis
+  tsaxis->SetRangeUser(0.0, (double) nevents);
+  auto temphisto = hsparse->Projection(0, 1, 2);
+  auto histo     = std::make_shared<TH3D>(*temphisto);
+  delete temphisto;
+  histo->SetName("calibbasehisto");
+  histo->SetTitle("calibbasehisto");
+  // We need the uniqueIDs for the correct handling within TAF
+  histo->GetXaxis()->SetUniqueID(hsparse->GetAxis(0)->GetUniqueID());
+  histo->GetYaxis()->SetUniqueID(hsparse->GetAxis(1)->GetUniqueID());
+  histo->GetZaxis()->SetUniqueID(hsparse->GetAxis(2)->GetUniqueID());
+
+  return histo;
+}
+
+// ---- GetCalibHisto ----
+std::shared_ptr<TH2>
+CbmMcbm2020TrdTshiftPar::GetCalibHisto(std::shared_ptr<TH3> calibbasehisto) {
+  ///< Extract the timeshiftHisto from the calibbase histogram. The calibbase histogram is a TH3* with the tsIdx, the module channels and the timeshifts on the axes.
+
+  // Get the x-axis definitions
+  size_t nevents    = calibbasehisto->GetNbinsX();
+  size_t firstTsIdx = calibbasehisto->GetXaxis()->GetBinLowEdge(
+    calibbasehisto->GetXaxis()->GetFirst());
+  size_t lastTsIdx = calibbasehisto->GetXaxis()->GetBinUpEdge(
+    calibbasehisto->GetXaxis()->GetLast());
+
+  // Get the y-axis definitions
+  size_t nchannels    = calibbasehisto->GetNbinsY();
+  size_t firstChannel = calibbasehisto->GetYaxis()->GetBinLowEdge(
+    calibbasehisto->GetYaxis()->GetFirst());
+  size_t lastChannel = calibbasehisto->GetYaxis()->GetBinUpEdge(
+    calibbasehisto->GetYaxis()->GetLast());
+
+  std::shared_ptr<TH2I> calibhisto = std::make_shared<TH2I>("calibhisto",
+                                                            "calibhisto",
+                                                            nevents,
+                                                            firstTsIdx,
+                                                            lastTsIdx,
+                                                            nchannels,
+                                                            firstChannel,
+                                                            lastChannel);
+
+
+  for (size_t itsidx = 1; itsidx < nevents; itsidx++) {
+    for (size_t ichannel = 1; ichannel < nchannels; ichannel++) {
+      auto dominantshift =
+        GetDominantShift(calibbasehisto, itsidx, ichannel) < 255
+          ? GetDominantShift(calibbasehisto, itsidx, ichannel)
+          : calibhisto->GetBinContent(itsidx - 1, ichannel);
+      if (itsidx - 1 == 0) dominantshift = 0;
+      calibhisto->SetBinContent(itsidx, ichannel, dominantshift);
+    }
+  }
+  return calibhisto;
+}
+
+// ---- GetDominantShift ----
+Int_t CbmMcbm2020TrdTshiftPar::GetDominantShift(
+  std::shared_ptr<TH3> calibbasehisto,
+  size_t itsidx,
+  size_t ichannel) {
+  auto hdomshift =
+    calibbasehisto->ProjectionZ("domshift", itsidx, itsidx, ichannel, ichannel);
+
+  // Scale histo to one
+  hdomshift->Scale(1. / hdomshift->Integral());
+
+  auto domshift  = hdomshift->GetBinCenter(hdomshift->GetMaximumBin());
+  auto dominance = hdomshift->GetMaximum();
+
+  delete hdomshift;
+
+  // Only add dominant slices with at least 50% of the entries in the dominant slice
+
+  if (dominance > 0.5)
+    // if (dominance > 0.25)
+    return domshift;
+  else
+    return 255;
+}
+
+// ---- FillTimeshiftArray ----
+bool CbmMcbm2020TrdTshiftPar::FillTimeshiftArray(
+  std::shared_ptr<TFile> mcbmanafile) {
+  ///< Extract the timeshift values from a TAF output file that contains the required histograms and write them to fTimeshifts.
+
+  auto calibbasehisto = GetCalibHisto(mcbmanafile);
+  auto calibhisto     = GetCalibHisto(calibbasehisto);
+  return GetTimeshifts(calibhisto);
+}
+
+ClassImp(CbmMcbm2020TrdTshiftPar)

fles/mcbm2018/parameter/CbmMcbm2020TrdTshiftPar.h

+#ifndef CbmMcbm2020TrdTshiftPar_H
+#define CbmMcbm2020TrdTshiftPar_H
+
+#include "FairParGenericSet.h"  // mother class
+
+// STD
+#include <map>
+
+// ROOT
+#include <TArrayD.h>
+#include <TArrayI.h>
+#include <TFile.h>
+#include <TH2.h>
+#include <TH3.h>
+
+class CbmMcbm2020TrdTshiftPar : public FairParGenericSet {
+public:
+  CbmMcbm2020TrdTshiftPar(
+    const char* name    = "CbmMcbm2020TrdTshiftPar",
+    const char* title   = "TRD timeshift unpacker parameters mCbm 2020",
+    const char* context = "Default");
+
+  virtual ~CbmMcbm2020TrdTshiftPar();
+
+  virtual void printparams();
+
+  /** Reset all parameters **/
+  virtual void clear();
+
+  void putParams(FairParamList*);
+  Bool_t getParams(FairParamList*);
+
+  bool GetTimeshifts(std::shared_ptr<TH2> timeshiftHisto);
+  ///< Extract the timeshift values from a histo containing the information and write them to fTimeshifts.
+
+  double GetNEvents(std::shared_ptr<TFile> mcbmanafile);
+  ///< Extract the number of events from a mcbmana task output file, which has to contain the required histogram.
+
+
+  std::vector<Int_t> GetTimeshiftsVec(size_t tsidx);
+  ///< Return the timeshift vector for a given Timeslice Index. Works only if getParams() was run before
+
+  std::map<size_t, std::vector<Int_t>>* GetTimeshiftsMap() {
+    return &fmapTimeshifts;
+  }
+  ///< Return the timeshift map.
+
+
+  std::shared_ptr<TH3> GetCalibHisto(std::shared_ptr<TFile> mcbmanafile);
+  ///< Extract the required base histogram from a mcbmana task output file.
+
+  std::shared_ptr<TH2> GetCalibHisto(std::shared_ptr<TH3> calibbasehisto);
+  ///< Extract the timeshiftHisto from the calibbase histogram. The calibbase histogram is a TH3* with the tsIdx, the module channels and the timeshifts on the axes.
+
+  Int_t GetDominantShift(std::shared_ptr<TH3> calibbasehisto,
+                         size_t itsidx,
+                         size_t ichannel);
+  ///< Extract the dominant shift value from the calibbase histogram. For a give tsIdx and channel.
+
+  bool FillTimeshiftArray(std::shared_ptr<TFile> mcbmanafile);
+  ///< Extract the timeshift values from a TAF output file that contains the required histograms and write them to fTimeshifts.
+
+private:
+  const std::vector<std::string> pararraynames = {"nTimeslices",
+                                                  "TrdTimeshifts"};
+  ///< Names of the parameter arrays
+
+  const size_t fgNchannels = 768;
+  ///< Number of channels on the single Trd module used during mCbm2020
+  // In principle this could be retrieved from the standard parameters. But since we hopefully need these ugly parameters only for mCbm2020 it is hardcoded, because a handling via the parameters would require full linking of the standard Trd parameter classes
+
+  const UInt_t fgMicroSliceLength = 1.280e6;
+  ///< Length of a single microslice during mCbm 2020 data taking
+
+  TArrayI fTimeshifts = {};
+  ///< Array contains the timeshifts. Everytime a timeshift appears the tsIdx from the TS MetaData is followed by the correct shift value for each single channel
+
+  TArrayI fNtimeslices;
+  ///< Number of timeslices in the given run. This is required to fill a complete fmapTimeshifts, when reading the compressed information from the ascii parameter file
+
+  std::map<size_t, std::vector<Int_t>> fmapTimeshifts = {};  //!
+  ///< Keys of the map are the tsIdx when a shift value for any channel changes, the vector contains the shift values for the given period for all channels. It is represented by fvecCurrentTimeshifts
+
+  std::vector<Int_t> fvecCurrentTimeshifts = {};  //!
+  ///< Vector containing the "current" timeshift values for all Trd channles of the mCbm2020 module
+
+public:
+  ClassDef(CbmMcbm2020TrdTshiftPar, 2);
+};
+
+#endif  // CbmMcbm2020TrdTshiftPar_H
\ No newline at end of file