From a384dd811b2e998b320bd49a3575a10934c066fa Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Fri, 21 Apr 2023 14:12:18 +0200 Subject: [PATCH 1/5] Remove dependency from the track propagation --- PWGLF/Tasks/CMakeLists.txt | 3 +- PWGLF/Tasks/NucleiSpectraTask.cxx | 132 ++++++++++++++++++++++++------ 2 files changed, 107 insertions(+), 28 deletions(-) diff --git a/PWGLF/Tasks/CMakeLists.txt b/PWGLF/Tasks/CMakeLists.txt index 3fe1ec0c692..2be2d30ef53 100644 --- a/PWGLF/Tasks/CMakeLists.txt +++ b/PWGLF/Tasks/CMakeLists.txt @@ -63,7 +63,7 @@ o2physics_add_dpl_workflow(spectra-tpc-tiny-pikapr o2physics_add_dpl_workflow(nuclei-spectra SOURCES NucleiSpectraTask.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(nuclei-efficiency @@ -160,4 +160,3 @@ o2physics_add_dpl_workflow(stqa SOURCES stqa.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::ReconstructionDataFormats O2Physics::AnalysisCore COMPONENT_NAME Analysis) - diff --git a/PWGLF/Tasks/NucleiSpectraTask.cxx b/PWGLF/Tasks/NucleiSpectraTask.cxx index 03dda0ef262..1cb5b11678d 100644 --- a/PWGLF/Tasks/NucleiSpectraTask.cxx +++ b/PWGLF/Tasks/NucleiSpectraTask.cxx @@ -15,15 +15,19 @@ // Executable + dependencies: // // Data (run3): -// o2-analysis-lf-nuclei-spectra, o2-analysis-track-propagation, o2-analysis-timestamp -// o2-analysis-trackselection, o2-analysis-pid-tof-base, o2-analysis-pid-tof-full -// o2-analysis-pid-tpc-base -// o2-analysis-pid-tpc-full, o2-analysis-multiplicity-table, o2-analysis-event-selection +// o2-analysis-lf-nuclei-spectra, o2-analysis-timestamp +// o2-analysis-pid-tof-base, o2-analysis-pid-tof-full +// o2-analysis-pid-tpc-base, o2-analysis-pid-tpc-full, +// o2-analysis-multiplicity-table, o2-analysis-event-selection #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/EventSelection.h" #include "Common/DataModel/PIDResponse.h" @@ -31,7 +35,11 @@ #include "Common/Core/PID/PIDTOF.h" #include "Common/TableProducer/PID/pidTOFBase.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" #include "DataFormatsTPC/BetheBlochAleph.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" @@ -213,6 +221,7 @@ struct NucleiSpectraTask { Produces nucleiTable; Produces nucleiTableMC; + Service ccdb; Configurable cfgCentralityEstimator{"cfgCentralityEstimator", "V0A", "Centrality estimator name"}; Configurable cfgCMrapidity{"cfgCMrapidity", 0.f, "Rapidity of the center of mass (only for p-Pb)"}; @@ -220,7 +229,7 @@ struct NucleiSpectraTask { Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; Configurable cfgCutRapidityMin{"cfgCutRapidityMin", -0.5, "Minimum rapidity for tracks"}; Configurable cfgCutRapidityMax{"cfgCutRapidityMax", 0.5, "Maximum rapidity for tracks"}; - Configurable cfgCutNclusITS{"cfgCutNclusITS", 4, "Minimum number of ITS clusters"}; + Configurable cfgCutNclusITS{"cfgCutNclusITS", 5, "Minimum number of ITS clusters"}; Configurable cfgCutNclusTPC{"cfgCutNclusTPC", 70, "Minimum number of TPC clusters"}; Configurable> cfgMomentumScalingBetheBloch{"cfgMomentumScalingBetheBloch", {nuclei::bbMomScalingDefault[0], 4, 2, nuclei::names, nuclei::chargeLabelNames}, "TPC Bethe-Bloch momentum scaling for light nuclei"}; @@ -247,15 +256,66 @@ struct NucleiSpectraTask { ConfigurableAxis cfgNsigmaTOFbins{"cfgNsigmaTOFbins", {100, -5., 5.}, "nsigma_TOF binning"}; ConfigurableAxis cfgTOFmassBins{"cfgTOFmassBins", {200, -5., 5.}, "TOF mass binning"}; + // CCDB options + Configurable cfgBz{"cfgBz", -999, "bz field, -999 is automatic"}; + 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"}; + Configurable cfgGRPpath{"cfgGRPpath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable cfgGRPmagPath{"cfgGRPmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable cfgLUTpath{"cfgLUTpath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; + Configurable cfgGeoPath{"cfgGeoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + int mRunNumber = 0; + float mBz = 0.f; + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; - Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (requireGlobalTrackWoDCAInFilter()); + Filter trackFilter = nabs(aod::track::eta) < cfgCutEta; + + using TrackCandidates = soa::Filtered>; - using TrackCandidates = soa::Filtered>; HistogramRegistry spectra{"spectra", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; o2::pid::tof::Beta responseBeta; + void initCCDB(aod::BCsWithTimestamps::iterator const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + auto run3grp_timestamp = bc.timestamp(); + + o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(cfgGRPpath, run3grp_timestamp); + o2::parameters::GRPMagField* grpmag = 0x0; + if (grpo) { + o2::base::Propagator::initFieldFromGRP(grpo); + if (cfgBz < -990) { + // Fetch magnetic field from ccdb for current collision + mBz = grpo->getNominalL3Field(); + LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << mBz << " kZG"; + } else { + mBz = cfgBz; + } + } else { + grpmag = ccdb->getForTimeStamp(cfgGRPmagPath, run3grp_timestamp); + if (!grpmag) { + LOG(fatal) << "Got nullptr from CCDB for path " << cfgGRPmagPath << " of object GRPMagField and " << cfgGRPpath << " of object GRPObject for timestamp " << run3grp_timestamp; + } + o2::base::Propagator::initFieldFromGRP(grpmag); + if (cfgBz < -990) { + // Fetch magnetic field from ccdb for current collision + mBz = std::lround(5.f * grpmag->getL3Current() / 30000.f); + LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << mBz << " kZG"; + } else { + mBz = cfgBz; + } + } + mRunNumber = bc.runNumber(); + } + void init(o2::framework::InitContext&) { + ccdb->setURL(cfgCCDBurl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(true); const AxisSpec centAxis{cfgCentralityBins, fmt::format("{} percentile", (std::string)cfgCentralityEstimator)}; const AxisSpec nSigmaAxes[2]{{cfgNsigmaTPCbins, "n#sigma_{TPC}"}, {cfgNsigmaTOFbins, "n#sigma_{TOF}"}}; @@ -298,6 +358,9 @@ struct NucleiSpectraTask { void fillDataInfo(soa::Filtered>::iterator const& collision, TrackCandidates const& tracks) { + auto bc = collision.bc_as(); + initCCDB(bc); + nuclei::candidates.clear(); // collision process loop if (!collision.sel8()) { @@ -305,6 +368,8 @@ struct NucleiSpectraTask { } spectra.fill(HIST("hRecVtxZData"), collision.posZ()); + const o2::math_utils::Point3D collVtx{collision.posX(), collision.posY(), collision.posZ()}; + const double bgScalings[4][2]{ {nuclei::charges[0] * cfgMomentumScalingBetheBloch->get(0u, 0u) / nuclei::masses[0], nuclei::charges[0] * cfgMomentumScalingBetheBloch->get(0u, 1u) / nuclei::masses[0]}, {nuclei::charges[1] * cfgMomentumScalingBetheBloch->get(1u, 0u) / nuclei::masses[1], nuclei::charges[1] * cfgMomentumScalingBetheBloch->get(1u, 1u) / nuclei::masses[1]}, @@ -314,42 +379,57 @@ struct NucleiSpectraTask { for (auto& track : tracks) { // start loop over tracks if (track.itsNCls() < cfgCutNclusITS || track.tpcNClsFound() < cfgCutNclusTPC || - std::abs(track.eta()) > cfgCutEta) { + track.tpcNClsCrossedRows() < 70 || + track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() || + track.tpcChi2NCl() > 4.f || + track.itsChi2NCl() > 36.f) { continue; } - const int iC{track.sign() < 0}; float nSigma[2][4]{ {track.tpcNSigmaDe(), track.tpcNSigmaTr(), track.tpcNSigmaHe(), track.tpcNSigmaAl()}, {track.tofNSigmaDe(), track.tofNSigmaTr(), track.tofNSigmaHe(), track.tofNSigmaAl()}}; + const int iC{track.sign() < 0}; + bool selectedTPC[4]{false}, goodToAnalyse{false}; + for (int iS{0}; iS < nuclei::species; ++iS) { + if (cfgBetheBlochParams->get(iS, 5u) > 0.f) { + double expBethe{tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * bgScalings[iS][iC]), cfgBetheBlochParams->get(iS, 0u), cfgBetheBlochParams->get(iS, 1u), cfgBetheBlochParams->get(iS, 2u), cfgBetheBlochParams->get(iS, 3u), cfgBetheBlochParams->get(iS, 4u))}; + double expSigma{expBethe * cfgBetheBlochParams->get(iS, 5u)}; + nSigma[0][iS] = static_cast((track.tpcSignal() - expBethe) / expSigma); + } + selectedTPC[iS] = (nSigma[0][iS] > nuclei::pidCuts[0][iS][0] && nSigma[0][iS] < nuclei::pidCuts[0][iS][1]); + goodToAnalyse = goodToAnalyse || selectedTPC[iS]; + } + if (!goodToAnalyse) { + continue; + } + + auto trackParCov = getTrackParCov(track); // should we set the charge according to the nucleus? + gpu::gpustd::array dcaInfo; + o2::base::Propagator::Instance()->propagateToDCABxByBz(collVtx, trackParCov, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfo); + float beta{responseBeta.GetBeta(track)}; - spectra.fill(HIST("hTofSignalData"), track.p(), beta); spectra.fill(HIST("hTpcSignalData"), track.tpcInnerParam() * track.sign(), track.tpcSignal()); + spectra.fill(HIST("hTofSignalData"), track.tpcInnerParam(), beta); for (int iS{0}; iS < nuclei::species; ++iS) { - bool selectedTPC{false}, selectedTOF{false}; - if (std::abs(track.dcaZ()) > cfgDCAcut->get(iS, 1)) { + bool selectedTOF{false}; + if (std::abs(dcaInfo[1]) > cfgDCAcut->get(iS, 1)) { continue; } - ROOT::Math::LorentzVector> fvector{track.pt() * nuclei::charges[iS], track.eta(), track.phi(), nuclei::masses[iS]}; + ROOT::Math::LorentzVector> fvector{trackParCov.getPt() * nuclei::charges[iS], trackParCov.getEta(), trackParCov.getPhi(), nuclei::masses[iS]}; float y{fvector.Rapidity() + cfgCMrapidity}; if (y < cfgCutRapidityMin || y > cfgCutRapidityMax) { continue; } - if (cfgBetheBlochParams->get(iS, 5u) > 0.f) { - double expBethe{tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * bgScalings[iS][iC]), cfgBetheBlochParams->get(iS, 0u), cfgBetheBlochParams->get(iS, 1u), cfgBetheBlochParams->get(iS, 2u), cfgBetheBlochParams->get(iS, 3u), cfgBetheBlochParams->get(iS, 4u))}; - double expSigma{expBethe * cfgBetheBlochParams->get(iS, 5u)}; - nSigma[0][iS] = static_cast((track.tpcSignal() - expBethe) / expSigma); - } for (int iPID{0}; iPID < 2; ++iPID) { - if (nSigma[0][iS] > nuclei::pidCuts[0][iS][0] && nSigma[0][iS] < nuclei::pidCuts[0][iS][1]) { - selectedTPC = true; + if (selectedTPC[iS]) { if (iPID && (!track.hasTOF() || nSigma[1][iS] < nuclei::pidCuts[1][iS][0] || nSigma[1][iS] > nuclei::pidCuts[1][iS][1])) { continue; } else if (iPID) { selectedTOF = true; } - nuclei::hDCAxy[iPID][iS][iC]->Fill(1., fvector.pt(), track.dcaXY()); - if (std::abs(track.dcaXY()) < cfgDCAcut->get(iS, 0u)) { + nuclei::hDCAxy[iPID][iS][iC]->Fill(1., fvector.pt(), dcaInfo[0]); + if (std::abs(dcaInfo[0]) < cfgDCAcut->get(iS, 0u)) { nuclei::hNsigma[iPID][iS][iC]->Fill(1., fvector.pt(), nSigma[iPID][iS]); if (iPID) { float mass{track.tpcInnerParam() * nuclei::charges[iS] * std::sqrt(1.f / (beta * beta) - 1.f) - nuclei::masses[iS]}; @@ -359,7 +439,7 @@ struct NucleiSpectraTask { } } uint16_t flag{kIsReconstructed}; - if (cfgTreeConfig->get(iS, 0u) && selectedTPC) { + if (cfgTreeConfig->get(iS, 0u) && selectedTPC[iS]) { int8_t massTOF{0u}; if (cfgTreeConfig->get(iS, 1u) && !selectedTOF) { continue; @@ -369,12 +449,12 @@ struct NucleiSpectraTask { massTOF = getBinnedValue(beta > 1.e-6f ? track.tpcInnerParam() * nuclei::charges[iS] * std::sqrt(1.f / (beta * beta) - 1.f) - nuclei::masses[iS] : -999.f, cfgBinnedVariables->get(4u, 1u)); } flag |= BIT(iS); - int8_t dcaxy = getBinnedValue(track.dcaXY(), cfgBinnedVariables->get(0u, 1u)); - int8_t dcaz = getBinnedValue(track.dcaZ(), cfgBinnedVariables->get(1u, 1u)); + int8_t dcaxy = getBinnedValue(dcaInfo[0], cfgBinnedVariables->get(0u, 1u)); + int8_t dcaz = getBinnedValue(dcaInfo[1], cfgBinnedVariables->get(1u, 1u)); int8_t nsigmaTPC = getBinnedValue(nSigma[0][iS], cfgBinnedVariables->get(2u, 1u)); int8_t nsigmaTOF = getBinnedValue(nSigma[1][iS], cfgBinnedVariables->get(3u, 1u)); - nuclei::candidates.emplace_back(track.globalIndex(), track.sign() * track.pt() * nuclei::charges[iS], track.eta(), track.itsClusterMap(), track.tpcNClsFound(), dcaxy, dcaz, flag, nsigmaTPC, nsigmaTOF, massTOF); + nuclei::candidates.emplace_back(track.globalIndex(), fvector.pt(), fvector.eta(), track.itsClusterMap(), track.tpcNClsFound(), dcaxy, dcaz, flag, nsigmaTPC, nsigmaTOF, massTOF); } } } // end loop over tracks From c49f2d24789f6f97770c71a13292bee3866219df Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Fri, 21 Apr 2023 15:00:44 +0200 Subject: [PATCH 2/5] Remove TOF pid full dependency --- PWGLF/Tasks/NucleiSpectraTask.cxx | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/PWGLF/Tasks/NucleiSpectraTask.cxx b/PWGLF/Tasks/NucleiSpectraTask.cxx index 1cb5b11678d..a379a6f8d08 100644 --- a/PWGLF/Tasks/NucleiSpectraTask.cxx +++ b/PWGLF/Tasks/NucleiSpectraTask.cxx @@ -77,7 +77,7 @@ float getBinCenter(uint8_t bin, double max) } struct NucleusCandidate { - NucleusCandidate(int idx, float aPt, float aEta, uint8_t aITSclsMap, uint8_t aTPCnCls, int8_t aDCAxy, int8_t aDCAz, uint16_t aFlags, uint8_t aTPCnsigma, uint8_t aTOFnsigma, uint8_t aTOFmass) : globalIndex(idx), pt(aPt), eta(aEta), ITSclsMap(aITSclsMap), TPCnCls(aTPCnCls), DCAxy(aDCAxy), DCAz(aDCAz), flags(aFlags), TPCnsigma(aTPCnsigma), TOFnsigma(aTOFnsigma), TOFmass(aTOFmass) + NucleusCandidate(int idx, float aPt, float aEta, uint8_t aITSclsMap, uint8_t aTPCnCls, int8_t aDCAxy, int8_t aDCAz, uint16_t aFlags, uint8_t aTPCnsigma, uint8_t aTOFmass) : globalIndex(idx), pt(aPt), eta(aEta), ITSclsMap(aITSclsMap), TPCnCls(aTPCnCls), DCAxy(aDCAxy), DCAz(aDCAz), flags(aFlags), TPCnsigma(aTPCnsigma), TOFmass(aTOFmass) { } int globalIndex; @@ -89,7 +89,6 @@ struct NucleusCandidate { int8_t DCAz; uint16_t flags; uint8_t TPCnsigma; - uint8_t TOFnsigma; uint8_t TOFmass; }; @@ -169,7 +168,6 @@ DECLARE_SOA_COLUMN(DCAxy, dcaxy, int8_t); DECLARE_SOA_COLUMN(DCAz, dcaz, int8_t); DECLARE_SOA_COLUMN(Flags, flags, uint16_t); DECLARE_SOA_COLUMN(TPCnsigma, tpcnsigma, uint8_t); -DECLARE_SOA_COLUMN(TOFnsigma, tofnsigma, uint8_t); DECLARE_SOA_COLUMN(TOFmass, tofmass, uint8_t); DECLARE_SOA_COLUMN(gPt, genPt, float); DECLARE_SOA_COLUMN(gEta, genEta, float); @@ -185,7 +183,6 @@ DECLARE_SOA_TABLE(NucleiTable, "AOD", "NUCLEITABLE", NucleiTableNS::DCAz, NucleiTableNS::Flags, NucleiTableNS::TPCnsigma, - NucleiTableNS::TOFnsigma, NucleiTableNS::TOFmass) DECLARE_SOA_TABLE(NucleiTableMC, "AOD", "NUCLEITABLEMC", @@ -197,7 +194,6 @@ DECLARE_SOA_TABLE(NucleiTableMC, "AOD", "NUCLEITABLEMC", NucleiTableNS::DCAz, NucleiTableNS::Flags, NucleiTableNS::TPCnsigma, - NucleiTableNS::TOFnsigma, NucleiTableNS::TOFmass, NucleiTableNS::gPt, NucleiTableNS::gEta, @@ -235,7 +231,6 @@ struct NucleiSpectraTask { Configurable> cfgMomentumScalingBetheBloch{"cfgMomentumScalingBetheBloch", {nuclei::bbMomScalingDefault[0], 4, 2, nuclei::names, nuclei::chargeLabelNames}, "TPC Bethe-Bloch momentum scaling for light nuclei"}; Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {nuclei::betheBlochDefault[0], 4, 6, nuclei::names, nuclei::betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; Configurable> cfgNsigmaTPC{"cfgNsigmaTPC", {nuclei::nSigmaTPCdefault[0], 4, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"}; - Configurable> cfgNsigmaTOF{"cfgNsigmaTOF", {nuclei::nSigmaTOFdefault[0], 4, 2, nuclei::names, nuclei::nSigmaConfigName}, "TOF nsigma selection for light nuclei"}; Configurable> cfgDCAcut{"cfgDCAcut", {nuclei::DCAcutDefault[0], 4, 2, nuclei::names, nuclei::nDCAConfigName}, "Max DCAxy and DCAz for light nuclei"}; Configurable> cfgTreeConfig{"cfgTreeConfig", {nuclei::TreeConfigDefault[0], 4, 2, nuclei::names, nuclei::treeConfigNames}, "Filtered trees configuration"}; Configurable> cfgBinnedVariables{"cfgBinnedVariables", {nuclei::BinnedVariablesDefaultMax[0], 5, 1, nuclei::binnedVariableNames, nuclei::binnedLabelNames}, "Maximum value for the binned variables"}; @@ -270,7 +265,7 @@ struct NucleiSpectraTask { Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; Filter trackFilter = nabs(aod::track::eta) < cfgCutEta; - using TrackCandidates = soa::Filtered>; + using TrackCandidates = soa::Filtered>; HistogramRegistry spectra{"spectra", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; o2::pid::tof::Beta responseBeta; @@ -338,8 +333,9 @@ struct NucleiSpectraTask { spectra.add("hTofSignalData", "TOF beta", HistType::kTH2F, {{500, 0., 5., "#it{p} (GeV/#it{c})"}, {750, 0, 1.5, "TOF #beta"}}); for (int iC{0}; iC < 2; ++iC) { for (int iS{0}; iS < nuclei::species; ++iS) { + nuclei::hNsigma[0][iS][iC] = spectra.add(fmt::format("h{}nsigma{}_{}", nuclei::pidName[0], nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("n#sigma_{{}} {} {}", nuclei::pidName[0], nuclei::matter[iC], nuclei::names[iS]).data(), HistType::kTH3D, {centAxis, ptAxes[iS], nSigmaAxes[0]}); + for (int iPID{0}; iPID < 2; ++iPID) { - nuclei::hNsigma[iPID][iS][iC] = spectra.add(fmt::format("h{}nsigma{}_{}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("n#sigma_{{}} {} {}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), HistType::kTH3D, {centAxis, ptAxes[iS], nSigmaAxes[iPID]}); nuclei::hDCAxy[iPID][iS][iC] = spectra.add(fmt::format("hDCAxy{}_{}_{}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("DCAxy {} {} {}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), HistType::kTH3D, {centAxis, ptAxes[iS], dcaAxes[iS]}); } nuclei::hTOFmass[iS][iC] = spectra.add(fmt::format("h{}TOFmass{}", nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("TOF mass - {} PDG mass", nuclei::names[iS]).data(), HistType::kTH3D, {centAxis, ptAxes[iS], tofMassAxis}); @@ -351,7 +347,6 @@ struct NucleiSpectraTask { for (int iS{0}; iS < 4; ++iS) { for (int iMax{0}; iMax < 2; ++iMax) { nuclei::pidCuts[0][iS][iMax] = cfgNsigmaTPC->get(iS, iMax); - nuclei::pidCuts[1][iS][iMax] = cfgNsigmaTOF->get(iS, iMax); } } } @@ -387,7 +382,7 @@ struct NucleiSpectraTask { } float nSigma[2][4]{ {track.tpcNSigmaDe(), track.tpcNSigmaTr(), track.tpcNSigmaHe(), track.tpcNSigmaAl()}, - {track.tofNSigmaDe(), track.tofNSigmaTr(), track.tofNSigmaHe(), track.tofNSigmaAl()}}; + {0.f, 0.f, 0.f, 0.f}}; /// then we will calibrate the TOF mass for the He3 and Alpha const int iC{track.sign() < 0}; bool selectedTPC[4]{false}, goodToAnalyse{false}; for (int iS{0}; iS < nuclei::species; ++iS) { @@ -423,7 +418,7 @@ struct NucleiSpectraTask { for (int iPID{0}; iPID < 2; ++iPID) { if (selectedTPC[iS]) { - if (iPID && (!track.hasTOF() || nSigma[1][iS] < nuclei::pidCuts[1][iS][0] || nSigma[1][iS] > nuclei::pidCuts[1][iS][1])) { + if (iPID && !track.hasTOF()) { continue; } else if (iPID) { selectedTOF = true; @@ -452,9 +447,8 @@ struct NucleiSpectraTask { int8_t dcaxy = getBinnedValue(dcaInfo[0], cfgBinnedVariables->get(0u, 1u)); int8_t dcaz = getBinnedValue(dcaInfo[1], cfgBinnedVariables->get(1u, 1u)); int8_t nsigmaTPC = getBinnedValue(nSigma[0][iS], cfgBinnedVariables->get(2u, 1u)); - int8_t nsigmaTOF = getBinnedValue(nSigma[1][iS], cfgBinnedVariables->get(3u, 1u)); - nuclei::candidates.emplace_back(track.globalIndex(), fvector.pt(), fvector.eta(), track.itsClusterMap(), track.tpcNClsFound(), dcaxy, dcaz, flag, nsigmaTPC, nsigmaTOF, massTOF); + nuclei::candidates.emplace_back(track.globalIndex(), fvector.pt(), fvector.eta(), track.itsClusterMap(), track.tpcNClsFound(), dcaxy, dcaz, flag, nsigmaTPC, massTOF); } } } // end loop over tracks @@ -464,7 +458,7 @@ struct NucleiSpectraTask { { fillDataInfo(collision, tracks); for (auto& c : nuclei::candidates) { - nucleiTable(c.pt, c.eta, c.ITSclsMap, c.TPCnCls, c.DCAxy, c.DCAz, c.flags, c.TPCnsigma, c.TOFnsigma, c.TOFmass); + nucleiTable(c.pt, c.eta, c.ITSclsMap, c.TPCnCls, c.DCAxy, c.DCAz, c.flags, c.TPCnsigma, c.TOFmass); } } PROCESS_SWITCH(NucleiSpectraTask, processData, "Data analysis", true); @@ -488,7 +482,7 @@ struct NucleiSpectraTask { c.flags |= kIsSecondaryFromMaterial; } - nucleiTableMC(c.pt, c.eta, c.ITSclsMap, c.TPCnCls, c.DCAxy, c.DCAz, c.flags, c.TPCnsigma, c.TOFnsigma, c.TOFmass, particle.pt(), particle.eta(), particle.pdgCode()); + nucleiTableMC(c.pt, c.eta, c.ITSclsMap, c.TPCnCls, c.DCAxy, c.DCAz, c.flags, c.TPCnsigma, c.TOFmass, particle.pt(), particle.eta(), particle.pdgCode()); for (int iS{0}; iS < nuclei::species; ++iS) { if (std::abs(particle.pdgCode()) == nuclei::codes[iS]) { nuclei::hMomRes[iS][particle.pdgCode() < 0]->Fill(1., std::abs(c.pt), 1. - std::abs(c.pt) / particle.pt()); @@ -517,7 +511,7 @@ struct NucleiSpectraTask { } if (!isReconstructed[index] && (cfgTreeConfig->get(iS, 0u) || cfgTreeConfig->get(iS, 1u))) { - nucleiTableMC(0, 0, 0, 0, 0, 0, flags, 0, 0, 0, particle.pt(), particle.eta(), particle.pdgCode()); + nucleiTableMC(0, 0, 0, 0, 0, 0, flags, 0, 0, particle.pt(), particle.eta(), particle.pdgCode()); } break; } From 36681b069ec5122def7b65e47bf324a288d3e779 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Fri, 21 Apr 2023 15:32:06 +0200 Subject: [PATCH 3/5] Final touches --- PWGLF/Tasks/NucleiSpectraTask.cxx | 32 +++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/PWGLF/Tasks/NucleiSpectraTask.cxx b/PWGLF/Tasks/NucleiSpectraTask.cxx index a379a6f8d08..503f86459d3 100644 --- a/PWGLF/Tasks/NucleiSpectraTask.cxx +++ b/PWGLF/Tasks/NucleiSpectraTask.cxx @@ -16,9 +16,7 @@ // // Data (run3): // o2-analysis-lf-nuclei-spectra, o2-analysis-timestamp -// o2-analysis-pid-tof-base, o2-analysis-pid-tof-full -// o2-analysis-pid-tpc-base, o2-analysis-pid-tpc-full, -// o2-analysis-multiplicity-table, o2-analysis-event-selection +// o2-analysis-pid-tof-base, o2-analysis-multiplicity-table, o2-analysis-event-selection #include @@ -152,6 +150,7 @@ std::shared_ptr hTOFmass[4][2]; std::shared_ptr hGenNuclei[4][2]; std::shared_ptr hMomRes[4][2]; std::shared_ptr hDCAxy[2][4][2]; +o2::base::MatLayerCylSet* lut = nullptr; std::vector candidates; } // namespace nuclei @@ -265,7 +264,7 @@ struct NucleiSpectraTask { Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; Filter trackFilter = nabs(aod::track::eta) < cfgCutEta; - using TrackCandidates = soa::Filtered>; + using TrackCandidates = soa::Filtered>; HistogramRegistry spectra{"spectra", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; o2::pid::tof::Beta responseBeta; @@ -303,6 +302,8 @@ struct NucleiSpectraTask { } } mRunNumber = bc.runNumber(); + o2::base::Propagator::initFieldFromGRP(grpmag); + o2::base::Propagator::Instance()->setMatLUT(nuclei::lut); } void init(o2::framework::InitContext&) @@ -310,7 +311,7 @@ struct NucleiSpectraTask { ccdb->setURL(cfgCCDBurl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(true); + ccdb->setFatalWhenNull(false); const AxisSpec centAxis{cfgCentralityBins, fmt::format("{} percentile", (std::string)cfgCentralityEstimator)}; const AxisSpec nSigmaAxes[2]{{cfgNsigmaTPCbins, "n#sigma_{TPC}"}, {cfgNsigmaTOFbins, "n#sigma_{TOF}"}}; @@ -349,6 +350,9 @@ struct NucleiSpectraTask { nuclei::pidCuts[0][iS][iMax] = cfgNsigmaTPC->get(iS, iMax); } } + + nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get("GLO/Param/MatLUT")); + } void fillDataInfo(soa::Filtered>::iterator const& collision, TrackCandidates const& tracks) @@ -381,16 +385,14 @@ struct NucleiSpectraTask { continue; } float nSigma[2][4]{ - {track.tpcNSigmaDe(), track.tpcNSigmaTr(), track.tpcNSigmaHe(), track.tpcNSigmaAl()}, + {-10., -10., -10., -10.}, {0.f, 0.f, 0.f, 0.f}}; /// then we will calibrate the TOF mass for the He3 and Alpha const int iC{track.sign() < 0}; bool selectedTPC[4]{false}, goodToAnalyse{false}; for (int iS{0}; iS < nuclei::species; ++iS) { - if (cfgBetheBlochParams->get(iS, 5u) > 0.f) { - double expBethe{tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * bgScalings[iS][iC]), cfgBetheBlochParams->get(iS, 0u), cfgBetheBlochParams->get(iS, 1u), cfgBetheBlochParams->get(iS, 2u), cfgBetheBlochParams->get(iS, 3u), cfgBetheBlochParams->get(iS, 4u))}; - double expSigma{expBethe * cfgBetheBlochParams->get(iS, 5u)}; - nSigma[0][iS] = static_cast((track.tpcSignal() - expBethe) / expSigma); - } + double expBethe{tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * bgScalings[iS][iC]), cfgBetheBlochParams->get(iS, 0u), cfgBetheBlochParams->get(iS, 1u), cfgBetheBlochParams->get(iS, 2u), cfgBetheBlochParams->get(iS, 3u), cfgBetheBlochParams->get(iS, 4u))}; + double expSigma{expBethe * cfgBetheBlochParams->get(iS, 5u)}; + nSigma[0][iS] = static_cast((track.tpcSignal() - expBethe) / expSigma); selectedTPC[iS] = (nSigma[0][iS] > nuclei::pidCuts[0][iS][0] && nSigma[0][iS] < nuclei::pidCuts[0][iS][1]); goodToAnalyse = goodToAnalyse || selectedTPC[iS]; } @@ -425,7 +427,9 @@ struct NucleiSpectraTask { } nuclei::hDCAxy[iPID][iS][iC]->Fill(1., fvector.pt(), dcaInfo[0]); if (std::abs(dcaInfo[0]) < cfgDCAcut->get(iS, 0u)) { - nuclei::hNsigma[iPID][iS][iC]->Fill(1., fvector.pt(), nSigma[iPID][iS]); + if (!iPID) { /// temporary exclusion of the TOF nsigma PID for the He3 and Alpha + nuclei::hNsigma[iPID][iS][iC]->Fill(1., fvector.pt(), nSigma[iPID][iS]); + } if (iPID) { float mass{track.tpcInnerParam() * nuclei::charges[iS] * std::sqrt(1.f / (beta * beta) - 1.f) - nuclei::masses[iS]}; nuclei::hTOFmass[iS][iC]->Fill(1., fvector.pt(), mass); @@ -454,7 +458,7 @@ struct NucleiSpectraTask { } // end loop over tracks } - void processData(soa::Filtered>::iterator const& collision, TrackCandidates const& tracks) + void processData(soa::Filtered>::iterator const& collision, TrackCandidates const& tracks, aod::BCsWithTimestamps const&) { fillDataInfo(collision, tracks); for (auto& c : nuclei::candidates) { @@ -463,7 +467,7 @@ struct NucleiSpectraTask { } PROCESS_SWITCH(NucleiSpectraTask, processData, "Data analysis", true); - void processMC(soa::Filtered>::iterator const& collision, TrackCandidates const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC) + void processMC(soa::Filtered>::iterator const& collision, TrackCandidates const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, aod::BCsWithTimestamps const&) { fillDataInfo(collision, tracks); std::vector isReconstructed(particlesMC.size(), false); From 4c357aa52cd5701caf1f97b395aa150eccd68f32 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 21 Apr 2023 13:36:51 +0000 Subject: [PATCH 4/5] Please consider the following formatting changes --- PWGLF/Tasks/NucleiSpectraTask.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/PWGLF/Tasks/NucleiSpectraTask.cxx b/PWGLF/Tasks/NucleiSpectraTask.cxx index 503f86459d3..1ca09e94d15 100644 --- a/PWGLF/Tasks/NucleiSpectraTask.cxx +++ b/PWGLF/Tasks/NucleiSpectraTask.cxx @@ -264,7 +264,7 @@ struct NucleiSpectraTask { Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; Filter trackFilter = nabs(aod::track::eta) < cfgCutEta; - using TrackCandidates = soa::Filtered>; + using TrackCandidates = soa::Filtered>; HistogramRegistry spectra{"spectra", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; o2::pid::tof::Beta responseBeta; @@ -352,7 +352,6 @@ struct NucleiSpectraTask { } nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get("GLO/Param/MatLUT")); - } void fillDataInfo(soa::Filtered>::iterator const& collision, TrackCandidates const& tracks) From 8c669ba5d8a043c3081d02621e066b2153332ab2 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Mon, 24 Apr 2023 16:06:04 +0200 Subject: [PATCH 5/5] Fix complain of osx CI --- PWGLF/Tasks/NucleiSpectraTask.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/PWGLF/Tasks/NucleiSpectraTask.cxx b/PWGLF/Tasks/NucleiSpectraTask.cxx index 1ca09e94d15..f0bb154d8fe 100644 --- a/PWGLF/Tasks/NucleiSpectraTask.cxx +++ b/PWGLF/Tasks/NucleiSpectraTask.cxx @@ -107,11 +107,11 @@ constexpr double nSigmaTPCdefault[4][2]{ {-5., 5.}, {-5., 5.}, {-5., 5.}}; -constexpr double nSigmaTOFdefault[4][2]{ - {-5., 5.}, - {-5., 5.}, - {-5., 5.}, - {-5., 5.}}; +// constexpr double nSigmaTOFdefault[4][2]{ +// {-5., 5.}, +// {-5., 5.}, +// {-5., 5.}, +// {-5., 5.}}; constexpr double DCAcutDefault[4][2]{ {1., 1.}, {1., 1.},