From 1df39f785300832701414ad09f0e55caa56ec1c3 Mon Sep 17 00:00:00 2001 From: Yash Patley <52608802+yashpatley@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:31:52 +0530 Subject: [PATCH 1/6] [PWGCF] Added chargeBalanceFunction.cxx Task to calculate R2, P2 and BF of charged hadrons --- .../Tasks/chargeBalanceFunction.cxx | 783 ++++++++++++++++++ 1 file changed, 783 insertions(+) create mode 100644 PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx diff --git a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx new file mode 100644 index 00000000000..77d959447ab --- /dev/null +++ b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx @@ -0,0 +1,783 @@ +// 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 lambdaR2Correlation.cxx +/// \brief R2 P2 and BF of charged hadrons. +/// \author Yash Patley + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::physics; +using namespace o2::constants::math; + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::physics; +using namespace o2::constants::math; + +enum ParticlePairType { + kPM = 0, + kPP, + kMM +}; + +enum RecGenType { + kRec = 0, + kGen +}; + +enum DMCType { + kData = 0, + kMC +}; + +struct ChargeBalanceFunction { + // Collisions + Configurable cZVtxCut{"cZVtxCut", 10.0, "Min VtxZ cut"}; + Configurable cMinCent{"cMinCent", 0., "Minumum Centrality"}; + Configurable cMaxCent{"cMaxCent", 100.0, "Maximum Centrality"}; + Configurable cSel8Trig{"cSel8Trig", true, "Sel8 (T0A + T0C) Selection Run3"}; + Configurable cPileupReject{"cPileupReject", true, "Pileup rejection"}; + Configurable cZVtxTimeDiff{"cZVtxTimeDiff", true, "z-vtx time diff selection"}; + Configurable cIsGoodITSLayers{"cIsGoodITSLayers", true, "Good ITS Layers All"}; + Configurable cMinOccupancy{"cMinOccupancy", 0, "Minimum FT0C Occupancy"}; + Configurable cMaxOccupancy{"cMaxOccupancy", 1e6, "Maximum FT0C Occupancy"}; + + // Tracks + Configurable cTrackNPtBins{"cTrackNPtBins", 20, "N pT bins"}; + Configurable cTrackMinPt{"cTrackMinPt", 0.2, "p_{T} minimum"}; + Configurable cTrackMaxPt{"cTrackMaxPt", 2.0, "p_{T} maximum"}; + Configurable cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"}; + Configurable cTrackGlobal{"cTrackGlobal", true, "Global Track"}; + Configurable cTrackDcaXYCut{"cTrackDcaXYCut", 0.1, "DcaXY Cut"}; + Configurable cTrackDcaZCut{"cTrackDcaZCut", 1., "DcaZ Cut"}; + Configurable cTpcElRejCutMin{"cTpcElRejCutMin", -3., "Electron Rejection Cut Minimum"}; + Configurable cTpcElRejCutMax{"cTpcElRejCutMax", 5., "Electron Rejection Cut Maximum"}; + Configurable cTpcRejCut{"cTpcRejCut", 3, "TPC Rej Cut"}; + + // Correlation + Configurable cNEtaBins{"cNEtaBins", 16, "N Eta Bins"}; + Configurable cNPhiBins{"cNPhiBins", 36, "N Phi Bins"}; + + // Femtoscopic correction + Configurable cApplyFemtoSel{"cApplyFemtoSel", false, "Femto qinv selection"}; + Configurable cFemtoCut{"cFemtoCut", 0.1, "Kaon--Lambda Femto qinv cut"}; + + // Two-track cuts + Configurable cApplyTwoTrackCut{"cApplyTwoTrackCut", false, "Flag for two track cut"}; + Configurable cDEtaCut{"cDEtaCut", 0.02, "DEta cut"}; + Configurable cDPhiStarCut{"cDPhiStarCut", 0.02, "DPhiStar cut"}; + + // Efficiency Correction + Configurable cGetCorrectionFlag{"cGetCorrectionFlag", false, "Apply correction flag"}; + Configurable cGetNuaCorrectionFlag{"cGetNuaCorrectionFlag", false, "Apply NUA correction flag"}; + + // CCDB + Configurable cUrlCCDB{"cUrlCCDB", "http://alice-ccdb.cern.ch", "ALICE CCDB URL"}; + Configurable cPathCCDBRecoEff{"cPathCCDBRecoEff", "Users/y/ypatley/CBF/Test/RecoEfficiency", "Path for ccdb-object for reco efficiency"}; + Configurable cPathCCDBNuaCorr{"cPathCCDBNuaCorr", "Users/y/ypatley/CBF/Test/Nua", "Path for ccdb-object for NUA correction"}; + Configurable nolaterthan{"nolaterthan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; + + // Configurable Axis + ConfigurableAxis cCentBins{"cCentBins", {VARIABLE_WIDTH, 0.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f}, "Variable Centrality Bins"}; + ConfigurableAxis cPosZBins{"cPosZBins", {VARIABLE_WIDTH, -10.f, -8.0f, -6.0f, -4.0f, -2.0f, 0.f, 2.0f, 4.0f, 6.0f, 8.0f, 10.f}, "Variable Vz Bins"}; + + // Ccdb service + Service ccdb{}; + o2::parameters::GRPMagField* grpo = nullptr; + + // Pdg service + Service pdg{}; + + // Histogram Registry. + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // Initialize global variables + float nrapbins = 0.; + float kminrap = 0.; + float kmaxrap = 0.; + float nphibins = 0.; + float kminphi = 0.; + float kmaxphi = TwoPI; + float rapbinwidth = 0.; + float phibinwidth = 0.; + float q = 0., e = 0., qinv = 0.; + float posz = 0., cent = 0., mult = 0.; + float magField = 0.; + std::array vTpcRadii = {0.8, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24}; + + // Efficiency object + TList* ccdbObjRecoEff = nullptr; + TList* ccdbObjNuaCorr = nullptr; + struct CorrHist { + THnSparseF *hRecEffP, *hRecEffM; + TH3F *hNuaP, *hNuaM; + } corrHist; + int nHistDim = 0; + + void init(InitContext const&) + { + // Set CCDB url + ccdb->setURL(cUrlCCDB.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(nolaterthan.value); + + // Set Density Histogram Attributes + nrapbins = static_cast(cNEtaBins); + kminrap = static_cast(-cTrackEtaCut); + kmaxrap = static_cast(cTrackEtaCut); + nphibins = static_cast(cNPhiBins); + + rapbinwidth = (kmaxrap - kminrap) / nrapbins; + phibinwidth = (kmaxphi - kminphi) / nphibins; + + int knrapphibins = static_cast(cNEtaBins) * static_cast(cNPhiBins); + float kminrapphi = 0.; + float kmaxrapphi = knrapphibins; + + // Initialize axis specifications + const AxisSpec axisCheck(1, 0, 1, ""); + const AxisSpec axisChMult(200, 0, 200, "N_{ch}"); + const AxisSpec axisVz(220, -11, 11, "V_{z} (cm)"); + + const AxisSpec axisCent(cCentBins, "FT0C (%)"); + const AxisSpec axisPosZ(cPosZBins, "V_{z} (cm)"); + + const AxisSpec axisDEta(320, -1.6, 1.6, "#Delta#eta"); + const AxisSpec axisDPhi(640, -PIHalf, 3. * PIHalf, "#Delta#varphi"); + + const AxisSpec axisTrackDcaXY{60, -0.15, 0.15, "DCA_{XY}"}; + const AxisSpec axisTrackDcaZ{230, -1.15, 1.15, "DCA_{Z}"}; + const AxisSpec axisdEdx(360, 20, 200, "#frac{dE}{dx}"); + + const AxisSpec axisTrackPt(cTrackNPtBins, cTrackMinPt, cTrackMaxPt, "p_{T} (GeV/#it{c})"); + const AxisSpec axisTrackEta(cNEtaBins, -cTrackEtaCut, cTrackEtaCut, "#eta"); + const AxisSpec axisTrackPhi(cNPhiBins, 0., TwoPI, "#varphi (rad)"); + + const AxisSpec axisRapPhi(knrapphibins, kminrapphi, kmaxrapphi, "#eta#varphi"); + + // Create Histograms + // Event histograms + histos.add("Event/h1f_collision_cent", "FT0C(%)", kTH1F, {axisCent}); + histos.add("Event/h1f_collision_Vz", "V_{z}-distribution", kTH1F, {axisVz}); + + // Track QA + histos.add("TrackQA/hPtDcaXY", "DCA_{XY} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaXY}); + histos.add("TrackQA/hPtDcaZ", "DCA_{Z} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaZ}); + histos.add("TrackQA/hTrackTPCdEdX", "hTrackTPCdEdX", kTH2F, {axisTrackPt, axisdEdx}); + + // Two track cut + histos.add("QA/TwoTrackCut/Before/h2d_n2_detadphi", "#rho_{2}", kTH2D, {axisDEta, axisDPhi}); + histos.add("QA/TwoTrackCut/After/h2d_n2_detadphi", "#rho_{2}", kTH2D, {axisDEta, axisDPhi}); + histos.add("QA/FemtoCut/Before/h2d_n2_detadphi", "#rho_{2}", kTH2D, {axisDEta, axisDPhi}); + histos.add("QA/FemtoCut/After/h2d_n2_detadphi", "#rho_{2}", kTH2D, {axisDEta, axisDPhi}); + + // Efficiency Histograms + // Single Particle Efficiencies + histos.add("Reco/Efficiency/h1f_n1_pt_P", "#rho_{1}^{#plus}", kTH1F, {axisTrackPt}); + histos.add("Reco/Efficiency/h1f_n1_pt_M", "#rho_{1}^{#minus}", kTH1F, {axisTrackPt}); + histos.add("Reco/Efficiency/h4f_n1_centvzptrap_P", "#rho_{1}^{#plus}", kTHnSparseF, {axisCent, axisPosZ, axisTrackPt, axisTrackEta}); + histos.add("Reco/Efficiency/h4f_n1_centvzptrap_M", "#rho_{1}^{#minus}", kTHnSparseF, {axisCent, axisPosZ, axisTrackPt, axisTrackEta}); + + // NUA phi + histos.add("Reco/NUA/h3f_n1_vzrapphi_P", "#rho_{1}^{#plus}", kTH3F, {axisPosZ, axisTrackEta, axisTrackPhi}); + histos.add("Reco/NUA/h3f_n1_vzrapphi_M", "#rho_{1}^{#minus}", kTH3F, {axisPosZ, axisTrackEta, axisTrackPhi}); + + // Correction checks + histos.add("Reco/h1f_n1_pt_P", "#rho_{1}^{#plus}", kTH1F, {axisTrackPt}); + histos.add("Reco/h1f_n1_pt_M", "#rho_{1}^{#minus}", kTH1F, {axisTrackPt}); + histos.add("Reco/h1f_n1_rap_P", "#rho_{1}^{#plus}", kTH1F, {axisTrackEta}); + histos.add("Reco/h1f_n1_rap_M", "#rho_{1}^{#minus}", kTH1F, {axisTrackEta}); + histos.add("Reco/h1f_n1_phi_P", "#rho_{1}^{#plus}", kTH1F, {axisTrackPhi}); + histos.add("Reco/h1f_n1_phi_M", "#rho_{1}^{#minus}", kTH1F, {axisTrackPhi}); + + // Single and Two Particle Densities + // Rho1 for R2 RapPhi + histos.add("Reco/h3f_n1_rapphi_P", "#rho_{1}^{#plus}", kTH3F, {axisCent, axisTrackEta, axisTrackPhi}); + histos.add("Reco/h3f_n1_rapphi_M", "#rho_{1}^{#minus}", kTH3F, {axisCent, axisTrackEta, axisTrackPhi}); + + // Rho1 for P2 RapPhi + histos.add("Reco/h3f_pt_rapphi_P", "#rho_{1}^{#plus}", kTH3F, {axisCent, axisTrackEta, axisTrackPhi}); + histos.add("Reco/h3f_pt_rapphi_M", "#rho_{1}^{#minus}", kTH3F, {axisCent, axisTrackEta, axisTrackPhi}); + + // Rho2 for R2 RapPhi + histos.add("Reco/h3f_n2_rapphi_PM", "#rho_{2}^{#plus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_n2_rapphi_PP", "#rho_{2}^{#plus#plus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_n2_rapphi_MM", "#rho_{2}^{#minus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + + // Rho2 for P2 RapPhi + histos.add("Reco/h3f_ptpt_rapphi_PM", "#rho_{2}^{#plus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_ptpt_rapphi_PP", "#rho_{2}^{#plus#plus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_ptpt_rapphi_MM", "#rho_{2}^{#minus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_npt_rapphi_PM", "#rho_{2}^{#plus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_npt_rapphi_PP", "#rho_{2}^{#plus#plus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_npt_rapphi_MM", "#rho_{2}^{#minus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_ptn_rapphi_PM", "#rho_{2}^{#plus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_ptn_rapphi_PP", "#rho_{2}^{#plus#plus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + histos.add("Reco/h3f_ptn_rapphi_MM", "#rho_{2}^{#minus#minus}", kTH3F, {axisCent, axisRapPhi, axisRapPhi}); + + histos.addClone("Reco/", "McGen/"); + + // MC Generated Histograms + if (doprocessMCRecoGen) { + // McGen Histos + histos.add("McGen/h1f_collision_recgen", "# of Reco Collision Associated to One Mc Generator Collision", kTH1F, {axisChMult}); + histos.add("McGen/h2f_collision_posZ", "V_{z}-distribution", kTH2F, {axisVz, axisVz}); + histos.add("McGen/h2f_collision_cent", "FT0M Centrality", kTH2F, {axisCent, axisCent}); + } + + // Load correction factor + if (cGetCorrectionFlag) { + // Set CCDB url + ccdb->setURL(cUrlCCDB.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(nolaterthan.value); + + // Get CCDB object + ccdbObjRecoEff = ccdb->getForTimeStamp(cPathCCDBRecoEff.value, nolaterthan.value); + ccdbObjNuaCorr = ccdb->getForTimeStamp(cPathCCDBNuaCorr.value, nolaterthan.value); + + // Load reco eff corrections + LoadRecoEfficiencyHistograms(); + } + } + + // Load reco efficiency histograms + void LoadRecoEfficiencyHistograms() + { + // Efficiency correction histograms + corrHist.hRecEffP = dynamic_cast(ccdbObjRecoEff->FindObject("h_RecEff_P")); + corrHist.hRecEffM = dynamic_cast(ccdbObjRecoEff->FindObject("h_RecEff_M")); + if (!corrHist.hRecEffP || !corrHist.hRecEffM) { + LOGF(fatal, "CCDB efficiency object doesn't exist !"); + } + nHistDim = corrHist.hRecEffP->GetNdimensions(); + LOGF(info, "Efficiency correction histogram dimensions: %d", nHistDim); + + // Nua correction histograms + corrHist.hNuaP = dynamic_cast(ccdbObjNuaCorr->FindObject("h_Nua_P")); + corrHist.hNuaM = dynamic_cast(ccdbObjNuaCorr->FindObject("h_Nua_M")); + if (!corrHist.hNuaP || !corrHist.hNuaM) { + LOGF(fatal, "CCDB NUA object doesn't exist !"); + } + } + + // Get magnetic field + float getMagneticField(int64_t const& timestamp) + { + if (grpo == nullptr) { + grpo = ccdb->getForTimeStamp("GLO/Config/GRPMagField", timestamp); + if (grpo == nullptr) { + LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); + return 1; + } + } + auto field = std::lround(5.f * grpo->getL3Current() / 30000.f); + return 0.1 * field; + } + + template + bool selCollision(C const& col) + { + posz = col.posZ(); + if (std::abs(posz) >= cZVtxCut) { // VtxZ selection + return false; + } + + if (cSel8Trig && !col.sel8()) { // Sel8 selection + return false; + } + + cent = col.centFT0C(); + if (cent <= cMinCent || cent >= cMaxCent) { // Centrality selection + return false; + } + + if (col.ft0cOccupancyInTimeRange() < cMinOccupancy || col.ft0cOccupancyInTimeRange() > cMaxOccupancy) { // Occupancy cut + return false; + } + + if (cPileupReject && !col.selection_bit(aod::evsel::kNoSameBunchPileup)) { // Pile-up rejection + return false; + } + + if (cZVtxTimeDiff && !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { // ZvtxFT0 vs PV + return false; + } + + if (cIsGoodITSLayers && !col.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { // All ITS layer active + return false; + } + + // // Set Multiplicity + // mult = col.multTPC(); + + return true; + } + + // Track Selection + template + bool selectTrack(T const& track) + { + // Kinematic selection + if (track.pt() <= cTrackMinPt || track.pt() >= cTrackMaxPt || std::abs(track.eta()) >= cTrackEtaCut) { + return false; + } + + // Global track selection + if (cTrackGlobal && !track.isGlobalTrackWoDCA()) { + return false; + } + + // DCA selection + if (std::abs(track.dcaXY()) >= cTrackDcaXYCut || std::abs(track.dcaZ()) >= cTrackDcaZCut) { + return false; + } + + // Electron rejection + if (std::abs(track.tpcNSigmaPi()) > cTpcRejCut && std::abs(track.tpcNSigmaKa()) > cTpcRejCut && std::abs(track.tpcNSigmaPr()) > cTpcRejCut && track.tpcNSigmaEl() > cTpcElRejCutMin && track.tpcNSigmaEl() < cTpcElRejCutMax) { + return false; + } + + // All selection passed + return true; + } + + template + float GetCorrectionFactor(T const& track, S const& sign) + { + if (!cGetCorrectionFlag) { + return 1.; + } + + if constexpr (rec_gen == kGen) { + return 1.; + } + + int binarray[4]; + std::array v = {cent, posz, track.pt(), track.eta()}; + if (sign > 0) { + for (int i = 0; i < nHistDim; ++i) { + binarray[i] = corrHist.hRecEffP->GetAxis(i)->FindBin(v[i]); + } + return corrHist.hRecEffP->GetBinContent(corrHist.hRecEffP->GetBin(binarray)); + } else if (sign < 0) { + for (int i = 0; i < nHistDim; ++i) { + binarray[i] = corrHist.hRecEffM->GetAxis(i)->FindBin(v[i]); + } + return corrHist.hRecEffM->GetBinContent(corrHist.hRecEffM->GetBin(binarray)); + } + + return 1.; + } + + template + float GetNuaCorrectionFactor(T const& track, S const& sign) + { + if (!cGetNuaCorrectionFlag) { + return 1.; + } + + if constexpr (rec_gen == kGen) { + return 1.; + } + + if (sign > 0) { + return corrHist.hNuaP->GetBinContent(corrHist.hNuaP->FindBin(posz, track.eta(), track.phi())); + } else if (sign < 0) { + return corrHist.hNuaM->GetBinContent(corrHist.hNuaM->FindBin(posz, track.eta(), track.phi())); + } + + return 1.; + } + + template + bool isClosePair(T const& p1, T const& p2) + { + // Before + histos.fill(HIST("QA/TwoTrackCut/Before/h2d_n2_detadphi"), p1.eta() - p2.eta(), RecoDecay::constrainAngle(p1.phi() - p2.phi(), -PIHalf)); + + // DPhiStar average over TPC + float dphistar = 0., n = 0.; + for (auto const& radii : vTpcRadii) { + float arg1 = 0.15 * magField * radii / p1.pt(); + float arg2 = 0.15 * magField * radii / p2.pt(); + if (std::abs(arg1) < 1.0 && std::abs(arg2) < 1.0) { + dphistar += (p1.phi() - p2.phi() - (p1.sign() * std::abs(std::asin(arg1))) + (p2.sign() * std::abs(std::asin(arg2)))); + ++n; + } + } + + // Nan check + if (n == 0) { + return false; + } + + // DPhistar + dphistar = RecoDecay::constrainAngle(dphistar / n, -PIHalf); + + // DEta + float deta = p1.eta() - p2.eta(); + + // Return flag + bool retFlag = (std::abs(deta) < cDEtaCut && std::abs(dphistar) < cDPhiStarCut); + + // Before + if (!retFlag) { + histos.fill(HIST("QA/TwoTrackCut/After/h2d_n2_detadphi"), deta, RecoDecay::constrainAngle(p1.phi() - p2.phi(), -PIHalf)); + } + + return retFlag; + } + + template + bool isCloseQinv(T const& p1, T const& p2) + { + // Before + histos.fill(HIST("QA/FemtoCut/Before/h2d_n2_detadphi"), p1.eta() - p2.eta(), RecoDecay::constrainAngle(p1.phi() - p2.phi(), -PIHalf)); + float dpt = std::abs(p1.pt() - p2.pt()); + bool retFlag = (dpt < cFemtoCut); + + if (!retFlag) { + histos.fill(HIST("QA/FemtoCut/After/h2d_n2_detadphi"), p1.eta() - p2.eta(), RecoDecay::constrainAngle(p1.phi() - p2.phi(), -PIHalf)); + } + + return retFlag; + } + + template + void fillPairHist(T const& trk_1, T const& trk_2, S const& sign_1, S const& sign_2) + { + // Check for same index + if (trk_1.index() == trk_2.index()) { + return; + } + + // Close pair rejection + if constexpr (rec_gen == kRec) { + if (trk_1.sign() * trk_2.sign() < 0) { + // Close pair + if (cApplyTwoTrackCut && isClosePair(trk_1, trk_2)) { + return; + } + // Femto selection + if (cApplyFemtoSel && isCloseQinv(trk_1, trk_2)) { + return; + } + } + } + + // Reco/Gen Dir + static constexpr auto SubDirRecGen = std::array{"Reco/", "McGen/"}; + + // Rap-phi bins + const auto rapbin1 = static_cast((trk_1.eta() - kminrap) / rapbinwidth); + const auto rapbin2 = static_cast((trk_2.eta() - kminrap) / rapbinwidth); + const auto phibin1 = static_cast(trk_1.phi() / phibinwidth); + const auto phibin2 = static_cast(trk_2.phi() / phibinwidth); + + float corfac = GetCorrectionFactor(trk_1, sign_1) * GetCorrectionFactor(trk_2, sign_2); + + if (rapbin1 >= 0 && rapbin2 >= 0 && phibin1 >= 0 && phibin2 >= 0 && rapbin1 < nrapbins && rapbin2 < nrapbins && phibin1 < nphibins && phibin2 < nphibins) { + + int rapphix = rapbin1 * nphibins + phibin1; + int rapphiy = rapbin2 * nphibins + phibin2; + + if ((sign_1 > 0 && sign_2 < 0) || (sign_1 < 0 && sign_2 > 0)) { + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_n2_rapphi_PM"), cent, rapphix + 0.5, rapphiy + 0.5, corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_ptpt_rapphi_PM"), cent, rapphix + 0.5, rapphiy + 0.5, trk_1.pt() * trk_2.pt() * corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_npt_rapphi_PM"), cent, rapphix + 0.5, rapphiy + 0.5, trk_1.pt() * corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_ptn_rapphi_PM"), cent, rapphix + 0.5, rapphiy + 0.5, trk_2.pt() * corfac); + } else { + if (sign_1 > 0 && sign_2 > 0) { + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_n2_rapphi_PP"), cent, rapphix + 0.5, rapphiy + 0.5, corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_ptpt_rapphi_PP"), cent, rapphix + 0.5, rapphiy + 0.5, trk_1.pt() * trk_2.pt() * corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_npt_rapphi_PP"), cent, rapphix + 0.5, rapphiy + 0.5, trk_1.pt() * corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_ptn_rapphi_PP"), cent, rapphix + 0.5, rapphiy + 0.5, trk_2.pt() * corfac); + } else if (sign_1 < 0 && sign_2 < 0) { + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_n2_rapphi_MM"), cent, rapphix + 0.5, rapphiy + 0.5, corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_ptpt_rapphi_MM"), cent, rapphix + 0.5, rapphiy + 0.5, trk_1.pt() * trk_2.pt() * corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_npt_rapphi_MM"), cent, rapphix + 0.5, rapphiy + 0.5, trk_1.pt() * corfac); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_ptn_rapphi_MM"), cent, rapphix + 0.5, rapphiy + 0.5, trk_2.pt() * corfac); + } + } + } + } + + template + void fillTrackQA(T const& track) + { + histos.fill(HIST("TrackQA/hPtDcaZ"), track.pt(), track.dcaZ()); + histos.fill(HIST("TrackQA/hPtDcaXY"), track.pt(), track.dcaXY()); + histos.fill(HIST("TrackQA/hTrackTPCdEdX"), track.pt(), track.tpcSignal()); + } + + template + void fillSingleHist(T const& track, S const& sign) + { + // Hist array + static constexpr auto SubDirRecGen = std::array{"Reco/", "McGen/"}; + + // Correction factor + float corrFact = GetCorrectionFactor(track, sign); + float nuaCorr = GetNuaCorrectionFactor(track, sign); + + // Histograms + if (sign > 0) { + // Corrections + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("Efficiency/h1f_n1_pt_P"), track.pt()); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("Efficiency/h4f_n1_centvzptrap_P"), cent, posz, track.pt(), track.eta()); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("NUA/h3f_n1_vzrapphi_P"), posz, track.eta(), track.phi()); + + // Checks + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h1f_n1_pt_P"), track.pt(), corrFact); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h1f_n1_rap_P"), track.eta(), corrFact); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h1f_n1_phi_P"), track.phi(), nuaCorr); + + // R2 Rho1 (Eta,Phi) + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_n1_rapphi_P"), cent, track.eta(), track.phi(), corrFact); + + // P2 Rho1 (Eta,Phi) + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_pt_rapphi_P"), cent, track.eta(), track.phi(), track.pt() * corrFact); + } else if (sign < 0) { + // Corrections + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("Efficiency/h1f_n1_pt_M"), track.pt()); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("Efficiency/h4f_n1_centvzptrap_M"), cent, posz, track.pt(), track.eta()); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("NUA/h3f_n1_vzrapphi_M"), posz, track.eta(), track.phi()); + + // Checks + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h1f_n1_pt_M"), track.pt(), corrFact); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h1f_n1_rap_M"), track.eta(), corrFact); + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h1f_n1_phi_M"), track.phi(), nuaCorr); + + // R2 Rho1 (Eta,Phi) + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_n1_rapphi_M"), cent, track.eta(), track.phi(), corrFact); + + // P2 Rho1 (Eta,Phi) + histos.fill(HIST(SubDirRecGen[rec_gen]) + HIST("h3f_pt_rapphi_M"), cent, track.eta(), track.phi(), track.pt() * corrFact); + } + } + + template + void analyzeCollision(B const&, C const& collision, T const& tracks) + { + // Select Collision (Only for Data... McRec has been selected already !!!) + if constexpr (dmc == kData) { + if (!selCollision(collision)) { + return; + } + } + + // Get magnetic field + magField = getMagneticField(collision.template foundBC_as().timestamp()); + + // Fill Event QA + histos.fill(HIST("Event/h1f_collision_cent"), cent); + histos.fill(HIST("Event/h1f_collision_Vz"), collision.posZ()); + + // Loop over tracks + for (auto const& track_1 : tracks) { + // Check for Mc matched particle + if constexpr (dmc == kMC) { + if (!track_1.has_mcParticle()) { + continue; + } + } + + // Select track + if (!selectTrack(track_1)) + continue; + + // Fill QA + fillTrackQA(track_1); + + // Fill singles + fillSingleHist(track_1, track_1.sign()); + + // Pair + for (auto const& track_2 : tracks) { + // Check for Mc matched particle + if constexpr (dmc == kMC) { + if (!track_2.has_mcParticle()) { + continue; + } + } + + // Select track + if (!selectTrack(track_2)) + continue; + + // Fill pair histograms + fillPairHist(track_1, track_2, track_1.sign(), track_2.sign()); + } + } + } + + template + bool selectMcParticle(M const& mcpart) + { + // Check for Primary Charged Particle + if (!mcpart.isPhysicalPrimary()) { + return false; + } + + // Check pdg info + auto pdgCode = mcpart.pdgCode(); + auto pdgInfo = pdg->GetParticle(pdgCode); + if (pdgInfo == nullptr) { // particle with unknown pdg code + return false; + } + + // Remove electron + if (pdgCode == kElectron || pdgCode == kPositron) { + return false; + } + + // Apply kinematic selection + if (mcpart.pt() <= cTrackMinPt || mcpart.pt() >= cTrackMaxPt || std::abs(mcpart.eta()) >= cTrackEtaCut) { + return false; + } + + // All selection passed + return true; + } + + // MC Gen analysis + template + void analyzeMcGenCollision(M const&, P const& mcParticles) + { + // Loop 1 over MC particles + for (auto const& mcpart_1 : mcParticles) { + // Select mc particle + if (!selectMcParticle(mcpart_1)) { + continue; + } + + // Fill single particle densities + int sign_1 = pdg->GetParticle(mcpart_1.pdgCode())->Charge() / 3; + fillSingleHist(mcpart_1, sign_1); + + // Loop 2 + for (auto const& mcpart_2 : mcParticles) { + // Select mc particle + if (!selectMcParticle(mcpart_2)) { + continue; + } + + // Fill pair densities + int sign_2 = pdg->GetParticle(mcpart_2.pdgCode())->Charge() / 3; + fillPairHist(mcpart_1, mcpart_2, sign_1, sign_2); + } + } + } + + // MC Reco-Gen analysis + template + void analyzeMcRecoGen(M const& mcCollision, C const& collisions, B const& bc, T const& tracks, P const& mcParticles) + { + // Number of Rec Collisions Associated to the McGen Collision + int nRecCols = collisions.size(); + if (nRecCols != 0) { + histos.fill(HIST("McGen/h1f_collision_recgen"), nRecCols); + } + // Do not analyze if more than one reco collision is accociated to one mc gen collision + if (nRecCols != 1) { + return; + } + // Check the reco collision + if (!collisions.begin().has_mcCollision() || !selCollision(collisions.begin()) || collisions.begin().mcCollisionId() != mcCollision.globalIndex()) { + return; + } + histos.fill(HIST("McGen/h2f_collision_posZ"), mcCollision.posZ(), collisions.begin().posZ()); + auto tracksThisCollision = tracks.sliceBy(tracksPerCollision, collisions.begin().globalIndex()); + analyzeCollision(bc, collisions.begin(), tracksThisCollision); + analyzeMcGenCollision(mcCollision, mcParticles); + } + + // BF, Collision and Track Table + using BCsRun3 = soa::Join; + using Collisions = soa::Join; + using Tracks = soa::Join; + using TracksMC = soa::Join; + + SliceCache cache; + Preslice tracksPerCollision = aod::track::collisionId; + + void processDummy(aod::Collisions const&) {} + + PROCESS_SWITCH(ChargeBalanceFunction, processDummy, "Dummy process", true); + + void processChargedData(Collisions::iterator const& collision, BCsRun3 const& bc, Tracks const& tracks) + { + analyzeCollision(bc, collision, tracks); + } + + PROCESS_SWITCH(ChargeBalanceFunction, processChargedData, "Charged particle BF process", false); + + void processMCRecoGen(aod::McCollisions::iterator const& mcCollision, + soa::SmallGroups> const& collisions, BCsRun3 const& bc, + TracksMC const& tracks, + aod::McParticles const& mcParticles) + { + analyzeMcRecoGen(mcCollision, collisions, bc, tracks, mcParticles); + } + + PROCESS_SWITCH(ChargeBalanceFunction, processMCRecoGen, "Process for MC RecoGen", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} From c347844ccbf50cbca42c47b69d834d78c1c7d647 Mon Sep 17 00:00:00 2001 From: Yash Patley <52608802+yashpatley@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:33:04 +0530 Subject: [PATCH 2/6] Add charge balance function workflow to CMakeLists --- PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt index 2001bcdde26..ea65ecf1162 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt +++ b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt @@ -122,3 +122,8 @@ o2physics_add_dpl_workflow(two-particle-correlations-mpi SOURCES twoParticleCorrelationsMpi.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(charge-balance-function + SOURCES chargeBalanceFunction.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) From d4bcf77600f03036171c51f4934182db63f57cf2 Mon Sep 17 00:00:00 2001 From: Yash Patley <52608802+yashpatley@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:38:50 +0530 Subject: [PATCH 3/6] Rename functions to follow camelCase convention --- .../Tasks/chargeBalanceFunction.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx index 77d959447ab..3add53b023a 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file lambdaR2Correlation.cxx +/// \file chargeBalanceFunction.cxx /// \brief R2 P2 and BF of charged hadrons. /// \author Yash Patley @@ -289,12 +289,12 @@ struct ChargeBalanceFunction { ccdbObjNuaCorr = ccdb->getForTimeStamp(cPathCCDBNuaCorr.value, nolaterthan.value); // Load reco eff corrections - LoadRecoEfficiencyHistograms(); + loadRecoEfficiencyHistograms(); } } // Load reco efficiency histograms - void LoadRecoEfficiencyHistograms() + void loadRecoEfficiencyHistograms() { // Efficiency correction histograms corrHist.hRecEffP = dynamic_cast(ccdbObjRecoEff->FindObject("h_RecEff_P")); @@ -395,7 +395,7 @@ struct ChargeBalanceFunction { } template - float GetCorrectionFactor(T const& track, S const& sign) + float getCorrectionFactor(T const& track, S const& sign) { if (!cGetCorrectionFlag) { return 1.; @@ -527,7 +527,7 @@ struct ChargeBalanceFunction { const auto phibin1 = static_cast(trk_1.phi() / phibinwidth); const auto phibin2 = static_cast(trk_2.phi() / phibinwidth); - float corfac = GetCorrectionFactor(trk_1, sign_1) * GetCorrectionFactor(trk_2, sign_2); + float corfac = getCorrectionFactor(trk_1, sign_1) * getCorrectionFactor(trk_2, sign_2); if (rapbin1 >= 0 && rapbin2 >= 0 && phibin1 >= 0 && phibin2 >= 0 && rapbin1 < nrapbins && rapbin2 < nrapbins && phibin1 < nphibins && phibin2 < nphibins) { @@ -570,7 +570,7 @@ struct ChargeBalanceFunction { static constexpr auto SubDirRecGen = std::array{"Reco/", "McGen/"}; // Correction factor - float corrFact = GetCorrectionFactor(track, sign); + float corrFact = getCorrectionFactor(track, sign); float nuaCorr = GetNuaCorrectionFactor(track, sign); // Histograms From 5027608aae75864d4f74dc5aef719ae214918583 Mon Sep 17 00:00:00 2001 From: Yash Patley <52608802+yashpatley@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:46:18 +0530 Subject: [PATCH 4/6] Update chargeBalanceFunction.cxx --- .../TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx index 3add53b023a..9bbeb4be8de 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx @@ -159,8 +159,10 @@ struct ChargeBalanceFunction { TList* ccdbObjRecoEff = nullptr; TList* ccdbObjNuaCorr = nullptr; struct CorrHist { - THnSparseF *hRecEffP, *hRecEffM; - TH3F *hNuaP, *hNuaM; + THnSparseF* hRecEffP = nullptr; + THnSparseF* hRecEffM = nullptr; + TH3F* hNuaP = nullptr; + TH3F* hNuaM = nullptr; } corrHist; int nHistDim = 0; From 4f08d8ceaa790df00c794c35ab577a5333d97e82 Mon Sep 17 00:00:00 2001 From: Yash Patley <52608802+yashpatley@users.noreply.github.com> Date: Thu, 6 Aug 2026 16:23:33 +0530 Subject: [PATCH 5/6] [PWGCF] Adding charge balance function workflow --- .../Tasks/chargeBalanceFunction.cxx | 20 ++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx index 9bbeb4be8de..a840eb1e294 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx @@ -407,18 +407,20 @@ struct ChargeBalanceFunction { return 1.; } - int binarray[4]; + std::array binarray; std::array v = {cent, posz, track.pt(), track.eta()}; if (sign > 0) { for (int i = 0; i < nHistDim; ++i) { binarray[i] = corrHist.hRecEffP->GetAxis(i)->FindBin(v[i]); } - return corrHist.hRecEffP->GetBinContent(corrHist.hRecEffP->GetBin(binarray)); - } else if (sign < 0) { + return corrHist.hRecEffP->GetBinContent(corrHist.hRecEffP->GetBin(binarray.data())); + } + + if (sign < 0) { for (int i = 0; i < nHistDim; ++i) { binarray[i] = corrHist.hRecEffM->GetAxis(i)->FindBin(v[i]); } - return corrHist.hRecEffM->GetBinContent(corrHist.hRecEffM->GetBin(binarray)); + return corrHist.hRecEffM->GetBinContent(corrHist.hRecEffM->GetBin(binarray.data())); } return 1.; @@ -437,7 +439,9 @@ struct ChargeBalanceFunction { if (sign > 0) { return corrHist.hNuaP->GetBinContent(corrHist.hNuaP->FindBin(posz, track.eta(), track.phi())); - } else if (sign < 0) { + } + + if (sign < 0) { return corrHist.hNuaM->GetBinContent(corrHist.hNuaM->FindBin(posz, track.eta(), track.phi())); } @@ -638,8 +642,9 @@ struct ChargeBalanceFunction { } // Select track - if (!selectTrack(track_1)) + if (!selectTrack(track_1)) { continue; + } // Fill QA fillTrackQA(track_1); @@ -657,8 +662,9 @@ struct ChargeBalanceFunction { } // Select track - if (!selectTrack(track_2)) + if (!selectTrack(track_2)) { continue; + } // Fill pair histograms fillPairHist(track_1, track_2, track_1.sign(), track_2.sign()); From c98072aa903201df699bc8201263aa5c9d65dbe4 Mon Sep 17 00:00:00 2001 From: Yash Patley <52608802+yashpatley@users.noreply.github.com> Date: Thu, 6 Aug 2026 19:51:13 +0530 Subject: [PATCH 6/6] Update chargeBalanceFunction.cxx --- PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx index a840eb1e294..c8b3095fe3f 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/chargeBalanceFunction.cxx @@ -407,7 +407,7 @@ struct ChargeBalanceFunction { return 1.; } - std::array binarray; + std::array binarray{}; std::array v = {cent, posz, track.pt(), track.eta()}; if (sign > 0) { for (int i = 0; i < nHistDim; ++i) {