From 2f5e710db243ea4090edbc2cd684d8b3bfec263f Mon Sep 17 00:00:00 2001 From: Maxim Virta Date: Wed, 28 Oct 2026 13:51:43 +0100 Subject: [PATCH 1/3] [PWGCF] Bit masks added for derived data --- .../TableProducer/CMakeLists.txt | 5 + PWGCF/GenericFramework/Tasks/flowGfwV02.cxx | 326 +++++++++++------- 2 files changed, 200 insertions(+), 131 deletions(-) diff --git a/PWGCF/GenericFramework/TableProducer/CMakeLists.txt b/PWGCF/GenericFramework/TableProducer/CMakeLists.txt index bbfd7adac2b..ab012ec889e 100644 --- a/PWGCF/GenericFramework/TableProducer/CMakeLists.txt +++ b/PWGCF/GenericFramework/TableProducer/CMakeLists.txt @@ -8,3 +8,8 @@ # 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. + +o2physics_add_dpl_workflow(flow-gfw-derived-producer + SOURCES flowGfwDerivedProducer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx b/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx index 5768bb877f7..ad8dcfef573 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx @@ -17,6 +17,7 @@ #include "PWGCF/GenericFramework/Core/GFW.h" #include "PWGCF/GenericFramework/Core/GFWConfig.h" #include "PWGCF/GenericFramework/Core/GFWWeights.h" +#include "PWGCF/GenericFramework/DataModel/GfwDerived.h" #include "PWGCF/JCorran/DataModel/JCatalyst.h" #include "Common/CCDB/EventSelectionParams.h" @@ -213,8 +214,14 @@ struct FlowGfwV02 { // // Filters to be applied to the received data. // // The analysis assumes the data has been subjected to a QA of its selection, // // and thus only the final distributions of the data for analysis are saved. + // Derived-data bit masks. Every listed track bit is required; see gfwcollision::EventSel and gfwtrack::TrackSelection. + O2_DEFINE_CONFIGURABLE(cfgEventBitMask, uint16_t, o2::aod::gfwcollision::DefaultEventMask, "Required event-selection bits for GFW derived data"); + Configurable> cfgBitmaskSelections{"cfgBitmaskSelections", {0, 3, 6, 9, 12}, "Enabled TrackSelection bit indices for GFW derived data. DCA 0-5, TPC track 6-13, TPC PID 14-22, TOF PID 23-31. The default is the nominal quality selection"}; + O2_DEFINE_CONFIGURABLE(cfgPidWindow, int, o2::aod::gfwtrack::PidWindowNominal, "PID nSigma window for GFW derived data: 0 nominal, 1 tight, 2 loose"); + o2::framework::expressions::Filter collFilter = (nabs(aod::collision::posZ) < cfgEventCuts.cfgZvtxMax); o2::framework::expressions::Filter trackFilter = (aod::track::pt > cfgTrackCuts.cfgPtMin) && (aod::track::pt < cfgTrackCuts.cfgPtMax) && nabs(aod::track::eta) < cfgTrackCuts.cfgEtaMax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && (aod::track::itsChi2NCl < cfgTrackCuts.cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgTrackCuts.cfgDCAz; + o2::framework::expressions::Filter gfwTrackFilter = (aod::gfwtrack::pt > cfgTrackCuts.cfgPtMin) && (aod::gfwtrack::pt < cfgTrackCuts.cfgPtMax) && (nabs(aod::gfwtrack::eta) < cfgTrackCuts.cfgEtaMax); // Connect to ccdb Service ccdb{}; @@ -292,6 +299,7 @@ struct FlowGfwV02 { std::unique_ptr fMultPVGlobalCutHigh; std::unique_ptr fPtDepDCAxy; + uint32_t trackBitMask = o2::aod::gfwtrack::NominalMask; using GFWTracks = soa::Filtered>; @@ -313,6 +321,12 @@ struct FlowGfwV02 { void init(InitContext const&) { + trackBitMask = 0; + for (const int bit : cfgBitmaskSelections.value) { + if (bit >= 0 && bit < 32) { + trackBitMask |= static_cast(1u << bit); + } + } projectMatrix(nSigmas->getData(), pidStates.tpcNsigmaCut, pidStates.tofNsigmaCut, pidStates.itsNsigmaCut); projectMatrix(nSigmasHighPurity->getData(), pidStates.highPurityTpcNsigmaCut, pidStates.highPurityTofNsigmaCut, pidStates.highPurityItsNsigmaCut); @@ -605,39 +619,61 @@ struct FlowGfwV02 { template int getNsigmaPID(TTrack const& track) { - // Computing Nsigma arrays for pion, kaon, and protons - std::array nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; - std::array nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; - std::array nSigmaITS = {pidStates.itsResponse.nSigmaITS(track), pidStates.itsResponse.nSigmaITS(track), pidStates.itsResponse.nSigmaITS(track)}; int pid = -1; // -1 = not identified, 1 = pion, 2 = kaon, 3 = proton - if (cfgUseHighPurityPID && track.pt() > cfgHighPurityPtMin) { - bool isHighPurityPion = nSigmaTPC[IndPionUp] < pidStates.highPurityTpcNsigmaCut[IndPionUp] && nSigmaTPC[IndPionUp] > pidStates.highPurityTpcNsigmaCut[IndPionLow] && nSigmaTOF[IndPionUp] < pidStates.highPurityTofNsigmaCut[IndPionUp] && nSigmaTOF[IndPionUp] > pidStates.highPurityTofNsigmaCut[IndPionLow]; - bool isHighPurityKaon = nSigmaTPC[IndKaonUp] < pidStates.highPurityTpcNsigmaCut[IndKaonUp] && nSigmaTPC[IndKaonUp] > pidStates.highPurityTpcNsigmaCut[IndKaonLow] && nSigmaTOF[IndKaonUp] < pidStates.highPurityTofNsigmaCut[IndKaonUp] && nSigmaTOF[IndKaonUp] > pidStates.highPurityTofNsigmaCut[IndKaonLow]; - bool isHighPurityProton = nSigmaTPC[IndProtonUp] < pidStates.highPurityTpcNsigmaCut[IndProtonUp] && nSigmaTPC[IndProtonUp] > pidStates.highPurityTpcNsigmaCut[IndProtonLow] && nSigmaTOF[IndProtonUp] < pidStates.highPurityTofNsigmaCut[IndProtonUp] && nSigmaTOF[IndProtonUp] > pidStates.highPurityTofNsigmaCut[IndProtonLow]; - if (isHighPurityPion) { - pid = PidPions; - } else if (isHighPurityKaon) { - pid = PidKaons; - } else if (isHighPurityProton) { - pid = PidProtons; - } else { + bool isPion = false, isKaon = false, isProton = false; + bool isDetectedPion = false; + bool isDetectedKaon = false; + bool isDetectedProton = false; + bool isTofPion = false; + bool isTofKaon = false; + bool isTofProton = false; + if constexpr (requires { track.trackSel(); }) { + const int window = cfgPidWindow; + const int shift = (window == aod::gfwtrack::PidWindowTight || window == aod::gfwtrack::PidWindowLoose) ? window * aod::gfwtrack::PidWindowBits : aod::gfwtrack::PidWindowNominal; + const uint32_t trackSel = track.trackSel(); + isDetectedPion = (trackSel & (static_cast(aod::gfwtrack::kTpcPi) << shift)) != 0; + isDetectedKaon = (trackSel & (static_cast(aod::gfwtrack::kTpcKa) << shift)) != 0; + isDetectedProton = (trackSel & (static_cast(aod::gfwtrack::kTpcPr) << shift)) != 0; + isTofPion = (trackSel & (static_cast(aod::gfwtrack::kTofPi) << shift)) != 0; + isTofKaon = (trackSel & (static_cast(aod::gfwtrack::kTofKa) << shift)) != 0; + isTofProton = (trackSel & (static_cast(aod::gfwtrack::kTofPr) << shift)) != 0; + } else { + const std::array nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; + const std::array nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; + constexpr bool HasItsResponse = requires { + track.itsClusterSizes(); + track.p(); + }; + if (cfgUseHighPurityPID && track.pt() > cfgHighPurityPtMin) { + const bool isHighPurityPion = nSigmaTPC[IndPionUp] < pidStates.highPurityTpcNsigmaCut[IndPionUp] && nSigmaTPC[IndPionUp] > pidStates.highPurityTpcNsigmaCut[IndPionLow] && nSigmaTOF[IndPionUp] < pidStates.highPurityTofNsigmaCut[IndPionUp] && nSigmaTOF[IndPionUp] > pidStates.highPurityTofNsigmaCut[IndPionLow]; + const bool isHighPurityKaon = nSigmaTPC[IndKaonUp] < pidStates.highPurityTpcNsigmaCut[IndKaonUp] && nSigmaTPC[IndKaonUp] > pidStates.highPurityTpcNsigmaCut[IndKaonLow] && nSigmaTOF[IndKaonUp] < pidStates.highPurityTofNsigmaCut[IndKaonUp] && nSigmaTOF[IndKaonUp] > pidStates.highPurityTofNsigmaCut[IndKaonLow]; + const bool isHighPurityProton = nSigmaTPC[IndProtonUp] < pidStates.highPurityTpcNsigmaCut[IndProtonUp] && nSigmaTPC[IndProtonUp] > pidStates.highPurityTpcNsigmaCut[IndProtonLow] && nSigmaTOF[IndProtonUp] < pidStates.highPurityTofNsigmaCut[IndProtonUp] && nSigmaTOF[IndProtonUp] > pidStates.highPurityTofNsigmaCut[IndProtonLow]; + if (isHighPurityPion) { + return PidPions; + } + if (isHighPurityKaon) { + return PidKaons; + } + if (isHighPurityProton) { + return PidProtons; + } return -1; } - return pid; + std::array nSigmaITS{}; + if constexpr (HasItsResponse) { + nSigmaITS = {pidStates.itsResponse.nSigmaITS(track), pidStates.itsResponse.nSigmaITS(track), pidStates.itsResponse.nSigmaITS(track)}; + } + const bool useItsPid = cfgUseItsPID && HasItsResponse; + const std::array nSigmaToUse = useItsPid ? nSigmaITS : nSigmaTPC; + const std::array detectorNsigmaCut = useItsPid ? pidStates.itsNsigmaCut : pidStates.tpcNsigmaCut; + isDetectedPion = nSigmaToUse[IndPionUp] < detectorNsigmaCut[IndPionUp] && nSigmaToUse[IndPionUp] > detectorNsigmaCut[IndPionLow]; + isDetectedKaon = nSigmaToUse[IndKaonUp] < detectorNsigmaCut[IndKaonUp] && nSigmaToUse[IndKaonUp] > detectorNsigmaCut[IndKaonLow]; + isDetectedProton = nSigmaToUse[IndProtonUp] < detectorNsigmaCut[IndProtonUp] && nSigmaToUse[IndProtonUp] > detectorNsigmaCut[IndProtonLow]; + isTofPion = nSigmaTOF[IndPionUp] < pidStates.tofNsigmaCut[IndPionUp] && nSigmaTOF[IndPionUp] > pidStates.tofNsigmaCut[IndPionLow]; + isTofKaon = nSigmaTOF[IndKaonUp] < pidStates.tofNsigmaCut[IndKaonUp] && nSigmaTOF[IndKaonUp] > pidStates.tofNsigmaCut[IndKaonLow]; + isTofProton = nSigmaTOF[IndProtonUp] < pidStates.tofNsigmaCut[IndProtonUp] && nSigmaTOF[IndProtonUp] > pidStates.tofNsigmaCut[IndProtonLow]; } - std::array nSigmaToUse = cfgUseItsPID ? nSigmaITS : nSigmaTPC; // Choose which nSigma to use: TPC or ITS - std::array detectorNsigmaCut = cfgUseItsPID ? pidStates.itsNsigmaCut : pidStates.tpcNsigmaCut; // Choose which nSigma to use: TPC or ITS - - bool isPion = false, isKaon = false, isProton = false; - bool isDetectedPion = nSigmaToUse[IndPionUp] < detectorNsigmaCut[IndPionUp] && nSigmaToUse[IndPionUp] > detectorNsigmaCut[IndPionLow]; - bool isDetectedKaon = nSigmaToUse[IndKaonUp] < detectorNsigmaCut[IndKaonUp] && nSigmaToUse[IndKaonUp] > detectorNsigmaCut[IndKaonLow]; - bool isDetectedProton = nSigmaToUse[IndProtonUp] < detectorNsigmaCut[IndProtonUp] && nSigmaToUse[IndProtonUp] > detectorNsigmaCut[IndProtonLow]; - - bool isTofPion = nSigmaTOF[IndPionUp] < pidStates.tofNsigmaCut[IndPionUp] && nSigmaTOF[IndPionUp] > pidStates.tofNsigmaCut[IndPionLow]; - bool isTofKaon = nSigmaTOF[IndKaonUp] < pidStates.tofNsigmaCut[IndKaonUp] && nSigmaTOF[IndKaonUp] > pidStates.tofNsigmaCut[IndKaonLow]; - bool isTofProton = nSigmaTOF[IndProtonUp] < pidStates.tofNsigmaCut[IndProtonUp] && nSigmaTOF[IndProtonUp] > pidStates.tofNsigmaCut[IndProtonLow]; - if (track.pt() > cfgTofPtCut && !track.hasTOF()) { return -1; } @@ -905,60 +941,65 @@ struct FlowGfwV02 { double ptMeanMid = pidStates.hPtMid[PidCharged]->GetMean(); double ptFractionMid = 0.; double dnxAB = fGFW->Calculate(corrconfigs.at(0), 0, kTRUE).real(); // V22 weight for AB - auto valAB = fGFW->Calculate(corrconfigs.at(0), 0, kFALSE).real() / dnxAB; double dnxABAB = fGFW->Calculate(corrconfigs.at(1), 0, kTRUE).real(); // V22 weight for ABAB - auto valABAB = fGFW->Calculate(corrconfigs.at(1), 0, kFALSE).real() / dnxABAB; - if (std::abs(valAB) > Threshold || std::isnan(valAB) || std::isinf(valAB) || dnxAB < MinDnxAB || std::abs(valABAB) > Threshold || std::isnan(valABAB) || std::isinf(valABAB) || dnxABAB < MinDnxAB) { - return; - } - double v22pt = valAB * ptMeanMid; - double weightAB = (cfgUseMultiplicityFlowWeights) ? dnxAB : 1.0; - double weightABAB = (cfgUseMultiplicityFlowWeights) ? dnxABAB : 1.0; - double weightC = 1.0; - // Calculate V02 for each particle type - for (int pid = 0; pid < PidTotal; pid++) { - int normIndex = (cfgNormalizeByCharged) ? PidCharged : pid; - if (!(pidStates.hPtMid[normIndex]->Integral() > 0)) { - continue; // Mid pT distribution is not defined - } - if (cfgUseMultiplicityFracWeights) { - weightC = pidStates.hPtMid[PidCharged]->Integral(); // Mean pt/npt C weight - } - for (int i = 1; i <= fSecondAxis->GetNbins(); i++) { - ptFractionMid = pidStates.hPtMid[pid]->GetBinContent(i) / pidStates.hPtMid[normIndex]->Integral(); - v22npt = valAB * ptFractionMid; - switch (pid) { - case PidCharged: - registry.fill(HIST("analysis/charged/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); - registry.fill(HIST("analysis/charged/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); - if (cfgUseMultiplicityFracWeights) { - registry.fill(HIST("analysis/charged/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); - } - break; - case PidPions: - registry.fill(HIST("analysis/pion/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); - registry.fill(HIST("analysis/pion/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); - if (cfgUseMultiplicityFracWeights) { - registry.fill(HIST("analysis/pion/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); - } - break; - case PidKaons: - registry.fill(HIST("analysis/kaon/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); - registry.fill(HIST("analysis/kaon/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); - if (cfgUseMultiplicityFracWeights) { - registry.fill(HIST("analysis/kaon/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); - } - break; - case PidProtons: - registry.fill(HIST("analysis/proton/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); - registry.fill(HIST("analysis/proton/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); - if (cfgUseMultiplicityFracWeights) { - registry.fill(HIST("analysis/proton/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); - } - break; - default: - break; + // Check if the V22 weight is valid + const bool validAB = dnxAB >= MinDnxAB; + const bool validABAB = dnxABAB >= MinDnxAB; + const double valAB = validAB ? fGFW->Calculate(corrconfigs.at(0), 0, kFALSE).real() / dnxAB : 0.; + const double valABAB = validABAB ? fGFW->Calculate(corrconfigs.at(1), 0, kFALSE).real() / dnxABAB : 0.; + const bool acceptAB = validAB && std::isfinite(valAB) && std::abs(valAB) <= Threshold; + const bool acceptABAB = validABAB && std::isfinite(valABAB) && std::abs(valABAB) <= Threshold; + if (acceptAB && acceptABAB) { + double v22pt = valAB * ptMeanMid; + double weightAB = (cfgUseMultiplicityFlowWeights) ? dnxAB : 1.0; + double weightABAB = (cfgUseMultiplicityFlowWeights) ? dnxABAB : 1.0; + double weightC = 1.0; + // Calculate V02 for each particle type + for (int pid = 0; pid < PidTotal; pid++) { + int normIndex = (cfgNormalizeByCharged) ? PidCharged : pid; + if (!(pidStates.hPtMid[normIndex]->Integral() > 0)) { + continue; // Mid pT distribution is not defined + } + if (cfgUseMultiplicityFracWeights) { + weightC = pidStates.hPtMid[PidCharged]->Integral(); // Mean pt/npt C weight + } + for (int i = 1; i <= fSecondAxis->GetNbins(); i++) { + ptFractionMid = pidStates.hPtMid[pid]->GetBinContent(i) / pidStates.hPtMid[normIndex]->Integral(); + v22npt = valAB * ptFractionMid; + + switch (pid) { + case PidCharged: + registry.fill(HIST("analysis/charged/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); + registry.fill(HIST("analysis/charged/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); + if (cfgUseMultiplicityFracWeights) { + registry.fill(HIST("analysis/charged/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); + } + break; + case PidPions: + registry.fill(HIST("analysis/pion/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); + registry.fill(HIST("analysis/pion/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); + if (cfgUseMultiplicityFracWeights) { + registry.fill(HIST("analysis/pion/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); + } + break; + case PidKaons: + registry.fill(HIST("analysis/kaon/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); + registry.fill(HIST("analysis/kaon/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); + if (cfgUseMultiplicityFracWeights) { + registry.fill(HIST("analysis/kaon/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); + } + break; + case PidProtons: + registry.fill(HIST("analysis/proton/v22npt"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, v22npt, weightC * weightAB); + registry.fill(HIST("analysis/proton/nchC"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC); + if (cfgUseMultiplicityFracWeights) { + registry.fill(HIST("analysis/proton/nchC3pc"), bootstrap, fSecondAxis->GetBinCenter(i), centmult, ptFractionMid, weightC * weightAB); + } + break; + default: + break; + } } } } @@ -1186,10 +1227,15 @@ struct FlowGfwV02 { template bool trackSelected(TTrack const& track) { - if (cfgTrackCuts.cfgDCAxyNSigma && (std::fabs(track.dcaXY()) > fPtDepDCAxy->Eval(track.pt()))) { - return false; + if constexpr (requires { track.trackSel(); }) { + const uint32_t mask = trackBitMask; + return (track.trackSel() & mask) == mask; + } else { + if (cfgTrackCuts.cfgDCAxyNSigma && (std::fabs(track.dcaXY()) > fPtDepDCAxy->Eval(track.pt()))) { + return false; + } + return ((track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgNTPCXrows) && (track.tpcNClsFound() >= cfgTrackCuts.cfgNTPCCls) && (track.itsNCls() >= cfgTrackCuts.cfgMinNITSCls)); } - return ((track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgNTPCXrows) && (track.tpcNClsFound() >= cfgTrackCuts.cfgNTPCCls) && (track.itsNCls() >= cfgTrackCuts.cfgMinNITSCls)); } template @@ -1278,67 +1324,57 @@ struct FlowGfwV02 { } } + template + void fillPidDedxQA(TTrack const& track, const int& pid) + { + if constexpr (requires { track.tpcSignal(); track.tpcExpSigmaPi(); track.tpcExpSigmaKa(); track.tpcExpSigmaPr(); }) { + if (!cfgGetdEdx) { + return; + } + if (ft == kBefore || pid == PidPions) { + const double tpcExpSignalPi = track.tpcSignal() - (track.tpcNSigmaPi() * track.tpcExpSigmaPi()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_pions"), track.pt(), track.tpcSignal(), track.tofNSigmaPi()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_pions"), track.pt(), tpcExpSignalPi, track.tofNSigmaPi()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_pions"), track.pt(), track.tpcExpSigmaPi(), track.tofNSigmaPi()); + } + if (ft == kBefore || pid == PidKaons) { + const double tpcExpSignalKa = track.tpcSignal() - (track.tpcNSigmaKa() * track.tpcExpSigmaKa()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_kaons"), track.pt(), track.tpcSignal(), track.tofNSigmaKa()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_kaons"), track.pt(), tpcExpSignalKa, track.tofNSigmaKa()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_kaons"), track.pt(), track.tpcExpSigmaKa(), track.tofNSigmaKa()); + } + if (ft == kBefore || pid == PidProtons) { + const double tpcExpSignalPr = track.tpcSignal() - (track.tpcNSigmaPr() * track.tpcExpSigmaPr()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_protons"), track.pt(), track.tpcSignal(), track.tofNSigmaPr()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_protons"), track.pt(), tpcExpSignalPr, track.tofNSigmaPr()); + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_protons"), track.pt(), track.tpcExpSigmaPr(), track.tofNSigmaPr()); + } + } + } + template inline void fillPidQA(TTrack const& track, const int& pid) { - // Fill Nsigma QA + // Fill Nsigma QA. The derived tables store PID bits and have no nSigma columns. + if constexpr (requires { track.tpcNSigmaPi(); track.tofNSigmaPi(); }) { if (!cfgUseItsPID) { if (ft == kBefore) { registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_pions"), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt()); - if (cfgGetdEdx) { - double tpcExpSignalPi = track.tpcSignal() - (track.tpcNSigmaPi() * track.tpcExpSigmaPi()); - - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_pions"), track.pt(), track.tpcSignal(), track.tofNSigmaPi()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_pions"), track.pt(), tpcExpSignalPi, track.tofNSigmaPi()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_pions"), track.pt(), track.tpcExpSigmaPi(), track.tofNSigmaPi()); - } registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_kaons"), track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt()); - if (cfgGetdEdx) { - double tpcExpSignalKa = track.tpcSignal() - (track.tpcNSigmaKa() * track.tpcExpSigmaKa()); - - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_kaons"), track.pt(), track.tpcSignal(), track.tofNSigmaKa()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_kaons"), track.pt(), tpcExpSignalKa, track.tofNSigmaKa()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_kaons"), track.pt(), track.tpcExpSigmaKa(), track.tofNSigmaKa()); - } registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_protons"), track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt()); - if (cfgGetdEdx) { - double tpcExpSignalPr = track.tpcSignal() - (track.tpcNSigmaPr() * track.tpcExpSigmaPr()); - - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_protons"), track.pt(), track.tpcSignal(), track.tofNSigmaPr()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_protons"), track.pt(), tpcExpSignalPr, track.tofNSigmaPr()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_protons"), track.pt(), track.tpcExpSigmaPr(), track.tofNSigmaPr()); - } } else if (ft == kAfter) { if (pid == PidPions) { registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_pions"), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt()); - if (cfgGetdEdx) { - double tpcExpSignalPi = track.tpcSignal() - (track.tpcNSigmaPi() * track.tpcExpSigmaPi()); - - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_pions"), track.pt(), track.tpcSignal(), track.tofNSigmaPi()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_pions"), track.pt(), tpcExpSignalPi, track.tofNSigmaPi()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_pions"), track.pt(), track.tpcExpSigmaPi(), track.tofNSigmaPi()); - } } if (pid == PidKaons) { registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_kaons"), track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt()); - if (cfgGetdEdx) { - double tpcExpSignalKa = track.tpcSignal() - (track.tpcNSigmaKa() * track.tpcExpSigmaKa()); - - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_kaons"), track.pt(), track.tpcSignal(), track.tofNSigmaKa()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_kaons"), track.pt(), tpcExpSignalKa, track.tofNSigmaKa()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_kaons"), track.pt(), track.tpcExpSigmaKa(), track.tofNSigmaKa()); - } } if (pid == PidProtons) { registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_protons"), track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt()); - if (cfgGetdEdx) { - double tpcExpSignalPr = track.tpcSignal() - (track.tpcNSigmaPr() * track.tpcExpSigmaPr()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_protons"), track.pt(), track.tpcSignal(), track.tofNSigmaPr()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_protons"), track.pt(), tpcExpSignalPr, track.tofNSigmaPr()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_protons"), track.pt(), track.tpcExpSigmaPr(), track.tofNSigmaPr()); - } } } + fillPidDedxQA(track, pid); + } } } @@ -1348,13 +1384,17 @@ struct FlowGfwV02 { double wacc = getAcceptance(track, vtxz); registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("phi_eta_vtxZ"), track.phi(), track.eta(), vtxz, (ft == kAfter) ? wacc : 1.0); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_dcaXY_dcaZ"), track.pt(), track.dcaXY(), track.dcaZ()); + if constexpr (requires { track.dcaXY(); track.dcaZ(); }) { + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_dcaXY_dcaZ"), track.pt(), track.dcaXY(), track.dcaZ()); + } - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prTPCcls"), track.tpcChi2NCl()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prITScls"), track.itsChi2NCl()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCClusters"), track.tpcNClsFound()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nITSClusters"), track.itsNCls()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCCrossedRows"), track.tpcNClsCrossedRows()); + if constexpr (requires { track.tpcChi2NCl(); track.itsChi2NCl(); track.tpcNClsFound(); track.itsNCls(); track.tpcNClsCrossedRows(); }) { + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prTPCcls"), track.tpcChi2NCl()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prITScls"), track.itsChi2NCl()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCClusters"), track.tpcNClsFound()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nITSClusters"), track.itsNCls()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCCrossedRows"), track.tpcNClsCrossedRows()); + } if (ft == kAfter) { registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_ref"), track.pt()); @@ -1447,6 +1487,30 @@ struct FlowGfwV02 { // processCollision(collision, tracks, xaxis); } PROCESS_SWITCH(FlowGfwV02, processCFDerivedCorrected, "Process analysis for CF derived data with corrections", false); + + void processDerived(aod::GfwCollision const& collision, soa::Filtered const& tracks) + { + int run = collision.runNumber(); + if (run != lastRun) { + lastRun = run; + LOGF(info, "run = %d", run); + } + loadCorrections(run); + + const uint16_t eventMask = cfgEventBitMask; + if ((collision.eventSel() & eventMask) != eventMask) { + return; + } + if (std::abs(collision.posZ()) > cfgEventCuts.cfgZvtxMax) { + return; + } + + const XAxis xaxis{.centrality = static_cast(cfgCentralityFactor * collision.centrality()), .multiplicity = static_cast(collision.multiplicity())}; + registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); + registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); + processCollision(collision, tracks, xaxis); + } + PROCESS_SWITCH(FlowGfwV02, processDerived, "Process analysis for GFW derived data", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From aab791a7891e452a58a9d43efbcf171be2fb2d93 Mon Sep 17 00:00:00 2001 From: Maxim Virta Date: Wed, 28 Oct 2026 16:39:58 +0100 Subject: [PATCH 2/3] [PWGCF] Derived data for Generic Framework analyses created --- PWGCF/GenericFramework/DataModel/GfwDerived.h | 145 ++++++ .../TableProducer/flowGfwDerivedProducer.cxx | 443 ++++++++++++++++++ PWGCF/GenericFramework/Tasks/flowGfwV02.cxx | 128 +++-- PWGCF/JCorran/Tasks/jflucWeightsLoader.cxx | 19 +- 4 files changed, 689 insertions(+), 46 deletions(-) create mode 100644 PWGCF/GenericFramework/DataModel/GfwDerived.h create mode 100644 PWGCF/GenericFramework/TableProducer/flowGfwDerivedProducer.cxx diff --git a/PWGCF/GenericFramework/DataModel/GfwDerived.h b/PWGCF/GenericFramework/DataModel/GfwDerived.h new file mode 100644 index 00000000000..5e490aa4e10 --- /dev/null +++ b/PWGCF/GenericFramework/DataModel/GfwDerived.h @@ -0,0 +1,145 @@ +// 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. + +/// \file GfwDerived.h +/// \brief Derived collision and track tables for the GFW v02 flow analysis. +/// +/// Selection columns are bit masks. A later analysis keeps a row when every bit +/// set in its mask is also set on the row (`ncheckbit` / `(value & mask) == mask`). +/// PID is stored only as bits. The nominal, tight, and loose windows are fixed +/// at production. +/// +/// \author Maxim Virta, AU, maxim.virta@cern.ch + +#ifndef PWGCF_GENERICFRAMEWORK_DATAMODEL_GFWDERIVED_H_ +#define PWGCF_GENERICFRAMEWORK_DATAMODEL_GFWDERIVED_H_ + +#include +#include +#include + +#include + +namespace o2::aod +{ +namespace gfwcollision +{ +DECLARE_SOA_COLUMN(Centrality, centrality, float); //! Centrality estimator chosen at production +DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); //! Global-track multiplicity, multNTracksGlobal +DECLARE_SOA_COLUMN(EventSel, eventSel, uint16_t); //! Event-selection bitmask, see EventSelection + +/// One bit per event-level decision. Bits record that the decision passed. +/// sel8 is required by the producer and is not stored here. +enum EventSelection : uint16_t { + kNoSameBunchPileup = BIT(0), + kIsGoodZvtxFT0vsPV = BIT(1), + kNoCollInTimeRangeStandard = BIT(2), + kNoCollInRofStandard = BIT(3), + kNoHighMultCollInPrevRof = BIT(4), + kIsVertexITSTPC = BIT(5), + kIsGoodITSLayersAll = BIT(6), + kNoTimeFrameBorder = BIT(7), + kNoITSROFrameBorder = BIT(8), + kOccupancy = BIT(9), //! track occupancy is within the production threshold + kNotTVXinTRD = BIT(10), //! collision is not a TRD-triggered TVX event + kGoodVertex = BIT(11), //! PV contributor and z-resolution quality + kMultCorr = BIT(12) //! multiplicity-correlation cut +}; + +inline constexpr uint16_t DefaultEventMask = 0; +} // namespace gfwcollision + +DECLARE_SOA_TABLE(GfwCollisions, "AOD", "GFWCOLLISION", //! + o2::soa::Index<>, + bc::RunNumber, + collision::PosZ, + gfwcollision::Centrality, + gfwcollision::Multiplicity, + timestamp::Timestamp, + gfwcollision::EventSel); +using GfwCollision = GfwCollisions::iterator; + +namespace gfwtrack +{ +DECLARE_SOA_INDEX_COLUMN(GfwCollision, gfwCollision); //! Parent derived collision +DECLARE_SOA_COLUMN(Pt, pt, float); //! pT (GeV/c) +DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity +DECLARE_SOA_COLUMN(Phi, phi, float); //! Azimuthal angle +DECLARE_SOA_COLUMN(HasTOF, hasTOF, bool); //! Track has a TOF match +DECLARE_SOA_COLUMN(TrackSel, trackSel, uint32_t); //! Track and PID bitmask, see TrackSelection + +/// Track-quality bits occupy 0-13. TPC PID bits occupy 14-22 and TOF PID bits occupy 23-31. +/// Global-track, ITS-cluster, and ITS-chi2 requirements are applied by the producer and are not stored. +/// DCA bits are |DCA| <= nSigma times a pT-dependent shape. +/// A PID bit is set when low * scale < nSigma < up * scale. +/// Above the TOF pT cut, every PID bit stays unset when the track has no TOF match. +/// Within each detector, PidWindowBits covers π/K/p, followed by the tight window and then the loose window. +enum TrackSelection : uint32_t { + kDCAxy = BIT(0), //! |DCAxy| within the nominal nSigma times f_xy(pT) + kDCAxyTight = BIT(1), + kDCAxyLoose = BIT(2), + kDCAz = BIT(3), //! |DCAz| within the nominal nSigma times f_z(pT) + kDCAzTight = BIT(4), + kDCAzLoose = BIT(5), // DCA bit selections + + kTPCXrows = BIT(6), + kTPCXrowsTight = BIT(7), + kTPCXrowsLoose = BIT(8), + kTPCCls = BIT(9), + kTPCClsTight = BIT(10), + kTPCClsLoose = BIT(11), + kChi2TPC = BIT(12), + kChi2TPCLoose = BIT(13), // TPC bit selections + + // TPC PID selections: nominal, tight, then loose + kTpcPi = BIT(14), + kTpcKa = BIT(15), + kTpcPr = BIT(16), + kTpcPiTight = BIT(17), + kTpcKaTight = BIT(18), + kTpcPrTight = BIT(19), + kTpcPiLoose = BIT(20), + kTpcKaLoose = BIT(21), + kTpcPrLoose = BIT(22), + // TOF PID selections: nominal, tight, then loose + kTofPi = BIT(23), + kTofKa = BIT(24), + kTofPr = BIT(25), + kTofPiTight = BIT(26), + kTofKaTight = BIT(27), + kTofPrTight = BIT(28), + kTofPiLoose = BIT(29), + kTofKaLoose = BIT(30), + kTofPrLoose = BIT(31) +}; + +/// Quality bits corresponding to the default flowGfwV02 track filter plus trackSelected(). +inline constexpr uint32_t NominalMask = static_cast(kDCAxy | kDCAz | kTPCXrows | kTPCCls | kChi2TPC); +/// Loosest quality bits. Production can optionally require these while still leaving the nominal cuts free. +inline constexpr uint32_t LooseSkimMask = static_cast(kDCAxyLoose | kDCAzLoose | kTPCXrowsLoose | kTPCClsLoose | kChi2TPCLoose); +/// Tracks that enter the multiplicity-correlation count. Global track and ITS chi2 are already required by the producer. +inline constexpr uint32_t MultCorrTrackMask = static_cast(kDCAz | kChi2TPC); + +inline constexpr int PidWindowBits = 3; +inline constexpr int PidWindowNominal = 0; +inline constexpr int PidWindowTight = 1; +inline constexpr int PidWindowLoose = 2; +} // namespace gfwtrack + +DECLARE_SOA_TABLE(GfwTracks, "AOD", "GFWTRACK", //! + o2::soa::Index<>, + gfwtrack::GfwCollisionId, + gfwtrack::Pt, gfwtrack::Eta, gfwtrack::Phi, + gfwtrack::HasTOF, gfwtrack::TrackSel); +using GfwTrack = GfwTracks::iterator; +} // namespace o2::aod + +#endif // PWGCF_GENERICFRAMEWORK_DATAMODEL_GFWDERIVED_H_ diff --git a/PWGCF/GenericFramework/TableProducer/flowGfwDerivedProducer.cxx b/PWGCF/GenericFramework/TableProducer/flowGfwDerivedProducer.cxx new file mode 100644 index 00000000000..1bfed843783 --- /dev/null +++ b/PWGCF/GenericFramework/TableProducer/flowGfwDerivedProducer.cxx @@ -0,0 +1,443 @@ +// 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. + +/// \file flowGfwDerivedProducer.cxx +/// \brief Produce the GFW derived dataset used by flowGfwV02. +/// +/// Keeps sel8 collisions and global tracks that pass the ITS cluster and ITS chi2 +/// cuts. Stores the event mask and the combined track and PID bitmask. PID bits +/// stay unset above the TOF pT cut when the track has no TOF match. Run it after +/// event selection, track selection, multiplicity, centrality, and the TPC/TOF PID +/// helpers (pid-tpc, pid-tof-full, pid-tof-beta). +/// +/// \author Maxim Virta, NBI, maxim.virta@cern.ch + +#include "PWGCF/GenericFramework/DataModel/GfwDerived.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/CCDB/TriggerAliases.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace +{ +constexpr float VertexResMax = 0.25f; +constexpr int VertexContribMin = 20; +constexpr float DcaShapeMinPt = 0.001f; +constexpr float DcaShapeMaxPt = 100.f; +constexpr float DcaZReferenceNsigma = 7.f; +constexpr float NsigmaScaleNominal = 1.f; +constexpr float NsigmaScaleTight = 2.f / 3.f; +constexpr float NsigmaScaleLoose = 4.f / 3.f; +constexpr uint32_t NsigmaRows = 6; +constexpr uint32_t NsigmaCols = 3; +constexpr int NsigmaUpPi = 0; +constexpr int NsigmaUpKa = 1; +constexpr int NsigmaUpPr = 2; +constexpr int NsigmaLowPi = 3; +constexpr int NsigmaLowKa = 4; +constexpr int NsigmaLowPr = 5; +constexpr int NsigmaColTpc = 0; +constexpr int NsigmaColTof = 1; +constexpr std::array, NsigmaRows> NsigmaDefault = {{{{1.1f, 2.1f, 3.1f}}, {{1.2f, 2.2f, 3.2f}}, {{1.3f, 2.3f, 3.3f}}, {{-1.1f, -2.1f, -3.1f}}, {{-1.2f, -2.2f, -3.2f}}, {{-1.3f, -2.3f, -3.3f}}}}; +} // namespace + +struct FlowGfwDerivedProducer { + using CollisionsWithInfo = soa::Join; + using TracksWithPid = soa::Join; + + enum CentEstimator : int { + CentFT0C = 0, + CentFT0CVariant1, + CentFT0M, + CentFV0A, + CentNTPV, + CentNGlobal, + CentMFT + }; + + Configurable cfgCutVertex{"cfgCutVertex", 10.f, "Maximum |z| of collisions written to the derived dataset (cm)"}; + Configurable cfgPtMin{"cfgPtMin", 0.2f, "Minimum track pT written to the derived dataset (GeV/c)"}; + Configurable cfgPtMax{"cfgPtMax", 10.f, "Maximum track pT written to the derived dataset (GeV/c)"}; + Configurable cfgEtaMax{"cfgEtaMax", 0.8f, "Maximum |eta| of tracks written to the derived dataset"}; + Configurable cfgSkipEmptyEvents{"cfgSkipEmptyEvents", true, "Do not write collisions that keep no tracks"}; + Configurable cfgRequireLooseTrackBits{"cfgRequireLooseTrackBits", false, "Keep only tracks that pass the loose DCA and TPC-crossed-rows bits"}; + Configurable cfgCentEstimator{"cfgCentEstimator", 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FV0A; 4:NTPV; 5:NGlobal; 6:MFT"}; + Configurable cfgOccupancyMax{"cfgOccupancyMax", 2000, "Maximum track occupancy for the occupancy bit; negative disables the cut"}; + + Configurable cfgDCAxyNSigma{"cfgDCAxyNSigma", 7.f, "Nominal pT-dependent DCAxy nSigma. Non-positive disables the nominal bit"}; + Configurable cfgDCAxyNSigmaTight{"cfgDCAxyNSigmaTight", 5.f, "Tight pT-dependent DCAxy nSigma"}; + Configurable cfgDCAxyNSigmaLoose{"cfgDCAxyNSigmaLoose", 10.f, "Loose pT-dependent DCAxy nSigma"}; + Configurable cfgDCAxy{"cfgDCAxy", "(0.0026+0.005/(x^1.01))", "1-sigma shape of the pT-dependent DCAxy cut; a bit passes when |DCAxy| <= nSigma*shape(pT)"}; + Configurable cfgDCAzNSigma{"cfgDCAzNSigma", 7.f, "Nominal DCAz nSigma relative to the 7-sigma pT shape. Non-positive disables the nominal bit"}; + Configurable cfgDCAzNSigmaTight{"cfgDCAzNSigmaTight", 5.f, "Tight DCAz nSigma"}; + Configurable cfgDCAzNSigmaLoose{"cfgDCAzNSigmaLoose", 10.f, "Loose DCAz nSigma"}; + Configurable cfgDCAz{"cfgDCAz", "(0.0105+0.0350/(x^1.1))", "7-sigma pT-dependent DCAz shape; a bit passes when |DCAz| <= (nSigma/7)*shape(pT)"}; + Configurable cfgNTPCXrows{"cfgNTPCXrows", 70.f, "Nominal minimum number of TPC crossed rows"}; + Configurable cfgNTPCXrowsTight{"cfgNTPCXrowsTight", 80.f, "Tight minimum number of TPC crossed rows"}; + Configurable cfgNTPCXrowsLoose{"cfgNTPCXrowsLoose", 60.f, "Loose minimum number of TPC crossed rows"}; + Configurable cfgNTPCCls{"cfgNTPCCls", 50.f, "Nominal minimum number of found TPC clusters"}; + Configurable cfgNTPCClsTight{"cfgNTPCClsTight", 60.f, "Tight minimum number of found TPC clusters"}; + Configurable cfgNTPCClsLoose{"cfgNTPCClsLoose", 40.f, "Loose minimum number of found TPC clusters"}; + Configurable cfgMinNITSCls{"cfgMinNITSCls", 5.f, "Minimum number of ITS clusters required to keep a track"}; + Configurable cfgChi2PrTPCCls{"cfgChi2PrTPCCls", 2.5f, "Nominal maximum chi2 per TPC cluster"}; + Configurable cfgChi2PrTPCClsLoose{"cfgChi2PrTPCClsLoose", 4.f, "Loose maximum chi2 per TPC cluster"}; + Configurable cfgChi2PrITSCls{"cfgChi2PrITSCls", 36.f, "Maximum chi2 per ITS cluster required to keep a track"}; + + Configurable> cfgMultGlobalCutPars{"cfgMultGlobalCutPars", std::vector{2272.16, -76.6932, 1.01204, -0.00631545, 1.59868e-05, 136.336, -4.97006, 0.121199, -0.0015921, 7.66197e-06}, "Global vs FT0C multiplicity cut parameters"}; + Configurable> cfgMultPVCutPars{"cfgMultPVCutPars", std::vector{3074.43, -106.192, 1.46176, -0.00968364, 2.61923e-05, 182.128, -7.43492, 0.193901, -0.00256715, 1.22594e-05}, "PV vs FT0C multiplicity cut parameters"}; + Configurable> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector{-0.223013, 0.715849, 0.664242, 0.0829653, -0.000503733, 1.21185e-06}, "Global vs PV multiplicity cut parameters"}; + Configurable cfgMultCorrHighCutFunction{"cfgMultCorrHighCutFunction", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Upper multiplicity-correlation cut"}; + Configurable cfgMultCorrLowCutFunction{"cfgMultCorrLowCutFunction", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x - 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Lower multiplicity-correlation cut"}; + Configurable cfgMultGlobalPVCorrCutFunction{"cfgMultGlobalPVCorrCutFunction", "[0] + [1]*x + 3*([2] + [3]*x + [4]*x*x + [5]*x*x*x)", "Global vs PV multiplicity-correlation cut"}; + Configurable cfgTofPtCut{"cfgTofPtCut", 0.5f, "Above this pT, PID bits stay unset when the track has no TOF match"}; + Configurable> nSigmas{"nSigmas", {NsigmaDefault.front().data(), NsigmaRows, NsigmaCols, {"UpCut_pi", "UpCut_ka", "UpCut_pr", "LowCut_pi", "LowCut_ka", "LowCut_pr"}, {"TPC", "TOF", "ITS"}}, "Nominal nSigma windows for TPC, TOF, and ITS. ITS is accepted so the analysis JSON can be reused and is not stored as a bit"}; + Configurable cfgNsigmaScaleTight{"cfgNsigmaScaleTight", static_cast(NsigmaScaleTight), "Scale applied to the nominal nSigma bounds for the tight PID bits"}; + Configurable cfgNsigmaScaleLoose{"cfgNsigmaScaleLoose", static_cast(NsigmaScaleLoose), "Scale applied to the nominal nSigma bounds for the loose PID bits"}; + + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter trackFilter = (aod::track::pt > cfgPtMin) && (aod::track::pt < cfgPtMax) && (nabs(aod::track::eta) < cfgEtaMax); + + Produces outputCollisions; + Produces outputTracks; + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; + std::unique_ptr fDcaXyShape; + std::unique_ptr fDcaZShape; + std::array tpcNsigmaCut{}; + std::array tofNsigmaCut{}; + std::unique_ptr fMultPVCutLow; + std::unique_ptr fMultPVCutHigh; + std::unique_ptr fMultCutLow; + std::unique_ptr fMultCutHigh; + std::unique_ptr fMultPVGlobalCutHigh; + + void init(InitContext const&) + { + const AxisSpec axisCounter{1, 0., 1., ""}; + registry.add("events", "Written derived collisions", HistType::kTH1F, {axisCounter}); + registry.add("tracks", "Written derived tracks", HistType::kTH1F, {axisCounter}); + + const std::string dcaXyFormula = "[0]*" + cfgDCAxy.value; + fDcaXyShape = std::make_unique("gfwDerivedDcaXy", dcaXyFormula.c_str(), DcaShapeMinPt, DcaShapeMaxPt); + fDcaXyShape->SetParameter(0, 1.f); + fDcaZShape = std::make_unique("gfwDerivedDcaZ", cfgDCAz.value.c_str(), DcaShapeMinPt, DcaShapeMaxPt); + + const auto nsigma = nSigmas->getData(); + for (uint32_t row = 0; row < NsigmaRows; ++row) { + tpcNsigmaCut[row] = nsigma(row, NsigmaColTpc); + tofNsigmaCut[row] = nsigma(row, NsigmaColTof); + } + + constexpr double MultCorrMaxCent = 100.; + constexpr double MultCorrMaxNch = 3000.; + fMultPVCutLow = std::make_unique("gfwDerivedMultPVCutLow", cfgMultCorrLowCutFunction.value.c_str(), 0., MultCorrMaxCent); + fMultPVCutLow->SetParameters(cfgMultPVCutPars.value.data()); + fMultPVCutHigh = std::make_unique("gfwDerivedMultPVCutHigh", cfgMultCorrHighCutFunction.value.c_str(), 0., MultCorrMaxCent); + fMultPVCutHigh->SetParameters(cfgMultPVCutPars.value.data()); + fMultCutLow = std::make_unique("gfwDerivedMultCutLow", cfgMultCorrLowCutFunction.value.c_str(), 0., MultCorrMaxCent); + fMultCutLow->SetParameters(cfgMultGlobalCutPars.value.data()); + fMultCutHigh = std::make_unique("gfwDerivedMultCutHigh", cfgMultCorrHighCutFunction.value.c_str(), 0., MultCorrMaxCent); + fMultCutHigh->SetParameters(cfgMultGlobalCutPars.value.data()); + fMultPVGlobalCutHigh = std::make_unique("gfwDerivedMultPVGlobalCutHigh", cfgMultGlobalPVCorrCutFunction.value.c_str(), 0., MultCorrMaxNch); + fMultPVGlobalCutHigh->SetParameters(cfgMultGlobalPVCutPars.value.data()); + } + + template + float centrality(TCollision const& collision) const + { + switch (cfgCentEstimator.value) { + case CentFT0C: + return collision.centFT0C(); + case CentFT0CVariant1: + return collision.centFT0CVariant1(); + case CentFT0M: + return collision.centFT0M(); + case CentFV0A: + return collision.centFV0A(); + case CentNTPV: + return collision.centNTPV(); + case CentNGlobal: + return collision.centNGlobal(); + case CentMFT: + return collision.centMFT(); + default: + return collision.centFT0C(); + } + } + + template + bool vertexQuality(TCollision const& collision) const + { + if (collision.numContrib() <= 1) { + return false; + } + const float zRes = std::sqrt(collision.covZZ()); + return !(zRes > VertexResMax && collision.numContrib() < VertexContribMin); + } + + template + bool multiplicityCorrelated(TCollision const& collision, float cent, float multTrk) const + { + const float multPV = collision.multNTracksPV(); + if (multPV < fMultPVCutLow->Eval(cent) || multPV > fMultPVCutHigh->Eval(cent)) { + return false; + } + if (multTrk < fMultCutLow->Eval(cent) || multTrk > fMultCutHigh->Eval(cent)) { + return false; + } + return multTrk <= fMultPVGlobalCutHigh->Eval(multPV); + } + + template + uint16_t packEventSel(TCollision const& collision, float cent, float multTrk) const + { + uint16_t bits = 0; + if (collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + bits |= aod::gfwcollision::kNoSameBunchPileup; + } + if (collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + bits |= aod::gfwcollision::kIsGoodZvtxFT0vsPV; + } + if (collision.selection_bit(aod::evsel::kNoCollInTimeRangeStandard)) { + bits |= aod::gfwcollision::kNoCollInTimeRangeStandard; + } + if (collision.selection_bit(aod::evsel::kNoCollInRofStandard)) { + bits |= aod::gfwcollision::kNoCollInRofStandard; + } + if (collision.selection_bit(aod::evsel::kNoHighMultCollInPrevRof)) { + bits |= aod::gfwcollision::kNoHighMultCollInPrevRof; + } + if (collision.selection_bit(aod::evsel::kIsVertexITSTPC)) { + bits |= aod::gfwcollision::kIsVertexITSTPC; + } + if (collision.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { + bits |= aod::gfwcollision::kIsGoodITSLayersAll; + } + if (collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + bits |= aod::gfwcollision::kNoTimeFrameBorder; + } + if (collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { + bits |= aod::gfwcollision::kNoITSROFrameBorder; + } + if (cfgOccupancyMax.value < 0 || collision.trackOccupancyInTimeRange() <= cfgOccupancyMax.value) { + bits |= aod::gfwcollision::kOccupancy; + } + if (!collision.alias_bit(kTVXinTRD)) { + bits |= aod::gfwcollision::kNotTVXinTRD; + } + if (vertexQuality(collision)) { + bits |= aod::gfwcollision::kGoodVertex; + } + if (multiplicityCorrelated(collision, cent, multTrk)) { + bits |= aod::gfwcollision::kMultCorr; + } + return bits; + } + + template + bool withinDcaXy(TTrack const& track, float nSigma) const + { + return std::abs(track.dcaXY()) <= nSigma * fDcaXyShape->Eval(track.pt()); + } + + template + bool withinDcaZ(TTrack const& track, float nSigma) const + { + return std::abs(track.dcaZ()) <= (nSigma / DcaZReferenceNsigma) * fDcaZShape->Eval(track.pt()); + } + + template + bool passProducerTrackCuts(TTrack const& track) const + { + const bool global = track.isGlobalTrack() || track.isGlobalTrackSDD() != 0; + return global && track.itsNCls() >= cfgMinNITSCls.value && track.itsChi2NCl() <= cfgChi2PrITSCls.value; + } + + template + uint32_t packTrackSel(TTrack const& track) const + { + uint32_t bits = 0; + if (cfgDCAxyNSigma.value <= 0.f || withinDcaXy(track, cfgDCAxyNSigma.value)) { + bits |= aod::gfwtrack::kDCAxy; + } + if (withinDcaXy(track, cfgDCAxyNSigmaTight.value)) { + bits |= aod::gfwtrack::kDCAxyTight; + } + if (withinDcaXy(track, cfgDCAxyNSigmaLoose.value)) { + bits |= aod::gfwtrack::kDCAxyLoose; + } + if (cfgDCAzNSigma.value <= 0.f || withinDcaZ(track, cfgDCAzNSigma.value)) { + bits |= aod::gfwtrack::kDCAz; + } + if (withinDcaZ(track, cfgDCAzNSigmaTight.value)) { + bits |= aod::gfwtrack::kDCAzTight; + } + if (withinDcaZ(track, cfgDCAzNSigmaLoose.value)) { + bits |= aod::gfwtrack::kDCAzLoose; + } + if (track.tpcNClsCrossedRows() >= cfgNTPCXrows.value) { + bits |= aod::gfwtrack::kTPCXrows; + } + if (track.tpcNClsCrossedRows() >= cfgNTPCXrowsTight.value) { + bits |= aod::gfwtrack::kTPCXrowsTight; + } + if (track.tpcNClsCrossedRows() >= cfgNTPCXrowsLoose.value) { + bits |= aod::gfwtrack::kTPCXrowsLoose; + } + if (track.tpcNClsFound() >= cfgNTPCCls.value) { + bits |= aod::gfwtrack::kTPCCls; + } + if (track.tpcNClsFound() >= cfgNTPCClsTight.value) { + bits |= aod::gfwtrack::kTPCClsTight; + } + if (track.tpcNClsFound() >= cfgNTPCClsLoose.value) { + bits |= aod::gfwtrack::kTPCClsLoose; + } + if (track.tpcChi2NCl() <= cfgChi2PrTPCCls.value) { + bits |= aod::gfwtrack::kChi2TPC; + } + if (track.tpcChi2NCl() <= cfgChi2PrTPCClsLoose.value) { + bits |= aod::gfwtrack::kChi2TPCLoose; + } + return bits; + } + + static bool insideNsigma(float nsigma, float up, float low, float scale) + { + return nsigma < up * scale && nsigma > low * scale; + } + + template + uint32_t packPidWindow(TTrack const& track, float scale) const + { + if (track.pt() > cfgTofPtCut.value && !track.hasTOF()) { + return 0; + } + uint32_t bits = 0; + if (insideNsigma(track.tpcNSigmaPi(), tpcNsigmaCut[NsigmaUpPi], tpcNsigmaCut[NsigmaLowPi], scale)) { + bits |= aod::gfwtrack::kTpcPi; + } + if (insideNsigma(track.tpcNSigmaKa(), tpcNsigmaCut[NsigmaUpKa], tpcNsigmaCut[NsigmaLowKa], scale)) { + bits |= aod::gfwtrack::kTpcKa; + } + if (insideNsigma(track.tpcNSigmaPr(), tpcNsigmaCut[NsigmaUpPr], tpcNsigmaCut[NsigmaLowPr], scale)) { + bits |= aod::gfwtrack::kTpcPr; + } + if (insideNsigma(track.tofNSigmaPi(), tofNsigmaCut[NsigmaUpPi], tofNsigmaCut[NsigmaLowPi], scale)) { + bits |= aod::gfwtrack::kTofPi; + } + if (insideNsigma(track.tofNSigmaKa(), tofNsigmaCut[NsigmaUpKa], tofNsigmaCut[NsigmaLowKa], scale)) { + bits |= aod::gfwtrack::kTofKa; + } + if (insideNsigma(track.tofNSigmaPr(), tofNsigmaCut[NsigmaUpPr], tofNsigmaCut[NsigmaLowPr], scale)) { + bits |= aod::gfwtrack::kTofPr; + } + return bits; + } + + template + uint32_t packPidSel(TTrack const& track) const + { + const uint32_t nominal = packPidWindow(track, NsigmaScaleNominal); + const uint32_t tight = packPidWindow(track, cfgNsigmaScaleTight.value); + const uint32_t loose = packPidWindow(track, cfgNsigmaScaleLoose.value); + // Nominal TPC bits start at bit 14 and nominal TOF bits at bit 23. Shifting by one window moves each detector into its tight or loose block. + return nominal | (tight << (aod::gfwtrack::PidWindowTight * aod::gfwtrack::PidWindowBits)) | (loose << (aod::gfwtrack::PidWindowLoose * aod::gfwtrack::PidWindowBits)); + } + + template + void fillTrack(TTrack const& track, uint32_t bits) + { + outputTracks(outputCollisions.lastIndex(), track.pt(), track.eta(), track.phi(), + track.hasTOF(), bits | packPidSel(track)); + registry.fill(HIST("tracks"), 0.5); + } + + void process(soa::Filtered::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered const& tracks) + { + if (!collision.sel8()) { + return; + } + + int nKept = 0; + int nForMultCorr = 0; + for (const auto& track : tracks) { + if (!passProducerTrackCuts(track)) { + continue; + } + const uint32_t bits = packTrackSel(track); + const bool passLoose = !cfgRequireLooseTrackBits || ((bits & aod::gfwtrack::LooseSkimMask) == aod::gfwtrack::LooseSkimMask); + if (!passLoose) { + continue; + } + ++nKept; + if ((bits & aod::gfwtrack::MultCorrTrackMask) == aod::gfwtrack::MultCorrTrackMask) { + ++nForMultCorr; + } + } + if (nKept == 0 && cfgSkipEmptyEvents) { + return; + } + + const auto bc = collision.bc_as(); + const float cent = centrality(collision); + outputCollisions(bc.runNumber(), collision.posZ(), cent, static_cast(collision.multNTracksGlobal()), bc.timestamp(), packEventSel(collision, cent, static_cast(nForMultCorr))); + registry.fill(HIST("events"), 0.5); + + for (const auto& track : tracks) { + if (!passProducerTrackCuts(track)) { + continue; + } + const uint32_t bits = packTrackSel(track); + const bool passLoose = !cfgRequireLooseTrackBits || ((bits & aod::gfwtrack::LooseSkimMask) == aod::gfwtrack::LooseSkimMask); + if (!passLoose) { + continue; + } + fillTrack(track, bits); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + }; +} diff --git a/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx b/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx index ad8dcfef573..2829e2e2316 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx @@ -215,9 +215,9 @@ struct FlowGfwV02 { // // The analysis assumes the data has been subjected to a QA of its selection, // // and thus only the final distributions of the data for analysis are saved. // Derived-data bit masks. Every listed track bit is required; see gfwcollision::EventSel and gfwtrack::TrackSelection. - O2_DEFINE_CONFIGURABLE(cfgEventBitMask, uint16_t, o2::aod::gfwcollision::DefaultEventMask, "Required event-selection bits for GFW derived data"); + O2_DEFINE_CONFIGURABLE(cfgEventBitMask, uint16_t, static_cast(o2::aod::gfwcollision::DefaultEventMask), "Required event-selection bits for GFW derived data"); Configurable> cfgBitmaskSelections{"cfgBitmaskSelections", {0, 3, 6, 9, 12}, "Enabled TrackSelection bit indices for GFW derived data. DCA 0-5, TPC track 6-13, TPC PID 14-22, TOF PID 23-31. The default is the nominal quality selection"}; - O2_DEFINE_CONFIGURABLE(cfgPidWindow, int, o2::aod::gfwtrack::PidWindowNominal, "PID nSigma window for GFW derived data: 0 nominal, 1 tight, 2 loose"); + O2_DEFINE_CONFIGURABLE(cfgPidWindow, int, static_cast(o2::aod::gfwtrack::PidWindowNominal), "PID nSigma window for GFW derived data: 0 nominal, 1 tight, 2 loose"); o2::framework::expressions::Filter collFilter = (nabs(aod::collision::posZ) < cfgEventCuts.cfgZvtxMax); o2::framework::expressions::Filter trackFilter = (aod::track::pt > cfgTrackCuts.cfgPtMin) && (aod::track::pt < cfgTrackCuts.cfgPtMax) && nabs(aod::track::eta) < cfgTrackCuts.cfgEtaMax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && (aod::track::itsChi2NCl < cfgTrackCuts.cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgTrackCuts.cfgDCAz; @@ -637,7 +637,7 @@ struct FlowGfwV02 { isTofPion = (trackSel & (static_cast(aod::gfwtrack::kTofPi) << shift)) != 0; isTofKaon = (trackSel & (static_cast(aod::gfwtrack::kTofKa) << shift)) != 0; isTofProton = (trackSel & (static_cast(aod::gfwtrack::kTofPr) << shift)) != 0; - } else { + } else if constexpr (requires { track.tpcNSigmaPi(); track.tofNSigmaPi(); }) { const std::array nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; const std::array nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; constexpr bool HasItsResponse = requires { @@ -672,19 +672,25 @@ struct FlowGfwV02 { isTofPion = nSigmaTOF[IndPionUp] < pidStates.tofNsigmaCut[IndPionUp] && nSigmaTOF[IndPionUp] > pidStates.tofNsigmaCut[IndPionLow]; isTofKaon = nSigmaTOF[IndKaonUp] < pidStates.tofNsigmaCut[IndKaonUp] && nSigmaTOF[IndKaonUp] > pidStates.tofNsigmaCut[IndKaonLow]; isTofProton = nSigmaTOF[IndProtonUp] < pidStates.tofNsigmaCut[IndProtonUp] && nSigmaTOF[IndProtonUp] > pidStates.tofNsigmaCut[IndProtonLow]; - } - - if (track.pt() > cfgTofPtCut && !track.hasTOF()) { + } else { return -1; } - if (track.pt() > cfgTofPtCut && track.hasTOF()) { - isPion = isTofPion && isDetectedPion; - isKaon = isTofKaon && isDetectedKaon; - isProton = isTofProton && isDetectedProton; + + if constexpr (requires { track.hasTOF(); }) { + if (track.pt() > cfgTofPtCut && !track.hasTOF()) { + return -1; + } + if (track.pt() > cfgTofPtCut && track.hasTOF()) { + isPion = isTofPion && isDetectedPion; + isKaon = isTofKaon && isDetectedKaon; + isProton = isTofProton && isDetectedProton; + } else { + isPion = isDetectedPion; + isKaon = isDetectedKaon; + isProton = isDetectedProton; + } } else { - isPion = isDetectedPion; - isKaon = isDetectedKaon; - isProton = isDetectedProton; + return -1; } if ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton)) { @@ -748,12 +754,15 @@ struct FlowGfwV02 { cfg.correctionsLoaded = true; } + // JWeights, filled by jflucWeightsLoader from a local file or CCDB, stores the same + // correction used by flowJSPCAnalysis: acceptance is 1/weightNUA and efficiency is weightEff. template double getJTrackAcceptance(TTrack const& track) { double wacc = 1; if constexpr (requires { track.weightNUA(); }) { - wacc = 1. / track.weightNUA(); + const double nua = track.weightNUA(); + wacc = (nua == 0.) ? 1. : 1. / nua; } return wacc; } @@ -771,24 +780,36 @@ struct FlowGfwV02 { template double getAcceptance(TTrack const& track, const double& vtxz) { - double wacc = 1; - if (cfg.mAcceptance) { - wacc = cfg.mAcceptance->getNUA(track.phi(), track.eta(), vtxz); + if constexpr (requires { track.weightNUA(); }) { + return getJTrackAcceptance(track); + } else { + double wacc = 1; + if (cfg.mAcceptance) { + wacc = cfg.mAcceptance->getNUA(track.phi(), track.eta(), vtxz); + } + return wacc; } - return wacc; } template double getEfficiency(TTrack const& track, const int& pid = PidCharged) { - double eff = 1.; - if (cfg.mEfficiency[pid]) { - eff = cfg.mEfficiency[pid]->GetBinContent(cfg.mEfficiency[pid]->FindBin(track.pt())); - } - if (eff == 0) { - return -1.; + if constexpr (requires { track.weightEff(); }) { + const double eff = getJTrackEfficiency(track); + if (eff == 0.) { + return -1.; + } + return eff; + } else { + double eff = 1.; + if (cfg.mEfficiency[pid]) { + eff = cfg.mEfficiency[pid]->GetBinContent(cfg.mEfficiency[pid]->FindBin(track.pt())); + } + if (eff == 0) { + return -1.; + } + return 1. / eff; } - return 1. / eff; } template @@ -1002,21 +1023,21 @@ struct FlowGfwV02 { } } } - } - // Calculate NCv22pt - if (cfgUseMultiplicityFracWeights) { - weightC = pidStates.hPtMid[PidCharged]->Integral(); - } + // Calculate NCv22pt + if (cfgUseMultiplicityFracWeights) { + weightC = pidStates.hPtMid[PidCharged]->Integral(); + } - registry.fill(HIST("analysis/charged/v22pt"), bootstrap, centmult, v22pt, weightC * weightAB); - registry.fill(HIST("analysis/charged/v22"), bootstrap, centmult, valAB, weightAB); - registry.fill(HIST("analysis/charged/v24"), bootstrap, centmult, valABAB, weightABAB); - registry.fill(HIST("analysis/charged/ptC"), bootstrap, centmult, ptMeanMid, weightC); - if (cfgUseMultiplicityFracWeights) { - registry.fill(HIST("analysis/charged/v223pc"), bootstrap, centmult, valAB, weightC * weightAB); - registry.fill(HIST("analysis/charged/v224pc"), bootstrap, centmult, valAB, weightABAB); - registry.fill(HIST("analysis/charged/ptC3pc"), bootstrap, centmult, ptMeanMid, weightC * weightAB); + registry.fill(HIST("analysis/charged/v22pt"), bootstrap, centmult, v22pt, weightC * weightAB); + registry.fill(HIST("analysis/charged/v22"), bootstrap, centmult, valAB, weightAB); + registry.fill(HIST("analysis/charged/v24"), bootstrap, centmult, valABAB, weightABAB); + registry.fill(HIST("analysis/charged/ptC"), bootstrap, centmult, ptMeanMid, weightC); + if (cfgUseMultiplicityFracWeights) { + registry.fill(HIST("analysis/charged/v223pc"), bootstrap, centmult, valAB, weightC * weightAB); + registry.fill(HIST("analysis/charged/v224pc"), bootstrap, centmult, valAB, weightABAB); + registry.fill(HIST("analysis/charged/ptC3pc"), bootstrap, centmult, ptMeanMid, weightC * weightAB); + } } } // End of V02 @@ -1230,11 +1251,13 @@ struct FlowGfwV02 { if constexpr (requires { track.trackSel(); }) { const uint32_t mask = trackBitMask; return (track.trackSel() & mask) == mask; - } else { + } else if constexpr (requires { track.dcaXY(); track.tpcNClsCrossedRows(); track.tpcNClsFound(); track.itsNCls(); }) { if (cfgTrackCuts.cfgDCAxyNSigma && (std::fabs(track.dcaXY()) > fPtDepDCAxy->Eval(track.pt()))) { return false; } return ((track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgNTPCXrows) && (track.tpcNClsFound() >= cfgTrackCuts.cfgNTPCCls) && (track.itsNCls() >= cfgTrackCuts.cfgMinNITSCls)); + } else { + return true; } } @@ -1484,9 +1507,9 @@ struct FlowGfwV02 { const XAxis xaxis{.centrality = collision.multiplicity(), .multiplicity = tracks.size()}; registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); - // processCollision(collision, tracks, xaxis); + processCollision(collision, tracks, xaxis); } - PROCESS_SWITCH(FlowGfwV02, processCFDerivedCorrected, "Process analysis for CF derived data with corrections", false); + PROCESS_SWITCH(FlowGfwV02, processCFDerivedCorrected, "Process CF derived data with JWeights acceptance and efficiency corrections", false); void processDerived(aod::GfwCollision const& collision, soa::Filtered const& tracks) { @@ -1511,6 +1534,29 @@ struct FlowGfwV02 { processCollision(collision, tracks, xaxis); } PROCESS_SWITCH(FlowGfwV02, processDerived, "Process analysis for GFW derived data", false); + + void processDerivedCorrected(aod::GfwCollision const& collision, soa::Filtered> const& tracks) + { + int run = collision.runNumber(); + if (run != lastRun) { + lastRun = run; + LOGF(info, "run = %d", run); + } + + const uint16_t eventMask = cfgEventBitMask; + if ((collision.eventSel() & eventMask) != eventMask) { + return; + } + if (std::abs(collision.posZ()) > cfgEventCuts.cfgZvtxMax) { + return; + } + + const XAxis xaxis{.centrality = static_cast(cfgCentralityFactor * collision.centrality()), .multiplicity = static_cast(collision.multiplicity())}; + registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); + registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); + processCollision(collision, tracks, xaxis); + } + PROCESS_SWITCH(FlowGfwV02, processDerivedCorrected, "Process GFW derived data with JWeights acceptance and efficiency corrections", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/JCorran/Tasks/jflucWeightsLoader.cxx b/PWGCF/JCorran/Tasks/jflucWeightsLoader.cxx index 02088eab148..bbc5d528faf 100644 --- a/PWGCF/JCorran/Tasks/jflucWeightsLoader.cxx +++ b/PWGCF/JCorran/Tasks/jflucWeightsLoader.cxx @@ -16,6 +16,7 @@ /// The weights are loaded from the local files or CCDB and stored in the JWeights table. #include "PWGCF/DataModel/CorrelationsDerived.h" +#include "PWGCF/GenericFramework/DataModel/GfwDerived.h" #include "PWGCF/JCorran/DataModel/JCatalyst.h" #include @@ -116,14 +117,16 @@ struct JflucWeightsLoader { void init(InitContext const&) { - if (!doprocessLoadWeights && !doprocessLoadWeightsCF && !doprocessJNUACreatorCFDerived) { + if (!doprocessLoadWeights && !doprocessLoadWeightsCF && !doprocessLoadWeightsGFW && !doprocessJNUACreatorCFDerived) { return; } - if (doprocessLoadWeights && doprocessLoadWeightsCF) - LOGF(fatal, "Only one of JTracks or CFTracks processing can be enabled at a time."); + const int nWeightProcesses = static_cast(doprocessLoadWeights) + static_cast(doprocessLoadWeightsCF) + static_cast(doprocessLoadWeightsGFW); + if (nWeightProcesses > 1) { + LOGF(fatal, "Only one of JTracks, CFTracks, or GFWTracks weight loading can be enabled at a time."); + } - if (doprocessLoadWeights || doprocessLoadWeightsCF) { + if (doprocessLoadWeights || doprocessLoadWeightsCF || doprocessLoadWeightsGFW) { // NUA corrections from local file or CCDB if (cfgPathPhiWeights.value.substr(0, 8) == "local://") { LOGF(info, "Using NUA corrections locally from: %s", cfgPathPhiWeights.value.substr(8).c_str()); @@ -172,7 +175,7 @@ struct JflucWeightsLoader { ccdb->setFatalWhenNull(false); } } - } // doprocessLoadWeights || doprocessLoadWeightsCF + } // doprocessLoadWeights || doprocessLoadWeightsCF || doprocessLoadWeightsGFW if (doprocessJNUACreatorCFDerived) { qaHistRegistry.add("trackType", "trackType;trackType;counts", HistType::kTH1F, {{65, -0.5, 64.5}}); @@ -276,6 +279,12 @@ struct JflucWeightsLoader { } PROCESS_SWITCH(JflucWeightsLoader, processLoadWeightsCF, "Load weights histograms for CF derived data table", true); + void processLoadWeightsGFW(aod::GfwCollision const& collision, aod::GfwTracks const& tracks) + { + loadWeights(output, collision, tracks); + } + PROCESS_SWITCH(JflucWeightsLoader, processLoadWeightsGFW, "Load weights histograms for GFW derived data table", false); + Produces output2p; void processLoadWeightsCF2Prong(aod::CFCollision const& collision, aod::CF2ProngTracks const& tracks2p) { From 97f7d864e5c2ad2ffbbc46e504ba1a7d9dad1d7a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 28 Oct 2026 15:41:47 +0000 Subject: [PATCH 3/3] Please consider the following formatting changes --- PWGCF/GenericFramework/DataModel/GfwDerived.h | 11 ++++---- PWGCF/GenericFramework/Tasks/flowGfwV02.cxx | 28 +++++++++---------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/PWGCF/GenericFramework/DataModel/GfwDerived.h b/PWGCF/GenericFramework/DataModel/GfwDerived.h index 5e490aa4e10..9fc0edf7f5f 100644 --- a/PWGCF/GenericFramework/DataModel/GfwDerived.h +++ b/PWGCF/GenericFramework/DataModel/GfwDerived.h @@ -24,6 +24,7 @@ #include #include + #include #include @@ -70,11 +71,11 @@ using GfwCollision = GfwCollisions::iterator; namespace gfwtrack { DECLARE_SOA_INDEX_COLUMN(GfwCollision, gfwCollision); //! Parent derived collision -DECLARE_SOA_COLUMN(Pt, pt, float); //! pT (GeV/c) -DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity -DECLARE_SOA_COLUMN(Phi, phi, float); //! Azimuthal angle -DECLARE_SOA_COLUMN(HasTOF, hasTOF, bool); //! Track has a TOF match -DECLARE_SOA_COLUMN(TrackSel, trackSel, uint32_t); //! Track and PID bitmask, see TrackSelection +DECLARE_SOA_COLUMN(Pt, pt, float); //! pT (GeV/c) +DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity +DECLARE_SOA_COLUMN(Phi, phi, float); //! Azimuthal angle +DECLARE_SOA_COLUMN(HasTOF, hasTOF, bool); //! Track has a TOF match +DECLARE_SOA_COLUMN(TrackSel, trackSel, uint32_t); //! Track and PID bitmask, see TrackSelection /// Track-quality bits occupy 0-13. TPC PID bits occupy 14-22 and TOF PID bits occupy 23-31. /// Global-track, ITS-cluster, and ITS-chi2 requirements are applied by the producer and are not stored. diff --git a/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx b/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx index 2829e2e2316..327372e0bc3 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwV02.cxx @@ -961,7 +961,7 @@ struct FlowGfwV02 { double v22npt = 0; double ptMeanMid = pidStates.hPtMid[PidCharged]->GetMean(); double ptFractionMid = 0.; - double dnxAB = fGFW->Calculate(corrconfigs.at(0), 0, kTRUE).real(); // V22 weight for AB + double dnxAB = fGFW->Calculate(corrconfigs.at(0), 0, kTRUE).real(); // V22 weight for AB double dnxABAB = fGFW->Calculate(corrconfigs.at(1), 0, kTRUE).real(); // V22 weight for ABAB // Check if the V22 weight is valid @@ -1380,24 +1380,24 @@ struct FlowGfwV02 { { // Fill Nsigma QA. The derived tables store PID bits and have no nSigma columns. if constexpr (requires { track.tpcNSigmaPi(); track.tofNSigmaPi(); }) { - if (!cfgUseItsPID) { - if (ft == kBefore) { - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_pions"), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_kaons"), track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt()); - registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_protons"), track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt()); - } else if (ft == kAfter) { - if (pid == PidPions) { + if (!cfgUseItsPID) { + if (ft == kBefore) { registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_pions"), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt()); - } - if (pid == PidKaons) { registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_kaons"), track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt()); - } - if (pid == PidProtons) { registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_protons"), track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt()); + } else if (ft == kAfter) { + if (pid == PidPions) { + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_pions"), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt()); + } + if (pid == PidKaons) { + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_kaons"), track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt()); + } + if (pid == PidProtons) { + registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_protons"), track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt()); + } } + fillPidDedxQA(track, pid); } - fillPidDedxQA(track, pid); - } } }