Skip to content

Commit 730d050

Browse files
skundu692skundu692
committed
Add recentering part of ZDC Q vector calibration
1 parent e76d66c commit 730d050

1 file changed

Lines changed: 161 additions & 30 deletions

File tree

‎PWGLF/TableProducer/Common/zdc2stagecalibration.cxx‎

Lines changed: 161 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,9 @@ struct zdc2stagecalibration {
7676
Configurable<bool> useSpatialCalib{"useSpatialCalib", false, "use ZDC spatial calibration"};
7777
Configurable<std::string> confSpatialPath{"confSpatialPath", "", "CCDB path to ZDC spatial calibration"};
7878
Configurable<bool> deriveSpatialCalib{"deriveSpatialCalib", false, "store spatial regression moments after gain correction"};
79+
Configurable<bool> deriveQRecentering{"deriveQRecentering", false, "store Q-vector recentering regression moments after gain and spatial correction"};
80+
Configurable<bool> useQRecentering{"useQRecentering", false, "apply Q-vector recentering calibration"};
81+
Configurable<std::string> confQRecenteringPath{"confQRecenteringPath", "", "CCDB path to the run-wise Q-vector recentering coefficients"};
7982

8083
Configurable<int> cfgSpatialCalibBins{
8184
"cfgSpatialCalibBins", 30,
@@ -94,6 +97,12 @@ struct zdc2stagecalibration {
9497
uint64_t eorTimestamp = 0;
9598
TH2D* gainprofile = nullptr;
9699
TH3D* spatialprofile = nullptr;
100+
TH2D* qRecenteringProfile = nullptr;
101+
102+
static constexpr int kQRecenteringNFeatures = 21;
103+
static constexpr int kQRecenteringNMatrixMoments = kQRecenteringNFeatures * (kQRecenteringNFeatures + 1) / 2;
104+
static constexpr int kQRecenteringNComponents = 4;
105+
static constexpr int kQRecenteringNMoments = kQRecenteringNMatrixMoments + kQRecenteringNComponents * kQRecenteringNFeatures;
97106

98107
void init(o2::framework::InitContext&)
99108
{
@@ -104,10 +113,16 @@ struct zdc2stagecalibration {
104113
AxisSpec towerAxis = {4, 0.0, 4.0, "tower"};
105114
AxisSpec crossAxis = {120, -1.0, 1.0, "cross asymmetry"};
106115
AxisSpec momentAxis = {15, 0.0, 15.0, "regression moment"};
116+
AxisSpec spatialMomentAxis = {6, 0.0, 6.0, "normalized pair regression moment"};
107117
AxisSpec spatialCalibAxis = {cfgSpatialCalibBins.value, -1.0, 1.0, "cross asymmetry"};
108118
AxisSpec phiAxis = {72, -3.14159265358979323846, 3.14159265358979323846, "#phi"};
109119
AxisSpec crossCoordAxis = {2, 0.0, 2.0, "coordinate"};
110120
AxisSpec ratioAxis = {240, 0.0, 2.4, "#Sigma tower/common"};
121+
AxisSpec qRecenteringMomentAxis = {kQRecenteringNMoments, 0.0, static_cast<double>(kQRecenteringNMoments), "regression moment"};
122+
AxisSpec centralityAxis = {80, 0.0, 80.0, "centrality (%)"};
123+
AxisSpec vertexXYAxis = {100, -0.5, 0.5, "vertex x/y (cm)"};
124+
AxisSpec vertexZAxis = {100, -10.0, 10.0, "vertex z (cm)"};
125+
AxisSpec qComponentAxis = {4, 0.0, 4.0, "Q component"};
111126

112127
histos.add("hEvtSelInfo", "hEvtSelInfo", kTH1F, {{10, 0.0, 10.0}});
113128
auto hEvtSelInfo = histos.get<TH1>(HIST("hEvtSelInfo"));
@@ -150,8 +165,31 @@ struct zdc2stagecalibration {
150165
// 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.
151166
histos.add("GainCalibration/hGainMomentsZNA", "ZNA linear-regression moments;time from SOR (h);moment index", kTH2D, {timeAxis, momentAxis});
152167
histos.add("GainCalibration/hGainMomentsZNC", "ZNC linear-regression moments;time from SOR (h);moment index", kTH2D, {timeAxis, momentAxis});
153-
histos.add("SpatialCalibration/hSpatialMomentsZNA", "ZNA spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, momentAxis});
154-
histos.add("SpatialCalibration/hSpatialMomentsZNC", "ZNC spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, momentAxis});
168+
// Normalized two-pair spatial regression moments. For each (Xcross,Ycross) cell:
169+
// u=(T1+T4)/C, v=(T2+T3)/C, and moments are u^2, u*v, v^2, u, v, N.
170+
histos.add("SpatialCalibration/hSpatialMomentsZNA", "ZNA normalized pair spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, spatialMomentAxis});
171+
histos.add("SpatialCalibration/hSpatialMomentsZNC", "ZNC normalized pair spatial regression moments;X cross;Y cross;moment index", kTH3D, {spatialCalibAxis, spatialCalibAxis, spatialMomentAxis});
172+
173+
// Q-recentering regression moments. Feature order:
174+
// 0:1, 1:C, 2:T, 3:Vx, 4:Vy, 5:Vz, 6:C2, 7:T2, 8:Vx2, 9:Vy2, 10:Vz2,
175+
// 11:C*T, 12:C*Vx, 13:C*Vy, 14:C*Vz, 15:T*Vx, 16:T*Vy, 17:T*Vz,
176+
// 18:Vx*Vy, 19:Vx*Vz, 20:Vy*Vz.
177+
// Moments 0..230 are the upper triangle of sum(X_i X_j).
178+
// Then 4 blocks of 21 sum(X_i Q): QxA, QyA, QxC, QyC.
179+
histos.add("QRecenteringCalibration/hRegressionMoments", "Q recentering regression moments;moment index;sum", kTH1D, {qRecenteringMomentAxis});
180+
181+
histos.add("QRecenteringQA/pQBeforeVsCentrality", "Q before recentering vs centrality;centrality (%);Q component;<Q>", kTProfile2D, {centralityAxis, qComponentAxis});
182+
histos.add("QRecenteringQA/pQBeforeVsTime", "Q before recentering vs time;time from SOR (h);Q component;<Q>", kTProfile2D, {timeAxis, qComponentAxis});
183+
histos.add("QRecenteringQA/pQBeforeVsVx", "Q before recentering vs v_{x};v_{x} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
184+
histos.add("QRecenteringQA/pQBeforeVsVy", "Q before recentering vs v_{y};v_{y} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
185+
histos.add("QRecenteringQA/pQBeforeVsVz", "Q before recentering vs v_{z};v_{z} (cm);Q component;<Q>", kTProfile2D, {vertexZAxis, qComponentAxis});
186+
histos.add("QRecenteringQA/pQAfterVsCentrality", "Q after recentering vs centrality;centrality (%);Q component;<Q>", kTProfile2D, {centralityAxis, qComponentAxis});
187+
histos.add("QRecenteringQA/pQAfterVsTime", "Q after recentering vs time;time from SOR (h);Q component;<Q>", kTProfile2D, {timeAxis, qComponentAxis});
188+
histos.add("QRecenteringQA/pQAfterVsVx", "Q after recentering vs v_{x};v_{x} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
189+
histos.add("QRecenteringQA/pQAfterVsVy", "Q after recentering vs v_{y};v_{y} (cm);Q component;<Q>", kTProfile2D, {vertexXYAxis, qComponentAxis});
190+
histos.add("QRecenteringQA/pQAfterVsVz", "Q after recentering vs v_{z};v_{z} (cm);Q component;<Q>", kTProfile2D, {vertexZAxis, qComponentAxis});
191+
histos.add("PhiQA/hPhiAfterRecenteringZNA", "ZNA #phi after Q recentering;#phi;events", kTH1F, {phiAxis});
192+
histos.add("PhiQA/hPhiAfterRecenteringZNC", "ZNC #phi after Q recentering;#phi;events", kTH1F, {phiAxis});
155193

156194
ccdb->setURL(cfgCcdbParam.cfgURL);
157195
ccdb->setCaching(true);
@@ -172,6 +210,7 @@ struct zdc2stagecalibration {
172210
bcSOR = static_cast<int64_t>(runInfo.orbitSOR) * static_cast<int64_t>(o2::constants::lhc::LHCMaxBunches);
173211
gainprofile = nullptr;
174212
spatialprofile = nullptr;
213+
qRecenteringProfile = nullptr;
175214
if (calibrationStage.value == 2 && useGainCallib.value && !confGainPath.value.empty()) {
176215
gainprofile = ccdb->getForTimeStamp<TH2D>(confGainPath.value, timestamp);
177216
if (!gainprofile) {
@@ -185,6 +224,12 @@ struct zdc2stagecalibration {
185224
LOGF(warn, "No ZDC spatial calibration found for run %d at timestamp %llu", runNumber, static_cast<unsigned long long>(timestamp));
186225
}
187226
}
227+
if (calibrationStage.value == 2 && useQRecentering.value && !confQRecenteringPath.value.empty()) {
228+
qRecenteringProfile = ccdb->getForTimeStamp<TH2D>(confQRecenteringPath.value, timestamp);
229+
if (!qRecenteringProfile) {
230+
LOGF(warn, "No ZDC Q-recentering calibration found for run %d at timestamp %llu", runNumber, static_cast<unsigned long long>(timestamp));
231+
}
232+
}
188233
}
189234

190235
float getTimeFromSOR(uint64_t globalBC) const
@@ -213,6 +258,24 @@ struct zdc2stagecalibration {
213258
return gain;
214259
}
215260

261+
std::array<double, kQRecenteringNFeatures> makeQRecenteringFeatures(float centrality, float timeFromSOR, float vx, float vy, float vz) const
262+
{
263+
// Fixed scaling keeps the polynomial basis numerically well behaved and must be
264+
// used identically when solving/applying the CCDB coefficients.
265+
const double c = (static_cast<double>(centrality) - 40.0) / 40.0;
266+
const double runHours = (eorTimestamp > sorTimestamp) ? static_cast<double>(eorTimestamp - sorTimestamp) * 1.e-3 / 3600.0 : static_cast<double>(cfgMaxRunHours.value);
267+
const double t = (runHours > 0.0) ? (2.0 * static_cast<double>(timeFromSOR) / runHours - 1.0) : 0.0;
268+
const double x = static_cast<double>(vx) / 0.1;
269+
const double y = static_cast<double>(vy) / 0.1;
270+
const double z = static_cast<double>(vz) / 10.0;
271+
272+
return {1.0, c, t, x, y, z,
273+
c * c, t * t, x * x, y * y, z * z,
274+
c * t, c * x, c * y, c * z,
275+
t * x, t * y, t * z,
276+
x * y, x * z, y * z};
277+
}
278+
216279
using MyCollisions = o2::soa::Join<o2::aod::Collisions, o2::aod::EvSels, o2::aod::Mults, o2::aod::FT0sCorrected, o2::aod::CentFT0Cs>;
217280

218281
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 {
373436
const double crossLookupYC = (zncCorr[2] - zncCorr[1]) / (zncCorr[2] + zncCorr[1]);
374437
if (calibrationStage.value == 2 && useGainCallib.value && deriveSpatialCalib.value && !useSpatialCalib.value) {
375438

376-
std::array<double, 15> spatialMomentsA{};
377-
std::array<double, 15> spatialMomentsC{};
378-
379-
int moment = 0;
380-
381-
for (int i = 0; i < 4; ++i) {
382-
for (int j = i; j < 4; ++j) {
383-
spatialMomentsA[moment] = znaCorr[i] * znaCorr[j];
384-
spatialMomentsC[moment] = zncCorr[i] * zncCorr[j];
385-
++moment;
439+
const double uA = (znaCorr[0] + znaCorr[3]) / static_cast<double>(znaEnergycommon);
440+
const double vA = (znaCorr[1] + znaCorr[2]) / static_cast<double>(znaEnergycommon);
441+
const double uC = (zncCorr[0] + zncCorr[3]) / static_cast<double>(zncEnergycommon);
442+
const double vC = (zncCorr[1] + zncCorr[2]) / static_cast<double>(zncEnergycommon);
443+
444+
if (std::isfinite(uA) && std::isfinite(vA) && std::isfinite(uC) && std::isfinite(vC)) {
445+
const std::array<double, 6> spatialMomentsA = {
446+
uA * uA,
447+
uA * vA,
448+
vA * vA,
449+
uA,
450+
vA,
451+
1.0};
452+
const std::array<double, 6> spatialMomentsC = {
453+
uC * uC,
454+
uC * vC,
455+
vC * vC,
456+
uC,
457+
vC,
458+
1.0};
459+
460+
for (int i = 0; i < 6; ++i) {
461+
histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNA"), crossLookupXA, crossLookupYA, i + 0.5, spatialMomentsA[i]);
462+
histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNC"), crossLookupXC, crossLookupYC, i + 0.5, spatialMomentsC[i]);
386463
}
387464
}
388-
389-
for (int i = 0; i < 4; ++i) {
390-
spatialMomentsA[10 + i] =
391-
static_cast<double>(znaEnergycommon) * znaCorr[i];
392-
393-
spatialMomentsC[10 + i] =
394-
static_cast<double>(zncEnergycommon) * zncCorr[i];
395-
}
396-
397-
spatialMomentsA[14] = 1.0;
398-
spatialMomentsC[14] = 1.0;
399-
400-
for (int i = 0; i < 15; ++i) {
401-
histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNA"), crossLookupXA, crossLookupYA, i + 0.5, spatialMomentsA[i]);
402-
histos.fill(HIST("SpatialCalibration/hSpatialMomentsZNC"), crossLookupXC, crossLookupYC, i + 0.5, spatialMomentsC[i]);
403-
}
404465
}
405466

406467
bool spatialOK = true;
@@ -533,9 +594,79 @@ struct zdc2stagecalibration {
533594
fillTable(false, 0.f, 0.f, 0.f, 0.f);
534595
return;
535596
}
536-
histos.fill(HIST("hEvtSelInfo"), 7.5);
537597

538-
fillTable(true, static_cast<float>(qxZDCA), static_cast<float>(qxZDCC), static_cast<float>(qyZDCA), static_cast<float>(qyZDCC));
598+
const auto qFeatures = makeQRecenteringFeatures(centrality, timeFromSOR, vx, vy, vz);
599+
std::array<double, 4> qValues = {qxZDCA, qyZDCA, qxZDCC, qyZDCC};
600+
601+
for (int i = 0; i < 4; ++i) {
602+
const double component = i + 0.5;
603+
histos.fill(HIST("QRecenteringQA/pQBeforeVsCentrality"), centrality, component, qValues[i]);
604+
histos.fill(HIST("QRecenteringQA/pQBeforeVsTime"), timeFromSOR, component, qValues[i]);
605+
histos.fill(HIST("QRecenteringQA/pQBeforeVsVx"), vx, component, qValues[i]);
606+
histos.fill(HIST("QRecenteringQA/pQBeforeVsVy"), vy, component, qValues[i]);
607+
histos.fill(HIST("QRecenteringQA/pQBeforeVsVz"), vz, component, qValues[i]);
608+
}
609+
610+
if (calibrationStage.value == 2 && useGainCallib.value && useSpatialCalib.value && deriveQRecentering.value && !useQRecentering.value) {
611+
int moment = 0;
612+
for (int i = 0; i < kQRecenteringNFeatures; ++i) {
613+
for (int j = i; j < kQRecenteringNFeatures; ++j) {
614+
histos.fill(HIST("QRecenteringCalibration/hRegressionMoments"), moment + 0.5, qFeatures[i] * qFeatures[j]);
615+
++moment;
616+
}
617+
}
618+
for (int component = 0; component < kQRecenteringNComponents; ++component) {
619+
for (int i = 0; i < kQRecenteringNFeatures; ++i) {
620+
const int index = kQRecenteringNMatrixMoments + component * kQRecenteringNFeatures + i;
621+
histos.fill(HIST("QRecenteringCalibration/hRegressionMoments"), index + 0.5, qFeatures[i] * qValues[component]);
622+
}
623+
}
624+
}
625+
626+
bool qRecenteringOK = true;
627+
if (calibrationStage.value == 2 && useQRecentering.value) {
628+
if (!qRecenteringProfile) {
629+
qRecenteringOK = false;
630+
} else {
631+
for (int component = 0; component < kQRecenteringNComponents; ++component) {
632+
double predictedBias = 0.0;
633+
for (int i = 0; i < kQRecenteringNFeatures; ++i) {
634+
const double coefficient = qRecenteringProfile->GetBinContent(qRecenteringProfile->FindBin(i + 0.5, component + 0.5));
635+
if (!std::isfinite(coefficient)) {
636+
qRecenteringOK = false;
637+
break;
638+
}
639+
predictedBias += coefficient * qFeatures[i];
640+
}
641+
if (!qRecenteringOK || !std::isfinite(predictedBias)) {
642+
qRecenteringOK = false;
643+
break;
644+
}
645+
qValues[component] -= predictedBias;
646+
}
647+
}
648+
}
649+
if (!qRecenteringOK) {
650+
fillTable(false, 0.f, 0.f, 0.f, 0.f);
651+
return;
652+
}
653+
654+
for (int i = 0; i < 4; ++i) {
655+
const double component = i + 0.5;
656+
histos.fill(HIST("QRecenteringQA/pQAfterVsCentrality"), centrality, component, qValues[i]);
657+
histos.fill(HIST("QRecenteringQA/pQAfterVsTime"), timeFromSOR, component, qValues[i]);
658+
histos.fill(HIST("QRecenteringQA/pQAfterVsVx"), vx, component, qValues[i]);
659+
histos.fill(HIST("QRecenteringQA/pQAfterVsVy"), vy, component, qValues[i]);
660+
histos.fill(HIST("QRecenteringQA/pQAfterVsVz"), vz, component, qValues[i]);
661+
}
662+
663+
const double phiRecenteringA = std::atan2(qValues[1], qValues[0]);
664+
const double phiRecenteringC = std::atan2(qValues[3], qValues[2]);
665+
histos.fill(HIST("PhiQA/hPhiAfterRecenteringZNA"), phiRecenteringA);
666+
histos.fill(HIST("PhiQA/hPhiAfterRecenteringZNC"), phiRecenteringC);
667+
668+
histos.fill(HIST("hEvtSelInfo"), 7.5);
669+
fillTable(true, static_cast<float>(qValues[0]), static_cast<float>(qValues[2]), static_cast<float>(qValues[1]), static_cast<float>(qValues[3]));
539670
histos.fill(HIST("hEvtSelInfo"), 8.5);
540671
}
541672
};

0 commit comments

Comments
 (0)