From 6e64f20ae1731086f9f05a6f42740b0e30b8cbb5 Mon Sep 17 00:00:00 2001 From: mjkim525 Date: Wed, 20 Aug 2025 22:58:17 +0900 Subject: [PATCH 1/6] added initializerevetqa.cxx and updated xi1530Analysisqa.cxx --- PWGLF/Tasks/Resonances/CMakeLists.txt | 5 + PWGLF/Tasks/Resonances/initializereventqa.cxx | 477 ++++++++++++++++++ PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx | 161 ++++-- 3 files changed, 601 insertions(+), 42 deletions(-) create mode 100644 PWGLF/Tasks/Resonances/initializereventqa.cxx diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index de23fa67237..27aad7c24c3 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -174,6 +174,11 @@ o2physics_add_dpl_workflow(xi1530analysisqa PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(initializereventqa + SOURCES initializereventqa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(kaonkaonanalysis SOURCES kaonkaonanalysis.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Resonances/initializereventqa.cxx b/PWGLF/Tasks/Resonances/initializereventqa.cxx new file mode 100644 index 00000000000..d3ddae10cbb --- /dev/null +++ b/PWGLF/Tasks/Resonances/initializereventqa.cxx @@ -0,0 +1,477 @@ +// Copyright 2019-2025 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +/// \file initializereventqa.cxx +/// \brief QA for the event-loss and signal-loss correction at the generator level for the ResonanceInitializer in pp collisions (referred to TableProducer/Strangeness/cascqaanalysis.cxx) +/// +/// Following the discussions at the two PAG meetings (https://indico.cern.ch/event/1518979, https://indico.cern.ch/event/1575984) +/// we have introduced an auxiliary task that, when the resonanceInitializer.cxx is used, +/// computes the event-loss and signal-loss correction factors at the generator level. +/// With minor configuration tuning for a truth-tagging, +/// we expect it to be applicable to most analyses that rely on the initializer. +/// +/// \author Minjae Kim (minjae.kim@cern.ch) + +#include +#include +#include + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/Centrality.h" +#include "PWGLF/DataModel/cascqaanalysis.h" +#include "TRandom2.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "PWGLF/Utils/inelGt.h" +#include "PWGLF/DataModel/mcCentrality.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + + +using TrkPidInfo = soa::Join; +using DauTracks = soa::Join; + +struct Initializereventqa { + + // Templates used, new hedder wil be added + Produces mycascades; + Produces myMCcascades; + + HistogramRegistry registry{"registry"}; + + // Axes + ConfigurableAxis ptAxis{"ptAxis", {400, 0.0f, 20.0f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis rapidityAxis{"rapidityAxis", {200, -2.0f, 2.0f}, "y"}; + ConfigurableAxis centFT0MAxis{"centFT0MAxis", + {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0, 110.0}, + "FT0M (%)"}; + ConfigurableAxis eventTypeAxis{"eventTypeAxis", {2, -0.5f, 1.5f}, "Event Type"}; + + ConfigurableAxis nAssocCollAxis{"nAssocCollAxis", {5, -0.5f, 4.5f}, "N_{assoc.}"}; + ConfigurableAxis nChargedFT0MGenAxis{"nChargedFT0MGenAxis", {300, 0, 300}, "N_{FT0M, gen.}"}; + ConfigurableAxis multNTracksAxis{"multNTracksAxis", {500, 0, 500}, "N_{tracks}"}; + ConfigurableAxis signalFT0MAxis{"signalFT0MAxis", {4000, 0, 40000}, "FT0M amplitude"}; + ConfigurableAxis signalFV0AAxis{"signalFV0AAxis", {4000, 0, 40000}, "FV0A amplitude"}; + ConfigurableAxis nCandidates{"nCandidates", {30, -0.5, 29.5}, "N_{cand.}"}; + + // Event selection criteria + Configurable cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"}; + Configurable isZvtxcut{"isZvtxcut", 1, "Select collisions with Accepted z-vertex"}; + Configurable isVertexITSTPC{"isVertexITSTPC", 0, "Select collisions with at least one ITS-TPC track"}; + Configurable isNoSameBunchPileup{"isNoSameBunchPileup", 0, "Same found-by-T0 bunch crossing rejection"}; + Configurable isGoodZvtxFT0vsPV{"isGoodZvtxFT0vsPV", 0, "z of PV by tracks and z of PV from FT0 A-C time difference cut"}; + Configurable isVertexTOFmatched{"isVertexTOFmatched", 0, "Is Vertex TOF matched"}; + + Configurable isTriggerTVX{"isTriggerTVX", 1, "TVX trigger"}; + Configurable isNoTimeFrameBorder{"isNoTimeFrameBorder", 1, "TF border cut"}; + Configurable isNoITSROFrameBorder{"isNoITSROFrameBorder", 1, "ITS ROF border cut"}; + Configurable isNoCollInTimeRangeNarrow{"isNoCollInTimeRangeNarrow", 0, "No collisions in +-2us window"}; + + // QA histograms for the multiplicity estimation + Configurable multQA{"multQA", 1, "0 - not to do QA, 1 - do the QA"}; + + // Selection for signal-loss corrections + Configurable isDaughterCheck{"isDaughterCheck", 1, "Check if the candidate has the correct daughters when it is considered"}; + + Configurable cfgRapidityCut{"cfgRapidityCut", 0.5, "Rapidity cut for the truth particle"}; + + Configurable pdgTruthMother{"pdgTruthMother", 3324, "pdgcode for the truth mother particle, e.g. Xi(1530) (3324)"}; + Configurable pdgTruthDaughter1{"pdgTruthDaughter1", 3312, "pdgcode for the first daughter particle, e.g. Xi-3312"}; + Configurable pdgTruthDaughter2{"pdgTruthDaughter2", 211, "pdgcode for the second daughter particle, e.g. Xi-3312"}; + + // Necessary for particle charges + Service pdgDB; + + SliceCache cache; + + // Struct to select on event type + typedef struct CollisionIndexAndType { + int64_t index; + uint8_t typeFlag; + } CollisionIndexAndType; + + void init(InitContext const&) + { + TString hNEventsMCLabels[5] = {"All", "z vrtx", "INEL", "INEL>0", "Associated with rec. collision"}; + TString hNEventsLabels[12] = {"All", "kIsTriggerTVX", "kNoTimeFrameBorder", "kNoITSROFrameBorder", "kIsVertexITSTPC", "kNoSameBunchPileup", "kIsGoodZvtxFT0vsPV", "isVertexTOFmatched", "kNoCollInTimeRangeNarrow", "z vrtx", "INEL", "INEL>0"}; + + registry.add("hNEvents", "hNEvents", {HistType::kTH1D, {{12, 0.f, 12.f}}}); + + for (int n = 1; n <= registry.get(HIST("hNEvents"))->GetNbinsX(); n++) { + registry.get(HIST("hNEvents"))->GetXaxis()->SetBinLabel(n, hNEventsLabels[n - 1]); + } + registry.add("hZCollision", "hZCollision", {HistType::kTH1D, {{200, -20.f, 20.f}}}); + + + registry.add("fakeEvents", "fakeEvents", {HistType::kTH1F, {{1, -0.5f, 0.5f}}}); + + registry.add("hNEventsMC", "hNEventsMC", {HistType::kTH1D, {{5, 0.0f, 5.0f}}}); + for (int n = 1; n <= registry.get(HIST("hNEventsMC"))->GetNbinsX(); n++) { + registry.get(HIST("hNEventsMC"))->GetXaxis()->SetBinLabel(n, hNEventsMCLabels[n - 1]); + } + registry.add("hZCollisionGen", "hZCollisionGen", {HistType::kTH1D, {{200, -20.f, 20.f}}}); + registry.add("hCentFT0MNAssocMCCollisions", "hCentFT0MNAssocMCCollisions", {HistType::kTH3D, {centFT0MAxis, nAssocCollAxis, eventTypeAxis}}); + registry.add("hCentFT0MNAssocMCCollisionsSameType", "hCentFT0MNAssocMCCollisionsSameType", {HistType::kTH3D, {centFT0MAxis, nAssocCollAxis, eventTypeAxis}}); + registry.add("hNchFT0MNAssocMCCollisions", "hNchFT0MNAssocMCCollisions", {HistType::kTH3D, {nChargedFT0MGenAxis, nAssocCollAxis, eventTypeAxis}}); + registry.add("hNchFT0MNAssocMCCollisionsSameType", "hNchFT0MNAssocMCCollisionsSameType", {HistType::kTH3D, {nChargedFT0MGenAxis, nAssocCollAxis, eventTypeAxis}}); + registry.add("hNContributorsCorrelation", "hNContributorsCorrelation", {HistType::kTH2F, {{250, -0.5f, 249.5f, "Secondary Contributor"}, {250, -0.5f, 249.5f, "Main Contributor"}}}); + registry.add("hNchFT0MGenEvType", "hNchFT0MGenEvType", {HistType::kTH2D, {nChargedFT0MGenAxis, eventTypeAxis}}); + registry.add("hCentFT0M_genMC", "hCentFT0M_genMC", {HistType::kTH2D, {centFT0MAxis, eventTypeAxis}}); + + registry.add("hCentFT0M_rec", "hCentFT0M_rec", {HistType::kTH2D, {centFT0MAxis, eventTypeAxis}}); + registry.add("hCentFT0M_corr", "hCentFT0M_Corr", {HistType::kTH2D, {centFT0MAxis, centFT0MAxis}}); + + + if (multQA) { + registry.add("hNchFT0Mglobal", "hNchFT0Mglobal", {HistType::kTH3D, {nChargedFT0MGenAxis, multNTracksAxis, eventTypeAxis}}); + registry.add("hNchFT0MPVContr", "hNchFT0MPVContr", {HistType::kTH3D, {nChargedFT0MGenAxis, multNTracksAxis, eventTypeAxis}}); + registry.add("hFT0MpvContr", "hFT0MpvContr", {HistType::kTH3D, {centFT0MAxis, multNTracksAxis, eventTypeAxis}}); + registry.add("hFT0Mglobal", "hFT0Mglobal", {HistType::kTH3D, {centFT0MAxis, multNTracksAxis, eventTypeAxis}}); + registry.add("hFT0MsignalPVContr", "hFT0MsignalPVContr", {HistType::kTH3D, {signalFT0MAxis, multNTracksAxis, eventTypeAxis}}); + } + + registry.add("h3ResonanceTruth", "pT distribution of True Resonance", kTHnSparseF, {eventTypeAxis, ptAxis, centFT0MAxis}); + registry.add("h3ResonanceTruthAnti", "pT distribution of True Resonance Anti", kTHnSparseF, {eventTypeAxis, ptAxis, centFT0MAxis}); + } +float pvEta1 = 1.0f; + float globalEta05 = 0.5f; + + Partition pvContribTracksIUEta1 = (nabs(aod::track::eta) < pvEta1) && ((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)); + Partition globalTracksIUEta05 = (nabs(aod::track::eta) < globalEta05) && (requireGlobalTrackInFilter()); + + + template + uint16_t getGenNchInFT0Mregion(TMcParticles particles) + { + float region1FT0 = -3.3f; + float region2FT0 = -2.1f; + float region3FT0 = 3.5f; + float region4FT0 = 4.9f; + // Particle counting in FITFT0: -3.3<η<-2.1; 3.5<η<4.9 + uint16_t nchFT0 = 0; + for (const auto& mcParticle : particles) { + if (!mcParticle.isPhysicalPrimary()) { + continue; + } + const auto& pdgInfo = pdgDB->GetParticle(mcParticle.pdgCode()); + if (!pdgInfo) { + continue; + } + if (pdgInfo->Charge() == 0) { + continue; + } + if (mcParticle.eta() < region1FT0 || mcParticle.eta() > region4FT0 || (mcParticle.eta() > region2FT0 && mcParticle.eta() < region3FT0)) { + continue; // select on T0M Nch region + } + nchFT0++; // increment + } + return nchFT0; + } + + template + bool acceptEvent(TCollision const& collision, bool isFillEventSelectionQA) + { + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 0.5); + } + + // kIsTriggerTVX selection + if (isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + return false; + } + + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 1.5); + } + + // kNoTimeFrameBorder selection + if (isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + return false; + } + + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 2.5); + } + + // kNoITSROFrameBorder selection + if (isNoITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { + return false; + } + + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 3.5); + } + + // kIsVertexITSTPC selection + if (isVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) { + return false; + } + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 4.5); + } + // kNoSameBunchPileup selection + if (isNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + return false; + } + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 5.5); + } + // kIsGoodZvtxFT0vsPV selection + if (isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; + } + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 6.5); + } + // isVertexTOFmatched selection + if (isVertexTOFmatched && !collision.selection_bit(aod::evsel::kIsVertexTOFmatched)) { + return false; + } + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 7.5); + } + // kNoCollInTimeRangeNarrow selection + if (isNoCollInTimeRangeNarrow && !collision.selection_bit(aod::evsel::kNoCollInTimeRangeNarrow)) { + return false; + } + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 8.5); + } + + // Z vertex selection + if (isZvtxcut && std::fabs(collision.posZ()) > cutzvertex) { + return false; + } + if (isFillEventSelectionQA) { + registry.fill(HIST("hNEvents"), 9.5); + registry.fill(HIST("hZCollision"), collision.posZ()); + } + + return true; + } + + + template + void fillMCParticles(TotalMCParts const& mcParticles, MultMCGen const& multiplicity, evtType const& eventType) + { + for (auto const& mcPart : mcParticles) { + + if (std::abs(mcPart.pdgCode()) != pdgTruthMother || std::abs(mcPart.y()) >= cfgRapidityCut) + continue; + std::vector daughterPDGs; + if (mcPart.has_daughters()) { + auto daughter01 = mcParticles.rawIteratorAt(mcPart.daughtersIds()[0] - mcParticles.offset()); + auto daughter02 = mcParticles.rawIteratorAt(mcPart.daughtersIds()[1] - mcParticles.offset()); + daughterPDGs = {daughter01.pdgCode(), daughter02.pdgCode()}; + } else { + daughterPDGs = {-1, -1}; + } + + if(isDaughterCheck){ + bool pass1 = std::abs(daughterPDGs[0]) == pdgTruthDaughter1 || std::abs(daughterPDGs[1]) == pdgTruthDaughter1; + bool pass2 = std::abs(daughterPDGs[0]) == pdgTruthDaughter2 || std::abs(daughterPDGs[1]) == pdgTruthDaughter2; + if (!pass1 || !pass2) + continue; + } + if (mcPart.pdgCode() > 0) // Consider INELt0 or INEL + registry.fill(HIST("h3ResonanceTruth"), eventType, mcPart.pt(), multiplicity); + else + registry.fill(HIST("h3ResonanceTruthAnti"), eventType, mcPart.pt(), multiplicity); + + daughterPDGs.clear(); + } + } + void processData(soa::Join::iterator const& collision, + DauTracks const&) + { + if (!acceptEvent(collision, 1)) { + return; + } + + int evType = 0; + registry.fill(HIST("hNEvents"), 10.5); // INEL + if (collision.isInelGt0()) { + evType += 1; + registry.fill(HIST("hNEvents"), 11.5); // INEL>0 + } + + auto tracksGroupedPVcontr = pvContribTracksIUEta1->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + int nTracksPVcontr = tracksGroupedPVcontr.size(); + + auto tracksGroupedGlobal = globalTracksIUEta05->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + int nTracksGlobal = tracksGroupedGlobal.size(); + + registry.fill(HIST("hCentFT0M_rec"), collision.centFT0M(), evType); + + if (multQA) { + registry.fill(HIST("hFT0MpvContr"), collision.centFT0M(), nTracksPVcontr, evType); + registry.fill(HIST("hFT0Mglobal"), collision.centFT0M(), nTracksGlobal, evType); + registry.fill(HIST("hFT0MsignalPVContr"), collision.multFT0A() + collision.multFT0C(), nTracksPVcontr, evType); + } + + } + + Preslice perMcCollision = aod::mcparticle::mcCollisionId; + void processMCrec(soa::Join::iterator const& collision, + soa::Join const&, + DauTracks const&, + aod::McParticles const& mcParticles) + { + if (!acceptEvent(collision, 1)) { + return; + } + + if (!collision.has_mcCollision()) { + registry.fill(HIST("fakeEvents"), 0); // no assoc. MC collisions + return; + } + + const auto& mcCollision = collision.mcCollision_as>(); + + auto tracksGroupedPVcontr = pvContribTracksIUEta1->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + int nTracksPVcontr = tracksGroupedPVcontr.size(); + + auto tracksGroupedGlobal = globalTracksIUEta05->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + int nTracksGlobal = tracksGroupedGlobal.size(); + + // N charged in FT0M region in corresponding gen. MC collision + auto mcPartSlice = mcParticles.sliceBy(perMcCollision, collision.mcCollision_as>().globalIndex()); + uint16_t nchFT0 = getGenNchInFT0Mregion(mcPartSlice); + + int evType = 0; + registry.fill(HIST("hNEvents"), 10.5); // reco INEL + if (collision.isInelGt0()) { + evType += 1; + registry.fill(HIST("hNEvents"), 11.5); // reco INEL>0 + } + + registry.fill(HIST("hCentFT0M_rec"), mcCollision.centFT0M(), evType); // correction only reco level in this stage + registry.fill(HIST("hCentFT0M_corr"), mcCollision.centFT0M(), mcCollision.centFT0M(), evType); + + if (multQA) { + registry.fill(HIST("hNchFT0MPVContr"), nchFT0, nTracksPVcontr, evType); + registry.fill(HIST("hFT0MpvContr"), mcCollision.centFT0M(), nTracksPVcontr, evType); + registry.fill(HIST("hFT0Mglobal"), mcCollision.centFT0M(), nTracksGlobal, evType); + registry.fill(HIST("hNchFT0Mglobal"), nchFT0, nTracksGlobal, evType); + registry.fill(HIST("hFT0MsignalPVContr"), collision.multFT0A() + collision.multFT0C(), nTracksPVcontr, evType); + } + + } + + void processMCgen(soa::Join::iterator const& mcCollision, + aod::McParticles const& mcParticles, + const soa::SmallGroups>& collisions) + { + auto cent = mcCollision.centFT0M(); + + registry.fill(HIST("hNEventsMC"), 0.5); + + if (isZvtxcut && std::fabs(mcCollision.posZ()) > cutzvertex) { + return; + } + registry.fill(HIST("hZCollisionGen"), mcCollision.posZ()); + registry.fill(HIST("hNEventsMC"), 1.5); + + int evType = 0; + registry.fill(HIST("hNEventsMC"), 2.5); + if (pwglf::isINELgtNmc(mcParticles, 0, pdgDB)) { // Truth INEL>0 + evType++; + registry.fill(HIST("hNEventsMC"), 3.5); + } + + fillMCParticles(mcParticles,cent, evType); + + registry.fill(HIST("hCentFT0M_genMC"), cent, evType); + + uint16_t nchFT0 = getGenNchInFT0Mregion(mcParticles); + registry.fill(HIST("hNchFT0MGenEvType"), nchFT0, evType); + + std::vector selectedEvents(collisions.size()); + std::vector numberOfContributors; + int nevts = 0; + int nAssocColl = 0; + const int nContSize = 2; + for (const auto& collision : collisions) { + CollisionIndexAndType collWithType = {0, 0x0}; + if (!acceptEvent(collision, 0)) { + continue; + } + collWithType.index = collision.mcCollision_as>().globalIndex(); + collWithType.typeFlag |= o2::aod::myMCcascades::EvFlags::EvINEL; + + if (collision.isInelGt0()) { // reco INEL>0 + collWithType.typeFlag |= o2::aod::myMCcascades::EvFlags::EvINELgt0; + } + selectedEvents[nevts++] = collWithType; + if (collision.mcCollision_as>().globalIndex() == mcCollision.globalIndex()) { + nAssocColl++; + numberOfContributors.push_back(collision.numContrib()); + } + } + selectedEvents.resize(nevts); + + registry.fill(HIST("hCentFT0MNAssocMCCollisions"),cent, nAssocColl, evType); + registry.fill(HIST("hNchFT0MNAssocMCCollisions"), nchFT0, nAssocColl, evType); + + if (numberOfContributors.size() == nContSize) { + std::sort(numberOfContributors.begin(), numberOfContributors.end()); + registry.fill(HIST("hNContributorsCorrelation"), numberOfContributors[0], numberOfContributors[1]); + } + + auto isAssocToINEL = [&mcCollision](CollisionIndexAndType i) { return (i.index == mcCollision.globalIndex()) && ((i.typeFlag & o2::aod::myMCcascades::EvFlags::EvINEL) == o2::aod::myMCcascades::EvFlags::EvINEL); }; + auto isAssocToINELgt0 = [&mcCollision](CollisionIndexAndType i) { return (i.index == mcCollision.globalIndex()) && ((i.typeFlag & o2::aod::myMCcascades::EvFlags::EvINELgt0) == o2::aod::myMCcascades::EvFlags::EvINELgt0); }; + // number of reconstructed INEL events that have the same global index as mcCollision + const auto evtReconstructedAndINEL = std::count_if(selectedEvents.begin(), selectedEvents.end(), isAssocToINEL); + // number of reconstructed INEL > 0 events that have the same global index as mcCollision + const auto evtReconstructedAndINELgt0 = std::count_if(selectedEvents.begin(), selectedEvents.end(), isAssocToINELgt0); + switch (evType) { + case 0: { + registry.fill(HIST("hCentFT0MNAssocMCCollisionsSameType"), cent, evtReconstructedAndINEL, evType); + registry.fill(HIST("hNchFT0MNAssocMCCollisionsSameType"), nchFT0, evtReconstructedAndINEL, evType); + break; + } + case 1: { + registry.fill(HIST("hCentFT0MNAssocMCCollisionsSameType"), cent, evtReconstructedAndINELgt0, evType); + registry.fill(HIST("hNchFT0MNAssocMCCollisionsSameType"), nchFT0, evtReconstructedAndINELgt0, evType); + break; + } + default: + LOGF(fatal, "incorrect evType in event task"); + break; + } + + if (evtReconstructedAndINELgt0) { // N INEL>0 reconstructed events associated with the MC collision + registry.fill(HIST("hNEventsMC"), 4.5); + } + } + PROCESS_SWITCH(Initializereventqa, processData, "Process Run 3 data", false); + PROCESS_SWITCH(Initializereventqa, processMCrec, "Process Run 3 mc, Reconstructed", true); + PROCESS_SWITCH(Initializereventqa, processMCgen, "Process Run 3 mc, genereated", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + }; +} diff --git a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx index 42773fa2947..ebfb994b255 100644 --- a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx +++ b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx @@ -29,6 +29,9 @@ #include "CommonConstants/PhysicsConstants.h" #include "Common/Core/RecoDecay.h" #include "Framework/O2DatabasePDGPlugin.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/inelGt.h" + using namespace o2; using namespace o2::framework; @@ -36,7 +39,7 @@ using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::constants::physics; using LorentzVectorPtEtaPhiMass = ROOT::Math::PtEtaPhiMVector; -// Service pdgDB; +Service pdgDB; enum { kData = 0, @@ -57,6 +60,10 @@ struct Xi1530Analysisqa { SliceCache cache; Preslice perRCol = aod::resodaughter::resoCollisionId; Preslice perCollision = aod::track::collisionId; + Preslice perResoCollision = + aod::resodaughter::resoCollisionId; + Preslice perResoCollisionCasc = + aod::resodaughter::resoCollisionId; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; using ResoMCCols = soa::Join; @@ -65,7 +72,7 @@ struct Xi1530Analysisqa { // Associated with histograms ConfigurableAxis binsPt{"binsPt", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10.0, 10.1, 10.2, 10.3, 10.4, 10.5, 10.6, 10.7, 10.8, 10.9, 11.0, 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7, 11.8, 11.9, 12.0, 12.1, 12.2, 12.3, 12.4, 12.5, 12.6, 12.7, 12.8, 12.9, 13.0, 13.1, 13.2, 13.3, 13.4, 13.5, 13.6, 13.7, 13.8, 13.9, 14.0, 14.1, 14.2, 14.3, 14.4, 14.5, 14.6, 14.7, 14.8, 14.9, 15.0}, "Binning of the pT axis"}; ConfigurableAxis binsPtQA{"binsPtQA", {VARIABLE_WIDTH, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2, 7.4, 7.6, 7.8, 8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.8, 10.0}, "Binning of the pT axis"}; - ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0}, "Binning of the centrality axis"}; + ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0, 110.0}, "Binning of the centrality axis"}; Configurable cInvMassStart{"cInvMassStart", 1.4, "Invariant mass start"}; Configurable cInvMassEnd{"cInvMassEnd", 3.0, "Invariant mass end"}; @@ -98,8 +105,8 @@ struct Xi1530Analysisqa { Configurable cfgTPCRows{"cfgTPCRows", 80, "Minimum Number of TPC Crossed Rows "}; Configurable cfgRatioTPCRowsOverFindableCls{"cfgRatioTPCRowsOverFindableCls", 0.8, "Minimum of TPC Crossed Rows to Findable Clusters"}; // Minmimum - Configurable cfgUseTPCRefit{"cfgUseTPCRefit", true, "Require TPC Refit"}; - Configurable cfgUseITSRefit{"cfgUseITSRefit", true, "Require ITS Refit"}; + // Configurable cfgUseTPCRefit{"cfgUseTPCRefit", true, "Require TPC Refit"}; //refit is included in global track selection + // Configurable cfgUseITSRefit{"cfgUseITSRefit", true, "Require ITS Refit"}; Configurable cfgHasTOF{"cfgHasTOF", false, "Require TOF"}; @@ -200,7 +207,6 @@ struct Xi1530Analysisqa { // MC Event selection // Configurable cZvertCutMC{"cZvertCutMC", 10.0, "MC Z-vertex cut"}; - Configurable cIsPhysicalPrimaryMC{"cIsPhysicalPrimaryMC", true, "Physical primary selection for a MC Parent"}; //*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*// @@ -209,8 +215,10 @@ struct Xi1530Analysisqa { Configurable cMaxPtMotherCut{"cMaxPtMotherCut", 9.0, "Maximum pt of mother cut"}; Configurable cMaxMinvMotherCut{"cMaxMinvMotherCut", 3.0, "Maximum Minv of mother cut"}; - Configurable cMicroTrack{"cMicroTrack", false, "Using Micro track for first pion"}; - Configurable studyStableXi{"studyStableXi", true, "Study stable Xi"}; + Configurable studyStableXi{"studyStableXi", false, "Study stable Xi"}; + + Configurable cMCCent{"cMCCent", true, "Using calibrated MC centrality (for FT0M)"}; + Configurable cRecoINELgt0{"cRecoINELgt0", true, "check if INEL>0 for reco events"}; TRandom* rn = new TRandom(); //*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*//*// @@ -251,11 +259,11 @@ struct Xi1530Analysisqa { // event histograms histos.add("QAevent/hEvtCounterSameE", "Number of analyzed Same Events", HistType::kTH1F, {{1, 0.5, 1.5}}); histos.add("QAevent/hVertexZSameE", "Collision Vertex Z position", HistType::kTH1F, {{100, -15., 15.}}); - histos.add("QAevent/hMultiplicityPercentSameE", "Multiplicity percentile of collision", HistType::kTH1F, {{120, 0.0f, 120.0f}}); + histos.add("QAevent/hMultiplicityPercentSameE", "Multiplicity percentile of collision", HistType::kTH1F, {centAxis}); histos.add("QAevent/hEvtCounterMixedE", "Number of analyzed Mixed Events", HistType::kTH1F, {{1, 0.5, 1.5}}); histos.add("QAevent/hVertexZMixedE", "Collision Vertex Z position", HistType::kTH1F, {{100, -15., 15.}}); - histos.add("QAevent/hMultiplicityPercentMixedE", "Multiplicity percentile of collision", HistType::kTH1F, {{120, 0.0f, 120.0f}}); + histos.add("QAevent/hMultiplicityPercentMixedE", "Multiplicity percentile of collision", HistType::kTH1F, {centAxis}); } if (invMass1D) { @@ -456,10 +464,10 @@ struct Xi1530Analysisqa { } if (cfgHasTOF && !track.hasTOF()) return false; - if (cfgUseITSRefit && !track.passedITSRefit()) - return false; - if (cfgUseTPCRefit && !track.passedTPCRefit()) - return false; + // if (cfgUseITSRefit && !track.passedITSRefit()) + // return false; + // if (cfgUseTPCRefit && !track.passedTPCRefit()) + // return false; if (cfgPVContributor && !track.isPVContributor()) return false; if (cfgPrimaryTrack && !track.isPrimaryTrack()) @@ -467,7 +475,7 @@ struct Xi1530Analysisqa { return true; } - bool hasSubsystemInfo(float Nsigma) // this will be replaced // .hasXX() was not appied in resocascade yet + bool hasSubsystemInfo(float Nsigma) { return std::abs(Nsigma) < cPIDBound; } @@ -492,7 +500,7 @@ struct Xi1530Analysisqa { return true; } - // Secondary track selection for cascades // need to more information, + // Secondary track selection for cascades // // Topological cuts for cascades template @@ -663,20 +671,20 @@ struct Xi1530Analysisqa { return lConsistentWithXi && lConsistentWithLambda; } - template - void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksTypeCasc& dTracks2) // Order: ResoColl, ResoTrack, ResoCascTrack + template + void fillHistograms(const CollisionType& collision, const CenMult& multiplicity, const TracksType& dTracks1, const TracksTypeCasc& dTracks2) // Order: ResoColl, ResoTrack, ResoCascTrack { - auto multiplicity = collision.cent(); + //auto multiplicity = collision.cent(); { if constexpr (!IsMix) { histos.fill(HIST("QAevent/hVertexZSameE"), collision.posZ()); - histos.fill(HIST("QAevent/hMultiplicityPercentSameE"), collision.cent()); + histos.fill(HIST("QAevent/hMultiplicityPercentSameE"), multiplicity); histos.fill(HIST("TestME/hCollisionIndexSameE"), collision.globalIndex()); histos.fill(HIST("TestME/hnTrksSameE"), dTracks1.size()); } else { histos.fill(HIST("QAevent/hVertexZMixedE"), collision.posZ()); - histos.fill(HIST("QAevent/hMultiplicityPercentMixedE"), collision.cent()); + histos.fill(HIST("QAevent/hMultiplicityPercentMixedE"), multiplicity); histos.fill(HIST("TestME/hCollisionIndexMixedE"), collision.globalIndex()); histos.fill(HIST("TestME/hnTrksMixedE"), dTracks1.size()); } @@ -742,7 +750,7 @@ struct Xi1530Analysisqa { float trk2NSigmaPiNegTOF = trk2.daughterTOFNSigmaNegPi(); if constexpr (!IsMix) { - //// QA plots before the selection // need to pt for cascade tracks + //// QA plots before the selection // // --- PID QA if (pidPlots) { histos.fill(HIST("QAbefore/TPC_Nsigma_pi_first_all"), multiplicity, trk1ptPi, trk1NSigmaPiTPC); @@ -750,7 +758,6 @@ struct Xi1530Analysisqa { histos.fill(HIST("QAbefore/TOF_Nsigma_pi_first_all"), multiplicity, trk1ptPi, trk1NSigmaPiTOF); histos.fill(HIST("QAbefore/TOF_TPC_Map_pi_first_all"), trk1NSigmaPiTOF, trk1NSigmaPiTPC); } - // hasSubsystemInfo is Temporary, it will be replaced. histos.fill(HIST("QAbefore/TPC_Nsigma_pi_bachelor_all"), multiplicity, 0, trk2NSigmaPiBachelorTPC); // can't take pt information for the cascade secondary if (hasSubsystemInfo(trk2NSigmaPiBachelorTOF)) { @@ -819,10 +826,6 @@ struct Xi1530Analysisqa { if (!casctopCut(trk2)) continue; - // TPCncluster distributions - // histos.fill(HIST("TPCncluster/TPCnclusterpifirst"), trk1.tpcNClsFound()); - // histos.fill(HIST("TPCncluster/TPCnclusterPhipifirst"), trk1.tpcNClsFound(), trk1.phi()); - if constexpr (!IsMix) { //// QA plots after the selection // --- PID QA @@ -1048,25 +1051,75 @@ struct Xi1530Analysisqa { } } - void processData(aod::ResoCollision const& resoCollision, aod::ResoTracks const& resoTracks, aod::ResoCascades const& cascTracks) + void processData(aod::ResoCollision const& resoCollision, + aod::ResoCollisionColls const& collisionIndex, + soa::Join const& collisions, + aod::ResoTracks const& resoTracks, + aod::ResoCascades const& cascTracks) { + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision + + auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision + + if(cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; histos.fill(HIST("QAevent/hEvtCounterSameE"), 1.0); - fillHistograms(resoCollision, resoTracks, cascTracks); + auto multiplicity = resoCollision.cent(); + fillHistograms(resoCollision,multiplicity,resoTracks, cascTracks); } + // Reconstructed level MC for the track void processMC(ResoMCCols::iterator const& resoCollision, + aod::ResoCollisionColls const& resoCollisionIndex, + soa::Join const& collisionsMC, soa::Join const& cascTracks, - soa::Join const& resoTracks) + soa::Join const& resoTracks, + soa::Join const&) { if (!resoCollision.isInAfterAllCuts() || (std::abs(resoCollision.posZ()) > cZvertCutMC)) // MC event selection, all cuts missing vtx cut return; - fillHistograms(resoCollision, resoTracks, cascTracks); + + + auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); + const int collId = linkRow.collisionId(); + + auto coll = collisionsMC.iteratorAt(collId); + + if(cRecoINELgt0 && !coll.isInelGt0()) + return; + + auto mcColl = coll.mcCollision_as>(); + + auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. + auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. + + float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; + + fillHistograms(resoCollision, multiplicity, resoTracks, cascTracks); } - void processMCTrue(ResoMCCols::iterator const& resoCollision, aod::ResoMCParents const& resoParents) + // Truth level MC for the track with reco event + void processMCTrue(ResoMCCols::iterator const& resoCollision, + aod::ResoCollisionColls const& resoCollisionIndex, + aod::ResoMCParents const& resoParents, + aod::ResoCollisionCandidatesMC const& collisionsMC, + soa::Join const&) { - auto multiplicity = resoCollision.cent(); + + auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); + const int collId = linkRow.collisionId(); + + auto coll = collisionsMC.iteratorAt(collId); + + auto mcColl = coll.mcCollision_as>(); + + auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. + auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. + + float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; + for (const auto& part : resoParents) { // loop over all pre-filtered MC particles if (std::abs(part.pdgCode()) != kXiStar || std::abs(part.y()) >= cfgRapidityCut) continue; @@ -1076,9 +1129,6 @@ struct Xi1530Analysisqa { if (!pass1 || !pass2) continue; - if (cIsPhysicalPrimaryMC && !part.isPhysicalPrimary()) - continue; - if (part.pdgCode() > 0) // INELt0 or INEL histos.fill(HIST("h3Xi1530Gen"), -1, part.pt(), multiplicity); else @@ -1115,11 +1165,23 @@ struct Xi1530Analysisqa { } } - void processDataMicro(aod::ResoCollision const& resoCollision, aod::ResoMicroTracks const& resomicrotracks, aod::ResoCascades const& cascTracks) + void processDataMicro(aod::ResoCollision const& resoCollision, + aod::ResoCollisionColls const& collisionIndex, + soa::Join const& collisions, + aod::ResoMicroTracks const& resomicrotracks, + aod::ResoCascades const& cascTracks) { + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision + + auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision + + if(cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; histos.fill(HIST("QAevent/hEvtCounterSameE"), 1.0); - fillHistograms(resoCollision, resomicrotracks, cascTracks); + auto multiplicity = resoCollision.cent(); + fillHistograms(resoCollision, multiplicity, resomicrotracks, cascTracks); } using BinningTypeVtxZT0M = ColumnBinningPolicy; @@ -1137,16 +1199,21 @@ struct Xi1530Analysisqa { for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { histos.fill(HIST("QAevent/hEvtCounterMixedE"), 1.0); - fillHistograms(collision1, tracks1, tracks2); + auto multiplicity = collision1.cent(); + fillHistograms(collision1,multiplicity, tracks1, tracks2); } } void processDataDF(aod::ResoCollisionDF const& resoCollision, aod::ResoTrackDFs const& resotracks, aod::ResoCascadeDFs const& cascTracks) { - - fillHistograms(resoCollision, resotracks, cascTracks); + auto multiplicity = resoCollision.cent(); + fillHistograms(resoCollision, multiplicity, resotracks, cascTracks); } - void processMEMicro(aod::ResoCollisions const& resoCollisions, aod::ResoMicroTracks const& resomicrotracks, aod::ResoCascades const& cascTracks) + void processMEMicro(aod::ResoCollisions const& resoCollisions, + aod::ResoCollisionColls const& collisionIndex, + soa::Join const& collisions, + aod::ResoMicroTracks const& resomicrotracks, + aod::ResoCascades const& cascTracks) { auto tracksTuple = std::make_tuple(resomicrotracks, cascTracks); @@ -1155,8 +1222,18 @@ struct Xi1530Analysisqa { for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { + const auto rcIdx = collision1.globalIndex(); + + const auto linkRow = collisionIndex.iteratorAt(rcIdx); + const auto collId = linkRow.collisionId(); + + auto coll = collisions.iteratorAt(collId); + if(cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + continue; + histos.fill(HIST("QAevent/hEvtCounterMixedE"), 1.0); - fillHistograms(collision1, tracks1, tracks2); + auto multiplicity = collision1.cent(); + fillHistograms(collision1, multiplicity,tracks1, tracks2); } } From a4821a95cf2e1c31a0a743fe3a5f82becb45ea40 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 20 Aug 2025 14:17:25 +0000 Subject: [PATCH 2/6] Please consider the following formatting changes --- PWGLF/Tasks/Resonances/initializereventqa.cxx | 90 ++++++------ PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx | 130 +++++++++--------- 2 files changed, 108 insertions(+), 112 deletions(-) diff --git a/PWGLF/Tasks/Resonances/initializereventqa.cxx b/PWGLF/Tasks/Resonances/initializereventqa.cxx index d3ddae10cbb..076163e6ef0 100644 --- a/PWGLF/Tasks/Resonances/initializereventqa.cxx +++ b/PWGLF/Tasks/Resonances/initializereventqa.cxx @@ -13,36 +13,38 @@ /// \brief QA for the event-loss and signal-loss correction at the generator level for the ResonanceInitializer in pp collisions (referred to TableProducer/Strangeness/cascqaanalysis.cxx) /// /// Following the discussions at the two PAG meetings (https://indico.cern.ch/event/1518979, https://indico.cern.ch/event/1575984) -/// we have introduced an auxiliary task that, when the resonanceInitializer.cxx is used, +/// we have introduced an auxiliary task that, when the resonanceInitializer.cxx is used, /// computes the event-loss and signal-loss correction factors at the generator level. -/// With minor configuration tuning for a truth-tagging, +/// With minor configuration tuning for a truth-tagging, /// we expect it to be applicable to most analyses that rely on the initializer. /// /// \author Minjae Kim (minjae.kim@cern.ch) -#include -#include -#include - -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Common/DataModel/TrackSelectionTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/cascqaanalysis.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/inelGt.h" + +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/Centrality.h" -#include "PWGLF/DataModel/cascqaanalysis.h" -#include "TRandom2.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "Framework/AnalysisTask.h" #include "Framework/O2DatabasePDGPlugin.h" -#include "PWGLF/Utils/inelGt.h" -#include "PWGLF/DataModel/mcCentrality.h" +#include "Framework/runDataProcessing.h" + +#include "TRandom2.h" +#include + +#include +#include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; - using TrkPidInfo = soa::Join; using DauTracks = soa::Join; @@ -58,7 +60,7 @@ struct Initializereventqa { ConfigurableAxis ptAxis{"ptAxis", {400, 0.0f, 20.0f}, "#it{p}_{T} (GeV/#it{c})"}; ConfigurableAxis rapidityAxis{"rapidityAxis", {200, -2.0f, 2.0f}, "y"}; ConfigurableAxis centFT0MAxis{"centFT0MAxis", - {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0, 110.0}, + {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0, 110.0}, "FT0M (%)"}; ConfigurableAxis eventTypeAxis{"eventTypeAxis", {2, -0.5f, 1.5f}, "Event Type"}; @@ -117,7 +119,6 @@ struct Initializereventqa { } registry.add("hZCollision", "hZCollision", {HistType::kTH1D, {{200, -20.f, 20.f}}}); - registry.add("fakeEvents", "fakeEvents", {HistType::kTH1F, {{1, -0.5f, 0.5f}}}); registry.add("hNEventsMC", "hNEventsMC", {HistType::kTH1D, {{5, 0.0f, 5.0f}}}); @@ -125,18 +126,17 @@ struct Initializereventqa { registry.get(HIST("hNEventsMC"))->GetXaxis()->SetBinLabel(n, hNEventsMCLabels[n - 1]); } registry.add("hZCollisionGen", "hZCollisionGen", {HistType::kTH1D, {{200, -20.f, 20.f}}}); - registry.add("hCentFT0MNAssocMCCollisions", "hCentFT0MNAssocMCCollisions", {HistType::kTH3D, {centFT0MAxis, nAssocCollAxis, eventTypeAxis}}); + registry.add("hCentFT0MNAssocMCCollisions", "hCentFT0MNAssocMCCollisions", {HistType::kTH3D, {centFT0MAxis, nAssocCollAxis, eventTypeAxis}}); registry.add("hCentFT0MNAssocMCCollisionsSameType", "hCentFT0MNAssocMCCollisionsSameType", {HistType::kTH3D, {centFT0MAxis, nAssocCollAxis, eventTypeAxis}}); registry.add("hNchFT0MNAssocMCCollisions", "hNchFT0MNAssocMCCollisions", {HistType::kTH3D, {nChargedFT0MGenAxis, nAssocCollAxis, eventTypeAxis}}); registry.add("hNchFT0MNAssocMCCollisionsSameType", "hNchFT0MNAssocMCCollisionsSameType", {HistType::kTH3D, {nChargedFT0MGenAxis, nAssocCollAxis, eventTypeAxis}}); registry.add("hNContributorsCorrelation", "hNContributorsCorrelation", {HistType::kTH2F, {{250, -0.5f, 249.5f, "Secondary Contributor"}, {250, -0.5f, 249.5f, "Main Contributor"}}}); registry.add("hNchFT0MGenEvType", "hNchFT0MGenEvType", {HistType::kTH2D, {nChargedFT0MGenAxis, eventTypeAxis}}); registry.add("hCentFT0M_genMC", "hCentFT0M_genMC", {HistType::kTH2D, {centFT0MAxis, eventTypeAxis}}); - + registry.add("hCentFT0M_rec", "hCentFT0M_rec", {HistType::kTH2D, {centFT0MAxis, eventTypeAxis}}); registry.add("hCentFT0M_corr", "hCentFT0M_Corr", {HistType::kTH2D, {centFT0MAxis, centFT0MAxis}}); - if (multQA) { registry.add("hNchFT0Mglobal", "hNchFT0Mglobal", {HistType::kTH3D, {nChargedFT0MGenAxis, multNTracksAxis, eventTypeAxis}}); registry.add("hNchFT0MPVContr", "hNchFT0MPVContr", {HistType::kTH3D, {nChargedFT0MGenAxis, multNTracksAxis, eventTypeAxis}}); @@ -145,23 +145,22 @@ struct Initializereventqa { registry.add("hFT0MsignalPVContr", "hFT0MsignalPVContr", {HistType::kTH3D, {signalFT0MAxis, multNTracksAxis, eventTypeAxis}}); } - registry.add("h3ResonanceTruth", "pT distribution of True Resonance", kTHnSparseF, {eventTypeAxis, ptAxis, centFT0MAxis}); - registry.add("h3ResonanceTruthAnti", "pT distribution of True Resonance Anti", kTHnSparseF, {eventTypeAxis, ptAxis, centFT0MAxis}); + registry.add("h3ResonanceTruth", "pT distribution of True Resonance", kTHnSparseF, {eventTypeAxis, ptAxis, centFT0MAxis}); + registry.add("h3ResonanceTruthAnti", "pT distribution of True Resonance Anti", kTHnSparseF, {eventTypeAxis, ptAxis, centFT0MAxis}); } -float pvEta1 = 1.0f; + float pvEta1 = 1.0f; float globalEta05 = 0.5f; Partition pvContribTracksIUEta1 = (nabs(aod::track::eta) < pvEta1) && ((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)); Partition globalTracksIUEta05 = (nabs(aod::track::eta) < globalEta05) && (requireGlobalTrackInFilter()); - template uint16_t getGenNchInFT0Mregion(TMcParticles particles) { float region1FT0 = -3.3f; float region2FT0 = -2.1f; float region3FT0 = 3.5f; - float region4FT0 = 4.9f; + float region4FT0 = 4.9f; // Particle counting in FITFT0: -3.3<η<-2.1; 3.5<η<4.9 uint16_t nchFT0 = 0; for (const auto& mcParticle : particles) { @@ -264,7 +263,6 @@ float pvEta1 = 1.0f; return true; } - template void fillMCParticles(TotalMCParts const& mcParticles, MultMCGen const& multiplicity, evtType const& eventType) @@ -282,24 +280,24 @@ float pvEta1 = 1.0f; daughterPDGs = {-1, -1}; } - if(isDaughterCheck){ - bool pass1 = std::abs(daughterPDGs[0]) == pdgTruthDaughter1 || std::abs(daughterPDGs[1]) == pdgTruthDaughter1; - bool pass2 = std::abs(daughterPDGs[0]) == pdgTruthDaughter2 || std::abs(daughterPDGs[1]) == pdgTruthDaughter2; - if (!pass1 || !pass2) - continue; + if (isDaughterCheck) { + bool pass1 = std::abs(daughterPDGs[0]) == pdgTruthDaughter1 || std::abs(daughterPDGs[1]) == pdgTruthDaughter1; + bool pass2 = std::abs(daughterPDGs[0]) == pdgTruthDaughter2 || std::abs(daughterPDGs[1]) == pdgTruthDaughter2; + if (!pass1 || !pass2) + continue; } if (mcPart.pdgCode() > 0) // Consider INELt0 or INEL registry.fill(HIST("h3ResonanceTruth"), eventType, mcPart.pt(), multiplicity); else registry.fill(HIST("h3ResonanceTruthAnti"), eventType, mcPart.pt(), multiplicity); - daughterPDGs.clear(); + daughterPDGs.clear(); } } - void processData(soa::Join::iterator const& collision, - DauTracks const&) + DauTracks const&) { if (!acceptEvent(collision, 1)) { return; @@ -325,7 +323,6 @@ float pvEta1 = 1.0f; registry.fill(HIST("hFT0Mglobal"), collision.centFT0M(), nTracksGlobal, evType); registry.fill(HIST("hFT0MsignalPVContr"), collision.multFT0A() + collision.multFT0C(), nTracksPVcontr, evType); } - } Preslice perMcCollision = aod::mcparticle::mcCollisionId; @@ -354,8 +351,8 @@ float pvEta1 = 1.0f; int nTracksGlobal = tracksGroupedGlobal.size(); // N charged in FT0M region in corresponding gen. MC collision - auto mcPartSlice = mcParticles.sliceBy(perMcCollision, collision.mcCollision_as>().globalIndex()); - uint16_t nchFT0 = getGenNchInFT0Mregion(mcPartSlice); + auto mcPartSlice = mcParticles.sliceBy(perMcCollision, collision.mcCollision_as>().globalIndex()); + uint16_t nchFT0 = getGenNchInFT0Mregion(mcPartSlice); int evType = 0; registry.fill(HIST("hNEvents"), 10.5); // reco INEL @@ -374,19 +371,18 @@ float pvEta1 = 1.0f; registry.fill(HIST("hNchFT0Mglobal"), nchFT0, nTracksGlobal, evType); registry.fill(HIST("hFT0MsignalPVContr"), collision.multFT0A() + collision.multFT0C(), nTracksPVcontr, evType); } - } void processMCgen(soa::Join::iterator const& mcCollision, aod::McParticles const& mcParticles, - const soa::SmallGroups>& collisions) + const soa::SmallGroups>& collisions) { auto cent = mcCollision.centFT0M(); registry.fill(HIST("hNEventsMC"), 0.5); - if (isZvtxcut && std::fabs(mcCollision.posZ()) > cutzvertex) { + if (isZvtxcut && std::fabs(mcCollision.posZ()) > cutzvertex) { return; } registry.fill(HIST("hZCollisionGen"), mcCollision.posZ()); @@ -399,7 +395,7 @@ float pvEta1 = 1.0f; registry.fill(HIST("hNEventsMC"), 3.5); } - fillMCParticles(mcParticles,cent, evType); + fillMCParticles(mcParticles, cent, evType); registry.fill(HIST("hCentFT0M_genMC"), cent, evType); @@ -423,16 +419,16 @@ float pvEta1 = 1.0f; collWithType.typeFlag |= o2::aod::myMCcascades::EvFlags::EvINELgt0; } selectedEvents[nevts++] = collWithType; - if (collision.mcCollision_as>().globalIndex() == mcCollision.globalIndex()) { + if (collision.mcCollision_as>().globalIndex() == mcCollision.globalIndex()) { nAssocColl++; numberOfContributors.push_back(collision.numContrib()); } } selectedEvents.resize(nevts); - registry.fill(HIST("hCentFT0MNAssocMCCollisions"),cent, nAssocColl, evType); + registry.fill(HIST("hCentFT0MNAssocMCCollisions"), cent, nAssocColl, evType); registry.fill(HIST("hNchFT0MNAssocMCCollisions"), nchFT0, nAssocColl, evType); - + if (numberOfContributors.size() == nContSize) { std::sort(numberOfContributors.begin(), numberOfContributors.end()); registry.fill(HIST("hNContributorsCorrelation"), numberOfContributors[0], numberOfContributors[1]); @@ -473,5 +469,5 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ adaptAnalysisTask(cfgc), - }; + }; } diff --git a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx index ebfb994b255..07ed10a9460 100644 --- a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx +++ b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx @@ -14,24 +14,25 @@ /// /// \author Min-jae Kim , Bong-Hwi Lim // #include -#include "Math/Vector4D.h" -#include "TF1.h" -#include "TRandom3.h" +#include "PWGLF/DataModel/LFResonanceTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/inelGt.h" -#include "Common/DataModel/PIDResponse.h" +#include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" -#include "Framework/AnalysisTask.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/runDataProcessing.h" -#include "PWGLF/DataModel/LFResonanceTables.h" -#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "Common/DataModel/PIDResponse.h" + #include "CommonConstants/PhysicsConstants.h" -#include "Common/Core/RecoDecay.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" #include "Framework/O2DatabasePDGPlugin.h" -#include "PWGLF/DataModel/mcCentrality.h" -#include "PWGLF/Utils/inelGt.h" +#include "Framework/runDataProcessing.h" +#include "Math/Vector4D.h" +#include "TF1.h" +#include "TRandom3.h" using namespace o2; using namespace o2::framework; @@ -475,7 +476,7 @@ struct Xi1530Analysisqa { return true; } - bool hasSubsystemInfo(float Nsigma) + bool hasSubsystemInfo(float Nsigma) { return std::abs(Nsigma) < cPIDBound; } @@ -674,7 +675,7 @@ struct Xi1530Analysisqa { template void fillHistograms(const CollisionType& collision, const CenMult& multiplicity, const TracksType& dTracks1, const TracksTypeCasc& dTracks2) // Order: ResoColl, ResoTrack, ResoCascTrack { - //auto multiplicity = collision.cent(); + // auto multiplicity = collision.cent(); { if constexpr (!IsMix) { @@ -750,7 +751,7 @@ struct Xi1530Analysisqa { float trk2NSigmaPiNegTOF = trk2.daughterTOFNSigmaNegPi(); if constexpr (!IsMix) { - //// QA plots before the selection // + //// QA plots before the selection // // --- PID QA if (pidPlots) { histos.fill(HIST("QAbefore/TPC_Nsigma_pi_first_all"), multiplicity, trk1ptPi, trk1NSigmaPiTPC); @@ -1051,23 +1052,23 @@ struct Xi1530Analysisqa { } } - void processData(aod::ResoCollision const& resoCollision, - aod::ResoCollisionColls const& collisionIndex, - soa::Join const& collisions, - aod::ResoTracks const& resoTracks, - aod::ResoCascades const& cascTracks) + void processData(aod::ResoCollision const& resoCollision, + aod::ResoCollisionColls const& collisionIndex, + soa::Join const& collisions, + aod::ResoTracks const& resoTracks, + aod::ResoCascades const& cascTracks) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision - if(cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + if (cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision return; histos.fill(HIST("QAevent/hEvtCounterSameE"), 1.0); auto multiplicity = resoCollision.cent(); - fillHistograms(resoCollision,multiplicity,resoTracks, cascTracks); + fillHistograms(resoCollision, multiplicity, resoTracks, cascTracks); } // Reconstructed level MC for the track @@ -1080,46 +1081,45 @@ struct Xi1530Analysisqa { { if (!resoCollision.isInAfterAllCuts() || (std::abs(resoCollision.posZ()) > cZvertCutMC)) // MC event selection, all cuts missing vtx cut return; - - auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); - const int collId = linkRow.collisionId(); + auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); + const int collId = linkRow.collisionId(); - auto coll = collisionsMC.iteratorAt(collId); + auto coll = collisionsMC.iteratorAt(collId); - if(cRecoINELgt0 && !coll.isInelGt0()) - return; + if (cRecoINELgt0 && !coll.isInelGt0()) + return; - auto mcColl = coll.mcCollision_as>(); + auto mcColl = coll.mcCollision_as>(); - auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. - auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. + auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. + auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. - float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; + float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; fillHistograms(resoCollision, multiplicity, resoTracks, cascTracks); } // Truth level MC for the track with reco event - void processMCTrue(ResoMCCols::iterator const& resoCollision, - aod::ResoCollisionColls const& resoCollisionIndex, - aod::ResoMCParents const& resoParents, - aod::ResoCollisionCandidatesMC const& collisionsMC, - soa::Join const&) + void processMCTrue(ResoMCCols::iterator const& resoCollision, + aod::ResoCollisionColls const& resoCollisionIndex, + aod::ResoMCParents const& resoParents, + aod::ResoCollisionCandidatesMC const& collisionsMC, + soa::Join const&) { - - auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); - const int collId = linkRow.collisionId(); - auto coll = collisionsMC.iteratorAt(collId); + auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); + const int collId = linkRow.collisionId(); - auto mcColl = coll.mcCollision_as>(); + auto coll = collisionsMC.iteratorAt(collId); - auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. - auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. + auto mcColl = coll.mcCollision_as>(); + + auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. + auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. + + float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; - float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; - for (const auto& part : resoParents) { // loop over all pre-filtered MC particles if (std::abs(part.pdgCode()) != kXiStar || std::abs(part.y()) >= cfgRapidityCut) continue; @@ -1165,18 +1165,18 @@ struct Xi1530Analysisqa { } } - void processDataMicro(aod::ResoCollision const& resoCollision, - aod::ResoCollisionColls const& collisionIndex, - soa::Join const& collisions, - aod::ResoMicroTracks const& resomicrotracks, - aod::ResoCascades const& cascTracks) + void processDataMicro(aod::ResoCollision const& resoCollision, + aod::ResoCollisionColls const& collisionIndex, + soa::Join const& collisions, + aod::ResoMicroTracks const& resomicrotracks, + aod::ResoCascades const& cascTracks) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision - if(cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + if (cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision return; histos.fill(HIST("QAevent/hEvtCounterSameE"), 1.0); @@ -1200,7 +1200,7 @@ struct Xi1530Analysisqa { histos.fill(HIST("QAevent/hEvtCounterMixedE"), 1.0); auto multiplicity = collision1.cent(); - fillHistograms(collision1,multiplicity, tracks1, tracks2); + fillHistograms(collision1, multiplicity, tracks1, tracks2); } } void processDataDF(aod::ResoCollisionDF const& resoCollision, aod::ResoTrackDFs const& resotracks, aod::ResoCascadeDFs const& cascTracks) @@ -1209,11 +1209,11 @@ struct Xi1530Analysisqa { fillHistograms(resoCollision, multiplicity, resotracks, cascTracks); } - void processMEMicro(aod::ResoCollisions const& resoCollisions, - aod::ResoCollisionColls const& collisionIndex, - soa::Join const& collisions, - aod::ResoMicroTracks const& resomicrotracks, - aod::ResoCascades const& cascTracks) + void processMEMicro(aod::ResoCollisions const& resoCollisions, + aod::ResoCollisionColls const& collisionIndex, + soa::Join const& collisions, + aod::ResoMicroTracks const& resomicrotracks, + aod::ResoCascades const& cascTracks) { auto tracksTuple = std::make_tuple(resomicrotracks, cascTracks); @@ -1222,18 +1222,18 @@ struct Xi1530Analysisqa { for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { - const auto rcIdx = collision1.globalIndex(); + const auto rcIdx = collision1.globalIndex(); - const auto linkRow = collisionIndex.iteratorAt(rcIdx); - const auto collId = linkRow.collisionId(); + const auto linkRow = collisionIndex.iteratorAt(rcIdx); + const auto collId = linkRow.collisionId(); - auto coll = collisions.iteratorAt(collId); - if(cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - continue; + auto coll = collisions.iteratorAt(collId); + if (cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + continue; histos.fill(HIST("QAevent/hEvtCounterMixedE"), 1.0); auto multiplicity = collision1.cent(); - fillHistograms(collision1, multiplicity,tracks1, tracks2); + fillHistograms(collision1, multiplicity, tracks1, tracks2); } } From aeba5a805128e778e5d1785de9a384c6511f6539 Mon Sep 17 00:00:00 2001 From: mjkim525 Date: Thu, 21 Aug 2025 21:23:42 +0900 Subject: [PATCH 3/6] Fix xi1530Analysisqa.cxx --- PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx | 121 ++++++++++++++------ 1 file changed, 86 insertions(+), 35 deletions(-) diff --git a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx index 07ed10a9460..ecf631fb9e1 100644 --- a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx +++ b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx @@ -1058,13 +1058,15 @@ struct Xi1530Analysisqa { aod::ResoTracks const& resoTracks, aod::ResoCascades const& cascTracks) { - auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); - auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision + if(cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision + auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision - if (cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - return; + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; + } histos.fill(HIST("QAevent/hEvtCounterSameE"), 1.0); auto multiplicity = resoCollision.cent(); @@ -1073,52 +1075,98 @@ struct Xi1530Analysisqa { // Reconstructed level MC for the track void processMC(ResoMCCols::iterator const& resoCollision, - aod::ResoCollisionColls const& resoCollisionIndex, + aod::ResoCollisionColls const& collisionIndex, soa::Join const& collisionsMC, soa::Join const& cascTracks, soa::Join const& resoTracks, soa::Join const&) { - if (!resoCollision.isInAfterAllCuts() || (std::abs(resoCollision.posZ()) > cZvertCutMC)) // MC event selection, all cuts missing vtx cut - return; + float multiplicity; + if(cMCCent && cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision + + auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision + + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; + + auto mcColl = coll.mcCollision_as>(); + multiplicity = mcColl.centFT0M(); + } + else if(!cMCCent && cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); - const int collId = linkRow.collisionId(); + auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision - auto coll = collisionsMC.iteratorAt(collId); + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; - if (cRecoINELgt0 && !coll.isInelGt0()) - return; + multiplicity = resoCollision.cent(); + } + else if(cMCCent && !cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - auto mcColl = coll.mcCollision_as>(); + auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision - auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. - auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. + auto mcColl = coll.mcCollision_as>(); + multiplicity = mcColl.centFT0M(); + } + else { + multiplicity = resoCollision.cent(); + } - float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; + if (!resoCollision.isInAfterAllCuts() || (std::abs(resoCollision.posZ()) > cZvertCutMC)) // MC event selection, all cuts missing vtx cut + return; fillHistograms(resoCollision, multiplicity, resoTracks, cascTracks); } // Truth level MC for the track with reco event void processMCTrue(ResoMCCols::iterator const& resoCollision, - aod::ResoCollisionColls const& resoCollisionIndex, + aod::ResoCollisionColls const& collisionIndex, aod::ResoMCParents const& resoParents, aod::ResoCollisionCandidatesMC const& collisionsMC, soa::Join const&) { + float multiplicity; + if(cMCCent && cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - auto linkRow = resoCollisionIndex.iteratorAt(resoCollision.globalIndex()); - const int collId = linkRow.collisionId(); + auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision - auto coll = collisionsMC.iteratorAt(collId); + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; - auto mcColl = coll.mcCollision_as>(); + auto mcColl = coll.mcCollision_as>(); + multiplicity = mcColl.centFT0M(); + } + else if(!cMCCent && cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision + + auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision + + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; + + multiplicity = resoCollision.cent(); + } + else if(cMCCent && !cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - auto multiplicityReco = resoCollision.cent(); // Reco level multiplicity per. - auto multiplicityGen = mcColl.centFT0M(); // Gen level multiplicity per. + auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision - float multiplicity = cMCCent ? multiplicityGen : multiplicityReco; + auto mcColl = coll.mcCollision_as>(); + multiplicity = mcColl.centFT0M(); + } + else { + multiplicity = resoCollision.cent(); + } for (const auto& part : resoParents) { // loop over all pre-filtered MC particles if (std::abs(part.pdgCode()) != kXiStar || std::abs(part.y()) >= cfgRapidityCut) @@ -1171,13 +1219,15 @@ struct Xi1530Analysisqa { aod::ResoMicroTracks const& resomicrotracks, aod::ResoCascades const& cascTracks) { - auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); - auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision + if(cRecoINELgt0){ + auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision + auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision - if (cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - return; + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + return; + } histos.fill(HIST("QAevent/hEvtCounterSameE"), 1.0); auto multiplicity = resoCollision.cent(); @@ -1222,15 +1272,16 @@ struct Xi1530Analysisqa { for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { + if(cRecoINELgt0){ const auto rcIdx = collision1.globalIndex(); + auto linkRow = collisionIndex.iteratorAt(rcIdx); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - const auto linkRow = collisionIndex.iteratorAt(rcIdx); - const auto collId = linkRow.collisionId(); + auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision - auto coll = collisions.iteratorAt(collId); - if (cRecoINELgt0 && !coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision continue; - + } histos.fill(HIST("QAevent/hEvtCounterMixedE"), 1.0); auto multiplicity = collision1.cent(); fillHistograms(collision1, multiplicity, tracks1, tracks2); From f90f3e233e047d0bb5f608da82549f7b591938ea Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 21 Aug 2025 12:28:45 +0000 Subject: [PATCH 4/6] Please consider the following formatting changes --- PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx | 44 +++++++++------------ 1 file changed, 19 insertions(+), 25 deletions(-) diff --git a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx index ecf631fb9e1..d33b0e589cf 100644 --- a/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx +++ b/PWGLF/Tasks/Resonances/xi1530Analysisqa.cxx @@ -1058,7 +1058,7 @@ struct Xi1530Analysisqa { aod::ResoTracks const& resoTracks, aod::ResoCascades const& cascTracks) { - if(cRecoINELgt0){ + if (cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1082,7 +1082,7 @@ struct Xi1530Analysisqa { soa::Join const&) { float multiplicity; - if(cMCCent && cRecoINELgt0){ + if (cMCCent && cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1093,8 +1093,7 @@ struct Xi1530Analysisqa { auto mcColl = coll.mcCollision_as>(); multiplicity = mcColl.centFT0M(); - } - else if(!cMCCent && cRecoINELgt0){ + } else if (!cMCCent && cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1104,8 +1103,7 @@ struct Xi1530Analysisqa { return; multiplicity = resoCollision.cent(); - } - else if(cMCCent && !cRecoINELgt0){ + } else if (cMCCent && !cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1113,13 +1111,12 @@ struct Xi1530Analysisqa { auto mcColl = coll.mcCollision_as>(); multiplicity = mcColl.centFT0M(); - } - else { + } else { multiplicity = resoCollision.cent(); } if (!resoCollision.isInAfterAllCuts() || (std::abs(resoCollision.posZ()) > cZvertCutMC)) // MC event selection, all cuts missing vtx cut - return; + return; fillHistograms(resoCollision, multiplicity, resoTracks, cascTracks); } @@ -1132,7 +1129,7 @@ struct Xi1530Analysisqa { soa::Join const&) { float multiplicity; - if(cMCCent && cRecoINELgt0){ + if (cMCCent && cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1143,8 +1140,7 @@ struct Xi1530Analysisqa { auto mcColl = coll.mcCollision_as>(); multiplicity = mcColl.centFT0M(); - } - else if(!cMCCent && cRecoINELgt0){ + } else if (!cMCCent && cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1154,8 +1150,7 @@ struct Xi1530Analysisqa { return; multiplicity = resoCollision.cent(); - } - else if(cMCCent && !cRecoINELgt0){ + } else if (cMCCent && !cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1163,8 +1158,7 @@ struct Xi1530Analysisqa { auto mcColl = coll.mcCollision_as>(); multiplicity = mcColl.centFT0M(); - } - else { + } else { multiplicity = resoCollision.cent(); } @@ -1219,7 +1213,7 @@ struct Xi1530Analysisqa { aod::ResoMicroTracks const& resomicrotracks, aod::ResoCascades const& cascTracks) { - if(cRecoINELgt0){ + if (cRecoINELgt0) { auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex()); auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision @@ -1272,16 +1266,16 @@ struct Xi1530Analysisqa { for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { - if(cRecoINELgt0){ - const auto rcIdx = collision1.globalIndex(); - auto linkRow = collisionIndex.iteratorAt(rcIdx); - auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision + if (cRecoINELgt0) { + const auto rcIdx = collision1.globalIndex(); + auto linkRow = collisionIndex.iteratorAt(rcIdx); + auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision - auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision + auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision - if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - continue; - } + if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision + continue; + } histos.fill(HIST("QAevent/hEvtCounterMixedE"), 1.0); auto multiplicity = collision1.cent(); fillHistograms(collision1, multiplicity, tracks1, tracks2); From a03e4da236d74b6d8670ca9ec1d1032d663e1c51 Mon Sep 17 00:00:00 2001 From: mjkim525 Date: Tue, 3 Feb 2026 22:14:48 +0900 Subject: [PATCH 5/6] Reso-Framework Update: add gen-level info.,MCCent. calibration --- PWGLF/DataModel/LFResonanceTables.h | 16 +- .../Resonances/resonanceInitializer.cxx | 284 ++++++++++++++++-- .../Resonances/resonanceModuleInitializer.cxx | 14 +- 3 files changed, 284 insertions(+), 30 deletions(-) diff --git a/PWGLF/DataModel/LFResonanceTables.h b/PWGLF/DataModel/LFResonanceTables.h index fcbe3b6378a..060fac0810f 100644 --- a/PWGLF/DataModel/LFResonanceTables.h +++ b/PWGLF/DataModel/LFResonanceTables.h @@ -68,6 +68,7 @@ DECLARE_SOA_COLUMN(EvtPlResAB, evtPlResAB, float); // DECLARE_SOA_COLUMN(EvtPlResAC, evtPlResAC, float); //! Second harmonic event plane resolution of A-C sub events DECLARE_SOA_COLUMN(EvtPlResBC, evtPlResBC, float); //! Second harmonic event plane resolution of B-C sub events DECLARE_SOA_COLUMN(BMagField, bMagField, float); //! Magnetic field +DECLARE_SOA_COLUMN(IsRecINELgt0, isRecINELgt0, bool); //! Is reconstructed INEL>0 event // MC DECLARE_SOA_COLUMN(IsVtxIn10, isVtxIn10, bool); //! Vtx10 DECLARE_SOA_COLUMN(IsINELgt0, isINELgt0, bool); //! INEL>0 @@ -75,16 +76,20 @@ DECLARE_SOA_COLUMN(IsTriggerTVX, isTriggerTVX, bool); //! TriggerTVX DECLARE_SOA_COLUMN(IsInSel8, isInSel8, bool); //! InSel8 DECLARE_SOA_COLUMN(IsInAfterAllCuts, isInAfterAllCuts, bool); //! InAfterAllCuts DECLARE_SOA_COLUMN(ImpactParameter, impactParameter, float); //! ImpactParameter +DECLARE_SOA_COLUMN(MCMultiplicity, mcMultiplicity, float); //! MC Multiplicity } // namespace resocollision DECLARE_SOA_TABLE(ResoCollisions, "AOD", "RESOCOLLISION", o2::soa::Index<>, o2::aod::mult::MultNTracksPV, + o2::aod::mult::MultNTracksPVeta1, + o2::aod::mult::MultNTracksPVetaHalf, collision::PosX, collision::PosY, collision::PosZ, resocollision::Cent, - resocollision::BMagField); + resocollision::BMagField, + resocollision::IsRecINELgt0); using ResoCollision = ResoCollisions::iterator; DECLARE_SOA_TABLE(ResoCollisionColls, "AOD", "RESOCOLLISIONCOL", @@ -98,7 +103,8 @@ DECLARE_SOA_TABLE(ResoMCCollisions, "AOD", "RESOMCCOLLISION", resocollision::IsTriggerTVX, resocollision::IsInSel8, resocollision::IsInAfterAllCuts, - resocollision::ImpactParameter); + resocollision::ImpactParameter, + resocollision::MCMultiplicity); using ResoMCCollision = ResoMCCollisions::iterator; DECLARE_SOA_TABLE(ResoSpheroCollisions, "AOD", "RESOSPHEROCOLLISION", @@ -229,6 +235,8 @@ DECLARE_SOA_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, bool); DECLARE_SOA_COLUMN(ProducedByGenerator, producedByGenerator, bool); DECLARE_SOA_COLUMN(MotherId, motherId, int); //! Id of the mother particle DECLARE_SOA_COLUMN(MotherPDG, motherPDG, int); //! PDG code of the mother particle +DECLARE_SOA_COLUMN(MotherPt, motherPt, float); //! Pt of the mother particle +DECLARE_SOA_COLUMN(MotherRap, motherRap, float); //! Rapidity of the mother particle DECLARE_SOA_COLUMN(DaughterPDG1, daughterPDG1, int); //! PDG code of the first Daughter particle DECLARE_SOA_COLUMN(DaughterPDG2, daughterPDG2, int); //! PDG code of the second Daughter particle DECLARE_SOA_COLUMN(DaughterID1, daughterID1, int); //! Id of the first Daughter particle @@ -812,6 +820,8 @@ DECLARE_SOA_TABLE(ResoMCV0s, "AOD", "RESOMCV0", mcparticle::PdgCode, resodaughter::MotherId, resodaughter::MotherPDG, + resodaughter::MotherPt, + resodaughter::MotherRap, resodaughter::DaughterID1, resodaughter::DaughterID2, resodaughter::DaughterPDG1, @@ -824,6 +834,8 @@ DECLARE_SOA_TABLE(ResoMCCascades, "AOD", "RESOMCCASCADE", mcparticle::PdgCode, resodaughter::MotherId, resodaughter::MotherPDG, + resodaughter::MotherPt, + resodaughter::MotherRap, resodaughter::BachTrkID, resodaughter::V0ID, resodaughter::DaughterPDG1, diff --git a/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx b/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx index 79514728d10..dfaf46f0e39 100644 --- a/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx +++ b/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx @@ -11,12 +11,13 @@ /// /// \file resonanceInitializer.cxx /// \brief Initializes variables for the resonance candidate producers -/// \author Bong-Hwi Lim +/// \author Bong-Hwi Lim , Minjae Kim /// #include "PWGLF/DataModel/LFResonanceTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/Utils/collisionCuts.h" +#include "PWGLF/DataModel/mcCentrality.h" #include "Common/Core/EventPlaneHelper.h" #include "Common/Core/RecoDecay.h" @@ -24,6 +25,7 @@ #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/McCollisionExtra.h" #include "Common/DataModel/Qvectors.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -205,6 +207,27 @@ struct ResonanceInitializer { Configurable cfgSecondaryCrossMassCutWindow{"cfgSecondaryCrossMassCutWindow", 0.05, "Secondary inv mass selection window with (anti)lambda hypothesis"}; } SecondaryCuts; + struct : ConfigurableGroup { + Configurable cfgGenMult05{"cfgGenMult05",true, "GenEvent: multiplicity in |eta| < 0.5"}; + Configurable cfgGenMult10{"cfgGenMult10",false, "GenEvent: multiplicity in |eta| < 1.0"}; + Configurable cfgGenMultPercentile{"cfgGenMultPercentile",false, "Inherit Centrality(Multiplicity) percentile from MC collision only using LF-mc-centrality task"}; + + Configurable isZvtxcutGen{"isZvtxcutGen", true, "z-vertex cut for the GenCollision"}; + Configurable cutzvertexGen{"cutzvertexGen", 10.0f, "z-vertex cut for the GenCollision"}; + Configurable checkIsTrueINELgt0{"checkIsTrueINELgt0", true, "Check true INEL>0 for the Gen. Collision"}; + + ConfigurableAxis ptAxisGen{"ptAxisGen", {400, 0.0f, 20.0f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis multNTracksAxis{"multNTracksAxis", {500, 0.0f, +5000.0f}, "Number of charged particles"}; + ConfigurableAxis impactParameterAxis{"impactParameterAxis", {500, 0, 50}, "IP (fm)"}; + + Configurable isDaughterCheck{"isDaughterCheck", 1, "Check if the mother has the correct daughters when it is considered"}; + Configurable cfgRapidityCutGen{"cfgRapidityCutGen", 0.5, "Rapidity cut for the truth particle"}; + Configurable pdgTruthMother{"pdgTruthMother", 3324, "pdgcode for the truth mother e.g. Xi(1530) (3324)"}; + Configurable pdgTruthDaughter1{"pdgTruthDaughter1", 3312, "pdgcode for the daughter 1, e.g. Xi- 3312"}; + Configurable pdgTruthDaughter2{"pdgTruthDaughter2", 211, "pdgcode for the daughter 2, e.g. pi+ 211"}; + } GenCuts; + Configurable checkIsRecINELgt0{"checkIsRecINELgt0", true, "Check rec INEL>0 for the Rec. Collision"}; + HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject}; // Pre-filters for efficient process @@ -245,6 +268,7 @@ struct ResonanceInitializer { || (nabs(aod::mcparticle::pdgCode) == 123324); // Xi(1820)-0 using ResoEvents = soa::Join; + using ResoEvents001 = soa::Join; using ResoRun2Events = soa::Join; using ResoEventsMC = soa::Join; using ResoRun2EventsMC = soa::Join; @@ -931,6 +955,22 @@ struct ResonanceInitializer { } return lMothersPDGs; }; + auto getMothersPt = [&](auto const& theMcParticle) { + std::vector lMothersPts{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pdgcode lMother: %f", lMother.pt()); + lMothersPts.push_back(lMother.pt()); + } + return lMothersPts; + }; + auto getMothersRap = [&](auto const& theMcParticle) { + std::vector lMothersRaps{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother rap lMother: %f", lMother.y()); + lMothersRaps.push_back(lMother.y()); + } + return lMothersRaps; + }; auto getDaughtersIndeces = [&](auto const& theMcParticle) { std::vector lDaughtersIndeces{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { @@ -954,6 +994,8 @@ struct ResonanceInitializer { // ------ std::vector mothers = {-1, -1}; std::vector motherPDGs = {-1, -1}; + std::vector mothersPts = {-1.0f,-1.0f}; + std::vector mothersRaps = {-1.0f,-1.0f}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; if (v0.has_mcParticle()) { @@ -961,10 +1003,13 @@ struct ResonanceInitializer { if (v0mc.has_mothers()) { mothers = getMothersIndeces(v0mc); motherPDGs = getMothersPDGCodes(v0mc); + mothersPts = getMothersPt(v0mc); + mothersRaps = getMothersRap(v0mc); } while (mothers.size() > 2) { mothers.pop_back(); motherPDGs.pop_back(); + mothersPts.pop_back(); } if (v0mc.has_daughters()) { daughters = getDaughtersIndeces(v0mc); @@ -978,6 +1023,8 @@ struct ResonanceInitializer { reso2mcv0s(v0mc.pdgCode(), mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -988,6 +1035,8 @@ struct ResonanceInitializer { reso2mcv0s(0, mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1017,6 +1066,22 @@ struct ResonanceInitializer { } return lMothersPDGs; }; + auto getMothersPt = [&](auto const& theMcParticle) { + std::vector lMothersPts{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pdgcode lMother: %f", lMother.pt()); + lMothersPts.push_back(lMother.pt()); + } + return lMothersPts; + }; + auto getMothersRap = [&](auto const& theMcParticle) { + std::vector lMothersRaps{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother rap lMother: %f", lMother.y()); + lMothersRaps.push_back(lMother.y()); + } + return lMothersRaps; + }; auto getDaughtersIndeces = [&](auto const& theMcParticle) { std::vector lDaughtersIndeces{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { @@ -1042,15 +1107,20 @@ struct ResonanceInitializer { std::vector motherPDGs = {-1, -1}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; + std::vector mothersPts = {-1.0f,-1.0f}; + std::vector mothersRaps = {-1.0f,-1.0f}; if (casc.has_mcParticle()) { auto cascmc = casc.mcParticle(); if (cascmc.has_mothers()) { mothers = getMothersIndeces(cascmc); + mothersPts = getMothersPt(cascmc); + mothersRaps = getMothersRap(cascmc); motherPDGs = getMothersPDGCodes(cascmc); } while (mothers.size() > 2) { mothers.pop_back(); motherPDGs.pop_back(); + mothersPts.pop_back(); } if (cascmc.has_daughters()) { daughters = getDaughtersIndeces(cascmc); @@ -1064,6 +1134,8 @@ struct ResonanceInitializer { reso2mccascades(cascmc.pdgCode(), mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1074,6 +1146,8 @@ struct ResonanceInitializer { reso2mccascades(0, mothers[0], motherPDGs[0], + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1112,20 +1186,53 @@ struct ResonanceInitializer { } } + template + void fillMCGenParticles(TotalMCParts const& mcParticles, MCCentGen const& Cent, MCMultGen const& MCMult, MCIPGen const& IP, evtType const& eventType) + { + for (auto const& mcPart : mcParticles) { + + if (std::abs(mcPart.pdgCode()) != GenCuts.pdgTruthMother || std::abs(mcPart.y()) >= GenCuts.cfgRapidityCutGen) + continue; + std::vector daughterPDGs; + if (mcPart.has_daughters()) { + auto daughter01 = mcParticles.rawIteratorAt(mcPart.daughtersIds()[0] - mcParticles.offset()); + auto daughter02 = mcParticles.rawIteratorAt(mcPart.daughtersIds()[1] - mcParticles.offset()); + daughterPDGs = {daughter01.pdgCode(), daughter02.pdgCode()}; + } else { + daughterPDGs = {-1, -1}; + } + + if (GenCuts.isDaughterCheck) { + bool pass1 = std::abs(daughterPDGs[0]) == GenCuts.pdgTruthDaughter1 || std::abs(daughterPDGs[1]) == GenCuts.pdgTruthDaughter1; + bool pass2 = std::abs(daughterPDGs[0]) == GenCuts.pdgTruthDaughter2 || std::abs(daughterPDGs[1]) == GenCuts.pdgTruthDaughter2; + if (!pass1 || !pass2) + continue; + } + if (mcPart.pdgCode() > 0) // Consider INELt0 or INEL + qaRegistry.fill(HIST("EventGen/h5ResonanceTruth"), eventType, mcPart.pt(), Cent, MCMult,IP); + else + qaRegistry.fill(HIST("EventGen/h5ResonanceTruthAnti"), eventType, mcPart.pt(), Cent, MCMult, IP); + + daughterPDGs.clear(); + } + } + template - void fillMCCollision(MCCol const& mccol, MCPart const& mcparts, float impactpar = -999.0) + void fillMCCollision(MCCol const& mccol, MCPart const& mcparts, float impactpar = -999.0, float mult = -1.0) { auto centrality = 0.0; if constexpr (!isRun2) centrality = centEst(mccol); else centrality = mccol.centRun2V0M(); - bool inVtx10 = (std::abs(mccol.mcCollision().posZ()) > 10.) ? false : true; + + // bool inVtx10 = (std::abs(mccol.mcCollision().posZ()) > 10.) ? false : true; -> Gen. level informations will be processed in processMCGen + bool inVtx10 = (std::abs(mccol.posZ()) > 10.) ? false : true; bool isTrueINELgt0 = isTrueINEL0(mccol, mcparts); bool isTriggerTVX = mccol.selection_bit(aod::evsel::kIsTriggerTVX); bool isSel8 = mccol.sel8(); bool isSelected = colCuts.isSelected(mccol); - resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, impactpar); + resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, impactpar, mult); // QA for Trigger efficiency qaRegistry.fill(HIST("Event/hMCEventIndices"), centrality, aod::resocollision::kINEL); @@ -1239,6 +1346,23 @@ struct ResonanceInitializer { qaRegistry.add("hGoodMCCascIndices", "hGoodMCCascIndices", kTH1F, {idxAxis}); qaRegistry.add("Phi", "#phi distribution", kTH1F, {{65, -0.1, 6.4}}); } + + + TString hNEventsMCLabels[4] = {"All", "z vrtx", "INEL", "INEL>0"}; + if (doprocessMCgen){ + AxisSpec centAxisGen = {binsCent, "Centrality (%)"}; + qaRegistry.add("EventGen/hNEventsMC", "EventGen/hNEventsMC", kTH1D, {{4, 0.0f, 4.0f}}); + for (int n = 1; n <= qaRegistry.get(HIST("EventGen/hNEventsMC"))->GetNbinsX(); n++) { + qaRegistry.get(HIST("EventGen/hNEventsMC"))->GetXaxis()->SetBinLabel(n, hNEventsMCLabels[n - 1]); + } + qaRegistry.add("EventGen/h5ResonanceTruth", "EventGen/h5ResonanceTruth", kTHnSparseD, {{2, 0.0f, 2.0f}, GenCuts.ptAxisGen,centAxisGen,GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/h5ResonanceTruthAnti", "EventGen/h5ResonanceTruthAnti", kTHnSparseD, {{2, 0.0f, 2.0f}, GenCuts.ptAxisGen, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/hZCollisionGen", "EventGen/hZCollisionGen", kTH1D, {{100, -20.0f, 20.0f}}); + + qaRegistry.add("EventGen/h4MultCent_genMC", "EventGen/h4MultCent_genMC", kTHnSparseD, {{2, 0.0f, 2.0f}, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/h4MultCent_recMC", "EventGen/h4MultCent_recMC", kTHnSparseD, {{2, 0.0f, 2.0f}, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/h2CentralityVsMultMC", "EventGen/h2CentralityVsMultMC", kTH2D, {centAxisGen, GenCuts.multNTracksAxis}); + } } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) // Simple copy from LambdaKzeroFinder.cxx @@ -1303,7 +1427,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1329,7 +1453,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1356,7 +1480,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1383,7 +1507,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1409,7 +1533,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1424,7 +1548,7 @@ struct ResonanceInitializer { } PROCESS_SWITCH(ResonanceInitializer, processTrackV0DataRun2, "Process for data", false); - void processTrackV0CascData(ResoEvents::iterator const& collision, + void processTrackV0CascData(ResoEvents001::iterator const& collision, soa::Filtered const& tracks, ResoV0s const& V0s, ResoCascades const& Cascades, @@ -1438,8 +1562,11 @@ struct ResonanceInitializer { if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) return; colCuts.fillQA(collision); + bool isRecINELgt0 = 0; + if(checkIsRecINELgt0) + isRecINELgt0 = collision.isInelGt0(); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(collision.multNTracksPV(),collision.multNTracksPVeta1(),collision.multNTracksPVetaHalf(),collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, isRecINELgt0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1466,7 +1593,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1493,7 +1620,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1525,7 +1652,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1552,7 +1679,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1582,7 +1709,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1611,7 +1738,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1632,25 +1759,45 @@ struct ResonanceInitializer { } PROCESS_SWITCH(ResonanceInitializer, processTrackV0MCRun2, "Process for MC", false); - void processTrackV0CascMC(soa::Join::iterator const& collision, - aod::McCollisions const&, soa::Filtered const& tracks, + void processTrackV0CascMC(soa::Join::iterator const& collision, + soa::Join const& mcCollisions, soa::Filtered const& tracks, ResoV0sMC const& V0s, ResoCascadesMC const& Cascades, aod::McParticles const& mcParticles, aod::BCsWithTimestamps const&) { + auto bc = collision.bc_as(); /// adding timestamp to access magnetic field later initCCDB(bc); if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) return; colCuts.fillQA(collision); - - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz); + + float Cent = 100.5f; + + const auto mcId = collision.mcCollisionId(); + auto mcCollision = mcCollisions.iteratorAt(mcId); + float impactpar = mcCollision.impactParameter(); + float mult = -1.0f; + if(GenCuts.cfgGenMultPercentile) + Cent = mcCollision.centFT0M(); + else + Cent = centEst(collision); + + bool isRecINELgt0 = 0; + if(checkIsRecINELgt0) + isRecINELgt0 = collision.isInelGt0(); + + resoCollisions(collision.multNTracksPV(),collision.multNTracksPVeta1(),collision.multNTracksPVetaHalf(), collision.posX(), collision.posY(), collision.posZ(), Cent, dBz, isRecINELgt0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } resoSpheroCollisions(computeSpherocity(tracks, trackSphMin, trackSphDef)); resoEvtPlCollisions(0, 0, 0, 0); - fillMCCollision(collision, mcParticles); + + if(GenCuts.cfgGenMult05) mult = mcCollision.multMCNParticlesEta05(); + else if (GenCuts.cfgGenMult10) mult = mcCollision.multMCNParticlesEta10(); + + fillMCCollision(collision, mcParticles, impactpar, mult); // Loop over tracks fillTracks(collision, tracks); @@ -1661,13 +1808,100 @@ struct ResonanceInitializer { fillCascades(collision, Cascades, tracks); // Loop over all MC particles - auto mcParts = selectedMCParticles->sliceBy(perMcCollision, collision.mcCollision().globalIndex()); + auto mcParts = selectedMCParticles->sliceBy(perMcCollision, mcId); fillMCParticles(mcParts, mcParticles); } PROCESS_SWITCH(ResonanceInitializer, processTrackV0CascMC, "Process for MC", false); + +// Following the discussions at the PAG meeting (https://indico.cern.ch/event/1583408/) +// we have introduced an auxiliary task that, when the resonanceInitializer.cxx is used, +// Only consider N_rec / N_gen i.e. not consider level of N_gen at least once + void processMCgen(soa::Join::iterator const& mcCollision, + aod::McParticles const& mcParticles, + const soa::SmallGroups>& collisions, + aod::BCsWithTimestamps const&) + { + auto bc = mcCollision.bc_as(); /// adding timestamp to access magnetic field later + initCCDB(bc); + + auto getCentGen = [&]() { + if (cfgMultName.value == "FT0M") { // FT0A,C results wiill be updated later when CCDB is available + return mcCollision.centFT0M(); + } + return 100.5f; + }; + auto getCentReco = [&](auto const& col) { + if (cfgMultName.value == "FT0M") { + return col.centFT0M(); + } else if (cfgMultName.value == "FT0C") { + return col.centFT0C(); + } else if (cfgMultName.value == "FT0A") { + return col.centFT0A(); + } + return 100.5f; + }; + + float cent = getCentGen(); + float IP = mcCollision.impactParameter(); + float mult = -1; + if (GenCuts.cfgGenMult05) { + mult = mcCollision.multMCNParticlesEta05(); + } else if (GenCuts.cfgGenMult10) { + mult = mcCollision.multMCNParticlesEta10(); + } + + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 0.5); + + if (GenCuts.isZvtxcutGen && std::fabs(mcCollision.posZ()) > GenCuts.cutzvertexGen) { + return; + } + qaRegistry.fill(HIST("EventGen/hZCollisionGen"), mcCollision.posZ()); + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 1.5); + + int evType = 0; + + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 2.5); + if (GenCuts.checkIsTrueINELgt0 && mcCollision.isInelGt0()) { + evType++; + qaRegistry.fill(HIST("EventGen/hNEventsMC"), 3.5); + } + + bool atLeastOne = false; + int biggestNContribs = -1; + + float centReco = 100.5f; + for (const auto& collision : collisions) { + if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) + continue; + if(!colCuts.isSelected(collision)) // Bug is appeared in colCuts-> double counting in event QA histo, will be fixed later + continue; + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + centReco = getCentReco(collision); + } + + atLeastOne = true; + } + + if (GenCuts.cfgGenMultPercentile) { + fillMCGenParticles(mcParticles, cent, mult, IP, evType); + qaRegistry.fill(HIST("EventGen/h4MultCent_genMC"), evType, cent, mult, IP); + } else { + fillMCGenParticles(mcParticles, centReco, mult, IP, evType); + qaRegistry.fill(HIST("EventGen/h4MultCent_genMC"), evType, centReco, mult, IP); + qaRegistry.fill(HIST("EventGen/h2CentralityVsMultMC"),centReco, mult); + } + + if (atLeastOne) { + qaRegistry.fill(HIST("EventGen/h4MultCent_recMC"), evType, centReco, mult, IP); + } + + } + PROCESS_SWITCH(ResonanceInitializer, processMCgen, "Process for MCGen", true); + void processTrackV0CascMCRun2(soa::Join::iterator const& collision, - aod::McCollisions const&, soa::Filtered const& tracks, + aod::McCollisions const&, ResoTracksMC const& tracks, ResoV0sMC const& V0s, ResoCascadesMC const& Cascades, aod::McParticles const& mcParticles, BCsWithRun2Info const&) @@ -1675,7 +1909,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } diff --git a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx index a1ade00c627..857ea2180f6 100644 --- a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx +++ b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx @@ -475,7 +475,7 @@ struct ResonanceModuleInitializer { bool isTriggerTVX = mccol.selection_bit(aod::evsel::kIsTriggerTVX); bool isSel8 = mccol.sel8(); bool isSelected = colCuts.isSelected(mccol); - resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, mcCent); + resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, mcCent, -1.0f); // QA for Trigger efficiency qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINEL); @@ -546,7 +546,7 @@ struct ResonanceModuleInitializer { colCuts.fillQA(collision); centrality = centEst(collision); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz,0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -569,7 +569,7 @@ struct ResonanceModuleInitializer { colCuts.fillQARun2(collision); centrality = collision.centRun2V0M(); - resoCollisions(0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz); + resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz,0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1123,6 +1123,8 @@ struct ResonanceDaughterInitializer { reso2mcv0s(v0mc.pdgCode(), mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0], @@ -1133,6 +1135,8 @@ struct ResonanceDaughterInitializer { reso2mcv0s(0, mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0], @@ -1289,6 +1293,8 @@ struct ResonanceDaughterInitializer { reso2mccascades(cascmc.pdgCode(), mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0], @@ -1299,6 +1305,8 @@ struct ResonanceDaughterInitializer { reso2mccascades(0, mothers[0], motherPDGs[0], + 0, + 0, daughters[0], daughters[1], daughterPDGs[0], From eba3850a50cd8ec5ff85239f0d246201b021e624 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 3 Feb 2026 13:36:42 +0000 Subject: [PATCH 6/6] Please consider the following formatting changes --- PWGLF/DataModel/LFResonanceTables.h | 2 +- .../Resonances/resonanceInitializer.cxx | 141 +++++++++--------- .../Resonances/resonanceModuleInitializer.cxx | 4 +- 3 files changed, 73 insertions(+), 74 deletions(-) diff --git a/PWGLF/DataModel/LFResonanceTables.h b/PWGLF/DataModel/LFResonanceTables.h index 060fac0810f..8a04eb4ffdd 100644 --- a/PWGLF/DataModel/LFResonanceTables.h +++ b/PWGLF/DataModel/LFResonanceTables.h @@ -76,7 +76,7 @@ DECLARE_SOA_COLUMN(IsTriggerTVX, isTriggerTVX, bool); //! TriggerTVX DECLARE_SOA_COLUMN(IsInSel8, isInSel8, bool); //! InSel8 DECLARE_SOA_COLUMN(IsInAfterAllCuts, isInAfterAllCuts, bool); //! InAfterAllCuts DECLARE_SOA_COLUMN(ImpactParameter, impactParameter, float); //! ImpactParameter -DECLARE_SOA_COLUMN(MCMultiplicity, mcMultiplicity, float); //! MC Multiplicity +DECLARE_SOA_COLUMN(MCMultiplicity, mcMultiplicity, float); //! MC Multiplicity } // namespace resocollision DECLARE_SOA_TABLE(ResoCollisions, "AOD", "RESOCOLLISION", diff --git a/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx b/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx index dfaf46f0e39..e1c4b3cc5d8 100644 --- a/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx +++ b/PWGLF/TableProducer/Resonances/resonanceInitializer.cxx @@ -16,8 +16,8 @@ #include "PWGLF/DataModel/LFResonanceTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" -#include "PWGLF/Utils/collisionCuts.h" #include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/collisionCuts.h" #include "Common/Core/EventPlaneHelper.h" #include "Common/Core/RecoDecay.h" @@ -208,23 +208,23 @@ struct ResonanceInitializer { } SecondaryCuts; struct : ConfigurableGroup { - Configurable cfgGenMult05{"cfgGenMult05",true, "GenEvent: multiplicity in |eta| < 0.5"}; - Configurable cfgGenMult10{"cfgGenMult10",false, "GenEvent: multiplicity in |eta| < 1.0"}; - Configurable cfgGenMultPercentile{"cfgGenMultPercentile",false, "Inherit Centrality(Multiplicity) percentile from MC collision only using LF-mc-centrality task"}; - - Configurable isZvtxcutGen{"isZvtxcutGen", true, "z-vertex cut for the GenCollision"}; - Configurable cutzvertexGen{"cutzvertexGen", 10.0f, "z-vertex cut for the GenCollision"}; - Configurable checkIsTrueINELgt0{"checkIsTrueINELgt0", true, "Check true INEL>0 for the Gen. Collision"}; - - ConfigurableAxis ptAxisGen{"ptAxisGen", {400, 0.0f, 20.0f}, "#it{p}_{T} (GeV/#it{c})"}; - ConfigurableAxis multNTracksAxis{"multNTracksAxis", {500, 0.0f, +5000.0f}, "Number of charged particles"}; - ConfigurableAxis impactParameterAxis{"impactParameterAxis", {500, 0, 50}, "IP (fm)"}; - - Configurable isDaughterCheck{"isDaughterCheck", 1, "Check if the mother has the correct daughters when it is considered"}; - Configurable cfgRapidityCutGen{"cfgRapidityCutGen", 0.5, "Rapidity cut for the truth particle"}; - Configurable pdgTruthMother{"pdgTruthMother", 3324, "pdgcode for the truth mother e.g. Xi(1530) (3324)"}; - Configurable pdgTruthDaughter1{"pdgTruthDaughter1", 3312, "pdgcode for the daughter 1, e.g. Xi- 3312"}; - Configurable pdgTruthDaughter2{"pdgTruthDaughter2", 211, "pdgcode for the daughter 2, e.g. pi+ 211"}; + Configurable cfgGenMult05{"cfgGenMult05", true, "GenEvent: multiplicity in |eta| < 0.5"}; + Configurable cfgGenMult10{"cfgGenMult10", false, "GenEvent: multiplicity in |eta| < 1.0"}; + Configurable cfgGenMultPercentile{"cfgGenMultPercentile", false, "Inherit Centrality(Multiplicity) percentile from MC collision only using LF-mc-centrality task"}; + + Configurable isZvtxcutGen{"isZvtxcutGen", true, "z-vertex cut for the GenCollision"}; + Configurable cutzvertexGen{"cutzvertexGen", 10.0f, "z-vertex cut for the GenCollision"}; + Configurable checkIsTrueINELgt0{"checkIsTrueINELgt0", true, "Check true INEL>0 for the Gen. Collision"}; + + ConfigurableAxis ptAxisGen{"ptAxisGen", {400, 0.0f, 20.0f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis multNTracksAxis{"multNTracksAxis", {500, 0.0f, +5000.0f}, "Number of charged particles"}; + ConfigurableAxis impactParameterAxis{"impactParameterAxis", {500, 0, 50}, "IP (fm)"}; + + Configurable isDaughterCheck{"isDaughterCheck", 1, "Check if the mother has the correct daughters when it is considered"}; + Configurable cfgRapidityCutGen{"cfgRapidityCutGen", 0.5, "Rapidity cut for the truth particle"}; + Configurable pdgTruthMother{"pdgTruthMother", 3324, "pdgcode for the truth mother e.g. Xi(1530) (3324)"}; + Configurable pdgTruthDaughter1{"pdgTruthDaughter1", 3312, "pdgcode for the daughter 1, e.g. Xi- 3312"}; + Configurable pdgTruthDaughter2{"pdgTruthDaughter2", 211, "pdgcode for the daughter 2, e.g. pi+ 211"}; } GenCuts; Configurable checkIsRecINELgt0{"checkIsRecINELgt0", true, "Check rec INEL>0 for the Rec. Collision"}; @@ -268,7 +268,7 @@ struct ResonanceInitializer { || (nabs(aod::mcparticle::pdgCode) == 123324); // Xi(1820)-0 using ResoEvents = soa::Join; - using ResoEvents001 = soa::Join; + using ResoEvents001 = soa::Join; using ResoRun2Events = soa::Join; using ResoEventsMC = soa::Join; using ResoRun2EventsMC = soa::Join; @@ -994,8 +994,8 @@ struct ResonanceInitializer { // ------ std::vector mothers = {-1, -1}; std::vector motherPDGs = {-1, -1}; - std::vector mothersPts = {-1.0f,-1.0f}; - std::vector mothersRaps = {-1.0f,-1.0f}; + std::vector mothersPts = {-1.0f, -1.0f}; + std::vector mothersRaps = {-1.0f, -1.0f}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; if (v0.has_mcParticle()) { @@ -1107,13 +1107,13 @@ struct ResonanceInitializer { std::vector motherPDGs = {-1, -1}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; - std::vector mothersPts = {-1.0f,-1.0f}; - std::vector mothersRaps = {-1.0f,-1.0f}; + std::vector mothersPts = {-1.0f, -1.0f}; + std::vector mothersRaps = {-1.0f, -1.0f}; if (casc.has_mcParticle()) { auto cascmc = casc.mcParticle(); if (cascmc.has_mothers()) { mothers = getMothersIndeces(cascmc); - mothersPts = getMothersPt(cascmc); + mothersPts = getMothersPt(cascmc); mothersRaps = getMothersRap(cascmc); motherPDGs = getMothersPDGCodes(cascmc); } @@ -1186,7 +1186,7 @@ struct ResonanceInitializer { } } - template + template void fillMCGenParticles(TotalMCParts const& mcParticles, MCCentGen const& Cent, MCMultGen const& MCMult, MCIPGen const& IP, evtType const& eventType) { for (auto const& mcPart : mcParticles) { @@ -1209,7 +1209,7 @@ struct ResonanceInitializer { continue; } if (mcPart.pdgCode() > 0) // Consider INELt0 or INEL - qaRegistry.fill(HIST("EventGen/h5ResonanceTruth"), eventType, mcPart.pt(), Cent, MCMult,IP); + qaRegistry.fill(HIST("EventGen/h5ResonanceTruth"), eventType, mcPart.pt(), Cent, MCMult, IP); else qaRegistry.fill(HIST("EventGen/h5ResonanceTruthAnti"), eventType, mcPart.pt(), Cent, MCMult, IP); @@ -1225,7 +1225,7 @@ struct ResonanceInitializer { centrality = centEst(mccol); else centrality = mccol.centRun2V0M(); - + // bool inVtx10 = (std::abs(mccol.mcCollision().posZ()) > 10.) ? false : true; -> Gen. level informations will be processed in processMCGen bool inVtx10 = (std::abs(mccol.posZ()) > 10.) ? false : true; bool isTrueINELgt0 = isTrueINEL0(mccol, mcparts); @@ -1346,23 +1346,22 @@ struct ResonanceInitializer { qaRegistry.add("hGoodMCCascIndices", "hGoodMCCascIndices", kTH1F, {idxAxis}); qaRegistry.add("Phi", "#phi distribution", kTH1F, {{65, -0.1, 6.4}}); } - - + TString hNEventsMCLabels[4] = {"All", "z vrtx", "INEL", "INEL>0"}; - if (doprocessMCgen){ + if (doprocessMCgen) { AxisSpec centAxisGen = {binsCent, "Centrality (%)"}; qaRegistry.add("EventGen/hNEventsMC", "EventGen/hNEventsMC", kTH1D, {{4, 0.0f, 4.0f}}); for (int n = 1; n <= qaRegistry.get(HIST("EventGen/hNEventsMC"))->GetNbinsX(); n++) { - qaRegistry.get(HIST("EventGen/hNEventsMC"))->GetXaxis()->SetBinLabel(n, hNEventsMCLabels[n - 1]); + qaRegistry.get(HIST("EventGen/hNEventsMC"))->GetXaxis()->SetBinLabel(n, hNEventsMCLabels[n - 1]); } - qaRegistry.add("EventGen/h5ResonanceTruth", "EventGen/h5ResonanceTruth", kTHnSparseD, {{2, 0.0f, 2.0f}, GenCuts.ptAxisGen,centAxisGen,GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); + qaRegistry.add("EventGen/h5ResonanceTruth", "EventGen/h5ResonanceTruth", kTHnSparseD, {{2, 0.0f, 2.0f}, GenCuts.ptAxisGen, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); qaRegistry.add("EventGen/h5ResonanceTruthAnti", "EventGen/h5ResonanceTruthAnti", kTHnSparseD, {{2, 0.0f, 2.0f}, GenCuts.ptAxisGen, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); qaRegistry.add("EventGen/hZCollisionGen", "EventGen/hZCollisionGen", kTH1D, {{100, -20.0f, 20.0f}}); qaRegistry.add("EventGen/h4MultCent_genMC", "EventGen/h4MultCent_genMC", kTHnSparseD, {{2, 0.0f, 2.0f}, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); qaRegistry.add("EventGen/h4MultCent_recMC", "EventGen/h4MultCent_recMC", kTHnSparseD, {{2, 0.0f, 2.0f}, centAxisGen, GenCuts.multNTracksAxis, GenCuts.impactParameterAxis}); qaRegistry.add("EventGen/h2CentralityVsMultMC", "EventGen/h2CentralityVsMultMC", kTH2D, {centAxisGen, GenCuts.multNTracksAxis}); - } + } } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) // Simple copy from LambdaKzeroFinder.cxx @@ -1427,7 +1426,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1453,7 +1452,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1480,7 +1479,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1507,7 +1506,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1533,7 +1532,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1563,10 +1562,10 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); bool isRecINELgt0 = 0; - if(checkIsRecINELgt0) + if (checkIsRecINELgt0) isRecINELgt0 = collision.isInelGt0(); - resoCollisions(collision.multNTracksPV(),collision.multNTracksPVeta1(),collision.multNTracksPVetaHalf(),collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, isRecINELgt0); + resoCollisions(collision.multNTracksPV(), collision.multNTracksPVeta1(), collision.multNTracksPVetaHalf(), collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, isRecINELgt0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1593,7 +1592,7 @@ struct ResonanceInitializer { return; colCuts.fillQARun2(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1620,7 +1619,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1652,7 +1651,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1679,7 +1678,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1709,7 +1708,7 @@ struct ResonanceInitializer { return; colCuts.fillQA(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1738,7 +1737,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -1771,31 +1770,33 @@ struct ResonanceInitializer { if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) return; colCuts.fillQA(collision); - + float Cent = 100.5f; const auto mcId = collision.mcCollisionId(); auto mcCollision = mcCollisions.iteratorAt(mcId); float impactpar = mcCollision.impactParameter(); float mult = -1.0f; - if(GenCuts.cfgGenMultPercentile) + if (GenCuts.cfgGenMultPercentile) Cent = mcCollision.centFT0M(); else Cent = centEst(collision); - + bool isRecINELgt0 = 0; - if(checkIsRecINELgt0) + if (checkIsRecINELgt0) isRecINELgt0 = collision.isInelGt0(); - - resoCollisions(collision.multNTracksPV(),collision.multNTracksPVeta1(),collision.multNTracksPVetaHalf(), collision.posX(), collision.posY(), collision.posZ(), Cent, dBz, isRecINELgt0); + + resoCollisions(collision.multNTracksPV(), collision.multNTracksPVeta1(), collision.multNTracksPVetaHalf(), collision.posX(), collision.posY(), collision.posZ(), Cent, dBz, isRecINELgt0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } resoSpheroCollisions(computeSpherocity(tracks, trackSphMin, trackSphDef)); resoEvtPlCollisions(0, 0, 0, 0); - if(GenCuts.cfgGenMult05) mult = mcCollision.multMCNParticlesEta05(); - else if (GenCuts.cfgGenMult10) mult = mcCollision.multMCNParticlesEta10(); + if (GenCuts.cfgGenMult05) + mult = mcCollision.multMCNParticlesEta05(); + else if (GenCuts.cfgGenMult10) + mult = mcCollision.multMCNParticlesEta10(); fillMCCollision(collision, mcParticles, impactpar, mult); @@ -1812,16 +1813,16 @@ struct ResonanceInitializer { fillMCParticles(mcParts, mcParticles); } PROCESS_SWITCH(ResonanceInitializer, processTrackV0CascMC, "Process for MC", false); - -// Following the discussions at the PAG meeting (https://indico.cern.ch/event/1583408/) -// we have introduced an auxiliary task that, when the resonanceInitializer.cxx is used, -// Only consider N_rec / N_gen i.e. not consider level of N_gen at least once + + // Following the discussions at the PAG meeting (https://indico.cern.ch/event/1583408/) + // we have introduced an auxiliary task that, when the resonanceInitializer.cxx is used, + // Only consider N_rec / N_gen i.e. not consider level of N_gen at least once void processMCgen(soa::Join::iterator const& mcCollision, - aod::McParticles const& mcParticles, + aod::McParticles const& mcParticles, const soa::SmallGroups>& collisions, aod::BCsWithTimestamps const&) { - auto bc = mcCollision.bc_as(); /// adding timestamp to access magnetic field later + auto bc = mcCollision.bc_as(); /// adding timestamp to access magnetic field later initCCDB(bc); auto getCentGen = [&]() { @@ -1866,19 +1867,19 @@ struct ResonanceInitializer { qaRegistry.fill(HIST("EventGen/hNEventsMC"), 3.5); } - bool atLeastOne = false; - int biggestNContribs = -1; + bool atLeastOne = false; + int biggestNContribs = -1; float centReco = 100.5f; for (const auto& collision : collisions) { if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) continue; - if(!colCuts.isSelected(collision)) // Bug is appeared in colCuts-> double counting in event QA histo, will be fixed later + if (!colCuts.isSelected(collision)) // Bug is appeared in colCuts-> double counting in event QA histo, will be fixed later continue; - if (biggestNContribs < collision.multPVTotalContributors()) { - biggestNContribs = collision.multPVTotalContributors(); - centReco = getCentReco(collision); - } + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + centReco = getCentReco(collision); + } atLeastOne = true; } @@ -1889,17 +1890,15 @@ struct ResonanceInitializer { } else { fillMCGenParticles(mcParticles, centReco, mult, IP, evType); qaRegistry.fill(HIST("EventGen/h4MultCent_genMC"), evType, centReco, mult, IP); - qaRegistry.fill(HIST("EventGen/h2CentralityVsMultMC"),centReco, mult); + qaRegistry.fill(HIST("EventGen/h2CentralityVsMultMC"), centReco, mult); } if (atLeastOne) { qaRegistry.fill(HIST("EventGen/h4MultCent_recMC"), evType, centReco, mult, IP); } - } PROCESS_SWITCH(ResonanceInitializer, processMCgen, "Process for MCGen", true); - void processTrackV0CascMCRun2(soa::Join::iterator const& collision, aod::McCollisions const&, ResoTracksMC const& tracks, ResoV0sMC const& V0s, @@ -1909,7 +1908,7 @@ struct ResonanceInitializer { // auto bc = collision.bc_as(); colCuts.fillQARun2(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), collision.centRun2V0M(), dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } diff --git a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx index 857ea2180f6..5a4bc9f2be1 100644 --- a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx +++ b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx @@ -546,7 +546,7 @@ struct ResonanceModuleInitializer { colCuts.fillQA(collision); centrality = centEst(collision); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz,0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); } @@ -569,7 +569,7 @@ struct ResonanceModuleInitializer { colCuts.fillQARun2(collision); centrality = collision.centRun2V0M(); - resoCollisions(0,0,0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz,0); + resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz, 0); if (!cfgBypassCollIndexFill) { resoCollisionColls(collision.globalIndex()); }