From eef74740e956e3c804d852a8650a4e4253b7e4d5 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 25 Sep 2025 13:56:45 +0200 Subject: [PATCH 1/9] uploading nucleiQC, moving nucleiUtils.h and updating dependancies --- PWGLF/DataModel/LFSlimNucleiTables.h | 15 + PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx | 8 +- PWGLF/TableProducer/Nuspex/nucleiUtils.h | 187 ------ PWGLF/TableProducer/QC/CMakeLists.txt | 5 + PWGLF/TableProducer/QC/nucleiQC.cxx | 427 ++++++++++++++ PWGLF/Utils/nucleiUtils.h | 539 ++++++++++++++++++ 6 files changed, 990 insertions(+), 191 deletions(-) delete mode 100644 PWGLF/TableProducer/Nuspex/nucleiUtils.h create mode 100644 PWGLF/TableProducer/QC/nucleiQC.cxx create mode 100644 PWGLF/Utils/nucleiUtils.h diff --git a/PWGLF/DataModel/LFSlimNucleiTables.h b/PWGLF/DataModel/LFSlimNucleiTables.h index e1c0e0eb364..5f14cf5ab46 100644 --- a/PWGLF/DataModel/LFSlimNucleiTables.h +++ b/PWGLF/DataModel/LFSlimNucleiTables.h @@ -185,6 +185,21 @@ DECLARE_SOA_TABLE(NucleiPairTable, "AOD", "NUCLEIPAIRTABLE", NucleiPairTableNS::ClusterSizesITS2, NucleiPairTableNS::Flags2); +// Reduced table +DECLARE_SOA_TABLE(NucleiTableRed, "AOD", "NUCLEITABLERED", + NucleiTableNS::Pt, + NucleiTableNS::Eta, + NucleiTableNS::Phi, + NucleiTableNS::TPCInnerParam, + NucleiTableNS::ITSclusterSizes, + NucleiTableNS::TPCsignal, + NucleiTableNS::Beta, + NucleiTableNS::DCAxy, + NucleiTableNS::DCAz, + NucleiTableNS::Flags, + NucleiTableNS::PDGcode, + NucleiTableNS::MotherPDGcode); + } // namespace o2::aod #endif // PWGLF_DATAMODEL_LFSLIMNUCLEITABLES_H_ diff --git a/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx b/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx index 14192166fae..c336bd8a74c 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx @@ -66,7 +66,7 @@ #include "TRandom3.h" -#include "nucleiUtils.h" +#include "PWGLF/Utils/nucleiUtils.h" using namespace o2; using namespace o2::framework; @@ -257,7 +257,7 @@ struct nucleiFlowTree { spectra.add("hTpcSignalDataSelected", "Specific energy loss for selected particles", HistType::kTH2F, {{600, -6., 6., "#it{p} (GeV/#it{c})"}, {1400, 0, 1400, "d#it{E} / d#it{X} (a. u.)"}}); spectra.add("hTofSignalData", "TOF beta", HistType::kTH2F, {{500, 0., 5., "#it{p} (GeV/#it{c})"}, {750, 0, 1.5, "TOF #beta"}}); - for (int iS{0}; iS < nuclei::species; ++iS) { + for (int iS{0}; iS < nuclei::Species::kNspecies; ++iS) { for (int iMax{0}; iMax < 2; ++iMax) { nuclei::pidCutTPC[iS][iMax] = cfgNsigmaTPC->get(iS, iMax); // changed pidCut to pidCutTPC so that it compiles TODO: check if it is correct } @@ -328,7 +328,7 @@ struct nucleiFlowTree { bool selectedTPC[5]{false}, goodToAnalyse{false}; std::array nSigmaTPC; - for (int iS{0}; iS < nuclei::species; ++iS) { + for (int iS{0}; iS < nuclei::Species::kNspecies; ++iS) { double expBethe{tpc::BetheBlochAleph(static_cast(correctedTpcInnerParam * bgScalings[iS][iC]), cfgBetheBlochParams->get(iS, 0u), cfgBetheBlochParams->get(iS, 1u), cfgBetheBlochParams->get(iS, 2u), cfgBetheBlochParams->get(iS, 3u), cfgBetheBlochParams->get(iS, 4u))}; @@ -371,7 +371,7 @@ struct nucleiFlowTree { if (!collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { flag |= kITSrof; } - for (int iS{0}; iS < nuclei::species; ++iS) { + for (int iS{0}; iS < nuclei::Species::kNspecies; ++iS) { bool selectedTOF{false}; if (std::abs(dcaInfo[1]) > cfgDCAcut->get(iS, 1)) { continue; diff --git a/PWGLF/TableProducer/Nuspex/nucleiUtils.h b/PWGLF/TableProducer/Nuspex/nucleiUtils.h deleted file mode 100644 index 29cd489768d..00000000000 --- a/PWGLF/TableProducer/Nuspex/nucleiUtils.h +++ /dev/null @@ -1,187 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -#ifndef PWGLF_TABLEPRODUCER_NUSPEX_NUCLEIUTILS_H_ -#define PWGLF_TABLEPRODUCER_NUSPEX_NUCLEIUTILS_H_ - -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace o2::constants::physics; - -struct NucleusCandidate { - int globalIndex; - int collTrackIndex; - float pt; - float eta; - float phi; - float tpcInnerParam; - float beta; - float zVertex; - int nContrib; - float DCAxy; - float DCAz; - float TPCsignal; - float ITSchi2; - float TPCchi2; - float TOFchi2; - std::array nSigmaTPC; - std::array tofMasses; - bool fillTree; - bool fillDCAHist; - bool correctPV; - bool isSecondary; - bool fromWeakDecay; - uint16_t flags; - uint8_t TPCfindableCls; - uint8_t TPCcrossedRows; - uint8_t ITSclsMap; - uint8_t TPCnCls; - uint8_t TPCnClsShared; - uint8_t ITSnCls; - uint32_t clusterSizesITS; -}; - -struct NucleusCandidateFlow { - float centFV0A; - float centFT0M; - float centFT0A; - float centFT0C; - float psiFT0A; - float psiFT0C; - float psiTPC; - float psiTPCl; - float psiTPCr; - float qFT0A; - float qFT0C; - float qTPC; - float qTPCl; - float qTPCr; -}; - -namespace nuclei -{ -constexpr double bbMomScalingDefault[5][2]{ - {1., 1.}, - {1., 1.}, - {1., 1.}, - {1., 1.}, - {1., 1.}}; -constexpr double betheBlochDefault[5][6]{ - {-136.71, 0.441, 0.2269, 1.347, 0.8035, 0.09}, - {-136.71, 0.441, 0.2269, 1.347, 0.8035, 0.09}, - {-239.99, 1.155, 1.099, 1.137, 1.006, 0.09}, - {-321.34, 0.6539, 1.591, 0.8225, 2.363, 0.09}, - {-586.66, 1.859, 4.435, 0.282, 3.201, 0.09}}; -constexpr double nSigmaTPCdefault[5][2]{ - {-5., 5.}, - {-5., 5.}, - {-5., 5.}, - {-5., 5.}, - {-5., 5.}}; -// constexpr double nSigmaTOFdefault[5][2]{ -// {-5., 5.}, -// {-5., 5.}, -// {-5., 5.}, -// {-5., 5.}, -// {-5., 5.}}; -constexpr double DCAcutDefault[5][2]{ - {1., 1.}, - {1., 1.}, - {1., 1.}, - {1., 1.}, - {1., 1.}}; -constexpr int TreeConfigDefault[5][2]{ - {0, 0}, - {0, 0}, - {0, 0}, - {0, 0}, - {0, 0}}; -constexpr int FlowHistDefault[5][1]{ - {0}, - {0}, - {0}, - {0}, - {0}}; -constexpr int DCAHistDefault[5][2]{ - {0, 0}, - {0, 0}, - {0, 0}, - {0, 0}, - {0, 0}}; -constexpr double DownscalingDefault[5][1]{ - {1.}, - {1.}, - {1.}, - {1.}, - {1.}}; -// constexpr bool storeTreesDefault[5]{false, false, false, false, false}; -constexpr int species{5}; -constexpr float charges[5]{1.f, 1.f, 1.f, 2.f, 2.f}; -constexpr float masses[5]{MassProton, MassDeuteron, MassTriton, MassHelium3, MassAlpha}; -static const std::vector matter{"M", "A"}; -static const std::vector pidName{"TPC", "TOF"}; -static const std::vector names{"proton", "deuteron", "triton", "He3", "alpha"}; -static const std::vector treeConfigNames{"Filter trees", "Use TOF selection"}; -static const std::vector flowConfigNames{"Save flow hists"}; -static const std::vector DCAConfigNames{"Save DCA hist", "Matter/Antimatter"}; -static const std::vector nSigmaConfigName{"nsigma_min", "nsigma_max"}; -static const std::vector nDCAConfigName{"max DCAxy", "max DCAz"}; -static const std::vector DownscalingConfigName{"Fraction of kept candidates"}; -static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; -static const std::vector chargeLabelNames{"Positive", "Negative"}; - -float pidCutTPC[5][2]; //[species][lower/upper limit] - -o2::base::MatLayerCylSet* lut = nullptr; - -std::vector candidates; -std::vector candidates_flow; - -enum centDetectors { - kFV0A = 0, - kFT0M = 1, - kFT0A = 2, - kFT0C = 3 -}; - -static const std::vector centDetectorNames{"FV0A", "FT0M", "FT0A", "FT0C"}; - -enum evSel { - kTVX = 0, - kZvtx, - kTFborder, - kITSROFborder, - kNoSameBunchPileup, - kIsGoodZvtxFT0vsPV, - kIsGoodITSLayersAll, - kIsEPtriggered, - kNevSels -}; - -static const std::vector eventSelectionTitle{"Event selections"}; -static const std::vector eventSelectionLabels{"TVX", "Z vtx", "TF border", "ITS ROF border", "No same-bunch pile-up", "kIsGoodZvtxFT0vsPV", "isGoodITSLayersAll", "isEPtriggered"}; - -constexpr int EvSelDefault[8][1]{ - {1}, - {1}, - {0}, - {0}, - {0}, - {0}, - {0}, - {0}}; -} // namespace nuclei - -#endif // PWGLF_TABLEPRODUCER_NUSPEX_NUCLEIUTILS_H_ diff --git a/PWGLF/TableProducer/QC/CMakeLists.txt b/PWGLF/TableProducer/QC/CMakeLists.txt index db3fdb4876b..79501158034 100644 --- a/PWGLF/TableProducer/QC/CMakeLists.txt +++ b/PWGLF/TableProducer/QC/CMakeLists.txt @@ -18,3 +18,8 @@ o2physics_add_dpl_workflow(flow-qc SOURCES flowQC.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(nucleiqc + SOURCES nucleiQC.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/QC/nucleiQC.cxx b/PWGLF/TableProducer/QC/nucleiQC.cxx new file mode 100644 index 00000000000..f0cfbdf52f8 --- /dev/null +++ b/PWGLF/TableProducer/QC/nucleiQC.cxx @@ -0,0 +1,427 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// Nuclei spectra analysis task +// ======================== +// +// Executable + dependencies: +// +// Data (run3): +// o2-analysis-lf-nuclei-spectra, o2-analysis-timestamp +// o2-analysis-pid-tof-base, o2-analysis-multiplicity-table, o2-analysis-event-selection +// (to add flow: o2-analysis-qvector-table, o2-analysis-centrality-table) + +#include +#include +#include +#include +#include +#include + +#include "Math/Vector4D.h" + +#include "CCDB/BasicCCDBManager.h" + +#include "Common/Core/TrackSelection.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/Core/PID/PIDTOF.h" +#include "Common/TableProducer/PID/pidTOFBase.h" +#include "Common/Core/EventPlaneHelper.h" +#include "Common/DataModel/Qvectors.h" +#include "Common/Tools/TrackTuner.h" +#include "Common/Core/RecoDecay.h" + +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsTPC/BetheBlochAleph.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" + +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" + +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" + +#include "ReconstructionDataFormats/Track.h" + +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "PWGLF/DataModel/LFSlimNucleiTables.h" + +#include "TRandom3.h" + +#include "PWGLF/Utils/nucleiUtils.h" + +using namespace o2; +using namespace o2::framework; + + +struct nucleiQC { + +//using TrackCandidates = soa::Join; + using TrackCandidates = soa::Join; +//using TrackCandidatesMC = soa::Join; + using TrackCandidatesMC = soa::Join; + using Collision = soa::Join::iterator; + + Configurable cfgFillTable{"cfgFillTable", true, "Fill output tree"}; + Configurable> cfgSpeciesToProcess{"cfgSpeciesToProcess", {nuclei::speciesToProcessDefault[0], nuclei::Species::kNspecies, 1, nuclei::names, {"processNucleus"}}, "Nuclei to process"}; + Configurable> cfgEventSelections{"cfgEventSelections", {nuclei::EvSelDefault[0], 8, 1, nuclei::eventSelectionLabels, nuclei::eventSelectionTitle}, "Event selections"}; + Configurable cfgCentralityEstimator{"cfgCentralityEstimator", 0, "Centrality estimator (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3)"}; + Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {nuclei::betheBlochDefault[0], nuclei::Species::kNspecies, 6, nuclei::names, nuclei::betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; + Configurable> cfgUseCentralTpcCalibration{"cfgUseCentralTpcCalibration", {nuclei::useCentralTpcCalibrationDefault[0], nuclei::Species::kNspecies, 1, nuclei::names, {"UseCentralTpcCalibration"}}, "Use central TPC calibration"}; + + Configurable cfgRapidityMin{"cfgRapidityMin", -1., "Minimum rapidity value"}; + Configurable cfgRapidityMax{"cfgRapidityMax", 1., "Maximum rapidity value"}; + Configurable cfgRapidityCenterMass{"cfgRapidityCenterMass", 0.0f, "Center of mass rapidity"}; + + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable cfgCutTpcMom{"cfgCutTpcMom", 0.2f, "Minimum TPC momentum for tracks"}; + Configurable cfgCutNclusITS{"cfgCutNclusITS", 5, "Minimum number of ITS clusters"}; + Configurable cfgCutNclusTPC{"cfgCutNclusTPC", 70, "Minimum number of TPC clusters"}; + + Configurable> cfgNsigmaTPC{"cfgNsigmaTPC", {nuclei::nSigmaTPCdefault[0], nuclei::Species::kNspecies, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"}; + Configurable> cfgNsigmaTOF{"cfgNsigmaTOF", {nuclei::nSigmaTOFdefault[0], nuclei::Species::kNspecies, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"}; + + HistogramRegistry mQaHistograms { + "QA", + { + {"hEventSelections", "Event selections; Selection step; Counts", {HistType::kTH1D, {{nuclei::evSel::kNevSels + 1, -0.5f, static_cast(nuclei::evSel::kNevSels) + 0.5f}}}}, + {"hVtxZBefore", "Vertex distribution in Z before selections;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}}, + {"hVtxZ", "Vertex distribution in Z;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}}, + }, + OutputObjHandlingPolicy::AnalysisObject, + false, + true + }; + std::array(nuclei::Species::kNspecies)> mNucleiHistogramRegistries; + std::vector mSpeciesToProcess; + Produces mNucleiTableRed; + + std::vector mNucleiCandidates; + std::vector mFilledMcParticleIds; + + std::array(nuclei::Species::kNspecies)> mPidManagers; + + void init(o2::framework::InitContext&) { + + LOG(info) << "Check init"; + + for (int iSel = 0; iSel < nuclei::evSel::kNevSels; iSel++) { + mQaHistograms.get(HIST("hEventSelections"))->GetXaxis()->SetBinLabel(iSel + 1, nuclei::eventSelectionLabels[iSel].c_str()); + } + + LOG(info) << " check after hist ev sel"; + + for (int iSpecies = 0; iSpecies < static_cast(nuclei::Species::kNspecies); iSpecies++) { + if (cfgSpeciesToProcess->get(iSpecies) == 1) { + mSpeciesToProcess.emplace_back(iSpecies); + } + } + + LOG(info) << " check after species to be processed were saved"; + + for (auto iSpecies: mSpeciesToProcess) { + + float tpcBetheBlochParams[6]; + for (int iParam = 0; iParam < 6; iParam++) { + tpcBetheBlochParams[iParam] = cfgBetheBlochParams->get(iSpecies, iParam); + } + + LOG(info) << " check after bb params, ispecies: " << iSpecies; + + mNucleiHistogramRegistries[iSpecies] = nuclei::createHistogramRegistryNucleus(iSpecies); + + if (cfgUseCentralTpcCalibration->get(uint32_t(iSpecies), uint32_t(0)) == 0) { + mPidManagers[iSpecies] = nuclei::PidManager(nuclei::Species(iSpecies), tpcBetheBlochParams); + } else { + mPidManagers[iSpecies] = nuclei::PidManager(nuclei::Species(iSpecies)); + } + + LOG(info) << " check after pidManagers, ispecies: " << iSpecies; + } + + LOG(info) << " check end init"; + } + + template + bool trackSelection(const Ttrack& track) { + if (std::abs(track.eta()) > cfgCutEta || + track.tpcInnerParam() < cfgCutTpcMom || + track.itsNCls() < cfgCutNclusITS || + track.tpcNClsFound() < cfgCutNclusTPC || + track.tpcNClsCrossedRows() < 70 || + track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() || + track.tpcChi2NCl() > 4.f || + track.itsChi2NCl() > 36.f) { + return false; + } + return true; + } + + template + bool pidSelection(const int iSpecies, const Ttrack& track, const Tcollision& collision) + { + if (!nuclei::checkSpeciesValidity(iSpecies)) { + std::runtime_error("species contains invalid nucleus index"); + } + + const float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator); + const float nsigmaTPC = mPidManagers[iSpecies].getNSigmaTPC(track); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NsigmaTPC_preselectionVsCentrality"), track.pt(), nsigmaTPC, centrality); + if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(iSpecies, 1)) { + return false; + } + mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NsigmaTPCVsCentrality"), track.pt(), nsigmaTPC, centrality); + + const float nsigmaITS = mPidManagers[iSpecies].getNSigmaITS(track); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaITS_preselectionVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); + // add nsigmaITS cut ? + mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaITSVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); + + //const float nsigmaTOF = mPidManagers[iSpecies].getNSigmaTOF(track); + //mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaTOF_preselectionVsCentrality"), track.pt(), nsigmaTOF, centrality); + //if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(iSpecies, 1)) { + // return false; + //} + //mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaTOF"), track.pt(), nsigmaTOF); + + return true; + } + + template + void fillNucleusFlagsPdgsMc(const Tparticle& particle, nuclei::SlimCandidate& candidate) + { + candidate.pdgCode = particle.pdgCode(); + + if (particle.isPhysicalPrimary()) { + candidate.flags |= nuclei::Flags::kIsPhysicalPrimary; + + // heavy flavour mother + //if (particle.has_mothers()) { + // for (auto& motherparticle : particle.mothers_as()) { + // if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { + // flags |= kIsSecondaryFromWeakDecay; + // motherPdgCode = motherparticle.pdgCode(); + // break; + // } + // } + //} + + } else if (particle.has_mothers()) { + candidate.flags |= nuclei::Flags::kIsSecondaryFromWeakDecay; + for (auto& motherparticle : particle.template mothers_as()) { + candidate.motherPdgCode = motherparticle.pdgCode(); + } + + } else { + candidate.flags |= nuclei::Flags::kIsSecondaryFromMaterial; + } + + mFilledMcParticleIds.emplace_back(particle.globalIndex()); + } + + template + void fillNucleusFlagsPdgs(const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) + { + candidate.flags = static_cast((track.pidForTracking() & 0xF) << 12); + candidate.flags |= iSpecies == nuclei::Species::kPr ? nuclei::Flags::kProton : + iSpecies == nuclei::Species::kDe ? nuclei::Flags::kDeuteron : + iSpecies == nuclei::Species::kTr ? nuclei::Flags::kTriton : + iSpecies == nuclei::Species::kHe ? nuclei::Flags::kHe3 : + iSpecies == nuclei::Species::kAl ? nuclei::Flags::kHe4 : 0; + + if (track.hasTOF()) { + candidate.flags |= nuclei::Flags::kHasTOF; + } + if (track.hasTRD()) { + candidate.flags |= nuclei::Flags::kHasTRD; + } + if (!collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + candidate.flags |= nuclei::Flags::kITSrof; + } + } + + template + void fillNucleusGeneratedVariables(const Tparticle& particle, nuclei::SlimCandidate& candidate) + { + candidate.ptGenerated = particle.pt() * (particle.pdgCode() > 0 ? 1.f : -1.f); + candidate.etaGenerated = particle.eta(); + candidate.phiGenerated = particle.phi(); + } + + template + nuclei::SlimCandidate fillCandidate(const int iSpecies, Tcollision const& collision, Ttrack const& track) + { + if (!nuclei::checkSpeciesValidity(iSpecies)) { + std::runtime_error("species contains invalid nucleus index"); + } + nuclei::SlimCandidate candidate = {.pt = track.pt() * track.sign(), + .eta = track.eta(), + .phi = track.phi(), + .tpcInnerParam = track.tpcInnerParam(), + .clusterSizesITS = track.itsClusterSizes(), + .TPCsignal = track.tpcSignal(), + //.beta = mPidManagers[iSpecies].getBetaTOF(track), + .beta = 0.f, + .DCAxy = track.dcaXY(), + .DCAz = track.dcaZ(), + .flags = 0, + .pdgCode = 0, + .motherPdgCode = 0, + .ptGenerated = 0.f, // to be filled for mc + .etaGenerated = 0.f, + .phiGenerated = 0.f, + .centrality = nuclei::getCentrality(collision, cfgCentralityEstimator) + }; + + fillNucleusFlagsPdgs(iSpecies, collision, track, candidate); + + if constexpr (isMc) { + if (track.has_mcParticle()) { + + const auto& particle = track.mcParticle(); + fillNucleusFlagsPdgsMc(particle, candidate); + fillNucleusGeneratedVariables(particle, candidate); + } + } + + mNucleiCandidates.emplace_back(candidate); + return candidate; + + } + + template + void fillHistograms(const int iSpecies, const nuclei::SlimCandidate& candidate) + { + if (!nuclei::checkSpeciesValidity(iSpecies)) { + std::runtime_error("species contains invalid nucleus index"); + } + + if (isGenerated) { + const float ptGenerated = (iSpecies == nuclei::Species::kPr || iSpecies == nuclei::Species::kDe || iSpecies == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f; + mNucleiHistogramRegistries[iSpecies].fill(HIST("hPtGenerated"), ptGenerated); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h3PtVsEtaVsCentralityGenerated"), ptGenerated, candidate.etaGenerated, candidate.centrality); + mNucleiHistogramRegistries[iSpecies].fill(HIST("hPhiPhiVsEtaVsCentralityGenerated"), candidate.phiGenerated, candidate.etaGenerated, candidate.centrality); + } else { + mNucleiHistogramRegistries[iSpecies].fill(HIST("hPt"), candidate.pt); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h3PtVsEtaVsCentralityReconstructed"), candidate.pt, candidate.eta, candidate.centrality); + mNucleiHistogramRegistries[iSpecies].fill(HIST("hPhiPhiVsEtaVsCentralityReconstructed"), candidate.phi, candidate.eta, candidate.centrality); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h3DCAxyVsPtVsCentrality"), candidate.DCAxy, candidate.pt, candidate.centrality); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h3DCAzVsPtVsCentrality"), candidate.DCAz, candidate.pt, candidate.centrality); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h3BetaVsPtVsCentrality"), candidate.beta, candidate.pt, candidate.centrality); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h3dEdxVsPVsCentrality"), candidate.TPCsignal, candidate.pt, candidate.centrality); + mNucleiHistogramRegistries[iSpecies].fill(HIST("h3ClusterSizeVsPtVsCentrality"), mPidManagers[iSpecies].getClusterSizeCosLambdaITS(candidate.clusterSizesITS, candidate.eta), candidate.pt, candidate.centrality); + } + } + + void processMc(const Collision& collision, const TrackCandidatesMC& tracks, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles) + { + mNucleiCandidates.clear(); + mFilledMcParticleIds.clear(); + + LOG(info) << "Running"; + + for (const auto iSpecies: mSpeciesToProcess) { + mPidManagers[iSpecies].initMcResponseITS(); + } + + if (!nuclei::eventSelection(collision, mQaHistograms, cfgEventSelections, cfgCutVertex)) { + return; + } + + for (const auto& track: tracks) + { + for (const auto iSpecies: mSpeciesToProcess) { + + mNucleiHistogramRegistries[iSpecies].fill(HIST("hTrackSelections"), nuclei::trackSelection::kNoCuts); + if (!trackSelection(track)) { + continue; + } + + mNucleiHistogramRegistries[iSpecies].fill(HIST("hTrackSelections"), nuclei::trackSelection::kTrackCuts); + if (!pidSelection(iSpecies, track, collision)) { + continue; + } + + mNucleiHistogramRegistries[iSpecies].fill(HIST("hTrackSelections"), nuclei::trackSelection::kPidCuts); + + nuclei::SlimCandidate candidate; + if (track.has_mcParticle()) { + const auto& particle = track.mcParticle(); + if (particle.y() - cfgRapidityCenterMass < cfgRapidityMin || particle.y() - cfgRapidityCenterMass > cfgRapidityMax) { + continue; + } + candidate = fillCandidate(iSpecies, collision, track); + fillHistograms(iSpecies, candidate); + } else { + candidate = fillCandidate(iSpecies, collision, track); + } + + fillHistograms(iSpecies, candidate); + } + } + + for (const auto& particle: mcParticles) + { + if (std::find(mFilledMcParticleIds.begin(), mFilledMcParticleIds.end(), particle.globalIndex()) != mFilledMcParticleIds.end()) { + continue; + } + + nuclei::SlimCandidate candidate; + fillNucleusFlagsPdgsMc(particle, candidate); + fillNucleusGeneratedVariables(particle, candidate); + mNucleiCandidates.emplace_back(candidate); + + const int iSpecies = nuclei::getSpeciesFromPdg(particle.pdgCode()); + fillHistograms(iSpecies, candidate); + } + + if (!cfgFillTable) { + return; + } + + for (const auto& candidate: mNucleiCandidates) + { + mNucleiTableRed( + candidate.pt, + candidate.eta, + candidate.phi, + candidate.tpcInnerParam, + candidate.clusterSizesITS, + candidate.TPCsignal, + candidate.beta, + candidate.DCAxy, + candidate.DCAz, + candidate.flags, + candidate.pdgCode, + candidate.motherPdgCode + ); + } + } + PROCESS_SWITCH(nucleiQC, processMc, "Mc analysis", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h new file mode 100644 index 00000000000..2e483926e39 --- /dev/null +++ b/PWGLF/Utils/nucleiUtils.h @@ -0,0 +1,539 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef PWGLF_UTILS_NUCLEIUTILS_H_ +#define PWGLF_UTILS_NUCLEIUTILS_H_ + +#include +#include + +#include "Framework/HistogramRegistry.h" +#include "Framework/HistogramSpec.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::physics; + +struct NucleusCandidate { + int globalIndex; + int collTrackIndex; + float pt; + float eta; + float phi; + float tpcInnerParam; + float beta; + float zVertex; + int nContrib; + float DCAxy; + float DCAz; + float TPCsignal; + float ITSchi2; + float TPCchi2; + float TOFchi2; + std::array nSigmaTPC; + std::array tofMasses; + bool fillTree; + bool fillDCAHist; + bool correctPV; + bool isSecondary; + bool fromWeakDecay; + uint16_t flags; + uint8_t TPCfindableCls; + uint8_t TPCcrossedRows; + uint8_t ITSclsMap; + uint8_t TPCnCls; + uint8_t TPCnClsShared; + uint8_t ITSnCls; + uint32_t clusterSizesITS; +}; + +struct NucleusCandidateFlow { + float centFV0A; + float centFT0M; + float centFT0A; + float centFT0C; + float psiFT0A; + float psiFT0C; + float psiTPC; + float psiTPCl; + float psiTPCr; + float qFT0A; + float qFT0C; + float qTPC; + float qTPCl; + float qTPCr; +}; + +namespace nuclei +{ + +struct SlimCandidate { + float pt = -999.f; + float eta = -999.f; + float phi = -999.f; + float tpcInnerParam = -999.f; + uint32_t clusterSizesITS = 0.f; + float TPCsignal = -999.f; + float beta = -999.f; + float DCAxy = -999.f; + float DCAz = -999.f; + uint16_t flags = 0; + int pdgCode = 0; + int motherPdgCode = 0; + float ptGenerated = -999.f; + float etaGenerated = -999.f; + float phiGenerated = -999.f; + float centrality = -1.f; +}; + +enum Species { + kPr = 0, + kDe = 1, + kTr = 2, + kHe = 3, + kAl = 4, + kNspecies = 5 +}; + +enum Flags { + kProton = BIT(0), + kDeuteron = BIT(1), + kTriton = BIT(2), + kHe3 = BIT(3), + kHe4 = BIT(4), + kHasTOF = BIT(5), + kHasTRD = BIT(6), + kIsAmbiguous = BIT(7), /// just a placeholder now + kITSrof = BIT(8), + kIsPhysicalPrimary = BIT(9), /// MC flags starting from the second half of the short + kIsSecondaryFromMaterial = BIT(10), + kIsSecondaryFromWeakDecay = BIT(11) /// the last 4 bits are reserved for the PID in tracking +}; + +int getSpeciesFromPdg(int pdg) { + switch (std::abs(pdg)) { + case PDG_t::kProton: + return Species::kPr; + case o2::constants::physics::Pdg::kDeuteron: + return Species::kDe; + case o2::constants::physics::Pdg::kTriton: + return Species::kTr; + case o2::constants::physics::Pdg::kHelium3: + return Species::kHe; + case o2::constants::physics::Pdg::kAlpha: + return Species::kAl; + default: + return -1; + } +} + +bool checkSpeciesValidity(const int species) { + if (species < 0 || species > Species::kNspecies) { + return false; + } + return true; +} + +constexpr int speciesToProcessDefault[Species::kNspecies][1]{ + {0}, + {0}, + {0}, + {0}, + {0}}; +constexpr int useCentralTpcCalibrationDefault[Species::kNspecies][1]{ + {1}, + {1}, + {1}, + {1}, + {1}}; +constexpr double bbMomScalingDefault[Species::kNspecies][2]{ + {1., 1.}, + {1., 1.}, + {1., 1.}, + {1., 1.}, + {1., 1.}}; +constexpr double betheBlochDefault[Species::kNspecies][6]{ + {-136.71, 0.441, 0.2269, 1.347, 0.8035, 0.09}, + {-136.71, 0.441, 0.2269, 1.347, 0.8035, 0.09}, + {-239.99, 1.155, 1.099, 1.137, 1.006, 0.09}, + {-321.34, 0.6539, 1.591, 0.8225, 2.363, 0.09}, + {-586.66, 1.859, 4.435, 0.282, 3.201, 0.09}}; +constexpr double nSigmaTPCdefault[Species::kNspecies][2]{ + {-5., 5.}, + {-5., 5.}, + {-5., 5.}, + {-5., 5.}, + {-5., 5.}}; +constexpr double nSigmaTOFdefault[Species::kNspecies][2]{ + {-5., 5.}, + {-5., 5.}, + {-5., 5.}, + {-5., 5.}, + {-5., 5.}}; +constexpr double DCAcutDefault[Species::kNspecies][2]{ + {1., 1.}, + {1., 1.}, + {1., 1.}, + {1., 1.}, + {1., 1.}}; +constexpr int TreeConfigDefault[Species::kNspecies][2]{ + {0, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {0, 0}}; +constexpr int FlowHistDefault[Species::kNspecies][1]{ + {0}, + {0}, + {0}, + {0}, + {0}}; +constexpr int DCAHistDefault[Species::kNspecies][2]{ + {0, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {0, 0}}; +constexpr double DownscalingDefault[Species::kNspecies][1]{ + {1.}, + {1.}, + {1.}, + {1.}, + {1.}}; +// constexpr bool storeTreesDefault[Species::kNspecies]{false, false, false, false, false}; +constexpr float charges[Species::kNspecies]{1.f, 1.f, 1.f, 2.f, 2.f}; +constexpr float masses[Species::kNspecies]{MassProton, MassDeuteron, MassTriton, MassHelium3, MassAlpha}; +static const std::vector matter{"M", "A"}; +static const std::vector pidName{"TPC", "TOF"}; +static const std::vector names{"proton", "deuteron", "triton", "He3", "alpha"}; +static const std::vector treeConfigNames{"Filter trees", "Use TOF selection"}; +static const std::vector flowConfigNames{"Save flow hists"}; +static const std::vector DCAConfigNames{"Save DCA hist", "Matter/Antimatter"}; +static const std::vector nSigmaConfigName{"nsigma_min", "nsigma_max"}; +static const std::vector nDCAConfigName{"max DCAxy", "max DCAz"}; +static const std::vector DownscalingConfigName{"Fraction of kept candidates"}; +static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; +static const std::vector chargeLabelNames{"Positive", "Negative"}; + +float pidCutTPC[Species::kNspecies][2]; //[species][lower/upper limit] + +o2::base::MatLayerCylSet* lut = nullptr; + +std::vector candidates; +std::vector candidates_flow; + +enum centDetectors { + kFV0A = 0, + kFT0M = 1, + kFT0A = 2, + kFT0C = 3 +}; + +static const std::vector centDetectorNames{"FV0A", "FT0M", "FT0A", "FT0C"}; + +// Event selections + +enum evSel { + kTVX = 0, + kZvtx, + kTFborder, + kITSROFborder, + kNoSameBunchPileup, + kIsGoodZvtxFT0vsPV, + kIsGoodITSLayersAll, + kIsEPtriggered, + kNevSels +}; + +static const std::vector eventSelectionTitle{"Event selections"}; +static const std::vector eventSelectionLabels{"TVX", "Z vtx", "TF border", "ITS ROF border", "No same-bunch pile-up", "kIsGoodZvtxFT0vsPV", "isGoodITSLayersAll", "isEPtriggered"}; + +constexpr int EvSelDefault[8][1]{ + {1}, + {1}, + {0}, + {0}, + {0}, + {0}, + {0}, + {0}}; + +template // move to nucleiUtils +bool eventSelection(const Tcollision& collision, HistogramRegistry& registry, LabeledArray eventSelections, const float cutVertex) +{ + if (!registry.contains(HIST("hEventSelections"))) { + registry.add("hEventSelections", "hEventSelections", {HistType::kTH1D, {{evSel::kNevSels + 1, -0.5f, static_cast(evSel::kNevSels) + 0.5f}}}); + } + registry.fill(HIST("hEventSelections"), 0); + + if (eventSelections.get(evSel::kTVX) && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kTVX + 1); + + if (eventSelections.get(evSel::kZvtx) && std::abs(collision.posZ()) > cutVertex) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kZvtx + 1); + + if (eventSelections.get(evSel::kTFborder) && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kTFborder + 1); + + if (eventSelections.get(evSel::kITSROFborder) && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kITSROFborder + 1); + + if (eventSelections.get(evSel::kNoSameBunchPileup) && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kNoSameBunchPileup + 1); + + if (eventSelections.get(evSel::kIsGoodZvtxFT0vsPV) && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kIsGoodZvtxFT0vsPV + 1); + + if (eventSelections.get(evSel::kIsGoodITSLayersAll) && !collision.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kIsGoodITSLayersAll + 1); + + if constexpr ( + requires { + collision.triggereventep(); + }) { + if (eventSelections.get(evSel::kIsEPtriggered) && !collision.triggereventep()) { + return false; + } + registry.fill(HIST("hEventSelections"), evSel::kIsEPtriggered + 1); + } + + return true; +} + +template +float getCentrality(Tcollision const& collision, const int centralityEstimator) +{ + if constexpr (o2::aod::HasCentrality) { + if (centralityEstimator == centDetectors::kFV0A) { + return collision.centFV0A(); + } else if (centralityEstimator == centDetectors::kFT0M) { + return collision.centFT0M(); + } else if (centralityEstimator == centDetectors::kFT0A) { + return collision.centFT0A(); + } else if (centralityEstimator == centDetectors::kFT0C) { + return collision.centFT0C(); + } else { + LOG(warning) << "Centrality estimator not valid. Possible values: (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3). Centrality set to -1."; + } + } + return -1.f; +} + +// Track selections + +enum trackSelection { + kNoCuts = 0, + kTrackCuts = 1, + kPidCuts = 2, + kNtrackSelections = 3 +}; +static const std::array(trackSelection::kNtrackSelections)> trackSelectionLabels{"All", "Track cuts", "PID cuts"}; + +HistogramRegistry createHistogramRegistryNucleus(const int iSpecies) { + + if (!checkSpeciesValidity(iSpecies)) { + std::runtime_error("species contains invalid nucleus index"); + } + + HistogramRegistry registry{ + names[iSpecies].c_str(), + { + {"hTrackSelections", *Form("%s", names[iSpecies].c_str()) + " track selections; Selection step; Counts", {HistType::kTH1D, {{trackSelection::kNtrackSelections, -0.5f, static_cast(trackSelection::kNtrackSelections) - 0.5f}}}}, + {"hPtReconstructed", *Form("%s", names[iSpecies].c_str()) + " - reconstructed variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); Counts", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, + {"h3PtVsEtaVsCentralityReconstructed", *Form("%s", names[iSpecies].c_str()) + " - reconstructed variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}}}, + {"h3PhiVsEtaVsCentralityReconstructed", *Form("%s", names[iSpecies].c_str()) + " - reconstructed variables; #phi (radians); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}}}, + {"h3DCAxyVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{xy} (cm); CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, + {"h3DCAzVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{z} (cm); CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, + {"h3NsigmaTPC_preselectionVsCentrality", *Form("Nsigma%s TPC distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TPC}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}}}, + {"h3NsigmaTPCVsCentrality", *Form("Nsigma%s TPC distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TPC}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, + {"h3NSigmaITS_preselectionVsCentrality", *Form("Nsigma%s ITS distribution; signed #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{ITS}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}}}, + {"h3NSigmaITSVsCentrality", *Form("Nsigma%s ITS distribution; signed #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{ITS}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}}}, + {"h3NSigmaTOF_preselectionVsCentrality", *Form("Nsigma%s TOF distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TOF}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}}}, + {"h3NSigmaTOFVsCentrality", *Form("Nsigma%s TOF distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TOF}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, + {"h3BetaVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); #beta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {100, 0.0f, 1.0f}, {20, 0.0f, 100.0f}}}}, + {"h3dEdxVsPVsCentrality", *Form("dEdx distribution for %s; #it{p} (GeV/#it{c}); d#it{E}/d#it{x} (a.u.);", names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{200, -6.0f, 6.0f}, {100, 0.0f, 2000.0f}, {20, 0.0f, 100.0f}}}}, + {"h3ClusterSizeVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); Cluster size ITS; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {90, 0.f, 15.f}, {20, 0.0f, 100.0f}}}}, + {"hPtGenerated", *Form("%s", names[iSpecies].c_str()) + " - generated variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); Counts", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, + {"h3PtVsEtaVsCentralityGenerated", *Form("%s", names[iSpecies].c_str()) + " - generated variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}}}, + {"h3PhiVsEtaVsCentralityGenerated", *Form("%s", names[iSpecies].c_str()) + " - generated variables; #phi (radians); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}}}, + }, + OutputObjHandlingPolicy::AnalysisObject, + false, + true}; + + for (size_t iSel = 0; iSel < trackSelection::kNtrackSelections; iSel++) { + registry.get(HIST("hTrackSelections"))->GetXaxis()->SetBinLabel(iSel + 1, trackSelectionLabels[iSel].c_str()); + } + + return registry; +}; + +// PID manager class + +class PidManager { + + public: + PidManager(const int species, const float* tpcBetheBlochParams = nullptr) + : mSpecies(species) + { + if (!checkSpeciesValidity(species)) { + std::runtime_error("species contains invalid nucleus index"); + } + + if (!tpcBetheBlochParams) { + mUseTpcCentralCalibration = true; + return; + } + + for (int i = 0; i < 6; i++) { + mTpcBetheBlochParams[i] = tpcBetheBlochParams[i]; + } + } + PidManager() = default; + ~PidManager() = default; + + + // TOF + template + float getBetaTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; + } + float beta = o2::pid::tof::Beta::GetBeta(track); + return std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + } + + template + float getMassTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; + } + const float charge{1.f + static_cast(mSpecies == Species::kHe || mSpecies == Species::kAl)}; + const float beta = getBetaTOF(track); + return track.tpcInnerParam() * charge * std::sqrt(1.f / (beta * beta) - 1.f); + } + + template + float getNSigmaTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; + } + + switch (mSpecies) { + case Species::kPr: return track.tofNSigmaPr(); + case Species::kDe: return track.tofNSigmaDe(); + case Species::kTr: return track.tofNSigmaTr(); + case Species::kHe: return track.tofNSigmaHe(); + case Species::kAl: return track.tofNSigmaAl(); + default: return -999.f; + } + } + + // ITS + void initMcResponseITS() { + mResponseITS.setMCDefaultParameters(); + } + + template + float getClusterSizeCosLambdaITS(const Ttrack& track) + { + return mResponseITS.averageClusterSize(track.itsClusterSizes()) / std::cosh(track.eta()); + } + + float getClusterSizeCosLambdaITS(const u_int32_t clusterSizesITS, const float eta) + { + return mResponseITS.averageClusterSize(clusterSizesITS) / std::cosh(eta); + } + + template + float getNSigmaITS(const Ttrack& track) + { + switch (mSpecies) { + case Species::kPr: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kDe: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kTr: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kHe: return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); + case Species::kAl: return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); + default: return -999.f; + } + } + + // TPC + float getExpectedTPCsignal(const float p) + { + if (!mUseTpcCentralCalibration) { + return -999.f; + } + float pScaled = p * mMomScaling[0] + mMomScaling[1]; + float betaGamma = pScaled / masses[mSpecies]; + return tpc::BetheBlochAleph(betaGamma, + mTpcBetheBlochParams[0], + mTpcBetheBlochParams[1], + mTpcBetheBlochParams[2], + mTpcBetheBlochParams[3], + mTpcBetheBlochParams[4]); + } + + template + float getNSigmaTPC(const Ttrack& track) + { + if (!mUseTpcCentralCalibration) { + return getNSigmaTPCcentral(track); + } + float expectedSignal = getExpectedTPCsignal(track.tpcInnerParam()); + float resolution = mTpcBetheBlochParams[5]; + return (track.tpcSignal() - expectedSignal) / (expectedSignal * resolution); + } + + protected: + // TPC + template + float getNSigmaTPCcentral(const Ttrack& track) { + switch (mSpecies) { + case Species::kPr: return track.tpcNSigmaPr(); + case Species::kDe: return track.tpcNSigmaDe(); + case Species::kTr: return track.tpcNSigmaTr(); + case Species::kHe: return track.tpcNSigmaHe(); + case Species::kAl: return track.tpcNSigmaAl(); + default: return -999.f; + } + } + + private: + float mTpcBetheBlochParams[6]; + bool mUseTpcCentralCalibration = true; // this just becomes a check for the null pointer in the parameters + o2::aod::ITSResponse mResponseITS; + float mMomScaling[2]{1., 0.}; + int mSpecies; +}; + + +} // namespace nuclei + +#endif // PWGLF_UTILS_NUCLEIUTILS_H_ From ea32ba0e246806a687a3c9f154827e2b92677cea Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 2 Oct 2025 18:02:20 +0200 Subject: [PATCH 2/9] moved to the appropriate directory, added pid manager, central collision selection and histogram registry generator --- PWGLF/Utils/nucleiUtils.h | 90 ++++++++++++++++++--------------------- 1 file changed, 41 insertions(+), 49 deletions(-) diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h index 2e483926e39..eeb70aa51b2 100644 --- a/PWGLF/Utils/nucleiUtils.h +++ b/PWGLF/Utils/nucleiUtils.h @@ -119,20 +119,14 @@ enum Flags { kIsSecondaryFromWeakDecay = BIT(11) /// the last 4 bits are reserved for the PID in tracking }; -int getSpeciesFromPdg(int pdg) { +constexpr int getSpeciesFromPdg(int pdg) { switch (std::abs(pdg)) { - case PDG_t::kProton: - return Species::kPr; - case o2::constants::physics::Pdg::kDeuteron: - return Species::kDe; - case o2::constants::physics::Pdg::kTriton: - return Species::kTr; - case o2::constants::physics::Pdg::kHelium3: - return Species::kHe; - case o2::constants::physics::Pdg::kAlpha: - return Species::kAl; - default: - return -1; + case PDG_t::kProton: return Species::kPr; + case o2::constants::physics::Pdg::kDeuteron: return Species::kDe; + case o2::constants::physics::Pdg::kTriton: return Species::kTr; + case o2::constants::physics::Pdg::kHelium3: return Species::kHe; + case o2::constants::physics::Pdg::kAlpha: return Species::kAl; + default: return -1; } } @@ -212,9 +206,11 @@ constexpr double DownscalingDefault[Species::kNspecies][1]{ // constexpr bool storeTreesDefault[Species::kNspecies]{false, false, false, false, false}; constexpr float charges[Species::kNspecies]{1.f, 1.f, 1.f, 2.f, 2.f}; constexpr float masses[Species::kNspecies]{MassProton, MassDeuteron, MassTriton, MassHelium3, MassAlpha}; +constexpr int pdgCodes[Species::kNspecies]{PDG_t::kProton, o2::constants::physics::Pdg::kDeuteron, o2::constants::physics::Pdg::kTriton, o2::constants::physics::Pdg::kHelium3, o2::constants::physics::Pdg::kAlpha}; +static constexpr std::string_view cNames[] = {"proton", "deuteron", "triton", "He3", "alpha"}; +static const std::vector names{"proton", "deuteron", "triton", "He3", "alpha"}; static const std::vector matter{"M", "A"}; static const std::vector pidName{"TPC", "TOF"}; -static const std::vector names{"proton", "deuteron", "triton", "He3", "alpha"}; static const std::vector treeConfigNames{"Filter trees", "Use TOF selection"}; static const std::vector flowConfigNames{"Save flow hists"}; static const std::vector DCAConfigNames{"Save DCA hist", "Matter/Antimatter"}; @@ -270,10 +266,17 @@ constexpr int EvSelDefault[8][1]{ template // move to nucleiUtils bool eventSelection(const Tcollision& collision, HistogramRegistry& registry, LabeledArray eventSelections, const float cutVertex) { + if (!registry.contains(HIST("hVtxZBefore"))) { + registry.add("hVtxZBefore", "Vertex distribution in Z before selections;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}); + } + if (!registry.contains(HIST("hVtxZ"))) { + registry.add("hVtxZ", "Vertex distribution in Z;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}); + } if (!registry.contains(HIST("hEventSelections"))) { registry.add("hEventSelections", "hEventSelections", {HistType::kTH1D, {{evSel::kNevSels + 1, -0.5f, static_cast(evSel::kNevSels) + 0.5f}}}); } registry.fill(HIST("hEventSelections"), 0); + registry.fill(HIST("hVtxZBefore"), collision.posZ()); if (eventSelections.get(evSel::kTVX) && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { return false; @@ -319,6 +322,7 @@ bool eventSelection(const Tcollision& collision, HistogramRegistry& registry, La } registry.fill(HIST("hEventSelections"), evSel::kIsEPtriggered + 1); } + registry.fill(HIST("hVtxZ"), collision.posZ()); return true; } @@ -352,43 +356,36 @@ enum trackSelection { }; static const std::array(trackSelection::kNtrackSelections)> trackSelectionLabels{"All", "Track cuts", "PID cuts"}; -HistogramRegistry createHistogramRegistryNucleus(const int iSpecies) { +template +void createHistogramRegistryNucleus(HistogramRegistry& registry) { - if (!checkSpeciesValidity(iSpecies)) { + constexpr int index = iSpecies; + if (!checkSpeciesValidity(index)) { std::runtime_error("species contains invalid nucleus index"); } - HistogramRegistry registry{ - names[iSpecies].c_str(), - { - {"hTrackSelections", *Form("%s", names[iSpecies].c_str()) + " track selections; Selection step; Counts", {HistType::kTH1D, {{trackSelection::kNtrackSelections, -0.5f, static_cast(trackSelection::kNtrackSelections) - 0.5f}}}}, - {"hPtReconstructed", *Form("%s", names[iSpecies].c_str()) + " - reconstructed variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); Counts", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, - {"h3PtVsEtaVsCentralityReconstructed", *Form("%s", names[iSpecies].c_str()) + " - reconstructed variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}}}, - {"h3PhiVsEtaVsCentralityReconstructed", *Form("%s", names[iSpecies].c_str()) + " - reconstructed variables; #phi (radians); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}}}, - {"h3DCAxyVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{xy} (cm); CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, - {"h3DCAzVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{z} (cm); CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, - {"h3NsigmaTPC_preselectionVsCentrality", *Form("Nsigma%s TPC distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TPC}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}}}, - {"h3NsigmaTPCVsCentrality", *Form("Nsigma%s TPC distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TPC}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, - {"h3NSigmaITS_preselectionVsCentrality", *Form("Nsigma%s ITS distribution; signed #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{ITS}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}}}, - {"h3NSigmaITSVsCentrality", *Form("Nsigma%s ITS distribution; signed #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{ITS}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}}}, - {"h3NSigmaTOF_preselectionVsCentrality", *Form("Nsigma%s TOF distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TOF}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}}}, - {"h3NSigmaTOFVsCentrality", *Form("Nsigma%s TOF distribution; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); n#sigma_{TOF}(%s);", names[iSpecies].c_str(), names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}}}, - {"h3BetaVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); #beta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {100, 0.0f, 1.0f}, {20, 0.0f, 100.0f}}}}, - {"h3dEdxVsPVsCentrality", *Form("dEdx distribution for %s; #it{p} (GeV/#it{c}); d#it{E}/d#it{x} (a.u.);", names[iSpecies].c_str()) + "CentralityFT0C (%)", {HistType::kTH3F, {{200, -6.0f, 6.0f}, {100, 0.0f, 2000.0f}, {20, 0.0f, 100.0f}}}}, - {"h3ClusterSizeVsPtVsCentrality", *Form("%s", names[iSpecies].c_str()) + "; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); Cluster size ITS; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {90, 0.f, 15.f}, {20, 0.0f, 100.0f}}}}, - {"hPtGenerated", *Form("%s", names[iSpecies].c_str()) + " - generated variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); Counts", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, - {"h3PtVsEtaVsCentralityGenerated", *Form("%s", names[iSpecies].c_str()) + " - generated variables; #it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}}}, - {"h3PhiVsEtaVsCentralityGenerated", *Form("%s", names[iSpecies].c_str()) + " - generated variables; #phi (radians); #eta; CentralityFT0C (%)", {HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}}}, - }, - OutputObjHandlingPolicy::AnalysisObject, - false, - true}; + registry.add(fmt::format("{}/hTrackSelections", cNames[index]).c_str(), (fmt::format("{} track selections;", cNames[index]) + std::string("Selection step; Counts")).c_str(), HistType::kTH1D, {{trackSelection::kNtrackSelections, -0.5f, static_cast(trackSelection::kNtrackSelections) - 0.5f}}); + registry.add(fmt::format("{}/hPtReconstructed", cNames[index]).c_str(), (fmt::format("{} - reconstructed variables;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); Counts")).c_str(), HistType::kTH1F, {{240, -6.0f, 6.0f}}); + registry.add(fmt::format("{}/h3PtVsEtaVsCentralityReconstructed", cNames[index]).c_str(), (fmt::format("{} - reconstructed variables;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3PhiVsEtaVsCentralityReconstructed", cNames[index]).c_str(), (fmt::format("{} - reconstructed variables;", cNames[index]) + std::string("#phi (radians); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3DCAxyVsPtVsCentrality", cNames[index]).c_str(), (fmt::format(";", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{xy} (cm); CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3DCAzVsPtVsCentrality", cNames[index]).c_str(), (fmt::format("{};", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{z} (cm); CentralityFT0C (%)")).c_str() , HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3NsigmaTPC_preselectionVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TPC distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TPC}}({}); CentralityFT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3NsigmaTPCVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TPC distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TPC}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3NsigmaITS_preselectionVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} ITS distribution;", cNames[index]) + std::string("signed #it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{ITS}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3NsigmaITSVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} ITS distribution;", cNames[index]) + std::string("signed #it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{ITS}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3NsigmaTOF_preselectionVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TOF distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TOF}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3NsigmaTOFVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TOF distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TOF}}({}); Centrality FT0C (%)", cNames[index], cNames[index])).c_str(), HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3BetaVsPtVsCentrality", cNames[index]).c_str(), (fmt::format("{};", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); #beta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {100, 0.0f, 1.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3dEdxVsPVsCentrality", cNames[index]).c_str(), (fmt::format("dEdx distribution for {};", cNames[index]) + std::string("#it{p} (GeV/#it{c}); d#it{E}/d#it{x} (a.u.); Centrality FT0C (%)")).c_str(), HistType::kTH3F, {{200, -6.0f, 6.0f}, {100, 0.0f, 2000.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3ClusterSizeVsPtVsCentrality", cNames[index]).c_str(), (fmt::format("{};", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); Cluster size ITS; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {90, 0.f, 15.f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/hPtGenerated", cNames[index]).c_str(), (fmt::format("{} - generated variables;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); Counts")).c_str(), HistType::kTH1F, {{240, -6.0f, 6.0f}}); + registry.add(fmt::format("{}/h3PtVsEtaVsCentralityGenerated", cNames[index]).c_str(), (fmt::format("{} - generated variables;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3PhiVsEtaVsCentralityGenerated", cNames[index]).c_str(), (fmt::format("{} - generated variables;", cNames[index]) + std::string("#phi (radians); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}); for (size_t iSel = 0; iSel < trackSelection::kNtrackSelections; iSel++) { - registry.get(HIST("hTrackSelections"))->GetXaxis()->SetBinLabel(iSel + 1, trackSelectionLabels[iSel].c_str()); + registry.get(HIST(cNames[index]) + HIST("/hTrackSelections"))->GetXaxis()->SetBinLabel(iSel + 1, trackSelectionLabels[iSel].c_str()); } - - return registry; }; // PID manager class @@ -455,11 +452,6 @@ class PidManager { } } - // ITS - void initMcResponseITS() { - mResponseITS.setMCDefaultParameters(); - } - template float getClusterSizeCosLambdaITS(const Ttrack& track) { From ac8847f44d6d1535265985dbfbbb2518773de247 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Thu, 2 Oct 2025 18:03:13 +0200 Subject: [PATCH 3/9] added new qc task for nuclei, with tablefor dca studies --- PWGLF/TableProducer/QC/nucleiQC.cxx | 251 +++++++++++++++++----------- 1 file changed, 156 insertions(+), 95 deletions(-) diff --git a/PWGLF/TableProducer/QC/nucleiQC.cxx b/PWGLF/TableProducer/QC/nucleiQC.cxx index f0cfbdf52f8..afd71da2204 100644 --- a/PWGLF/TableProducer/QC/nucleiQC.cxx +++ b/PWGLF/TableProducer/QC/nucleiQC.cxx @@ -59,6 +59,7 @@ #include "Framework/ASoAHelpers.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" +#include "Framework/StaticFor.h" #include "ReconstructionDataFormats/Track.h" @@ -75,13 +76,13 @@ using namespace o2::framework; struct nucleiQC { -//using TrackCandidates = soa::Join; - using TrackCandidates = soa::Join; -//using TrackCandidatesMC = soa::Join; - using TrackCandidatesMC = soa::Join; - using Collision = soa::Join::iterator; + using TrackCandidates = soa::Join; + using TrackCandidatesMC = soa::Join; + using Collision = soa::Join::iterator; + Preslice mMcParticlesPerCollision = o2::aod::mcparticle::mcCollisionId; Configurable cfgFillTable{"cfgFillTable", true, "Fill output tree"}; + Configurable cfgDoCheckPdgCode{"cfgDoCheckPdgCode", true, "Should you only select tracks associated to a mc particle with the correct PDG code?"}; Configurable> cfgSpeciesToProcess{"cfgSpeciesToProcess", {nuclei::speciesToProcessDefault[0], nuclei::Species::kNspecies, 1, nuclei::names, {"processNucleus"}}, "Nuclei to process"}; Configurable> cfgEventSelections{"cfgEventSelections", {nuclei::EvSelDefault[0], 8, 1, nuclei::eventSelectionLabels, nuclei::eventSelectionTitle}, "Event selections"}; Configurable cfgCentralityEstimator{"cfgCentralityEstimator", 0, "Centrality estimator (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3)"}; @@ -101,8 +102,15 @@ struct nucleiQC { Configurable> cfgNsigmaTPC{"cfgNsigmaTPC", {nuclei::nSigmaTPCdefault[0], nuclei::Species::kNspecies, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"}; Configurable> cfgNsigmaTOF{"cfgNsigmaTOF", {nuclei::nSigmaTOFdefault[0], nuclei::Species::kNspecies, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"}; - HistogramRegistry mQaHistograms { - "QA", + // CCDB options + Configurable cfgMaterialCorrection{"cfgMaterialCorrection", static_cast(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Type of material correction"}; + Configurable cfgCCDBurl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Service mCcdb; + int mRunNumber = 0; + float mBz = 0.f; + + HistogramRegistry mHistograms { + "histos", { {"hEventSelections", "Event selections; Selection step; Counts", {HistType::kTH1D, {{nuclei::evSel::kNevSels + 1, -0.5f, static_cast(nuclei::evSel::kNevSels) + 0.5f}}}}, {"hVtxZBefore", "Vertex distribution in Z before selections;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}}, @@ -112,54 +120,71 @@ struct nucleiQC { false, true }; - std::array(nuclei::Species::kNspecies)> mNucleiHistogramRegistries; std::vector mSpeciesToProcess; Produces mNucleiTableRed; std::vector mNucleiCandidates; std::vector mFilledMcParticleIds; + o2::dataformats::DCA mDcaInfoCov; + o2::track::TrackParametrizationWithError mTrackParCov; std::array(nuclei::Species::kNspecies)> mPidManagers; void init(o2::framework::InitContext&) { - LOG(info) << "Check init"; + mCcdb->setURL(cfgCCDBurl); + mCcdb->setCaching(true); + mCcdb->setLocalObjectValidityChecking(); + mCcdb->setFatalWhenNull(false); + nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get("GLO/Param/MatLUT")); for (int iSel = 0; iSel < nuclei::evSel::kNevSels; iSel++) { - mQaHistograms.get(HIST("hEventSelections"))->GetXaxis()->SetBinLabel(iSel + 1, nuclei::eventSelectionLabels[iSel].c_str()); + mHistograms.get(HIST("hEventSelections"))->GetXaxis()->SetBinLabel(iSel + 1, nuclei::eventSelectionLabels[iSel].c_str()); } - LOG(info) << " check after hist ev sel"; - for (int iSpecies = 0; iSpecies < static_cast(nuclei::Species::kNspecies); iSpecies++) { if (cfgSpeciesToProcess->get(iSpecies) == 1) { mSpeciesToProcess.emplace_back(iSpecies); } } - LOG(info) << " check after species to be processed were saved"; + static_for<0, nuclei::kNspecies-1> ([&] (auto iSpecies) { + + constexpr int iSpeciesCt = decltype(iSpecies)::value; + const int iSpeciesRt = iSpeciesCt; - for (auto iSpecies: mSpeciesToProcess) { + if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpeciesCt) == mSpeciesToProcess.end()) { + return; + } float tpcBetheBlochParams[6]; for (int iParam = 0; iParam < 6; iParam++) { - tpcBetheBlochParams[iParam] = cfgBetheBlochParams->get(iSpecies, iParam); + tpcBetheBlochParams[iParam] = cfgBetheBlochParams->get(iSpeciesRt, iParam); } - LOG(info) << " check after bb params, ispecies: " << iSpecies; - - mNucleiHistogramRegistries[iSpecies] = nuclei::createHistogramRegistryNucleus(iSpecies); + nuclei::createHistogramRegistryNucleus(mHistograms); - if (cfgUseCentralTpcCalibration->get(uint32_t(iSpecies), uint32_t(0)) == 0) { - mPidManagers[iSpecies] = nuclei::PidManager(nuclei::Species(iSpecies), tpcBetheBlochParams); + if (cfgUseCentralTpcCalibration->get(uint32_t(iSpeciesRt), uint32_t(0)) == 0) { + mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt, tpcBetheBlochParams); } else { - mPidManagers[iSpecies] = nuclei::PidManager(nuclei::Species(iSpecies)); + mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt); } + }); + } - LOG(info) << " check after pidManagers, ispecies: " << iSpecies; + void initCCDB(const aod::BCsWithTimestamps::iterator& bc) + { + if (mRunNumber == bc.runNumber()) { + return; } - - LOG(info) << " check end init"; + auto timestamp = bc.timestamp(); + mRunNumber = bc.runNumber(); + + o2::parameters::GRPMagField* grpmag = mCcdb->getForTimeStamp("GLO/Config/GRPMagField", timestamp); + o2::base::Propagator::initFieldFromGRP(grpmag); + o2::base::Propagator::Instance()->setMatLUT(nuclei::lut); + mBz = static_cast(grpmag->getNominalL3Field()); + LOGF(info, "Retrieved GRP for timestamp %ull (%i) with magnetic field of %1.2f kZG", timestamp, mRunNumber, mBz); } template @@ -177,32 +202,29 @@ struct nucleiQC { return true; } - template - bool pidSelection(const int iSpecies, const Ttrack& track, const Tcollision& collision) + template + bool pidSelection(const Ttrack& track, const Tcollision& collision) { - if (!nuclei::checkSpeciesValidity(iSpecies)) { + constexpr int index = iSpecies; + if (!nuclei::checkSpeciesValidity(index)) { std::runtime_error("species contains invalid nucleus index"); } const float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator); - const float nsigmaTPC = mPidManagers[iSpecies].getNSigmaTPC(track); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NsigmaTPC_preselectionVsCentrality"), track.pt(), nsigmaTPC, centrality); - if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(iSpecies, 1)) { - return false; - } - mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NsigmaTPCVsCentrality"), track.pt(), nsigmaTPC, centrality); + const float nsigmaTPC = mPidManagers[index].getNSigmaTPC(track); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPC_preselectionVsCentrality"), track.pt(), nsigmaTPC, centrality); + if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(index, 1)) return false; + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPCVsCentrality"), track.pt(), nsigmaTPC, centrality); - const float nsigmaITS = mPidManagers[iSpecies].getNSigmaITS(track); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaITS_preselectionVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); + const float nsigmaITS = mPidManagers[index].getNSigmaITS(track); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaITS_preselectionVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); // add nsigmaITS cut ? - mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaITSVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaITSVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); - //const float nsigmaTOF = mPidManagers[iSpecies].getNSigmaTOF(track); - //mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaTOF_preselectionVsCentrality"), track.pt(), nsigmaTOF, centrality); - //if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(iSpecies, 1)) { - // return false; - //} - //mNucleiHistogramRegistries[iSpecies].fill(HIST("h2NSigmaTOF"), track.pt(), nsigmaTOF); + const float nsigmaTOF = mPidManagers[index].getNSigmaTOF(track); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOF_preselectionVsCentrality"), track.pt(), nsigmaTOF, centrality); + if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(index, 1)) return false; + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOFVsCentrality"), track.pt(), nsigmaTOF, centrality); return true; } @@ -268,6 +290,22 @@ struct nucleiQC { candidate.phiGenerated = particle.phi(); } + template + void fillDcaInformation(const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) { + + const o2::math_utils::Point3D collisionVertex{collision.posX(), collision.posY(), collision.posZ()}; + + mDcaInfoCov.set(999, 999, 999, 999, 999); + setTrackParCov(track, mTrackParCov); + mTrackParCov.setPID(track.pidForTracking()); + std::array dcaInfo; + o2::base::Propagator::Instance()->propagateToDCA(collisionVertex, mTrackParCov, mBz, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfo); + + candidate.DCAxy = dcaInfo[0]; + candidate.DCAz = dcaInfo[1]; + + } + template nuclei::SlimCandidate fillCandidate(const int iSpecies, Tcollision const& collision, Ttrack const& track) { @@ -280,10 +318,9 @@ struct nucleiQC { .tpcInnerParam = track.tpcInnerParam(), .clusterSizesITS = track.itsClusterSizes(), .TPCsignal = track.tpcSignal(), - //.beta = mPidManagers[iSpecies].getBetaTOF(track), - .beta = 0.f, - .DCAxy = track.dcaXY(), - .DCAz = track.dcaZ(), + .beta = mPidManagers[iSpecies].getBetaTOF(track), + .DCAxy = 0.f, + .DCAz = 0.f, .flags = 0, .pdgCode = 0, .motherPdgCode = 0, @@ -293,6 +330,7 @@ struct nucleiQC { .centrality = nuclei::getCentrality(collision, cfgCentralityEstimator) }; + fillDcaInformation(collision, track, candidate); fillNucleusFlagsPdgs(iSpecies, collision, track, candidate); if constexpr (isMc) { @@ -309,27 +347,40 @@ struct nucleiQC { } - template - void fillHistograms(const int iSpecies, const nuclei::SlimCandidate& candidate) + template + void dispatchFillHistograms(const int iSpecies, const nuclei::SlimCandidate& candidate) { + switch (iSpecies) { + case nuclei::Species::kPr: return fillHistograms(candidate); + case nuclei::Species::kDe: return fillHistograms(candidate); + case nuclei::Species::kTr: return fillHistograms(candidate); + case nuclei::Species::kHe: return fillHistograms(candidate); + case nuclei::Species::kAl: return fillHistograms(candidate); + default: return; + } + } + + template + void fillHistograms(const nuclei::SlimCandidate& candidate) { - if (!nuclei::checkSpeciesValidity(iSpecies)) { + constexpr int index = iSpecies; + if (!nuclei::checkSpeciesValidity(index)) { std::runtime_error("species contains invalid nucleus index"); } - + if (isGenerated) { - const float ptGenerated = (iSpecies == nuclei::Species::kPr || iSpecies == nuclei::Species::kDe || iSpecies == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f; - mNucleiHistogramRegistries[iSpecies].fill(HIST("hPtGenerated"), ptGenerated); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h3PtVsEtaVsCentralityGenerated"), ptGenerated, candidate.etaGenerated, candidate.centrality); - mNucleiHistogramRegistries[iSpecies].fill(HIST("hPhiPhiVsEtaVsCentralityGenerated"), candidate.phiGenerated, candidate.etaGenerated, candidate.centrality); + const float ptGenerated = (index == nuclei::Species::kPr || index == nuclei::Species::kDe || index == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f; + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/hPtGenerated"), ptGenerated); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PtVsEtaVsCentralityGenerated"), ptGenerated, candidate.etaGenerated, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PhiVsEtaVsCentralityGenerated"), candidate.phiGenerated, candidate.etaGenerated, candidate.centrality); } else { - mNucleiHistogramRegistries[iSpecies].fill(HIST("hPt"), candidate.pt); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h3PtVsEtaVsCentralityReconstructed"), candidate.pt, candidate.eta, candidate.centrality); - mNucleiHistogramRegistries[iSpecies].fill(HIST("hPhiPhiVsEtaVsCentralityReconstructed"), candidate.phi, candidate.eta, candidate.centrality); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h3DCAxyVsPtVsCentrality"), candidate.DCAxy, candidate.pt, candidate.centrality); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h3DCAzVsPtVsCentrality"), candidate.DCAz, candidate.pt, candidate.centrality); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h3BetaVsPtVsCentrality"), candidate.beta, candidate.pt, candidate.centrality); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h3dEdxVsPVsCentrality"), candidate.TPCsignal, candidate.pt, candidate.centrality); - mNucleiHistogramRegistries[iSpecies].fill(HIST("h3ClusterSizeVsPtVsCentrality"), mPidManagers[iSpecies].getClusterSizeCosLambdaITS(candidate.clusterSizesITS, candidate.eta), candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/hPtReconstructed"), candidate.pt); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PtVsEtaVsCentralityReconstructed"), candidate.pt, candidate.eta, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PhiVsEtaVsCentralityReconstructed"), candidate.phi, candidate.eta, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3DCAxyVsPtVsCentrality"), candidate.DCAxy, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3DCAzVsPtVsCentrality"), candidate.DCAz, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3BetaVsPtVsCentrality"), candidate.beta, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3dEdxVsPVsCentrality"), candidate.TPCsignal, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3ClusterSizeVsPtVsCentrality"), mPidManagers[index].getClusterSizeCosLambdaITS(candidate.clusterSizesITS, candidate.eta), candidate.pt, candidate.centrality); } } @@ -337,67 +388,77 @@ struct nucleiQC { { mNucleiCandidates.clear(); mFilledMcParticleIds.clear(); + + auto bc = collision.template bc_as(); + initCCDB(bc); - LOG(info) << "Running"; - - for (const auto iSpecies: mSpeciesToProcess) { - mPidManagers[iSpecies].initMcResponseITS(); - } - - if (!nuclei::eventSelection(collision, mQaHistograms, cfgEventSelections, cfgCutVertex)) { + if (!nuclei::eventSelection(collision, mHistograms, cfgEventSelections, cfgCutVertex)) { return; } for (const auto& track: tracks) { - for (const auto iSpecies: mSpeciesToProcess) { - - mNucleiHistogramRegistries[iSpecies].fill(HIST("hTrackSelections"), nuclei::trackSelection::kNoCuts); - if (!trackSelection(track)) { - continue; - } + + static_for<0, nuclei::kNspecies-1> ([&] (auto iSpecies) { - mNucleiHistogramRegistries[iSpecies].fill(HIST("hTrackSelections"), nuclei::trackSelection::kTrackCuts); - if (!pidSelection(iSpecies, track, collision)) { - continue; + constexpr int iSpeciesCt = decltype(iSpecies)::value; + const int iSpeciesRt = iSpeciesCt; + + if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpeciesRt) == mSpeciesToProcess.end()) { + return; } - mNucleiHistogramRegistries[iSpecies].fill(HIST("hTrackSelections"), nuclei::trackSelection::kPidCuts); - - nuclei::SlimCandidate candidate; if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); + if (cfgDoCheckPdgCode) { + if (particle.pdgCode() != nuclei::pdgCodes[iSpeciesRt]) return; + } if (particle.y() - cfgRapidityCenterMass < cfgRapidityMin || particle.y() - cfgRapidityCenterMass > cfgRapidityMax) { - continue; + return; } - candidate = fillCandidate(iSpecies, collision, track); - fillHistograms(iSpecies, candidate); + } + + mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts); + if (!trackSelection(track)) return; + mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kTrackCuts); + + if (!pidSelection(track, collision)) return; + mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kPidCuts); + + nuclei::SlimCandidate candidate; + if (track.has_mcParticle()) { + candidate = fillCandidate(iSpeciesCt, collision, track); + dispatchFillHistograms(iSpeciesRt, candidate); } else { - candidate = fillCandidate(iSpecies, collision, track); + candidate = fillCandidate(iSpeciesCt, collision, track); } - - fillHistograms(iSpecies, candidate); - } - } - for (const auto& particle: mcParticles) + dispatchFillHistograms(iSpeciesRt, candidate); + }); + } + + const int mcCollisionId = collision.mcCollisionId(); + auto mcParticlesThisCollision = mcParticles.sliceBy(mMcParticlesPerCollision, mcCollisionId); + mcParticlesThisCollision.bindExternalIndices(&mcParticles); + for (const auto& particle: mcParticlesThisCollision) { if (std::find(mFilledMcParticleIds.begin(), mFilledMcParticleIds.end(), particle.globalIndex()) != mFilledMcParticleIds.end()) { continue; } + int iSpecies = nuclei::getSpeciesFromPdg(particle.pdgCode()); + if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpecies) == mSpeciesToProcess.end()) { + continue; + } nuclei::SlimCandidate candidate; fillNucleusFlagsPdgsMc(particle, candidate); fillNucleusGeneratedVariables(particle, candidate); mNucleiCandidates.emplace_back(candidate); - const int iSpecies = nuclei::getSpeciesFromPdg(particle.pdgCode()); - fillHistograms(iSpecies, candidate); + dispatchFillHistograms(iSpecies, candidate); } - if (!cfgFillTable) { - return; - } + if (!cfgFillTable) return; for (const auto& candidate: mNucleiCandidates) { From f436d6ec8b1436dfd771ee0f964dd28a34ca1715 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 2 Oct 2025 16:08:55 +0000 Subject: [PATCH 4/9] Please consider the following formatting changes --- PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx | 42 ++- PWGLF/TableProducer/QC/nucleiQC.cxx | 225 +++++++------- PWGLF/Utils/nucleiUtils.h | 283 ++++++++++-------- 3 files changed, 285 insertions(+), 265 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx b/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx index c336bd8a74c..beb4da99806 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx @@ -19,54 +19,48 @@ // o2-analysis-pid-tof-base, o2-analysis-multiplicity-table, o2-analysis-event-selection // (to add flow: o2-analysis-qvector-table, o2-analysis-centrality-table) -#include -#include -#include -#include -#include - -#include "Math/Vector4D.h" - -#include "CCDB/BasicCCDBManager.h" +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "PWGLF/DataModel/LFSlimNucleiTables.h" +#include "PWGLF/Utils/nucleiUtils.h" +#include "Common/Core/EventPlaneHelper.h" +#include "Common/Core/PID/PIDTOF.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" -#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/Qvectors.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/PID/PIDTOF.h" #include "Common/TableProducer/PID/pidTOFBase.h" -#include "Common/Core/EventPlaneHelper.h" -#include "Common/DataModel/Qvectors.h" #include "Common/Tools/TrackTuner.h" -#include "Common/Core/RecoDecay.h" +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" +#include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPMagField.h" #include "DataFormatsParameters/GRPObject.h" #include "DataFormatsTPC/BetheBlochAleph.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" - -#include "EventFiltering/Zorro.h" -#include "EventFiltering/ZorroSummary.h" - +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" -#include "Framework/ASoAHelpers.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" - #include "ReconstructionDataFormats/Track.h" -#include "PWGLF/DataModel/EPCalibrationTables.h" -#include "PWGLF/DataModel/LFSlimNucleiTables.h" - +#include "Math/Vector4D.h" #include "TRandom3.h" -#include "PWGLF/Utils/nucleiUtils.h" +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; diff --git a/PWGLF/TableProducer/QC/nucleiQC.cxx b/PWGLF/TableProducer/QC/nucleiQC.cxx index afd71da2204..2672ef1029d 100644 --- a/PWGLF/TableProducer/QC/nucleiQC.cxx +++ b/PWGLF/TableProducer/QC/nucleiQC.cxx @@ -19,61 +19,54 @@ // o2-analysis-pid-tof-base, o2-analysis-multiplicity-table, o2-analysis-event-selection // (to add flow: o2-analysis-qvector-table, o2-analysis-centrality-table) -#include -#include -#include -#include -#include -#include - -#include "Math/Vector4D.h" - -#include "CCDB/BasicCCDBManager.h" +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "PWGLF/DataModel/LFSlimNucleiTables.h" +#include "PWGLF/Utils/nucleiUtils.h" +#include "Common/Core/EventPlaneHelper.h" +#include "Common/Core/PID/PIDTOF.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" -#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/Qvectors.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/PID/PIDTOF.h" #include "Common/TableProducer/PID/pidTOFBase.h" -#include "Common/Core/EventPlaneHelper.h" -#include "Common/DataModel/Qvectors.h" #include "Common/Tools/TrackTuner.h" -#include "Common/Core/RecoDecay.h" +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" +#include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPMagField.h" #include "DataFormatsParameters/GRPObject.h" #include "DataFormatsTPC/BetheBlochAleph.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" - -#include "EventFiltering/Zorro.h" -#include "EventFiltering/ZorroSummary.h" - +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" -#include "Framework/ASoAHelpers.h" #include "Framework/HistogramRegistry.h" -#include "Framework/runDataProcessing.h" #include "Framework/StaticFor.h" - +#include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/Track.h" -#include "PWGLF/DataModel/EPCalibrationTables.h" -#include "PWGLF/DataModel/LFSlimNucleiTables.h" - +#include "Math/Vector4D.h" #include "TRandom3.h" -#include "PWGLF/Utils/nucleiUtils.h" +#include +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; - struct nucleiQC { using TrackCandidates = soa::Join; @@ -109,7 +102,7 @@ struct nucleiQC { int mRunNumber = 0; float mBz = 0.f; - HistogramRegistry mHistograms { + HistogramRegistry mHistograms{ "histos", { {"hEventSelections", "Event selections; Selection step; Counts", {HistType::kTH1D, {{nuclei::evSel::kNevSels + 1, -0.5f, static_cast(nuclei::evSel::kNevSels) + 0.5f}}}}, @@ -118,19 +111,19 @@ struct nucleiQC { }, OutputObjHandlingPolicy::AnalysisObject, false, - true - }; + true}; std::vector mSpeciesToProcess; Produces mNucleiTableRed; - + std::vector mNucleiCandidates; std::vector mFilledMcParticleIds; - + o2::dataformats::DCA mDcaInfoCov; o2::track::TrackParametrizationWithError mTrackParCov; std::array(nuclei::Species::kNspecies)> mPidManagers; - void init(o2::framework::InitContext&) { + void init(o2::framework::InitContext&) + { mCcdb->setURL(cfgCCDBurl); mCcdb->setCaching(true); @@ -148,8 +141,7 @@ struct nucleiQC { } } - static_for<0, nuclei::kNspecies-1> ([&] (auto iSpecies) { - + static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) { constexpr int iSpeciesCt = decltype(iSpecies)::value; const int iSpeciesRt = iSpeciesCt; @@ -163,7 +155,7 @@ struct nucleiQC { } nuclei::createHistogramRegistryNucleus(mHistograms); - + if (cfgUseCentralTpcCalibration->get(uint32_t(iSpeciesRt), uint32_t(0)) == 0) { mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt, tpcBetheBlochParams); } else { @@ -188,8 +180,9 @@ struct nucleiQC { } template - bool trackSelection(const Ttrack& track) { - if (std::abs(track.eta()) > cfgCutEta || + bool trackSelection(const Ttrack& track) + { + if (std::abs(track.eta()) > cfgCutEta || track.tpcInnerParam() < cfgCutTpcMom || track.itsNCls() < cfgCutNclusITS || track.tpcNClsFound() < cfgCutNclusTPC || @@ -201,10 +194,10 @@ struct nucleiQC { } return true; } - + template - bool pidSelection(const Ttrack& track, const Tcollision& collision) - { + bool pidSelection(const Ttrack& track, const Tcollision& collision) + { constexpr int index = iSpecies; if (!nuclei::checkSpeciesValidity(index)) { std::runtime_error("species contains invalid nucleus index"); @@ -213,7 +206,8 @@ struct nucleiQC { const float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator); const float nsigmaTPC = mPidManagers[index].getNSigmaTPC(track); mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPC_preselectionVsCentrality"), track.pt(), nsigmaTPC, centrality); - if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(index, 1)) return false; + if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(index, 1)) + return false; mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPCVsCentrality"), track.pt(), nsigmaTPC, centrality); const float nsigmaITS = mPidManagers[index].getNSigmaITS(track); @@ -223,54 +217,55 @@ struct nucleiQC { const float nsigmaTOF = mPidManagers[index].getNSigmaTOF(track); mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOF_preselectionVsCentrality"), track.pt(), nsigmaTOF, centrality); - if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(index, 1)) return false; + if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(index, 1)) + return false; mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOFVsCentrality"), track.pt(), nsigmaTOF, centrality); return true; } template - void fillNucleusFlagsPdgsMc(const Tparticle& particle, nuclei::SlimCandidate& candidate) + void fillNucleusFlagsPdgsMc(const Tparticle& particle, nuclei::SlimCandidate& candidate) { candidate.pdgCode = particle.pdgCode(); if (particle.isPhysicalPrimary()) { - candidate.flags |= nuclei::Flags::kIsPhysicalPrimary; - - // heavy flavour mother - //if (particle.has_mothers()) { - // for (auto& motherparticle : particle.mothers_as()) { - // if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { - // flags |= kIsSecondaryFromWeakDecay; - // motherPdgCode = motherparticle.pdgCode(); - // break; - // } - // } - //} - - } else if (particle.has_mothers()) { - candidate.flags |= nuclei::Flags::kIsSecondaryFromWeakDecay; - for (auto& motherparticle : particle.template mothers_as()) { - candidate.motherPdgCode = motherparticle.pdgCode(); - } - - } else { - candidate.flags |= nuclei::Flags::kIsSecondaryFromMaterial; + candidate.flags |= nuclei::Flags::kIsPhysicalPrimary; + + // heavy flavour mother + // if (particle.has_mothers()) { + // for (auto& motherparticle : particle.mothers_as()) { + // if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { + // flags |= kIsSecondaryFromWeakDecay; + // motherPdgCode = motherparticle.pdgCode(); + // break; + // } + // } + //} + + } else if (particle.has_mothers()) { + candidate.flags |= nuclei::Flags::kIsSecondaryFromWeakDecay; + for (auto& motherparticle : particle.template mothers_as()) { + candidate.motherPdgCode = motherparticle.pdgCode(); } - mFilledMcParticleIds.emplace_back(particle.globalIndex()); + } else { + candidate.flags |= nuclei::Flags::kIsSecondaryFromMaterial; + } + + mFilledMcParticleIds.emplace_back(particle.globalIndex()); } template - void fillNucleusFlagsPdgs(const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) + void fillNucleusFlagsPdgs(const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) { candidate.flags = static_cast((track.pidForTracking() & 0xF) << 12); - candidate.flags |= iSpecies == nuclei::Species::kPr ? nuclei::Flags::kProton : - iSpecies == nuclei::Species::kDe ? nuclei::Flags::kDeuteron : - iSpecies == nuclei::Species::kTr ? nuclei::Flags::kTriton : - iSpecies == nuclei::Species::kHe ? nuclei::Flags::kHe3 : - iSpecies == nuclei::Species::kAl ? nuclei::Flags::kHe4 : 0; - + candidate.flags |= iSpecies == nuclei::Species::kPr ? nuclei::Flags::kProton : iSpecies == nuclei::Species::kDe ? nuclei::Flags::kDeuteron + : iSpecies == nuclei::Species::kTr ? nuclei::Flags::kTriton + : iSpecies == nuclei::Species::kHe ? nuclei::Flags::kHe3 + : iSpecies == nuclei::Species::kAl ? nuclei::Flags::kHe4 + : 0; + if (track.hasTOF()) { candidate.flags |= nuclei::Flags::kHasTOF; } @@ -291,7 +286,8 @@ struct nucleiQC { } template - void fillDcaInformation(const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) { + void fillDcaInformation(const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) + { const o2::math_utils::Point3D collisionVertex{collision.posX(), collision.posY(), collision.posZ()}; @@ -303,7 +299,6 @@ struct nucleiQC { candidate.DCAxy = dcaInfo[0]; candidate.DCAz = dcaInfo[1]; - } template @@ -327,35 +322,40 @@ struct nucleiQC { .ptGenerated = 0.f, // to be filled for mc .etaGenerated = 0.f, .phiGenerated = 0.f, - .centrality = nuclei::getCentrality(collision, cfgCentralityEstimator) - }; - + .centrality = nuclei::getCentrality(collision, cfgCentralityEstimator)}; + fillDcaInformation(collision, track, candidate); fillNucleusFlagsPdgs(iSpecies, collision, track, candidate); - + if constexpr (isMc) { if (track.has_mcParticle()) { - + const auto& particle = track.mcParticle(); fillNucleusFlagsPdgsMc(particle, candidate); fillNucleusGeneratedVariables(particle, candidate); } } - + mNucleiCandidates.emplace_back(candidate); return candidate; - } - template - void dispatchFillHistograms(const int iSpecies, const nuclei::SlimCandidate& candidate) { + template + void dispatchFillHistograms(const int iSpecies, const nuclei::SlimCandidate& candidate) + { switch (iSpecies) { - case nuclei::Species::kPr: return fillHistograms(candidate); - case nuclei::Species::kDe: return fillHistograms(candidate); - case nuclei::Species::kTr: return fillHistograms(candidate); - case nuclei::Species::kHe: return fillHistograms(candidate); - case nuclei::Species::kAl: return fillHistograms(candidate); - default: return; + case nuclei::Species::kPr: + return fillHistograms(candidate); + case nuclei::Species::kDe: + return fillHistograms(candidate); + case nuclei::Species::kTr: + return fillHistograms(candidate); + case nuclei::Species::kHe: + return fillHistograms(candidate); + case nuclei::Species::kAl: + return fillHistograms(candidate); + default: + return; } } @@ -366,7 +366,7 @@ struct nucleiQC { if (!nuclei::checkSpeciesValidity(index)) { std::runtime_error("species contains invalid nucleus index"); } - + if (isGenerated) { const float ptGenerated = (index == nuclei::Species::kPr || index == nuclei::Species::kDe || index == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f; mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/hPtGenerated"), ptGenerated); @@ -388,7 +388,7 @@ struct nucleiQC { { mNucleiCandidates.clear(); mFilledMcParticleIds.clear(); - + auto bc = collision.template bc_as(); initCCDB(bc); @@ -396,11 +396,9 @@ struct nucleiQC { return; } - for (const auto& track: tracks) - { - - static_for<0, nuclei::kNspecies-1> ([&] (auto iSpecies) { + for (const auto& track : tracks) { + static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) { constexpr int iSpeciesCt = decltype(iSpecies)::value; const int iSpeciesRt = iSpeciesCt; @@ -411,37 +409,39 @@ struct nucleiQC { if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); if (cfgDoCheckPdgCode) { - if (particle.pdgCode() != nuclei::pdgCodes[iSpeciesRt]) return; + if (particle.pdgCode() != nuclei::pdgCodes[iSpeciesRt]) + return; } if (particle.y() - cfgRapidityCenterMass < cfgRapidityMin || particle.y() - cfgRapidityCenterMass > cfgRapidityMax) { return; } - } + } mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts); - if (!trackSelection(track)) return; + if (!trackSelection(track)) + return; mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kTrackCuts); - if (!pidSelection(track, collision)) return; + if (!pidSelection(track, collision)) + return; mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kPidCuts); nuclei::SlimCandidate candidate; if (track.has_mcParticle()) { - candidate = fillCandidate(iSpeciesCt, collision, track); - dispatchFillHistograms(iSpeciesRt, candidate); + candidate = fillCandidate(iSpeciesCt, collision, track); + dispatchFillHistograms(iSpeciesRt, candidate); } else { - candidate = fillCandidate(iSpeciesCt, collision, track); + candidate = fillCandidate(iSpeciesCt, collision, track); } - dispatchFillHistograms(iSpeciesRt, candidate); + dispatchFillHistograms(iSpeciesRt, candidate); }); } - + const int mcCollisionId = collision.mcCollisionId(); auto mcParticlesThisCollision = mcParticles.sliceBy(mMcParticlesPerCollision, mcCollisionId); mcParticlesThisCollision.bindExternalIndices(&mcParticles); - for (const auto& particle: mcParticlesThisCollision) - { + for (const auto& particle : mcParticlesThisCollision) { if (std::find(mFilledMcParticleIds.begin(), mFilledMcParticleIds.end(), particle.globalIndex()) != mFilledMcParticleIds.end()) { continue; } @@ -455,13 +455,13 @@ struct nucleiQC { fillNucleusGeneratedVariables(particle, candidate); mNucleiCandidates.emplace_back(candidate); - dispatchFillHistograms(iSpecies, candidate); + dispatchFillHistograms(iSpecies, candidate); } - if (!cfgFillTable) return; + if (!cfgFillTable) + return; - for (const auto& candidate: mNucleiCandidates) - { + for (const auto& candidate : mNucleiCandidates) { mNucleiTableRed( candidate.pt, candidate.eta, @@ -474,8 +474,7 @@ struct nucleiQC { candidate.DCAz, candidate.flags, candidate.pdgCode, - candidate.motherPdgCode - ); + candidate.motherPdgCode); } } PROCESS_SWITCH(nucleiQC, processMc, "Mc analysis", false); diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h index eeb70aa51b2..e954ef76c1e 100644 --- a/PWGLF/Utils/nucleiUtils.h +++ b/PWGLF/Utils/nucleiUtils.h @@ -12,12 +12,12 @@ #ifndef PWGLF_UTILS_NUCLEIUTILS_H_ #define PWGLF_UTILS_NUCLEIUTILS_H_ -#include -#include - #include "Framework/HistogramRegistry.h" #include "Framework/HistogramSpec.h" +#include +#include + using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -119,18 +119,26 @@ enum Flags { kIsSecondaryFromWeakDecay = BIT(11) /// the last 4 bits are reserved for the PID in tracking }; -constexpr int getSpeciesFromPdg(int pdg) { +constexpr int getSpeciesFromPdg(int pdg) +{ switch (std::abs(pdg)) { - case PDG_t::kProton: return Species::kPr; - case o2::constants::physics::Pdg::kDeuteron: return Species::kDe; - case o2::constants::physics::Pdg::kTriton: return Species::kTr; - case o2::constants::physics::Pdg::kHelium3: return Species::kHe; - case o2::constants::physics::Pdg::kAlpha: return Species::kAl; - default: return -1; + case PDG_t::kProton: + return Species::kPr; + case o2::constants::physics::Pdg::kDeuteron: + return Species::kDe; + case o2::constants::physics::Pdg::kTriton: + return Species::kTr; + case o2::constants::physics::Pdg::kHelium3: + return Species::kHe; + case o2::constants::physics::Pdg::kAlpha: + return Species::kAl; + default: + return -1; } } -bool checkSpeciesValidity(const int species) { +bool checkSpeciesValidity(const int species) +{ if (species < 0 || species > Species::kNspecies) { return false; } @@ -356,8 +364,9 @@ enum trackSelection { }; static const std::array(trackSelection::kNtrackSelections)> trackSelectionLabels{"All", "Track cuts", "PID cuts"}; -template -void createHistogramRegistryNucleus(HistogramRegistry& registry) { +template +void createHistogramRegistryNucleus(HistogramRegistry& registry) +{ constexpr int index = iSpecies; if (!checkSpeciesValidity(index)) { @@ -369,7 +378,7 @@ void createHistogramRegistryNucleus(HistogramRegistry& registry) { registry.add(fmt::format("{}/h3PtVsEtaVsCentralityReconstructed", cNames[index]).c_str(), (fmt::format("{} - reconstructed variables;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3PhiVsEtaVsCentralityReconstructed", cNames[index]).c_str(), (fmt::format("{} - reconstructed variables;", cNames[index]) + std::string("#phi (radians); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3DCAxyVsPtVsCentrality", cNames[index]).c_str(), (fmt::format(";", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{xy} (cm); CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); - registry.add(fmt::format("{}/h3DCAzVsPtVsCentrality", cNames[index]).c_str(), (fmt::format("{};", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{z} (cm); CentralityFT0C (%)")).c_str() , HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3DCAzVsPtVsCentrality", cNames[index]).c_str(), (fmt::format("{};", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{z} (cm); CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3NsigmaTPC_preselectionVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TPC distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TPC}}({}); CentralityFT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3NsigmaTPCVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TPC distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TPC}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3NsigmaITS_preselectionVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} ITS distribution;", cNames[index]) + std::string("signed #it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{ITS}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}); @@ -390,141 +399,159 @@ void createHistogramRegistryNucleus(HistogramRegistry& registry) { // PID manager class -class PidManager { - - public: - PidManager(const int species, const float* tpcBetheBlochParams = nullptr) +class PidManager +{ + + public: + PidManager(const int species, const float* tpcBetheBlochParams = nullptr) : mSpecies(species) - { - if (!checkSpeciesValidity(species)) { - std::runtime_error("species contains invalid nucleus index"); - } - - if (!tpcBetheBlochParams) { - mUseTpcCentralCalibration = true; - return; - } - - for (int i = 0; i < 6; i++) { - mTpcBetheBlochParams[i] = tpcBetheBlochParams[i]; - } + { + if (!checkSpeciesValidity(species)) { + std::runtime_error("species contains invalid nucleus index"); } - PidManager() = default; - ~PidManager() = default; - - // TOF - template - float getBetaTOF(const Ttrack& track) - { - if (!track.hasTOF()) { - return -999.f; - } - float beta = o2::pid::tof::Beta::GetBeta(track); - return std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + if (!tpcBetheBlochParams) { + mUseTpcCentralCalibration = true; + return; } - template - float getMassTOF(const Ttrack& track) - { - if (!track.hasTOF()) { - return -999.f; - } - const float charge{1.f + static_cast(mSpecies == Species::kHe || mSpecies == Species::kAl)}; - const float beta = getBetaTOF(track); - return track.tpcInnerParam() * charge * std::sqrt(1.f / (beta * beta) - 1.f); + for (int i = 0; i < 6; i++) { + mTpcBetheBlochParams[i] = tpcBetheBlochParams[i]; } - - template - float getNSigmaTOF(const Ttrack& track) - { - if (!track.hasTOF()) { - return -999.f; - } - - switch (mSpecies) { - case Species::kPr: return track.tofNSigmaPr(); - case Species::kDe: return track.tofNSigmaDe(); - case Species::kTr: return track.tofNSigmaTr(); - case Species::kHe: return track.tofNSigmaHe(); - case Species::kAl: return track.tofNSigmaAl(); - default: return -999.f; - } + } + PidManager() = default; + ~PidManager() = default; + + // TOF + template + float getBetaTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; } + float beta = o2::pid::tof::Beta::GetBeta(track); + return std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + } - template - float getClusterSizeCosLambdaITS(const Ttrack& track) - { - return mResponseITS.averageClusterSize(track.itsClusterSizes()) / std::cosh(track.eta()); + template + float getMassTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; } + const float charge{1.f + static_cast(mSpecies == Species::kHe || mSpecies == Species::kAl)}; + const float beta = getBetaTOF(track); + return track.tpcInnerParam() * charge * std::sqrt(1.f / (beta * beta) - 1.f); + } - float getClusterSizeCosLambdaITS(const u_int32_t clusterSizesITS, const float eta) - { - return mResponseITS.averageClusterSize(clusterSizesITS) / std::cosh(eta); + template + float getNSigmaTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; } - template - float getNSigmaITS(const Ttrack& track) - { - switch (mSpecies) { - case Species::kPr: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); - case Species::kDe: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); - case Species::kTr: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); - case Species::kHe: return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); - case Species::kAl: return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); - default: return -999.f; - } + switch (mSpecies) { + case Species::kPr: + return track.tofNSigmaPr(); + case Species::kDe: + return track.tofNSigmaDe(); + case Species::kTr: + return track.tofNSigmaTr(); + case Species::kHe: + return track.tofNSigmaHe(); + case Species::kAl: + return track.tofNSigmaAl(); + default: + return -999.f; } + } - // TPC - float getExpectedTPCsignal(const float p) - { - if (!mUseTpcCentralCalibration) { + template + float getClusterSizeCosLambdaITS(const Ttrack& track) + { + return mResponseITS.averageClusterSize(track.itsClusterSizes()) / std::cosh(track.eta()); + } + + float getClusterSizeCosLambdaITS(const u_int32_t clusterSizesITS, const float eta) + { + return mResponseITS.averageClusterSize(clusterSizesITS) / std::cosh(eta); + } + + template + float getNSigmaITS(const Ttrack& track) + { + switch (mSpecies) { + case Species::kPr: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kDe: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kTr: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kHe: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); + case Species::kAl: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); + default: return -999.f; - } - float pScaled = p * mMomScaling[0] + mMomScaling[1]; - float betaGamma = pScaled / masses[mSpecies]; - return tpc::BetheBlochAleph(betaGamma, - mTpcBetheBlochParams[0], - mTpcBetheBlochParams[1], - mTpcBetheBlochParams[2], - mTpcBetheBlochParams[3], - mTpcBetheBlochParams[4]); } - - template - float getNSigmaTPC(const Ttrack& track) - { - if (!mUseTpcCentralCalibration) { - return getNSigmaTPCcentral(track); - } - float expectedSignal = getExpectedTPCsignal(track.tpcInnerParam()); - float resolution = mTpcBetheBlochParams[5]; - return (track.tpcSignal() - expectedSignal) / (expectedSignal * resolution); + } + + // TPC + float getExpectedTPCsignal(const float p) + { + if (!mUseTpcCentralCalibration) { + return -999.f; } + float pScaled = p * mMomScaling[0] + mMomScaling[1]; + float betaGamma = pScaled / masses[mSpecies]; + return tpc::BetheBlochAleph(betaGamma, + mTpcBetheBlochParams[0], + mTpcBetheBlochParams[1], + mTpcBetheBlochParams[2], + mTpcBetheBlochParams[3], + mTpcBetheBlochParams[4]); + } - protected: - // TPC - template - float getNSigmaTPCcentral(const Ttrack& track) { - switch (mSpecies) { - case Species::kPr: return track.tpcNSigmaPr(); - case Species::kDe: return track.tpcNSigmaDe(); - case Species::kTr: return track.tpcNSigmaTr(); - case Species::kHe: return track.tpcNSigmaHe(); - case Species::kAl: return track.tpcNSigmaAl(); - default: return -999.f; - } + template + float getNSigmaTPC(const Ttrack& track) + { + if (!mUseTpcCentralCalibration) { + return getNSigmaTPCcentral(track); } + float expectedSignal = getExpectedTPCsignal(track.tpcInnerParam()); + float resolution = mTpcBetheBlochParams[5]; + return (track.tpcSignal() - expectedSignal) / (expectedSignal * resolution); + } - private: - float mTpcBetheBlochParams[6]; - bool mUseTpcCentralCalibration = true; // this just becomes a check for the null pointer in the parameters - o2::aod::ITSResponse mResponseITS; - float mMomScaling[2]{1., 0.}; - int mSpecies; -}; + protected: + // TPC + template + float getNSigmaTPCcentral(const Ttrack& track) + { + switch (mSpecies) { + case Species::kPr: + return track.tpcNSigmaPr(); + case Species::kDe: + return track.tpcNSigmaDe(); + case Species::kTr: + return track.tpcNSigmaTr(); + case Species::kHe: + return track.tpcNSigmaHe(); + case Species::kAl: + return track.tpcNSigmaAl(); + default: + return -999.f; + } + } + private: + float mTpcBetheBlochParams[6]; + bool mUseTpcCentralCalibration = true; // this just becomes a check for the null pointer in the parameters + o2::aod::ITSResponse mResponseITS; + float mMomScaling[2]{1., 0.}; + int mSpecies; +}; } // namespace nuclei From 1f13e261c3d6a52ece4fac14158ad53df6583c6b Mon Sep 17 00:00:00 2001 From: Giorgio Alberto Lucia <87222843+GiorgioAlbertoLucia@users.noreply.github.com> Date: Thu, 2 Oct 2025 18:18:46 +0200 Subject: [PATCH 5/9] MegaLinter fix --- PWGLF/TableProducer/QC/nucleiQC.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/QC/nucleiQC.cxx b/PWGLF/TableProducer/QC/nucleiQC.cxx index 2672ef1029d..ee1e0f9c8ab 100644 --- a/PWGLF/TableProducer/QC/nucleiQC.cxx +++ b/PWGLF/TableProducer/QC/nucleiQC.cxx @@ -156,7 +156,7 @@ struct nucleiQC { nuclei::createHistogramRegistryNucleus(mHistograms); - if (cfgUseCentralTpcCalibration->get(uint32_t(iSpeciesRt), uint32_t(0)) == 0) { + if (cfgUseCentralTpcCalibration->get(static_cast(iSpeciesRt), static_cast(0)) == 0) { mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt, tpcBetheBlochParams); } else { mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt); From e530ae827c5f4273934a668d56e9cb32d11da845 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 3 Oct 2025 10:01:10 +0200 Subject: [PATCH 6/9] fixed build errors spotted during the PR --- PWGLF/TableProducer/QC/nucleiQC.cxx | 78 ++++++++++++++--------------- PWGLF/Utils/nucleiUtils.h | 64 ++++++++++------------- 2 files changed, 66 insertions(+), 76 deletions(-) diff --git a/PWGLF/TableProducer/QC/nucleiQC.cxx b/PWGLF/TableProducer/QC/nucleiQC.cxx index ee1e0f9c8ab..30323d4f6ab 100644 --- a/PWGLF/TableProducer/QC/nucleiQC.cxx +++ b/PWGLF/TableProducer/QC/nucleiQC.cxx @@ -142,24 +142,24 @@ struct nucleiQC { } static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) { - constexpr int iSpeciesCt = decltype(iSpecies)::value; - const int iSpeciesRt = iSpeciesCt; + constexpr int kSpeciesCt = decltype(iSpecies)::value; + const int kSpeciesRt = kSpeciesCt; - if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpeciesCt) == mSpeciesToProcess.end()) { + if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesCt) == mSpeciesToProcess.end()) { return; } float tpcBetheBlochParams[6]; for (int iParam = 0; iParam < 6; iParam++) { - tpcBetheBlochParams[iParam] = cfgBetheBlochParams->get(iSpeciesRt, iParam); + tpcBetheBlochParams[iParam] = cfgBetheBlochParams->get(kSpeciesRt, iParam); } nuclei::createHistogramRegistryNucleus(mHistograms); - if (cfgUseCentralTpcCalibration->get(static_cast(iSpeciesRt), static_cast(0)) == 0) { - mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt, tpcBetheBlochParams); + if (cfgUseCentralTpcCalibration->get(static_cast(kSpeciesRt), static_cast(0)) == 0) { + mPidManagers[kSpeciesRt] = nuclei::PidManager(kSpeciesRt, tpcBetheBlochParams); } else { - mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt); + mPidManagers[kSpeciesRt] = nuclei::PidManager(kSpeciesRt); } }); } @@ -198,28 +198,28 @@ struct nucleiQC { template bool pidSelection(const Ttrack& track, const Tcollision& collision) { - constexpr int index = iSpecies; - if (!nuclei::checkSpeciesValidity(index)) { - std::runtime_error("species contains invalid nucleus index"); + constexpr int kIndex = iSpecies; + if (!nuclei::checkSpeciesValidity(kIndex)) { + std::runtime_error("species contains invalid nucleus kIndex"); } - const float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator); - const float nsigmaTPC = mPidManagers[index].getNSigmaTPC(track); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPC_preselectionVsCentrality"), track.pt(), nsigmaTPC, centrality); - if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(index, 1)) + float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator); + float nsigmaTPC = mPidManagers[kIndex].getNSigmaTPC(track); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3NsigmaTPC_preselectionVsCentrality"), track.pt(), nsigmaTPC, centrality); + if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(kIndex, 1)) return false; - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPCVsCentrality"), track.pt(), nsigmaTPC, centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3NsigmaTPCVsCentrality"), track.pt(), nsigmaTPC, centrality); - const float nsigmaITS = mPidManagers[index].getNSigmaITS(track); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaITS_preselectionVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); + float nsigmaITS = mPidManagers[kIndex].getNSigmaITS(track); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3NsigmaITS_preselectionVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); // add nsigmaITS cut ? - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaITSVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3NsigmaITSVsCentrality"), track.sign() * track.pt(), nsigmaITS, centrality); - const float nsigmaTOF = mPidManagers[index].getNSigmaTOF(track); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOF_preselectionVsCentrality"), track.pt(), nsigmaTOF, centrality); - if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(index, 1)) + float nsigmaTOF = mPidManagers[kIndex].getNSigmaTOF(track); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3NsigmaTOF_preselectionVsCentrality"), track.pt(), nsigmaTOF, centrality); + if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(kIndex, 1)) return false; - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOFVsCentrality"), track.pt(), nsigmaTOF, centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3NsigmaTOFVsCentrality"), track.pt(), nsigmaTOF, centrality); return true; } @@ -234,7 +234,7 @@ struct nucleiQC { // heavy flavour mother // if (particle.has_mothers()) { - // for (auto& motherparticle : particle.mothers_as()) { + // for (const auto& motherparticle : particle.mothers_as()) { // if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { // flags |= kIsSecondaryFromWeakDecay; // motherPdgCode = motherparticle.pdgCode(); @@ -245,7 +245,7 @@ struct nucleiQC { } else if (particle.has_mothers()) { candidate.flags |= nuclei::Flags::kIsSecondaryFromWeakDecay; - for (auto& motherparticle : particle.template mothers_as()) { + for (const auto& motherparticle : particle.template mothers_as()) { candidate.motherPdgCode = motherparticle.pdgCode(); } @@ -362,25 +362,25 @@ struct nucleiQC { template void fillHistograms(const nuclei::SlimCandidate& candidate) { - constexpr int index = iSpecies; - if (!nuclei::checkSpeciesValidity(index)) { - std::runtime_error("species contains invalid nucleus index"); + constexpr int kIndex = iSpecies; + if (!nuclei::checkSpeciesValidity(kIndex)) { + std::runtime_error("species contains invalid nucleus kIndex"); } if (isGenerated) { - const float ptGenerated = (index == nuclei::Species::kPr || index == nuclei::Species::kDe || index == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f; - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/hPtGenerated"), ptGenerated); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PtVsEtaVsCentralityGenerated"), ptGenerated, candidate.etaGenerated, candidate.centrality); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PhiVsEtaVsCentralityGenerated"), candidate.phiGenerated, candidate.etaGenerated, candidate.centrality); + const float ptGenerated = (kIndex == nuclei::Species::kPr || kIndex == nuclei::Species::kDe || kIndex == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f; + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/hPtGenerated"), ptGenerated); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3PtVsEtaVsCentralityGenerated"), ptGenerated, candidate.etaGenerated, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3PhiVsEtaVsCentralityGenerated"), candidate.phiGenerated, candidate.etaGenerated, candidate.centrality); } else { - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/hPtReconstructed"), candidate.pt); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PtVsEtaVsCentralityReconstructed"), candidate.pt, candidate.eta, candidate.centrality); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3PhiVsEtaVsCentralityReconstructed"), candidate.phi, candidate.eta, candidate.centrality); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3DCAxyVsPtVsCentrality"), candidate.DCAxy, candidate.pt, candidate.centrality); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3DCAzVsPtVsCentrality"), candidate.DCAz, candidate.pt, candidate.centrality); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3BetaVsPtVsCentrality"), candidate.beta, candidate.pt, candidate.centrality); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3dEdxVsPVsCentrality"), candidate.TPCsignal, candidate.pt, candidate.centrality); - mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3ClusterSizeVsPtVsCentrality"), mPidManagers[index].getClusterSizeCosLambdaITS(candidate.clusterSizesITS, candidate.eta), candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/hPtReconstructed"), candidate.pt); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3PtVsEtaVsCentralityReconstructed"), candidate.pt, candidate.eta, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3PhiVsEtaVsCentralityReconstructed"), candidate.phi, candidate.eta, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3DCAxyVsPtVsCentrality"), candidate.DCAxy, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3DCAzVsPtVsCentrality"), candidate.DCAz, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3BetaVsPtVsCentrality"), candidate.beta, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3dEdxVsPVsCentrality"), candidate.TPCsignal, candidate.pt, candidate.centrality); + mHistograms.fill(HIST(nuclei::cNames[kIndex]) + HIST("/h3ClusterSizeVsPtVsCentrality"), mPidManagers[kIndex].getClusterSizeCosLambdaITS(candidate.clusterSizesITS, candidate.eta), candidate.pt, candidate.centrality); } } diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h index e954ef76c1e..10720befcb4 100644 --- a/PWGLF/Utils/nucleiUtils.h +++ b/PWGLF/Utils/nucleiUtils.h @@ -15,8 +15,17 @@ #include "Framework/HistogramRegistry.h" #include "Framework/HistogramSpec.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/PIDResponseITS.h" +#include "Common/TableProducer/PID/pidTOFBase.h" + +#include "DataFormatsTPC/BetheBlochAleph.h" + #include #include +#include using namespace o2; using namespace o2::framework; @@ -403,7 +412,7 @@ class PidManager { public: - PidManager(const int species, const float* tpcBetheBlochParams = nullptr) + explicit PidManager(const int species, const float* tpcBetheBlochParams = nullptr) : mSpecies(species) { if (!checkSpeciesValidity(species)) { @@ -426,9 +435,7 @@ class PidManager template float getBetaTOF(const Ttrack& track) { - if (!track.hasTOF()) { - return -999.f; - } + if (!track.hasTOF()) return -999.f; float beta = o2::pid::tof::Beta::GetBeta(track); return std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked } @@ -436,9 +443,7 @@ class PidManager template float getMassTOF(const Ttrack& track) { - if (!track.hasTOF()) { - return -999.f; - } + if (!track.hasTOF()) return -999.f; const float charge{1.f + static_cast(mSpecies == Species::kHe || mSpecies == Species::kAl)}; const float beta = getBetaTOF(track); return track.tpcInnerParam() * charge * std::sqrt(1.f / (beta * beta) - 1.f); @@ -447,23 +452,15 @@ class PidManager template float getNSigmaTOF(const Ttrack& track) { - if (!track.hasTOF()) { - return -999.f; - } + if (!track.hasTOF()) return -999.f; switch (mSpecies) { - case Species::kPr: - return track.tofNSigmaPr(); - case Species::kDe: - return track.tofNSigmaDe(); - case Species::kTr: - return track.tofNSigmaTr(); - case Species::kHe: - return track.tofNSigmaHe(); - case Species::kAl: - return track.tofNSigmaAl(); - default: - return -999.f; + case Species::kPr: return track.tofNSigmaPr(); + case Species::kDe: return track.tofNSigmaDe(); + case Species::kTr: return track.tofNSigmaTr(); + case Species::kHe: return track.tofNSigmaHe(); + case Species::kAl: return track.tofNSigmaAl(); + default: return -999.f; } } @@ -500,9 +497,8 @@ class PidManager // TPC float getExpectedTPCsignal(const float p) { - if (!mUseTpcCentralCalibration) { - return -999.f; - } + if (!mUseTpcCentralCalibration) return -999.f; + float pScaled = p * mMomScaling[0] + mMomScaling[1]; float betaGamma = pScaled / masses[mSpecies]; return tpc::BetheBlochAleph(betaGamma, @@ -530,18 +526,12 @@ class PidManager float getNSigmaTPCcentral(const Ttrack& track) { switch (mSpecies) { - case Species::kPr: - return track.tpcNSigmaPr(); - case Species::kDe: - return track.tpcNSigmaDe(); - case Species::kTr: - return track.tpcNSigmaTr(); - case Species::kHe: - return track.tpcNSigmaHe(); - case Species::kAl: - return track.tpcNSigmaAl(); - default: - return -999.f; + case Species::kPr: return track.tpcNSigmaPr(); + case Species::kDe: return track.tpcNSigmaDe(); + case Species::kTr: return track.tpcNSigmaTr(); + case Species::kHe: return track.tpcNSigmaHe(); + case Species::kAl: return track.tpcNSigmaAl(); + default: return -999.f; } } From 3a5b78c0ba5dc0959760c664416972f483adb4a2 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 3 Oct 2025 08:03:13 +0000 Subject: [PATCH 7/9] Please consider the following formatting changes --- PWGLF/Utils/nucleiUtils.h | 55 +++++++++++++++++++++++++-------------- 1 file changed, 35 insertions(+), 20 deletions(-) diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h index 10720befcb4..f75dab12a74 100644 --- a/PWGLF/Utils/nucleiUtils.h +++ b/PWGLF/Utils/nucleiUtils.h @@ -12,9 +12,6 @@ #ifndef PWGLF_UTILS_NUCLEIUTILS_H_ #define PWGLF_UTILS_NUCLEIUTILS_H_ -#include "Framework/HistogramRegistry.h" -#include "Framework/HistogramSpec.h" - #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" @@ -22,10 +19,12 @@ #include "Common/TableProducer/PID/pidTOFBase.h" #include "DataFormatsTPC/BetheBlochAleph.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/HistogramSpec.h" +#include #include #include -#include using namespace o2; using namespace o2::framework; @@ -435,7 +434,8 @@ class PidManager template float getBetaTOF(const Ttrack& track) { - if (!track.hasTOF()) return -999.f; + if (!track.hasTOF()) + return -999.f; float beta = o2::pid::tof::Beta::GetBeta(track); return std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked } @@ -443,7 +443,8 @@ class PidManager template float getMassTOF(const Ttrack& track) { - if (!track.hasTOF()) return -999.f; + if (!track.hasTOF()) + return -999.f; const float charge{1.f + static_cast(mSpecies == Species::kHe || mSpecies == Species::kAl)}; const float beta = getBetaTOF(track); return track.tpcInnerParam() * charge * std::sqrt(1.f / (beta * beta) - 1.f); @@ -452,15 +453,22 @@ class PidManager template float getNSigmaTOF(const Ttrack& track) { - if (!track.hasTOF()) return -999.f; + if (!track.hasTOF()) + return -999.f; switch (mSpecies) { - case Species::kPr: return track.tofNSigmaPr(); - case Species::kDe: return track.tofNSigmaDe(); - case Species::kTr: return track.tofNSigmaTr(); - case Species::kHe: return track.tofNSigmaHe(); - case Species::kAl: return track.tofNSigmaAl(); - default: return -999.f; + case Species::kPr: + return track.tofNSigmaPr(); + case Species::kDe: + return track.tofNSigmaDe(); + case Species::kTr: + return track.tofNSigmaTr(); + case Species::kHe: + return track.tofNSigmaHe(); + case Species::kAl: + return track.tofNSigmaAl(); + default: + return -999.f; } } @@ -497,7 +505,8 @@ class PidManager // TPC float getExpectedTPCsignal(const float p) { - if (!mUseTpcCentralCalibration) return -999.f; + if (!mUseTpcCentralCalibration) + return -999.f; float pScaled = p * mMomScaling[0] + mMomScaling[1]; float betaGamma = pScaled / masses[mSpecies]; @@ -526,12 +535,18 @@ class PidManager float getNSigmaTPCcentral(const Ttrack& track) { switch (mSpecies) { - case Species::kPr: return track.tpcNSigmaPr(); - case Species::kDe: return track.tpcNSigmaDe(); - case Species::kTr: return track.tpcNSigmaTr(); - case Species::kHe: return track.tpcNSigmaHe(); - case Species::kAl: return track.tpcNSigmaAl(); - default: return -999.f; + case Species::kPr: + return track.tpcNSigmaPr(); + case Species::kDe: + return track.tpcNSigmaDe(); + case Species::kTr: + return track.tpcNSigmaTr(); + case Species::kHe: + return track.tpcNSigmaHe(); + case Species::kAl: + return track.tpcNSigmaAl(); + default: + return -999.f; } } From c0ee2c452d6e4914438c620cc11d5ab1c9989af3 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 3 Oct 2025 10:46:11 +0200 Subject: [PATCH 8/9] fixed compilation errors spotted during PR --- PWGLF/TableProducer/QC/nucleiQC.cxx | 26 +++++++++++++------------- PWGLF/Utils/nucleiUtils.h | 2 ++ 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/PWGLF/TableProducer/QC/nucleiQC.cxx b/PWGLF/TableProducer/QC/nucleiQC.cxx index 30323d4f6ab..739890077b7 100644 --- a/PWGLF/TableProducer/QC/nucleiQC.cxx +++ b/PWGLF/TableProducer/QC/nucleiQC.cxx @@ -154,7 +154,7 @@ struct nucleiQC { tpcBetheBlochParams[iParam] = cfgBetheBlochParams->get(kSpeciesRt, iParam); } - nuclei::createHistogramRegistryNucleus(mHistograms); + nuclei::createHistogramRegistryNucleus(mHistograms); if (cfgUseCentralTpcCalibration->get(static_cast(kSpeciesRt), static_cast(0)) == 0) { mPidManagers[kSpeciesRt] = nuclei::PidManager(kSpeciesRt, tpcBetheBlochParams); @@ -399,17 +399,17 @@ struct nucleiQC { for (const auto& track : tracks) { static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) { - constexpr int iSpeciesCt = decltype(iSpecies)::value; - const int iSpeciesRt = iSpeciesCt; + constexpr int kSpeciesCt = decltype(iSpecies)::value; + const int kSpeciesRt = kSpeciesCt; - if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpeciesRt) == mSpeciesToProcess.end()) { + if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesRt) == mSpeciesToProcess.end()) { return; } if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); if (cfgDoCheckPdgCode) { - if (particle.pdgCode() != nuclei::pdgCodes[iSpeciesRt]) + if (particle.pdgCode() != nuclei::pdgCodes[kSpeciesRt]) return; } if (particle.y() - cfgRapidityCenterMass < cfgRapidityMin || particle.y() - cfgRapidityCenterMass > cfgRapidityMax) { @@ -417,24 +417,24 @@ struct nucleiQC { } } - mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts); + mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts); if (!trackSelection(track)) return; - mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kTrackCuts); + mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kTrackCuts); - if (!pidSelection(track, collision)) + if (!pidSelection(track, collision)) return; - mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kPidCuts); + mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kPidCuts); nuclei::SlimCandidate candidate; if (track.has_mcParticle()) { - candidate = fillCandidate(iSpeciesCt, collision, track); - dispatchFillHistograms(iSpeciesRt, candidate); + candidate = fillCandidate(kSpeciesCt, collision, track); + dispatchFillHistograms(kSpeciesRt, candidate); } else { - candidate = fillCandidate(iSpeciesCt, collision, track); + candidate = fillCandidate(kSpeciesCt, collision, track); } - dispatchFillHistograms(iSpeciesRt, candidate); + dispatchFillHistograms(kSpeciesRt, candidate); }); } diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h index f75dab12a74..01aec18c286 100644 --- a/PWGLF/Utils/nucleiUtils.h +++ b/PWGLF/Utils/nucleiUtils.h @@ -21,6 +21,8 @@ #include "DataFormatsTPC/BetheBlochAleph.h" #include "Framework/HistogramRegistry.h" #include "Framework/HistogramSpec.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" #include #include From 60d9f7befee3d11700bcafac667433b4f8c96c60 Mon Sep 17 00:00:00 2001 From: GiorgioAlbertoLucia Date: Fri, 3 Oct 2025 18:53:47 +0200 Subject: [PATCH 9/9] clang format --- PWGLF/Utils/nucleiUtils.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h index 01aec18c286..967f017cbe7 100644 --- a/PWGLF/Utils/nucleiUtils.h +++ b/PWGLF/Utils/nucleiUtils.h @@ -19,10 +19,10 @@ #include "Common/TableProducer/PID/pidTOFBase.h" #include "DataFormatsTPC/BetheBlochAleph.h" -#include "Framework/HistogramRegistry.h" -#include "Framework/HistogramSpec.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/HistogramSpec.h" #include #include