diff --git a/PWGHF/D2H/Tasks/CMakeLists.txt b/PWGHF/D2H/Tasks/CMakeLists.txt index 76caf26673f..a6bab371fde 100644 --- a/PWGHF/D2H/Tasks/CMakeLists.txt +++ b/PWGHF/D2H/Tasks/CMakeLists.txt @@ -71,7 +71,7 @@ o2physics_add_dpl_workflow(task-charm-reso-to-d-trk-reduced o2physics_add_dpl_workflow(task-d0 SOURCES taskD0.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::SGCutParHolder COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(task-directed-flow-charm-hadrons @@ -81,7 +81,7 @@ o2physics_add_dpl_workflow(task-directed-flow-charm-hadrons o2physics_add_dpl_workflow(task-dplus SOURCES taskDplus.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::SGCutParHolder COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(task-ds diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index 3f74a0b4efa..219142f3d78 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -14,6 +14,7 @@ /// /// \author Gian Michele Innocenti , CERN /// \author Vít Kučera , CERN +/// \author Minjung Kim , CERN #include "PWGHF/Core/CentralityEstimation.h" #include "PWGHF/Core/DecayChannels.h" @@ -24,6 +25,8 @@ #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGHF/DataModel/TrackIndexSkimmingTables.h" #include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGHF/Utils/utilsUpcHf.h" +#include "PWGUD/Core/UPCHelpers.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" @@ -58,6 +61,8 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::hf_centrality; using namespace o2::hf_occupancy; +using namespace o2::hf_evsel; +using namespace o2::analysis::hf_upc; /// D0 analysis task namespace @@ -91,26 +96,32 @@ struct HfTaskD0 { Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable irSource{"irSource", "ZNC hadronic", "Estimator of the interaction rate (Recommended: pp --> T0VTX, Pb-Pb --> ZNC hadronic)"}; + HfEventSelection hfEvSel; // event selection and monitoring + HfUpcGapThresholds upcThresholds; // UPC gap determination thresholds ctpRateFetcher mRateFetcher; SliceCache cache; Service ccdb; - using D0Candidates = soa::Join; - using D0CandidatesMc = soa::Join; - using D0CandidatesKF = soa::Join; - using D0CandidatesMcKF = soa::Join; + using D0Candidates = soa::Filtered>; + using D0CandidatesMc = soa::Filtered>; + using D0CandidatesKF = soa::Filtered>; + using D0CandidatesMcKF = soa::Filtered>; - using D0CandidatesMl = soa::Join; - using D0CandidatesMlMc = soa::Join; - using D0CandidatesMlKF = soa::Join; - using D0CandidatesMlMcKF = soa::Join; + using D0CandidatesMl = soa::Filtered>; + using D0CandidatesMlMc = soa::Filtered>; + using D0CandidatesMlKF = soa::Filtered>; + using D0CandidatesMlMcKF = soa::Filtered>; using Collisions = soa::Join; using CollisionsCent = soa::Join; using CollisionsWithMcLabels = soa::Join; using CollisionsWithMcLabelsCent = soa::Join; using TracksSelQuality = soa::Join; + + Filter filterD0Flag = (o2::aod::hf_track_index::hfflag & static_cast(BIT(aod::hf_cand_2prong::DecayType::D0ToPiK))) != static_cast(0); + + Preslice candD0PerCollision = aod::hf_cand::collisionId; PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; PresliceUnsorted colPerMcCollisionCent = aod::mccollisionlabel::mcCollisionId; @@ -142,6 +153,11 @@ struct HfTaskD0 { ConfigurableAxis thnConfigAxisMinItsNCls{"thnConfigAxisMinItsNCls", {5, 3, 8}, "axis for minimum ITS NCls of candidate prongs"}; ConfigurableAxis thnConfigAxisMinTpcNCrossedRows{"thnConfigAxisMinTpcNCrossedRows", {10, 70, 180}, "axis for minimum TPC NCls crossed rows of candidate prongs"}; ConfigurableAxis thnConfigAxisIR{"thnConfigAxisIR", {5000, 0, 500}, "Interaction rate (kHz)"}; + ConfigurableAxis thnConfigAxisGapType{"thnConfigAxisGapType", {7, -1.5, 5.5}, "axis for UPC gap type (see TrueGap enum in o2::aod::sgselector)"}; + ConfigurableAxis thnConfigAxisFT0{"thnConfigAxisFT0", {1001, -1.5, 999.5}, "axis for FT0 amplitude (a.u.)"}; + ConfigurableAxis thnConfigAxisFV0A{"thnConfigAxisFV0A", {2001, -1.5, 1999.5}, "axis for FV0-A amplitude (a.u.)"}; + ConfigurableAxis thnConfigAxisFDD{"thnConfigAxisFDD", {200, 0., 4000.}, "axis for FDD amplitude (a.u.)"}; + ConfigurableAxis thnConfigAxisZN{"thnConfigAxisZN", {510, -1.5, 49.5}, "axis for ZN energy (a.u.)"}; HistogramRegistry registry{ "registry", @@ -221,15 +237,15 @@ struct HfTaskD0 { void init(InitContext&) { - std::array doprocess{doprocessDataWithDCAFitterN, doprocessDataWithDCAFitterNCent, doprocessDataWithKFParticle, doprocessMcWithDCAFitterN, doprocessMcWithDCAFitterNCent, doprocessMcWithKFParticle, doprocessDataWithDCAFitterNMl, doprocessDataWithDCAFitterNMlCent, doprocessDataWithKFParticleMl, doprocessMcWithDCAFitterNMl, doprocessMcWithDCAFitterNMlCent, doprocessMcWithKFParticleMl}; + std::array doprocess{doprocessDataWithDCAFitterN, doprocessDataWithDCAFitterNCent, doprocessDataWithKFParticle, doprocessMcWithDCAFitterN, doprocessMcWithDCAFitterNCent, doprocessMcWithKFParticle, doprocessDataWithDCAFitterNMl, doprocessDataWithDCAFitterNMlCent, doprocessDataWithKFParticleMl, doprocessMcWithDCAFitterNMl, doprocessMcWithDCAFitterNMlCent, doprocessMcWithKFParticleMl, doprocessDataWithDCAFitterNWithUpc, doprocessDataWithDCAFitterNMlWithUpc}; if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) == 0) { LOGP(fatal, "At least one process function should be enabled at a time."); } if ((doprocessDataWithDCAFitterN || doprocessDataWithDCAFitterNCent || doprocessMcWithDCAFitterN || doprocessMcWithDCAFitterNCent || doprocessDataWithDCAFitterNMl || doprocessDataWithDCAFitterNMlCent || doprocessMcWithDCAFitterNMl || doprocessMcWithDCAFitterNMlCent) && (doprocessDataWithKFParticle || doprocessMcWithKFParticle || doprocessDataWithKFParticleMl || doprocessMcWithKFParticleMl)) { LOGP(fatal, "DCAFitterN and KFParticle can not be enabled at a time."); } - if ((storeCentrality || storeOccupancyAndIR) && !(doprocessDataWithDCAFitterNCent || doprocessMcWithDCAFitterNCent || doprocessDataWithDCAFitterNMlCent || doprocessMcWithDCAFitterNMlCent)) { - LOGP(fatal, "Can't enable the storeCentrality and storeOccupancu without cent process"); + if ((storeCentrality || storeOccupancyAndIR) && !(doprocessDataWithDCAFitterNCent || doprocessMcWithDCAFitterNCent || doprocessDataWithDCAFitterNMlCent || doprocessMcWithDCAFitterNMlCent || doprocessDataWithDCAFitterNWithUpc || doprocessDataWithDCAFitterNMlWithUpc)) { + LOGP(fatal, "Can't enable the storeCentrality and storeOccupancu without cent process or UPC process"); } auto vbins = (std::vector)binsPt; registry.add("hMass", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{500, 0., 5.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); @@ -284,6 +300,14 @@ struct HfTaskD0 { const AxisSpec thnAxisMinItsNCls{thnConfigAxisMinItsNCls, "Minimum ITS cluster found"}; const AxisSpec thnAxisMinTpcNCrossedRows{thnConfigAxisMinTpcNCrossedRows, "Minimum TPC crossed rows"}; const AxisSpec thnAxisIR{thnConfigAxisIR, "Interaction rate"}; + const AxisSpec thnAxisGapType{thnConfigAxisGapType, "Gap type"}; + const AxisSpec thnAxisFT0A{thnConfigAxisFT0, "FT0-A amplitude"}; + const AxisSpec thnAxisFT0C{thnConfigAxisFT0, "FT0-C amplitude"}; + const AxisSpec thnAxisFV0A{thnConfigAxisFV0A, "FV0-A amplitude"}; + const AxisSpec thnAxisFDDA{thnConfigAxisFDD, "FDD-A amplitude"}; + const AxisSpec thnAxisFDDC{thnConfigAxisFDD, "FDD-C amplitude"}; + const AxisSpec thnAxisZNA{thnConfigAxisZN, "ZNA energy"}; + const AxisSpec thnAxisZNC{thnConfigAxisZN, "ZNC energy"}; if (doprocessMcWithDCAFitterN || doprocessMcWithDCAFitterNCent || doprocessMcWithKFParticle || doprocessMcWithDCAFitterNMl || doprocessMcWithDCAFitterNMlCent || doprocessMcWithKFParticleMl) { std::vector axesAcc = {thnAxisGenPtD, thnAxisGenPtB, thnAxisY, thnAxisOrigin, thnAxisNumPvContr}; @@ -322,10 +346,29 @@ struct HfTaskD0 { const AxisSpec thnAxisNonPromptScore{thnConfigAxisNonPromptScore, "BDT score non-prompt."}; const AxisSpec thnAxisPromptScore{thnConfigAxisPromptScore, "BDT score prompt."}; - axes.insert(axes.begin(), thnAxisPromptScore); - axes.insert(axes.begin(), thnAxisNonPromptScore); - axes.insert(axes.begin(), thnAxisBkgScore); + // Insert ML scores after pt (position 2) to match taskDplus structure: [mass, pt, mlScores, ...] + if (doprocessDataWithDCAFitterNMlWithUpc) { + axes.insert(axes.begin() + 2, thnAxisPromptScore); + axes.insert(axes.begin() + 2, thnAxisNonPromptScore); + axes.insert(axes.begin() + 2, thnAxisBkgScore); + } else { + axes.insert(axes.begin(), thnAxisPromptScore); + axes.insert(axes.begin(), thnAxisNonPromptScore); + axes.insert(axes.begin(), thnAxisBkgScore); + } + } + if (doprocessDataWithDCAFitterNMlWithUpc || doprocessDataWithDCAFitterNWithUpc) { + axes.push_back(thnAxisGapType); + axes.push_back(thnAxisFT0A); + axes.push_back(thnAxisFT0C); + axes.push_back(thnAxisFV0A); + axes.push_back(thnAxisFDDA); + axes.push_back(thnAxisFDDC); + axes.push_back(thnAxisZNA); + axes.push_back(thnAxisZNC); + } + if (applyMl) { registry.add("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type", "Thn for D0 candidates", HistType::kTHnSparseD, axes); registry.get(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"))->Sumw2(); } else { @@ -333,6 +376,12 @@ struct HfTaskD0 { registry.get(HIST("hMassVsPtVsPtBVsYVsOriginVsD0Type"))->Sumw2(); } + registry.add("Data/fitInfo/ampFT0A_vs_ampFT0C", "FT0-A vs FT0-C amplitude;FT0-A amplitude (a.u.);FT0-C amplitude (a.u.)", {HistType::kTH2F, {{2500, 0., 250}, {2500, 0., 250}}}); + registry.add("Data/zdc/energyZNA_vs_energyZNC", "ZNA vs ZNC common energy;E_{ZNA}^{common} (a.u.);E_{ZNC}^{common} (a.u.)", {HistType::kTH2F, {{200, 0., 20}, {200, 0., 20}}}); + registry.add("Data/hUpcGapAfterSelection", "UPC gap type after selection;Gap type;Counts", {HistType::kTH1F, {{7, -1.5, 5.5}}}); + + hfEvSel.addHistograms(registry); + ccdb->setURL(ccdbUrl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -517,6 +566,144 @@ struct HfTaskD0 { } } } + + template + void runAnalysisPerCollisionDataWithUpc(CollType const& collisions, + CandType const& candidates, + BCsType const& bcs, + aod::FT0s const& ft0s, + aod::FV0As const& fv0as, + aod::FDDs const& fdds) + { + for (const auto& collision : collisions) { + float centrality{-1.f}; + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMaskWithUpc(collision, centrality, ccdb, registry, bcs); + if (rejectionMask != 0) { + continue; + } + const auto& bc = collision.template bc_as(); + + // Determine gap type using SGSelector with BC range checking + const auto gapResult = hf_upc::determineGapType(collision, bcs, upcThresholds); + const int gap = gapResult.value; + + // Use the BC with FIT activity if available from SGSelector + auto bcForUPC = bc; + if (gapResult.bc) { + bcForUPC = *(gapResult.bc); + } + + // Get FIT information from the UPC BC + upchelpers::FITInfo fitInfo{}; + udhelpers::getFITinfo(fitInfo, bcForUPC, bcs, ft0s, fv0as, fdds); + + // Get ZDC energies if available (extract once and reuse) + const bool hasZdc = bcForUPC.has_zdc(); + float zdcEnergyZNA = -1.f; + float zdcEnergyZNC = -1.f; + if (hasZdc) { + const auto& zdc = bcForUPC.zdc(); + zdcEnergyZNA = zdc.energyCommonZNA(); + zdcEnergyZNC = zdc.energyCommonZNC(); + } + + // Fill QA histograms using the UPC BC for both FIT and ZDC + if (hasZdc) { + registry.fill(HIST("Data/fitInfo/ampFT0A_vs_ampFT0C"), fitInfo.ampFT0A, fitInfo.ampFT0C); + registry.fill(HIST("Data/zdc/energyZNA_vs_energyZNC"), zdcEnergyZNA, zdcEnergyZNC); + registry.fill(HIST("Data/hUpcGapAfterSelection"), gap); + } + + if (hf_upc::isSingleSidedGap(gap)) { + const auto thisCollId = collision.globalIndex(); + const auto& groupedD0Candidates = candidates.sliceBy(candD0PerCollision, thisCollId); + + // Calculate occupancy and interaction rate if needed + float occ{-1.f}; + float ir{-1.f}; + if (storeOccupancyAndIR && occEstimator != OccupancyEstimator::None) { + occ = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + ir = mRateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), irSource, true) * 1.e-3; // kHz + } + + for (const auto& candidate : groupedD0Candidates) { + if (yCandRecoMax >= 0. && std::abs(HfHelper::yD0(candidate)) > yCandRecoMax) { + continue; + } + + const float massD0 = HfHelper::invMassD0ToPiK(candidate); + const float massD0bar = HfHelper::invMassD0barToKPi(candidate); + const auto ptCandidate = candidate.pt(); + + if (candidate.isSelD0() >= selectionFlagD0) { + registry.fill(HIST("hMass"), massD0, ptCandidate); + registry.fill(HIST("hMassFinerBinning"), massD0, ptCandidate); + registry.fill(HIST("hMassVsPhi"), massD0, ptCandidate, candidate.phi()); + } + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + registry.fill(HIST("hMass"), massD0bar, ptCandidate); + registry.fill(HIST("hMassFinerBinning"), massD0bar, ptCandidate); + registry.fill(HIST("hMassVsPhi"), massD0bar, ptCandidate, candidate.phi()); + } + + // Fill THnSparse with structure matching histogram axes: [mass, pt, (mlScores if FillMl), rapidity, d0Type, (cent if storeCentrality), (occ, ir if storeOccupancyAndIR), gapType, FT0A, FT0C, FV0A, FDDA, FDDC, ZNA, ZNC] + auto fillTHnData = [&](float mass, int d0Type) { + // Pre-calculate vector size to avoid reallocations + constexpr int NAxesBase = 12; // mass, pt, rapidity, d0Type, gapType, FT0A, FT0C, FV0A, FDDA, FDDC, ZNA, ZNC + constexpr int NAxesMl = FillMl ? 3 : 0; // 3 ML scores if FillMl + int const nAxesCent = storeCentrality ? 1 : 0; // centrality if storeCentrality + int const nAxesOccIR = storeOccupancyAndIR ? 2 : 0; // occupancy and IR if storeOccupancyAndIR + int const nAxesTotal = NAxesBase + NAxesMl + nAxesCent + nAxesOccIR; + + std::vector valuesToFill; + valuesToFill.reserve(nAxesTotal); + + // Fill values in order matching histogram axes + valuesToFill.push_back(static_cast(mass)); + valuesToFill.push_back(static_cast(ptCandidate)); + if constexpr (FillMl) { + valuesToFill.push_back(candidate.mlProbD0()[0]); + valuesToFill.push_back(candidate.mlProbD0()[1]); + valuesToFill.push_back(candidate.mlProbD0()[2]); + } + valuesToFill.push_back(static_cast(HfHelper::yD0(candidate))); + valuesToFill.push_back(static_cast(d0Type)); + if (storeCentrality) { + valuesToFill.push_back(centrality); + } + if (storeOccupancyAndIR) { + valuesToFill.push_back(occ); + valuesToFill.push_back(ir); + } + valuesToFill.push_back(static_cast(gap)); + valuesToFill.push_back(static_cast(fitInfo.ampFT0A)); + valuesToFill.push_back(static_cast(fitInfo.ampFT0C)); + valuesToFill.push_back(static_cast(fitInfo.ampFV0A)); + valuesToFill.push_back(static_cast(fitInfo.ampFDDA)); + valuesToFill.push_back(static_cast(fitInfo.ampFDDC)); + valuesToFill.push_back(static_cast(zdcEnergyZNA)); + valuesToFill.push_back(static_cast(zdcEnergyZNC)); + + if constexpr (FillMl) { + registry.get(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"))->Fill(valuesToFill.data()); + } else { + registry.get(HIST("hMassVsPtVsPtBVsYVsOriginVsD0Type"))->Fill(valuesToFill.data()); + } + }; + + if (candidate.isSelD0() >= selectionFlagD0) { + fillTHnData(massD0, SigD0); + fillTHnData(massD0, candidate.isSelD0bar() ? ReflectedD0 : PureSigD0); + } + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + fillTHnData(massD0bar, SigD0bar); + fillTHnData(massD0bar, candidate.isSelD0() ? ReflectedD0bar : PureSigD0bar); + } + } + } + } + } + void processDataWithDCAFitterN(D0Candidates const&, Collisions const& collisions, aod::TracksWExtra const& tracks, aod::BcFullInfos const& bcs) { processData(selectedD0Candidates, collisions, tracks, bcs); @@ -978,6 +1165,32 @@ struct HfTaskD0 { } PROCESS_SWITCH(HfTaskD0, processMcWithKFParticleMl, "Process MC with KFParticle and ML selections", false); // TODO: add the processMcWithKFParticleMlCent + + void processDataWithDCAFitterNWithUpc(soa::Join const& collisions, + aod::BcFullInfos const& bcs, + D0Candidates const&, + aod::TracksWExtra const&, + aod::FT0s const& ft0s, + aod::FV0As const& fv0as, + aod::FDDs const& fdds, + aod::Zdcs const& /*zdcs*/) + { + runAnalysisPerCollisionDataWithUpc(collisions, selectedD0Candidates, bcs, ft0s, fv0as, fdds); + } + PROCESS_SWITCH(HfTaskD0, processDataWithDCAFitterNWithUpc, "Process real data with DCAFitterN w/o ML with UPC", false); + + void processDataWithDCAFitterNMlWithUpc(soa::Join const& collisions, + aod::BcFullInfos const& bcs, + D0CandidatesMl const&, + aod::TracksWExtra const&, + aod::FT0s const& ft0s, + aod::FV0As const& fv0as, + aod::FDDs const& fdds, + aod::Zdcs const& /*zdcs*/) + { + runAnalysisPerCollisionDataWithUpc(collisions, selectedD0CandidatesMl, bcs, ft0s, fv0as, fdds); + } + PROCESS_SWITCH(HfTaskD0, processDataWithDCAFitterNMlWithUpc, "Process real data with DCAFitterN and ML with UPC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/D2H/Tasks/taskDplus.cxx b/PWGHF/D2H/Tasks/taskDplus.cxx index 5948f00d1f1..2e3ceec0884 100644 --- a/PWGHF/D2H/Tasks/taskDplus.cxx +++ b/PWGHF/D2H/Tasks/taskDplus.cxx @@ -16,21 +16,27 @@ /// \author Fabio Catalano , Politecnico and INFN Torino /// \author Vít Kučera , CERN /// \author Luca Aglietta , University and INFN Torino +/// \author Minjung Kim , CERN #include "PWGHF/Core/CentralityEstimation.h" #include "PWGHF/Core/DecayChannels.h" #include "PWGHF/Core/HfHelper.h" #include "PWGHF/Core/SelectorCuts.h" +#include "PWGHF/DataModel/AliasTables.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGHF/DataModel/TrackIndexSkimmingTables.h" #include "PWGHF/Utils/utilsAnalysis.h" #include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGHF/Utils/utilsUpcHf.h" +#include "PWGUD/Core/UPCHelpers.h" +#include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" +#include #include #include #include @@ -50,6 +56,7 @@ #include #include #include +#include #include using namespace o2; @@ -58,6 +65,8 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::hf_centrality; using namespace o2::hf_occupancy; +using namespace o2::hf_evsel; +using namespace o2::analysis::hf_upc; /// D± analysis task struct HfTaskDplus { @@ -71,8 +80,17 @@ struct HfTaskDplus { Configurable> classMl{"classMl", {0, 1, 2}, "Indexes of ML scores to be stored. Three indexes max."}; Configurable storeCentrality{"storeCentrality", false, "Flag to store centrality information"}; Configurable storeOccupancy{"storeOccupancy", false, "Flag to store occupancy information"}; + Configurable storeIR{"storeIR", false, "Flag to store interaction rate information"}; Configurable storePvContributors{"storePvContributors", false, "Flag to store number of PV contributors information"}; Configurable fillMcBkgHistos{"fillMcBkgHistos", false, "Flag to fill and store histograms for MC background"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable irSource{"irSource", "ZNC hadronic", "Estimator of the interaction rate (Recommended: pp --> T0VTX, Pb-Pb --> ZNC hadronic)"}; + + HfEventSelection hfEvSel; // event selection and monitoring + HfUpcGapThresholds upcThresholds; // UPC gap determination thresholds + ctpRateFetcher mRateFetcher; // interaction rate fetcher + + Service ccdb; using CandDplusData = soa::Filtered>; using CandDplusDataWithMl = soa::Filtered>; @@ -85,6 +103,7 @@ struct HfTaskDplus { Filter filterDplusFlag = (o2::aod::hf_track_index::hfflag & static_cast(BIT(aod::hf_cand_3prong::DecayType::DplusToPiKPi))) != static_cast(0); + Preslice candDplusPerCollision = aod::hf_cand::collisionId; Preslice mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId; PresliceUnsorted recoColPerMcCollision = aod::mccollisionlabel::mcCollisionId; @@ -103,12 +122,18 @@ struct HfTaskDplus { ConfigurableAxis thnConfigAxisY{"thnConfigAxisY", {40, -1, 1}, "Cand. rapidity bins"}; ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {110, 0., 110.}, "axis for centrality"}; ConfigurableAxis thnConfigAxisOccupancy{"thnConfigAxisOccupancy", {14, 0, 14000}, "axis for occupancy"}; + ConfigurableAxis thnConfigAxisIR{"thnConfigAxisIR", {5000, 0, 500}, "Interaction rate (kHz)"}; ConfigurableAxis thnConfigAxisPvContributors{"thnConfigAxisPvContributors", {100, 0., 100.}, "axis for PV contributors"}; ConfigurableAxis thnConfigAxisPtBHad{"thnConfigAxisPtBHad", {25, 0., 50}, "axis for pt of B hadron decayed into D candidate"}; ConfigurableAxis thnConfigAxisFlagBHad{"thnConfigAxisFlagBHad", {5, 0., 5}, "axis for PDG of B hadron"}; ConfigurableAxis thnConfigAxisMlScore0{"thnConfigAxisMlScore0", {100, 0., 1.}, "axis for ML output score 0"}; ConfigurableAxis thnConfigAxisMlScore1{"thnConfigAxisMlScore1", {100, 0., 1.}, "axis for ML output score 1"}; ConfigurableAxis thnConfigAxisMlScore2{"thnConfigAxisMlScore2", {100, 0., 1.}, "axis for ML output score 2"}; + ConfigurableAxis thnConfigAxisGapType{"thnConfigAxisGapType", {7, -1.5, 5.5}, "axis for UPC gap type (see TrueGap enum in o2::aod::sgselector)"}; + ConfigurableAxis thnConfigAxisFT0{"thnConfigAxisFT0", {1001, -1.5, 999.5}, "axis for FT0 amplitude (a.u.)"}; + ConfigurableAxis thnConfigAxisFV0A{"thnConfigAxisFV0A", {2001, -1.5, 1999.5}, "axis for FV0-A amplitude (a.u.)"}; + ConfigurableAxis thnConfigAxisFDD{"thnConfigAxisFDD", {200, 0., 4000.}, "axis for FDD amplitude (a.u.)"}; + ConfigurableAxis thnConfigAxisZN{"thnConfigAxisZN", {510, -1.5, 49.5}, "axis for ZN energy (a.u.)"}; HistogramRegistry registry{ "registry", @@ -124,7 +149,7 @@ struct HfTaskDplus { void init(InitContext&) { - std::array doprocess{doprocessData, doprocessDataWithMl, doprocessMc, doprocessMcWithMl}; + std::array doprocess{doprocessData, doprocessDataWithMl, doprocessMc, doprocessMcWithMl, doprocessDataWithUpc, doprocessDataWithMlWithUpc}; if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) { LOGP(fatal, "Only one process function should be enabled! Please check your configuration!"); } @@ -139,7 +164,16 @@ struct HfTaskDplus { AxisSpec const thnAxisFlagBHad{thnConfigAxisFlagBHad, "B Hadron flag"}; AxisSpec const thnAxisCent{thnConfigAxisCent, "Centrality"}; AxisSpec const thnAxisOccupancy{thnConfigAxisOccupancy, "Occupancy"}; + AxisSpec const thnAxisIR{thnConfigAxisIR, "Interaction rate (kHz)"}; AxisSpec const thnAxisPvContributors{thnConfigAxisPvContributors, "PV contributors"}; + AxisSpec const thnAxisGapType{thnConfigAxisGapType, "Gap type"}; + AxisSpec const thnAxisFT0A{thnConfigAxisFT0, "FT0-A amplitude"}; + AxisSpec const thnAxisFT0C{thnConfigAxisFT0, "FT0-C amplitude"}; + AxisSpec const thnAxisFV0A{thnConfigAxisFV0A, "FV0-A amplitude"}; + AxisSpec const thnAxisFDDA{thnConfigAxisFDD, "FDD-A amplitude"}; + AxisSpec const thnAxisFDDC{thnConfigAxisFDD, "FDD-C amplitude"}; + AxisSpec const thnAxisZNA{thnConfigAxisZN, "ZNA energy"}; + AxisSpec const thnAxisZNC{thnConfigAxisZN, "ZNC energy"}; registry.add("hMass", "3-prong candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{350, 1.7, 2.05}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); registry.add("hEta", "3-prong candidates;candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); @@ -179,10 +213,10 @@ struct HfTaskDplus { registry.add("hPtVsYGenPrompt", "MC particles (matched, prompt);#it{p}_{T}^{gen.}; #it{y}", {HistType::kTH2F, {{vbins, "#it{p}_{T} (GeV/#it{c})"}, {100, -5., 5.}}}); registry.add("hPtVsYGenNonPrompt", "MC particles (matched, non-prompt);#it{p}_{T}^{gen.}; #it{y}", {HistType::kTH2F, {{vbins, "#it{p}_{T} (GeV/#it{c})"}, {100, -5., 5.}}}); - if (doprocessDataWithMl || doprocessData) { + if (doprocessDataWithMl || doprocessData || doprocessDataWithMlWithUpc) { std::vector axes = {thnAxisMass, thnAxisPt}; - if (doprocessDataWithMl) { + if (doprocessDataWithMl || doprocessDataWithMlWithUpc) { axes.push_back(thnAxisMlScore0); axes.push_back(thnAxisMlScore1); axes.push_back(thnAxisMlScore2); @@ -193,6 +227,19 @@ struct HfTaskDplus { if (storeOccupancy) { axes.push_back(thnAxisOccupancy); } + if (storeIR) { + axes.push_back(thnAxisIR); + } + if (doprocessDataWithMlWithUpc) { + axes.push_back(thnAxisGapType); + axes.push_back(thnAxisFT0A); + axes.push_back(thnAxisFT0C); + axes.push_back(thnAxisFV0A); + axes.push_back(thnAxisFDDA); + axes.push_back(thnAxisFDDC); + axes.push_back(thnAxisZNA); + axes.push_back(thnAxisZNC); + } registry.add("hSparseMass", "THn for Dplus", HistType::kTHnSparseF, axes); } @@ -239,6 +286,16 @@ struct HfTaskDplus { registry.add("hSparseMassGenPrompt", "THn for gen Prompt Dplus", HistType::kTHnSparseF, axesGenPrompt); registry.add("hSparseMassGenFD", "THn for gen FD Dplus", HistType::kTHnSparseF, axesGenFD); } + + registry.add("Data/fitInfo/ampFT0A_vs_ampFT0C", "FT0-A vs FT0-C amplitude;FT0-A amplitude (a.u.);FT0-C amplitude (a.u.)", {HistType::kTH2F, {{2500, 0., 250}, {2500, 0., 250}}}); + registry.add("Data/zdc/energyZNA_vs_energyZNC", "ZNA vs ZNC common energy;E_{ZNA}^{common} (a.u.);E_{ZNC}^{common} (a.u.)", {HistType::kTH2F, {{200, 0., 20}, {200, 0., 20}}}); + registry.add("Data/hUpcGapAfterSelection", "UPC gap type after selection;Gap type;Counts", {HistType::kTH1F, {{7, -1.5, 5.5}}}); + + hfEvSel.addHistograms(registry); // collision monitoring + + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); } // Fill histograms of quantities for the reconstructed Dplus candidates @@ -646,6 +703,124 @@ struct HfTaskDplus { } } + template + void runAnalysisPerCollisionDataWithUpc(CollType const& collisions, + CandType const& candidates, + BCsType const& bcs, + aod::FT0s const& ft0s, + aod::FV0As const& fv0as, + aod::FDDs const& fdds) + { + for (const auto& collision : collisions) { + float centrality{-1.f}; + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMaskWithUpc(collision, centrality, ccdb, registry, bcs); + if (rejectionMask != 0) { + /// at least one event selection not satisfied --> reject the candidate + continue; + } + const auto& bc = collision.template bc_as(); + + // Determine gap type using SGSelector with BC range checking + const auto gapResult = hf_upc::determineGapType(collision, bcs, upcThresholds); + const int gap = gapResult.value; + + // Use the BC with FIT activity if available from SGSelector + auto bcForUPC = bc; + if (gapResult.bc) { + bcForUPC = *(gapResult.bc); + } + + // Get FIT information from the UPC BC + upchelpers::FITInfo fitInfo{}; + udhelpers::getFITinfo(fitInfo, bcForUPC, bcs, ft0s, fv0as, fdds); + + // Get ZDC energies if available (extract once and reuse) + const bool hasZdc = bcForUPC.has_zdc(); + float zdcEnergyZNA = -1.f; + float zdcEnergyZNC = -1.f; + if (hasZdc) { + const auto& zdc = bcForUPC.zdc(); + zdcEnergyZNA = zdc.energyCommonZNA(); + zdcEnergyZNC = zdc.energyCommonZNC(); + } + + // Fill QA histograms using the UPC BC for both FIT and ZDC + if (hasZdc) { + registry.fill(HIST("Data/fitInfo/ampFT0A_vs_ampFT0C"), fitInfo.ampFT0A, fitInfo.ampFT0C); + registry.fill(HIST("Data/zdc/energyZNA_vs_energyZNC"), zdcEnergyZNA, zdcEnergyZNC); + registry.fill(HIST("Data/hUpcGapAfterSelection"), gap); + } + + if (hf_upc::isSingleSidedGap(gap)) { + const auto thisCollId = collision.globalIndex(); + const auto& groupedDplusCandidates = candidates.sliceBy(candDplusPerCollision, thisCollId); + + float cent{-1.f}; + float occ{-1.f}; + float ir{-1.f}; + if (storeOccupancy && occEstimator != OccupancyEstimator::None) { + occ = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + } + if (storeIR) { + ir = mRateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), irSource, true) * 1.e-3; // kHz + } + static constexpr auto HSparseMass = HIST("hSparseMass"); + // Lambda function to fill THn - handles both ML and non-ML cases + auto fillTHnData = [&](const auto& candidate) { + // Pre-calculate vector size to avoid reallocations + constexpr int NAxesBase = 10; // mass, pt, gapType, FT0A, FT0C, FV0A, FDDA, FDDC, ZNA, ZNC + constexpr int NAxesMl = FillMl ? 3 : 0; // 3 ML scores if FillMl + int const nAxesCent = storeCentrality ? 1 : 0; // centrality if storeCentrality + int const nAxesOcc = storeOccupancy ? 1 : 0; // occupancy if storeOccupancy + int const nAxesIR = storeIR ? 1 : 0; // IR if storeIR + int const nAxesTotal = NAxesBase + NAxesMl + nAxesCent + nAxesOcc + nAxesIR; + + std::vector valuesToFill; + valuesToFill.reserve(nAxesTotal); + + // Fill values in order matching histogram axes + valuesToFill.push_back(HfHelper::invMassDplusToPiKPi(candidate)); + valuesToFill.push_back(static_cast(candidate.pt())); + if constexpr (FillMl) { + std::vector outputMl = {-999., -999., -999.}; + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMl[iclass] = candidate.mlProbDplusToPiKPi()[classMl->at(iclass)]; + } + valuesToFill.push_back(outputMl[0]); + valuesToFill.push_back(outputMl[1]); + valuesToFill.push_back(outputMl[2]); + } + if (storeCentrality) { + valuesToFill.push_back(cent); + } + if (storeOccupancy) { + valuesToFill.push_back(occ); + } + if (storeIR) { + valuesToFill.push_back(ir); + } + valuesToFill.push_back(static_cast(gap)); + valuesToFill.push_back(static_cast(fitInfo.ampFT0A)); + valuesToFill.push_back(static_cast(fitInfo.ampFT0C)); + valuesToFill.push_back(static_cast(fitInfo.ampFV0A)); + valuesToFill.push_back(static_cast(fitInfo.ampFDDA)); + valuesToFill.push_back(static_cast(fitInfo.ampFDDC)); + valuesToFill.push_back(static_cast(zdcEnergyZNA)); + valuesToFill.push_back(static_cast(zdcEnergyZNC)); + registry.get(HSparseMass)->Fill(valuesToFill.data()); + }; + + for (const auto& candidate : groupedDplusCandidates) { + if ((yCandRecoMax >= 0. && std::abs(HfHelper::yDplus(candidate)) > yCandRecoMax)) { + continue; + } + fillHisto(candidate); + fillTHnData(candidate); + } + } + } + } + // process functions void processData(CandDplusData const& candidates, CollisionsCent const& collisions) { @@ -678,6 +853,32 @@ struct HfTaskDplus { runAnalysisMcGen(mcGenCollisions, mcRecoCollisions, mcGenParticles); } PROCESS_SWITCH(HfTaskDplus, processMcWithMl, "Process MC with ML", false); + + void processDataWithUpc(soa::Join const& collisions, + aod::BcFullInfos const& bcs, + CandDplusData const& selectedDplusCandidates, + aod::Tracks const&, + aod::FT0s const& ft0s, + aod::FV0As const& fv0as, + aod::FDDs const& fdds, + aod::Zdcs const& /*zdcs*/) + { + runAnalysisPerCollisionDataWithUpc(collisions, selectedDplusCandidates, bcs, ft0s, fv0as, fdds); + } + PROCESS_SWITCH(HfTaskDplus, processDataWithUpc, "Process real data w/o ML with UPC", false); + + void processDataWithMlWithUpc(soa::Join const& collisions, + aod::BcFullInfos const& bcs, + CandDplusDataWithMl const& selectedDplusCandidatesMl, + aod::Tracks const&, + aod::FT0s const& ft0s, + aod::FV0As const& fv0as, + aod::FDDs const& fdds, + aod::Zdcs const& /*zdcs*/) + { + runAnalysisPerCollisionDataWithUpc(collisions, selectedDplusCandidatesMl, bcs, ft0s, fv0as, fdds); + } + PROCESS_SWITCH(HfTaskDplus, processDataWithMlWithUpc, "Process real data with the ML method with UPC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/Utils/utilsUpcHf.h b/PWGHF/Utils/utilsUpcHf.h new file mode 100644 index 00000000000..45190185c74 --- /dev/null +++ b/PWGHF/Utils/utilsUpcHf.h @@ -0,0 +1,150 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file utilsUpcHf.h +/// \brief Utility functions for Ultra-Peripheral Collision (UPC) analysis in Heavy Flavor physics +/// +/// \author Minjung Kim , CERN + +#ifndef PWGHF_UTILS_UTILSUPCHF_H_ +#define PWGHF_UTILS_UTILSUPCHF_H_ + +#include "PWGUD/Core/SGCutParHolder.h" +#include "PWGUD/Core/SGSelector.h" +#include "PWGUD/Core/UDHelpers.h" + +#include + +#include + +namespace o2::analysis::hf_upc +{ + +/// \brief Use TrueGap enum from SGSelector for gap type classification +using o2::aod::sgselector::TrueGap; + +/// \brief Configurable group for UPC gap determination thresholds +struct HfUpcGapThresholds : o2::framework::ConfigurableGroup { + std::string prefix = "upc"; // JSON group name + o2::framework::Configurable fv0aThreshold{"fv0aThreshold", 100.0f, "FV0-A amplitude threshold for UPC gap determination (a.u.)"}; + o2::framework::Configurable ft0aThreshold{"ft0aThreshold", 100.0f, "FT0-A amplitude threshold for UPC gap determination (a.u.)"}; + o2::framework::Configurable ft0cThreshold{"ft0cThreshold", 50.0f, "FT0-C amplitude threshold for UPC gap determination (a.u.)"}; + o2::framework::Configurable zdcThreshold{"zdcThreshold", 1.0f, "ZDC energy threshold for UPC gap determination (a.u.)"}; +}; + +/// \brief Default thresholds for gap determination +namespace defaults +{ +constexpr float AmplitudeThresholdFV0A = 100.0f; ///< Amplitude threshold for FV0-A (a.u.) +constexpr float AmplitudeThresholdFT0A = 100.0f; ///< Amplitude threshold for FT0-A (a.u.) +constexpr float AmplitudeThresholdFT0C = 50.0f; ///< Amplitude threshold for FT0-C (a.u.) +constexpr float MaxFITTime = 4.0f; ///< Maximum FIT time (ns) +constexpr int NDtColl = 1000; ///< Time window for BC range (ns) +constexpr int MinNBCs = 7; ///< Minimum number of BCs to check +constexpr int MinNTracks = 0; ///< Minimum number of tracks +constexpr int MaxNTracks = 100; ///< Maximum number of tracks +} // namespace defaults + +/// \brief Determine gap type using SGSelector with BC range checking +/// \tparam TCollision Collision type +/// \tparam TBCs BC table type +/// \param collision Collision object +/// \param bcs BC table +/// \param amplitudeThresholdFV0A Threshold for FV0-A (default: 100.0) +/// \param amplitudeThresholdFT0A Threshold for FT0-A (default: 100.0) +/// \param amplitudeThresholdFT0C Threshold for FT0-C (default: 50.0) +/// \return SelectionResult with gap type value and BC pointer +template +inline auto determineGapType(TCollision const& collision, + TBCs const& bcs, + float amplitudeThresholdFV0A = defaults::AmplitudeThresholdFV0A, + float amplitudeThresholdFT0A = defaults::AmplitudeThresholdFT0A, + float amplitudeThresholdFT0C = defaults::AmplitudeThresholdFT0C) +{ + using BCType = std::decay_t())>; + + // Configure SGSelector thresholds + SGCutParHolder sgCuts; + sgCuts.SetNDtcoll(defaults::NDtColl); + sgCuts.SetMinNBCs(defaults::MinNBCs); + sgCuts.SetNTracks(defaults::MinNTracks, defaults::MaxNTracks); + sgCuts.SetMaxFITtime(defaults::MaxFITTime); + sgCuts.SetFITAmpLimits({amplitudeThresholdFV0A, amplitudeThresholdFT0A, amplitudeThresholdFT0C}); + + // Get BC and BC range + if (!collision.has_foundBC()) { + return SelectionResult{TrueGap::NoGap, nullptr}; + } + + const auto bc = collision.template foundBC_as(); + const auto bcRange = udhelpers::compatibleBCs(collision, sgCuts.NDtcoll(), bcs, sgCuts.minNBCs()); + + // Create SGSelector instance and determine gap type with BC range checking + SGSelector sgSelector; + const auto sgResult = sgSelector.IsSelected(sgCuts, collision, bcRange, bc); + + return sgResult; +} + +/// \brief Determine gap type using SGSelector with BC range checking and HfUpcGapThresholds +/// \tparam TCollision Collision type +/// \tparam TBCs BC table type +/// \param collision Collision object +/// \param bcs BC table +/// \param thresholds HfUpcGapThresholds object containing all UPC thresholds +/// \return SelectionResult with gap type value and BC pointer +template +inline auto determineGapType(TCollision const& collision, + TBCs const& bcs, + HfUpcGapThresholds const& thresholds) +{ + return determineGapType(collision, bcs, + thresholds.fv0aThreshold.value, + thresholds.ft0aThreshold.value, + thresholds.ft0cThreshold.value); +} + +/// \brief Check if the gap type is a single-sided gap (SingleGapA or SingleGapC) +/// \param gap TrueGap enum value +/// \return true if single-sided gap, false otherwise +constexpr bool isSingleSidedGap(int gap) noexcept +{ + return (gap == TrueGap::SingleGapA || gap == TrueGap::SingleGapC); +} + +/// \brief Get gap type name as string +/// \param gap TrueGap enum value +/// \return String representation of gap type +constexpr const char* getGapTypeName(int gap) noexcept +{ + switch (gap) { + case TrueGap::NoGap: + return "NoGap"; + case TrueGap::SingleGapA: + return "SingleGapA"; + case TrueGap::SingleGapC: + return "SingleGapC"; + case TrueGap::DoubleGap: + return "DoubleGap"; + case TrueGap::NoUpc: + return "NoUpc"; + case TrueGap::TrkOutOfRange: + return "TrkOutOfRange"; + case TrueGap::BadDoubleGap: + return "BadDoubleGap"; + default: + return "Unknown"; + } +} + +} // namespace o2::analysis::hf_upc + +#endif // PWGHF_UTILS_UTILSUPCHF_H_