diff --git a/core/base/CbmTrackingDetectorInterfaceBase.cxx b/core/base/CbmTrackingDetectorInterfaceBase.cxx index 0548f85f72e7f44db3dc14ee66f26a0f2873f25b..1cc11e8fc38f60fea3f8ef7df46cfca1f3162c46 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.cxx +++ b/core/base/CbmTrackingDetectorInterfaceBase.cxx @@ -1,11 +1,11 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ -/// \file CbmTrackingDetectorInterfaceBase.h -/// \brief Base abstract class for tracking detector interface to L1 (implementation of Checker) -/// \since 31.05.2022 -/// \author S.Zharko (s.zharko@gsi.de) +/// @file CbmTrackingDetectorInterfaceBase.h +/// @brief Base abstract class for tracking detector interface to L1 (implementation of Checker) +/// @since 31.05.2022 +/// @author S.Zharko (s.zharko@gsi.de) #include "CbmTrackingDetectorInterfaceBase.h" @@ -37,7 +37,7 @@ bool CbmTrackingDetectorInterfaceBase::Check() const { // Position along Z-axis double z0 = this->GetZmin(iSt); - double z1 = this->GetZref(iSt); + double z1 = this->GetZ(iSt); double z2 = this->GetZmax(iSt); if (!std::isfinite(z0) || !std::isfinite(z1) || !std::isfinite(z2) || !(z0 < z1 && z1 < z2)) { msg << prefix << " wrong Z position (" << z0 << " < " << z1 << " < " << z2 << " cm)\n"; @@ -63,7 +63,7 @@ bool CbmTrackingDetectorInterfaceBase::Check() const // Position along beam axis std::vector<double> zPositions(this->GetNtrackingStations()); for (int iSt = 0; iSt < this->GetNtrackingStations(); ++iSt) { - zPositions[iSt] = this->GetZref(iSt); + zPositions[iSt] = this->GetZ(iSt); } std::set<double> zPositionSet(zPositions.begin(), zPositions.end()); if (zPositions.size() != zPositionSet.size()) { @@ -91,23 +91,28 @@ std::string CbmTrackingDetectorInterfaceBase::ToString() const using std::setfill; using std::setw; std::stringstream table; - table << "\nTracking detector interface: " << setw(5) << setfill(' ') << GetDetectorName() << '\n'; - table << setw(5) << setfill(' ') << "st.No" << ' '; - table << setw(10) << setfill(' ') << "z[cm]" << ' '; - table << setw(10) << setfill(' ') << "z_min[cm]" << ' '; - table << setw(10) << setfill(' ') << "z_max[cm]" << ' '; - table << setw(10) << setfill(' ') << "x_max[cm]" << ' '; - table << setw(10) << setfill(' ') << "y_max[cm]" << ' '; - table << setw(11) << setfill(' ') << "angleF[rad]" << ' '; - table << setw(11) << setfill(' ') << "angleB[rad]" << ' '; + table << std::boolalpha; + table << "\nTracking detector interface: " << setw(5) << GetDetectorName() << '\n'; + table << setw(5) << "st.No" << ' '; + table << setw(10) << "z[cm]" << ' '; + table << setw(10) << "z_min[cm]" << ' '; + table << setw(10) << "z_max[cm]" << ' '; + table << setw(10) << "x_min[cm]" << ' '; + table << setw(10) << "x_max[cm]" << ' '; + table << setw(10) << "y_min[cm]" << ' '; + table << setw(10) << "y_max[cm]" << ' '; + table << setw(10) << "time info" << ' '; table << '\n'; for (int iSt = 0; iSt < GetNtrackingStations(); ++iSt) { - table << setw(5) << setfill(' ') << iSt << ' '; - table << setw(10) << setfill(' ') << GetZref(iSt) << ' '; - table << setw(10) << setfill(' ') << GetZmin(iSt) << ' '; - table << setw(10) << setfill(' ') << GetZmax(iSt) << ' '; - table << setw(10) << setfill(' ') << GetXmax(iSt) << ' '; - table << setw(10) << setfill(' ') << GetYmax(iSt) << ' '; + table << setw(5) << iSt << ' '; + table << setw(10) << GetZ(iSt) << ' '; + table << setw(10) << GetZmin(iSt) << ' '; + table << setw(10) << GetZmax(iSt) << ' '; + table << setw(10) << GetXmin(iSt) << ' '; + table << setw(10) << GetXmax(iSt) << ' '; + table << setw(10) << GetYmin(iSt) << ' '; + table << setw(10) << GetYmax(iSt) << ' '; + table << setw(10) << IsTimeInfoProvided(iSt) << ' '; table << '\n'; } return table.str(); diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h index ab62e284aaa025d2bba71d8bd772cdfb359d7652..b532e2deb39c73c370d13d64d41fae68927d905a 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.h +++ b/core/base/CbmTrackingDetectorInterfaceBase.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ @@ -22,80 +22,91 @@ class CbmPixelHit; -/// Abstract class, which should be inherited by every detecting subsystem tracking interface class +/// @class CbmTrackingDetectorInterfaceBase +/// @brief Abstract class, which should be inherited by every detecting subsystem tracking interface class /// class CbmTrackingDetectorInterfaceBase { public: - /// Virtual destructor + /// @brief Virtual destructor virtual ~CbmTrackingDetectorInterfaceBase() {} + + /// @brief Checks detector interface: boundary conditions of the parameters + bool Check() const; - /// Returns the name of the detector subsystem + /// @brief Returns the name of the detector subsystem virtual std::string GetDetectorName() const = 0; - /// Gets actual number of stations, provided by the current geometry setup + /// @brief Gets actual number of stations, provided by the current geometry setup virtual int GetNtrackingStations() const = 0; - /// Gets a tracking station of a CbmPixelHit - /// \param hit A pointer to CbmPixelHit - /// \return Local index of the tracking station + /// @brief Gets stereo angles of the two independent measured coordinates + /// @note 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; + + /// @brief Gets a tracking station of a CbmPixelHit + /// @param hit A pointer to CbmPixelHit + /// @return Local index of the tracking station virtual int GetTrackingStationIndex(const CbmPixelHit* hit) const = 0; - /// Gets a tracking station by the address of element - /// \param address Unique element address - /// \return Local index of the tracking station + /// @brief Gets a tracking station by the address of element + /// @param address Unique element address + /// @return Local index of the tracking station virtual int GetTrackingStationIndex(int address) const = 0; - /// Gets a tracking station index of MC point + /// @brief Gets a tracking station index of MC point /// @param zPos z-coordinate of hit/point position [cm] /// - /// Function searches the closest station in beam axis direction + /// The function searches for the closest station in beam axis direction int GetTrackingStationIndex(double zPos) const; - - /// Gets max size of a station along the X-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the X-axis [cm] + + /// @brief Gets upper bound of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] virtual double GetXmax(int stationId) const = 0; - /// Gets max size of a station along the Y-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the Y-axis [cm] - virtual double GetYmax(int stationId) const = 0; + /// @brief Gets lower bound of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] + virtual double GetXmin(int stationId) const = 0; - /// Gets reference z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of the station [cm] - virtual double GetZref(int stationId) const = 0; + /// @brief Gets upper bound of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + virtual double GetYmax(int stationId) const = 0; - /// Gets min z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return min Z of the station [cm] + /// @brief Gets lower bound of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + virtual double GetYmin(int stationId) const = 0; + + /// @brief Gets reference z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Z position of the station [cm] + virtual double GetZ(int stationId) const = 0; + + /// @brief Gets min z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return min Z of the station [cm] virtual double GetZmin(int stationId) const = 0; - /// Gets max z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return max Z of the station [cm] + /// @brief Gets max z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return max Z of the station [cm] virtual double GetZmax(int stationId) const = 0; - /// Check if station provides time measurements - /// \param stationId Tracking station ID in the setup - /// \return Flag: true - station provides time measurements, false - station does not provide time measurements + /// @brief Check if station provides time measurements + /// @param stationId Tracking station ID in the setup + /// @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) + /// 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}; - /// Checks detector interface: boundary conditions of the parameters - bool Check() const; - - /// Prints all the parameters into table and saves the table as a string + /// @brief Prints all the parameters into table and saves the table as a string std::string ToString() const; }; @@ -111,7 +122,7 @@ inline int CbmTrackingDetectorInterfaceBase::GetTrackingStationIndex(double zPos int nStations = GetNtrackingStations(); int iStSelected = -1; for (int iSt = 0; iSt < nStations; ++iSt) { - auto dist = std::fabs(zPos - GetZref(iSt)); + auto dist = std::fabs(zPos - GetZ(iSt)); if (dist < bestDist) { bestDist = dist; iStSelected = iSt; diff --git a/core/detectors/much/CbmMuchTrackingInterface.h b/core/detectors/much/CbmMuchTrackingInterface.h index c5c78cf7bde686ba9e6078356306042d0bb5fdc4..7f15c29783f27188e8bcd3531e6def9093c31d5e 100644 --- a/core/detectors/much/CbmMuchTrackingInterface.h +++ b/core/detectors/much/CbmMuchTrackingInterface.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ @@ -29,99 +29,121 @@ class CbmMuchLayer; -/// Class CbmMuchTrackingInterface is a CbmL1 subtask, which provides necessary methods for L1 tracker -/// to access the geometry and dataflow settings. +/// @class CbmMuchTrackingInterface +/// @brief A CbmL1 subtask, which provides necessary methods for L1 tracker to access the geometry and dataflow +/// settings. /// -/// NOTE: For MuCh one tracking station is a MuCh layer! +/// @note For MuCh one tracking station is a MuCh layer! /// class CbmMuchTrackingInterface : public FairTask, public CbmTrackingDetectorInterfaceBase { public: - /// Default constructor + /// @brief Default constructor CbmMuchTrackingInterface(); - /// Destructor + /// @brief Destructor ~CbmMuchTrackingInterface(); + + /// @brief Copy constructor + CbmMuchTrackingInterface(const CbmMuchTrackingInterface&) = delete; + + /// @brief Move constructor + CbmMuchTrackingInterface(CbmMuchTrackingInterface&&) = delete; + + /// @brief Copy assignment operator + CbmMuchTrackingInterface& operator=(const CbmMuchTrackingInterface&) = delete; + + /// @brief Move assignment operator + CbmMuchTrackingInterface& operator=(CbmMuchTrackingInterface&&) = delete; - /// FairTask: Init method + /// @brief FairTask: Init method InitStatus Init(); - /// FairTask: ReInit method + /// @brief FairTask: ReInit method InitStatus ReInit(); - /// Gets pointer to the instance of the CbmMuchTrackingInterface + /// @brief Gets pointer to the instance of the CbmMuchTrackingInterface static CbmMuchTrackingInterface* Instance() { return fpInstance; } - /// Gets name of this subsystem + /// @brief Gets name of this subsystem std::string GetDetectorName() const { return "MuCh"; } - /// Gets actual number of tracking stations, provided by the current geometry setup + /// @brief Gets actual number of tracking stations, provided by the current geometry setup int GetNtrackingStations() const; - /// Gets reference z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of the station [cm] - double GetZref(int stationId) const { return GetMuchLayer(stationId)->GetZ(); } + /// @brief Gets reference z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Z position of the station [cm] + double GetZ(int stationId) const { return GetMuchLayer(stationId)->GetZ(); } - /// Gets min z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return min Z of the station [cm] - double GetZmin(int stationId) const { return GetZref(stationId) - 0.5 * GetMuchLayer(stationId)->GetDz(); } + /// @brief Gets min z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return min Z of the station [cm] + double GetZmin(int stationId) const { return GetZ(stationId) - 0.5 * GetMuchLayer(stationId)->GetDz(); } - /// Gets max z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return max Z of the station [cm] - double GetZmax(int stationId) const { return GetZref(stationId) + 0.5 * GetMuchLayer(stationId)->GetDz(); } + /// @brief Gets max z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return max Z of the station [cm] + double GetZmax(int stationId) const { return GetZ(stationId) + 0.5 * GetMuchLayer(stationId)->GetDz(); } - /// Gets a tracking station of a CbmPixelHit - /// \param hit A pointer to CbmPixelHit - /// \return Local index of the tracking station + /// @brief Gets a tracking station of a CbmPixelHit + /// @param hit A pointer to CbmPixelHit + /// @return Local index of the tracking station int GetTrackingStationIndex(const CbmPixelHit* hit) const { return GetTrackingStationIndex(hit->GetAddress()); } - /// Gets a tracking station by the address of element - /// \param address Unique element address - /// \return Local index of the tracking station + /// @brief Gets a tracking station by the address of element + /// @param address Unique element address + /// @return Local index of the tracking station int GetTrackingStationIndex(int address) const { return CbmMuchGeoScheme::GetStationIndex(address); } - /// Gets max size of a station along the X-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the X-axis [cm] + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets upper bound of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] double GetXmax(int /*stationId*/) const { return 100.; } - /// Gets max size of a station along the Y-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the Y-axis [cm] + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets lower bound of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] + double GetXmin(int /*stationId*/) const { return -100.; } + + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets upper bound of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] double GetYmax(int /*stationId*/) const { return 100.; } - /// Check if station provides time measurements - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Flag: true - station provides time measurements, false - station does not provide time measurements + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets lower bound of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + double GetYmin(int /*stationId*/) const { return -100.; } + + /// @brief Check if station provides time measurements + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @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] + /// @brief Gets stereo angles of the two independent measured coordinates + /// @note 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 + /// @brief FairTask: sets parameter containers up void SetParContainers(); - /// Copy and move constructers and assign operators are prohibited - CbmMuchTrackingInterface(const CbmMuchTrackingInterface&) = delete; - CbmMuchTrackingInterface(CbmMuchTrackingInterface&&) = delete; - CbmMuchTrackingInterface& operator=(const CbmMuchTrackingInterface&) = delete; - CbmMuchTrackingInterface& operator=(CbmMuchTrackingInterface&&) = delete; - private: - /// Gets pointer to the TRD module - /// \param stationId Tracking staton ID - /// \return Pointer to the particular CbmTrdParModDigi object + /// @brief Gets pointer to the TRD module + /// @param stationId Tracking staton ID + /// @return Pointer to the particular CbmTrdParModDigi object __attribute__((always_inline)) CbmMuchLayer* GetMuchLayer(int stationId) const { return fGeoScheme->GetLayer(GetMuchStationId(stationId), GetMuchLayerId(stationId)); } - /// Calculates MuCh layer ID from tracker station ID + /// @brief Calculates MuCh layer ID from tracker station ID + /// @param stationId Index of the tracking station __attribute__((always_inline)) int GetMuchLayerId(int stationId) const { int layersCount = 0; @@ -133,7 +155,8 @@ private: return stationId - layersCount; } - /// Calculates MuCh station ID from tracker station ID + /// @brief Calculates MuCh station ID from tracker station ID + /// @param stationId Index of the tracking station __attribute__((always_inline)) int GetMuchStationId(int stationId) const { int layersCount = 0; @@ -146,7 +169,6 @@ private: return stationsCount; } - static CbmMuchTrackingInterface* fpInstance; ///< Instance of the class CbmMuchGeoScheme* fGeoScheme {nullptr}; ///< MuCh geometry scheme instance diff --git a/core/detectors/mvd/CbmMvdTrackingInterface.h b/core/detectors/mvd/CbmMvdTrackingInterface.h index 47e5201530c0a9c883d7aef9a327d86a30a1b550..d8410cb07392eeb882d7db265989f6e675c3e90c 100644 --- a/core/detectors/mvd/CbmMvdTrackingInterface.h +++ b/core/detectors/mvd/CbmMvdTrackingInterface.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ @@ -30,50 +30,60 @@ class TBuffer; class TClass; class TMemberInspector; -/// Class CbmMvdTrackingInterface is a CbmL1 subtask, which provides necessary methods for L1 tracker -/// to access the geometry and dataflow settings. +/// @brief CbmMvdTrackingInterface +/// @brief A CbmL1 subtask, which provides necessary methods for L1 tracker to access the geometry and dataflow +/// settings. /// class CbmMvdTrackingInterface : public FairTask, public CbmTrackingDetectorInterfaceBase, public CbmMvdDetectorId { public: - /// Default constructor + /// @brief Default constructor CbmMvdTrackingInterface(); - /// Destructor + /// @brief Destructor ~CbmMvdTrackingInterface(); - /// FairTask: Init method - InitStatus Init(); - - /// FairTask: ReInit method - InitStatus ReInit(); - - /// Gets pointer to the instance of the CbmMvdTrackingInterface - static CbmMvdTrackingInterface* Instance() { return fpInstance; } + /// @brief Copy constructor + CbmMvdTrackingInterface(const CbmMvdTrackingInterface&) = delete; + + /// @brief Move constructor + CbmMvdTrackingInterface(CbmMvdTrackingInterface&&) = delete; + + /// @brief Copy assignment operator + CbmMvdTrackingInterface& operator=(const CbmMvdTrackingInterface&) = delete; + + /// @brief Move assignment operator + CbmMvdTrackingInterface& operator=(CbmMvdTrackingInterface&&) = delete; - /// Gets name of this subsystem + /// @brief Gets name of this subsystem std::string GetDetectorName() const { return "MVD"; } - /// Gets actual number of tracking stations, provided by the current geometry setup + /// @brief Gets actual number of tracking stations, provided by the current geometry setup int GetNtrackingStations() const { return fMvdStationPar->GetStationCount(); } - /// Gets reference z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of the station [cm] - double GetZref(int stationId) const { return fMvdStationPar->GetZPosition(stationId); } - - /// Gets min z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return min Z of the station [cm] - double GetZmin(int stationId) const { return GetZref(stationId) - 0.5 * fMvdStationPar->GetZThickness(stationId); } + /// --- to be removed --- Gets the tracking station radiation length + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Radiation length [cm] + [[deprecated]] + double GetRadLength(int stationId) const + { + return fMvdStationPar->GetZThickness(stationId) / (10. * fMvdStationPar->GetZRadThickness(stationId)); + } - /// Gets max z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return max Z of the station [cm] - double GetZmax(int stationId) const { return GetZref(stationId) + 0.5 * fMvdStationPar->GetZThickness(stationId); } + /// @brief Gets stereo angles of the two independent measured coordinates + /// @note 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.); } - /// Gets a tracking station of a CbmPixelHit - /// \param hit A pointer to CbmPixelHit - /// \return Local index of the tracking station + /// --- to be removed --- 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] + [[deprecated]] + double GetSensorThickness(int stationId) const { return fMvdStationPar->GetZThickness(stationId); } + + /// @brief Gets a tracking station of a CbmPixelHit + /// @param hit A pointer to CbmPixelHit + /// @return Local index of the tracking station int GetTrackingStationIndex(const CbmPixelHit* hit) const { auto hitMvd = [&] { @@ -85,53 +95,62 @@ public: return hitMvd->GetStationNr(); } - /// Gets a tracking station by the address of element (detectorID in terms of MVD) - /// \param detectorId Unique element address (detectorID in terms of MVD) - /// \return Local index of the tracking station + /// @brief Gets a tracking station by the address of element (detectorID in terms of MVD) + /// @param detectorId Unique element address (detectorID in terms of MVD) + /// @return Local index of the tracking station int GetTrackingStationIndex(int detectorId) const { return StationNr(detectorId); } - /// Gets max size of a tracking station along the X-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the X-axis [cm] + /// @brief Gets upper bound of a tracking station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] double GetXmax(int stationId) const { return fMvdStationPar->GetWidth(stationId); } - - /// Gets max size of a tracking station along the Y-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the Y-axis [cm] + + /// @brief Gets lower bound of a tracking station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] + double GetXmin(int stationId) const { return -fMvdStationPar->GetWidth(stationId); } + + /// @brief Gets upper bound of a tracking station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] double GetYmax(int stationId) const { return fMvdStationPar->GetHeight(stationId); } - /// Check if the detector provides time measurements - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Flag: true - station provides time measurements, false - station does not provide time measurements - bool IsTimeInfoProvided(int /*stationId*/) const { return false; } + /// @brief Gets lower bound of a tracking station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + double GetYmin(int stationId) const { return -fMvdStationPar->GetHeight(stationId); } - /// 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.); } + /// @brief Gets reference z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Z position of the station [cm] + double GetZ(int stationId) const { return fMvdStationPar->GetZPosition(stationId); } - /// FairTask: sets parameter containers up - void SetParContainers(); + /// @brief Gets max z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return max Z of the station [cm] + double GetZmax(int stationId) const { return GetZ(stationId) + 0.5 * fMvdStationPar->GetZThickness(stationId); } - /// --- to be removed --- Gets the tracking station radiation length - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Radiation length [cm] - double GetRadLength(int stationId) const - { - return fMvdStationPar->GetZThickness(stationId) / (10. * fMvdStationPar->GetZRadThickness(stationId)); - } + /// @brief Gets min z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return min Z of the station [cm] + double GetZmin(int stationId) const { return GetZ(stationId) - 0.5 * fMvdStationPar->GetZThickness(stationId); } - /// --- to be removed --- 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] - double GetSensorThickness(int stationId) const { return fMvdStationPar->GetZThickness(stationId); } + /// @brief FairTask: Init method + InitStatus Init(); - /// Copy and move constructers and assign operators are prohibited - CbmMvdTrackingInterface(const CbmMvdTrackingInterface&) = delete; - CbmMvdTrackingInterface(CbmMvdTrackingInterface&&) = delete; - CbmMvdTrackingInterface& operator=(const CbmMvdTrackingInterface&) = delete; - CbmMvdTrackingInterface& operator=(CbmMvdTrackingInterface&&) = delete; + /// @brief Gets pointer to the instance of the CbmMvdTrackingInterface + static CbmMvdTrackingInterface* Instance() { return fpInstance; } + + /// @brief Check if the detector provides time measurements + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Flag: true - station provides time measurements, false - station does not provide time measurements + bool IsTimeInfoProvided(int /*stationId*/) const { return false; } + + /// @brief FairTask: ReInit method + InitStatus ReInit(); + + /// @brief FairTask: sets parameter containers up + void SetParContainers(); private: static CbmMvdTrackingInterface* fpInstance; ///< Instance of the class diff --git a/core/detectors/sts/CbmStsTrackingInterface.h b/core/detectors/sts/CbmStsTrackingInterface.h index f78b2fd8165acad3189b02261a483d5d55c79e2d..9b2484618fb759edc58f61f52f778c4aa1ad2f9d 100644 --- a/core/detectors/sts/CbmStsTrackingInterface.h +++ b/core/detectors/sts/CbmStsTrackingInterface.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ @@ -27,118 +27,123 @@ #include <iostream> -/// Class CbmStsTrackingInterface is a CbmL1 subtask, which provides necessary methods for L1 tracker -/// to access the geometry and dataflow settings. +/// @class CbmStsTrackingInterface +/// @brief A CbmL1 subtask, which provides necessary methods for CA tracker to access the geometry and dataflow +/// settings /// class CbmStsTrackingInterface : public FairTask, public CbmTrackingDetectorInterfaceBase { public: - /// Default constructor + /// @brief Default constructor CbmStsTrackingInterface(); - /// Destructor + /// @brief Destructor ~CbmStsTrackingInterface(); - /// FairTask: Init method + /// @brief Copy constructor + CbmStsTrackingInterface(const CbmStsTrackingInterface&) = delete; + + /// @brief Move constructor + CbmStsTrackingInterface(CbmStsTrackingInterface&&) = delete; + + /// @brief Copy assignment operator + CbmStsTrackingInterface& operator=(const CbmStsTrackingInterface&) = delete; + + /// @brief Move assignment operator + CbmStsTrackingInterface& operator=(CbmStsTrackingInterface&&) = delete; + + + /// @brief FairTask: Init method InitStatus Init(); - /// FairTask: ReInit method + /// @brief FairTask: ReInit method InitStatus ReInit(); - /// Gets pointer to the instance of the CbmStsTrackingInterface class + /// @brief Gets pointer to the instance of the CbmStsTrackingInterface class static CbmStsTrackingInterface* Instance() { return fpInstance; } - /// Gets name of this subsystem + /// @brief Gets name of this subsystem std::string GetDetectorName() const { return "STS"; } - /// Gets actual number of the tracking stations, provided by the current geometry setup + /// @brief Gets actual number of the tracking stations, provided by the current geometry setup int GetNtrackingStations() const { return CbmStsSetup::Instance()->GetNofStations(); } - /// 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] + /// @brief Gets stereo angles of the two independent measured coordinates + /// @note 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 a tracking station of a CbmPixelHit - /// \param hit A pointer to CbmPixelHit - /// \return Local index of the tracking station + /// @brief Gets a tracking station of a CbmPixelHit + /// @param hit A pointer to CbmPixelHit + /// @return Local index of the tracking station int GetTrackingStationIndex(const CbmPixelHit* hit) const { return GetTrackingStationIndex(hit->GetAddress()); } - /// Gets a tracking station by the address of element - /// \param address Unique element address - /// \return Local index of the tracking station + /// @brief Gets a tracking station by the address of element + /// @param address Unique element address + /// @return Local index of the tracking station int GetTrackingStationIndex(int address) const { return CbmStsSetup::Instance()->GetStationNumber(address); } - /// Gets min size of a station along the X-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the X-axis [cm] - double GetXmin(int stationId) const { return GetStsStation(stationId)->GetXmin(); } - - /// Gets max size of a station along the X-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the X-axis [cm] + /// @brief Gets upper edge of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] double GetXmax(int stationId) const { return GetStsStation(stationId)->GetXmax(); } - /// Gets min size of a station along the Y-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the Y-axis [cm] - double GetYmin(int stationId) const { return GetStsStation(stationId)->GetYmin(); } + /// @brief Gets lower edge of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] + double GetXmin(int stationId) const { return GetStsStation(stationId)->GetXmin(); } - /// Gets max size of a station along the Y-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the Y-axis [cm] + /// @brief Gets max size of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] double GetYmax(int stationId) const { return GetStsStation(stationId)->GetYmax(); } - /// Gets reference Z position of the station: approximately == mean Z of all its component - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of station [cm] - double GetZref(int stationId) const { return GetStsStation(stationId)->GetZ(); } + /// @brief Gets lower edge of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + double GetYmin(int stationId) const { return GetStsStation(stationId)->GetYmin(); } - /// Gets reference Z position of the station: approximately == mean Z of all its component - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of station [cm] - double GetZ(int stationId) const { return GetZref(stationId); } + /// @brief Gets reference Z position of the station: approximately == mean Z of all its component + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Z position of station [cm] + double GetZ(int stationId) const { return GetStsStation(stationId)->GetZ(); } - /// Gets minimal z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return min z position of station [cm] + /// @brief Gets maximal z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return max z position of station [cm] // TODO: calculate from components - double GetZmin(int stationId) const { return GetStsStation(stationId)->GetZ() - 2.; } + double GetZmax(int stationId) const { return GetStsStation(stationId)->GetZ() + 2.; } - /// Gets maximal z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return max z position of station [cm] + /// @brief Gets minimal z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return min z position of station [cm] // TODO: calculate from components - double GetZmax(int stationId) const { return GetStsStation(stationId)->GetZ() + 2.; } + double GetZmin(int stationId) const { return GetStsStation(stationId)->GetZ() - 2.; } - /// Check if station provides time measurements - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Flag: true - station provides time measurements, false - station does not provide time measurements + /// @brief Check if station provides time measurements + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Flag: true - station provides time measurements, false - station does not provide time measurements bool IsTimeInfoProvided(int /*stationId*/) const { return true; } /// -- to be removed -- Gets station radiation length /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) /// \return Radiation length [cm] + [[deprecated]] double GetRadLength(int stationId) const { return GetStsStation(stationId)->GetRadLength(); } /// -- to be removed -- 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] + [[deprecated]] double GetSensorThickness(int stationId) const { return GetStsStation(stationId)->GetSensorD(); } - /// FairTask: sets parameter containers up + /// @brief FairTask: sets parameter containers up void SetParContainers(); - /// Copy and move constructors and assign operators are forbidden - CbmStsTrackingInterface(const CbmStsTrackingInterface&) = delete; - CbmStsTrackingInterface(CbmStsTrackingInterface&&) = delete; - CbmStsTrackingInterface& operator=(const CbmStsTrackingInterface&) = delete; - CbmStsTrackingInterface& operator=(CbmStsTrackingInterface&&) = delete; - private: - /// Gets pointer to the STS station object - /// \param stationId Tracking staton ID - /// \return Pointer to the particular CbmStsStation object + /// @brief Gets pointer to the STS station object + /// @param stationId Tracking staton ID + /// @return Pointer to the particular CbmStsStation object __attribute__((always_inline)) CbmStsStation* GetStsStation(int stationId) const { return CbmStsSetup::Instance()->GetStation(stationId); diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h index 662f4f52d8308bde17c57f5ff1a25edd93425643..81833eab0ff3f9a422e3169f34eddbb872fb99ee 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.h +++ b/core/detectors/tof/CbmTofTrackingInterface.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ @@ -35,50 +35,44 @@ class CbmTofDigiPar; /// class CbmTofTrackingInterface : public FairTask, public CbmTrackingDetectorInterfaceBase { public: - /// Default constructor + /// @brief Default constructor CbmTofTrackingInterface(); - /// Destructor + /// @brief Destructor ~CbmTofTrackingInterface(); + + /// @brief Copy constructor + CbmTofTrackingInterface(const CbmTofTrackingInterface&) = delete; + + /// @brief Move constructor + CbmTofTrackingInterface(CbmTofTrackingInterface&&) = delete; + + /// @brief Copy assignment operator + CbmTofTrackingInterface& operator=(const CbmTofTrackingInterface&) = delete; + + /// @brief Move assignment operator + CbmTofTrackingInterface& operator=(CbmTofTrackingInterface&&) = delete; - /// FairTask: Init method - InitStatus Init(); - - /// FairTask: ReInit method - InitStatus ReInit(); - - /// Gets pointer to the instance of the CbmTofTrackingInterface - static CbmTofTrackingInterface* Instance() { return fpInstance; } - - /// Gets name of this subsystem + /// @brief Gets name of this subsystem std::string GetDetectorName() const { return "TOF"; } - /// Gets actual number of tracking stations, provided by the current geometry setup + /// @brief Gets actual number of tracking stations, provided by the current geometry setup int GetNtrackingStations() const { return fDigiBdfPar->GetNbTrackingStations(); } - /// Gets reference z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of the station [cm] - double GetZref(int stationId) const { return fTofStationZ[stationId]; } - - /// Gets min z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return min Z of the station [cm] - double GetZmin(int stationId) const { return GetZref(stationId) - 5.; } - - /// Gets max z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return max Z of the station [cm] - double GetZmax(int stationId) const { return GetZref(stationId) + 5.; } + /// @brief Gets stereo angles of the two independent measured coordinates + /// @note 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.); } - /// Gets a tracking station of a ToF hit - /// \param hit A pointer to ToF hit - /// \return Local index of the tracking station + /// @brief Gets a tracking station of a ToF hit + /// @param hit A pointer to ToF hit + /// @return Local index of the tracking station int GetTrackingStationIndex(const CbmPixelHit* hit) const { return GetTrackingStationIndex(hit->GetAddress()); } - /// Gets a tracking station by the address of element - /// \param address Unique element address - /// \return Local index of the tracking station + /// @brief Gets a tracking station by the address of element + /// @param address Unique element address + /// @return Local index of the tracking station int GetTrackingStationIndex(int address) const { int iSt = fDigiBdfPar->GetTrackingStation(CbmTofAddress::GetSmType(address), CbmTofAddress::GetSmId(address), @@ -86,36 +80,61 @@ public: return iSt; } + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets upper bound of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] + double GetXmax(int /*stationId*/) const { return 100.; } + + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets lower bound of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] + double GetXmin(int /*stationId*/) const { return -100.; } + + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets max size of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + double GetYmax(int /*stationId*/) const { return 100.; } + + // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + /// @brief Gets max size of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + double GetYmin(int /*stationId*/) const { return -100.; } + + /// @brief Gets reference z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Z position of the station [cm] + double GetZ(int stationId) const { return fTofStationZ[stationId]; } - /// Gets max size of a station along the X-axis + /// Gets max z of the station /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the X-axis [cm] - double GetXmax(int /*stationId*/) const { return 20.; } + /// \return max Z of the station [cm] + double GetZmax(int stationId) const { return GetZ(stationId) + 5.; } - /// Gets max size of a station along the Y-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the Y-axis [cm] - double GetYmax(int /*stationId*/) const { return 20.; } + /// @brief Gets min z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return min Z of the station [cm] + double GetZmin(int stationId) const { return GetZ(stationId) - 5.; } - /// Check if station provides time measurements - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Flag: true - station provides time measurements, false - station does not provide time measurements + /// @brief Check if station provides time measurements + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @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.); } + /// @brief FairTask: Init method + InitStatus Init(); - /// FairTask: sets parameter containers up - void SetParContainers(); + /// @brief Gets pointer to the instance of the CbmTofTrackingInterface + static CbmTofTrackingInterface* Instance() { return fpInstance; } - /// Copy and move constructers and assign operators are prohibited - CbmTofTrackingInterface(const CbmTofTrackingInterface&) = delete; - CbmTofTrackingInterface(CbmTofTrackingInterface&&) = delete; - CbmTofTrackingInterface& operator=(const CbmTofTrackingInterface&) = delete; - CbmTofTrackingInterface& operator=(CbmTofTrackingInterface&&) = delete; + /// @brief FairTask: ReInit method + InitStatus ReInit(); + + /// @brief FairTask: sets parameter containers up + void SetParContainers(); private: static CbmTofTrackingInterface* fpInstance; ///< Instance of the class diff --git a/core/detectors/trd/CbmTrdTrackingInterface.h b/core/detectors/trd/CbmTrdTrackingInterface.h index 590c3fa4a160e2ff818f2cedab9f2ce436087906..5f03e1e20668a80bd8b7e4c173f5ccaefcd8dd18 100644 --- a/core/detectors/trd/CbmTrdTrackingInterface.h +++ b/core/detectors/trd/CbmTrdTrackingInterface.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ @@ -25,93 +25,110 @@ #include <iostream> #include <vector> -/// Class CbmTrdTrackingInterface is a CbmL1 subtask, which provides necessary methods for L1 tracker -/// to access the geometry and dataflow settings. +/// @class CbmTrdTrackingInterface +/// @brief A CbmL1 subtask, which provides necessary methods for CA tracker to access the geometry and dataflow s +/// ettings. /// -/// NOTE: For TRD one tracking station is a TRD module! +/// @note For TRD one tracking station is a TRD module! /// class CbmTrdTrackingInterface : public FairTask, public CbmTrackingDetectorInterfaceBase { public: - /// Default constructor + /// @brief Default constructor CbmTrdTrackingInterface(); - /// Destructor + /// @brief Destructor ~CbmTrdTrackingInterface(); - /// FairTask: Init method - InitStatus Init(); - - /// FairTask: ReInit method - InitStatus ReInit(); - - /// Gets pointer to the instance of the CbmTrdTrackingInterface - static CbmTrdTrackingInterface* Instance() { return fpInstance; } + /// @brief Copy constructor + CbmTrdTrackingInterface(const CbmTrdTrackingInterface&) = delete; + + /// @brief Move constructor + CbmTrdTrackingInterface(CbmTrdTrackingInterface&&) = delete; + + /// @brief Copy assignment operator + CbmTrdTrackingInterface& operator=(const CbmTrdTrackingInterface&) = delete; + + /// @brief Move assignment operator + CbmTrdTrackingInterface& operator=(CbmTrdTrackingInterface&&) = delete; - /// Gets name of this subsystem + /// @brief Gets name of this subsystem std::string GetDetectorName() const { return "TRD"; } - /// Gets actual number of tracking stations, provided by the current geometry setup + /// @brief Gets actual number of tracking stations, provided by the current geometry setup int GetNtrackingStations() const; - /// Gets reference z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of the station [cm] - double GetZref(int stationId) const { return GetTrdModulePar(stationId)->GetZ(); } - - /// Gets min z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return min Z of the station [cm] - double GetZmin(int stationId) const { return GetZref(stationId) - GetTrdModulePar(stationId)->GetSizeZ(); } - - /// Gets max z of the station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return max Z of the station [cm] - double GetZmax(int stationId) const { return GetZref(stationId) + GetTrdModulePar(stationId)->GetSizeZ(); } + /// @brief Gets stereo angles of the two independent measured coordinates + /// @note 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.); } - /// Gets a tracking station of a CbmPixelHit - /// \param hit A pointer to CbmPixelHit - /// \return Local index of the tracking station + /// @brief Gets a tracking station of a CbmPixelHit + /// @param hit A pointer to CbmPixelHit + /// @return Local index of the tracking station int GetTrackingStationIndex(const CbmPixelHit* hit) const { return GetTrackingStationIndex(hit->GetAddress()); } - /// Gets a tracking station by the address - /// \param address Unique element address - /// \return Local index of the tracking station + /// @brief Gets a tracking station by the address + /// @param address Unique element address + /// @return Local index of the tracking station int GetTrackingStationIndex(int address) const { return CbmTrdAddress::GetLayerId(address); } - /// Gets max size of a station along the X-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the X-axis [cm] + /// @brief Gets upper edge of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] double GetXmax(int stationId) const { return GetTrdModulePar(stationId)->GetSizeX(); } - /// Gets max size of a station along the Y-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Size of station along the Y-axis [cm] + /// @brief Gets lower edge coordinate of a station along the X-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the X-axis [cm] + double GetXmin(int stationId) const { return -GetTrdModulePar(stationId)->GetSizeX(); } + + /// @brief Gets upper edge of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] double GetYmax(int stationId) const { return GetTrdModulePar(stationId)->GetSizeY(); } - /// Check if station provides time measurements - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Flag: true - station provides time measurements, false - station does not provide time measurements + /// @brief Gets lower edge of a station along the Y-axis + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Size of station along the Y-axis [cm] + double GetYmin(int stationId) const { return -GetTrdModulePar(stationId)->GetSizeY(); } + + /// @brief Gets reference z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return Z position of the station [cm] + double GetZ(int stationId) const { return GetTrdModulePar(stationId)->GetZ(); } + + /// @brief Gets max z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return max Z of the station [cm] + double GetZmax(int stationId) const { return GetZ(stationId) + GetTrdModulePar(stationId)->GetSizeZ(); } + + /// @brief Gets min z of the station + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @return min Z of the station [cm] + double GetZmin(int stationId) const { return GetZ(stationId) - GetTrdModulePar(stationId)->GetSizeZ(); } + + /// @brief Check if station provides time measurements + /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @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.); } + /// @brief FairTask: Init method + InitStatus Init(); - /// FairTask: sets parameter containers up - void SetParContainers(); + /// @brief Gets pointer to the instance of the CbmTrdTrackingInterface + static CbmTrdTrackingInterface* Instance() { return fpInstance; } - /// Copy and move constructers and assign operators are prohibited - CbmTrdTrackingInterface(const CbmTrdTrackingInterface&) = delete; - CbmTrdTrackingInterface(CbmTrdTrackingInterface&&) = delete; - CbmTrdTrackingInterface& operator=(const CbmTrdTrackingInterface&) = delete; - CbmTrdTrackingInterface& operator=(CbmTrdTrackingInterface&&) = delete; + /// @brief FairTask: ReInit method + InitStatus ReInit(); + + /// @brief FairTask: sets parameter containers up + void SetParContainers(); private: - /// Gets pointer to the TRD module - /// \param moduleId Id of the Trd module - /// \return Pointer to the particular CbmTrdParModDigi object + /// @brief Gets pointer to the TRD module + /// @param moduleId Id of the Trd module + /// @return Pointer to the particular CbmTrdParModDigi object __attribute__((always_inline)) CbmTrdParModDigi* GetTrdModulePar(int moduleId) const { return static_cast<CbmTrdParModDigi*>(fTrdDigiPar->GetModulePar(fTrdDigiPar->GetModuleId(moduleId))); diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C index d4ec8fcedd1af27ced867db1bddedf52e6386d94..52709d434d06653a34f011fe97ccbfde35765698 100644 --- a/macro/mcbm/mcbm_qa.C +++ b/macro/mcbm/mcbm_qa.C @@ -212,28 +212,28 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test", run->AddTask(muchHitFinderQa); // CA Input QA - auto* pCaInputMuch = new CbmCaInputQaMuch(verbose, bUseMC); - pCaInputMuch->SetEfficiencyThrsh(0.5, 0, 100); - run->AddTask(pCaInputMuch); + //auto* pCaInputMuch = new CbmCaInputQaMuch(verbose, bUseMC); + //pCaInputMuch->SetEfficiencyThrsh(0.5, 0, 100); + //run->AddTask(pCaInputMuch); } // ------------------------------------------------------------------------ // ----- TRD QA ----------------------------------------------------------- if (bUseTrd) { // CA Input QA - auto* pCaInputTrd = new CbmCaInputQaTrd(verbose, bUseMC); - pCaInputTrd->SetEfficiencyThrsh(0.5, 0, 100); - //pCaInputTrd->SetConfigName(qaConfig); - run->AddTask(pCaInputTrd); + //auto* pCaInputTrd = new CbmCaInputQaTrd(verbose, bUseMC); + //pCaInputTrd->SetEfficiencyThrsh(0.5, 0, 100); + ////pCaInputTrd->SetConfigName(qaConfig); + //run->AddTask(pCaInputTrd); } // ------------------------------------------------------------------------ // ----- TOF QA ----------------------------------------------------------- if (bUseTof) { // CA Input QA - auto* pCaInputTof = new CbmCaInputQaTof(verbose, bUseMC); - pCaInputTof->SetEfficiencyThrsh(0.5, 0, 100); - run->AddTask(pCaInputTof); + //auto* pCaInputTof = new CbmCaInputQaTof(verbose, bUseMC); + //pCaInputTof->SetEfficiencyThrsh(0.5, 0, 100); + //run->AddTask(pCaInputTof); } // ------------------------------------------------------------------------ diff --git a/reco/KF/CbmKF.cxx b/reco/KF/CbmKF.cxx index f0ed83263a13979d1a5131faa60c71c71f058c2b..102359af8bef87a71432210df9c85e2bb2f6dbee 100644 --- a/reco/KF/CbmKF.cxx +++ b/reco/KF/CbmKF.cxx @@ -157,7 +157,7 @@ InitStatus CbmKF::Init() tube.ID = 1101 + ist; // tube.F = 1.; - tube.z = mvdInterface->GetZref(ist); + tube.z = mvdInterface->GetZ(ist); tube.dz = mvdInterface->GetSensorThickness(ist); // TODO: verify the thickness of MVD stations tube.RadLength = mvdInterface->GetRadLength(ist); @@ -193,7 +193,7 @@ InitStatus CbmKF::Init() tube.ID = 1000 + ist; tube.F = 1.; - tube.z = stsInterface->GetZref(ist); + tube.z = stsInterface->GetZ(ist); tube.dz = stsInterface->GetSensorThickness(ist); tube.RadLength = stsInterface->GetRadLength(ist); tube.r = 0.; diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index c5e6357083ee21e3a2e6b5d5bdeba752cb43d593..0b74f0cc066cefc9c6186ef2c8a820efc2809676 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -80,10 +80,11 @@ set(SRCS qa/CbmTrackerInputQaTrd.cxx qa/CbmTrackerInputQaTof.cxx + qa/CbmCaInputQaBase.cxx qa/CbmCaInputQaSts.cxx - qa/CbmCaInputQaMuch.cxx - qa/CbmCaInputQaTrd.cxx - qa/CbmCaInputQaTof.cxx + #qa/CbmCaInputQaMuch.cxx + #qa/CbmCaInputQaTrd.cxx + #qa/CbmCaInputQaTof.cxx qa/CbmCaInputQaSetup.cxx qa/CbmCaOutputQa.cxx qa/CbmCaTrackTypeQa.cxx @@ -119,6 +120,7 @@ set(HEADERS #utils/CbmCaIdealHitProducer.h catools/CaToolsMaterialHelper.h qa/CbmCaInputQaBase.h + qa/CbmCaHitQaData.h ) If(CMAKE_CXX_COMPILER_ID MATCHES "Clang") diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index bf0a7c1cdd5f357226fbb5d2b540c6cc263a6fae..01f8f836b13f559560fc35515e34f5634231d9df 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -469,7 +469,7 @@ InitStatus CbmL1::Init() stationInfo.SetStationType(1); // MVD stationInfo.SetTimeInfo(mvdInterface->IsTimeInfoProvided(iSt)); stationInfo.SetFieldStatus(fTrackingMode == L1Algo::TrackingMode::kMcbm ? 0 : 1); - stationInfo.SetZref(mvdInterface->GetZref(iSt)); + stationInfo.SetZref(mvdInterface->GetZ(iSt)); stationInfo.SetZmin(mvdInterface->GetZmin(iSt)); stationInfo.SetZmax(mvdInterface->GetZmax(iSt)); stationInfo.SetXmax(mvdInterface->GetXmax(iSt)); @@ -493,7 +493,7 @@ InitStatus CbmL1::Init() stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false); stationInfo.SetFieldStatus(L1Algo::TrackingMode::kMcbm == fTrackingMode ? 0 : 1); - stationInfo.SetZref(stsInterface->GetZref(iSt)); + stationInfo.SetZref(stsInterface->GetZ(iSt)); stationInfo.SetZmin(stsInterface->GetZmin(iSt)); stationInfo.SetZmax(stsInterface->GetZmax(iSt)); stationInfo.SetXmax(stsInterface->GetXmax(iSt)); @@ -518,7 +518,7 @@ InitStatus CbmL1::Init() stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false); stationInfo.SetFieldStatus(0); - stationInfo.SetZref(muchInterface->GetZref(iSt)); + stationInfo.SetZref(muchInterface->GetZ(iSt)); stationInfo.SetZmin(muchInterface->GetZmin(iSt)); stationInfo.SetZmax(muchInterface->GetZmax(iSt)); stationInfo.SetXmax(muchInterface->GetXmax(iSt)); @@ -543,7 +543,7 @@ InitStatus CbmL1::Init() stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false); stationInfo.SetFieldStatus(0); - stationInfo.SetZref(trdInterface->GetZref(iSt)); + stationInfo.SetZref(trdInterface->GetZ(iSt)); stationInfo.SetZmin(trdInterface->GetZmin(iSt)); stationInfo.SetZmax(trdInterface->GetZmax(iSt)); stationInfo.SetXmax(trdInterface->GetXmax(iSt)); @@ -569,7 +569,7 @@ InitStatus CbmL1::Init() stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false); stationInfo.SetFieldStatus(0); - stationInfo.SetZref(tofInterface->GetZref(iSt)); + stationInfo.SetZref(tofInterface->GetZ(iSt)); stationInfo.SetZmin(tofInterface->GetZmin(iSt)); stationInfo.SetZmax(tofInterface->GetZmax(iSt)); stationInfo.SetXmax(tofInterface->GetXmax(iSt)); diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h index 6890d9bd0806031c4beabe3f3f5afb930c443dbb..792d212511f1f4f7fe4644f7294eb0fef3d40622 100644 --- a/reco/L1/CbmL1DetectorID.h +++ b/reco/L1/CbmL1DetectorID.h @@ -76,6 +76,14 @@ namespace cbm::ca /// read-out corruption. constexpr DetIdArr_t<const char*> kDetName = {{"MVD", "STS", "MUCH", "TRD", "TOF"}}; + /// @brief Name of hit branches for each detector + constexpr DetIdArr_t<const char*> kDetHitBrName = {{"MvdHit", "StsHit", "MuchPixelHit", "TrdHit", "TofHit"}}; + + /// @brief Name of point branches for each detector + constexpr DetIdArr_t<const char*> kDetPointBrName = {{"MvdPoint", "StsPoint", "MuchPoint", "TrdPoint", "TofPoint"}}; + + /// @brief Name + /// @brief Types of MC point objects for each detector using PointTypes_t = DetIdTypeArr_t<CbmMvdPoint, CbmStsPoint, CbmMuchPoint, CbmTrdPoint, CbmTofPoint>; diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h index 3bd3e5527734bc2679fb6dede587e2597f438662..c1dc866aae060386d505547df0245713a05101c7 100644 --- a/reco/L1/L1LinkDef.h +++ b/reco/L1/L1LinkDef.h @@ -24,10 +24,11 @@ //#pragma link C++ class CbmL1SttTrack+; #pragma link C++ class CbmTrackerInputQaTrd + ; #pragma link C++ class CbmTrackerInputQaTof + ; -#pragma link C++ class CbmCaInputQaMuch + ; +#pragma link C++ class CbmCaInputQaBase + ; +//#pragma link C++ class CbmCaInputQaMuch + ; #pragma link C++ class CbmCaInputQaSts + ; -#pragma link C++ class CbmCaInputQaTrd + ; -#pragma link C++ class CbmCaInputQaTof + ; +//#pragma link C++ class CbmCaInputQaTrd + ; +//#pragma link C++ class CbmCaInputQaTof + ; #pragma link C++ enum cbm::ca::ETrackType; #pragma link C++ enum class L1DetectorID; #pragma link C++ class cbm::ca::OutputQa + ; diff --git a/reco/L1/qa/CbmCaInputQaBase.h b/reco/L1/qa/CbmCaInputQaBase.h index 0169b0f68cc2275e83b7ed5df4d124456d194d8a..f775b093ebc485f04bac68a579569ab852965eba 100644 --- a/reco/L1/qa/CbmCaInputQaBase.h +++ b/reco/L1/qa/CbmCaInputQaBase.h @@ -2,44 +2,76 @@ SPDX-License-Identifier: GPL-3.0-only Authors: Sergei Zharko [committer] */ -/// \file CbmCaInputQaBase.h -/// \date 13.01.2023 -/// \brief Class providing basic functionality for CA input QA-tasks -/// \author S.Zharko <s.zharko@lx-pool.gsi.de> +/// @file CbmCaInputQaBase.h +/// @date 13.01.2023 +/// @brief Class providing basic functionality for CA input QA-tasks +/// @author S.Zharko <s.zharko@gsi.de> #ifndef CbmCaInputQaBase_h #define CbmCaInputQaBase_h 1 -#include "Logger.h" +#include "CbmL1DetectorID.h" +#include "CbmCaInputQaBase.h" +#include "CbmCaHitQaData.h" +#include "CbmMCDataManager.h" +#include "CbmQaTask.h" -#include <array> +#include "TMath.h" + +#include <set> +#include <unordered_map> +#include <vector> + +class CbmMatch; +class CbmMCEventList; +class CbmMCDataArray; +class CbmMCDataManager; +class CbmMCTrack; +class CbmStsHit; +class CbmStsPoint; +class CbmTrackingDetectorInterfaceBase; +class CbmTimeSlice; +class FairMCPoint; +class TClonesArray; +class TH1F; +class TH2F; +class TProfile; +class TProfile2D; + + +/// A QA-task class, which provides assurance of MuCh hits and geometry +template <L1DetectorID DetID> +class CbmCaInputQaBase : public CbmQaTask { + using Point_t = cbm::ca::PointTypes_t::at<DetID>; ///< Point type for detector ID + using Hit_t = cbm::ca::HitTypes_t::at<DetID>; ///< Hit type for detector ID -/// Class CbmCaInputQaBase contains common functionality for all of the CA input QA-tasks -/// -class CbmCaInputQaBase { public: - /// Default constructor - CbmCaInputQaBase() = default; + /// @brief Constructor from parameters + /// @param name Name of the task + /// @param verbose Verbose level + /// @param isMCUsed Flag, whether MC information is available for this task + CbmCaInputQaBase(const char* name, int verbose, bool isMCUsed); - /// Default destructor - ~CbmCaInputQaBase() = default; + /// @brief Destructor + virtual ~CbmCaInputQaBase() = default; - /// Sets maximum allowed deviation of pull distribution mean from 0 - /// \param devThrsh Deviation threshold (- devThrsh, + devThrsh) + + /// @brief Sets maximum allowed deviation of pull distribution mean from 0 + /// @param devThrsh Deviation threshold (- devThrsh, + devThrsh) void SetPullMeanDeviationThrsh(double devThrsh) { fPullMeanThrsh = devThrsh; } - /// Sets maximum allowed deviation of pull distribution width from unity - /// \param devThrsh Deviation threshold (1 - devThrsh, 1 + devThrsh) + /// @brief Sets maximum allowed deviation of pull distribution width from unity + /// @param devThrsh Deviation threshold (1 - devThrsh, 1 + devThrsh) void SetPullWidthDeviationThrsh(double devThrsh) { fPullWidthThrsh = devThrsh; } - /// Sets maximum allowed deviation of residual mean from zero (-devThrsh * rms, +devThrsh * rms) - /// \param devThrsh Deviation threshold in sigmas of residual distribution + /// @brief Sets maximum allowed deviation of residual mean from zero (-devThrsh * rms, +devThrsh * rms) + /// @param devThrsh Deviation threshold in sigmas of residual distribution void SetResidualMeanDeviationThrsh(double devThrsh) { fResMeanThrsh = devThrsh; } - /// Sets hit efficiency analysis parameters: efficiency fit range - /// \param thrsh Threshold for efficiency in a given range - /// \param loRange Lower bound of the integration range - /// \param upRange Upper limit of the integration range + /// @brief Sets hit efficiency analysis parameters: efficiency fit range + /// @param thrsh Threshold for efficiency in a given range + /// @param loRange Lower bound of the integration range + /// @param upRange Upper limit of the integration range void SetEfficiencyThrsh(double thrsh, double loRange, double upRange) { LOG(info) << "Call: SetEfficiencyThrsh"; @@ -49,36 +81,209 @@ public: LOG(info) << fEffRange[0] << ' ' << fEffRange[1]; } - /// Sets maximum allowed difference between z-position of hit and station - /// \param dz Difference between station and hit position z components [cm] + /// @brief Sets maximum allowed difference between z-position of hit and station + /// @param dz Difference between station and hit position z components [cm] void SetMaxDiffZStationHit(double dz) { fMaxDiffZStHit = dz; } - /// Sets lower momentum threshold - /// \param mom Minimum absolute value of momentum + /// @brief Sets lower momentum threshold + /// @param mom Minimum absolute value of momentum void SetMinMomentum(double mom) { fMinMomentum = mom; } protected: - // ********************* - // ** Data structures ** - // ********************* + // ******************************************** + // ** Virtual method override from CbmQaTask ** + // ******************************************** + + /// @brief Checks results of the QA and returns some flag + bool Check(); + + /// @brief Initializes data branches + InitStatus InitDataBranches(); + + /// @brief Initializes canvases + InitStatus InitCanvases(); + + /// @brief Initializes histograms + InitStatus InitHistograms(); + + /// @brief Fills histograms per event or time-slice + void FillHistograms(); + + /// @brief Fills histograms per hit + void FillHistogramsPerHit() {} + + /// @brief Fills histograms per MC point + void FillHistogramsPerPoint() {} + + /// @brief De-initializes histograms + void DeInit(); - /// Structure to store fit residuals result + /// @brief Defines parameter for a derived class + void DefineParameters(); + + /// @brief Checks ranges for mean and standard deviation + /// @return False, if variable exits the range + bool CheckRangePull(TH1* h) + { + return CbmQaTask::CheckRange(h, fPullMeanThrsh, 1. - fPullWidthThrsh, 1. + fPullWidthThrsh); + } + + /// @brief Returns a pointer to current hit QA data object + /// + /// The hit QA data object is filled on each step for hit reconstruction and available inside + /// functions FillHit() and FillMCPoint() + const cbm::ca::HitQaData& GetHitQaData() const { return fHitQaData; } + + /// @brief checks if the track at mc point passes the track selection cuts + /// @return true when selected + bool IsTrackSelected(const CbmMCTrack* track, const Point_t* point) const; + + /// @struct ResidualFitResult + /// @brief Stores fit residuals result struct ResidualFitResult { double mean = 0; ///< mean of the distribution double lo = 0; ///< lower limit for the mean double hi = 0; ///< higher limit for the mean }; - // ----- QA constants (thresholds, ranges etc.) + // ****************** + // ** Parameters ** + // ****************** + + static constexpr double kNAN = std::numeric_limits<double>::signaling_NaN(); + //TODO: remove check of residuals - double fResMeanThrsh = .5; ///< Maximum allowed deviation of residual mean from zero [sigma] - double fPullMeanThrsh = .05; ///< Maximum allowed deviation of pull mean from zero - double fPullWidthThrsh = .1; ///< Maximum allowed deviation of pull width from unity - double fEffThrsh = .5; ///< Threshold for hit efficiency in the selected range - double fMaxDiffZStHit = 1.; ///< Maximum allowed difference between z-position of hit and station [cm] - double fMinMomentum = .05; ///< Minimum momentum of particle [GeV/c] - - std::array<double, 2> fEffRange = {0., 100.}; ///< Range for hit efficiency approximation + // ----- QA constants (thresholds, ranges etc.) + double fResMeanThrsh = kNAN; ///< Maximum allowed deviation of residual mean from zero [sigma] + double fPullMeanThrsh = kNAN; ///< Maximum allowed deviation of pull mean from zero + double fPullWidthThrsh = kNAN; ///< Maximum allowed deviation of pull width from unity + double fEffThrsh = kNAN; ///< Threshold for hit efficiency in the selected range + double fMaxDiffZStHit = kNAN; ///< Maximum allowed difference between z-position of hit and station [cm] + double fMinMomentum = kNAN; ///< Minimum momentum of particle [GeV/c] + + std::array<double, 2> fEffRange = {kNAN, kNAN}; ///< Range for hit efficiency approximation + + // ----- Track selection for efficiencies and pulls + bool fTrackSelectionPrimary = true; ///< must the track come from the primary vertex + double fTrackSelectionMomentum = kNAN; ///< track momentum GeV when it exits the tracking station + double fTrackSelectionAngle = kNAN; ///< track incl. angle [grad] when it exits the station + + // ----- Histogram binning parameters + int fNbins = 200; ///< General number of bins + int fNbinsZ = 800; ///< Number of bins for z axis in overall views + + std::array<double, 2> fRHitDx = {kNAN, kNAN}; ///< Range for hit x coordinate error [cm] + std::array<double, 2> fRHitDy = {kNAN, kNAN}; ///< Range for hit y coordinate error [cm] + std::array<double, 2> fRHitDu = {kNAN, kNAN}; ///< Range for hit u coordinate error [cm] + std::array<double, 2> fRHitDv = {kNAN, kNAN}; ///< Range for hit v coordinate error [cm] + std::array<double, 2> fRHitDt = {kNAN, kNAN}; ///< Range for hit time error [ns] + + std::array<double, 2> fRResX = {kNAN, kNAN}; ///< Range for residual of x coordinate [cm] + std::array<double, 2> fRResY = {kNAN, kNAN}; ///< Range for residual of y coordinate [cm] + std::array<double, 2> fRResU = {kNAN, kNAN}; ///< Range for residual of u coordinate [cm] + std::array<double, 2> fRResV = {kNAN, kNAN}; ///< Range for residual of v coordinate [cm] + std::array<double, 2> fRResT = {kNAN, kNAN}; ///< Range for residual of time [ns] + + // NOTE: Pull binning is fixed by convention, since it is used for hit finder calibrations. Please, + // do not modify! + static constexpr int kNbinsPull = 200; + static constexpr double kRPull[2] = {-10., 10.}; ///< Range for pull histograms coordinate + + std::vector<double> frXmin; ///< Range for hit x coordinate [cm] + std::vector<double> frXmax; ///< Range for hit x coordinate [cm] + + std::vector<double> frYmin; ///< Range for hit y coordinate [cm] + std::vector<double> frYmax; ///< Range for hit y coordinate [cm] + + std::vector<double> frZmin; ///< Range for hit z coordinate [cm] + std::vector<double> frZmax; ///< Range for hit z coordinate [cm] + + // ********************* + // ** Data branches ** + // ********************* + + CbmTrackingDetectorInterfaceBase* fpDetInterface = nullptr; ///< Instance of detector interface + + CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice + + TClonesArray* fpHits = nullptr; ///< Array of hits + TClonesArray* fpClusters = nullptr; ///< Array of hit clusters + + CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager + CbmMCEventList* fpMCEventList = nullptr; ///< MC event list + CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks + CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points + + TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches + + + + // ----- Histograms + + + + + // NOTE: the last element of each vector stands for integral distribution over all stations + + // ----- Additional histograms + // Hit occupancy + std::vector<TH2F*> fvph_hit_xy; ///< hit y vs x in different stations + std::vector<TH2F*> fvph_hit_zx; ///< hit x vs z in different stations + std::vector<TH2F*> fvph_hit_zy; ///< hit y vs z in different stations + + std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm] + + // Hit errors + std::vector<TH1F*> fvph_hit_dx; + std::vector<TH1F*> fvph_hit_dy; + std::vector<TH1F*> fvph_hit_du; + std::vector<TH1F*> fvph_hit_dv; + std::vector<TH1F*> fvph_hit_kuv; + std::vector<TH1F*> fvph_hit_dt; + + // MC points occupancy + std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations + + std::vector<TH2F*> fvph_point_xy; ///< point y vs x in different stations + std::vector<TH2F*> fvph_point_zx; ///< point x vs z in different stations + std::vector<TH2F*> fvph_point_zy; ///< point y vs z in different stations + + std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations + + // Residuals + std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations + std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations + std::vector<TH1F*> fvph_res_u; ///< residual for u coordinate in different stations + std::vector<TH1F*> fvph_res_v; ///< residual for v coordinate in different stations + std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations + + std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations + std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations + std::vector<TH2F*> fvph_res_u_vs_u; ///< residual for u coordinate in different stations + std::vector<TH2F*> fvph_res_v_vs_v; ///< residual for v coordinate in different stations + std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations + + // Pulls + std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations + std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations + std::vector<TH1F*> fvph_pull_u; ///< pull for u coordinate in different stations + std::vector<TH1F*> fvph_pull_v; ///< pull for v coordinate in different stations + std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations + + std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations + std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations + std::vector<TH2F*> fvph_pull_u_vs_u; ///< pull for u coordinate in different stations + std::vector<TH2F*> fvph_pull_v_vs_v; ///< pull for v coordinate in different stations + std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations + + // Hit efficiencies + std::vector<TProfile2D*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC + std::vector<TH1F*> fvph_reco_eff; ///< Distribution of hit reconstruction efficiency in bins of fvpe_reco_eff_vs_xy + +// FIXME: change to private +protected: + cbm::ca::HitQaData fHitQaData {}; ///< Current hit QA data object + }; #endif // CbmCaInputQaBase_h diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx index 57711228296b698b31eb5e3b5c6a4f430a91f7e8..d98c2148b7473d67d30be7ce34859cf2afc7ac13 100644 --- a/reco/L1/qa/CbmCaInputQaSts.cxx +++ b/reco/L1/qa/CbmCaInputQaSts.cxx @@ -53,7 +53,8 @@ namespace phys = L1Constants::phys; // from L1Constants.h // --------------------------------------------------------------------------------------------------------------------- // -CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaSts", verbose, isMCUsed) +CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed) +: CbmCaInputQaBase("CbmCaInputQaSts", verbose, isMCUsed) { // TODO: remove an enlarged threshold when the pulls are fixed SetPullMeanDeviationThrsh(0.1); @@ -65,163 +66,9 @@ CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed) : CbmQaTask("CbmCaI // bool CbmCaInputQaSts::Check() { - bool res = true; - - int nSt = fpDetInterface->GetNtrackingStations(); - - - // ************************************************************** - // ** Basic checks, available both for real and simulated data ** - // ************************************************************** - - // Fill efficiency distributions - - for (int iSt = 0; iSt < nSt; ++iSt) { - TProfile2D* effxy = fvpe_reco_eff_vs_xy[iSt]; - for (int i = 1; i < effxy->GetNbinsX() - 1; i++) { - for (int j = 1; j < effxy->GetNbinsY() - 1; j++) { - int bin = effxy->GetBin(i, j); - if (effxy->GetBinEntries(bin) >= 1) { - fvph_reco_eff[iSt]->Fill(effxy->GetBinContent(bin)); - fvph_reco_eff[nSt]->Fill(effxy->GetBinContent(bin)); - } - } - } - } - - // ----- Checks for mismatches in the ordering of the stations - // - std::vector<double> vStationPos(nSt, 0.); - for (int iSt = 0; iSt < nSt; ++iSt) { - vStationPos[iSt] = fpDetInterface->GetZ(iSt); - } - - if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) { - if (fVerbose > 0) { - LOG(error) << fName << ": stations are ordered improperly along the beam axis:"; - for (auto z : vStationPos) { - LOG(error) << "\t- " << z; - } - } - res = false; - } - - // ----- Checks for mismatch between station and hit z positions - // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking - // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of - // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the - // positions of the hits along z-axis are distributed relatively to the position of the station in some range. - // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range - // one should select large enough values to cover the difference described above and in the same time small enough - // to avoid overlaps with the neighboring stations. - // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the - // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to - // another station. - // - for (int iSt = 0; iSt < nSt; ++iSt) { - int nHits = fvph_hit_station_delta_z[iSt]->GetEntries(); - if (!nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits"; - res = false; - continue; - } - int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit); - int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit); - - if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions"; - res = false; - } - } - - - // ******************************************************* - // ** Additional checks, if MC information is available ** - // ******************************************************* - + bool res = CbmCaInputQaBase::Check(); + if (IsMCUsed()) { - // ----- Check efficiencies - // Error is raised, if any station has integrated efficiency lower then a selected threshold. - // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform - // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant - // - // Fit efficiency curves - LOG(info) << "-- Hit efficiency integrated over hit distance from station center"; - - auto* pEffTable = MakeQaObject<CbmQaTable>("vs Station/eff_table", "Efficiency table", nSt, 2); - pEffTable->SetNamesOfCols({"Station ID", "Efficiency"}); - pEffTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt; ++iSt) { - auto eff = fvph_reco_eff[iSt]->GetMean(); - pEffTable->SetCell(iSt, 0, iSt); - pEffTable->SetCell(iSt, 1, eff); - res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000); - } - - LOG(info) << '\n' << pEffTable->ToString(3); - - // ----- Checks for residuals - // Check hits for potential biases from central values - - auto* pResidualsTable = - MakeQaObject<CbmQaTable>("vs Station/residuals_mean", "Residual mean values in different stations", nSt, 4); - pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"}); - pResidualsTable->SetColWidth(20); - - // Fit residuals - for (int iSt = 0; iSt <= nSt; ++iSt) { - - cbm::qa::util::SetLargeStats(fvph_res_x[iSt]); - cbm::qa::util::SetLargeStats(fvph_res_y[iSt]); - cbm::qa::util::SetLargeStats(fvph_res_t[iSt]); - cbm::qa::util::SetLargeStats(fvph_res_u[iSt]); - cbm::qa::util::SetLargeStats(fvph_res_v[iSt]); - - pResidualsTable->SetCell(iSt, 0, iSt); - pResidualsTable->SetCell(iSt, 1, fvph_res_x[iSt]->GetStdDev()); - pResidualsTable->SetCell(iSt, 2, fvph_res_y[iSt]->GetStdDev()); - pResidualsTable->SetCell(iSt, 3, fvph_res_t[iSt]->GetStdDev()); - } - - LOG(info) << '\n' << pResidualsTable->ToString(8); - - // ----- Checks for pulls - // - // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component - // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case, - // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal - // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out - // this range, QA task fails. - - auto* pPullsTable = - MakeQaObject<CbmQaTable>("vs Station/pulls_rms", "Pulls std. dev. values in different stations", nSt, 4); - pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(t) sigm"}); - pPullsTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt + 1; ++iSt) { - - // Fit pull distributions for nicer representation. Fit results are not used in further checks. - - cbm::qa::util::SetLargeStats(fvph_pull_x[iSt]); - cbm::qa::util::SetLargeStats(fvph_pull_y[iSt]); - cbm::qa::util::SetLargeStats(fvph_pull_t[iSt]); - cbm::qa::util::SetLargeStats(fvph_pull_u[iSt]); - cbm::qa::util::SetLargeStats(fvph_pull_v[iSt]); - - // Check the pull quality - res = CheckRangePull(fvph_pull_x[iSt]) && res; - res = CheckRangePull(fvph_pull_y[iSt]) && res; - res = CheckRangePull(fvph_pull_t[iSt]) && res; - res = CheckRangePull(fvph_pull_u[iSt]) && res; - res = CheckRangePull(fvph_pull_v[iSt]) && res; - - pPullsTable->SetCell(iSt, 0, iSt); - pPullsTable->SetCell(iSt, 1, fvph_pull_x[iSt]->GetStdDev()); - pPullsTable->SetCell(iSt, 2, fvph_pull_y[iSt]->GetStdDev()); - pPullsTable->SetCell(iSt, 3, fvph_pull_t[iSt]->GetStdDev()); - } - for (int idig = 0; idig <= fkMaxDigisInClusterForPulls; idig++) { cbm::qa::util::SetLargeStats(fvph_pull_u_Ndig[idig]); res = CheckRangePull(fvph_pull_u_Ndig[idig]) && res; @@ -229,8 +76,6 @@ bool CbmCaInputQaSts::Check() cbm::qa::util::SetLargeStats(fvph_pull_v_Ndig[idig]); res = CheckRangePull(fvph_pull_v_Ndig[idig]) && res; } - - LOG(info) << '\n' << pPullsTable->ToString(3); } // McUsed return res; @@ -240,413 +85,72 @@ bool CbmCaInputQaSts::Check() // void CbmCaInputQaSts::DeInit() { - // Vectors with pointers to histograms - fvph_hit_xy.clear(); - fvph_hit_zx.clear(); - fvph_hit_zy.clear(); - - fvph_hit_station_delta_z.clear(); - - fvph_hit_dx.clear(); - fvph_hit_dy.clear(); - fvph_hit_du.clear(); - fvph_hit_dv.clear(); - fvph_hit_kuv.clear(); - fvph_hit_dt.clear(); - - fvph_n_points_per_hit.clear(); - - fvph_point_xy.clear(); - fvph_point_zx.clear(); - fvph_point_zy.clear(); - - fvph_point_hit_delta_z.clear(); - - fvph_res_x.clear(); - fvph_res_y.clear(); - fvph_res_u.clear(); - fvph_res_v.clear(); - fvph_res_t.clear(); - - fvph_pull_x.clear(); - fvph_pull_y.clear(); - fvph_pull_u.clear(); - fvph_pull_v.clear(); - fvph_pull_t.clear(); + CbmCaInputQaBase::DeInit(); + // Vectors with pointers to histograms fvph_pull_u_Ndig.clear(); fvph_pull_v_Ndig.clear(); - - fvph_res_x_vs_x.clear(); - fvph_res_y_vs_y.clear(); - fvph_res_u_vs_u.clear(); - fvph_res_v_vs_v.clear(); - fvph_res_t_vs_t.clear(); - - fvph_pull_x_vs_x.clear(); - fvph_pull_y_vs_y.clear(); - fvph_pull_u_vs_u.clear(); - fvph_pull_v_vs_v.clear(); - fvph_pull_t_vs_t.clear(); - - fvpe_reco_eff_vs_xy.clear(); - fvph_reco_eff.clear(); } // --------------------------------------------------------------------------------------------------------------------- // -void CbmCaInputQaSts::FillHistograms() +void CbmCaInputQaSts::DefineParameters() { - int nSt = fpDetInterface->GetNtrackingStations(); - int nHits = fpHits->GetEntriesFast(); - int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1; + auto SetRange = [] (std::array<double, 2>& range, double min, double max) { range[0] = min; range[1] = max; }; + // Hit errors + SetRange(fRHitDx, 0.0000, 0.0050); // [cm] + SetRange(fRHitDy, 0.0000, 0.0200); // [cm] + SetRange(fRHitDu, 0.0000, 0.0050); // [cm] + SetRange(fRHitDv, 0.0000, 0.0050); // [cm] + SetRange(fRHitDt, 0.0000, 10.000); // [ns] + // Residuals + SetRange(fRResX, -0.02, 0.02); + SetRange(fRResY, -0.10, 0.10); + SetRange(fRResU, -0.02, 0.02); + SetRange(fRResV, -0.02, 0.02); + SetRange(fRResT, -25.0, 25.0); + // QA result selection criteria + SetRange(fEffRange, 0.0, 100.0); ///< Range for hit efficiency approximation + fResMeanThrsh = 0.50; ///< Maximum allowed deviation of residual mean from zero [sigma] + fPullMeanThrsh = 0.05; ///< Maximum allowed deviation of pull mean from zero + fPullWidthThrsh = 0.10; ///< Maximum allowed deviation of pull width from unity + fEffThrsh = 0.50; ///< Threshold for hit efficiency in the selected range + fMaxDiffZStHit = 1.00; ///< Maximum allowed difference between z-position of hit and station [cm] + fMinMomentum = 0.05; ///< Minimum momentum of particle [GeV/c] + // Track selection criteria + fTrackSelectionPrimary = true; ///< must the track come from the primary vertex + fTrackSelectionMomentum = 0.1; ///< track momentum GeV when it exits the tracking station + fTrackSelectionAngle = 60.; ///< track incl. angle [grad] when it exits the station - std::vector<std::vector<std::vector<int>>> - vNofHitsPerMcTrack; // Number of hits per MC track per station in different MC events +} - if (IsMCUsed()) { - vNofHitsPerMcTrack.resize(nMCevents); - for (int iE = 0; iE < nMCevents; ++iE) { - int iFile = fpMCEventList->GetFileIdByIndex(iE); - int iEvent = fpMCEventList->GetEventIdByIndex(iE); - int nMcTracks = fpMCTracks->Size(iFile, iEvent); - vNofHitsPerMcTrack[iE].resize(nSt); - for (int iSt = 0; iSt < nSt; iSt++) { - vNofHitsPerMcTrack[iE][iSt].resize(nMcTracks, 0); - } - } +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaInputQaSts::FillHistogramsPerHit() +{ + const auto* pHit = dynamic_cast<const CbmStsHit*>(fpHits->At(fHitQaData.GetHitIndex())); + + { // u-coordinate residual per number of digis + const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetFrontClusterId())); + assert(pCluster); + int nDigis = pCluster->GetNofDigis(); + if (nDigis > fkMaxDigisInClusterForPulls) { nDigis = 0; } + fvph_pull_u_Ndig[nDigis]->Fill(fHitQaData.GetPullU()); + } + + { // v-coordinate residual per number of digis + const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetBackClusterId())); + assert(pCluster); + int nDigis = pCluster->GetNofDigis(); + if (nDigis > fkMaxDigisInClusterForPulls) { nDigis = 0; } + fvph_pull_v_Ndig[nDigis]->Fill(fHitQaData.GetPullV()); } - - for (int iH = 0; iH < nHits; ++iH) { - const auto* pHit = dynamic_cast<const CbmStsHit*>(fpHits->At(iH)); - LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmStsHit (dynamic cast failed)"; - - // ************************* - // ** Reconstructed hit QA - - // Hit station geometry info - int iSt = fpDetInterface->GetTrackingStationIndex(pHit); - LOG_IF(fatal, iSt < 0 || iSt >= nSt) << fName << ": index of station (" << iSt << ") is out of range " - << "for hit with id = " << iH; - double stPhiU; // STS front strips stereo angle - double stPhiV; // STS back strips stereo angle - - std::tie(stPhiU, stPhiV) = fpDetInterface->GetStereoAnglesSensor(pHit->GetAddress()); - - double sinU = sin(stPhiU); - double cosU = cos(stPhiU); - double sinV = sin(stPhiV); - double cosV = cos(stPhiV); - - // hit cov matrix - double Cxx = pHit->GetDx() * pHit->GetDx(); - double Cxy = pHit->GetDxy(); - double Cyy = pHit->GetDy() * pHit->GetDy(); - - // Hit position - double xHit = pHit->GetX(); - double yHit = pHit->GetY(); - double zHit = pHit->GetZ(); - double uHit = xHit * cosU + yHit * sinU; - double vHit = xHit * cosV + yHit * sinV; - - // Hit errors - //double duHit = pHit->GetDu(); - //double dvHit = pHit->GetDv(); - double duHit = sqrt(cosU * cosU * Cxx + 2. * sinU * cosU * Cxy + sinU * sinU * Cyy); - double dvHit = sqrt(cosV * cosV * Cxx + 2. * sinV * cosV * Cxy + sinV * sinV * Cyy); - double dxHit = pHit->GetDx(); - double dyHit = pHit->GetDy(); - - double duvHit = Cxx * cosU * cosV + Cxy * (sinU * cosV + cosU * sinV) + Cyy * sinU * sinV; - double kuvHit = duvHit / duHit / dvHit; - - // Hit time - double tHit = pHit->GetTime(); - double dtHit = pHit->GetTimeError(); - - fvph_hit_xy[iSt]->Fill(xHit, yHit); - fvph_hit_zx[iSt]->Fill(zHit, xHit); - fvph_hit_zy[iSt]->Fill(zHit, yHit); - - fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZref(iSt)); - - fvph_hit_dx[iSt]->Fill(dxHit); - fvph_hit_dy[iSt]->Fill(dyHit); - fvph_hit_du[iSt]->Fill(duHit); - fvph_hit_dv[iSt]->Fill(dvHit); - fvph_hit_kuv[iSt]->Fill(kuvHit); - fvph_hit_dt[iSt]->Fill(dtHit); - - fvph_hit_xy[nSt]->Fill(xHit, yHit); - fvph_hit_zx[nSt]->Fill(zHit, xHit); - fvph_hit_zy[nSt]->Fill(zHit, yHit); - - fvph_hit_station_delta_z[nSt]->Fill(zHit - fpDetInterface->GetZref(iSt)); - - fvph_hit_dx[nSt]->Fill(dxHit); - fvph_hit_dy[nSt]->Fill(dyHit); - fvph_hit_du[nSt]->Fill(duHit); - fvph_hit_dv[nSt]->Fill(dvHit); - fvph_hit_kuv[nSt]->Fill(kuvHit); - fvph_hit_dt[nSt]->Fill(dtHit); - - // ********************** - // ** MC information QA - - if (IsMCUsed()) { - const auto* pHitMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH)); - assert(pHitMatch); - - // Evaluate number of hits per one MC point - int nMCpoints = 0; // Number of MC points for this hit - for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) { - const auto& link = pHitMatch->GetLink(iLink); - - int iP = link.GetIndex(); // Index of MC point - - // Skip noisy links - if (iP < 0) { continue; } - - ++nMCpoints; - - int iE = fpMCEventList->GetEventIndex(link); - - LOG_IF(fatal, iE < 0 || iE >= nMCevents) << fName << ": id of MC event is out of range (hit id = " << iH - << ", link id = " << iLink << ", event id = " << iE << ')'; - - // matched point - const auto* pMCPoint = dynamic_cast<const CbmStsPoint*>(fpMCPoints->Get(link)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH; - - int iTr = pMCPoint->GetTrackID(); - - LOG_IF(fatal, iTr >= (int) vNofHitsPerMcTrack[iE][iSt].size()) - << fName << ": index of MC track is out of range (hit id = " << iH << ", link id = " << iLink - << ", event id = " << iE << ", track id = " << iTr << ')'; - - vNofHitsPerMcTrack[iE][iSt][iTr]++; - } - - fvph_n_points_per_hit[iSt]->Fill(nMCpoints); - fvph_n_points_per_hit[nSt]->Fill(nMCpoints); - - // - // Fill the following histograms exclusively for isolated hits with a one-to-one correspondence to the mc track - // The mc track must satisfy the cuts - // - - if (pHitMatch->GetNofLinks() != 1) { continue; } - - // The best link to in the match (probably, the cut on nMCpoints is meaningless) - assert(pHitMatch->GetNofLinks() > 0); // Should be always true due to the cut above - const auto& bestPointLink = pHitMatch->GetMatchedLink(); - - // Skip noisy links - if (bestPointLink.GetIndex() < 0) { continue; } - - // Point matched by the best link - const auto* pMCPoint = dynamic_cast<const CbmStsPoint*>(fpMCPoints->Get(bestPointLink)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH; - - // MC track for this point - CbmLink bestTrackLink = bestPointLink; - bestTrackLink.SetIndex(pMCPoint->GetTrackID()); - const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(bestTrackLink)); - LOG_IF(fatal, !pMCTrack) << fName << ": MC track object does not exist for hit " << iH << " and link: "; - - double t0MC = fpMCEventList->GetEventTime(bestPointLink); - LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC; - - // cut on the mc track quality - if (!IsTrackSelected(pMCTrack, pMCPoint)) { continue; } - - { // skip the case when the mc point participates to several hits - int iE = fpMCEventList->GetEventIndex(bestTrackLink); - if (vNofHitsPerMcTrack[iE][iSt][pMCPoint->GetTrackID()] != 1) { continue; } - } - - // ----- MC point properties - // - double mass = pMCTrack->GetMass(); - //double charge = pMCTrack->GetCharge(); - //double pdg = pMCTrack->GetPdgCode(); - - // Entrance position and time - double xMC = pMCPoint->GetXIn(); - double yMC = pMCPoint->GetYIn(); - double zMC = pMCPoint->GetZIn(); - double tMC = pMCPoint->GetTime() + t0MC; - - // MC point entrance momenta - double pxMC = pMCPoint->GetPx(); - double pyMC = pMCPoint->GetPy(); - double pzMC = pMCPoint->GetPz(); - double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC); - - // MC point exit momenta - // double pxMCo = pMCPoint->GetPxOut(); - // double pyMCo = pMCPoint->GetPyOut(); - // double pzMCo = pMCPoint->GetPzOut(); - // double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo); - - // Position and time shifted to z-coordinate of the hit - double shiftZ = zHit - zMC; // Difference between hit and point z positions - double xMCs = xMC + shiftZ * pxMC / pzMC; - double yMCs = yMC + shiftZ * pyMC / pzMC; - double uMCs = xMCs * cosU + yMCs * sinU; - double vMCs = xMCs * cosV + yMCs * sinV; - double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC); - - // Residuals - double xRes = xHit - xMCs; - double yRes = yHit - yMCs; - double uRes = uHit - uMCs; - double vRes = vHit - vMCs; - double tRes = tHit - tMCs; - double zRes = zHit - 0.5 * (pMCPoint->GetZIn() + pMCPoint->GetZOut()); - - fvph_point_hit_delta_z[iSt]->Fill(zRes); - - fvph_res_x[iSt]->Fill(xRes); - fvph_res_y[iSt]->Fill(yRes); - fvph_res_u[iSt]->Fill(uRes); - fvph_res_v[iSt]->Fill(vRes); - fvph_res_t[iSt]->Fill(tRes); - - fvph_pull_x[iSt]->Fill(xRes / dxHit); - fvph_pull_y[iSt]->Fill(yRes / dyHit); - fvph_pull_u[iSt]->Fill(uRes / duHit); - fvph_pull_v[iSt]->Fill(vRes / dvHit); - fvph_pull_t[iSt]->Fill(tRes / dtHit); - - fvph_res_x_vs_x[iSt]->Fill(xHit, xRes); - fvph_res_y_vs_y[iSt]->Fill(yHit, yRes); - fvph_res_u_vs_u[iSt]->Fill(uHit, uRes); - fvph_res_v_vs_v[iSt]->Fill(vHit, vRes); - fvph_res_t_vs_t[iSt]->Fill(tHit, tRes); - - fvph_pull_x_vs_x[iSt]->Fill(xMCs, xRes / dxHit); - fvph_pull_y_vs_y[iSt]->Fill(yMCs, yRes / dyHit); - fvph_pull_u_vs_u[iSt]->Fill(uMCs, uRes / duHit); - fvph_pull_v_vs_v[iSt]->Fill(vMCs, vRes / dvHit); - fvph_pull_t_vs_t[iSt]->Fill(tMCs, tRes / dtHit); - - // fill summary histos - - fvph_point_hit_delta_z[nSt]->Fill(zRes); - - fvph_res_x[nSt]->Fill(xRes); - fvph_res_y[nSt]->Fill(yRes); - fvph_res_u[nSt]->Fill(uRes); - fvph_res_v[nSt]->Fill(vRes); - fvph_res_t[nSt]->Fill(tRes); - - fvph_pull_x[nSt]->Fill(xRes / dxHit); - fvph_pull_y[nSt]->Fill(yRes / dyHit); - fvph_pull_u[nSt]->Fill(uRes / duHit); - fvph_pull_v[nSt]->Fill(vRes / dvHit); - fvph_pull_t[nSt]->Fill(tRes / dtHit); - - fvph_res_x_vs_x[nSt]->Fill(xMCs, xRes); - fvph_res_y_vs_y[nSt]->Fill(yMCs, yRes); - fvph_res_u_vs_u[nSt]->Fill(uMCs, uRes); - fvph_res_v_vs_v[nSt]->Fill(vMCs, vRes); - fvph_res_t_vs_t[nSt]->Fill(tMCs, tRes); - - fvph_pull_x_vs_x[nSt]->Fill(xMCs, xRes / dxHit); - fvph_pull_y_vs_y[nSt]->Fill(yMCs, yRes / dyHit); - fvph_pull_u_vs_u[nSt]->Fill(uMCs, uRes / duHit); - fvph_pull_v_vs_v[nSt]->Fill(vMCs, vRes / dvHit); - fvph_pull_t_vs_t[nSt]->Fill(tMCs, tRes / dtHit); - - { - const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetFrontClusterId())); - assert(pCluster); - int ndigis = pCluster->GetNofDigis(); - if (ndigis > fkMaxDigisInClusterForPulls) { ndigis = 0; } - fvph_pull_u_Ndig[ndigis]->Fill(uRes / duHit); - } - - { - const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetBackClusterId())); - assert(pCluster); - int ndigis = pCluster->GetNofDigis(); - if (ndigis > fkMaxDigisInClusterForPulls) { ndigis = 0; } - fvph_pull_v_Ndig[ndigis]->Fill(vRes / dvHit); - } - } - } // loop over hits: end - - // Fill hit reconstruction efficiencies - if (IsMCUsed()) { - for (int iE = 0; iE < nMCevents; ++iE) { - int iFile = fpMCEventList->GetFileIdByIndex(iE); - int iEvent = fpMCEventList->GetEventIdByIndex(iE); - int nPoints = fpMCPoints->Size(iFile, iEvent); - int nTracks = fpMCTracks->Size(iFile, iEvent); - - // If efficiency for the track at the station is already evaluated - std::vector<std::vector<bool>> vIsTrackProcessed(nSt); - for (int iSt = 0; iSt < nSt; iSt++) { - vIsTrackProcessed[iSt].resize(nTracks, 0); - } - - for (int iP = 0; iP < nPoints; ++iP) { - const auto* pMCPoint = dynamic_cast<const CbmStsPoint*>(fpMCPoints->Get(iFile, iEvent, iP)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile - << ", iEvent = " << iEvent << ", iP = " << iP; - - int address = pMCPoint->GetDetectorID(); - int iSt = fpDetInterface->GetTrackingStationIndex(address); - LOG_IF(fatal, iSt < 0 || iSt >= nSt) - << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address " << address - << " has wrong station index: iSt = " << iSt; - - - double xMC = pMCPoint->GetXIn(); - double yMC = pMCPoint->GetYIn(); - double zMC = pMCPoint->GetZIn(); - - { - fvph_point_xy[iSt]->Fill(xMC, yMC); - fvph_point_zx[iSt]->Fill(zMC, xMC); - fvph_point_zy[iSt]->Fill(zMC, yMC); - - fvph_point_xy[nSt]->Fill(xMC, yMC); - fvph_point_zx[nSt]->Fill(zMC, xMC); - fvph_point_zy[nSt]->Fill(zMC, yMC); - } - - int iTr = pMCPoint->GetTrackID(); - - LOG_IF(fatal, iTr >= nTracks) << fName << ": index of MC track is out of range (point id = " << iP - << ", event id = " << iE << ", track id = " << iTr << ')'; - - const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTr)); - - LOG_IF(fatal, pMCTrack == nullptr) << fName << ": null MC track pointer for file id = " << iFile - << ", event id = " << iEvent << ", track id = " << iTr; - - // check efficiency only once per mc track per station - if (vIsTrackProcessed[iSt][iTr]) { continue; } - - vIsTrackProcessed[iSt][iTr] = true; - - // cut on the mc track quality - if (!IsTrackSelected(pMCTrack, pMCPoint)) { continue; } - - // Conditions under which point is accounted as reconstructed: point - bool ifTrackHasHits = (vNofHitsPerMcTrack[iE][iSt][iTr] > 0); - - fvpe_reco_eff_vs_xy[iSt]->Fill(xMC, yMC, ifTrackHasHits); - fvpe_reco_eff_vs_xy[nSt]->Fill(xMC, yMC, ifTrackHasHits); - - } // loop over MC-points: end - } // loop over MC-events: end - } // MC is used: end } +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaInputQaSts::FillHistogramsPerPoint() {} + // --------------------------------------------------------------------------------------------------------------------- // InitStatus CbmCaInputQaSts::InitDataBranches() @@ -654,107 +158,13 @@ InitStatus CbmCaInputQaSts::InitDataBranches() // STS tracking detector interface fpDetInterface = CbmStsTrackingInterface::Instance(); - LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m"; - - // FairRootManager - auto* pFairRootManager = FairRootManager::Instance(); - LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m"; - - // Time-slice - fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice.")); - LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m"; - - // ----- Reconstructed data-branches initialization - - // Hits container - fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("StsHit")); - LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in STS is not found\033[0m"; - + auto baseInitStatus = CbmCaInputQaBase::InitDataBranches(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + // Clusters container - fpClusters = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("StsCluster")); + fpClusters = dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsCluster")); LOG_IF(fatal, !fpClusters) << "\033[1;31m" << fName << ": container of hit clusters in STS is not found\033[0m"; - - // ----- MC-information branches initialization - if (IsMCUsed()) { - // MC manager - fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager")); - LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m"; - - // MC event list - fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList.")); - LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m"; - - // MC tracks - fpMCTracks = fpMCDataManager->InitBranch("MCTrack"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m"; - - // MC points - fpMCPoints = fpMCDataManager->InitBranch("StsPoint"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for STS\033[0m"; - - // Hit matches - fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("StsHitMatch")); - LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for STS\033[0m"; - } - - - // get hit distribution ranges from the geometry - - int nSt = fpDetInterface->GetNtrackingStations(); - - frXmin.resize(nSt + 1, 0.); - frXmax.resize(nSt + 1, 0.); - - frYmin.resize(nSt + 1, 0.); - frYmax.resize(nSt + 1, 0.); - - frZmin.resize(nSt + 1, 0.); - frZmax.resize(nSt + 1, 0.); - - for (int i = nSt; i >= 0; --i) { - { // read the ranges for station i - int j = (i == nSt ? 0 : i); - - frXmin[i] = fpDetInterface->GetXmin(j); - frXmax[i] = fpDetInterface->GetXmax(j); - - frYmin[i] = fpDetInterface->GetYmin(j); - frYmax[i] = fpDetInterface->GetYmax(j); - - frZmin[i] = fpDetInterface->GetZmin(j); - frZmax[i] = fpDetInterface->GetZmax(j); - } - - // update overall ranges - - frXmin[nSt] = std::min(frXmin[nSt], frXmin[i]); - frXmax[nSt] = std::max(frXmax[nSt], frXmax[i]); - - frYmin[nSt] = std::min(frYmin[nSt], frYmin[i]); - frYmax[nSt] = std::max(frYmax[nSt], frYmax[i]); - - frZmin[nSt] = std::min(frZmin[nSt], frZmin[i]); - frZmax[nSt] = std::max(frZmax[nSt], frZmax[i]); - } - - // add 5% margins for a better view - - for (int i = 0; i <= nSt; ++i) { - double dx = 0.05 * fabs(frXmax[i] - frXmin[i]); - frXmin[i] -= dx; - frXmax[i] += dx; - - double dy = 0.05 * fabs(frYmax[i] - frYmin[i]); - frYmin[i] -= dy; - frYmax[i] += dy; - - double dz = 0.05 * fabs(frZmax[i] - frZmin[i]); - frZmin[i] -= dz; - frZmax[i] += dz; - } - - return kSUCCESS; } @@ -764,344 +174,12 @@ InitStatus CbmCaInputQaSts::InitDataBranches() // InitStatus CbmCaInputQaSts::InitCanvases() { - gStyle->SetOptFit(1); - int nSt = fpDetInterface->GetNtrackingStations(); - - // ******************** - // ** Drawing options - - // Contours - constexpr auto contColor = kOrange + 7; - constexpr auto contWidth = 2; // Line width - constexpr auto contStyle = 2; // Line style - constexpr auto contFill = 0; // Fill style - - // ********************** - // ** summary for all stations - - { - auto* canv = MakeQaObject<CbmQaCanvas>("occ_hit", "Hit Occupancy", 1600, 800); - canv->Divide2D(3); - canv->cd(1); - fvph_hit_xy[nSt]->DrawCopy("colz", ""); - canv->cd(2); - fvph_hit_zx[nSt]->DrawCopy("colz", ""); - canv->cd(3); - fvph_hit_zy[nSt]->DrawCopy("colz", ""); - } - - { - auto* canv = MakeQaObject<CbmQaCanvas>("occ_point", "Point Occupancy", 1600, 800); - canv->Divide2D(3); - - canv->cd(1); - fvph_point_xy[nSt]->DrawCopy("colz", ""); - canv->cd(2); - fvph_point_zx[nSt]->DrawCopy("colz", ""); - canv->cd(3); - fvph_point_zy[nSt]->DrawCopy("colz", ""); - } - - { - auto* canv = MakeQaObject<CbmQaCanvas>("residual", "Hit Residuals", 1600, 800); - canv->Divide2D(6); - canv->cd(1); - fvph_res_x[nSt]->DrawCopy("colz", ""); - canv->cd(2); - fvph_res_y[nSt]->DrawCopy("colz", ""); - canv->cd(3); - fvph_res_t[nSt]->DrawCopy("colz", ""); - canv->cd(4); - fvph_res_u[nSt]->DrawCopy("colz", ""); - canv->cd(5); - fvph_res_v[nSt]->DrawCopy("colz", ""); - } - - { - auto* canv = MakeQaObject<CbmQaCanvas>("pull", "Hit Pulls", 1600, 800); - canv->Divide2D(6); - canv->cd(1); - fvph_pull_x[nSt]->DrawCopy("colz", ""); - canv->cd(2); - fvph_pull_y[nSt]->DrawCopy("colz", ""); - canv->cd(3); - fvph_pull_t[nSt]->DrawCopy("colz", ""); - canv->cd(4); - fvph_pull_u[nSt]->DrawCopy("colz", ""); - canv->cd(5); - fvph_pull_v[nSt]->DrawCopy("colz", ""); - } - - { - auto* canv = MakeQaObject<CbmQaCanvas>("eff", "Hit Reconstruction Efficiency", 1600, 800); - canv->Divide2D(2); - canv->cd(1); - fvpe_reco_eff_vs_xy[nSt]->DrawCopy("colz", ""); - canv->cd(2); - fvph_reco_eff[nSt]->DrawCopy("colz", ""); - } - - { - auto* canv = MakeQaObject<CbmQaCanvas>("err", "Hit Errors", 1600, 800); - canv->Divide2D(6); - canv->cd(1); - fvph_hit_dx[nSt]->DrawCopy(); - canv->cd(2); - fvph_hit_dy[nSt]->DrawCopy(); - canv->cd(3); - fvph_hit_dt[nSt]->DrawCopy(); - canv->cd(4); - fvph_hit_du[nSt]->DrawCopy(); - canv->cd(5); - fvph_hit_dv[nSt]->DrawCopy(); - canv->cd(6); - fvph_hit_kuv[nSt]->DrawCopy(); - } - - { - auto* canv = MakeQaObject<CbmQaCanvas>("other", "Other histograms", 1600, 800); - canv->Divide2D(3); - canv->cd(1); - fvph_n_points_per_hit[nSt]->DrawCopy("colz", ""); - canv->cd(2); - fvph_hit_station_delta_z[nSt]->DrawCopy("colz", ""); - canv->cd(3); - fvph_point_hit_delta_z[nSt]->DrawCopy("colz", ""); - } - - - // ********************************* - // ** Hit occupancies - - { // ** Occupancy in XY plane - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_xy", "", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_hit_xy[iSt]->DrawCopy("colz", ""); - - // Square boarders of the active area of the station - double stXmin = fpDetInterface->GetXmin(iSt); - double stXmax = fpDetInterface->GetXmax(iSt); - double stYmin = fpDetInterface->GetYmin(iSt); - double stYmax = fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - } + auto baseInitStatus = CbmCaInputQaBase::InitCanvases(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } - { // ** Occupancy in XZ plane - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_zx", "", 1600, 800); - canv->cd(); - fvph_hit_zx[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = fpDetInterface->GetXmin(iSt); - double stHmax = fpDetInterface->GetXmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - } - - { // ** Occupancy in YZ plane - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_zy", "", 1600, 800); - canv->cd(); - fvph_hit_zy[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = fpDetInterface->GetYmin(iSt); - double stHmax = fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - } - - // ************ - // ************ + gStyle->SetOptFit(1); if (IsMCUsed()) { - - // ********************** - // ** Point occupancies - - { // ** Occupancy in XY plane - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_xy", "", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_point_xy[iSt]->DrawCopy("colz", ""); - - // Square boarders of the active area of the station - double stXmin = fpDetInterface->GetXmin(iSt); - double stXmax = fpDetInterface->GetXmax(iSt); - double stYmin = fpDetInterface->GetYmin(iSt); - double stYmax = fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - } - - { // ** Occupancy in XZ plane - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_zx", "", 1600, 800); - canv->cd(); - fvph_point_zx[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = fpDetInterface->GetXmin(iSt); - double stHmax = fpDetInterface->GetXmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - } - - { // ** Occupancy in YZ plane - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_zy", "", 1600, 800); - canv->cd(); - fvph_point_zy[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = fpDetInterface->GetYmin(iSt); - double stHmax = fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - } - - // ************** - // ** Residuals - - { // x-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_x", "Residuals for x coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_res_x[iSt]->DrawCopy("", ""); - } - } - - { // y-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_y", "Residuals for y coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_res_y[iSt]->DrawCopy("", ""); - } - } - - { // u-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_u", "Residuals for u coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_res_u[iSt]->DrawCopy("", ""); - } - } - - { // v-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_v", "Residuals for v coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_res_v[iSt]->DrawCopy("", ""); - } - } - - { // t-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_t", "Residuals for time", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_res_t[iSt]->DrawCopy("", ""); - } - } - - - // ********** - // ** Pulls - - { // x-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_x", "Pulls for x coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_pull_x[iSt]->DrawCopy("", ""); - } - } - - { // y-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_y", "Pulls for y coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_pull_y[iSt]->DrawCopy("", ""); - } - } - - { // u-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_u", "Pulls for u coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_pull_u[iSt]->DrawCopy("", ""); - } - } - - { // v-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_v", "Pulls for v coordinate", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_pull_v[iSt]->DrawCopy("", ""); - } - } - - { // t-coordinate - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_t", "Pulls for time", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvph_pull_t[iSt]->DrawCopy("", ""); - } - } - { // u-coordinate vs N digis auto* canv = MakeQaObject<CbmQaCanvas>("vs N digi/pull_u", "Pulls for u coordinate vs N digis", 1600, 800); canv->Divide2D(fkMaxDigisInClusterForPulls + 1); @@ -1119,28 +197,6 @@ InitStatus CbmCaInputQaSts::InitCanvases() fvph_pull_v_Ndig[i]->DrawCopy("", ""); } } - - // ************************************ - // ** Hit reconstruction efficiencies - - { - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/reco_eff", "Hit efficiencies in xy bins", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - //fvph_reco_eff[iSt]->SetMarkerStyle(20); - fvph_reco_eff[iSt]->DrawCopy(); - } - } - - { - auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800); - canv->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - canv->cd(iSt + 1); - fvpe_reco_eff_vs_xy[iSt]->DrawCopy("colz", ""); - } - } } return kSUCCESS; @@ -1150,265 +206,16 @@ InitStatus CbmCaInputQaSts::InitCanvases() // InitStatus CbmCaInputQaSts::InitHistograms() { - int nSt = fpDetInterface->GetNtrackingStations(); - - //SetStoringMode(EStoringMode::kSAMEDIR); - SetStoringMode(EStoringMode::kSUBDIR); - - // ----- Subdirectories - MakeQaDirectory("Summary"); - MakeQaDirectory("Summary/vs Station"); - MakeQaDirectory("Summary/vs N digi"); - MakeQaDirectory("All stations"); - - // ----- Histograms initialization - fvph_hit_xy.resize(nSt + 1, nullptr); - fvph_hit_zy.resize(nSt + 1, nullptr); - fvph_hit_zx.resize(nSt + 1, nullptr); - - fvph_hit_station_delta_z.resize(nSt + 1, nullptr); - - fvph_hit_dx.resize(nSt + 1, nullptr); - fvph_hit_dy.resize(nSt + 1, nullptr); - fvph_hit_dt.resize(nSt + 1, nullptr); - fvph_hit_dv.resize(nSt + 1, nullptr); - fvph_hit_kuv.resize(nSt + 1, nullptr); - fvph_hit_du.resize(nSt + 1, nullptr); - - for (int iSt = 0; iSt <= nSt; ++iSt) { - TString sF = ""; // Subfolder - sF += (iSt == nSt) ? "All stations/" : Form("Station %d/", iSt); - TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix - TString tsuff = (iSt == nSt) ? "" : Form(" in STS station %d", iSt); // Histogram title suffix - TString sN = ""; - TString sT = ""; - - // place directories in a prefered order - MakeQaDirectory(sF + "occup/"); - if (IsMCUsed()) { - MakeQaDirectory(sF + "res/"); - MakeQaDirectory(sF + "pull/"); - MakeQaDirectory(sF + "eff/"); - } - MakeQaDirectory(sF + "err/"); - - int nBinsZ = (iSt < nSt) ? kNbins : kNbinsZ; - - // Hit occupancy - sN = (TString) "hit_xy" + nsuff; - sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]"; - fvph_hit_xy[iSt] = - MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]); - - sN = (TString) "hit_zx" + nsuff; - sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]"; - fvph_hit_zx[iSt] = - MakeQaObject<TH2F>(sF + "occup/" + sN, sT, nBinsZ, frZmin[iSt], frZmax[iSt], kNbins, frXmin[iSt], frXmax[iSt]); - - sN = (TString) "hit_zy" + nsuff; - sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]"; - fvph_hit_zy[iSt] = - MakeQaObject<TH2F>(sF + "occup/" + sN, sT, nBinsZ, frZmin[iSt], frZmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]); - - // Hit errors - sN = (TString) "hit_dx" + nsuff; - sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]"; - fvph_hit_dx[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDx[0], kRHitDx[1]); - - sN = (TString) "hit_dy" + nsuff; - sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]"; - fvph_hit_dy[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDy[0], kRHitDy[1]); - - sN = (TString) "hit_du" + nsuff; - sT = (TString) "Hit position error across front strips" + tsuff + ";du_{hit} [cm]"; - fvph_hit_du[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDu[0], kRHitDu[1]); - - sN = (TString) "hit_dv" + nsuff; - sT = (TString) "Hit position error across back strips" + tsuff + ";dv_{hit} [cm]"; - fvph_hit_dv[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDv[0], kRHitDv[1]); - - sN = (TString) "hit_kuv" + nsuff; - sT = (TString) "Hit error correlation between front and back strips" + tsuff + ";kuv_{hit} [unitless]"; - fvph_hit_kuv[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins / 2 * 2 + 1, -1.1, 1.1); - - sN = (TString) "hit_dt" + nsuff; - sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]"; - fvph_hit_dt[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, kNbins, kRHitDt[0], kRHitDt[1]); + auto baseInitStatus = CbmCaInputQaBase::InitHistograms(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } - sN = (TString) "hit_station_delta_z" + nsuff; - sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]"; - fvph_hit_station_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -5., 5.); - - } // loop over station index: end // ----- Initialize histograms, which are use MC-information if (IsMCUsed()) { // Resize histogram vectors - fvph_n_points_per_hit.resize(nSt + 1, nullptr); - fvph_point_xy.resize(nSt + 1, nullptr); - fvph_point_zx.resize(nSt + 1, nullptr); - fvph_point_zy.resize(nSt + 1, nullptr); - fvph_point_hit_delta_z.resize(nSt + 1, nullptr); - fvph_res_x.resize(nSt + 1, nullptr); - fvph_res_y.resize(nSt + 1, nullptr); - fvph_res_u.resize(nSt + 1, nullptr); - fvph_res_v.resize(nSt + 1, nullptr); - fvph_res_t.resize(nSt + 1, nullptr); - fvph_pull_x.resize(nSt + 1, nullptr); - fvph_pull_y.resize(nSt + 1, nullptr); - fvph_pull_u.resize(nSt + 1, nullptr); - fvph_pull_u_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr); fvph_pull_v_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr); - fvph_pull_v.resize(nSt + 1, nullptr); - fvph_pull_t.resize(nSt + 1, nullptr); - fvph_res_x_vs_x.resize(nSt + 1, nullptr); - fvph_res_y_vs_y.resize(nSt + 1, nullptr); - fvph_res_u_vs_u.resize(nSt + 1, nullptr); - fvph_res_v_vs_v.resize(nSt + 1, nullptr); - fvph_res_t_vs_t.resize(nSt + 1, nullptr); - fvph_pull_x_vs_x.resize(nSt + 1, nullptr); - fvph_pull_y_vs_y.resize(nSt + 1, nullptr); - fvph_pull_u_vs_u.resize(nSt + 1, nullptr); - fvph_pull_v_vs_v.resize(nSt + 1, nullptr); - fvph_pull_t_vs_t.resize(nSt + 1, nullptr); - fvpe_reco_eff_vs_xy.resize(nSt + 1, nullptr); - fvph_reco_eff.resize(nSt + 1, nullptr); - - for (int iSt = 0; iSt <= nSt; ++iSt) { - TString sF = ""; // Subfolder - sF += (iSt == nSt) ? "All stations/" : Form("Station %d/", iSt); - TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix - TString tsuff = (iSt == nSt) ? "" : Form(" in STS station %d", iSt); // Histogram title suffix - TString sN = ""; - TString sT = ""; - - sN = (TString) "n_points_per_hit" + nsuff; - sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit"; - fvph_n_points_per_hit[iSt] = MakeQaObject<TH1F>(sF + sN, sT, 10, -0.5, 9.5); - - // Point occupancy - sN = (TString) "point_xy" + nsuff; - sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]"; - fvph_point_xy[iSt] = - MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]); - - sN = (TString) "point_zx" + nsuff; - sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]"; - fvph_point_zx[iSt] = - MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, frZmin[iSt], frZmax[iSt], kNbins, frXmin[iSt], frXmax[iSt]); - - sN = (TString) "point_zy" + nsuff; - sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]"; - fvph_point_zy[iSt] = - MakeQaObject<TH2F>(sF + "occup/" + sN, sT, kNbinsZ, frZmin[iSt], frZmax[iSt], kNbins, frYmin[iSt], frYmax[iSt]); - - // Difference between MC input z and hit z coordinates - sN = (TString) "point_hit_delta_z" + nsuff; - sT = (TString) "Distance between STS point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]"; - fvph_point_hit_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, kNbins, -0.05, 0.05); - - sN = (TString) "res_x" + nsuff; - sT = (TString) "Residuals for X" + tsuff + ";x_{reco} - x_{MC} [cm]"; - fvph_res_x[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResX[0], kRResX[1]); - - sN = (TString) "res_y" + nsuff; - sT = (TString) "Residuals for Y" + tsuff + ";y_{reco} - y_{MC} [cm]"; - fvph_res_y[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResY[0], kRResY[1]); - - sN = (TString) "res_u" + nsuff; - sT = (TString) "Residuals for Front strip coordinate" + tsuff + ";u_{reco} - u_{MC} [cm]"; - fvph_res_u[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResU[0], kRResU[1]); - - sN = (TString) "res_v" + nsuff; - sT = (TString) "Residuals for Back strip coordinate" + tsuff + ";v_{reco} - v_{MC} [cm]"; - fvph_res_v[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResV[0], kRResV[1]); - - sN = (TString) "res_t" + nsuff; - sT = (TString) "Residuals for Time" + tsuff + ";t_{reco} - t_{MC} [ns]"; - fvph_res_t[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, kNbins, kRResT[0], kRResT[1]); - - sN = (TString) "pull_x" + nsuff; - sT = (TString) "Pulls for X" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_y" + nsuff; - sT = (TString) "Pulls for Y" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_u" + nsuff; - sT = (TString) "Pulls for Front strip coordinate" + tsuff + ";(u_{reco} - u_{MC}) / #sigma_{u}^{reco}"; - fvph_pull_u[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_v" + nsuff; - sT = (TString) "Pulls for Back strip coordinate" + tsuff + ";(v_{reco} - v_{MC}) / #sigma_{v}^{reco}"; - fvph_pull_v[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_t" + nsuff; - sT = (TString) "Pulls for Time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "res_x_vs_x" + nsuff; - sT = (TString) "Residuals for X" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]"; - fvph_res_x_vs_x[iSt] = - MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, kRResX[0], kRResX[1]); - - sN = (TString) "res_y_vs_y" + nsuff; - sT = (TString) "Residuals for Y" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]"; - fvph_res_y_vs_y[iSt] = - MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frYmin[iSt], frYmax[iSt], kNbins, kRResY[0], kRResY[1]); - - sN = (TString) "res_u_vs_u" + nsuff; - sT = (TString) "Residuals for Front strip coordinate" + tsuff + ";u_{reco} [cm];u_{reco} - u_{MC} [cm]"; - fvph_res_u_vs_u[iSt] = - MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, kRResU[0], kRResU[1]); - - sN = (TString) "res_v_vs_v" + nsuff; - sT = (TString) "Residuals for Back strip coordinate" + tsuff + ";v_{reco} [cm];v_{reco} - v_{MC} [cm]"; - fvph_res_v_vs_v[iSt] = - MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, kRResV[0], kRResV[1]); - - sN = (TString) "res_t_vs_t" + nsuff; - sT = (TString) "Residuals for Time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]"; - fvph_res_t_vs_t[iSt] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]); - - sN = (TString) "pull_x_vs_x" + nsuff; - sT = (TString) "Pulls for X" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x_vs_x[iSt] = - MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_y_vs_y" + nsuff; - sT = (TString) "Pulls for Y" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y_vs_y[iSt] = - MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frYmin[iSt], frYmax[iSt], kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_u_vs_u" + nsuff; - sT = - (TString) "Pulls for Front strip coordinate" + tsuff + ";u_{reco} [cm];(u_{reco} - u_{MC}) / #sigma_{u}^{reco}"; - fvph_pull_u_vs_u[iSt] = - MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_v_vs_v" + nsuff; - sT = - (TString) "Pulls for Back strip coordinate" + tsuff + ";v_{reco} [cm];(v_{reco} - v_{MC}) / #sigma_{v}^{reco}"; - fvph_pull_v_vs_v[iSt] = - MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "pull_t_vs_t" + nsuff; - sT = (TString) "Pulls for Time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t_vs_t[iSt] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, kNbins, 0, 0, kNbins, kRPull[0], kRPull[1]); - - sN = (TString) "reco_eff_vs_xy" + nsuff; - sT = (TString) "Hit rec. efficiency in XY" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon"; - fvpe_reco_eff_vs_xy[iSt] = MakeQaObject<TProfile2D>(sF + "eff/" + sN, sT, kNbins, frXmin[iSt], frXmax[iSt], - kNbins, frYmin[iSt], frYmax[iSt]); - - sN = (TString) "reco_eff" + nsuff; - sT = (TString) "Hit rec. efficiency in XY bins" + tsuff + ";eff"; - fvph_reco_eff[iSt] = MakeQaObject<TH1F>(sF + "eff/" + sN, sT, 130, -0.005, 1.30 - 0.005); - } // stations - for (int idig = 0; idig <= fkMaxDigisInClusterForPulls; idig++) { TString sN = "All stations/pull/Ndigi/pull_u_"; TString sT = "Pulls for Front strip coordinate, hits with "; @@ -1421,7 +228,7 @@ InitStatus CbmCaInputQaSts::InitHistograms() sT += Form("%d digi", idig); } sT += "; (u_{reco} - u_{MC}) / #sigma_{u}^{reco}"; - fvph_pull_u_Ndig[idig] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPull[0], kRPull[1]); + fvph_pull_u_Ndig[idig] = MakeQaObject<TH1F>(sN, sT, kNbinsPull, kRPull[0], kRPull[1]); } for (int idig = 0; idig <= fkMaxDigisInClusterForPulls; idig++) { @@ -1436,30 +243,10 @@ InitStatus CbmCaInputQaSts::InitHistograms() sT += Form("%d digi", idig); } sT += "; (v_{reco} - v_{MC}) / #sigma_{v}^{reco}"; - fvph_pull_v_Ndig[idig] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPull[0], kRPull[1]); + fvph_pull_v_Ndig[idig] = MakeQaObject<TH1F>(sN, sT, kNbinsPull, kRPull[0], kRPull[1]); } } return kSUCCESS; } -// --------------------------------------------------------------------------------------------------------------------- -// -Bool_t CbmCaInputQaSts::IsTrackSelected(const CbmMCTrack* track, const CbmStsPoint* point) const -{ - - if (kTrackSelectionPrimary && track->GetMotherId() >= 0) { return false; } - - double px = point->GetPxOut(); - double py = point->GetPyOut(); - double pz = point->GetPzOut(); - double p = sqrt(px * px + py * py + pz * pz); - - if (pz <= 0.) { return false; } - - if (p < kTrackSelectionMomentum) { return false; } - - if (TMath::ATan2(sqrt(px * px + py * py), pz) * TMath::RadToDeg() > kTrackSelectionAngle) { return false; } - - return true; -} diff --git a/reco/L1/qa/CbmCaInputQaSts.h b/reco/L1/qa/CbmCaInputQaSts.h index 651965088989293be4cb29a2aa5b6cc4c2e97a53..6fef1d4d6c26b2b1afbd229629137554b2d7ff5e 100644 --- a/reco/L1/qa/CbmCaInputQaSts.h +++ b/reco/L1/qa/CbmCaInputQaSts.h @@ -37,7 +37,7 @@ class TProfile; class TProfile2D; /// A QA-task class, which provides assurance of MuCh hits and geometry -class CbmCaInputQaSts : public CbmQaTask, public CbmCaInputQaBase { +class CbmCaInputQaSts : public CbmCaInputQaBase<L1DetectorID::kSts> { public: /// Constructor from parameters /// \param verbose Verbose level @@ -62,11 +62,15 @@ protected: InitStatus InitHistograms(); /// Fills histograms per event or time-slice - void FillHistograms(); + void FillHistograms() { CbmCaInputQaBase::FillHistograms(); } /// De-initializes histograms void DeInit(); + /// @brief Defines parameters of the task + void DefineParameters(); + + /// Checks ranges for mean and standard deviation /// \return False, if variable exits the range bool CheckRangePull(TH1* h) @@ -74,121 +78,25 @@ protected: return CbmQaTask::CheckRange(h, fPullMeanThrsh, 1. - fPullWidthThrsh, 1. + fPullWidthThrsh); } - /// \brief checks if the track at mc point passes the track selection cuts - /// \return true when selected - Bool_t IsTrackSelected(const CbmMCTrack* track, const CbmStsPoint* point) const; + /// @brief Fills histograms per hit + void FillHistogramsPerHit(); -private: - // ----- Data branches - CbmStsTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface + /// @brief Fills histograms per MC point + void FillHistogramsPerPoint(); - CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice - TClonesArray* fpHits = nullptr; ///< Array of hits +private: + // ----- Data branches TClonesArray* fpClusters = nullptr; ///< Array of hit clusters - CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager - CbmMCEventList* fpMCEventList = nullptr; ///< MC event list - CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks - CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points - - TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches - - // ----- Track selection for efficiencies and pulls - - static constexpr bool kTrackSelectionPrimary {true}; // must the track come from the primary vertex - static constexpr double kTrackSelectionMomentum {0.1}; // track momentum GeV when it exits the tracking station - static constexpr double kTrackSelectionAngle {60}; // track incl. angle [grad] when it exits the tracking station - - // ----- Histograms - - static constexpr int kNbins = 200; ///< General number of bins - static constexpr int kNbinsZ = 800; ///< Number of bins for z axis in overall views - - static constexpr double kRHitDx[2] = {0., .0050}; ///< Range for hit x coordinate error [cm] - static constexpr double kRHitDy[2] = {0., .0200}; ///< Range for hit y coordinate error [cm] - static constexpr double kRHitDu[2] = {0., .0050}; ///< Range for hit u coordinate error [cm] - static constexpr double kRHitDv[2] = {0., .0050}; ///< Range for hit v coordinate error [cm] - static constexpr double kRHitDt[2] = {0., 10.}; ///< Range for hit time error [ns] - - static constexpr double kRResX[2] = {-.02, .02}; ///< Range for residual of x coordinate [cm] - static constexpr double kRResY[2] = {-.1, .1}; ///< Range for residual of y coordinate [cm] - static constexpr double kRResU[2] = {-.02, .02}; ///< Range for residual of u coordinate [cm] - static constexpr double kRResV[2] = {-.02, .02}; ///< Range for residual of v coordinate [cm] - static constexpr double kRResT[2] = {-25., 25.}; ///< Range for residual of time [ns] - - static constexpr double kRPull[2] = {-10., 10.}; ///< Range for pull histograms coordinate - - std::vector<double> frXmin; ///< Range for hit x coordinate [cm] - std::vector<double> frXmax; ///< Range for hit x coordinate [cm] - - std::vector<double> frYmin; ///< Range for hit y coordinate [cm] - std::vector<double> frYmax; ///< Range for hit y coordinate [cm] - - std::vector<double> frZmin; ///< Range for hit z coordinate [cm] - std::vector<double> frZmax; ///< Range for hit z coordinate [cm] - + // ----- Histograms (additional to ones from the base class) // NOTE: the last element of each vector stands for integral distribution over all stations - // Hit occupancy - std::vector<TH2F*> fvph_hit_xy; ///< hit y vs x in different stations - std::vector<TH2F*> fvph_hit_zx; ///< hit x vs z in different stations - std::vector<TH2F*> fvph_hit_zy; ///< hit y vs z in different stations - - std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm] - - // Hit errors - std::vector<TH1F*> fvph_hit_dx; - std::vector<TH1F*> fvph_hit_dy; - std::vector<TH1F*> fvph_hit_du; - std::vector<TH1F*> fvph_hit_dv; - std::vector<TH1F*> fvph_hit_kuv; - std::vector<TH1F*> fvph_hit_dt; - - // MC points occupancy - std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations - - std::vector<TH2F*> fvph_point_xy; ///< point y vs x in different stations - std::vector<TH2F*> fvph_point_zx; ///< point x vs z in different stations - std::vector<TH2F*> fvph_point_zy; ///< point y vs z in different stations - - std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations - - // Residuals - std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations - std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations - std::vector<TH1F*> fvph_res_u; ///< residual for u coordinate in different stations - std::vector<TH1F*> fvph_res_v; ///< residual for v coordinate in different stations - std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations - - std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations - std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations - std::vector<TH2F*> fvph_res_u_vs_u; ///< residual for u coordinate in different stations - std::vector<TH2F*> fvph_res_v_vs_v; ///< residual for v coordinate in different stations - std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations - - // Pulls - std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations - std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations - std::vector<TH1F*> fvph_pull_u; ///< pull for u coordinate in different stations - std::vector<TH1F*> fvph_pull_v; ///< pull for v coordinate in different stations - std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations - const int fkMaxDigisInClusterForPulls {5}; ///< max digis in cluster for separate histogramming of puls std::vector<TH1F*> fvph_pull_u_Ndig; ///< pull for u coordinate, depending on N digis in the cluster std::vector<TH1F*> fvph_pull_v_Ndig; ///< pull for v coordinate, depending on N digis in the cluster - std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations - std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations - std::vector<TH2F*> fvph_pull_u_vs_u; ///< pull for u coordinate in different stations - std::vector<TH2F*> fvph_pull_v_vs_v; ///< pull for v coordinate in different stations - std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations - - // Hit efficiencies - std::vector<TProfile2D*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC - std::vector<TH1F*> fvph_reco_eff; ///< Distribution of hit reconstruction efficiency in bins of fvpe_reco_eff_vs_xy - ClassDef(CbmCaInputQaSts, 0); }; diff --git a/reco/L1/qa/CbmTrackerInputQaTof.cxx b/reco/L1/qa/CbmTrackerInputQaTof.cxx index 62a99819c8b0a928a504058e5a26198f80567b19..5f374d598320bd84558f1902e48d0da08f17b2a7 100644 --- a/reco/L1/qa/CbmTrackerInputQaTof.cxx +++ b/reco/L1/qa/CbmTrackerInputQaTof.cxx @@ -375,8 +375,8 @@ int CbmTrackerInputQaTof::GetStationIndex(CbmTofPoint* p) int iSta = -1; auto tofInterface = CbmTofTrackingInterface::Instance(); for (int iSt = 0; iSt < fNtrackingStations; iSt++) { - if (fabs(p->GetZ() - tofInterface->GetZref(iSt)) < dist) { - dist = fabs(p->GetZ() - tofInterface->GetZref(iSt)); + if (fabs(p->GetZ() - tofInterface->GetZ(iSt)) < dist) { + dist = fabs(p->GetZ() - tofInterface->GetZ(iSt)); iSta = iSt; } } @@ -400,7 +400,7 @@ InitStatus CbmTrackerInputQaTof::GeometryQa() // tofInterface->GetTimeResolution(iStation); // tofInterface->IsTimeInfoProvided(iStation); - double staZ = tofInterface->GetZref(iStation); + double staZ = tofInterface->GetZ(iStation); // check that the stations are properly ordered in Z if (((iStation > 0) && (staZ <= lastZ)) || ((staZ != staZ))) {