Ca: track fit: bugfix in masked SIMD operation

algo/ca/core/tracking/CaTrackFit.cxx

@@ -40,13 +40,16 @@ namespace cbm::algo::ca
     // correction to HCH is needed for the case when sigma2 is so small
     // with respect to HCH that it disappears due to the roundoff error
     //
-    DataT wi     = doProtect ? (fMaskF / (m.Du2() + DataT(1.0000001) * HCH)) : (fMaskF / (m.Du2() + HCH));
-    DataT zetawi = fMaskF * zeta / (kf::utils::iif(maskDoFilter, m.Du2(), DataT(0.)) + HCH);
+    DataT w  = m.Du2() + (doProtect ? (DataT(1.0000001) * HCH) : HCH);
+    DataT wi = kf::utils::iif(fMask, DataT(1.) / w, DataT(0.));
+
+    DataT zetawi = zeta / (kf::utils::iif(maskDoFilter, m.Du2(), DataT(0.)) + HCH);
+    zetawi       = kf::utils::iif(fMask, zetawi, DataT(0.));
 
     wi = kf::utils::iif(m.Du2() > DataT(0.), wi, DataT(0.));
 
     fTr.ChiSq() += m.Ndf() * zeta * zeta * wi;
-    fTr.Ndf() += m.Ndf() * fMaskF;
+    fTr.Ndf() += kf::utils::iif(fMask, m.Ndf(), DataT(0.));
 
     K1 = F1 * wi;
     K2 = F2 * wi;
@@ -115,21 +118,21 @@ namespace cbm::algo::ca
 
     DataT HCH = fTr.C55();
 
-    DataT w = kf::utils::iif(timeInfo, fMaskF, DataT(0.));
+    DataTmask mask = fMask && timeInfo;
 
     // when dt0 is much smaller than current time error,
     // set track time exactly to the measurement value without filtering
     // it helps to keep the initial time errors reasonably small
     // the calculations in the covariance matrix are not affected
 
-    const DataTmask maskDoFilter = (HCH < dt2 * 16.f);
+    const DataTmask maskDoFilter = mask && (HCH < dt2 * 16.f);
 
-    DataT wi     = w / (dt2 + DataT(1.0000001) * HCH);
-    DataT zeta   = fTr.Time() - t;
-    DataT zetawi = w * zeta / (kf::utils::iif(maskDoFilter, dt2, DataT(0.)) + HCH);
+    DataT wi     = kf::utils::iif(mask, DataT(1.) / (dt2 + DataT(1.0000001) * HCH), DataT(0.));
+    DataT zeta   = kf::utils::iif(mask, fTr.Time() - t, DataT(0.));
+    DataT zetawi = kf::utils::iif(mask, zeta / (kf::utils::iif(maskDoFilter, dt2, DataT(0.)) + HCH), DataT(0.));
 
     fTr.ChiSqTime() += kf::utils::iif(maskDoFilter, zeta * zeta * wi, DataT(0.));
-    fTr.NdfTime() += w;
+    fTr.NdfTime() += kf::utils::iif(mask, DataT(1.), DataT(0.));
 
     DataT K1 = F1 * wi;
     DataT K2 = F2 * wi;
@@ -461,10 +464,10 @@ namespace cbm::algo::ca
 
     HCH = F4 * h - F6;
 
-    DataT wi     = fMaskF / HCH;
-    DataT zetawi = fMaskF * zeta / HCH;
-    fTr.ChiSqTime() += zeta * zeta * wi;
-    fTr.NdfTime() += fMaskF;
+    DataT wi     = kf::utils::iif(fMask, DataT(1.) / HCH, DataT(0.));
+    DataT zetawi = kf::utils::iif(fMask, zeta / HCH, DataT(0.));
+    fTr.ChiSqTime() += kf::utils::iif(fMask, zeta * zeta * wi, DataT(0.));
+    fTr.NdfTime() += kf::utils::iif(fMask, DataT(1.), DataT(0.));
 
     K1 = F1 * wi;
     K2 = F2 * wi;
@@ -544,10 +547,10 @@ namespace cbm::algo::ca
 
     HCH = F6;
 
-    DataT wi     = fMaskF / HCH;
-    DataT zetawi = fMaskF * zeta / HCH;
-    fTr.ChiSqTime() += zeta * zeta * wi;
-    fTr.NdfTime() += fMaskF;
+    DataT wi     = kf::utils::iif(fMask, DataT(1.) / HCH, DataT(0.));
+    DataT zetawi = kf::utils::iif(fMask, zeta / HCH, DataT(0.));
+    fTr.ChiSqTime() += kf::utils::iif(fMask, zeta * zeta * wi, DataT(0.));
+    fTr.NdfTime() += kf::utils::iif(fMask, DataT(1.), DataT(0.));
 
     K1 = F1 * wi;
     K2 = F2 * wi;
@@ -621,7 +624,8 @@ namespace cbm::algo::ca
     }
     else {
       DataT sgn = kf::utils::iif(fTr.GetZ() < z_out, DataT(1.), DataT(-1.));
-      while (!kf::utils::isFull(fMaskF * kf::utils::fabs(z_out - fTr.GetZ()) <= DataT(1.e-6))) {
+      while (
+        !kf::utils::isFull(kf::utils::iif(fMask, kf::utils::fabs(z_out - fTr.GetZ()), DataT(0.)) <= DataT(1.e-6))) {
         DataT zNew = fTr.GetZ() + sgn * fMaxExtraplationStep;  // max. 50 cm step
         zNew       = kf::utils::iif(sgn * (z_out - zNew) <= DataT(0.), z_out, zNew);
         ExtrapolateStep(zNew, F);
@@ -1007,9 +1011,9 @@ namespace cbm::algo::ca
     //DataT a = ( (kONE+mass2*qp0*qp0t)*radThick*s0*s0 );
     DataT a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
 
-    fTr.C22() += fMaskF * txtx1 * a;
-    fTr.C32() += fMaskF * tx * ty * a;
-    fTr.C33() += fMaskF * (kONE + tyty);
+    fTr.C22() += kf::utils::iif(fMask, txtx1 * a, DataT(0.));
+    fTr.C32() += kf::utils::iif(fMask, tx * ty * a, DataT(0.));
+    fTr.C33() += kf::utils::iif(fMask, (kONE + tyty) * a, DataT(0.));
   }
 
   template<typename DataT>
@@ -1038,18 +1042,18 @@ namespace cbm::algo::ca
     DataT T23 = (thickness * thickness) / DataT(3.0);
     DataT T2  = thickness / DataT(2.0);
 
-    fTr.C00() += fMaskF * txtx1 * a * T23;
-    fTr.C10() += fMaskF * tx * ty * a * T23;
-    fTr.C20() += fMaskF * txtx1 * a * D * T2;
-    fTr.C30() += fMaskF * tx * ty * a * D * T2;
+    fTr.C00() += kf::utils::iif(fMask, txtx1 * a * T23, DataT(0.));
+    fTr.C10() += kf::utils::iif(fMask, tx * ty * a * T23, DataT(0.));
+    fTr.C20() += kf::utils::iif(fMask, txtx1 * a * D * T2, DataT(0.));
+    fTr.C30() += kf::utils::iif(fMask, tx * ty * a * D * T2, DataT(0.));
 
-    fTr.C11() += fMaskF * (kONE + tyty) * a * T23;
-    fTr.C21() += fMaskF * tx * ty * a * D * T2;
-    fTr.C31() += fMaskF * (kONE + tyty) * a * D * T2;
+    fTr.C11() += kf::utils::iif(fMask, (kONE + tyty) * a * T23, DataT(0.));
+    fTr.C21() += kf::utils::iif(fMask, tx * ty * a * D * T2, DataT(0.));
+    fTr.C31() += kf::utils::iif(fMask, (kONE + tyty) * a * D * T2, DataT(0.));
 
-    fTr.C22() += fMaskF * txtx1 * a;
-    fTr.C32() += fMaskF * tx * ty * a;
-    fTr.C33() += fMaskF * (kONE + tyty) * a;
+    fTr.C22() += kf::utils::iif(fMask, txtx1 * a, DataT(0.));
+    fTr.C32() += kf::utils::iif(fMask, tx * ty * a, DataT(0.));
+    fTr.C33() += kf::utils::iif(fMask, (kONE + tyty) * a, DataT(0.));
   }
 
 

algo/ca/core/tracking/CaTrackFit.h

@@ -54,16 +54,7 @@ namespace cbm::algo::ca
       SetTrack(t);
     }
 
-    void SetMask(const DataTmask& m)
-    {
-      fMask = m;
-      if constexpr (std::is_same<DataT, fvec>::value) {
-        fMaskF = iif(m, fvec::One(), fvec::Zero());
-      }
-      else {
-        fMaskF = (m ? 1. : 0.);
-      }
-    }
+    void SetMask(const DataTmask& m) { fMask = m; }
 
     template<typename T>
     void SetTrack(const TrackParamBase<T>& t)
@@ -253,7 +244,6 @@ namespace cbm::algo::ca
     /// Data members
 
     DataTmask fMask{true};  ///< mask of active elements in simd vectors
-    DataT fMaskF{1.};       ///< floating point representation of fMask
 
     TrackParamBase<DataT> fTr{};  ///< track parameters
     DataT fQp0{0.};