diff --git a/core/base/CbmTrackingDetectorInterfaceBase.cxx b/core/base/CbmTrackingDetectorInterfaceBase.cxx
index a30cb8aa62c6735d0c692a921aba287c46dab67c..8a46391ae24ce7376a8f65604a704f8dd1e5ceb0 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.cxx
+++ b/core/base/CbmTrackingDetectorInterfaceBase.cxx
@@ -67,20 +67,6 @@ bool CbmTrackingDetectorInterfaceBase::Check() const
msg << prefix << " negative or NaN inner radius (" << rMin << " cm)\n";
res = false && res;
}
-
- // Front strips stereo angle
- auto angleF = this->GetStripsStereoAngleFront(iSt);
- if (std::isnan(angleF)) {
- msg << prefix << " NaN front strips stereo angle (" << angleF << " rad)\n";
- res = false && res;
- }
-
- // Back strips stereo angle
- auto angleB = this->GetStripsStereoAngleBack(iSt);
- if (std::isnan(angleF)) {
- msg << prefix << " NaN back strips stereo angle (" << angleB << " rad)\n";
- res = false && res;
- }
}
// Position along beam axis
@@ -132,40 +118,8 @@ std::string CbmTrackingDetectorInterfaceBase::ToString() const
table << setw(10) << setfill(' ') << GetRmax(iSt) << ' ';
table << setw(10) << setfill(' ') << GetXmax(iSt) << ' ';
table << setw(10) << setfill(' ') << GetYmax(iSt) << ' ';
- table << setw(11) << setfill(' ') << GetStripsStereoAngleFront(iSt) << ' ';
- table << setw(11) << setfill(' ') << GetStripsStereoAngleBack(iSt) << ' ';
table << setw(10) << setfill(' ') << GetThickness(iSt) << ' ';
table << setw(10) << setfill(' ') << GetRadLength(iSt) << '\n';
}
return table.str();
}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void CbmTrackingDetectorInterfaceBase::InitConversionMatrices()
-{
- int nStations = GetNtrackingStations();
- fvConvUVtoXY.clear();
- fvConvUVtoXY.resize(nStations);
- fvConvXYtoUV.clear();
- fvConvXYtoUV.resize(nStations);
- for (int iS = 0; iS < nStations; ++iS) {
- double phiF = GetStripsStereoAngleFront(iS);
- double phiB = GetStripsStereoAngleBack(iS);
- double cF = std::cos(phiF);
- double sF = std::sin(phiF);
- double cB = std::cos(phiB);
- double sB = std::sin(phiB);
- double det = cF * sB - sF * cB;
- // (x,y) -> (u,v)
- fvConvUVtoXY[iS][0] = cF;
- fvConvUVtoXY[iS][1] = sF;
- fvConvUVtoXY[iS][2] = cB;
- fvConvUVtoXY[iS][3] = sB;
- // (u,v) -> (x,y)
- fvConvXYtoUV[iS][0] = sB / det;
- fvConvXYtoUV[iS][1] = -sF / det;
- fvConvXYtoUV[iS][2] = -cB / det;
- fvConvXYtoUV[iS][3] = cF / det;
- }
-}
diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h
index 2f53e4a81cb1e750ccf3f846eaa6ea4bddd438cf..c0928d96a1b45df8f36980c7d5034ba281530072 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.h
+++ b/core/base/CbmTrackingDetectorInterfaceBase.h
@@ -54,16 +54,6 @@ public:
/// \return Size of station outer radius [cm]
virtual double GetRmax(int stationId) const = 0;
- /// Gets front strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for front strips [rad]
- virtual double GetStripsStereoAngleFront(int stationId) const = 0;
-
- /// Gets back strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for back strips [rad]
- virtual double GetStripsStereoAngleBack(int stationId) const = 0;
-
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
@@ -105,71 +95,22 @@ public:
/// \return Flag: true - station provides time measurements, false - station does not provide time measurements
virtual bool IsTimeInfoProvided(int stationId) const = 0;
+ /// Gets stereo angles of the two independent measured coordinates
+ /// The tracking does not use this method. It is only used by the QA task.
+ /// \param address detector unique identifier
+ /// \return [phiU, phiV] - Stereo angles [rad]
+ virtual std::tuple<double, double> GetStereoAnglesSensor(int address) const = 0;
+
/// Technical flag: true - all downcasts are done with dynamic_cast followed by checks for nullptr (increases
// computing time, better for debug), false - all downcasts are done with static_cast without sanity checks
// (decreases computing time, but can cause bugs)
static constexpr bool kUseDynamicCast {true};
- /// @brief Conversion function (x,y) -> (u,v)
- /// @param iSt Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// @param x X-coordinate
- /// @param y Y-coordinate
- /// @return pair [u, v]
- std::pair<double, double> ConvertXYtoUV(int iSt, double x, double y) const
- {
- return std::make_pair(fvConvXYtoUV[iSt][0] * x + fvConvXYtoUV[iSt][1] * y,
- fvConvXYtoUV[iSt][2] * x + fvConvXYtoUV[iSt][3] * y);
- }
-
- /// @brief Conversion function (x,y) -> (u,v)
- /// @param iSt Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// @param u U-coordinate
- /// @param v V-coordinate
- /// @return pair [x, y]
- std::pair<double, double> ConvertUVtoXY(int iSt, double u, double v) const
- {
- return std::make_pair(fvConvUVtoXY[iSt][0] * u + fvConvUVtoXY[iSt][1] * v,
- fvConvUVtoXY[iSt][2] * u + fvConvUVtoXY[iSt][3] * v);
- }
-
- /// @brief Conversion function (du2, dv2) -> (dx2, dxy, dy2)
- /// @param iSt Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// @param du2 Variance of U-coordinate measurement
- /// @param dv2 Variance of V-coordinate measurement
- /// @return tuple [dx2, dxy, dy2]
- std::tuple<double, double, double> ConvertCovMatrixUVtoXY(int iSt, double du2, double dv2) const
- {
- return std::make_tuple(
- fvConvUVtoXY[iSt][0] * fvConvUVtoXY[iSt][0] * du2 + fvConvUVtoXY[iSt][1] * fvConvUVtoXY[iSt][1] * dv2,
- fvConvUVtoXY[iSt][0] * fvConvUVtoXY[iSt][2] * du2 + fvConvUVtoXY[iSt][1] * fvConvUVtoXY[iSt][3] * dv2,
- fvConvUVtoXY[iSt][2] * fvConvUVtoXY[iSt][2] * du2 + fvConvUVtoXY[iSt][3] * fvConvUVtoXY[iSt][3] * dv2);
- }
-
/// Checks detector interface: boundary conditions of the parameters
bool Check() const;
/// Prints all the parameters into table and saves the table as a string
std::string ToString() const;
-
-protected:
- /// @brief Initialized conversion matrices between (u,v) and (x,y) coordinate systems
- void InitConversionMatrices();
-
- /// @brief Vector of arrays of conversion matrices from (u,v) to (x,y) coordinate systems in different stations
- ///
- /// @note Indeces specification:
- /// |x| = |fvConvUVtoXY[0] fvConvUVtoXY[1]| * |u|
- /// |y| |fvConvUVtoXY[2] fvConvUVtoXY[2]| |v|
- ///
- std::vector<std::array<double, 4>> fvConvUVtoXY;
-
- /// @brief Vector of arrays of conversion matrices from (x,y) to (u,v) coordinate systems in different stations
- ///
- /// @note Indeces specification:
- /// |u| = |fvConvXYtoUV[0] fvConvXYtoUV[1]| * |x|
- /// |v| |fvConvXYtoUV[2] fvConvXYtoUV[2]| |y|
- ///
- std::vector<std::array<double, 4>> fvConvXYtoUV;
};
// **********************************************************
diff --git a/core/detectors/much/CbmMuchTrackingInterface.cxx b/core/detectors/much/CbmMuchTrackingInterface.cxx
index 3f58aadf976d672e5c7e628897bb25b9f2a8c9aa..465a4b898a1ea7193c03413441c50b454069eea1 100644
--- a/core/detectors/much/CbmMuchTrackingInterface.cxx
+++ b/core/detectors/much/CbmMuchTrackingInterface.cxx
@@ -77,9 +77,6 @@ InitStatus CbmMuchTrackingInterface::Init()
return kFATAL;
}
- // Init conversion matrices
- InitConversionMatrices();
-
// Keep Track of how many layers are present in each station (may be heterogenous)
for (int iMuchSt = 0; iMuchSt < fGeoScheme->GetNStations(); ++iMuchSt) {
fNbLayersPerStation.push_back(fGeoScheme->GetStation(iMuchSt)->GetNLayers());
diff --git a/core/detectors/much/CbmMuchTrackingInterface.h b/core/detectors/much/CbmMuchTrackingInterface.h
index ec82756563560582aa386bb0535e5063b561031e..bb84f4942a21777efed1fd501762fd8830b418a3 100644
--- a/core/detectors/much/CbmMuchTrackingInterface.h
+++ b/core/detectors/much/CbmMuchTrackingInterface.h
@@ -72,16 +72,6 @@ public:
/// \return Size of station inner radius [cm]
double GetRmin(int /*stationId*/) const { return 10.; }
- /// Gets back strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for back strips [rad]
- double GetStripsStereoAngleBack(int /*stationId*/) const { return TMath::Pi() / 2.; }
-
- /// Gets front strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for front strips [rad]
- double GetStripsStereoAngleFront(int /*stationId*/) const { return 0.; }
-
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
@@ -117,6 +107,12 @@ public:
/// \return Flag: true - station provides time measurements, false - station does not provide time measurements
bool IsTimeInfoProvided(int /*stationId*/) const { return true; }
+ /// Gets stereo angles of the two independent measured coordinates
+ /// The tracking does not use this method. It is only used by the QA task.
+ /// \param address detector unique identifier
+ /// \return [phiU, phiV] - Stereo angles [rad]
+ std::tuple<double, double> GetStereoAnglesSensor(int /*address*/) const { return std::tuple(0., TMath::Pi() / 2.); }
+
/// FairTask: sets parameter containers up
void SetParContainers();
diff --git a/core/detectors/mvd/CbmMvdTrackingInterface.cxx b/core/detectors/mvd/CbmMvdTrackingInterface.cxx
index af92e9908715eb66fc812d79f139fe0f98f6be03..a4db0e0483c77791d52ad0e6e31f9beaeb34b11a 100644
--- a/core/detectors/mvd/CbmMvdTrackingInterface.cxx
+++ b/core/detectors/mvd/CbmMvdTrackingInterface.cxx
@@ -41,9 +41,6 @@ InitStatus CbmMvdTrackingInterface::Init()
if (!fMvdStationPar) { return kFATAL; }
- // Init conversion matrices
- InitConversionMatrices();
-
// Check the validity of the parameters
if (!this->Check()) {
LOG(error)
diff --git a/core/detectors/mvd/CbmMvdTrackingInterface.h b/core/detectors/mvd/CbmMvdTrackingInterface.h
index 7cdefa87849197361f6f67c8c64fa4ec963cd9d1..e26f52c7987c7e8d9c4f9cde9d9f43b842e2e2e8 100644
--- a/core/detectors/mvd/CbmMvdTrackingInterface.h
+++ b/core/detectors/mvd/CbmMvdTrackingInterface.h
@@ -80,26 +80,6 @@ public:
return std::min(fMvdStationPar->GetBeamHeight(stationId), fMvdStationPar->GetBeamWidth(stationId));
}
- /// Gets spatial resolution (RMS) for back strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsBack(int stationId) const { return fMvdStationPar->GetYRes(stationId) / 10000.; }
-
- /// Gets spatial resolution (RMS) for front strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsFront(int stationId) const { return fMvdStationPar->GetXRes(stationId) / 10000.; }
-
- /// Gets back strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for back strips [rad]
- double GetStripsStereoAngleBack(int /*stationId*/) const { return TMath::Pi() / 2.; }
-
- /// Gets front strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for front strips [rad]
- double GetStripsStereoAngleFront(int /*stationId*/) const { return 0.; }
-
/// Gets the tracking station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
@@ -144,6 +124,12 @@ public:
/// \return Flag: true - station provides time measurements, false - station does not provide time measurements
bool IsTimeInfoProvided(int /*stationId*/) const { return false; }
+ /// Gets stereo angles of the two independent measured coordinates
+ /// The tracking does not use this method. It is only used by the QA task.
+ /// \param address detector unique identifier
+ /// \return [phiU, phiV] - Stereo angles [rad]
+ std::tuple<double, double> GetStereoAnglesSensor(int /*address*/) const { return std::tuple(0., TMath::Pi() / 2.); }
+
/// FairTask: sets parameter containers up
void SetParContainers();
diff --git a/core/detectors/sts/CbmStsTrackingInterface.cxx b/core/detectors/sts/CbmStsTrackingInterface.cxx
index 035a3a2a0e4e985c6f9fde870fe427624cb1404a..f80ee34dcba3e48e1c934361adbac837544e0bf7 100644
--- a/core/detectors/sts/CbmStsTrackingInterface.cxx
+++ b/core/detectors/sts/CbmStsTrackingInterface.cxx
@@ -37,23 +37,7 @@ CbmStsTrackingInterface::~CbmStsTrackingInterface()
//-------------------------------------------------------------------------------------------------------------------------------------
//
-double CbmStsTrackingInterface::GetStripsStereoAngleFront(int stationId) const
-{
- auto station = GetStsStation(stationId);
- return station->GetSensorRotation() + station->GetSensorStereoAngle(0) * TMath::Pi() / 180.;
-}
-
-//-------------------------------------------------------------------------------------------------------------------------------------
-//
-double CbmStsTrackingInterface::GetStripsStereoAngleBack(int stationId) const
-{
- auto station = GetStsStation(stationId);
- return station->GetSensorRotation() + station->GetSensorStereoAngle(1) * TMath::Pi() / 180.;
-}
-
-//-------------------------------------------------------------------------------------------------------------------------------------
-//
-std::tuple<double, double> CbmStsTrackingInterface::GetStripStereoAnglesSensor(int address) const
+std::tuple<double, double> CbmStsTrackingInterface::GetStereoAnglesSensor(int address) const
{
/// Gets front & back strip angle
@@ -107,9 +91,6 @@ InitStatus CbmStsTrackingInterface::Init()
if (!stsSetup->IsSensorParsInit()) { stsSetup->SetSensorParameters(fStsParSetSensor); }
if (!stsSetup->IsSensorCondInit()) { stsSetup->SetSensorConditions(fStsParSetSensorCond); }
- // Init conversion matrices
- InitConversionMatrices();
-
// Check the validity of the parameters
if (!this->Check()) {
LOG(error)
diff --git a/core/detectors/sts/CbmStsTrackingInterface.h b/core/detectors/sts/CbmStsTrackingInterface.h
index 26694a34a2c1cb23e6c188bf45cf2f7e38372f4f..a629218908e76a3f9935eb99bdc9b847c807b4c9 100644
--- a/core/detectors/sts/CbmStsTrackingInterface.h
+++ b/core/detectors/sts/CbmStsTrackingInterface.h
@@ -71,20 +71,11 @@ public:
/// \return Size of station inner radius [cm]
double GetRmin(int /*stationId*/) const { return 0.; }
- /// Gets back strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for back strips [rad]
- double GetStripsStereoAngleBack(int stationId) const;
-
- /// Gets front strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for front strips [rad]
- double GetStripsStereoAngleFront(int stationId) const;
-
- /// Gets front and back strip stereo angle
- /// \param address Unique element address
- /// \return std::tie(front,back) - Stereo angle for the front and back strips [rad]
- std::tuple<double, double> GetStripStereoAnglesSensor(int address) const;
+ /// Gets stereo angles of the two independent measured coordinates
+ /// The tracking does not use this method. It is only used by the QA task.
+ /// \param address detector unique identifier
+ /// \return [phiU, phiV] - Stereo angles [rad]
+ std::tuple<double, double> GetStereoAnglesSensor(int address) const;
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
diff --git a/core/detectors/tof/CbmTofTrackingInterface.cxx b/core/detectors/tof/CbmTofTrackingInterface.cxx
index 9cccf72f1ec931ce3009e5b1a62977881a4b46ce..cbb3c55691a6224638d0993f8ba2d435d81b30a1 100644
--- a/core/detectors/tof/CbmTofTrackingInterface.cxx
+++ b/core/detectors/tof/CbmTofTrackingInterface.cxx
@@ -83,9 +83,6 @@ InitStatus CbmTofTrackingInterface::Init()
fTofStationZ[iSt] = fTofStationZ[iSt] / nTofStationModules[iSt];
}
- // Init conversion matrices
- InitConversionMatrices();
-
// Check the validity of the parameters
if (!this->Check()) {
LOG(error)
diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h
index 6167bf824178155f43a89735e5a74e0a680fa586..923ea541892b3ad997e017d0e4f56a635b16b5eb 100644
--- a/core/detectors/tof/CbmTofTrackingInterface.h
+++ b/core/detectors/tof/CbmTofTrackingInterface.h
@@ -71,16 +71,6 @@ public:
/// \return Size of station inner radius [cm]
double GetRmin(int /*stationId*/) const { return 0.; }
- /// Gets back strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for back strips [rad]
- double GetStripsStereoAngleBack(int /*stationId*/) const { return 0.5 * TMath::Pi(); }
-
- /// Gets front strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for front strips [rad]
- double GetStripsStereoAngleFront(int /*stationId*/) const { return 0.; }
-
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
@@ -128,6 +118,12 @@ public:
/// \return Flag: true - station provides time measurements, false - station does not provide time measurements
bool IsTimeInfoProvided(int /*stationId*/) const { return true; }
+ /// Gets stereo angles of the two independent measured coordinates
+ /// The tracking does not use this method. It is only used by the QA task.
+ /// \param address detector unique identifier
+ /// \return [phiU, phiV] - Stereo angles [rad]
+ std::tuple<double, double> GetStereoAnglesSensor(int /*address*/) const { return std::tuple(0., TMath::Pi() / 2.); }
+
/// FairTask: sets parameter containers up
void SetParContainers();
diff --git a/core/detectors/trd/CbmTrdTrackingInterface.cxx b/core/detectors/trd/CbmTrdTrackingInterface.cxx
index 6e0ad828afe28865ba7455fd0edcbc786f3f9c40..00e6a8c0c9e9dc1027f9f3de05f1497b164c7fbe 100644
--- a/core/detectors/trd/CbmTrdTrackingInterface.cxx
+++ b/core/detectors/trd/CbmTrdTrackingInterface.cxx
@@ -68,9 +68,6 @@ InitStatus CbmTrdTrackingInterface::Init()
}
}
- // Init conversion matrices
- InitConversionMatrices();
-
// Check the validity of the parameters
if (!this->Check()) {
LOG(error)
diff --git a/core/detectors/trd/CbmTrdTrackingInterface.h b/core/detectors/trd/CbmTrdTrackingInterface.h
index fd732b9282149fe9b6fe6ba8d7bc07a054edb83a..be23a9472bba68a7ef614d900ae322b3103bc900 100644
--- a/core/detectors/trd/CbmTrdTrackingInterface.h
+++ b/core/detectors/trd/CbmTrdTrackingInterface.h
@@ -68,16 +68,6 @@ public:
/// \return Size of station inner radius [cm]
double GetRmin(int /*stationId*/) const { return 0.; }
- /// Gets back strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for back strips [rad]
- double GetStripsStereoAngleBack(int /*stationId*/) const { return 0.5 * TMath::Pi(); }
-
- /// Gets front strips stereo angle
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Absolute stereo angle for front strips [rad]
- double GetStripsStereoAngleFront(int /*stationId*/) const { return 0.; }
-
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
@@ -113,6 +103,12 @@ public:
/// \return Flag: true - station provides time measurements, false - station does not provide time measurements
bool IsTimeInfoProvided(int /*stationId*/) const { return true; }
+ /// Gets stereo angles of the two independent measured coordinates
+ /// The tracking does not use this method. It is only used by the QA task.
+ /// \param address detector unique identifier
+ /// \return [phiU, phiV] - Stereo angles [rad]
+ std::tuple<double, double> GetStereoAnglesSensor(int /*address*/) const { return std::tuple(0., TMath::Pi() / 2.); }
+
/// FairTask: sets parameter containers up
void SetParContainers();
diff --git a/macro/L1/configs/config_ideal_hits.yaml b/macro/L1/configs/config_ideal_hits.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..46b1533380a573eae73d9bfdf1cf0c5eb3c3aab4
--- /dev/null
+++ b/macro/L1/configs/config_ideal_hits.yaml
@@ -0,0 +1,87 @@
+# Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+# SPDX-License-Identifier: GPL-3.0-only
+# Authors: Sergei Zharko [committer]
+#
+## @file config_ideal_hits_mcbm.yaml
+## @brief Configuration file for ideal hit producers
+## @since 28.06.2023
+## @author Sergei Zharko <s.zharko@gsi.de>
+
+ideal_hit_producer:
+ #
+ # Flags are defined for each station
+ # 0 - real hits
+ # 1 - real hit position and time are adjusted with MC (exact quantity values are used)
+ # 2 - real hit position and time are adjusted with MC (true values are smeared within corresponding error)
+ # 3 - ideal hits are created (exact quantity values are used)
+ # 4 - ideal hits are created (true values are smeared withing errors, defined in section parameters)
+ flags:
+ sts:
+ 0: 4
+ 1: 4
+ 2: 4
+ 3: 4
+ 4: 4
+ 5: 4
+ 6: 4
+ 7: 4
+
+ #
+ # Notes:
+ # 1) For MVD, MuCh, TRD and TOF du corresponds to dx, and dv corresponds to dy
+ # 2) Spatial errors are expressed in cm, time errors are expressed in ns
+ # 3) pdf stands for probability density function
+ #
+ # Available PDFs:
+ # g: bounded gaussian distribution
+ # u: uniform distribution
+ parameters:
+ sts:
+ 0:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
+ 1:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
+ 2:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
+ 3:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
+ 4:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
+ 5:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
+ 6:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
+ 7:
+ {
+ du: { angle: 0.0, value: 15.e-4, pdf: g },
+ dv: { angle: 7.5, value: 15.e-4, pdf: g },
+ dt: { value: 2.7, pdf: g },
+ }
diff --git a/macro/L1/configs/config_ideal_hits_mcbm.yaml b/macro/L1/configs/config_ideal_hits_mcbm.yaml
index 81e896ff43643efda2377db57ed6518bd363ece6..d1abdd1f8a61c5587b59c12d9cc015b3cd8c5cca 100644
--- a/macro/L1/configs/config_ideal_hits_mcbm.yaml
+++ b/macro/L1/configs/config_ideal_hits_mcbm.yaml
@@ -34,9 +34,9 @@ ideal_hit_producer:
# u: uniform distribution
parameters:
sts:
- 0: { du: { value: 10.e-4, pdf: g }, dv: { value: 10.e-4, pdf: g }, dt: { value: 2.9, pdf: g }}
- 1: { du: { value: 10.e-4, pdf: g }, dv: { value: 10.e-4, pdf: g }, dt: { value: 2.9, pdf: g }}
+ 0: { du: { angle: 0, value: 10.e-4, pdf: g }, dv: { angle: 7.5, value: 10.e-4, pdf: g }, dt: { value: 2.9, pdf: g }}
+ 1: { du: { angle: 0, value: 10.e-4, pdf: g }, dv: { angle: 7.5, value: 10.e-4, pdf: g }, dt: { value: 2.9, pdf: g }}
trd:
- 0: { du: { value: .5, pdf: g }, dv: { value: .5, pdf: g }, dt: { value: 10., pdf: g }}
- 1: { du: { value: .05, pdf: g }, dv: { value: 4.0, pdf: u }, dt: { value: 10., pdf: g }}
- 2: { du: { value: 4.0, pdf: u }, dv: { value: .05, pdf: g }, dt: { value: 10., pdf: g }}
+ 0: { du: { angle: 0, value: .50, pdf: g }, dv: { angle: 90, value: .50, pdf: g }, dt: { value: 10., pdf: g }}
+ 1: { du: { angle: 0, value: .05, pdf: g }, dv: { angle: 90, value: 4.0, pdf: u }, dt: { value: 10., pdf: g }}
+ 2: { du: { angle: 0, value: 4.0, pdf: u }, dv: { angle: 90, value: .05, pdf: g }, dt: { value: 10., pdf: g }}
diff --git a/macro/L1/run_reco_L1global.C b/macro/L1/run_reco_L1global.C
index 22380619443e8ec3661f190596b646a017ed07b1..a1f6a0780812a2d85a681dfe7e647b8dd8825f9b 100644
--- a/macro/L1/run_reco_L1global.C
+++ b/macro/L1/run_reco_L1global.C
@@ -371,7 +371,6 @@ void run_reco_L1global(TString input = "", Int_t nTimeSlices = -1, Int_t firstTi
l1 = new CbmL1("L1", 0);
}
l1->SetGlobalMode();
- //l1->SetUseMcHit(0, 0, 0, 1, 0);
//l1->SetInputConfigName("");
run->AddTask(l1);
diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C
index 373050ec88cbc58c947d324516052816127cd67e..6806562132a51a9dc6b7947cefd549c7f4204dd7 100644
--- a/macro/run/run_reco.C
+++ b/macro/run/run_reco.C
@@ -387,6 +387,12 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice =
auto kalman = new CbmKF();
run->AddTask(kalman);
+ if (debugWithMC) {
+ auto* pIdealHitProducer = new cbm::ca::IdealHitProducer("IdealHitProducer", 3);
+ pIdealHitProducer->SetConfigName(srcDir + "/macro/L1/configs/config_ideal_hits.yaml");
+ //run->AddTask(pIdealHitProducer);
+ }
+
// L1 tracking
auto l1 = (debugWithMC) ? new CbmL1("CA", 2, 3) : new CbmL1("CA");
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index ad21d0b959463683dc5a6e702095c0f87bcf0f41..37f26dfaa2e0049ee4a86d5671a45cf02b4ebe4e 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -136,14 +136,14 @@ inline int CbmL1PFFitter::GetStsStationIndex(const CbmStsHit* hit)
}
-void FilterFirst(L1Fit& fit, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
+void FilterFirst(L1Fit& fit, fvec& x, fvec& y, fvec& t, L1XYMeasurementInfo& cov, fvec& dt2)
{
L1TrackPar& tr = fit.Tr();
- tr.ResetErrors(dxx, dyy, 1., 1., 1., 1.e6, 1.e2);
- tr.C10 = dxy;
+ tr.ResetErrors(cov.C00, cov.C11, 1., 1., 1., dt2, 1.e2);
+ tr.C10 = cov.C10;
tr.x = x;
tr.y = y;
- tr.t = 0.;
+ tr.t = t;
tr.vi = 0.;
tr.NDF = -3.0;
}
@@ -172,24 +172,21 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
int ista;
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
- fvec* x = new fvec[nHits];
- fvec* u = new fvec[nHits];
- fvec* v = new fvec[nHits];
- fvec* t = new fvec[nHits];
- fvec* du2 = new fvec[nHits];
- fvec* dv2 = new fvec[nHits];
- fvec* dt2 = new fvec[nHits];
- fvec* y = new fvec[nHits];
- fvec* z = new fvec[nHits];
- fmask* w = new fmask[nHits];
+ fvec x[nHits];
+ fvec y[nHits];
+ fvec z[nHits];
+ fvec t[nHits];
+ L1XYMeasurementInfo cov[nHits];
+ fvec dt2[nHits];
+ fmask w[nHits];
// fvec y_temp;
- fvec x_first, y_first, x_last, y_last;
- fvec fstC00, fstC10, fstC11;
- fvec lstC00, lstC10, lstC11;
+ fvec x_first, y_first, t_first, x_last, y_last, t_last;
+ L1XYMeasurementInfo cov_first, cov_last;
+ fvec dt2_first, dt2_last;
fvec z0, z1, z2, dz, z_start, z_end;
- L1FieldValue* fB = new L1FieldValue[nHits];
+ L1FieldValue fB[nHits];
L1FieldValue fB_temp _fvecalignment;
unsigned short N_vTracks = Tracks.size();
@@ -253,72 +250,59 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
int nHitsTrackSts = tr[iVec]->GetNofStsHits();
int nHitsTrack = nHitsTrackMvd + nHitsTrackSts;
for (i = 0; i < nHitsTrack; i++) {
- float posx = 0.f, posy = 0.f, posz = 0.f, time = 0.;
+ const CbmPixelHit* hit;
if (i < nHitsTrackMvd) {
- int hitIndex = tr[iVec]->GetMvdHitIndex(i);
- const CbmMvdHit* hit = &(vMvdHits[hitIndex]);
-
- posx = hit->GetX();
- posy = hit->GetY();
- posz = hit->GetZ();
- time = hit->GetTime();
- ista = GetMvdStationIndex(hit);
+ int hitIndex = tr[iVec]->GetMvdHitIndex(i);
+ const CbmMvdHit* mvdHit = &(vMvdHits[hitIndex]);
+ ista = GetMvdStationIndex(mvdHit);
if (ista < 0) continue;
-
- du2[ista][iVec] = hit->GetDx() * hit->GetDx();
- dv2[ista][iVec] = hit->GetDy() * hit->GetDy();
- dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
+ hit = mvdHit;
}
else {
- int hitIndex = tr[iVec]->GetHitIndex(i - nHitsTrackMvd);
- const CbmStsHit* hit = &(vStsHits[hitIndex]);
-
- posx = hit->GetX();
- posy = hit->GetY();
- posz = hit->GetZ();
- time = hit->GetTime();
- ista = GetStsStationIndex(hit);
+ int hitIndex = tr[iVec]->GetHitIndex(i - nHitsTrackMvd);
+ const CbmStsHit* stsHit = &(vStsHits[hitIndex]);
+ ista = GetStsStationIndex(stsHit);
if (ista < 0) continue;
-
- du2[ista][iVec] = hit->GetDu() * hit->GetDu();
- dv2[ista][iVec] = hit->GetDv() * hit->GetDv();
- dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
+ hit = stsHit;
}
+
+ x[ista][iVec] = hit->GetX();
+ y[ista][iVec] = hit->GetY();
+ z[ista][iVec] = hit->GetZ();
+ t[ista][iVec] = hit->GetTime();
+
+ cov[ista].C00[iVec] = hit->GetDx() * hit->GetDx();
+ cov[ista].C10[iVec] = hit->GetDxy();
+ cov[ista].C11[iVec] = hit->GetDy() * hit->GetDy();
+ dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
+
w[ista][iVec] = true;
- x[ista][iVec] = posx;
- y[ista][iVec] = posy;
- u[ista][iVec] = posx * sta[ista].frontInfo.cos_phi[0] + posy * sta[ista].frontInfo.sin_phi[0];
- v[ista][iVec] = posx * sta[ista].backInfo.cos_phi[0] + posy * sta[ista].backInfo.sin_phi[0];
- z[ista][iVec] = posz;
- t[ista][iVec] = time;
sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
if (i == 0) {
- z_start[iVec] = posz;
- x_first[iVec] = x[ista][iVec];
- y_first[iVec] = y[ista][iVec];
-
- auto [dxx, dxy, dyy] = sta[ista].FormXYCovarianceMatrix((fscal) du2[ista][iVec], (fscal) dv2[ista][iVec]);
-
- fstC00[iVec] = dxx;
- fstC10[iVec] = dxy;
- fstC11[iVec] = dyy;
+ z_start[iVec] = z[ista][iVec];
+ x_first[iVec] = x[ista][iVec];
+ y_first[iVec] = y[ista][iVec];
+ t_first[iVec] = t[ista][iVec];
+ cov_first.C00[iVec] = cov[ista].C00[iVec];
+ cov_first.C10[iVec] = cov[ista].C10[iVec];
+ cov_first.C11[iVec] = cov[ista].C11[iVec];
+ dt2_first[iVec] = dt2[ista][iVec];
}
if (i == nHitsTrack - 1) {
- z_end[iVec] = posz;
- x_last[iVec] = x[ista][iVec];
- y_last[iVec] = y[ista][iVec];
-
- auto [dxx, dxy, dyy] = sta[ista].FormXYCovarianceMatrix((fscal) du2[ista][iVec], (fscal) dv2[ista][iVec]);
-
- lstC00[iVec] = dxx;
- lstC10[iVec] = dxy;
- lstC11[iVec] = dyy;
+ z_end[iVec] = z[ista][iVec];
+ x_last[iVec] = x[ista][iVec];
+ y_last[iVec] = y[ista][iVec];
+ t_last[iVec] = t[ista][iVec];
+ cov_last.C00[iVec] = cov[ista].C00[iVec];
+ cov_last.C10[iVec] = cov[ista].C10[iVec];
+ cov_last.C11[iVec] = cov[ista].C11[iVec];
+ dt2_last[iVec] = dt2[ista][iVec];
}
}
}
@@ -327,7 +311,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
i = 0;
- FilterFirst(fit, x_first, y_first, fstC00, fstC10, fstC11);
+ FilterFirst(fit, x_first, y_first, t_first, cov_first, dt2_first);
fit.SetQp0(fit.Tr().qp);
z1 = z[i];
@@ -354,8 +338,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.EnergyLossCorrection(radThick, -fvec::One());
fit.SetMask(initialised && w[i]);
- fit.Filter(sta[i].frontInfo, u[i], du2[i]);
- fit.Filter(sta[i].backInfo, v[i], dv2[i]);
+ fit.FilterXY(cov[i], x[i], y[i]);
fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
b2 = b1;
@@ -398,7 +381,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
i = nHits - 1;
- FilterFirst(fit, x_last, y_last, lstC00, lstC10, lstC11);
+ FilterFirst(fit, x_last, y_last, t_last, cov_last, dt2_last);
z1 = z[i];
sta[i].fieldSlice.GetFieldValue(T.x, T.y, b1);
@@ -425,8 +408,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.EnergyLossCorrection(radThick, fvec::One());
fit.SetMask(initialised && w[i]);
- fit.Filter(sta[i].frontInfo, u[i], du2[i]);
- fit.Filter(sta[i].backInfo, v[i], dv2[i]);
+ fit.FilterXY(cov[i], x[i], y[i]);
fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
b2 = b1;
@@ -465,14 +447,6 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
tr[iVec]->SetNDF(static_cast<int>(T.NDF[iVec]));
}
}
-
- delete[] x;
- delete[] u;
- delete[] v;
- delete[] y;
- delete[] z;
- delete[] w;
- delete[] fB;
}
void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, const vector<int>& pidHypo)
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index b47941903b5a7a80ac354852e4af053b66d20530..6c3b9e75aaf183faa2ae81c642a04d52e0898ec7 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -66,6 +66,8 @@ set(SRCS
L1Algo/utils/L1AlgoDraw.cxx
L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
L1Algo/utils/L1AlgoPulls.cxx
+ L1Algo/utils/CaUvConverter.cxx
+
catools/CaToolsMCData.cxx
catools/CaToolsMCPoint.cxx
catools/CaToolsMCTrack.cxx
@@ -108,6 +110,7 @@ set(HEADERS
L1Algo/L1Vector.h
L1Algo/L1EArray.h
L1Algo/L1Undef.h
+ L1Algo/utils/CaUvConverter.h
catools/CaToolsWindowFinder.h
catools/CaToolsLinkKey.h
catools/CaToolsHitRecord.h
@@ -215,6 +218,7 @@ install(FILES CbmL1Counters.h
L1Algo/L1Track.h
utils/CbmCaIdealHitProducer.h
utils/CbmCaIdealHitProducerDet.h
+ L1Algo/utils/CaUvConverter.cxx
qa/CbmCaInputQaBase.h
DESTINATION include
)
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index 63cdb5f2ab9ac6795223fc3139ccb8b76df8a56c..b8c99253113c01ea16e37a9e47c34bb74e56b6ff 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -390,17 +390,18 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
// Save the algo hit
if (fpIODataManager.get()) {
L1Hit aHit;
- aHit.iSt = hitRecord.fStaId;
aHit.f = hitRecord.fStripF;
aHit.b = hitRecord.fStripB;
- aHit.ID = static_cast<int>(fpIODataManager->GetNofHits());
+ aHit.x = hitRecord.fX;
+ aHit.y = hitRecord.fY;
aHit.z = hitRecord.fZ;
- aHit.u = hitRecord.fU;
- aHit.v = hitRecord.fV;
aHit.t = hitRecord.fT;
+ aHit.dx = hitRecord.fDx;
+ aHit.dxy = hitRecord.fDxy;
+ aHit.dy = hitRecord.fDy;
aHit.dt2 = hitRecord.fDt * hitRecord.fDt;
- aHit.du2 = hitRecord.fDu * hitRecord.fDu;
- aHit.dv2 = hitRecord.fDv * hitRecord.fDv;
+ aHit.ID = static_cast<int>(fpIODataManager->GetNofHits());
+ aHit.iSt = hitRecord.fStaId;
fpIODataManager->PushBackHit(aHit, hitRecord.fDataStream);
}
diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h
index d8a67b8324fee986ca656f2821b3bb4c260e99dc..06c2c374dbd2547c3789d87bc313a2a5c00f98d6 100644
--- a/reco/L1/CbmCaTimeSliceReader.h
+++ b/reco/L1/CbmCaTimeSliceReader.h
@@ -240,8 +240,6 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector()
iStGeom = pDetInterface->GetTrackingStationIndex(pStsHit->GetAddress());
hitRecord.fStripF = fFirstHitKey + pStsHit->GetFrontClusterId();
hitRecord.fStripB = fFirstHitKey + pStsHit->GetBackClusterId();
- hitRecord.fDu = pStsHit->GetDu();
- hitRecord.fDv = pStsHit->GetDv();
}
else if constexpr (L1DetectorID::kMuch == DetID) {
CbmMuchPixelHit* pMuchHit = static_cast<CbmMuchPixelHit*>(fvpBrHits[DetID]->At(iHext));
@@ -284,10 +282,6 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector()
hitRecord.fDt = pPixelHit->GetTimeError();
hitRecord.fDet = static_cast<int>(DetID);
hitRecord.fDataStream = (static_cast<int64_t>(hitRecord.fDet) << 60) | pPixelHit->GetAddress();
- float phiF = pDetInterface->GetStripsStereoAngleFront(iStGeom);
- float phiB = pDetInterface->GetStripsStereoAngleBack(iStGeom);
- hitRecord.fU = hitRecord.fX * cos(phiF) + hitRecord.fY * sin(phiF);
- hitRecord.fV = hitRecord.fX * cos(phiB) + hitRecord.fY * sin(phiB);
hitRecord.fExtId = iHext;
// Update number of hit keys
@@ -296,8 +290,6 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector()
if (fNofHitKeys <= hitRecord.fStripB) { fNofHitKeys += hitRecord.fStripB; }
}
else {
- hitRecord.fDu = hitRecord.fDx;
- hitRecord.fDv = hitRecord.fDy;
hitRecord.fStripF = fFirstHitKey + iHext;
hitRecord.fStripB = hitRecord.fStripF;
if (fNofHitKeys <= hitRecord.fStripF) { fNofHitKeys += hitRecord.fStripF; }
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 4cf54dc1b3add62c736fe771400da8b0de133f83..30a908259f2ea5087c37e6f812a73d3fa0b08bc4 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -473,9 +473,6 @@ InitStatus CbmL1::Init()
stationInfo.SetRmin(mvdInterface->GetRmin(iSt));
stationInfo.SetRmax(mvdInterface->GetRmax(iSt));
stationInfo.SetZthickness(mvdInterface->GetThickness(iSt));
- // TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
- stationInfo.SetFrontBackStripsGeometry((fscal) mvdInterface->GetStripsStereoAngleFront(iSt),
- (fscal) mvdInterface->GetStripsStereoAngleBack(iSt));
stationInfo.SetTrackingStatus(fTargetZ < stationInfo.GetZdouble() ? true : false);
if (fvmDisabledStationIDs[L1DetectorID::kMvd].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kMvd].cend()) {
stationInfo.SetTrackingStatus(false);
@@ -502,9 +499,6 @@ InitStatus CbmL1::Init()
stationInfo.SetRmax(stsInterface->GetRmax(iSt));
stationInfo.SetZthickness(stsInterface->GetThickness(iSt));
- // TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
- stationInfo.SetFrontBackStripsGeometry((fscal) stsInterface->GetStripsStereoAngleFront(iSt),
- (fscal) stsInterface->GetStripsStereoAngleBack(iSt));
stationInfo.SetTrackingStatus(fTargetZ < stationInfo.GetZdouble() ? true : false);
if (fvmDisabledStationIDs[L1DetectorID::kSts].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kSts].cend()) {
stationInfo.SetTrackingStatus(false);
@@ -530,9 +524,6 @@ InitStatus CbmL1::Init()
stationInfo.SetRmin(muchInterface->GetRmin(iSt));
stationInfo.SetRmax(muchInterface->GetRmax(iSt));
stationInfo.SetZthickness(muchInterface->GetThickness(iSt));
- // TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
- stationInfo.SetFrontBackStripsGeometry((fscal) muchInterface->GetStripsStereoAngleFront(iSt),
- (fscal) muchInterface->GetStripsStereoAngleBack(iSt));
stationInfo.SetTrackingStatus(fTargetZ < stationInfo.GetZdouble() ? true : false);
if (fvmDisabledStationIDs[L1DetectorID::kMuch].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kMuch].cend()) {
stationInfo.SetTrackingStatus(false);
@@ -558,10 +549,7 @@ InitStatus CbmL1::Init()
stationInfo.SetRmin(trdInterface->GetRmin(iSt));
stationInfo.SetRmax(trdInterface->GetRmax(iSt));
stationInfo.SetZthickness(trdInterface->GetThickness(iSt));
- fscal trdFrontPhi = trdInterface->GetStripsStereoAngleFront(iSt);
- fscal trdBackPhi = trdInterface->GetStripsStereoAngleBack(iSt);
if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { stationInfo.SetTimeInfo(false); }
- stationInfo.SetFrontBackStripsGeometry(trdFrontPhi, trdBackPhi);
stationInfo.SetTrackingStatus(fTargetZ < stationInfo.GetZdouble() ? true : false);
if (iSt == 1 && L1Algo::TrackingMode::kMcbm == fTrackingMode && fMissingHits) {
stationInfo.SetTrackingStatus(false);
@@ -591,9 +579,6 @@ InitStatus CbmL1::Init()
stationInfo.SetYmax(tofInterface->GetYmax(iSt));
stationInfo.SetRmin(tofInterface->GetRmin(iSt));
stationInfo.SetRmax(tofInterface->GetRmax(iSt));
- fscal tofFrontPhi = tofInterface->GetStripsStereoAngleFront(iSt);
- fscal tofBackPhi = tofInterface->GetStripsStereoAngleBack(iSt);
- stationInfo.SetFrontBackStripsGeometry(tofFrontPhi, tofBackPhi);
stationInfo.SetTrackingStatus(fTargetZ < stationInfo.GetZdouble() ? true : false);
if (fvmDisabledStationIDs[L1DetectorID::kTof].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kTof].cend()) {
stationInfo.SetTrackingStatus(false);
@@ -976,10 +961,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
{
for (L1HitIndex_t i = 0; i < fpAlgo->GetInputData().GetNhits(); i++) {
const L1Hit& h = fpAlgo->GetInputData().GetHit(i);
- const L1Station& sta = fpAlgo->GetParameters()->GetStation(h.iSt);
- fscal x, y;
- std::tie(x, y) = sta.ConvUVtoXY<fscal>(h.u, h.v);
- fMaterialMonitor[h.iSt].MarkActiveBin(x, y);
+ fMaterialMonitor[h.iSt].MarkActiveBin(h.x, h.y);
}
}
//
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index e2c580c1a732982ab45b3a3069fdd090e6aea527..81e48c78590cc3fc111c3f257b257e860875605d 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -280,16 +280,6 @@ public:
const L1Vector<int>& GetHitMCRefs() const { return fvHitPointIndexes; }
- void SetUseMcHit(int MvdUseMcHit = 0, int StsUseMcHit = 0, int MuchUseMcHit = 0, int TrdUseMcHit = 0,
- int TofUseMcHit = 0)
- {
- fMvdUseMcHit = MvdUseMcHit;
- fStsUseMcHit = StsUseMcHit;
- fMuchUseMcHit = MuchUseMcHit;
- fTrdUseMcHit = TrdUseMcHit;
- fTofUseMcHit = TofUseMcHit;
- }
-
static double boundedGaus(double sigma);
/// Obsolete setters to be removed
@@ -570,16 +560,6 @@ private:
Double_t fMomentumCutOff = 0.1; // currently not used
Bool_t fGhostSuppression = true; // currently not used
- /// Flags to adjust hits with MC information
- /// 1: Position of a reconstructed hit is taken from matched MC point and smeared (optionally) with the position
- /// resolution
- /// 2: A set of hits is created from the set of MC points
- Int_t fMvdUseMcHit = -1; ///< MC info flag for MVD hits
- Int_t fStsUseMcHit = -1; ///< MC info flag for STS hits
- Int_t fMuchUseMcHit = -1; ///< MC info flag for MuCh hits
- Int_t fTrdUseMcHit = -1; ///< MC info flag for TRD hits
- Int_t fTofUseMcHit = -1; ///< MC info flag for TOF hits
-
bool fUseMVD = false; ///< if Mvd data should be processed
bool fUseSTS = false; ///< if Mvd data should be processed
bool fUseMUCH = false; ///< if Much data should be processed
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 1453dc42cd170c110f6098e807d4017baa51e336..5d90973d45bac9c63823c56978e99909531998ae 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -69,73 +69,16 @@ struct TmpHit {
int iStation; ///< index of station
int64_t fDataStream; ///< global index of detector module
int ExtIndex; ///< index of hit in the external TClonesArray array (NOTE: negative for MVD)
- double u; ///< position of hit along axis orthogonal to front strips [cm]
- double v; ///< position of hit along axis orthogonal to back strips [cm]
double x; ///< position of hit along x axis [cm]
double y; ///< position of hit along y axis [cm]
double z; ///< position of hit along z axis [cm]
double dx; ///< hit position uncertainty along x axis [cm]
double dy; ///< hit position uncertainty along y axis [cm]
double dxy; ///< hit position covariance along x and y axes [cm2]
- double du; ///< hit position uncertainty along axis orthogonal to front strips [cm]
- double dv; ///< hit position uncertainty along axis orthogonal to back strips [cm]
int iMC; ///< index of MCPoint in the fvMCPoints array
double time = 0.; ///< time of the hit [ns]
double dt = 1.e4; ///< time error of the hit [ns]
int Det;
-
- /// Creates a hit from the CbmL1MCPoint object
- /// \param point constant reference to the input MC-point
- /// \param det
- /// \param nStripF
- /// \param ip
- /// \param NStrips
- /// \param st reference to the station info object
- void CreateHitFromPoint(const CbmL1MCPoint& point, int ip, int det, int& NStrips, const L1Station& st, double du_,
- double dv_, double dt_, bool doSmear)
- {
- ExtIndex = 0;
- Det = det;
- iStation = point.iStation;
- fDataStream = (static_cast<int64_t>(Det) << 60);
- dt = dt_;
- time = point.time;
-
- iStripF = NStrips;
- iStripB = iStripF;
- NStrips++;
-
- du = du_;
- dv = dv_;
-
- double c00, c10, c11;
- std::tie(c00, c10, c11) = st.FormXYCovarianceMatrix(du_ * du_, dv_ * dv_);
-
- dx = sqrt(c00);
- dy = sqrt(c11);
- dxy = c10;
-
- SetHitFromPoint(point, ip, st, doSmear);
- }
-
- /// Sets randomized position and time of the hit
- /// The positions are smeared within predefined errors dx, dy, dt; z coordinate
- /// of the hit is known precisely
- /// \param point constant reference to the input MC-point
- /// \param st reference to the station info object
- ///
- void SetHitFromPoint(const CbmL1MCPoint& point, int ip, const L1Station& st, bool doSmear = 1)
- {
- std::tie(u, v) = st.ConvXYtoUV<double>(point.x, point.y);
- if (doSmear) {
- u += CbmL1::boundedGaus(du);
- v += CbmL1::boundedGaus(dv);
- time += CbmL1::boundedGaus(dt);
- }
- std::tie(x, y) = st.ConvUVtoXY<double>(u, v);
- z = point.z;
- iMC = ip;
- }
};
@@ -198,12 +141,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
int nTrdHits = 0;
int nTofHits = 0;
- int firstMvdPoint = 0;
- int firstStsPoint = 0;
- int firstMuchPoint = 0;
- int firstTrdPoint = 0;
- int firstTofPoint = 0;
-
/*
* MC hits and tracks gathering: START
*
@@ -238,7 +175,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
for (DFSET::iterator set_it = fvSelectedMcEvents.begin(); set_it != fvSelectedMcEvents.end(); ++set_it) {
Int_t iFile = set_it->first;
Int_t iEvent = set_it->second;
- firstMvdPoint = fvMCPoints.size();
if (fUseMVD && fpMvdPoints) {
Int_t fNpointsMvdInEvent = fpMvdPoints->Size(iFile, iEvent);
double maxDeviation = 0;
@@ -276,7 +212,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
// assert(fabs(maxDeviation) < 1.);
} // fpMvdPoints
- firstStsPoint = fvMCPoints.size();
if (fUseSTS && fpStsPoints) {
Int_t nMC = fpStsPoints->Size(iFile, iEvent);
@@ -314,7 +249,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
// assert(fabs(maxDeviation) < 1.);
} // fpStsPoints
- firstMuchPoint = fvMCPoints.size();
if (fUseMUCH && fpMuchPoints) {
Int_t nMC = fpMuchPoints->Size(iFile, iEvent);
for (Int_t iMC = 0; iMC < nMC; iMC++) {
@@ -335,8 +269,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
}
} // fpMuchPoints
- firstTrdPoint = fvMCPoints.size();
-
if (fUseTRD && fpTrdPoints) {
for (Int_t iMC = 0; iMC < fpTrdPoints->Size(iFile, iEvent); iMC++) {
CbmL1MCPoint MC;
@@ -362,7 +294,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
} // fpTrdPoints
- firstTofPoint = fvMCPoints.size();
if (fUseTOF && fpTofPoints) {
vector<float> vTofPointToTrackdZ[fNTofStationsGeom]; // active stations!
@@ -491,7 +422,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
//
// get MVD hits
- if (fUseMVD && fpMvdHits && (2 != fMvdUseMcHit)) {
+ if (fUseMVD && fpMvdHits) {
int firstDetStrip = NStrips;
@@ -525,27 +456,18 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.dx = h->GetDx();
th.dy = h->GetDy();
- th.du = h->GetDx();
- th.dv = h->GetDy();
th.dxy = h->GetDxy();
th.x = pos.X();
th.y = pos.Y();
th.z = pos.Z();
- const L1Station& st = fpAlgo->GetParameters()->GetStation(th.iStation);
- th.u = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.v = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
-
// Get time
th.time = h->GetTime(); // currently ignored by the tracking
th.dt = h->GetTimeError(); // currently ignored by the tracking
}
- th.Det = 0;
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMvd>(hitIndex) : -1;
- if (1 == fMvdUseMcHit && th.iMC >= 0) {
- th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation));
- }
+ th.Det = 0;
+ th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMvd>(hitIndex) : -1;
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
//if( th.iMC >=0 ) // DEBUG !!!!
{
@@ -555,24 +477,10 @@ void CbmL1::ReadEvent(CbmEvent* event)
} // for j
} // if fpMvdHits
- if (fUseMVD && (2 == fMvdUseMcHit)) { // create hits from points
-
- for (int ip = firstMvdPoint; ip < firstMvdPoint + fNpointsMvd; ip++) {
- const CbmL1MCPoint& p = fvMCPoints[ip];
- TmpHit th;
- int DetId = 0;
- double du = 5.e-4;
- double dt = 1000.;
- th.CreateHitFromPoint(p, ip, DetId, NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du, du, dt, true);
- tmpHits.push_back(th);
- nMvdHits++;
- }
- }
-
//
// Get STS hits
//
- if (fUseSTS && fpStsHits && (2 != fStsUseMcHit)) {
+ if (fUseSTS && fpStsHits) {
Int_t nEntSts = (event ? event->GetNofData(ECbmDataType::kStsHit) : fpStsHits->GetEntriesFast());
@@ -621,19 +529,8 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.dx = h->GetDx();
th.dy = h->GetDy();
th.dxy = h->GetDxy();
-
- th.du = h->GetDu();
- th.dv = h->GetDv();
-
- const L1Station& st = fpAlgo->GetParameters()->GetStation(th.iStation);
- th.u = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.v = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
- }
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kSts>(hitIndex) : -1;
-
- if (1 == fStsUseMcHit && th.iMC >= 0) {
- th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation));
}
+ th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kSts>(hitIndex) : -1;
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
tmpHits.push_back(th);
nStsHits++;
@@ -641,25 +538,12 @@ void CbmL1::ReadEvent(CbmEvent* event)
} // for j
} // if fpStsHits
- if (fUseSTS && (2 == fStsUseMcHit)) { // create hits from points
- for (int ip = firstStsPoint; ip < firstStsPoint + fNpointsSts; ip++) {
- const CbmL1MCPoint& p = fvMCPoints[ip];
- TmpHit th;
- int DetId = 1;
- double du = 10.e-4;
- double dt = 5.;
- th.CreateHitFromPoint(p, ip, DetId, NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du, du, dt, true);
- tmpHits.push_back(th);
- nStsHits++;
- }
- }
-
//
// Get MuCh hits
//
- if (fUseMUCH && fpMuchPixelHits && (2 != fMuchUseMcHit)) {
+ if (fUseMUCH && fpMuchPixelHits) {
Int_t nEnt = (event ? event->GetNofData(ECbmDataType::kMuchPixelHit) : fpMuchPixelHits->GetEntriesFast());
@@ -699,42 +583,8 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.dx = h->GetDx();
th.dy = h->GetDy();
th.dxy = 0; /// FIXME: ???
-
- th.du = h->GetDx();
- th.dv = h->GetDy();
-
-
- const L1Station& st = fpAlgo->GetParameters()->GetStation(th.iStation);
- th.u = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.v = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
}
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMuch>(th.ExtIndex) : -1;
- if (1 == fMuchUseMcHit && th.iMC > -1) {
- th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation));
- }
- //int iMC = -1;
- //if (fPerformance) {
- // if (fpMuchHitMatches) {
- // CbmMatch* matchHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fpMuchHitMatches->At(j));
-
-
- // for (Int_t iLink = 0; iLink < matchHitMatch->GetNofLinks(); iLink++) {
- // if (matchHitMatch->GetLink(iLink).GetIndex() < fNpointsMuch) {
- // iMC = matchHitMatch->GetLink(iLink).GetIndex();
- // Int_t iIndex = iMC + fNpointsMvd + fNpointsSts;
-
- // Int_t iFile = matchHitMatch->GetLink(0).GetFile();
- // Int_t iEvent = matchHitMatch->GetLink(0).GetEntry();
-
- // auto itPoint = fmMCPointsLinksMap.find(CbmL1LinkKey(iIndex, iEvent, iFile));
- // if (itPoint == fmMCPointsLinksMap.cend()) continue;
- // th.iMC = itPoint->second;
- // if ((1 == fMuchUseMcHit) && (th.iMC > -1))
- // th.SetHitFromPoint(fvMCPoints[th.iMC], fpAlgo->GetParameters()->GetStation(th.iStation));
- // }
- // }
- // }
- //}
+ th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMuch>(th.ExtIndex) : -1;
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
tmpHits.push_back(th);
nMuchHits++;
@@ -742,33 +592,12 @@ void CbmL1::ReadEvent(CbmEvent* event)
} // for j
} // if listMuchHits
- if ((2 == fMuchUseMcHit) && fUseMUCH) { // create hits from points
-
- for (int ip = firstMuchPoint; ip < firstMuchPoint + fNpointsMuch; ip++) {
- const CbmL1MCPoint& p = fvMCPoints[ip];
-
- // int mcTrack = p.ID;
- // if (mcTrack < 0) continue;
- // const CbmL1MCTrack& t = fvMCTracks[mcTrack];
- //if (t.p < 1) continue;
- // if (t.Points.size() > 4) continue;
-
- TmpHit th;
- int DetId = 2;
- double du = 100.e-4;
- double dt = 3.9;
- th.CreateHitFromPoint(p, ip, DetId, NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du, du, dt, true);
- tmpHits.push_back(th);
- nMuchHits++;
- }
- }
-
//
// Get TRD hits
//
- if (fUseTRD && fpTrdHits && (2 != fTrdUseMcHit)) {
+ if (fUseTRD && fpTrdHits) {
assert(fpTrdHits);
@@ -821,70 +650,19 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.dy = fabs(h->GetDy());
th.dxy = 0;
- th.du = fabs(h->GetDx());
- th.dv = fabs(h->GetDy());
-
- const L1Station& st = fpAlgo->GetParameters()->GetStation(th.iStation);
-
- std::tie(th.u, th.v) = st.ConvXYtoUV<double>(th.x, th.y);
-
th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTrd>(hitIndex) : -1;
- if (1 == fTrdUseMcHit) { // replace hit by MC point
-
- assert(fPerformance && fpTrdHitMatches);
-
- if (th.iMC < 0) continue; // skip noise hits
- int iMcTrd = th.iMC - fNpointsMvd - fNpointsSts - fNpointsMuch;
- assert(iMcTrd >= 0 && iMcTrd < fNpointsTrd);
- if (mcUsed[iMcTrd]) continue; // one hit per MC point
- bool doSmear = true;
- if (1) { // SGtrd2d!! set artificial errors when the input errors are unrealistically large
- if (th.dx > 0.3) { th.dx = 0.3; }
- if (th.dy > 0.2) { th.dy = 0.2; }
- th.du = th.dx;
- th.dv = th.dy;
- }
- th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation), doSmear);
- mcUsed[iMcTrd] = 1;
- } // replace hit by MC points
-
- if (0) { // select specific tracks
- int mcTrack = -1;
- float mcP = 0;
- if (th.iMC >= 0) {
- mcTrack = fvMCPoints[th.iMC].ID;
- if (mcTrack >= 0) { mcP = fvMCTracks[mcTrack].p; }
- }
- if (mcP < 1.) continue;
- }
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
tmpHits.push_back(th);
nTrdHits++;
} // for fpTrdHits
} // read TRD hits
- if (fUseTRD && (2 == fTrdUseMcHit)) { // create hits from MC points
- for (int ip = firstTrdPoint; ip < firstTrdPoint + fNpointsTrd; ip++) {
- const CbmL1MCPoint& p = fvMCPoints[ip];
- // int mcTrack = p.ID;
- // if (mcTrack < 0) continue;
- // const CbmL1MCTrack& t = fvMCTracks[mcTrack];
- TmpHit th;
- int DetId = 3;
- double du = 0.1;
- double dt = 10.;
- th.CreateHitFromPoint(p, ip, DetId, NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du, du, dt, true);
- tmpHits.push_back(th);
- nTrdHits++;
- }
- }
-
//
// Get ToF hits
//
- if (fUseTOF && fpTofHits && (2 != fTofUseMcHit)) {
+ if (fUseTOF && fpTofHits) {
int firstDetStrip = NStrips;
Int_t nEntTof = (event ? event->GetNofData(ECbmDataType::kTofHit) : fpTofHits->GetEntriesFast());
@@ -914,10 +692,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.dx = h->GetDx();
th.dy = h->GetDy();
- th.dxy = 0;
-
- th.du = h->GetDx();
- th.dv = h->GetDy();
+ th.dxy = h->GetDxy();
// th.iSector = 0;
th.iStripF = firstDetStrip + j;
@@ -933,21 +708,15 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.z = pos.Z();
- if (th.z > 400) continue; // is it still needed here? (S.Zharko)
+ //TODO: is it still needed here? (S.Zharko)
+ if (L1Algo::TrackingMode::kMcbm == fTrackingMode && th.z > 400) continue;
int stIdx = fpAlgo->GetParameters()->GetStationIndexActive(sttof, L1DetectorID::kTof);
if (stIdx == -1) continue;
th.iStation = stIdx;
- const L1Station& st = fpAlgo->GetParameters()->GetStation(th.iStation);
- th.u = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.v = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
-
th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTof>(hitIndex) : -1;
- if (1 == fTofUseMcHit && th.iMC > -1) {
- th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation));
- }
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
tmpHits.push_back(th);
nTofHits++;
@@ -955,26 +724,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
} // for j
} // if listTofHits
- if ((2 == fTofUseMcHit) && fUseTOF) { // create hits from points
- for (int ip = firstTofPoint; ip < firstTofPoint + fNpointsTof; ip++) {
- const CbmL1MCPoint& p = fvMCPoints[ip];
-
- // int mcTrack = p.ID;
- // if (mcTrack < 0) continue;
- //const CbmL1MCTrack& t = fvMCTracks[mcTrack];
- //if (t.p < 1) continue;
- //if (t.Points.size() > 4) continue;
-
- TmpHit th;
- int DetId = 4;
- double du = 0.1;
- double dt = 0.075;
- th.CreateHitFromPoint(p, ip, DetId, NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du, du, dt, true);
- tmpHits.push_back(th);
- nTofHits++;
- }
- }
-
if (fVerbose > 1) {
LOG(info) << "L1 ReadEvent: nhits mvd " << nMvdHits << " sts " << nStsHits << " much " << nMuchHits << " trd "
@@ -1018,18 +767,19 @@ void CbmL1::ReadEvent(CbmEvent* event)
assert(th.iStripB >= 0 || th.iStripB < NStrips);
L1Hit h;
- h.iSt = th.iStation;
+
h.f = th.iStripF;
h.b = th.iStripB;
- h.ID = iHit;
+ h.x = th.x;
+ h.y = th.y;
h.z = th.z;
- h.u = th.u;
- h.v = th.v;
h.t = th.time;
+ h.dx = th.dx;
+ h.dy = th.dy;
+ h.dxy = th.dxy;
h.dt2 = th.dt * th.dt;
- h.du2 = th.du * th.du;
- h.dv2 = th.dv * th.dv;
-
+ h.ID = iHit;
+ h.iSt = th.iStation;
// save hit
fvExternalHits.push_back(CbmL1HitId(th.Det, th.ExtIndex));
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index d454ddad23e597bb97d73b4e184b77917902c73d..3a7bb63f035ac9a8d3e4c55e5b196744f85da42b 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -53,6 +53,12 @@ void L1Algo::SetNThreads(unsigned int n)
fTripletHitR_dUfront[i].reserve(L1Constants::size::kMaxPortionTripletsP);
fTripletHitR_dUback[i].reserve(L1Constants::size::kMaxPortionTripletsP);
+ fTripletHitR_X[i].reserve(L1Constants::size::kMaxPortionTripletsP);
+ fTripletHitR_Y[i].reserve(L1Constants::size::kMaxPortionTripletsP);
+ fTripletHitR_dX[i].reserve(L1Constants::size::kMaxPortionTripletsP);
+ fTripletHitR_dXY[i].reserve(L1Constants::size::kMaxPortionTripletsP);
+ fTripletHitR_dY[i].reserve(L1Constants::size::kMaxPortionTripletsP);
+
for (unsigned int j = 0; j < L1Constants::size::kMaxNstations; j++) {
fTriplets[j][i].SetName(std::stringstream() << "L1Algo::fTriplets[" << i << "][" << j << "]");
}
@@ -169,33 +175,12 @@ void L1Algo::ReceiveParameters(L1Parameters&& parameters)
L1FieldRegion::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
}
-/// TODO: Move to L1Hit
-std::pair<fscal, fscal> L1Algo::GetHitCoorOnGrid(const L1Hit& h, char iS)
-{
- const L1Station& sta = fParameters.GetStation(int(iS));
- fscal x, y;
- std::tie(x, y) = sta.ConvUVtoXY<fscal>(h.u, h.v);
- float dz = h.z - fParameters.GetTargetPositionZ()[0];
- return {x / dz, y / dz};
-}
-
-/// TODO: Move to L1Hit
-void L1Algo::GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta)
-{
- std::tie(_x, _y) = sta.ConvUVtoXY<fscal>(_h.u, _h.v);
- _z = _h.z;
-}
-
-L1HitPoint L1Algo::CreateHitPoint(const L1Hit& hit)
-{
- /// full the hit point by hit information: takes hit as input (2 strips)
- /// and creates hit_point with all coordinates (x,y,z,u,v,t);
- return L1HitPoint(hit.z, hit.u, hit.v, hit.t, hit.du2, hit.dv2, hit.dt2);
-}
-
-void L1Algo::CreateHitPoint(const L1Hit& hit, L1HitPoint& point)
+std::pair<fscal, fscal> L1Algo::GetHitCoorOnGrid(const L1Hit& h)
{
- point.Set(hit.z, hit.u, hit.v, hit.t, hit.du2, hit.dv2, hit.dt2);
+ float dx = h.x - fParameters.GetTargetPositionX()[0];
+ float dy = h.y - fParameters.GetTargetPositionY()[0];
+ float dz = h.z - fParameters.GetTargetPositionZ()[0];
+ return {dx / dz, dy / dz};
}
int L1Algo::GetMcTrackIdForCaHit(int iHit) const
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 2500a4425738567fe36ec26af22759a1959c8761..88ac58a8724a4f34969c0bb1cd5a9950ddef1be5 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -170,16 +170,10 @@ public:
/// Gets pointer to input data object for external access
const L1InputData& GetInputData() const { return fInputData; }
- /// ----- Hit-point-strips conversion routines ------
+ /// Hit coordinates on the grid
+ std::pair<fscal, fscal> GetHitCoorOnGrid(const L1Hit& h);
- void GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta);
- std::pair<fscal, fscal> GetHitCoorOnGrid(const L1Hit& h, char iS);
-
-
- L1HitPoint CreateHitPoint(const L1Hit& hit); // full the hit point by hit information.
-
- void CreateHitPoint(const L1Hit& hit, L1HitPoint& point);
inline int PackIndex(const int& a, const int& b, const int& c);
inline int UnPackIndex(const int& i, int& a, int& b, int& c);
@@ -271,16 +265,14 @@ public:
void findSingletsStep0( // input
Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
// output
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* t_l, fvec* dt2_l, fvec* du2_l,
- fvec* dv2_l);
+ L1HitIndex_t* hitsl, fvec* x_l, fvec* y_l, fvec* z_l, fvec* t_l, fvec* dx_l, fvec* dxy_l, fvec* dy_l, fvec* dt2_l);
/// Get the field approximation. Add the target to parameters estimation. Propagate to middle station.
void findSingletsStep1( // input
- int istal, int istam, Tindex n1_V,
-
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* t_l, fvec* dt2_l,
+ int istal, int istam, Tindex n1_V, fvec* x_l, fvec* y_l, fvec* zPos_l, fvec* t_l, fvec* dx_l, fvec* dxy_l,
+ fvec* dy_l, fvec* dt2_l,
// output
- L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* du2_l, fvec* dv2_l);
+ L1TrackPar* T_1, L1FieldRegion* fld_1);
/// Find the doublets. Reformat data in the portion of doublets.
void findDoubletsStep0( // input
@@ -303,13 +295,15 @@ public:
// output
Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
- L1Vector<L1HitIndex_t>& hitsr_3, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3,
- L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3);
+ L1Vector<L1HitIndex_t>& hitsr_3,
+
+ L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3,
+ L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3);
/// Add the right hits to parameters estimation.
void findTripletsStep1( // input
- Tindex n3_V, const L1Station& star, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3,
- L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3,
+ Tindex n3_V, const L1Station& star, L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3,
+ L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3, L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3,
// output
L1Vector<L1TrackPar>& T_3);
@@ -501,6 +495,12 @@ public:
L1Vector<fvec> fTripletHitR_dUfront[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dUfront"};
L1Vector<fvec> fTripletHitR_dUback[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dUback"};
+ L1Vector<fvec> fTripletHitR_X[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_X"};
+ L1Vector<fvec> fTripletHitR_Y[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Y"};
+ L1Vector<fvec> fTripletHitR_dX[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dX"};
+ L1Vector<fvec> fTripletHitR_dXY[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dXY"};
+ L1Vector<fvec> fTripletHitR_dY[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dY"};
+
// Tindex NHits_l[L1Constants::size::kMaxNstations];
// Tindex NHits_l_P[L1Constants::size::kMaxNstations];
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
index 667badd8e45c9bc2ebc8ce188cd01abbca751d3f..18c9073d431b83234ebb32f0b1ea31aa133398fd 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.cxx
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -291,37 +291,6 @@ void L1BaseStationInfo::SetFieldStatus(int fieldStatus)
fInitController.SetFlag(EInitKey::kFieldStatus);
}
-//------------------------------------------------------------------------------------------------------------------------
-//
-void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double backPhi)
-{
- //----- Original code from L1Algo -----------------------------------------------------------------------
- double cFront = cos(frontPhi);
- double sFront = sin(frontPhi);
- double cBack = cos(backPhi);
- double sBack = sin(backPhi);
-
- // NOTE: Here additional double variables are used to save the precission
- fL1Station.frontInfo.cos_phi = cFront;
- fL1Station.frontInfo.sin_phi = sFront;
-
- fL1Station.backInfo.cos_phi = cBack;
- fL1Station.backInfo.sin_phi = sBack;
-
- double det = cFront * sBack - sFront * cBack;
-
- assert(fabs(det) > 1.e-2);
-
- fL1Station.xInfo.cos_phi = sBack / det;
- fL1Station.xInfo.sin_phi = -sFront / det;
-
- fL1Station.yInfo.cos_phi = -cBack / det;
- fL1Station.yInfo.sin_phi = cFront / det;
- //-------------------------------------------------------------------------------------------------------
-
- fInitController.SetFlag(EInitKey::kStripsFrontPhi);
- fInitController.SetFlag(EInitKey::kStripsBackPhi);
-}
//------------------------------------------------------------------------------------------------------------------------
//
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index 4045f81ac5b0037767b9a7bef5abae37ca06079d..906ebafde32d13f8f633df0b929c9796a732b007 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -55,8 +55,6 @@ public:
kZthickness, ///< Z thickness of the station
kThicknessMap, ///< thickness map of the station (optional?)
kFieldSlice, ///< L1Station.L1FieldSlice object initialization
- kStripsFrontPhi, ///< strips geometry initialization
- kStripsBackPhi, ///<
// The last item is equal to the number of bits in fInitFlags
kEnd
};
@@ -176,11 +174,6 @@ public:
/// of magnetic field components in position
void SetFieldFunction(const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue);
- /// Sets stereo angles for front and back strips
- /// \param f_phi Stereoangle of front strips [rad]
- /// \param b_phi Stereoangle of back strips [rad]
- void SetFrontBackStripsGeometry(double fPhi, double bPhi);
-
/// Sets station thickness and radiation length
/// \param thickness Thickness of station [arb. units]
void SetZthickness(double thickness);
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index ba15b1fcd71fe43bd4e7b25a776c1f67b06ad01b..d7d7e2fd893595b1a058814bbe265b51e1821ed8 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -19,10 +19,9 @@ using std::endl;
/// Fit track
/// t - track with hits
/// T - track params
-/// dir - 0 - forward, 1 - backward
/// qp0 - momentum for extrapolation
/// initialize - should be params ititialized. 1 - yes.
-void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool dir, const fvec qp0,
+void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool upstream, const fvec qp0,
const bool initParams)
{
L1_assert(t.NHits >= 3);
@@ -37,9 +36,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
const L1Vector<L1HitIndex_t>& hits = t.fHits; // array of indeses of hits of current track
const int nHits = t.NHits;
- const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
- const int iFirstHit = (dir) ? nHits - 1 : 0;
- const int iLastHit = (dir) ? 0 : nHits - 1;
+ const signed short int step = -2 * static_cast<int>(upstream) + 1; // increment for station index
+ const int iFirstHit = (upstream) ? nHits - 1 : 0;
+ const int iLastHit = (upstream) ? 0 : nHits - 1;
const L1Hit& hit0 = fInputData.GetHit(hits[iFirstHit]);
const L1Hit& hit1 = fInputData.GetHit(hits[iFirstHit + step]);
@@ -53,29 +52,24 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
const L1Station& sta1 = fParameters.GetStation(ista1);
const L1Station& sta2 = fParameters.GetStation(ista2);
- fvec u0 = hit0.u;
- fvec v0 = hit0.v;
- auto [x0, y0] = sta0.ConvUVtoXY<fvec>(u0, v0);
- fvec z0 = hit0.z;
+ fvec x0 = hit0.x;
+ fvec y0 = hit0.y;
+ fvec z0 = hit0.z;
- fvec u1 = hit1.u;
- fvec v1 = hit1.v;
- auto [x1, y1] = sta1.ConvUVtoXY<fvec>(u1, v1);
- fvec z1 = hit1.z;
+ fvec x1 = hit1.x;
+ fvec y1 = hit1.y;
+ fvec z1 = hit1.z;
- fvec u2 = hit2.u;
- fvec v2 = hit2.v;
- auto [x2, y2] = sta2.ConvUVtoXY<fvec>(u2, v2);
- // fvec z2 = hit2.z;
-
- fvec dzi = fvec(1.) / (z1 - z0);
+ fvec x2 = hit2.x;
+ fvec y2 = hit2.y;
T.x = x0;
T.y = y0;
if (initParams) {
- T.tx = (x1 - x0) * dzi;
- T.ty = (y1 - y0) * dzi;
- T.qp = qp0;
+ fvec dzi = fvec(1.) / (z1 - z0);
+ T.tx = (x1 - x0) * dzi;
+ T.ty = (y1 - y0) * dzi;
+ T.qp = qp0;
}
fit.SetQp0(qp0);
@@ -86,7 +80,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
T.ResetErrors(1., 1., .1, .1, 1., (sta0.timeInfo ? hit0.dt2 : 1.e6), 1.e6);
T.NDF = fvec(2.);
- std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2);
+ T.C00 = hit0.dx * hit0.dx;
+ T.C10 = hit0.dxy;
+ T.C11 = hit0.dy * hit0.dy;
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -101,31 +97,22 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
- //int ista_prev = ista1;
- int ista = ista2;
-
for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
- const L1Hit& hit = fInputData.GetHit(hits[i]);
- //ista_prev = ista;
- ista = hit.iSt;
-
+ const L1Hit& hit = fInputData.GetHit(hits[i]);
+ int ista = hit.iSt;
const L1Station& sta = fParameters.GetStation(ista);
fit.Extrapolate(hit.z, fld);
- fit.Filter(sta.frontInfo, hit.u, hit.du2);
- fit.Filter(sta.backInfo, hit.v, hit.dv2);
- fit.FilterTime(hit.t, hit.dt2, sta.timeInfo);
+ fit.FilterHit(sta, hit);
auto radThick = fParameters.GetMaterialThickness(ista, T.x, T.y);
fit.AddMsInMaterial(radThick);
- fit.EnergyLossCorrection(radThick, fvec(-1.f));
-
- fldB0 = fldB1;
- fldB1 = fldB2;
- fldZ0 = fldZ1;
- fldZ1 = fldZ2;
- auto [x, y] = sta.ConvUVtoXY<fvec>(hit.u, hit.v);
- sta.fieldSlice.GetFieldValue(x, y, fldB2);
+ fit.EnergyLossCorrection(radThick, fvec(upstream ? 1. : -1.));
+ fldB0 = fldB1;
+ fldB1 = fldB2;
+ fldZ0 = fldZ1;
+ fldZ1 = fldZ2;
+ sta.fieldSlice.GetFieldValue(hit.x, hit.y, fldB2);
fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
} // i
@@ -134,23 +121,21 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
} // void L1Algo::BranchFitterFast
/// like BranchFitterFast but more precise
-void L1Algo::BranchFitter(const L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0, const bool initParams)
+void L1Algo::BranchFitter(const L1Branch& t, L1TrackPar& T, const bool upstream, const fvec qp0, const bool initParams)
{
- BranchFitterFast(t, T, dir, qp0, initParams);
+ BranchFitterFast(t, T, upstream, qp0, initParams);
for (int i = 0; i < 1; i++) {
- BranchFitterFast(t, T, !dir, T.qp, false);
- BranchFitterFast(t, T, dir, T.qp, false);
+ BranchFitterFast(t, T, !upstream, T.qp, false);
+ BranchFitterFast(t, T, upstream, T.qp, false);
}
} // void L1Algo::BranchFitter
/// Find additional hits for existing track
/// t - track with hits
/// T - track params
-/// dir - 0 - forward, 1 - backward
/// qp0 - momentum for extrapolation
/// initialize - should be params ititialized. 1 - yes.
-void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
- const fvec qp0) // TODO take into account pipe
+void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool upstream, const fvec qp0)
{
L1Vector<L1HitIndex_t> newHits {"L1TrackExtender::newHits"};
newHits.reserve(fParameters.GetNstationsActive());
@@ -161,8 +146,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
fit.SetTrack(Tout);
fit.SetQp0(qp0);
- const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
- const int iFirstHit = (dir) ? 2 : t.NHits - 3;
+ const signed short int step = -2 * static_cast<int>(upstream) + 1; // increment for station index
+ const int iFirstHit = (upstream) ? 2 : t.NHits - 3;
// int ista = fInputData.GetHit(t.Hits[iFirstHit]).iSt + 2 * step; // current station. set to the end of track
const L1Hit& hit0 = fInputData.GetHit(t.fHits[iFirstHit]); // optimize
@@ -177,17 +162,14 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
const L1Station& sta1 = fParameters.GetStation(ista1);
const L1Station& sta2 = fParameters.GetStation(ista2);
- fvec u0 = hit0.u;
- fvec v0 = hit0.v;
- auto [x0, y0] = sta0.ConvUVtoXY<fvec>(u0, v0);
+ fvec x0 = hit0.x;
+ fvec y0 = hit0.y;
- fvec u1 = hit1.u;
- fvec v1 = hit1.v;
- auto [x1, y1] = sta1.ConvUVtoXY<fvec>(u1, v1);
+ fvec x1 = hit1.x;
+ fvec y1 = hit1.y;
- fvec u2 = hit2.u;
- fvec v2 = hit2.v;
- auto [x2, y2] = sta2.ConvUVtoXY<fvec>(u2, v2);
+ fvec x2 = hit2.x;
+ fvec y2 = hit2.y;
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -253,21 +235,18 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
//L1_SHOW(hit.b);
if (fvHitKeyFlags[hit.f] || fvHitKeyFlags[hit.b]) { continue; }
- fscal xh, yh, zh;
- GetHitCoor(hit, xh, yh, zh, sta); // faster
-
fvec y, C11;
- fit.ExtrapolateYC11Line(zh, y, C11);
+ fit.ExtrapolateYC11Line(hit.z, y, C11);
// fscal dym_est = ( fPickGather * sqrt(fabs(C11[0]+sta.XYInfo.C11[0])) );
// fscal y_minus_new = y[0] - dym_est;
// if (yh < y_minus_new) continue; // CHECKME take into account overlaping?
fvec x, C00;
- fit.ExtrapolateXC00Line(zh, x, C00);
+ fit.ExtrapolateXC00Line(hit.z, x, C00);
- fscal d_x = xh - x[0];
- fscal d_y = yh - y[0];
+ fscal d_x = hit.x - x[0];
+ fscal d_y = hit.y - y[0];
fscal d2 = d_x * d_x + d_y * d_y;
if (d2 > r2_best) continue;
@@ -277,27 +256,24 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
r2_best = d2;
iHit_best = globalInd;
}
+
if (iHit_best < 0) break;
newHits.push_back(fGridHitIds[iHit_best]);
const L1Hit& hit = fGridHits[iHit_best];
- auto [x, y] = sta.ConvUVtoXY<fvec>(hit.u, hit.v);
-
fit.Extrapolate(hit.z, fld);
- fit.Filter(sta.frontInfo, hit.u, hit.du2);
- fit.Filter(sta.backInfo, hit.v, hit.dv2);
- fit.FilterTime(hit.t, hit.dt2, sta.timeInfo);
+ fit.FilterHit(sta, hit);
auto radThick = fParameters.GetMaterialThickness(ista, tr.x, tr.y);
fit.AddMsInMaterial(radThick);
- fit.EnergyLossCorrection(radThick, dir ? fvec(1.f) : fvec(-1.f));
+ fit.EnergyLossCorrection(radThick, upstream ? fvec(1.f) : fvec(-1.f));
fldB0 = fldB1;
fldB1 = fldB2;
fldZ0 = fldZ1;
fldZ1 = fldZ2;
- sta.fieldSlice.GetFieldValue(x, y, fldB2);
+ sta.fieldSlice.GetFieldValue(hit.x, hit.y, fldB2);
fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
}
@@ -308,7 +284,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
t.fHits.enlarge(NNewHits + NOldHits);
t.NHits = (NNewHits + NOldHits);
- if (dir) { // backward
+ if (upstream) { // backward
for (int i = NOldHits - 1; i >= 0; i--) {
t.fHits[NNewHits + i] = t.fHits[i];
}
@@ -333,21 +309,21 @@ fscal L1Algo::BranchExtender(L1Branch& t) // TODO Simdize
L1TrackPar T;
- // forward
- bool dir = 0;
+ // downstream
+ bool upstream = 0;
- BranchFitter(t, T, dir);
- // BranchFitterFast (t, T, dir, 0, true);
+ BranchFitter(t, T, upstream);
+ // BranchFitterFast (t, T, upstream, 0, true);
// if (t.NHits < minNHits) return T.chi2[0];
- FindMoreHits(t, T, dir, T.qp);
+ FindMoreHits(t, T, upstream, T.qp);
- // backward
- dir = 1;
- BranchFitterFast(t, T, dir, T.qp, false); // 577
+ // upstream
+ upstream = 1;
+ BranchFitterFast(t, T, upstream, T.qp, false);
- FindMoreHits(t, T, dir, T.qp);
+ FindMoreHits(t, T, upstream, T.qp);
return T.chi2[0];
}
diff --git a/reco/L1/L1Algo/L1CaTrackFinder.cxx b/reco/L1/L1Algo/L1CaTrackFinder.cxx
index 113a0cc7e468c4a802bdd8240d1bbaf070af97f5..589b6b3f5efc789b8b9cbebbe32c113dbecb2290 100644
--- a/reco/L1/L1Algo/L1CaTrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinder.cxx
@@ -76,9 +76,8 @@ void L1Algo::CaTrackFinder()
const L1Hit& h = fInputData.GetHit(caHitId);
const L1Station& st = fParameters.GetStation(h.iSt);
- auto [x, y] = st.ConvUVtoXY(h.u, h.v);
- fscal dx = x - targX;
- fscal dy = y - targY;
+ fscal dx = h.x - targX;
+ fscal dy = h.y - targY;
fscal dz = h.z - targZ;
fscal l = sqrt(dx * dx + dy * dy + dz * dz);
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index 09002f5a78a5d361c3abe3a088f4ae32869993b6..3bc126a16b0864f297f73c0bd7ead6a930a4e5bc 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -125,7 +125,7 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
inline void L1Algo::findSingletsStep0( // input
Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
// output
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* t_l, fvec* dt2_l, fvec* du2_l, fvec* dv2_l)
+ L1HitIndex_t* hitsl, fvec* x_l, fvec* y_l, fvec* z_l, fvec* t_l, fvec* dx_l, fvec* dxy_l, fvec* dy_l, fvec* dt2_l)
{
/// Prepare the portion of data of left hits of a triplet:
@@ -137,32 +137,29 @@ inline void L1Algo::findSingletsStep0( // input
if (lastV >= 0) {
// set some positive errors to unfilled part of vectors in order to avoid nans
L1HitPoint& hitl = Hits_l[0];
- du2_l[lastV] = hitl.dU2();
- dv2_l[lastV] = hitl.dV2();
- dt2_l[lastV] = hitl.dT2();
- u_front_l[lastV] = hitl.U();
- u_back_l[lastV] = hitl.V();
+ x_l[lastV] = hitl.X();
+ y_l[lastV] = hitl.Y();
+ z_l[lastV] = hitl.Z();
t_l[lastV] = hitl.T();
- zPos_l[lastV] = hitl.Z();
+ dx_l[lastV] = hitl.dX();
+ dxy_l[lastV] = hitl.dXY();
+ dy_l[lastV] = hitl.dY();
+ dt2_l[lastV] = hitl.dT2();
}
for (Tindex ilh = start_lh, i1 = 0; ilh < end_lh; ++ilh, ++i1) {
- Tindex i1_V = i1 / fvec::size();
- Tindex i1_4 = i1 % fvec::size();
- L1HitPoint& hitl = Hits_l[ilh];
-
-
+ Tindex i1_V = i1 / fvec::size();
+ Tindex i1_4 = i1 % fvec::size();
+ L1HitPoint& hitl = Hits_l[ilh];
+ hitsl[i1] = ilh;
+ x_l[i1_V][i1_4] = hitl.X();
+ y_l[i1_V][i1_4] = hitl.Y();
+ z_l[i1_V][i1_4] = hitl.Z();
t_l[i1_V][i1_4] = hitl.T();
+ dx_l[i1_V][i1_4] = hitl.dX();
+ dxy_l[i1_V][i1_4] = hitl.dXY();
+ dy_l[i1_V][i1_4] = hitl.dY();
dt2_l[i1_V][i1_4] = hitl.dT2();
-
- hitsl[i1] = ilh;
- u_front_l[i1_V][i1_4] = hitl.U();
- u_back_l[i1_V][i1_4] = hitl.V();
-
- du2_l[i1_V][i1_4] = hitl.dU2();
- dv2_l[i1_V][i1_4] = hitl.dV2();
-
- zPos_l[i1_V][i1_4] = hitl.Z();
}
}
@@ -171,10 +168,10 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
int istal,
int istam, /// indexes of left and middle stations of a triplet
Tindex n1_V, ///
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* t_l, fvec* dt2_l,
+ fvec* x_l, fvec* y_l, fvec* z_l, fvec* t_l, fvec* dx_l, fvec* dxy_l, fvec* dy_l, fvec* dt2_l,
// output
// L1TrackPar *T_1,
- L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* du2_l, fvec* dv2_l)
+ L1TrackPar* T_1, L1FieldRegion* fld_1)
{
/// Get the field approximation. Add the target to parameters estimation.
@@ -236,10 +233,9 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
// get the magnetic field along the track.
// (suppose track is straight line with origin in the target)
- fvec& u = u_front_l[i1_V];
- fvec& v = u_back_l[i1_V];
- auto [xl, yl] = stal.ConvUVtoXY<fvec>(u, v);
- fvec zl = zPos_l[i1_V];
+ fvec& xl = x_l[i1_V];
+ fvec& yl = y_l[i1_V];
+ fvec& zl = z_l[i1_V];
fvec& time = t_l[i1_V];
fvec& timeEr2 = dt2_l[i1_V];
const fvec dzli = 1.f / (zl - fTargZ);
@@ -272,6 +268,9 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
fvec qpErr2 = fMaxInvMom * fMaxInvMom / fvec(9.);
T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (stal.timeInfo ? timeEr2 : 1.e6), 1.e2);
+ T.C00 = dx_l[i1_V] * dx_l[i1_V];
+ T.C10 = dxy_l[i1_V];
+ T.C11 = dy_l[i1_V] * dy_l[i1_V];
T.chi2 = 0.;
T.NDF = (fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
@@ -279,12 +278,9 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
// TODO: iteration parameter: "Starting NDF of track parameters"
// NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target
- { // add the target constraint
-
- std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(du2_l[i1_V], dv2_l[i1_V]);
+ // add the target constraint
- fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fld0);
- }
+ fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fld0);
fit.AddMsInMaterial(fParameters.GetMaterialThickness(istal, T.x, T.y));
@@ -391,14 +387,11 @@ inline void L1Algo::findDoubletsStep0(
fvec y, C11;
fit.ExtrapolateYC11Line(zm, y, C11);
- //L1ExtrapolateYC11Line(T1, zm, y, C11);
-
- /// Covariation matrix of the hit
- auto [dxxScalMhit, dxyScalMhit, dyyScalMhit] = stam.FormXYCovarianceMatrix(hitm.dU2(), hitm.dV2());
- fscal dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + dyyScalMhit);
+ fscal dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + hitm.dY2());
- auto [xm, ym] = stam.ConvUVtoXY<fscal>(hitm.U(), hitm.V());
+ fscal xm = hitm.X();
+ fscal ym = hitm.Y();
const fscal dY = ym - y[i1_4];
@@ -408,38 +401,34 @@ inline void L1Algo::findDoubletsStep0(
fvec x, C00;
fit.ExtrapolateXC00Line(zm, x, C00);
- fscal dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + dxxScalMhit);
+ fscal dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + hitm.dX2());
const fscal dX = xm - x[i1_4];
if (dX * dX > dx_est2) continue;
// check chi2
+
fvec C10;
fit.ExtrapolateC10Line(zm, C10);
- fvec chi2 = T1.chi2;
-
- L1Fit::FilterChi2XYC00C10C11(stam.frontInfo, x, y, C00, C10, C11, chi2, hitm.U(), hitm.dU2());
+ auto [chi2x, chi2u] =
+ L1Fit::GetChi2XChi2U(hitm.X(), hitm.Y(), hitm.dX2(), hitm.dXY(), hitm.dY2(), x, y, C00, C10, C11);
+ // TODO: adjust the cut, cut on chi2x & chi2u separately
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2[i1_4] > fDoubletChi2Cut) continue;
+ if (chi2x[i1_4] > fDoubletChi2Cut) continue;
+ if (chi2x[i1_4] + chi2u[i1_4] > fDoubletChi2Cut) continue;
}
- L1Fit::FilterChi2(stam.backInfo, x, y, C00, C10, C11, chi2, hitm.V(), hitm.dV2());
-
// FilterTime(T1, hitm.T(), hitm.dT2());
- if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2[i1_4] > fDoubletChi2Cut) continue;
- }
-
// check if there is a second hit on the same station
{
bool isOtherHit = 0;
fscal dz = hitm.Z() - hitl.Z();
- fscal tu = (hitm.U() - hitl.U()) / dz;
- fscal tv = (hitm.V() - hitl.V()) / dz;
+ fscal tx = (hitm.X() - hitl.X()) / dz;
+ fscal ty = (hitm.Y() - hitl.Y()) / dz;
fscal tt = (hitm.T() - hitl.T()) / dz;
for (unsigned int imh1 = indFirstDoublet; imh1 < hitsm_2.size(); imh1++) {
@@ -450,11 +439,11 @@ inline void L1Algo::findDoubletsStep0(
if (dt * dt > 30. * (hitm.dT2() + hitm1.dT2())) { continue; }
}
- fscal du = hitm.U() + tu * (hitm1.Z() - hitm.Z()) - hitm1.U();
- if (du * du > 20. * (hitm.dU2() + hitm1.dU2())) { continue; }
+ fscal dx = hitm.X() + tx * (hitm1.Z() - hitm.Z()) - hitm1.X();
+ if (dx * dx > 20. * (hitm.dX2() + hitm1.dX2())) { continue; }
- fscal dv = hitm.V() + tv * (hitm1.Z() - hitm.Z()) - hitm1.V();
- if (dv * dv > 30. * (hitm.dV2() + hitm1.dV2())) { continue; }
+ fscal dy = hitm.Y() + ty * (hitm1.Z() - hitm.Z()) - hitm1.Y();
+ if (dy * dy > 30. * (hitm.dY2() + hitm1.dY2())) { continue; }
if (fParameters.DevIsSuppressOverlapHitsViaMc()) {
int indL = fGridHitStartIndex[iStaL] + hitsl_1[i1];
@@ -509,8 +498,8 @@ inline void L1Algo::findTripletsStep0( // input
L1HitIndex_t* hitsl_1, Tindex n2, L1Vector<L1HitIndex_t>& hitsm_2, L1Vector<L1HitIndex_t>& i1_2,
// output
Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
- L1Vector<L1HitIndex_t>& hitsr_3, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3,
- L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3)
+ L1Vector<L1HitIndex_t>& hitsr_3, L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3, L1Vector<fvec>& t_3,
+ L1Vector<fvec>& dx_3, L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3)
{
const L1Station& stal = fParameters.GetStation(iStaL);
const L1Station& stam = fParameters.GetStation(iStaM);
@@ -546,12 +535,15 @@ inline void L1Algo::findTripletsStep0( // input
Tindex n3_V = 0, n3_4 = 0;
T_3.reset(1, L1TrackPar_0);
- u_front_3.reset(1, fvec::Zero());
- u_back_3.reset(1, fvec::Zero());
- z_Pos_3.reset(1, fvec::Zero());
- du2_3.reset(1, fvec::One());
- dv2_3.reset(1, fvec::One());
+
+ x_3.reset(1, fvec::Zero());
+ y_3.reset(1, fvec::Zero());
+ z_3.reset(1, fvec::Zero());
t_3.reset(1, fvec::Zero());
+
+ dx_3.reset(1, fvec::One());
+ dxy_3.reset(1, fvec::Zero());
+ dy_3.reset(1, fvec::One());
dt2_3.reset(1, fvec::One());
assert(iStaR < fParameters.GetNstationsActive()); //TODO SG!!! check if it is needed
@@ -567,13 +559,15 @@ inline void L1Algo::findTripletsStep0( // input
L1FieldRegion f2;
// pack the data
- fvec u_front_2 = 0.f;
- fvec u_back_2 = 0.f;
- fvec du2_2 = 1.f;
- fvec dv2_2 = 1.f;
- fvec zPos_2 = 0.f;
- fvec t_2 = 0.f;
- fvec dt2_2 = 1.f;
+ fvec x_2 = 0.f;
+ fvec y_2 = 0.f;
+ fvec z_2 = 0.f;
+ fvec t_2 = 0.f;
+ L1XYMeasurementInfo cov_2;
+ cov_2.C00 = 1.f;
+ cov_2.C10 = 0.f;
+ cov_2.C11 = 1.f;
+ fvec dt2_2 = 1.f;
size_t n2_4 = 0;
for (; n2_4 < fvec::size() && i2 < n2; i2++, n2_4++) {
@@ -593,14 +587,15 @@ inline void L1Algo::findTripletsStep0( // input
const Tindex& imh = hitsm_2[i2];
const L1HitPoint& hitm = vHits_m[imh];
- u_front_2[n2_4] = hitm.U();
- u_back_2[n2_4] = hitm.V();
- zPos_2[n2_4] = hitm.Z();
- t_2[n2_4] = hitm.T();
- dt2_2[n2_4] = hitm.dT2();
- // num[n2_4] = hitm.track;
- du2_2[n2_4] = hitm.dU2();
- dv2_2[n2_4] = hitm.dV2();
+
+ x_2[n2_4] = hitm.X();
+ y_2[n2_4] = hitm.Y();
+ z_2[n2_4] = hitm.Z();
+ t_2[n2_4] = hitm.T();
+ cov_2.C00[n2_4] = hitm.dX2();
+ cov_2.C10[n2_4] = hitm.dXY();
+ cov_2.C11[n2_4] = hitm.dY2();
+ dt2_2[n2_4] = hitm.dT2();
hitsl_2[n2_4] = hitsl_1[i1];
hitsm_2_tmp[n2_4] = hitsm_2[i2];
@@ -608,10 +603,10 @@ inline void L1Algo::findTripletsStep0( // input
// add the middle hit
- fit.ExtrapolateLineNoField(zPos_2);
+ fit.ExtrapolateLineNoField(z_2);
+
+ fit.FilterXY(cov_2, x_2, y_2);
- fit.Filter(stam.frontInfo, u_front_2, du2_2);
- fit.Filter(stam.backInfo, u_back_2, dv2_2);
fit.FilterTime(t_2, dt2_2, stam.timeInfo);
fit.SetQp0(isMomentumFitted ? fit.Tr().qp : fMaxInvMom);
@@ -675,9 +670,11 @@ inline void L1Algo::findTripletsStep0( // input
}
const fscal zr = hitr.Z();
- // const fscal yr = hitr.Y();
- auto [xr, yr] = star.ConvUVtoXY<fscal>(hitr.U(), hitr.V());
+ fscal xr = hitr.X();
+ fscal yr = hitr.Y();
+
+ //auto [xr, yr] = star.ConvUVtoXY<fscal>(hitr.U(), hitr.V());
L1Fit fit3;
fit3.SetParticleMass(fDefaultMass);
@@ -702,10 +699,7 @@ inline void L1Algo::findTripletsStep0( // input
fvec y, C11;
fit3.ExtrapolateYC11Line(zr, y, C11);
- /// Covariation matrix of the hit
- auto [dxxScalRhit, dxyScalRhit, dyyScalRhit] = star.FormXYCovarianceMatrix(hitr.dU2(), hitr.dV2());
-
- fscal dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + dyyScalRhit)));
+ fscal dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + hitr.dY2())));
const fscal dY = yr - y[i2_4];
const fscal dY2 = dY * dY;
@@ -716,7 +710,7 @@ inline void L1Algo::findTripletsStep0( // input
fvec x, C00;
fit3.ExtrapolateXC00Line(zr, x, C00);
- fscal dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + dxxScalRhit)));
+ fscal dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + hitr.dX2())));
const fscal dX = xr - x[i2_4];
if (dX * dX > dx_est2) continue;
@@ -724,16 +718,13 @@ inline void L1Algo::findTripletsStep0( // input
fvec C10;
fit3.ExtrapolateC10Line(zr, C10);
- fvec chi2 = T2.chi2;
-
- L1Fit::FilterChi2XYC00C10C11(star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U(), hitr.dU2());
-
- L1Fit::FilterChi2(star.backInfo, x, y, C00, C10, C11, chi2, hitr.V(), hitr.dV2());
+ auto [chi2x, chi2u] =
+ L1Fit::GetChi2XChi2U(hitr.X(), hitr.Y(), hitr.dX2(), hitr.dXY(), hitr.dY2(), x, y, C00, C10, C11);
- fit3.FilterTime(hitr.T(), hitr.dT2(), star.timeInfo);
+ //fit3.FilterTime(hitr.T(), hitr.dT2(), star.timeInfo);
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2[i2_4] > fTripletChi2Cut || C00[i2_4] < 0 || C11[i2_4] < 0 || T_cur.C55[i2_4] < 0) {
+ if (chi2x[i2_4] + chi2u[i2_4] > fTripletChi2Cut || T_cur.C55[i2_4] < 0) {
continue; // chi2_triplet < CHI2_CUT
}
}
@@ -768,13 +759,15 @@ inline void L1Algo::findTripletsStep0( // input
T3.SetOneEntry(n3_4, T2, i2_4);
- u_front_3[n3_V][n3_4] = hitr.U();
- u_back_3[n3_V][n3_4] = hitr.V();
- du2_3[n3_V][n3_4] = hitr.dU2();
- dv2_3[n3_V][n3_4] = hitr.dV2();
- z_Pos_3[n3_V][n3_4] = zr;
- t_3[n3_V][n3_4] = hitr.T();
- dt2_3[n3_V][n3_4] = hitr.dT2();
+
+ x_3[n3_V][n3_4] = hitr.X();
+ y_3[n3_V][n3_4] = hitr.Y();
+ z_3[n3_V][n3_4] = zr;
+ t_3[n3_V][n3_4] = hitr.T();
+ dx_3[n3_V][n3_4] = hitr.dX();
+ dxy_3[n3_V][n3_4] = hitr.dXY();
+ dy_3[n3_V][n3_4] = hitr.dY();
+ dt2_3[n3_V][n3_4] = hitr.dT2();
n3++;
n3_V = n3 / fvec::size();
@@ -786,12 +779,13 @@ inline void L1Algo::findTripletsStep0( // input
if (0 == n3_4) {
T_3.push_back(L1TrackPar_0);
- u_front_3.push_back(fvec::Zero());
- u_back_3.push_back(fvec::Zero());
- z_Pos_3.push_back(fvec::Zero());
- du2_3.push_back(fvec::One());
- dv2_3.push_back(fvec::One());
+ x_3.push_back(fvec::Zero());
+ y_3.push_back(fvec::Zero());
+ z_3.push_back(fvec::Zero());
t_3.push_back(fvec::Zero());
+ dx_3.push_back(fvec::One());
+ dxy_3.push_back(fvec::Zero());
+ dy_3.push_back(fvec::One());
dt2_3.push_back(fvec::One());
}
} // search area
@@ -802,10 +796,9 @@ inline void L1Algo::findTripletsStep0( // input
/// Add the right hits to parameters estimation.
inline void L1Algo::findTripletsStep1( // input
- Tindex n3_V, const L1Station& star, L1Vector<fvec>& u_front_, L1Vector<fvec>& u_back_, L1Vector<fvec>& z_Pos,
- L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3,
+ Tindex n3_V, const L1Station& star, L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3,
+ L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3, L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3,
// output
- // L1TrackPar *T_3
L1Vector<L1TrackPar>& T_3)
{
L1Fit fit;
@@ -813,9 +806,12 @@ inline void L1Algo::findTripletsStep1( // input
fit.SetMask(fmask::One());
for (Tindex i3_V = 0; i3_V < n3_V; ++i3_V) {
fit.SetTrack(T_3[i3_V]);
- fit.ExtrapolateLineNoField(z_Pos[i3_V]);
- fit.Filter(star.frontInfo, u_front_[i3_V], du2_3[i3_V]);
- fit.Filter(star.backInfo, u_back_[i3_V], dv2_3[i3_V]);
+ fit.ExtrapolateLineNoField(z_3[i3_V]);
+ L1XYMeasurementInfo info;
+ info.C00 = dx_3[i3_V] * dx_3[i3_V];
+ info.C10 = dxy_3[i3_V];
+ info.C11 = dy_3[i3_V] * dy_3[i3_V];
+ fit.FilterXY(info, x_3[i3_V], y_3[i3_V]);
if (kMcbm != fTrackingMode) { fit.FilterTime(t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
T_3[i3_V] = fit.Tr();
}
@@ -880,20 +876,20 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
if ((mc1 != mc2) || (mc1 != mc3)) { continue; }
}
- fscal u[NHits], v[NHits], t[NHits], x[NHits], y[NHits], z[NHits];
- fscal du2[NHits], dv2[NHits], dt2[NHits];
+ fscal x[NHits], y[NHits], z[NHits], t[NHits];
+ fscal dt2[NHits];
+ L1XYMeasurementInfo cov[NHits];
for (int ih = 0; ih < NHits; ++ih) {
- const L1Hit& hit = fInputData.GetHit(ihit[ih]);
- u[ih] = hit.u;
- v[ih] = hit.v;
- t[ih] = hit.t;
- std::tie(x[ih], y[ih]) = sta[ih].ConvUVtoXY<fscal>(u[ih], v[ih]);
- z[ih] = hit.z;
-
- du2[ih] = hit.du2;
- dv2[ih] = hit.dv2;
- dt2[ih] = hit.dt2;
+ const L1Hit& hit = fInputData.GetHit(ihit[ih]);
+ x[ih] = hit.x;
+ y[ih] = hit.y;
+ z[ih] = hit.z;
+ t[ih] = hit.t;
+ cov[ih].C00 = hit.dx * hit.dx;
+ cov[ih].C10 = hit.dxy;
+ cov[ih].C11 = hit.dy * hit.dy;
+ dt2[ih] = hit.dt2;
};
// find the field along the track
@@ -940,7 +936,9 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.t = t[ih0];
T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
- std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]);
+ T.C00 = cov[ih0].C00;
+ T.C10 = cov[ih0].C10;
+ T.C11 = cov[ih0].C11;
T.NDF = ndf;
@@ -952,8 +950,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
auto radThick = fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
fit.AddMsInMaterial(radThick);
fit.EnergyLossCorrection(radThick, fvec(-1.f));
- fit.Filter(sta[ih].frontInfo, u[ih], du2[ih]);
- fit.Filter(sta[ih].backInfo, v[ih], dv2[ih]);
+ fit.FilterXY(cov[ih], x[ih], y[ih]);
fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
}
}
@@ -976,7 +973,9 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.t = t[ih0];
T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
- std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]);
+ T.C00 = cov[ih0].C00;
+ T.C10 = cov[ih0].C10;
+ T.C11 = cov[ih0].C11;
T.NDF = ndf;
@@ -985,8 +984,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
auto radThick = fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
fit.AddMsInMaterial(radThick);
fit.EnergyLossCorrection(radThick, fvec(1.f));
- fit.Filter(sta[ih].frontInfo, u[ih], du2[ih]);
- fit.Filter(sta[ih].backInfo, v[ih], dv2[ih]);
+ fit.FilterXY(cov[ih], x[ih], y[ih]);
fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
}
}
@@ -1257,12 +1255,13 @@ inline void L1Algo::CreatePortionOfDoublets(
L1HitPoint* vHits_l = &(fGridPoints[0]) + fGridHitStartIndex[istal];
L1HitPoint* vHits_m = &(fGridPoints[0]) + fGridHitStartIndex[istam];
- fvec u_front[L1Constants::size::kSingletPortionSizeVec];
- fvec u_back[L1Constants::size::kSingletPortionSizeVec];
- fvec dv2[L1Constants::size::kSingletPortionSizeVec];
- fvec du2[L1Constants::size::kSingletPortionSizeVec];
- fvec zPos[L1Constants::size::kSingletPortionSizeVec];
+ fvec x[L1Constants::size::kSingletPortionSizeVec];
+ fvec y[L1Constants::size::kSingletPortionSizeVec];
+ fvec z[L1Constants::size::kSingletPortionSizeVec];
fvec t[L1Constants::size::kSingletPortionSizeVec];
+ fvec dx[L1Constants::size::kSingletPortionSizeVec];
+ fvec dxy[L1Constants::size::kSingletPortionSizeVec];
+ fvec dy[L1Constants::size::kSingletPortionSizeVec];
fvec dt2[L1Constants::size::kSingletPortionSizeVec];
/// prepare the portion of left hits data
@@ -1270,7 +1269,7 @@ inline void L1Algo::CreatePortionOfDoublets(
findSingletsStep0( // input
iSingletPortion * L1Constants::size::kSingletPortionSize, singletPortionSize, vHits_l,
// output
- u_front, u_back, zPos, hitsl_1, t, dt2, du2, dv2);
+ hitsl_1, x, y, z, t, dx, dxy, dy, dt2);
for (Tindex i = 0; i < singletPortionSize; ++i)
L1_ASSERT(hitsl_1[i] < fGridHitStopIndex[istal] - fGridHitStartIndex[istal],
@@ -1280,9 +1279,9 @@ inline void L1Algo::CreatePortionOfDoublets(
/// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station.
- findSingletsStep1(istal, istam, portionSize_V, u_front, u_back, zPos, t, dt2,
+ findSingletsStep1(istal, istam, portionSize_V, x, y, z, t, dx, dxy, dy, dt2,
// output
- T_1, fld_1, du2, dv2);
+ T_1, fld_1);
/// Find the doublets. Reformat data in the portion of doublets.
@@ -1373,24 +1372,30 @@ inline void L1Algo::CreatePortionOfTriplets(
L1Vector<L1HitIndex_t>& hitsl_3 = fTripletHitsL[Thread];
L1Vector<L1HitIndex_t>& hitsm_3 = fTripletHitsM[Thread];
L1Vector<L1HitIndex_t>& hitsr_3 = fTripletHitsR[Thread];
- L1Vector<fvec>& u_front3 = fTripletHitR_Ufront[Thread];
- L1Vector<fvec>& u_back3 = fTripletHitR_Uback[Thread];
- L1Vector<fvec>& z_pos3 = fTripletHitR_Z[Thread];
- L1Vector<fvec>& t_3 = fTripletHitR_Time[Thread];
- L1Vector<fvec>& dt2_3 = fTripletHitR_TimeErr[Thread];
- L1Vector<fvec>& du2_3 = fTripletHitR_dUfront[Thread];
- L1Vector<fvec>& dv2_3 = fTripletHitR_dUback[Thread];
+
+ L1Vector<fvec>& x_3 = fTripletHitR_X[Thread];
+ L1Vector<fvec>& y_3 = fTripletHitR_Y[Thread];
+ L1Vector<fvec>& z_3 = fTripletHitR_Z[Thread];
+ L1Vector<fvec>& t_3 = fTripletHitR_Time[Thread];
+
+ L1Vector<fvec>& dx_3 = fTripletHitR_dX[Thread];
+ L1Vector<fvec>& dxy_3 = fTripletHitR_dXY[Thread];
+ L1Vector<fvec>& dy_3 = fTripletHitR_dY[Thread];
+ L1Vector<fvec>& dt2_3 = fTripletHitR_TimeErr[Thread];
T_3.clear();
hitsl_3.clear();
hitsm_3.clear();
hitsr_3.clear();
- u_front3.clear();
- u_back3.clear();
- z_pos3.clear();
- du2_3.clear();
- dv2_3.clear();
+
+ x_3.clear();
+ y_3.clear();
+ z_3.clear();
t_3.clear();
+
+ dx_3.clear();
+ dxy_3.clear();
+ dy_3.clear();
dt2_3.clear();
/// Find the triplets(right hit). Reformat data in the portion of triplets.
@@ -1403,7 +1408,7 @@ inline void L1Algo::CreatePortionOfTriplets(
n_2, hitsm_2, i1_2,
// output
- n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, du2_3, dv2_3, t_3, dt2_3);
+ n3, T_3, hitsl_3, hitsm_3, hitsr_3, x_3, y_3, z_3, t_3, dx_3, dxy_3, dy_3, dt2_3);
for (Tindex i = 0; i < static_cast<Tindex>(hitsl_3.size()); ++i)
@@ -1418,7 +1423,7 @@ inline void L1Algo::CreatePortionOfTriplets(
Tindex n3_V = (n3 + fvec::size() - 1) / fvec::size();
findTripletsStep1( // input
- n3_V, star, u_front3, u_back3, z_pos3, du2_3, dv2_3, t_3, dt2_3,
+ n3_V, star, x_3, y_3, z_3, t_3, dx_3, dxy_3, dy_3, dt2_3,
// output
T_3);
@@ -1574,7 +1579,7 @@ void L1Algo::CaTrackFinderSlice()
L1HitIndex_t nGridHitsFilled = 0;
for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
- const L1Station& st = fParameters.GetStation(iS);
+
fMaxDx[iS] = 0.;
fMaxDy[iS] = 0.;
fMaxDt[iS] = 0.;
@@ -1591,16 +1596,13 @@ void L1Algo::CaTrackFinderSlice()
for (L1HitIndex_t ih = 0; ih < fSliceHitIds[iS].size(); ++ih) {
const L1Hit& h = fInputData.GetHit(fSliceHitIds[iS][ih]);
- auto [dxx, dxy, dyy] = st.FormXYCovarianceMatrix(h.du2, h.dv2);
- fscal dx = sqrt(dxx);
- fscal dy = sqrt(dyy);
fscal dt = sqrt(h.dt2);
- if (fMaxDx[iS] < dx) { fMaxDx[iS] = dx; }
- if (fMaxDy[iS] < dy) { fMaxDy[iS] = dy; }
+ if (fMaxDx[iS] < h.dx) { fMaxDx[iS] = h.dx; }
+ if (fMaxDy[iS] < h.dy) { fMaxDy[iS] = h.dy; }
if (fMaxDt[iS] < dt) { fMaxDt[iS] = dt; }
- auto [x, y] = GetHitCoorOnGrid(h, iS);
+ auto [x, y] = GetHitCoorOnGrid(h);
if (x < gridMinX) { gridMinX = x; }
if (x > gridMaxX) { gridMaxX = x; }
if (y < gridMinY) { gridMinY = y; }
@@ -1654,8 +1656,7 @@ void L1Algo::CaTrackFinderSlice()
#pragma omp parallel for schedule(dynamic, 5)
#endif
for (unsigned int ih = 0; ih < fGridHits.size(); ++ih) {
- const L1Hit& h = fGridHits[ih];
- CreateHitPoint(h, fGridPoints[ih]);
+ fGridPoints[ih].Set(fGridHits[ih]);
}
#ifdef XXX
@@ -2175,14 +2176,11 @@ void L1Algo::CaTrackFinderSlice()
fvHitKeyFlags[hit.f] = 1;
fvHitKeyFlags[hit.b] = 1;
fSliceRecoHits_thread[num_thread].push_back(*phIt);
- L1HitPoint tempPoint = CreateHitPoint(hit); //TODO take number of station from hit
-
- const L1Station& stah = fParameters.GetStation(hit.iSt);
- auto [xcoor, ycoor] = stah.ConvUVtoXY<fscal>(tempPoint.U(), tempPoint.V());
-
- fscal zcoor = tempPoint.Z() - fParameters.GetTargetPositionZ()[0];
+ fscal dx = hit.x - fParameters.GetTargetPositionX()[0];
+ fscal dy = hit.y - fParameters.GetTargetPositionY()[0];
+ fscal dz = hit.z - fParameters.GetTargetPositionZ()[0];
- fscal timeFlight = sqrt(xcoor * xcoor + ycoor * ycoor + zcoor * zcoor) / 30.f; // c = 30[cm/ns]
+ fscal timeFlight = sqrt(dx * dx + dy * dy + dz * dz) / 30.f; // c = 30[cm/ns]
sumTime += (hit.t - timeFlight);
}
@@ -2273,7 +2271,7 @@ void L1Algo::CaTrackFinderSlice()
int NHitsOnStationTmp = NHitsOnStation;
vGridTime[ista].UpdateIterGrid(Nelements, staHits, fGridHitIdsBuf, staHitIndices, fGridHitsBuf, fGridPointsBuf,
- staHitPoints, NHitsOnStation, ista, *this, fvHitKeyFlags);
+ staHitPoints, NHitsOnStation, *this, fvHitKeyFlags);
fGridHitStartIndex[ista] = NHitsOnStationTmp;
fGridHitStopIndex[ista] = NHitsOnStation;
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index 8fe46a0534673e32b0b472161869f44ca4546da6..b8ef95582e130e78f823074c72634b6f28ace13b 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -4,9 +4,12 @@
#include "L1Fit.h"
+#include "L1Hit.h"
+#include "L1Station.h"
+
#define cnst const fvec
-void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
+void L1Fit::FilterU(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
{
fvec zeta, HCH;
fvec F0, F1, F2, F3, F4, F5, F6;
@@ -35,6 +38,9 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
fvec wi = fMaskF / (sigma2 + fvec(1.0000001) * HCH);
fvec zetawi = fMaskF * zeta / (iif(maskDoFilter, sigma2, fvec::Zero()) + HCH);
+ //wi = iif(wi > fvec::Zero(), wi, fvec::Zero());
+ wi = iif(sigma2 > fvec::Zero(), wi, fvec::Zero());
+
// fTr.chi2 += iif( maskDoFilter, zeta * zetawi, fvec::Zero() );
fTr.chi2 += zeta * zeta * wi;
fTr.NDF += fMaskF;
@@ -90,7 +96,6 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
fTr.C66 -= K6 * F6;
}
-
void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
{
// filter track with a time measurement
@@ -180,6 +185,22 @@ void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
{
+ {
+ L1UMeasurementInfo xInfo;
+ xInfo.cos_phi = fvec::One();
+ xInfo.sin_phi = fvec::Zero();
+
+ L1UMeasurementInfo uInfo;
+ uInfo.cos_phi = -info.C10 / info.C00;
+ uInfo.sin_phi = fvec::One();
+ fvec u = uInfo.cos_phi * x + y;
+ fvec u2 = info.C11 - info.C10 * info.C10 / info.C00;
+
+ FilterU(xInfo, x, info.C00);
+ FilterU(uInfo, u, u2);
+ return;
+ }
+
cnst TWO(2.);
fvec zeta0, zeta1, S00, S10, S11, si;
@@ -286,6 +307,16 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
fTr.C66 -= K60 * F60 + K61 * F61;
}
+void L1Fit::FilterHit(const L1Station& sta, const L1Hit& hit)
+{
+ L1XYMeasurementInfo info;
+ info.C00 = hit.dx * hit.dx;
+ info.C10 = hit.dxy;
+ info.C11 = hit.dy * hit.dy;
+ FilterXY(info, hit.x, hit.y);
+ FilterTime(hit.t, hit.dt2, sta.timeInfo);
+}
+
void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY,
cnst Jx[6], cnst Jy[6])
{
@@ -1822,6 +1853,28 @@ void L1Fit::FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C
chi2 += zeta * zeta / (du2 + HCH);
}
+std::tuple<fvec, fvec> L1Fit::GetChi2XChi2U(fvec xm, fvec ym, fvec V00, fvec V10, fvec V11, fvec xt, fvec yt, fvec C00,
+ fvec C10, fvec C11)
+{
+
+ L1UMeasurementInfo xInfo;
+ xInfo.cos_phi = fvec::One();
+ xInfo.sin_phi = fvec::Zero();
+
+ L1UMeasurementInfo uInfo;
+ uInfo.cos_phi = -V10 / V00;
+ uInfo.sin_phi = fvec::One();
+ fvec u = uInfo.cos_phi * xm + ym;
+ fvec du2 = V11 - V10 * V10 / V00;
+
+ fvec chi2x = 0.f;
+ fvec chi2u = 0.f;
+
+ FilterChi2XYC00C10C11(xInfo, xt, yt, C00, C10, C11, chi2x, xm, V00);
+ FilterChi2(uInfo, xt, yt, C00, C10, C11, chi2u, u, du2);
+ return std::tuple<fvec, fvec>(chi2x, chi2u);
+}
+
void L1Fit::GuessTrack(const fvec& trackZ, const fvec hitX[], const fvec hitY[], const fvec hitZ[], const fvec hitT[],
const fvec By[], const fmask hitW[], const fmask hitWtime[], int NHits)
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index 2be29385f56a804b759cdca1fafb23192e90b09f..a3540cf85cd69b68386ee566406a579e10da003d 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -17,6 +17,9 @@
#include "L1UMeasurementInfo.h"
#include "L1XYMeasurementInfo.h"
+class L1Hit;
+class L1Station;
+
#define cnst const fvec
/// Kalman filter fitter for L1 tracking
@@ -89,13 +92,15 @@ public:
/// get the particle mass
fvec GetMaxExtrapolationStep() const { return fMaxExtraplationStep; }
-
- void Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2);
+ void FilterU(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2);
void FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y);
void FilterTime(cnst& t, cnst& dt2, cnst& timeInfo);
+ void FilterHit(const L1Station& s, const L1Hit& h);
+
+ void FilterVi(fvec vi);
+
void FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY, cnst Jx[6],
cnst Jy[6]);
- void FilterVi(fvec vi);
void MeasureVelocityWithQp();
@@ -140,6 +145,9 @@ public:
static void FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11, fvec& chi2,
cnst& u, cnst& du2);
+ static std::tuple<fvec, fvec> GetChi2XChi2U(fvec xm, fvec ym, fvec V00, fvec V10, fvec V11, fvec xt, fvec yt,
+ fvec C00, fvec C10, fvec C11);
+
void GuessTrack(const fvec& trackZ, const fvec hitX[], const fvec hitY[], const fvec hitZ[], const fvec hitT[],
const fvec By[], const fmask hitW[], const fmask hitWtime[], int NHits);
diff --git a/reco/L1/L1Algo/L1Grid.cxx b/reco/L1/L1Algo/L1Grid.cxx
index 67b2208f984d5ca2e5cd3951b1f45483334be97c..9bade04c188cdd21d31d456efa5173e547f81bbc 100644
--- a/reco/L1/L1Algo/L1Grid.cxx
+++ b/reco/L1/L1Algo/L1Grid.cxx
@@ -39,7 +39,7 @@ inline void memset(T* dest, T i, size_t num)
void L1Grid::UpdateIterGrid(unsigned int Nelements, L1Hit* hits, L1Vector<L1HitIndex_t>& indicesBuf,
L1HitIndex_t* indices, L1Vector<L1Hit>& hitsBuf, L1Vector<L1HitPoint>& pointsBuf,
- L1HitPoint* points, int& NHitsOnStation, char iS, L1Algo& Algo,
+ L1HitPoint* points, int& NHitsOnStation, L1Algo& Algo,
const L1Vector<unsigned char>& vSFlag)
{
//L1_SHOW(vSFlag.size());
@@ -56,7 +56,7 @@ void L1Grid::UpdateIterGrid(unsigned int Nelements, L1Hit* hits, L1Vector<L1HitI
const L1Hit& hit = hits[x];
if (!(vSFlag[hit.f] || vSFlag[hit.b])) {
- std::tie(xs, ys) = Algo.GetHitCoorOnGrid(hit, iS);
+ std::tie(xs, ys) = Algo.GetHitCoorOnGrid(hit);
const L1HitIndex_t& bin = GetBinBounded(xs, ys, hit.t);
fHitsInBin[x] = fFirstHitInBin[bin + 1];
@@ -97,7 +97,7 @@ void L1Grid::UpdateIterGrid(unsigned int Nelements, L1Hit* hits, L1Vector<L1HitI
const L1Hit& hit = hits[x];
if (!(vSFlag[hit.f] || vSFlag[hit.b])) {
- std::tie(xs, ys) = Algo.GetHitCoorOnGrid(hit, iS);
+ std::tie(xs, ys) = Algo.GetHitCoorOnGrid(hit);
const L1HitIndex_t& bin = GetBinBounded(xs, ys, hit.t);
@@ -179,10 +179,9 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, L1HitIndex_t& nGridHitsFilled)
for (L1HitIndex_t ih = 0; ih < nHits; ih++) {
L1HitIndex_t caHitId = algo.fSliceHitIds[iS][ih];
const L1Hit& h = algo.GetInputData().GetHit(caHitId);
- fscal x, y;
- std::tie(x, y) = algo.GetHitCoorOnGrid(h, iS);
- auto bin = GetBinBounded(x, y, h.t);
- fHitsInBin[ih] = fFirstHitInBin[bin + 1];
+ auto [x, y] = algo.GetHitCoorOnGrid(h);
+ auto bin = GetBinBounded(x, y, h.t);
+ fHitsInBin[ih] = fFirstHitInBin[bin + 1];
#ifdef _OPENMP
#pragma omp atomic
#endif
@@ -218,9 +217,8 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, L1HitIndex_t& nGridHitsFilled)
for (L1HitIndex_t ih = 0; ih < nHits; ih++) {
L1HitIndex_t caHitId = algo.fSliceHitIds[iS][ih];
const L1Hit& h = algo.GetInputData().GetHit(caHitId);
- fscal x, y;
- std::tie(x, y) = algo.GetHitCoorOnGrid(h, iS);
- auto bin = GetBinBounded(x, y, h.t);
+ auto [x, y] = algo.GetHitCoorOnGrid(h);
+ auto bin = GetBinBounded(x, y, h.t);
{
const L1HitIndex_t& index1 = fFirstHitInBin[bin] + fHitsInBin[ih];
algo.fGridHits[nGridHitsFilled + index1] = h;
diff --git a/reco/L1/L1Algo/L1Grid.h b/reco/L1/L1Algo/L1Grid.h
index f6a8653fa0bad11d4600e3507f09f15667b6dd0d..9945ab3441321a36d44fc09a7a18756f3759ba52 100644
--- a/reco/L1/L1Algo/L1Grid.h
+++ b/reco/L1/L1Algo/L1Grid.h
@@ -98,7 +98,7 @@ public:
void UpdateIterGrid(unsigned int Nelements, L1Hit* hits, L1Vector<L1HitIndex_t>& indicesBuf, L1HitIndex_t* indices,
L1Vector<L1Hit>& hitsBuf, L1Vector<L1HitPoint>& pointsBuf, L1HitPoint* points,
- int& NHitsOnStation, char iS, L1Algo& Algo, const L1Vector<unsigned char>& vSFlag);
+ int& NHitsOnStation, L1Algo& Algo, const L1Vector<unsigned char>& vSFlag);
private:
diff --git a/reco/L1/L1Algo/L1Hit.h b/reco/L1/L1Algo/L1Hit.h
index 983b0462a980e44b317772f5de7e7b1b32072ef4..b3c375a84ad3717f9aec1d9323dd07d6e7483400 100644
--- a/reco/L1/L1Algo/L1Hit.h
+++ b/reco/L1/L1Algo/L1Hit.h
@@ -34,16 +34,16 @@ public:
/// tracking detectors (MVD, MuCh, TRD, TOF) f == b and corresponds to the index of the hit. Indexes f and b
/// do not intersect between different detector stations.
- float u = 0.f; ///< measured U coordinate [cm]
- float v = 0.f; ///< measured V coordinate [cm]
- float t = 0.f; ///< measured time [ns]
- float z = 0.f; ///< fixed Z coordinate [cm]
- float du2 = 0.f; ///< measured uncertainty of U coordinate [cm]
- float dv2 = 0.f; ///< measured uncertainty of V coordinate [cm]
+ float x {0.f}; ///< measured X coordinate [cm]
+ float y {0.f}; ///< measured Y coordinate [cm]
+ float z = 0.f; ///< fixed Z coordinate [cm]
+ float t = 0.f; ///< measured time [ns]
+ float dx {0.f}; ///< uncertainty of X coordinate [cm]
+ float dy {0.f}; ///< uncertainty of Y coordinate [cm]
+ float dxy {0.f}; ///< X/Y covariance [cm2]
float dt2 = 0.f; ///< measured uncertainty of time [ns]
- int ID = 0; ///< index of hit before hits sorting
- int iSt = 0; ///< index of station in the active stations array
- // TODO: Test speed penalty of using iSt index
+ int ID = 0; ///< index of hit before hit sorting
+ int iSt = 0; ///< index of station in the active stations array
private:
/// Serialization method, used to save L1Hit objects into binary or text file in a defined order
@@ -52,12 +52,13 @@ private:
{
ar& f;
ar& b;
- ar& u;
- ar& v;
- ar& t;
+ ar& x;
+ ar& y;
ar& z;
- ar& du2;
- ar& dv2;
+ ar& t;
+ ar& dx;
+ ar& dxy;
+ ar& dy;
ar& dt2;
ar& ID;
ar& iSt;
diff --git a/reco/L1/L1Algo/L1HitPoint.h b/reco/L1/L1Algo/L1HitPoint.h
index a7a367c7e6a76fe9fd8fdd4d474d8f17820f8663..5ded6ff21997b52ed6c606b8f505f50ebbdb3de7 100644
--- a/reco/L1/L1Algo/L1HitPoint.h
+++ b/reco/L1/L1Algo/L1HitPoint.h
@@ -5,55 +5,47 @@
#ifndef _L1HitPoint_h_
#define _L1HitPoint_h_
+#include "L1Hit.h"
+
/// contain strips positions and coordinates of the hit
struct L1HitPoint {
L1HitPoint() = default;
- L1HitPoint(fscal z_, fscal v_, fscal u_, fscal t_, fscal du2_, fscal dv2_, fscal dt2_)
- : z(z_)
- , u(u_)
- , v(v_)
- , t(t_)
- , du2(du2_)
- , dv2(dv2_)
- , dt2(dt2_) {};
+ L1HitPoint(const L1Hit& hit) { Set(hit); };
+
+ void Set(const L1Hit& hit)
+ {
+ x = hit.x;
+ y = hit.y;
+ z = hit.z;
+ t = hit.t;
+ dx = hit.dx;
+ dxy = hit.dxy;
+ dy = hit.dy;
+ dt2 = hit.dt2;
+
+ fIsSuppressed = 0;
+ }
fscal Z() const { return z; }
- fscal U() const { return u; }
- fscal V() const { return v; }
fscal T() const { return t; }
-
- fscal dU2() const { return du2; }
- fscal dV2() const { return dv2; }
+ fscal X() const { return x; }
+ fscal Y() const { return y; }
+ fscal dX() const { return dx; }
+ fscal dXY() const { return dxy; }
+ fscal dY() const { return dy; }
+ fscal dX2() const { return dx * dx; }
+ fscal dY2() const { return dy * dy; }
fscal dT2() const { return dt2; }
bool IsSuppressed() const { return fIsSuppressed; }
- void SetZ(fscal z_) { z = z_; }
- void SetU(fscal u_) { u = u_; }
- void SetV(fscal v_) { v = v_; }
- void SetT(fscal t_) { t = t_; }
-
- void Set(const fscal& z_, const fscal& u_, const fscal& v_, const fscal& t_, const fscal& du2_, const fscal& dv2_,
- const fscal& dt2_)
- {
- z = z_;
- u = u_;
- v = v_;
- t = t_;
- du2 = du2_;
- dv2 = dv2_;
- dt2 = dt2_;
-
- fIsSuppressed = 0;
- }
-
void SetIsSuppresed(bool val) { fIsSuppressed = val; }
private:
- fscal z {0.}, u {0.}, v {0.}, t {0.}, du2 {0.}, dv2 {0.}, dt2 {0.};
+ fscal x {0.}, y {0.}, z {0.}, t {0.}, dx {0.}, dxy {0.}, dy {0.}, dt2 {0.};
// fIsSuppressed flag is used to suppress duplicated hits at the module overlaps
// TODO: collect those hits on the track instead of suppressing them
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 4629b5090af0ef52fdea1afa7ec5a713c4362ddf..e4ad499dcc3c4c67c1045a4f6ebf55903c4ec48d 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -149,16 +149,16 @@ public:
const InitController_t& GetInitController() const { return fInitController; }
/// Gets total number of active stations
- [[deprecated]] int GetNstationsActive() const;
+ int GetNstationsActive() const;
/// Gets number of active stations for given detector ID
- [[deprecated]] int GetNstationsActive(L1DetectorID detectorID) const;
+ int GetNstationsActive(L1DetectorID detectorID) const;
/// Gets total number of stations, provided by setup geometry
- [[deprecated]] int GetNstationsGeometry() const;
+ int GetNstationsGeometry() const;
/// Gets number of stations, provided by setup geometry for given detector ID
- [[deprecated]] int GetNstationsGeometry(L1DetectorID detectorID) const;
+ int GetNstationsGeometry(L1DetectorID detectorID) const;
/// Gets a name of the output configuration file
const std::string& GetOutputConfigName() const { return fConfigOutputName; }
diff --git a/reco/L1/L1Algo/L1Station.cxx b/reco/L1/L1Algo/L1Station.cxx
index 85be4deca5a28ede207faa67c9e20ae9f1faf978..2a8594642feae07453ca8632191a75301ec5bfc5 100644
--- a/reco/L1/L1Algo/L1Station.cxx
+++ b/reco/L1/L1Algo/L1Station.cxx
@@ -67,26 +67,6 @@ void L1Station::CheckConsistency() const
// TODO: Temporarily switched off, because Much has RL = 0, which is actually incorrect, but really is not used.
// One should provide average radiation length values for each Much layer (S.Zharko)
fieldSlice.CheckConsistency();
- frontInfo.CheckConsistency();
- backInfo.CheckConsistency();
- xInfo.CheckConsistency();
- yInfo.CheckConsistency();
-
- // Check consistency of coordinate transformations
- for (float x = -1.f; x < 1.1f; x += 0.2) {
- for (float y = -1.f; y < 1.1f; y += 0.2) {
- auto [u, v] = ConvXYtoUV<float>(x, y);
- auto [x1, y1] = ConvUVtoXY<float>(u, v);
-
- if (fabs(x1 - x) > 1.e-6 || fabs(y1 - y) > 1.e-6) {
- std::stringstream msg;
- msg << "L1Station: " << this->ToString() << " makes wrong coordinate convertion: "
- << "x,y " << x << "," << y << " -> u,v " << u << "," << v << " -> x,y " << x1 << "," << y1
- << " are not the same: dx,dy " << x1 - x << "," << y1 - y;
- throw std::logic_error(msg.str());
- }
- }
- }
}
// TODO: Improve log style (S.Zh.)
@@ -116,17 +96,6 @@ std::string L1Station::ToString(int verbosityLevel, int indentLevel) const
aStream << indent << "Field approcimation coefficients:\n";
aStream << fieldSlice.ToString(indentLevel + 1) << '\n';
}
- if (verbosityLevel > 2) {
- aStream << indent << "Strips geometry:\n";
- aStream << indent << indentChar << "Front:\n";
- aStream << frontInfo.ToString(indentLevel + 2) << '\n';
- aStream << indent << indentChar << "Back:\n";
- aStream << backInfo.ToString(indentLevel + 2) << '\n';
- aStream << indent << indentChar << "X layer:\n";
- aStream << xInfo.ToString(indentLevel + 2) << '\n';
- aStream << indent << indentChar << "Y layer:\n";
- aStream << yInfo.ToString(indentLevel + 2) << '\n';
- }
}
return aStream.str();
}
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index 0993b435e35446879cf700ecbedea640e8063947..110394c535fbc8f4317734331c97a3f4edf7691b 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -33,10 +33,6 @@ public:
fvec Rmax = undef::kFvc; ///< max radius of the station [cm]
L1FieldSlice fieldSlice {};
- L1UMeasurementInfo frontInfo {};
- L1UMeasurementInfo backInfo {};
- L1UMeasurementInfo xInfo {}; ///< x axis in front,back coordinates
- L1UMeasurementInfo yInfo {}; ///< y axis in front,back coordinates
/// Serialization function
friend class boost::serialization::access;
@@ -53,10 +49,6 @@ public:
ar& Rmax;
ar& fieldSlice;
- ar& frontInfo;
- ar& backInfo;
- ar& xInfo;
- ar& yInfo;
}
@@ -74,78 +66,7 @@ public:
/// initialization
void CheckConsistency() const;
- /// Converts Cartesian coordinates X and Y into strip coordinates U and V in transverse plane
- /// \param x X coordinate [length unit]
- /// \param y Y coordinate [length unit]
- /// \return Pair of strip coordinates [u, v]
- template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool> = true>
- std::pair<T, T> ConvXYtoUV(T x, T y) const;
-
- /// Converts strip coordinates U and V into Cartesian coordinates X and Y in transverse plane
- /// \param u U coordinate [length unit]
- /// \param v V coordinate [length unit]
- /// \return Pair of Cartesian coordinates [x, y]
- template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool> = true>
- std::pair<T, T> ConvUVtoXY(T u, T v) const;
-
- /// Converts variances of strip coordinates U and V into covariance matrix of Cartesian coordinates X and Y
- /// \param du U coordinate uncertainty [length unit]
- /// \param dv V coordinate uncertainty [length unit]
- /// \return Covariance matrix of hit position in Cartesian coordinates: [dxx, dxy, dyy] [length unit squared]
- template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool> = true>
- std::tuple<T, T, T> FormXYCovarianceMatrix(T du2, T dv2) const;
-
} _fvecalignment;
-//
-// Template member functions implementation
-//
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool>>
-std::pair<T, T> L1Station::ConvXYtoUV(T x, T y) const
-{
- if constexpr (std::is_same<T, fvec>::value) {
- return std::make_pair(x * frontInfo.cos_phi + y * frontInfo.sin_phi, x * backInfo.cos_phi + y * backInfo.sin_phi);
- }
- else {
- return std::make_pair(x * frontInfo.cos_phi[0] + y * frontInfo.sin_phi[0],
- x * backInfo.cos_phi[0] + y * backInfo.sin_phi[0]);
- }
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool>>
-std::pair<T, T> L1Station::ConvUVtoXY(T u, T v) const
-{
- if constexpr (std::is_same<T, fvec>::value) {
- return std::make_pair(u * xInfo.cos_phi + v * xInfo.sin_phi, u * yInfo.cos_phi + v * yInfo.sin_phi);
- }
- else {
- return std::make_pair(u * xInfo.cos_phi[0] + v * xInfo.sin_phi[0], u * yInfo.cos_phi[0] + v * yInfo.sin_phi[0]);
- }
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool>>
-std::tuple<T, T, T> L1Station::FormXYCovarianceMatrix(T du2, T dv2) const
-{
- if constexpr (std::is_same<T, fvec>::value) {
- return std::make_tuple(xInfo.cos_phi * xInfo.cos_phi * du2 + xInfo.sin_phi * xInfo.sin_phi * dv2, // dx2
- xInfo.cos_phi * yInfo.cos_phi * du2 + xInfo.sin_phi * yInfo.sin_phi * dv2, // dxy
- yInfo.cos_phi * yInfo.cos_phi * du2 + yInfo.sin_phi * yInfo.sin_phi * dv2 // dy2
- );
- }
- else {
- return std::make_tuple(
- xInfo.cos_phi[0] * xInfo.cos_phi[0] * du2 + xInfo.sin_phi[0] * xInfo.sin_phi[0] * dv2, // dx2
- xInfo.cos_phi[0] * yInfo.cos_phi[0] * du2 + xInfo.sin_phi[0] * yInfo.sin_phi[0] * dv2, // dxy
- yInfo.cos_phi[0] * yInfo.cos_phi[0] * du2 + yInfo.sin_phi[0] * yInfo.sin_phi[0] * dv2 // dy2
- );
- }
-}
#endif // L1Station_h
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 1e51f2762e6aab01e7dc34861c29e03cfeef57f9..d5058c45a8fb61a1313ec1daff5649be302f0ae7 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -48,17 +48,10 @@ void L1Algo::L1KFTrackFitter()
const L1Station* sta = fParameters.GetStations().begin();
// Spatial-time position of a hit vs. station and track in the portion
- // NOTE: u- and v-axes are axes, orthogonal to front and back strips of the station, respectively.
- fvec u[L1Constants::size::kMaxNstations]; // Hit position along the u-axis [cm]
- fvec v[L1Constants::size::kMaxNstations]; // Hit position along the v-axis [cm]
- fvec du2[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the u-axis [cm]
- fvec dv2[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the v-axis [cm]
fvec x[L1Constants::size::kMaxNstations]; // Hit position along the x-axis [cm]
fvec y[L1Constants::size::kMaxNstations]; // Hit position along the y-axis [cm]
- fvec d_xx[L1Constants::size::kMaxNstations]; // Variance of the x hit position coordinate [cm2]
- fvec d_yy[L1Constants::size::kMaxNstations]; // Variance of the y hit position coordinate [cm2]
- fvec d_xy[L1Constants::size::kMaxNstations]; // Covariance between the x and y hit position coordinates [cm2]
+ L1XYMeasurementInfo cov_xy[L1Constants::size::kMaxNstations]; // Covariance matrix for x,y
fvec z[L1Constants::size::kMaxNstations]; // Hit position along the z-axis (precised) [cm]
@@ -67,9 +60,7 @@ void L1Algo::L1KFTrackFitter()
fvec x_first;
fvec y_first;
- fvec d_xx_fst;
- fvec d_yy_fst;
- fvec d_xy_fst;
+ L1XYMeasurementInfo cov_xy_fst;
fvec time_first;
fvec wtime_first;
@@ -77,9 +68,7 @@ void L1Algo::L1KFTrackFitter()
fvec x_last;
fvec y_last;
- fvec d_xx_lst;
- fvec d_yy_lst;
- fvec d_xy_lst;
+ L1XYMeasurementInfo cov_xy_lst;
fvec time_last;
fvec wtime_last;
@@ -103,7 +92,10 @@ void L1Algo::L1KFTrackFitter()
fvec ZSta[L1Constants::size::kMaxNstations];
for (int ista = 0; ista < nStations; ista++) {
- ZSta[ista] = sta[ista].fZ;
+ ZSta[ista] = sta[ista].fZ;
+ cov_xy[ista].C00 = 1.f;
+ cov_xy[ista].C10 = 0.f;
+ cov_xy[ista].C11 = 1.f;
}
unsigned short N_vTracks = fSliceRecoTracks.size(); // number of tracks processed per one SSE register
@@ -146,45 +138,41 @@ void L1Algo::L1KFTrackFitter()
w[ista][iVec] = true;
if (sta[ista].timeInfo) { w_time[ista][iVec] = true; }
- u[ista][iVec] = hit.u;
- v[ista][iVec] = hit.v;
- du2[ista][iVec] = hit.du2;
- dv2[ista][iVec] = hit.dv2;
- std::tie(x_temp, y_temp) = sta[ista].ConvUVtoXY<fvec>(u[ista], v[ista]);
- x[ista][iVec] = x_temp[iVec];
- y[ista][iVec] = y_temp[iVec];
+ x[ista][iVec] = hit.x; //x_temp[iVec];
+ y[ista][iVec] = hit.y; //y_temp[iVec];
time[ista][iVec] = hit.t;
dt2[ista][iVec] = hit.dt2;
if (!sta[ista].timeInfo) { dt2[ista][iVec] = 1.e4; }
z[ista][iVec] = hit.z;
sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
- std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(du2[ista], dv2[ista]);
+ cov_xy[ista].C00[iVec] = hit.dx * hit.dx;
+ cov_xy[ista].C10[iVec] = hit.dxy;
+ cov_xy[ista].C11[iVec] = hit.dy * hit.dy;
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
if (ih == 0) {
- z_start[iVec] = z[ista][iVec];
- x_first[iVec] = x[ista][iVec];
- y_first[iVec] = y[ista][iVec];
- time_first[iVec] = time[ista][iVec];
- wtime_first[iVec] = sta[ista].timeInfo ? 1. : 0.;
- dt2_first[iVec] = dt2[ista][iVec];
- d_xx_fst[iVec] = d_xx[ista][iVec];
- d_yy_fst[iVec] = d_yy[ista][iVec];
- d_xy_fst[iVec] = d_xy[ista][iVec];
+ z_start[iVec] = z[ista][iVec];
+ x_first[iVec] = x[ista][iVec];
+ y_first[iVec] = y[ista][iVec];
+ time_first[iVec] = time[ista][iVec];
+ wtime_first[iVec] = sta[ista].timeInfo ? 1. : 0.;
+ dt2_first[iVec] = dt2[ista][iVec];
+ cov_xy_fst.C00[iVec] = cov_xy[ista].C00[iVec];
+ cov_xy_fst.C10[iVec] = cov_xy[ista].C10[iVec];
+ cov_xy_fst.C11[iVec] = cov_xy[ista].C11[iVec];
}
-
else if (ih == nHitsTrack - 1) {
- z_end[iVec] = z[ista][iVec];
- x_last[iVec] = x[ista][iVec];
- y_last[iVec] = y[ista][iVec];
- d_xx_lst[iVec] = d_xx[ista][iVec];
- d_yy_lst[iVec] = d_yy[ista][iVec];
- d_xy_lst[iVec] = d_xy[ista][iVec];
- time_last[iVec] = time[ista][iVec];
- dt2_last[iVec] = dt2[ista][iVec];
- wtime_last[iVec] = sta[ista].timeInfo ? 1. : 0.;
+ z_end[iVec] = z[ista][iVec];
+ x_last[iVec] = x[ista][iVec];
+ y_last[iVec] = y[ista][iVec];
+ cov_xy_lst.C00[iVec] = cov_xy[ista].C00[iVec];
+ cov_xy_lst.C10[iVec] = cov_xy[ista].C10[iVec];
+ cov_xy_lst.C11[iVec] = cov_xy[ista].C11[iVec];
+ time_last[iVec] = time[ista][iVec];
+ dt2_last[iVec] = dt2[ista][iVec];
+ wtime_last[iVec] = sta[ista].timeInfo ? 1. : 0.;
}
}
@@ -212,8 +200,8 @@ void L1Algo::L1KFTrackFitter()
time_last = iif(w_time[ista], time_last, fvec::Zero());
dt2_last = iif(w_time[ista], dt2_last, fvec(1.e6));
- tr.ResetErrors(d_xx_lst, d_yy_lst, 0.1, 0.1, 1.0, dt2_last, 1.e-2);
- tr.C10 = d_xy_lst;
+ tr.ResetErrors(cov_xy_lst.C00, cov_xy_lst.C11, 0.1, 0.1, 1.0, dt2_last, 1.e-2);
+ tr.C10 = cov_xy_lst.C10;
tr.x = x_last;
tr.y = y_last;
tr.t = time_last;
@@ -253,8 +241,7 @@ void L1Algo::L1KFTrackFitter()
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), fvec(1.f));
fit.SetMask(initialised && w[ista]);
- fit.Filter(sta[ista].frontInfo, u[ista], du2[ista]);
- fit.Filter(sta[ista].backInfo, v[ista], dv2[ista]);
+ fit.FilterXY(cov_xy[ista], x[ista], y[ista]);
fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
@@ -270,13 +257,10 @@ void L1Algo::L1KFTrackFitter()
{
fitpv.SetMask(fmask::One());
- L1UMeasurementInfo vtxInfoX;
- vtxInfoX.cos_phi = 1.;
- vtxInfoX.sin_phi = 0.;
-
- L1UMeasurementInfo vtxInfoY;
- vtxInfoY.cos_phi = 0;
- vtxInfoY.sin_phi = 1.;
+ L1XYMeasurementInfo vtxInfo;
+ vtxInfo.C00 = fvec(1.e-8);
+ vtxInfo.C10 = fvec::Zero();
+ vtxInfo.C11 = fvec(1.e-8);
if (kGlobal == fTrackingMode) {
for (int vtxIter = 0; vtxIter < 2; vtxIter++) {
@@ -284,8 +268,7 @@ void L1Algo::L1KFTrackFitter()
fitpv.Tr() = fit.Tr();
fitpv.Tr().qp = fitpv.Qp0();
fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
- fitpv.Filter(vtxInfoX, fParameters.GetTargetPositionX(), fvec(1.e-8));
- fitpv.Filter(vtxInfoY, fParameters.GetTargetPositionY(), fvec(1.e-8));
+ fitpv.FilterXY(vtxInfo, fParameters.GetTargetPositionX(), fParameters.GetTargetPositionY());
}
}
else {
@@ -318,8 +301,8 @@ void L1Algo::L1KFTrackFitter()
ista = 0;
- tr.ResetErrors(d_xx_fst, d_yy_fst, 0.1, 0.1, 1., dt2_first, 1.e-2);
- tr.C10 = d_xy_fst;
+ tr.ResetErrors(cov_xy_fst.C00, cov_xy_fst.C11, 0.1, 0.1, 1., dt2_first, 1.e-2);
+ tr.C10 = cov_xy_fst.C10;
tr.x = x_first;
tr.y = y_first;
@@ -357,8 +340,7 @@ void L1Algo::L1KFTrackFitter()
fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y));
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), fvec(-1.f));
fit.SetMask(initialised && w[ista]);
- fit.Filter(sta[ista].frontInfo, u[ista], du2[ista]);
- fit.Filter(sta[ista].backInfo, v[ista], dv2[ista]);
+ fit.FilterXY(cov_xy[ista], x[ista], y[ista]);
fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
fldB2 = fldB1;
diff --git a/reco/L1/L1Algo/utils/CaUvConverter.cxx b/reco/L1/L1Algo/utils/CaUvConverter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4eac026d647c7c0b4ff98eaf216898da46590b5f
--- /dev/null
+++ b/reco/L1/L1Algo/utils/CaUvConverter.cxx
@@ -0,0 +1,68 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+#include "CaUvConverter.h"
+
+#include <Logger.h>
+
+#include <iostream>
+
+#include <cmath>
+
+using namespace cbm::ca;
+
+//----------------------------------------------------------------------------------------------------------------------
+//
+void CaUvConverter::SetFromUV(double phiU, double phiV)
+{
+ fcosU = cos(phiU);
+ fsinU = sin(phiU);
+ fcosV = cos(phiV);
+ fsinV = sin(phiV);
+
+ double det = fcosU * fsinV - fsinU * fcosV;
+ if (fabs(det) < 1.e-2) {
+ LOG(error) << "CaUvConverter: U & V coordinates are too close: " << phiU << " " << phiV;
+ phiV = phiU + 10. / 180. * M_PI;
+ fcosV = cos(phiV);
+ fsinV = sin(phiV);
+ det = fcosU * fsinV - fsinU * fcosV;
+ }
+ fcosX = fsinV / det;
+ fsinX = -fsinU / det;
+ fcosY = -fcosV / det;
+ fsinY = fcosU / det;
+}
+
+
+//----------------------------------------------------------------------------------------------------------------------
+//
+void CaUvConverter::SetFromXYCovMatrix(double phiU, double dx2, double dxy, double dy2)
+{
+ // take U coordinate from fPhiU
+ double cosU = cos(phiU);
+ double sinU = sin(phiU);
+ double du2 = cosU * cosU * dx2 + 2 * sinU * cosU * dxy + sinU * sinU * dy2;
+
+ // take V coordinate from the hit covariance matrix, making V uncorrelated with U
+
+ // rotate X,Y coordinates on angle fPhiU (X,Y) -> (U,W)
+
+ //double du2 = cosU * cosU * dx2 + 2. * sinU * cosU * dxy + sinU * sinU * dy2;
+ double duw = sinU * cosU * (dy2 - dx2) + (cosU * cosU - sinU * sinU) * dxy;
+ //double dw2 = cosU * cosU * dy2 - 2. * sinU * cosU * dxy + sinU * sinU * dx2;
+
+ // find an angle in rotaded coordinates, that has 0 covariance with U
+ double phiV = phiU + atan2(du2, -duw);
+
+ SetFromUV(phiU, phiV);
+
+ // check that u/v covariance is 0.
+ auto [tmp1, duv, dv2] = ConvertCovMatrixXYtoUV(dx2, dxy, dy2);
+
+ if (fabs(duv) > 1.e-8) { LOG(error) << "can not define V coordinate from XY covariance matrix"; }
+
+ //std::cout << " u " << phiU / M_PI * 180. << "v " << phiV / M_PI * 180. << " du " << sqrt(du2) << " dv " << sqrt(dv2)
+ // << " duv " << duv << std::endl;
+}
diff --git a/reco/L1/L1Algo/utils/CaUvConverter.h b/reco/L1/L1Algo/utils/CaUvConverter.h
new file mode 100644
index 0000000000000000000000000000000000000000..442a60b98734f4f017dcb2c0598172ce58a53313
--- /dev/null
+++ b/reco/L1/L1Algo/utils/CaUvConverter.h
@@ -0,0 +1,91 @@
+/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer], Igor Kulakov, Sergei Zharko */
+
+#ifndef CaUvConverter_h
+#define CaUvConverter_h 1
+
+
+#include "L1Def.h"
+#include "L1Undef.h"
+
+namespace cbm::ca
+{
+ /// @brief A class to convert XY coordinates to UV coordinates
+
+ class CaUvConverter {
+
+ public:
+ /// @brief construct from U,V angles
+ CaUvConverter(double phiU, double phiV) { SetFromUV(phiU, phiV); }
+
+ /// @brief construct from U angle and XY covariance matrix
+ CaUvConverter(double phiU, double dx2, double dxy, double dy2) { SetFromXYCovMatrix(phiU, dx2, dxy, dy2); }
+
+ /// @brief construct from U,V angles
+ void SetFromUV(double phiU, double phiV);
+
+ /// @brief construct from U angle and XY covariance matrix
+ void SetFromXYCovMatrix(double phiU, double dx2, double dxy, double dy2);
+
+ /// @brief Conversion function (x,y) -> (u,v)
+ /// @param x X-coordinate
+ /// @param y Y-coordinate
+ /// @return pair [u, v] coordinates
+ std::pair<double, double> ConvertXYtoUV(double x, double y) const
+ {
+ return std::make_pair(fcosU * x + fsinU * y, fcosV * x + fsinV * y);
+ }
+
+ /// @brief Conversion function (x,y) -> (u,v)
+ /// @param u U-coordinate
+ /// @param v V-coordinate
+ /// @return pair [x, y] coordinates
+ std::pair<double, double> ConvertUVtoXY(double u, double v) const
+ {
+ return std::make_pair(fcosX * u + fsinX * v, fcosY * u + fsinY * v);
+ }
+
+ /// @brief Conversion function (du2, duv, dv2) -> (dx2, dxy, dy2)
+ /// @param du2 Variance of U-coordinate measurement
+ /// @param duv Covariance of U & V - coordinate measurement
+ /// @param dv2 Variance of V-coordinate measurement
+ /// @return tuple [dx2, dxy, dy2]
+ std::tuple<double, double, double> ConvertCovMatrixUVtoXY(double du2, double duv, double dv2) const
+ {
+ return std::make_tuple(fcosX * fcosX * du2 + 2. * fsinX * fcosX * duv + fsinX * fsinX * dv2,
+ fcosX * fcosY * du2 + (fcosX * fsinY + fsinX * fcosY) * duv + fsinX * fsinY * dv2,
+ fcosY * fcosY * du2 + 2. * fsinY * fcosY * duv + fsinY * fsinY * dv2);
+ }
+
+ /// @brief Conversion function (dx2, dxy, dy2) -> (du2, duv, dv2)
+ /// @param dx2 Variance of X-coordinate measurement
+ /// @param dxy Covariance of X & Y - coordinate measurement
+ /// @param dy2 Variance of Y-coordinate measurement
+ /// @return tuple [du2, duv, dv2]
+ std::tuple<double, double, double> ConvertCovMatrixXYtoUV(double dx2, double dxy, double dy2) const
+ {
+ return std::make_tuple(fcosU * fcosU * dx2 + 2. * fsinU * fcosU * dxy + fsinU * fsinU * dy2,
+ fcosU * fcosV * dx2 + (fcosU * fsinV + fsinU * fcosV) * dxy + fsinU * fsinV * dy2,
+ fcosV * fcosV * dx2 + 2. * fsinV * fcosV * dxy + fsinV * fsinV * dy2);
+ }
+
+
+ private:
+ double fcosU {undef::kD}; ///< U coordinate in XY
+ double fsinU {undef::kD};
+
+ double fcosV {undef::kD}; ///< V coordinate in XY
+ double fsinV {undef::kD};
+
+ double fcosX {undef::kD}; ///< X coordinate in UV
+ double fsinX {undef::kD};
+
+ double fcosY {undef::kD}; ///< Y coordinate in UV
+ double fsinY {undef::kD};
+ };
+
+} // namespace cbm::ca
+
+
+#endif
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index 5d4098beaf76faff879f3c95466e7385677475af..f24ec91ca557f931b72fb0bc684aae120e29d3c3 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -561,9 +561,8 @@ void L1AlgoDraw::DrawInputHits()
int iMC = CbmL1::Instance()->fvHitPointIndexes[ih];
//if( (vSFlag[h.f] | vSFlagB[h.b] )&0x02 ) continue; // if used
- fscal x, y, z;
+ fscal x = h.x, y = h.y, z = h.z;
fscal x_t, z_t;
- algo->GetHitCoor(h, x, y, z, st);
TVector3 v3(x, y, z);
v3.RotateX(TMath::Pi() / 5);
@@ -695,8 +694,8 @@ void L1AlgoDraw::DrawRestHits(L1HitIndex_t* StsRestHitsStartIndex, L1HitIndex_t*
int iMC = CbmL1::Instance()->fvHitPointIndexes[ih];
//if( (vSFlag[h.f] | vSFlagB[h.b] )&0x02 ) continue; // if used
- fscal x, y, z;
- algo->GetHitCoor(h, x, y, z, st);
+ fscal x = h.x, y = h.y, z = h.z;
+
if (iMC >= 0) {
x_poly[n_poly] = x;
y_poly[n_poly] = y;
@@ -792,18 +791,7 @@ void L1AlgoDraw::ClearVeiw()
L1AlgoDraw::Point L1AlgoDraw::GetHitCoor(int ih)
{
L1Hit& hit = vHits[ih];
- // find station
- int ista = 0;
- for (int i = 0; i < NStations; i++) {
- if ((HitsStartIndex[i] <= ih) && (HitsStopIndex[i] > ih)) {
- ista = i;
- break;
- }
- }
- L1Station& sta = vStations[ista];
- fscal x, y, z;
- algo->GetHitCoor(hit, x, y, z, sta);
- return Point(x, y, z);
+ return Point(hit.x, hit.y, hit.z);
};
void L1AlgoDraw::SaveCanvas(TString name)
diff --git a/reco/L1/catools/CaToolsHitRecord.h b/reco/L1/catools/CaToolsHitRecord.h
index 9915a0661b8c7cef1330971d8fc3cee0de335b17..e7cba94cb9720b18ec04e134123e01d9780c7a94 100644
--- a/reco/L1/catools/CaToolsHitRecord.h
+++ b/reco/L1/catools/CaToolsHitRecord.h
@@ -22,15 +22,11 @@ namespace ca::tools
struct HitRecord {
double fX = -1.; ///< x component of hit position [cm]
double fY = -1.; ///< y component of hit position [cm]
+ double fZ = -1.; ///< z component of hit position [cm]
+ double fT = -1.; ///< time of hit [ns]
double fDx = -1.; ///< error of x component of hit position [cm]
double fDy = -1.; ///< error of y component of hit position [cm]
double fDxy = -1.; ///< correlation between x and y components [cm]
- double fU = -1.; ///< hit position in direction of front strips [cm]
- double fV = -1.; ///< hit position in direction of back strips [cm]
- double fDu = -1.; ///< hit position error in direction of front strips [cm]
- double fDv = -1.; ///< hit position error in direction of back strips [cm]
- double fZ = -1.; ///< z component of hit position [cm]
- double fT = -1.; ///< time of hit [ns]
double fDt = -1.; ///< time error of hit [ns]
int64_t fDataStream = -1; ///< Global index of detector module
int fPointId = -1; ///< index of MC point
diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx
index b8dd85b0dbf0479a998a2be8491b42a044da7ef0..be0dddbd415abce04e77cd9a07e33fc5885a4b7c 100644
--- a/reco/L1/qa/CbmCaInputQaSts.cxx
+++ b/reco/L1/qa/CbmCaInputQaSts.cxx
@@ -330,7 +330,7 @@ void CbmCaInputQaSts::FillHistograms()
double stPhiU; // STS front strips stereo angle
double stPhiV; // STS back strips stereo angle
- std::tie(stPhiU, stPhiV) = fpDetInterface->GetStripStereoAnglesSensor(pHit->GetAddress());
+ std::tie(stPhiU, stPhiV) = fpDetInterface->GetStereoAnglesSensor(pHit->GetAddress());
double sinU = sin(stPhiU);
double cosU = cos(stPhiU);
diff --git a/reco/L1/utils/CbmCaIdealHitProducerDet.h b/reco/L1/utils/CbmCaIdealHitProducerDet.h
index b4f28a8be345b27513f5069c6aadcdb2c0ff92a0..51e467223562a356eac42f6d8acdd4d56662a2dc 100644
--- a/reco/L1/utils/CbmCaIdealHitProducerDet.h
+++ b/reco/L1/utils/CbmCaIdealHitProducerDet.h
@@ -53,6 +53,7 @@
#include <yaml-cpp/yaml.h>
#include "CaAlgoRandom.h"
+#include "L1Algo/utils/CaUvConverter.h"
#include "L1Constants.h"
#include "L1Def.h"
#include "L1Undef.h"
@@ -126,6 +127,10 @@ namespace cbm::ca
/// @brief A structure to keep hit adjustment/creation settings for each station
struct HitParameters {
+
+ double fPhiU {undef::kD}; ///< Angle of the first independent measured coordinate [rad]
+ double fPhiV {undef::kD}; ///< Angle of the second independent measured coordinate [rad]
+
double fDu = undef::kD; ///< Error of u-coordinate measurement [cm]
double fDv = undef::kD; ///< Error of v-coordinate measurement [cm]
double fDt = undef::kD; ///< Error of time measurement [ns]
@@ -332,13 +337,15 @@ namespace cbm::ca
auto& node = parNode[iSt];
short mode = flagNode[iSt].as<int>(0);
if (mode > 0) {
- par.fDu = node["du"]["value"].as<double>();
+ par.fPhiU = node["du"]["angle"].as<double>() * 180. / M_PI;
+ par.fDu = node["du"]["value"].as<double>();
switch (node["du"]["pdf"].as<char>()) {
case 'g': par.fPdfU = 0; break;
case 'u': par.fPdfU = 1; break;
default: par.fPdfU = -1; break;
}
- par.fDv = node["dv"]["value"].as<double>();
+ par.fPhiV = node["dv"]["angle"].as<double>() * 180. / M_PI;
+ par.fDv = node["dv"]["value"].as<double>();
switch (node["dv"]["pdf"].as<char>()) {
case 'g': par.fPdfV = 0; break;
case 'u': par.fPdfV = 1; break;
@@ -536,19 +543,21 @@ namespace cbm::ca
double t = pPoint->GetTime() + fpMCEventList->GetEventTime(link);
if (setting.fbSmear) {
+ double dx2 = pHit->GetDx() * pHit->GetDx();
+ double dxy = pHit->GetDxy();
+ double dy2 = pHit->GetDy() * pHit->GetDy();
+
+ CaUvConverter conv(setting.fPhiU, dx2, dxy, dy2);
+
+ auto [u, v] = conv.ConvertXYtoUV(x, y);
+ auto [du2, duv, dv2] = conv.ConvertCovMatrixXYtoUV(dx2, dxy, dy2);
+
double dt = pHit->GetTimeError();
- double du = pHit->GetDx();
- double dv = pHit->GetDy();
- if constexpr (L1DetectorID::kSts == DetID) {
- du = pHit->GetDu();
- dv = pHit->GetDv();
- }
- auto [u, v] = fpDetInterface->ConvertXYtoUV(iSt, x, y);
- this->SmearValue(u, du, setting.fPdfU);
- this->SmearValue(v, dv, setting.fPdfV);
+ this->SmearValue(u, sqrt(fabs(du2)), setting.fPdfU);
+ this->SmearValue(v, sqrt(fabs(dv2)), setting.fPdfV);
this->SmearValue(t, dt, setting.fPdfT);
- std::tie(x, y) = fpDetInterface->ConvertUVtoXY(iSt, u, v);
+ std::tie(x, y) = conv.ConvertUVtoXY(u, v);
}
pHit->SetX(x);
pHit->SetY(y);
@@ -582,11 +591,13 @@ namespace cbm::ca
int iSt = fpDetInterface->GetTrackingStationIndex(pPoint->GetDetectorID());
const auto& setting = fvStationPars[iSt];
if (setting.fMode != 2) { continue; }
- double du = setting.fDu;
- double dv = setting.fDv;
- double dt = setting.fDt;
- double dz = 0.;
- auto [dx2, dxy, dy2] = fpDetInterface->ConvertCovMatrixUVtoXY(iSt, du * du, dv * dv);
+ double du = setting.fDu;
+ double dv = setting.fDv;
+ double dt = setting.fDt;
+ double dz = 0.;
+
+ CaUvConverter conv(setting.fPhiU, setting.fPhiV);
+ auto [dx2, dxy, dy2] = conv.ConvertCovMatrixUVtoXY(du * du, 0., dv * dv);
// Get point position and time
auto posIn = TVector3 {};
@@ -606,13 +617,12 @@ namespace cbm::ca
double z = pos.Z();
double t = pPoint->GetTime() + fpMCEventList->GetEventTime(eventId, fileId);
- if (setting.fbSmear) {
- auto [u, v] = fpDetInterface->ConvertXYtoUV(iSt, x, y);
- // TODO: Provide more realistic profiles for TRD
+ if (setting.fbSmear) { // TODO: Provide more realistic profiles for TRD
+ auto [u, v] = conv.ConvertXYtoUV(x, y);
this->SmearValue(u, du, setting.fPdfU);
this->SmearValue(v, dv, setting.fPdfV);
this->SmearValue(t, dt, setting.fPdfT);
- std::tie(x, y) = fpDetInterface->ConvertUVtoXY(iSt, u, v);
+ std::tie(x, y) = conv.ConvertUVtoXY(u, v);
}
pos.SetX(x);
pos.SetY(y);