diff --git a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx index 811639344e5..c0211f56dbc 100644 --- a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx +++ b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx @@ -1,4 +1,4 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2024 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. // @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// /// \author Veronika Barbasova (veronika.barbasova@cern.ch) -/// \since October 12, 2023 +/// \since April 3, 2024 #include @@ -35,7 +35,12 @@ struct phianalysisTHnSparse { SliceCache cache; + Configurable produceQA{"produce-qa", false, "Produce QA histograms."}; Configurable produceTrue{"produce-true", false, "Produce True and Gen histograms."}; + Configurable produceLikesign{"produce-likesign", false, "Produce Like sign histograms."}; + Configurable eventMixing{"event-mixing", false, "Produce Event Mixing histograms."}; + + Configurable numberofMixedEvents{"number-of-mixed-events", 5, "Number of events that should be mixed."}; Configurable verboselevel{"verbose-level", 0, "Verbose level"}; Configurable refresh{"print-refresh", 0, "Freqency of print event information."}; Configurable refresh_index{"print-refresh-index", 0, "Freqency of print event information index."}; @@ -44,12 +49,17 @@ struct phianalysisTHnSparse { Configurable tpcnSigmaNeg{"tpc-ns-neg", 3.0f, "TPC NSigma cut of the negative particle."}; Configurable dautherPos{"dauther-type-pos", 3, "Particle type of the positive dauther according to ReconstructionDataFormats/PID.h (Default = Kaon)"}; Configurable dautherNeg{"dauther-type-neg", 3, "Particle type of the negative dauther according to ReconstructionDataFormats/PID.h (Default = Kaon)"}; - Configurable zVertexCut{"zvertex-cut", 10.0f, "Z vertex range."}; - Configurable rapidityCut{"rapidity-max", 0.5, "Rapidity cut."}; Configurable motherPDG{"mother-pdg", 333, "PDG code of mother particle."}; Configurable dautherPosPDG{"dauther-pdg-pos", 321, "PDG code of positive dauther particle."}; Configurable dautherNegPDG{"dauther-pdg-neg", 321, "PDG code of negative dauther particle."}; + // Cuts + Configurable zVertexCut{"zvertex-cut", 10.0f, "Z vertex range."}; + Configurable rapidityCut{"rapidity-max", 0.5, "Rapidity cut."}; + Configurable ptCut{"ptCut", 0.15f, "Cut: Minimal value of tracks pt."}; + Configurable dcaXYCut{"dcaXYCut", 1.0f, "Cut: Maximal value of tracks DCA XY."}; + Configurable dcaZCut{"dcaZCut", 1.0f, "Cut: Maximal value of tracks DCA Z."}; + Configurable> sparseAxes{"sparse-axes", std::vector{o2::analysis::rsn::PariAxis::names}, "Axes."}; ConfigurableAxis invaxis{"invAxis", {130, 0.97, 1.1}, "Invariant mass axis binning."}; @@ -60,39 +70,57 @@ struct phianalysisTHnSparse { ConfigurableAxis nsigmaaxisPos{"nsigmaAxisPos", {1, 0., tpcnSigmaPos}, "NSigma of positive particle axis binning in THnSparse."}; ConfigurableAxis nsigmaaxisNeg{"nsigmaAxisNeg", {1, 0., tpcnSigmaNeg}, "NSigma of negative particle axis binning in THnSparse."}; + // mixing + using BinningType = ColumnBinningPolicy>; + ConfigurableAxis axisVertexMixing{"axisVertexMixing", {1, -10, 10}, "vertex axis for bin"}; + ConfigurableAxis axisMultiplicityMixing{"axisMultiplicityMixing", {10, 0, 5000}, "TPC multiplicity for bin"}; + BinningType binning{{axisVertexMixing, axisMultiplicityMixing}, true}; + // defined in DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h float massPos = o2::track::PID::getMass(dautherPos); float massNeg = o2::track::PID::getMass(dautherNeg); + // Axes specifications + AxisSpec posZAxis = {400, -20., 20., "Z vx (cm)"}; + AxisSpec dcaxyAxis = {120, -3.0, 3.0, "DCA_{xy} (cm)"}; + AxisSpec dcazAxis = {120, -3.0, 3.0, "DCA_{z} (cm)"}; + + // Sparse axes AxisSpec invAxis = {invaxis, "Inv. mass (GeV/c^{2})", "im"}; AxisSpec ptAxis = {ptaxis, "p_{T} (GeV/c)", "pt"}; - AxisSpec mAxis = {multiplicityaxis, "N", "mu"}; + AxisSpec mu1Axis = {multiplicityaxis, "N", "mu1"}; + AxisSpec mu2Axis = {multiplicityaxis, "N", "mu2"}; AxisSpec yAxis = {rapidityaxis, "y", "y"}; AxisSpec nsigmatrackaxisPos = {nsigmaaxisPos, fmt::format("nSigma of positive particle ({})", massPos), "ns1"}; AxisSpec nsigmatrackaxisNeg = {nsigmaaxisNeg, fmt::format("nSigma of negative particle ({})", massNeg), "ns2"}; + AxisSpec zv1Axis = {posZ, "Z vx (cm)", "zv1"}; + AxisSpec zv2Axis = {posZ, "Z vx (cm)", "zv2"}; // All axes has to have same order as defined enum o2::analysis::rsn::PairAxisType (name from AxisSpec is taken to compare in o2::analysis::rsn::Output::init()) - std::vector allAxes = {invAxis, ptAxis, mAxis, nsigmatrackaxisPos, nsigmatrackaxisNeg, yAxis}; + std::vector allAxes = {invAxis, ptAxis, mu1Axis, nsigmatrackaxisPos, nsigmatrackaxisNeg, yAxis, zv1Axis, mu2Axis, zv2Axis}; HistogramRegistry registry{"registry"}; o2::analysis::rsn::Output* rsnOutput = nullptr; Service pdg; - float multiplicity; - float multiplicityMC; + int n = 0; double* pointPair = nullptr; double* pointEvent = nullptr; TLorentzVector d1, d2, mother; - Filter eventFilter = (o2::aod::evsel::sel8 == true); + Filter triggerFilter = (o2::aod::evsel::sel8 == true); Filter vtxFilter = (nabs(o2::aod::collision::posZ) < zVertexCut); + Filter ptFilter = nabs(aod::track::pt) > ptCut; + Filter etaFilter = nabs(aod::track::eta) < rapidityCut; + Filter dcaFilter = (nabs(o2::aod::track::dcaXY) < dcaXYCut) && (nabs(o2::aod::track::dcaZ) < dcaZCut); + using EventCandidates = soa::Filtered>; using EventCandidate = EventCandidates::iterator; - using TrackCandidates = soa::Join; + using TrackCandidates = soa::Filtered>; using EventCandidatesMC = soa::Filtered>; - using TrackCandidatesMC = soa::Join; + using TrackCandidatesMC = soa::Filtered>; Partition positive = (aod::track::signed1Pt > 0.0f) && (nabs(o2::aod::pidtpc::tpcNSigmaKa) < tpcnSigmaPos); Partition negative = (aod::track::signed1Pt < 0.0f) && (nabs(o2::aod::pidtpc::tpcNSigmaKa) < tpcnSigmaNeg); @@ -105,7 +133,27 @@ struct phianalysisTHnSparse { pointPair = new double[static_cast(o2::analysis::rsn::PairAxisType::unknown)]; pointEvent = new double[static_cast(o2::analysis::rsn::EventType::all)]; rsnOutput = new o2::analysis::rsn::OutputSparse(); - rsnOutput->init(sparseAxes, allAxes, produceTrue, ®istry); + rsnOutput->init(sparseAxes, allAxes, produceTrue, eventMixing, produceLikesign, ®istry); + + if (produceQA) { + // Event QA + registry.add("QAEvent/hVtxZ", "", kTH1F, {posZAxis}); + // Track QA + registry.add("QATrack/unlikepm/beforeSelection/hTrack1pt", "", kTH1F, {ptAxis}); + registry.add("QATrack/unlikepm/beforeSelection/hTrackDCAxy", "", kTH1F, {dcaxyAxis}); + registry.add("QATrack/unlikepm/beforeSelection/hTrackDCAz", "", kTH1F, {dcazAxis}); + registry.add("QATrack/unlikepm/beforeSelection/hTrack1eta", "", kTH1F, {{100, -1.0, 1.0}}); + + registry.add("QATrack/unlikepm/afterSelection/hTrack1pt", "", kTH1F, {ptAxis}); + registry.add("QATrack/unlikepm/afterSelection/hTrackDCAxy", "", kTH1F, {dcaxyAxis}); + registry.add("QATrack/unlikepm/afterSelection/hTrackDCAz", "", kTH1F, {dcazAxis}); + registry.add("QATrack/unlikepm/afterSelection/hTrack1eta", "", kTH1F, {{100, -1.0, 1.0}}); + + registry.add("QATrack/unlikepm/TPCPID/h2TracknSigma", "", kTH2F, {{120, -tpcnSigmaPos, tpcnSigmaPos}, {120, -tpcnSigmaNeg, tpcnSigmaNeg}}); + + // Mixing QA + registry.add("QAMixing/s4Multiplicity", "", kTHnSparseF, {axisMultiplicityMixing, axisMultiplicityMixing, axisVertexMixing, axisVertexMixing}); + } } template @@ -113,6 +161,8 @@ struct phianalysisTHnSparse { { if (!track.isPrimaryTrack()) return false; + if (!track.isPVContributor()) + return false; return true; } @@ -128,101 +178,134 @@ struct phianalysisTHnSparse { return true; } + template + float GetMultiplicity(const T& collision) + { + float multiplicity = collision.multFT0C() + collision.multFT0A(); + return multiplicity; + } + void processData(EventCandidate const& collision, TrackCandidates const& tracks) { + float multiplicity; + auto posDauthers = positive->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); auto negDauthers = negative->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); if (ignorezeroevent && collision.globalIndex() == 0) { if (verboselevel > 0) LOGF(info, "BAD pos=%lld neg=%lld, Z vertex position: %f [cm], %d, mult:%f.0", posDauthers.size(), negDauthers.size(), collision.posZ(), - collision.globalIndex(), multiplicity); + collision.globalIndex(), GetMultiplicity(collision)); return; } - multiplicity = collision.multFT0A() + collision.multFT0C(); + multiplicity = GetMultiplicity(collision); if (verboselevel > 0 && refresh > 0 && collision.globalIndex() % refresh == refresh_index) LOGF(info, "pos=%lld neg=%lld, Z vertex position: %f [cm], %d, mult:%f.0", posDauthers.size(), negDauthers.size(), collision.posZ(), collision.globalIndex(), multiplicity); - if (std::abs(collision.posZ()) > zVertexCut) - return; + if (produceQA) + registry.fill(HIST("QAEvent/hVtxZ"), collision.posZ()); pointEvent[0] = collision.posZ(); rsnOutput->fill(o2::analysis::rsn::EventType::zvertex, pointEvent); for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthers, negDauthers))) { - + if (produceQA) { + registry.fill(HIST("QATrack/unlikepm/beforeSelection/hTrack1pt"), track1.pt()); + registry.fill(HIST("QATrack/unlikepm/beforeSelection/hTrackDCAxy"), track1.dcaXY()); + registry.fill(HIST("QATrack/unlikepm/beforeSelection/hTrackDCAz"), track1.dcaZ()); + registry.fill(HIST("QATrack/unlikepm/beforeSelection/hTrack1eta"), track1.eta()); + } if (!selectedTrack(track1)) - continue; + if (!selectedTrack(track2)) continue; + if (produceQA) { + registry.fill(HIST("QATrack/unlikepm/afterSelection/hTrack1pt"), track1.pt()); + registry.fill(HIST("QATrack/unlikepm/afterSelection/hTrackDCAxy"), track1.dcaXY()); + registry.fill(HIST("QATrack/unlikepm/afterSelection/hTrackDCAz"), track1.dcaZ()); + registry.fill(HIST("QATrack/unlikepm/afterSelection/hTrack1eta"), track1.eta()); + } + if (!selectedPair(mother, track1, track2)) continue; + if (produceQA) + registry.fill(HIST("QATrack/unlikepm/TPCPID/h2TracknSigma"), track1.tpcNSigmaKa(), track2.tpcNSigmaKa()); + if (verboselevel > 1) LOGF(info, "Unlike-sign: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.Mag()); pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = multiplicity; pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = collision.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = collision.posZ(); - rsnOutput->fillUnlike(pointPair); + rsnOutput->fillUnlikepm(pointPair); } - for (auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posDauthers, posDauthers))) { - - if (!selectedTrack(track1)) - continue; - if (!selectedTrack(track2)) - continue; + if (produceLikesign) { - if (!selectedPair(mother, track1, track2)) - continue; + for (auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posDauthers, posDauthers))) { + if (!selectedTrack(track1)) + continue; + if (!selectedTrack(track2)) + continue; - if (verboselevel > 1) - LOGF(info, "Like-sign positive: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.Mag()); + if (!selectedPair(mother, track1, track2)) + continue; - pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu)] = multiplicity; - pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + if (verboselevel > 1) + LOGF(info, "Like-sign positive: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.Mag()); - rsnOutput->fillLikepp(pointPair); - } + pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = collision.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = collision.posZ(); - for (auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negDauthers, negDauthers))) { + rsnOutput->fillLikepp(pointPair); + } - if (!selectedTrack(track1)) - continue; - if (!selectedTrack(track2)) - continue; + for (auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negDauthers, negDauthers))) { + if (!selectedTrack(track1)) + continue; + if (!selectedTrack(track2)) + continue; - if (!selectedPair(mother, track1, track2)) - continue; + if (!selectedPair(mother, track1, track2)) + continue; - if (verboselevel > 1) - LOGF(info, "Like-sign negative: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.Mag()); + if (verboselevel > 1) + LOGF(info, "Like-sign negative: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.Mag()); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu)] = multiplicity; - pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = collision.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = collision.posZ(); - rsnOutput->fillLikemm(pointPair); + rsnOutput->fillLikemm(pointPair); + } } } - PROCESS_SWITCH(phianalysisTHnSparse, processData, "Process Event for Data", true); void processTrue(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions) @@ -230,6 +313,8 @@ struct phianalysisTHnSparse { if (!produceTrue) return; + float multiplicity; + auto posDauthersMC = positiveMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); auto negDauthersMC = negativeMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); @@ -241,9 +326,9 @@ struct phianalysisTHnSparse { if (std::abs(collision.posZ()) > zVertexCut) return; - multiplicityMC = collision.multFT0A() + collision.multFT0C(); + multiplicity = GetMultiplicity(collision); - for (auto& [track1, track2] : combinations(o2::soa::CombinationsUpperIndexPolicy(posDauthersMC, negDauthersMC))) { + for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersMC, negDauthersMC))) { if (!track1.has_mcParticle()) { LOGF(warning, "No MC particle for track, skip..."); @@ -295,10 +380,13 @@ struct phianalysisTHnSparse { pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = multiplicity; pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = collision.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = collision.posZ(); rsnOutput->fillUnliketrue(pointPair); } @@ -312,6 +400,8 @@ struct phianalysisTHnSparse { void processGen(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles) { + float multiplicity = 0; + if (!produceTrue) return; @@ -351,10 +441,13 @@ struct phianalysisTHnSparse { pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); - pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu)] = multiplicityMC; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = multiplicity; pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(tpcnSigmaPos / 2.0); pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(tpcnSigmaNeg / 2.0); pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = mcCollision.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = multiplicity; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = mcCollision.posZ(); rsnOutput->fillUnlikegen(pointPair); @@ -368,6 +461,135 @@ struct phianalysisTHnSparse { } PROCESS_SWITCH(phianalysisTHnSparse, processGen, "Process generated.", false); + + double posz1; + double posz2; + + void processMixed(EventCandidates const& collisions, TrackCandidates const& tracks) + { + if (!eventMixing) + return; + + auto tracksTuple = std::make_tuple(tracks); + + SameKindPair pair{binning, numberofMixedEvents, -1, collisions, tracksTuple, &cache}; + for (auto& [c1, tracks1, c2, tracks2] : pair) { + if (!c1.sel8()) { + continue; + } + if (!c2.sel8()) { + continue; + } + + auto posDauthersc1 = positive->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); + auto posDauthersc2 = positive->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); + auto negDauthersc1 = negative->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); + auto negDauthersc2 = negative->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); + + for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc1, negDauthersc2))) { + + if (!selectedTrack(track1)) + + continue; + if (!selectedTrack(track2)) + continue; + + if (!selectedPair(mother, track1, track2)) + continue; + + pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = GetMultiplicity(c1); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = posz1; + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = GetMultiplicity(c2); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = posz2; + + rsnOutput->fillMixingpm(pointPair); + + if (produceQA) + registry.fill(HIST("QAMixing/s4Multiplicity"), GetMultiplicity(c1), GetMultiplicity(c2), c1.posZ(), c2.posZ()); + } + + if (produceLikesign) { + + for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc1, posDauthersc2))) { + + if (!selectedTrack(track1)) + + continue; + if (!selectedTrack(track2)) + continue; + + if (!selectedPair(mother, track1, track2)) + continue; + + pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = GetMultiplicity(c1); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = c1.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = GetMultiplicity(c2); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = c2.posZ(); + + rsnOutput->fillMixingpp(pointPair); + } + + for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(negDauthersc1, negDauthersc2))) { + + if (!selectedTrack(track1)) + + continue; + if (!selectedTrack(track2)) + continue; + + if (!selectedPair(mother, track1, track2)) + continue; + + pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = GetMultiplicity(c1); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = c1.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = GetMultiplicity(c2); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = c2.posZ(); + + rsnOutput->fillMixingmm(pointPair); + } + } + + for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc2, negDauthersc1))) { + + if (!selectedTrack(track1)) + + continue; + if (!selectedTrack(track2)) + continue; + + if (!selectedPair(mother, track1, track2)) + continue; + + pointPair[static_cast(o2::analysis::rsn::PairAxisType::im)] = mother.Mag(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::pt)] = mother.Pt(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu1)] = GetMultiplicity(c1); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns1)] = std::abs(track1.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::ns2)] = std::abs(track2.tpcNSigmaKa()); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::y)] = mother.Rapidity(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv1)] = c1.posZ(); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::mu2)] = GetMultiplicity(c2); + pointPair[static_cast(o2::analysis::rsn::PairAxisType::zv2)] = c2.posZ(); + + rsnOutput->fillMixingmp(pointPair); + } + } + } + PROCESS_SWITCH(phianalysisTHnSparse, processMixed, "Process Mixing Event.", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGLF/Utils/rsnOutput.h b/PWGLF/Utils/rsnOutput.h index e2fd9845bea..d04c97ae099 100644 --- a/PWGLF/Utils/rsnOutput.h +++ b/PWGLF/Utils/rsnOutput.h @@ -1,4 +1,4 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2024 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. // @@ -8,20 +8,15 @@ // 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 rsnOutput.h -/// \brief Dynamic sparse output for resonances for O2 analysis -/// -/// Simply copied from FemtoDreamCollisionSelection.h -/// original author: Laura Serksnyte, TU München /// -/// \author Veronika Barbasova +/// \author Veronika Barbasova (veronika.barbasova@cern.ch) +/// \since April 3, 2024 #ifndef PWGLF_UTILS_RSNOUTPUT_H_ #define PWGLF_UTILS_RSNOUTPUT_H_ -#include #include +#include #include #include "Framework/HistogramRegistry.h" @@ -46,33 +41,41 @@ enum class TrackType { }; enum class PairType { - unlike, + unlikepm, + unlikemp, likepp, likemm, unliketrue, unlikegen, + mixingpm, + mixingpp, + mixingmm, + mixingmp, all }; enum class PairAxisType { - im = 0, - pt = 1, - mu = 2, - ns1 = 3, - ns2 = 4, - y = 5, + im, + pt, + mu1, + ns1, + ns2, + y, + zv1, + mu2, + zv2, unknown }; namespace PariAxis { -std::vector names{"im", "pt", "mu", "ns1", "ns2", "y"}; +std::vector names{"im", "pt", "mu1", "ns1", "ns2", "y", "zv1", "mu2", "zv2"}; } class Output { public: virtual ~Output() = default; - virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, bool produceTrue = false, HistogramRegistry* registry = nullptr) + virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, bool produceTrue = false, bool eventMixing = false, bool produceLikesign = false, HistogramRegistry* registry = nullptr) { mHistogramRegistry = registry; if (mHistogramRegistry == nullptr) @@ -104,7 +107,7 @@ class Output mFillPoint = new double[mCurrentAxisTypes.size()]; LOGF(info, "Number of axis added: %d", mCurrentAxes.size()); - mPairHisto = new HistogramConfigSpec(HistType::kTHnSparseF, mCurrentAxes); + mPairHisto = new HistogramConfigSpec({HistType::kTHnSparseF}, mCurrentAxes); }; template @@ -130,11 +133,16 @@ class Output LOGF(warning, "Abstract method : 'virtual void rsn::Output::fill(PairType t, double* point)' !!! Please implement it first."); }; - virtual void fillUnlike(double* point) = 0; + virtual void fillUnlikepm(double* point) = 0; + virtual void fillUnlikemp(double* point) = 0; virtual void fillLikepp(double* point) = 0; virtual void fillLikemm(double* point) = 0; virtual void fillUnliketrue(double* point) = 0; virtual void fillUnlikegen(double* point) = 0; + virtual void fillMixingpm(double* point) = 0; + virtual void fillMixingpp(double* point) = 0; + virtual void fillMixingmm(double* point) = 0; + virtual void fillMixingmp(double* point) = 0; PairAxisType type(std::string name) { @@ -169,19 +177,29 @@ class Output class OutputSparse : public Output { public: - virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, bool produceTrue = false, HistogramRegistry* registry = nullptr) + virtual void init(std::vector const& sparseAxes, std::vector const& allAxes, bool produceTrue = false, bool eventMixing = false, bool produceLikesign = false, HistogramRegistry* registry = nullptr) { - Output::init(sparseAxes, allAxes, produceTrue, registry); + Output::init(sparseAxes, allAxes, produceTrue, eventMixing, produceLikesign, registry); mHistogramRegistry->add("hVz", "; vtx_{z} (cm); Entries", kTH1F, {{40, -20., 20.}}); - mHistogramRegistry->add("unlike", "Unlike", *mPairHisto); - - mHistogramRegistry->add("likepp", "Like PP", *mPairHisto); - mHistogramRegistry->add("likemm", "Like MM", *mPairHisto); + mHistogramRegistry->add("unlikepm", "Unlike pm", *mPairHisto); + mHistogramRegistry->add("unlikemp", "Unlike mp", *mPairHisto); + if (produceLikesign) { + mHistogramRegistry->add("likepp", "Like PP", *mPairHisto); + mHistogramRegistry->add("likemm", "Like MM", *mPairHisto); + } if (produceTrue) { mHistogramRegistry->add("unliketrue", "Unlike True", *mPairHisto); mHistogramRegistry->add("unlikegen", "Unlike Gen", *mPairHisto); } + if (eventMixing) { + mHistogramRegistry->add("mixingpm", "Event Mixing pm", *mPairHisto); + if (produceLikesign) { + mHistogramRegistry->add("mixingpp", "Event Mixing pp", *mPairHisto); + mHistogramRegistry->add("mixingmm", "Event Mixing mm", *mPairHisto); + } + mHistogramRegistry->add("mixingmp", "Event Mixing mp", *mPairHisto); + } } virtual void @@ -199,8 +217,11 @@ class OutputSparse : public Output virtual void fill(PairType t, double* point) { switch (t) { - case PairType::unlike: - fillUnlike(point); + case PairType::unlikepm: + fillUnlikepm(point); + break; + case PairType::unlikemp: + fillUnlikemp(point); break; case PairType::likepp: fillLikepp(point); @@ -214,16 +235,31 @@ class OutputSparse : public Output case PairType::unlikegen: fillUnlikegen(point); break; + case PairType::mixingpm: + fillMixingpm(point); + break; + case PairType::mixingpp: + fillMixingpp(point); + break; + case PairType::mixingmm: + fillMixingmm(point); + break; + case PairType::mixingmp: + fillMixingmp(point); + break; default: break; } } - virtual void fillUnlike(double* point) + virtual void fillUnlikepm(double* point) { - fillSparse(HIST("unlike"), point); + fillSparse(HIST("unlikepm"), point); + } + virtual void fillUnlikemp(double* point) + { + fillSparse(HIST("unlikemp"), point); } - virtual void fillLikepp(double* point) { fillSparse(HIST("likepp"), point); @@ -243,6 +279,22 @@ class OutputSparse : public Output { fillSparse(HIST("unlikegen"), point); } + virtual void fillMixingpm(double* point) + { + fillSparse(HIST("mixingpm"), point); + } + virtual void fillMixingpp(double* point) + { + fillSparse(HIST("mixingpp"), point); + } + virtual void fillMixingmm(double* point) + { + fillSparse(HIST("mixingmm"), point); + } + virtual void fillMixingmp(double* point) + { + fillSparse(HIST("mixingmp"), point); + } }; } // namespace rsn } // namespace o2::analysis