From 082e4a4361f8f6710e25070ccff38db923b0dfd2 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Fri, 27 Mar 2026 18:02:43 +0100 Subject: [PATCH 01/12] Add reduced q-vectors --- Common/DataModel/Qvectors.h | 10 ++ Common/TableProducer/qVectorsTable.cxx | 24 ++++ PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 136 +++++++++++++++++------ 3 files changed, 139 insertions(+), 31 deletions(-) diff --git a/Common/DataModel/Qvectors.h b/Common/DataModel/Qvectors.h index d723d659bc5..78bd0ef2181 100644 --- a/Common/DataModel/Qvectors.h +++ b/Common/DataModel/Qvectors.h @@ -77,6 +77,12 @@ DECLARE_SOA_COLUMN(LabelsTPCpos, labelsTPCpos, std::vector); DECLARE_SOA_COLUMN(LabelsTPCneg, labelsTPCneg, std::vector); DECLARE_SOA_COLUMN(LabelsTPCall, labelsTPCall, std::vector); +// Normalized amplitudes for Event-shape Engineering +DECLARE_SOA_COLUMN(QVecRedFT0C, qVecRedFT0C, float); +DECLARE_SOA_COLUMN(QVecRedTpcPos, qVecRedTpcPos, float); +DECLARE_SOA_COLUMN(QVecRedTpcNeg, qVecRedTpcNeg, float); +DECLARE_SOA_COLUMN(QVecRedTpcAll, qVecRedTpcAll, float); + // Deprecated, will be removed in future after transition time // DECLARE_SOA_COLUMN(QvecBPosReVec, qvecBPosReVec, std::vector); DECLARE_SOA_COLUMN(QvecBPosImVec, qvecBPosImVec, std::vector); @@ -121,6 +127,8 @@ DECLARE_SOA_TABLE(QvectorTPCposVecs, "AOD", "QVECTORSTPCPVEC", qvec::IsCalibrate DECLARE_SOA_TABLE(QvectorTPCnegVecs, "AOD", "QVECTORSTPCNVEC", qvec::IsCalibrated, qvec::QvecTPCnegReVec, qvec::QvecTPCnegImVec, qvec::NTrkTPCneg, qvec::LabelsTPCneg); DECLARE_SOA_TABLE(QvectorTPCallVecs, "AOD", "QVECTORSTPCAVEC", qvec::IsCalibrated, qvec::QvecTPCallReVec, qvec::QvecTPCallImVec, qvec::NTrkTPCall, qvec::LabelsTPCall); +DECLARE_SOA_TABLE(QvectorsReds, "AOD", "QVECTORSRED", qvec::QVecRedFT0C, qvec::QVecRedTpcPos, qvec::QVecRedTpcNeg, qvec::QVecRedTpcAll); + using QvectorFT0C = QvectorFT0Cs::iterator; using QvectorFT0A = QvectorFT0As::iterator; using QvectorFT0M = QvectorFT0Ms::iterator; @@ -137,6 +145,8 @@ using QvectorTPCposVec = QvectorTPCposVecs::iterator; using QvectorTPCnegVec = QvectorTPCnegVecs::iterator; using QvectorTPCallVec = QvectorTPCallVecs::iterator; +using QvectorRed = QvectorsReds::iterator; + // Deprecated, will be removed in future after transition time // DECLARE_SOA_TABLE(QvectorBPoss, "AOD", "QVECTORSBPOS", qvec::IsCalibrated, qvec::QvecBPosRe, qvec::QvecBPosIm, qvec::NTrkBPos, qvec::LabelsBPos); DECLARE_SOA_TABLE(QvectorBNegs, "AOD", "QVECTORSBNEG", qvec::IsCalibrated, qvec::QvecBNegRe, qvec::QvecBNegIm, qvec::NTrkBNeg, qvec::LabelsBNeg); diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index 273394a8e67..e3e676a2e5c 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -112,6 +112,7 @@ struct qVectorsTable { Configurable cfgUseTPCpos{"cfgUseTPCpos", false, "Initial value for using TPCpos. By default obtained from DataModel."}; Configurable cfgUseTPCneg{"cfgUseTPCneg", false, "Initial value for using TPCneg. By default obtained from DataModel."}; Configurable cfgUseTPCall{"cfgUseTPCall", false, "Initial value for using TPCall. By default obtained from DataModel."}; + Configurable cfgProduceRedQVecs{"cfgProduceRedQVecs", false, "Produce reduced Q-vectors for Event-Shape Engineering"}; // Table. Produces qVector; @@ -131,6 +132,8 @@ struct qVectorsTable { Produces qVectorTPCnegVec; Produces qVectorTPCallVec; + Produces qVectorRed; + std::vector FT0RelGainConst{}; std::vector FV0RelGainConst{}; @@ -730,6 +733,27 @@ struct qVectorsTable { if (useDetector["QvectorTPCalls"]) qVectorTPCall(IsCalibrated, qvecReTPCall.at(0), qvecImTPCall.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + + double qVecRedFT0C{-999.}, qVecRedTpcPos{-999.}, qVecRedTpcNeg{-999.}, qVecRedTpcAll{-999.}; + if (cfgProduceRedQVecs) { + // Correct normalization to remove multiplicity dependence, + // taking into account that the Q-vector is normalized by 1/M + // and the EsE reduced Q-vector must be normalized to 1/sqrt(M) + if (useDetector["QvectorFT0Cs"]) { + qVecRedFT0C = TMath::Sqrt(qvecReFT0C.at(0) * qvecReFT0C.at(0) + qvecImFT0C.at(0) * qvecImFT0C.at(0)) * TMath::Sqrt(qvecAmp[kFT0C]); + } + if (useDetector["QvectorTPCposs"]) { + qVecRedTpcPos = TMath::Sqrt(qvecReTPCpos.at(0) * qvecReTPCpos.at(0) + qvecImTPCpos.at(0) * qvecImTPCpos.at(0)) * TMath::Sqrt(qvecAmp[kTPCpos]); + } + if (useDetector["QvectorTPCnegs"]) { + qVecRedTpcNeg = TMath::Sqrt(qvecReTPCneg.at(0) * qvecReTPCneg.at(0) + qvecImTPCneg.at(0) * qvecImTPCneg.at(0)) * TMath::Sqrt(qvecAmp[kTPCneg]); + } + if (useDetector["QvectorTPCalls"]) { + qVecRedTpcAll = TMath::Sqrt(qvecReTPCall.at(0) * qvecReTPCall.at(0) + qvecImTPCall.at(0) * qvecImTPCall.at(0)) * TMath::Sqrt(qvecAmp[kTPCall]); + } + } + qVectorRed(qVecRedFT0C, qVecRedTpcPos, qVecRedTpcNeg, qVecRedTpcAll); + qVectorFT0CVec(IsCalibrated, qvecReFT0C, qvecImFT0C, qvecAmp[kFT0C]); qVectorFT0AVec(IsCalibrated, qvecReFT0A, qvecImFT0A, qvecAmp[kFT0A]); qVectorFT0MVec(IsCalibrated, qvecReFT0M, qvecImFT0M, qvecAmp[kFT0M]); diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 78462e10a73..7f54543aef0 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -75,6 +75,7 @@ DECLARE_SOA_COLUMN(MlScore0, mlScore0, float); //! ML score of the first con DECLARE_SOA_COLUMN(MlScore1, mlScore1, float); //! ML score of the second configured index DECLARE_SOA_COLUMN(ScalarProd, scalarProd, float); //! Scalar product DECLARE_SOA_COLUMN(Cent, cent, float); //! Centrality +DECLARE_SOA_COLUMN(RedQVec, redQVec, float); //! Reduced Q-vector } // namespace full DECLARE_SOA_TABLE(HfCandMPtInfos, "AOD", "HFCANDMPTINFO", full::M, @@ -89,6 +90,7 @@ DECLARE_SOA_TABLE(HfCandFlowInfos, "AOD", "HFCANDFLOWINFO", full::MlScore1, full::ScalarProd, full::Cent); +DECLARE_SOA_TABLE(HfRedQVecEsEs, "AOD", "HFREDQVECESE", full::RedQVec); } // namespace o2::aod enum DecayChannel { DplusToPiKPi = 0, @@ -106,15 +108,18 @@ enum DecayChannel { DplusToPiKPi = 0, struct HfTaskFlowCharmHadrons { Produces rowCandMassPtMl; Produces rowCandMassPtMlSpCent; + Produces rowRedQVecEsE; Configurable harmonic{"harmonic", 2, "harmonic number"}; Configurable qVecDetector{"qVecDetector", 3, "Detector for Q vector estimation (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3, TPC Pos: 4, TPC Neg: 5, TPC Tot: 6)"}; + Configurable qVecRedDetector{"qVecRedDetector", 6, "Detector for Q vector estimation (FT0C: 3, TPC Pos: 4, TPC Neg: 5, TPC Tot: 6)"}; Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4, NTracksPV: 5, FT0CVariant2: 6)"}; Configurable selectionFlag{"selectionFlag", 1, "Selection Flag for hadron (e.g. 1 for skimming, 3 for topo. and kine., 7 for PID)"}; Configurable centralityMin{"centralityMin", 0., "Minimum centrality accepted in SP/EP computation (not applied in resolution process)"}; Configurable centralityMax{"centralityMax", 100., "Maximum centrality accepted in SP/EP computation (not applied in resolution process)"}; Configurable storeEP{"storeEP", false, "Flag to store EP-related axis"}; Configurable storeMl{"storeMl", false, "Flag to store ML scores"}; + Configurable storeRedQVec{"storeRedQVec", false, "Flag to store reduced Q-vectors for ESE"}; Configurable fillMassPtMlTree{"fillMassPtMlTree", false, "Flag to fill mass, pt and ML scores tree"}; Configurable fillMassPtMlSpCentTree{"fillMassPtMlSpCentTree", false, "Flag to fill mass, pt, ML scores, SP and centrality tree"}; Configurable fillSparse{"fillSparse", true, "Flag to fill sparse"}; @@ -148,7 +153,7 @@ struct HfTaskFlowCharmHadrons { using CandXic0DataWMl = soa::Filtered>; using CandD0DataWMl = soa::Filtered>; using CandD0Data = soa::Filtered>; - using CollsWithQvecs = soa::Join; + using CollsWithQvecs = soa::Join; using TracksWithExtra = soa::Join; Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; @@ -198,6 +203,7 @@ struct HfTaskFlowCharmHadrons { ConfigurableAxis thnConfigAxisResoFV0aTPCtot{"thnConfigAxisResoFV0aTPCtot", {160, -8, 8}, ""}; ConfigurableAxis thnConfigAxisCandidateEta{"thnConfigAxisCandidateEta", {100, -5, 5}, ""}; ConfigurableAxis thnConfigAxisSign{"thnConfigAxisSign", {6, -3.0, 3.0}, ""}; + ConfigurableAxis thnConfigAxisRedQVec = ConfigurableAxis{"thnConfigAxisRedQVec", {1000, 0, 100}, ""}; HistogramRegistry registry{"registry", {}}; @@ -225,6 +231,7 @@ struct HfTaskFlowCharmHadrons { const AxisSpec thnAxisNoCollInTimeRangeNarrow{thnConfigAxisNoCollInTimeRangeNarrow, "NoCollInTimeRangeNarrow"}; const AxisSpec thnAxisNoCollInTimeRangeStandard{thnConfigAxisNoCollInTimeRangeStandard, "NoCollInTimeRangeStandard"}; const AxisSpec thnAxisNoCollInRofStandard{thnConfigAxisNoCollInRofStandard, "NoCollInRofStandard"}; + const AxisSpec thnAxisRedQVec{thnConfigAxisRedQVec, "Reduced Q-vector"}; // TODO: currently only the Q vector of FT0c FV0a and TPCtot are considered const AxisSpec thnAxisResoFT0cFV0a{thnConfigAxisResoFT0cFV0a, "Q_{FT0c} #bullet Q_{FV0a}"}; const AxisSpec thnAxisResoFT0cTPCtot{thnConfigAxisResoFT0cTPCtot, "Q_{FT0c} #bullet Q_{TPCtot}"}; @@ -252,6 +259,9 @@ struct HfTaskFlowCharmHadrons { thnAxisNoCollInTimeRangeNarrow, thnAxisNoCollInTimeRangeStandard, thnAxisNoCollInRofStandard}); } } + if (storeRedQVec) { + axes.insert(axes.end(), {thnAxisRedQVec}); + } registry.add("hSparseFlowCharm", "THn for SP", HistType::kTHnSparseF, axes); registry.add("hCentEventWithCand", "Centrality distributions with charm candidates;Cent;entries", HistType::kTH1F, {{100, 0.f, 100.f}}); registry.add("hCentEventWithCandInSigRegion", "Centrality distributions with charm candidates in signal range;Cent;entries", HistType::kTH1F, {{100, 0.f, 100.f}}); @@ -288,6 +298,14 @@ struct HfTaskFlowCharmHadrons { registry.add("spReso/hSpResoFV0aTPCtot", "hSpResoFV0aTPCtot; centrality; Q_{FV0a} #bullet Q_{TPCtot}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}}); registry.add("spReso/hSpResoTPCposTPCneg", "hSpResoTPCposTPCneg; centrality; Q_{TPCpos} #bullet Q_{TPCneg}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}}); + if (storeRedQVec) { + registry.add("redQVecs/hSparseRedQVecs", "hSpResoRedQVec; centrality; Q_{red} #bullet Q_{TPCtot}", {HistType::kTHnSparseF, {thnAxisCent, thnAxisRedQVec, thnAxisRedQVec, thnAxisRedQVec, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecFT0C", "hRedQVecFT0C; centrality; q_{FT0c}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecTpcPos", "hRedQVecTpcPos; centrality; q_{TPCpos}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecTpcNeg", "hRedQVecTpcNeg; centrality; q_{TPCneg}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecTpcAll", "hRedQVecTpcAll; centrality; q_{TPCall}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + } + if (saveEpResoHisto) { registry.add("epReso/hEpResoFT0cFT0a", "hEpResoFT0cFT0a; centrality; #Delta#Psi_{sub}", {HistType::kTH2F, {thnAxisCent, thnAxisCosNPhi}}); registry.add("epReso/hEpResoFT0cFV0a", "hEpResoFT0cFV0a; centrality; #Delta#Psi_{sub}", {HistType::kTH2F, {thnAxisCent, thnAxisCosNPhi}}); @@ -430,6 +448,7 @@ struct HfTaskFlowCharmHadrons { /// \param outputMl are the ML scores /// \param occupancy is the occupancy of the collision using the track estimator /// \param hfevselflag flag of the collision associated to utilsEvSelHf.h + /// \param redQVec is the reduced Q vector for EsE void fillThn(const float mass, const float pt, const float eta, @@ -441,7 +460,8 @@ struct HfTaskFlowCharmHadrons { const float sp, const std::vector& outputMl, const float occupancy, - const o2::hf_evsel::HfCollisionRejectionMask hfevselflag) + const o2::hf_evsel::HfCollisionRejectionMask hfevselflag, + const float redQVec) { auto hSparse = registry.get(HIST("hSparseFlowCharm")); const int ndim = hSparse->GetNdimensions(); @@ -480,6 +500,9 @@ struct HfTaskFlowCharmHadrons { values.push_back(evtSelFlags[3]); values.push_back(evtSelFlags[4]); } + if (storeRedQVec) { + values.push_back(redQVec); + } if (static_cast(values.size()) != ndim) { LOGF(fatal, @@ -487,7 +510,6 @@ struct HfTaskFlowCharmHadrons { "does not match THnSparse dimensionality (%d).", static_cast(values.size()), ndim); } - hSparse->Fill(values.data()); } @@ -543,6 +565,16 @@ struct HfTaskFlowCharmHadrons { float yQVec = qVecs[1]; float const amplQVec = qVecs[2]; float const evtPl = epHelper.GetEventPlane(xQVec, yQVec, harmonic); + + float redQVec{-999.f}; + std::array qVecRedComps{-999.f, -999.f, -999.f}; + if (storeRedQVec) { + qVecRedComps = getQvec(collision, qVecRedDetector.value); + } + float xRedQVec = qVecRedComps[0]; + float yRedQVec = qVecRedComps[1]; + float const amplRedQVec = qVecRedComps[2]; + for (const auto& candidate : candidates) { float massCand = 0.; float signCand = 0.; @@ -682,9 +714,10 @@ struct HfTaskFlowCharmHadrons { } // If TPC is used for the SP estimation, the tracks of the hadron candidate must be removed from the corresponding TPC Q vector to avoid self-correlations - if (qVecDetector == QvecEstimator::TPCNeg || - qVecDetector == QvecEstimator::TPCPos || - qVecDetector == QvecEstimator::TPCTot) { + bool subtractDaugsFromQVec = (qVecDetector == QvecEstimator::TPCNeg || + qVecDetector == QvecEstimator::TPCPos || + qVecDetector == QvecEstimator::TPCTot); + if (subtractDaugsFromQVec) { std::vector tracksQx; std::vector tracksQy; @@ -722,8 +755,41 @@ struct HfTaskFlowCharmHadrons { rowCandMassPtMlSpCent(massCand, ptCand, outputMl[0], outputMl[1], scalprodCand, cent); } } + + bool subtractDaugsFromRedQVec = storeRedQVec && (qVecRedDetector == QvecEstimator::TPCNeg || + qVecRedDetector == QvecEstimator::TPCPos || + qVecRedDetector == QvecEstimator::TPCTot); + if (subtractDaugsFromRedQVec) { + std::vector tracksRedQx; + std::vector tracksRedQy; + + // IMPORTANT: use the amplitude of the reduced Q-vector to build this Q-vector + if constexpr (std::is_same_v || std::is_same_v) { + getQvecXic0Tracks(candidate, tracksRedQx, tracksRedQy, amplRedQVec, static_cast(qVecRedDetector.value)); + } else { + getQvecDtracks(candidate, tracksRedQx, tracksRedQy, amplRedQVec, static_cast(qVecRedDetector.value)); + } + + // subtract daughters' contribution from the (normalized) Q-vector + for (std::size_t iTrack = 0; iTrack < tracksRedQx.size(); ++iTrack) { + xRedQVec -= tracksRedQx[iTrack]; + yRedQVec -= tracksRedQy[iTrack]; + } + + if (qVecRedDetector.value == QvecEstimator::TPCTot || qVecRedDetector.value == QvecEstimator::TPCPos || qVecRedDetector.value == QvecEstimator::TPCNeg) { + // Correct for track multiplicity + redQVec = TMath::Sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * amplRedQVec / TMath::Sqrt(amplRedQVec - tracksRedQx.size()); + } + else { + redQVec = TMath::Sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * TMath::Sqrt(amplRedQVec); + } + } + if (storeRedQVec) { + rowRedQVecEsE(redQVec); + } if (fillSparse) { - fillThn(massCand, ptCand, etaCand, signCand, cent, cosNPhi, sinNPhi, cosDeltaPhi, scalprodCand, outputMl, occupancy, hfevflag); + fillThn(massCand, ptCand, etaCand, signCand, cent, cosNPhi, sinNPhi, + cosDeltaPhi, scalprodCand, outputMl, occupancy, hfevflag, redQVec); } } if (hasCandInMassWin) { @@ -878,12 +944,12 @@ struct HfTaskFlowCharmHadrons { float const yQVecFT0m = collision.qvecFT0MIm(); float const xQVecFV0a = collision.qvecFV0ARe(); float const yQVecFV0a = collision.qvecFV0AIm(); - float const xQVecBPos = collision.qvecBPosRe(); - float const yQVecBPos = collision.qvecBPosIm(); - float const xQVecBNeg = collision.qvecBNegRe(); - float const yQVecBNeg = collision.qvecBNegIm(); - float const xQVecBTot = collision.qvecBTotRe(); - float const yQVecBTot = collision.qvecBTotIm(); + float const xQVecTPCPos = collision.qvecTPCposRe(); + float const yQVecTPCPos = collision.qvecTPCposIm(); + float const xQVecTPCNeg = collision.qvecTPCnegRe(); + float const yQVecTPCNeg = collision.qvecTPCnegIm(); + float const xQVecTPCAll = collision.qvecTPCallRe(); + float const yQVecTPCAll = collision.qvecTPCallIm(); centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); if (storeResoOccu) { @@ -892,8 +958,8 @@ struct HfTaskFlowCharmHadrons { const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); std::vector evtSelFlags = getEventSelectionFlags(rejectionMask); registry.fill(HIST("spReso/hSparseReso"), centrality, xQVecFT0c * xQVecFV0a + yQVecFT0c * yQVecFV0a, - xQVecFT0c * xQVecBTot + yQVecFT0c * yQVecBTot, - xQVecFV0a * xQVecBTot + yQVecFV0a * yQVecBTot, + xQVecFT0c * xQVecTPCAll + yQVecFT0c * yQVecTPCAll, + xQVecFV0a * xQVecTPCAll + yQVecFV0a * yQVecTPCAll, occupancy, evtSelFlags[0], evtSelFlags[1], evtSelFlags[2], evtSelFlags[3], evtSelFlags[4]); } @@ -901,6 +967,9 @@ struct HfTaskFlowCharmHadrons { registry.fill(HIST("hSparseCentEstimators"), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(0)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(1)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(2)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(3))); } + if (storeRedQVec) { + registry.fill(HIST("redQVecs/hSparseRedQVecs"), centrality, collision.qVecRedFT0C(), collision.qVecRedTpcPos(), collision.qVecRedTpcNeg(), collision.qVecRedTpcAll()); + } if (!isCollSelected(collision, bcs, centrality)) { // no selection on the centrality is applied, but on event selection flags @@ -909,30 +978,35 @@ struct HfTaskFlowCharmHadrons { registry.fill(HIST("spReso/hSpResoFT0cFT0a"), centrality, xQVecFT0c * xQVecFT0a + yQVecFT0c * yQVecFT0a); registry.fill(HIST("spReso/hSpResoFT0cFV0a"), centrality, xQVecFT0c * xQVecFV0a + yQVecFT0c * yQVecFV0a); - registry.fill(HIST("spReso/hSpResoFT0cTPCpos"), centrality, xQVecFT0c * xQVecBPos + yQVecFT0c * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0cTPCneg"), centrality, xQVecFT0c * xQVecBNeg + yQVecFT0c * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0cTPCtot"), centrality, xQVecFT0c * xQVecBTot + yQVecFT0c * yQVecBTot); + registry.fill(HIST("spReso/hSpResoFT0cTPCpos"), centrality, xQVecFT0c * xQVecTPCPos + yQVecFT0c * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFT0cTPCneg"), centrality, xQVecFT0c * xQVecTPCNeg + yQVecFT0c * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFT0cTPCtot"), centrality, xQVecFT0c * xQVecTPCAll + yQVecFT0c * yQVecTPCAll); registry.fill(HIST("spReso/hSpResoFT0aFV0a"), centrality, xQVecFT0a * xQVecFV0a + yQVecFT0a * yQVecFV0a); - registry.fill(HIST("spReso/hSpResoFT0aTPCpos"), centrality, xQVecFT0a * xQVecBPos + yQVecFT0a * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0aTPCneg"), centrality, xQVecFT0a * xQVecBNeg + yQVecFT0a * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0aTPCtot"), centrality, xQVecFT0a * xQVecBTot + yQVecFT0a * yQVecBTot); + registry.fill(HIST("spReso/hSpResoFT0aTPCpos"), centrality, xQVecFT0a * xQVecTPCPos + yQVecFT0a * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFT0aTPCneg"), centrality, xQVecFT0a * xQVecTPCNeg + yQVecFT0a * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFT0aTPCtot"), centrality, xQVecFT0a * xQVecTPCAll + yQVecFT0a * yQVecTPCAll); registry.fill(HIST("spReso/hSpResoFT0mFV0a"), centrality, xQVecFT0m * xQVecFV0a + yQVecFT0m * yQVecFV0a); - registry.fill(HIST("spReso/hSpResoFT0mTPCpos"), centrality, xQVecFT0m * xQVecBPos + yQVecFT0m * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0mTPCneg"), centrality, xQVecFT0m * xQVecBNeg + yQVecFT0m * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0mTPCtot"), centrality, xQVecFT0m * xQVecBTot + yQVecFT0m * yQVecBTot); - registry.fill(HIST("spReso/hSpResoFV0aTPCpos"), centrality, xQVecFV0a * xQVecBPos + yQVecFV0a * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFV0aTPCneg"), centrality, xQVecFV0a * xQVecBNeg + yQVecFV0a * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecBTot + yQVecFV0a * yQVecBTot); - registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecBPos * xQVecBNeg + yQVecBPos * yQVecBNeg); + registry.fill(HIST("spReso/hSpResoFT0mTPCpos"), centrality, xQVecFT0m * xQVecTPCPos + yQVecFT0m * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFT0mTPCneg"), centrality, xQVecFT0m * xQVecTPCNeg + yQVecFT0m * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFT0mTPCtot"), centrality, xQVecFT0m * xQVecTPCAll + yQVecFT0m * yQVecTPCAll); + registry.fill(HIST("spReso/hSpResoFV0aTPCpos"), centrality, xQVecFV0a * xQVecTPCPos + yQVecFV0a * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFV0aTPCneg"), centrality, xQVecFV0a * xQVecTPCNeg + yQVecFV0a * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecTPCAll + yQVecFV0a * yQVecTPCAll); + registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecTPCPos * xQVecTPCNeg + yQVecTPCPos * yQVecTPCNeg); + + registry.fill(HIST("redQVecs/hRedQVecFT0C"), centrality, collision.qVecRedFT0C()); + registry.fill(HIST("redQVecs/hRedQVecTpcPos"), centrality, collision.qVecRedTpcPos()); + registry.fill(HIST("redQVecs/hRedQVecTpcNeg"), centrality, collision.qVecRedTpcNeg()); + registry.fill(HIST("redQVecs/hRedQVecTpcAll"), centrality, collision.qVecRedTpcAll()); if (saveEpResoHisto) { float const epFT0a = epHelper.GetEventPlane(xQVecFT0a, yQVecFT0a, harmonic); float const epFT0c = epHelper.GetEventPlane(xQVecFT0c, yQVecFT0c, harmonic); float const epFT0m = epHelper.GetEventPlane(xQVecFT0m, yQVecFT0m, harmonic); float const epFV0a = epHelper.GetEventPlane(xQVecFV0a, yQVecFV0a, harmonic); - float const epBPoss = epHelper.GetEventPlane(xQVecBPos, yQVecBPos, harmonic); - float const epBNegs = epHelper.GetEventPlane(xQVecBNeg, yQVecBNeg, harmonic); - float const epBTots = epHelper.GetEventPlane(xQVecBTot, yQVecBTot, harmonic); + float const epBPoss = epHelper.GetEventPlane(xQVecTPCPos, yQVecTPCPos, harmonic); + float const epBNegs = epHelper.GetEventPlane(xQVecTPCNeg, yQVecTPCNeg, harmonic); + float const epBTots = epHelper.GetEventPlane(xQVecTPCAll, yQVecTPCAll, harmonic); registry.fill(HIST("epReso/hEpResoFT0cFT0a"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epFT0a, harmonic))); registry.fill(HIST("epReso/hEpResoFT0cFV0a"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epFV0a, harmonic))); From 74a9512c44055e29e67a4161705d000b0a75947a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 2 Apr 2026 19:05:19 +0000 Subject: [PATCH 02/12] Please consider the following formatting changes --- Common/TableProducer/qVectorsTable.cxx | 1 - PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 5 ++--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index e3e676a2e5c..489fd0dc2c9 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -733,7 +733,6 @@ struct qVectorsTable { if (useDetector["QvectorTPCalls"]) qVectorTPCall(IsCalibrated, qvecReTPCall.at(0), qvecImTPCall.at(0), qvecAmp[kTPCall], TrkTPCallLabel); - double qVecRedFT0C{-999.}, qVecRedTpcPos{-999.}, qVecRedTpcNeg{-999.}, qVecRedTpcAll{-999.}; if (cfgProduceRedQVecs) { // Correct normalization to remove multiplicity dependence, diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 7f54543aef0..aa482e9ad6f 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -779,8 +779,7 @@ struct HfTaskFlowCharmHadrons { if (qVecRedDetector.value == QvecEstimator::TPCTot || qVecRedDetector.value == QvecEstimator::TPCPos || qVecRedDetector.value == QvecEstimator::TPCNeg) { // Correct for track multiplicity redQVec = TMath::Sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * amplRedQVec / TMath::Sqrt(amplRedQVec - tracksRedQx.size()); - } - else { + } else { redQVec = TMath::Sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * TMath::Sqrt(amplRedQVec); } } @@ -788,7 +787,7 @@ struct HfTaskFlowCharmHadrons { rowRedQVecEsE(redQVec); } if (fillSparse) { - fillThn(massCand, ptCand, etaCand, signCand, cent, cosNPhi, sinNPhi, + fillThn(massCand, ptCand, etaCand, signCand, cent, cosNPhi, sinNPhi, cosDeltaPhi, scalprodCand, outputMl, occupancy, hfevflag, redQVec); } } From e605c551f5418c951ba51142872052456337b140 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Thu, 2 Apr 2026 21:19:38 +0200 Subject: [PATCH 03/12] Fix linter --- PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index aa482e9ad6f..a86f0a8752a 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -203,7 +203,7 @@ struct HfTaskFlowCharmHadrons { ConfigurableAxis thnConfigAxisResoFV0aTPCtot{"thnConfigAxisResoFV0aTPCtot", {160, -8, 8}, ""}; ConfigurableAxis thnConfigAxisCandidateEta{"thnConfigAxisCandidateEta", {100, -5, 5}, ""}; ConfigurableAxis thnConfigAxisSign{"thnConfigAxisSign", {6, -3.0, 3.0}, ""}; - ConfigurableAxis thnConfigAxisRedQVec = ConfigurableAxis{"thnConfigAxisRedQVec", {1000, 0, 100}, ""}; + ConfigurableAxis thnConfigAxisRedQVec{"thnConfigAxisRedQVec", {1000, 0, 100}, ""}; HistogramRegistry registry{"registry", {}}; From e796572d328b8a3ebb0dbe6eb3136306fa472501 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Thu, 2 Apr 2026 21:24:30 +0200 Subject: [PATCH 04/12] Fix linter 2 --- PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index a86f0a8752a..8f2b2b07734 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -778,9 +778,9 @@ struct HfTaskFlowCharmHadrons { if (qVecRedDetector.value == QvecEstimator::TPCTot || qVecRedDetector.value == QvecEstimator::TPCPos || qVecRedDetector.value == QvecEstimator::TPCNeg) { // Correct for track multiplicity - redQVec = TMath::Sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * amplRedQVec / TMath::Sqrt(amplRedQVec - tracksRedQx.size()); + redQVec = std::sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * amplRedQVec / std::sqrt(amplRedQVec - tracksRedQx.size()); } else { - redQVec = TMath::Sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * TMath::Sqrt(amplRedQVec); + redQVec = std::sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * std::sqrt(amplRedQVec); } } if (storeRedQVec) { From 1dca5af3f73b6ab42bb952538a5ee776a668ccba Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 7 Apr 2026 11:11:50 +0200 Subject: [PATCH 05/12] Implement Vit comments --- Common/TableProducer/qVectorsTable.cxx | 9 ++--- PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 44 +++++++++++++----------- 2 files changed, 29 insertions(+), 24 deletions(-) diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index 489fd0dc2c9..cd8ac007be0 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -50,6 +50,7 @@ #include #include +#include #include #include #include @@ -739,16 +740,16 @@ struct qVectorsTable { // taking into account that the Q-vector is normalized by 1/M // and the EsE reduced Q-vector must be normalized to 1/sqrt(M) if (useDetector["QvectorFT0Cs"]) { - qVecRedFT0C = TMath::Sqrt(qvecReFT0C.at(0) * qvecReFT0C.at(0) + qvecImFT0C.at(0) * qvecImFT0C.at(0)) * TMath::Sqrt(qvecAmp[kFT0C]); + qVecRedFT0C = std::hypot(qvecReFT0C.at(0), qvecImFT0C.at(0)) * std::sqrt(qvecAmp[kFT0C]); } if (useDetector["QvectorTPCposs"]) { - qVecRedTpcPos = TMath::Sqrt(qvecReTPCpos.at(0) * qvecReTPCpos.at(0) + qvecImTPCpos.at(0) * qvecImTPCpos.at(0)) * TMath::Sqrt(qvecAmp[kTPCpos]); + qVecRedTpcPos = std::hypot(qvecReTPCpos.at(0), qvecImTPCpos.at(0)) * std::sqrt(qvecAmp[kTPCpos]); } if (useDetector["QvectorTPCnegs"]) { - qVecRedTpcNeg = TMath::Sqrt(qvecReTPCneg.at(0) * qvecReTPCneg.at(0) + qvecImTPCneg.at(0) * qvecImTPCneg.at(0)) * TMath::Sqrt(qvecAmp[kTPCneg]); + qVecRedTpcNeg = std::hypot(qvecReTPCneg.at(0), qvecImTPCneg.at(0)) * std::sqrt(qvecAmp[kTPCneg]); } if (useDetector["QvectorTPCalls"]) { - qVecRedTpcAll = TMath::Sqrt(qvecReTPCall.at(0) * qvecReTPCall.at(0) + qvecImTPCall.at(0) * qvecImTPCall.at(0)) * TMath::Sqrt(qvecAmp[kTPCall]); + qVecRedTpcAll = std::hypot(qvecReTPCall.at(0), qvecImTPCall.at(0)) * std::sqrt(qvecAmp[kTPCall]); } } qVectorRed(qVecRedFT0C, qVecRedTpcPos, qVecRedTpcNeg, qVecRedTpcAll); diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 8f2b2b07734..0afa8d2927b 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -523,13 +523,16 @@ struct HfTaskFlowCharmHadrons { aod::BCsWithTimestamps const&, float& centrality) { - const auto occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + float occupancy{-999.f}; + if (occEstimator != 0) { + occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); + } const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); centrality = o2::hf_centrality::getCentralityColl(collision, CentEstimator); /// monitor the satisfied event selections hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); - registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); return rejectionMask == 0; } @@ -713,10 +716,15 @@ struct HfTaskFlowCharmHadrons { etaCand = candidate.eta(); } + float const cosNPhi = std::cos(harmonic * phiCand); + float const sinNPhi = std::sin(harmonic * phiCand); + float const cosDeltaPhi = std::cos(harmonic * (phiCand - evtPl)); + // If TPC is used for the SP estimation, the tracks of the hadron candidate must be removed from the corresponding TPC Q vector to avoid self-correlations bool subtractDaugsFromQVec = (qVecDetector == QvecEstimator::TPCNeg || qVecDetector == QvecEstimator::TPCPos || qVecDetector == QvecEstimator::TPCTot); + float scalprodCand{-999.f}; if (subtractDaugsFromQVec) { std::vector tracksQx; @@ -730,17 +738,13 @@ struct HfTaskFlowCharmHadrons { } // subtract daughters' contribution from the (normalized) Q-vector - for (std::size_t iTrack = 0; iTrack < tracksQx.size(); ++iTrack) { - xQVec -= tracksQx[iTrack]; - yQVec -= tracksQy[iTrack]; - } + float xQVecDaugSubtr = xQVec - std::accumulate(tracksQx.begin(), tracksQx.end(), 0.0); + float yQVecDaugSubtr = yQVec - std::accumulate(tracksQy.begin(), tracksQy.end(), 0.0); + scalprodCand = cosNPhi * xQVecDaugSubtr + sinNPhi * yQVecDaugSubtr; + } else { + scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; } - float const cosNPhi = std::cos(harmonic * phiCand); - float const sinNPhi = std::sin(harmonic * phiCand); - float const scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; - float const cosDeltaPhi = std::cos(harmonic * (phiCand - evtPl)); - if (fillMassPtMlTree || fillMassPtMlSpCentTree) { if (downSampleFactor < 1.) { float const pseudoRndm = ptCand * 1000. - static_cast(ptCand * 1000); @@ -771,16 +775,13 @@ struct HfTaskFlowCharmHadrons { } // subtract daughters' contribution from the (normalized) Q-vector - for (std::size_t iTrack = 0; iTrack < tracksRedQx.size(); ++iTrack) { - xRedQVec -= tracksRedQx[iTrack]; - yRedQVec -= tracksRedQy[iTrack]; - } - + const float redQVecXDaugSubtr = xRedQVec - std::accumulate(tracksRedQx.begin(), tracksRedQx.end(), 0.0); + const float redQVecYDaugSubtr = yRedQVec - std::accumulate(tracksRedQy.begin(), tracksRedQy.end(), 0.0); if (qVecRedDetector.value == QvecEstimator::TPCTot || qVecRedDetector.value == QvecEstimator::TPCPos || qVecRedDetector.value == QvecEstimator::TPCNeg) { // Correct for track multiplicity - redQVec = std::sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * amplRedQVec / std::sqrt(amplRedQVec - tracksRedQx.size()); + redQVec = std::hypot(redQVecXDaugSubtr, redQVecYDaugSubtr) * amplRedQVec / std::sqrt(amplRedQVec - tracksRedQx.size()); } else { - redQVec = std::sqrt(xRedQVec * xRedQVec + yRedQVec * yRedQVec) * std::sqrt(amplRedQVec); + redQVec = std::hypot(redQVecXDaugSubtr, redQVecYDaugSubtr) * std::sqrt(amplRedQVec); } } if (storeRedQVec) { @@ -952,8 +953,11 @@ struct HfTaskFlowCharmHadrons { centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); if (storeResoOccu) { - const auto occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); - registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); + float occupancy{-999.f}; + if (occEstimator != 0) { + occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); + } const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); std::vector evtSelFlags = getEventSelectionFlags(rejectionMask); registry.fill(HIST("spReso/hSparseReso"), centrality, xQVecFT0c * xQVecFV0a + yQVecFT0c * yQVecFV0a, From 2419970a482650c8b1ab1f6d8b1dea5c058c656c Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 7 Apr 2026 09:12:39 +0000 Subject: [PATCH 06/12] Please consider the following formatting changes --- PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 0afa8d2927b..a218cedb276 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -724,7 +724,7 @@ struct HfTaskFlowCharmHadrons { bool subtractDaugsFromQVec = (qVecDetector == QvecEstimator::TPCNeg || qVecDetector == QvecEstimator::TPCPos || qVecDetector == QvecEstimator::TPCTot); - float scalprodCand{-999.f}; + float scalprodCand{-999.f}; if (subtractDaugsFromQVec) { std::vector tracksQx; From 45c7404b55d5fcdad0140dcb996c1d8c6598398c Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 7 Apr 2026 12:32:09 +0200 Subject: [PATCH 07/12] Dummy commit --- PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index a218cedb276..522162f4275 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -742,7 +742,7 @@ struct HfTaskFlowCharmHadrons { float yQVecDaugSubtr = yQVec - std::accumulate(tracksQy.begin(), tracksQy.end(), 0.0); scalprodCand = cosNPhi * xQVecDaugSubtr + sinNPhi * yQVecDaugSubtr; } else { - scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; + scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; } if (fillMassPtMlTree || fillMassPtMlSpCentTree) { From 2c61a23401450064aaecd34d6ff0bccbdb1c6b6f Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 7 Apr 2026 10:32:56 +0000 Subject: [PATCH 08/12] Please consider the following formatting changes --- PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 522162f4275..a218cedb276 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -742,7 +742,7 @@ struct HfTaskFlowCharmHadrons { float yQVecDaugSubtr = yQVec - std::accumulate(tracksQy.begin(), tracksQy.end(), 0.0); scalprodCand = cosNPhi * xQVecDaugSubtr + sinNPhi * yQVecDaugSubtr; } else { - scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; + scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; } if (fillMassPtMlTree || fillMassPtMlSpCentTree) { From c42d6c344ac95c4829d65d88a500bc24e2d89e44 Mon Sep 17 00:00:00 2001 From: Marcello Di Costanzo Date: Fri, 24 Apr 2026 14:40:53 +0200 Subject: [PATCH 09/12] Independent calibrations for EsE q-vectors --- Common/DataModel/Qvectors.h | 99 +++- Common/TableProducer/qVectorsTable.cxx | 660 +++++++++++++---------- PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 26 +- PWGHF/D2H/Utils/utilsFlow.h | 50 ++ 4 files changed, 518 insertions(+), 317 deletions(-) diff --git a/Common/DataModel/Qvectors.h b/Common/DataModel/Qvectors.h index 78bd0ef2181..71b4a496ff6 100644 --- a/Common/DataModel/Qvectors.h +++ b/Common/DataModel/Qvectors.h @@ -77,12 +77,6 @@ DECLARE_SOA_COLUMN(LabelsTPCpos, labelsTPCpos, std::vector); DECLARE_SOA_COLUMN(LabelsTPCneg, labelsTPCneg, std::vector); DECLARE_SOA_COLUMN(LabelsTPCall, labelsTPCall, std::vector); -// Normalized amplitudes for Event-shape Engineering -DECLARE_SOA_COLUMN(QVecRedFT0C, qVecRedFT0C, float); -DECLARE_SOA_COLUMN(QVecRedTpcPos, qVecRedTpcPos, float); -DECLARE_SOA_COLUMN(QVecRedTpcNeg, qVecRedTpcNeg, float); -DECLARE_SOA_COLUMN(QVecRedTpcAll, qVecRedTpcAll, float); - // Deprecated, will be removed in future after transition time // DECLARE_SOA_COLUMN(QvecBPosReVec, qvecBPosReVec, std::vector); DECLARE_SOA_COLUMN(QvecBPosImVec, qvecBPosImVec, std::vector); @@ -127,8 +121,6 @@ DECLARE_SOA_TABLE(QvectorTPCposVecs, "AOD", "QVECTORSTPCPVEC", qvec::IsCalibrate DECLARE_SOA_TABLE(QvectorTPCnegVecs, "AOD", "QVECTORSTPCNVEC", qvec::IsCalibrated, qvec::QvecTPCnegReVec, qvec::QvecTPCnegImVec, qvec::NTrkTPCneg, qvec::LabelsTPCneg); DECLARE_SOA_TABLE(QvectorTPCallVecs, "AOD", "QVECTORSTPCAVEC", qvec::IsCalibrated, qvec::QvecTPCallReVec, qvec::QvecTPCallImVec, qvec::NTrkTPCall, qvec::LabelsTPCall); -DECLARE_SOA_TABLE(QvectorsReds, "AOD", "QVECTORSRED", qvec::QVecRedFT0C, qvec::QVecRedTpcPos, qvec::QVecRedTpcNeg, qvec::QVecRedTpcAll); - using QvectorFT0C = QvectorFT0Cs::iterator; using QvectorFT0A = QvectorFT0As::iterator; using QvectorFT0M = QvectorFT0Ms::iterator; @@ -145,8 +137,6 @@ using QvectorTPCposVec = QvectorTPCposVecs::iterator; using QvectorTPCnegVec = QvectorTPCnegVecs::iterator; using QvectorTPCallVec = QvectorTPCallVecs::iterator; -using QvectorRed = QvectorsReds::iterator; - // Deprecated, will be removed in future after transition time // DECLARE_SOA_TABLE(QvectorBPoss, "AOD", "QVECTORSBPOS", qvec::IsCalibrated, qvec::QvecBPosRe, qvec::QvecBPosIm, qvec::NTrkBPos, qvec::LabelsBPos); DECLARE_SOA_TABLE(QvectorBNegs, "AOD", "QVECTORSBNEG", qvec::IsCalibrated, qvec::QvecBNegRe, qvec::QvecBNegIm, qvec::NTrkBNeg, qvec::LabelsBNeg); @@ -165,6 +155,95 @@ using QvectorBNegVec = QvectorBNegVecs::iterator; using QvectorBTotVec = QvectorBTotVecs::iterator; ///////////////////////////////////////////////////////////////// +namespace ese_qvec +{ +DECLARE_SOA_COLUMN(IsCalibrated, isCalibrated, bool); + +DECLARE_SOA_COLUMN(EseQvecRe, eseQvecRe, std::vector); +DECLARE_SOA_COLUMN(EseQvecIm, eseQvecIm, std::vector); +DECLARE_SOA_COLUMN(EseQvecAmp, eseQvecAmp, std::vector); + +DECLARE_SOA_COLUMN(EseQvecFT0CReVec, eseQvecFT0CReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0CImVec, eseQvecFT0CImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0AReVec, eseQvecFT0AReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0AImVec, eseQvecFT0AImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0MReVec, eseQvecFT0MReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0MImVec, eseQvecFT0MImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFV0AReVec, eseQvecFV0AReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFV0AImVec, eseQvecFV0AImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCposReVec, eseQvecTPCposReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCposImVec, eseQvecTPCposImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCnegReVec, eseQvecTPCnegReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCnegImVec, eseQvecTPCnegImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCallReVec, eseQvecTPCallReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCallImVec, eseQvecTPCallImVec, std::vector); + +DECLARE_SOA_COLUMN(EseQvecFT0CRe, eseQvecFT0CRe, float); +DECLARE_SOA_COLUMN(EseQvecFT0CIm, eseQvecFT0CIm, float); +DECLARE_SOA_COLUMN(EseQvecFT0ARe, eseQvecFT0ARe, float); +DECLARE_SOA_COLUMN(EseQvecFT0AIm, eseQvecFT0AIm, float); +DECLARE_SOA_COLUMN(EseQvecFT0MRe, eseQvecFT0MRe, float); +DECLARE_SOA_COLUMN(EseQvecFT0MIm, eseQvecFT0MIm, float); +DECLARE_SOA_COLUMN(EseQvecFV0ARe, eseQvecFV0ARe, float); +DECLARE_SOA_COLUMN(EseQvecFV0AIm, eseQvecFV0AIm, float); +DECLARE_SOA_COLUMN(EseQvecTPCposRe, eseQvecTPCposRe, float); +DECLARE_SOA_COLUMN(EseQvecTPCposIm, eseQvecTPCposIm, float); +DECLARE_SOA_COLUMN(EseQvecTPCnegRe, eseQvecTPCnegRe, float); +DECLARE_SOA_COLUMN(EseQvecTPCnegIm, eseQvecTPCnegIm, float); +DECLARE_SOA_COLUMN(EseQvecTPCallRe, eseQvecTPCallRe, float); +DECLARE_SOA_COLUMN(EseQvecTPCallIm, eseQvecTPCallIm, float); + +DECLARE_SOA_COLUMN(EseRedQFT0C, eseRedQFT0C, float); +DECLARE_SOA_COLUMN(EseRedQFT0A, eseRedQFT0A, float); +DECLARE_SOA_COLUMN(EseRedQFT0M, eseRedQFT0M, float); +DECLARE_SOA_COLUMN(EseRedQFV0A, eseRedQFV0A, float); +///////////////////////////////////////////////////////////////// +} // namespace ese_qvec + +DECLARE_SOA_TABLE(EseQvectors, "AOD", "ESEQVECTORDEVS", //! Table with all Qvectors. + qvec::Cent, ese_qvec::IsCalibrated, ese_qvec::EseQvecRe, ese_qvec::EseQvecIm, ese_qvec::EseQvecAmp); +using Qvector = Qvectors::iterator; + +DECLARE_SOA_TABLE(EseQvecFT0Cs, "AOD", "ESEQVECFT0C", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0CRe, ese_qvec::EseQvecFT0CIm, qvec::SumAmplFT0C); +DECLARE_SOA_TABLE(EseQvecFT0As, "AOD", "ESEQVECFT0A", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0ARe, ese_qvec::EseQvecFT0AIm, qvec::SumAmplFT0A); +DECLARE_SOA_TABLE(EseQvecFT0Ms, "AOD", "ESEQVECFT0M", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0MRe, ese_qvec::EseQvecFT0MIm, qvec::SumAmplFT0M); +DECLARE_SOA_TABLE(EseQvecFV0As, "AOD", "ESEQVECFV0A", ese_qvec::IsCalibrated, ese_qvec::EseQvecFV0ARe, ese_qvec::EseQvecFV0AIm, qvec::SumAmplFV0A); +DECLARE_SOA_TABLE(EseQvecTPCposs, "AOD", "ESEQVECTPCPOS", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCposRe, ese_qvec::EseQvecTPCposIm, qvec::NTrkTPCpos, qvec::LabelsTPCpos); +DECLARE_SOA_TABLE(EseQvecTPCnegs, "AOD", "ESEQVECTPCNEG", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCnegRe, ese_qvec::EseQvecTPCnegIm, qvec::NTrkTPCneg, qvec::LabelsTPCneg); +DECLARE_SOA_TABLE(EseQvecTPCalls, "AOD", "ESEQVECTPCALL", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCallRe, ese_qvec::EseQvecTPCallIm, qvec::NTrkTPCall, qvec::LabelsTPCall); + +DECLARE_SOA_TABLE(EseQvecFT0CVecs, "AOD", "ESEQVECFT0CVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0CReVec, ese_qvec::EseQvecFT0CImVec, qvec::SumAmplFT0C); +DECLARE_SOA_TABLE(EseQvecFT0AVecs, "AOD", "ESEQVECFT0AVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0AReVec, ese_qvec::EseQvecFT0AImVec, qvec::SumAmplFT0A); +DECLARE_SOA_TABLE(EseQvecFT0MVecs, "AOD", "ESEQVECFT0MVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0MReVec, ese_qvec::EseQvecFT0MImVec, qvec::SumAmplFT0M); +DECLARE_SOA_TABLE(EseQvecFV0AVecs, "AOD", "ESEQVECFV0AVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFV0AReVec, ese_qvec::EseQvecFV0AImVec, qvec::SumAmplFV0A); +DECLARE_SOA_TABLE(EseQvecTPCposVecs, "AOD", "ESEQVECTPCPVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCposReVec, ese_qvec::EseQvecTPCposImVec, qvec::NTrkTPCpos, qvec::LabelsTPCpos); +DECLARE_SOA_TABLE(EseQvecTPCnegVecs, "AOD", "ESEQVECTPCNVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCnegReVec, ese_qvec::EseQvecTPCnegImVec, qvec::NTrkTPCneg, qvec::LabelsTPCneg); +DECLARE_SOA_TABLE(EseQvecTPCallVecs, "AOD", "ESEQVECTPCAVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCallReVec, ese_qvec::EseQvecTPCallImVec, qvec::NTrkTPCall, qvec::LabelsTPCall); + +DECLARE_SOA_TABLE(EseQvecPercs, "AOD", "ESEQVECPERC", ese_qvec::EseQvecFT0CRe, ese_qvec::EseQvecFT0CIm, qvec::SumAmplFT0C, + ese_qvec::EseQvecFT0ARe, ese_qvec::EseQvecFT0AIm, qvec::SumAmplFT0A, + ese_qvec::EseQvecFT0MRe, ese_qvec::EseQvecFT0MIm, qvec::SumAmplFT0M, + ese_qvec::EseQvecFV0ARe, ese_qvec::EseQvecFV0AIm, qvec::SumAmplFV0A, + ese_qvec::EseQvecTPCposRe, ese_qvec::EseQvecTPCposIm, qvec::NTrkTPCpos, + ese_qvec::EseQvecTPCnegRe, ese_qvec::EseQvecTPCnegIm, qvec::NTrkTPCneg, + ese_qvec::EseQvecTPCallRe, ese_qvec::EseQvecTPCallIm, qvec::NTrkTPCall); + +using EseQvectorFT0C = EseQvecFT0Cs::iterator; +using EseQvectorFT0A = EseQvecFT0As::iterator; +using EseQvectorFT0M = EseQvecFT0Ms::iterator; +using EseQvectorFV0A = EseQvecFV0As::iterator; +using EseQvectorTPCpos = EseQvecTPCposs::iterator; +using EseQvectorTPCneg = EseQvecTPCnegs::iterator; +using EseQvectorTPCall = EseQvecTPCalls::iterator; + +using EseQvectorFT0CVec = EseQvecFT0CVecs::iterator; +using EseQvectorFT0AVec = EseQvecFT0AVecs::iterator; +using EseQvectorFT0MVec = EseQvecFT0MVecs::iterator; +using EseQvectorFV0AVec = EseQvecFV0AVecs::iterator; +using EseQvectorTPCposVec = EseQvecTPCposVecs::iterator; +using EseQvectorTPCnegVec = EseQvecTPCnegVecs::iterator; +using EseQvectorTPCallVec = EseQvecTPCallVecs::iterator; +///////////////////////////////////////////////////////////////// } // namespace o2::aod #endif // COMMON_DATAMODEL_QVECTORS_H_ diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index cd8ac007be0..06ee7c6321a 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -73,7 +73,15 @@ struct qVectorsTable { kFV0A, kTPCpos, kTPCneg, - kTPCall + kTPCall, + kNDetectors + }; + enum Corrections { + kNoCorr = 0, + kRecenter, + kTwist, + kRescale, + kNCorrections }; // Configurables. @@ -99,7 +107,8 @@ struct qVectorsTable { Configurable useCorrectionForRun{"useCorrectionForRun", true, "Get Qvector corrections based on run number instead of timestamp"}; Configurable cfgGainEqPath{"cfgGainEqPath", "Users/j/junlee/Qvector/GainEq", "CCDB path for gain equalization constants"}; - Configurable cfgQvecCalibPath{"cfgQvecCalibPath", "Analysis/EventPlane/QVecCorrections", "CCDB pasth for Q-vecteor calibration constants"}; + Configurable cfgQvecCalibPath{"cfgQvecCalibPath", "Analysis/EventPlane/QVecCorrections", "CCDB path for Q-vector calibration constants"}; + Configurable cfgQvecEseCalibPath{"cfgQvecEseCalibPath", "Analysis/EventPlane/QVecEseCorrections", "CCDB path for EsE Q-vector calibration constants"}; Configurable cfgShiftCorr{"cfgShiftCorr", false, "configurable flag for shift correction"}; Configurable cfgShiftPath{"cfgShiftPath", "", "CCDB path for shift correction"}; @@ -133,7 +142,23 @@ struct qVectorsTable { Produces qVectorTPCnegVec; Produces qVectorTPCallVec; - Produces qVectorRed; + Produces eseQVector; + Produces eseQVectorFT0C; + Produces eseQVectorFT0A; + Produces eseQVectorFT0M; + Produces eseQVectorFV0A; + Produces eseQVectorTPCpos; + Produces eseQVectorTPCneg; + Produces eseQVectorTPCall; + + Produces eseQVectorFT0CVec; + Produces eseQVectorFT0AVec; + Produces eseQVectorFT0MVec; + Produces eseQVectorFV0AVec; + Produces eseQVectorTPCposVec; + Produces eseQVectorTPCnegVec; + Produces eseQVectorTPCallVec; + Produces eseQVectorPerc; std::vector FT0RelGainConst{}; std::vector FV0RelGainConst{}; @@ -154,7 +179,8 @@ struct qVectorsTable { int runNumber{-1}; float cent; - std::vector objQvec{}; + std::vector corrsQvecSp{}; + std::vector corrsQvecEse{}; std::vector shiftprofile{}; // Deprecated, will be removed in future after transition time // @@ -261,19 +287,34 @@ struct qVectorsTable { LOGF(fatal, "Could not get the alignment parameters for FV0."); } - objQvec.clear(); + corrsQvecSp.clear(); for (std::size_t i = 0; i < cfgnMods->size(); i++) { int ind = cfgnMods->at(i); fullPath = cfgQvecCalibPath; fullPath += "/v"; fullPath += std::to_string(ind); - auto objqvec = getForTsOrRun(fullPath, timestamp, runnumber); - if (!objqvec) { + auto modeCorrQvecSp = getForTsOrRun(fullPath, timestamp, runnumber); + if (!modeCorrQvecSp) { fullPath = cfgQvecCalibPath; fullPath += "/v2"; - objqvec = getForTsOrRun(fullPath, timestamp, runnumber); + modeCorrQvecSp = getForTsOrRun(fullPath, timestamp, runnumber); + } + corrsQvecSp.push_back(modeCorrQvecSp); + } + + corrsQvecEse.clear(); + for (std::size_t i = 0; i < cfgnMods->size(); i++) { + int ind = cfgnMods->at(i); + fullPath = cfgQvecEseCalibPath; + fullPath += "/eseq"; + fullPath += std::to_string(ind); + auto modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); + if (!modeCorrQvecEse) { + fullPath = cfgQvecCalibPath; // cfgQvecEseCalibPath; + fullPath += "/v2"; // "/eseq2"; + modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); } - objQvec.push_back(objqvec); + corrsQvecEse.push_back(modeCorrQvecEse); } if (cfgShiftCorr) { @@ -351,16 +392,139 @@ struct qVectorsTable { } } + void NormalizeQvec(std::vector& QvecReNorm, std::vector& QvecImNorm, std::vector QvecReRaw, std::vector QvecImRaw, std::vector& QvecAmp, bool useSqrt = false) { + + for (std::size_t i = 0; i < kNDetectors; i++) { + float qVecDetReNorm{999.}, qVecDetImNorm{999.}; + if (QvecAmp[i] > 1e-8) { + if (useSqrt) { + qVecDetReNorm = QvecReRaw[i] / std::sqrt(QvecAmp[i]); + qVecDetImNorm = QvecImRaw[i] / std::sqrt(QvecAmp[i]); + } else { + qVecDetReNorm = QvecReRaw[i] / QvecAmp[i]; + qVecDetImNorm = QvecImRaw[i] / QvecAmp[i]; + } + } + for (int iCorr=0; iCorr < Corrections::kNCorrections; iCorr++) { + QvecReNorm.push_back(qVecDetReNorm); + QvecImNorm.push_back(qVecDetImNorm); + } + } + } + + void CorrectQvec(float cent, std::vector& qvecRe, std::vector& qvecIm, TH3F* histsCorrs, int nMode) { + int nCorrections = static_cast(Corrections::kNCorrections); + if (cent < cfgMaxCentrality) { + for (auto i{0u}; i < kTPCall + 1; i++) { + int idxDet = i * Corrections::kNCorrections; + helperEP.DoRecenter(qvecRe[idxDet + Corrections::kRecenter], qvecIm[idxDet + Corrections::kRecenter], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + + helperEP.DoRecenter(qvecRe[idxDet + Corrections::kTwist], qvecIm[idxDet + Corrections::kTwist], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + helperEP.DoTwist(qvecRe[idxDet + Corrections::kTwist], qvecIm[idxDet + Corrections::kTwist], + histsCorrs->GetBinContent(static_cast(cent) + 1, 3, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 4, i + 1)); + + helperEP.DoRecenter(qvecRe[idxDet + Corrections::kRescale], qvecIm[idxDet + Corrections::kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + helperEP.DoTwist(qvecRe[idxDet + Corrections::kRescale], qvecIm[idxDet + Corrections::kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 3, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 4, i + 1)); + helperEP.DoRescale(qvecRe[idxDet + Corrections::kRescale], qvecIm[idxDet + Corrections::kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 5, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 6, i + 1)); + } + if (cfgShiftCorr) { + auto deltapsiFT0C = 0.0; + auto deltapsiFT0A = 0.0; + auto deltapsiFT0M = 0.0; + auto deltapsiFV0A = 0.0; + auto deltapsiTPCpos = 0.0; + auto deltapsiTPCneg = 0.0; + auto deltapsiTPCall = 0.0; + + auto psidefFT0C = TMath::ATan2(qvecIm[kFT0C * nCorrections + Corrections::kRescale], qvecRe[kFT0C * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefFT0A = TMath::ATan2(qvecIm[kFT0A * nCorrections + Corrections::kRescale], qvecRe[kFT0A * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefFT0M = TMath::ATan2(qvecIm[kFT0M * nCorrections + Corrections::kRescale], qvecRe[kFT0M * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefFV0A = TMath::ATan2(qvecIm[kFV0A * nCorrections + Corrections::kRescale], qvecRe[kFV0A * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefTPCpos = TMath::ATan2(qvecIm[kTPCpos * nCorrections + Corrections::kRescale], qvecRe[kTPCpos * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefTPCneg = TMath::ATan2(qvecIm[kTPCneg * nCorrections + Corrections::kRescale], qvecRe[kTPCneg * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefTPCall = TMath::ATan2(qvecIm[kTPCall * nCorrections + Corrections::kRescale], qvecRe[kTPCall * nCorrections + Corrections::kRescale]) / static_cast(nMode); + + for (int ishift = 1; ishift <= 10; ishift++) { + auto coeffshiftxFT0C = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFT0C, ishift - 0.5)); + auto coeffshiftyFT0C = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFT0C + 1, ishift - 0.5)); + auto coeffshiftxFT0A = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFT0A, ishift - 0.5)); + auto coeffshiftyFT0A = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFT0A + 1, ishift - 0.5)); + auto coeffshiftxFT0M = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFT0M, ishift - 0.5)); + auto coeffshiftyFT0M = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFT0M + 1, ishift - 0.5)); + auto coeffshiftxFV0A = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFV0A, ishift - 0.5)); + auto coeffshiftyFV0A = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kFV0A + 1, ishift - 0.5)); + auto coeffshiftxTPCpos = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kTPCpos, ishift - 0.5)); + auto coeffshiftyTPCpos = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kTPCpos + 1, ishift - 0.5)); + auto coeffshiftxTPCneg = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kTPCneg, ishift - 0.5)); + auto coeffshiftyTPCneg = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kTPCneg + 1, ishift - 0.5)); + auto coeffshiftxTPCall = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kTPCall, ishift - 0.5)); + auto coeffshiftyTPCall = shiftprofile.at(nMode - 2)->GetBinContent(shiftprofile.at(nMode - 2)->FindBin(cent, 2 * kTPCall + 1, ishift - 0.5)); + + deltapsiFT0C += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0C * TMath::Cos(ishift * static_cast(nMode) * psidefFT0C) + coeffshiftyFT0C * TMath::Sin(ishift * static_cast(nMode) * psidefFT0C))) / static_cast(nMode); + deltapsiFT0A += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0A * TMath::Cos(ishift * static_cast(nMode) * psidefFT0A) + coeffshiftyFT0A * TMath::Sin(ishift * static_cast(nMode) * psidefFT0A))) / static_cast(nMode); + deltapsiFT0M += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0M * TMath::Cos(ishift * static_cast(nMode) * psidefFT0M) + coeffshiftyFT0M * TMath::Sin(ishift * static_cast(nMode) * psidefFT0M))) / static_cast(nMode); + deltapsiFV0A += ((2. / (1.0 * ishift)) * (-coeffshiftxFV0A * TMath::Cos(ishift * static_cast(nMode) * psidefFV0A) + coeffshiftyFV0A * TMath::Sin(ishift * static_cast(nMode) * psidefFV0A))) / static_cast(nMode); + deltapsiTPCpos += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCpos * TMath::Cos(ishift * static_cast(nMode) * psidefTPCpos) + coeffshiftyTPCpos * TMath::Sin(ishift * static_cast(nMode) * psidefTPCpos))) / static_cast(nMode); + deltapsiTPCneg += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCneg * TMath::Cos(ishift * static_cast(nMode) * psidefTPCneg) + coeffshiftyTPCneg * TMath::Sin(ishift * static_cast(nMode) * psidefTPCneg))) / static_cast(nMode); + deltapsiTPCall += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCall * TMath::Cos(ishift * static_cast(nMode) * psidefTPCall) + coeffshiftyTPCall * TMath::Sin(ishift * static_cast(nMode) * psidefTPCall))) / static_cast(nMode); + } + + deltapsiFT0C *= static_cast(nMode); + deltapsiFT0A *= static_cast(nMode); + deltapsiFT0M *= static_cast(nMode); + deltapsiFV0A *= static_cast(nMode); + deltapsiTPCpos *= static_cast(nMode); + deltapsiTPCneg *= static_cast(nMode); + deltapsiTPCall *= static_cast(nMode); + + float qvecReShiftedFT0C = qvecRe[kFT0C * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0C) - qvecIm[kFT0C * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0C); + float qvecImShiftedFT0C = qvecRe[kFT0C * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0C) + qvecIm[kFT0C * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0C); + float qvecReShiftedFT0A = qvecRe[kFT0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0A) - qvecIm[kFT0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0A); + float qvecImShiftedFT0A = qvecRe[kFT0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0A) + qvecIm[kFT0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0A); + float qvecReShiftedFT0M = qvecRe[kFT0M * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0M) - qvecIm[kFT0M * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0M); + float qvecImShiftedFT0M = qvecRe[kFT0M * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0M) + qvecIm[kFT0M * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0M); + float qvecReShiftedFV0A = qvecRe[kFV0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFV0A) - qvecIm[kFV0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFV0A); + float qvecImShiftedFV0A = qvecRe[kFV0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFV0A) + qvecIm[kFV0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFV0A); + float qvecReShiftedTPCpos = qvecRe[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCpos) - qvecIm[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCpos); + float qvecImShiftedTPCpos = qvecRe[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCpos) + qvecIm[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCpos); + float qvecReShiftedTPCneg = qvecRe[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCneg) - qvecIm[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCneg); + float qvecImShiftedTPCneg = qvecRe[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCneg) + qvecIm[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCneg); + float qvecReShiftedTPCall = qvecRe[kTPCall * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCall) - qvecIm[kTPCall * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCall); + float qvecImShiftedTPCall = qvecRe[kTPCall * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCall) + qvecIm[kTPCall * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCall); + + qvecRe[kFT0C * nCorrections + Corrections::kRescale] = qvecReShiftedFT0C; + qvecIm[kFT0C * nCorrections + Corrections::kRescale] = qvecImShiftedFT0C; + qvecRe[kFT0A * nCorrections + Corrections::kRescale] = qvecReShiftedFT0A; + qvecIm[kFT0A * nCorrections + Corrections::kRescale] = qvecImShiftedFT0A; + qvecRe[kFT0M * nCorrections + Corrections::kRescale] = qvecReShiftedFT0M; + qvecIm[kFT0M * nCorrections + Corrections::kRescale] = qvecImShiftedFT0M; + qvecRe[kFV0A * nCorrections + Corrections::kRescale] = qvecReShiftedFV0A; + qvecIm[kFV0A * nCorrections + Corrections::kRescale] = qvecImShiftedFV0A; + qvecRe[kTPCpos * nCorrections + Corrections::kRescale] = qvecReShiftedTPCpos; + qvecIm[kTPCpos * nCorrections + Corrections::kRescale] = qvecImShiftedTPCpos; + qvecRe[kTPCneg * nCorrections + Corrections::kRescale] = qvecReShiftedTPCneg; + qvecIm[kTPCneg * nCorrections + Corrections::kRescale] = qvecImShiftedTPCneg; + qvecRe[kTPCall * nCorrections + Corrections::kRescale] = qvecReShiftedTPCall; + qvecIm[kTPCall * nCorrections + Corrections::kRescale] = qvecImShiftedTPCall; + } + } + } + template - void CalQvec(const Nmode nmode, const CollType& coll, const TrackType& track, std::vector& QvecRe, std::vector& QvecIm, std::vector& QvecAmp, std::vector& TrkTPCposLabel, std::vector& TrkTPCnegLabel, std::vector& TrkTPCallLabel) + void CalQvec(const Nmode nMode, const CollType& coll, const TrackType& track, std::vector& QvecRe, std::vector& QvecIm, std::vector& QvecAmp, std::vector& TrkTPCposLabel, std::vector& TrkTPCnegLabel, std::vector& TrkTPCallLabel) { - float qVectFT0A[2] = {0.}; - float qVectFT0C[2] = {0.}; - float qVectFT0M[2] = {0.}; - float qVectFV0A[2] = {0.}; - float qVectTPCpos[2] = {0.}; - float qVectTPCneg[2] = {0.}; - float qVectTPCall[2] = {0.}; + float qVectFT0A[2] = {-999., -999.}; + float qVectFT0C[2] = {-999., -999.}; + float qVectFT0M[2] = {-999., -999.}; + float qVectFV0A[2] = {-999., -999.}; + float qVectTPCpos[2] = {0., 0.}; // Always computed + float qVectTPCneg[2] = {0., 0.}; // Always computed + float qVectTPCall[2] = {0., 0.}; // Always computed TComplex QvecDet(0); TComplex QvecFT0M(0); @@ -380,17 +544,14 @@ struct qVectorsTable { histosQA.fill(HIST("FT0Amp"), ampl, FT0AchId); histosQA.fill(HIST("FT0AmpCor"), ampl / FT0RelGainConst[FT0AchId], FT0AchId); - helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nmode, QvecDet, sumAmplFT0A, ft0geom, fv0geom); - helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nmode, QvecFT0M, sumAmplFT0M, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nMode, QvecDet, sumAmplFT0A, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nMode, QvecFT0M, sumAmplFT0M, ft0geom, fv0geom); } if (sumAmplFT0A > 1e-8) { QvecDet /= sumAmplFT0A; qVectFT0A[0] = QvecDet.Re(); qVectFT0A[1] = QvecDet.Im(); } - } else { - qVectFT0A[0] = 999.; - qVectFT0A[1] = 999.; } if (useDetector["QvectorFT0Cs"]) { @@ -402,65 +563,42 @@ struct qVectorsTable { histosQA.fill(HIST("FT0Amp"), ampl, FT0CchId); histosQA.fill(HIST("FT0AmpCor"), ampl / FT0RelGainConst[FT0CchId], FT0CchId); - helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nmode, QvecDet, sumAmplFT0C, ft0geom, fv0geom); - helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nmode, QvecFT0M, sumAmplFT0M, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nMode, QvecDet, sumAmplFT0C, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nMode, QvecFT0M, sumAmplFT0M, ft0geom, fv0geom); } if (sumAmplFT0C > 1e-8) { QvecDet /= sumAmplFT0C; qVectFT0C[0] = QvecDet.Re(); qVectFT0C[1] = QvecDet.Im(); - } else { - qVectFT0C[0] = 999.; - qVectFT0C[1] = 999.; } - } else { - qVectFT0C[0] = -999.; - qVectFT0C[1] = -999.; - } - - if (sumAmplFT0M > 1e-8 && useDetector["QvectorFT0Ms"]) { - QvecFT0M /= sumAmplFT0M; - qVectFT0M[0] = QvecFT0M.Re(); - qVectFT0M[1] = QvecFT0M.Im(); - } else { - qVectFT0M[0] = 999.; - qVectFT0M[1] = 999.; + if (sumAmplFT0M > 1e-8 && useDetector["QvectorFT0Ms"]) { + QvecFT0M /= sumAmplFT0M; + qVectFT0M[0] = QvecFT0M.Re(); + qVectFT0M[1] = QvecFT0M.Im(); + } } - } else { - qVectFT0A[0] = -999.; - qVectFT0A[1] = -999.; - qVectFT0C[0] = -999.; - qVectFT0C[1] = -999.; - qVectFT0M[0] = -999.; - qVectFT0M[1] = -999.; - } - QvecDet = TComplex(0., 0.); - sumAmplFV0A = 0; - if (coll.has_foundFV0() && useDetector["QvectorFV0As"]) { - auto fv0 = coll.foundFV0(); + QvecDet = TComplex(0., 0.); + sumAmplFV0A = 0; + if (coll.has_foundFV0() && useDetector["QvectorFV0As"]) { + auto fv0 = coll.foundFV0(); - for (std::size_t iCh = 0; iCh < fv0.channel().size(); iCh++) { - float ampl = fv0.amplitude()[iCh]; - int FV0AchId = fv0.channel()[iCh]; - histosQA.fill(HIST("FV0Amp"), ampl, FV0AchId); - histosQA.fill(HIST("FV0AmpCor"), ampl / FV0RelGainConst[FV0AchId], FV0AchId); + for (std::size_t iCh = 0; iCh < fv0.channel().size(); iCh++) { + float ampl = fv0.amplitude()[iCh]; + int FV0AchId = fv0.channel()[iCh]; + histosQA.fill(HIST("FV0Amp"), ampl, FV0AchId); + histosQA.fill(HIST("FV0AmpCor"), ampl / FV0RelGainConst[FV0AchId], FV0AchId); - helperEP.SumQvectors(1, FV0AchId, ampl / FV0RelGainConst[FV0AchId], nmode, QvecDet, sumAmplFV0A, ft0geom, fv0geom); - } + helperEP.SumQvectors(1, FV0AchId, ampl / FV0RelGainConst[FV0AchId], nMode, QvecDet, sumAmplFV0A, ft0geom, fv0geom); + } - if (sumAmplFV0A > 1e-8) { - QvecDet /= sumAmplFV0A; - qVectFV0A[0] = QvecDet.Re(); - qVectFV0A[1] = QvecDet.Im(); - } else { - qVectFV0A[0] = 999.; - qVectFV0A[1] = 999.; + if (sumAmplFV0A > 1e-8) { + QvecDet /= sumAmplFV0A; + qVectFV0A[0] = QvecDet.Re(); + qVectFV0A[1] = QvecDet.Im(); + } } - } else { - qVectFV0A[0] = -999.; - qVectFV0A[1] = -999.; } int nTrkTPCpos = 0; @@ -478,77 +616,40 @@ struct qVectorsTable { if (trk.eta() < cfgEtaMin) { continue; } - qVectTPCall[0] += trk.pt() * std::cos(trk.phi() * nmode); - qVectTPCall[1] += trk.pt() * std::sin(trk.phi() * nmode); + qVectTPCall[0] += trk.pt() * std::cos(trk.phi() * nMode); + qVectTPCall[1] += trk.pt() * std::sin(trk.phi() * nMode); TrkTPCallLabel.push_back(trk.globalIndex()); nTrkTPCall++; if (std::abs(trk.eta()) < 0.1) { continue; } if (trk.eta() > 0 && (useDetector["QvectorTPCposs"] || useDetector["QvectorBPoss"])) { - qVectTPCpos[0] += trk.pt() * std::cos(trk.phi() * nmode); - qVectTPCpos[1] += trk.pt() * std::sin(trk.phi() * nmode); + qVectTPCpos[0] += trk.pt() * std::cos(trk.phi() * nMode); + qVectTPCpos[1] += trk.pt() * std::sin(trk.phi() * nMode); TrkTPCposLabel.push_back(trk.globalIndex()); nTrkTPCpos++; } else if (trk.eta() < 0 && (useDetector["QvectorTPCnegs"] || useDetector["QvectorBNegs"])) { - qVectTPCneg[0] += trk.pt() * std::cos(trk.phi() * nmode); - qVectTPCneg[1] += trk.pt() * std::sin(trk.phi() * nmode); + qVectTPCneg[0] += trk.pt() * std::cos(trk.phi() * nMode); + qVectTPCneg[1] += trk.pt() * std::sin(trk.phi() * nMode); TrkTPCnegLabel.push_back(trk.globalIndex()); nTrkTPCneg++; } } - if (nTrkTPCpos > 0) { - qVectTPCpos[0] /= nTrkTPCpos; - qVectTPCpos[1] /= nTrkTPCpos; - } else { - qVectTPCpos[0] = 999.; - qVectTPCpos[1] = 999.; - } - - if (nTrkTPCneg > 0) { - qVectTPCneg[0] /= nTrkTPCneg; - qVectTPCneg[1] /= nTrkTPCneg; - } else { - qVectTPCneg[0] = 999.; - qVectTPCneg[1] = 999.; - } - - if (nTrkTPCall > 0) { - qVectTPCall[0] /= nTrkTPCall; - qVectTPCall[1] /= nTrkTPCall; - } else { - qVectTPCall[0] = 999.; - qVectTPCall[1] = 999.; - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFT0C[0]); - QvecIm.push_back(qVectFT0C[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFT0A[0]); - QvecIm.push_back(qVectFT0A[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFT0M[0]); - QvecIm.push_back(qVectFT0M[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFV0A[0]); - QvecIm.push_back(qVectFV0A[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectTPCpos[0]); - QvecIm.push_back(qVectTPCpos[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectTPCneg[0]); - QvecIm.push_back(qVectTPCneg[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectTPCall[0]); - QvecIm.push_back(qVectTPCall[1]); - } + QvecRe.push_back(qVectFT0C[0]); + QvecIm.push_back(qVectFT0C[1]); + QvecRe.push_back(qVectFT0A[0]); + QvecIm.push_back(qVectFT0A[1]); + QvecRe.push_back(qVectFT0M[0]); + QvecIm.push_back(qVectFT0M[1]); + QvecRe.push_back(qVectFV0A[0]); + QvecIm.push_back(qVectFV0A[1]); + QvecRe.push_back(qVectTPCpos[0]); + QvecIm.push_back(qVectTPCpos[1]); + QvecRe.push_back(qVectTPCneg[0]); + QvecIm.push_back(qVectTPCneg[1]); + QvecRe.push_back(qVectTPCall[0]); + QvecIm.push_back(qVectTPCall[1]); QvecAmp.push_back(sumAmplFT0C); QvecAmp.push_back(sumAmplFT0A); @@ -564,24 +665,41 @@ struct qVectorsTable { std::vector TrkTPCposLabel{}; std::vector TrkTPCnegLabel{}; std::vector TrkTPCallLabel{}; - std::vector qvecRe{}; - std::vector qvecIm{}; std::vector qvecAmp{}; - std::vector qvecReFT0C{}; - std::vector qvecImFT0C{}; - std::vector qvecReFT0A{}; - std::vector qvecImFT0A{}; - std::vector qvecReFT0M{}; - std::vector qvecImFT0M{}; - std::vector qvecReFV0A{}; - std::vector qvecImFV0A{}; - std::vector qvecReTPCpos{}; - std::vector qvecImTPCpos{}; - std::vector qvecReTPCneg{}; - std::vector qvecImTPCneg{}; - std::vector qvecReTPCall{}; - std::vector qvecImTPCall{}; + std::vector qvecReSp{}; + std::vector qvecImSp{}; + std::vector qvecReFT0CSp{}; + std::vector qvecImFT0CSp{}; + std::vector qvecReFT0ASp{}; + std::vector qvecImFT0ASp{}; + std::vector qvecReFT0MSp{}; + std::vector qvecImFT0MSp{}; + std::vector qvecReFV0ASp{}; + std::vector qvecImFV0ASp{}; + std::vector qvecReTPCposSp{}; + std::vector qvecImTPCposSp{}; + std::vector qvecReTPCnegSp{}; + std::vector qvecImTPCnegSp{}; + std::vector qvecReTPCallSp{}; + std::vector qvecImTPCallSp{}; + + std::vector qvecReEse{}; + std::vector qvecImEse{}; + std::vector qvecReFT0CEse{}; + std::vector qvecImFT0CEse{}; + std::vector qvecReFT0AEse{}; + std::vector qvecImFT0AEse{}; + std::vector qvecReFT0MEse{}; + std::vector qvecImFT0MEse{}; + std::vector qvecReFV0AEse{}; + std::vector qvecImFV0AEse{}; + std::vector qvecReTPCposEse{}; + std::vector qvecImTPCposEse{}; + std::vector qvecReTPCnegEse{}; + std::vector qvecImTPCnegEse{}; + std::vector qvecReTPCallEse{}; + std::vector qvecImTPCallEse{}; auto bc = coll.bc_as(); int currentRun = bc.runNumber(); @@ -599,182 +717,132 @@ struct qVectorsTable { cent = 110.; IsCalibrated = false; } - for (std::size_t id = 0; id < cfgnMods->size(); id++) { - int nmode = cfgnMods->at(id); - CalQvec(nmode, coll, tracks, qvecRe, qvecIm, qvecAmp, TrkTPCposLabel, TrkTPCnegLabel, TrkTPCallLabel); - if (cent < cfgMaxCentrality) { - for (auto i{0u}; i < kTPCall + 1; i++) { - helperEP.DoRecenter(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 1], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 1], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 1, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 2, i + 1)); - - helperEP.DoRecenter(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 2], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 2], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 1, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 2, i + 1)); - helperEP.DoTwist(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 2], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 2], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 3, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 4, i + 1)); - - helperEP.DoRecenter(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 3], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 3], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 1, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 2, i + 1)); - helperEP.DoTwist(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 3], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 3], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 3, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 4, i + 1)); - helperEP.DoRescale(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 3], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 3], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 5, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 6, i + 1)); - } - if (cfgShiftCorr) { - auto deltapsiFT0C = 0.0; - auto deltapsiFT0A = 0.0; - auto deltapsiFT0M = 0.0; - auto deltapsiFV0A = 0.0; - auto deltapsiTPCpos = 0.0; - auto deltapsiTPCneg = 0.0; - auto deltapsiTPCall = 0.0; - - auto psidefFT0C = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3]) / static_cast(nmode); - auto psidefFT0A = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3]) / static_cast(nmode); - auto psidefFT0M = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3]) / static_cast(nmode); - auto psidefFV0A = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3]) / static_cast(nmode); - auto psidefTPCpos = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3]) / static_cast(nmode); - auto psidefTPCneg = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3]) / static_cast(nmode); - auto psidefTPCall = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3]) / static_cast(nmode); - - for (int ishift = 1; ishift <= 10; ishift++) { - auto coeffshiftxFT0C = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0C, ishift - 0.5)); - auto coeffshiftyFT0C = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0C + 1, ishift - 0.5)); - auto coeffshiftxFT0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0A, ishift - 0.5)); - auto coeffshiftyFT0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0A + 1, ishift - 0.5)); - auto coeffshiftxFT0M = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0M, ishift - 0.5)); - auto coeffshiftyFT0M = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0M + 1, ishift - 0.5)); - auto coeffshiftxFV0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFV0A, ishift - 0.5)); - auto coeffshiftyFV0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFV0A + 1, ishift - 0.5)); - auto coeffshiftxTPCpos = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCpos, ishift - 0.5)); - auto coeffshiftyTPCpos = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCpos + 1, ishift - 0.5)); - auto coeffshiftxTPCneg = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCneg, ishift - 0.5)); - auto coeffshiftyTPCneg = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCneg + 1, ishift - 0.5)); - auto coeffshiftxTPCall = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCall, ishift - 0.5)); - auto coeffshiftyTPCall = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCall + 1, ishift - 0.5)); - - deltapsiFT0C += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0C * TMath::Cos(ishift * static_cast(nmode) * psidefFT0C) + coeffshiftyFT0C * TMath::Sin(ishift * static_cast(nmode) * psidefFT0C))) / static_cast(nmode); - deltapsiFT0A += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0A * TMath::Cos(ishift * static_cast(nmode) * psidefFT0A) + coeffshiftyFT0A * TMath::Sin(ishift * static_cast(nmode) * psidefFT0A))) / static_cast(nmode); - deltapsiFT0M += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0M * TMath::Cos(ishift * static_cast(nmode) * psidefFT0M) + coeffshiftyFT0M * TMath::Sin(ishift * static_cast(nmode) * psidefFT0M))) / static_cast(nmode); - deltapsiFV0A += ((2. / (1.0 * ishift)) * (-coeffshiftxFV0A * TMath::Cos(ishift * static_cast(nmode) * psidefFV0A) + coeffshiftyFV0A * TMath::Sin(ishift * static_cast(nmode) * psidefFV0A))) / static_cast(nmode); - deltapsiTPCpos += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCpos * TMath::Cos(ishift * static_cast(nmode) * psidefTPCpos) + coeffshiftyTPCpos * TMath::Sin(ishift * static_cast(nmode) * psidefTPCpos))) / static_cast(nmode); - deltapsiTPCneg += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCneg * TMath::Cos(ishift * static_cast(nmode) * psidefTPCneg) + coeffshiftyTPCneg * TMath::Sin(ishift * static_cast(nmode) * psidefTPCneg))) / static_cast(nmode); - deltapsiTPCall += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCall * TMath::Cos(ishift * static_cast(nmode) * psidefTPCall) + coeffshiftyTPCall * TMath::Sin(ishift * static_cast(nmode) * psidefTPCall))) / static_cast(nmode); - } - deltapsiFT0C *= static_cast(nmode); - deltapsiFT0A *= static_cast(nmode); - deltapsiFT0M *= static_cast(nmode); - deltapsiFV0A *= static_cast(nmode); - deltapsiTPCpos *= static_cast(nmode); - deltapsiTPCneg *= static_cast(nmode); - deltapsiTPCall *= static_cast(nmode); - - float qvecReShiftedFT0C = qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Cos(deltapsiFT0C) - qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Sin(deltapsiFT0C); - float qvecImShiftedFT0C = qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Sin(deltapsiFT0C) + qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Cos(deltapsiFT0C); - float qvecReShiftedFT0A = qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Cos(deltapsiFT0A) - qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Sin(deltapsiFT0A); - float qvecImShiftedFT0A = qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Sin(deltapsiFT0A) + qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Cos(deltapsiFT0A); - float qvecReShiftedFT0M = qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Cos(deltapsiFT0M) - qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Sin(deltapsiFT0M); - float qvecImShiftedFT0M = qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Sin(deltapsiFT0M) + qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Cos(deltapsiFT0M); - float qvecReShiftedFV0A = qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Cos(deltapsiFV0A) - qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Sin(deltapsiFV0A); - float qvecImShiftedFV0A = qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Sin(deltapsiFV0A) + qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Cos(deltapsiFV0A); - float qvecReShiftedTPCpos = qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Cos(deltapsiTPCpos) - qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Sin(deltapsiTPCpos); - float qvecImShiftedTPCpos = qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Sin(deltapsiTPCpos) + qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Cos(deltapsiTPCpos); - float qvecReShiftedTPCneg = qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Cos(deltapsiTPCneg) - qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Sin(deltapsiTPCneg); - float qvecImShiftedTPCneg = qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Sin(deltapsiTPCneg) + qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Cos(deltapsiTPCneg); - float qvecReShiftedTPCall = qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Cos(deltapsiTPCall) - qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Sin(deltapsiTPCall); - float qvecImShiftedTPCall = qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Sin(deltapsiTPCall) + qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Cos(deltapsiTPCall); - - qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] = qvecReShiftedFT0C; - qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] = qvecImShiftedFT0C; - qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] = qvecReShiftedFT0A; - qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] = qvecImShiftedFT0A; - qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] = qvecReShiftedFT0M; - qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] = qvecImShiftedFT0M; - qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] = qvecReShiftedFV0A; - qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] = qvecImShiftedFV0A; - qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] = qvecReShiftedTPCpos; - qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] = qvecImShiftedTPCpos; - qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] = qvecReShiftedTPCneg; - qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] = qvecImShiftedTPCneg; - qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] = qvecReShiftedTPCall; - qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] = qvecImShiftedTPCall; - } - } + for (std::size_t id = 0; id < cfgnMods->size(); id++) { + int nMode = cfgnMods->at(id); + + // Raw Q-vectors, no multiplicity normalization and no corrections + std::vector qvecReRaw{}; + std::vector qvecImRaw{}; + CalQvec(nMode, coll, tracks, qvecReRaw, qvecImRaw, qvecAmp, TrkTPCposLabel, TrkTPCnegLabel, TrkTPCallLabel); + + // Scalar Product Q-vectors, normalization by multiplicity/amplitude + std::vector nModeQvecReSp{}; + std::vector nModeQvecImSp{}; + NormalizeQvec(nModeQvecReSp, nModeQvecImSp, qvecReRaw, qvecImRaw, qvecAmp, false); + CorrectQvec(cent, nModeQvecReSp, nModeQvecImSp, corrsQvecSp[id], nMode); + + // Add to summary vector + qvecReSp.insert(qvecReSp.end(), nModeQvecReSp.begin(), nModeQvecReSp.end()); + qvecImSp.insert(qvecImSp.end(), nModeQvecImSp.begin(), nModeQvecImSp.end()); + + // Ese Q-vectors, normalization by sqrt(multiplicity/amplitude) + std::vector nModeQvecReEse{}; + std::vector nModeQvecImEse{}; + NormalizeQvec(nModeQvecReEse, nModeQvecImEse, qvecReRaw, qvecImRaw, qvecAmp, true); + CorrectQvec(cent, nModeQvecReEse, nModeQvecImEse, corrsQvecEse[id], nMode); + + // Add to summary vector + qvecReEse.insert(qvecReEse.end(), nModeQvecReEse.begin(), nModeQvecReEse.end()); + qvecImEse.insert(qvecImEse.end(), nModeQvecImEse.begin(), nModeQvecImEse.end()); + + // Pick the desired correction level for the Q-vectors to be stored in the analysis table int CorrLevel = cfgCorrLevel == 0 ? 0 : cfgCorrLevel - 1; - qvecReFT0C.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + CorrLevel]); - qvecImFT0C.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + CorrLevel]); - qvecReFT0A.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + CorrLevel]); - qvecImFT0A.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + CorrLevel]); - qvecReFT0M.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + CorrLevel]); - qvecImFT0M.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + CorrLevel]); - qvecReFV0A.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + CorrLevel]); - qvecImFV0A.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + CorrLevel]); - qvecReTPCpos.push_back(qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + CorrLevel]); - qvecImTPCpos.push_back(qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + CorrLevel]); - qvecReTPCneg.push_back(qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + CorrLevel]); - qvecImTPCneg.push_back(qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + CorrLevel]); - qvecReTPCall.push_back(qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + CorrLevel]); - qvecImTPCall.push_back(qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + CorrLevel]); + int nCorrections = static_cast(Corrections::kNCorrections); + + qvecReFT0CSp.push_back(nModeQvecReSp[kFT0C * nCorrections + CorrLevel]); + qvecImFT0CSp.push_back(nModeQvecImSp[kFT0C * nCorrections + CorrLevel]); + qvecReFT0ASp.push_back(nModeQvecReSp[kFT0A * nCorrections + CorrLevel]); + qvecImFT0ASp.push_back(nModeQvecImSp[kFT0A * nCorrections + CorrLevel]); + qvecReFT0MSp.push_back(nModeQvecReSp[kFT0M * nCorrections + CorrLevel]); + qvecImFT0MSp.push_back(nModeQvecImSp[kFT0M * nCorrections + CorrLevel]); + qvecReFV0ASp.push_back(nModeQvecReSp[kFV0A * nCorrections + CorrLevel]); + qvecImFV0ASp.push_back(nModeQvecImSp[kFV0A * nCorrections + CorrLevel]); + qvecReTPCposSp.push_back(nModeQvecReSp[kTPCpos * nCorrections + CorrLevel]); + qvecImTPCposSp.push_back(nModeQvecImSp[kTPCpos * nCorrections + CorrLevel]); + qvecReTPCnegSp.push_back(nModeQvecReSp[kTPCneg * nCorrections + CorrLevel]); + qvecImTPCnegSp.push_back(nModeQvecImSp[kTPCneg * nCorrections + CorrLevel]); + qvecReTPCallSp.push_back(nModeQvecReSp[kTPCall * nCorrections + CorrLevel]); + qvecImTPCallSp.push_back(nModeQvecImSp[kTPCall * nCorrections + CorrLevel]); + + qvecReFT0CEse.push_back(nModeQvecReEse[kFT0C * nCorrections + CorrLevel]); + qvecImFT0CEse.push_back(nModeQvecImEse[kFT0C * nCorrections + CorrLevel]); + qvecReFT0AEse.push_back(nModeQvecReEse[kFT0A * nCorrections + CorrLevel]); + qvecImFT0AEse.push_back(nModeQvecImEse[kFT0A * nCorrections + CorrLevel]); + qvecReFT0MEse.push_back(nModeQvecReEse[kFT0M * nCorrections + CorrLevel]); + qvecImFT0MEse.push_back(nModeQvecImEse[kFT0M * nCorrections + CorrLevel]); + qvecReFV0AEse.push_back(nModeQvecReEse[kFV0A * nCorrections + CorrLevel]); + qvecImFV0AEse.push_back(nModeQvecImEse[kFV0A * nCorrections + CorrLevel]); + qvecReTPCposEse.push_back(nModeQvecReEse[kTPCpos * nCorrections + CorrLevel]); + qvecImTPCposEse.push_back(nModeQvecImEse[kTPCpos * nCorrections + CorrLevel]); + qvecReTPCnegEse.push_back(nModeQvecReEse[kTPCneg * nCorrections + CorrLevel]); + qvecImTPCnegEse.push_back(nModeQvecImEse[kTPCneg * nCorrections + CorrLevel]); + qvecReTPCallEse.push_back(nModeQvecReEse[kTPCall * nCorrections + CorrLevel]); + qvecImTPCallEse.push_back(nModeQvecImEse[kTPCall * nCorrections + CorrLevel]); } // Fill the columns of the Qvectors table. - qVector(cent, IsCalibrated, qvecRe, qvecIm, qvecAmp); + qVector(cent, IsCalibrated, qvecReSp, qvecImSp, qvecAmp); if (useDetector["QvectorFT0Cs"]) - qVectorFT0C(IsCalibrated, qvecReFT0C.at(0), qvecImFT0C.at(0), qvecAmp[kFT0C]); + qVectorFT0C(IsCalibrated, qvecReFT0CSp.at(0), qvecImFT0CSp.at(0), qvecAmp[kFT0C]); if (useDetector["QvectorFT0As"]) - qVectorFT0A(IsCalibrated, qvecReFT0A.at(0), qvecImFT0A.at(0), qvecAmp[kFT0A]); + qVectorFT0A(IsCalibrated, qvecReFT0ASp.at(0), qvecImFT0ASp.at(0), qvecAmp[kFT0A]); if (useDetector["QvectorFT0Ms"]) - qVectorFT0M(IsCalibrated, qvecReFT0M.at(0), qvecImFT0M.at(0), qvecAmp[kFT0M]); + qVectorFT0M(IsCalibrated, qvecReFT0MSp.at(0), qvecImFT0MSp.at(0), qvecAmp[kFT0M]); if (useDetector["QvectorFV0As"]) - qVectorFV0A(IsCalibrated, qvecReFV0A.at(0), qvecImFV0A.at(0), qvecAmp[kFV0A]); + qVectorFV0A(IsCalibrated, qvecReFV0ASp.at(0), qvecImFV0ASp.at(0), qvecAmp[kFV0A]); if (useDetector["QvectorTPCposs"]) - qVectorTPCpos(IsCalibrated, qvecReTPCpos.at(0), qvecImTPCpos.at(0), qvecAmp[kTPCpos], TrkTPCposLabel); + qVectorTPCpos(IsCalibrated, qvecReTPCposSp.at(0), qvecImTPCposSp.at(0), qvecAmp[kTPCpos], TrkTPCposLabel); if (useDetector["QvectorTPCnegs"]) - qVectorTPCneg(IsCalibrated, qvecReTPCneg.at(0), qvecImTPCneg.at(0), qvecAmp[kTPCneg], TrkTPCnegLabel); + qVectorTPCneg(IsCalibrated, qvecReTPCnegSp.at(0), qvecImTPCnegSp.at(0), qvecAmp[kTPCneg], TrkTPCnegLabel); if (useDetector["QvectorTPCalls"]) - qVectorTPCall(IsCalibrated, qvecReTPCall.at(0), qvecImTPCall.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + qVectorTPCall(IsCalibrated, qvecReTPCallSp.at(0), qvecImTPCallSp.at(0), qvecAmp[kTPCall], TrkTPCallLabel); - double qVecRedFT0C{-999.}, qVecRedTpcPos{-999.}, qVecRedTpcNeg{-999.}, qVecRedTpcAll{-999.}; - if (cfgProduceRedQVecs) { - // Correct normalization to remove multiplicity dependence, - // taking into account that the Q-vector is normalized by 1/M - // and the EsE reduced Q-vector must be normalized to 1/sqrt(M) - if (useDetector["QvectorFT0Cs"]) { - qVecRedFT0C = std::hypot(qvecReFT0C.at(0), qvecImFT0C.at(0)) * std::sqrt(qvecAmp[kFT0C]); - } - if (useDetector["QvectorTPCposs"]) { - qVecRedTpcPos = std::hypot(qvecReTPCpos.at(0), qvecImTPCpos.at(0)) * std::sqrt(qvecAmp[kTPCpos]); - } - if (useDetector["QvectorTPCnegs"]) { - qVecRedTpcNeg = std::hypot(qvecReTPCneg.at(0), qvecImTPCneg.at(0)) * std::sqrt(qvecAmp[kTPCneg]); - } - if (useDetector["QvectorTPCalls"]) { - qVecRedTpcAll = std::hypot(qvecReTPCall.at(0), qvecImTPCall.at(0)) * std::sqrt(qvecAmp[kTPCall]); - } - } - qVectorRed(qVecRedFT0C, qVecRedTpcPos, qVecRedTpcNeg, qVecRedTpcAll); - - qVectorFT0CVec(IsCalibrated, qvecReFT0C, qvecImFT0C, qvecAmp[kFT0C]); - qVectorFT0AVec(IsCalibrated, qvecReFT0A, qvecImFT0A, qvecAmp[kFT0A]); - qVectorFT0MVec(IsCalibrated, qvecReFT0M, qvecImFT0M, qvecAmp[kFT0M]); - qVectorFV0AVec(IsCalibrated, qvecReFV0A, qvecImFV0A, qvecAmp[kFV0A]); - qVectorTPCposVec(IsCalibrated, qvecReTPCpos, qvecImTPCpos, qvecAmp[kTPCpos], TrkTPCposLabel); - qVectorTPCnegVec(IsCalibrated, qvecReTPCneg, qvecImTPCneg, qvecAmp[kTPCneg], TrkTPCnegLabel); - qVectorTPCallVec(IsCalibrated, qvecReTPCall, qvecImTPCall, qvecAmp[kTPCall], TrkTPCallLabel); + qVectorFT0CVec(IsCalibrated, qvecReFT0CSp, qvecImFT0CSp, qvecAmp[kFT0C]); + qVectorFT0AVec(IsCalibrated, qvecReFT0ASp, qvecImFT0ASp, qvecAmp[kFT0A]); + qVectorFT0MVec(IsCalibrated, qvecReFT0MSp, qvecImFT0MSp, qvecAmp[kFT0M]); + qVectorFV0AVec(IsCalibrated, qvecReFV0ASp, qvecImFV0ASp, qvecAmp[kFV0A]); + qVectorTPCposVec(IsCalibrated, qvecReTPCposSp, qvecImTPCposSp, qvecAmp[kTPCpos], TrkTPCposLabel); + qVectorTPCnegVec(IsCalibrated, qvecReTPCnegSp, qvecImTPCnegSp, qvecAmp[kTPCneg], TrkTPCnegLabel); + qVectorTPCallVec(IsCalibrated, qvecReTPCallSp, qvecImTPCallSp, qvecAmp[kTPCall], TrkTPCallLabel); // Deprecated, will be removed in future after transition time // if (useDetector["QvectorBPoss"]) - qVectorBPos(IsCalibrated, qvecReTPCpos.at(0), qvecImTPCpos.at(0), qvecAmp[kTPCpos], TrkTPCposLabel); + qVectorBPos(IsCalibrated, qvecReTPCposSp.at(0), qvecImTPCposSp.at(0), qvecAmp[kTPCpos], TrkTPCposLabel); if (useDetector["QvectorBNegs"]) - qVectorBNeg(IsCalibrated, qvecReTPCneg.at(0), qvecImTPCneg.at(0), qvecAmp[kTPCneg], TrkTPCnegLabel); + qVectorBNeg(IsCalibrated, qvecReTPCnegSp.at(0), qvecImTPCnegSp.at(0), qvecAmp[kTPCneg], TrkTPCnegLabel); if (useDetector["QvectorBTots"]) - qVectorBTot(IsCalibrated, qvecReTPCall.at(0), qvecImTPCall.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + qVectorBTot(IsCalibrated, qvecReTPCallSp.at(0), qvecImTPCallSp.at(0), qvecAmp[kTPCall], TrkTPCallLabel); - qVectorBPosVec(IsCalibrated, qvecReTPCpos, qvecImTPCpos, qvecAmp[kTPCpos], TrkTPCposLabel); - qVectorBNegVec(IsCalibrated, qvecReTPCneg, qvecImTPCneg, qvecAmp[kTPCneg], TrkTPCnegLabel); - qVectorBTotVec(IsCalibrated, qvecReTPCall, qvecImTPCall, qvecAmp[kTPCall], TrkTPCallLabel); + qVectorBPosVec(IsCalibrated, qvecReTPCposSp, qvecImTPCposSp, qvecAmp[kTPCpos], TrkTPCposLabel); + qVectorBNegVec(IsCalibrated, qvecReTPCnegSp, qvecImTPCnegSp, qvecAmp[kTPCneg], TrkTPCnegLabel); + qVectorBTotVec(IsCalibrated, qvecReTPCallSp, qvecImTPCallSp, qvecAmp[kTPCall], TrkTPCallLabel); ///////////////////////////////////////////////////////////////// + if (cfgProduceRedQVecs) { + eseQVector(cent, IsCalibrated, qvecReEse, qvecImEse, qvecAmp); + eseQVectorFT0C(IsCalibrated, qvecReFT0CEse.at(0), qvecImFT0CEse.at(0), qvecAmp[kFT0C]); + eseQVectorFT0A(IsCalibrated, qvecReFT0AEse.at(0), qvecImFT0AEse.at(0), qvecAmp[kFT0A]); + eseQVectorFT0M(IsCalibrated, qvecReFT0MEse.at(0), qvecImFT0MEse.at(0), qvecAmp[kFT0M]); + eseQVectorFV0A(IsCalibrated, qvecReFV0AEse.at(0), qvecImFV0AEse.at(0), qvecAmp[kFV0A]); + eseQVectorTPCpos(IsCalibrated, qvecReTPCposEse.at(0), qvecImTPCposEse.at(0), qvecAmp[kTPCpos], TrkTPCposLabel); + eseQVectorTPCneg(IsCalibrated, qvecReTPCnegEse.at(0), qvecImTPCnegEse.at(0), qvecAmp[kTPCneg], TrkTPCnegLabel); + eseQVectorTPCall(IsCalibrated, qvecReTPCallEse.at(0), qvecImTPCallEse.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + eseQVectorFT0CVec(IsCalibrated, qvecReFT0CEse, qvecImFT0CEse, qvecAmp[kFT0C]); + eseQVectorFT0AVec(IsCalibrated, qvecReFT0AEse, qvecImFT0AEse, qvecAmp[kFT0A]); + eseQVectorFT0MVec(IsCalibrated, qvecReFT0MEse, qvecImFT0MEse, qvecAmp[kFT0M]); + eseQVectorFV0AVec(IsCalibrated, qvecReFV0AEse, qvecImFV0AEse, qvecAmp[kFV0A]); + eseQVectorTPCposVec(IsCalibrated, qvecReTPCposEse, qvecImTPCposEse, qvecAmp[kTPCpos], TrkTPCposLabel); + eseQVectorTPCnegVec(IsCalibrated, qvecReTPCnegEse, qvecImTPCnegEse, qvecAmp[kTPCneg], TrkTPCnegLabel); + eseQVectorTPCallVec(IsCalibrated, qvecReTPCallEse, qvecImTPCallEse, qvecAmp[kTPCall], TrkTPCallLabel); + eseQVectorPerc(qvecReFT0CEse.at(0), qvecImFT0CEse.at(0), qvecAmp[kFT0C], + qvecReFT0AEse.at(0), qvecImFT0AEse.at(0), qvecAmp[kFT0A], + qvecReFT0MEse.at(0), qvecImFT0MEse.at(0), qvecAmp[kFT0M], + qvecReFV0AEse.at(0), qvecImFV0AEse.at(0), qvecAmp[kFV0A], + qvecReTPCposEse.at(0), qvecImTPCposEse.at(0), qvecAmp[kTPCpos], + qvecReTPCnegEse.at(0), qvecImTPCnegEse.at(0), qvecAmp[kTPCneg], + qvecReTPCallEse.at(0), qvecImTPCallEse.at(0), qvecAmp[kTPCall]); + } } // End process. }; diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index a218cedb276..797f0504a58 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -153,7 +153,7 @@ struct HfTaskFlowCharmHadrons { using CandXic0DataWMl = soa::Filtered>; using CandD0DataWMl = soa::Filtered>; using CandD0Data = soa::Filtered>; - using CollsWithQvecs = soa::Join; + using CollsWithQvecs = soa::Join; using TracksWithExtra = soa::Join; Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; @@ -572,7 +572,7 @@ struct HfTaskFlowCharmHadrons { float redQVec{-999.f}; std::array qVecRedComps{-999.f, -999.f, -999.f}; if (storeRedQVec) { - qVecRedComps = getQvec(collision, qVecRedDetector.value); + qVecRedComps = getEseQvec(collision, qVecRedDetector.value); } float xRedQVec = qVecRedComps[0]; float yRedQVec = qVecRedComps[1]; @@ -769,9 +769,9 @@ struct HfTaskFlowCharmHadrons { // IMPORTANT: use the amplitude of the reduced Q-vector to build this Q-vector if constexpr (std::is_same_v || std::is_same_v) { - getQvecXic0Tracks(candidate, tracksRedQx, tracksRedQy, amplRedQVec, static_cast(qVecRedDetector.value)); + getQvecXic0Tracks(candidate, tracksRedQx, tracksRedQy, std::sqrt(amplRedQVec), static_cast(qVecRedDetector.value)); } else { - getQvecDtracks(candidate, tracksRedQx, tracksRedQy, amplRedQVec, static_cast(qVecRedDetector.value)); + getQvecDtracks(candidate, tracksRedQx, tracksRedQy, std::sqrt(amplRedQVec), static_cast(qVecRedDetector.value)); } // subtract daughters' contribution from the (normalized) Q-vector @@ -779,9 +779,9 @@ struct HfTaskFlowCharmHadrons { const float redQVecYDaugSubtr = yRedQVec - std::accumulate(tracksRedQy.begin(), tracksRedQy.end(), 0.0); if (qVecRedDetector.value == QvecEstimator::TPCTot || qVecRedDetector.value == QvecEstimator::TPCPos || qVecRedDetector.value == QvecEstimator::TPCNeg) { // Correct for track multiplicity - redQVec = std::hypot(redQVecXDaugSubtr, redQVecYDaugSubtr) * amplRedQVec / std::sqrt(amplRedQVec - tracksRedQx.size()); + redQVec = std::hypot(redQVecXDaugSubtr, redQVecYDaugSubtr) * std::sqrt(amplRedQVec) / std::sqrt(amplRedQVec - tracksRedQx.size()); } else { - redQVec = std::hypot(redQVecXDaugSubtr, redQVecYDaugSubtr) * std::sqrt(amplRedQVec); + redQVec = std::hypot(xRedQVec, yRedQVec); } } if (storeRedQVec) { @@ -971,7 +971,11 @@ struct HfTaskFlowCharmHadrons { o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(2)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(3))); } if (storeRedQVec) { - registry.fill(HIST("redQVecs/hSparseRedQVecs"), centrality, collision.qVecRedFT0C(), collision.qVecRedTpcPos(), collision.qVecRedTpcNeg(), collision.qVecRedTpcAll()); + registry.fill(HIST("redQVecs/hSparseRedQVecs"), centrality, + std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm()), + std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm()), + std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm()), + std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm())); } if (!isCollSelected(collision, bcs, centrality)) { @@ -997,10 +1001,10 @@ struct HfTaskFlowCharmHadrons { registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecTPCAll + yQVecFV0a * yQVecTPCAll); registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecTPCPos * xQVecTPCNeg + yQVecTPCPos * yQVecTPCNeg); - registry.fill(HIST("redQVecs/hRedQVecFT0C"), centrality, collision.qVecRedFT0C()); - registry.fill(HIST("redQVecs/hRedQVecTpcPos"), centrality, collision.qVecRedTpcPos()); - registry.fill(HIST("redQVecs/hRedQVecTpcNeg"), centrality, collision.qVecRedTpcNeg()); - registry.fill(HIST("redQVecs/hRedQVecTpcAll"), centrality, collision.qVecRedTpcAll()); + registry.fill(HIST("redQVecs/hRedQVecFT0C"), centrality, std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm())); + registry.fill(HIST("redQVecs/hRedQVecTpcPos"), centrality, std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm())); + registry.fill(HIST("redQVecs/hRedQVecTpcNeg"), centrality, std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm())); + registry.fill(HIST("redQVecs/hRedQVecTpcAll"), centrality, std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm())); if (saveEpResoHisto) { float const epFT0a = epHelper.GetEventPlane(xQVecFT0a, yQVecFT0a, harmonic); diff --git a/PWGHF/D2H/Utils/utilsFlow.h b/PWGHF/D2H/Utils/utilsFlow.h index 7efaedcdda3..911fc60dd2c 100644 --- a/PWGHF/D2H/Utils/utilsFlow.h +++ b/PWGHF/D2H/Utils/utilsFlow.h @@ -201,6 +201,56 @@ std::array getQvec(TCollision const& collision, const int qvecEst) } return std::array{-999.f, -999.f, -999.f}; } + +/// Get the Ese Q vector choosing your favourite estimator +/// \param collision is the collision with the Q vector information +/// \param qvecEst is the chosen Q-vector estimator +template +std::array getEseQvec(TCollision const& collision, const int qvecEst) +{ + switch (qvecEst) { + case QvecEstimator::FV0A: + if constexpr (HasQvecFV0A) { + return std::array{collision.eseQvecFV0ARe(), collision.eseQvecFV0AIm(), collision.sumAmplFV0A()}; + } + break; + case QvecEstimator::FT0A: + if constexpr (HasQvecFT0A) { + return std::array{collision.eseQvecFT0ARe(), collision.eseQvecFT0AIm(), collision.sumAmplFT0A()}; + } + break; + case QvecEstimator::FT0C: + if constexpr (HasQvecFT0C) { + return std::array{collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm(), collision.sumAmplFT0C()}; + } + break; + case QvecEstimator::FT0M: + if constexpr (HasQvecFT0M) { + return std::array{collision.eseQvecFT0MRe(), collision.eseQvecFT0MIm(), collision.sumAmplFT0M()}; + } + break; + case QvecEstimator::TPCPos: + if constexpr (HasQvecTPCpos) { + return std::array{collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm(), static_cast(collision.nTrkTPCpos())}; + } + break; + case QvecEstimator::TPCNeg: + if constexpr (HasQvecTPCneg) { + return std::array{collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm(), static_cast(collision.nTrkTPCneg())}; + } + break; + case QvecEstimator::TPCTot: + if constexpr (HasQvecTPCtot) { + return std::array{collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm(), static_cast(collision.nTrkTPCall())}; + } + break; + default: + LOGP(fatal, "Q-vector estimator not valid. Please choose between FV0A, FT0M, FT0A, FT0C, TPCPos, TPCNeg, TPCTot"); + break; + } + return std::array{-999.f, -999.f, -999.f}; +} + } // namespace hf_flow_utils } // namespace o2::analysis From 49f17349966878bc8eb48c09ea93137b9bf30f1a Mon Sep 17 00:00:00 2001 From: Marcello Di Costanzo Date: Tue, 28 Apr 2026 08:22:52 +0200 Subject: [PATCH 10/12] Normalization enum + debug prints --- Common/TableProducer/qVectorsTable.cxx | 64 +++++++++++++++++++------- 1 file changed, 48 insertions(+), 16 deletions(-) diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index 06ee7c6321a..cd795a1ba63 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -83,6 +83,12 @@ struct qVectorsTable { kRescale, kNCorrections }; + enum MultNorms { + kNoNorm = 0, + kScalarProd, + kEsE, + kMultNormTypes + }; // Configurables. struct : ConfigurableGroup { @@ -392,18 +398,30 @@ struct qVectorsTable { } } - void NormalizeQvec(std::vector& QvecReNorm, std::vector& QvecImNorm, std::vector QvecReRaw, std::vector QvecImRaw, std::vector& QvecAmp, bool useSqrt = false) { - + void NormalizeQvec(std::vector& QvecReNorm, + std::vector& QvecImNorm, + std::vector QvecReRaw, + std::vector QvecImRaw, + std::vector& QvecAmp, + MultNorms normType) + { for (std::size_t i = 0; i < kNDetectors; i++) { float qVecDetReNorm{999.}, qVecDetImNorm{999.}; if (QvecAmp[i] > 1e-8) { - if (useSqrt) { - qVecDetReNorm = QvecReRaw[i] / std::sqrt(QvecAmp[i]); - qVecDetImNorm = QvecImRaw[i] / std::sqrt(QvecAmp[i]); - } else { - qVecDetReNorm = QvecReRaw[i] / QvecAmp[i]; - qVecDetImNorm = QvecImRaw[i] / QvecAmp[i]; + switch (normType) { + case MultNorms::kScalarProd: + qVecDetReNorm = QvecReRaw[i] / QvecAmp[i]; + qVecDetImNorm = QvecImRaw[i] / QvecAmp[i]; + break; + case MultNorms::kEsE: + qVecDetReNorm = QvecReRaw[i] / std::sqrt(QvecAmp[i]); + qVecDetImNorm = QvecImRaw[i] / std::sqrt(QvecAmp[i]); + break; + default: + LOGP(fatal, "Undefined normalization type for Q-vector amplitude. Check the configuration."); + break; } + std::cout << "[NORMALIZED] " << i << " Re: " << qVecDetReNorm << ", Im: " << qVecDetImNorm << ", amp: " << QvecAmp[i] << std::endl; } for (int iCorr=0; iCorr < Corrections::kNCorrections; iCorr++) { QvecReNorm.push_back(qVecDetReNorm); @@ -548,7 +566,6 @@ struct qVectorsTable { helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nMode, QvecFT0M, sumAmplFT0M, ft0geom, fv0geom); } if (sumAmplFT0A > 1e-8) { - QvecDet /= sumAmplFT0A; qVectFT0A[0] = QvecDet.Re(); qVectFT0A[1] = QvecDet.Im(); } @@ -568,12 +585,10 @@ struct qVectorsTable { } if (sumAmplFT0C > 1e-8) { - QvecDet /= sumAmplFT0C; qVectFT0C[0] = QvecDet.Re(); qVectFT0C[1] = QvecDet.Im(); } if (sumAmplFT0M > 1e-8 && useDetector["QvectorFT0Ms"]) { - QvecFT0M /= sumAmplFT0M; qVectFT0M[0] = QvecFT0M.Re(); qVectFT0M[1] = QvecFT0M.Im(); } @@ -594,7 +609,6 @@ struct qVectorsTable { } if (sumAmplFV0A > 1e-8) { - QvecDet /= sumAmplFV0A; qVectFV0A[0] = QvecDet.Re(); qVectFV0A[1] = QvecDet.Im(); } @@ -658,10 +672,20 @@ struct qVectorsTable { QvecAmp.push_back(static_cast(nTrkTPCpos)); QvecAmp.push_back(static_cast(nTrkTPCneg)); QvecAmp.push_back(static_cast(nTrkTPCall)); + + LOG(info) << "[RAW] qVectFT0A: " << qVectFT0A[0] << ", " << qVectFT0A[1] << ", ampl: " << sumAmplFT0A; + LOG(info) << "[RAW] qVectFT0C: " << qVectFT0C[0] << ", " << qVectFT0C[1] << ", ampl: " << sumAmplFT0C; + LOG(info) << "[RAW] qVectFT0M: " << qVectFT0M[0] << ", " << qVectFT0M[1] << ", ampl: " << sumAmplFT0M; + LOG(info) << "[RAW] qVectFV0A: " << qVectFV0A[0] << ", " << qVectFV0A[1] << ", ampl: " << sumAmplFV0A; + LOG(info) << "[RAW] qVectTPCpos: " << qVectTPCpos[0] << ", " << qVectTPCpos[1] << ", nTrk: " << nTrkTPCpos; + LOG(info) << "[RAW] qVectTPCneg: " << qVectTPCneg[0] << ", " << qVectTPCneg[1] << ", nTrk: " << nTrkTPCneg; + LOG(info) << "[RAW] qVectTPCall: " << qVectTPCall[0] << ", " << qVectTPCall[1] << ", nTrk: " << nTrkTPCall; + } void process(MyCollisions::iterator const& coll, aod::BCsWithTimestamps const&, aod::FT0s const&, aod::FV0As const&, MyTracks const& tracks) { + LOG(info) << "---------------------------- Processing Event ---------------------------"; std::vector TrkTPCposLabel{}; std::vector TrkTPCnegLabel{}; std::vector TrkTPCallLabel{}; @@ -729,9 +753,8 @@ struct qVectorsTable { // Scalar Product Q-vectors, normalization by multiplicity/amplitude std::vector nModeQvecReSp{}; std::vector nModeQvecImSp{}; - NormalizeQvec(nModeQvecReSp, nModeQvecImSp, qvecReRaw, qvecImRaw, qvecAmp, false); + NormalizeQvec(nModeQvecReSp, nModeQvecImSp, qvecReRaw, qvecImRaw, qvecAmp, MultNorms::kScalarProd); CorrectQvec(cent, nModeQvecReSp, nModeQvecImSp, corrsQvecSp[id], nMode); - // Add to summary vector qvecReSp.insert(qvecReSp.end(), nModeQvecReSp.begin(), nModeQvecReSp.end()); qvecImSp.insert(qvecImSp.end(), nModeQvecImSp.begin(), nModeQvecImSp.end()); @@ -739,9 +762,8 @@ struct qVectorsTable { // Ese Q-vectors, normalization by sqrt(multiplicity/amplitude) std::vector nModeQvecReEse{}; std::vector nModeQvecImEse{}; - NormalizeQvec(nModeQvecReEse, nModeQvecImEse, qvecReRaw, qvecImRaw, qvecAmp, true); + NormalizeQvec(nModeQvecReEse, nModeQvecImEse, qvecReRaw, qvecImRaw, qvecAmp, MultNorms::kEsE); CorrectQvec(cent, nModeQvecReEse, nModeQvecImEse, corrsQvecEse[id], nMode); - // Add to summary vector qvecReEse.insert(qvecReEse.end(), nModeQvecReEse.begin(), nModeQvecReEse.end()); qvecImEse.insert(qvecImEse.end(), nModeQvecImEse.begin(), nModeQvecImEse.end()); @@ -798,6 +820,16 @@ struct qVectorsTable { if (useDetector["QvectorTPCalls"]) qVectorTPCall(IsCalibrated, qvecReTPCallSp.at(0), qvecImTPCallSp.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + // Debug prints of values after corrections + std::cout << "[CORRECTED] FT0C, Re: " << qvecReFT0CSp.at(0) << ", Im: " << qvecImFT0CSp.at(0) << std::endl; + std::cout << "[CORRECTED] FT0A, Re: " << qvecReFT0ASp.at(0) << ", Im: " << qvecImFT0ASp.at(0) << std::endl; + std::cout << "[CORRECTED] FT0M, Re: " << qvecReFT0MSp.at(0) << ", Im: " << qvecImFT0MSp.at(0) << std::endl; + std::cout << "[CORRECTED] FV0A, Re: " << qvecReFV0ASp.at(0) << ", Im: " << qvecImFV0ASp.at(0) << std::endl; + std::cout << "[CORRECTED] TPCpos, Re: " << qvecReTPCposSp.at(0) << ", Im: " << qvecImTPCposSp.at(0) << std::endl; + std::cout << "[CORRECTED] TPCneg, Re: " << qvecReTPCnegSp.at(0) << ", Im: " << qvecImTPCnegSp.at(0) << std::endl; + std::cout << "[CORRECTED] TPCall, Re: " << qvecReTPCallSp.at(0) << ", Im: " << qvecImTPCallSp.at(0) << std::endl; + + qVectorFT0CVec(IsCalibrated, qvecReFT0CSp, qvecImFT0CSp, qvecAmp[kFT0C]); qVectorFT0AVec(IsCalibrated, qvecReFT0ASp, qvecImFT0ASp, qvecAmp[kFT0A]); qVectorFT0MVec(IsCalibrated, qvecReFT0MSp, qvecImFT0MSp, qvecAmp[kFT0M]); From 42a7fe0ed5a3b53590271525af4eeb7dd9c8f1df Mon Sep 17 00:00:00 2001 From: Marcello Di Costanzo Date: Tue, 28 Apr 2026 09:14:15 +0200 Subject: [PATCH 11/12] Add documentation --- Common/TableProducer/qVectorsTable.cxx | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index cd795a1ba63..2fcdc885eb2 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -398,6 +398,13 @@ struct qVectorsTable { } } + /// Function to normalize the q-vectors + /// \param QvecReNorm is the vector with the normalized real part of the q-vector for each detector and correction step + /// \param QvecImNorm is the vector with the normalized imaginary part of the q-vector for each detector and correction step + /// \param QvecReRaw is the vector with the raw real part of the q-vector for each detector and correction step + /// \param QvecImRaw is the vector with the raw imaginary part of the q-vector for each detector and correction step + /// \param QvecAmp is the vector with the amplitude of the q-vector for each detector and correction step + /// \param normType is the type of normalization to apply to the q-vectors void NormalizeQvec(std::vector& QvecReNorm, std::vector& QvecImNorm, std::vector QvecReRaw, @@ -430,6 +437,12 @@ struct qVectorsTable { } } + /// Function to calculate the un-normalized q-vectors + /// \param cent is the collision centrality + /// \param qvecRe is the vector with the real part of the q-vector for each detector and correction step + /// \param qvecIm is the vector with the imaginary part of the q-vector for each detector and correction step + /// \param histsCorrs is the vector with the histograms with the correction constants for each detector and correction step + /// \param nMode is the modulation of interest void CorrectQvec(float cent, std::vector& qvecRe, std::vector& qvecIm, TH3F* histsCorrs, int nMode) { int nCorrections = static_cast(Corrections::kNCorrections); if (cent < cfgMaxCentrality) { @@ -533,6 +546,16 @@ struct qVectorsTable { } } + /// Function to calculate the un-normalized q-vectors + /// \param nMode is the harmonic number of the q-vector + /// \param coll is the collision object + /// \param track are the tracks associated to the collision + /// \param QvecRe is the vector with the real part of the q-vector for each detector + /// \param QvecIm is the vector with the imaginary part of the q-vector for each detector + /// \param QvecAmp is the vector with the amplitude of the signal in each detector + /// \param TrkTPCposLabel is the vector with the number of TPC tracks with positive eta + /// \param TrkTPCnegLabel is the vector with the number of TPC tracks with negative eta + /// \param TrkTPCallLabel is the vector with the number of TPC tracks with any eta template void CalQvec(const Nmode nMode, const CollType& coll, const TrackType& track, std::vector& QvecRe, std::vector& QvecIm, std::vector& QvecAmp, std::vector& TrkTPCposLabel, std::vector& TrkTPCnegLabel, std::vector& TrkTPCallLabel) { From be545a48622ddecaae43b58cc503998d8dca8c84 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 28 Apr 2026 07:15:36 +0000 Subject: [PATCH 12/12] Please consider the following formatting changes --- Common/DataModel/Qvectors.h | 12 +++++------ Common/TableProducer/qVectorsTable.cxx | 27 ++++++++++++------------ PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 10 ++++----- 3 files changed, 24 insertions(+), 25 deletions(-) diff --git a/Common/DataModel/Qvectors.h b/Common/DataModel/Qvectors.h index 71b4a496ff6..f62e8851a2b 100644 --- a/Common/DataModel/Qvectors.h +++ b/Common/DataModel/Qvectors.h @@ -221,12 +221,12 @@ DECLARE_SOA_TABLE(EseQvecTPCnegVecs, "AOD", "ESEQVECTPCNVEC", ese_qvec::IsCalibr DECLARE_SOA_TABLE(EseQvecTPCallVecs, "AOD", "ESEQVECTPCAVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCallReVec, ese_qvec::EseQvecTPCallImVec, qvec::NTrkTPCall, qvec::LabelsTPCall); DECLARE_SOA_TABLE(EseQvecPercs, "AOD", "ESEQVECPERC", ese_qvec::EseQvecFT0CRe, ese_qvec::EseQvecFT0CIm, qvec::SumAmplFT0C, - ese_qvec::EseQvecFT0ARe, ese_qvec::EseQvecFT0AIm, qvec::SumAmplFT0A, - ese_qvec::EseQvecFT0MRe, ese_qvec::EseQvecFT0MIm, qvec::SumAmplFT0M, - ese_qvec::EseQvecFV0ARe, ese_qvec::EseQvecFV0AIm, qvec::SumAmplFV0A, - ese_qvec::EseQvecTPCposRe, ese_qvec::EseQvecTPCposIm, qvec::NTrkTPCpos, - ese_qvec::EseQvecTPCnegRe, ese_qvec::EseQvecTPCnegIm, qvec::NTrkTPCneg, - ese_qvec::EseQvecTPCallRe, ese_qvec::EseQvecTPCallIm, qvec::NTrkTPCall); + ese_qvec::EseQvecFT0ARe, ese_qvec::EseQvecFT0AIm, qvec::SumAmplFT0A, + ese_qvec::EseQvecFT0MRe, ese_qvec::EseQvecFT0MIm, qvec::SumAmplFT0M, + ese_qvec::EseQvecFV0ARe, ese_qvec::EseQvecFV0AIm, qvec::SumAmplFV0A, + ese_qvec::EseQvecTPCposRe, ese_qvec::EseQvecTPCposIm, qvec::NTrkTPCpos, + ese_qvec::EseQvecTPCnegRe, ese_qvec::EseQvecTPCnegIm, qvec::NTrkTPCneg, + ese_qvec::EseQvecTPCallRe, ese_qvec::EseQvecTPCallIm, qvec::NTrkTPCall); using EseQvectorFT0C = EseQvecFT0Cs::iterator; using EseQvectorFT0A = EseQvecFT0As::iterator; diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index 2fcdc885eb2..c801d4c4ab3 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -317,7 +317,7 @@ struct qVectorsTable { auto modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); if (!modeCorrQvecEse) { fullPath = cfgQvecCalibPath; // cfgQvecEseCalibPath; - fullPath += "/v2"; // "/eseq2"; + fullPath += "/v2"; // "/eseq2"; modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); } corrsQvecEse.push_back(modeCorrQvecEse); @@ -405,12 +405,12 @@ struct qVectorsTable { /// \param QvecImRaw is the vector with the raw imaginary part of the q-vector for each detector and correction step /// \param QvecAmp is the vector with the amplitude of the q-vector for each detector and correction step /// \param normType is the type of normalization to apply to the q-vectors - void NormalizeQvec(std::vector& QvecReNorm, - std::vector& QvecImNorm, - std::vector QvecReRaw, - std::vector QvecImRaw, - std::vector& QvecAmp, - MultNorms normType) + void NormalizeQvec(std::vector& QvecReNorm, + std::vector& QvecImNorm, + std::vector QvecReRaw, + std::vector QvecImRaw, + std::vector& QvecAmp, + MultNorms normType) { for (std::size_t i = 0; i < kNDetectors; i++) { float qVecDetReNorm{999.}, qVecDetImNorm{999.}; @@ -430,7 +430,7 @@ struct qVectorsTable { } std::cout << "[NORMALIZED] " << i << " Re: " << qVecDetReNorm << ", Im: " << qVecDetImNorm << ", amp: " << QvecAmp[i] << std::endl; } - for (int iCorr=0; iCorr < Corrections::kNCorrections; iCorr++) { + for (int iCorr = 0; iCorr < Corrections::kNCorrections; iCorr++) { QvecReNorm.push_back(qVecDetReNorm); QvecImNorm.push_back(qVecDetImNorm); } @@ -443,7 +443,8 @@ struct qVectorsTable { /// \param qvecIm is the vector with the imaginary part of the q-vector for each detector and correction step /// \param histsCorrs is the vector with the histograms with the correction constants for each detector and correction step /// \param nMode is the modulation of interest - void CorrectQvec(float cent, std::vector& qvecRe, std::vector& qvecIm, TH3F* histsCorrs, int nMode) { + void CorrectQvec(float cent, std::vector& qvecRe, std::vector& qvecIm, TH3F* histsCorrs, int nMode) + { int nCorrections = static_cast(Corrections::kNCorrections); if (cent < cfgMaxCentrality) { for (auto i{0u}; i < kTPCall + 1; i++) { @@ -563,9 +564,9 @@ struct qVectorsTable { float qVectFT0C[2] = {-999., -999.}; float qVectFT0M[2] = {-999., -999.}; float qVectFV0A[2] = {-999., -999.}; - float qVectTPCpos[2] = {0., 0.}; // Always computed - float qVectTPCneg[2] = {0., 0.}; // Always computed - float qVectTPCall[2] = {0., 0.}; // Always computed + float qVectTPCpos[2] = {0., 0.}; // Always computed + float qVectTPCneg[2] = {0., 0.}; // Always computed + float qVectTPCall[2] = {0., 0.}; // Always computed TComplex QvecDet(0); TComplex QvecFT0M(0); @@ -703,7 +704,6 @@ struct qVectorsTable { LOG(info) << "[RAW] qVectTPCpos: " << qVectTPCpos[0] << ", " << qVectTPCpos[1] << ", nTrk: " << nTrkTPCpos; LOG(info) << "[RAW] qVectTPCneg: " << qVectTPCneg[0] << ", " << qVectTPCneg[1] << ", nTrk: " << nTrkTPCneg; LOG(info) << "[RAW] qVectTPCall: " << qVectTPCall[0] << ", " << qVectTPCall[1] << ", nTrk: " << nTrkTPCall; - } void process(MyCollisions::iterator const& coll, aod::BCsWithTimestamps const&, aod::FT0s const&, aod::FV0As const&, MyTracks const& tracks) @@ -852,7 +852,6 @@ struct qVectorsTable { std::cout << "[CORRECTED] TPCneg, Re: " << qvecReTPCnegSp.at(0) << ", Im: " << qvecImTPCnegSp.at(0) << std::endl; std::cout << "[CORRECTED] TPCall, Re: " << qvecReTPCallSp.at(0) << ", Im: " << qvecImTPCallSp.at(0) << std::endl; - qVectorFT0CVec(IsCalibrated, qvecReFT0CSp, qvecImFT0CSp, qvecAmp[kFT0C]); qVectorFT0AVec(IsCalibrated, qvecReFT0ASp, qvecImFT0ASp, qvecAmp[kFT0A]); qVectorFT0MVec(IsCalibrated, qvecReFT0MSp, qvecImFT0MSp, qvecAmp[kFT0M]); diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 797f0504a58..14d9c2c190f 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -971,11 +971,11 @@ struct HfTaskFlowCharmHadrons { o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(2)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(3))); } if (storeRedQVec) { - registry.fill(HIST("redQVecs/hSparseRedQVecs"), centrality, - std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm()), - std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm()), - std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm()), - std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm())); + registry.fill(HIST("redQVecs/hSparseRedQVecs"), centrality, + std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm()), + std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm()), + std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm()), + std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm())); } if (!isCollSelected(collision, bcs, centrality)) {