Skip to content

[PWGLF] Add recentering part of ZDC Q vector calibration #17895

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
skundu692 merged 1 commit into from
Sep 13, 2026
+161 −30
Merged

[PWGLF] Add recentering part of ZDC Q vector calibration #17895

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
191 changes: 161 additions & 30 deletions PWGLF/TableProducer/Common/zdc2stagecalibration.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,9 @@
Configurable<bool> useSpatialCalib{"useSpatialCalib", false, "use ZDC spatial calibration"};
Configurable<std::string> confSpatialPath{"confSpatialPath", "", "CCDB path to ZDC spatial calibration"};
Configurable<bool> deriveSpatialCalib{"deriveSpatialCalib", false, "store spatial regression moments after gain correction"};
Configurable<bool> deriveQRecentering{"deriveQRecentering", false, "store Q-vector recentering regression moments after gain and spatial correction"};
Configurable<bool> useQRecentering{"useQRecentering", false, "apply Q-vector recentering calibration"};
Configurable<std::string> confQRecenteringPath{"confQRecenteringPath", "", "CCDB path to the run-wise Q-vector recentering coefficients"};

Configurable<int> cfgSpatialCalibBins{
"cfgSpatialCalibBins", 30,
Expand All @@ -94,6 +97,12 @@
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&)
{
Expand All @@ -104,10 +113,16 @@
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<double>(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<TH1>(HIST("hEvtSelInfo"));
Expand Down Expand Up @@ -150,8 +165,31 @@
// 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;<Q>", kTProfile2D, {centralityAxis, qComponentAxis});
histos.add("QRecenteringQA/pQBeforeVsTime", "Q before recentering vs time;time from SOR (h);Q component;<Q>", kTProfile2D, {timeAxis, qComponentAxis});
histos.add("QRecenteringQA/pQBeforeVsVx", "Q before recentering vs v_{x};v_{x} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
histos.add("QRecenteringQA/pQBeforeVsVy", "Q before recentering vs v_{y};v_{y} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
histos.add("QRecenteringQA/pQBeforeVsVz", "Q before recentering vs v_{z};v_{z} (cm);Q component;<Q>", kTProfile2D, {vertexZAxis, qComponentAxis});
histos.add("QRecenteringQA/pQAfterVsCentrality", "Q after recentering vs centrality;centrality (%);Q component;<Q>", kTProfile2D, {centralityAxis, qComponentAxis});
histos.add("QRecenteringQA/pQAfterVsTime", "Q after recentering vs time;time from SOR (h);Q component;<Q>", kTProfile2D, {timeAxis, qComponentAxis});
histos.add("QRecenteringQA/pQAfterVsVx", "Q after recentering vs v_{x};v_{x} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
histos.add("QRecenteringQA/pQAfterVsVy", "Q after recentering vs v_{y};v_{y} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
histos.add("QRecenteringQA/pQAfterVsVz", "Q after recentering vs v_{z};v_{z} (cm);Q component;<Q>", 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);
Expand All @@ -172,19 +210,26 @@
bcSOR = static_cast<int64_t>(runInfo.orbitSOR) * static_cast<int64_t>(o2::constants::lhc::LHCMaxBunches);
gainprofile = nullptr;
spatialprofile = nullptr;
qRecenteringProfile = nullptr;
if (calibrationStage.value == 2 && useGainCallib.value && !confGainPath.value.empty()) {

Check failure on line 214 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
gainprofile = ccdb->getForTimeStamp<TH2D>(confGainPath.value, timestamp);
if (!gainprofile) {
LOGF(warn, "No ZDC stage-1 gain map found for run %d at timestamp %llu", runNumber, static_cast<unsigned long long>(timestamp));
}
}
if (calibrationStage.value == 2 && useSpatialCalib.value && !confSpatialPath.value.empty()) {

Check failure on line 220 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
spatialprofile = ccdb->getForTimeStamp<TH3D>(confSpatialPath.value, timestamp);

if (!spatialprofile) {
LOGF(warn, "No ZDC spatial calibration found for run %d at timestamp %llu", runNumber, static_cast<unsigned long long>(timestamp));
}
}
if (calibrationStage.value == 2 && useQRecentering.value && !confQRecenteringPath.value.empty()) {

Check failure on line 227 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
qRecenteringProfile = ccdb->getForTimeStamp<TH2D>(confQRecenteringPath.value, timestamp);
if (!qRecenteringProfile) {
LOGF(warn, "No ZDC Q-recentering calibration found for run %d at timestamp %llu", runNumber, static_cast<unsigned long long>(timestamp));
}
}
}

float getTimeFromSOR(uint64_t globalBC) const
Expand All @@ -198,7 +243,7 @@

double getGain(float timeFromSOR, int channel, bool& gainOK) const
{
if (calibrationStage.value != 2 || !useGainCallib.value) {

Check failure on line 246 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
return 1.0;
}
if (!gainprofile || timeFromSOR < 0.f) {
Expand All @@ -213,6 +258,24 @@
return gain;
}

std::array<double, kQRecenteringNFeatures> 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<double>(centrality) - 40.0) / 40.0;
const double runHours = (eorTimestamp > sorTimestamp) ? static_cast<double>(eorTimestamp - sorTimestamp) * 1.e-3 / 3600.0 : static_cast<double>(cfgMaxRunHours.value);
const double t = (runHours > 0.0) ? (2.0 * static_cast<double>(timeFromSOR) / runHours - 1.0) : 0.0;
const double x = static_cast<double>(vx) / 0.1;
const double y = static_cast<double>(vy) / 0.1;
const double z = static_cast<double>(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<o2::aod::Collisions, o2::aod::EvSels, o2::aod::Mults, o2::aod::FT0sCorrected, o2::aod::CentFT0Cs>;

void process(MyCollisions::iterator const& collision, o2::aod::FT0s const& /*ft0s*/, o2::aod::FV0As const& /*fv0s*/, BCsRun3 const& /*bcs*/, o2::aod::Zdcs const& /*zdcs*/)
Expand Down Expand Up @@ -254,7 +317,7 @@
histos.fill(HIST("hEvtSelInfo"), 2.5);

bool validTowers = true;
for (int i = 0; i < 4; ++i) {

Check failure on line 320 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
if (!std::isfinite(znaEnergy[i]) || !std::isfinite(zncEnergy[i]) || znaEnergy[i] <= 0.f || zncEnergy[i] <= 0.f) {
validTowers = false;
}
Expand Down Expand Up @@ -291,7 +354,7 @@
double rawYC = 0.0;
double rawSumA = 0.0;
double rawSumC = 0.0;
for (int i = 0; i < 4; ++i) {

Check failure on line 357 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
rawXA -= znaEnergy[i] * x[i];
rawYA += znaEnergy[i] * y[i];
rawXC += zncEnergy[i] * x[i];
Expand All @@ -312,20 +375,20 @@
std::array<double, 15> momentsA{};
std::array<double, 15> momentsC{};
int moment = 0;
for (int i = 0; i < 4; ++i) {

Check failure on line 378 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
for (int j = i; j < 4; ++j) {

Check failure on line 379 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
momentsA[moment] = static_cast<double>(znaEnergy[i]) * static_cast<double>(znaEnergy[j]);
momentsC[moment] = static_cast<double>(zncEnergy[i]) * static_cast<double>(zncEnergy[j]);
++moment;
}
}
for (int i = 0; i < 4; ++i) {

Check failure on line 385 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
momentsA[10 + i] = static_cast<double>(znaEnergycommon) * static_cast<double>(znaEnergy[i]);
momentsC[10 + i] = static_cast<double>(zncEnergycommon) * static_cast<double>(zncEnergy[i]);
}
momentsA[14] = 1.0;
momentsC[14] = 1.0;
for (int i = 0; i < 15; ++i) {

Check failure on line 391 in PWGLF/TableProducer/Common/zdc2stagecalibration.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
histos.fill(HIST("GainCalibration/hGainMomentsZNA"), timeFromSOR, i + 0.5, momentsA[i]);
histos.fill(HIST("GainCalibration/hGainMomentsZNC"), timeFromSOR, i + 0.5, momentsC[i]);
}
Expand Down Expand Up @@ -373,34 +436,32 @@
const double crossLookupYC = (zncCorr[2] - zncCorr[1]) / (zncCorr[2] + zncCorr[1]);
if (calibrationStage.value == 2 && useGainCallib.value && deriveSpatialCalib.value && !useSpatialCalib.value) {

std::array<double, 15> spatialMomentsA{};
std::array<double, 15> 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<double>(znaEnergycommon);
const double vA = (znaCorr[1] + znaCorr[2]) / static_cast<double>(znaEnergycommon);
const double uC = (zncCorr[0] + zncCorr[3]) / static_cast<double>(zncEnergycommon);
const double vC = (zncCorr[1] + zncCorr[2]) / static_cast<double>(zncEnergycommon);

if (std::isfinite(uA) && std::isfinite(vA) && std::isfinite(uC) && std::isfinite(vC)) {
const std::array<double, 6> spatialMomentsA = {
uA * uA,
uA * vA,
vA * vA,
uA,
vA,
1.0};
const std::array<double, 6> 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<double>(znaEnergycommon) * znaCorr[i];

spatialMomentsC[10 + i] =
static_cast<double>(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;
Expand Down Expand Up @@ -533,9 +594,79 @@
fillTable(false, 0.f, 0.f, 0.f, 0.f);
return;
}
histos.fill(HIST("hEvtSelInfo"), 7.5);

fillTable(true, static_cast<float>(qxZDCA), static_cast<float>(qxZDCC), static_cast<float>(qyZDCA), static_cast<float>(qyZDCC));
const auto qFeatures = makeQRecenteringFeatures(centrality, timeFromSOR, vx, vy, vz);
std::array<double, 4> 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<float>(qValues[0]), static_cast<float>(qValues[2]), static_cast<float>(qValues[1]), static_cast<float>(qValues[3]));
histos.fill(HIST("hEvtSelInfo"), 8.5);
}
};
Expand Down
Loading