From 86ff13083644294f29aa10c7a0e2604ed4d4d1b8 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Tue, 21 Jan 2025 10:26:31 +0100
Subject: [PATCH] Introduction of mCBM V0 finder task
---
algo/ca/core/tracking/CaTrackFinderWindow.cxx | 4 +-
core/qa/CbmQaIO.cxx | 2 +-
macro/beamtime/mcbm2022/mcbm_event_reco_L1.C | 41 +-
macro/beamtime/mcbm2024/reco_mcbm.sh | 9 +-
macro/mcbm/CMakeLists.txt | 2 +-
macro/mcbm/mcbm_hadron_kfp_ana.C | 281 +++++++
macro/run/run_qa.C | 7 +-
reco/KF/CbmKFParticleFinder.cxx | 2 +-
reco/KF/CbmKFParticleFinder.h | 3 +-
reco/KF/CbmKFV0FinderQa.cxx | 61 ++
reco/KF/CbmKFV0FinderQa.h | 66 ++
reco/KF/CbmKFV0FinderTask.cxx | 774 ++++++++++++++++++
reco/KF/CbmKFV0FinderTask.h | 461 +++++++++++
reco/KF/KFParticleInterface.cmake | 2 +
reco/KF/KFParticleInterfaceLinkDef.h | 2 +
reco/L1/CbmCaTimeSliceReader.cxx | 1 +
.../digis/CbmAlgoBuildRawEvents.cxx | 3 +-
reco/tasks/CbmTaskInspectDigiTimeslice.cxx | 2 +
18 files changed, 1697 insertions(+), 26 deletions(-)
create mode 100644 macro/mcbm/mcbm_hadron_kfp_ana.C
create mode 100644 reco/KF/CbmKFV0FinderQa.cxx
create mode 100644 reco/KF/CbmKFV0FinderQa.h
create mode 100644 reco/KF/CbmKFV0FinderTask.cxx
create mode 100644 reco/KF/CbmKFV0FinderTask.h
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
index 1dd20f4040..cce227803a 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
@@ -498,7 +498,9 @@ namespace cbm::algo::ca
if (wData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
}
}
- fvTrackCandidates.push_back(best_tr);
+ // NOTE: Temporarely disable the warnings, one has to provide more specific coefficient in the capacity
+ // reservation
+ fvTrackCandidates.push_back_no_warning(best_tr);
ca::Branch& tr = fvTrackCandidates.back();
tr.SetStation(istaF);
tr.SetId(fvTrackCandidates.size() - 1);
diff --git a/core/qa/CbmQaIO.cxx b/core/qa/CbmQaIO.cxx
index 426e3e57f3..926cfb4bfa 100644
--- a/core/qa/CbmQaIO.cxx
+++ b/core/qa/CbmQaIO.cxx
@@ -36,7 +36,7 @@ CbmQaIO::~CbmQaIO() {}
void CbmQaIO::SetTH1Properties(TH1* pHist) const
{
// Set default histo properties
- pHist->GetYaxis()->SetLabelOffset(0.95);
+ //pHist->GetYaxis()->SetLabelOffset(0.95);
TPaveStats* stats = cbm::qa::util::GetHistStats(pHist);
assert(stats);
diff --git a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
index b9bd63107f..35ba48cbc8 100644
--- a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
+++ b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
@@ -36,13 +36,11 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
Bool_t bFSD = kFALSE,
TString sInpFile = ""
)
+/* clang-format on */
{
- /// FIXME: Re-enable clang formatting after parameters initial values setting
- /* clang-format on */
-
// --- Logger settings ----------------------------------------------------
TString logLevel = "INFO"; //"INFO";
- TString logVerbosity = "VERYLOW"; //"VERYLOW";
+ TString logVerbosity = "LOW"; //"VERYLOW";
// ------------------------------------------------------------------------
// ----- Environment --------------------------------------------------
@@ -421,17 +419,28 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
Int_t stsAddress11 = CbmStsAddress::GetAddress(1, 2, 1, 2, 0, 0); // U1 L2 M2 6 cm
- std::cout << "STS address01 " << std::dec << stsAddress01 << " " << std::hex << stsAddress01 << std::endl;
- std::cout << "STS address02 " << std::dec << stsAddress02 << " " << std::hex << stsAddress02 << std::endl;
- std::cout << "STS address03 " << std::dec << stsAddress03 << " " << std::hex << stsAddress03 << std::endl;
- std::cout << "STS address04 " << std::dec << stsAddress04 << " " << std::hex << stsAddress04 << std::endl;
- std::cout << "STS address05 " << std::dec << stsAddress05 << " " << std::hex << stsAddress05 << std::endl;
- std::cout << "STS address06 " << std::dec << stsAddress06 << " " << std::hex << stsAddress06 << std::endl;
- std::cout << "STS address07 " << std::dec << stsAddress07 << " " << std::hex << stsAddress07 << std::endl;
- std::cout << "STS address08 " << std::dec << stsAddress08 << " " << std::hex << stsAddress08 << std::endl;
- std::cout << "STS address09 " << std::dec << stsAddress09 << " " << std::hex << stsAddress09 << std::endl;
- std::cout << "STS address10 " << std::dec << stsAddress10 << " " << std::hex << stsAddress10 << std::endl;
- std::cout << "STS address11 " << std::dec << stsAddress11 << " " << std::hex << stsAddress11 << std::endl;
+ std::cout << "STS address01 " << std::dec << stsAddress01 << " " << std::hex << stsAddress01 << " "
+ << CbmStsAddress::ToString(stsAddress01) << std::endl;
+ std::cout << "STS address02 " << std::dec << stsAddress02 << " " << std::hex << stsAddress02 << " "
+ << CbmStsAddress::ToString(stsAddress02) << std::endl;
+ std::cout << "STS address03 " << std::dec << stsAddress03 << " " << std::hex << stsAddress03 << " "
+ << CbmStsAddress::ToString(stsAddress03) << std::endl;
+ std::cout << "STS address04 " << std::dec << stsAddress04 << " " << std::hex << stsAddress04 << " "
+ << CbmStsAddress::ToString(stsAddress04) << std::endl;
+ std::cout << "STS address05 " << std::dec << stsAddress05 << " " << std::hex << stsAddress05 << " "
+ << CbmStsAddress::ToString(stsAddress05) << std::endl;
+ std::cout << "STS address06 " << std::dec << stsAddress06 << " " << std::hex << stsAddress06 << " "
+ << CbmStsAddress::ToString(stsAddress06) << std::endl;
+ std::cout << "STS address07 " << std::dec << stsAddress07 << " " << std::hex << stsAddress07 << " "
+ << CbmStsAddress::ToString(stsAddress07) << std::endl;
+ std::cout << "STS address08 " << std::dec << stsAddress08 << " " << std::hex << stsAddress08 << " "
+ << CbmStsAddress::ToString(stsAddress08) << std::endl;
+ std::cout << "STS address09 " << std::dec << stsAddress09 << " " << std::hex << stsAddress09 << " "
+ << CbmStsAddress::ToString(stsAddress09) << std::endl;
+ std::cout << "STS address10 " << std::dec << stsAddress10 << " " << std::hex << stsAddress10 << " "
+ << CbmStsAddress::ToString(stsAddress10) << std::endl;
+ std::cout << "STS address11 " << std::dec << stsAddress11 << " " << std::hex << stsAddress11 << " "
+ << CbmStsAddress::ToString(stsAddress11) << std::endl;
// --- Now we can define the sensor parameter set and tell recoSts to use it
auto sensorParSet = new CbmStsParSetSensor("CbmStsParSetSensor", "STS sensor parameters"
@@ -712,7 +721,7 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
// ----- Function needed for CTest runtime dependency -----------------
- // RemoveGeoManager();
+ RemoveGeoManager();
// ------------------------------------------------------------------------
/// --- Screen output for automatic tests
diff --git a/macro/beamtime/mcbm2024/reco_mcbm.sh b/macro/beamtime/mcbm2024/reco_mcbm.sh
index 8f09705466..51c9e57176 100755
--- a/macro/beamtime/mcbm2024/reco_mcbm.sh
+++ b/macro/beamtime/mcbm2024/reco_mcbm.sh
@@ -36,6 +36,7 @@
# Auxiliary options:
# -n [--nts] <N_TS> Number of timeslices to procede
# --setup, --unpack, --reco, --qa, --qa-module
+# --data-dir Path to output directory
# --param-online <path> Path to the online parameters
# NOTE: should be ${VMCWORKDIR}/parameters/online, but local user may have copy for testing
#
@@ -143,6 +144,9 @@ while [[ $# -gt 0 ]]; do
-n | --nts )
N_TS=${2}
;;
+ --data-dir )
+ DATA_TOP_DIR=${2}
+ ;;
esac
shift
done
@@ -272,7 +276,10 @@ if [[ ${DO_RECO} -eq 1 ]]; then
DATA_PREF="${DATA_TOP_DIR}"
PARS="${RUN},${N_TS},\"${DATA_PREF}\",\"${DATA_PREF}\",${UNP_FILE_ID},${RECO_MVD},${RECO_STS},${RECO_TRD}"
PARS="${PARS},${RECO_TRD2d},${RECO_RICH},${RECO_MUCH},${RECO_TOF},${RECO_TOFtr},${RECO_PSD},${RECO_ALI},${RECO_EvB}"
- PARS="${PARS},${RECO_CA},${RECO_QA},${RECO_FSD},\"${INP_FILE}\",${RECO_PV}"
+ PARS="${PARS},${RECO_CA},${RECO_QA},${RECO_FSD},\"${INP_FILE}\""
+ if [[ ${MACRO_RECO} == "*.mcbm2024*" ]]; then
+ PARS="${PARS},${RECO_PV}"
+ fi
root -b -l -q ${MACRO_RECO}"(${PARS})" &> ${RECO_LOG}
cat ${RECO_LOG}
diff --git a/macro/mcbm/CMakeLists.txt b/macro/mcbm/CMakeLists.txt
index 8056b9eeed..7313634988 100644
--- a/macro/mcbm/CMakeLists.txt
+++ b/macro/mcbm/CMakeLists.txt
@@ -138,7 +138,7 @@ ForEach(setup IN LISTS cbm_setup)
EndForEach(setup IN LISTS cbm_setup)
Install(FILES .rootrc mcbm_transport.C mcbm_digi.C mcbm_reco_event.C mcbm_match_check.C mcbm_hadron_analysis.C mcbm_qa.C
- mcbm_transport_boxgen.C
+ mcbm_transport_boxgen.C mcbm_hadron_kfp_ana.C
DESTINATION share/cbmroot/macro/mcbm
)
Install(DIRECTORY modules DESTINATION share/cbmroot/macro/mcbm)
diff --git a/macro/mcbm/mcbm_hadron_kfp_ana.C b/macro/mcbm/mcbm_hadron_kfp_ana.C
new file mode 100644
index 0000000000..b48dd99436
--- /dev/null
+++ b/macro/mcbm/mcbm_hadron_kfp_ana.C
@@ -0,0 +1,281 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file mcbm_hadron_kfp_ana.C
+/// \brief Macro for hadron analysis with the KfParticleFinder in mCBM
+/// \since 08.01.2025
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+
+/// \brief Main function of the macro
+/// \param uRunId Run identifier
+/// \param firstTimeslice First timeslice
+/// \param nTimeslices Number of timeslices to proceed
+/// \param sRecoFile Input with reconstructed data tree
+/// \param sCollTraFile Input with transported collision MC data tree [optional]
+/// \param sSignTraFile Input with transported signal MC data tree [optional]
+/// \param sBeamTraFile Input with transported beam MC data tree [optional]
+/// \param sSinkFile Output file
+/// \param sGeoFile Geometry file
+/// \param sParInFile Input parameters file
+/// \param bUseAlignment Flag: use alignment
+/// \param bMixEvent Flag: perform mixed-event analysis
+/// \param bUseMc Flag: use MC-information (if available)
+///
+/// \note File extention conventions used:
+/// 1. Input:
+/// - geometry: <geoSetupTag>.geo.root
+/// - reco: <dataLabel>.rec.root
+/// - transport:
+/// - collisions: <genLabel>.tra.root
+/// - signal: <genLabel>.tra.root
+/// - beam: <genLabel>.tra.root
+/// 2. Output:
+/// - sink: <dataLabel>.kfp.ana.root
+/// - monitor: <dataLabel>.kfp.mon.root
+/// - parameters: <dataLabel>.kfp.par.root
+///
+/// <geoSetupTag> -- geometry setup tag provided by cbm::mcbm::GetSetupFromRunId(uRunId)
+/// <dataLabel> -- label of data files reflecting the name prefix of the original TSA file
+/// E.g.: <dataLabel> = 2391_node8_0_0055
+/// <genLabel> -- label of the simulated files reflecting the name prefix of the original generated file:
+/// E.g.: <genLabel> = urqmd.auau.1.24gev.mbias.00109
+///
+/// TODO: Probably, use a combined label for output: <dataLabel>.<genLabel>.kfp.ana.root
+///
+/* clang-format off */
+Bool_t mcbm_hadron_kfp_ana(UInt_t uRunId = 2391,
+ Int_t firstTimeslice = 0,
+ Int_t nTimeslices = 10,
+ TString sRecoFile = "./data/2391_first20Ts.rec.root",
+ TString sCollTraFile = "",
+ TString sSignTraFile = "",
+ TString sBeamTraFile = "",
+ TString sSinkFile = "./data/2391_first20Ts.kfp.ana.root",
+ TString sGeoFile = "./data/mcbm_beam_2022_05_23_nickel.geo.root",
+ TString sParInFile = "./data/2391_first20Ts.par.root",
+ Bool_t bUseAlignment = true,
+ Bool_t bMixEvent = false,
+ Bool_t bUseMc = false
+)
+/* clang-format on */
+{
+ // -- Logger settings ------------------------------------------------------------------------------------------------
+ const TString logLevel{"INFO"}; // Logger level
+ const TString logVerbosity{"LOW"}; // Logger verbosity
+ const bool logColored{true}; // Colored logger
+ const TString myName{"mcbm_hadron_kfp_ana"}; // Macro name
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Environment settings -------------------------------------------------------------------------------------------
+ TString sOutDir{};
+ if (auto lastSlashPos = sSinkFile.Last('/'); lastSlashPos != -1) {
+ sOutDir = sSinkFile(0, lastSlashPos);
+ }
+ else {
+ sOutDir = "./";
+ }
+ const TString srcDir{gSystem->Getenv("VMCWORKDIR")}; // Top source directory
+ gSystem->Exec("mkdir " + sOutDir); // Output directory
+ //gSystem->Exec("cp $VMCWORKDIR/macro/run/.rootrc ."); // Loading environment
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Auxilary common variables --------------------------------------------------------------------------------------
+ const TString sRunId{TString::Format("%04u", uRunId)}; // String representation of the run id
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Geometry setup tag ---------------------------------------------------------------------------------------------
+ cbm::mcbm::ToForceLibLoad dummy; /// Needed to trigger loading of the library as no fct dict in ROOT6 and CLING
+ TString geoSetupTag{""};
+ try {
+ geoSetupTag = cbm::mcbm::GetSetupFromRunId(uRunId);
+ }
+ catch (const std::invalid_argument& err) {
+ std::cout << "-E- " << myName << ": mapping from runID to geometry setup tag failed: " << err.what() << '\n';
+ return false;
+ }
+ auto* geoSetup = CbmSetup::Instance();
+ geoSetup->LoadSetup(geoSetupTag);
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Start timer ----------------------------------------------------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- FairRun initialisation -----------------------------------------------------------------------------------------
+ FairRunAna* run = new FairRunAna();
+ run->SetGeomFile(sGeoFile);
+ FairFileSource* inputSource = new FairFileSource(sRecoFile);
+ if (bUseMc) {
+ int nTraFiles{0};
+ if (!sCollTraFile.IsNull()) {
+ inputSource->AddFriend(sCollTraFile);
+ nTraFiles++;
+ }
+ if (!sSignTraFile.IsNull()) {
+ inputSource->AddFriend(sSignTraFile);
+ nTraFiles++;
+ }
+ if (!sBeamTraFile.IsNull()) {
+ inputSource->AddFriend(sBeamTraFile);
+ nTraFiles++;
+ }
+ if (nTraFiles == 0) {
+ std::cout << "-E- " << myName << ": MC information was required, but no transport files were provided";
+ return false;
+ }
+ }
+ run->SetSource(inputSource);
+
+ FairRootFileSink* outputSink = new FairRootFileSink(sSinkFile);
+ run->SetSink(outputSink);
+
+ TString sMonitorFile{sSinkFile};
+ sMonitorFile.ReplaceAll(".ana.", ".mon.");
+ FairMonitor::GetMonitor()->EnableMonitor(true, sMonitorFile);
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Logger settings ------------------------------------------------------------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ FairLogger::GetLogger()->SetColoredLog(logColored);
+ //FairLogger::GetLogger()->SetLogScreenLevel("DEBUG");
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Alignment corrections ------------------------------------------------------------------------------------------
+ if (bUseAlignment) {
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TString sAlignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geoSetupTag + ".root";
+ std::cout << "-I- " << myName << ": loading alignment matrices from file " << sAlignmentMatrixFileName << " ... ";
+ TFile* misalignmentMatrixRootfile = TFile::Open(sAlignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile && misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ std::cout << "done\n";
+ }
+ else {
+ std::cout << "failed: running with misaligned geometry\n";
+ }
+ if (matrices) run->AddAlignmentMatrices(*matrices);
+ }
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- KF initialization ----------------------------------------------------------------------------------------------
+ run->AddTask(new CbmTrackingDetectorInterfaceInit());
+ run->AddTask(new CbmKF());
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- V0 Finder ------------------------------------------------------------------------------------------------------
+ auto* pV0 = new cbm::kfp::V0FinderTask(/*verbose = */ 3);
+ // TODO: define all the setters initialization within a config file
+ // The task parameters:
+ pV0->SetMinPionDca(1.5);
+ pV0->SetMinProtonDca(0.5);
+ pV0->SetTzeroOffset(0.12);
+ pV0->SetQpAssignedUncertainty(0.10); // 10% uncertainty for track q/p
+ pV0->SetMixedEventMode(bMixEvent);
+ pV0->SetProcessingMode(cbm::kfp::V0FinderTask::EProcessingMode::TimeBased);
+ pV0->SetPidApproach(cbm::kfp::V0FinderTask::EPidApproach::Topo);
+ pV0->SetPvFindingMode(cbm::kfp::V0FinderTask::EPvUsageMode::Target);
+
+ // Decays to be reconstructed:
+ pV0->AddDecayToReconstructionList(3122); // Lambda
+
+ // KFParticleFinder cuts:
+ pV0->SetChiPrimaryCut2D(0.f);
+ pV0->SetLdLCut2D(1.f);
+ pV0->SetLCut(0.1);
+ pV0->SetChi2Cut2D(50.f);
+
+ // QA output setting
+ {
+ TString sQaFile{sSinkFile};
+ sQaFile.ReplaceAll(".ana.", "_qa.ana.");
+ pV0->SetQaOutputFileName(sQaFile);
+ pV0->SetRunQa(true);
+ }
+ run->AddTask(pV0);
+ // -------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Parameter database ---------------------------------------------------------------------------------------------
+ std::cout << "\n\n";
+ std::cout << "-I- " << myName << ": setting runtime DB\n";
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parInp = new FairParRootFileIo();
+ FairParRootFileIo* parOut = new FairParRootFileIo();
+ parInp->open(sParInFile.Data(), "READ");
+ rtdb->setFirstInput(parInp);
+ TString sParOutFile{sSinkFile};
+ sParOutFile.ReplaceAll(".ana.", ".par.");
+ parOut->open(sParOutFile.Data(), "RECREATE");
+ rtdb->setOutput(parOut);
+ // ------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Run initialisation --------------------------------------------------------------------------------------------
+ std::cout << "\n\n";
+ std::cout << "-I- " << myName << ": initialising run\n";
+ run->Init();
+ rtdb->print();
+ // ------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Run execution -------------------------------------------------------------------------------------------------
+ std::cout << "\n\n";
+ std::cout << "-I- " << myName << ": starting run\n";
+ run->Run(firstTimeslice, nTimeslices);
+ // ------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Finish --------------------------------------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ std::cout << "-I- " << myName << "Macro finished successfully.\n";
+ std::cout << "-I- " << myName << "Output file is " << sSinkFile << '\n';
+ std::cout << "-I- " << myName << "Parameter file is " << sParOutFile << '\n';
+ std::cout << "-I- " << myName << "Real time " << rtime << " s, CPU time " << ctime << " s\n\n";
+ // ------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Resourse monitoring -------------------------------------------------------------------------------------------
+ // Extract the maximal used memory an add is as Dart measurement
+ // This line is filtered by CTest and the value send to CDash
+ FairSystemInfo sysInfo;
+ Float_t maxMemory = sysInfo.GetMaxMemory();
+ std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ std::cout << maxMemory;
+ std::cout << "</DartMeasurement>" << std::endl;
+
+ Float_t cpuUsage = ctime / rtime;
+ std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ std::cout << cpuUsage;
+ std::cout << "</DartMeasurement>" << std::endl;
+ // ------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Memory clean-up -----------------------------------------------------------------------------------------------
+ RemoveGeoManager();
+ // ------------------------------------------------------------------------------------------------------------------
+
+
+ // -- Screen output for automatic tests -----------------------------------------------------------------------------
+ std::cout << " Test passed" << std::endl;
+ std::cout << " All ok " << std::endl;
+ // ------------------------------------------------------------------------------------------------------------------
+
+
+ return true;
+}
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index 16880b48a3..07f9da0f08 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -218,10 +218,11 @@ void run_qa(TString dataTraColl,
// ------------------------------------------------------------------------
// ----- MCDataManager -----------------------------------
+ // TODO: provide to mcbm_qa.C, run_reco_.....C, mcbm_reco_......C
CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0);
- mcManager->AddFile(traFileColl);
- if (!dataTraSign.IsNull()) mcManager->AddFile(traFileSign);
- if (!dataTraBeam.IsNull()) mcManager->AddFile(traFileBeam);
+ mcManager->AddFileToChain(traFileColl, 0);
+ if (!dataTraSign.IsNull()) mcManager->AddFileToChain(traFileSign, 1);
+ if (!dataTraBeam.IsNull()) mcManager->AddFileToChain(traFileBeam, 2);
qaManager->AddTask(mcManager);
// ------------------------------------------------------------------------
diff --git a/reco/KF/CbmKFParticleFinder.cxx b/reco/KF/CbmKFParticleFinder.cxx
index f7da8967f2..7ffc2dd17b 100644
--- a/reco/KF/CbmKFParticleFinder.cxx
+++ b/reco/KF/CbmKFParticleFinder.cxx
@@ -581,7 +581,7 @@ void CbmKFParticleFinder::FillKFPTrackVector(KFPTrackVector* tracks, const vecto
}
}
-double CbmKFParticleFinder::InversedChi2Prob(double p, int ndf) const
+double CbmKFParticleFinder::InversedChi2Prob(double p, int ndf)
{
double epsilon = 1.e-14;
double chi2Left = 0.f;
diff --git a/reco/KF/CbmKFParticleFinder.h b/reco/KF/CbmKFParticleFinder.h
index 9e656426b0..5a38b57db5 100644
--- a/reco/KF/CbmKFParticleFinder.h
+++ b/reco/KF/CbmKFParticleFinder.h
@@ -83,8 +83,9 @@ class CbmKFParticleFinder : public FairTask {
//Add decay to the reconstruction list. By default, if list is empty - all decays are reconstructed.
void AddDecayToReconstructionList(int pdg);
+ static double InversedChi2Prob(double p, int ndf);
+
private:
- double InversedChi2Prob(double p, int ndf) const;
void FillKFPTrackVector(KFPTrackVector* tracks, const std::vector<CbmStsTrack>& vRTracks,
const std::vector<KFFieldVector>& vField, const std::vector<int>& pdg,
const std::vector<int>& trackId, const std::vector<float>& vChiToPrimVtx,
diff --git a/reco/KF/CbmKFV0FinderQa.cxx b/reco/KF/CbmKFV0FinderQa.cxx
new file mode 100644
index 0000000000..6b6b1e5463
--- /dev/null
+++ b/reco/KF/CbmKFV0FinderQa.cxx
@@ -0,0 +1,61 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmKFV0FinderQa.cxx
+/// \brief A simple QA for the V0 finder class (runs together with the finder task) (implementation)
+/// \since 16.01.2025
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#include "CbmKFV0FinderQa.h"
+
+#include <Logger.h>
+
+#include <TFile.h>
+
+using cbm::kfp::V0FinderQa;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderQa::InitHistograms()
+{
+ TString n{""}; // name of histogram
+ TString t{""}; // title of histogram
+
+ n = "tof_lst_hit_time";
+ t = "Time of last TOF hit in track;t [ns];Counts";
+ fph_tof_lst_hit_time = MakeQaObject<TH1D>(n, t, 200, -100., 100.);
+
+ n = "beta_all";
+ t = "Speed of tracks (all);#beta;Counts";
+ fph_beta_all = MakeQaObject<TH1D>(n, t, 150, 0., 1.5);
+
+ {
+ n = "dca";
+ t = "DCA to the origin of the selected tracks;dca [cm];Counts";
+ fph_dca = MakeQaObject<TH1D>(n, t, 120, 0., 6.);
+
+ n = "dca2D";
+ t = "DCA to the origin of the selected tracks;dca_{x} [cm];dca_{y} [cm];Counts";
+ fph_dca2D = MakeQaObject<TH2D>(n, t, 240, -6., 6., 240, -6., 6.);
+
+ n = "dca_projectionX";
+ t = "DCA to the origin of the selected tracks;dca_{x} [cm];Counts";
+ fph_dca_projectionX = MakeQaObject<TH1D>(n, t, 240, -6., 6.);
+
+ n = "dca_projectionY";
+ t = "DCA to the origin of the selected tracks;dca_{y} [cm];Counts";
+ fph_dca_projectionY = MakeQaObject<TH1D>(n, t, 240, -6., 6.);
+ }
+}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderQa::WriteHistograms(const TString& fileName)
+{
+ LOG(info) << "V0FinderQa: writing histograms to file: " << fileName;
+ TFile fileOut(fileName, "RECREATE");
+ this->WriteToFile(&fileOut);
+ fileOut.Close();
+}
diff --git a/reco/KF/CbmKFV0FinderQa.h b/reco/KF/CbmKFV0FinderQa.h
new file mode 100644
index 0000000000..53091e0112
--- /dev/null
+++ b/reco/KF/CbmKFV0FinderQa.h
@@ -0,0 +1,66 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmKFV0FinderQa.h
+/// \brief A simple QA for the V0 finder class (runs together with the finder task)
+/// \since 15.01.2025
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+// FIXME: Later move this class into a separate QA-task
+
+#pragma once
+
+#include "CbmQaIO.h"
+
+#include <TString.h>
+
+#include <array>
+#include <string>
+#include <vector>
+
+class TH1D;
+class TH2D;
+
+namespace cbm::kfp
+{
+ /// \class V0FinderQa
+ /// \brief A simple QA for the V0 finder
+ class V0FinderQa : public CbmQaIO {
+ public:
+ /// \brief Constructor from parameters
+ /// \param bUseMc If MC truth should be used
+ V0FinderQa(bool bUseMc) : CbmQaIO("v0Finder"), fbUseMc(bUseMc) {}
+
+ /// \brief Destructor
+ ~V0FinderQa() = default;
+
+ V0FinderQa(const V0FinderQa&) = delete;
+ V0FinderQa(V0FinderQa&&) = delete;
+ V0FinderQa& operator=(const V0FinderQa&) = delete;
+ V0FinderQa& operator=(V0FinderQa&&) = delete;
+
+ /// \brief Initializes histograms
+ void InitHistograms();
+
+ /// \brief Writes histograms into file
+ void WriteHistograms(const TString&);
+
+ //* Histograms (public for easier access)
+ TH1D* fph_tof_lst_hit_time{nullptr}; ///< Time of tof hits, used for track momentum estimation
+ TH1D* fph_beta_all{nullptr}; ///< Beta of tracks
+ TH1D* fph_dca{nullptr}; ///< Track DCA to origin
+ TH2D* fph_dca2D{nullptr}; ///< Track DCA to origin (2D)
+ TH1D* fph_dca_projectionX{nullptr}; ///< Track DCA to origin (x-component)
+ TH1D* fph_dca_projectionY{nullptr}; ///< Track DCA to origin (y-component)
+ TH2D* fph_track_rapidity_vs_pt_all{nullptr}; ///< Phase space of all accepted tracks
+ TH2D* fph_track_rapidity_vs_pt_pion{nullptr}; ///< Phase space of pion candidates
+ TH2D* fph_track_rapidity_vs_pt_proton{nullptr}; ///< Phase space of proton candidates
+
+
+ private:
+ bool fbUseMc{false};
+
+ ClassDef(V0FinderQa, 0);
+ };
+} // namespace cbm::kfp
diff --git a/reco/KF/CbmKFV0FinderTask.cxx b/reco/KF/CbmKFV0FinderTask.cxx
new file mode 100644
index 0000000000..2ec69f809c
--- /dev/null
+++ b/reco/KF/CbmKFV0FinderTask.cxx
@@ -0,0 +1,774 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmKFV0FinderTask.cxx
+/// \brief A FairTask for V0 candidates finding in mCBM (implementation)
+/// \since 10.01.2025
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#include "CbmKFV0FinderTask.h"
+
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmKfTarget.h"
+#include "CbmStsHit.h"
+#include "CbmStsTrack.h"
+#include "CbmTofAddress.h"
+#include "CbmTofHit.h"
+#include "CbmTofTrack.h"
+#include "CbmTrdHit.h"
+#include "CbmTrdTrack.h"
+#include "FairRunAna.h"
+#include "KFPTrackVector.h"
+#include "KFVertex.h"
+#include "Logger.h"
+#include "TClonesArray.h"
+#include "TMatrixTSym.h"
+
+#include <algorithm>
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <sstream>
+
+// Log macros:
+#define ERR_ LOG(error) << fName << ": "
+#define LOG_(SEVERITY, VERBOSITY) LOG_IF(SEVERITY, fVerbose >= VERBOSITY) << fName << ": "
+
+using cbm::kfp::V0FinderTask;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool V0FinderTask::AssignMomentum(CbmGlobalTrack* pTrack, int pdg)
+{
+ const int iTofTrk = pTrack->GetTofTrackIndex();
+ if (iTofTrk < 0) {
+ ++fCounters[ECounter::TracksWoTofHits];
+ return false; // Skip tracks, which do not have hits in TOF
+ }
+
+ const auto* pTofTrack{static_cast<const CbmTofTrack*>(fpBrTofTracks->At(iTofTrk))};
+ assert(pTofTrack);
+
+ int nTofHits = pTofTrack->GetNofTofHits();
+ if (nTofHits < 1) {
+ // Must not be called
+ return false;
+ }
+ ++fCounters[ECounter::TracksWAtLeastOneTofHit];
+
+ //int iFstTofHit = pTofTrack->GetTofHitIndex(0);
+ int iLstTofHit = pTofTrack->GetTofHitIndex(nTofHits - 1); // last hit in TOF
+
+ const auto* pLstTofHit{static_cast<CbmTofHit*>(fpBrTofHits->At(iLstTofHit))};
+ if (fbRunQa) {
+ fpQa->fph_tof_lst_hit_time->Fill(pLstTofHit->GetTime());
+ }
+ assert(pLstTofHit);
+
+ if (pLstTofHit->GetTime() < 0) { //< wrongly calibrated T0
+ ++fCounters[ECounter::TracksWNegativeTofHitTime];
+ return false;
+ }
+
+ auto qpAndBeta = EstimateQp(pLstTofHit, pdg);
+ if (fbRunQa) {
+ fpQa->fph_beta_all->Fill(qpAndBeta.fBeta);
+ }
+
+ if (qpAndBeta.fBeta > 1) {
+ ++fCounters[ECounter::TracksWithUnphysicalBeta];
+ return false; // unphysical track
+ }
+
+ // Update the track with the PID hypothesis and momentum
+ FairTrackParam parFst(*pTrack->GetParamFirst());
+ FairTrackParam parLst(*pTrack->GetParamLast());
+
+ parFst.SetQp(qpAndBeta.fQp);
+ parLst.SetQp(qpAndBeta.fQp);
+ parFst.SetCovariance(4, 4, qpAndBeta.fQpVar);
+ parFst.SetCovariance(4, 4, qpAndBeta.fQpVar);
+
+ pTrack->SetParamFirst(&parFst);
+ pTrack->SetParamLast(&parLst);
+ return true;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<bool UseEvent>
+bool V0FinderTask::ProcessEvent(const CbmEvent* pEvent)
+{
+ int nTracks = UseEvent ? pEvent->GetNofData(ECbmDataType::kGlobalTrack) : fpBrGlobalTracks->GetEntriesFast();
+ fCounters[ECounter::TracksTotal] += nTracks;
+
+
+ // Local event counters
+ int nProtons{0};
+ int nPions{0};
+ int nTracksSelected{0};
+
+ // Local containers (per event)
+ std::vector<const CbmGlobalTrack*> vpSelectedTracks; // Selected tracks [n selected tracks]
+ vpSelectedTracks.reserve(nTracks);
+ std::vector<int> vSelectedTrackIds; // Indices of selected tracks [n selected tracks]
+ vSelectedTrackIds.reserve(nTracks);
+
+ // ----- Shift TOF hit times to t0
+ if constexpr (UseEvent) {
+ ShiftTofHitsToTzero(pEvent);
+ }
+
+ // ----- Select tracks
+ for (int iTrkEvt{0}; iTrkEvt < nTracks; ++iTrkEvt) {
+ const int iTrk = UseEvent ? pEvent->GetIndex(ECbmDataType::kGlobalTrack, iTrkEvt) : iTrkEvt;
+ auto* pGlobalTrack{static_cast<CbmGlobalTrack*>(fpBrGlobalTracks->At(iTrk))};
+ if (!SelectTrack(pGlobalTrack, iTrk)) {
+ continue;
+ }
+ ++nTracksSelected;
+ switch (pGlobalTrack->GetPidHypo()) {
+ case -211: // pions
+ ++nPions;
+ break;
+ case 2212: // protons
+ ++nProtons;
+ break;
+ }
+ vpSelectedTracks.push_back(pGlobalTrack);
+ vSelectedTrackIds.push_back(iTrk);
+ } // END LOOP: tracks in event
+
+ if (nProtons < 1 || nPions < 1) {
+ return false; // Nothing to search for in the event
+ }
+
+ fCounters[ECounter::TracksSelected] += nTracksSelected;
+ fCounters[ECounter::Pions] += nPions;
+ fCounters[ECounter::Protons] += nProtons;
+
+ ++fCounters[ECounter::EventsLambdaCand];
+
+ // NOTE: The variable is used to assign a KfParticle to a corresponding primary vertex:
+ // a) pvChi2 < 3 => primary particle
+ // b) pvChi2 > 3 => secondary particle
+ // Here we assign all particles as secondary (TODO: test and try to separate primary and secondary particles)
+ std::vector<float> vChi2ToPv; // Chi-square to primary vertex (NOTE: not used)
+ vChi2ToPv.resize(nTracksSelected, 999.); // Assign all tracks as secondary
+
+ KFPTrackVector kfpTracksFst = MakeKfpTrackVector(vpSelectedTracks, vSelectedTrackIds, vChi2ToPv, true);
+ KFPTrackVector kfpTracksLst = MakeKfpTrackVector(vpSelectedTracks, vSelectedTrackIds, vChi2ToPv, false);
+
+ // ----- Reconstruct topology in the event
+ fpTopoReconstructorEvent->Clear();
+ fpTopoReconstructorEvent->Init(kfpTracksFst, kfpTracksLst);
+ if (EPvUsageMode::Target == fPvUsageMode) {
+ float x{float(fpOrigin->GetX())};
+ float y{float(fpOrigin->GetY())};
+ float z{float(fpOrigin->GetZ())};
+ // NOTE: the target size can be applied as a covariance matrix, now it has all nulls
+ fpTopoReconstructorEvent->AddPV(MakeKfpPrimaryVertex(x, y, z), std::vector<int>{});
+ }
+ fpTopoReconstructorEvent->SortTracks();
+ fpTopoReconstructorEvent->ReconstructParticles();
+
+ // ----- Count number of found lambda-candidates
+ const auto& particles = fpTopoReconstructorEvent->GetParticles();
+ int nLambda = std::count_if(particles.begin(), particles.end(), [](const auto& p) { return p.GetPDG() == 3122; });
+ fCounters[ECounter::KfpLambdaCandidates] += nLambda;
+ if (nLambda > 0) {
+ ++fCounters[ECounter::KfpEventsLambdaCand];
+ }
+
+ // ----- Store particles to the run topology
+ StoreParticles();
+
+ return true;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool V0FinderTask::CheckTrackParam(const FairTrackParam* pParam)
+{
+ bool bOk{true};
+ bOk &= std::isfinite(pParam->GetX());
+ bOk &= std::isfinite(pParam->GetY());
+ bOk &= std::isfinite(pParam->GetZ());
+ bOk &= std::isfinite(pParam->GetTx());
+ bOk &= std::isfinite(pParam->GetTy());
+ bOk &= std::isfinite(pParam->GetQp());
+ std::array<double, 15> covMatrix = {0.};
+ if (bOk) {
+ for (int i = 0, iCov = 0; i < 5; ++i) {
+ for (int j = i; j < 5; j++, iCov++) {
+ // NOTE: In the FairTrackParam::GetCovariance(i, j) one can exchange the i and j indices, but for more
+ // efficient calculation it is recommended to keep i < j.
+ covMatrix[iCov] = pParam->GetCovariance(i, j);
+ bOk &= std::isfinite(covMatrix[iCov]);
+ }
+ }
+ }
+ return bOk;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+V0FinderTask::DcaVector V0FinderTask::EstimateDcaToOrigin(const CbmStsTrack* pStsTrack) const
+{
+ DcaVector res;
+ if (pStsTrack->GetNofStsHits() < 2) {
+ // Too few STS hits
+ return res;
+ }
+ const auto* pHitFst{static_cast<const CbmStsHit*>(fpBrStsHits->At(pStsTrack->GetStsHitIndex(0)))}; // first hit
+ const auto* pHitSnd{static_cast<const CbmStsHit*>(fpBrStsHits->At(pStsTrack->GetStsHitIndex(1)))}; // second hit
+ assert(pHitFst);
+ assert(pHitSnd);
+ double factor{(pHitFst->GetZ() - fpOrigin->GetZ()) / (pHitSnd->GetZ() - pHitFst->GetZ())};
+ res.fX = pHitFst->GetX() - fpOrigin->GetX() - factor * (pHitSnd->GetX() - pHitFst->GetX());
+ res.fY = pHitFst->GetY() - fpOrigin->GetY() - factor * (pHitSnd->GetY() - pHitFst->GetY());
+ res.fAbs = std::sqrt(res.fX * res.fX + res.fY * res.fY);
+ res.fX /= res.fAbs;
+ res.fY /= res.fAbs;
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+V0FinderTask::QpAndBeta V0FinderTask::EstimateQp(const CbmTofHit* pTofHit, int pdg) const
+{
+ // FIXME: For now we consider the origin as a precise measurement, which is not correct. The corresponding
+ // uncertainty of the origin must be accounted here.
+ // Also: redefine origin with reconstructed or MC PV (in other function).
+ // Also: the t0 uncertainty should be estimated and accounted in the momentum uncertainty
+ QpAndBeta res;
+ double x{pTofHit->GetX() - fpOrigin->GetX()};
+ double y{pTofHit->GetY() - fpOrigin->GetY()};
+ double z{pTofHit->GetZ() - fpOrigin->GetZ()};
+ double x2{x * x};
+ double y2{y * y};
+ double z2{z * z};
+ double t2{pTofHit->GetTime() * pTofHit->GetTime()};
+ double r2{x2 + y2 + z2};
+ double r4{r2 * r2};
+ // FIXME: In general case, the x, y, z coordinates of the origin can have non-zero covariances, which should be
+ // accounted in the calculation.
+ double xVar{pTofHit->GetDx() * pTofHit->GetDx() + fpOrigin->GetCovariance(0, 0)};
+ double yVar{pTofHit->GetDy() * pTofHit->GetDy() + fpOrigin->GetCovariance(1, 1)};
+ double zVar{pTofHit->GetDz() * pTofHit->GetDz() + fpOrigin->GetCovariance(2, 2)};
+ double tVar{pTofHit->GetTimeError() * pTofHit->GetTimeError()};
+ double factor{(x2 * xVar + y2 * yVar + z2 * zVar) / r4 + tVar / t2};
+ res.fBeta = std::sqrt(r2) / (pTofHit->GetTime() * kSpeedOfLight);
+ res.fBetaVar = res.fBeta * factor;
+ if (res.fBeta >= 1) {
+ return res; //< Non-physical beta
+ // FIXME: I would set beta to fixed large value, like 0.999999 for this case
+ }
+ double gamma{1. / std::sqrt(1. - res.fBeta * res.fBeta)};
+ double m{0};
+ double q{0};
+ switch (pdg) {
+ case -211:
+ m = kPionMass;
+ q = -1;
+ break;
+ case 2212:
+ m = kProtonMass;
+ q = 1;
+ break;
+ default:
+ res.fQp = -1. / 99999.; //< Give small momentum to primary tracks
+ res.fQpVar = res.fQp * res.fQp * fQpAssignedUncertainty * fQpAssignedUncertainty;
+ return res;
+ // NOTE: it would be more straight forward to assign the primary tracks a mass of pion and keep it
+ }
+ res.fQp = q / (gamma * res.fBeta * m);
+ if (fQpAssignedUncertainty > 0) {
+ res.fQpVar = res.fQp * res.fQp * fQpAssignedUncertainty * fQpAssignedUncertainty;
+ }
+ else {
+ res.fQpVar = q * q / (m * m * res.fBeta * res.fBeta) * factor;
+ }
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderTask::Exec(Option_t*)
+{
+ fvTrackDca.clear();
+ fvTrackDca.resize(fpBrGlobalTracks->GetEntriesFast());
+
+ // ----- Reset data containers per timeslice
+ fpTopoReconstructorRun->Clear();
+
+ // ----- Process timeslice
+ if (EProcessingMode::EventBased == fProcessingMode) {
+ ++fCounters[ECounter::EventsTotal];
+ ProcessEvent</*UseEvent = */ false>(nullptr);
+ }
+ else {
+ const int nEvents{fpBrRecoEvents->GetEntriesFast()};
+ fCounters[ECounter::EventsTotal] += nEvents;
+ for (int iEvent = 0; iEvent < nEvents; ++iEvent) {
+ const auto* pEvent = static_cast<CbmEvent*>(fpBrRecoEvents->At(iEvent));
+ ProcessEvent</*UseEvent = */ true>(pEvent);
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderTask::Finish()
+{
+ // ----- Store QA histograms
+ if (fbRunQa) {
+ fpQa->WriteHistograms(fsQaOutputName);
+ }
+
+ // ----- Counters summary log
+ std::stringstream msg;
+ msg << "*** Counters summary ***\n";
+ msg << "\n tracks total: " << fCounters[ECounter::TracksTotal];
+ msg << "\n tracks selected: " << fCounters[ECounter::TracksSelected];
+ msg << "\n tracks w/ infinit parameters: " << fCounters[ECounter::TracksInfiniteParam];
+ msg << "\n tracks w/o TOF hits: " << fCounters[ECounter::TracksWoTofHits];
+ msg << "\n tracks w/ last TOF hit having t<0: " << fCounters[ECounter::TracksWNegativeTofHitTime];
+ msg << "\n tracks w/o STS hits: " << fCounters[ECounter::TracksWoStsHits];
+ msg << "\n tracks w/o PID: " << fCounters[ECounter::TracksWoPid];
+ msg << "\n tracks w/o momentum: " << fCounters[ECounter::TracksWoMomentum];
+ msg << "\n tracks w/ at least one TOF hit: " << fCounters[ECounter::TracksWAtLeastOneTofHit];
+ msg << "\n tracks w/ at least two TOF hits: " << fCounters[ECounter::TracksWAtLeastTwoTofHits];
+ msg << "\n tracks w/ beta > 1 " << fCounters[ECounter::TracksWithUnphysicalBeta];
+ msg << "\n pion candidates: " << fCounters[ECounter::Pions];
+ msg << "\n proton candidates: " << fCounters[ECounter::Protons];
+ msg << "\n events total: " << fCounters[ECounter::EventsTotal];
+ msg << "\n events w/ potential Lambda-candidate: " << fCounters[ECounter::EventsLambdaCand];
+ msg << "\n events w/ Lambda-candidate from KFPF: " << fCounters[ECounter::KfpEventsLambdaCand];
+ msg << "\n Lambda-candidates: " << fCounters[ECounter::KfpLambdaCandidates];
+ LOG_(info, 1) << msg.str();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+int V0FinderTask::InferTrackPidTopo(double dca) const
+{
+ if (std::isnan(dca)) {
+ return -2; //< Track PID cannot be infered
+ }
+ if (dca > fMinPionDca) {
+ return -211; // pi-
+ }
+ else if (dca > fMinProtonDca) {
+ return 2212; // proton
+ }
+ return kPrimaryPdg;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus V0FinderTask::Init()
+{
+ LOG_(info, 1) << "initializing the task ... ";
+
+ LOG_(info, 1) << "using the following configuration:\n processing mode: " << ToString(fProcessingMode)
+ << "\n PID approach: " << ToString(fPidApproach) << "\n " << ToString(fPvUsageMode);
+
+ // Temporary not implemented:
+ if (EPidApproach::Mc == fPidApproach || EPvUsageMode::Mc == fPvUsageMode
+ || EPvUsageMode::Reconstructed == fPvUsageMode || EPvUsageMode::ReconstructSingle == fPvUsageMode
+ || EPvUsageMode::ReconstructMultiple == fPvUsageMode) {
+ ERR_ << "Desirable configuration of PID or PV handling was not implemented yet";
+ return kFATAL;
+ }
+
+ // ----- Input data branches initialization
+ const auto* pTarget = kf::Target::Instance(); // CBM target info
+
+ auto* pFairManager = FairRootManager::Instance();
+ if (!pFairManager) {
+ ERR_ << "FairRootManager not found";
+ return kERROR;
+ }
+
+ auto InitBranch = [&](TClonesArray*& branch, const char* name) -> bool {
+ branch = dynamic_cast<TClonesArray*>(pFairManager->GetObject(name));
+ if (!branch) {
+ ERR_ << "branch \"" << name << "\" not found";
+ return false;
+ }
+ return true;
+ };
+
+ bool bBranchesInitialized{true};
+ if (EProcessingMode::TimeBased == fProcessingMode) {
+ bBranchesInitialized = InitBranch(fpBrRecoEvents, "CbmEvent") && bBranchesInitialized;
+ }
+ bBranchesInitialized = InitBranch(fpBrGlobalTracks, "GlobalTrack") && bBranchesInitialized;
+ bBranchesInitialized = InitBranch(fpBrStsTracks, "StsTrack") && bBranchesInitialized;
+ bBranchesInitialized = InitBranch(fpBrTrdTracks, "TrdTrack") && bBranchesInitialized;
+ bBranchesInitialized = InitBranch(fpBrTofTracks, "TofTrack") && bBranchesInitialized;
+ bBranchesInitialized = InitBranch(fpBrStsHits, "StsHit") && bBranchesInitialized;
+ bBranchesInitialized = InitBranch(fpBrTrdHits, "TrdHit") && bBranchesInitialized;
+ bBranchesInitialized = InitBranch(fpBrTofHits, "TofHit") && bBranchesInitialized;
+
+ if (EPvUsageMode::Reconstructed == fPvUsageMode) {
+ fpBrPrimaryVertex = dynamic_cast<CbmVertex*>(pFairManager->GetObject("PrimaryVertex."));
+ if (!fpBrPrimaryVertex) {
+ ERR_ << "branch \"PrimaryVertex.\" not found";
+ bBranchesInitialized = false;
+ }
+ }
+ else if (EPvUsageMode::Target == fPvUsageMode) {
+ double x{pTarget->GetX()};
+ double y{pTarget->GetY()};
+ double z{pTarget->GetZ()};
+ LOG_(info, 1) << "using target center as origin: r = (" << x << ", " << y << ", " << z << ") [cm]";
+ TMatrixFSym covMatrix(3);
+ covMatrix(1, 0) = 0.;
+ covMatrix(2, 0) = 0.;
+ covMatrix(2, 1) = 0.;
+ //double transverseVariance{pTarget->GetRmax() * pTarget->GetRmax() / 12.};
+ //covMatrix(0, 0) = transverseVariance;
+ //covMatrix(1, 1) = transverseVariance;
+ //covMatrix(2, 2) = pTarget->GetDz() * pTarget->GetDz() / 3.;
+ covMatrix(0, 0) = 0.;
+ covMatrix(1, 1) = 0.;
+ covMatrix(2, 2) = 0.;
+ fpOrigin->SetVertex(x, y, z, 0., 1, 0, covMatrix);
+ }
+
+ if (!bBranchesInitialized) {
+ return kERROR;
+ }
+
+ // ----- Topology reconstructor initialization
+ fpTopoReconstructorRun->SetTarget({float(pTarget->GetX()), float(pTarget->GetY()), float(pTarget->GetZ())});
+ fpTopoReconstructorEvent->SetTarget(fpTopoReconstructorRun->GetTargetPosition());
+ fpTopoReconstructorEvent->CopyCuts(fpTopoReconstructorRun.get());
+ fpTopoReconstructorEvent->GetKFParticleFinder()->SetReconstructionList(
+ fpTopoReconstructorRun->GetKFParticleFinder()->GetReconstructionList());
+
+ // ----- Auxilary variables initialization
+ std::fill(fCounters.begin(), fCounters.end(), 0);
+
+ // ----- QA initialization
+ if (fbRunQa) {
+ fpQa = std::make_unique<V0FinderQa>(fbUseMc);
+ fpQa->InitHistograms();
+ }
+
+ LOG_(info, 1) << "initializing the task ... done";
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+KFPTrackVector V0FinderTask::MakeKfpTrackVector(const std::vector<const CbmGlobalTrack*>& vpTracks,
+ const std::vector<int>& vTrackIds, const std::vector<float>& vChi2ToPv,
+ bool bAtFirstPoint) const
+{
+ // TODO: Cross-check with the CbmKFParticleFinder results
+ // NOTE: Actually tracks themselves are not needed here, so one could fill and pass a vector of FairTrackParam*,
+ // put in KFPTrackVector there is an option to set number of Pixel hits, so we will keep the signature of
+ // the function as it is.
+ KFPTrackVector trackVector;
+ int nTracks = vpTracks.size();
+ trackVector.Resize(nTracks);
+ for (int iTrk = 0; iTrk < nTracks; ++iTrk) {
+ const auto* pTrack{vpTracks[iTrk]};
+ const auto* pParam{bAtFirstPoint ? pTrack->GetParamFirst() : pTrack->GetParamLast()};
+
+ // ----- Parameters definition
+ int pdg{pTrack->GetPidHypo()};
+ const double& tx{pParam->GetTx()};
+ const double& ty{pParam->GetTy()};
+ const double& qp{pParam->GetQp()};
+ int q{qp > 0. ? 1 : -1};
+ if (std::fabs(pdg) == 10000020030 || std::fabs(pdg == 1000020040)) { // He3 and He4
+ // NOTE: at the moment not called, but never the less leave it here for the future pid procedures
+ q *= 2;
+ }
+ double p{q / qp};
+ double p2{p * p};
+ double t2inv{1. / (1. + tx * tx + ty * ty)};
+ double pz{std::sqrt(t2inv * p2)};
+ double px{tx * pz};
+ double py{ty * pz};
+
+ trackVector.SetParameter(pParam->GetX(), 0, iTrk);
+ trackVector.SetParameter(pParam->GetY(), 1, iTrk);
+ trackVector.SetParameter(pParam->GetZ(), 2, iTrk);
+ trackVector.SetParameter(px, 3, iTrk);
+ trackVector.SetParameter(py, 4, iTrk);
+ trackVector.SetParameter(pz, 5, iTrk);
+
+ // Jacobian matrix for (tx, ty, qp) -> (px, py, pz)
+ std::array<std::array<double, 3>, 3> Jp;
+ Jp[2][0] = -t2inv * px; // d(pz)/d(tx)
+ Jp[2][1] = -t2inv * py; // d(pz)/d(ty)
+ Jp[2][2] = -pz / qp; // d(pz)/d(qp)
+ Jp[0][0] = tx * Jp[2][0] + pz; // d(px)/d(tx)
+ Jp[0][1] = tx * Jp[2][1]; // d(px)/d(ty)
+ Jp[0][2] = tx * Jp[2][2]; // d(px)/d(qp)
+ Jp[1][0] = ty * Jp[2][0]; // d(py)/d(tx)
+ Jp[1][1] = ty * Jp[2][1] + pz; // d(py)/d(ty)
+ Jp[1][2] = ty * Jp[2][2]; // d(py)/d(qp)
+
+ // ----- Covariance matrix definition
+ // Position covariances
+ trackVector.SetCovariance(pParam->GetCovariance(0, 0), 0, iTrk); // var(x)
+ trackVector.SetCovariance(pParam->GetCovariance(0, 1), 1, iTrk); // cov(x, y)
+ trackVector.SetCovariance(pParam->GetCovariance(1, 1), 2, iTrk); // var(y)
+
+ // Momentum-position covariances
+ auto MomPosCovariance = [&](const int k, const int l) constexpr->double
+ {
+ double val{0.};
+ const auto& JpA = Jp[k];
+ for (int i = 0; i < 3; ++i) {
+ val += JpA[i] * pParam->GetCovariance(i + 2, l);
+ }
+ return val;
+ };
+ trackVector.SetCovariance(MomPosCovariance(0, 0), 6, iTrk); // cov(x, px)
+ trackVector.SetCovariance(MomPosCovariance(0, 1), 7, iTrk); // cov(y, px)
+ trackVector.SetCovariance(MomPosCovariance(1, 0), 10, iTrk); // cov(x, py)
+ trackVector.SetCovariance(MomPosCovariance(1, 1), 11, iTrk); // cov(y, py)
+ trackVector.SetCovariance(MomPosCovariance(2, 0), 15, iTrk); // cov(x, pz)
+ trackVector.SetCovariance(MomPosCovariance(2, 1), 16, iTrk); // cov(y, pz)
+
+ // Momentum covariances
+ auto MomentumCovariance = [&](const int k, const int l) constexpr->double
+ {
+ double val{0.};
+ const auto& JpA = Jp[k];
+ const auto& JpB = Jp[l];
+ for (int i = 0; i < 3; ++i) {
+ double factor{0.};
+ for (int j = 0; j < 3; ++j) {
+ factor += JpB[j] * pParam->GetCovariance(i + 2, j + 2);
+ }
+ val += JpA[i] * factor;
+ }
+ return val;
+ };
+ trackVector.SetCovariance(MomentumCovariance(0, 0), 9, iTrk); // var(px)
+ trackVector.SetCovariance(MomentumCovariance(1, 0), 13, iTrk); // cov(px, py)
+ trackVector.SetCovariance(MomentumCovariance(1, 1), 14, iTrk); // var(py)
+ trackVector.SetCovariance(MomentumCovariance(2, 0), 18, iTrk); // cov(px, pz)
+ trackVector.SetCovariance(MomentumCovariance(2, 1), 19, iTrk); // cov(py, pz)
+ trackVector.SetCovariance(MomentumCovariance(2, 2), 20, iTrk); // var(pz)
+
+ // Zero covariances
+ // NOTE: from the tracking point of view a z-coordinate in known precisely, so all the corresponding covariances
+ // should be set to null.
+ trackVector.SetCovariance(0.f, 3, iTrk); // cov(x,z)
+ trackVector.SetCovariance(0.f, 4, iTrk); // cov(y,z)
+ trackVector.SetCovariance(0.f, 5, iTrk); // var(z)
+ trackVector.SetCovariance(0.f, 8, iTrk); // cov(z,px)
+ trackVector.SetCovariance(0.f, 12, iTrk); // cov(z,py)
+ trackVector.SetCovariance(0.f, 17, iTrk); // var(z,pz)
+
+ // ----- Other parameters
+ // Magnetic field
+ // NOTE: In mCBM there is no magnetic field, so here we define the coefficients with zeros
+ for (int iF = 0; iF < 10; ++iF) {
+ trackVector.SetFieldCoefficient(0.f, iF, iTrk);
+ }
+
+ trackVector.SetId(vTrackIds[iTrk], iTrk);
+ trackVector.SetPDG(pdg, iTrk);
+ trackVector.SetQ(q, iTrk);
+ trackVector.SetNPixelHits(0, iTrk); // Number of MVD hits in track (used to define PV inside the KFPFinder)
+
+ if (EPvUsageMode::Reconstructed == fPvUsageMode || EPvUsageMode::Mc == fPvUsageMode) {
+ if (vChi2ToPv[iTrk] < kChi2PvPrimThrsh) {
+ trackVector.SetPVIndex(0, iTrk); // Primary track
+ }
+ else {
+ trackVector.SetPVIndex(-1, iTrk); // Secondary track
+ }
+ }
+ else {
+ trackVector.SetPVIndex(-1, iTrk); // Secondary track
+ }
+ }
+ return trackVector;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+KFVertex V0FinderTask::MakeKfpPrimaryVertex(float x, float y, float z) const
+{
+ KFVertex kfVertex;
+ kfVertex.X() = x;
+ kfVertex.Y() = y;
+ kfVertex.Z() = z;
+ for (int iC = 0; iC < 6; ++iC) {
+ kfVertex.Covariance(iC) = 0.f;
+ }
+ kfVertex.Chi2() = -100.f;
+ return kfVertex;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool V0FinderTask::SelectTrack(CbmGlobalTrack* pTrack, int iTrk)
+{
+ const int iStsTrk = pTrack->GetStsTrackIndex();
+ if (iStsTrk < 0) {
+ ++fCounters[ECounter::TracksWoStsHits];
+ return false; // Skip tracks, which do not have hits in TOF
+ }
+ const auto* pStsTrack{static_cast<const CbmStsTrack*>(fpBrStsTracks->At(iStsTrk))};
+
+ // Get PID hypothesis for the track
+ int pdgHypo{kUndefPdg};
+ if (EPidApproach::Topo == fPidApproach) {
+ auto& trkDca{fvTrackDca[iTrk]};
+ trkDca = EstimateDcaToOrigin(pStsTrack);
+ pdgHypo = InferTrackPidTopo(trkDca.fAbs);
+ if (fbRunQa) {
+ fpQa->fph_dca->Fill(trkDca.fAbs);
+ fpQa->fph_dca2D->Fill(trkDca.fAbs * trkDca.fX, trkDca.fAbs * trkDca.fY);
+ fpQa->fph_dca_projectionX->Fill(trkDca.fAbs * trkDca.fX);
+ fpQa->fph_dca_projectionY->Fill(trkDca.fAbs * trkDca.fY);
+ }
+ }
+
+ pTrack->SetPidHypo(pdgHypo);
+ if (pdgHypo == kUndefPdg) {
+ ++fCounters[ECounter::TracksWoPid];
+ return false; // Tracks, which do not obey the PID criteria
+ }
+
+ // Assign a momentum to the track using the last hit in TOF, reject the track, if it has no TOF hits
+ if (!AssignMomentum(pTrack, pdgHypo)) {
+ ++fCounters[ECounter::TracksWoMomentum];
+ return false;
+ }
+
+ // Check track parameters for validity
+ // FIXME: understand the reason, why the parameters are infinite
+ if (!CheckTrackParam(pTrack->GetParamFirst())) {
+ ++fCounters[ECounter::TracksInfiniteParam];
+ return false;
+ }
+
+ return true;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+double V0FinderTask::ShiftTofHitsToTzero(const CbmEvent* pEvent)
+{
+ // Define t0 from the Bmon hits
+ double t0{std::numeric_limits<double>::signaling_NaN()};
+ int nTofHits = pEvent->GetNofData(ECbmDataType::kTofHit);
+ for (int iTofHitEvt{0}; iTofHitEvt < nTofHits; ++iTofHitEvt) {
+ int iTofHit = pEvent->GetIndex(ECbmDataType::kTofHit, iTofHitEvt);
+ const auto* pTofHit{static_cast<const CbmTofHit*>(fpBrTofHits->At(iTofHit))};
+ //if (5 == CbmTofAddress::GetSmType(pTofHit->GetAddress())) { // selection by the supermodule type
+ if (pTofHit->GetZ() == 0) { // Take some small z-coordinate to identify BMon (TODO: provide a flag by some task)
+ t0 = pTofHit->GetTime();
+ }
+ }
+ // NOTE: t0 must be defined for each event, since the Bmon digis are used to seed the digi event builder. Basically,
+ // the tZero must be defined as a field of CbmEvent, e.g. in a separate FairTask, or directly
+ // in CbmAlgoBuildRawEvent
+ assert(!std::isnan(t0));
+
+ // Shift TOF times to found t0:
+ for (int iTofHitEvt{0}; iTofHitEvt < nTofHits; ++iTofHitEvt) {
+ int iTofHit = pEvent->GetIndex(ECbmDataType::kTofHit, iTofHitEvt);
+ auto* pTofHit{static_cast<CbmTofHit*>(fpBrTofHits->At(iTofHit))};
+ pTofHit->SetTime(pTofHit->GetTime() - t0 - fTzeroOffset);
+ }
+ return t0;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderTask::StoreParticles()
+{
+ int indexOffset = fpTopoReconstructorRun->GetParticles().size();
+ for (int iPv = 0; iPv < fpTopoReconstructorEvent->NPrimaryVertices(); ++iPv) {
+ std::vector<int> vPVTrackIndexArray = fpTopoReconstructorEvent->GetPVTrackIndexArray(iPv);
+ std::for_each(vPVTrackIndexArray.begin(), vPVTrackIndexArray.end(), [&](auto& item) { item += indexOffset; });
+ fpTopoReconstructorRun->AddPV(fpTopoReconstructorEvent->GetPrimKFVertex(iPv), vPVTrackIndexArray);
+ }
+
+ for (size_t iP = 0; iP < fpTopoReconstructorEvent->GetParticles().size(); ++iP) {
+ const KFParticle& particleEvt{fpTopoReconstructorEvent->GetParticles()[iP]};
+ KFParticle particle{particleEvt};
+ particle.CleanDaughtersId();
+ int nDaughters{particleEvt.NDaughters()};
+ if (nDaughters == 1) {
+ particle.AddDaughterId(particleEvt.DaughterIds()[0]);
+ }
+ else if (nDaughters > 1) {
+ for (int iD = 0; iD < particleEvt.NDaughters(); ++iD) {
+ particle.AddDaughterId(particleEvt.DaughterIds()[iD] + indexOffset);
+ }
+ }
+ particle.SetId(particleEvt.Id() + indexOffset);
+ fpTopoReconstructorRun->AddParticle(particle);
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string V0FinderTask::ToString(EProcessingMode mode)
+{
+ switch (mode) {
+ case EProcessingMode::EventBased: return "event-based";
+ case EProcessingMode::TimeBased: return "time-based";
+ default: return "";
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string V0FinderTask::ToString(EPidApproach pid)
+{
+ switch (pid) {
+ case EPidApproach::Topo: return "track topology";
+ case EPidApproach::Mc: return "MC-truth";
+ default: return "";
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string V0FinderTask::ToString(EPvUsageMode pvMode)
+{
+ switch (pvMode) {
+ case EPvUsageMode::Reconstructed: return "reconstructed primary vertex is used";
+ case EPvUsageMode::ReconstructSingle: return "a single primary vertex will be reconstructed";
+ case EPvUsageMode::ReconstructMultiple: return "multiple primary vertices will be reconstructed";
+ case EPvUsageMode::Target: return "target center is used for primary vertex";
+ case EPvUsageMode::Mc: return "MC-true primary vertex is used";
+ default: return "";
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string V0FinderTask::ToString(const FairTrackParam* pParam)
+{
+ std::stringstream msg;
+ msg << "Track parameters: r=(" << pParam->GetX() << ", " << pParam->GetY() << ", " << pParam->GetZ()
+ << "), tx=" << pParam->GetTx() << ", ty=" << pParam->GetTy() << ", q/p=" << pParam->GetQp()
+ << ", Covariance matrix:";
+ for (int i = 0; i < 5; ++i) {
+ msg << '\n';
+ for (int j = 0; j <= i; ++j) {
+ msg << std::setw(15) << pParam->GetCovariance(i, j) << ' ';
+ }
+ }
+ return msg.str();
+}
diff --git a/reco/KF/CbmKFV0FinderTask.h b/reco/KF/CbmKFV0FinderTask.h
new file mode 100644
index 0000000000..567c0e721f
--- /dev/null
+++ b/reco/KF/CbmKFV0FinderTask.h
@@ -0,0 +1,461 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmV0FinderTask.h
+/// \brief A FairTask for V0 candidates finding in mCBM
+/// \since 10.01.2025
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+// TODO: Provide a configurator to the class
+
+#pragma once
+
+#include "CbmEnumArray.h"
+#include "CbmKFParticleFinder.h" // for the CbmKFParticleFinder::InverseChi2Prob()
+#include "CbmKFV0FinderQa.h"
+#include "CbmKFVertex.h"
+#include "CbmVertex.h"
+#include "FairTask.h"
+#include "KFParticleTopoReconstructor.h"
+
+#include <memory>
+
+class CbmEvent;
+class CbmGlobalTrack;
+class CbmStsTrack;
+class CbmTofHit;
+class TClonesArray;
+class FairTrackParam;
+class KFPTrackVector;
+class KFVertex;
+class KFParticleFinder;
+
+namespace cbm::kfp
+{
+ /// \class V0FinderTask
+ /// \brief A class to find V0 candidates in mCBM
+ class V0FinderTask : public FairTask {
+ public:
+ /// \struct QpAndBeta
+ /// \brief Qp and beta container
+ struct QpAndBeta {
+ double fQp{std::numeric_limits<double>::signaling_NaN()};
+ double fQpVar{std::numeric_limits<double>::signaling_NaN()};
+ double fBeta{std::numeric_limits<double>::signaling_NaN()};
+ double fBetaVar{std::numeric_limits<double>::signaling_NaN()};
+ };
+
+ /// \struct DcaVector
+ /// \brief A vector representation of DCA to target
+ struct DcaVector {
+ double fAbs{std::numeric_limits<double>::signaling_NaN()}; ///< Absolute value
+ double fX{std::numeric_limits<double>::signaling_NaN()}; ///< X-component of the unit-vector
+ double fY{std::numeric_limits<double>::signaling_NaN()}; ///< Y-component of the unit-vector
+ };
+
+ /// \enum Counter
+ /// \brief Enumeration of counters (for the entire run)
+ enum class ECounter : uint8_t
+ {
+ TracksTotal, //< Total number of tracks
+ TracksSelected, //< Tracks, which satisfy topology PID applicability
+ TracksInfiniteParam, //< Tracks, which have infinite parameters
+ TracksWoTofHits, //< Tracks, which have no TOF hits
+ TracksWNegativeTofHitTime, //< Tracks, the last TOF hit of which has a negative time (it's time is less then the t0)
+ TracksWoStsHits, //< Tracks, which have no STS hits
+ TracksWoPid, //< Tracks, which has undefined PID
+ TracksWoMomentum, //< Tracks, which has undefined PID
+ TracksWAtLeastOneTofHit, //< Tracks with at least one tof hit
+ TracksWAtLeastTwoTofHits, //< Tracks with at least two tof hits
+ TracksWithUnphysicalBeta, //< Tracks with beta > 1
+ TracksW2Sts2Tof, //< Tracks with at least two STS hits and two TOF hits
+ Pions, //< Number of pion-candidates
+ Protons, //< Number of proton-candidates
+ EventsTotal, //< Total number of events
+ EventsLambdaCand, //< Events with at least one pion and one proton candidate
+ KfpEventsLambdaCand, //< Events with lambda-candidates in KF-particle
+ KfpLambdaCandidates, //< Number of lambda-candidates
+ END
+ };
+
+ /// \enum EPidApproach
+ /// \brief PID approach used
+ enum class EPidApproach : uint8_t
+ {
+ Topo, //< Topological PID (DCA calculation for a track)
+ Mc //< MC-truth PID
+ };
+
+ /// \enum EProcessingMode
+ /// \brief Data processing mode
+ enum class EProcessingMode : uint8_t
+ {
+ TimeBased, //< Time based simulation or real data
+ EventBased //< Event based simulation
+ };
+
+ /// \enum EPvUsageMode
+ /// \brief Primary vertex finding/handling mode
+ enum class EPvUsageMode : uint8_t
+ {
+ Target, //< Use target center as a primary vertex
+ Reconstructed, //< Use reconstructed primary vertex
+ ReconstructSingle, //< Reconstruct a single PV from given track vectors
+ ReconstructMultiple, //< Reconstruct multiple PVs (in case of pile-up)
+ Mc //< Use MC-true primary vertex
+ };
+
+ /// \brief Constructor
+ /// \param verbose Verbosity of the task
+ explicit V0FinderTask(int verbose = 1) : FairTask("V0FinderTask", verbose){};
+
+ /// \brief Destructor
+ ~V0FinderTask() = default;
+
+ // Eliminate copy and move
+ V0FinderTask(const V0FinderTask&) = delete;
+ V0FinderTask(V0FinderTask&&) = delete;
+ V0FinderTask& operator=(const V0FinderTask&) = delete;
+ V0FinderTask& operator=(V0FinderTask&&) = delete;
+
+
+ //* FairTask overriden methods
+
+ /// \brief Executes the task
+ void Exec(Option_t*) override;
+
+ /// \brief Action on the end of the run
+ void Finish() override;
+
+ /// \brief Initializes the task
+ InitStatus Init() override;
+
+
+ //* Accessors
+
+ /// \brief Mutable access to the KfParticleFinder of the run topology reconstructor
+ KFParticleFinder* GetKFParticleFinder() { return fpTopoReconstructorRun->GetKFParticleFinder(); }
+
+ /// \brief Constant access to the KfParticleFinder of the run topology reconstructor
+ const KFParticleFinder* GetKFParticleFinder() const { return fpTopoReconstructorRun->GetKFParticleFinder(); }
+
+ /// \brief Gets DCA to origin
+ const auto& GetTrackDcaToOrigin() const { return fvTrackDca; }
+
+
+ //* Task logic and preliminary selection cut setters
+
+ /// \brief Sets processing mode (time-based/event-based)
+ void SetProcessingMode(EProcessingMode mode) { fProcessingMode = mode; }
+
+ /// \brief Sets PID approach
+ void SetPidApproach(EPidApproach pid) { fPidApproach = pid; }
+
+ /// \brief Sets PV finding mode
+ void SetPvFindingMode(EPvUsageMode mode) { fPvUsageMode = mode; }
+
+ /// \brief Sets the MC flag (if MC information required)
+ void SetUseMc(bool bUseMc) { fbUseMc = bUseMc; }
+
+ /// \brief Sets the flag: if the QA should be executed
+ void SetRunQa(bool bRunQa) { fbRunQa = bRunQa; }
+
+ /// \brief Sets minimal pion DCA to primary vertex
+ /// \param dca DCA [cm]
+ void SetMinPionDca(double dca) { fMinPionDca = dca; }
+
+ /// \brief Sets minimal proton DCA to primary vertex
+ /// \param dca DCA [cm]
+ void SetMinProtonDca(double dca) { fMinProtonDca = dca; }
+
+ /// \brief Sets a file name for the QA
+ /// \param fileName Name of file
+ void SetQaOutputFileName(const TString& fileName) { fsQaOutputName = fileName; }
+
+ /// \brief Assignes an uncertainty to the momentum measurement
+ /// \param uncertainty Relative uncertainty ( = sqrt(var(q/p)) / (q/p))
+ void SetQpAssignedUncertainty(double uncertainty) { fQpAssignedUncertainty = uncertainty; }
+
+ /// \brief Sets an offset to t0
+ /// \param offset An offset [ns]
+ void SetTzeroOffset(double offset) { fTzeroOffset = offset; }
+
+ /// \brief Special settings in the mixed-event analysis mode
+ /// \param bMixedEvent Flag: true - run in the mixed-event mode
+ // NOTE: The KFParticleFinder has its internal mixed-event analysis mode, but was not used in mCBM Lambda-analysis
+ // the internal mixed-event mode can be selected using the following method:
+ // fpTopoReconstructorRun->SetMixedEventAnalysis()
+ void SetMixedEventMode(bool bMixedEvent) { fbMixedEventMode = bMixedEvent; }
+
+
+ //* KFParticleFinder setters
+
+ /// \brief Adds particle to reconstruction list
+ /// \param pdg A PDG code of the particle to be reconstructed
+ void AddDecayToReconstructionList(int pdg) { GetKFParticleFinder()->AddDecayToReconstructionList(pdg); }
+
+ /// \brief Sets P-value cut to the topology reconstructor
+ /// \param pVal P-value cut
+ void SetPrimaryProbCut(double pVal)
+ {
+ constexpr int Ndf{2};
+ fpTopoReconstructorRun->SetChi2PrimaryCut(CbmKFParticleFinder::InversedChi2Prob(pVal, Ndf));
+ }
+
+ /// \brief Sets cut on the distance between secondary tracks at the DCA point
+ /// \param cut Cut value [cm]
+ void SetMaxDistanceBetweenParticlesCut(float cut) { GetKFParticleFinder()->SetMaxDistanceBetweenParticlesCut(cut); }
+
+ /// \brief Sets cut on the distance to the primary vertex from the decay vertex
+ /// \param cut Cut value [cm]
+ void SetLCut(float cut) { GetKFParticleFinder()->SetLCut(cut); }
+
+ /// \brief Sets cut on \f$\chi^2_{prim}\f$ of each track for 2-daughter decays
+ /// \param cut Cut value
+ void SetChiPrimaryCut2D(float cut) { GetKFParticleFinder()->SetChiPrimaryCut2D(cut); }
+
+ /// \brief Sets cut on \f$\chi^2_{geo}\f$ for 2-daughter decays
+ /// \param cut Cut value
+ void SetChi2Cut2D(float cut) { GetKFParticleFinder()->SetChi2Cut2D(cut); }
+
+ /// \brief Sets cut on \f$l/\Delta l\f$ for 2-daughter decays
+ /// \param cut Cut value
+ void SetLdLCut2D(float cut) { GetKFParticleFinder()->SetLdLCut2D(cut); }
+
+ /// \brief Sets cuts on selection of secondary and primary candidates
+ /// \param sigmaMass \f$\sigma_{M}\f$ [GeV/c2]
+ /// \param chi2Topo \f$\chi^2_{topo}\f$
+ /// \param ldl \f$l/\Delta l\f$
+ void SetSecondaryCuts(float sigmaMass, float chi2Topo, float ldl)
+ {
+ GetKFParticleFinder()->SetSecondaryCuts(sigmaMass, chi2Topo, ldl);
+ }
+
+ /// \brief Sets \f$l/\Delta l\f$ cut for \f$\Xi\f$ and \f$\Omega\f$
+ /// \param cut Cut value
+ void SetLdLCutXiOmega(float cut) { GetKFParticleFinder()->SetLdLCutXiOmega(cut); }
+
+ /// \brief Sets \f$\chi^2_{topo}\f$ cut for \f$\Xi\f$ and \f$\Omega\f$
+ /// \param cut Cut value
+ void SetChi2TopoCutXiOmega(float cut) { GetKFParticleFinder()->SetChi2TopoCutXiOmega(cut); }
+
+ /// \brief Sets \f$\chi^2_{geo}\f$ cut for \f$\Xi\f$ and \f$\Omega\f$
+ /// \param cut Cut value
+ void SetChi2CutXiOmega(float cut) { GetKFParticleFinder()->SetChi2CutXiOmega(cut); }
+
+ /// \brief Sets \f$\chi^2_{topo}\f$ cut for resonances
+ /// \param cut Cut value
+ void SetChi2TopoCutResonances(float cut) { GetKFParticleFinder()->SetChi2TopoCutResonances(cut); }
+
+ /// \brief Sets \f$\chi^2_{geo}\f$ cut for resonances
+ /// \param cut Cut value
+ void SetChi2CutResonances(float cut) { GetKFParticleFinder()->SetChi2CutResonances(cut); }
+
+ /// \brief Sets the cut on transverse momentum of each daughter track of low mass vector mesons
+ /// \param cut Cut value [GeV/c]
+ void SetPtCutLMVM(float cut) { GetKFParticleFinder()->SetPtCutLMVM(cut); }
+
+ /// \brief Sets the cut on momentum of each daughter track of low mass vector mesons in dimuon channel
+ /// \param cut Cut value [GeV/c]
+ void SetPCutLMVM(float cut) { GetKFParticleFinder()->SetPCutLMVM(cut); }
+
+ /// \brief Sets the cut on transverse momentum of each daughter track of \f$J/\psi\f$
+ /// \param cut Cut value [GeV/c]
+ void SetPtCutJPsi(float cut) { GetKFParticleFinder()->SetPtCutJPsi(cut); }
+
+ /// \brief Sets the cut on transverse momentum of each daughter track of open charm particles
+ /// \param cut Cut value [GeV/c]
+ void SetPtCutCharm(float cut) { GetKFParticleFinder()->SetPtCutCharm(cut); }
+
+ /// \brief Sets cut on \f$\chi^2_{prim}\f$ of each track for open charm particles
+ /// \param cut Cut value [GeV/c]
+ void SetChiPrimaryCutCharm(float cut) { GetKFParticleFinder()->SetChiPrimaryCutCharm(cut); }
+
+ /// \brief Sets \f$l/\Delta l\f$ cut for open charm with >=3 daughters
+ /// \param cut Cut value
+ void SetLdLCutCharmManybodyDecays(float cut) { GetKFParticleFinder()->SetLdLCutCharmManybodyDecays(cut); }
+
+ /// \brief Sets \f$\chi^2_{topo}\f$ cut for open charm with >=3 daughters
+ /// \param cut Cut value
+ void SetChi2TopoCutCharmManybodyDecays(float cut) { GetKFParticleFinder()->SetChi2TopoCutCharmManybodyDecays(cut); }
+
+ /// \brief Sets \f$\chi^2_{geo}\f$ cut for open charm with >=3 daughters
+ /// \param cut Cut value
+ void SetChi2CutCharmManybodyDecays(float cut) { GetKFParticleFinder()->SetChi2CutCharmManybodyDecays(cut); }
+
+ /// \brief Sets \f$l/\Delta l\f$ cut for open charm with 2 daughters
+ /// \param cut Cut value
+ void SetLdLCutCharm2D(float cut) { GetKFParticleFinder()->SetLdLCutCharm2D(cut); }
+
+ /// \brief Sets \f$\chi^2_{topo}\f$ cut for open charm with 2 daughters
+ /// \param cut Cut value
+ void SetChi2TopoCutCharm2D(float cut) { GetKFParticleFinder()->SetChi2TopoCutCharm2D(cut); }
+
+ /// \brief Sets \f$\chi^2_{geo}\f$ cut for open charm with 2 daughters
+ /// \param cut Cut value
+ void SetChi2CutCharm2D(float cut) { GetKFParticleFinder()->SetChi2CutCharm2D(cut); }
+
+
+ //* Static helper methods
+
+ /// \brief String representation of processing mode
+ /// \param mode Data processing mode
+ static std::string ToString(EProcessingMode mode);
+
+ /// \brief String representation of PID approach
+ /// \param pid PID approach
+ static std::string ToString(EPidApproach pid);
+
+ /// \brief String representation of the PV finding mode
+ /// \param pvMode PV finding mode
+ static std::string ToString(EPvUsageMode pvMode);
+
+ /// \brief String representation of
+ /// \param pParam Parameters of the track
+ static std::string ToString(const FairTrackParam* pParam);
+
+ /// \brief Checks track parameter validity
+ /// \param pParam Pointer to a track parameter
+ /// \return true The parameters are ok
+ /// \return false The parameters are invalid
+ static bool CheckTrackParam(const FairTrackParam* pParam);
+
+ private:
+ //* Auxilary internal methods
+
+ /// \brief Assigns momentum to a global track
+ /// \param pTrack Pointer to a global track
+ /// \param pdg PID hypothesis
+ /// \return true Momentum is set
+ /// \return false Momentum is not set (the track must be rejected)
+ bool AssignMomentum(CbmGlobalTrack* pTrack, int pdg);
+
+ /// \brief Estimates distance to the origin
+ /// \param pStsTrack STS track
+ /// \return Vector of the DCA {abs, ex, ey}
+ DcaVector EstimateDcaToOrigin(const CbmStsTrack* pTrack) const;
+
+ /// \brief Estimates q/p of the track, using one TOF hit (Norbert's method)
+ /// \param pTofHit Pointer to a TOF hit
+ /// \param pdg PID hypothesis
+ /// \note Gives smaller statistical uncertainty, but has a systematic effect (to be studied with MC)
+ QpAndBeta EstimateQp(const CbmTofHit* pTofHit, int pdg) const;
+
+ /// \brief Makes a KF-particle track vector
+ /// \param vpTracks Vector of pointers to global tracks
+ /// \param vTrackIds Vector of track indices
+ /// \param vChi2ToPv Vector of chi2 of assigning a track as primary
+ /// \param bAtFirstPoint true: the first track point is used, false: the last track point is used
+ /// \note No track selection is performed
+ KFPTrackVector MakeKfpTrackVector(const std::vector<const CbmGlobalTrack*>& vpTracks,
+ const std::vector<int>& vTrackIds, const std::vector<float>& vChi2ToPv,
+ bool bAtFirstPoint) const;
+
+ /// \brief Makes a KF vertex
+ /// \param x x-coordinate of PV [cm]
+ /// \param y y-coordinate of PV [cm]
+ /// \param z z-coordinate of PV [cm]
+ KFVertex MakeKfpPrimaryVertex(float x, float y, float z) const;
+
+ /// \brief Infers PID hypothesis for a track using track topology
+ /// \param dca An absolute value of DCA to the PV estimation
+ /// \return PDG code of the particle
+ /// \return -2 If the PID analysis cannot be performed for a track of a given quality
+ int InferTrackPidTopo(double dca) const;
+
+ /// \brief Processes one event
+ /// \tparam UseEvent Explicit flag of using the CbmEvent pointer to access data
+ /// \param pEvent Pointer to event
+ /// \return true Event is accepted for V0 reconstruction
+ /// \return false Event is rejected for V0 reconstruction
+ /// \note In the event-based mode, the ProcessEvent(nullptr) must be called explicitly
+ template<bool UseEvent>
+ bool ProcessEvent(const CbmEvent* pEvent);
+
+ /// \brief Selects/rejects a track
+ /// \param pTrack A pointer to global track
+ /// \param iTrk Global index of the track
+ /// \return true: Track is selected
+ /// \return false: Track is rejected
+ /// \note Modifies the track, assigning momentum and PID hypothesis
+ ///
+ /// In this function the preliminary selection of tracks is performed. The selection criteria:
+ /// 1) the track has at least two STS hits for estimation of the PDG codes (pi-, proton, or primary (kaon PDG))
+ /// 2) the track has at least one TOF hit for estimation of its momentum
+ /// 2.1) tracks with the TOF hit times (t - t0) < 0 are rejected
+ /// 2.2) tracks with the beta > 1 are rejected
+ /// 3) the track has defined parameters in the first hit
+ /// NOTE: this cut was used, because the tracks with negative TOF times were not rejected. Now it seems,
+ /// that the track parameters are always defined.
+ // TODO: Replace with a functor (a function from a separate FairTask)
+ bool SelectTrack(CbmGlobalTrack* pTrack, int iTrk);
+
+ /// \brief Shifts TOF hits to the t0 estimation
+ /// \param pEvent Pointer to the event
+ /// \return value of estimated t0 [ns]
+ double ShiftTofHitsToTzero(const CbmEvent* pEvent);
+
+ /// \brief Stores particles, reconstructed in event to the run topology reconstructor
+ void StoreParticles();
+
+
+ //* Task constants (To be moved to a configuration)
+ static constexpr int kPrimaryPdg{321}; ///< PID hypothesis of primary tracks (kaons?)
+ static constexpr int kUndefPdg{-2}; ///< Undefined value, such tracks will be skipped
+ static constexpr float kChi2PvPrimThrsh{3.}; ///< Chi2 threshold of assigning tracks as primaries in PV reco
+
+ //* Physical constants
+ static constexpr double kPionMass{0.13957039}; ///< Pion mass [GeV/c2]
+ static constexpr double kProtonMass{0.938272088}; ///< Proton mass [GeV/c2]
+ static constexpr double kSpeedOfLight{29.9792458}; ///< Speed of light [cm/ns]
+
+ //* Different run-time flags;
+ double fTzeroOffset{0.}; ///< Offset for T0
+ double fMinPionDca{1.5}; ///< Minimum DCA to PV for pions
+ double fMinProtonDca{0.5}; ///< Minimum DCA to PV for protons
+ double fQpAssignedUncertainty{0.1}; ///< Assigned relative uncertainty for q/p estimation
+ int fPrimaryAssignedPdg{321}; ///< Assigned PDG hypothesis for primary particles
+
+ //* Input data
+ TClonesArray* fpBrRecoEvents{nullptr};
+ TClonesArray* fpBrGlobalTracks{nullptr};
+ TClonesArray* fpBrStsTracks{nullptr};
+ TClonesArray* fpBrTrdTracks{nullptr};
+ TClonesArray* fpBrTofTracks{nullptr};
+ TClonesArray* fpBrStsHits{nullptr};
+ TClonesArray* fpBrTrdHits{nullptr};
+ TClonesArray* fpBrTofHits{nullptr};
+ CbmVertex* fpBrPrimaryVertex{nullptr};
+
+ //* Temporary data and auxilary variables (per FairTask::Exec() call)
+ std::vector<DcaVector> fvTrackDca; ///< Track DCA vector [n selected tracks]
+ std::unique_ptr<CbmVertex> fpOrigin{
+ std::make_unique<CbmVertex>()}; ///< Origin (e.g., can be either reconstructed PV or target)
+
+ //* KFParticleFinder utilities
+ /// \brief Main topology reconstructor
+ // Main topology reconstructor
+ std::unique_ptr<KFParticleTopoReconstructor> fpTopoReconstructorRun{
+ std::make_unique<KFParticleTopoReconstructor>()};
+ std::unique_ptr<KFParticleTopoReconstructor> fpTopoReconstructorEvent{
+ std::make_unique<KFParticleTopoReconstructor>()};
+
+ //* Control flags
+ EProcessingMode fProcessingMode{EProcessingMode::TimeBased}; ///< Processing mode
+ EPidApproach fPidApproach{EPidApproach::Topo}; ///< PID approach used
+ EPvUsageMode fPvUsageMode{EPvUsageMode::Target}; ///< Primary vertex mode
+ bool fbMixedEventMode{false}; ///< Run in a mix-event mode
+ bool fbUseMc{false}; ///< Run using MC-information
+ bool fbRunQa{false}; ///< Run QA
+ bool fbUsePvChi2Selection{false}; ///< Select
+
+ //* QA and monitoring variables
+ cbm::core::EnumArray<ECounter, size_t> fCounters{{0}}; ///< Counters per run
+ std::unique_ptr<V0FinderQa> fpQa{nullptr}; ///< If QA is processed
+ TString fsQaOutputName{"./V0FinderQa.root"}; ///< Output QA name
+
+
+ ClassDefOverride(V0FinderTask, 0);
+ };
+} // namespace cbm::kfp
diff --git a/reco/KF/KFParticleInterface.cmake b/reco/KF/KFParticleInterface.cmake
index d5cd174857..7f594ec0a9 100644
--- a/reco/KF/KFParticleInterface.cmake
+++ b/reco/KF/KFParticleInterface.cmake
@@ -13,6 +13,8 @@ CbmKFParticleFinder.cxx
CbmKFParticleFinderPID.cxx
CbmKFParticleFinderQa.cxx
CbmKFParticleInterface.cxx
+CbmKFV0FinderTask.cxx
+CbmKFV0FinderQa.cxx
)
If(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
diff --git a/reco/KF/KFParticleInterfaceLinkDef.h b/reco/KF/KFParticleInterfaceLinkDef.h
index dedd55f91f..f04b91c6bf 100644
--- a/reco/KF/KFParticleInterfaceLinkDef.h
+++ b/reco/KF/KFParticleInterfaceLinkDef.h
@@ -12,5 +12,7 @@
#pragma link C++ class CbmKFParticleFinderPID + ;
#pragma link C++ class CbmKFParticleFinderQa + ;
#pragma link C++ class CbmKFParticleInterface + ;
+#pragma link C++ class cbm::kfp::V0FinderTask + ;
+#pragma link C++ class cbm::kfp::V0FinderQa + ;
#endif
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index 0655131131..0a41dd11d3 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -145,6 +145,7 @@ try {
InitHitBranch(ca::EDetectorID::kMuch, "MuchPixelHit");
InitHitBranch(ca::EDetectorID::kTrd, "TrdHit");
if (!InitHitBranch(ca::EDetectorID::kTof, "TofCalHit")) {
+ LOG(warn) << "TimeSliceReader: TofCalHit branch not found, trying to get TofHit branch of uncalibrated hits";
InitHitBranch(ca::EDetectorID::kTof, "TofHit");
}
diff --git a/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx b/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx
index 38b304fa67..b84a2d8a3e 100644
--- a/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx
+++ b/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx
@@ -732,8 +732,9 @@ bool CbmAlgoBuildRawEvents::SetBmonEventTime(CbmEvent* event)
eventTime = std::min(pDigi->GetTime(), eventTime);
}
- if (timeSet)
+ if (timeSet) {
event->SetStartTime(eventTime);
+ }
else
LOG(warning) << "CbmAlgoBuildRawEvents::SetBmonEventTime : failed";
diff --git a/reco/tasks/CbmTaskInspectDigiTimeslice.cxx b/reco/tasks/CbmTaskInspectDigiTimeslice.cxx
index 169d3cbf0e..bc8d30d69a 100644
--- a/reco/tasks/CbmTaskInspectDigiTimeslice.cxx
+++ b/reco/tasks/CbmTaskInspectDigiTimeslice.cxx
@@ -4,6 +4,8 @@
#include "CbmTaskInspectDigiTimeslice.h"
+#include "CbmStsAddress.h"
+
#include <FairRootManager.h>
#include <Logger.h>
--
GitLab