KfpV0FinderCutter.cxx

Go to the documentation of this file.

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

 
#include "algo/kfp/KfpV0FinderCutter.h"
 
#include "KFParticleFinder.h"
 
#include <iomanip>
#include <set>
#include <sstream>
 
using cbm::algo::kfp::ETrackPid;
using cbm::algo::kfp::EV0Type;
using cbm::algo::kfp::V0FinderCutter;
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
V0FinderCutter<V0Type>::V0FinderCutter(const std::string& filename)
  : fCuts(cbm::algo::yaml::ReadFromFile<CutValues>(filename))
{
  fCuts.Init();
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
void V0FinderCutter<V0Type>::ApplyCutsToKFParticleFinder(KFParticleFinder* pFinder) const
{
  const auto& cuts = fCuts.kfpFinderCuts;
  pFinder->SetChiPrimaryCut2D(cuts.maxChi2NdfPrim);
  pFinder->SetLdLCut2D(cuts.minDecayLDL);
  pFinder->SetLCut(cuts.minDecayLength);
  pFinder->SetChi2Cut2D(cuts.maxChi2NdfGeo);
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
void V0FinderCutter<V0Type>::ReadCutsFromYaml(const std::string& filename)
{
  fCuts = cbm::algo::yaml::ReadFromFile<CutValues>(filename);
  fCuts.Init();
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
bool V0FinderCutter<V0Type>::SelectParticle(const ParticleProperty& property) const
{
  constexpr auto SelectDaughter = [](const DaughterProperty& pr, const DaughterCuts& cut) constexpr->bool
  {
    if (pr.beta < cut.beta[0] || pr.beta > cut.beta[1]) return false;
    if (pr.dcaRelToOrigin < cut.dcaRelToOrigin[0] || pr.dcaRelToOrigin > cut.dcaRelToOrigin[1]) return false;
    return true;
  };
 
  const auto& cuts = fCuts.particleCuts;
  if (!SelectDaughter(property.daughterA, cuts.daughterA)) return false;
  if (!SelectDaughter(property.daughterB, cuts.daughterB)) return false;
  if (property.dcaRelToPrimVertex > cuts.maxDcaRelToPrimVertex) return false;
  if (property.dcaBetweenDaughters > cuts.maxDcaBetweenDaughters) return false;
  if (property.decayLength < cuts.decayLength[0] || property.decayLength > cuts.decayLength[1]) return false;
  if (property.daughterVtxZ < cuts.daughterVtxZ[0] || property.daughterVtxZ > cuts.daughterVtxZ[1]) return false;
  if (property.openingAngle < cuts.minOpeningAngle) return false;
  if (property.chi2NdfTopo < cuts.chi2NdfTopo[0] || property.chi2NdfTopo > cuts.chi2NdfTopo[1]) return false;
  if (property.chi2NdfGeo < cuts.chi2NdfGeo[0] || property.chi2NdfGeo > cuts.chi2NdfGeo[1]) return false;
 
  return true;
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
bool V0FinderCutter<V0Type>::SelectTrack(const TrackProperty& property) const
{
  // Check, if track has a defined PID
  if (property.pid == ETrackPid::Undef) return false;
 
  // Check, if DCA was defined (NOTE: equivalent to the previous cut, but only with TopoPid)
  // (a) undefined => less then two STS hits
  if (std::isnan(property.dcaRelToOrigin)) return false;
 
  // Beta physical cuts
  // (a) undefined   => no TOF hits
  // (b) < 0 or > 1  => miscalibration of TOF/BMON
  if (std::isnan(property.beta) || property.beta <= 0 || property.beta >= 1) return false;
 
  // Other cuts (depending on PID)
  const auto& cuts = fCuts.trackCuts;
 
  // Track selection depending on the PDG
  // (a) track beta
  auto Check = [&property, &cuts](ETrackPid trackPid) constexpr->bool
  {
    const auto& pidCuts = cuts.pidCuts[static_cast<size_t>(trackPid)];
    if (property.beta < pidCuts.beta[0] || property.beta > pidCuts.beta[1]) return false;
    return true;
  };
 
  if (!Check(property.pid)) return false;
 
  return true;
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
void V0FinderCutter<V0Type>::SetCuts(const CutValues& cuts)
{
  // Apply validity check here (maybe)
  fCuts = cuts;
  fCuts.Init();
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
std::string V0FinderCutter<V0Type>::ToString() const
{
  using std::setw;
  std::stringstream msg;
 
  std::string sV0{};
  std::string sDaughterA{};
  std::string sDaughterB{};
  switch (V0Type) {
    case EV0Type::Lambda:
      sV0        = "Lambda";
      sDaughterA = "pi-";
      sDaughterB = "p";
      break;
  }
 
  auto PrintTrackPidCuts = [&](int iPid) {
    const auto& cuts = fCuts.trackCuts.pidCuts[iPid];
    msg << "\n\t" << cbm::algo::kfp::ToString(static_cast<ETrackPid>(iPid)) << ':';
    msg << "\n\t\tbeta: " << setw(12) << cuts.beta[0] << " -- " << setw(12) << cuts.beta[1];
  };
 
  msg << "======= V0-finder selector for " << sV0 << " -> " << sDaughterA << " + " << sDaughterB;
 
  msg << "\nTRACK PRE-SELECTION CUTS:";
  for (int iPid = 0; iPid < static_cast<int>(ETrackPid::Undef) + 1; ++iPid) {
    PrintTrackPidCuts(iPid);
  }
 
  const auto& kfpFinderCuts = fCuts.kfpFinderCuts;
  msg << "\nKFP-FINDER INTRINSIC CUTS:";
  msg << "\n\tmin. decay length [cm]: " << setw(12) << kfpFinderCuts.minDecayLength;
  msg << "\n\tmin. l / dl:            " << setw(12) << kfpFinderCuts.minDecayLDL;
  msg << "\n\tmax chi2/NDF primary:   " << setw(12) << kfpFinderCuts.maxChi2NdfPrim;
  msg << "\n\tmax chi2/NDF geo:       " << setw(12) << kfpFinderCuts.maxChi2NdfGeo;
 
  auto PrintDaughterCuts = [&](const DaughterCuts& cuts) {
    msg << "\n\t\tDCA rel. to origin: [cm]" << setw(12) << cuts.dcaRelToOrigin[0] << " -- " << setw(12)
        << cuts.dcaRelToOrigin[1];
    msg << "\n\t\tbeta:                   " << setw(12) << cuts.beta[0] << " -- " << setw(12) << cuts.beta[1];
  };
 
  const auto& prtCuts = fCuts.particleCuts;
  msg << "\nV0 SELECTION CUTS:";
  msg << "\n\t" << sDaughterA << ':';
  PrintDaughterCuts(prtCuts.daughterA);
  msg << "\n\t" << sDaughterB << ':';
  PrintDaughterCuts(prtCuts.daughterB);
  msg << "\n\tmax DCA rel. to primary vertex [cm]: " << prtCuts.maxDcaRelToPrimVertex;
  msg << "\n\tmax DCA between daughters [cm]:      " << prtCuts.maxDcaBetweenDaughters;
  msg << "\n\tdecay length [cm]:                   " << prtCuts.decayLength[0] << " -- " << prtCuts.decayLength[1];
  msg << "\n\tdaughter z-vertex [cm]:              " << prtCuts.daughterVtxZ[0] << " -- " << prtCuts.daughterVtxZ[1];
  msg << "\n\tmin opening angle [radians]:         " << prtCuts.minOpeningAngle;
  msg << "\n\tchi2/NDF topology:                   " << prtCuts.chi2NdfTopo[0] << " -- " << prtCuts.chi2NdfTopo[1];
  msg << "\n\tchi2/NDF geo:                        " << prtCuts.chi2NdfGeo[0] << " -- " << prtCuts.chi2NdfGeo[1];
 
  return msg.str();
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<EV0Type V0Type>
void cbm::algo::kfp::V0FinderCutter<V0Type>::CutValues::Init()
{
  auto& trackPidCuts = trackCuts.pidCuts;
  std::set<ETrackPid> knownPids;
  for (const auto& entry : trackPidCuts) {
    knownPids.insert(entry.pid);
  }
 
  if (knownPids.size() != static_cast<size_t>(ETrackPid::Undef) + 1) {
    std::stringstream msg;
    msg << "V0FinderCutter::CutValues::Init(): the node \"trackCuts/pidCuts\" must contain exactly the same amount ";
    msg
      << "of elements, as the ETrackPid enumerator. Each element must have a unique pid value. The pid values, which ";
    msg << "were read from the parameter file, are: ";
    for (const auto& p : knownPids) {
      msg << cbm::algo::kfp::ToString(p) << ' ';
    }
    throw std::runtime_error(msg.str());
  }
 
  std::sort(trackPidCuts.begin(), trackPidCuts.end(), [](const auto& lhs, const auto& rhs) {
    return static_cast<size_t>(lhs.pid) < static_cast<size_t>(rhs.pid);
  });
}
 
template class cbm::algo::kfp::V0FinderCutter<EV0Type::Lambda>;