@@ -61,6 +61,32 @@ using VTIndex = o2::dataformats::VtxTrackIndex;
6161using GTrackID = o2::dataformats::GlobalTrackID;
6262using TrackD = o2::track::TrackParCovD;
6363
64+ // compute normalized (u,v) in [-1,1] from global position on a sensor
65+ static std::pair<double , double > computeUV (double gloX, double gloY, double gloZ, int sensorID, double radius)
66+ {
67+ const bool isTop = sensorID % 2 == 0 ;
68+ const double phi = o2::math_utils::to02Pid (std::atan2 (gloY, gloX));
69+ const double phiBorder1 = o2::math_utils::to02Pid (((isTop ? 0 . : 1 .) * TMath::Pi ()) + std::asin (constants::equatorialGap / 2 . / radius));
70+ const double phiBorder2 = o2::math_utils::to02Pid (((isTop ? 1 . : 2 .) * TMath::Pi ()) - std::asin (constants::equatorialGap / 2 . / radius));
71+ const double u = (((phi - phiBorder1) * 2 .) / (phiBorder2 - phiBorder1)) - 1 .;
72+ const double v = ((2 . * gloZ + constants::segment::lengthSensitive) / constants::segment::lengthSensitive) - 1 .;
73+ return {u, v};
74+ }
75+
76+ // evaluate Legendre polynomials P_0(x) through P_order(x) via recurrence
77+ static std::vector<double > legendrePols (int order, double x)
78+ {
79+ std::vector<double > p (order + 1 );
80+ p[0 ] = 1 .;
81+ if (order > 0 ) {
82+ p[1 ] = x;
83+ }
84+ for (int n = 1 ; n < order; ++n) {
85+ p[n + 1 ] = ((2 * n + 1 ) * x * p[n] - n * p[n - 1 ]) / (n + 1 );
86+ }
87+ return p;
88+ }
89+
6490class AlignmentSpec final : public Task
6591{
6692 public:
@@ -129,8 +155,8 @@ class AlignmentSpec final : public Task
129155 GTrackID::mask_t mTracksSrc ;
130156 int mNThreads {1 };
131157 const AlignmentParams* mParams {nullptr };
132- std::array<o2::math_utils::Legendre2DPolynominal, 6 > mDeformations ; // one per sensorID (0-5)
133- std::array<Eigen::Matrix<double , 6 , 1 >, 6 > mRigidBodyParams ; // (dx,dy,dz,rx,ry,rz) in LOC per sensorID
158+ std::array<o2::math_utils::Legendre2DPolynominal, 6 > mDeformations ; // one per sensorID (0-5)
159+ std::array<Eigen::Matrix<double , 6 , 1 >, 6 > mRigidBodyParams ; // (dx,dy,dz,rx,ry,rz) in LOC per sensorID
134160};
135161
136162void AlignmentSpec::init (InitContext& ic)
@@ -295,37 +321,45 @@ void AlignmentSpec::process()
295321 // derivatives for all sensitive volumes and their parents
296322 // this is the derivative in TRK but we want to align in LOC
297323 // so dr/da_(LOC) = dr/da_(TRK) * da_(TRK)/da_(LOC)
298- const auto * sensorVol = mChip2Hiearchy .at (lbl);
324+ const auto * tileVol = mChip2Hiearchy .at (lbl);
299325 Matrix36 der = getRigidBodyDerivatives (wTrk);
300326
301- // count columns: only volumes with real DOFs (not DOFPseudo)
302- int nCol {0 };
303- for (const auto * v = sensorVol ; v && !v->isRoot (); v = v->getParent ()) {
327+ // count rigid body columns: only volumes with real DOFs (not DOFPseudo)
328+ int nColRB {0 };
329+ for (const auto * v = tileVol ; v && !v->isRoot (); v = v->getParent ()) {
304330 if (v->getRigidBodyDOF () != RigidBodyDOFPseudo) {
305- nCol += RigidBodyDOF::NDOF;
331+ nColRB += RigidBodyDOF::NDOF;
306332 }
307333 }
334+
335+ // count Legendre columns for ITS3 sensors (applied at sensor = tile's parent)
336+ const bool isITS3 = constants::detID::isDetITS3 (frame.sens );
337+ const int N = mParams ->legendreOrder ;
338+ const int nLeg = isITS3 ? (N + 1 ) * (N + 2 ) / 2 : 0 ;
339+ const int nCol = nColRB + nLeg;
340+
308341 std::vector<int > gLabels ;
309342 gLabels .reserve (nCol);
310343 Eigen::MatrixXd gDer (3 , nCol);
344+ gDer .setZero ();
311345 Eigen::Index curCol{0 };
312346
313- // 1) sensor : TRK -> LOC via precomputed T2L and J_L2T
347+ // 1) tile : TRK -> LOC via precomputed T2L and J_L2T
314348 const double posTrk[3 ] = {frame.x , 0 ., 0 .};
315349 double posLoc[3 ];
316- sensorVol ->getT2L ().LocalToMaster (posTrk, posLoc);
350+ tileVol ->getT2L ().LocalToMaster (posTrk, posLoc);
317351 Matrix66 jacL2T;
318- sensorVol ->computeJacobianL2T (posLoc, jacL2T);
352+ tileVol ->computeJacobianL2T (posLoc, jacL2T);
319353 der *= jacL2T;
320- if (sensorVol ->getRigidBodyDOF () != RigidBodyDOFPseudo) {
354+ if (tileVol ->getRigidBodyDOF () != RigidBodyDOFPseudo) {
321355 for (uint8_t iDOF{0 }; iDOF < RigidBodyDOF::NDOF; ++iDOF) {
322- gLabels .push_back (sensorVol ->getLabel ().rawGBL (iDOF));
356+ gLabels .push_back (tileVol ->getLabel ().rawGBL (iDOF));
323357 }
324358 gDer .middleCols (curCol++ * RigidBodyDOF::NDOF, RigidBodyDOF::NDOF) = der;
325359 }
326360
327361 // 2) chain through parents: child's J_L2P
328- for (const auto * child = sensorVol ; child->getParent () && !child->getParent ()->isRoot (); child = child->getParent ()) {
362+ for (const auto * child = tileVol ; child->getParent () && !child->getParent ()->isRoot (); child = child->getParent ()) {
329363 der *= child->getJL2P ();
330364 const auto * parent = child->getParent ();
331365 if (parent->getRigidBodyDOF () != RigidBodyDOFPseudo) {
@@ -335,6 +369,46 @@ void AlignmentSpec::process()
335369 gDer .middleCols (curCol++ * RigidBodyDOF::NDOF, RigidBodyDOF::NDOF) = der;
336370 }
337371 }
372+
373+ // 3) Legendre global derivatives for ITS3 sensors
374+ // h(u,v) = sum_{i=0}^{N} sum_{j=0}^{i} c_{ij} P_{i-j}(u) P_j(v)
375+ // dres/dc_{ij} = slope * P_{i-j}(u) * P_j(v)
376+ // Legendre DOFs live on the sensor volume (tile's direct parent)
377+ if (isITS3) {
378+ const int sensorID = constants::detID::getSensorID (frame.sens );
379+ const int layerID = constants::detID::getDetID2Layer (frame.sens );
380+ const auto * sensorVol = tileVol->getParent ();
381+
382+ // compute (u,v) from global position
383+ const double r = frame.x ;
384+ const double gX = r * std::cos (frame.alpha );
385+ const double gY = r * std::sin (frame.alpha );
386+ const double gZ = frame.positionTrackingFrame [1 ];
387+ auto [u, v] = computeUV (gX , gY , gZ , sensorID, constants::radii[layerID]);
388+
389+ // track slopes at this cluster (reconstructed track)
390+ const double snp = wTrk.getSnp ();
391+ const double tgl = wTrk.getTgl ();
392+ const double csci = 1 . / std::sqrt (1 . - (snp * snp));
393+ const double dydx = snp * csci;
394+ const double dzdx = tgl * csci;
395+
396+ // evaluate Legendre polynomials P_0..P_N for both u and v
397+ auto pu = legendrePols (N, u);
398+ auto pv = legendrePols (N, v);
399+
400+ int legIdx = 0 ;
401+ const int legColStart = nColRB;
402+ for (int i = 0 ; i <= N; ++i) {
403+ for (int j = 0 ; j <= i; ++j) {
404+ const double basis = pu[i - j] * pv[j];
405+ gLabels .push_back (sensorVol->getLabel ().asCalib ().rawGBL (legIdx));
406+ gDer (0 , legColStart + legIdx) = dydx * basis; // dres_y / dc_{ij}
407+ gDer (1 , legColStart + legIdx) = dzdx * basis; // dres_z / dc_{ij}
408+ ++legIdx;
409+ }
410+ }
411+ }
338412 point.addGlobals (gLabels , gDer );
339413 }
340414
@@ -771,16 +845,6 @@ bool AlignmentSpec::applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt&
771845 // --- Legendre deformation (non-rigid-body) ---
772846 if (mParams ->doMisalignment && mIT3Dict && mUseMC ) {
773847 const auto prop = o2::base::PropagatorD::Instance ();
774- // compute normalized (u,v) in [-1,1] from global position
775- auto computeUV = [](double gloX, double gloY, double gloZ, int sensorID, double radius) -> std::pair<double , double > {
776- const bool isTop = sensorID % 2 == 0 ;
777- const double phi = o2::math_utils::to02Pid (std::atan2 (gloY, gloX));
778- const double phiBorder1 = o2::math_utils::to02Pid (((isTop ? 0 . : 1 .) * TMath::Pi ()) + std::asin (constants::equatorialGap / 2 . / radius));
779- const double phiBorder2 = o2::math_utils::to02Pid (((isTop ? 1 . : 2 .) * TMath::Pi ()) - std::asin (constants::equatorialGap / 2 . / radius));
780- const double u = (((phi - phiBorder1) * 2 .) / (phiBorder2 - phiBorder1)) - 1 .;
781- const double v = ((2 . * gloZ + constants::segment::lengthSensitive) / constants::segment::lengthSensitive) - 1 .;
782- return {u, v};
783- };
784848
785849 const auto lbl = mRecoData ->getITSTracksMCLabels ()[iTrk];
786850 const auto mcTrk = mcReader->getTrack (lbl);
0 commit comments