diff --git a/algo/ca/core/data/CaTrackParam.h b/algo/ca/core/data/CaTrackParam.h index 25fc57b716971ab499cd03cc7507e8479a4a60b2..8dbc71b78d008f4e08ec3cc1c555a5952582afb4 100644 --- a/algo/ca/core/data/CaTrackParam.h +++ b/algo/ca/core/data/CaTrackParam.h @@ -139,26 +139,53 @@ namespace cbm::algo::ca /// \return covariance matrix element /// DataT C00() const { return C<0, 0>(); } + DataT C01() const { return C<0, 1>(); } + DataT C02() const { return C<0, 2>(); } + DataT C03() const { return C<0, 3>(); } + DataT C04() const { return C<0, 4>(); } + DataT C05() const { return C<0, 5>(); } + DataT C06() const { return C<0, 6>(); } + DataT C10() const { return C<1, 0>(); } DataT C11() const { return C<1, 1>(); } + DataT C12() const { return C<1, 2>(); } + DataT C13() const { return C<1, 3>(); } + DataT C14() const { return C<1, 4>(); } + DataT C15() const { return C<1, 5>(); } + DataT C16() const { return C<1, 6>(); } + DataT C20() const { return C<2, 0>(); } DataT C21() const { return C<2, 1>(); } DataT C22() const { return C<2, 2>(); } + DataT C23() const { return C<2, 3>(); } + DataT C24() const { return C<2, 4>(); } + DataT C25() const { return C<2, 5>(); } + DataT C26() const { return C<2, 6>(); } + DataT C30() const { return C<3, 0>(); } DataT C31() const { return C<3, 1>(); } DataT C32() const { return C<3, 2>(); } DataT C33() const { return C<3, 3>(); } + DataT C34() const { return C<3, 4>(); } + DataT C35() const { return C<3, 5>(); } + DataT C36() const { return C<3, 6>(); } + DataT C40() const { return C<4, 0>(); } DataT C41() const { return C<4, 1>(); } DataT C42() const { return C<4, 2>(); } DataT C43() const { return C<4, 3>(); } DataT C44() const { return C<4, 4>(); } + DataT C45() const { return C<4, 5>(); } + DataT C46() const { return C<4, 6>(); } + DataT C50() const { return C<5, 0>(); } DataT C51() const { return C<5, 1>(); } DataT C52() const { return C<5, 2>(); } DataT C53() const { return C<5, 3>(); } DataT C54() const { return C<5, 4>(); } DataT C55() const { return C<5, 5>(); } + DataT C56() const { return C<5, 6>(); } + DataT C60() const { return C<6, 0>(); } DataT C61() const { return C<6, 1>(); } DataT C62() const { return C<6, 2>(); } @@ -437,26 +464,53 @@ namespace cbm::algo::ca /// \brief Individual references to covariance matrix elements /// DataT& C00() { return C<0, 0>(); } + DataT& C01() { return C<0, 1>(); } + DataT& C02() { return C<0, 2>(); } + DataT& C03() { return C<0, 3>(); } + DataT& C04() { return C<0, 4>(); } + DataT& C05() { return C<0, 5>(); } + DataT& C06() { return C<0, 6>(); } + DataT& C10() { return C<1, 0>(); } DataT& C11() { return C<1, 1>(); } + DataT& C12() { return C<1, 2>(); } + DataT& C13() { return C<1, 3>(); } + DataT& C14() { return C<1, 4>(); } + DataT& C15() { return C<1, 5>(); } + DataT& C16() { return C<1, 6>(); } + DataT& C20() { return C<2, 0>(); } DataT& C21() { return C<2, 1>(); } DataT& C22() { return C<2, 2>(); } + DataT& C23() { return C<2, 3>(); } + DataT& C24() { return C<2, 4>(); } + DataT& C25() { return C<2, 5>(); } + DataT& C26() { return C<2, 6>(); } + DataT& C30() { return C<3, 0>(); } DataT& C31() { return C<3, 1>(); } DataT& C32() { return C<3, 2>(); } DataT& C33() { return C<3, 3>(); } + DataT& C34() { return C<3, 4>(); } + DataT& C35() { return C<3, 5>(); } + DataT& C36() { return C<3, 6>(); } + DataT& C40() { return C<4, 0>(); } DataT& C41() { return C<4, 1>(); } DataT& C42() { return C<4, 2>(); } DataT& C43() { return C<4, 3>(); } DataT& C44() { return C<4, 4>(); } + DataT& C45() { return C<4, 5>(); } + DataT& C46() { return C<4, 6>(); } + DataT& C50() { return C<5, 0>(); } DataT& C51() { return C<5, 1>(); } DataT& C52() { return C<5, 2>(); } DataT& C53() { return C<5, 3>(); } DataT& C54() { return C<5, 4>(); } DataT& C55() { return C<5, 5>(); } + DataT& C56() { return C<5, 6>(); } + DataT& C60() { return C<6, 0>(); } DataT& C61() { return C<6, 1>(); } DataT& C62() { return C<6, 2>(); } diff --git a/algo/ca/core/pars/CaConstants.h b/algo/ca/core/pars/CaConstants.h index 38f0f7b1df820efc97f99a4602d28edefe746b3c..4fd8d9bf84eea25c78622441c03a4a2952522e8b 100644 --- a/algo/ca/core/pars/CaConstants.h +++ b/algo/ca/core/pars/CaConstants.h @@ -91,9 +91,11 @@ namespace cbm::algo::ca::constants /// Miscellaneous constants namespace misc { - constexpr int AssertionLevel = 0; ///< Assertion level - constexpr int Alignment = 16; ///< Default alignment of data (bytes) - } // namespace misc + constexpr int AssertionLevel = 0; ///< Assertion level + constexpr int Alignment = 16; ///< Default alignment of data (bytes) + constexpr fscal NegligibleFieldkG = 1.e-4; ///< Negligible field [kG] + + } // namespace misc /// \brief Undefined values template<typename T1, typename T2 = T1> diff --git a/algo/ca/core/pars/CaField.cxx b/algo/ca/core/pars/CaField.cxx index de8f60531808f477d064ab056f01da358c1bc4e1..fc9dfe59e1433b80f423e364dc9683fd720adb2b 100644 --- a/algo/ca/core/pars/CaField.cxx +++ b/algo/ca/core/pars/CaField.cxx @@ -240,6 +240,8 @@ template<typename DataT> void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z, const FieldValue<DataT>& b2, const DataT b2z) { + fIsZeroField = (b0.IsZeroField() && b1.IsZeroField() && b2.IsZeroField()); + z0 = b0z; DataT dz1 = b1z - b0z; DataT dz2 = b2z - b0z; @@ -274,6 +276,8 @@ void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const template<typename DataT> void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z) { + fIsZeroField = (b0.IsZeroField() && b1.IsZeroField()); + z0 = b0z; DataT dzi = utils::simd::One<DataT>() / (b1z - b0z); cx0 = b0.x; diff --git a/algo/ca/core/pars/CaField.h b/algo/ca/core/pars/CaField.h index b00b4508da66fbc90942fd54433d84801b9f45a2..57290fe35fa4bbb213fadb10a85f6bfc18de702f 100644 --- a/algo/ca/core/pars/CaField.h +++ b/algo/ca/core/pars/CaField.h @@ -45,6 +45,19 @@ namespace cbm::algo::ca void Combine(const FieldValue<fvec>& B, const fmask w); + /// Is the field value zero? + bool IsZeroField() const + { + constexpr auto e = constants::misc::NegligibleFieldkG; + auto isZero = (x * x + y * y + z * z <= e * e); + if constexpr (std::is_same<DataT, fvec>::value) { + return isZero.isFull(); + } + else { + return isZero; + } + } + /// Consistency checker void CheckConsistency() const; @@ -205,6 +218,9 @@ namespace cbm::algo::ca /// \return the magnetic field value FieldValue<DataT> Get(const DataT& x, const DataT& y, const DataT& z) const; + /// Is the field region empty? + bool IsZeroField() const { return fIsZeroField; } + /// Interpolates the magnetic field by three nodal points and sets the result to this FieldRegion object /// The result is a quadratic interpolation of the field as a function of z /// \param b0 field value in the first nodal point @@ -298,7 +314,9 @@ namespace cbm::algo::ca DataT z0{0.f}; ///< z-coordinate of the field region central point - bool fUseOriginalField{false}; + bool fUseOriginalField{false}; ///< Use original field instead of the field approximation + + bool fIsZeroField{false}; ///< Is the field region empty? /// Serialization function friend class boost::serialization::access; @@ -314,6 +332,9 @@ namespace cbm::algo::ca ar& cz0; ar& cz1; ar& cz2; + ar& z0; + ar& fUseOriginalField; + ar& fIsZeroField; } } _fvecalignment; @@ -340,6 +361,10 @@ namespace cbm::algo::ca copy(cz2, Tb.cz2); copy(z0, Tb.z0); + + fUseOriginalField = Tb.fUseOriginalField; + fIsZeroField = Tb.fIsZeroField; + } // CopyBase /// \brief Explicit instantiation declarations for FieldRegion class with specific template types to ensure proper fgOdiginalField definition diff --git a/algo/ca/core/pars/CaParameters.cxx b/algo/ca/core/pars/CaParameters.cxx index 333b286fe4a5223675726adca499568fb2e8721f..43db55d2ae635805f51061c5edaeead5735d334b 100644 --- a/algo/ca/core/pars/CaParameters.cxx +++ b/algo/ca/core/pars/CaParameters.cxx @@ -263,6 +263,16 @@ std::string Parameters<DataT>::ToString(int verbosity, int indentLevel) const msg << indent << indentCh << indentCh << indentCh << char(120 + dim) << " = " << utils::simd::Cast<DataT, float>(fTargetPos[dim], 0) << " cm\n"; } + msg << indent << indentCh << clrs::CLb << "FIELD:\n" << clrs::CL; + msg << indent << indentCh << indentCh << "Target field:\n"; + msg << indent << indentCh << indentCh << indentCh + << "Bx = " << utils::simd::Cast<DataT, float>(fVertexFieldValue.x, 0) << " Kg\n"; + msg << indent << indentCh << indentCh << indentCh + << "By = " << utils::simd::Cast<DataT, float>(fVertexFieldValue.y, 0) << " Kg\n"; + msg << indent << indentCh << indentCh << indentCh + << "Bz = " << utils::simd::Cast<DataT, float>(fVertexFieldValue.z, 0) << " Kg\n"; + + msg << indent << indentCh << clrs::CLb << "NUMBER OF STATIONS:\n" << clrs::CL; msg << indent << indentCh << indentCh << "Number of stations (Geometry): "; for (int iDet = 0; iDet < constants::size::MaxNdetectors; ++iDet) { diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx index c997f2fa1a669d69f2e1d5a0efdf87e993c29a5e..9e82e6c46d8ba4d1f647603899a802765135e583 100644 --- a/algo/ca/core/tracking/CaTrackFinder.cxx +++ b/algo/ca/core/tracking/CaTrackFinder.cxx @@ -103,6 +103,7 @@ void TrackFinder::FindTracks() ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih; const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); + const ca::Station<fvec>& st = frAlgo.GetParameters().GetStation(h.Station()); fscal dx = h.X() - targX; @@ -201,6 +202,9 @@ void TrackFinder::FindTracks() // cut data into sub-timeslices and process them one by one + std::vector<int> statNwindows(frAlgo.GetNofThreads(), 0); + std::vector<int> statNhitsProcessed(frAlgo.GetNofThreads(), 0); + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { bool areUntouchedDataLeft = true; // is the whole TS processed ca::HitIndex_t sliceFirstHit[nDataStreams]; @@ -209,8 +213,9 @@ void TrackFinder::FindTracks() sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream); } - int statNwindows = 0; - int statNhitsProcessed = 0; + statNwindows[iThread] = 0; + statNhitsProcessed[iThread] = 0; + int statLastLogTimeChunk = -1; while (areUntouchedDataLeft) { @@ -224,7 +229,7 @@ void TrackFinder::FindTracks() // TODO: SG: skip empty regions and start the subslice with the earliest hit - statNwindows++; + statNwindows[iThread]++; int statNwindowHits = 0; for (int iStream = 0; iStream < nDataStreams; ++iStream) { @@ -259,7 +264,7 @@ void TrackFinder::FindTracks() } // else the hit has been alread processed in previous sub-slices } } - statNhitsProcessed += statNwindowHits; + statNhitsProcessed[iThread] += statNwindowHits; // print the LOG for every 10 ms of data processed { @@ -270,10 +275,10 @@ void TrackFinder::FindTracks() if (dataRead > 100.) { dataRead = 100.; } - LOG(info) << "CA tracker process sliding window N " << statNwindows << ": time " << tsStart / 1.e6 << " ms + " - << fWindowLength / 1.e3 << " us) with " << statNwindowHits << " hits. " + LOG(info) << "CA tracker process sliding window N " << statNwindows[iThread] << ": time " << tsStart / 1.e6 + << " ms + " << fWindowLength / 1.e3 << " us) with " << statNwindowHits << " hits. " << " Processing " << dataRead << " % of the TS time and " - << 100. * statNhitsProcessed / statNhitsTotal << " % of TS hits." + << 100. * statNhitsProcessed[iThread] / statNhitsTotal << " % of TS hits." << " Already reconstructed " << frAlgo.fRecoTracks.size() << " tracks "; } } @@ -354,6 +359,17 @@ void TrackFinder::FindTracks() auto timerEnd = std::chrono::high_resolution_clock::now(); frAlgo.fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count()); + + int statNhitsProcessedTotal = 0; + int statNwindowsTotal = 0; + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { + statNhitsProcessedTotal += statNhitsProcessed[iThread]; + statNwindowsTotal += statNwindows[iThread]; + } + + LOG(info) << "CA tracker: time slice finished. Reconstructed " << frAlgo.fRecoTracks.size() << " tracks with " + << frAlgo.fRecoHits.size() << " hits. Processed " << statNhitsProcessedTotal << " hits in " + << statNwindowsTotal << " time windows. Reco time " << frAlgo.fCaRecoTime / 1.e9 << " s"; } // --------------------------------------------------------------------------------------------------------------------- diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx index d0b8b74b5a76af8691fec5bcde2b7312de93a07d..f1a33a5c5b894a9f09d23dee044fe3439fc6db09 100644 --- a/algo/ca/core/tracking/CaTrackFit.cxx +++ b/algo/ca/core/tracking/CaTrackFit.cxx @@ -603,11 +603,16 @@ namespace cbm::algo::ca { // use Q/p linearisation at fQp0 - fvec sgn = iif(fTr.GetZ() < z_out, fvec(1.), fvec(-1.)); - while (!(fMaskF * abs(z_out - fTr.GetZ()) <= fvec(1.e-6)).isFull()) { - fvec zNew = fTr.GetZ() + sgn * fMaxExtraplationStep; // max. 50 cm step - zNew(sgn * (z_out - zNew) <= fvec(0.)) = z_out; - ExtrapolateStep(zNew, F); + if (F.IsZeroField()) { + ExtrapolateLineNoField(z_out); + } + else { + fvec sgn = iif(fTr.GetZ() < z_out, fvec(1.), fvec(-1.)); + while (!(fMaskF * abs(z_out - fTr.GetZ()) <= fvec(1.e-6)).isFull()) { + fvec zNew = fTr.GetZ() + sgn * fMaxExtraplationStep; // max. 50 cm step + zNew(sgn * (z_out - zNew) <= fvec(0.)) = z_out; + ExtrapolateStep(zNew, F); + } } } @@ -832,62 +837,134 @@ namespace cbm::algo::ca fQp0 = qp0; } - void TrackFit::ExtrapolateLineNoField(fvec z_out) + void TrackFit::ExtrapolateLineNoField(fvec zOut) { // extrapolate the track assuming no field - // - // it is a copy of a sequence two routines - // L1ExtrapolateTime() and L1ExtrapolateLine() - // TODO: make it right - // - cnst c_light(29.9792458); - fvec dz = (z_out - fTr.GetZ()); - dz.setZero(!fMask); + // x += dz * tx + // y += dz * ty + // t += dz * sqrt ( 1 + tx*tx + ty*ty ) * vi - TrackParamV& T = fTr; + if (0) { // debug: full Runge-Kutta extrapolation with zero field values + FieldRegion<fvec> F; + ExtrapolateStep(zOut, F); + return; + } - // extrapolate time - fvec L = sqrt(T.Tx() * T.Tx() + T.Ty() * T.Ty() + fvec(1.)); - T.Time() += dz * L / c_light; + auto& t = fTr; // use reference to shorten the text - cnst k1 = dz * T.Tx() / (c_light * L); - cnst k2 = dz * T.Ty() / (c_light * L); + cnst zMasked = iif(fMask, zOut, t.GetZ()); - fvec c52 = T.C52(); - fvec c53 = T.C53(); + cnst dz = (zMasked - t.GetZ()); - T.C50() += k1 * T.C20() + k2 * T.C30(); - T.C51() += k1 * T.C21() + k2 * T.C31(); - T.C52() += k1 * T.C22() + k2 * T.C32(); - T.C53() += k1 * T.C32() + k2 * T.C33(); - T.C54() += k1 * T.C42() + k2 * T.C43(); - T.C55() += k1 * (T.C52() + c52) + k2 * (T.C53() + c53); + fvec tx = t.GetTx(); + fvec ty = t.GetTy(); + fvec vi = t.GetVi(); - // extrapolate trajectory + fvec L = sqrt(fvec(1.) + tx * tx + ty * ty); - T.X() += T.Tx() * dz; - T.Y() += T.Ty() * dz; - T.Z() += dz; + fvec j52 = dz * tx * vi / L; + fvec j53 = dz * ty * vi / L; + fvec j56 = dz * L; - cnst dzC32_in = dz * T.C32(); + // transport parameters - T.C21() += dzC32_in; - T.C10() += dz * (T.C21() + T.C30()); + t.X() += tx * dz; + t.Y() += ty * dz; + t.Time() += L * vi * dz; + t.Z() = zMasked; - cnst C20_in = T.C20(); + // transport covariance matrix - T.C20() += dz * T.C22(); - T.C00() += dz * (T.C20() + C20_in); + // ( 1 0 dz 0 0 0 0 ) + // ( 0 1 0 dz 0 0 0 ) + // ( 0 0 1 0 0 0 0 ) + // J = ( 0 0 0 1 0 0 0 ) + // ( 0 0 0 0 1 0 0 ) + // ( 0 0 j52 j53 0 1 j56 ) + // ( 0 0 0 0 0 0 1 ) - cnst C31_in = T.C31(); - T.C31() += dz * T.C33(); - T.C11() += dz * (T.C31() + C31_in); - T.C30() += dzC32_in; + // JC = J * C - T.C40() += dz * T.C42(); - T.C41() += dz * T.C43(); + fvec jc00 = t.C00() + dz * t.C20(); + //fvec jc01 = t.C01() + dz * t.C21(); + fvec jc02 = t.C02() + dz * t.C22(); + //fvec jc03 = t.C03() + dz * t.C23(); + //fvec jc04 = t.C04() + dz * t.C24(); + //fvec jc05 = t.C05() + dz * t.C25(); + //fvec jc06 = t.C06() + dz * t.C26(); + + + fvec jc10 = t.C10() + dz * t.C30(); + fvec jc11 = t.C11() + dz * t.C31(); + fvec jc12 = t.C12() + dz * t.C32(); + fvec jc13 = t.C13() + dz * t.C33(); + //fvec jc14 = t.C14() + dz * t.C34(); + //fvec jc15 = t.C15() + dz * t.C35(); + //fvec jc16 = t.C16() + dz * t.C36(); + + // jc2? = t.C2? + // jc3? = t.C3? + // jc4? = t.C4? + + fvec jc50 = t.C50() + j52 * t.C20() + j53 * t.C30() + j56 * t.C60(); + fvec jc51 = t.C51() + j52 * t.C21() + j53 * t.C31() + j56 * t.C61(); + fvec jc52 = t.C52() + j52 * t.C22() + j53 * t.C32() + j56 * t.C62(); + fvec jc53 = t.C53() + j52 * t.C23() + j53 * t.C33() + j56 * t.C63(); + fvec jc54 = t.C54() + j52 * t.C24() + j53 * t.C34() + j56 * t.C64(); + fvec jc55 = t.C55() + j52 * t.C25() + j53 * t.C35() + j56 * t.C65(); + fvec jc56 = t.C56() + j52 * t.C26() + j53 * t.C36() + j56 * t.C66(); + + // jc6? = t.C6? + + // transpose J + // + // ( 1 0 0 0 0 0 0 ) + // ( 0 1 0 0 0 0 0 ) + // ( dz 0 1 0 0 j52 0 ) + // J' = ( 0 dz 0 1 0 j53 0 ) + // ( 0 0 0 0 1 0 0 ) + // ( 0 0 0 0 0 1 0 ) + // ( 0 0 0 0 0 j56 1 ) + + + // C = JC * J' + + t.C00() = jc00 + jc02 * dz; + t.C10() = jc10 + jc12 * dz; + t.C20() = t.C20() + t.C22() * dz; + t.C30() = t.C30() + t.C32() * dz; + t.C40() = t.C40() + t.C42() * dz; + t.C50() = jc50 + jc52 * dz; + t.C60() = t.C60() + t.C62() * dz; + + t.C11() = jc11 + jc13 * dz; + t.C21() = t.C21() + t.C23() * dz; + t.C31() = t.C31() + t.C33() * dz; + t.C41() = t.C41() + t.C43() * dz; + t.C51() = jc51 + jc53 * dz; + t.C61() = t.C61() + t.C63() * dz; + + // t.C22 = jc22 == t.C22 -> unchanged + // t.C32 = jc32 == t.C32 -> unchanged + // t.C42 = jc42 == t.C42 -> unchanged + t.C52() = jc52; + // t.C62 = jc62 == t.C62 -> unchanged + + // t.C33 = jc33 == t.C33 -> unchanged + // t.C43 = jc43 == t.C43 -> unchanged + t.C53() = jc53; + // t.C63 = jc63 == t.C63 -> unchanged + + // t.C44 = jc44 == t.C44 -> unchanged + t.C54() = jc54; + // t.C64 = jc64 == t.C64 -> unchanged + + t.C55() = jc55 + jc52 * j52 + jc53 * j53 + jc56 * j56; + t.C65() = t.C65() + t.C62() * j52 + t.C63() * j53 + t.C66() * j56; + + // t.C66 = jc66 = t.C66 -> unchanged }