ITS alignment monitoring (check with QCv20240903), revised 20240919

Modules/ITS/include/ITS/ITSTrackTask.h

@@ -25,6 +25,13 @@
 #include <TTree.h>
 #include "Common/TH1Ratio.h"
 #include "Common/TH2Ratio.h"
+#include <Fit/Fitter.h>
+#include <Math/Functor.h>
+#include <TMatrixD.h>
+#include <TVector3.h>
+#include "TMath.h"
+#include <cassert>
+#include "ITS/ITSHelpers.h"
 
 class TH1;
 class TH2;
@@ -118,6 +125,152 @@ class ITSTrackTask : public TaskInterface
   double ptBins[141];     // pt bins
 
   o2::itsmft::TopologyDictionary* mDict;
+
+ private:
+  // analysis for its-only residual
+  o2::its::GeometryTGeo* mGeom;
+
+  std::vector<double> FitStepSize{ 0.3, 1.0e-5, 1.0e-5, 1.0e-5 };
+
+  double FitTolerance = 1.0e-8;
+  double ITS_AbsBz = 0.5; // T
+
+  double inputG_C[3 * NLayer];
+  double FitparXY[6];
+  double FitparDZ[2];
+  ROOT::Fit::Fitter fitterA;
+  ROOT::Fit::Fitter fitterB;
+  Int_t mAlignmentMonitor = 0;
+  Int_t mDefaultMomResPar = 0;
+  Int_t mResMonNclMin = 0;
+  float mResMonTrackMinPt = 0;
+
+  std::array<std::unique_ptr<TH2F>, NLayer> hResidualXY{}; //[NLayer];
+  std::array<std::unique_ptr<TH2F>, NLayer> hResidualZD{}; //[NLayer];
+
+  void circleFitXY(double* input, double* par, double& MSEvalue, std::vector<bool> hitUpdate, int step = 0);
+
+  // default setting
+  //  function Object to be minimized
+  struct se_circlefitXY {
+    // the TGraph is a data member of the object
+    std::vector<TVector3> fHits;
+    double thetaR;
+    double Bz;
+    double sigma_meas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, { 40, 40, 40, 40, 40, 40, 40 } };    // um unit
+    double sigma_msc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, { 25, 25, 25, 75, 75, 75, 75 } }; // um unit
+
+    se_circlefitXY() = default;
+    se_circlefitXY(std::vector<TVector3>& h, double tR, double bz)
+    {
+      fHits = h;
+      thetaR = tR;
+      Bz = bz;
+    };
+
+    void loadparameters(double arrpar_meas[][NLayer], double arrpar_msc[][NLayer])
+    {
+      for (int a = 0; a < 2; a++) {
+        for (int l = 0; l < NLayer; l++) {
+          sigma_meas[a][l] = arrpar_meas[a][l];
+          sigma_msc[a][l] = arrpar_msc[a][l];
+        }
+      }
+    };
+
+    void init()
+    {
+      fHits.clear();
+      thetaR = 0;
+      Bz = 0;
+    };
+
+    void set(std::vector<TVector3>& h, double tR, double bz)
+    {
+      fHits = h;
+      thetaR = tR;
+      Bz = bz;
+    };
+
+    double getsigma(double R, int L, double B, int axis)
+    {
+      // R : cm
+      // B : T
+      if (L < 0 || L >= NLayer)
+        return 1;
+      double aL = sigma_meas[axis][L] * 1e-4; // um -> cm
+      double bL = sigma_msc[axis][L] * 1e-4;  // um -> cm
+      double Beff = 0.299792458 * B;
+      double sigma = std::sqrt(std::pow(aL, 2) + std::pow(bL, 2) / (std::pow(Beff, 2) * std::pow(R * 1e-2, 2)));
+
+      return sigma;
+    };
+
+    // calculate distance line-point
+    double distance2(double x, double y, const double* p, double tR, int charge)
+    {
+
+      double R = std::abs(1 / p[0]);
+
+      double Xc = p[0] > 0 ? R * std::cos(p[1] + tR + 0.5 * TMath::Pi()) : R * std::cos(p[1] + tR - 0.5 * TMath::Pi());
+      double Yc = p[0] > 0 ? R * std::sin(p[1] + tR + 0.5 * TMath::Pi()) : R * std::sin(p[1] + tR - 0.5 * TMath::Pi());
+
+      double dx = x - (Xc + p[2]);
+      double dy = y - (Yc + p[3]);
+      double dxy = R - std::sqrt(dx * dx + dy * dy);
+
+      double d2 = dxy * dxy;
+      return d2;
+    };
+
+    // implementation of the function to be minimized
+    double operator()(const double* par)
+    { // const double -> double
+      assert(fHits != 0);
+
+      int nhits = fHits.size();
+      double sum = 0;
+
+      double charge = par[0] > 0 ? +1 : -1;
+      double RecRadius = std::abs(1 / par[0]);
+
+      double Sigma_tot[NLayer];
+      double sum_Sigma_tot = 0;
+      for (int l = 0; l < nhits; l++) {
+        Sigma_tot[l] = getsigma(RecRadius, fHits[l].Z(), Bz, 1);
+        sum_Sigma_tot += 1 / (std::pow(Sigma_tot[l], 2));
+      }
+
+      for (int l = 0; l < nhits; l++) {
+        double d = distance2(fHits[l].X(), fHits[l].Y(), par, thetaR, charge) / (std::pow(Sigma_tot[l], 2));
+        sum += d;
+      }
+      sum = sum / sum_Sigma_tot;
+
+      return sum;
+    };
+  };
+
+  se_circlefitXY fitfuncXY;
+
+  void lineFitDZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector<bool> hitUpdate);
+
+  double mSigmaMeas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, { 40, 40, 40, 40, 40, 40, 40 } };    // um unit
+  double mSigmaMsc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, { 25, 25, 25, 75, 75, 75, 75 } }; // um unit
+
+  double getsigma(double R, int L, double B, int axis)
+  {
+    // R : cm
+    // B : T
+    if (L < 0 || L >= NLayer)
+      return 1;
+    double aL = mSigmaMeas[axis][L] * 1e-4; // um -> cm
+    double bL = mSigmaMsc[axis][L] * 1e-4;  // um -> cm
+    double Beff = 0.299792458 * B;
+    double sigma = std::sqrt(std::pow(aL, 2) + std::pow(bL, 2) / (std::pow(Beff, 2) * std::pow(R * 1e-2, 2)));
+
+    return sigma;
+  };
 };
 } // namespace o2::quality_control_modules::its
 

Modules/ITS/include/ITS/README.md

@@ -0,0 +1,31 @@
+# Analysis for ITS-only residual
+
+We use the Weighted Least Square(WLS) fit and calculate the residuals from each cluster positions to the fitted trajectory.
+
+## Circular Fit in the XY plane : circleFitXY
+
+In the XY plane, the trajectory of particle under the magnetic field in z-axis can be fitted with circular model.
+There are 4 fit parameters : radius(R), direction at origin($$\theta_R$$), estimated vertex X(vx), and estimated vertex Y(vy).
+
+$$ x_c = R \cdot cos(\theta_R + \frac{R}{|R|} \cdot \pi/2) $$
+
+$$ y_c = R \cdot sin(\theta_R + \frac{R}{|R|} \cdot \pi/2) $$
+
+$$ (x - (x_c + v_x))^2 + (y - (y_c + v_y))^2 = |R|^2 $$
+
+For the fit stability, 
+1. A simple algorithm to determine seeds of fit parameters(radius and direction at origin) is applied.
+2. Circular fit is performed for both (+) and (-) signs of the circular trajectory (even if we have approximated seeds from 1.), the one with the lower chi2(sum of weighted residuals) is finally adopted. This procedure ensures the stability in case of hard(high pT) tracks. 
+
+## Linear Fit in the DZ plane : lineFitDZ
+
+After 4 fit parameters of circular trajectory is determined, we can describe the fit positions of associated clusters as well as the global positions.
+This makes it possible to represent each cluster position in the poloar coordinates $$(r, \beta)$$, whose origin is $$(x_c, y_c)$$.
+
+One of the simplest residual in z direction can be calculated in the RZ plane, where $$r = \sqrt(x^2 + y^2)$$.
+However, in this case, the fitting stability cannot be guaranteed if the collision position is not located in the same quadrant.
+
+On the other hand, another parameter beta determined above the equation of the circle determined (in the XY plane) corresponds to the distance of the particle moving(=D) in the XY plane.
+In this case, the fitting stability can be guaranteed even if the collision position is not located in the same quadrant.
+
+So we proceed with the linear fit in the DZ plane.

Modules/ITS/itsTrack.json

@@ -45,9 +45,16 @@
           "nBCbins": "103",
           "dicttimestamp" : "0",
           "doNorm" : "1",
-          "InvMasses" : "0"
+          "InvMasses" : "0",
+          "doAlignmentMonitor" : "1",
+          "UseDefaultMomResPar" : "0",
+          "MomResParMEAS1": "45, 45, 45, 55, 55, 55, 55",
+          "MomResParMEAS2": "40, 40, 40, 40, 40, 40, 40",
+          "MomResParMSC1": "30, 30, 30, 110, 110, 110, 110",
+          "MomResParMSC2": "25, 25, 25, 75, 75, 75, 75",
+          "ResidualMonitorNclMin" : "5",
+          "ResidualMonitorTrackMinPt" : "0.00"
         }
-
       }
     },
     "checks" : {