diff --git a/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx b/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx index 3dabb700b77..3d39e10fbdc 100644 --- a/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx +++ b/PWGLF/TableProducer/Common/zdc2stagecalibration.cxx @@ -76,6 +76,9 @@ struct zdc2stagecalibration { Configurable useSpatialCalib{"useSpatialCalib", false, "use ZDC spatial calibration"}; Configurable confSpatialPath{"confSpatialPath", "", "CCDB path to ZDC spatial calibration"}; 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 confQRecenteringPath{"confQRecenteringPath", "", "CCDB path to the run-wise Q-vector recentering coefficients"}; Configurable cfgSpatialCalibBins{ "cfgSpatialCalibBins", 30, @@ -94,6 +97,12 @@ struct zdc2stagecalibration { uint64_t eorTimestamp = 0; TH2D* gainprofile = nullptr; TH3D* spatialprofile = nullptr; + TH2D* qRecenteringProfile = nullptr; + + static constexpr int kQRecenteringNFeatures = 21; + static constexpr int kQRecenteringNMatrixMoments = kQRecenteringNFeatures * (kQRecenteringNFeatures + 1) / 2; + static constexpr int kQRecenteringNComponents = 4; + static constexpr int kQRecenteringNMoments = kQRecenteringNMatrixMoments + kQRecenteringNComponents * kQRecenteringNFeatures; void init(o2::framework::InitContext&) { @@ -104,10 +113,16 @@ struct zdc2stagecalibration { AxisSpec towerAxis = {4, 0.0, 4.0, "tower"}; AxisSpec crossAxis = {120, -1.0, 1.0, "cross asymmetry"}; AxisSpec momentAxis = {15, 0.0, 15.0, "regression moment"}; + AxisSpec spatialMomentAxis = {6, 0.0, 6.0, "normalized pair regression moment"}; AxisSpec spatialCalibAxis = {cfgSpatialCalibBins.value, -1.0, 1.0, "cross asymmetry"}; AxisSpec phiAxis = {72, -3.14159265358979323846, 3.14159265358979323846, "#phi"}; AxisSpec crossCoordAxis = {2, 0.0, 2.0, "coordinate"}; AxisSpec ratioAxis = {240, 0.0, 2.4, "#Sigma tower/common"}; + AxisSpec qRecenteringMomentAxis = {kQRecenteringNMoments, 0.0, static_cast(kQRecenteringNMoments), "regression moment"}; + AxisSpec centralityAxis = {80, 0.0, 80.0, "centrality (%)"}; + AxisSpec vertexXYAxis = {100, -0.5, 0.5, "vertex x/y (cm)"}; + AxisSpec vertexZAxis = {100, -10.0, 10.0, "vertex z (cm)"}; + AxisSpec qComponentAxis = {4, 0.0, 4.0, "Q component"}; histos.add("hEvtSelInfo", "hEvtSelInfo", kTH1F, {{10, 0.0, 10.0}}); auto hEvtSelInfo = histos.get(HIST("hEvtSelInfo")); @@ -150,8 +165,31 @@ struct zdc2stagecalibration { // Stage-1 correction input. For each time bin the first 10 y bins contain sum(T_i*T_j), the next 4 contain sum(C*T_i), and the last contains the event count. histos.add("GainCalibration/hGainMomentsZNA", "ZNA linear-regression moments;time from SOR (h);moment index", kTH2D, {timeAxis, momentAxis}); histos.add("GainCalibration/hGainMomentsZNC", "ZNC linear-regression moments;time from SOR (h);moment index", kTH2D, {timeAxis, momentAxis}); - histos.add("SpatialCalibration/hSpatialMomentsZNA", "ZNA spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, momentAxis}); - histos.add("SpatialCalibration/hSpatialMomentsZNC", "ZNC spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, momentAxis}); + // Normalized two-pair spatial regression moments. For each (Xcross,Ycross) cell: + // u=(T1+T4)/C, v=(T2+T3)/C, and moments are u^2, u*v, v^2, u, v, N. + histos.add("SpatialCalibration/hSpatialMomentsZNA", "ZNA normalized pair spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, spatialMomentAxis}); + histos.add("SpatialCalibration/hSpatialMomentsZNC", "ZNC normalized pair spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, spatialMomentAxis}); + + // Q-recentering regression moments. Feature order: + // 0:1, 1:C, 2:T, 3:Vx, 4:Vy, 5:Vz, 6:C2, 7:T2, 8:Vx2, 9:Vy2, 10:Vz2, + // 11:C*T, 12:C*Vx, 13:C*Vy, 14:C*Vz, 15:T*Vx, 16:T*Vy, 17:T*Vz, + // 18:Vx*Vy, 19:Vx*Vz, 20:Vy*Vz. + // Moments 0..230 are the upper triangle of sum(X_i X_j). + // Then 4 blocks of 21 sum(X_i Q): QxA, QyA, QxC, QyC. + histos.add("QRecenteringCalibration/hRegressionMoments", "Q recentering regression moments;moment index;sum", kTH1D, {qRecenteringMomentAxis}); + + histos.add("QRecenteringQA/pQBeforeVsCentrality", "Q before recentering vs centrality;centrality (%);Q component;", kTProfile2D, {centralityAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQBeforeVsTime", "Q before recentering vs time;time from SOR (h);Q component;", kTProfile2D, {timeAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQBeforeVsVx", "Q before recentering vs v_{x};v_{x} (cm);Q component;", kTProfile2D, {vertexXYAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQBeforeVsVy", "Q before recentering vs v_{y};v_{y} (cm);Q component;", kTProfile2D, {vertexXYAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQBeforeVsVz", "Q before recentering vs v_{z};v_{z} (cm);Q component;", kTProfile2D, {vertexZAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQAfterVsCentrality", "Q after recentering vs centrality;centrality (%);Q component;", kTProfile2D, {centralityAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQAfterVsTime", "Q after recentering vs time;time from SOR (h);Q component;", kTProfile2D, {timeAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQAfterVsVx", "Q after recentering vs v_{x};v_{x} (cm);Q component;", kTProfile2D, {vertexXYAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQAfterVsVy", "Q after recentering vs v_{y};v_{y} (cm);Q component;", kTProfile2D, {vertexXYAxis, qComponentAxis}); + histos.add("QRecenteringQA/pQAfterVsVz", "Q after recentering vs v_{z};v_{z} (cm);Q component;", kTProfile2D, {vertexZAxis, qComponentAxis}); + histos.add("PhiQA/hPhiAfterRecenteringZNA", "ZNA #phi after Q recentering;#phi;events", kTH1F, {phiAxis}); + histos.add("PhiQA/hPhiAfterRecenteringZNC", "ZNC #phi after Q recentering;#phi;events", kTH1F, {phiAxis}); ccdb->setURL(cfgCcdbParam.cfgURL); ccdb->setCaching(true); @@ -172,6 +210,7 @@ struct zdc2stagecalibration { bcSOR = static_cast(runInfo.orbitSOR) * static_cast(o2::constants::lhc::LHCMaxBunches); gainprofile = nullptr; spatialprofile = nullptr; + qRecenteringProfile = nullptr; if (calibrationStage.value == 2 && useGainCallib.value && !confGainPath.value.empty()) { gainprofile = ccdb->getForTimeStamp(confGainPath.value, timestamp); if (!gainprofile) { @@ -185,6 +224,12 @@ struct zdc2stagecalibration { LOGF(warn, "No ZDC spatial calibration found for run %d at timestamp %llu", runNumber, static_cast(timestamp)); } } + if (calibrationStage.value == 2 && useQRecentering.value && !confQRecenteringPath.value.empty()) { + qRecenteringProfile = ccdb->getForTimeStamp(confQRecenteringPath.value, timestamp); + if (!qRecenteringProfile) { + LOGF(warn, "No ZDC Q-recentering calibration found for run %d at timestamp %llu", runNumber, static_cast(timestamp)); + } + } } float getTimeFromSOR(uint64_t globalBC) const @@ -213,6 +258,24 @@ struct zdc2stagecalibration { return gain; } + std::array makeQRecenteringFeatures(float centrality, float timeFromSOR, float vx, float vy, float vz) const + { + // Fixed scaling keeps the polynomial basis numerically well behaved and must be + // used identically when solving/applying the CCDB coefficients. + const double c = (static_cast(centrality) - 40.0) / 40.0; + const double runHours = (eorTimestamp > sorTimestamp) ? static_cast(eorTimestamp - sorTimestamp) * 1.e-3 / 3600.0 : static_cast(cfgMaxRunHours.value); + const double t = (runHours > 0.0) ? (2.0 * static_cast(timeFromSOR) / runHours - 1.0) : 0.0; + const double x = static_cast(vx) / 0.1; + const double y = static_cast(vy) / 0.1; + const double z = static_cast(vz) / 10.0; + + return {1.0, c, t, x, y, z, + c * c, t * t, x * x, y * y, z * z, + c * t, c * x, c * y, c * z, + t * x, t * y, t * z, + x * y, x * z, y * z}; + } + 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*/) @@ -373,34 +436,32 @@ struct zdc2stagecalibration { const double crossLookupYC = (zncCorr[2] - zncCorr[1]) / (zncCorr[2] + zncCorr[1]); if (calibrationStage.value == 2 && useGainCallib.value && deriveSpatialCalib.value && !useSpatialCalib.value) { - std::array spatialMomentsA{}; - std::array spatialMomentsC{}; - - int moment = 0; - - for (int i = 0; i < 4; ++i) { - for (int j = i; j < 4; ++j) { - spatialMomentsA[moment] = znaCorr[i] * znaCorr[j]; - spatialMomentsC[moment] = zncCorr[i] * zncCorr[j]; - ++moment; + const double uA = (znaCorr[0] + znaCorr[3]) / static_cast(znaEnergycommon); + const double vA = (znaCorr[1] + znaCorr[2]) / static_cast(znaEnergycommon); + const double uC = (zncCorr[0] + zncCorr[3]) / static_cast(zncEnergycommon); + const double vC = (zncCorr[1] + zncCorr[2]) / static_cast(zncEnergycommon); + + if (std::isfinite(uA) && std::isfinite(vA) && std::isfinite(uC) && std::isfinite(vC)) { + const std::array spatialMomentsA = { + uA * uA, + uA * vA, + vA * vA, + uA, + vA, + 1.0}; + const std::array spatialMomentsC = { + uC * uC, + uC * vC, + vC * vC, + uC, + vC, + 1.0}; + + for (int i = 0; i < 6; ++i) { + histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNA"), crossLookupXA, crossLookupYA, i + 0.5, spatialMomentsA[i]); + histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNC"), crossLookupXC, crossLookupYC, i + 0.5, spatialMomentsC[i]); } } - - for (int i = 0; i < 4; ++i) { - spatialMomentsA[10 + i] = - static_cast(znaEnergycommon) * znaCorr[i]; - - spatialMomentsC[10 + i] = - static_cast(zncEnergycommon) * zncCorr[i]; - } - - spatialMomentsA[14] = 1.0; - spatialMomentsC[14] = 1.0; - - for (int i = 0; i < 15; ++i) { - histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNA"), crossLookupXA, crossLookupYA, i + 0.5, spatialMomentsA[i]); - histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNC"), crossLookupXC, crossLookupYC, i + 0.5, spatialMomentsC[i]); - } } bool spatialOK = true; @@ -533,9 +594,79 @@ struct zdc2stagecalibration { fillTable(false, 0.f, 0.f, 0.f, 0.f); return; } - histos.fill(HIST("hEvtSelInfo"), 7.5); - fillTable(true, static_cast(qxZDCA), static_cast(qxZDCC), static_cast(qyZDCA), static_cast(qyZDCC)); + const auto qFeatures = makeQRecenteringFeatures(centrality, timeFromSOR, vx, vy, vz); + std::array qValues = {qxZDCA, qyZDCA, qxZDCC, qyZDCC}; + + for (int i = 0; i < 4; ++i) { + const double component = i + 0.5; + histos.fill(HIST("QRecenteringQA/pQBeforeVsCentrality"), centrality, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQBeforeVsTime"), timeFromSOR, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQBeforeVsVx"), vx, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQBeforeVsVy"), vy, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQBeforeVsVz"), vz, component, qValues[i]); + } + + if (calibrationStage.value == 2 && useGainCallib.value && useSpatialCalib.value && deriveQRecentering.value && !useQRecentering.value) { + int moment = 0; + for (int i = 0; i < kQRecenteringNFeatures; ++i) { + for (int j = i; j < kQRecenteringNFeatures; ++j) { + histos.fill(HIST("QRecenteringCalibration/hRegressionMoments"), moment + 0.5, qFeatures[i] * qFeatures[j]); + ++moment; + } + } + for (int component = 0; component < kQRecenteringNComponents; ++component) { + for (int i = 0; i < kQRecenteringNFeatures; ++i) { + const int index = kQRecenteringNMatrixMoments + component * kQRecenteringNFeatures + i; + histos.fill(HIST("QRecenteringCalibration/hRegressionMoments"), index + 0.5, qFeatures[i] * qValues[component]); + } + } + } + + bool qRecenteringOK = true; + if (calibrationStage.value == 2 && useQRecentering.value) { + if (!qRecenteringProfile) { + qRecenteringOK = false; + } else { + for (int component = 0; component < kQRecenteringNComponents; ++component) { + double predictedBias = 0.0; + for (int i = 0; i < kQRecenteringNFeatures; ++i) { + const double coefficient = qRecenteringProfile->GetBinContent(qRecenteringProfile->FindBin(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; + } + } + } + if (!qRecenteringOK) { + fillTable(false, 0.f, 0.f, 0.f, 0.f); + return; + } + + for (int i = 0; i < 4; ++i) { + const double component = i + 0.5; + histos.fill(HIST("QRecenteringQA/pQAfterVsCentrality"), centrality, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQAfterVsTime"), timeFromSOR, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQAfterVsVx"), vx, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQAfterVsVy"), vy, component, qValues[i]); + histos.fill(HIST("QRecenteringQA/pQAfterVsVz"), vz, component, qValues[i]); + } + + const double phiRecenteringA = std::atan2(qValues[1], qValues[0]); + const double phiRecenteringC = std::atan2(qValues[3], qValues[2]); + histos.fill(HIST("PhiQA/hPhiAfterRecenteringZNA"), phiRecenteringA); + histos.fill(HIST("PhiQA/hPhiAfterRecenteringZNC"), phiRecenteringC); + + histos.fill(HIST("hEvtSelInfo"), 7.5); + fillTable(true, static_cast(qValues[0]), static_cast(qValues[2]), static_cast(qValues[1]), static_cast(qValues[3])); histos.fill(HIST("hEvtSelInfo"), 8.5); } };