diff --git a/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx b/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx index e9aa45fd7c8..36726e80d5d 100644 --- a/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx +++ b/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx @@ -80,73 +80,8 @@ enum AssociatedParticleType { kAssocPartSize }; -/*enum ParticleOfInterest { - Phi = 0, - K0S, - Pion, - ParticleOfInterestSize -};*/ - -/* -#define LIST_OF_PARTICLES_OF_INTEREST \ - X(Phi) \ - X(K0S) \ - X(Pion) \ - //X(PionTPC) \ - //X(PionTPCTOF) - -enum ParticleOfInterest { -#define X(name) name, - LIST_OF_PARTICLES_OF_INTEREST -#undef X - ParticleOfInterestSize -}; - -static constexpr std::array particleOfInterestLabels{ -#define X(name) #name, - LIST_OF_PARTICLES_OF_INTEREST -#undef X -}; - -static constexpr auto particleOfInterestLabels = std::to_array({ -#define X(name) #name, - LIST_OF_PARTICLES_OF_INTEREST -#undef X -}); -*/ - using EffMapPtr = std::variant, std::shared_ptr>; -/*struct BoundEfficiencyMap { - using CoordsTuple = std::tuple; - - const TH3* effMap; - CoordsTuple coords; - - BoundEfficiencyMap(const std::shared_ptr& effMap, float x, float y, float z) : effMap(effMap.get()), coords(x, y, z) {} - BoundEfficiencyMap(const std::shared_ptr& effMap, const CoordsTuple& coords) : effMap(effMap.get()), coords(coords) {} - - float getBinEfficiency() const - { - if (!effMap) { - return 1.0f; - } - - const auto& [x, y, z] = coords; - return effMap->GetBinContent(effMap->FindFixBin(x, y, z)); - } - - float interpolateEfficiency() const - { - if (!effMap) { - return 1.0f; - } - - const auto& [x, y, z] = coords; - return effMap->Interpolate(x, y, z); - } -};*/ - struct BoundEfficiencyMap { using CoordsTuple = std::tuple; @@ -156,12 +91,12 @@ struct BoundEfficiencyMap { BoundEfficiencyMap(const EffMapPtr& effMap, float x, float y, float z) : effMap(effMap), coords(x, y, z) {} BoundEfficiencyMap(const EffMapPtr& effMap, const CoordsTuple& coords) : effMap(effMap), coords(coords) {} - float getBinEfficiency() const + std::pair getBinEfficiencyAndError() const { return std::visit( - [this](auto&& mapPtr) -> float { + [this](auto&& mapPtr) -> std::pair { if (!mapPtr) - return 1.0f; + return {1.0f, 0.0f}; const auto& [x, y, z] = coords; @@ -170,9 +105,11 @@ struct BoundEfficiencyMap { // Compile-time branching: generates the exact correct function call if constexpr (std::is_same_v) { - return mapPtr->GetBinContent(mapPtr->FindFixBin(y, z)); // 2D case only + int bin = mapPtr->FindFixBin(y, z); // 2D case only + return {mapPtr->GetBinContent(bin), mapPtr->GetBinError(bin)}; } else { - return mapPtr->GetBinContent(mapPtr->FindFixBin(x, y, z)); // Full 3D case + int bin = mapPtr->FindFixBin(x, y, z); // Full 3D case + return {mapPtr->GetBinContent(bin), mapPtr->GetBinError(bin)}; } }, effMap); @@ -264,7 +201,8 @@ struct PhiStrangenessCorrelation { struct : ConfigurableGroup { Configurable applyEfficiency{"applyEfficiency", false, "Use efficiency for filling histograms"}; Configurable useEffInterpolation{"useEffInterpolation", false, "If true, interpolates efficiency map, else uses bin center"}; - Configurable applyPhiEfficiency{"applyPhiEfficiency", true, "Apply efficiency for Phi candidates"}; + Configurable applyPhiEfficiency{"applyPhiEfficiency", false, "Apply efficiency for Phi candidates"}; + Configurable propagateEffError{"propagateEffError", false, "Propagate efficiency error"}; } efficiencyConfigs; // Configurable for event mixing @@ -324,9 +262,8 @@ struct PhiStrangenessCorrelation { using BinningTypeVertexCent = ColumnBinningPolicy; BinningTypeVertexCent binningOnVertexAndCent{{axisVertexMixing, axisCentralityMixing}, true}; - static constexpr std::array phiMassRegionLabels{"Signal", "Sideband"}; - // static constexpr std::array particleOfInterestLabels{"Phi", "K0S", "Pion" /*"PionTPC", "PionTPCTOF"*/}; - static constexpr std::array assocParticleLabels{"K0S", "Xi", "Pi"}; + static constexpr std::array PhiMassRegionLabels{"Signal", "Sideband"}; + static constexpr std::array AssocParticleLabels{"K0S", "Xi", "Pi"}; // Light structures to store only the necessary information for the correlation analysis at MCGen level struct MiniParticle { @@ -377,7 +314,7 @@ struct PhiStrangenessCorrelation { histos.add("phiPi/h6PhiPiTPCDataME", "Invariant mass of Phi vs nSigmaTPC Pion in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, massPhiAxis, nSigmaPiAxis}); histos.add("phiPi/h6PhiPiTOFDataME", "Invariant mass of Phi vs nSigmaTOF Pion in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, massPhiAxis, nSigmaPiAxis}); - for (const auto& label : phiMassRegionLabels) { + for (const auto& label : PhiMassRegionLabels) { histos.add(fmt::format("phiK0S/h5PhiK0SData{}", label).c_str(), "Deltay vs deltaphi for Phi and K0Short in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis}); histos.add(fmt::format("phiPi/h5PhiPiData{}", label).c_str(), "Deltay vs deltaphi for Phi and Pion in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis}); @@ -459,7 +396,7 @@ struct PhiStrangenessCorrelation { // Only load the associated maps that are explicitly enabled for (size_t i = 0; i < kAssocPartSize; ++i) { if (doAssocCorrelations[i]) - fetchSingleEfficiencyMapFromCCDB(effMapsAssoc[i], assocParticleLabels[i]); + fetchSingleEfficiencyMapFromCCDB(effMapsAssoc[i], AssocParticleLabels[i]); } } @@ -473,14 +410,39 @@ struct PhiStrangenessCorrelation { // Compute weight based on efficiencies template - float computeWeight(const BoundEffMaps&... boundEffMaps) + std::pair computeWeightAndError(const BoundEffMaps&... boundEffMaps) { if (!efficiencyConfigs.applyEfficiency) - return 1.0f; + return {1.0f, 0.0f}; + + float effTot = 1.0f; + float relErrSqSum = 0.0f; + + auto processMap = [&](const auto& boundMap) { + auto [eff, err] = boundMap.getBinEfficiencyAndError(); // Unpack efficiency and error from the map + + if (efficiencyConfigs.useEffInterpolation) { + eff = boundMap.interpolateEfficiency(); + // For simplicity, we keep the error from the bin content even when interpolating, but this can be refined if needed + } + + effTot *= eff; - float totalEfficiency = ((efficiencyConfigs.useEffInterpolation ? boundEffMaps.interpolateEfficiency() : boundEffMaps.getBinEfficiency()) * ...); + if (eff > 0.0f) { + float mapErr = efficiencyConfigs.propagateEffError ? err : 0.0f; // Optionally propagate error, otherwise treat as zero + relErrSqSum += (mapErr / eff) * (mapErr / eff); + } + }; + + (processMap(boundEffMaps), ...); // Fold expression to process all bound efficiency maps + + if (effTot <= 0.0f) + return {1.0f, 0.0f}; - return totalEfficiency <= 0.0f ? 1.0f : 1.0f / totalEfficiency; + float weight = 1.0f / effTot; + float weightErr = weight * std::sqrt(relErrSqSum); // Propagate relative error to the weight + + return {weight, weightErr}; } float getDeltaPhi(float phiTrigger, float phiAssociated) @@ -497,6 +459,72 @@ struct PhiStrangenessCorrelation { return std::distance(binsMult->begin(), it) - 1; } + template + void customFillHist(auto histId, const std::pair& weightPair, Args... coords) + { + auto hist = histos.get(histId); + if (!hist) { + return; + } + + // Extract the weight and its propagated uncertainty + auto [w, wErr] = weightPair; + + // Find the global bin number for the given physical coordinates + int bin = hist->FindFixBin(coords...); + + // Retrieve the previous bin content and its absolute error + double prevContent = hist->GetBinContent(bin); + double prevErr = hist->GetBinError(bin); + + // Calculate the new content by adding the current weight + double newContent = prevContent + w; + + // Error propagation in quadrature: + // prevErr^2 : previous accumulated variance + // (w * w) : statistical fluctuation added by this specific particle count + // (wErr * wErr) : systematic uncertainty added by the efficiency map inaccuracy + double newErr = std::sqrt(prevErr * prevErr + (w * w) + (wErr * wErr)); + + // Update the histogram bin + hist->SetBinContent(bin, newContent); + hist->SetBinError(bin, newErr); + } + + template + void customFillTHn(auto histId, const std::pair& weightPair, Args... coords) + { + auto hist = histos.get(histId); + if (!hist) { + return; + } + + // Extract the weight and its propagated uncertainty + auto [w, wErr] = weightPair; + + // Create an array of floats for the coordinates to pass to GetBin and SetBinContent + double coordArray[] = {static_cast(coords)...}; + + // Find the bin number. + // The 'true' flag is mandatory: it allocates the bin in memory if it doesn't exist yet. + auto bin = hist->GetBin(coordArray, true); + + // Retrieve the previous content and the squared error (variance) + double prevContent = hist->GetBinContent(bin); + double prevErr2 = hist->GetBinError2(bin); + + // Calculate the new content + double newContent = prevContent + w; + + // Add the new variances to the accumulated variance + // (w * w) is the statistical term, (wErr * wErr) is the map uncertainty term + double newErr2 = prevErr2 + (w * w) + (wErr * wErr); + + // Update the THnSparse bin + hist->SetBinContent(bin, newContent); + hist->SetBinError2(bin, newErr2); // Note: SetBinError2 takes the squared error directly + } + template void processPhiK0SPionSE(TCollision const& collision, TPhiCands const& phiCandidates, TK0SCands const& k0sReduced, TPionCands const& pionTracks) { @@ -527,9 +555,11 @@ struct PhiStrangenessCorrelation { if (efficiencyConfigs.applyEfficiency && efficiencyConfigs.applyPhiEfficiency && phiCand.pt() >= binspTPhi->back()) continue; - float weightPhi = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); + // float weightPhi = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); + auto weightPhi = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); - histos.fill(HIST("phi/h3PhiData"), multiplicity, phiCand.pt(), phiCand.m(), weightPhi); + // histos.fill(HIST("phi/h3PhiData"), multiplicity, phiCand.pt(), phiCand.m(), weightPhi); + customFillHist(HIST("phi/h3PhiData"), weightPhi, multiplicity, phiCand.pt(), phiCand.m()); auto processCorrelations = [&](auto fillK0S, auto fillPion) { if (doAssocCorrelations[kK0S]) { @@ -538,8 +568,10 @@ struct PhiStrangenessCorrelation { if (!isK0sValid(k0s)) continue; - float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); + // float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); + auto weightPhiK0S = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); fillK0S(k0s, weightPhiK0S); } } @@ -550,8 +582,10 @@ struct PhiStrangenessCorrelation { if (!isPionValid(pionTrack)) continue; - float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, pionTrack.pt(), pionTrack.y())); + // float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, pionTrack.pt(), pionTrack.y())); + auto weightPhiPion = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, pionTrack.pt(), pionTrack.y())); fillPion(pionTrack, weightPhiPion); } } @@ -563,41 +597,50 @@ struct PhiStrangenessCorrelation { auto piTOFHistID = HIST("phiPi/h6PhiPiTOFData"); processCorrelations( - [&](const auto& k0s, float w) { - histos.fill(k0sHistID, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m(), w); + //[&](const auto& k0s, float w) { + // histos.fill(k0sHistID, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m(), w); + [&](const auto& k0s, const std::pair& w) { + customFillTHn(k0sHistID, w, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m()); }, - [&](const auto& pion, float w) { - histos.fill(piTPCHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC(), w); - histos.fill(piTOFHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF(), w); + //[&](const auto& pion, float w) { + // histos.fill(piTPCHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC(), w); + // histos.fill(piTOFHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF(), w); + [&](const auto& pion, const std::pair& w) { + customFillTHn(piTPCHistID, w, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC()); + customFillTHn(piTOFHistID, w, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF()); }); } else if (analysisMode == kDeltaYvsDeltaPhi) { auto k0sHistID = std::make_tuple(HIST("phiK0S/h5PhiK0SDataSignal"), HIST("phiK0S/h5PhiK0SDataSideband")); auto piHistID = std::make_tuple(HIST("phiPi/h5PhiPiDataSignal"), HIST("phiPi/h5PhiPiDataSideband")); - static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; + static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[i]; + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) return; - // auto k0sHistID = HIST("phiK0S/h5PhiK0SData") + HIST(phiMassRegionLabels[i]); - // auto piHistID = HIST("phiPi/h5PhiPiData") + HIST(phiMassRegionLabels[i]); + // auto k0sHistID = HIST("phiK0S/h5PhiK0SData") + HIST(PhiMassRegionLabels[Idx]); + // auto piHistID = HIST("phiPi/h5PhiPiData") + HIST(PhiMassRegionLabels[Idx]); processCorrelations( - [&](const auto& k0s, float w) { - histos.fill(std::get(k0sHistID), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), w); + //[&](const auto& k0s, float w) { + // histos.fill(std::get(k0sHistID), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), w); + [&](const auto& k0s, const std::pair& w) { + customFillTHn(std::get(k0sHistID), w, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi())); }, - [&](const auto& pion, float w) { - histos.fill(std::get(piHistID), multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi()), w); + //[&](const auto& pion, float w) { + // histos.fill(std::get(piHistID), multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi()), w); + [&](const auto& pion, const std::pair& w) { + customFillTHn(std::get(piHistID), w, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi())); }); }); } - /*static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; + /*static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[i]; + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) return; @@ -612,7 +655,7 @@ struct PhiStrangenessCorrelation { float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), BoundEfficiencyMap(effMaps[K0S], multiplicity, k0s.pt(), k0s.y())); - histos.fill(HIST("phiK0S/h5PhiK0SData") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), weightPhiK0S); + histos.fill(HIST("phiK0S/h5PhiK0SData") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), weightPhiK0S); } // Loop over all primary pion candidates @@ -620,7 +663,7 @@ struct PhiStrangenessCorrelation { float weightPhiPion = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), BoundEfficiencyMap(effMaps[Pion], multiplicity, pionTrack.pt(), pionTrack.y())); - histos.fill(HIST("phiPi/h5PhiPiData") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), pionTrack.pt(), phiCand.y() - pionTrack.y(), getDeltaPhi(phiCand.phi(), pionTrack.phi()), weightPhiPion); + histos.fill(HIST("phiPi/h5PhiPiData") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), pionTrack.pt(), phiCand.y() - pionTrack.y(), getDeltaPhi(phiCand.phi(), pionTrack.phi()), weightPhiPion); } });*/ } @@ -668,8 +711,10 @@ struct PhiStrangenessCorrelation { continue; auto processCorrelations = [&](auto fillK0S) { - float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); + // float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); + auto weightPhiK0S = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); fillK0S(k0s, weightPhiK0S); }; @@ -677,30 +722,34 @@ struct PhiStrangenessCorrelation { auto k0sHistID = HIST("phiK0S/h6PhiK0SDataME"); processCorrelations( - [&](const auto& k0s, float w) { - histos.fill(k0sHistID, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m(), w); + //[&](const auto& k0s, float w) { + // histos.fill(k0sHistID, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m(), w); + [&](const auto& k0s, const std::pair& w) { + customFillTHn(k0sHistID, w, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m()); }); } else if (analysisMode == kDeltaYvsDeltaPhi) { auto k0sHistID = std::make_tuple(HIST("phiK0S/h5PhiK0SDataMESignal"), HIST("phiK0S/h5PhiK0SDataMESideband")); - static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; + static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[i]; + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) return; processCorrelations( - [&](const auto& k0s, float w) { - histos.fill(std::get(k0sHistID), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), w); + //[&](const auto& k0s, float w) { + // histos.fill(std::get(k0sHistID), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), w); + [&](const auto& k0s, const std::pair& w) { + customFillTHn(std::get(k0sHistID), w, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi())); }); }); } - /*static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; + /*static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[i]; + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) return; @@ -713,7 +762,7 @@ struct PhiStrangenessCorrelation { float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), BoundEfficiencyMap(effMaps[K0S], multiplicity, k0s.pt(), k0s.y())); - histos.fill(HIST("phiK0S/h5PhiK0SDataME") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), weightPhiK0S); + histos.fill(HIST("phiK0S/h5PhiK0SDataME") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), weightPhiK0S); });*/ } } @@ -763,8 +812,10 @@ struct PhiStrangenessCorrelation { continue; auto processCorrelations = [&](auto fillPion) { - float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, piTrack.pt(), piTrack.y())); + // float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, piTrack.pt(), piTrack.y())); + auto weightPhiPion = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, piTrack.pt(), piTrack.y())); fillPion(piTrack, weightPhiPion); }; @@ -773,38 +824,43 @@ struct PhiStrangenessCorrelation { auto piTOFHistID = HIST("phiPi/h6PhiPiTOFDataME"); processCorrelations( - [&](const auto& pion, float w) { - histos.fill(piTPCHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC(), w); - histos.fill(piTOFHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF(), w); + //[&](const auto& pion, float w) { + // histos.fill(piTPCHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC(), w); + // histos.fill(piTOFHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF(), w); + [&](const auto& pion, const std::pair& w) { + customFillTHn(piTPCHistID, w, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC()); + customFillTHn(piTOFHistID, w, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF()); }); } else if (analysisMode == kDeltaYvsDeltaPhi) { auto piHistID = std::make_tuple(HIST("phiPi/h5PhiPiDataMESignal"), HIST("phiPi/h5PhiPiDataMESideband")); - static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; + static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[i]; + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) return; processCorrelations( - [&](const auto& pion, float w) { - histos.fill(std::get(piHistID), multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi()), w); + //[&](const auto& pion, float w) { + // histos.fill(std::get(piHistID), multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi()), w); + [&](const auto& pion, const std::pair& w) { + customFillTHn(std::get(piHistID), w, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi())); }); }); } - /*static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; + /*static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[i]; + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) return; float weightPhiPion = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), BoundEfficiencyMap(effMaps[Pion], multiplicity, piTrack.pt(), piTrack.y())); - histos.fill(HIST("phiPi/h5PhiPiDataME") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), piTrack.pt(), phiCand.y() - piTrack.y(), getDeltaPhi(phiCand.phi(), piTrack.phi()), weightPhiPion); + histos.fill(HIST("phiPi/h5PhiPiDataME") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), piTrack.pt(), phiCand.y() - piTrack.y(), getDeltaPhi(phiCand.phi(), piTrack.phi()), weightPhiPion); });*/ } } @@ -999,12 +1055,12 @@ struct PhiStrangenessCorrelation { auto& phiParticle = mcParticles.rawIteratorAt(phiIndices[iTrigg]); static_for<0, assocIndices.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; + constexpr unsigned int Idx = i_idx.value; - for (std::size_t iAssoc{0}; iAssoc < assocIndices[i].size(); ++iAssoc) { - auto& assocParticle = mcParticles.rawIteratorAt(assocIndices[i][iAssoc]); + for (std::size_t iAssoc{0}; iAssoc < assocIndices[Idx].size(); ++iAssoc) { + auto& assocParticle = mcParticles.rawIteratorAt(assocIndices[Idx][iAssoc]); - histos.fill(HIST("mcGenClosure/h5Phi") + HIST(assocParticleLabels[i]) + HIST("ClosureGenSE"), multiplicity, phiParticle.pt(), assocParticle.pt(), phiParticle.y() - assocParticle.y(), getDeltaPhi(phiParticle.phi(), assocParticle.phi())); + histos.fill(HIST("mcGenClosure/h5Phi") + HIST(AssocParticleLabels[Idx]) + HIST("ClosureGenSE"), multiplicity, phiParticle.pt(), assocParticle.pt(), phiParticle.y() - assocParticle.y(), getDeltaPhi(phiParticle.phi(), assocParticle.phi())); } }); } @@ -1130,13 +1186,13 @@ struct PhiStrangenessCorrelation { for (const auto& phiParticle : phiParticles) { histos.fill(HIST("phi/h3PhiMCClosureGen"), multiplicity, phiParticle.pt, phiParticle.y); - static_for<0, assocParticleLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; - if (!doAssocCorrelations[i]) + static_for<0, AssocParticleLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; + if (!doAssocCorrelations[Idx]) return; - for (const auto& assocParticle : *(currentAssocParticles[i])) { - histos.fill(HIST("phi") + HIST(assocParticleLabels[i]) + HIST("/h5Phi") + HIST(assocParticleLabels[i]) + HIST("ClosureMCGen"), + for (const auto& assocParticle : *(currentAssocParticles[Idx])) { + histos.fill(HIST("phi") + HIST(AssocParticleLabels[Idx]) + HIST("/h5Phi") + HIST(AssocParticleLabels[Idx]) + HIST("ClosureMCGen"), multiplicity, phiParticle.pt, assocParticle.pt, phiParticle.y - assocParticle.y, getDeltaPhi(phiParticle.phi, assocParticle.phi)); @@ -1150,13 +1206,13 @@ struct PhiStrangenessCorrelation { // Loop over past events in the same multiplicity bin and fill histograms with all combinations of current phi particles and past associated particles for (const auto& phiParticle : phiParticles) { - static_for<0, assocParticleLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int i = i_idx.value; - if (!doAssocCorrelations[i]) + static_for<0, AssocParticleLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; + if (!doAssocCorrelations[Idx]) return; - for (const auto& assocParticle : *(pastAssocParticles[i])) { - histos.fill(HIST("phi") + HIST(assocParticleLabels[i]) + HIST("/h5Phi") + HIST(assocParticleLabels[i]) + HIST("ClosureMCGenME"), + for (const auto& assocParticle : *(pastAssocParticles[Idx])) { + histos.fill(HIST("phi") + HIST(AssocParticleLabels[Idx]) + HIST("/h5Phi") + HIST(AssocParticleLabels[Idx]) + HIST("ClosureMCGenME"), multiplicity, phiParticle.pt, assocParticle.pt, phiParticle.y - assocParticle.y, getDeltaPhi(phiParticle.phi, assocParticle.phi));