[MFT] Use residual instead of measurement to make Mille records (#10726)

arakotoz · web-flow · commit b503777a722e · 2023-02-08T18:12:09.000+01:00
* test with residual instead of measurement used to make records

* also store Pede errors and pulls on alignment corrections

* Use ClassDefOverride to avoid countless warnings in ROOT interpreter

* Forgot to use ClassDefOverride here
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/Aligner.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/Aligner.h
@@ -54,6 +54,9 @@ class Aligner
   void setAllowedVariationDeltaRz(const double value) { mAllowVar[2] = value; }
   void setChi2CutFactor(const double value) { mStartFac = value; }
 
+  /// \brief return the number of DOF per sensor
+  int getNDofPerSensor() const { return mNDofPerSensor; }
+
  protected:
   static constexpr int mNumberOfTrackParam = 4;                                  ///< Number of track (= local) parameters (X0, Tx, Y0, Ty)
   static constexpr int mNDofPerSensor = 4;                                       ///< translation in global x, y, z, and rotation Rz around global z-axis
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSparse.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSparse.h
@@ -77,7 +77,7 @@ class MatrixSparse : public MatrixSq
  protected:
   VectorSparse** fVecs = nullptr;
 
-  ClassDef(MatrixSparse, 0);
+  ClassDefOverride(MatrixSparse, 0);
 };
 
 //___________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSq.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSq.h
@@ -144,7 +144,7 @@ class MatrixSq : public TMatrixDBase
  protected:
   Bool_t fSymmetric; ///< is the matrix symmetric? Only lower triangle is filled
 
-  ClassDef(MatrixSq, 1);
+  ClassDefOverride(MatrixSq, 1);
 };
 
 //___________________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePede2.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePede2.h
@@ -181,7 +181,9 @@ class MillePede2
   void SetSigmaPar(int i, double par);
 
   /// \brief performs a requested number of global iterations
-  int GlobalFit(double* par = nullptr, double* error = nullptr, double* pull = nullptr);
+  int GlobalFit(std::vector<double>& par,
+                std::vector<double>& error,
+                std::vector<double>& pull);
 
   /// \brief perform global parameters fit once all the local equations have been fitted
   int GlobalFitIteration();
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePedeRecord.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePedeRecord.h
@@ -127,7 +127,7 @@ class MillePedeRecord : public TObject
   Double32_t* fValue; ///< [fSize] array of values: derivs,residuals
   Double32_t fWeight; ///< global weight for the record
 
-  ClassDef(MillePedeRecord, 3);
+  ClassDefOverride(MillePedeRecord, 3);
 };
 
 //_____________________________________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MinResSolve.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MinResSolve.h
@@ -132,7 +132,7 @@ class MinResSolve : public TObject
   MatrixSparse* fMatU; // aux. space
   SymBDMatrix* fMatBD; // aux. space
 
-  ClassDef(MinResSolve, 0);
+  ClassDefOverride(MinResSolve, 0);
 };
 
 } // namespace mft
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RecordsToAlignParams.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RecordsToAlignParams.h
@@ -16,6 +16,7 @@
 #ifndef ALICEO2_MFT_RECORDS_TO_ALIGN_PARAMS_H
 #define ALICEO2_MFT_RECORDS_TO_ALIGN_PARAMS_H
 
+#include <vector>
 #include <TChain.h>
 
 #include "MFTAlignment/MillePede2.h"
@@ -53,6 +54,15 @@ class RecordsToAlignParams : public Aligner
   /// \brief provide access to the AlignParam vector
   void getAlignParams(std::vector<o2::detectors::AlignParam>& alignParams) { alignParams = mAlignParams; }
 
+  /// \brief provide access to the vector of alignment corrections
+  void getPedeOutParams(std::vector<double>& output) { output = mPedeOutParams; }
+
+  /// \brief provide access to the vector of errors on the alignement corrections
+  void getPedeOutParamsErrors(std::vector<double>& output) { output = mPedeOutParamsErrors; }
+
+  /// \brief provide access to the vector of pulls on the alignement corrections
+  void getPedeOutParamsPulls(std::vector<double>& output) { output = mPedeOutParamsPulls; }
+
   /// \brief connect data record reader to input TChain of records
   void connectRecordReaderToChain(TChain* ch);
 
@@ -67,8 +77,11 @@ class RecordsToAlignParams : public Aligner
   bool mWithConstraintsRecReader;                      ///< boolean to set to true if one wants to also read constraints records
   o2::mft::MilleRecordReader* mConstraintsRecReader;   ///< utility that handles the reading of the constraints records
   o2::mft::MillePede2* mMillepede;                     ///< Millepede2 implementation copied from AliROOT
+  std::vector<double> mPedeOutParams;                  ///< Vector to store the outputs (alignment corrections) of the MillePede simulatenous fit
+  std::vector<double> mPedeOutParamsErrors;            ///< Vector to store the outputs (errors on the alignement corrections) of the MillePede simulatenous fit
+  std::vector<double> mPedeOutParamsPulls;             ///< Vector to store the outputs (pulls on the alignement corrections) of the MillePede simulatenous fit
 
-  ClassDef(RecordsToAlignParams, 0);
+  ClassDefOverride(RecordsToAlignParams, 0);
 };
 
 } // namespace mft
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RectMatrix.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RectMatrix.h
@@ -67,7 +67,7 @@ class RectMatrix : public TObject
   Int_t fNCols;     ///< Number of columns
   Double_t** fRows; ///< pointers on rows
 
-  ClassDef(RectMatrix, 0);
+  ClassDefOverride(RectMatrix, 0);
 };
 
 //___________________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymBDMatrix.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymBDMatrix.h
@@ -106,7 +106,7 @@ class SymBDMatrix : public MatrixSq
  protected:
   Double_t* fElems; ///< Elements booked by constructor
 
-  ClassDef(SymBDMatrix, 0);
+  ClassDefOverride(SymBDMatrix, 0);
 };
 
 //___________________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymMatrix.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymMatrix.h
@@ -215,7 +215,7 @@ class SymMatrix : public MatrixSq
   static SymMatrix* fgBuffer; ///< buffer for fast solution
   static Int_t fgCopyCnt;     ///< matrix copy counter
 
-  ClassDef(SymMatrix, 0);
+  ClassDefOverride(SymMatrix, 0);
 };
 
 //___________________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/TracksToRecords.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/TracksToRecords.h
@@ -144,7 +144,7 @@ class TracksToRecords : public Aligner
   /// \brief set the last component of the local equation vector for a given alignment point
   bool setLocalEquationZ();
 
-  ClassDef(TracksToRecords, 0);
+  ClassDefOverride(TracksToRecords, 0);
 };
 
 } // namespace mft
diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/VectorSparse.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/VectorSparse.h
@@ -82,7 +82,7 @@ class VectorSparse : public TObject
   UShort_t* fIndex; ///< Index of stored elems
   Double_t* fElems; ///< pointer on elements
 
-  ClassDef(VectorSparse, 0);
+  ClassDefOverride(VectorSparse, 0);
 };
 
 //___________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/src/MillePede2.cxx b/Detectors/ITSMFT/MFT/alignment/src/MillePede2.cxx
@@ -875,7 +875,9 @@ int MillePede2::LocalFit(std::vector<double>& localParams)
 }
 
 //_____________________________________________________________________________
-int MillePede2::GlobalFit(double* par, double* error, double* pull)
+int MillePede2::GlobalFit(std::vector<double>& par,
+                          std::vector<double>& error,
+                          std::vector<double>& pull)
 {
   if (fRecordReader == nullptr) {
     LOG(fatal) << "MillePede2::GlobalFit() - aborted, input record reader is a null pointer";
@@ -911,19 +913,19 @@ int MillePede2::GlobalFit(double* par, double* error, double* pull)
     return 0;
   }
 
-  if (par) {
+  if (par.size()) {
     for (int i = fNGloParIni; i--;) {
       par[i] = GetFinalParam(i);
     }
   }
 
   if (fGloSolveStatus == kInvert) { // errors on params are available
-    if (error) {
+    if (error.size()) {
       for (int i = fNGloParIni; i--;) {
         error[i] = GetFinalError(i);
       }
     }
-    if (pull) {
+    if (pull.size()) {
       for (int i = fNGloParIni; i--;) {
         pull[i] = GetPull(i);
       }
diff --git a/Detectors/ITSMFT/MFT/alignment/src/RecordsToAlignParams.cxx b/Detectors/ITSMFT/MFT/alignment/src/RecordsToAlignParams.cxx
@@ -36,6 +36,14 @@ RecordsToAlignParams::RecordsToAlignParams()
   if (mWithConstraintsRecReader) {
     mConstraintsRecReader = new MilleRecordReader();
   }
+
+  // allocate memory in vectors to store the results of the global fit
+
+  double default_value = 0;
+  mPedeOutParams = std::vector<double>(mNumberOfGlobalParam, default_value);
+  mPedeOutParamsErrors = std::vector<double>(mNumberOfGlobalParam, default_value);
+  mPedeOutParamsPulls = std::vector<double>(mNumberOfGlobalParam, default_value);
+
   LOGF(debug, "RecordsToAlignParams instantiated");
 }
 
@@ -117,23 +125,9 @@ void RecordsToAlignParams::globalFit()
     mMillepede->InitChi2Storage(mNEntriesAutoSave);
   }
 
-  // allocate memory in arrays to temporarily store the results of the global fit
-
-  double* params = (double*)malloc(sizeof(double) * mNumberOfGlobalParam);
-  double* paramsErrors = (double*)malloc(sizeof(double) * mNumberOfGlobalParam);
-  double* paramsPulls = (double*)malloc(sizeof(double) * mNumberOfGlobalParam);
-
-  // initialise the content of each array
-
-  for (int ii = 0; ii < mNumberOfGlobalParam; ii++) {
-    params[ii] = 0.;
-    paramsErrors[ii] = 0.;
-    paramsPulls[ii] = 0.;
-  }
-
   // perform the simultaneous fit of track and alignement parameters
 
-  mMillepede->GlobalFit(params, paramsErrors, paramsPulls);
+  mMillepede->GlobalFit(mPedeOutParams, mPedeOutParamsErrors, mPedeOutParamsPulls);
 
   if (mWithControl) {
     mMillepede->EndChi2Storage();
@@ -154,10 +148,10 @@ void RecordsToAlignParams::globalFit()
   for (int chipId = 0; chipId < mNumberOfSensors; chipId++) {
     chipHelper.setSensorOnlyInfo(chipId);
     std::stringstream name = chipHelper.getSensorFullName(withSymName);
-    dx = params[chipId * mNDofPerSensor + 0];
-    dy = params[chipId * mNDofPerSensor + 1];
-    dz = params[chipId * mNDofPerSensor + 3];
-    dRz = params[chipId * mNDofPerSensor + 2];
+    dx = mPedeOutParams[chipId * mNDofPerSensor + 0];
+    dy = mPedeOutParams[chipId * mNDofPerSensor + 1];
+    dz = mPedeOutParams[chipId * mNDofPerSensor + 3];
+    dRz = mPedeOutParams[chipId * mNDofPerSensor + 2];
     LOGF(info,
          "%s, %.3e, %.3e, %.3e, %.3e",
          name.str().c_str(), dx, dy, dz, dRz);
@@ -168,12 +162,6 @@ void RecordsToAlignParams::globalFit()
       dRx, dRy, dRz,
       global);
   }
-
-  // free allocated memory
-
-  free(params);
-  free(paramsErrors);
-  free(paramsPulls);
 }
 
 //__________________________________________________________________________
diff --git a/Detectors/ITSMFT/MFT/alignment/src/TracksToRecords.cxx b/Detectors/ITSMFT/MFT/alignment/src/TracksToRecords.cxx
@@ -568,13 +568,13 @@ bool TracksToRecords::setLocalEquationX()
            mGlobalDerivatives[chipId * mNDofPerSensor + 1],
            mGlobalDerivatives[chipId * mNDofPerSensor + 2],
            mGlobalDerivatives[chipId * mNDofPerSensor + 3],
-           mAlignPoint->getLocalMeasuredPosition().X(),
+           mAlignPoint->getLocalResidual().X(),
            mAlignPoint->getLocalMeasuredPositionSigma().X());
     }
     mMillepede->SetLocalEquation(
       mGlobalDerivatives,
       mLocalDerivatives,
-      mAlignPoint->getLocalMeasuredPosition().X(),
+      mAlignPoint->getLocalResidual().X(),
       mAlignPoint->getLocalMeasuredPositionSigma().X());
   } else {
     mCounterLocalEquationFailed++;
@@ -632,13 +632,13 @@ bool TracksToRecords::setLocalEquationY()
            mGlobalDerivatives[chipId * mNDofPerSensor + 1],
            mGlobalDerivatives[chipId * mNDofPerSensor + 2],
            mGlobalDerivatives[chipId * mNDofPerSensor + 3],
-           mAlignPoint->getLocalMeasuredPosition().Y(),
+           mAlignPoint->getLocalResidual().Y(),
            mAlignPoint->getLocalMeasuredPositionSigma().Y());
     }
     mMillepede->SetLocalEquation(
       mGlobalDerivatives,
       mLocalDerivatives,
-      mAlignPoint->getLocalMeasuredPosition().Y(),
+      mAlignPoint->getLocalResidual().Y(),
       mAlignPoint->getLocalMeasuredPositionSigma().Y());
   } else {
     mCounterLocalEquationFailed++;
@@ -696,13 +696,13 @@ bool TracksToRecords::setLocalEquationZ()
            mGlobalDerivatives[chipId * mNDofPerSensor + 1],
            mGlobalDerivatives[chipId * mNDofPerSensor + 2],
            mGlobalDerivatives[chipId * mNDofPerSensor + 3],
-           mAlignPoint->getLocalMeasuredPosition().Z(),
+           mAlignPoint->getLocalResidual().Z(),
            mAlignPoint->getLocalMeasuredPositionSigma().Z());
     }
     mMillepede->SetLocalEquation(
       mGlobalDerivatives,
       mLocalDerivatives,
-      mAlignPoint->getLocalMeasuredPosition().Z(),
+      mAlignPoint->getLocalResidual().Z(),
       mAlignPoint->getLocalMeasuredPositionSigma().Z());
   } else {
     mCounterLocalEquationFailed++;
