CbmMvdMcQaTask.cxx

Go to the documentation of this file.

/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Ajit Kumar [committer] */

 
#include "CbmMvdMcQaTask.h"
 
#include "CbmMCDataArray.h"
#include "CbmMCDataManager.h"
#include "CbmMCEventList.h"
#include "CbmMCTrack.h"
#include "CbmMvdPoint.h"
 
#include <FairRootManager.h>
 
using cbm::mvd::CbmMvdMcQaTask;
 
CbmMvdMcQaTask::CbmMvdMcQaTask(int verbose, bool isMcUsed) : CbmQaTask("CbmMvdMcQaTask", verbose, isMcUsed) {}
 
// --------------------------------------------------------------------------------------------------------------------- //
InitStatus CbmMvdMcQaTask::InitQa()
{
  fManager = FairRootManager::Instance();
 
  if (!fManager) {
    LOG(error) << "No FairRootManager found";
    return kERROR;
  }
  if (IsMCUsed()) {
    fMcManager = dynamic_cast<CbmMCDataManager*>(fManager->GetObject("MCDataManager"));
    if (!fMcManager) LOG(fatal) << GetName() << ": No MCDataManager!";
 
    fMcEventList = dynamic_cast<CbmMCEventList*>(fManager->GetObject("MCEventList."));
    if (!fMcEventList) {
      fMcEventList = dynamic_cast<CbmMCEventList*>(fManager->GetObject("EventList."));
      if (!fMcEventList)
        LOG(fatal) << GetName()
                   << ": No MCEventList.! Add digitization and reconstruction files in the macro. MCEventList comes in "
                      "digitzation or reconstruction.";
    }
 
    // Get MCTrack array
    fMcTracks = fMcManager->InitBranch("MCTrack");
    if (!fMcTracks) {
      LOG(error) << "No MC tracks array found in CbmMCDataManager!";
      return kERROR;
    }
    // Get MvdPoint array
    fMvdPoints = fMcManager->InitBranch("MvdPoint");
    if (!fMvdPoints) {
      LOG(error) << "No MC points array found in CbmMCDataManager!";
      return kERROR;
    }
  }
  for (int oi = 0; oi < kNumOrigins; ++oi) {
    fLists[oi] = new TList();
    fLists[oi]->SetOwner(kTRUE);
    fLists[oi]->SetName(Form("MvdMcQa_%s", kOriginName[oi]));
    fManager->Register(Form("MvdMcQaHists_%s", kOriginName[oi]), "MVD", fLists[oi], kTRUE);
  }
 
  BookHists_();
 
  return kSUCCESS;
}
// ----- Book histograms ----------------------------------------------------- //
void CbmMvdMcQaTask::BookHists_()
{
  // MCTrack level histograms
  TString name  = Form("mc_tracks_per_event");
  TString title = Form("MC tracks per event;# MC tracks;Counts");
  fh_nTracksAll = MakeQaObject<TH1D>(name, title, 20000, 0, 20000);
 
  // --- MCTrack-level QA ---
  const int nSpeciesBins = 18;
 
  fh_mc_pdg_all  = MakeQaObject<TH1D>("h_mc_pdg_all", "MCTrack particle species (all);particle;counts", nSpeciesBins,
                                     0.5, nSpeciesBins + 0.5);
  fh_mc_pdg_prim = MakeQaObject<TH1D>("h_mc_pdg_prim", "MCTrack particle species (primaries);particle;counts",
                                      nSpeciesBins, 0.5, nSpeciesBins + 0.5);
  fh_mc_pdg_sec  = MakeQaObject<TH1D>("h_mc_pdg_sec", "MCTrack particle species (secondaries);particle;counts",
                                     nSpeciesBins, 0.5, nSpeciesBins + 0.5);
 
  auto labelSpecies = [](TH1* h) {
    auto* ax = h->GetXaxis();
    ax->SetBinLabel(1, "#gamma");
    ax->SetBinLabel(2, "e^{-}");
    ax->SetBinLabel(3, "e^{+}");
    ax->SetBinLabel(4, "#mu^{-}");
    ax->SetBinLabel(5, "#mu^{+}");
    ax->SetBinLabel(6, "#pi^{+}");
    ax->SetBinLabel(7, "#pi^{-}");
    ax->SetBinLabel(8, "#pi^{0}");
    ax->SetBinLabel(9, "K^{+}");
    ax->SetBinLabel(10, "K^{-}");
    ax->SetBinLabel(11, "K^{0}_{S}");
    ax->SetBinLabel(12, "K^{0}_{L}");
    ax->SetBinLabel(13, "p");
    ax->SetBinLabel(14, "#bar{p}");
    ax->SetBinLabel(15, "n");
    ax->SetBinLabel(16, "#bar{n}");
    ax->SetBinLabel(17, "ions");
    ax->SetBinLabel(18, "other");
  };
 
  labelSpecies(fh_mc_pdg_all);
  labelSpecies(fh_mc_pdg_prim);
  labelSpecies(fh_mc_pdg_sec);
 
  // momentum distributions (tune max if you like)
  fh_mc_p_prim = MakeQaObject<TH1D>("h_mc_p_prim", "MCTrack momentum (primaries);p [GeV/c];counts", 200, 0., 10.);
  fh_mc_p_sec  = MakeQaObject<TH1D>("h_mc_p_sec", "MCTrack momentum (secondaries);p [GeV/c];counts", 200, 0., 10.);
 
  // start vertex distributions
  fh_mc_startZ = MakeQaObject<TH1D>("h_mc_startZ", "MCTrack start z;z_{start} [cm];counts", 400, -200., 200.);
 
  fh_mc_startX_vs_Z = MakeQaObject<TH2D>("h_mc_startX_vs_Z", "MCTrack start x vs z;x_{start} [cm];z_{start} [cm]", 400,
                                         -50., 50., 400, -200., 200.);
 
  fh_mc_startY_vs_Z = MakeQaObject<TH2D>("h_mc_startY_vs_Z", "MCTrack start y vs z;y_{start} [cm];z_{start} [cm]", 400,
                                         -50., 50., 400, -200., 200.);
 
  fh_mc_startY_vs_X = MakeQaObject<TH2D>("h_mc_startY_vs_X", "MCTrack start y vs x;x_{start} [cm];y_{start} [cm]", 400,
                                         -50., 50., 400, -50., 50.);
 
  // MvdPoint histograms
  // per-event per-station
  fh_mvd_nPointsEvt_station.resize(fNStations);
  for (int iSt = 0; iSt < fNStations; ++iSt) {
    TString namest                 = Form("h_mvd_nPointsEvt_station_%d", iSt);
    TString titlest                = Form("MVD points per event, station %d;N_{MVD points};Events", iSt);
    fh_mvd_nPointsEvt_station[iSt] = MakeQaObject<TH1D>(namest, titlest, 5000, 0., 5000.);
  }
 
  // XY maps per station, split by origin
  fh_mvd_xy_prim_station.resize(fNStations);
  fh_mvd_xy_sec_station.resize(fNStations);
 
  for (int iSt = 0; iSt < fNStations; ++iSt) {
    TString nameP  = Form("h_mvd_xy_prim_station_%d", iSt);
    TString titleP = Form("MVD primary points X-Y, station %d;x [cm];y [cm]", iSt);
    fh_mvd_xy_prim_station[iSt] =
      MakeQaObject<TH2D>(nameP, titleP, 400, -5., 5., 400, -5., 5.);  // adjust ranges to your detector
 
    TString nameS              = Form("h_mvd_xy_sec_station_%d", iSt);
    TString titleS             = Form("MVD secondary points X-Y, station %d;x [cm];y [cm]", iSt);
    fh_mvd_xy_sec_station[iSt] = MakeQaObject<TH2D>(nameS, titleS, 400, -5., 5., 400, -5., 5.);
  }
 
  fh_mvd_pdg_all = MakeQaObject<TH1D>("h_mvd_pdg_all", "MVD point particle species (all);particle;counts", nSpeciesBins,
                                      0.5, nSpeciesBins + 0.5);
 
  fh_mvd_pdg_prim = MakeQaObject<TH1D>("h_mvd_pdg_prim", "MVD point particle species (primaries);particle;counts",
                                       nSpeciesBins, 0.5, nSpeciesBins + 0.5);
 
  fh_mvd_pdg_sec = MakeQaObject<TH1D>("h_mvd_pdg_sec", "MVD point particle species (secondaries);particle;counts",
                                      nSpeciesBins, 0.5, nSpeciesBins + 0.5);
 
  labelSpecies(fh_mvd_pdg_all);
  labelSpecies(fh_mvd_pdg_prim);
  labelSpecies(fh_mvd_pdg_sec);
 
  // Momentum distributions from parent MCTrack
  fh_mvd_p_prim =
    MakeQaObject<TH1D>("h_mvd_p_prim", "Momentum at MVD point (primaries);p [GeV/c];counts", 200, 0., 10.);
 
  fh_mvd_p_sec =
    MakeQaObject<TH1D>("h_mvd_p_sec", "Momentum at MVD point (secondaries);p [GeV/c];counts", 200, 0., 10.);
 
  // dE vs momentum and dE vs PDG species
  fh_mvd_dEkeV_vs_p = MakeQaObject<TH2D>("h_mvd_dEkeV_vs_p", "MVD dE vs momentum; p [GeV/c]; dE [keV]", 200, 0., 10.,
                                         200, 0., 500.);  // tune dE range
 
  fh_mvd_dEkeV_vs_pdg = MakeQaObject<TH2D>("h_mvd_dEkeV_vs_pdg", "MVD dE vs particle species;particle;dE [keV]",
                                           nSpeciesBins, 0.5, nSpeciesBins + 0.5, 200, 0., 500.);
 
  labelSpecies(fh_mvd_dEkeV_vs_pdg);
 
 
  fh_nMvdPointsEvt = MakeQaObject<TH1D>("h_nMvdPointsEvt", "MVD points per event;N points;Events", 5000, 0, 5000);
  fh_nDeltaEvt =
    MakeQaObject<TH1D>("h_nDeltaEvt", "Secondary e^{#pm} with MVD hits per event;N tracks;Events", 1000, 0, 1000);
 
  const int nbin = std::max(1, fKeVMax / fKeVBin);
  fh_sec_count_bySpecies =
    MakeQaObject<TH1D>("h_sec_count_bySpecies", "Secondary tracks with MVD hits by species;Class;Counts", 9, 0.5, 9.5);
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(1, "#gamma");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(2, "e^{#pm}");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(3, "#mu^{#pm}");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(4, "#pi^{#pm}");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(5, "K^{#pm}");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(6, "p/p#bar");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(7, "n/n#bar");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(8, "ions");
  fh_sec_count_bySpecies->GetXaxis()->SetBinLabel(9, "other");
 
  fh_dEkeV_vs_pidClass = MakeQaObject<TH2D>(
    "h_dEkeV_vs_pidClass", "#DeltaE per MVD point (keV) vs PID class (non-primaries with hits);PID class;#DeltaE [keV]",
    9, 0.5, 9.5, nbin, 0., fKeVMax);
  for (int b = 1; b <= 9; ++b) {
    const char* lbl = nullptr;
    switch (b) {
      case 1: lbl = "#gamma"; break;
      case 2: lbl = "e^{#pm}"; break;
      case 3: lbl = "#mu^{#pm}"; break;
      case 4: lbl = "#pi^{#pm}"; break;
      case 5: lbl = "K^{#pm}"; break;
      case 6: lbl = "p/p#bar"; break;
      case 7: lbl = "n/n#bar"; break;
      case 8: lbl = "ions"; break;
      default: lbl = "other"; break;
    }
    fh_dEkeV_vs_pidClass->GetXaxis()->SetBinLabel(b, lbl);
  }
 
  fhZdistMvd = MakeQaObject<TH1D>("hZdistMvd", "MVD hit z, z[cm]", 200, -100, 100);
 
  // dE (keV) vs momentum (GeV/c) at the point (non-primaries)
  fh_dEkeV_vs_p =
    MakeQaObject<TH2D>("h_dEkeV_vs_p", "#DeltaE (keV) vs p at point (non-primaries);p [GeV/c];#DeltaE [keV]", 200, 0.,
                       10.,  // tune to beam
                       nbin, 0., fKeVMax);
}
// --------------------------------------------------------------------------------------------------------------------- //
void CbmMvdMcQaTask::ExecQa()
{
  if (!IsMCUsed()) return;
  int nMcEvents = fMcEventList->GetNofEvents();
  for (int iEv = 0; iEv < nMcEvents; iEv++) {
    const auto fileId  = fMcEventList->GetFileIdByIndex(iEv);
    const auto eventId = fMcEventList->GetEventIdByIndex(iEv);
    ProcessMcTracks(fileId, eventId);
    ProcessMvdPoints(fileId, eventId);
  }
}
 
void CbmMvdMcQaTask::ProcessMcTracks(const int& fileId, const int& eventId)
{
  const int nMcTracks = fMcTracks->Size(fileId, eventId);
  fh_nTracksAll->Fill(nMcTracks);  
 
  fPdgBuf.assign(nMcTracks, 0);
  fOriginBuf.assign(nMcTracks, Origin::Other);
 
  for (int iTrack = 0; iTrack < nMcTracks; iTrack++) {
    const auto* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, iTrack));
    if (!mcTrack) continue;
    // const auto* motherTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, mcTrack->GetMotherId()));
    // if (!motherTrack) continue;
    const int pdg    = mcTrack->GetPdgCode();
    const Origin org = ClassifyTrack_(fileId, eventId, iTrack);
    const double p   = mcTrack->GetP();
    const double x0  = mcTrack->GetStartX();
    const double y0  = mcTrack->GetStartY();
    const double z0  = mcTrack->GetStartZ();
 
    fPdgBuf[iTrack]    = pdg;
    fOriginBuf[iTrack] = org;
    // map PDG to species bin for labelled histograms
    const int bin = MapPdgToSpeciesBin_(pdg);
    fh_mc_pdg_all->Fill(bin);
    fh_mc_startZ->Fill(z0);
    fh_mc_startX_vs_Z->Fill(x0, z0);
    fh_mc_startY_vs_Z->Fill(y0, z0);
    fh_mc_startY_vs_X->Fill(x0, y0);
 
    if (org == Origin::Primary) {
      ++fAllNPrim;
      ++fAllCountsPrim[pdg];
 
      fh_mc_pdg_prim->Fill(bin);
      fh_mc_p_prim->Fill(p);
    }
    else if (org == Origin::Secondary || org == Origin::DeltaSec) {
      ++fAllNSec;
      ++fAllCountsSec[pdg];
 
      fh_mc_pdg_sec->Fill(bin);
      fh_mc_p_sec->Fill(p);
    }
    else {
      // Origin::Other: ignore for now
    }
  }
}
// --------------------------------------------------------------------------------------------------------------------- //
void CbmMvdMcQaTask::ProcessMvdPoints(const int& fileId, const int& eventId)
{
  std::vector<int> nPointsStation(fNStations, 0);
  int nDeltaThisEvt = 0;
  TVector3 vIn;   
  TVector3 vOut;  
 
  const int nMvdPoints = fMvdPoints->Size(fileId, eventId);
  const int nMcTracks  = fMcTracks->Size(fileId, eventId);
 
  fh_nMvdPointsEvt->Fill(nMvdPoints);  
  for (int iPoint = 0; iPoint < nMvdPoints; iPoint++) {
    auto* mvdPoint = static_cast<CbmMvdPoint*>(fMvdPoints->Get(fileId, eventId, iPoint));
    if (!mvdPoint) continue;
 
    const auto* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, mvdPoint->GetTrackID()));
    if (!mcTrack) continue;
    // const auto* motherTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, mcTrack->GetMotherId()));
    // if (!motherTrack) continue;
    const int trackId = mvdPoint->GetTrackID();
 
    if (trackId < 0 || trackId >= nMcTracks) {
      LOG(fatal) << "Mvd point " << iPoint << ": trackId " << trackId << " doesn't belong to [0," << nMcTracks - 1
                 << "]";
      break;
    }
 
    const int pdg        = fPdgBuf[trackId];
    const auto origin    = fOriginBuf[trackId];
    const int speciesBin = MapPdgToSpeciesBin_(pdg);
    Int_t motherId       = mcTrack->GetMotherId();
    Int_t pdgCode        = mcTrack->GetPdgCode();
 
    mvdPoint->Position(vIn);
    mvdPoint->PositionOut(vOut);
    Double_t xMC = 0.5 * (vIn.X() + vOut.X());
    Double_t yMC = 0.5 * (vIn.Y() + vOut.Y());
    Double_t zMC = 0.5 * (vIn.Z() + vOut.Z());
 
    fhZdistMvd->Fill(zMC);
    const int stationId = GetStationIdFromZ_(zMC);
    if (stationId < 0 || stationId >= fNStations) continue;
 
    const int cls = PidClass_(pdgCode);
 
    double p            = 0.;
    p                   = mcTrack->GetP();
    const double dE_keV = PointELossGeV_(mvdPoint) * 1e6;
    fh_mvd_pdg_all->Fill(speciesBin);
 
    nPointsStation[stationId]++;
 
    if (origin == Origin::Primary) {
      fh_mvd_xy_prim_station[stationId]->Fill(xMC, yMC);
      fh_mvd_pdg_prim->Fill(speciesBin);
      fh_mvd_p_prim->Fill(p);
    }
    else {
      fh_mvd_xy_sec_station[stationId]->Fill(xMC, yMC);
      fh_mvd_pdg_sec->Fill(speciesBin);
      fh_mvd_p_sec->Fill(p);
    }
    fh_mvd_dEkeV_vs_p->Fill(p, dE_keV);
    fh_mvd_dEkeV_vs_pdg->Fill(speciesBin, dE_keV);
 
    if (motherId < 0) {
      GetOrBookElossHist(pdgCode, Origin::Primary)->Fill(dE_keV);
    }
    else {
      GetOrBookElossHist(pdgCode, Origin::Secondary)->Fill(dE_keV);
      fh_sec_count_bySpecies->Fill(cls);
    }
    fh_dEkeV_vs_pidClass->Fill(cls, dE_keV);
 
    if (fh_dEkeV_vs_p) {
      const double px  = mvdPoint->GetPx();
      const double py  = mvdPoint->GetPy();
      const double pz  = mvdPoint->GetPz();
      const double mom = std::sqrt(px * px + py * py + pz * pz);
      fh_dEkeV_vs_p->Fill(mom, dE_keV);
    }
 
    if (std::abs(pdgCode) == 11) ++nDeltaThisEvt;
  }
  for (int iSt = 0; iSt < fNStations; ++iSt) {
    fh_mvd_nPointsEvt_station[iSt]->Fill(nPointsStation[iSt]);
  }
 
  fh_nDeltaEvt->Fill(nDeltaThisEvt);
}
int CbmMvdMcQaTask::GetStationIdFromZ_(double z) const
{
  if (z > fZMinStation0 && z < fZMaxStation0) return 0;
  if (z > fZMinStation1 && z < fZMaxStation1) return 1;
  if (z > fZMinStation2 && z < fZMaxStation2) return 2;
  if (z > fZMinStation3 && z < fZMaxStation3) return 3;
  return -1;
}
 
// ----- helpers: PDG naming ------------------------------------------------- //
bool CbmMvdMcQaTask::DecodeIonPDG_(int pdg, int& Z, int& A)
{
  // PDG for ions: 10LZZZAAAI  -> commonly used as 1000000000 + Z*10000 + A*10 + I
  // where L is strangeness, Z is charge,
  // A is number of nucleons, and I is isomer level.
 
  int apdg = std::abs(pdg);
  if (apdg < 1000000000) return false;
  int code = apdg - 1000000000;
  Z        = (code / 10000) % 1000;
  A        = (code / 10) % 1000;
  return (Z > 0 && A >= Z);
}
 
// --------------------------------------------------------------------------------------------------------------------- //
void CbmMvdMcQaTask::Check()
{
  if (!IsMCUsed()) {
    return;  
  }
}
 
// --------------------------------------------------------------------------------------------------------------------- //
void CbmMvdMcQaTask::CreateSummary()
{
  LOG(info) << "=== MC primary species counts ===";
  for (auto& kv : fAllCountsPrim)
    LOG(info) << "PDG " << kv.first << " : " << kv.second;
 
  LOG(info) << "=== MC secondary species counts ===";
  for (auto& kv : fAllCountsSec)
    LOG(info) << "PDG " << kv.first << " : " << kv.second;
}
// --------------------------------------------------------------------------------------------------------------------- //
TH1F* CbmMvdMcQaTask::GetOrBookElossHist(int pdg, Origin origin)
{
  int oi = static_cast<int>(origin);
  if (oi < 0 || oi >= kNumOrigins) oi = static_cast<int>(Origin::Other);
 
  auto& map  = fHistsElossKeVByPDG[oi];
  auto* list = fLists[oi];
 
  if (auto it = map.find(pdg); it != map.end()) return it->second;
 
  const int nbins   = 400;
  const double xmin = 0.0;
  const double xmax = 4000.0;  // keV
 
  std::string name  = Form("hEloss_keV_%s_pdg_%s%d", kOriginName[oi], pdg < 0 ? "m" : "", std::abs(pdg));
  std::string title = Form("dE (keV) [%s], PDG %d;dE [keV];counts", kOriginName[oi], pdg);
 
  TH1F* h = new TH1F(name.c_str(), title.c_str(), nbins, xmin, xmax);
  h->SetDirectory(nullptr);
 
  map[pdg] = h;
  list->Add(h);
  return h;
}
// ----- per-event classification helpers ----------------------------------- //
int CbmMvdMcQaTask::PidClass_(int pdg) const
{
  if (pdg == 22) return 1;
  const int a = std::abs(pdg);
  if (a == 11) return 2;
  if (a == 13) return 3;
  if (a == 211) return 4;
  if (a == 321) return 5;
  if (a == 2212) return 6;
  if (a == 2112) return 7;
  if (a > 1000000000) return 8;  
  return 9;
}
// --------------------------------------------------------------------------------------------------------------------- //
int CbmMvdMcQaTask::MapPdgToSpeciesBin_(int pdg) const
{
  // bin layout:
  //  1: Gamma
  //  2: e-
  //  3: e+
  //  4: mu-
  //  5: mu+
  //  6: pi+
  //  7: pi-
  //  8: pi0
  //  9: K+
  // 10: K-
  // 11: K0_S
  // 12: K0_L
  // 13: p
  // 14: pbar
  // 15: n
  // 16: nbar
  // 17: ions (PDG > 1e9)
  // 18: other
 
  if (pdg == 22) return 1;
 
  if (pdg == 11) return 2;
  if (pdg == -11) return 3;
 
  if (pdg == 13) return 4;
  if (pdg == -13) return 5;
 
  if (pdg == 211) return 6;
  if (pdg == -211) return 7;
 
  if (pdg == 111) return 8;
 
  if (pdg == 321) return 9;
  if (pdg == -321) return 10;
 
  if (pdg == 310) return 11;
  if (pdg == 130) return 12;
 
  if (pdg == 2212) return 13;
  if (pdg == -2212) return 14;
 
  if (pdg == 2112) return 15;
  if (pdg == -2112) return 16;
 
  if (std::abs(pdg) > 1000000000) return 17;  // ions
 
  return 18;  // everything else
}
// --------------------------------------------------------------------------------------------------------------------- //
CbmMvdMcQaTask::Origin CbmMvdMcQaTask::ClassifyTrack_(const int& fileId, const int& eventId, Int_t trackindex) const
{
  auto* trk = static_cast<const CbmMCTrack*>(fMcTracks->Get(fileId, eventId, trackindex));
  if (!trk) {
    LOG(info) << GetName() << ": MCTrack not found for file " << fileId << " event " << eventId << " track "
              << trackindex;
    return Origin::Other;
  }
  const int pdg      = std::abs(trk->GetPdgCode());
  const int motherId = trk->GetMotherId();
 
  if (motherId < 0) {
    return Origin::Primary;
  }
 
  if (pdg == 11) {
    return Origin::DeltaSec;
  }
 
  return Origin::Secondary;
}
// --------------------------------------------------------------------------------------------------------------------- //
double CbmMvdMcQaTask::PointELossGeV_(const void* vp) const
{
  const auto* p = static_cast<const FairMCPoint*>(vp);  
  return p ? p->GetEnergyLoss() : 0.0;                  
}