Use stable lin.ref. point for alignment track initial fit

Author: shahor02

Detectors/Align/include/Align/AlignConfig.h

@@ -85,6 +85,7 @@ struct AlignConfig : public o2::conf::ConfigurableParamHelper<AlignConfig> {
   float controlFraction = -1.;     // fraction for which control output is requested, if negative - only 1st instance of device will write them
   float MPRecOutFraction = -1.;    // compact Millepede2Record fraction, if negative - only 1st instance of device will write them
 
+  bool useLinRef = true;      // use initial track for lienarization reference point
   bool MilleOut = true;       // Mille output
   bool KalmanResid = true;    // Kalman residuals
   bool MilleOutBin = true;    // text vs binary output for mille data

Detectors/Align/include/Align/AlignmentTrack.h

@@ -39,6 +39,7 @@ class AlignmentTrack : public trackParam_t, public TObject
 {
  public:
   using trackParam_t = o2::track::TrackParametrizationWithError<double>;
+  using trackPar_t = o2::track::TrackParametrization<double>;
   using PropagatorD = o2::base::PropagatorD;
   using MatCorrType = PropagatorD::MatCorrType;
   using GTrackID = o2::dataformats::GlobalTrackID;
@@ -83,9 +84,9 @@ class AlignmentTrack : public trackParam_t, public TObject
   //
   template <typename P>
   void copyFrom(const o2::track::TrackParametrizationWithError<P>& trc);
-  bool propagateToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tLT = nullptr, int signCorr = 0);
-  bool propagateParamToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0);             // param only
-  bool propagateParamToPoint(trackParam_t* trSet, int nTr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0); // params only
+  bool propagateToPoint(trackParam_t& tr, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tLT = nullptr, int signCorr = 0);
+  bool propagateParamToPoint(trackPar_t& tr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0);             // param only
+  bool propagateParamToPoint(trackPar_t* trSet, int nTr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0); // params only
   //
   bool calcResiduals(const double* params = nullptr);
   bool calcResidDeriv(double* params = nullptr);
@@ -119,23 +120,23 @@ class AlignmentTrack : public trackParam_t, public TObject
   void imposePtBOff(double pt) { setQ2Pt(1. / pt); }
   // propagation methods
   void copyFrom(const trackParam_t* etp);
-  bool applyMatCorr(trackParam_t& trPar, const double* corrDiag, const AlignmentPoint* pnt);
-  bool applyMatCorr(trackParam_t* trSet, int ntr, const double* corrDiaf, const AlignmentPoint* pnt);
-  bool applyMatCorr(trackParam_t& trPar, const double* corrpar);
+  bool applyMatCorr(trackPar_t& trPar, const double* corrDiag, const AlignmentPoint* pnt);
+  bool applyMatCorr(trackPar_t* trSet, int ntr, const double* corrDiaf, const AlignmentPoint* pnt);
+  bool applyMatCorr(trackPar_t& trPar, const double* corrpar);
   //
   double getResidual(int dim, int pntID) const { return mResid[dim][pntID]; }
   const double* getDResDLoc(int dim, int pntID) const { return mDResDLoc[dim].data() + (pntID * mNLocPar); }
   const double* getDResDGlo(int dim, int id) const { return mDResDGlo[dim].data() + id; }
   const int* getGloParID() const { return mGloParID.data(); }
   //
-  void setParams(trackParam_t& tr, double x, double alp, const double* par, bool add);
-  void setParams(trackParam_t* trSet, int ntr, double x, double alp, const double* par, bool add);
-  void setParam(trackParam_t& tr, int par, double val);
-  void setParam(trackParam_t* trSet, int ntr, int par, double val);
-  void modParam(trackParam_t& tr, int par, double delta);
-  void modParam(trackParam_t* trSet, int ntr, int par, double delta);
+  void setParams(trackPar_t& tr, double x, double alp, const double* par, bool add);
+  void setParams(trackPar_t* trSet, int ntr, double x, double alp, const double* par, bool add);
+  void setParam(trackPar_t& tr, int par, double val);
+  void setParam(trackPar_t* trSet, int ntr, int par, double val);
+  void modParam(trackPar_t& tr, int par, double delta);
+  void modParam(trackPar_t* trSet, int ntr, int par, double delta);
   //
-  void richardsonDeriv(const trackParam_t* trSet, const double* delta,
+  void richardsonDeriv(const trackPar_t* trSet, const double* delta,
                        const AlignmentPoint* pnt, double& derY, double& derZ);
   //
   const double* getLocPars() const { return mLocPar.data(); }
@@ -179,13 +180,14 @@ class AlignmentTrack : public trackParam_t, public TObject
   std::vector<double> mLocPar;              // local parameters array
   std::vector<int> mGloParID;               // IDs of relevant global params
  private:
-  bool propagate(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr = 0);
+  bool propagate(trackParam_t& tr, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr = 0);
+  bool propagate(trackPar_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr = 0);
   //
   ClassDefOverride(AlignmentTrack, 2)
 };
 
 //____________________________________________________________________________________________
-inline void AlignmentTrack::setParams(trackParam_t& tr, double x, double alp, const double* par, bool add)
+inline void AlignmentTrack::setParams(trackPar_t& tr, double x, double alp, const double* par, bool add)
 {
   // set track params
   const double kDefQ2PtCosm = 1;
@@ -205,7 +207,7 @@ inline void AlignmentTrack::setParams(trackParam_t& tr, double x, double alp, co
 }
 
 //____________________________________________________________________________________________
-inline void AlignmentTrack::setParams(trackParam_t* trSet, int ntr, double x, double alp, const double* par, bool add)
+inline void AlignmentTrack::setParams(trackPar_t* trSet, int ntr, double x, double alp, const double* par, bool add)
 {
   // set parames for multiple tracks (VECTORIZE THIS)
   if (!add) { // full parameter supplied
@@ -224,14 +226,14 @@ inline void AlignmentTrack::setParams(trackParam_t* trSet, int ntr, double x, do
 }
 
 //____________________________________________________________________________________________
-inline void AlignmentTrack::setParam(trackParam_t& tr, int par, double val)
+inline void AlignmentTrack::setParam(trackPar_t& tr, int par, double val)
 {
   // set track parameter
   tr.setParam(val, par);
 }
 
 //____________________________________________________________________________________________
-inline void AlignmentTrack::setParam(trackParam_t* trSet, int ntr, int par, double val)
+inline void AlignmentTrack::setParam(trackPar_t* trSet, int ntr, int par, double val)
 {
   // set parames for multiple tracks (VECTORIZE THIS)
   for (int i = 0; i < ntr; ++i) {
@@ -240,15 +242,15 @@ inline void AlignmentTrack::setParam(trackParam_t* trSet, int ntr, int par, doub
 }
 
 //____________________________________________________________________________________________
-inline void AlignmentTrack::modParam(trackParam_t& tr, int par, double delta)
+inline void AlignmentTrack::modParam(trackPar_t& tr, int par, double delta)
 {
   // modify track parameter
   const auto val = tr.getParam(par) + delta;
   setParam(tr, par, val);
 }
 
 //____________________________________________________________________________________________
-inline void AlignmentTrack::modParam(trackParam_t* trSet, int ntr, int par, double delta)
+inline void AlignmentTrack::modParam(trackPar_t* trSet, int ntr, int par, double delta)
 {
   // modify track parameter (VECTORIZE THOS)
   for (int i = 0; i < ntr; ++i) {

Detectors/Align/src/AlignmentTrack.cxx

@@ -168,7 +168,7 @@ bool AlignmentTrack::calcResidDeriv(double* extendedParams, bool invert, int pFr
   // (like http://root.cern.ch/root/html/ROOT__Math__RichardsonDerivator.html)
   //
   const auto& algConf = AlignConfig::Instance();
-  trackParam_t probD[kNRDClones]; // use this to vary supplied param for derivative calculation
+  trackPar_t probD[kNRDClones]; // use this to vary supplied param for derivative calculation
   double varDelta[kRichardsonN];
   const int kInvElem[kNKinParBON] = {-1, 1, 1, -1, -1};
   //
@@ -511,7 +511,7 @@ bool AlignmentTrack::calcResiduals(const double* extendedParams, bool invert, in
 }
 
 //______________________________________________________
-bool AlignmentTrack::propagateParamToPoint(trackParam_t* tr, int nTr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr)
+bool AlignmentTrack::propagateParamToPoint(trackPar_t* tr, int nTr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr)
 {
   // Propagate set of tracks to the point  (only parameters, no error matrix)
   // VECTORIZE this
@@ -521,7 +521,7 @@ bool AlignmentTrack::propagateParamToPoint(trackParam_t* tr, int nTr, const Alig
     if (!propagateParamToPoint(tr[itr], pnt, maxStep, maxSnp, mt, signCorr)) {
       if (algConf.verbose > 2) {
         LOG(error) << "Failed on clone " << itr << " propagation ";
-        tr[itr].print();
+        tr[itr].printParam();
         pnt->print(AlignmentPoint::kMeasurementBit | AlignmentPoint::kMaterialBit);
       }
       return false;
@@ -531,21 +531,33 @@ bool AlignmentTrack::propagateParamToPoint(trackParam_t* tr, int nTr, const Alig
 }
 
 //______________________________________________________
-bool AlignmentTrack::propagateParamToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr)
+bool AlignmentTrack::propagateParamToPoint(trackPar_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr)
 {
   // propagate tracks to the point (only parameters, no error matrix)
   return propagate(tr, pnt, maxStep, maxSnp, mt, nullptr, signCorr);
 }
 
 //______________________________________________________
-bool AlignmentTrack::propagateToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr)
+bool AlignmentTrack::propagateToPoint(trackParam_t& tr, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr)
 {
   // propagate tracks to the point. If matCor is true, then material corrections will be applied.
   // if matPar pointer is provided, it will be filled by total x2x0 and signed xrho
-  return propagate(tr, pnt, maxStep, maxSnp, mt, tLT, signCorr);
+  return propagate(tr, linRef, pnt, maxStep, maxSnp, mt, tLT, signCorr);
 }
 
-bool AlignmentTrack::propagate(trackParam_t& track, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr)
+bool AlignmentTrack::propagate(trackParam_t& track, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr)
+{
+  if (signCorr == 0) { // auto
+    // calculate the sign of the energy loss correction and ensure the upper leg of cosmics is calculated correctly.
+    double dx = pnt->getXTracking() - track.getX();
+    int dir = dx > 0.f ? 1 : -1;
+    signCorr = pnt->isInvDir() ? dir : -dir; // propagation along the track direction should have signCorr=-1
+  }
+  // do propagation in at least 2 step to reveal eventual effect of MS on the position
+  return PropagatorD::Instance()->propagateToAlphaX(track, linRef, pnt->getAlphaSens(), pnt->getXTracking(), pnt->getUseBzOnly(), maxSnp, maxStep, 2, mt, tLT, signCorr);
+}
+
+bool AlignmentTrack::propagate(trackPar_t& track, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr)
 {
   if (signCorr == 0) { // auto
     // calculate the sign of the energy loss correction and ensure the upper leg of cosmics is calculated correctly.
@@ -603,7 +615,7 @@ bool AlignmentTrack::ApplyMS(trackParam_t& trPar, double tms,double pms)
 */
 
 //______________________________________________________
-bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corrPar, const AlignmentPoint* pnt)
+bool AlignmentTrack::applyMatCorr(trackPar_t& trPar, const double* corrPar, const AlignmentPoint* pnt)
 {
   // Modify track param (e.g. trackParam_t) in the tracking frame
   // by delta accounting for material effects
@@ -630,7 +642,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corrPar, co
 }
 
 //______________________________________________________
-bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corr)
+bool AlignmentTrack::applyMatCorr(trackPar_t& trPar, const double* corr)
 {
   // Modify track param (e.g. trackParam_t) in the tracking frame
   // by delta accounting for material effects
@@ -645,7 +657,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corr)
         printf("%+.3e ", corr[i]);
       }
       printf("\n");
-      trPar.print();
+      trPar.printParam();
     }
     return false;
   }
@@ -656,7 +668,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corr)
 }
 
 //______________________________________________________
-bool AlignmentTrack::applyMatCorr(trackParam_t* trSet, int ntr, const double* corrDiag, const AlignmentPoint* pnt)
+bool AlignmentTrack::applyMatCorr(trackPar_t* trSet, int ntr, const double* corrDiag, const AlignmentPoint* pnt)
 {
   // Modify set of track params (e.g. trackParam_t) in the tracking frame
   // by delta accounting for material effects
@@ -683,7 +695,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t* trSet, int ntr, const double* co
     if (!applyMatCorr(trSet[itr], corr)) {
       if (algConf.verbose > 2) {
         LOGP(error, "Failed on clone {} materials", itr);
-        trSet[itr].print();
+        trSet[itr].printParam();
       }
       return false;
     }
@@ -732,7 +744,7 @@ double AlignmentTrack::richardsonExtrap(const double* val, int ord)
 }
 
 //______________________________________________
-void AlignmentTrack::richardsonDeriv(const trackParam_t* trSet, const double* delta, const AlignmentPoint* pnt, double& derY, double& derZ)
+void AlignmentTrack::richardsonDeriv(const trackPar_t* trSet, const double* delta, const AlignmentPoint* pnt, double& derY, double& derZ)
 {
   // Calculate Richardson derivatives for diagonalized Y and Z from a set of kRichardsonN pairs
   // of tracks with same parameter of i-th pair varied by +-delta[i]
@@ -882,7 +894,7 @@ bool AlignmentTrack::iniFit()
     //
     // propagate to reference point, which is the inner point of lower leg
     const AlignmentPoint* refP = getPoint(getInnerPointID());
-    if (!propagateToPoint(trcU, refP, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, -1)) { // moving along the track: energy is lost
+    if (!propagateToPoint(trcU, nullptr, refP, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, -1)) { // moving along the track: energy is lost
       return false;
     }
     //
@@ -1024,6 +1036,7 @@ bool AlignmentTrack::fitLeg(trackParam_t& trc, int pFrom, int pTo, bool& inv)
     }
     return false;
   }
+  trackPar_t linRef(trc), *linRefP = algConf.useLinRef ? &linRef : nullptr;
   trc.setCov(kIniErr);
   trc.setCov(16 * trc.getQ2Pt() * trc.getQ2Pt(), 4, 4); // lowest diagonal element (Q2Pt2)
   //
@@ -1042,7 +1055,7 @@ bool AlignmentTrack::fitLeg(trackParam_t& trc, int pFrom, int pTo, bool& inv)
   int pntCnt = 0;
   for (int ip = pFrom; ip != pTo; ip += pinc) { // inward fit from outer point
     AlignmentPoint* pnt = getPoint(ip);
-    if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // against track direction : e.loss is compensated
+    if (!propagateToPoint(trc, linRefP, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // against track direction : e.loss is compensated
       if (algConf.verbose > 2) {
         LOGF(warn, "Failed on propagateToPoint %d (%d : %d) %f", ip, pFrom, pTo, pnt->getXTracking());
         trc.print();
@@ -1139,7 +1152,7 @@ bool AlignmentTrack::residKalman()
       trc.invert();
       inv = !inv;
     }
-    if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) {
+    if (!propagateToPoint(trc, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) {
       return false;
     }
     if (!pnt->containsMeasurement()) {
@@ -1178,7 +1191,7 @@ bool AlignmentTrack::residKalman()
       trc.invert();
       inv = !inv;
     }
-    if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // we are going along track direction, e.loss is applied
+    if (!propagateToPoint(trc, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // we are going along track direction, e.loss is applied
       return false;
     }
     if (!pnt->containsMeasurement()) {
@@ -1335,7 +1348,7 @@ bool AlignmentTrack::processMaterials(trackParam_t& trc, int pFrom, int pTo)
     //
     matTL.clearFast();
     //    printf("-> ProcMat %d (%d->%d)\n",ip,pFrom,pTo);
-    if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), &matTL, signELoss)) { // with material corrections
+    if (!propagateToPoint(trc, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), &matTL, signELoss)) { // with material corrections
       if (algConf.verbose > 2) {
         LOG(error) << "Failed to take track to point" << ip << " (dir: " << pFrom << "->" << pTo << ") with mat.corr.";
         trc.print();
@@ -1346,7 +1359,7 @@ bool AlignmentTrack::processMaterials(trackParam_t& trc, int pFrom, int pTo)
     //
     // is there enough material to consider the point as a scatterer?
     bool hasMaterial = matTL.getX2X0() > minX2X0;
-    if (!propagateToPoint(tr0, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType::USEMatCorrNONE, nullptr, signELoss)) { // no material corrections
+    if (!propagateToPoint(tr0, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType::USEMatCorrNONE, nullptr, signELoss)) { // no material corrections
       if (algConf.verbose > 2) {
         LOG(error) << "Failed to take track to point" << ip << " (dir: " << pFrom << "->" << pTo << ") with mat.corr.";
         tr0.print();

Detectors/Base/include/DetectorsBase/Propagator.h

@@ -127,6 +127,9 @@ class PropagatorImpl
   GPUd() bool propagateToAlphaX(track_T& track, value_type alpha, value_type x, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, int minSteps = 1,
                                 MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const;
 
+  GPUd() bool propagateToAlphaX(TrackParCov_t& track, TrackPar_t* linRef, value_type alpha, value_type x, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, int minSteps = 1,
+                                MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const;
+
   template <typename track_T>
   GPUd() bool propagateToR(track_T& track, value_type r, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP,
                            MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const;