From 589f4ccb388ee2e985a58e64ab8e04c26dad2047 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Sun, 15 Nov 2026 11:35:57 +0100 Subject: [PATCH] Add correlation with TPC to estimate transverse response --- .../Common/zdc2stagecalibration.cxx | 274 ++++++++++++++---- 1 file changed, 215 insertions(+), 59 deletions(-) diff --git a/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx b/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx index 7b87fecc482..16b85a16ffe 100644 --- a/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx +++ b/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx @@ -20,6 +20,7 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/FT0Corrected.h" #include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" #include #include @@ -42,6 +43,7 @@ #include #include #include +#include #include #include @@ -78,7 +80,7 @@ struct zdc2stagecalibration { Configurable deriveSpatialCalib{"deriveSpatialCalib", false, "store spatial regression moments after gain correction"}; Configurable deriveQRecentering{"deriveQRecentering", false, "store Q-vector recentering regression moments after gain and spatial correction"}; Configurable useQRecentering{"useQRecentering", false, "apply Q-vector recentering calibration"}; - Configurable qRecenteringMethod{"qRecenteringMethod", 0, "Q recentering method: 0=current regression, 1=multidimensional cubic regression, 4=structured71 regression"}; + Configurable qRecenteringMethod{"qRecenteringMethod", 4, "Q recentering: 4=Structured71 derive/apply (recommended); 0=legacy21; 1=legacy35 apply / legacy combined derivation"}; Configurable confQRecenteringPath{"confQRecenteringPath", "", "CCDB path to the current run-wise Q-vector recentering coefficients"}; Configurable confQRecenteringMultidimPath{"confQRecenteringMultidimPath", "", "CCDB path to the multidimensional cubic Q-vector recentering coefficients"}; Configurable confQRecenteringStructured71Path{"confQRecenteringStructured71Path", "", "CCDB path to the 71-feature structured Q-vector recentering coefficients"}; @@ -86,6 +88,14 @@ struct zdc2stagecalibration { Configurable confQTwistPath{"confQTwistPath", "", "CCDB path to the run-wise Q-twist map"}; Configurable useQShift{"useQShift", false, "apply ZDC azimuthal shift correction after twist"}; Configurable confQShiftPath{"confQShiftPath", "", "CCDB path to the run-wise Q-shift Fourier map"}; + Configurable deriveQResponse{"deriveQResponse", false, "collect TPC-ZDC response moments after recentering"}; + Configurable useQResponse{"useQResponse", false, "apply approved run-wise diagonal Q response after recentering"}; + Configurable confQResponsePath{"confQResponsePath", "Users/s/skundu/ZDC/QResponse_1percent", "CCDB path for the 1-percent Q response map"}; + Configurable cfgQResponseCentMin{"cfgQResponseCentMin", 10, "lower centrality edge of approved response application range"}; + Configurable cfgQResponseCentMax{"cfgQResponseCentMax", 40, "upper centrality edge of approved response application range"}; + Configurable cfgTPCResponsePtMin{"cfgTPCResponsePtMin", 0.5f, "TPC response QA minimum pT (strict)"}; + Configurable cfgTPCResponsePtMax{"cfgTPCResponsePtMax", 3.0f, "TPC response QA maximum pT (strict)"}; + Configurable cfgTPCResponseEtaMax{"cfgTPCResponseEtaMax", 0.8f, "TPC response QA maximum absolute eta (strict)"}; Configurable cfgSpatialCalibBins{ "cfgSpatialCalibBins", 30, "number of bins per cross coordinate for spatial calibration"}; @@ -106,6 +116,7 @@ struct zdc2stagecalibration { TH3D* qRecenteringProfile = nullptr; TH2D* qTwistProfile = nullptr; TH2D* qShiftProfile = nullptr; + TH2D* qResponseMap = nullptr; static constexpr int kQRecenteringNFeatures = 21; static constexpr int kQRecenteringNMatrixMoments = kQRecenteringNFeatures * (kQRecenteringNFeatures + 1) / 2; @@ -142,6 +153,17 @@ struct zdc2stagecalibration { if (useQShift.value && (!useQTwist.value || confQShiftPath.value.empty())) { LOGF(fatal, "Shift correction requires twist correction and confQShiftPath"); } + if (deriveQResponse.value && (calibrationStage.value != 2 || !useQRecentering.value || qRecenteringMethod.value != 4 || + useQResponse.value || useQTwist.value || useQShift.value || + !(cfgTPCResponsePtMin.value >= 0.f && cfgTPCResponsePtMax.value > cfgTPCResponsePtMin.value && cfgTPCResponseEtaMax.value > 0.f))) { + LOGF(fatal, "Response derivation requires stage 2, Structured71 recentering, response/twist/shift disabled and valid TPC track selection"); + } + if (useQResponse.value && (calibrationStage.value != 2 || !useQRecentering.value || qRecenteringMethod.value != 4 || + useQTwist.value || useQShift.value || confQResponsePath.value.empty() || + cfgQResponseCentMin.value < 0 || cfgQResponseCentMax.value > 80 || + cfgQResponseCentMin.value >= cfgQResponseCentMax.value)) { + LOGF(fatal, "Q response correction requires stage 2, Structured71 recentering, twist/shift disabled and valid response centrality range/path"); + } const int nTimeBins = static_cast(std::ceil(cfgMaxRunHours.value * 60.f / cfgTimeSliceMinutes.value)); AxisSpec timeAxis = {nTimeBins, 0.0, cfgMaxRunHours.value, "time from SOR (h)"}; AxisSpec towerAxis = {4, 0.0, 4.0, "tower"}; @@ -159,6 +181,14 @@ struct zdc2stagecalibration { AxisSpec qRecenteringStructuredMomentAxis = {kQRecenteringStructuredNMoments, 0.0, static_cast(kQRecenteringStructuredNMoments), "structured regression moment"}; AxisSpec qRecenteringStructuredFlatnessMomentAxis = {kQRecenteringStructuredFlatnessNMoments, 0.0, static_cast(kQRecenteringStructuredFlatnessNMoments), "structured flatness sufficient-statistic index"}; AxisSpec centralityAxis = {80, 0.0, 80.0, "centrality (%)"}; + if (deriveQResponse.value) { + const AxisSpec responseCentAxis = {80, 0., 80., "centrality (%)"}; + const AxisSpec responseMomentAxis = {47, 0., 47., "response moment index"}; + histos.add("TPCResponseQA/pMomentsFull", "Raw TPC Q1/Q2 vs recentered ZDC (full eta);centrality (%);moment index;", kTProfile2D, {responseCentAxis, responseMomentAxis}); + histos.add("TPCResponseQA/pMomentsEtaPos", "Raw TPC Q1/Q2 vs recentered ZDC (positive eta);centrality (%);moment index;", kTProfile2D, {responseCentAxis, responseMomentAxis}); + histos.add("TPCResponseQA/pMomentsEtaNeg", "Raw TPC Q1/Q2 vs recentered ZDC (negative eta);centrality (%);moment index;", kTProfile2D, {responseCentAxis, responseMomentAxis}); + } + AxisSpec resolutionCentralityAxis = {50, 0.0, 50.0, "centrality (%)"}; AxisSpec resolutionVzAxis = {20, -10.0, 10.0, "v_{z} (cm)"}; AxisSpec vertexXYAxis = {100, -0.1, 0.1, "vertex x/y (cm)"}; @@ -236,7 +266,7 @@ struct zdc2stagecalibration { // an additional flatness penalty. They do not alter the 35-feature model or // the application-side CCDB format. For each (centrality,conditioning) bin: // indices 0..34 = sum(X_i), indices 35..38 = sum(QxA,QyA,QxC,QyC). - if (deriveQRecentering.value && qRecenteringMethod.value == 1) { + if (deriveQRecentering.value && (qRecenteringMethod.value == 1 || qRecenteringMethod.value == 4)) { histos.add("QRecenteringCalibration/hFlatnessMomentsTime", "Flatness sufficient statistics vs time;centrality (%);time from SOR (h);statistic index", kTH3D, {centralityAxis, timeAxis, qRecenteringFlatnessMomentAxis}); histos.add("QRecenteringCalibration/hFlatnessMomentsVx", "Flatness sufficient statistics vs v_{x};centrality (%);v_{x} (cm);statistic index", kTH3D, {centralityAxis, vertexXYAxis, qRecenteringFlatnessMomentAxis}); histos.add("QRecenteringCalibration/hFlatnessMomentsVy", "Flatness sufficient statistics vs v_{y};centrality (%);v_{y} (cm);statistic index", kTH3D, {centralityAxis, vertexXYAxis, qRecenteringFlatnessMomentAxis}); @@ -319,7 +349,7 @@ struct zdc2stagecalibration { histos.add("QCorrelationQA/pQACorrelationBeforeVsCentrality", "A-C Q correlations before recentering;centrality (%);correlation component;", kTProfile2D, {centralityAxis, qCorrelationComponentAxis}); histos.add("QCorrelationQA/pQACorrelationAfterVsCentrality", "A-C Q correlations after recentering;centrality (%);correlation component;", kTProfile2D, {centralityAxis, qCorrelationComponentAxis}); - // Final-calibration QA: after recentering, correlation-only twist and optional shift. + // Final-calibration QA: the last enabled stage, including response-only or recentering-only. histos.add("final_calibration/hQxAVsCentrality", "Final QxA vs centrality;centrality (%);Q_{x}^{A}", kTH2F, {centralityAxis, qValueAxis}); histos.add("final_calibration/hQyAVsCentrality", "Final QyA vs centrality;centrality (%);Q_{y}^{A}", kTH2F, {centralityAxis, qValueAxis}); histos.add("final_calibration/hQxCVsCentrality", "Final QxC vs centrality;centrality (%);Q_{x}^{C}", kTH2F, {centralityAxis, qValueAxis}); @@ -346,6 +376,11 @@ struct zdc2stagecalibration { histos.add("final_calibration/pQMagAVsCentrality", "Final |Q_{A}| vs centrality;centrality (%);<|Q_{A}|>", kTProfile, {centralityAxis}); histos.add("final_calibration/pQMagCVsCentrality", "Final |Q_{C}| vs centrality;centrality (%);<|Q_{C}|>", kTProfile, {centralityAxis}); + if (useQResponse.value) { + const AxisSpec scaleAxis = {4, 0., 4., "Ax,Ay,Cx,Cy"}; + histos.add("QResponseQA/pAppliedScalesVsCentrality", "Actual applied response factors;centrality (%);component;", kTProfile2D, {centralityAxis, scaleAxis}); + histos.add("QResponseQA/hAppliedVsCentrality", "Events with response correction;centrality (%)", kTH1F, {centralityAxis}); + } ccdb->setURL(cfgCcdbParam.cfgURL); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -368,6 +403,7 @@ struct zdc2stagecalibration { qRecenteringProfile = nullptr; qTwistProfile = nullptr; qShiftProfile = nullptr; + qResponseMap = nullptr; if (calibrationStage.value == 2 && useGainCallib.value && !confGainPath.value.empty()) { gainprofile = ccdb->getForTimeStamp(confGainPath.value, timestamp); if (!gainprofile) { @@ -397,6 +433,39 @@ struct zdc2stagecalibration { } } } + if (useQResponse.value) { + qResponseMap = ccdb->getForTimeStamp(confQResponsePath.value, timestamp); + if (!qResponseMap || qResponseMap->GetNbinsX() != 80 || qResponseMap->GetNbinsY() != 5 || + qResponseMap->GetXaxis()->GetXmin() != 0. || qResponseMap->GetXaxis()->GetXmax() != 80.) { + LOGF(fatal, "Missing or invalid 1-percent response map for run %d", runNumber); + } + const std::array labels{"ScaleAx", "ScaleAy", "ScaleCx", "ScaleCy", "Approved"}; + for (int j = 0; j < 5; ++j) { + if (std::string(qResponseMap->GetYaxis()->GetBinLabel(j + 1)) != labels[j]) { + LOGF(fatal, "Wrong Q response map component ordering for run %d", runNumber); + } + } + if (std::string(qResponseMap->GetName()) != std::string("hZDCQResponse_run") + std::to_string(runNumber)) { + LOGF(fatal, "Wrong run identity in Q response map for run %d", runNumber); + } + for (int i = 1; i <= 80; ++i) { + if (qResponseMap->GetXaxis()->GetBinLowEdge(i) != i - 1 || qResponseMap->GetXaxis()->GetBinUpEdge(i) != i) { + LOGF(fatal, "Nonuniform centrality axis in response map for run %d", runNumber); + } + if (i <= cfgQResponseCentMin.value || i > cfgQResponseCentMax.value) + continue; + if (qResponseMap->GetBinContent(i, 5) != 1.) { + LOGF(fatal, "Unapproved response map bin %d for run %d", i, runNumber); + } + for (int j = 1; j <= 4; ++j) { + const double scale = qResponseMap->GetBinContent(i, j); + if (!std::isfinite(scale) || scale <= 0.) { + LOGF(fatal, "Invalid response factor in bin %d for run %d", i, runNumber); + } + } + } + LOGF(info, "Applying run %d response map in [%d,%d) percent only; outside this range Q remains recentered", runNumber, cfgQResponseCentMin.value, cfgQResponseCentMax.value); + } if (calibrationStage.value == 2 && useQTwist.value) { qTwistProfile = ccdb->getForTimeStamp(confQTwistPath.value, timestamp); if (!qTwistProfile || qTwistProfile->GetNbinsX() != 80 || qTwistProfile->GetNbinsY() != 6) { @@ -529,7 +598,92 @@ struct zdc2stagecalibration { using MyCollisions = o2::soa::Join; - void process(MyCollisions::iterator const& collision, o2::aod::FT0s const& /*ft0s*/, o2::aod::FV0As const& /*fv0s*/, BCsRun3 const& /*bcs*/, o2::aod::Zdcs const& /*zdcs*/) + using ResponseTracks = o2::soa::Join; + + void fillTPCResponseMoments(ResponseTracks const& tracks, const std::array& qValues, + float centrality) + { + // QA only. No modification of qValues or table. + if (std::isfinite(qValues[0]) && std::isfinite(qValues[1]) && + std::isfinite(qValues[2]) && std::isfinite(qValues[3])) { + struct TPCSum { + std::array q{}; // X1,Y1,X2,Y2 + int n = 0; + }; + std::array sums{}; // full, positive eta, negative eta + for (const auto& track : tracks) { + if (!track.isGlobalTrack() || !std::isfinite(track.pt()) || !std::isfinite(track.eta()) || !std::isfinite(track.phi()) || + track.pt() <= cfgTPCResponsePtMin.value || track.pt() >= cfgTPCResponsePtMax.value || + std::abs(track.eta()) >= cfgTPCResponseEtaMax.value) { + continue; + } + const double c = std::cos(track.phi()); + const double t = std::sin(track.phi()); + const std::array u{c, t, c * c - t * t, 2. * c * t}; + auto add = [&](TPCSum& sum) { + ++sum.n; + for (int j = 0; j < 4; ++j) { + sum.q[j] += u[j]; + } + }; + add(sums[0]); + if (track.eta() > 0.f) { + add(sums[1]); + } else if (track.eta() < 0.f) { + add(sums[2]); + } + } + const double a = qValues[0], b = qValues[1], c = qValues[2], d = qValues[3]; + // Order: squares Ax,Ay,Cx,Cy; AC xx,yy,xy,yx; same-side xy A,C. + const std::array z2{a * a, b * b, c * c, d * d, a * c, b * d, a * d, b * c, a * b, c * d}; + for (int region = 0; region < 3; ++region) { + const auto& sum = sums[region]; + if (sum.n == 0) { + continue; // no TPC reference in this region; not an event/table cut + } + std::array tpc{}; + for (int j = 0; j < 4; ++j) { + tpc[j] = sum.q[j] / sum.n; + } + // All moments get the same event weight, so connected correlations can + // be constructed after merging. Indices below are ZERO-based. + std::array m{}; + m[0] = sum.n; + for (int j = 0; j < 4; ++j) { + m[1 + j] = tpc[j]; + m[5 + j] = qValues[j]; + } + for (int j = 0; j < 10; ++j) { + m[9 + j] = z2[j]; + } + // First harmonic diagnostic: X1 and Y1 times each ZDC component. + for (int h = 0; h < 2; ++h) { + for (int j = 0; j < 4; ++j) { + m[19 + 4 * h + j] = tpc[h] * qValues[j]; + } + } + // Second harmonic: X2 and Y2 times every ZDC quadratic component. + // Save all 20 products including null checks. Separate from Q1. + for (int h = 0; h < 2; ++h) { + for (int j = 0; j < 10; ++j) { + m[27 + 10 * h + j] = tpc[2 + h] * z2[j]; + } + } + for (int j = 0; j < static_cast(m.size()); ++j) { + const double idx = j + 0.5; + if (region == 0) { + histos.fill(HIST("TPCResponseQA/pMomentsFull"), centrality, idx, m[j]); + } else if (region == 1) { + histos.fill(HIST("TPCResponseQA/pMomentsEtaPos"), centrality, idx, m[j]); + } else { + histos.fill(HIST("TPCResponseQA/pMomentsEtaNeg"), centrality, idx, m[j]); + } + } + } + } + } + + void process(MyCollisions::iterator const& collision, ResponseTracks const& tracks, o2::aod::FT0s const&, o2::aod::FV0As const&, BCsRun3 const&, o2::aod::Zdcs const&) { histos.fill(HIST("hEvtSelInfo"), 0.5); @@ -1012,62 +1166,32 @@ struct zdc2stagecalibration { bool qRecenteringOK = true; if (calibrationStage.value == 2 && useQRecentering.value) { - if (!qRecenteringProfile) { - qRecenteringOK = false; - } else if (qRecenteringMethod.value == 0) { - const auto qFeatures = makeQRecenteringFeatures(centrality, timeFromSOR, vx, vy, vz); + // Same feature definitions, coefficient ordering and accumulation as the legacy implementation. + const auto subtractBias = [&](const auto& features) { for (int component = 0; component < kQRecenteringNComponents; ++component) { double predictedBias = 0.0; - for (int i = 0; i < kQRecenteringNFeatures; ++i) { + for (std::size_t i = 0; i < features.size(); ++i) { const double coefficient = qRecenteringProfile->GetBinContent(qRecenteringProfile->FindBin(centrality, i + 0.5, component + 0.5)); if (!std::isfinite(coefficient)) { - qRecenteringOK = false; - break; + return false; } - predictedBias += coefficient * qFeatures[i]; + predictedBias += coefficient * features[i]; } - if (!qRecenteringOK || !std::isfinite(predictedBias)) { - qRecenteringOK = false; - break; + if (!std::isfinite(predictedBias)) { + return false; } qValues[component] -= predictedBias; } + return true; + }; + if (!qRecenteringProfile) { + qRecenteringOK = false; } else if (qRecenteringMethod.value == 4) { - const auto qFeatures = makeQRecenteringStructuredFeatures(timeFromSOR, vx, vy, vz); - for (int component = 0; component < kQRecenteringNComponents; ++component) { - double predictedBias = 0.0; - for (int i = 0; i < kQRecenteringStructuredNFeatures; ++i) { - const double coefficient = qRecenteringProfile->GetBinContent(qRecenteringProfile->FindBin(centrality, i + 0.5, component + 0.5)); - if (!std::isfinite(coefficient)) { - qRecenteringOK = false; - break; - } - predictedBias += coefficient * qFeatures[i]; - } - if (!qRecenteringOK || !std::isfinite(predictedBias)) { - qRecenteringOK = false; - break; - } - qValues[component] -= predictedBias; - } + qRecenteringOK = subtractBias(makeQRecenteringStructuredFeatures(timeFromSOR, vx, vy, vz)); + } else if (qRecenteringMethod.value == 1) { + qRecenteringOK = subtractBias(makeQRecenteringMultidimFeatures(timeFromSOR, vx, vy, vz)); } else { - const auto qFeatures = makeQRecenteringMultidimFeatures(timeFromSOR, vx, vy, vz); - for (int component = 0; component < kQRecenteringNComponents; ++component) { - double predictedBias = 0.0; - for (int i = 0; i < kQRecenteringMultidimNFeatures; ++i) { - const double coefficient = qRecenteringProfile->GetBinContent(qRecenteringProfile->FindBin(centrality, i + 0.5, component + 0.5)); - if (!std::isfinite(coefficient)) { - qRecenteringOK = false; - break; - } - predictedBias += coefficient * qFeatures[i]; - } - if (!qRecenteringOK || !std::isfinite(predictedBias)) { - qRecenteringOK = false; - break; - } - qValues[component] -= predictedBias; - } + qRecenteringOK = subtractBias(makeQRecenteringFeatures(centrality, timeFromSOR, vx, vy, vz)); } } if (!qRecenteringOK) { @@ -1109,6 +1233,11 @@ struct zdc2stagecalibration { histos.fill(HIST("QCorrelationQA/pQACorrelationAfterVsCentrality"), centrality, 2.5, qValues[0] * qValues[3]); histos.fill(HIST("QCorrelationQA/pQACorrelationAfterVsCentrality"), centrality, 3.5, qValues[1] * qValues[2]); + // Response derivation always uses the recentered, unscaled Q vectors. + if (deriveQResponse.value) { + fillTPCResponseMoments(tracks, qValues, centrality); + } + const double phiRecenteringA = std::atan2(qValues[1], qValues[0]); const double phiRecenteringC = std::atan2(qValues[3], qValues[2]); histos.fill(HIST("PhiQA/hPhiAfterRecenteringZNA"), phiRecenteringA); @@ -1132,6 +1261,22 @@ struct zdc2stagecalibration { histos.fill(HIST("QTwistQA/pMomentsBeforeVsCentrality"), centrality, i + 0.5, momentsBeforeTwist[i]); } } + // Response QA above always observes recentered, unscaled Q. Correction is optional. + if (useQResponse.value && centrality >= cfgQResponseCentMin.value && centrality < cfgQResponseCentMax.value) { + const int centBin = qResponseMap->GetXaxis()->FindFixBin(centrality); + if (centBin < 1 || centBin > qResponseMap->GetNbinsX()) { + LOGF(fatal, "Invalid response centrality bin for run %d", runNumber); + } + for (int j = 0; j < 4; ++j) { + const double scale = qResponseMap->GetBinContent(centBin, j + 1); + qValues[j] *= scale; + if (!std::isfinite(qValues[j]) || !std::isfinite(static_cast(qValues[j]))) { + LOGF(fatal, "Nonfinite Q after response correction, run %d bin %d", runNumber, centBin); + } + histos.fill(HIST("QResponseQA/pAppliedScalesVsCentrality"), centrality, j + 0.5, scale); + } + histos.fill(HIST("QResponseQA/hAppliedVsCentrality"), centrality); + } if (useQTwist.value) { const int centBin = qTwistProfile->GetXaxis()->FindFixBin(centrality); if (centBin < 1 || centBin > qTwistProfile->GetNbinsX()) { @@ -1228,7 +1373,8 @@ struct zdc2stagecalibration { histos.fill(HIST("ResolutionQA/pCosPsiAPsiCVsCentralityVz"), centrality, vz, correlation); } - if (calibrationStage.value == 2 && useQRecentering.value && useQTwist.value && qMagA > 0. && qMagC > 0.) { + // Always describe the Q vectors actually written to the stage-2 table, regardless of enabled corrections. + if (calibrationStage.value == 2) { histos.fill(HIST("final_calibration/hQxAVsCentrality"), centrality, qValues[0]); histos.fill(HIST("final_calibration/hQyAVsCentrality"), centrality, qValues[1]); histos.fill(HIST("final_calibration/hQxCVsCentrality"), centrality, qValues[2]); @@ -1237,11 +1383,17 @@ struct zdc2stagecalibration { histos.fill(HIST("final_calibration/pQyAVsCentrality"), centrality, qValues[1]); histos.fill(HIST("final_calibration/pQxCVsCentrality"), centrality, qValues[2]); histos.fill(HIST("final_calibration/pQyCVsCentrality"), centrality, qValues[3]); - histos.fill(HIST("final_calibration/hPhiAVsCentrality"), centrality, phiFinalA); - histos.fill(HIST("final_calibration/hPhiCVsCentrality"), centrality, phiFinalC); - const double finalCorrelation = std::cos(phiFinalA - phiFinalC); - histos.fill(HIST("final_calibration/pCosPsiAPsiCVsCentrality"), centrality, finalCorrelation); - histos.fill(HIST("final_calibration/pCosPsiAPsiCVsCentralityVz"), centrality, vz, finalCorrelation); + if (qMagA > 0.) { + histos.fill(HIST("final_calibration/hPhiAVsCentrality"), centrality, phiFinalA); + } + if (qMagC > 0.) { + histos.fill(HIST("final_calibration/hPhiCVsCentrality"), centrality, phiFinalC); + } + if (qMagA > 0. && qMagC > 0.) { + const double finalCorrelation = std::cos(phiFinalA - phiFinalC); + histos.fill(HIST("final_calibration/pCosPsiAPsiCVsCentrality"), centrality, finalCorrelation); + histos.fill(HIST("final_calibration/pCosPsiAPsiCVsCentralityVz"), centrality, vz, finalCorrelation); + } histos.fill(HIST("final_calibration/pQxAQxCVsCentrality"), centrality, qValues[0] * qValues[2]); histos.fill(HIST("final_calibration/pQyAQyCVsCentrality"), centrality, qValues[1] * qValues[3]); histos.fill(HIST("final_calibration/pQxAQyCVsCentrality"), centrality, qValues[0] * qValues[3]); @@ -1255,10 +1407,14 @@ struct zdc2stagecalibration { histos.fill(HIST("final_calibration/pQMomentsVsCentrality"), centrality, i + 0.5, finalMoments[i]); } for (int n = 1; n <= 4; ++n) { - histos.fill(HIST("final_calibration/pPhiHarmonicsAVsCentrality"), centrality, 2 * n - 1.5, std::cos(n * phiFinalA)); - histos.fill(HIST("final_calibration/pPhiHarmonicsAVsCentrality"), centrality, 2 * n - 0.5, std::sin(n * phiFinalA)); - histos.fill(HIST("final_calibration/pPhiHarmonicsCVsCentrality"), centrality, 2 * n - 1.5, std::cos(n * phiFinalC)); - histos.fill(HIST("final_calibration/pPhiHarmonicsCVsCentrality"), centrality, 2 * n - 0.5, std::sin(n * phiFinalC)); + if (qMagA > 0.) { + histos.fill(HIST("final_calibration/pPhiHarmonicsAVsCentrality"), centrality, 2 * n - 1.5, std::cos(n * phiFinalA)); + histos.fill(HIST("final_calibration/pPhiHarmonicsAVsCentrality"), centrality, 2 * n - 0.5, std::sin(n * phiFinalA)); + } + if (qMagC > 0.) { + histos.fill(HIST("final_calibration/pPhiHarmonicsCVsCentrality"), centrality, 2 * n - 1.5, std::cos(n * phiFinalC)); + histos.fill(HIST("final_calibration/pPhiHarmonicsCVsCentrality"), centrality, 2 * n - 0.5, std::sin(n * phiFinalC)); + } } histos.fill(HIST("final_calibration/pQMagAVsCentrality"), centrality, qMagA); histos.fill(HIST("final_calibration/pQMagCVsCentrality"), centrality, qMagC);