CbmRecoQaTask.cxx

Go to the documentation of this file.

/* Copyright (C) 2024-2026 Hulubei National Institute of Physics and Nuclear Engineering - Horia Hulubei, Bucharest
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Alexandru Bercuci [committer], Omveer Singh */
 
//CBM headers
#include "CbmRecoQaTask.h"
 
#include "../../core/detectors/rich/utils/CbmRichUtil.h"
#include "CbmCluster.h"
#include "CbmDefs.h"
#include "CbmEvent.h"
#include "CbmFsdHit.h"
#include "CbmGlobalTrack.h"
#include "CbmMCDataArray.h"
#include "CbmMCDataManager.h"
#include "CbmMvdHit.h"
#include "CbmPsdHit.h"
#include "CbmRichHit.h"
#include "CbmSetup.h"
#include "CbmStsAddress.h"
#include "CbmStsCluster.h"
#include "CbmStsHit.h"
#include "CbmStsSetup.h"
#include "CbmTimeSlice.h"
#include "CbmTofAddress.h"
#include "CbmTofHit.h"
#include "CbmTrdAddress.h"
#include "CbmTrdCluster.h"
#include "CbmTrdHit.h"
// FAIR headers
#include "FairMCPoint.h"
 
#include <FairRootManager.h>
#include <FairSink.h>
// ROOT headers
#include <TClonesArray.h>
#include <TGeoManager.h>
#include <TGeoNode.h>
#include <TH2D.h>
#include <TH3D.h>
 
#include <iterator>
#include <regex>
 
using namespace std;
using namespace cbm::algo;
std::bitset<CbmRecoQaTask::eRecoConfig::kRecoQaNConfigs> CbmRecoQaTask::fuRecoConfig = {};
 
//_____________________________________________________________________
CbmRecoQaTask::CbmRecoQaTask() : FairTask("RecoQA", 0)
{
  fuRecoConfig.set(eRecoConfig::kUseSts);
  fuRecoConfig.set(eRecoConfig::kUseTrd);
  fuRecoConfig.set(eRecoConfig::kUseTof);
}
 
//_____________________________________________________________________
CbmRecoQaTask::Detector* CbmRecoQaTask::AddDetector(ECbmModuleId id)
{
  /* Interface to the user. Check that the det defined in the outside world
   * fulfills the quality criteria.
   */
  if (GetDetector(id)) {
    LOG(warn) << GetName() << "::AddDetector(" << id << ")."
              << " already registered. Using "
              << "\"CbmRecoQaTask::GetDetector(ECbmModuleId).\"";
    return GetDetector(id);
  }
  switch (id) {
    case ECbmModuleId::kMvd:
    case ECbmModuleId::kSts:
    case ECbmModuleId::kMuch:
    case ECbmModuleId::kRich:
    case ECbmModuleId::kTrd:
    case ECbmModuleId::kTrd2d:
    case ECbmModuleId::kTof:
    case ECbmModuleId::kFsd:
    case ECbmModuleId::kPsd:
      LOG(debug) << GetName() << "::AddDetector(" << id << ").";
      fDetQa.emplace(id, id);
      break;
    default: LOG(warn) << GetName() << "::AddDetector : unsupported det " << id; return nullptr;
  }
 
  return &fDetQa[id];
}
 
//_____________________________________________________________________
void CbmRecoQaTask::Activate(ECbmModuleId id, bool set)
{
  // translate to local parameterization
  eRecoConfig feature = eRecoConfig::kRecoQaNConfigs;
  switch (id) {
    case ECbmModuleId::kSts: feature = eRecoConfig::kUseSts; break;
    case ECbmModuleId::kTrd: feature = eRecoConfig::kUseTrd; break;
    case ECbmModuleId::kTof: feature = eRecoConfig::kUseTof; break;
    default: break;
  }
  if (feature == eRecoConfig::kRecoQaNConfigs) {
    LOG(warn) << GetName() << "::Activate : unsupported operation on det " << id;
    return;
  }
  if (set) {
    LOG(info) << GetName() << "::Activate " << id << " for track fit.";
    fuRecoConfig.set(feature);
  }
  else {
    LOG(info) << GetName() << "::Deactivate " << id << " for track fit.";
    fuRecoConfig.reset(feature);
  }
}
 
//_____________________________________________________________________
CbmRecoQaTask::Detector* CbmRecoQaTask::GetDetector(ECbmModuleId id)
{
  if (fDetQa.find(id) == fDetQa.end()) return nullptr;
  return &fDetQa[id];
}
 
//_____________________________________________________________________
template<>
std::map<ECbmModuleId, TClonesArray*> CbmRecoQaTask::MapDataList<CbmCluster>(std::string& s) const
{
  s = "cluster";
  return fClusters;
}
template<>
std::map<ECbmModuleId, TClonesArray*> CbmRecoQaTask::MapDataList<CbmHit>(std::string& s) const
{
  s = "hit";
  return fHits;
}
template<>
std::map<ECbmModuleId, TClonesArray*> CbmRecoQaTask::MapDataList<CbmTrack>(std::string& s) const
{
  s = "track";
  return fTracks;
}
template<>
std::map<ECbmModuleId, TClonesArray*> CbmRecoQaTask::MapDataList<CbmMatch>(std::string& s) const
{
  s = "match";
  return fHitMatch;
}
template<class Data>
const Data* CbmRecoQaTask::GetData(ECbmModuleId did, int id) const
{
  if (id < 0) {
    LOG(warning) << GetName() << " negative data idx for " << did;
    return nullptr;
  }
  std::string ss;
  auto datas = MapDataList<Data>(ss);
  if (datas.find(did) == datas.end()) {
    LOG(debug1) << GetName() << " missing " << ss << "s for " << did;
    return nullptr;
  }
  auto dd = (const Data*) (datas.at(did))->At(id);
  if (!dd) {
    LOG(debug1) << GetName() << " missing " << ss << "_" << id << " for " << did;
    return nullptr;
  }
  return dd;
}
 
//_____________________________________________________________________
InitStatus CbmRecoQaTask::Init()
{
  // INIT TRACKING INTERFACE
  fFitter.Init();
  fFitter.SetIgnorePoorCoordinates();
 
  // INIT DATA
  FairRootManager* ioman = FairRootManager::Instance();
  if (!ioman) {
    LOG(fatal) << GetName() << "::Init :: RootManager not instantiated!";
    return kERROR;
  }
  // read global tracks
  fGTracks = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrack"));
  if (!fGTracks)
    LOG(warn) << GetName() << "::Init: Global track array not found!";
  else
    fuRecoConfig.set(kRecoTracks);
  // read reco events
  fEvents = static_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
  if (!fEvents)
    LOG(warn) << GetName() << "::Init: No event found. Some results will be missing.";
  else
    fuRecoConfig.set(kRecoEvents);
  // read MC particles using the MC datat manager
  if (fuRecoConfig[kUseMC]) {
    cbm_mc_manager = dynamic_cast<CbmMCDataManager*>(ioman->GetObject("MCDataManager"));
    if (!cbm_mc_manager) {
      LOG(warn) << GetName() << "::Init: MC data manager not available even though asked by user !";
      fuRecoConfig.reset(kUseMC);
    }
    else {
      fTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrackMatch"));
      if (!fTrackMatches) {
        LOG(warn) << GetName() << "::Init: MC info for Global track not available !";
      }
    }
  }
  // initialize setup. For special cases use specific configurations
  // TODO should be moved to configuration
  switch (fSetupClass) {
    case kMcbm22:
      LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2022\".";
      InitMcbm22();
      break;
    case kMcbm24:
      LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2024\". ";
      InitMcbm24();
      break;
    case kMcbm25:
      LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2025\". ";
      InitMcbm25();
      break;
    case kDefault:
      LOG(info) << GetName() << "::Init: Setup config for \"Default\". ";
      InitDefault();
      break;
  }
  // init MC points, hits and clusters for all systems registered in the setup
  fOutFolder.Clear();
  for (auto& detp : fDetQa) {  // SYSTEM LOOP
    auto& det = detp.second;
    // read Hit -> MC linking
    if (cbm_mc_manager) {
      fPoints[det.id] = cbm_mc_manager->InitBranch(det.point.name.data());
      if (!fPoints[det.id]) LOG(warn) << GetName() << "::Init: MC Point array for " << det.id << " not found!";
      fHitMatch[det.id] = static_cast<TClonesArray*>(ioman->GetObject((det.hit.name + "Match").c_str()));
 
      if (!fHitMatch[det.id]) LOG(warn) << GetName() << "::Init: Hit Match array for " << det.id << " not found!";
    }
    // init reco hits
    fHits[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.hit.name.data()));
    if (!fHits[det.id]) {
      LOG(warn) << GetName() << "::Init: Hit array for " << cbm::util::ToString(det.id) << " not found!";
      continue;
    }
    switch (det.id) {
      case ECbmModuleId::kSts: fuRecoConfig.set(kStsHits); break;
      case ECbmModuleId::kTrd: fuRecoConfig.set(kTrdHits); break;
      case ECbmModuleId::kTof: fuRecoConfig.set(kTofHits); break;
      case ECbmModuleId::kRich: fuRecoConfig.set(kRichHits); break;
      case ECbmModuleId::kMuch: fuRecoConfig.set(kMuchHits); break;
      default: break;
    }
 
    // init reco clusters
    if (det.cl.id != ECbmDataType::kUnknown) {
      fClusters[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.cl.name.data()));
      if (!fClusters[det.id]) {
        LOG(warn) << GetName() << "::Init: Cluster array for " << det.id << " not found!";
        continue;
      }
    }
    // init reco track segments for each detector registered in the setup
    if (fGTracks && det.trk.id != ECbmDataType::kUnknown) {
      fTracks[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.trk.name.data()));
      if (!fTracks[det.id]) LOG(warn) << GetName() << "::Init: Track array for " << det.id << " not found!";
    }
    if (!det.Init(&fOutFolder, (bool) fuRecoConfig[kUseMC])) continue;
  }
 
  // INIT PROJECTION PLANES FOR TRACKS
  TDirectoryFile* trgDir = (fuRecoConfig[kRecoEvents] || fuRecoConfig[kUseMC] || GetNviews(eViewType::kTrkProj))
                             ? (TDirectoryFile*) fOutFolder.mkdir("TRG", "Target tomography with CA")
                             : nullptr;
  // register track projection on z planes
  if (fPrjPlanes.size() && trgDir) {
    fViews.emplace("prj", View("TrkProj", "", {}));
    fViews["prj"].fType = eViewType::kTrkProj;
    int i(0);
    for (auto plane : fPrjPlanes) {  // add projection planes according to user request
      if (i >= kNtrkProjections) {
        LOG(warn) << GetName() << "::Init: Only " << kNtrkProjections << " are supported. Skipping the rest.";
        fPrjPlanes.erase(fPrjPlanes.begin() + int(kNtrkProjections), fPrjPlanes.end());
        break;
      }
 
      fViews["prj"].AddProjection(CbmRecoQaTask::eProjectionType(int(CbmRecoQaTask::eProjectionType::kXYt0) + i),
                                  plane.Z());
      i++;
    }
    fViews["prj"].Init("");
    fViews["prj"].Register(trgDir);
  }
  // register also projection plans for track filtering if they are not already defined
  for (auto trkCut : fFilterTrk) {
    // special cut for selection on track composition. See FilterTrack(const CbmGlobalTrack*)
    if (trkCut.fDet != ECbmModuleId::kRef) continue;
    // read z coordinates of projection planes
    for (auto rcut : trkCut.GetProjCuts()) {
      // check if projection not already on the list
      bool projFound = false;
      for (auto plane : fPrjPlanes)
        if (abs(plane.Z() - rcut[0]) < 1.e-3) {
          projFound = true;
          break;
        }
      if (projFound)
        LOG(debug1) << GetName() << "::Init: Trk projection cut @ z=" << rcut[0]
                    << " already registered as projection plane";
      else
        fPrjPlanes.emplace_back(0., 0., rcut[0]);
    }
  }
 
  // INIT PRIMARY VERTEX QA
  if (fuRecoConfig[kRecoEvents] && trgDir) {
    fViews.emplace("vx", View("Vx", "", {}));
    fViews["vx"].fType = eViewType::kPV;
    fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVmult);
    fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVxy);
    fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVxz);
    fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVyz);
    fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVxsx);
    fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVysy);
    fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVzsz);
    fViews["vx"].Init("Prim");
    fViews["vx"].Register(trgDir);
  }
 
  // DUMP FILTERS SETTINGS
  for (auto evSelector : fFilterEv)
    LOG(info) << GetName() << evSelector.ToString();
  for (auto trkSelector : fFilterTrk)
    LOG(info) << GetName() << trkSelector.ToString();
  // add verbosity on active detectors
  LOG(info) << GetName() << "::Init: Using for tracking: " << (IsActive(eRecoConfig::kUseSts) ? "STS " : "")
            << (IsActive(eRecoConfig::kUseTrd) ? "TRD " : "") << (IsActive(eRecoConfig::kUseTof) ? "ToF " : "");
  return kSUCCESS;
}
 
//_________________________________________________________________
template<class Hit>
bool CbmRecoQaTask::View::HasAddress(const CbmHit* h, double&, double&, double&, double&) const
{
  LOG(error) << "Unprocessed hit in view " << name;
  cout << h->ToString();
 
  return false;
}
 
//_________________________________________________________________
template<ECbmModuleId d>
int CbmRecoQaTask::GetDebugInfo(const CbmHit*)
{
  return 0xab00;
}
template<>
int CbmRecoQaTask::GetDebugInfo<ECbmModuleId::kTrd>(const CbmHit* h)
{
 
  //auto htrd = dynamic_cast<const CbmTrdHit*>(h);
  auto ctrd = (const CbmTrdCluster*) GetData<CbmCluster>(ECbmModuleId::kTrd, h->GetRefId());
  if (!ctrd) return 0xab00;
  return ctrd->GetRow();
}
 
//_________________________________________________________________
template<>
size_t CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kSts>(const CbmHit* h, std::vector<int>& sel)
{
  int32_t a = h->GetAddress();
  sel.emplace_back(CbmStsAddress::GetElementId(a, EStsElementLevel::kStsUnit));
  sel.emplace_back(CbmStsAddress::GetElementId(a, EStsElementLevel::kStsLadder));
  sel.emplace_back(CbmStsAddress::GetElementId(a, EStsElementLevel::kStsModule));
  return sel.size();
}
 
//_________________________________________________________________
template<>
size_t CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kTrd>(const CbmHit* h, std::vector<int>& sel)
{
  int32_t a = h->GetAddress();
  sel.emplace_back(CbmTrdAddress::GetLayerId(CbmTrdAddress::GetModuleAddress(a)));
  sel.emplace_back(CbmTrdAddress::GetModuleId(a));
  return sel.size();
}
 
//_________________________________________________________________
template<>
size_t CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kTof>(const CbmHit* h, std::vector<int>& sel)
{
  int32_t a = h->GetAddress();
  sel.emplace_back(CbmTofAddress::GetSmId(a));
  sel.emplace_back(CbmTofAddress::GetSmType(a));
  sel.emplace_back(CbmTofAddress::GetRpcId(a));
  return sel.size();
}
 
//_________________________________________________________________
template<>
size_t CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kRich>(const CbmHit* h, std::vector<int>& sel)
{
  int32_t a   = h->GetAddress();
  int16_t uId = CbmRichUtil::GetDirichId(a);
  // TODO implement at RichUtil level
  sel.emplace_back((uId >> 8) & 0xF);
  return sel.size();
}
 
//_________________________________________________________________
template<ECbmModuleId d>
bool CbmRecoQaTask::TopoFilter::Accept(const CbmHit* h, const CbmCluster*, const CbmCluster*) const
{
 
  // check if there is any condition defined for current detector:
  if (!fReject.size() && !fAccept.size()) return true;
 
  // check individually the conditions
  std::vector<int> sel;
  auto nsel = CbmRecoQaTask::GetHitDetectorId<d>(h, sel);
 
  // Lambda function to filter specific location parameters
  auto filter = [&](vector<int> def) -> size_t {
    size_t accept(nsel);
    // check cut in topology [unit, ladder, module]
    for (size_t isel(0); isel < nsel; isel++) {
      if (def[isel] < 0 || def[isel] == sel[isel])
        accept--;
      else
        return accept;
    }
    return accept;
  };
  // check reject conditions
  for (auto def : fReject) {
    if (def.size() < nsel) {
      LOG(error) << "TopoFilter::Accept: Det[" << d << "] unsupported definition for [reject] cut of size["
                 << def.size() << "] expect " << nsel;
      break;
    }
    if (!filter(def)) return false;
  }
  // if not rejected check accepting conditions
  for (auto def : fAccept) {
    if (def.size() < nsel) {
      LOG(error) << "TopoFilter::Accept: Det[" << d << "] unsupported definition for [accept] cut of size["
                 << def.size() << "] expect " << nsel;
      break;
    }
    if (!filter(def)) return true;
  }
  return fAccept.size() ? false : true;
}
 
//_________________________________________________________________
template<>
bool CbmRecoQaTask::TopoFilter::Accept<ECbmModuleId::kTrd>(const CbmHit* h, const CbmCluster* c0,
                                                           const CbmCluster*) const
{
  if (!fReject.size() && !fAccept.size()) return true;
  std::vector<int> sel;
  auto nsel = CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kTrd>(h, sel);
  const CbmTrdCluster* c(nullptr);
  if (c0) c = dynamic_cast<const CbmTrdCluster*>(c0);
  auto htrd = dynamic_cast<const CbmTrdHit*>(h);
 
  // Lambda function to filter specific TRD(2D) parameters
  auto filter = [&](vector<int> def) -> size_t {
    size_t accept(trdSz);  // consider all extra cuts defined on the detector
    // check cut in topology [layer, module]
    size_t isel(0);
    for (; isel < nsel; isel++) {
      if (def[isel] < 0 || def[isel] == sel[isel])
        accept--;
      else
        return accept;
    }
    // check cut cluster size
    if (def[isel] < 0 || (c && c->GetSize() == def[isel]))
      accept--;
    else
      return accept;
    isel++;
    // check hit row cross
    auto rc = htrd->IsRowCross();
    if (def[isel] < 0 || (rc == bool(def[isel])))
      accept--;
    else
      return accept;
    isel++;
    // check hit maximum pairing
    auto mx = htrd->GetMaxType();
    if (def[isel] < 0 || (mx == bool(def[isel])))
      accept--;
    else
      return accept;
    isel++;
    // check hit maximum symmetry. Add 1 to shift it in the range {0, 1, 2} and skip the default [-1]
    auto symm = 0;  // htrd->GetMaxSymm() + 1; TODO need update of the CbmTrdHit class
    if (def[isel] < 0 || (symm == def[isel]))
      accept--;
    else
      return accept;
    isel++;
    // check cut cluster row
    if (def[isel] < 0 || (c && c->GetRow() == def[isel]))
      accept--;
    else
      return accept;
    isel++;
 
    return accept;
  };
 
  // check reject conditions
  for (auto def : fReject) {
    if (def.size() < nsel) {
      LOG(error) << "TopoFilter::Accept: Det[TRD] unsupported definition for [reject] cut of size[" << def.size()
                 << "] expect " << nsel;
      break;
    }
    if (!filter(def)) return false;
  }
  for (auto def : fAccept) {
    if (def.size() < nsel) {
      LOG(error) << "TopoFilter::Accept: Det[TRD] unsupported definition for [accept] cut of size[" << def.size()
                 << "] expect " << nsel;
      break;
    }
    if (!filter(def)) return true;
  }
  return fAccept.size() ? false : true;
}
 
//_________________________________________________________________
// HasAddress definitions
template<>
bool CbmRecoQaTask::View::HasAddress<CbmStsHit>(const CbmHit* h, double& x, double& y, double& dx, double& dy) const
{
  std::vector<int> sel;
  auto nsel = CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kSts>(h, sel);
  size_t isel(0);
  for (auto ii : fSelector) {
    if (ii != sel[isel]) break;
    isel++;
  }
  if (fSetup == eSetup::kDefault && sel[0] != fSelector[0]) {
    return false;
  }
  else if (fSetup != eSetup::kDefault && isel != nsel) {
    return false;
  }
  else
    LOG(debug4) << "Accept Sts hit for " << sel[0] << " " << sel[1] << " " << sel[2];
 
  const CbmStsHit* h0 = dynamic_cast<const CbmStsHit*>(h);
  if (!h0) {
    LOG(error) << "Failed loading STS hit in view " << name;
    return false;
  }
  x  = h0->GetX();
  dx = h0->GetDx();
  y  = h0->GetY();
  dy = h0->GetDy();
  return true;
}
 
//_________________________________________________________________
template<>
bool CbmRecoQaTask::View::HasAddress<CbmRichHit>(const CbmHit* h, double& x, double& y, double& dx, double& dy) const
{
  std::vector<int> sel;
  CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kRich>(h, sel);
 
  if (fSetup != eSetup::kDefault) {
    if (sel[0] >= fSelector[0]) return false;
  }
 
  const CbmRichHit* h0 = dynamic_cast<const CbmRichHit*>(h);
  if (!h0) {
    LOG(error) << "Failed loading RICH hit in view " << name;
    return false;
  }
  x  = h0->GetX();
  dx = h0->GetDx();
  y  = h0->GetY();
  dy = h0->GetDy();
  return true;
}
 
//_________________________________________________________________
template<>
bool CbmRecoQaTask::View::HasAddress<CbmTrdHit>(const CbmHit* h, double& x, double& y, double& dx, double& dy) const
{
  std::vector<int> sel;
  auto nsel = CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kTrd>(h, sel);
  /*
  if (fSetup == eSetup::kDefault) {
    sel = {static_cast<int>(CbmTrdAddress::GetLayerId(CbmTrdAddress::GetModuleAddress(h->GetAddress()))),
           static_cast<int>(CbmTrdAddress::GetModuleId(h->GetAddress())), -1};
  }
  else {
    sel = {static_cast<int>(CbmTrdAddress::GetModuleAddress(h->GetAddress())), -1, -1};
  }*/
 
  size_t isel(0);
  for (auto ii : fSelector) {
    if (ii != sel[isel]) break;
    isel++;
  }
 
  if (isel != nsel) return false;
 
  LOG(debug4) << "Accept TRD hit for " << sel[0] << " " << sel[1] << " " << sel[2];
 
  const CbmTrdHit* h0 = dynamic_cast<const CbmTrdHit*>(h);
  if (!h0) {
    LOG(error) << "Failed loading TRD hit in view " << name;
    return false;
  }
  x  = h0->GetX();
  dx = h0->GetDx();
  y  = h0->GetY();
  dy = h0->GetDy();
  return true;
}
 
//_________________________________________________________________
template<>
bool CbmRecoQaTask::View::HasAddress<CbmTofHit>(const CbmHit* h, double& x, double& y, double& dx, double& dy) const
{
  std::vector<int> sel;
  auto nsel = CbmRecoQaTask::GetHitDetectorId<ECbmModuleId::kTof>(h, sel);
  size_t isel(0);
  for (auto ii : fSelector) {
    if (ii != sel[isel]) break;
    isel++;
  }
  if (fSetup == eSetup::kDefault && nsel != 1) {
    return false;
  }
  else if (fSetup != eSetup::kDefault && isel != nsel) {
    return false;
  }
  else
    LOG(debug4) << "Accept Tof hit for " << sel[0] << " " << sel[1] << " " << sel[2];
 
  const CbmTofHit* h0 = dynamic_cast<const CbmTofHit*>(h);
  if (!h0) {
    LOG(error) << "Failed loading ToF hit in view " << name;
    return false;
  }
  x  = h0->GetX();
  dx = h0->GetDx();
  y  = h0->GetY();
  dy = h0->GetDy();
  return true;
}
 
//===============================================================================
//__________________ IMPLEMENTATIONS ____________________________________________
template<class Hit>
bool CbmRecoQaTask::View::Load(const CbmHit* h, const FairMCPoint* point, const CbmEvent* ev)
{
  double x(0), y(0), dx(0), dy(0);
  if (!HasAddress<Hit>(h, x, y, dx, dy)) {
    LOG(debug4) << "view " << name << " does not own hit " << h->ToString();
    return false;
  }
 
  int scale(0);
  TH2* hh(nullptr);
  auto fillView = [&](eProjectionType proj, double xx, double yy, bool scaleY = true) {
    auto it = fProjection.find(proj);
    if (it != fProjection.end()) {
      scale = std::get<0>(it->second);
      hh    = dynamic_cast<TH2*>(std::get<2>(it->second));
      if (hh) hh->Fill(xx, scaleY ? yy * scale : yy);
    }
  };
 
  if (ev) fillView(eProjectionType::kChdT, x, h->GetTime() - ev->GetStartTime());
 
  fillView(eProjectionType::kXYh, x, y, false);
 
  double event_time = ev ? ev->GetStartTime() : 1000;
 
  if (point) {
    // TODO dirty fix for MC point correction of Z coordinate. Assume straight
    // tracks from 0,0,0 !! Works for mCBM
    double prj = h->GetZ() / point->GetZ();
    fillView(eProjectionType::kXYhMC, prj * point->GetX(), prj * point->GetY(), false);
    fillView(eProjectionType::kResidualX, x, x - prj * point->GetX());
    fillView(eProjectionType::kResidualY, y, y - prj * point->GetY());
    fillView(eProjectionType::kResidualTX, x, h->GetTime() - event_time - point->GetTime());
    fillView(eProjectionType::kResidualTY, y, h->GetTime() - event_time - point->GetTime());
    fillView(eProjectionType::kPullX, x, (x - prj * point->GetX()) / dx);
    fillView(eProjectionType::kPullY, y, (y - prj * point->GetY()) / dy);
  }
  fMult++;  // register hit/point multiplicity
  return true;
}
 
bool CbmRecoQaTask::View::Load(const Trajectory::Node* n, const FairMCPoint* point, const int dinfo)
{
  const kf::TrackParamD& t = n->paramUp;
  // define aliases
  double xhit(n->measurementXY.X()), yhit(n->measurementXY.Y()), thit(n->measurementTime.T()),
    dx = xhit - t.X(), dy = yhit - t.Y(), dt = thit - t.Time(),
    pullx = dx / sqrt(n->measurementXY.Dx2() + t.GetCovariance(0, 0)),
    pully = dy / sqrt(n->measurementXY.Dy2() + t.GetCovariance(1, 1));
  // pullt = dt / sqrt(n->fMt.Dt2() + n->fTrack.GetCovariance(5, 5));
 
  for (auto& projection : fProjection) {
    int scale = get<0>(projection.second);
    auto h    = get<2>(projection.second);
    if (!h) continue;
 
    auto hh = dynamic_cast<TH2*>(h);
    switch (projection.first) {
      case eProjectionType::kXYa: hh->Fill(xhit, yhit); break;
      case eProjectionType::kXYp: hh->Fill(t.X(), t.Y()); break;
      case eProjectionType::kXdX: hh->Fill(xhit, scale * dx); break;
      case eProjectionType::kYdY: hh->Fill(yhit, scale * dy); break;
      case eProjectionType::kWdT: hh->Fill(xhit, scale * dt); break;
      case eProjectionType::kXpX: hh->Fill(xhit, pullx); break;
      case eProjectionType::kYpY: hh->Fill(yhit, pully); break;
      case eProjectionType::kXdXT: hh->Fill(xhit, scale * sqrt(t.GetCovariance(0, 0))); break;
      case eProjectionType::kYdYT: hh->Fill(yhit, scale * sqrt(t.GetCovariance(1, 1))); break;
 
      default: break;
    }
 
    if (point) {
      // TODO dirty fix for MC point correction of Z coordinate. Assume straight
      // tracks from 0,0,0 !! Works for mCBM
      double prj  = n->z / point->GetZ();
      double dxMC = prj * point->GetX() - t.X(), dyMC = prj * point->GetY() - t.Y();
 
      switch (projection.first) {
        case eProjectionType::kXdXMC: hh->Fill(xhit, scale * dxMC); break;
        case eProjectionType::kYdYMC: hh->Fill(yhit, scale * dyMC); break;
        default: break;
      }
    }
    if (dinfo != 0xab00) {
      auto hhh = dynamic_cast<TH3*>(h);
      switch (projection.first) {
        case eProjectionType::kDebug0: hhh->Fill(xhit, scale * dx, double(dinfo)); break;
        case eProjectionType::kDebug1: hhh->Fill(xhit, scale * dy, double(dinfo)); break;
        default: break;
      }
    }
  }
  return true;
}
 
bool CbmRecoQaTask::View::Load(TVector3* p, TVector3* pe)
{
  for (auto& projection : fProjection) {
    int scale = get<0>(projection.second);
    TH2* hh   = dynamic_cast<TH2*>(get<2>(projection.second));
    if (!hh) continue;
    switch (projection.first) {
      case eProjectionType::kPVmult:
        if (int(p->Z()) == -123) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      case eProjectionType::kDmult:
        if (int(p->Z()) == -124) hh->Fill(p->X(), p->Y());
        break;
      case eProjectionType::kXYt0:
        if (int(p->Z()) == 0) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      // fall through
      case eProjectionType::kXYt1:
        if (int(p->Z()) == 1) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      // fall through
      case eProjectionType::kXYt2:
        if (int(p->Z()) == 2) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      // fall through
      case eProjectionType::kXYt3:
        if (int(p->Z()) == 3) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      // fall through
      case eProjectionType::kXYt4:
        if (int(p->Z()) == 4) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      // fall through
      case eProjectionType::kXYt5:
        if (int(p->Z()) == 5) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      // fall through
      case eProjectionType::kPVxy:
        if (pe) hh->Fill(scale * p->X(), scale * p->Y());
        break;
      case eProjectionType::kPVxz:
        if (pe) hh->Fill(p->Z(), scale * p->X());
        break;
      case eProjectionType::kPVyz:
        if (pe) hh->Fill(p->Z(), scale * p->Y());
        break;
      case eProjectionType::kPVxsx:
        if (pe) hh->Fill(scale * p->X(), scale * pe->X());
        break;
      case eProjectionType::kPVysy:
        if (pe) hh->Fill(scale * p->Y(), scale * pe->Y());
        break;
      case eProjectionType::kPVzsz:
        if (pe) hh->Fill(scale * p->Z(), scale * pe->Z());
        break;
      default: break;
    }
  }
  return true;
}
 
//_____________________________________________________________________
 
void CbmRecoQaTask::Exec(Option_t*)
{
  LOG(info) << GetName() << "::Exec : Evs[" << (fEvents ? fEvents->GetEntriesFast() : 0) << "] Trks["
            << (fGTracks ? fGTracks->GetEntriesFast() : 0) << "] Hits["
            << (fHits[ECbmModuleId::kSts] ? fHits[ECbmModuleId::kSts]->GetEntriesFast() : 0) << " "
            << (fHits[ECbmModuleId::kTrd] ? fHits[ECbmModuleId::kTrd]->GetEntriesFast() : 0) << " "
            << (fHits[ECbmModuleId::kTof] ? fHits[ECbmModuleId::kTof]->GetEntriesFast() : 0) << "]";
 
  int iev = 0, itrack = 0;
  // match to MC lambda function
  auto matchMC = [&](const CbmHit* hit, int hid, ECbmModuleId d, const CbmEvent* ev = nullptr) -> const FairMCPoint* {
    const FairMCPoint* mcpoint = nullptr;
    if (fHitMatch[d] && fPoints[d]) {
      if (d == ECbmModuleId::kRich) {
        Int_t pointID = hit->GetRefId();
        if (pointID >= 0) {
          // mcpoint = dynamic_cast<FairMCPoint*>(fPoints[n.fHitSystemId]->At(pointID));
        }
      }
      else {
        auto match = dynamic_cast<const CbmMatch*>(fHitMatch[d]->At(hid));
        if (match && match->GetNofLinks() > 0) {
          const auto& link = match->GetMatchedLink();
          if (ev) {
            int file_id{0}, event_id{0};
            if (ev && ev->GetMatch() && ev->GetMatch()->GetNofLinks() > 0) {
              file_id  = ev->GetMatch()->GetMatchedLink().GetFile();
              event_id = ev->GetMatch()->GetMatchedLink().GetEntry();
            }
            else {
              event_id = FairRootManager::Instance()->GetEntryNr();
            }
            if (link.GetFile() != file_id || link.GetEntry() != event_id) {
              LOG(warn) << "match from different event";
            }
          }
          mcpoint = dynamic_cast<FairMCPoint*>(fPoints[d]->Get(link));
        }
      }
    }
    return mcpoint;
  };
  // hit processor lambda function
  auto processHits = [&](CbmEvent* ev) {
    for (auto& detp : fDetQa) {
      auto& det = detp.second;
      // skip detectors from setup for which reco was not activated
      switch (det.id) {
        case ECbmModuleId::kSts:
          if (!fuRecoConfig[kStsHits]) continue;
          break;
        case ECbmModuleId::kMuch:
          if (!fuRecoConfig[kMuchHits]) continue;
          break;
        case ECbmModuleId::kRich:
          if (!fuRecoConfig[kRichHits]) continue;
          break;
        case ECbmModuleId::kTrd2d:
        case ECbmModuleId::kTrd:
          if (!fuRecoConfig[kTrdHits]) continue;
          break;
        case ECbmModuleId::kTof:
          if (!fuRecoConfig[kTofHits]) continue;
          break;
        default: break;
      }
      if (!fHits[det.id]) {
        LOG(error) << GetName() << "::Exec() : Hits for " << cbm::util::ToString(det.id) << " not available. Skip.";
        continue;
      }
 
      const int nh = (ev) ? max(int(0), int(ev->GetNofData(det.hit.id))) : fHits[det.id]->GetEntriesFast();
      for (int ih = 0; ih < nh; ++ih) {
        const int jh = (ev) ? ev->GetIndex(det.hit.id, ih) : ih;
        // retrieve hit
        auto hit = dynamic_cast<const CbmHit*>(fHits[det.id]->At(jh));
        if (!hit) {
          LOG(warning) << GetName() << "::Exec() : Hit " << jh << " for " << det.id << " not available. Skip.";
          continue;
        }
        // filter hits for analysis for which specific properties are requiered
        if (!FilterHit(hit, det.id)) continue;
 
        // link MC info if available
        const FairMCPoint* mcpoint = matchMC(hit, jh, det.id, ev);
        // upload info [hit, point, event] to the corresponding view
        bool ret(false);
        for (auto& view : det.fViews) {
          switch (det.id) {
            case ECbmModuleId::kMvd: ret = view.Load<CbmMvdHit>(hit, mcpoint, ev); break;
            case ECbmModuleId::kSts: ret = view.Load<CbmStsHit>(hit, mcpoint, ev); break;
            case ECbmModuleId::kMuch: ret = view.Load<CbmMuchPixelHit>(hit, mcpoint, ev); break;
            case ECbmModuleId::kRich: ret = view.Load<CbmRichHit>(hit, mcpoint, ev); break;
            case ECbmModuleId::kTrd:
            case ECbmModuleId::kTrd2d: ret = view.Load<CbmTrdHit>(hit, mcpoint, ev); break;
            case ECbmModuleId::kTof: ret = view.Load<CbmTofHit>(hit, mcpoint, ev); break;
            case ECbmModuleId::kFsd: ret = view.Load<CbmFsdHit>(hit, mcpoint, ev); break;
            case ECbmModuleId::kPsd: ret = view.Load<CbmPsdHit>(hit, mcpoint, ev); break;
            default: LOG(fatal) << GetName() << "::Exec : unsupported det " << det.id; break;
          }
          if (ret) break;
        }
      }
      TVector3 mult;
      for (auto& v : det.fViews) {
        mult.SetXYZ(nh, double(v.fMult) / nh, -124);
        v.Load(&mult);
        v.fMult = 0;
      }
    }
  };
 
  auto processTracks = [&](CbmEvent* ev) {
    const int ntrk = (ev) ? ev->GetNofData(ECbmDataType::kGlobalTrack) : fGTracks->GetEntriesFast();
    // read in the vertex and the list of tracks used for its definition
    TVector3 pvx, evx;
    int ntrkDet[3] = {0};
 
    CbmVertex* vx(nullptr);
    int nTrkVx(0);
    if (ev && fViews.find("vx") != fViews.end()) {
      vx      = ev->GetVertex();
      nTrkVx  = vx->GetNTracks();
      auto& v = fViews["vx"];
      if (nTrkVx >= 2) {
        vx->Position(pvx);
        for (int i(0); i < 3; i++) {
          evx[i] = vx->GetCovariance(i, i);
          if (evx[i] > 0.) evx[i] = sqrt(evx[i]);
        }
        v.Load(&pvx, &evx);
      }
      pvx.SetXYZ(ntrk, nTrkVx, -123);
      v.Load(&pvx);
    }
 
    // find the track projection view if it is registered
    View* viewPrj = nullptr;
    if (fViews.find("prj") != fViews.end()) {
      viewPrj = &fViews["prj"];
      if (viewPrj->fType != eViewType::kTrkProj) {
        LOG(warning) << GetName() << "::Exec: view for track projection with wrong type. Skipping.";
        viewPrj = nullptr;
      }
    }
 
    for (int itrk = 0; itrk < ntrk; ++itrk) {
      int trkIdx = ev ? ev->GetIndex(ECbmDataType::kGlobalTrack, itrk) : itrk;
      auto track = dynamic_cast<const CbmGlobalTrack*>((*fGTracks)[trkIdx]);
 
      // FILTER on TRACK HIT COMPOSITION/TOPOLOGY
      if (!FilterTrack(track)) continue;
      if (nTrkVx >= 2) {
        if (vx->FindTrackByIndex(trkIdx)) {
          if (track->GetStsTrackIndex() >= 0) ntrkDet[0]++;
          if (track->GetTrdTrackIndex() >= 0) ntrkDet[1]++;
          if (track->GetTofTrackIndex() >= 0) ntrkDet[2]++;
        }
      }
 
      // CREATE THE TRAJECTORY FOR THE TRACK.
      // When some projection views are required, create trajectory nodes for all material layers.
      // otherwise, materials outside of the track hit area will be ignored.
      Trajectory trkKf;
      if (!fFitter.CreateFromGlobalTrack(trkKf, *track, (viewPrj ? true : false))) {
        LOG(fatal) << GetName() << "::Exec: can not create the track for the fit!";
        break;
      }
 
      // ADD REFERENCE TRAJECTORY NODES W/O MEASUREMENTS
      // for the track projections
      int nprj(0);
      if (viewPrj) {
        std::vector<Trajectory::Node> nodesPrj;
        for (auto plane : fPrjPlanes) {
          Trajectory::Node n;
          n.z = plane.Z();
          Trajectory::SetPlaneId(n, nprj++);
          nodesPrj.push_back(n);
        }
        trkKf.AddNodes(nodesPrj);
      }
 
      // SWITCH OFF DETECTORS
      // which are skipped for tracking by user
      for (unsigned int i = 0; i < trkKf.GetNofNodes(); ++i) {
        const auto& n = trkKf.GetNode(i);
        // special case TRD2D. remap det id to TRD
        auto nodeSystemId =
          (Trajectory::GetHitSystemId(n) == ECbmModuleId::kTrd2d ? ECbmModuleId::kTrd : Trajectory::GetHitSystemId(n));
        if (fDetQa.find(nodeSystemId) == fDetQa.end()) continue;
        switch (nodeSystemId) {
          case ECbmModuleId::kSts:
            if (!IsActive(eRecoConfig::kUseSts)) {
              trkKf.DisableMeasurementAtNode(i);
            }
            break;
          case ECbmModuleId::kTrd:
            if (!IsActive(eRecoConfig::kUseTrd)) {
              trkKf.DisableMeasurementAtNode(i);
            }
            break;
          case ECbmModuleId::kTof:
            if (!IsActive(eRecoConfig::kUseTof)) {
              trkKf.DisableMeasurementAtNode(i);
            }
            break;
          default: break;
        }
      }
 
      // FIT THE TRAJECTORY
      // first time
      if (!fFitter.FitTrajectory(trkKf)) {
        continue;
      }
      // second time using the previous results for the accurate linearization of equations
      if (!fFitter.FitTrajectory(trkKf)) {
        continue;
      }
 
      // FILTER ON TRACK PROJECTION
      if (!FilterTrack(&trkKf)) continue;
      LOG(debug) << GetName() << "::Exec: Accept track for QA.";
 
      // minimalistic QA tool for tracks used for target projection
      int nhit[(size_t) ECbmModuleId::kLastModule] = {0};
      for (int iNode = 0; iNode < (int) trkKf.GetNofNodes(); iNode++) {
        auto& n     = trkKf.GetNode(iNode);
        int planeId = Trajectory::GetPlaneId(n);
 
        if (viewPrj && planeId >= 0 && n.isFitted) {            // this is a projection node
          TVector3 xyz(n.paramUp.X(), n.paramUp.Y(), planeId);  // to communicate to the Load function
          viewPrj->Load(&xyz);
          continue;
        }
        // special case TRD2D. remap det id to TRD
        auto nodeSystemId =
          (Trajectory::GetHitSystemId(n) == ECbmModuleId::kTrd2d ? ECbmModuleId::kTrd : Trajectory::GetHitSystemId(n));
        int hitIndex     = Trajectory::GetHitIndex(n);
        // check if hit is well defined and detector is registered
        if (nodeSystemId == ECbmModuleId::kNotExist || fDetQa.find(nodeSystemId) == fDetQa.end()) continue;
        auto& det = fDetQa[nodeSystemId];
        // det.Print();
        auto view = det.FindView(n.measurementXY.X(), n.measurementXY.Y(), n.z);
        if (!view) {
          LOG(warning) << GetName() << "::Exec: " << cbm::util::ToString(det.id)
                       << " view not defined for hit @ x=" << n.measurementXY.X() << " y=" << n.measurementXY.Y()
                       << " z=" << n.z << ". Skipping.";
          LOG(debug1) << "This is symptomatic for using different alignment matrices during RECO and QA phases.";
          continue;
        }
        auto hit = dynamic_cast<CbmHit*>(fHits[nodeSystemId]->At(hitIndex));
        if (!hit) {
          LOG(error) << GetName() << "::Exec: " << cbm::util::ToString(det.id) << " hit @ id=" << hitIndex
                     << " missing. Skipping.";
          LOG(debug1) << "The hit index present in the track but missing in the hit array.";
          continue;
        }
        // match to MC
        auto mcpoint = matchMC(hit, hitIndex, nodeSystemId);
        // refit the trajectory without this hit to get unbiased residual.
        // use a temporary copy of the trajectory to keep the best linearization in trkKf
        auto trkTmp      = trkKf;
        const auto& nTmp = trkTmp.GetNode(iNode);
        trkTmp.DisableMeasurementAtNode(iNode);  // don't use the hit measurement in the fit
        if (!fFitter.FitTrajectory(trkTmp) || !nTmp.isFitted) continue;
        // nTmp contains now the trajectory at the hit Z position, but fitted without the hit
        if (view->HasDebug())
          view->Load(&nTmp, mcpoint, GetDebugInfo<ECbmModuleId::kTrd>(hit));
        else
          view->Load(&nTmp, mcpoint);
        nhit[(int) nodeSystemId]++;
      }
 
      itrack++;
    }
  };
 
  if (fEvents) {
    for (auto evObj : *fEvents) {
      auto ev = dynamic_cast<CbmEvent*>(evObj);
      if (!FilterEvent(ev)) continue;
      processHits(ev);
      if (fGTracks) processTracks(ev);
      iev++;
    }
    LOG(info) << GetName() << "::Exec : Evs(%)["
              << (fEvents->GetEntriesFast() ? 1.e2 * iev / fEvents->GetEntriesFast() : 0) << "] Trks(%)["
              << (fGTracks && fGTracks->GetEntriesFast() ? 1.e2 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
    return;  // END EbyE QA
  }
  else {
    processHits(nullptr);
 
    if (!fGTracks) {
      LOG(info) << GetName() << "::Exec : TS local reco only.";
      return;  // END TS QA (no tracking)
    }
 
    processTracks(nullptr);
  }
  LOG(info) << GetName() << "::Exec : Trks(%)["
            << (fGTracks && fGTracks->GetEntriesFast() ? 1.e2 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
  return;
}
 
//_____________________________________________________________________
int CbmRecoQaTask::GetNviews(eViewType type) const
{
  int n(0);
  for (auto v : fViews) {
    if (v.second.fType != type) continue;
    n++;
  }
  return n;
}
 
 
//_____________________________________________________________________
void CbmRecoQaTask::Finish()
{
  FairSink* sink = FairRootManager::Instance()->GetSink();
  sink->WriteObject(&fOutFolder, nullptr);
}
 
//_____________________________________________________________________
bool CbmRecoQaTask::FilterEvent(const CbmEvent* ptr) const
{
  // sanity checks
  if (ptr->GetStartTime() < 0) return false;  // this is found in MC
 
  for (auto evCut : fFilterEv) {
    if (!evCut.Accept(ptr, this)) return false;
  }
  return true;
}
 
//_____________________________________________________________________
bool CbmRecoQaTask::FilterHit(const CbmHit* hit, ECbmModuleId d) const
{
  bool ret(true);
  stringstream ss;
  for (auto tf : fFilterTrk) {
    if (tf.fDet != d) continue;
    auto filter = tf.GetHitsFilter();
    if (!filter) continue;
    const CbmCluster *cl0(nullptr), *cl1(nullptr);
    switch (d) {
      case ECbmModuleId::kSts:
        cl0 = (CbmCluster*) GetData<CbmCluster>(d, ((CbmStsHit*) hit)->GetFrontClusterId());
        cl1 = (CbmCluster*) GetData<CbmCluster>(d, ((CbmStsHit*) hit)->GetBackClusterId());
        ret = filter->Accept<ECbmModuleId::kSts>(hit, cl0, cl1);
        break;
      case ECbmModuleId::kTrd:
        cl0 = (CbmCluster*) GetData<CbmCluster>(d, hit->GetRefId());
        ret = filter->Accept<ECbmModuleId::kTrd>(hit, cl0);
        break;
      case ECbmModuleId::kTof: ret = filter->Accept<ECbmModuleId::kTof>(hit); break;
      default: break;
    }
    ss << cbm::util::ToString(d) << " topology hit :\n"
       << hit->ToString() << (d == ECbmModuleId::kSts ? "\n\t" : "\t") << (cl0 ? cl0->ToString() : "") << "\t"
       << (cl1 ? cl1->ToString() : "");
    if (!ret) break;
  }
  if (!ret) {
    LOG(debug) << "FilterHit::Reject " << cbm::util::ToString(d) << " hit.";
    LOG(debug1) << ss.str();
  }
  // else LOG(info) << "FilterHit::Accept " << cbm::util::ToString(d) << " hit.";
  return ret;
}
 
//_____________________________________________________________________
bool CbmRecoQaTask::FilterTrack(const CbmGlobalTrack* ptr) const
{
  for (auto trkCut : fFilterTrk) {
    // special cut for selection on track projection. See FilterTrack(const cbm::algo::kf::Trajectory<>*)
    if (trkCut.fDet == ECbmModuleId::kRef) continue;
    if (!trkCut.Accept(ptr, this)) return false;
  }
  LOG(debug) << GetName() << "::FilterTrack : Accept hit cuts Sts[" << ptr->GetStsTrackIndex() << "] Trd["
             << ptr->GetTrdTrackIndex() << "] Tof[" << ptr->GetTofTrackIndex() << "].";
  return true;
}
 
//_____________________________________________________________________
bool CbmRecoQaTask::FilterTrack(const Trajectory* ptr) const
{
  bool ret = true;
  // check projections
  for (auto& n : ptr->GetNodes()) {
    if (Trajectory::GetPlaneId(n) < 0) continue;
    for (auto trkCut : fFilterTrk) {
      if (trkCut.fDet != ECbmModuleId::kRef) continue;
      for (auto rcut : trkCut.GetProjCuts()) {
        if (abs(rcut[0] - n.z) >= 1.e-3) continue;
        double r = sqrt((n.paramUp.X() - rcut[1]) * (n.paramUp.X() - rcut[1])
                        + (n.paramUp.Y() - rcut[2]) * (n.paramUp.Y() - rcut[2]));
        if (r > rcut[3]) ret = false;
        if (!ret) {
          LOG(debug1) << GetName() << "::FilterTrack : Reject projection @ z[" << rcut[0] << "] r=" << r
                      << " R=" << rcut[3];
          break;
        }
      }
    }
  }
  if (ret) LOG(debug) << GetName() << "::FilterTrack : Accept projection cuts.";
  return ret;
}
 
//_____________________________________________________________________
TString CbmRecoQaTask::GetGeoTagForDetector(const TString& detector)
{
  gGeoManager->CdTop();
  TGeoNode* cave = gGeoManager->GetCurrentNode();
  if (!cave) {
    LOG(error) << "Error: Could not get the top node in the geometry manager." << std::endl;
    return "";
  }
 
  for (Int_t iNode = 0; iNode < cave->GetNdaughters(); iNode++) {
    TString name = cave->GetDaughter(iNode)->GetVolume()->GetName();
    if (name.Contains(detector, TString::kIgnoreCase)) {
      return name.Contains("mcbm", TString::kIgnoreCase) ? TString(name(5, name.Length()))
                                                         : TString(name(5, name.Length() - 5));
    }
  }
  return "";
}
 
//_____________________________________________________________________
std::vector<TString> CbmRecoQaTask::GetPath(TGeoNode* node, TString detector, TString activeNodeName, int depth,
                                            const TString& path)
{
  std::vector<TString> nodePaths;
 
  if (!node) {
    return nodePaths;
  }
 
  TString nodePath = path;
  if (!nodePath.IsNull()) {
    nodePath += "/";
  }
  nodePath += node->GetName();
 
  if (TString(node->GetName()).Contains(activeNodeName)) {
    if (nodePath.Contains(detector)) nodePaths.push_back(nodePath);
  }
 
  // Recursively traverse the daughters
  Int_t numDaughters = node->GetNdaughters();
  for (Int_t i = 0; i < numDaughters; ++i) {
    TGeoNode* daughterNode = node->GetDaughter(i);
 
    std::vector<TString> result = GetPath(daughterNode, detector, activeNodeName, depth + 1, nodePath);
    nodePaths.insert(nodePaths.end(), result.begin(), result.end());
  }
 
  return nodePaths;
}
//_____________________________________________________________________
void CbmRecoQaTask::InitMcbm22()
{
  CbmSetup* setup = CbmSetup::Instance();
  if (!setup) {
    LOG(fatal) << GetName() << "::InitMcbm22() : Missing setup definition.";
    return;
  }
  // fixed momentum no magentic field for mCBM
  fFitter.SetDefaultMomentumForMs(0.5);
  fFitter.SetEnforceDefaultMomentumForMs();
 
  View* v = nullptr;
 
  std::vector<TString> path;
  TGeoNode* topNode = gGeoManager->GetTopNode();
  if (!topNode) {
    LOG(error) << "Error: Top node not found.";
    return;
  }
 
  auto processDetector = [&](const std::string& detector, const std::string& component) {
    TString geoTag = GetGeoTagForDetector(detector);
 
    if (geoTag.Length() > 0) {
      LOG(info) << detector << ": geometry tag is " << geoTag;
    }
    else {
      LOG(warn) << "Warning: No geometry tag found for detector " << detector;
    }
 
    path = GetPath(topNode, detector, component, 0);
  };
 
  if (setup->IsActive(ECbmModuleId::kSts)) {
 
    processDetector("sts", "Sensor");
 
    std::regex pattern("/Station(\\d+)_(\\d+)/Ladder(\\d+)_(\\d+)/HalfLadder\\d+d_(\\d+)/"
                       "HalfLadder\\d+d_Module(\\d+)_(\\d+)/Sensor(\\d+)_(\\d+)");
 
    Detector* sts = AddDetector(ECbmModuleId::kSts);
    for (const auto& str : path) {
      std::string Str(str.Data());
      std::smatch match;
      if (std::regex_search(Str, match, pattern)) {
        int station = std::stoi(match[2]);
        int ladder  = std::stoi(match[4]);
        int module  = std::stoi(match[7]);
 
        const char* fType = Form("U%dL%dM%d", station - 1, ladder - 1, module - 1);
 
        v = sts->AddView(fType, str.Data(), {station - 1, ladder - 1, module - 1});
        v->SetProjection(eProjectionType::kXdX, (station == 1) ? 1 : 0.75, "mm");
        v->SetProjection(eProjectionType::kYdY, (station == 1) ? 4 : 3, "mm");
        if (fuRecoConfig[kUseMC]) {
          v->SetProjection(eProjectionType::kXdXMC, 1, "mm");
          v->SetProjection(eProjectionType::kYdYMC, 3, "mm");
          v->SetProjection(eProjectionType::kResidualX, 200, "um");
          v->SetProjection(eProjectionType::kResidualY, 500, "um");
        }
      }
      else {
        std::cout << "No match found in string: " << str << std::endl;
      }
    }
  }
 
  if (setup->IsActive(ECbmModuleId::kTrd)) {
    processDetector("trd", "module");
 
    // Trd2d
    Detector* trd = AddDetector(ECbmModuleId::kTrd);
    for (const auto& str : path) {
      if (!str.Contains("module9")) continue;
      v = trd->AddView("2D", str.Data(), {5});
      v->SetProjection(eProjectionType::kChdT, 400, "ns");
      v->SetProjection(eProjectionType::kXdX, 10, "mm");
      v->SetProjection(eProjectionType::kYdY, 20, "mm");
      if (fuRecoConfig[kUseMC]) {
        v->SetProjection(eProjectionType::kXdXMC, 10, "mm");
        v->SetProjection(eProjectionType::kYdYMC, 20, "mm");
        v->SetProjection(eProjectionType::kResidualX, 1.5, "mm");
        v->SetProjection(eProjectionType::kResidualY, 5.0, "mm");
      }
    }
    // Trd1DxDy
    for (const auto& str : path) {
      LOG(info) << str;
      if (str.Contains("layer02")) {
        v = trd->AddView("1Dx", str.Data(), {21});
        v->SetProjection(eProjectionType::kChdT, 400, "ns");
        v->SetProjection(eProjectionType::kXdX, 1.5, "cm");
        if (fuRecoConfig[kUseMC]) {
          v->SetProjection(eProjectionType::kXdXMC, 1.5, "cm");
          v->SetProjection(eProjectionType::kResidualX, 1.5, "mm");
        }
      }
      if (str.Contains("layer03")) {
        v = trd->AddView("1Dy", str.Data(), {37});
        v->SetProjection(eProjectionType::kChdT, 400, "ns");
        v->SetProjection(eProjectionType::kYdY, 1.5, "cm");
        if (fuRecoConfig[kUseMC]) {
          v->SetProjection(eProjectionType::kYdYMC, 1.5, "cm");
          v->SetProjection(eProjectionType::kResidualY, 1.5, "mm");
        }
      }
    }
  }
 
 
  if (setup->IsActive(ECbmModuleId::kTof)) {
    processDetector("tof", "counter");
 
    Detector* tof = AddDetector(ECbmModuleId::kTof);
    tof->hit.name = "TofHit";
    std::regex pattern("module_(\\d+)_(\\d+)/gas_box_(\\d+)/counter_(\\d+)");
 
    for (const auto& str : path) {
      std::string Str(str.Data());
      std::smatch match;
      if (std::regex_search(Str, match, pattern)) {
        int type         = std::stoi(match[1]);
        int smid         = std::stoi(match[2]);
        int rpc          = std::stoi(match[4]);
        const char* name = Form("Sm%d_%dRpc%d", type, smid, rpc);
        if (type == 0) {
          v = tof->AddView(name, str.Data(), {smid, type, rpc});
          v->SetProjection(eProjectionType::kChdT, 15, "ns");
        }
        else {
          v = tof->AddView(name, str.Data(), {smid, type, rpc});
        }
      }
      else {
        std::cout << "No match found in string: " << str << std::endl;
      }
    }
  }
  // ===============================================================================
  // TRG - upstream projections
  float angle = 25., L[] = {14.3, 0, -20, -38 /*-1*/, -50.5 /*-4.1*/};
  //R[]= {2.5, 3, 0.5, 0.6, 0.6}, dx[5], dz[5];
  for (int i(0); i < 5; i++) {
    fPrjPlanes.emplace_back(L[i] * TMath::Sin(angle * TMath::DegToRad()), 0.,
                            L[i] * TMath::Cos(angle * TMath::DegToRad()));
  }
  //     (cm); T_{trk} - T_{ev} (ns)", prj.first),
}
 
//_____________________________________________________________________
void CbmRecoQaTask::InitMcbm24()
{
  CbmSetup* setup = CbmSetup::Instance();
  if (!setup) {
    LOG(fatal) << GetName() << "::InitMcbm24() : Missing setup definition.";
    return;
  }
  // fixed momentum no magentic field for mCBM
  fFitter.SetDefaultMomentumForMs(1.5);
  fFitter.SetEnforceDefaultMomentumForMs();
 
  TString dtag;
  View* v(nullptr);
  if (setup->IsActive(ECbmModuleId::kSts)) {
    // setup->GetGeoTag(ECbmModuleId::kSts, dtag);
    // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
    dtag          = "v24c_mcbm";
    Detector* sts = AddDetector(ECbmModuleId::kSts);
    std::map<const char*, std::pair<const char*, std::vector<int>>> sensors = {
      {"U0L0M0", {"Station01_1/Ladder13_1/HalfLadder13u_1/HalfLadder13u_Module03_1/Sensor03_1", {0, 0, 0}}},
      {"U1L0M0", {"Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", {1, 0, 0}}},
      {"U1L0M1", {"Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", {1, 0, 1}}},
      {"U1L1M0", {"Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", {1, 1, 0}}},
      {"U1L1M1", {"Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", {1, 1, 1}}},
      {"U2L0M0", {"Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module03_1/Sensor03_1", {2, 0, 0}}},
      {"U2L0M1", {"Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module04_2/Sensor04_1", {2, 0, 1}}},
      {"U2L1M0", {"Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module03_1/Sensor03_1", {2, 1, 0}}},
      {"U2L1M1", {"Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module04_2/Sensor04_1", {2, 1, 1}}},
      {"U2L2M0", {"Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_1/Sensor03_1", {2, 2, 0}}},
      {"U2L2M1", {"Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1", {2, 2, 1}}},
      {"U2L2M2", {"Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1", {2, 2, 2}}}};
    for (auto& [name, def] : sensors) {
      v = sts->AddView(name, Form("/cave_1/sts_%s_0/%s", dtag.Data(), def.first), def.second);
      v->SetProjection(eProjectionType::kXdX, 1000, "um");
      v->SetProjection(eProjectionType::kYdY, 1500, "um");
      v->SetProjection(eProjectionType::kXdXT, 500, "um");
      v->SetProjection(eProjectionType::kYdYT, 500, "um");
    }
  }
  if (setup->IsActive(ECbmModuleId::kTrd)) {
    // setup->GetGeoTag(ECbmModuleId::kTrd, dtag);
    // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
    dtag          = "v24e_mcbm";
    Detector* trd = AddDetector(ECbmModuleId::kTrd);
    // Trd2D
    v = trd->AddView("2D", Form("/cave_1/trd_%s_0/layer01_20101/module9_101001001/gas_0", dtag.Data()), {0, 0});
    v->SetProjection(eProjectionType::kChdT, 300, "ns");
    // Trd1Dx
    v = trd->AddView("1Dx", Form("/cave_1/trd_%s_0/layer02_10202/module5_101002001", dtag.Data()), {1, 0});
    v->SetProjection(eProjectionType::kChdT, 300, "ns");
    v->SetProjection(eProjectionType::kXdX, 1.5, "cm");
    v->SetProjection(eProjectionType::kYdY, 3.0, "cm");
    // Trd1Dy
    v = trd->AddView("1Dy", Form("/cave_1/trd_%s_0/layer03_11303/module5_101103001", dtag.Data()), {2, 0});
    v->SetProjection(eProjectionType::kChdT, 300, "ns");
    v->SetProjection(eProjectionType::kXdX, 1.5, "cm");
    v->SetProjection(eProjectionType::kYdY, 3.0, "cm");
  }
  if (setup->IsActive(ECbmModuleId::kTof)) {
    // setup->GetGeoTag(ECbmModuleId::kTof, dtag);
    // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
    dtag          = "v24d_mcbm";
    Detector* tof = AddDetector(ECbmModuleId::kTof);
    vector<int> tofSelect(3);
    // add type 0
    for (int ism(0); ism < 6; ism++) {
      tofSelect[0] = ism;
      for (int irpc(0); irpc < 5; irpc++) {
        tofSelect[2] = irpc;
        v            = tof->AddView(Form("Sm%dRpc%d", ism, irpc),
                         Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/"
                              "gas_box_0/counter_%d",
                              dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
                         tofSelect);
        v->SetProjection(eProjectionType::kXdX, 3, "cm");
        v->SetProjection(eProjectionType::kYdY, 1, "cm");
      }
    }
    // add type 6 (Buch)
    tofSelect[0] = 0;
    tofSelect[1] = 6;
    for (int irpc(0); irpc < 2; irpc++) {
      tofSelect[2] = irpc;
      tof->AddView(Form("BuchRpc%d", irpc),
                   Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/"
                        "counter_%d",
                        dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
                   tofSelect);
    }
    // add type 9
    tofSelect[1] = 9;
    for (int ism(0); ism < 2; ism++) {
      tofSelect[0] = ism;
      for (int irpc(0); irpc < 2; irpc++) {
        tofSelect[2] = irpc;
        tof->AddView(Form("Test%dRpc%d", ism, irpc),
                     Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/"
                          "gas_box_0/counter_%d",
                          dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
                     tofSelect);
      }
    }
    // add type 2
    tofSelect[1] = 2;
    for (int ism(0); ism < 2; ism++) {
      tofSelect[0] = ism;
      for (int irpc(0); irpc < 5; irpc++) {
        tofSelect[2] = irpc;
        tof->AddView(Form("Sm2%dRpc%d", ism, irpc),
                     Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/"
                          "gas_box_0/counter_%d",
                          dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
                     tofSelect);
      }
    }
  }
  if (setup->IsActive(ECbmModuleId::kRich)) {
    // setup->GetGeoTag(ECbmModuleId::kRich, dtag);
    // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
    dtag           = "v24a_mcbm";
    Detector* rich = AddDetector(ECbmModuleId::kRich);
    // Aerogel 1
    rich->AddView("Aerogel", Form("/cave_1/rich_%s_0/box_1/Gas_1", dtag.Data()), {4});
  }
  // ===============================================================================
  // TRG - upstream projections
  float angle = 25., L[] = {14.3, 0, -20, -38, -50.5};
  //R[]= {2.5, 3, 0.5, 0.6, 0.6}, dx[5], dz[5];
  for (int i(0); i < 5; i++) {
    fPrjPlanes.emplace_back(L[i] * TMath::Sin(angle * TMath::DegToRad()), 0.,
                            L[i] * TMath::Cos(angle * TMath::DegToRad()));
  }
}
 
//_____________________________________________________________________
void CbmRecoQaTask::InitMcbm25()
{
  CbmSetup* setup = CbmSetup::Instance();
  if (!setup) {
    LOG(fatal) << GetName() << "::InitMcbm25() : Missing setup definition.";
    return;
  }
  // fixed momentum no magentic field for mCBM
  fFitter.SetDefaultMomentumForMs(1.5);
  fFitter.SetEnforceDefaultMomentumForMs();
 
  TString dtag;
  View* v(nullptr);
  if (setup->IsActive(ECbmModuleId::kSts)) {
    // setup->GetGeoTag(ECbmModuleId::kSts, dtag);
    // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
    dtag                                                                    = "v24c_mcbm";
    Detector* sts                                                           = AddDetector(ECbmModuleId::kSts);
    std::map<const char*, std::pair<const char*, std::vector<int>>> sensors = {
      {"U0L0M0", {"Station01_1/Ladder13_1/HalfLadder13u_1/HalfLadder13u_Module03_1/Sensor03_1", {0, 0, 0}}},
      {"U1L0M0", {"Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", {1, 0, 0}}},
      {"U1L0M1", {"Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", {1, 0, 1}}},
      {"U1L1M0", {"Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", {1, 1, 0}}},
      {"U1L1M1", {"Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", {1, 1, 1}}},
      {"U2L0M0", {"Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module03_1/Sensor03_1", {2, 0, 0}}},
      {"U2L0M1", {"Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module04_2/Sensor04_1", {2, 0, 1}}},
      {"U2L1M0", {"Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module03_1/Sensor03_1", {2, 1, 0}}},
      {"U2L1M1", {"Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module04_2/Sensor04_1", {2, 1, 1}}},
      {"U2L2M0", {"Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_1/Sensor03_1", {2, 2, 0}}},
      {"U2L2M1", {"Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1", {2, 2, 1}}},
      {"U2L2M2", {"Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1", {2, 2, 2}}}};
    for (auto& [name, def] : sensors) {
      v = sts->AddView(name, Form("/cave_1/sts_%s_0/%s", dtag.Data(), def.first), def.second);
      v->SetProjection(eProjectionType::kXdX, 1000, "um");
      if (fuRecoConfig[kUseMC]) v->SetProjection(eProjectionType::kXdXMC, 1000, "um");
      v->SetProjection(eProjectionType::kYdY, 1500, "um");
      if (fuRecoConfig[kUseMC]) v->SetProjection(eProjectionType::kYdYMC, 1500, "um");
      v->SetProjection(eProjectionType::kXdXT, 500, "um");
      v->SetProjection(eProjectionType::kYdYT, 500, "um");
      if (fuRecoConfig[kUseMC]) {
        v->SetProjection(eProjectionType::kXdXMC, 1000, "um");
        v->SetProjection(eProjectionType::kYdYMC, 1500, "um");
        v->SetProjection(eProjectionType::kResidualX, 200, "um");
        v->SetProjection(eProjectionType::kResidualY, 500, "um");
      }
    }
  }
  if (setup->IsActive(ECbmModuleId::kTrd)) {
    // setup->GetGeoTag(ECbmModuleId::kTrd, dtag);
    // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
    dtag            = "v24e_mcbm";
    Detector* trd                                                           = AddDetector(ECbmModuleId::kTrd);
    std::map<const char*, std::pair<const char*, std::vector<int>>> modules = {
      {"2D", {"layer01_20101/module9_101001001/gas_0", {0, 0}}},
      {"1Dx", {"layer02_10202/module5_101002001", {1, 0}}},
      {"1Dy", {"layer03_11303/module5_101103001", {2, 0}}},
    };
    for (auto& [name, def] : modules) {
      bool debug = true;
      if (strcmp("2D", name) == 0) {
        v = trd->AddView(name, Form("/cave_1/trd_%s_0/%s", dtag.Data(), def.first), def.second, debug);
        v->SetProjection(eProjectionType::kXdX, 2500, "um");
        v->SetProjection(eProjectionType::kYdY, 15, "mm");
        if (debug) {
          v->SetProjection(eProjectionType::kDebug0, 2500, "um");
          v->SetProjection(eProjectionType::kDebug1, 15, "mm");
        }
        if (fuRecoConfig[kUseMC]) {
          v->SetProjection(eProjectionType::kXdXMC, 2500, "um");
          v->SetProjection(eProjectionType::kYdYMC, 3000, "um");
          v->SetProjection(eProjectionType::kResidualX, 1000, "um");
          v->SetProjection(eProjectionType::kResidualY, 15, "mm");
        }
      }
      else if (strcmp("1Dx", name) == 0) {
        v = trd->AddView(name, Form("/cave_1/trd_%s_0/%s", dtag.Data(), def.first), def.second);
        v->SetProjection(eProjectionType::kYdY, 5, "cm");
        v->SetProjection(eProjectionType::kXdX, 15, "mm");
        if (fuRecoConfig[kUseMC]) {
          v->SetProjection(eProjectionType::kXdXMC, 5, "mm");
          v->SetProjection(eProjectionType::kYdYMC, 15, "mm");
          v->SetProjection(eProjectionType::kResidualX, 3000, "um");
          v->SetProjection(eProjectionType::kResidualY, 15, "mm");
        }
      }
      else if (strcmp("1Dy", name) == 0) {
        v = trd->AddView(name, Form("/cave_1/trd_%s_0/%s", dtag.Data(), def.first), def.second);
        v->SetProjection(eProjectionType::kYdY, 10, "mm");
        v->SetProjection(eProjectionType::kXdX, 5, "cm");
        if (fuRecoConfig[kUseMC]) {
          v->SetProjection(eProjectionType::kXdXMC, 15, "mm");
          v->SetProjection(eProjectionType::kYdYMC, 5, "mm");
          v->SetProjection(eProjectionType::kResidualX, 15, "mm");
          v->SetProjection(eProjectionType::kResidualY, 3000, "um");
        }
      }
      v->SetProjection(eProjectionType::kChdT, 300, "ns");
      v->SetProjection(eProjectionType::kXdXT, 2.5, "mm");
      v->SetProjection(eProjectionType::kYdYT, 2.5, "mm");
    }
  }
  if (setup->IsActive(ECbmModuleId::kTof)) {
    // setup->GetGeoTag(ECbmModuleId::kTof, dtag);
    // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
    dtag          = "v24f_mcbm";
    Detector* tof = AddDetector(ECbmModuleId::kTof);
    vector<int> tofSelect(3);
    // add type 0
    tofSelect[1] = 0;
    for (int ism(0); ism < 5; ism++) {
      tofSelect[0] = ism;
      for (int irpc(0); irpc < 5; irpc++) {
        tofSelect[2] = irpc;
        v            = tof->AddView(Form("Sm%dRpc%d", ism, irpc),
                         Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/"
                              "gas_box_0/counter_%d",
                              dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
                         tofSelect);
        v->SetProjection(eProjectionType::kXdX, 30, "mm");
        if (fuRecoConfig[kUseMC]) v->SetProjection(eProjectionType::kXdXMC, 30, "mm");
        v->SetProjection(eProjectionType::kYdY, 30, "mm");
        if (fuRecoConfig[kUseMC]) v->SetProjection(eProjectionType::kYdYMC, 30, "mm");
      }
    }
  }
 
  // ===============================================================================
  // TRG - upstream projections
  float angle = 25., sina = TMath::Sin(angle * TMath::DegToRad()), cosa = TMath::Cos(angle * TMath::DegToRad());
  std::vector<float> L = {14.3, 0, -20, -38, -50.5};
  for (auto l : L)
    fPrjPlanes.emplace_back(l * sina, 0., l * cosa);
}
//_____________________________________________________________________
void CbmRecoQaTask::InitDefault()
{
  LOG(info) << "Init Default ....";
  CbmSetup* setup = CbmSetup::Instance();
  if (!setup) {
    LOG(fatal) << GetName() << "::InitDefault() : Missing setup definition.";
    return;
  }
 
  View* v = nullptr;
  std::vector<TString> path;
  TGeoNode* topNode = gGeoManager->GetTopNode();
  if (!topNode) {
    LOG(error) << "Error: Top node not found.";
    return;
  }
 
  TString geoTag;
  auto processDetector = [&](const std::string& detector, const std::string& component) {
    geoTag = GetGeoTagForDetector(detector);
 
    if (geoTag.Length() > 0) {
      LOG(info) << detector << ": geometry tag is " << geoTag;
    }
    else {
      LOG(warn) << "Warning: No geometry tag found for detector " << detector;
    }
    path = GetPath(topNode, detector, component, 0);
  };
 
  if (setup->IsActive(ECbmModuleId::kSts)) {
    processDetector("sts", "Unit");
 
    // std::regex pattern("Unit\\d{2}[LR]_\\d+");
    std::regex pattern(R"(Unit(\d{2}[LR])_(\d+))");
 
    Detector* sts = AddDetector(ECbmModuleId::kSts);
 
    for (const auto& str : path) {
      std::string Str(str.Data());
      std::smatch match;
 
      if (std::regex_search(Str, match, pattern)) {
        int unitid           = std::stoi(match[2].str());
        std::string unitname = (Str.find("Unit") != std::string::npos) ? Str.substr(Str.find("Unit")) : "";
 
        v = sts->AddView(unitname.c_str(), str.Data(), {unitid - 1, -1, -1});
        v->SetSetup(eSetup::kDefault);
      }
      else {
        std::cout << "No match found in string: " << str << std::endl;
      }
    }
  }
 
  if (setup->IsActive(ECbmModuleId::kTrd)) {
    processDetector("trd", "module");
    std::regex pattern("layer(\\d+)_(\\d+)/module(\\d+)_(\\d+)");
    Detector* trd = AddDetector(ECbmModuleId::kTrd);
    char name[256];
    for (const auto& str : path) {
      std::string Str(str.Data());
      std::smatch match;
      if (std::regex_search(Str, match, pattern)) {
        int layer       = std::stoi(match[1]);
        int layercopyNr = std::stoi(match[2]);
 
        int module       = std::stoi(match[3]);
        int modulecopyNr = std::stoi(match[4]);
 
        int fLayer      = ((layercopyNr / 100) % 10);
        int fModuleCopy = (modulecopyNr % 100);
 
        sprintf(name, "layer%d_%d_module_%d_%d", layer, layercopyNr, module, modulecopyNr);
 
        v = trd->AddView(name, str.Data(), {fLayer - 1, fModuleCopy - 1, -1});
        v->SetSetup(eSetup::kDefault);
      }
      else {
        std::cout << "No match found in string: " << str << std::endl;
      }
    }
  }
 
  if (setup->IsActive(ECbmModuleId::kRich)) {
    processDetector("rich", "sens");
    Detector* rich = AddDetector(ECbmModuleId::kRich);
    for (const auto& str : path) {
      v = rich->AddView("rich", str.Data(), {-1, -1, -1});
      v->SetSetup(eSetup::kDefault);
    }
  }
 
  if (setup->IsActive(ECbmModuleId::kTof)) {
    processDetector("tof", "module");
 
    Detector* tof = AddDetector(ECbmModuleId::kTof);
    tof->hit.name = "TofHit";
    std::regex pattern("module_(\\d+)_(\\d+)");
    for (const auto& str : path) {
      std::string Str(str.Data());
      std::smatch match;
      if (std::regex_search(Str, match, pattern)) {
        int type = std::stoi(match[1]);
        int smid = std::stoi(match[2]);
 
        std::string modulename = (Str.find("module") != std::string::npos) ? Str.substr(Str.find("module")) : "";
 
        v = tof->AddView(modulename.c_str(), str.Data(), {smid, type, -1});
        v->SetSetup(eSetup::kDefault);
      }
      else {
        std::cout << "No match found in string: " << str << std::endl;
      }
    }
  }
}
 
//========= DETECTOR ======================
CbmRecoQaTask::Detector::Detector(ECbmModuleId did)
{
  id = did;
  switch (id) {
    case ECbmModuleId::kMvd:
      new (&hit) Data(ECbmDataType::kMvdHit, "MvdHit");
      new (&point) Data(ECbmDataType::kMvdPoint, "MvdPoint");
      break;
    case ECbmModuleId::kSts:
      new (&cl) Data(ECbmDataType::kStsCluster, "StsCluster");
      new (&hit) Data(ECbmDataType::kStsHit, "StsHit");
      new (&point) Data(ECbmDataType::kStsPoint, "StsPoint");
      new (&trk) Data(ECbmDataType::kStsTrack, "StsTrack");
      break;
    case ECbmModuleId::kMuch:
      new (&hit) Data(ECbmDataType::kMuchPixelHit, "MuchHit");
      new (&point) Data(ECbmDataType::kMuchPoint, "MuchPoint");
      break;
    case ECbmModuleId::kRich:
      new (&hit) Data(ECbmDataType::kRichHit, "RichHit");
      new (&point) Data(ECbmDataType::kRichPoint, "RichPoint");
      break;
    case ECbmModuleId::kTrd:
    case ECbmModuleId::kTrd2d:
      new (&cl) Data(ECbmDataType::kTrdCluster, "TrdCluster");
      new (&hit) Data(ECbmDataType::kTrdHit, "TrdHit");
      new (&point) Data(ECbmDataType::kTrdPoint, "TrdPoint");
      new (&trk) Data(ECbmDataType::kTrdTrack, "TrdTrack");
      break;
    case ECbmModuleId::kTof:
      // new (&cl) Data(ECbmDataType::kTofCalDigi, "TofCalDigi");
      new (&hit) Data(ECbmDataType::kTofHit, "TofHit");
      new (&point) Data(ECbmDataType::kTofPoint, "TofPoint");
      new (&trk) Data(ECbmDataType::kTofTrack, "TofTrack");
      break;
    case ECbmModuleId::kFsd:
      new (&hit) Data(ECbmDataType::kFsdHit, "FsdHit");
      new (&point) Data(ECbmDataType::kFsdHit, "FsdPoint");
      break;
    case ECbmModuleId::kPsd:
      new (&hit) Data(ECbmDataType::kPsdHit, "PsdHit");
      new (&point) Data(ECbmDataType::kPsdPoint, "PsdPoint");
      break;
    default:
      LOG(warn) << "QA unsupported for Detector=" << did;
      id = ECbmModuleId::kNotExist;
      break;
  }
}
 
//+++++++++++++++++++++  Detector  ++++++++++++++++++++++++++++++++
//_________________________________________________________________
CbmRecoQaTask::View* CbmRecoQaTask::Detector::AddView(const char* n, const char* p, std::vector<int> set, bool debug)
{
  View* v(nullptr);
  fViews.emplace_back(n, p, set);
  v = &fViews.back();
  v->AddProjection(CbmRecoQaTask::eProjectionType::kChdT);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXYh);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXYhMC);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kDmult);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kDmultMC);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kPullX);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kPullY);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualX);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualY);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualTX);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualTY);
 
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXYa);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXYp);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXdX);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXdXMC);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kYdY);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kYdYMC);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kWdT);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXpX);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kYpY);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kXdXT);
  v->AddProjection(CbmRecoQaTask::eProjectionType::kYdYT);
 
  if (debug) {
    v->AddProjection(CbmRecoQaTask::eProjectionType::kDebug0);
    v->AddProjection(CbmRecoQaTask::eProjectionType::kDebug1);
  }
  return v;
}
 
CbmRecoQaTask::View* CbmRecoQaTask::Detector::GetView(const char* n)
{
  for (auto& view : fViews)
    if (view.name.compare(n) == 0) return &view;
 
  return nullptr;
}
 
CbmRecoQaTask::View* CbmRecoQaTask::Detector::FindView(double x, double y, double z)
{
  for (auto& v : fViews) {
    // printf("Looking for p[%f %f %f] in V[%s-%s] @ %f %f %f (%f %f %f)\n", x,
    // y, z, ToString(id).data(), v.name.data(), v.pos[0], v.pos[1], v.pos[2],
    // v.size[0], v.size[1], v.size[2]);
    if (abs(v.pos[0] - x) > 0.5 * v.size[0]) continue;
    if (abs(v.pos[1] - y) > 0.5 * v.size[1]) continue;
    if (abs(v.pos[2] - z) > 0.5 * v.size[2]) continue;
    return &v;
  }
 
  return nullptr;
}
 
bool CbmRecoQaTask::Detector::Init(TDirectoryFile* fOut, bool mc)
{
  /* Query the geometry and trigger detector views init
   */
  if (!gGeoManager) {
    LOG(fatal) << "CbmRecoQaTask::Detector::Init() " << id << " missing geometry.";
    return false;
  }
  LOG(info) << "CbmRecoQaTask::Detector::Init() : " << cbm::util::ToString(id) << " MC = " << mc;
 
  TDirectoryFile* modDir = (TDirectoryFile*) fOut->mkdir(
    std::string(cbm::util::ToString(id)).c_str(), Form("Reco QA for %s", std::string(cbm::util::ToString(id)).c_str()));
  bool ret(true);
  for (auto& view : fViews) {
    bool vret = view.Init(std::string(cbm::util::ToString(id)).c_str(), mc);
    view.Register(modDir);
    ret = ret && vret;
  }
  return ret;
}
 
void CbmRecoQaTask::Detector::Print() const
{
  stringstream s;
  s << "D[" << id << "] views[" << fViews.size() << "]\n";
  for (auto v : fViews)
    s << v.ToString();
  cout << s.str();
}
 
//+++++++++++++++++++++  View  ++++++++++++++++++++++++++++++++
//_______________________________________________________________________
bool CbmRecoQaTask::View::AddProjection(eProjectionType prj, float range, const char* unit)
{
  if (fProjection.find(prj) != fProjection.end()) {
    LOG(warn) << "Projection " << ToString(prj) << " already registered";
    return false;
  }
  int scale(1);
  if (strcmp(unit, "mm") == 0)
    scale = 10;
  else if (strcmp(unit, "um") == 0)
    scale = 10000;
  else if (strcmp(unit, "ps") == 0)
    scale = 1000;
  else if (strcmp(unit, "eV") == 0)
    scale = 1000;
  else if (strcmp(unit, "cm") == 0)
    scale = 1;
  else if (strcmp(unit, "ns") == 0)
    scale = 1;
  else if (strcmp(unit, "keV") == 0)
    scale = 1;
  else if (strcmp(unit, "a.u.") == 0)
    scale = 1;
  else
    LOG(warn) << "Projection units " << unit << " not registered. Natural units will be used.";
 
  fProjection[prj] = make_tuple(scale, range, nullptr);
  return true;
}
 
bool CbmRecoQaTask::View::HasDebug() const
{
  if (fProjection.find(eProjectionType::kDebug0) != fProjection.end()) return true;
  if (fProjection.find(eProjectionType::kDebug1) != fProjection.end()) return true;
 
  return false;
}
 
bool CbmRecoQaTask::View::SetProjection(eProjectionType prj, float range, const char* unit)
{
  if (fProjection.find(prj) == fProjection.end()) {
    LOG(warn) << "Projection " << ToString(prj)
              << " not initialized. Calling "
                 "\"CbmRecoQaTask::View::AddProjection()\"";
    return AddProjection(prj, range, unit);
  }
  int scale(1);
  if (strcmp(unit, "mm") == 0)
    scale = 10;
  else if (strcmp(unit, "um") == 0)
    scale = 10000;
  else if (strcmp(unit, "ps") == 0)
    scale = 1000;
  else if (strcmp(unit, "eV") == 0)
    scale = 1000;
  else if (strcmp(unit, "cm") == 0)
    scale = 1;
  else if (strcmp(unit, "ns") == 0)
    scale = 1;
  else if (strcmp(unit, "keV") == 0)
    scale = 1;
  else if (strcmp(unit, "a.u.") == 0)
    scale = 1;
  else
    LOG(warn) << "Projection units " << unit << " not registered. Natural units will be used.";
 
  get<0>(fProjection[prj]) = scale;
  get<1>(fProjection[prj]) = range;
  return true;
}
 
bool CbmRecoQaTask::View::Init(const char* dname, bool mc)
{
  if (path.size() && fType == eViewType::kDetUnit) {  // applies only for detection units
    if (!gGeoManager->cd(path.data())) return false;
    TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
    const double* tr(m->GetTranslation());
    TGeoVolume* v = gGeoManager->GetCurrentVolume();
    TGeoShape* bb = v->GetShape();
    double w_lo, w_hi;
    for (int i(0); i < 3; i++) {  // loop over the spatial dimensions
      size[i] = bb->GetAxisRange(i + 1, w_lo, w_hi);
      pos[i]  = 0.5 * (w_lo + w_hi) + tr[i];
    }
 
    LOG(info) << "CbmRecoQaTask::Detector(" << dname << ")::View(" << name << ")::Init() : size [" << size[0] << " x "
              << size[1] << " x " << size[2] << "]. mc[" << mc << "].";
  }
  else
    LOG(info) << "CbmRecoQaTask::Detector(" << dname << ")::View(" << name << ")::Init() :  mc[" << mc << "].";
 
 
  // TODO short-hand for XY display resolution
  int dscale = 10;
  if (strcmp(dname, "Tof") == 0 || strcmp(dname, "Rich") == 0) dscale = 1;
 
  int nbinsx, nbinsy;
  string unit;
  for (auto& projection : fProjection) {
    int scale     = get<0>(projection.second);
    float yrange  = get<1>(projection.second);
    char xy_id    = 0;
    char mc_id[2] = {' ', ' '};
    double xlo = pos[0] - 0.5 * size[0], xhi = pos[0] + 0.5 * size[0], ylo = pos[1] - 0.5 * size[1],
           yhi = pos[1] + 0.5 * size[1];
    switch (projection.first) {
      case eProjectionType::kXdXMC:
        if (!mc) break;
        xy_id    = 'M';
        mc_id[0] = 'M';
        mc_id[1] = 'C';
      // fall through
      case eProjectionType::kXdX:
        if (!xy_id) xy_id = 't';
        nbinsx = 10 * ceil(size[0]);  // mm binning
        if (strcmp(dname, "Trd") == 0 && name.compare("2D") == 0) nbinsx *= 10;
        unit   = makeYrange(scale, yrange);
        nbinsy = 200;  //* ceil(yrange);
        get<2>(projection.second) =
          new TH2D(Form("hxx%s_%s%s", (xy_id == 'M' ? "MC" : ""), dname, name.data()),
                   Form("X resolution %s %s [%s]; X_{%s-%s} (cm); X_{%s} - X_{Trk} (%s)", name.data(), dname, mc_id,
                        dname, name.data(), (xy_id == 'M' ? "Pnt" : "Hit"), unit.data()),
                   nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kYdYMC:
        if (!mc) break;
        xy_id    = 'M';
        mc_id[0] = 'M';
        mc_id[1] = 'C';
      // fall through
      case eProjectionType::kYdY:
        if (!xy_id) xy_id = 't';
        nbinsx = 10 * ceil(size[1]);  // mm binning
        if (strcmp(dname, "Trd") == 0 && name.compare("2D") == 0) nbinsx *= 10;
        unit   = makeYrange(scale, yrange);
        nbinsy = 200;  //* ceil(yrange);
        get<2>(projection.second) =
          new TH2D(Form("hyy%s_%s%s", (xy_id == 'M' ? "MC" : ""), dname, name.data()),
                   Form("Y resolution %s %s [%s]; Y_{%s-%s} (cm); Y_{%s} - Y_{Trk} (%s)", name.data(), dname, mc_id,
                        dname, name.data(), (xy_id == 'M' ? "Pnt" : "Hit"), unit.data()),
                   nbinsx, ylo, yhi, nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kXdXT:
        nbinsx = 10 * ceil(size[0]);  // mm binning
        nbinsy = 200;                 //* ceil(yrange);
        unit   = makeYrange(scale, yrange);
        get<2>(projection.second) =
          new TH2D(Form("hdxT_%s%s", dname, name.data()),
                   Form("X_{TRK} resolution %s [%s]; X^{%s}_{%s-%s} (cm); #sigma_{Trk}(X) (%s)", name.data(), dname,
                        mc_id, dname, name.data(), unit.data()),
                   nbinsx, xlo, xhi, nbinsy, 0., yrange);
        break;
 
      case eProjectionType::kXpX:
        nbinsx = 10 * ceil(size[0]);  // mm binning
        if (yrange < 0) {
          LOG(debug) << "ProjectionP[" << name << "] using default range.";
          yrange = 5;  // +- 5 sigma range
        }
        nbinsy = 200;  //* ceil(yrange);
        get<2>(projection.second) =
          new TH2D(Form("hpx_%s%s", dname, name.data()),
                   Form("X pulls %s [%s]; X^{%s}_{%s-%s} (cm); pull(X)", name.data(), dname, mc_id, dname, name.data()),
                   nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kYdYT:
        nbinsx = 10 * ceil(size[0]);  // mm binning
        nbinsy = 200;                 //* ceil(yrange);
        unit   = makeYrange(scale, yrange);
        get<2>(projection.second) =
          new TH2D(Form("hdyT_%s%s", dname, name.data()),
                   Form("Y_{TRK} resolution %s [%s]; Y^{%s}_{%s-%s} (cm); #sigma_{Trk}(Y) (%s)", name.data(), dname,
                        mc_id, dname, name.data(), unit.data()),
                   nbinsx, xlo, xhi, nbinsy, 0., yrange);
        break;
 
      case eProjectionType::kYpY:
        nbinsx = 10 * ceil(size[1]);  // mm binning
        if (yrange < 0) {
          LOG(debug) << "ProjectionP[" << name << "] using default range.";
          yrange = 5;  // +- 5 sigma range
        }
        nbinsy = 200;  //* ceil(yrange);
        get<2>(projection.second) =
          new TH2D(Form("hpy_%s%s", dname, name.data()),
                   Form("Y pulls %s [%s]; Y^{%s}_{%s-%s} (cm); pull(Y)", name.data(), dname, mc_id, dname, name.data()),
                   nbinsx, ylo, yhi, nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kWdT:
        nbinsx                    = 10 * ceil(size[0]);  // mm binning
        unit                      = makeTrange(scale, yrange);
        nbinsy                    = 2 * ceil(yrange);
        get<2>(projection.second) = new TH2D(Form("hxt_%s%s", dname, name.data()),
                                             Form("Hit_Trk %s [%s]; X_{%s-%s} (cm); T_{Hit} - T_{Trk} (%s)",
                                                  name.data(), dname, dname, name.data(), unit.data()),
                                             nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kChdT:
        nbinsx                    = dscale * ceil(size[0]);  // mm binning
        unit                      = makeTrange(scale, yrange);
        nbinsy                    = 10 * ceil(yrange);
        get<2>(projection.second) = new TH2D(Form("hct_%s%s", dname, name.data()),
                                             Form("Hit Event %s [%s]; X_{%s-%s} (cm); T_{Hit} - T_{Evt} (%s)",
                                                  name.data(), dname, dname, name.data(), unit.data()),
                                             nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
        break;
      case eProjectionType::kDmultMC:
        if (!mc) break;
        xy_id = 'm';
        // fall through
      case eProjectionType::kDmult:
        nbinsx                    = 100;
        xlo                       = -0.5;
        xhi                       = xlo + nbinsx;
        get<2>(projection.second) = new TH2D(
          Form("hdm%s_%s%s", (xy_id == 'm' ? "MC" : ""), dname, name.data()),
          Form("%s multiplicity [EbyE] %s [%s]; N_{%s}^{%s}; N_{%s-%s}^{%s} (%%)", (xy_id ? "point" : "hit"),
               name.data(), dname, dname, (xy_id ? "point" : "hit"), dname, name.data(), (xy_id ? "point" : "hit")),
          nbinsx, xlo, xhi, 200, 0, 1);
        break;
 
      case eProjectionType::kXYa:
        xy_id = 'a';
        // fall through
      case eProjectionType::kXYh:
        if (!xy_id) xy_id = 'h';
        nbinsx = dscale * ceil(size[0]);  // mm binning
        nbinsy = dscale * ceil(size[1]);
        get<2>(projection.second) =
          new TH2D(Form("hxy%c_%s%s", xy_id, dname, name.data()),
                   Form("Hit_{%s} %s [%s]; X_{%s-%s} (cm);  Y_{%s-%s} (cm)", (xy_id == 'h' ? "reco" : "attach"),
                        name.data(), dname, dname, name.data(), dname, name.data()),
                   nbinsx, xlo, xhi, nbinsy, ylo, yhi);
        break;
 
      case eProjectionType::kXYhMC:
        if (!mc) break;
        nbinsx   = dscale * ceil(size[0]);  // mm binning
        nbinsy   = dscale * ceil(size[1]);
        get<2>(projection.second) = new TH2D(Form("hxymc_%s%s", dname, name.data()),
                                             Form("Point %s [%s]; X^{MC}_{%s-%s} (cm);  Y^{MC}_{%s-%s} (cm)",
                                                  name.data(), dname, dname, name.data(), dname, name.data()),
                                             nbinsx, xlo, xhi, nbinsy, ylo, yhi);
        break;
 
      case eProjectionType::kResidualX:
      // fall through
      case eProjectionType::kResidualY:
        if (!mc) break;
        xy_id  = (projection.first == eProjectionType::kResidualX ? 'X' : 'Y');
        nbinsx = dscale * ceil(size[0]);  // mm binning
        if (strcmp(dname, "Trd") == 0 && name.compare("2D") == 0) nbinsx *= 10;
        nbinsy = 200;  //dscale * ceil(size[1]);
        unit   = makeYrange(scale, yrange);
        get<2>(projection.second) =
          new TH2D(Form("h%c%cR_%s%s", xy_id, xy_id, dname, name.data()),
                   Form("%c resolution %s [%s]; %c_{%s-%s} (cm); %c_{Hit} - %c_{Pnt} (%s)", xy_id, name.data(), dname,
                        xy_id, dname, name.data(), xy_id, xy_id, unit.data()),
                   nbinsx, (xy_id == 'X' ? xlo : ylo), (xy_id == 'X' ? xhi : yhi), nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kResidualTX:
      // fall through
      case eProjectionType::kResidualTY:
        if (!mc) break;
        xy_id  = (projection.first == eProjectionType::kResidualTX ? 'X' : 'Y');
        nbinsx = dscale * ceil(size[0]);  // mm binning
        nbinsy = 200;                     // dscale * ceil(size[1]);
        unit   = makeYrange(scale, yrange);
        get<2>(projection.second) =
          new TH2D(Form("ht%cR_%s%s", xy_id, dname, name.data()),
                   Form("Time resolution %s [%s]; %c_{%s-%s} (cm); T_{Hit} - T_{Pnt} (ns)", name.data(), dname, xy_id,
                        dname, name.data()),
                   nbinsx, (xy_id == 'X' ? xlo : ylo), (xy_id == 'X' ? xhi : yhi), nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kPullX:
      // fall through
      case eProjectionType::kPullY:
        if (!mc) break;
        xy_id  = (projection.first == eProjectionType::kPullX ? 'X' : 'Y');
        nbinsx = dscale * ceil(size[0]);  // mm binning
        nbinsy = dscale * ceil(size[1]);
        unit   = makeYrange(scale, yrange);
 
        get<2>(projection.second) =
          new TH2D(Form("h%c%cP_%s%s", xy_id, xy_id, dname, name.data()),
                   Form("%c pull %s [%s]; %c_{%s-%s} (cm); pull(%c)", xy_id, name.data(), dname, xy_id, dname,
                        name.data(), xy_id),
                   nbinsx, (xy_id == 'X' ? xlo : ylo), (xy_id == 'X' ? xhi : yhi), nbinsy, -yrange, yrange);
        break;
 
      case eProjectionType::kXYp:
        nbinsx = dscale * ceil(size[0] * 2.);  // mm binning
        nbinsy = dscale * ceil(size[1] * 2.);
        xlo    = pos[0] - size[0];
        xhi    = pos[0] + size[0];
        ylo    = pos[1] - size[1];
        yhi    = pos[1] + size[1];
 
        get<2>(projection.second) = new TH2D(Form("hxyp_%s%s", dname, name.data()),
                                             Form("Trk_{proj} %s [%s]; X_{%s-%s} (cm);  Y_{%s-%s} (cm)", name.data(),
                                                  dname, dname, name.data(), dname, name.data()),
                                             nbinsx, xlo, xhi, nbinsy, ylo, yhi);
        break;
      case eProjectionType::kXYt0: xy_id = '0';
      // fall through
      case eProjectionType::kXYt1:
        if (!xy_id) xy_id = '1';
      // fall through
      case eProjectionType::kXYt2:
        if (!xy_id) xy_id = '2';
      // fall through
      case eProjectionType::kXYt3:
        if (!xy_id) xy_id = '3';
      // fall through
      case eProjectionType::kXYt4:
        if (!xy_id) xy_id = '4';
      // fall through
      case eProjectionType::kXYt5:
        if (!xy_id) xy_id = '5';
        nbinsx = 4500;  // mm binning
        nbinsy = 1000;
        xlo    = -25;
        xhi    = +20;
        ylo    = -5;
        yhi    = +5;
 
        get<2>(projection.second) = new TH2D(Form("hxyt%c_%s%s", xy_id, dname, name.data()),
                                             Form("Trk_{proj} z = %+.2f [%s]; X_{%s-%s} (cm);  Y_{%s-%s} (cm)", yrange,
                                                  dname, dname, name.data(), dname, name.data()),
                                             nbinsx, xlo, xhi, nbinsy, ylo, yhi);
        break;
 
      // Primary vertex plots
      case eProjectionType::kPVmult:  // define primary vertex multiplicity
        nbinsx = 50;
        nbinsy = 20;
        xlo    = -0.5;
        xhi    = xlo + nbinsx;
        ylo    = -0.5;
        yhi    = ylo + nbinsy;
 
        get<2>(projection.second) =
          new TH2D(Form("hMult_%s%s", dname, name.data()), "Vertex multiplicity; N_{trk}^{event};  N_{trk}^{PV}",
                   nbinsx, xlo, xhi, nbinsy, ylo, yhi);
        break;
      case eProjectionType::kPVxz:  // define X-Z vertex projection
        xy_id    = 'y';
        mc_id[0] = 'z';  // use to store x-axis title
        mc_id[1] = 'x';  // use to store y-axis title
        nbinsx   = 100;
        nbinsy   = 200;
        xlo      = -10;
        xhi      = +10;
        ylo      = -2;
        yhi      = +2;
        // fall through
      case eProjectionType::kPVyz:  // define Y-Z vertex projection
        if (!xy_id) {
          xy_id    = 'x';
          mc_id[0] = 'z';  // use to store x-axis title
          mc_id[1] = 'y';  // use to store y-axis title
          nbinsx   = 100;
          nbinsy   = 200;
          xlo      = -10;
          xhi      = +10;
          ylo      = -2;
          yhi      = +2;
        }
        // fall through
      case eProjectionType::kPVxy:  // define X-Y vertex projection
        if (!xy_id) {
          xy_id    = 'z';
          mc_id[0] = 'x';  // use to store x-axis title
          mc_id[1] = 'y';  // use to store y-axis title
          nbinsx   = 200;
          nbinsy   = 200;
          xlo      = -2;
          xhi      = +2;
          ylo      = -2;
          yhi      = +2;
        }
        get<2>(projection.second) =
          new TH2D(Form("h%c%c_%s%s", mc_id[0], mc_id[1], dname, name.data()),
                   Form("Vertex_{%c%c}; %c (cm);  %c (cm)", mc_id[0], mc_id[1], mc_id[0], mc_id[1]), nbinsx, xlo, xhi,
                   nbinsy, ylo, yhi);
        break;
 
      case eProjectionType::kPVzsz:  // define z-sz vertex projection
        xy_id = 'z';
        yhi   = +0.5;
        // fall through
      case eProjectionType::kPVysy:  // define y - sy vertex projection
        if (!xy_id) {
          xy_id = 'y';
          yhi   = +0.2;
        }
        // fall through
      case eProjectionType::kPVxsx:  // define x-sx vertex projection
        if (!xy_id) {
          xy_id = 'x';
          yhi   = +0.1;
        }
        nbinsx                    = 200;
        nbinsy                    = 200;
        xlo                       = -0.5;
        xhi                       = +0.5;
        ylo                       = 0;
        get<2>(projection.second) = new TH2D(Form("hs%c_%s%s", xy_id, dname, name.data()),
                                             Form("Vertex_{%c}; %c (cm);  #sigma %c (cm)", xy_id, xy_id, xy_id), nbinsx,
                                             xlo, xhi, nbinsy, ylo, yhi);
        break;
      case eProjectionType::kDebug0:  // expand kXdX:
        xy_id = 'x';
 
        // fall through
      case eProjectionType::kDebug1:  // expand kYdY:
        if (!xy_id) xy_id = 'y';
        nbinsx = 10 * ceil(size[0]);  // mm binning
        if (strcmp(dname, "Trd") == 0 && name.compare("2D") == 0) nbinsx *= 10;
        unit   = makeYrange(scale, yrange);
        nbinsy = 200;  //* ceil(yrange);
        get<2>(projection.second) =
          new TH3D(Form("h%c%cr_%s%s", xy_id, xy_id, dname, name.data()),
                   Form("%c residuals %s [%s]; x_{%s-%s} (cm); %c_{Hit} - %c_{Trk} (%s); row", xy_id, name.data(),
                        dname, dname, name.data(), xy_id, xy_id, unit.data()),
                   nbinsx, xlo, xhi, nbinsy, -yrange, yrange, 20, -0.5, 19.5);
        break;
    }
  }
  return true;
}
 
string CbmRecoQaTask::View::ToString() const
{
  stringstream s;
  s << "V[" << name << "] path[" << path << "] @ [" << pos[0] << ", " << pos[1] << ", " << pos[2] << "] size["
    << size[0] << ", " << size[1] << ", " << size[2] << "] projections[" << fProjection.size() << "] sel[";
  for (auto ii : fSelector)
    s << ii << " ";
  s << "]\n";
  return s.str();
}
 
uint CbmRecoQaTask::View::Register(TDirectoryFile* fOut)
{
  uint n(0);
  TDirectoryFile* lDir = (TDirectoryFile*) fOut->mkdir(name.data(), Form("QA for vol[%s]", path.data()));
  TDirectoryFile *sDir(nullptr), *tDir(nullptr);
  // run extra directory structure for the case of detector view
  if (fType == eViewType::kDetUnit) {
    if (CbmRecoQaTask::fuRecoConfig[kRecoEvents])
      sDir = (TDirectoryFile*) lDir->mkdir("event", "Local Reco QA");
    else
      sDir = (TDirectoryFile*) lDir->mkdir("TS", "TimeSlice QA");
    tDir = (CbmRecoQaTask::fuRecoConfig[kRecoTracks] ? (TDirectoryFile*) lDir->mkdir("track", "CA QA") : nullptr);
  }
 
  for (auto& projection : fProjection) {
    TH1* hh = get<2>(projection.second);
    if (!hh) continue;
    switch (projection.first) {
      case eProjectionType::kXYa:
      case eProjectionType::kXYp:
      case eProjectionType::kXdX:
      case eProjectionType::kXdXMC:
      case eProjectionType::kYdY:
      case eProjectionType::kYdYMC:
      case eProjectionType::kWdT:
      case eProjectionType::kXpX:
      case eProjectionType::kYpY:
      case eProjectionType::kXdXT:
      case eProjectionType::kYdYT:
      case eProjectionType::kDebug0:
      case eProjectionType::kDebug1:
        if (tDir) tDir->Add(hh);
        break;
      case eProjectionType::kChdT:
      case eProjectionType::kXYh:
      case eProjectionType::kXYhMC:
      case eProjectionType::kDmult:
      case eProjectionType::kDmultMC:
      case eProjectionType::kResidualX:
      case eProjectionType::kResidualY:
      case eProjectionType::kResidualTX:
      case eProjectionType::kResidualTY:
      case eProjectionType::kPullX:
      case eProjectionType::kPullY:
        if (sDir) sDir->Add(hh);
        break;
      case eProjectionType::kXYt0:
      case eProjectionType::kXYt1:
      case eProjectionType::kXYt2:
      case eProjectionType::kXYt3:
      case eProjectionType::kXYt4:
      case eProjectionType::kXYt5:
      case eProjectionType::kPVxz:
      case eProjectionType::kPVyz:
      case eProjectionType::kPVxy:
      case eProjectionType::kPVxsx:
      case eProjectionType::kPVysy:
      case eProjectionType::kPVzsz:
      case eProjectionType::kPVmult:
        if (lDir) lDir->Add(hh);
        break;
    }
    n++;
  }
  LOG(debug) << "CbmRecoQaTask::View[" << name << "]::Register() : " << n << " projections for " << fOut->GetName();
  return n;
}
 
string CbmRecoQaTask::View::ToString(eProjectionType prj)
{
  string s;
  switch (prj) {
    case eProjectionType::kXYa: s = "x-y (attach)"; break;
    case eProjectionType::kXYp:
    case eProjectionType::kXYt0:
    case eProjectionType::kXYt1:
    case eProjectionType::kXYt2:
    case eProjectionType::kXYt3:
    case eProjectionType::kXYt4:
    case eProjectionType::kXYt5: s = "x-y (track)"; break;
    case eProjectionType::kXdX: s = "dx-x"; break;
    case eProjectionType::kXdXMC: s = "dx-x (MC)"; break;
    case eProjectionType::kResidualX: s = "residual(x)-x:MC"; break;
    case eProjectionType::kResidualY: s = "residual(y)-y:MC"; break;
    case eProjectionType::kResidualTX: s = "residual(t)-x:MC"; break;
    case eProjectionType::kResidualTY: s = "residual(t)-y:MC"; break;
    case eProjectionType::kPullX: s = "pull(x)-x:MC"; break;
    case eProjectionType::kPullY: s = "pull(y)-y:MC"; break;
    case eProjectionType::kYdY: s = "dy-y"; break;
    case eProjectionType::kYdYMC: s = "dy-y (MC)"; break;
    case eProjectionType::kXpX: s = "pull(x)-x"; break;
    case eProjectionType::kYpY: s = "pull(y)-y"; break;
    case eProjectionType::kXdXT: s = "d_trk(x)-x"; break;
    case eProjectionType::kYdYT: s = "d_trk(y)-y"; break;
    case eProjectionType::kWdT: s = "dt(trk)-w"; break;
    case eProjectionType::kChdT: s = "dt(ev)-w"; break;
    case eProjectionType::kXYh: s = "x-y (hit)"; break;
    case eProjectionType::kXYhMC: s = "x-y (point)"; break;
    default: LOG(error) << "View::ToString() : Unknown projection " << int(prj); break;
  }
  return s;
}
 
string CbmRecoQaTask::View::makeYrange(const int scale, float& yrange)
{
  bool kDefaultRange = false;
  string unit        = "cm";
  if (yrange < 0) {
    LOG(debug) << "ProjectionY[" << name << "] using default range.";
    kDefaultRange = true;
    yrange        = 3;  // +- 6 cm default range
  }
  if (scale == 10) {  // mm resolution
    unit = "mm";
    if (kDefaultRange) yrange = 10;  // +- 20 mm default range
  }
  else if (scale == 10000) {  // micron resolution
    unit = "#mu m";
    if (kDefaultRange) yrange = 150;  // +- 300 um default range
  }
  return unit;
}
 
string CbmRecoQaTask::View::makeTrange(const int scale, float& yrange)
{
  bool kDefaultRange = false;
  string unit        = "ns";
  if (yrange < 0) {
    LOG(debug) << "ProjectionT[" << name << "] using default range.";
    kDefaultRange = true;
    yrange        = 40;  // +- 80 ns default range
  }
  if (scale == 1000) {  // ps resolution
    unit = "ps";
    if (kDefaultRange) yrange = 300;  // +- 600 ps default range
  }
  return unit;
}
 
CbmRecoQaTask::EventFilter* CbmRecoQaTask::AddEventFilter(CbmRecoQaTask::EventFilter::eEventCut type)
{
  for (auto cut : fFilterEv) {
    if (cut.fType == type) {
      LOG(warning) << GetName() << "::AddEventFilter event filter : " << cut.ToString() << " already on the list.";
      return nullptr;
    }
  }
  fFilterEv.emplace_back(type);
  return &fFilterEv.back();
}
 
//+++++++++++++++++++++  TrackFilter  ++++++++++++++++++++++++++++++++
//____________________________________________________________________
CbmRecoQaTask::TrackFilter* CbmRecoQaTask::AddTrackFilter(ECbmModuleId type)
{
  // check if not already registered
  for (auto cut : fFilterTrk) {
    if (cut.fDet == type) {
      LOG(warning) << GetName() << "::AddTrackFilter track filter : " << type << " already on the list.";
      return nullptr;
    }
  }
  fFilterTrk.emplace_back(type);
  return &fFilterTrk.back();
}
std::string CbmRecoQaTask::TrackFilter::ToString() const
{
  stringstream ss;
 
  if (fDet == ECbmModuleId::kRef) {  // take care of projection cuts
    for (auto cut : fProjFilter.fRef) {
      ss << "::TrackProjection(z = " << cut[0] << ") centered on [" << cut[1] << ", " << cut[2] << "] < " << cut[3]
         << "cm.";
    }
  }
  else {
    ss << "::TrackHitFilter(" << fDet << ") min hits[" << fNminHits << "]. ";
    auto detTopo = [&](vector<int> cut) {
      std::string label[3];
      switch (fDet) {
        case ECbmModuleId::kSts: ss << "U" << cut[0] << "L" << cut[1] << "M" << cut[2] << " / "; break;
        case ECbmModuleId::kTrd:
          switch (cut[3]) {
            case 0: label[0] = "n"; break;
            case 1: label[0] = "y"; break;
            default: label[0] = "any"; break;
          }
          switch (cut[4]) {
            case 0: label[1] = "rect"; break;
            case 1: label[1] = "tilt"; break;
            default: label[1] = "any"; break;
          }
          switch (cut[5]) {
            case 0: label[2] = "left"; break;
            case 1: label[2] = "center"; break;
            case 2: label[2] = "right"; break;
            default: label[2] = "any"; break;
          }
          ss << "Ly" << cut[0] << " Mod" << cut[1] << " ClSz[" << cut[2] << "] RowCross[" << label[0] << "] Max["
             << label[1] << "] Symm[" << label[2] << "] PadRow[" << cut[6] << "] / ";
          break;
        default: break;
      }
    };
    if (fHitsFilter.fAccept.size()) {
      ss << "accept topology : ";
      for (auto cut : fHitsFilter.fAccept)
        detTopo(cut);
    }
    if (fHitsFilter.fReject.size()) {
      ss << "reject topology : ";
      for (auto cut : fHitsFilter.fReject)
        detTopo(cut);
    }
  }
 
  return ss.str();
}
 
//+++++++++++++++++++++  EventFilter  ++++++++++++++++++++++++++++++++
//____________________________________________________________________
bool CbmRecoQaTask::EventFilter::Accept(const CbmEvent* ptr, const CbmRecoQaTask* /* lnk*/)
{
  stringstream ss;
  bool ret(true);
  size_t val(0);
  switch (fType) {
    case eEventCut::kMultTrk:
      val = ptr->GetNofData(ECbmDataType::kGlobalTrack);
      if (fMinTrack > 0 && val < size_t(fMinTrack)) {
        ss << "NofTrack[" << val << "] < min[" << fMinTrack << "].";
        ret = false;
        break;
      }
      if (ret && fMaxTrack > 0 && val > size_t(fMaxTrack)) {
        ss << "NofTrack[" << val << "] > max[" << fMaxTrack << "].";
        ret = false;
      }
      break;
    case eEventCut::kMultHit:
      val = ptr->GetNofData(ECbmDataType::kStsHit);
      if (fMultHit[0] > 0 && val > size_t(fMultHit[0])) {
        ss << "Sts hits [" << val << "] > max[" << fMultHit[0] << "].";
        ret = false;
        break;
      }
      val = ptr->GetNofData(ECbmDataType::kTrdHit);
      if (fMultHit[1] > 0 && val > size_t(fMultHit[1])) {
        ss << "Trd hits [" << val << "] > max[" << fMultHit[1] << "].";
        ret = false;
        break;
      }
      val = ptr->GetNofData(ECbmDataType::kTofHit);
      if (fMultHit[2] > 0 && val > size_t(fMultHit[2])) {
        ss << "Tof hits [" << val << "] > max[" << fMultHit[2] << "].";
        ret = false;
      }
      break;
    default: break;
  }
  if (!ret) LOG(debug2) << "EventFilter::Reject event for : " << ss.str();
  return ret;
}
bool CbmRecoQaTask::EventFilter::SetFilter(std::vector<float> cuts)
{
  switch (fType) {
    case eEventCut::kMultTrk:
      if (cuts.size() < 2 || cuts[1] < cuts[0]) {
        LOG(warning) << "Improper definition for event filter :\n\t" << ToString() << endl;
        HelpMess();
        return false;
      }
      fMinTrack = int(cuts[0]);
      fMaxTrack = int(cuts[1]);
      break;
    case eEventCut::kMultHit:
      if (cuts.size() < 3) {
        LOG(warning) << "Improper definition for event filter :\n\t" << ToString() << endl;
        HelpMess();
        return false;
      }
      for (int i(0); i < int(CbmRecoQaTask::EventFilter::eEventDef::kNofDetHit); i++)
        fMultHit[i] = int(cuts[i]);
      break;
    case eEventCut::kTrigger: break;
    case eEventCut::kVertex: break;
    default: break;
  }
  return true;
}
std::string CbmRecoQaTask::EventFilter::ToString() const
{
  stringstream ss;
  switch (fType) {
    case eEventCut::kMultTrk: ss << "kMultTrk : \"cut on track multiplicity"; break;
    case eEventCut::kMultHit: ss << "kMultHit : \"cut on hit multiplicity"; break;
    case eEventCut::kTrigger: ss << "kTrigger : \"cut on trigger conditions"; break;
    case eEventCut::kVertex: ss << "kVertex : \"cut on vertex definition"; break;
    default: break;
  }
  return ss.str();
}
void CbmRecoQaTask::EventFilter::HelpMess() const
{
  LOG(info) << "CbmRecoQaTask::EventFilter : Usage";
  switch (fType) {
    case eEventCut::kMultTrk:
      LOG(info) << ToString();
      LOG(info) << "\tDepends one two values : min and max of the no of global "
                   "tracks / event.";
      LOG(info) << "\tValue should follow the relation max >= min.";
      break;
    case eEventCut::kMultHit:
      LOG(info) << ToString();
      LOG(info) << "\tDepends one three values : total hits in the following "
                   "systems (in this order) : STS TRD TOF.";
      LOG(info) << "\tNegative values are interpreted as no-cut.";
      break;
    default:
      LOG(info) << ToString();
      LOG(info) << "\tNo user info.";
      break;
  }
}
//+++++++++++++++++++++  TrackFilter  ++++++++++++++++++++++++++++++++
//____________________________________________________________________
bool CbmRecoQaTask::TrackFilter::Accept(const CbmGlobalTrack* ptr, const CbmRecoQaTask* lnk) const
{
  bool ret(true);
  stringstream ss;
  int idx(-1);
  const CbmTrack* trk(nullptr);
  // check if we have any registered cut for detector
  if (fNminHits <= 0) return false;
  // check if we have any found track for detector
  switch (fDet) {
    case ECbmModuleId::kSts: idx = ptr->GetStsTrackIndex(); break;
    case ECbmModuleId::kTrd: idx = ptr->GetTrdTrackIndex(); break;
    case ECbmModuleId::kTof: idx = ptr->GetTofTrackIndex(); break;
    case ECbmModuleId::kMuch: idx = ptr->GetMuchTrackIndex(); break;
    default: break;
  }
  if (idx < 0) {
    ss << fDet << " trk missing";
    ret = false;
  }
  if (ret && !(trk = lnk->GetData<CbmTrack>(fDet, idx))) {
    ss << fDet << " trk[" << idx << "]=nullptr.";
    ret = false;
  }
  // cut on the total no of STS hits/trk
  if (ret && trk->GetNofHits() < fNminHits) {
    ss << fDet << " hits/trk[" << idx << "]=" << trk->GetNofHits() << " < min[" << fNminHits << "].";
    ret = false;
  }
  // check topology of hits on track if topological cuts are defined for detector
  if (ret && (fHitsFilter.fAccept.size() || fHitsFilter.fReject.size())) {
    for (int ihit(0), jhit(0); ihit < trk->GetNofHits(); ihit++) {
      jhit     = trk->GetHitIndex(ihit);
      auto hit = (CbmHit*) lnk->GetData<CbmHit>(fDet, jhit);
      if (!hit) {
        LOG(warning) << "TrackFilter::Accept : Missing " << ihit << "hit @ idx=" << jhit << "  for " << fDet;
        continue;
      }
      //printf("check hit %d\n", ihit);
      if (!(ret = lnk->FilterHit(hit, fDet))) break;
    }
  }
 
  if (!ret) LOG(debug) << "TrackFilter::Reject trk for : " << ss.str();
  return ret;
}
 
//____________________________________________________________________
bool CbmRecoQaTask::TrackFilter::AddTopoCut(std::vector<int> cuts, const char opt)
{
  size_t ncuts(0);
  switch (fDet) {
    case ECbmModuleId::kSts: ncuts = stsSz; break;
    case ECbmModuleId::kTof: ncuts = tofSz; break;
    case ECbmModuleId::kTrd: ncuts = trdSz; break;
    default:
      LOG(warning) << "CbmRecoQaTask::TrackFilter::AddTopoCut : k" << fDet << " : \"specific cut on hit topology for "
                   << fDet << "\" not implemented.";
      return false;
  }
  if (cuts.size() != ncuts) {
    LOG(warning) << "CbmRecoQaTask::TrackFilter::AddTopoCut : k" << fDet << " : \"specific cut on hit topology for "
                 << fDet << "\" incorrect size " << cuts.size() << ". Needed size " << ncuts;
    return false;
  }
  switch (opt) {
    case 'a': fHitsFilter.fAccept.push_back(cuts); break;
    case 'r': fHitsFilter.fReject.push_back(cuts); break;
  }
 
  return true;
}
 
//____________________________________________________________________
bool CbmRecoQaTask::TrackFilter::SetProjFilter(std::vector<float> cuts)
{
  if (cuts.size() < 1) {
    LOG(warning) << "::SetFilter : k" << fDet << " : \"cut on no of " << fDet << " hits / track\" not defined.";
    return false;
  }
  fProjFilter.fRef.emplace_back(std::array<float, projSz>{cuts[0], cuts[1], cuts[2], cuts[3]});
  return true;
}
 
/* clang-format off */
ClassImp(CbmRecoQaTask)
ClassImp(CbmRecoQaTask::Detector)
ClassImp(CbmRecoQaTask::View)
ClassImp(CbmRecoQaTask::TopoFilter)
ClassImp(CbmRecoQaTask::TrackFilter)
ClassImp(CbmRecoQaTask::EventFilter)
  /* clang-format on */