From fc28c5778a5f83e591b3881ffdb2096df56a0567 Mon Sep 17 00:00:00 2001 From: "s.zharko@gsi.de" <s.zharko@gsi.de> Date: Thu, 22 Jun 2023 13:14:06 +0200 Subject: [PATCH] CA: Ideal hit finder was introduced --- core/base/CbmMatchRecoToMC.cxx | 87 ++- core/base/CbmMatchRecoToMC.h | 15 +- .../base/CbmTrackingDetectorInterfaceBase.cxx | 30 ++ core/base/CbmTrackingDetectorInterfaceBase.h | 94 +++- core/data/CbmPixelHit.cxx | 5 +- core/data/much/CbmMuchPixelHit.cxx | 9 + core/data/much/CbmMuchPixelHit.h | 5 + core/data/mvd/CbmMvdHit.cxx | 13 + core/data/mvd/CbmMvdHit.h | 4 + core/data/sts/CbmStsHit.cxx | 10 +- core/data/tof/CbmTofHit.cxx | 7 +- core/data/trd/CbmTrdHit.cxx | 2 +- .../much/CbmMuchTrackingInterface.cxx | 5 +- .../detectors/mvd/CbmMvdTrackingInterface.cxx | 5 +- .../detectors/sts/CbmStsTrackingInterface.cxx | 2 + .../detectors/tof/CbmTofTrackingInterface.cxx | 7 +- core/detectors/tof/CbmTofTrackingInterface.h | 23 +- .../detectors/trd/CbmTrdTrackingInterface.cxx | 3 + macro/mcbm/mcbm_qa.C | 3 + macro/mcbm/mcbm_reco_event.C | 19 + macro/run/CMakeLists.txt | 2 +- macro/run/run_qa.C | 10 +- macro/run/run_reco.C | 32 ++ reco/KF/CbmTrackingDetectorInterfaceInit.h | 1 + reco/L1/CMakeLists.txt | 5 + reco/L1/CbmCaMCModule.h | 3 + reco/L1/CbmCaMCModule.h.orig | 503 ------------------ reco/L1/CbmCaTimeSliceReader.h | 2 +- reco/L1/CbmL1DetectorID.h | 45 +- reco/L1/L1Algo/L1Station.h | 8 +- reco/L1/L1LinkDef.h | 6 + 31 files changed, 388 insertions(+), 577 deletions(-) delete mode 100644 reco/L1/CbmCaMCModule.h.orig diff --git a/core/base/CbmMatchRecoToMC.cxx b/core/base/CbmMatchRecoToMC.cxx index a8d46e282d..2ad0c57ab8 100644 --- a/core/base/CbmMatchRecoToMC.cxx +++ b/core/base/CbmMatchRecoToMC.cxx @@ -117,49 +117,66 @@ InitStatus CbmMatchRecoToMC::Init() { ReadAndCreateDataBranches(); + // Notification to a user about matching suppression + if (fbSuppressMatching) { + std::stringstream msg; + msg << "\033[1;31mCbmMatchRecoToMC (\"" << fName << "\"): the cluster and hit matching routines for MVD, STS, "; + msg << "MuCh, TRD, TOF will be suppressed by a request from CbmMatchRecoToMC::SuppressHitMatching():\033[0m\n"; + LOG(warn) << msg.str(); + } + return kSUCCESS; } void CbmMatchRecoToMC::Exec(Option_t* /*opt*/) { - if (fStsClusterMatches != nullptr) fStsClusterMatches->Delete(); - if (fStsHitMatches != nullptr) fStsHitMatches->Delete(); + if (!fbSuppressMatching) { + if (fMvdHitMatches != nullptr) fMvdHitMatches->Delete(); + if (fMvdClusterMatches != nullptr) fMvdClusterMatches->Delete(); + if (fStsClusterMatches != nullptr) fStsClusterMatches->Delete(); + if (fStsHitMatches != nullptr) fStsHitMatches->Delete(); + if (fTrdClusterMatches != nullptr) fTrdClusterMatches->Delete(); + if (fTrdHitMatches != nullptr) fTrdHitMatches->Delete(); + if (fMuchClusterMatches != nullptr) fMuchClusterMatches->Delete(); + if (fMuchPixelHitMatches != nullptr) fMuchPixelHitMatches->Delete(); + if (fTofHitMatches != nullptr) fTofHitMatches->Delete(); + } + if (fTrdClusterMatches != nullptr) fTrdClusterMatches->Delete(); if (fStsTrackMatches != nullptr) fStsTrackMatches->Delete(); if (fRichTrackMatches != nullptr) fRichTrackMatches->Delete(); - if (fTrdClusterMatches != nullptr) fTrdClusterMatches->Delete(); - if (fTrdHitMatches != nullptr) fTrdHitMatches->Delete(); if (fTrdTrackMatches != nullptr) fTrdTrackMatches->Delete(); - if (fMuchClusterMatches != nullptr) fMuchClusterMatches->Delete(); - if (fMuchPixelHitMatches != nullptr) fMuchPixelHitMatches->Delete(); if (fMuchTrackMatches != nullptr) fMuchTrackMatches->Delete(); - if (fMvdHitMatches != nullptr) fMvdHitMatches->Delete(); - if (fMvdClusterMatches != nullptr) fMvdClusterMatches->Delete(); - if (fTofHitMatches != nullptr) fTofHitMatches->Delete(); // ----- MVD ----- - // MVD: digi->hit - if (fMvdHits && fMvdHitMatches && !fMvdCluster) { MatchHitsMvd(fMvdHits, fMvdHitMatches); } - else { - // MVD: digi->cluster - MatchClusters(ECbmModuleId::kMvd, fMvdCluster, fMvdClusterMatches); - // MVD: cluster->hit - MatchHits(fMvdClusterMatches, fMvdHits, fMvdHitMatches); + if (!fbSuppressMatching) { + // MVD: digi->hit + if (fMvdHits && fMvdHitMatches && !fMvdCluster) { MatchHitsMvd(fMvdHits, fMvdHitMatches); } + else { + // MVD: digi->cluster + MatchClusters(ECbmModuleId::kMvd, fMvdCluster, fMvdClusterMatches); + // MVD: cluster->hit + MatchHits(fMvdClusterMatches, fMvdHits, fMvdHitMatches); + } } // ----- STS ----- - // STS: digi->cluster - MatchClusters(ECbmModuleId::kSts, fStsClusters, fStsClusterMatches); - // STS: cluster->hit - MatchHitsSts(fStsClusterMatches, fStsHits, fStsHitMatches); + if (!fbSuppressMatching) { + // STS: digi->cluster + MatchClusters(ECbmModuleId::kSts, fStsClusters, fStsClusterMatches); + // STS: cluster->hit + MatchHitsSts(fStsClusterMatches, fStsHits, fStsHitMatches); + } // STS: hit->track MatchStsTracks(fMvdHitMatches, fStsHitMatches, fMvdPoints, fStsPoints, fStsTracks, fStsTrackMatches); // ----- MUCH ----- - // MUCH: digi->cluster - MatchClusters(ECbmModuleId::kMuch, fMuchClusters, fMuchClusterMatches); - // MUCH: cluster->hit - MatchHits(fMuchClusterMatches, fMuchPixelHits, fMuchPixelHitMatches); + if (!fbSuppressMatching) { + // MUCH: digi->cluster + MatchClusters(ECbmModuleId::kMuch, fMuchClusters, fMuchClusterMatches); + // MUCH: cluster->hit + MatchHits(fMuchClusterMatches, fMuchPixelHits, fMuchPixelHitMatches); + } // MUCH: hit->track MatchTracks(fMuchPixelHitMatches, fMuchPoints, fMuchTracks, fMuchTrackMatches); @@ -171,18 +188,23 @@ void CbmMatchRecoToMC::Exec(Option_t* /*opt*/) // TRD if (fTrdClusters && fTrdHits) { // MC->digi->cluster->hit->track - MatchClusters(ECbmModuleId::kTrd, fTrdClusters, fTrdClusterMatches); - MatchHits(fTrdClusterMatches, fTrdHits, fTrdHitMatches); + if (!fbSuppressMatching) { + MatchClusters(ECbmModuleId::kTrd, fTrdClusters, fTrdClusterMatches); + MatchHits(fTrdClusterMatches, fTrdHits, fTrdHitMatches); + } MatchTracks(fTrdHitMatches, fTrdPoints, fTrdTracks, fTrdTrackMatches); } else if (fTrdHits) { // MC->hit->track - MatchHitsToPoints(fTrdPoints, fTrdHits, fTrdHitMatches); + if (!fbSuppressMatching) { + MatchHitsToPoints(fTrdPoints, fTrdHits, fTrdHitMatches); + } MatchTracks(fTrdHitMatches, fTrdPoints, fTrdTracks, fTrdTrackMatches); } // TOF: (Digi->MC)+(Hit->Digi)=>(Hit->MC) - MatchHitsTof(fTofHitDigiMatches, fTofHits, fTofHitMatches); - + if (!fbSuppressMatching) { + MatchHitsTof(fTofHitDigiMatches, fTofHits, fTofHitMatches); + } //static Int_t eventNo = 0; LOG(info) << "CbmMatchRecoToMC::Exec eventNo=" << fEventNumber++; } @@ -218,6 +240,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches() if (nullptr != fStsClusters) { fStsClusterMatches = static_cast<TClonesArray*>(ioman->GetObject("StsClusterMatch")); if (nullptr == fStsClusterMatches) { + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fStsClusterMatches = new TClonesArray("CbmMatch", 100); ioman->Register("StsClusterMatch", "STS", fStsClusterMatches, IsOutputBranchPersistent("StsClusterMatch")); } @@ -227,6 +250,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches() if (nullptr != fStsHits) { fStsHitMatches = static_cast<TClonesArray*>(ioman->GetObject("StsHitMatch")); if (nullptr == fStsHitMatches) { + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fStsHitMatches = new TClonesArray("CbmMatch", 100); ioman->Register("StsHitMatch", "STS", fStsHitMatches, IsOutputBranchPersistent("StsHitMatch")); } @@ -260,6 +284,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches() if (nullptr != fTrdClusters) { fTrdClusterMatches = static_cast<TClonesArray*>(ioman->GetObject("TrdClusterMatch")); if (nullptr == fTrdClusterMatches) { + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fTrdClusterMatches = new TClonesArray("CbmMatch", 100); ioman->Register("TrdClusterMatch", "TRD", fTrdClusterMatches, IsOutputBranchPersistent("TrdClusterMatch")); } @@ -269,6 +294,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches() if (nullptr != fTrdHits) { fTrdHitMatches = static_cast<TClonesArray*>(ioman->GetObject("TrdHitMatch")); if (nullptr == fTrdHitMatches) { + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fTrdHitMatches = new TClonesArray("CbmMatch", 100); ioman->Register("TrdHitMatch", "TRD", fTrdHitMatches, IsOutputBranchPersistent("TrdHitMatch")); } @@ -292,6 +318,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches() if (nullptr != fMuchClusters) { fMuchClusterMatches = static_cast<TClonesArray*>(ioman->GetObject("MuchClusterMatch")); if (nullptr == fMuchClusterMatches) { + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fMuchClusterMatches = new TClonesArray("CbmMatch", 100); ioman->Register("MuchClusterMatch", "MUCH", fMuchClusterMatches, IsOutputBranchPersistent("MuchClusterMatch")); } @@ -301,6 +328,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches() if (nullptr != fMuchPixelHits) { fMuchPixelHitMatches = static_cast<TClonesArray*>(ioman->GetObject("MuchPixelHitMatch")); if (nullptr == fMuchPixelHitMatches) { + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fMuchPixelHitMatches = new TClonesArray("CbmMatch", 100); ioman->Register("MuchPixelHitMatch", "MUCH", fMuchPixelHitMatches, IsOutputBranchPersistent("MuchPixelHitMatch")); } @@ -876,4 +904,5 @@ vector<Int_t> CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(CbmDigiManager* di return result; } + ClassImp(CbmMatchRecoToMC); diff --git a/core/base/CbmMatchRecoToMC.h b/core/base/CbmMatchRecoToMC.h index a1a24fd3a7..6021b90299 100644 --- a/core/base/CbmMatchRecoToMC.h +++ b/core/base/CbmMatchRecoToMC.h @@ -127,11 +127,22 @@ public: static std::vector<Int_t> GetMcTrackMotherIdsForRichHit(CbmDigiManager* digiMan, const CbmRichHit* hit, const TClonesArray* richPoints, const TClonesArray* mcTracks); + /** + ** \brief Suppresses cluster and hit matching procedures + ** + ** The function sets a flag, which suppresses the cluster and hit matching procedures for MVD, STS, MuCh, TRD and + ** TOF, if and only if all the corresponding match branches are presented in the simulation tree. If at least one + ** of the branch is absent (and the corresponding cluster and hit branches are present), the flag will be set to + ** false, and the cluster/hit matching will be executed. + **/ + void SuppressReMatching() { fbSuppressMatching = true; } + private: static Int_t fEventNumber; - Bool_t fIsMvdActive = kTRUE; // is the Mvd module active - Bool_t fbDigiExpUsed = kTRUE; // Usage of CbmTofDigiExp instead of CbmTofDigi + Bool_t fIsMvdActive = kTRUE; // is the Mvd module active + Bool_t fbDigiExpUsed = kTRUE; // Usage of CbmTofDigiExp instead of CbmTofDigi + bool fbSuppressMatching = false; // Suppression of MC->hit matching, if the matches are already there CbmMCDataArray* fMCTracks = nullptr; //! Monte-Carlo tracks CbmDigiManager* fDigiManager = nullptr; //! Interface to digi branches diff --git a/core/base/CbmTrackingDetectorInterfaceBase.cxx b/core/base/CbmTrackingDetectorInterfaceBase.cxx index 28994558d2..6cc855074f 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.cxx +++ b/core/base/CbmTrackingDetectorInterfaceBase.cxx @@ -139,3 +139,33 @@ std::string CbmTrackingDetectorInterfaceBase::ToString() const } return table.str(); } + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmTrackingDetectorInterfaceBase::InitConversionMatrices() +{ + int nStations = GetNtrackingStations(); + fvConvUVtoXY.clear(); + fvConvUVtoXY.resize(nStations); + fvConvXYtoUV.clear(); + fvConvXYtoUV.resize(nStations); + for (int iS = 0; iS < nStations; ++iS) { + double phiF = GetStripsStereoAngleFront(iS); + double phiB = GetStripsStereoAngleBack(iS); + double cF = std::cos(phiF); + double sF = std::sin(phiF); + double cB = std::cos(phiB); + double sB = std::sin(phiB); + double det = cF * sB - sF * cB; + // (x,y) -> (u,v) + fvConvUVtoXY[iS][0] = cF; + fvConvUVtoXY[iS][1] = sF; + fvConvUVtoXY[iS][2] = cB; + fvConvUVtoXY[iS][3] = sB; + // (u,v) -> (x,y) + fvConvXYtoUV[iS][0] = sB / det; + fvConvXYtoUV[iS][1] = -sF / det; + fvConvXYtoUV[iS][2] = -cB / det; + fvConvXYtoUV[iS][3] = cF / det; + } +} diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h index e6cb0c9102..57926df434 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.h +++ b/core/base/CbmTrackingDetectorInterfaceBase.h @@ -13,6 +13,11 @@ #define CbmTrackingDetectorInterfaceBase_h 1 #include <string> +#include <array> +#include <utility> +#include <tuple> +#include <vector> +#include <cmath> class CbmPixelHit; @@ -73,6 +78,12 @@ public: /// \return Local index of the tracking station virtual int GetTrackingStationIndex(int address) const = 0; + /// 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 + 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] @@ -98,12 +109,93 @@ public: // (decreases computing time, but can cause bugs) static constexpr bool kUseDynamicCast {true}; + /// @brief Conversion function (x,y) -> (u,v) + /// @param iSt Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @param x X-coordinate + /// @param y Y-coordinate + /// @return pair [u, v] + std::pair<double, double> ConvertXYtoUV(int iSt, double x, double y) const + { + return std::make_pair( + fvConvXYtoUV[iSt][0] * x + fvConvXYtoUV[iSt][1] * y, + fvConvXYtoUV[iSt][2] * x + fvConvXYtoUV[iSt][3] * y + ); + } + + /// @brief Conversion function (x,y) -> (u,v) + /// @param iSt Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @param u U-coordinate + /// @param v V-coordinate + /// @return pair [x, y] + std::pair<double, double> ConvertUVtoXY(int iSt, double u, double v) const + { + return std::make_pair( + fvConvUVtoXY[iSt][0] * u + fvConvUVtoXY[iSt][1] * v, + fvConvUVtoXY[iSt][2] * u + fvConvUVtoXY[iSt][3] * v + ); + } + + /// @brief Conversion function (du2, dv2) -> (dx2, dxy, dy2) + /// @param iSt Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// @param du2 Variance of U-coordinate measurement + /// @param dv2 Variance of V-coordinate measurement + /// @return tuple [dx2, dxy, dy2] + std::tuple<double, double, double> ConvertCovMatrixUVtoXY(int iSt, double du2, double dv2) const + { + return std::make_tuple( + fvConvUVtoXY[iSt][0] * fvConvUVtoXY[iSt][0] * du2 + fvConvUVtoXY[iSt][1] * fvConvUVtoXY[iSt][1] * dv2, + fvConvUVtoXY[iSt][0] * fvConvUVtoXY[iSt][2] * du2 + fvConvUVtoXY[iSt][1] * fvConvUVtoXY[iSt][3] * dv2, + fvConvUVtoXY[iSt][2] * fvConvUVtoXY[iSt][2] * du2 + fvConvUVtoXY[iSt][3] * fvConvUVtoXY[iSt][3] * dv2 + ); + } /// Checks detector interface: boundary conditions of the parameters bool Check() const; - /// Prints all the parameters into table and saves the table as a string + /// Prints all the parameters into table and saves the table as a strinG std::string ToString() const; + +protected: + /// @brief Initialized conversion matrices between (u,v) and (x,y) coordinate systems + void InitConversionMatrices(); + + /// @brief Vector of arrays of conversion matrices from (u,v) to (x,y) coordinate systems in different stations + /// + /// @note Indeces specification: + /// |x| = |fvConvUVtoXY[0] fvConvUVtoXY[1]| * |u| + /// |y| |fvConvUVtoXY[2] fvConvUVtoXY[2]| |v| + /// + std::vector<std::array<double, 4>> fvConvUVtoXY; + + /// @brief Vector of arrays of conversion matrices from (x,y) to (u,v) coordinate systems in different stations + /// + /// @note Indeces specification: + /// |u| = |fvConvXYtoUV[0] fvConvXYtoUV[1]| * |x| + /// |v| |fvConvXYtoUV[2] fvConvXYtoUV[2]| |y| + /// + std::vector<std::array<double, 4>> fvConvXYtoUV; }; +// ********************************************************** +// ** Template and inline function implementations ** +// ********************************************************** + +// --------------------------------------------------------------------------------------------------------------------- +// +inline int CbmTrackingDetectorInterfaceBase::GetTrackingStationIndex(double zPos) const +{ + double bestDist = 1.e+20; // arbitrary large number + int nStations = GetNtrackingStations(); + int iStSelected = -1; + for (int iSt = 0; iSt < nStations; ++iSt) { + auto dist = std::fabs(zPos - GetZ(iSt)); + if (dist < bestDist) { + bestDist = dist; + iStSelected = iSt; + } + } + return iStSelected; +} + + #endif // CbmTrackingDetectorInterfaceBase_h diff --git a/core/data/CbmPixelHit.cxx b/core/data/CbmPixelHit.cxx index ea28bc257c..ef612099d5 100644 --- a/core/data/CbmPixelHit.cxx +++ b/core/data/CbmPixelHit.cxx @@ -40,8 +40,9 @@ CbmPixelHit::~CbmPixelHit() {} std::string CbmPixelHit::ToString() const { stringstream ss; - ss << "CbmPixelHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ") err=(" - << GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() << " refId=" << GetRefId() << endl; + ss << "CbmPixelHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ")[cm] err=(" + << GetDx() << "," << GetDy() << "," << GetDz() << ")[cm] dxy=" << GetDxy() << "[cm2] refId=" << GetRefId() + << "time=" << GetTime() << "+/-" << GetTimeError() << "[ns]\n"; return ss.str(); } diff --git a/core/data/much/CbmMuchPixelHit.cxx b/core/data/much/CbmMuchPixelHit.cxx index afb4ba292b..cc06144e11 100644 --- a/core/data/much/CbmMuchPixelHit.cxx +++ b/core/data/much/CbmMuchPixelHit.cxx @@ -12,6 +12,7 @@ #include "CbmHit.h" // for kMUCHPIXELHIT #include <TVector3.h> // for TVector3 +#include <sstream> CbmMuchPixelHit::CbmMuchPixelHit() : CbmPixelHit(), fPlaneId(-1), fFlag(0) { SetType(kMUCHPIXELHIT); } @@ -50,4 +51,12 @@ CbmMuchPixelHit::CbmMuchPixelHit(int32_t address, const TVector3& pos, const TVe SetTimeError(-1.); } +std::string CbmMuchPixelHit::ToString() const +{ + std::stringstream msg; + msg << CbmPixelHit::ToString() << "; CbmMuchPixelHit: planeId=" << fPlaneId << " flag=" << fFlag; + return msg.str(); +} + + ClassImp(CbmMuchPixelHit); diff --git a/core/data/much/CbmMuchPixelHit.h b/core/data/much/CbmMuchPixelHit.h index 04a31ad984..248a3e43d7 100644 --- a/core/data/much/CbmMuchPixelHit.h +++ b/core/data/much/CbmMuchPixelHit.h @@ -84,6 +84,11 @@ public: /** Modifiers **/ void SetFlag(int32_t flag) { fFlag = flag; } + /** + * \brief String representation of the MuchPixelHit + **/ + std::string ToString() const; + private: int32_t fPlaneId; // Plane number int32_t fFlag; // Flag diff --git a/core/data/mvd/CbmMvdHit.cxx b/core/data/mvd/CbmMvdHit.cxx index 53c965658d..ae6c3f893a 100644 --- a/core/data/mvd/CbmMvdHit.cxx +++ b/core/data/mvd/CbmMvdHit.cxx @@ -97,4 +97,17 @@ void CbmMvdHit::Print(const Option_t* /*opt*/) const // ------------------------------------------------------------------------- +// ----- String representation of the MVD hit ----------------------------- +std::string CbmMvdHit::ToString() const +{ + std::stringstream msg; + msg << CbmPixelHit::ToString() << "; CbmMvdHit: flag=" << GetFlag() << " detectorID=" << fDetectorID + << " fIndexCentralX=" << fIndexCentralX << " fIndexCentralY=" << fIndexCentralY + << " fClusterIndex=" << fClusterIndex; + return msg.str(); +} +// ------------------------------------------------------------------------- + + + ClassImp(CbmMvdHit) diff --git a/core/data/mvd/CbmMvdHit.h b/core/data/mvd/CbmMvdHit.h index a6d8efe863..95cd2a9319 100644 --- a/core/data/mvd/CbmMvdHit.h +++ b/core/data/mvd/CbmMvdHit.h @@ -66,6 +66,10 @@ public: // void GetDigiIndexVector(TClonesArray* cbmMvdClusterArray, std::vector<int32_t>* digiIndexVector); int32_t GetRefIndex() { return fFlag; } + /** + * \brief String representation of MVD hit + **/ + std::string ToString() const; protected: int32_t fFlag; // Hit flag; to be used later diff --git a/core/data/sts/CbmStsHit.cxx b/core/data/sts/CbmStsHit.cxx index 32f5f68fab..2fbcb69dcb 100644 --- a/core/data/sts/CbmStsHit.cxx +++ b/core/data/sts/CbmStsHit.cxx @@ -50,9 +50,13 @@ CbmStsHit::~CbmStsHit() {} string CbmStsHit::ToString() const { stringstream ss; - ss << "StsHit: address " << GetAddress() << " | time " << GetTime() << " +- " << GetTimeError() << " | Position (" - << std::setprecision(6) << GetX() << ", " << GetY() << ", " << GetZ() << ") cm | Error (" << GetDx() << ", " - << GetDy() << ", " << GetDz() << ") cm | Cluster (" << fFrontClusterId << ", " << fBackClusterId << ")"; + //ss << "StsHit: address " << GetAddress() << " | time " << GetTime() << " +- " << GetTimeError() << " | Position (" + // << std::setprecision(6) << GetX() << ", " << GetY() << ", " << GetZ() << ") cm | Error (" << GetDx() << ", " + // << GetDy() << ", " << GetDz() << ") cm | Cluster (" << fFrontClusterId << ", " << fBackClusterId << ")"; + ss << CbmPixelHit::ToString() << "; CbmStsHit: frontClusterId=" << fFrontClusterId + << " backClusterId=" << fBackClusterId << " du=" << fDu << "[cm] dv=" << fDv << "[cm]"; + + return ss.str(); } diff --git a/core/data/tof/CbmTofHit.cxx b/core/data/tof/CbmTofHit.cxx index fda9f18278..1f4553090a 100644 --- a/core/data/tof/CbmTofHit.cxx +++ b/core/data/tof/CbmTofHit.cxx @@ -66,10 +66,11 @@ CbmTofHit::~CbmTofHit() {} string CbmTofHit::ToString() const { stringstream ss; - ss << "CbmTofHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ") err=(" - << GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() //<< " refId=" << GetRefId() - << Form(" time=%8.2f", GetTime()) << " flag=" << GetFlag(); + //ss << "CbmTofHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ") err=(" + // << GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() //<< " refId=" << GetRefId() + // << Form(" time=%8.2f", GetTime()) << " flag=" << GetFlag(); // << " channel=" << GetCh(); // << endl; + ss << CbmPixelHit::ToString() << "; CbmTofHit: flag=" << fFlag << ", channel=" << fChannel; return ss.str(); } diff --git a/core/data/trd/CbmTrdHit.cxx b/core/data/trd/CbmTrdHit.cxx index 3b0d42252b..80eb152883 100644 --- a/core/data/trd/CbmTrdHit.cxx +++ b/core/data/trd/CbmTrdHit.cxx @@ -43,7 +43,7 @@ std::string CbmTrdHit::ToString() const { stringstream ss; ss << CbmPixelHit::ToString(); - ss << "CbmTrdHit" << (GetClassType() ? "2" : "1") << "D: time[ns]=" << GetTime() << "+-" << GetTimeError() + ss << "; CbmTrdHit" << (GetClassType() ? "2" : "1") << "D: time[ns]=" << GetTime() << "+-" << GetTimeError() << " eloss=" << GetELoss(); if (GetClassType()) ss << " Max=" << (GetMaxType() ? "T" : "R"); ss << " RC=" << (IsRowCross() ? 'y' : 'n') << " Ovf=" << (HasOverFlow() ? 'y' : 'n') << endl; diff --git a/core/detectors/much/CbmMuchTrackingInterface.cxx b/core/detectors/much/CbmMuchTrackingInterface.cxx index 0e43542eb5..6fadcb51fe 100644 --- a/core/detectors/much/CbmMuchTrackingInterface.cxx +++ b/core/detectors/much/CbmMuchTrackingInterface.cxx @@ -77,13 +77,16 @@ InitStatus CbmMuchTrackingInterface::Init() return kFATAL; } + // Init conversion matrices + InitConversionMatrices(); + // Check the validity of the parameters if (!this->Check()) { LOG(error) << "Some errors occurred in the tracking detector interface initialization for MuCh (see information above)"; return kFATAL; } - + return kSUCCESS; } diff --git a/core/detectors/mvd/CbmMvdTrackingInterface.cxx b/core/detectors/mvd/CbmMvdTrackingInterface.cxx index a4db0e0483..27c9d93271 100644 --- a/core/detectors/mvd/CbmMvdTrackingInterface.cxx +++ b/core/detectors/mvd/CbmMvdTrackingInterface.cxx @@ -40,6 +40,9 @@ InitStatus CbmMvdTrackingInterface::Init() fMvdStationPar = CbmMvdDetector::Instance()->GetParameterFile(); if (!fMvdStationPar) { return kFATAL; } + + // Init conversion matrices + InitConversionMatrices(); // Check the validity of the parameters if (!this->Check()) { @@ -47,7 +50,7 @@ InitStatus CbmMvdTrackingInterface::Init() << "Some errors occurred in the tracking detector interface initialization for MVD (see information above)"; return kFATAL; } - + return kSUCCESS; } diff --git a/core/detectors/sts/CbmStsTrackingInterface.cxx b/core/detectors/sts/CbmStsTrackingInterface.cxx index f57a2b520c..fdf43d62b0 100644 --- a/core/detectors/sts/CbmStsTrackingInterface.cxx +++ b/core/detectors/sts/CbmStsTrackingInterface.cxx @@ -65,6 +65,8 @@ InitStatus CbmStsTrackingInterface::Init() if (!stsSetup->IsSensorParsInit()) { stsSetup->SetSensorParameters(fStsParSetSensor); } if (!stsSetup->IsSensorCondInit()) { stsSetup->SetSensorConditions(fStsParSetSensorCond); } + // Init conversion matrices + InitConversionMatrices(); // Check the validity of the parameters if (!this->Check()) { diff --git a/core/detectors/tof/CbmTofTrackingInterface.cxx b/core/detectors/tof/CbmTofTrackingInterface.cxx index 61d5f2e744..df8f6fc323 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.cxx +++ b/core/detectors/tof/CbmTofTrackingInterface.cxx @@ -11,9 +11,7 @@ #include "CbmTofTrackingInterface.h" -#include "CbmTofCell.h" #include "CbmTofCreateDigiPar.h" -#include "CbmTofDigiPar.h" #include "FairDetector.h" #include "FairRunAna.h" @@ -85,13 +83,16 @@ InitStatus CbmTofTrackingInterface::Init() fTofStationZ[iSt] = fTofStationZ[iSt] / nTofStationModules[iSt]; } + // Init conversion matrices + InitConversionMatrices(); + // Check the validity of the parameters if (!this->Check()) { LOG(error) << "Some errors occurred in the tracking detector interface initialization for TOF (see information above)"; return kFATAL; } - + return kSUCCESS; } diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h index f6370fc789..10970068e2 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.h +++ b/core/detectors/tof/CbmTofTrackingInterface.h @@ -15,7 +15,9 @@ #include "CbmPixelHit.h" #include "CbmTofAddress.h" #include "CbmTofDigiBdfPar.h" +#include "CbmTofDigiPar.h" #include "CbmTofHit.h" +#include "CbmTofCell.h" #include "CbmTrackingDetectorInterfaceBase.h" #include "FairTask.h" @@ -89,15 +91,7 @@ public: /// \return Local index of the tracking station int GetTrackingStationIndex(const CbmPixelHit* hit) const { - int iSt = fDigiBdfPar->GetTrackingStation(CbmTofAddress::GetSmType(hit->GetAddress()), - CbmTofAddress::GetSmId(hit->GetAddress()), - CbmTofAddress::GetRpcId(hit->GetAddress())); - - if (fIfMissingHits) { - // Recalculate station index for inconsistent hits - if (hit->GetX() > 20. && hit->GetZ() > 270. && 1 == iSt) { iSt = 2; } - } - return iSt; + return GetTrackingStationIndex(hit->GetAddress()); } /// Gets a tracking station by the address of element @@ -105,17 +99,14 @@ public: /// \return Local index of the tracking station int GetTrackingStationIndex(int address) const { - LOG(fatal) << "CbmTofTrackingInterface::GetTrackingStationIndex(int address): Unfortunately this function cannot be" - << " currently used. Please, use its overloaded version " - << "CbmTofTrackingInterface::GetTrackingStationIndex(const CbmPixelHit* hit) instead"; int iSt = fDigiBdfPar->GetTrackingStation(CbmTofAddress::GetSmType(address), CbmTofAddress::GetSmId(address), CbmTofAddress::GetRpcId(address)); // NOTE: Implement, when the "mcbm_beam_2021_07_surveyed" parameters will be fixed // TODO: Invesitgate problem in mcbm_beam_2021_07_surveyed - //if (fIfMissingHits) { - // // Recalculate station index for inconsistent hits - // if (hit->GetX() > 20. && hit->GetZ() > 270. && 1 == iSt) { iSt = 2; } - //} + if (fIfMissingHits) { + auto* channelInfo = static_cast<CbmTofCell*>(fDigiPar->GetCell(address)); + if (channelInfo->GetX() > 20. && channelInfo->GetZ() > 270. && 1 == iSt) { iSt = 2; } + } return iSt; } diff --git a/core/detectors/trd/CbmTrdTrackingInterface.cxx b/core/detectors/trd/CbmTrdTrackingInterface.cxx index 00e6a8c0c9..6e0ad828af 100644 --- a/core/detectors/trd/CbmTrdTrackingInterface.cxx +++ b/core/detectors/trd/CbmTrdTrackingInterface.cxx @@ -68,6 +68,9 @@ InitStatus CbmTrdTrackingInterface::Init() } } + // Init conversion matrices + InitConversionMatrices(); + // Check the validity of the parameters if (!this->Check()) { LOG(error) diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C index 2d4d45fd99..d4f4641b3e 100644 --- a/macro/mcbm/mcbm_qa.C +++ b/macro/mcbm/mcbm_qa.C @@ -171,6 +171,9 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test", run->AddTask(mcManager); auto* matchRecoToMC = new CbmMatchRecoToMC(); + // NOTE: Matching is suppressed, if there are hit and cluster matches already in the tree. If there + // are no hit matches, they are produced on this stage. + matchRecoToMC->SuppressReMatching(); run->AddTask(matchRecoToMC); } // ------------------------------------------------------------------------ diff --git a/macro/mcbm/mcbm_reco_event.C b/macro/mcbm/mcbm_reco_event.C index c741e76129..6ccceba606 100644 --- a/macro/mcbm/mcbm_reco_event.C +++ b/macro/mcbm/mcbm_reco_event.C @@ -370,6 +370,25 @@ void mcbm_reco_event(Int_t nEvents = 10, TString dataset = "data/test", // Geometry interface initializer for tracker run->AddTask(new CbmTrackingDetectorInterfaceInit()); + if (debugWithMC) { + auto* pIdealHitProducer = new cbm::ca::IdealHitProducer("IdealHitProducer", 3); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 0, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 1, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 0, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 1, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 2, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 0, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 1, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 2, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 3, false); + //pIdealHitProducer->CreateNewHits(L1DetectorID::kSts, 1, false, 1.e-3, 1.e-3, 5.); + //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 1, true, .1, .1, 10.); + //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 1, false, .1, .1, 10.); + //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 2, false, .1, .1, 10.); + + run->AddTask(pIdealHitProducer); + } + // Kalman filter auto kalman = new CbmKF(); run->AddTask(kalman); diff --git a/macro/run/CMakeLists.txt b/macro/run/CMakeLists.txt index dd8fc5160f..7a382ad69d 100644 --- a/macro/run/CMakeLists.txt +++ b/macro/run/CMakeLists.txt @@ -578,7 +578,7 @@ EndIf() # If(DEFINED ENV{RAW_DATA_PATH} ) ##################### # ============================================================================ -Install(FILES .rootrc run_tra_file.C run_tra_beam.C run_digi.C run_reco.C run_unpack_online.C run_unpack_tsa.C qa_config.cbm.yaml +Install(FILES .rootrc run_tra_file.C run_tra_beam.C run_digi.C run_reco.C run_qa.C run_unpack_online.C run_unpack_tsa.C qa_config.cbm.yaml DESTINATION share/cbmroot/macro/run ) Install(PROGRAMS run_tests.sh diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C index 054bbb3aa7..b259bee06d 100644 --- a/macro/run/run_qa.C +++ b/macro/run/run_qa.C @@ -187,9 +187,13 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d run->AddTask(new CbmTrackingDetectorInterfaceInit()); // Geometry interface initializer for tracker // ----- Match reco to MC ------ - CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC(); - run->AddTask(matchTask); - + if (bUseMC) { + CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC(); + // NOTE: Matching is suppressed, if there are hit and cluster matches already in the tree. If there + // are no hit matches, they are produced on this stage. + matchTask->SuppressReMatching(); + run->AddTask(matchTask); + } // ----- MUCH QA --------------------------------- if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) { run->AddTask(new CbmMuchTransportQa()); diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C index 9f3be86520..1d732e5dcd 100644 --- a/macro/run/run_reco.C +++ b/macro/run/run_reco.C @@ -387,6 +387,38 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice = auto kalman = new CbmKF(); run->AddTask(kalman); + auto* pIdealHitProducer = new cbm::ca::IdealHitProducer("IdealHitProducer", 3); + for (int iSt = 0; iSt < 8; ++iSt) { + if (iSt == 2) continue; + pIdealHitProducer->CreateNewHits(L1DetectorID::kSts, iSt, false, 1.e-3, 1.e-3, 5.); + } + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 0, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 1, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 2, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 3, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 4, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 5, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 6, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 7, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 0, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 1, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 2, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 3, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 4, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 5, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 6, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 7, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 8, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 9, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 10, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 11, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 0, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 1, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 2, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 3, false); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 0, false); + run->AddTask(pIdealHitProducer); + // L1 tracking auto l1 = (debugWithMC) ? new CbmL1("CA", 2, 3) : new CbmL1("CA"); diff --git a/reco/KF/CbmTrackingDetectorInterfaceInit.h b/reco/KF/CbmTrackingDetectorInterfaceInit.h index 9e55e65145..ffc52c2c73 100644 --- a/reco/KF/CbmTrackingDetectorInterfaceInit.h +++ b/reco/KF/CbmTrackingDetectorInterfaceInit.h @@ -14,6 +14,7 @@ #include "FairTask.h" + class CbmMvdTrackingInterface; class CbmStsTrackingInterface; class CbmMuchTrackingInterface; diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 469c472599..2dca504946 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -18,6 +18,7 @@ set(INCLUDE_DIRECTORIES ${CMAKE_CURRENT_SOURCE_DIR}/qa ${CMAKE_CURRENT_SOURCE_DIR}/L1Algo/utils ${CMAKE_CURRENT_SOURCE_DIR}/catools + ${CMAKE_CURRENT_SOURCE_DIR}/utils ) set(SRCS @@ -92,6 +93,8 @@ set(SRCS qa/CbmCaTrackTypeQa.cxx qa/CbmCaTrackFitQa.cxx qa/CbmTofInteraction.cxx # Tests + + utils/CbmCaIdealHitProducer.cxx ) set(NO_DICT_SRCS @@ -116,6 +119,7 @@ set(HEADERS catools/CaToolsWindowFinder.h catools/CaToolsLinkKey.h catools/CaToolsHitRecord.h + #utils/CbmCaIdealHitProducer.h qa/CbmCaInputQaBase.h ) @@ -217,6 +221,7 @@ install(FILES CbmL1Counters.h L1Algo/L1SimdSerializer.h L1Algo/L1TrackPar.h L1Algo/L1Track.h + utils/CbmCaIdealHitProducer.h qa/CbmCaInputQaBase.h DESTINATION include ) diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h index b46bf53fc6..ab6cd543c4 100644 --- a/reco/L1/CbmCaMCModule.h +++ b/reco/L1/CbmCaMCModule.h @@ -408,6 +408,9 @@ std::optional<::ca::tools::MCPoint> cbm::ca::MCModule::FillMCPoint(int iExtId, i } } + + + // Update point time with event time time += fpMCEventList->GetEventTime(iEvent, iFile); diff --git a/reco/L1/CbmCaMCModule.h.orig b/reco/L1/CbmCaMCModule.h.orig deleted file mode 100644 index ebe5bbe25f..0000000000 --- a/reco/L1/CbmCaMCModule.h.orig +++ /dev/null @@ -1,503 +0,0 @@ -/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// @file CbmCaMCModule.h -/// @brief CA Tracking performance interface for CBM (header) -/// @since 23.09.2022 -/// @author S.Zharko <s.zharko@gsi.de> - -#ifndef CbmCaMCModule_h -#define CbmCaMCModule_h 1 - -#include "CbmL1DetectorID.h" -#include "CbmL1Hit.h" -#include "CbmL1Track.h" -#include "CbmLink.h" -#include "CbmMCDataArray.h" -#include "CbmMCEventList.h" -#include "CbmMCTrack.h" -#include "CbmMatch.h" -#include "CbmMuchPoint.h" -#include "CbmMvdPoint.h" -#include "CbmStsPoint.h" -#include "CbmTimeSlice.h" -#include "CbmTofPoint.h" -#include "CbmTrdPoint.h" - -#include "TClonesArray.h" -#include "TDatabasePDG.h" - -#include <numeric> -#include <string> -#include <string_view> -#include <type_traits> - -#include "CaToolsMCData.h" -#include "CaToolsMCPoint.h" -#include "L1Constants.h" -#include "L1InputData.h" -#include "L1Parameters.h" -#include "L1Undef.h" - -class CbmEvent; -class CbmMCDataObject; -class CbmL1HitDebugInfo; -class CbmL1Track; - -enum class L1DetectorID; - -namespace cbm::ca -{ -<<<<<<< HEAD - /// @brief Class CbmCaPerformance is an interface to communicate between -======= - /// Class CbmCaPerformance is an interface to communicate between ->>>>>>> CA: QA updates - /// - class MCModule { - public: - /// @brief Constructor - /// @param verbosity Verbosity level - /// @param perfMode Performance mode (defines cut on number of consecutive stations with hit or point) - MCModule(int verb = 0, int perfMode = 3) : fVerbose(verb), fPerformanceMode(perfMode) {} - - /// @brief Destructor - ~MCModule() = default; - - /// @brief Copy constructor - MCModule(const MCModule&) = delete; - - /// @brief Move constructor - MCModule(MCModule&&) = delete; - - /// @brief Copy assignment operator - MCModule& operator=(const MCModule&) = delete; - - /// @brief Move assignment operator - MCModule& operator=(MCModule&&) = delete; - - /// @brief Creates hits from MC points - /// - /// This function creates hits from MC points in selected tracking stations and re-fills hit arrays. - void CreateHitsFromPoints() {} - - /// @brief Sets hit parameters from the matched MC point - void AdjustHitsToPoints(); - - /// @brief Defines performance action in the end of the run - void Finish(); - - /// @brief Defines performance action in the beginning of each event or time slice - /// @note This function should be called before hits initialization - /// @param pEvent Pointer to a current CbmEvent - void InitEvent(CbmEvent* pEvent); - - /// @brief Defines action on the module in the beginning of the run - /// @return Success flag - bool InitRun(); - - /// @brief Match reconstructed and MC data - /// - /// Runs matching of hits with points and reconstructed tracks with MC ones. Reconstructed tracks are updated with - /// true information. - void MatchRecoAndMC(); - - /// @brief Processes event - /// - /// Fills histograms and tables, should be called after the tracking done - void ProcessEvent(CbmEvent* pEvent); - - /// @brief Sets first hit indexes container in different detectors - /// @param source Array of indexes - void RegisterFirstHitIndexes(const std::array<int, L1Constants::size::kMaxNdetectors + 1>& source) - { - fpvFstHitId = &source; - } - - /// @brief Sets used detector subsystems - /// @param detID Id of detector - /// @param flag Flag: true - detector is used - /// @note Should be called before this->Init() - void SetDetector(L1DetectorID detID, bool flag); - - /// @brief Sets legacy event mode: - /// @param flag Flag: - /// - true: runs on events base, - /// - false: runs on time-slice base - void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; } - - /// @brief Creates hits from MC points for a particular detector and station - /// @param detID DetectorID - /// @param iStLoc Local index of station, provided by geometry - /// @param du Error u-coordinate [cm] - /// @param dv Error v-coordinate [cm] - /// @param dt Error of time [ns] - /// @param ifSmear Flag: - /// - true: MC point quantities are smeared along the error by gaussian, - /// - false: Precise point quantities are used - ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, double du, double dv, double dt, bool ifSmear) {} - - /// @brief Adjusts existing reconstructed hits to MC points for a particular detector and station - /// @param detID DetectorID - /// @param iStLoc Local index of station, provided by geometry - /// @param ifSmear Flag: - /// - true: MC point quantities are smeared along the error by gaussian, - /// - false: Precise point quantities are used - ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, bool ifSmear) {} - - - /// @brief Registers MC data object - /// @param mcData Instance of MC data - void RegisterMCData(::ca::tools::MCData& mcData) { fpMCData = &mcData; } - - /// @brief Registers reconstructed track container - /// @param vRecoTrack Reference to reconstructed track container - void RegisterRecoTrackContainer(L1Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; } - - /// @brief Registers hit index container - /// @param vHitIds Reference to hit index container - void RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; } - - /// @brief Registers CA parameters object - /// @param pParameters A shared pointer to the parameters object - void RegisterParameters(std::shared_ptr<L1Parameters>& pParameters) { fpParameters = pParameters; } - - /// @brief Registers debug hit container - /// @param vQaHits Reference to debug hit container - void RegisterQaHitContainer(L1Vector<CbmL1HitDebugInfo>& vQaHits) { fpvQaHits = &vQaHits; } - - private: - /// @brief Check class initialization - /// @note The function throws std::logic_error, if initialization is incomplete - void CheckInit() const; - - /// @brief Initializes MC track - /// - /// Initializes information about arrangement of points and hits of MC tracks within stations and the status - /// of track ability to be reconstructed, calculates max number of points and hits on a station, number of - /// consecutive stations containing a hit or point and number of stations and points with hits. - void InitTrackInfo(); - - - /// @brief Matches hit with MC point - /// @tparam DetId Detector ID - /// @param iHitExt External index of hit - /// @return MC-point index in fvMCPoints array - /// - /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best - /// link in the match and returns the corresponding global ID of the MC points. - template<L1DetectorID DetId> - int MatchHitWithMc(int iHitExt) const; - - /// @brief Match MC points and reconstructed hits - /// - /// Writes indexes of MC points to each hit and indexes of hits to each MC point. - void MatchPointsAndHits(); - - /// @brief Matches reconstructed tracks with MC tracks - /// - /// In the procedure, the maps of associated MC track indexes to corresponding number of hits are filled out and the - /// max purity is calculated for each reconstructed track in the TS/event. If the value of purity for a given MC track - /// is larger then a threshold, the corresponding track index is saved to the reconstructed track object, in the same - /// time the index of the reconstructed track object is saved to this MC track object. If purity for the MC track does - /// not pass the threshold, its index will not be accounted as a matched track, and the index of reconstructed track - /// will be added as an index of "touched" track. - void MatchRecoAndMCTracks(); - - /// @brief Reads MC tracks from external trees and saves them to MCDataObject - void ReadMCTracks(); - - /// @brief Reads MC points from external trees and saves them to MCDataObject - void ReadMCPoints(); - - /// @brief Reads MC points in particular detector - template<L1DetectorID DetID> - void ReadMCPointsForDetector(CbmMCDataArray* pPoints); - - /// @brief Fills a single detector-specific MC point - /// @tparam DetID Detector subsystem ID - /// @param[in] iExtId Index of point in the external points container - /// @param[in] iEvent EventID of point in the external points container - /// @param[in] iFile FileID of point int the external points container - /// @param[out] intMCPoint Reference to the internal tracking MC point object (ca::tools::MCData) - /// @return true Point is correct and is to be saved to the MCData object - /// @return false Point is incorrect and will be ignored - template<L1DetectorID DetID> - bool FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point); - - std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object - - // ------ Flags - bool fbUseMvd = false; ///< is MVD used - bool fbUseSts = false; ///< is STS used - bool fbUseMuch = false; ///< is MuCh used - bool fbUseTrd = false; ///< is TRD used - bool fbUseTof = false; ///< is TOF used - - bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility) - int fVerbose = 0; ///< Verbosity level - int fPerformanceMode = -1; ///< Mode of performance - - // ------ Input data branches - const CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice - - // Mc-event - CbmMCEventList* fpMCEventList = nullptr; ///< MC event list - CbmMCDataObject* fpMCEventHeader = nullptr; ///< MC event header - CbmMCDataArray* fpMCTracks = nullptr; ///< MC tracks input - - // Mc-points - CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points input container - CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points input container - CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points input container - CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container - CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container - - // Matches - TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches - TClonesArray* fpStsHitMatches = nullptr; ///< Array of STS hit matches - TClonesArray* fpMuchHitMatches = nullptr; ///< Array of MuCh hit matches - TClonesArray* fpTrdHitMatches = nullptr; ///< Array of TRD hit matches - TClonesArray* fpTofHitMatches = nullptr; ///< Array of TOF hit matches - - TClonesArray* fpStsHits = nullptr; ///< Array of STS hits (currently needed for matching) - TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches - - // Matching information - std::set<std::pair<int, int>> fFileEventIDs; ///< Set of file and event indexes: first - iFile, second - iEvent - - // MC data - ::ca::tools::MCData* fpMCData = nullptr; ///< MC information (hits and tracks) instance - - // Reconstructed data - L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container - L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to hit index container - L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to QA hit container - - /// @brief Pointer to array of first hit indexes in the detector subsystem - /// - /// This array must be initialized in the run initialization function. - const std::array<int, L1Constants::size::kMaxNdetectors + 1>* fpvFstHitId = nullptr; - }; -} // namespace cbm::ca - - -// ********************************************** -// ** Template function implementation ** -// ********************************************** - -// --------------------------------------------------------------------------------------------------------------------- -// -template<L1DetectorID DetID> -bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point) -{ - // ----- Get positions, momenta, time and track ID - TVector3 posIn; // Position at entrance to station [cm] - TVector3 posOut; // Position at exist of station [cm] - TVector3 momIn; // 3-momentum at entrance to station [GeV/c] - TVector3 momOut; // 3-momentum at exit of station [GeV/c] - double time = undef::kD; // Time of MC point (after beginning of the event) [ns] - int iTmcExt = undef::kI32; // Track ID in the external container - // MVD - if constexpr (L1DetectorID::kMvd == DetID) { - auto* pExtPoint = dynamic_cast<CbmMvdPoint*>(fpMvdPoints->Get(iFile, iEvent, iExtId)); - if (!pExtPoint) { - LOG(warn) << "CbmCaMCModule: MVD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent - << ", iFile = " << iFile << " does not exist"; - return false; - } - pExtPoint->Position(posIn); - pExtPoint->PositionOut(posOut); - pExtPoint->Momentum(momIn); - pExtPoint->MomentumOut(momOut); - time = pExtPoint->GetTime(); - iTmcExt = pExtPoint->GetTrackID(); - } - // STS - else if constexpr (L1DetectorID::kSts == DetID) { - auto* pExtPoint = dynamic_cast<CbmStsPoint*>(fpStsPoints->Get(iFile, iEvent, iExtId)); - if (!pExtPoint) { - LOG(warn) << "CbmCaMCModule: STS MC point with iExtId = " << iExtId << ", iEvent = " << iEvent - << ", iFile = " << iFile << " does not exist"; - return false; - } - pExtPoint->Position(posIn); - pExtPoint->PositionOut(posOut); - pExtPoint->Momentum(momIn); - pExtPoint->MomentumOut(momOut); - time = pExtPoint->GetTime(); - iTmcExt = pExtPoint->GetTrackID(); - } - // MuCh - else if constexpr (L1DetectorID::kMuch == DetID) { - auto* pExtPoint = dynamic_cast<CbmMuchPoint*>(fpMuchPoints->Get(iFile, iEvent, iExtId)); - if (!pExtPoint) { - LOG(warn) << "CbmCaMCModule: MuCh MC point with iExtId = " << iExtId << ", iEvent = " << iEvent - << ", iFile = " << iFile << " does not exist"; - return false; - } - pExtPoint->Position(posIn); - pExtPoint->PositionOut(posOut); - pExtPoint->Momentum(momIn); - pExtPoint->Momentum(momOut); - time = pExtPoint->GetTime(); - iTmcExt = pExtPoint->GetTrackID(); - } - // TRD - else if constexpr (L1DetectorID::kTrd == DetID) { - auto* pExtPoint = dynamic_cast<CbmTrdPoint*>(fpTrdPoints->Get(iFile, iEvent, iExtId)); - if (!pExtPoint) { - LOG(warn) << "CbmCaMCModule: TRD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent - << ", iFile = " << iFile << " does not exist"; - return false; - } - pExtPoint->Position(posIn); - pExtPoint->PositionOut(posOut); - pExtPoint->Momentum(momIn); - pExtPoint->MomentumOut(momOut); - time = pExtPoint->GetTime(); - iTmcExt = pExtPoint->GetTrackID(); - } - // TOF - else if constexpr (L1DetectorID::kTof == DetID) { - auto* pExtPoint = dynamic_cast<CbmTofPoint*>(fpTofPoints->Get(iFile, iEvent, iExtId)); - if (!pExtPoint) { - LOG(warn) << "CbmCaMCModule: TOF MC point with iExtId = " << iExtId << ", iEvent = " << iEvent - << ", iFile = " << iFile << " does not exist"; - return false; - } - pExtPoint->Position(posIn); - pExtPoint->Position(posOut); - pExtPoint->Momentum(momIn); - pExtPoint->Momentum(momOut); - time = pExtPoint->GetTime(); - iTmcExt = pExtPoint->GetTrackID(); - } - if (iTmcExt < 0) { - LOG(warn) << "CbmCaMCModule: For MC point with iExtId = " << iExtId << ", iEvent = " << iEvent - << ", iFile = " << iFile << " MC track is undefined (has ID = " << iTmcExt << ')'; - return false; - } - TVector3 posMid = 0.5 * (posIn + posOut); - TVector3 momMid = 0.5 * (momIn + momOut); - - // // ----- Get station index - int stationID = undef::kI32; // Index of station in the array of active tracking stations - int nStationsGeo = fpParameters->GetNstationsGeometry(DetID); - // MVD, STS TODO: Try to extend for MuCh and TRD - if constexpr (L1DetectorID::kMvd == DetID || L1DetectorID::kSts == DetID) { - // Determine active station global index - double bestDist = 1.e+20; - for (int iStLocGeo = 0; iStLocGeo < nStationsGeo; ++iStLocGeo) { - int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, DetID); - if (iStActive < 0) { continue; } - const L1Station& station = fpParameters->GetStation(iStActive); - double dist = posIn.Z() - station.fZ[0]; - if (fabs(dist) < fabs(bestDist)) { - bestDist = dist; - stationID = iStActive; - } - } - } - // MuCh, TRD - else if constexpr (L1DetectorID::kMuch == DetID || L1DetectorID::kTrd == DetID) { - // Offset for MuCh - 2.5 cm, for TRD - 20. cm. TODO: Where did these values came from? - double offset = (L1DetectorID::kMuch == DetID) ? 2.5 : 20.; - for (int iStLocGeo = 0; iStLocGeo < nStationsGeo; ++iStLocGeo) { - int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, DetID); - if (iStActive < 0) { continue; } - const L1Station& station = fpParameters->GetStation(iStActive); - if (posMid.Z() > station.fZ[0] - offset) { stationID = iStActive; } - } - } - - // Update point time with event time - time += fpMCEventList->GetEventTime(iEvent, iFile); - - // ----- Reject MC points falling out of the time slice - // STS, MuCh, TRD, TOF - if constexpr (DetID != L1DetectorID::kMvd) { - double startT = fpTimeSlice->GetStartTime(); - double endT = fpTimeSlice->GetEndTime(); - - if ((startT > 0. && time < startT) || (endT > 0. && time > endT)) { - LOG(warn) << "CbmCaMCModule: MC point with iExtId = " << iExtId << ", iEvent = " << iEvent - << ", iFile = " << iFile << " and det id " << int(DetID) << " fell out of the TS duration [" << startT - << ", " << endT << "] with measured time = " << time << " [ns]"; - return false; - } - } - - // ----- Fill MC point - point.SetExternalId(iExtId); - point.SetEventId(iEvent); - point.SetFileId(iFile); - point.SetTime(time); - point.SetX(posMid.X()); - point.SetY(posMid.Y()); - point.SetZ(posMid.Z()); - point.SetXIn(posIn.X()); - point.SetYIn(posIn.Y()); - point.SetZIn(posIn.Z()); - point.SetXOut(posOut.X()); - point.SetYOut(posOut.Y()); - point.SetZOut(posOut.Z()); - point.SetPx(momMid.X()); - point.SetPy(momMid.Y()); - point.SetPz(momMid.Z()); - point.SetPxIn(momIn.X()); - point.SetPyIn(momIn.Y()); - point.SetPzIn(momIn.Z()); - point.SetPxOut(momOut.X()); - point.SetPyOut(momOut.Y()); - point.SetPzOut(momOut.Z()); - - int iTmcInt = fpMCData->FindInternalTrackIndex(iTmcExt, iEvent, iFile); - - point.SetId(fpMCData->GetNofPoints()); // select current number of points as a local id of points - point.SetTrackId(iTmcInt); - point.SetStationId(stationID); - point.SetDetectorId(DetID); - - auto* pExtTrk = L1_DYNAMIC_CAST<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTmcExt)); - if (!pExtTrk) { - LOG(warn) << "CbmCaMCModule: MC track with iTmcExt = " << iTmcExt << ", iEvent = " << iEvent - << ", iFile = " << iFile << " MC track is undefined (nullptr)"; - return false; - } - point.SetPdgCode(pExtTrk->GetPdgCode()); - point.SetMotherId(pExtTrk->GetMotherId()); - - auto* pPdgDB = TDatabasePDG::Instance()->GetParticle(point.GetPdgCode()); - point.SetMass(pPdgDB ? pPdgDB->Mass() : 0.); /// TODO: Set from track - point.SetCharge(pPdgDB ? pPdgDB->Charge() / 3. : 0.); - - return true; -} - -// --------------------------------------------------------------------------------------------------------------------- -// NOTE: template is used, because another template function FillMCPoint is used inside -template<L1DetectorID DetID> -void cbm::ca::MCModule::ReadMCPointsForDetector(CbmMCDataArray* pPoints) -{ - if constexpr (L1DetectorID::kTof != DetID) { - for (const auto& key : fFileEventIDs) { - int iFile = key.first; - int iEvent = key.second; - int nPointsEvent = pPoints->Size(iFile, iEvent); - for (int iP = 0; iP < nPointsEvent; ++iP) { - ::ca::tools::MCPoint aPoint; - if (FillMCPoint<DetID>(iP, iEvent, iFile, aPoint)) { - aPoint.SetExternalId(fpMCData->GetPointGlobExtIndex(DetID, iP)); - int iPInt = fpMCData->GetNofPoints(); // assign an index of point in internal container - if (aPoint.GetTrackId() > -1) { fpMCData->GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); } - fpMCData->AddPoint(aPoint); - } - } // iP: end - } // key: end - } -} - - -#endif // CbmCaMCModule_h diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h index 7cab89abfe..d8a67b8324 100644 --- a/reco/L1/CbmCaTimeSliceReader.h +++ b/reco/L1/CbmCaTimeSliceReader.h @@ -232,7 +232,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector() if constexpr (L1DetectorID::kMvd == DetID) { CbmMvdHit* pMvdHit = static_cast<CbmMvdHit*>(fvpBrHits[DetID]->At(iHext)); pPixelHit = static_cast<CbmPixelHit*>(pMvdHit); - iStGeom = pDetInterface->GetTrackingStationIndex(pMvdHit->GetStationNr()); + iStGeom = pDetInterface->GetTrackingStationIndex(pMvdHit); } else if constexpr (L1DetectorID::kSts == DetID) { CbmStsHit* pStsHit = static_cast<CbmStsHit*>(fvpBrHits[DetID]->At(iHext)); diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h index d4435388cc..a3f53f4181 100644 --- a/reco/L1/CbmL1DetectorID.h +++ b/reco/L1/CbmL1DetectorID.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: Sergei Zharko [committer] */ @@ -34,10 +34,49 @@ enum class L1DetectorID template<typename T> using CbmCaDetIdArr_t = L1EArray<L1DetectorID, T>; +/// @brief Array of types, indexed by L1DetectorID enum +/// +/// The array of types allows to treat different types of detector in a uniform manner +/// Example: +/// using HitTypes_t = CbmCaDetIdTypeArr_t<CbmMvdHit, CbmStsHit, CbmMuchPixelHit, CbmTrdHit, CbmTofHit>; +/// ... +/// HitTypes_t::at<L1DetectorID::kSts> hit; // Sts hit +template<class... Types> +struct CbmCaDetIdTypeArr_t { + template <L1DetectorID DetID> + using at = std::tuple_element_t<static_cast<std::size_t>(DetID), std::tuple<Types...>>; + static constexpr std::size_t size = sizeof...(Types); +}; + +/// ************************************************* +/// ** Detector-dependent common definitions ** +/// ************************************************* + +class CbmMvdPoint; +class CbmStsPoint; +class CbmMuchPoint; +class CbmTrdPoint; +class CbmTofPoint; + +class CbmMvdHit; +class CbmStsHit; +class CbmMuchPixelHit; +class CbmTrdHit; +class CbmTofHit; + namespace cbm::ca { - /// @brief Name of detector subsystems - constexpr CbmCaDetIdArr_t<const char*> kDetName = {{"MVD", "STS", "MuCh", "TRD", "TOF"}}; + /// @brief Names of detector subsystems + /// @note These names are used for data branches IO, thus any modification can lead to the + /// read-out corruption. + constexpr CbmCaDetIdArr_t<const char*> kDetName = {{"MVD", "STS", "MUCH", "TRD", "TOF"}}; + + + /// @brief Types of MC point objects for each detector + using PointTypes_t = CbmCaDetIdTypeArr_t<CbmMvdPoint, CbmStsPoint, CbmMuchPoint, CbmTrdPoint, CbmTofPoint>; + + /// @brief Types of hit objects for each detector + using HitTypes_t = CbmCaDetIdTypeArr_t<CbmMvdHit, CbmStsHit, CbmMuchPixelHit, CbmTrdHit, CbmTofHit>; } // namespace cbm::ca diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h index f22946a897..0993b435e3 100644 --- a/reco/L1/L1Algo/L1Station.h +++ b/reco/L1/L1Algo/L1Station.h @@ -77,14 +77,14 @@ public: /// Converts Cartesian coordinates X and Y into strip coordinates U and V in transverse plane /// \param x X coordinate [length unit] /// \param y Y coordinate [length unit] - /// \return Pair of U and V strip coordinates + /// \return Pair of strip coordinates [u, v] template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool> = true> std::pair<T, T> ConvXYtoUV(T x, T y) const; /// Converts strip coordinates U and V into Cartesian coordinates X and Y in transverse plane - /// \param x X coordinate [length unit] - /// \param y Y coordinate [length unit] - /// \return Pair of Y and X Cartesian coordinates + /// \param u U coordinate [length unit] + /// \param v V coordinate [length unit] + /// \return Pair of Cartesian coordinates [x, y] template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool> = true> std::pair<T, T> ConvUVtoXY(T u, T v) const; diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h index e5d1a2a947..5054f2b2b1 100644 --- a/reco/L1/L1LinkDef.h +++ b/reco/L1/L1LinkDef.h @@ -38,5 +38,11 @@ #pragma link C++ enum class L1DetectorID; #pragma link C++ class cbm::ca::OutputQa + ; #pragma link C++ class ca::tools::WindowFinder + ; +//#pragma link C++ class cbm::ca::IdealHitProducer < L1DetectorID::kMvd > + ; +//#pragma link C++ class cbm::ca::IdealHitProducer < L1DetectorID::kSts > + ; +//#pragma link C++ class cbm::ca::IdealHitProducer < L1DetectorID::kMuch > + ; +//#pragma link C++ class cbm::ca::IdealHitProducer < L1DetectorID::kTrd > + ; +//#pragma link C++ class cbm::ca::IdealHitProducer < L1DetectorID::kTof > + ; +#pragma link C++ class cbm::ca::IdealHitProducer + ; #endif -- GitLab