From e677746964c245df43ff49a8f26ed71885a380cb Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Fri, 20 Mar 2026 15:45:06 +0100 Subject: [PATCH] add resonances --- .../Tasks/flowGenericFramework.cxx | 712 ++++++++++++++++-- 1 file changed, 649 insertions(+), 63 deletions(-) diff --git a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx index de8731f61ca..e848d353761 100644 --- a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx @@ -22,19 +22,30 @@ #include "GFWWeights.h" #include "GFWWeightsList.h" +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" +#include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseITS.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "CommonConstants/PhysicsConstants.h" #include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/RunningWorkflowInfo.h" +#include "Framework/StepTHn.h" #include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/Track.h" #include #include #include @@ -46,6 +57,7 @@ #include #include +#include #include #include #include @@ -85,6 +97,31 @@ std::vector multGlobalV0ACutPars; std::vector multGlobalT0ACutPars; } // namespace o2::analysis::gfw +template +auto projectMatrix(Array2D const& mat, std::array& array1, std::array& array2, std::array& array3) +{ + for (auto j = 0; j < static_cast(mat.cols); ++j) { + array1[j] = mat(0, j); + array2[j] = mat(1, j); + array3[j] = mat(2, j); + } + return; +} +template +auto readMatrix(Array2D const& mat, P& array) +{ + for (auto i = 0; i < static_cast(mat.rows); ++i) { + for (auto j = 0; j < static_cast(mat.cols); ++j) { + array[i][j] = mat(i, j); + } + } + + return; +} + +static constexpr float LongArrayFloat[3][20] = {{1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {2.1, 2.2, 2.3, -2.1, -2.2, -2.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {3.1, 3.2, 3.3, -3.1, -3.2, -3.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}}; +static constexpr int LongArrayInt[3][20] = {{1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1}, {2, 2, 2, -2, -2, -2, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1}, {3, 3, 3, -3, -3, -3, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1}}; + struct FlowGenericFramework { O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples") @@ -105,6 +142,7 @@ struct FlowGenericFramework { O2_DEFINE_CONFIGURABLE(cfgPtmax, float, 10, "maximum pt (GeV/c)"); O2_DEFINE_CONFIGURABLE(cfgEta, float, 0.8, "eta cut"); O2_DEFINE_CONFIGURABLE(cfgEtaPtPt, float, 0.4, "eta cut for pt-pt correlations"); + O2_DEFINE_CONFIGURABLE(cfgUsePIDTotal, bool, false, "use fraction of PID total"); O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10, "vertex cut (cm)"); struct : ConfigurableGroup { O2_DEFINE_CONFIGURABLE(cfgDCAxyNSigma, float, 7, "Cut on number of sigma deviations from expected DCA in the transverse direction"); @@ -151,9 +189,22 @@ struct FlowGenericFramework { O2_DEFINE_CONFIGURABLE(cfgGlobalT0ALowSigma, float, -3., "Number of sigma deviations below expected value in global vs T0A correlation"); O2_DEFINE_CONFIGURABLE(cfgGlobalT0AHighSigma, float, 4, "Number of sigma deviations above expected value in global vs T0A correlation"); } cfgGlobalAsideCorrCuts; - struct ConfigurableGroup { + struct : ConfigurableGroup { + Configurable> nSigmas{"nSigmas", {LongArrayFloat[0], 3, 6, {"TPC", "TOF", "ITS"}, {"pos_pi", "pos_ka", "pos_pr", "neg_pi", "neg_ka", "neg_pr"}}, "Labeled array for n-sigma values for TPC, TOF, ITS for pions, kaons, protons (positive and negative)"}; + Configurable> resonanceCuts{"resonanceCuts", {LongArrayFloat[0], 3, 11, {"K0", "Lambda", "Phi"}, {"cos_PAs", "massMin", "massMax", "PosTrackPt", "NegTrackPt", "DCAPosToPVMin", "DCANegToPVMin", "Lifetime", "RadiusMin", "RadiusMax", "Rapidity"}}, "Labeled array (float) for various cuts on resonances"}; + Configurable> resonanceSwitches{"resonanceSwitches", {LongArrayInt[0], 3, 6, {"K0", "Lambda", "Phi"}, {"UseParticle", "UseCosPA", "NMassBins", "DCABetDaug", "UseProperLifetime", "UseV0Radius"}}, "Labeled array (int) for various cuts on resonances"}; + O2_DEFINE_CONFIGURABLE(cfgUseLsPhi, bool, true, "Use LikeSign for Phi v2") + O2_DEFINE_CONFIGURABLE(cfgUseOnlyTPC, bool, true, "Use only TPC PID for daughter selection") + O2_DEFINE_CONFIGURABLE(cfgFakeKaonCut, float, 0.1f, "Maximum difference in measured momentum and TPC inner ring momentum of particle") + O2_DEFINE_CONFIGURABLE(cfgUseAsymmetricPID, bool, false, "Use asymmetric PID cuts") + O2_DEFINE_CONFIGURABLE(cfgTPCNsigmaCut, float, 3.0f, "TPC N-sigma cut for pions, kaons, protons") O2_DEFINE_CONFIGURABLE(cfgUseStrictPID, bool, true, "Use strict PID cuts for TPC") O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 0.5, "pt cut on TOF for PID"); + O2_DEFINE_CONFIGURABLE(cfgUseItsPID, bool, true, "Use ITS PID for particle identification") + O2_DEFINE_CONFIGURABLE(cfgK0SignalMin, float, 0.48, "Minimum cut on K0 mT signal (upper limit of left sideband)"); + O2_DEFINE_CONFIGURABLE(cfgK0SignalMax, float, 0.51, "Minimum cut on K0 mT signal (lower limit of right sideband)"); + O2_DEFINE_CONFIGURABLE(cfgLambdaSignalMin, float, 1.1, "Minimum cut on Lambda mT signal (upper limit of left sideband)"); + O2_DEFINE_CONFIGURABLE(cfgLambdaSignalMax, float, 1.3, "Minimum cut on Lambda mT signal (lower limit of right sideband)"); } cfgPIDCuts; Configurable cfgGFWBinning{"cfgGFWBinning", {40, 16, 72, 300, 0, 3000, 0.2, 10.0, 0.2, 3.0, {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.5, 5, 5.5, 6, 7, 8, 9, 10}, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90}}, "Configuration for binning"}; @@ -162,10 +213,14 @@ struct FlowGenericFramework { Configurable cfgCorrConfig{"cfgCorrConfig", {{"refP {2} refN {-2}", "refP {3} refN {-3}", "refP {4} refN {-4}", "refFull {2 -2}", "refFull {2 2 -2 -2}"}, {"ChGap22", "ChGap32", "ChGap42", "ChFull22", "ChFull24"}, {0, 0, 0, 0, 0}, {15, 1, 1, 0, 0}}, "Configurations for each correlation to calculate"}; Configurable cfgCorrConfigRadial{"cfgCorrConfigRadial", {{"refP {2} refN {-2}", "refP {3} refN {-3}", "refP {4} refN {-4}"}, {"ChGap22", "ChGap32", "ChGap42"}, {1, 1, 1}, {0, 0, 0}}, "Configurations for each radial flow correlation to calculate"}; - // #include "PWGCF/TwoParticleCorrelations/TableProducer/Productions/skimmingconf_20221115.cxx" // NOLINT + ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {80, -5, 5}, "nsigmaTPC axis"}; + ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"}; + // Connect to ccdb Service ccdb; + o2::aod::ITSResponse itsResponse; + struct Config { TH1D* mEfficiency = nullptr; std::vector mAcceptance; @@ -181,6 +236,12 @@ struct FlowGenericFramework { OutputObj fFCgen{FlowContainer("FlowContainer_gen")}; HistogramRegistry registry{"registry"}; + std::array, 3> resoCutVals; + std::array, 3> resoSwitchVals; + std::array tofNsigmaCut; + std::array itsNsigmaCut; + std::array tpcNsigmaCut; + std::vector fFCpts = {&(*fFCpt_ch), &(*fFCpt_pi), &(*fFCpt_ka), &(*fFCpt_pr)}; // QA outputs @@ -249,6 +310,52 @@ struct FlowGenericFramework { {cfgEventCutFlags.cfgIsVertexITSTPC, kIsVertexITSTPC, o2::aod::evsel::kNoTimeFrameBorder}, {cfgEventCutFlags.cfgIsGoodITSLayersAll, kIsGoodITSLayersAll, o2::aod::evsel::kNoITSROFrameBorder}, }; + enum Particles { + PIONS, + KAONS, + PROTONS + }; + enum OutputSpecies { + K0 = 0, + LAMBDA = 1, + PHI = 2, + ANLAMBDA = 3, + REF = 4, + kCount_OutputSpecies + }; + enum ParticleCuts { + kCosPA = 0, + kMassMin, + kMassMax, + kPosTrackPt, + kNegTrackPt, + kDCAPosToPVMin, + kDCANegToPVMin, + kLifeTime, + kRadiusMin, + kRadiusMax, + kRapidity + }; + enum ParticleSwitches { + kUseParticle = 0, + kUseCosPA, + kMassBins, + kDCABetDaug, + kUseProperLifetime, + kUseV0Radius + }; + enum V0Selection { + kFillCandidate = 0, + kFillDaughterPt, + kFillMassCut, + kFillRapidityCut, + kFillDCAtoPV, + kFillDCAxDaughters, + kFillV0Radius, + kFillCosPA, + kFillProperLifetime, + kFillDaughterTrackSelection + }; // Define global variables // Generic Framework @@ -336,6 +443,11 @@ struct FlowGenericFramework { o2::analysis::gfw::multGlobalV0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalV0ACutPars; o2::analysis::gfw::multGlobalT0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalT0ACutPars; + projectMatrix(cfgPIDCuts.nSigmas->getData(), tpcNsigmaCut, tofNsigmaCut, itsNsigmaCut); + readMatrix(cfgPIDCuts.resonanceCuts->getData(), resoCutVals); + readMatrix(cfgPIDCuts.resonanceSwitches->getData(), resoSwitchVals); + PrintResoCuts(); + AxisSpec phiAxis = {o2::analysis::gfw::phibins, o2::analysis::gfw::philow, o2::analysis::gfw::phiup, "#phi"}; AxisSpec phiModAxis = {100, 0, constants::math::PI / 9, "fmod(#varphi,#pi/9)"}; AxisSpec etaAxis = {o2::analysis::gfw::etabins, -cfgEta, cfgEta, "#eta"}; @@ -359,6 +471,7 @@ struct FlowGenericFramework { : centAxis; AxisSpec dcaZAXis = {200, -2, 2, "DCA_{z} (cm)"}; AxisSpec dcaXYAXis = {200, -1, 1, "DCA_{xy} (cm)"}; + AxisSpec singleCount = {1, 0, 1}; ccdb->setURL("http://alice-ccdb.cern.ch"); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -389,6 +502,8 @@ struct FlowGenericFramework { registry.addClone("trackQA/before/", "trackQA/after/"); registry.add("trackQA/after/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfw::ptreflow, o2::analysis::gfw::ptrefup}}}); registry.add("trackQA/after/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfw::ptpoilow, o2::analysis::gfw::ptpoiup}}}); + registry.add("trackQA/after/Nch_corrected", "", {HistType::kTH1D, {nchAxis}}); + registry.add("trackQA/after/Nch_uncorrected", "", {HistType::kTH1D, {nchAxis}}); registry.add("eventQA/before/centrality", "; centrality (%); Counts", {HistType::kTH1D, {centAxis}}); registry.add("eventQA/before/multiplicity", "; N_{ch}; Counts", {HistType::kTH1D, {nchAxis}}); @@ -401,23 +516,33 @@ struct FlowGenericFramework { registry.add("eventQA/before/multT0C_centT0C", "; multT0C; FT0C centrality (%)", {HistType::kTH2D, {centAxis, t0cAxis}}); registry.add("eventQA/before/occ_mult_cent", "; occupancy; N_{ch}; centrality (%)", {HistType::kTH3D, {occAxis, nchAxis, centAxis}}); registry.addClone("eventQA/before/", "eventQA/after/"); - registry.add("eventQA/eventSel", "Number of Events;; Counts", {HistType::kTH1D, {{11, 0, 11}}}); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(1, "Filtered event"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(2, "sel8"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(3, "occupancy"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(4, "kTVXinTRD"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(5, "kNoSameBunchPileup"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(6, "kIsGoodZvtxFT0vsPV"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(7, "kNoCollInTimeRangeStandard"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(8, "kIsVertexITSTPC"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(9, "kIsGoodITSLayersAll"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(10, "after Mult cuts"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(11, "has track + within cent"); - - registry.add("npt_ch", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTH2D, {ptAxis, centAxis}}); - registry.add("npt_pi", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTH2D, {ptAxis, centAxis}}); - registry.add("npt_ka", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTH2D, {ptAxis, centAxis}}); - registry.add("npt_pr", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTH2D, {ptAxis, centAxis}}); + registry.add("eventQA/eventSel", "Number of Events;; Counts", {HistType::kTH1D, {{15, 0.5, 15.5}}}); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kFilteredEvent, "Filtered event"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kSel8, "sel8"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kOccupancy, "occupancy"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kTVXinTRD, "kTVXinTRD"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoSameBunchPileup, "kNoSameBunchPileup"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kIsGoodZvtxFT0vsPV, "kIsGoodZvtxFT0vsPV"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoCollInTimeRangeStandard, "kNoCollInTimeRangeStandard"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoCollInRofStandard, "kNoCollInRofStandard"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoHighMultCollInPrevRof, "kNoHighMultCollInPrevRof"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoTimeFrameBorder, "kNoTimeFrameBorder"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoITSROFrameBorder, "kNoITSROFrameBorder"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kIsVertexITSTPC, "kIsVertexITSTPC"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kIsGoodITSLayersAll, "kIsGoodITSLayersAll"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kMultCuts, "after Mult cuts"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kTrackCent, "has track + within cent"); + + registry.add("npt_ch", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_pi", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_ka", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_pr", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_K0_sig", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_K0_sb1", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_K0_sb2", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_Lambda_sig", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_Lambda_sb1", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); + registry.add("npt_Lambda_sb2", "; #it{p}_{T} (GeV/#it{c}; ; centrality (%); fraction)", {HistType::kTProfile2D, {ptAxis, centAxis}}); if (!cfgRunByRun) { if (cfgUsePID) { @@ -431,8 +556,62 @@ struct FlowGenericFramework { } } - registry.add("trackQA/after/Nch_corrected", "", {HistType::kTH1D, {nchAxis}}); - registry.add("trackQA/after/Nch_uncorrected", "", {HistType::kTH1D, {nchAxis}}); + AxisSpec axisK0Mass = {resoSwitchVals[K0][kMassBins], resoCutVals[K0][kMassMin], resoCutVals[K0][kMassMax]}; + AxisSpec axisLambdaMass = {resoSwitchVals[LAMBDA][kMassBins], resoCutVals[LAMBDA][kMassMin], resoCutVals[LAMBDA][kMassMax]}; + + // QA histograms for V0s + if (resoSwitchVals[K0][kUseParticle]) { + registry.add("K0/PiPlusTPC_K0", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}}); + registry.add("K0/PiMinusTPC_K0", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}}); + registry.add("K0/PiPlusTOF_K0", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTOF}}}); + registry.add("K0/PiMinusTOF_K0", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTOF}}}); + registry.add("K0/hK0Phi", "", {HistType::kTH1D, {phiAxis}}); + registry.add("K0/hK0Eta", "", {HistType::kTH1D, {etaAxis}}); + registry.add("K0/hK0Mass_sparse", "", {HistType::kTHnSparseF, {{axisK0Mass, ptAxis, nchAxis}}}); + registry.add("K0/hK0s", "", {HistType::kTH1D, {singleCount}}); + + registry.add("K0/hK0Count", "Number of K0;; Count", {HistType::kTH1D, {{10, 0.5, 10.5}}}); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillCandidate, "K0 candidates"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillDaughterPt, "Daughter pt"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillMassCut, "Mass cut"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillRapidityCut, "Rapidity cut"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillDCAtoPV, "DCA to PV"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillDCAxDaughters, "DCA between daughters"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillV0Radius, "V0radius"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillCosPA, "CosPA"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillProperLifetime, "Proper lifetime"); + registry.get(HIST("K0/hK0Count"))->GetXaxis()->SetBinLabel(kFillDaughterTrackSelection, "Daughter track selection"); + } + + if (resoSwitchVals[LAMBDA][kUseParticle]) { + registry.add("Lambda/PrPlusTPC_L", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}}); + registry.add("Lambda/PiMinusTPC_L", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}}); + registry.add("Lambda/PrPlusTOF_L", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTOF}}}); + registry.add("Lambda/PiMinusTOF_L", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTOF}}}); + registry.add("Lambda/hLambdaPhi", "", {HistType::kTH1D, {phiAxis}}); + registry.add("Lambda/hLambdaEta", "", {HistType::kTH1D, {etaAxis}}); + registry.add("Lambda/hLambdaMass_sparse", "", {HistType::kTHnSparseF, {{axisLambdaMass, ptAxis, nchAxis}}}); + registry.add("Lambda/PiPlusTPC_AL", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}}); + registry.add("Lambda/PrMinusTPC_AL", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}}); + registry.add("Lambda/PiPlusTOF_AL", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTOF}}}); + registry.add("Lambda/PrMinusTOF_AL", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTOF}}}); + registry.add("Lambda/hAntiLambdaPhi", "", {HistType::kTH1D, {phiAxis}}); + registry.add("Lambda/hAntiLambdaEta", "", {HistType::kTH1D, {etaAxis}}); + registry.add("Lambda/hAntiLambdaMass_sparse", "", {HistType::kTHnSparseF, {{axisLambdaMass, ptAxis, nchAxis}}}); + registry.add("Lambda/hLambdas", "", {HistType::kTH1D, {singleCount}}); + + registry.add("Lambda/hLambdaCount", "Number of Lambda;; Count", {HistType::kTH1D, {{10, 0.5, 10.5}}}); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillCandidate, "Lambda candidates"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillDaughterPt, "Daughter pt"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillMassCut, "Mass cut"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillRapidityCut, "Rapidity cut"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillDCAtoPV, "DCA to PV"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillDCAxDaughters, "DCA between daughters"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillV0Radius, "V0radius"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillCosPA, "CosPA"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillProperLifetime, "Proper lifetime"); + registry.get(HIST("Lambda/hLambdaCount"))->GetXaxis()->SetBinLabel(kFillDaughterTrackSelection, "Daughter track selection"); + } } if (o2::analysis::gfw::regions.GetSize() < 0) @@ -546,6 +725,70 @@ struct FlowGenericFramework { static constexpr std::string_view FillTimeName[] = {"before/", "after/"}; + void PrintResoCuts() + { + auto printTable = [](const auto& lbl, const auto& valuesMatrix, const std::string& title) { + LOGF(info, "===== %s =====", title.c_str()); + + size_t nRows = lbl.labels_rows.size(); + size_t nCols = lbl.labels_cols.size(); + + // Compute column widths: max of header or value string length + std::vector colWidths(nCols); + for (size_t c = 0; c < nCols; ++c) { + size_t maxValLen = 0; + for (size_t r = 0; r < nRows; ++r) { + std::ostringstream oss; + oss << std::fixed << std::setprecision(2) << valuesMatrix[r][c]; + maxValLen = std::max(maxValLen, oss.str().length()); + } + colWidths[c] = std::max(lbl.labels_cols[c].size(), maxValLen); + } + + // Determine row label width + size_t rowLabelWidth = 0; + for (const auto& rowName : lbl.labels_rows) + rowLabelWidth = std::max(rowLabelWidth, rowName.size()); + rowLabelWidth += 2; // for ": " + + // Print header + std::ostringstream header; + header << std::setw(rowLabelWidth) << " "; + for (size_t c = 0; c < nCols; ++c) + header << std::setw(colWidths[c]) << lbl.labels_cols[c] << " "; + LOGF(info, "%s", header.str().c_str()); + + // Print rows + for (size_t r = 0; r < nRows; ++r) { + std::ostringstream line; + line << std::setw(rowLabelWidth) << (lbl.labels_rows[r] + ":"); + for (size_t c = 0; c < nCols; ++c) { + line << std::setw(colWidths[c]) << std::fixed << std::setprecision(2) << valuesMatrix[r][c] << " "; + } + LOGF(info, "%s", line.str().c_str()); + } + }; + + // ----- nSigma PID ----- + // Map arrays into a 2D vector + std::vector> nSigmaVals = { + std::vector(tpcNsigmaCut.begin(), tpcNsigmaCut.end()), + std::vector(tofNsigmaCut.begin(), tofNsigmaCut.end()), + std::vector(itsNsigmaCut.begin(), itsNsigmaCut.end())}; + printTable(cfgPIDCuts.nSigmas.value, nSigmaVals, "nSigma PID Cuts"); + + // ----- Resonance Cuts ----- + std::vector> resoCutsVals(resoCutVals.size()); + for (size_t r = 0; r < resoCutVals.size(); ++r) + resoCutsVals[r] = std::vector(resoCutVals[r].begin(), resoCutVals[r].end()); + printTable(cfgPIDCuts.resonanceCuts.value, resoCutsVals, "Resonance Cuts"); + + // ----- Resonance Switches ----- + std::vector> resoSwitchValsF(resoSwitchVals.size()); + for (size_t r = 0; r < resoSwitchVals.size(); ++r) + resoSwitchValsF[r] = std::vector(resoSwitchVals[r].begin(), resoSwitchVals[r].end()); + printTable(cfgPIDCuts.resonanceSwitches.value, resoSwitchValsF, "Resonance Switches"); + } enum QAFillTime { kBefore, kAfter @@ -664,6 +907,56 @@ struct FlowGenericFramework { return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton } + template + int getNsigmaPIDAssymmetric(TTrack 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 = {itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track)}; + int pid = -1; + + std::array nSigmaToUse = cfgPIDCuts.cfgUseItsPID ? nSigmaITS : nSigmaTPC; // Choose which nSigma to use: TPC or ITS + std::array detectorNsigmaCut = cfgPIDCuts.cfgUseItsPID ? itsNsigmaCut : tpcNsigmaCut; // Choose which nSigma to use: TPC or ITS + + bool isPion, isKaon, isProton; + bool isDetectedPion = nSigmaToUse[0] < detectorNsigmaCut[0] && nSigmaToUse[0] > detectorNsigmaCut[0 + 3]; + bool isDetectedKaon = nSigmaToUse[1] < detectorNsigmaCut[1] && nSigmaToUse[1] > detectorNsigmaCut[1 + 3]; + bool isDetectedProton = nSigmaToUse[2] < detectorNsigmaCut[2] && nSigmaToUse[2] > detectorNsigmaCut[2 + 3]; + + bool isTofPion = nSigmaTOF[0] < tofNsigmaCut[0] && nSigmaTOF[0] > tofNsigmaCut[0 + 3]; + bool isTofKaon = nSigmaTOF[1] < tofNsigmaCut[1] && nSigmaTOF[1] > tofNsigmaCut[1 + 3]; + bool isTofProton = nSigmaTOF[2] < tofNsigmaCut[2] && nSigmaTOF[2] > tofNsigmaCut[2 + 3]; + + if (track.pt() > cfgPIDCuts.cfgTofPtCut && !track.hasTOF()) { + return 0; + } else if (track.pt() > cfgPIDCuts.cfgTofPtCut && track.hasTOF()) { + isPion = isTofPion && isDetectedPion; + isKaon = isTofKaon && isDetectedKaon; + isProton = isTofProton && isDetectedProton; + } else { + isPion = isDetectedPion; + isKaon = isDetectedKaon; + isProton = isDetectedProton; + } + + if ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton)) { + return 0; // more than one particle satisfy the criteria + } + + if (isPion) { + pid = PIONS; + } else if (isKaon) { + pid = KAONS; + } else if (isProton) { + pid = PROTONS; + } else { + return 0; // no particle satisfies the criteria + } + + return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton + } + template bool eventSelected(TCollision collision, const int multTrk, const float& centrality, const int run) { @@ -872,7 +1165,8 @@ struct FlowGenericFramework { struct AcceptedTracks { explicit AcceptedTracks(std::size_t nptbins) - : nch(nptbins, 0.f), + : pidtotal{0, 0, 0}, + nch(nptbins, 0.f), npi(nptbins, 0.f), nka(nptbins, 0.f), npr(nptbins, 0.f) @@ -882,6 +1176,7 @@ struct FlowGenericFramework { float total = 0; unsigned int total_uncorr = 0; + std::vector pidtotal; std::vector nch; std::vector npi; std::vector nka; @@ -898,7 +1193,6 @@ struct FlowGenericFramework { if (!cfgUseGapMethod) container->fillVnPtStdProfiles(centmult, rndm); } - for (uint l_ind = 0; l_ind < corrconfigs.size(); ++l_ind) { if (!corrconfigs.at(l_ind).pTDif) { auto dnx = fGFW->Calculate(corrconfigs.at(l_ind), 0, kTRUE).real(); @@ -924,35 +1218,39 @@ struct FlowGenericFramework { (dt == kGen) ? fFCgen->FillProfile(Form("%s_pt_%i", corrconfigs.at(l_ind).Head.c_str(), i), centmult, val, dnx, rndm) : fFC->FillProfile(Form("%s_pt_%i", corrconfigs.at(l_ind).Head.c_str(), i), centmult, val, dnx, rndm); } } + float chtotal = (cfgUseNchCorrection) ? acceptedtracks.total : acceptedtracks.total_uncorr; - int total = (cfgUseNchCorrection) ? acceptedtracks.total : acceptedtracks.total_uncorr; - - if (total == 0) - return; - - for (std::size_t i = 0; i < acceptedtracks.nch.size(); ++i) - registry.fill(HIST("npt_ch"), fPtAxis->GetBinCenter(i + 1), centmult, acceptedtracks.nch[i] / total); - for (std::size_t i = 0; i < acceptedtracks.npi.size(); ++i) - registry.fill(HIST("npt_pi"), fPtAxis->GetBinCenter(i + 1), centmult, acceptedtracks.npi[i] / total); - for (std::size_t i = 0; i < acceptedtracks.nka.size(); ++i) - registry.fill(HIST("npt_ka"), fPtAxis->GetBinCenter(i + 1), centmult, acceptedtracks.nka[i] / total); - for (std::size_t i = 0; i < acceptedtracks.npr.size(); ++i) - registry.fill(HIST("npt_pr"), fPtAxis->GetBinCenter(i + 1), centmult, acceptedtracks.npr[i] / total); - + // calculate fractions std::vector> inputs = {acceptedtracks.nch, acceptedtracks.npi, acceptedtracks.nka, acceptedtracks.npr}; std::vector> fractions; fractions.reserve(inputs.size()); - + int pidcounter = 0; for (const auto& vec : inputs) { fractions.emplace_back(); fractions.back().reserve(vec.size()); + float total = chtotal; + if (cfgUsePIDTotal) + (pidcounter) ? acceptedtracks.pidtotal[pidcounter] : chtotal; + + if (total == 0.) { + ++pidcounter; + continue; + } std::transform(vec.begin(), vec.end(), std::back_inserter(fractions.back()), [&](double x) { return x / total; }); + ++pidcounter; } - - for (uint l_ind = 0; l_ind < corrconfigsradial.size(); ++l_ind) { + for (std::size_t i = 0; i < fractions[0].size(); ++i) + registry.fill(HIST("npt_ch"), fPtAxis->GetBinCenter(i + 1), centmult, fractions[0][i]); + for (std::size_t i = 0; i < fractions[1].size(); ++i) + registry.fill(HIST("npt_pi"), fPtAxis->GetBinCenter(i + 1), centmult, fractions[1][i]); + for (std::size_t i = 0; i < fractions[2].size(); ++i) + registry.fill(HIST("npt_ka"), fPtAxis->GetBinCenter(i + 1), centmult, fractions[2][i]); + for (std::size_t i = 0; i < fractions[3].size(); ++i) + registry.fill(HIST("npt_pr"), fPtAxis->GetBinCenter(i + 1), centmult, fractions[3][i]); + for (uint l_ind = 0; l_ind < 4; ++l_ind) { for (int i = 1; i <= fPtAxis->GetNbins(); i++) { auto dnx = fGFW->Calculate(corrconfigsradial.at(l_ind), i - 1, kTRUE).real(); if (dnx == 0) @@ -962,21 +1260,20 @@ struct FlowGenericFramework { (dt == kGen) ? fFCgen->FillProfile(Form("%s_pt_%i", corrconfigsradial.at(l_ind).Head.c_str(), i), centmult, val * fractions[l_ind][i - 1], dnx, rndm) : fFC->FillProfile(Form("%s_pt_%i", corrconfigsradial.at(l_ind).Head.c_str(), i), centmult, val * fractions[l_ind][i - 1], dnx, rndm); } } - return; } - template - void processCollision(TCollision collision, TTracks tracks, const float& centrality, const int field, const int run) + template + void processCollision(TCollision collision, TTracks tracks, TV0s v0s, const float& centrality, const int field, const int run) { if (tracks.size() < 1) return; if (dt != kGen && (centrality < o2::analysis::gfw::centbinning.front() || centrality > o2::analysis::gfw::centbinning.back())) return; if (dt != kGen) { - registry.fill(HIST("eventQA/eventSel"), 10.5); + registry.fill(HIST("eventQA/eventSel"), kTrackCent); if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(10.5); + th1sList[run][hEventSel]->Fill(kTrackCent); } float vtxz = collision.posZ(); if (dt != kGen && cfgRunByRun) { @@ -989,7 +1286,6 @@ struct FlowGenericFramework { container->clearVector(); float lRandom = fRndm->Rndm(); - // be cautious, this only works for Pb-Pb // esimate the Event plane and vn for this event DensityCorr densitycorrections; @@ -1024,15 +1320,87 @@ struct FlowGenericFramework { densitycorrections.v4 = v4; densitycorrections.density = tracks.size(); } - + // process tracks AcceptedTracks acceptedTracks(o2::analysis::gfw::ptbinning.size() - 1); for (const auto& track : tracks) { processTrack(track, vtxz, field, run, densitycorrections, acceptedTracks); } - registry.fill(HIST("trackQA/after/Nch_corrected"), acceptedTracks.total); registry.fill(HIST("trackQA/after/Nch_uncorrected"), acceptedTracks.total_uncorr); + std::vector> npt_resonances(6, std::vector(o2::analysis::gfw::ptbinning.size())); + // Process V0s + for (const auto& v0 : v0s) { + if (resoSwitchVals[K0][kUseParticle]) { + if (selectK0(collision, v0, centrality)) { + int ptBinIndex = fPtAxis->FindBin(v0.pt()) - 1; + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { + if (v0.mK0Short() < cfgPIDCuts.cfgK0SignalMin) + npt_resonances[0][ptBinIndex] += getEfficiency(v0); + else if (v0.mK0Short() > cfgPIDCuts.cfgK0SignalMax) + npt_resonances[2][ptBinIndex] += getEfficiency(v0); + else + npt_resonances[1][ptBinIndex] += getEfficiency(v0); + registry.fill(HIST("K0/hK0s"), 1); + } + } + } + // Add lambdabar + if (resoSwitchVals[LAMBDA][kUseParticle]) { + if (selectLambda(collision, v0, centrality)) { + int ptBinIndex = fPtAxis->FindBin(v0.pt()) - 1; + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { + if (v0.mLambda() < cfgPIDCuts.cfgLambdaSignalMin) + npt_resonances[3][ptBinIndex] += getEfficiency(v0); + else if (v0.mLambda() > cfgPIDCuts.cfgLambdaSignalMax) + npt_resonances[5][ptBinIndex] += getEfficiency(v0); + else + npt_resonances[4][ptBinIndex] += getEfficiency(v0); + registry.fill(HIST("Lambda/hLambdas"), 1); + } + } + } + } + float chtotal = (cfgUseNchCorrection) ? acceptedTracks.total : acceptedTracks.total_uncorr; + // calculate fractions + std::vector> fractions_resonances = npt_resonances; + int pidcounter = 0; + for (const auto& vec : fractions_resonances) { + float total = chtotal; + if (cfgUsePIDTotal) + (pidcounter) ? std::accumulate(vec.begin(), vec.end(), 0) : chtotal; + + if (total == 0.) { + ++pidcounter; + continue; + } + std::transform(vec.begin(), vec.end(), + std::back_inserter(fractions_resonances.back()), + [&](double x) { return x / total; }); + ++pidcounter; + } + for (std::size_t i = 0; i < fractions_resonances[0].size(); ++i) + registry.fill(HIST("npt_K0_sb1"), fPtAxis->GetBinCenter(i + 1), centrality, fractions_resonances[0][i]); + for (std::size_t i = 0; i < fractions_resonances[2].size(); ++i) + registry.fill(HIST("npt_K0_sb2"), fPtAxis->GetBinCenter(i + 1), centrality, fractions_resonances[2][i]); + for (std::size_t i = 0; i < fractions_resonances[1].size(); ++i) + registry.fill(HIST("npt_K0_sig"), fPtAxis->GetBinCenter(i + 1), centrality, fractions_resonances[1][i]); + for (std::size_t i = 0; i < fractions_resonances[3].size(); ++i) + registry.fill(HIST("npt_Lambda_sb1"), fPtAxis->GetBinCenter(i + 1), centrality, fractions_resonances[3][i]); + for (std::size_t i = 0; i < fractions_resonances[5].size(); ++i) + registry.fill(HIST("npt_Lambda_sb2"), fPtAxis->GetBinCenter(i + 1), centrality, fractions_resonances[5][i]); + for (std::size_t i = 0; i < fractions_resonances[4].size(); ++i) + registry.fill(HIST("npt_Lambda_sig"), fPtAxis->GetBinCenter(i + 1), centrality, fractions_resonances[4][i]); + for (uint l_ind = 4; l_ind < corrconfigsradial.size(); ++l_ind) { + for (int i = 1; i <= fPtAxis->GetNbins(); i++) { + auto dnx = fGFW->Calculate(corrconfigsradial.at(l_ind), i - 1, kTRUE).real(); + if (dnx == 0) + continue; + auto val = fGFW->Calculate(corrconfigsradial.at(l_ind), i - 1, kFALSE).real() / dnx; + if (std::abs(val) < 1) + (dt == kGen) ? fFCgen->FillProfile(Form("%s_pt_%i", corrconfigsradial.at(l_ind).Head.c_str(), i), centrality, val * fractions_resonances[l_ind - 4][i - 1], dnx, lRandom) : fFC->FillProfile(Form("%s_pt_%i", corrconfigsradial.at(l_ind).Head.c_str(), i), centrality, val * fractions_resonances[l_ind - 4][i - 1], dnx, lRandom); + } + } int multiplicity = 0; switch (cfgUseNchCorrection) { case 0: @@ -1089,6 +1457,20 @@ struct FlowGenericFramework { pidIndex = 3; } + if (pidIndex) + acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track); + int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; + + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { + acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; + if (pidIndex == 1) + acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; + if (pidIndex == 2) + acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; + if (pidIndex == 3) + acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; + } + if (cfgFillWeights) { fillWeights(mcParticle, vtxz, 0, run); } else { @@ -1124,6 +1506,21 @@ struct FlowGenericFramework { } ++acceptedTracks.total; ++acceptedTracks.total_uncorr; + + if (pidIndex) + acceptedTracks.pidtotal[pidIndex - 1] += 1; + int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; + + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { + acceptedTracks.nch[ptBinIndex] += 1.0; + if (pidIndex == 1) + acceptedTracks.npi[ptBinIndex] += 1.0; + if (pidIndex == 2) + acceptedTracks.nka[ptBinIndex] += 1.0; + if (pidIndex == 3) + acceptedTracks.npr[ptBinIndex] += 1.0; + } + fillPtSums(track, vtxz, pidIndex); fillGFW(track, vtxz, pidIndex, densitycorrections); @@ -1144,15 +1541,24 @@ struct FlowGenericFramework { // if (cfgUsePID) Need PID for v02 int pidIndex = getNsigmaPID(track); - std::size_t ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; - if (!(ptBinIndex > o2::analysis::gfw::ptbinning.size())) { + if (pidIndex) + acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track); + + if (pidIndex) + acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track); + int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; + + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - if (pidIndex == 1) + if (pidIndex == 1) { acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - if (pidIndex == 2) + } + if (pidIndex == 2) { acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - if (pidIndex == 3) + } + if (pidIndex == 3) { acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; + } } if (cfgFillWeights) { @@ -1171,6 +1577,177 @@ struct FlowGenericFramework { } } + template + bool selectionV0Daughter(TTrack const& track, int pid) + { + if (!(track.itsNCls() > cfgTrackCuts.cfgMinNITSCls)) + return 0; + if (!track.hasTPC()) + return false; + if (track.tpcNClsFound() < cfgTrackCuts.cfgNTPCCls) + return false; + if (!(track.tpcNClsCrossedRows() > cfgTrackCuts.cfgNTPCXrows)) + return 0; + + if (cfgPIDCuts.cfgUseOnlyTPC) { + if (pid == PIONS && std::abs(track.tpcNSigmaPi()) > cfgPIDCuts.cfgTPCNsigmaCut) + return false; + if (pid == KAONS && std::abs(track.tpcNSigmaKa()) > cfgPIDCuts.cfgTPCNsigmaCut) + return false; + if (pid == PROTONS && std::abs(track.tpcNSigmaPr()) > cfgPIDCuts.cfgTPCNsigmaCut) + return false; + } else { + int partIndex = cfgPIDCuts.cfgUseAsymmetricPID ? getNsigmaPIDAssymmetric(track) : getNsigmaPID(track); + int pidIndex = partIndex - 1; // 0 = pion, 1 = kaon, 2 = proton + if (pidIndex != pid) + return false; + } + + return true; + } + + template + bool selectK0(TCollision const& collision, TV0 const& v0, const double& centrality) + { + + double mass_K0s = v0.mK0Short(); + + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + + registry.fill(HIST("K0/hK0Count"), kFillCandidate); + if (postrack.pt() < resoCutVals[K0][kPosTrackPt] || negtrack.pt() < resoCutVals[K0][kNegTrackPt]) + return false; + registry.fill(HIST("K0/hK0Count"), kFillDaughterPt); + if (mass_K0s < resoCutVals[K0][kMassMin] && mass_K0s > resoCutVals[K0][kMassMax]) + return false; + registry.fill(HIST("K0/hK0Count"), kFillMassCut); + // Rapidity correction + if (v0.yK0Short() > resoCutVals[K0][kRapidity]) + return false; + registry.fill(HIST("K0/hK0Count"), kFillRapidityCut); + // DCA cuts for K0short + if (std::abs(v0.dcapostopv()) < resoCutVals[K0][kDCAPosToPVMin] || std::abs(v0.dcanegtopv()) < resoCutVals[K0][kDCANegToPVMin]) + return false; + registry.fill(HIST("K0/hK0Count"), kFillDCAtoPV); + if (std::abs(v0.dcaV0daughters()) > resoSwitchVals[K0][kDCABetDaug]) + return false; + registry.fill(HIST("K0/hK0Count"), kFillDCAxDaughters); + // v0 radius cuts + if (resoSwitchVals[K0][kUseV0Radius] && (v0.v0radius() < resoCutVals[K0][kRadiusMin] || v0.v0radius() > resoCutVals[K0][kRadiusMax])) + return false; + registry.fill(HIST("K0/hK0Count"), kFillV0Radius); + // cosine pointing angle cuts + if (v0.v0cosPA() < resoCutVals[K0][kCosPA]) + return false; + registry.fill(HIST("K0/hK0Count"), kFillCosPA); + // Proper lifetime + if (resoSwitchVals[K0][kUseProperLifetime] && v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massK0Short > resoCutVals[K0][kLifeTime]) + return false; + registry.fill(HIST("K0/hK0Count"), kFillProperLifetime); + if (!selectionV0Daughter(postrack, PIONS) || !selectionV0Daughter(negtrack, PIONS)) + return false; + registry.fill(HIST("K0/hK0Count"), kFillDaughterTrackSelection); + + registry.fill(HIST("K0/hK0Mass_sparse"), mass_K0s, v0.pt(), centrality); + registry.fill(HIST("K0/hK0Phi"), v0.phi()); + registry.fill(HIST("K0/hK0Eta"), v0.eta()); + registry.fill(HIST("K0/PiPlusTPC_K0"), postrack.pt(), postrack.tpcNSigmaKa()); + registry.fill(HIST("K0/PiPlusTOF_K0"), postrack.pt(), postrack.tofNSigmaKa()); + registry.fill(HIST("K0/PiMinusTPC_K0"), negtrack.pt(), negtrack.tpcNSigmaKa()); + registry.fill(HIST("K0/PiMinusTOF_K0"), negtrack.pt(), negtrack.tofNSigmaKa()); + + return true; + } + + template + bool selectLambda(TCollision const& collision, TV0 const& v0, const double& centrality) + { + bool isL = false; // Is lambda candidate + bool isAL = false; // Is anti-lambda candidate + + double mlambda = v0.mLambda(); + double mantilambda = v0.mAntiLambda(); + + // separate the positive and negative V0 daughters + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + + registry.fill(HIST("Lambda/hLambdaCount"), kFillCandidate); + if (postrack.pt() < resoCutVals[LAMBDA][kPosTrackPt] || negtrack.pt() < resoCutVals[LAMBDA][kNegTrackPt]) + return false; + + registry.fill(HIST("Lambda/hLambdaCount"), kFillDaughterPt); + if (mlambda > resoCutVals[LAMBDA][kMassMin] && mlambda < resoCutVals[LAMBDA][kMassMax]) + isL = true; + if (mantilambda > resoCutVals[LAMBDA][kMassMin] && mantilambda < resoCutVals[LAMBDA][kMassMax]) + isAL = true; + + if (!isL && !isAL) { + return false; + } + registry.fill(HIST("Lambda/hLambdaCount"), kFillMassCut); + + // Rapidity correction + if (v0.yLambda() > resoCutVals[LAMBDA][kRapidity]) + return false; + registry.fill(HIST("Lambda/hLambdaCount"), kFillRapidityCut); + // DCA cuts for lambda and antilambda + if (isL) { + if (std::abs(v0.dcapostopv()) < resoCutVals[LAMBDA][kDCAPosToPVMin] || std::abs(v0.dcanegtopv()) < resoCutVals[LAMBDA][kDCANegToPVMin]) + return false; + } + if (isAL) { + if (std::abs(v0.dcapostopv()) < resoCutVals[LAMBDA][kDCANegToPVMin] || std::abs(v0.dcanegtopv()) < resoCutVals[LAMBDA][kDCAPosToPVMin]) + return false; + } + registry.fill(HIST("Lambda/hLambdaCount"), kFillDCAtoPV); + if (std::abs(v0.dcaV0daughters()) > resoSwitchVals[LAMBDA][kDCABetDaug]) + return false; + registry.fill(HIST("Lambda/hLambdaCount"), kFillDCAxDaughters); + // v0 radius cuts + if (resoSwitchVals[LAMBDA][kUseV0Radius] && (v0.v0radius() < resoCutVals[LAMBDA][kRadiusMin] || v0.v0radius() > resoCutVals[LAMBDA][kRadiusMax])) + return false; + registry.fill(HIST("Lambda/hLambdaCount"), kFillV0Radius); + // cosine pointing angle cuts + if (v0.v0cosPA() < resoCutVals[LAMBDA][kCosPA]) + return false; + registry.fill(HIST("Lambda/hLambdaCount"), kFillCosPA); + // Proper lifetime + if (resoSwitchVals[LAMBDA][kUseProperLifetime] && v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massLambda > resoCutVals[LAMBDA][kLifeTime]) + return false; + registry.fill(HIST("Lambda/hLambdaCount"), kFillProperLifetime); + if (isL) { + if (!selectionV0Daughter(postrack, PROTONS) || !selectionV0Daughter(negtrack, PIONS)) + return false; + } + if (isAL) { + if (!selectionV0Daughter(postrack, PIONS) || !selectionV0Daughter(negtrack, PROTONS)) + return false; + } + registry.fill(HIST("Lambda/hLambdaCount"), kFillDaughterTrackSelection); + + if (isL) { + registry.fill(HIST("Lambda/hLambdaMass_sparse"), mlambda, v0.pt(), centrality); + registry.fill(HIST("Lambda/hLambdaPhi"), v0.phi()); + registry.fill(HIST("Lambda/hLambdaEta"), v0.eta()); + registry.fill(HIST("Lambda/PrPlusTPC_L"), postrack.pt(), postrack.tpcNSigmaKa()); + registry.fill(HIST("Lambda/PrPlusTOF_L"), postrack.pt(), postrack.tofNSigmaKa()); + registry.fill(HIST("Lambda/PiMinusTPC_L"), negtrack.pt(), negtrack.tpcNSigmaKa()); + registry.fill(HIST("Lambda/PiMinusTOF_L"), negtrack.pt(), negtrack.tofNSigmaKa()); + } + if (isAL) { + registry.fill(HIST("Lambda/hAntiLambdaMass_sparse"), mantilambda, v0.pt(), centrality); + registry.fill(HIST("Lambda/hAntiLambdaPhi"), v0.phi()); + registry.fill(HIST("Lambda/hAntiLambdaEta"), v0.eta()); + registry.fill(HIST("Lambda/PiPlusTPC_AL"), postrack.pt(), postrack.tpcNSigmaKa()); + registry.fill(HIST("Lambda/PiPlusTOF_AL"), postrack.pt(), postrack.tofNSigmaKa()); + registry.fill(HIST("Lambda/PrMinusTPC_AL"), negtrack.pt(), negtrack.tpcNSigmaKa()); + registry.fill(HIST("Lambda/PrMinusTOF_AL"), negtrack.pt(), negtrack.tofNSigmaKa()); + } + return true; + } + template inline void fillPtSums(TTrack track, const double& vtxz, const int pidIndex) { @@ -1311,10 +1888,19 @@ struct FlowGenericFramework { o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ; o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::itsChi2NCl < cfgTrackCuts.cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgTrackCuts.cfgDCAz; + using GFWCollisions = soa::Filtered>; // using GFWTracks = soa::Filtered>; using GFWTracks = soa::Filtered>; - void processData(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) + SliceCache cache; + Partition posTracks = aod::track::signed1Pt > 0.0f; + Partition negTracks = aod::track::signed1Pt < 0.0f; + + double massKaPlus = o2::constants::physics::MassKPlus; + double massLambda = o2::constants::physics::MassLambda; + double massK0Short = o2::constants::physics::MassK0Short; + + void processData(GFWCollisions::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks, aod::V0Datas const& v0s) { auto bc = collision.bc_as(); int run = bc.runNumber(); @@ -1373,11 +1959,11 @@ struct FlowGenericFramework { // Get magnetic field polarity auto field = (cfgMagField == 99999) ? getMagneticField(bc.timestamp()) : cfgMagField; - processCollision(collision, tracks, centrality, field, run); + processCollision(collision, tracks, v0s, centrality, field, run); } PROCESS_SWITCH(FlowGenericFramework, processData, "Process analysis for non-derived data", true); - void processMCReco(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered> const& tracks, aod::McParticles const&) + void processMCReco(GFWCollisions::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered> const& tracks, aod::McParticles const&, aod::V0Datas const& v0s) { auto bc = collision.bc_as(); int run = bc.runNumber(); @@ -1422,12 +2008,12 @@ struct FlowGenericFramework { loadCorrections(bc); auto field = (cfgMagField == 99999) ? getMagneticField(bc.timestamp()) : cfgMagField; - processCollision(collision, tracks, centrality, field, run); + processCollision(collision, tracks, v0s, centrality, field, run); } PROCESS_SWITCH(FlowGenericFramework, processMCReco, "Process analysis for MC reconstructed events", false); o2::framework::expressions::Filter mcCollFilter = nabs(aod::mccollision::posZ) < cfgVtxZ; - void processMCGen(soa::Filtered::iterator const& mcCollision, soa::SmallGroups> const& collisions, aod::McParticles const& particles) + void processMCGen(soa::Filtered::iterator const& mcCollision, soa::SmallGroups> const& collisions, aod::McParticles const& particles, aod::V0Datas const& v0s) { if (collisions.size() != 1) return; @@ -1436,19 +2022,19 @@ struct FlowGenericFramework { centrality = getCentrality(collision); } int run = 0; - processCollision(mcCollision, particles, centrality, -999, run); + processCollision(mcCollision, particles, v0s, centrality, -999, run); } PROCESS_SWITCH(FlowGenericFramework, processMCGen, "Process analysis for MC generated events", false); - void processOnTheFly(soa::Filtered::iterator const& mcCollision, aod::McParticles const& mcParticles) + void processOnTheFly(soa::Filtered::iterator const& mcCollision, aod::McParticles const& mcParticles, aod::V0Datas const& v0s) { int run = 0; registry.fill(HIST("MCGen/impactParameter"), mcCollision.impactParameter(), mcParticles.size()); - processCollision(mcCollision, mcParticles, mcCollision.impactParameter(), -999, run); + processCollision(mcCollision, mcParticles, v0s, mcCollision.impactParameter(), -999, run); } PROCESS_SWITCH(FlowGenericFramework, processOnTheFly, "Process analysis for MC on-the-fly generated events", false); - void processRun2(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) + void processRun2(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks, aod::V0Datas const& v0s) { auto bc = collision.bc_as(); int run = bc.runNumber(); @@ -1464,7 +2050,7 @@ struct FlowGenericFramework { loadCorrections(bc); auto field = (cfgMagField == 99999) ? getMagneticField(bc.timestamp()) : cfgMagField; - processCollision(collision, tracks, centrality, field, run); + processCollision(collision, tracks, v0s, centrality, field, run); } PROCESS_SWITCH(FlowGenericFramework, processRun2, "Process analysis for Run 2 converted data", false); };