diff --git a/core/base/CbmMatchRecoToMC.cxx b/core/base/CbmMatchRecoToMC.cxx index 2ad0c57ab8a534ac89f6cff2175366b5b41ea68f..a01144afa24285c45f658d1e2511371c16d4cc9f 100644 --- a/core/base/CbmMatchRecoToMC.cxx +++ b/core/base/CbmMatchRecoToMC.cxx @@ -195,16 +195,12 @@ void CbmMatchRecoToMC::Exec(Option_t* /*opt*/) MatchTracks(fTrdHitMatches, fTrdPoints, fTrdTracks, fTrdTrackMatches); } else if (fTrdHits) { // MC->hit->track - if (!fbSuppressMatching) { - MatchHitsToPoints(fTrdPoints, fTrdHits, fTrdHitMatches); - } + if (!fbSuppressMatching) { MatchHitsToPoints(fTrdPoints, fTrdHits, fTrdHitMatches); } MatchTracks(fTrdHitMatches, fTrdPoints, fTrdTracks, fTrdTrackMatches); } // TOF: (Digi->MC)+(Hit->Digi)=>(Hit->MC) - if (!fbSuppressMatching) { - MatchHitsTof(fTofHitDigiMatches, fTofHits, fTofHitMatches); - } + if (!fbSuppressMatching) { MatchHitsTof(fTofHitDigiMatches, fTofHits, fTofHitMatches); } //static Int_t eventNo = 0; LOG(info) << "CbmMatchRecoToMC::Exec eventNo=" << fEventNumber++; } @@ -318,7 +314,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 + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fMuchClusterMatches = new TClonesArray("CbmMatch", 100); ioman->Register("MuchClusterMatch", "MUCH", fMuchClusterMatches, IsOutputBranchPersistent("MuchClusterMatch")); } @@ -328,7 +324,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 + fbSuppressMatching = false; // Never the less, do matching, because there are no matches fMuchPixelHitMatches = new TClonesArray("CbmMatch", 100); ioman->Register("MuchPixelHitMatch", "MUCH", fMuchPixelHitMatches, IsOutputBranchPersistent("MuchPixelHitMatch")); } diff --git a/core/base/CbmMatchRecoToMC.h b/core/base/CbmMatchRecoToMC.h index 6021b90299720d894fe1b7456d003aac22b7e9e4..f093b98ec29f61fdcf1d00171d615db5cd3026a5 100644 --- a/core/base/CbmMatchRecoToMC.h +++ b/core/base/CbmMatchRecoToMC.h @@ -127,7 +127,7 @@ 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 @@ -135,14 +135,14 @@ public: ** 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; } + 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 fbSuppressMatching = false; // Suppression of MC->hit matching, if the matches are already there + 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 6cc855074fd05546b5d81fe966174dba95963438..a30cb8aa62c6735d0c692a921aba287c46dab67c 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.cxx +++ b/core/base/CbmTrackingDetectorInterfaceBase.cxx @@ -152,11 +152,11 @@ void CbmTrackingDetectorInterfaceBase::InitConversionMatrices() 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; + 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; diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h index 57926df434a35801df46f808ff188d2685a2965a..c96a2457f7863b86c988e1a5949df2ecef4b100d 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.h +++ b/core/base/CbmTrackingDetectorInterfaceBase.h @@ -12,11 +12,12 @@ #ifndef CbmTrackingDetectorInterfaceBase_h #define CbmTrackingDetectorInterfaceBase_h 1 -#include <string> #include <array> -#include <utility> +#include <string> #include <tuple> +#include <utility> #include <vector> + #include <cmath> class CbmPixelHit; @@ -116,10 +117,8 @@ public: /// @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 - ); + 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) @@ -129,10 +128,8 @@ public: /// @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 - ); + 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) @@ -145,8 +142,7 @@ public: 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 - ); + fvConvUVtoXY[iSt][2] * fvConvUVtoXY[iSt][2] * du2 + fvConvUVtoXY[iSt][3] * fvConvUVtoXY[iSt][3] * dv2); } /// Checks detector interface: boundary conditions of the parameters @@ -164,15 +160,15 @@ protected: /// @note Indeces specification: /// |x| = |fvConvUVtoXY[0] fvConvUVtoXY[1]| * |u| /// |y| |fvConvUVtoXY[2] fvConvUVtoXY[2]| |v| - /// + /// std::vector<std::array<double, 4>> fvConvUVtoXY; - + /// @brief Vector of arrays of conversion matrices from (x,y) to (u,v) coordinate systems in different stations /// /// @note Indeces specification: /// |u| = |fvConvXYtoUV[0] fvConvXYtoUV[1]| * |x| /// |v| |fvConvXYtoUV[2] fvConvXYtoUV[2]| |y| - /// + /// std::vector<std::array<double, 4>> fvConvXYtoUV; }; diff --git a/core/data/CbmPixelHit.cxx b/core/data/CbmPixelHit.cxx index ef612099d516d158ce0f4a3af511c8b4279bab9e..52cbc53b760659eac8290f3dc47cc9eb748223b2 100644 --- a/core/data/CbmPixelHit.cxx +++ b/core/data/CbmPixelHit.cxx @@ -40,9 +40,9 @@ CbmPixelHit::~CbmPixelHit() {} std::string CbmPixelHit::ToString() const { stringstream ss; - ss << "CbmPixelHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ")[cm] err=(" - << GetDx() << "," << GetDy() << "," << GetDz() << ")[cm] dxy=" << GetDxy() << "[cm2] refId=" << GetRefId() - << "time=" << GetTime() << "+/-" << GetTimeError() << "[ns]\n"; + ss << "CbmPixelHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ") err=(" + << GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() << " refId=" << GetRefId() << endl; + return ss.str(); } diff --git a/core/data/much/CbmMuchPixelHit.cxx b/core/data/much/CbmMuchPixelHit.cxx index cc06144e1132cc9094fb90e72341b0224c27b151..afb4ba292b4f99b4fe300833663a82448c208206 100644 --- a/core/data/much/CbmMuchPixelHit.cxx +++ b/core/data/much/CbmMuchPixelHit.cxx @@ -12,7 +12,6 @@ #include "CbmHit.h" // for kMUCHPIXELHIT #include <TVector3.h> // for TVector3 -#include <sstream> CbmMuchPixelHit::CbmMuchPixelHit() : CbmPixelHit(), fPlaneId(-1), fFlag(0) { SetType(kMUCHPIXELHIT); } @@ -51,12 +50,4 @@ 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 248a3e43d77f58ba3df791fd5ae75ecf2d5feaf2..04a31ad984c1760150fca713d342106d166baa5d 100644 --- a/core/data/much/CbmMuchPixelHit.h +++ b/core/data/much/CbmMuchPixelHit.h @@ -84,11 +84,6 @@ 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 ae6c3f893ac49171fe7039360e863963346f3a10..53c965658d5935bb7f93f9a2120b4557d944033c 100644 --- a/core/data/mvd/CbmMvdHit.cxx +++ b/core/data/mvd/CbmMvdHit.cxx @@ -97,17 +97,4 @@ 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 95cd2a93196ae99368710cc76a9691294e1c72c9..a6d8efe863f8254a8ce8296072bf588e5cbc6778 100644 --- a/core/data/mvd/CbmMvdHit.h +++ b/core/data/mvd/CbmMvdHit.h @@ -66,10 +66,6 @@ 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 2fbcb69dcbb5b850dbb514683a3d9f48f0545b6f..32f5f68fabc2343cf8b75605989129202a439b1c 100644 --- a/core/data/sts/CbmStsHit.cxx +++ b/core/data/sts/CbmStsHit.cxx @@ -50,13 +50,9 @@ 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 << CbmPixelHit::ToString() << "; CbmStsHit: frontClusterId=" << fFrontClusterId - << " backClusterId=" << fBackClusterId << " du=" << fDu << "[cm] dv=" << fDv << "[cm]"; - - + ss << "StsHit: address " << GetAddress() << " | time " << GetTime() << " +- " << GetTimeError() << " | Position (" + << std::setprecision(6) << GetX() << ", " << GetY() << ", " << GetZ() << ") cm | Error (" << GetDx() << ", " + << GetDy() << ", " << GetDz() << ") cm | Cluster (" << fFrontClusterId << ", " << fBackClusterId << ")"; return ss.str(); } diff --git a/core/data/tof/CbmTofHit.cxx b/core/data/tof/CbmTofHit.cxx index 1f4553090ab9fd71c436d23cd984750da844e900..1d8cf9480343a88a47ccb7da2ff5583d6230ef43 100644 --- a/core/data/tof/CbmTofHit.cxx +++ b/core/data/tof/CbmTofHit.cxx @@ -66,11 +66,9 @@ 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(); - // << " channel=" << GetCh(); // << endl; - ss << CbmPixelHit::ToString() << "; CbmTofHit: flag=" << fFlag << ", channel=" << fChannel; + ss << "CbmTofHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ") err=(" + << GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() //<< " refId=" << GetRefId() + << Form(" time=%8.2f", GetTime()) << " flag=" << GetFlag(); return ss.str(); } diff --git a/core/data/trd/CbmTrdHit.cxx b/core/data/trd/CbmTrdHit.cxx index 80eb152883999caae50796169a69083c2a96b010..3b0d42252b9bfea683229fc6031b727b799c0f95 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 6fadcb51fe84902687028efefee555fa48d203d5..33446c0ffe93aebefc9708a5684cc6dcc83b4db7 100644 --- a/core/detectors/much/CbmMuchTrackingInterface.cxx +++ b/core/detectors/much/CbmMuchTrackingInterface.cxx @@ -86,7 +86,7 @@ InitStatus CbmMuchTrackingInterface::Init() << "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 27c9d9327124881a010091b96dd5c03872dffff6..af92e9908715eb66fc812d79f139fe0f98f6be03 100644 --- a/core/detectors/mvd/CbmMvdTrackingInterface.cxx +++ b/core/detectors/mvd/CbmMvdTrackingInterface.cxx @@ -40,7 +40,7 @@ InitStatus CbmMvdTrackingInterface::Init() fMvdStationPar = CbmMvdDetector::Instance()->GetParameterFile(); if (!fMvdStationPar) { return kFATAL; } - + // Init conversion matrices InitConversionMatrices(); @@ -50,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/tof/CbmTofTrackingInterface.cxx b/core/detectors/tof/CbmTofTrackingInterface.cxx index df8f6fc3234f11391faf350b2a2393ef7b685846..9cccf72f1ec931ce3009e5b1a62977881a4b46ce 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.cxx +++ b/core/detectors/tof/CbmTofTrackingInterface.cxx @@ -92,7 +92,7 @@ InitStatus CbmTofTrackingInterface::Init() << "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 10970068e2934247f25d9951e74dc20694a4bd88..6167bf824178155f43a89735e5a74e0a680fa586 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.h +++ b/core/detectors/tof/CbmTofTrackingInterface.h @@ -14,10 +14,10 @@ #include "CbmPixelHit.h" #include "CbmTofAddress.h" +#include "CbmTofCell.h" #include "CbmTofDigiBdfPar.h" #include "CbmTofDigiPar.h" #include "CbmTofHit.h" -#include "CbmTofCell.h" #include "CbmTrackingDetectorInterfaceBase.h" #include "FairTask.h" @@ -89,10 +89,7 @@ public: /// 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()); - } + int GetTrackingStationIndex(const CbmPixelHit* hit) const { return GetTrackingStationIndex(hit->GetAddress()); } /// Gets a tracking station by the address of element /// \param address Unique element address diff --git a/macro/mcbm/mcbm_reco_event.C b/macro/mcbm/mcbm_reco_event.C index 6ccceba60658b1e7e3650565546d4b4de3d0ea2a..b2ad563f1f6a7d2720b3a2ac31dc54123d8d76c1 100644 --- a/macro/mcbm/mcbm_reco_event.C +++ b/macro/mcbm/mcbm_reco_event.C @@ -372,8 +372,8 @@ void mcbm_reco_event(Int_t nEvents = 10, TString dataset = "data/test", 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::kSts, 0, true); + //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 1, true); //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 0, false); //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 1, false); //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 2, false); @@ -381,12 +381,14 @@ void mcbm_reco_event(Int_t nEvents = 10, TString dataset = "data/test", //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::kSts, 0, true, 1.e-3, 1.e-3, 5.); + //pIdealHitProducer->CreateNewHits(L1DetectorID::kSts, 1, true, 1.e-3, 1.e-3, 5.); + //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 0, true, .1, .1, 10.); //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 1, true, .1, .1, 10.); + //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 2, false, .1, .1, 10.); //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 1, false, .1, .1, 10.); //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 2, false, .1, .1, 10.); - - run->AddTask(pIdealHitProducer); + //run->AddTask(pIdealHitProducer); } // Kalman filter diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C index 1d732e5dcd108c9da41755f318ee2975d354dcda..71b5c867fbf67dd4c40e9dd7f3c485d54530fa3a 100644 --- a/macro/run/run_reco.C +++ b/macro/run/run_reco.C @@ -387,37 +387,39 @@ 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.); + if (debugWithMC) { + 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); } - //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/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 2dca5049464a69820de138ba2531521faf8d5e71..af613208b7afb827bfd34a5cee4966f77324a6bc 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -222,6 +222,7 @@ install(FILES CbmL1Counters.h L1Algo/L1TrackPar.h L1Algo/L1Track.h utils/CbmCaIdealHitProducer.h + utils/CbmCaIdealHitProducerDet.h qa/CbmCaInputQaBase.h DESTINATION include ) diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h index ab6cd543c4c1caedcd69f6bd0d4da49114188625..af9ea9ef2cc51c8f1c3ac806d58394f92b61e42b 100644 --- a/reco/L1/CbmCaMCModule.h +++ b/reco/L1/CbmCaMCModule.h @@ -408,8 +408,6 @@ 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/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h index a3f53f4181f686f1e3519afbb515608c642121e4..90c278fc4a1f6971d965801373b8743b12b94c96 100644 --- a/reco/L1/CbmL1DetectorID.h +++ b/reco/L1/CbmL1DetectorID.h @@ -35,16 +35,16 @@ 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: +/// 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...>>; + 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); }; @@ -67,16 +67,16 @@ class CbmTofHit; namespace cbm::ca { /// @brief Names of detector subsystems - /// @note These names are used for data branches IO, thus any modification can lead to the + /// @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>; + using HitTypes_t = CbmCaDetIdTypeArr_t<CbmMvdHit, CbmStsHit, CbmMuchPixelHit, CbmTrdHit, CbmTofHit>; } // namespace cbm::ca diff --git a/reco/L1/L1Algo/utils/CaAlgoRandom.h b/reco/L1/L1Algo/utils/CaAlgoRandom.h index a98d18e850e5f3a4237e1a4988807ba711561ebf..d506496c5ce371d517840ca0eb2d7c0ed1379e4b 100644 --- a/reco/L1/L1Algo/utils/CaAlgoRandom.h +++ b/reco/L1/L1Algo/utils/CaAlgoRandom.h @@ -10,6 +10,8 @@ #ifndef CaAlgoRandom_h #define CaAlgoRandom_h 1 +#include <Logger.h> + #include <cstdint> #include <random> #include <type_traits> @@ -80,8 +82,9 @@ namespace ca::algo template<typename T, std::enable_if_t<std::is_floating_point<T>::value, bool>> T ca::algo::Random::BoundedGaus(const T& mean, const T& sigma, const T& nSigmas) const { - assert(sigma > 0 && std::isfinite(sigma)); - assert(nSigmas > 0 && std::isfinite(nSigmas)); + LOG_IF(fatal, !(nSigmas > 0 && std::isfinite(nSigmas))) << "ca::algo::Random::BoundedGaus nSigmas = " << nSigmas; + LOG_IF(fatal, !(sigma > 0 && std::isfinite(sigma))) << "ca::algo::Random::BoundedGaus sigma = " << sigma; + std::normal_distribution rndGaus {mean, sigma}; double res = 0; do { diff --git a/reco/L1/qa/CbmCaInputQaTrd.cxx b/reco/L1/qa/CbmCaInputQaTrd.cxx index d9a85050ed369bd6e13c1f67eaf4e238b48bb211..ae7dea1a222346af574a100d96fcf11fcda3f704 100644 --- a/reco/L1/qa/CbmCaInputQaTrd.cxx +++ b/reco/L1/qa/CbmCaInputQaTrd.cxx @@ -312,7 +312,6 @@ void CbmCaInputQaTrd::FillHistograms() double dxHit = pHit->GetDx(); double dyHit = pHit->GetDy(); - //double dxyHit = pHit->GetDxy(); // Hit time double tHit = pHit->GetTime(); diff --git a/reco/L1/utils/CbmCaIdealHitProducer.cxx b/reco/L1/utils/CbmCaIdealHitProducer.cxx new file mode 100644 index 0000000000000000000000000000000000000000..873787954de4839c65a9799ff40c4221648945b8 --- /dev/null +++ b/reco/L1/utils/CbmCaIdealHitProducer.cxx @@ -0,0 +1,77 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaIdealHitProducerDetBase.h +/// @brief A FairTask to run ideal hit producers for CA tracking purposes (implementation) +/// @author S.Zharko<s.zharko@gsi.de> +/// @since 01.06.2023 + +#include "CbmCaIdealHitProducer.h" + +#include "CbmSetup.h" + +using cbm::ca::IdealHitProducer; + +ClassImp(cbm::ca::IdealHitProducer); + +// --------------------------------------------------------------------------------------------------------------------- +// +void IdealHitProducer::AdjustExistingHits(L1DetectorID detID, int iSt, bool ifSmear) +{ + switch (detID) { + case L1DetectorID::kMvd: fHitProducerMvd.AdjustExistingHits(iSt, ifSmear); break; + case L1DetectorID::kSts: fHitProducerSts.AdjustExistingHits(iSt, ifSmear); break; + case L1DetectorID::kMuch: fHitProducerMuch.AdjustExistingHits(iSt, ifSmear); break; + case L1DetectorID::kTrd: fHitProducerTrd.AdjustExistingHits(iSt, ifSmear); break; + case L1DetectorID::kTof: fHitProducerTof.AdjustExistingHits(iSt, ifSmear); break; + case L1DetectorID::kEND: + LOG(fatal) << fName << "::AdjustExistingHits(): L1DetectorID::kEND cannot be passed as a detID"; + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void IdealHitProducer::CreateNewHits(L1DetectorID detID, int iSt, bool ifSmear, double du, double dv, double dt) +{ + switch (detID) { + case L1DetectorID::kMvd: fHitProducerMvd.CreateNewHits(iSt, ifSmear, du, dv, dt); break; + case L1DetectorID::kSts: fHitProducerSts.CreateNewHits(iSt, ifSmear, du, dv, dt); break; + case L1DetectorID::kMuch: fHitProducerMuch.CreateNewHits(iSt, ifSmear, du, dv, dt); break; + case L1DetectorID::kTrd: fHitProducerTrd.CreateNewHits(iSt, ifSmear, du, dv, dt); break; + case L1DetectorID::kTof: fHitProducerTof.CreateNewHits(iSt, ifSmear, du, dv, dt); break; + case L1DetectorID::kEND: + LOG(fatal) << fName << "::AdjustExistingHits(): L1DetectorID::kEND cannot be passed as a detID"; + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus IdealHitProducer::Init() +{ + fbUseDet[L1DetectorID::kMvd] = false; //CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd); + fbUseDet[L1DetectorID::kSts] = CbmSetup::Instance()->IsActive(ECbmModuleId::kSts); + fbUseDet[L1DetectorID::kMuch] = CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch); + fbUseDet[L1DetectorID::kTrd] = CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd); + fbUseDet[L1DetectorID::kTof] = false; //CbmSetup::Instance()->IsActive(ECbmModuleId::kTof); + + + InitStatus ret = kSUCCESS; + if (fbUseDet[L1DetectorID::kMvd]) { ret = std::max(fHitProducerMvd.Init(), ret); } + if (fbUseDet[L1DetectorID::kSts]) { ret = std::max(fHitProducerSts.Init(), ret); } + if (fbUseDet[L1DetectorID::kMuch]) { ret = std::max(fHitProducerMuch.Init(), ret); } + if (fbUseDet[L1DetectorID::kTrd]) { ret = std::max(fHitProducerTrd.Init(), ret); } + if (fbUseDet[L1DetectorID::kTof]) { ret = std::max(fHitProducerTof.Init(), ret); } + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void IdealHitProducer::Exec(Option_t* option) +{ + if (fbUseDet[L1DetectorID::kMvd]) { fHitProducerMvd.Exec(option); } + if (fbUseDet[L1DetectorID::kSts]) { fHitProducerSts.Exec(option); } + if (fbUseDet[L1DetectorID::kMuch]) { fHitProducerMuch.Exec(option); } + if (fbUseDet[L1DetectorID::kTrd]) { fHitProducerTrd.Exec(option); } + if (fbUseDet[L1DetectorID::kTof]) { fHitProducerTof.Exec(option); } +} diff --git a/reco/L1/utils/CbmCaIdealHitProducer.h b/reco/L1/utils/CbmCaIdealHitProducer.h new file mode 100644 index 0000000000000000000000000000000000000000..6530c0de5ae5c83652d57a95175c8bfc72ae25c7 --- /dev/null +++ b/reco/L1/utils/CbmCaIdealHitProducer.h @@ -0,0 +1,82 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaIdealHitProducerDetBase.h +/// @brief A FairTask to run ideal hit producers for CA tracking purposes (header) +/// @author S.Zharko<s.zharko@gsi.de> +/// @since 01.06.2023 + +#ifndef CbmCaIdealHitProducer_h +#define CbmCaIdealHitProducer_h 1 + +#include "CbmCaIdealHitProducerDet.h" +#include "CbmL1DetectorID.h" + +#include "FairTask.h" + +namespace cbm::ca +{ + /// @brief Ideal hit producer task for CA tracking + class IdealHitProducer : public FairTask { + public: + /// @brief Constructor + /// @param name Name of the task + /// @param verbose Verbosity level + IdealHitProducer(const char* name, int verbose) : FairTask(name, verbose) {} + + /// @brief Destructor + ~IdealHitProducer() = default; + + /// @brief Copy constructor + IdealHitProducer(const IdealHitProducer&) = delete; + + /// @brief Move constructor + IdealHitProducer(IdealHitProducer&&) = delete; + + /// @brief Copy assignment operator + IdealHitProducer& operator=(const IdealHitProducer&) = delete; + + /// @brief Move assignment operator + IdealHitProducer& operator=(IdealHitProducer&&) = delete; + + /// @brief Sets a mode to adjust existing hits to MC points for a selected tracking station + /// @param detID Detector ID + /// @param iSt Index of station + /// @param ifSmear Flag: true - MC true position and time are smeared within error + void AdjustExistingHits(L1DetectorID detID, int iSt, bool ifSmear); + + /// @brief Sets a mode to create hits from MC points for a selected tracking station + /// @param detID Detector ID + /// @param iSt Index of station + /// @param ifSmear Flag: true - MC true position and time are smeared within error + /// @param du Error of u-coordinate measurement [cm] + /// @param dv Error of v-coordinate measurement [cm] + /// @param dt Error of time measurement [ns] + void CreateNewHits(L1DetectorID detID, int iSt, bool ifSmear, double du, double dv, double dt); + + /// @brief Initialization of the task + InitStatus Init(); + + /// @brief Re-initialization of the task + InitStatus ReInit() { return Init(); } + + /// @brief Execution of the task + void Exec(Option_t* option); + + ClassDef(IdealHitProducer, 1); + + private: + IdealHitProducerDet<L1DetectorID::kMvd> fHitProducerMvd; ///< Instance of hit producer for MVD + IdealHitProducerDet<L1DetectorID::kSts> fHitProducerSts; ///< Instance of hit producer for STS + IdealHitProducerDet<L1DetectorID::kMuch> fHitProducerMuch; ///< Instance of hit producer for MuCh + IdealHitProducerDet<L1DetectorID::kTrd> fHitProducerTrd; ///< Instance of hit producer for TRD + IdealHitProducerDet<L1DetectorID::kTof> fHitProducerTof; ///< Instance of hit producer for TOF + + CbmCaDetIdArr_t<bool> fbUseDet = {{false}}; ///< Usage flag of different detectors + }; + +} // namespace cbm::ca + + +#endif // CbmCaIdealhitProducer_h diff --git a/reco/L1/utils/CbmCaIdealHitProducerDet.cxx b/reco/L1/utils/CbmCaIdealHitProducerDet.cxx new file mode 100644 index 0000000000000000000000000000000000000000..15007488834accd7fc451d2536358ee00d6391a1 --- /dev/null +++ b/reco/L1/utils/CbmCaIdealHitProducerDet.cxx @@ -0,0 +1,18 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaIdealHitProducer.h +/// @brief Base class for the ideal hit producer task (implementation) +/// @author S.Zharko +/// @since 01.06.2023 + +#include "CbmCaIdealHitProducerDet.h" + +//templateClassImp(cbm::ca::IdealHitProducerDet); + +template class cbm::ca::IdealHitProducerDet<L1DetectorID::kMvd>; +template class cbm::ca::IdealHitProducerDet<L1DetectorID::kSts>; +template class cbm::ca::IdealHitProducerDet<L1DetectorID::kMuch>; +template class cbm::ca::IdealHitProducerDet<L1DetectorID::kTrd>; +template class cbm::ca::IdealHitProducerDet<L1DetectorID::kTof>; diff --git a/reco/L1/utils/CbmCaIdealHitProducerDet.h b/reco/L1/utils/CbmCaIdealHitProducerDet.h new file mode 100644 index 0000000000000000000000000000000000000000..2e6631c4ad26bd2ab5111568c17c16526d697edd --- /dev/null +++ b/reco/L1/utils/CbmCaIdealHitProducerDet.h @@ -0,0 +1,592 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaIdealHitProducerDetBase.h +/// @brief Base class for the ideal hit producer task (header) +/// @author S.Zharko +/// @since 30.05.2023 + +#ifndef CbmCaIdealHitProducerDet_h +#define CbmCaIdealHitProducerDet_h 1 + +#include "CbmL1DetectorID.h" +#include "CbmMCDataArray.h" +#include "CbmMCDataManager.h" +#include "CbmMCDataObject.h" +#include "CbmMCEventList.h" +#include "CbmMatch.h" +#include "CbmMuchAddress.h" +#include "CbmMuchPixelHit.h" +#include "CbmMuchPoint.h" +#include "CbmMuchTrackingInterface.h" +#include "CbmMvdHit.h" +#include "CbmMvdPoint.h" +#include "CbmMvdTrackingInterface.h" +#include "CbmStsHit.h" +#include "CbmStsPoint.h" +#include "CbmStsTrackingInterface.h" +#include "CbmTimeSlice.h" +#include "CbmTofAddress.h" +#include "CbmTofHit.h" +#include "CbmTofPoint.h" +#include "CbmTofTrackingInterface.h" +#include "CbmTrackingDetectorInterfaceBase.h" +#include "CbmTrdHit.h" +#include "CbmTrdPoint.h" +#include "CbmTrdTrackingInterface.h" + +#include "FairRootManager.h" +#include "FairTask.h" +#include "Logger.h" + +#include "TClonesArray.h" +#include "TVector3.h" + +#include <array> +#include <tuple> +#include <unordered_map> + +#include "CaAlgoRandom.h" +#include "L1Constants.h" +#include "L1Def.h" +#include "L1Undef.h" + +namespace cbm::ca +{ + /// @brief Ideal hit producer class + /// + /// + template<L1DetectorID DetID> + class IdealHitProducerDet { + public: + using Hit_t = HitTypes_t::at<DetID>; + using Point_t = PointTypes_t::at<DetID>; + + /// @brief Constructor + IdealHitProducerDet() = default; + + /// @brief Destructor + ~IdealHitProducerDet(); + + /// @brief Copy constructor + IdealHitProducerDet(const IdealHitProducerDet&) = delete; + + /// @brief Move constructor + IdealHitProducerDet(IdealHitProducerDet&&) = delete; + + /// @brief Copy assignment operator + IdealHitProducerDet& operator=(const IdealHitProducerDet&) = delete; + + /// @brief Move assignment operator + IdealHitProducerDet& operator=(IdealHitProducerDet&&) = delete; + + /// @brief Sets a mode to adjust existing hits to MC points for a selected tracking station + /// @param iSt Index of station + /// @param ifSmear Flag: true - MC true position and time are smeared within error + void AdjustExistingHits(int iSt, bool ifSmear); + + /// @brief Sets a mode to create hits from MC points for a selected tracking station + /// @param iSt Index of station + /// @param ifSmear Flag: true - MC true position and time are smeared within error + /// @param du Error of u-coordinate measurement [cm] + /// @param dv Error of v-coordinate measurement [cm] + /// @param dt Error of time measurement [ns] + void CreateNewHits(int iSt, bool ifSmear, double du, double dv, double dt); + + /// @brief Initialization of the task + InitStatus Init(); + + /// @brief Re-initialization of the task + InitStatus ReInit() { return Init(); } + + /// @brief Execution of the task + void Exec(Option_t* option); + + /// @brief Sets random seed (1 by default) + /// @param seed Random seed + void SetRandomSeed(int seed) { fRandom.SetSeed(seed); } + + private: + static constexpr double kHitWSize = 3.5; ///< half-width for the bounded gaussian [sigma] + + /// @brief Gets pointer to matched MC point + /// @param iH Index of hit + /// @return link to matched point + std::optional<CbmLink> GetMatchedPointLink(int iH); + + /// @brief Pushes back a hit into the hits branch + /// @param pHit Pointer to source hit + void PushBackHit(const Hit_t* pHit); + + + /// @brief A structure to keep hit adjustment/creation settings for each station + struct StationSetting { + double fDu = undef::kD; ///< Error of u-coordinate measurement [cm] + double fDv = undef::kD; ///< Error of v-coordinate measurement [cm] + double fDt = undef::kD; ///< Error of time measurement [ns] + + /// @brief Mode: + /// - 0: do nothing with hits + /// - 1: adjust measured hits with matched point info (hit errors used) + /// - 2: create new hits from points + short fMode = 0; + + /// @brief Flag: + /// - true: MC point quantities are smeared with gaussian using error + /// - false: Precise MC point quantities are used + bool fbSmear = false; + + /// @brief String representation of the structure object + /// @param verbose Verbosity level: + /// - 0: prints nothing + /// - 1: prints one-line log + /// - 2: prints verbose log + std::string ToString(int verbose = 1) const; + }; + + CbmTrackingDetectorInterfaceBase* fpDetInterface = nullptr; ///< Detector interface + + // ----- Input branches: + CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice + CbmMCEventList* fpMCEventList = nullptr; ///< MC event list + CbmMCDataObject* fpMCEventHeader = nullptr; ///< MC event header + CbmMCDataArray* fpBrPoints = nullptr; ///< Branch: array of MC points + + TClonesArray* fpBrHits = nullptr; ///< Branch: array of hits + TClonesArray* fpBrHitMatches = nullptr; ///< Branch: array of hit matches + + TClonesArray* fpBrHitsTmp = nullptr; ///< Temporary array of hits + TClonesArray* fpBrHitMatchesTmp = nullptr; ///< Temporary array of hit matches + + std::unordered_map<int, StationSetting> fmStationPars; ///< List of settings vs. tracking station + + ::ca::algo::Random fRandom {1}; ///< Random generator + + /// @brief Management flag, which does run the routine if there was a request + bool fbRunTheRoutine = true; + + int fHitCounter = 0; ///< Hit counter in a new branch + }; + + + // ********************************************************** + // ** Template and inline function implementations ** + // ********************************************************** + + // ------------------------------------------------------------------------------------------------------------------ + // + template<L1DetectorID DetID> + IdealHitProducerDet<DetID>::~IdealHitProducerDet() + { + if (fpBrHitsTmp) { + fpBrHitsTmp->Delete(); + delete fpBrHitsTmp; + fpBrHitsTmp = nullptr; + } + + if (fpBrHitMatchesTmp) { + fpBrHitMatchesTmp->Delete(); + delete fpBrHitMatchesTmp; + fpBrHitMatchesTmp = nullptr; + } + } + + // ------------------------------------------------------------------------------------------------------------------ + // + template<L1DetectorID DetID> + void IdealHitProducerDet<DetID>::AdjustExistingHits(int iSt, bool ifSmear) + { + StationSetting pars {}; + pars.fMode = 1; + pars.fbSmear = ifSmear; + fmStationPars[iSt] = pars; + } + + // ------------------------------------------------------------------------------------------------------------------ + // + template<L1DetectorID DetID> + void IdealHitProducerDet<DetID>::CreateNewHits(int iSt, bool ifSmear, double du, double dv, double dt) + { + StationSetting pars {}; + pars.fMode = 2; + pars.fbSmear = ifSmear; + pars.fDu = du; + pars.fDv = dv; + pars.fDt = dt; + fmStationPars[iSt] = pars; + } + + // ------------------------------------------------------------------------------------------------------------------ + // + template<L1DetectorID DetID> + InitStatus IdealHitProducerDet<DetID>::Init() + try { + auto CheckBranch = [](auto pBranch, const char* brName) { + if (!pBranch) { throw std::logic_error(Form("Branch %s is not available", brName)); } + }; + + // ------ Input branch initialization + auto* pFairManager = FairRootManager::Instance(); + auto* pMcManager = dynamic_cast<CbmMCDataManager*>(pFairManager->GetObject("MCDataManager")); + CheckBranch(pMcManager, "MCDataManager"); + + fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairManager->GetObject("TimeSlice.")); + fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairManager->GetObject("MCEventList.")); + CheckBranch(fpTimeSlice, "TimeSlice."); + CheckBranch(fpMCEventList, "MCEventList."); + fpMCEventHeader = pMcManager->GetObject("MCEventHeader."); + + switch (DetID) { + case L1DetectorID::kMvd: + fpDetInterface = CbmMvdTrackingInterface::Instance(); + fpBrPoints = pMcManager->InitBranch("MvdPoint"); + fpBrHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MvdHit")); + fpBrHitMatches = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MvdHitMatch")); + CheckBranch(fpBrPoints, "MvdPoint"); + CheckBranch(fpBrHits, "MvdHit"); + CheckBranch(fpBrHitMatches, "MvdHitMatch"); + break; + case L1DetectorID::kSts: + fpDetInterface = CbmStsTrackingInterface::Instance(); + fpBrPoints = pMcManager->InitBranch("StsPoint"); + fpBrHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsHit")); + fpBrHitMatches = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsHitMatch")); + CheckBranch(fpBrPoints, "StsPoint"); + CheckBranch(fpBrHits, "StsHit"); + CheckBranch(fpBrHitMatches, "StsHitMatch"); + break; + case L1DetectorID::kMuch: + fpDetInterface = CbmMuchTrackingInterface::Instance(); + fpBrPoints = pMcManager->InitBranch("MuchPoint"); + fpBrHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchPixelHit")); + fpBrHitMatches = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchPixelHitMatch")); + CheckBranch(fpBrPoints, "MuchPoint"); + CheckBranch(fpBrHits, "MuchPixelHit"); + CheckBranch(fpBrHitMatches, "MuchPixelHitMatch"); + break; + case L1DetectorID::kTrd: + fpDetInterface = CbmTrdTrackingInterface::Instance(); + fpBrPoints = pMcManager->InitBranch("TrdPoint"); + fpBrHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdHit")); + fpBrHitMatches = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdHitMatch")); + CheckBranch(fpBrPoints, "TrdPoint"); + CheckBranch(fpBrHits, "TrdHit"); + CheckBranch(fpBrHitMatches, "TrdHitMatch"); + break; + case L1DetectorID::kTof: + fpDetInterface = CbmTofTrackingInterface::Instance(); + fpBrPoints = pMcManager->InitBranch("TofPoint"); + fpBrHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofHit")); + fpBrHitMatches = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofHitMatch")); + CheckBranch(fpBrPoints, "TofPoint"); + CheckBranch(fpBrHits, "TofHit"); + CheckBranch(fpBrHitMatches, "TofHitMatch"); + break; + case L1DetectorID::kEND: + LOG(fatal) << "CA ideal hit producer for " << kDetName[DetID] + << ": Wrong template argument is used: L1DetectorID::kEND"; + break; + } + + // ------ Checks of the settings of hit creation/adjustment mode + int nStations = fpDetInterface->GetNtrackingStations(); + for (int iSt = 0; iSt < nStations; ++iSt) { + if (fmStationPars.find(iSt) == fmStationPars.end()) { + StationSetting pars {}; + pars.fMode = 0; + fmStationPars[iSt] = pars; + } + } + + // Skip event/time slice processing, if there are no requests for hit modification of the detector + // (fbRunTheRoutine = true, if there is at least one station with fMode > 0) + fbRunTheRoutine = + std::any_of(fmStationPars.begin(), fmStationPars.end(), [](const auto& p) { return p.second.fMode > 0; }); + + std::stringstream msg; + msg << "\033[1;31mATTENTION! IDEAL HIT PRODUCER IS USED FOR " << kDetName[DetID] << "\033[0m"; + + LOG_IF(info, fbRunTheRoutine) << msg.str(); + + + return kSUCCESS; + } + catch (const std::logic_error& error) { + LOG(error) << "CA ideal hit producer for " << kDetName[DetID] + << ": initialization failed. Reason: " << error.what(); + return kERROR; + } + + // ------------------------------------------------------------------------------------------------------------------- + // + template<L1DetectorID DetID> + void IdealHitProducerDet<DetID>::PushBackHit(const HitTypes_t::at<DetID>* pHit) + { + // Common fields for all the subsystems + auto pos = TVector3 {}; + auto dpos = TVector3 {}; + pHit->Position(pos); + pHit->PositionError(dpos); + + // Fields specific to each detector type + if constexpr (L1DetectorID::kMvd == DetID) { + auto statNr = pHit->GetStationNr(); + auto iCentrX = pHit->GetIndexCentralX(); + auto iCentrY = pHit->GetIndexCentralY(); + auto iClust = pHit->GetClusterIndex(); + auto flag = pHit->GetFlag(); + new ((*fpBrHits)[fHitCounter]) CbmMvdHit(statNr, pos, dpos, iCentrX, iCentrY, iClust, flag); + } + else { + auto address = pHit->GetAddress(); + auto time = pHit->GetTime(); + auto dtime = pHit->GetTimeError(); + if constexpr (L1DetectorID::kSts == DetID) { + auto dxy = pHit->GetDxy(); + auto iClustF = pHit->GetFrontClusterId(); + auto iClustB = pHit->GetBackClusterId(); + auto du = pHit->GetDu(); + auto dv = pHit->GetDv(); + new ((*fpBrHits)[fHitCounter]) CbmStsHit(address, pos, dpos, dxy, iClustF, iClustB, time, dtime, du, dv); + } + else if constexpr (L1DetectorID::kMuch == DetID) { + auto dxy = pHit->GetDxy(); + auto refId = pHit->GetRefId(); + auto planeId = pHit->GetPlaneId(); + LOG(info) << "MUCH: " << pHit->GetRefId() << ' ' << pHit->GetPlaneId(); + new ((*fpBrHits)[fHitCounter]) CbmMuchPixelHit(address, pos, dpos, dxy, refId, planeId, time, dtime); + } + else if constexpr (L1DetectorID::kTrd == DetID) { + auto dxy = pHit->GetDxy(); + auto refId = pHit->GetRefId(); + auto eLoss = pHit->GetELoss(); + new ((*fpBrHits)[fHitCounter]) CbmTrdHit(address, pos, dpos, dxy, refId, eLoss, time, dtime); + } + else if constexpr (L1DetectorID::kTof == DetID) { + auto refId = -1; // pHit->GetRefId(); // -> deprecated + auto flag = pHit->GetFlag(); + auto channel = pHit->GetCh(); + new ((*fpBrHits)[fHitCounter]) CbmTofHit(address, pos, dpos, refId, time, dtime, flag, channel); + } + } + } + + // ------------------------------------------------------------------------------------------------------------------- + // + template<L1DetectorID DetID> + std::optional<CbmLink> IdealHitProducerDet<DetID>::GetMatchedPointLink(int iH) + { + const auto* pHitMatch = dynamic_cast<CbmMatch*>(fpBrHitMatchesTmp->At(iH)); + if constexpr (L1DetectorID::kTof == DetID) { + LOG(fatal) << "IdealHitFinder: TOF cannot be used until the new MC points scheme will be introduced"; + } + + assert(pHitMatch); + if (pHitMatch->GetNofLinks() > 0) { + const auto& link = pHitMatch->GetMatchedLink(); + if (link.GetIndex() > -1) { return std::make_optional(link); } + } + return std::nullopt; + } + + // ------------------------------------------------------------------------------------------------------------------- + // + template<L1DetectorID DetID> + void IdealHitProducerDet<DetID>::Exec(Option_t*) + { + assert(L1DetectorID::kTof != DetID); + + // ------ Check, if there are any requirement to run the routine for this detector + if (!fbRunTheRoutine) { return; } + + // ------ Copy hit and match branches to the temporary array and clear the main arrays + if (fpBrHits) { + fpBrHitsTmp = static_cast<TClonesArray*>(fpBrHits->Clone("tmpHit")); + fpBrHits->Delete(); + } + if (fpBrHitMatches) { + fpBrHitMatchesTmp = static_cast<TClonesArray*>(fpBrHitMatches->Clone("tmpHitMatch")); + fpBrHitMatches->Delete(); + } + + assert(fpBrHits); + assert(fpBrHitMatches); + assert(fpBrHitsTmp); + assert(fpBrHitMatchesTmp); + + // ------ Fill main hits array + fHitCounter = 0; + int iCluster = 0; ///< Current cluster index + for (int iH = 0; iH < fpBrHitsTmp->GetEntriesFast(); ++iH) { + auto* pHit = static_cast<Hit_t*>(fpBrHitsTmp->At(iH)); + + // Get setting + int iSt = fpDetInterface->GetTrackingStationIndex(pHit); + const auto& setting = fmStationPars.at(iSt); + + // Skip hits from stations, where new hits are to be created from points + if (setting.fMode == 2) { continue; } + + else { + if (5 == CbmTofAddress::GetSmType(pHit->GetAddress())) { LOG(info) << "TOF is TZERO!"; } + } + + // ------ Adjust hit to matched MC point + if (setting.fMode == 1) { + // Access matched MC point + auto oLinkToPoint = GetMatchedPointLink(iH); + if (oLinkToPoint == std::nullopt) { continue; } // Hit had no links, so it's fake + auto& link = oLinkToPoint.value(); + auto* pPoint = static_cast<Point_t*>(fpBrPoints->Get(link)); + + // Get point position and time + auto posIn = TVector3 {}; + auto posOut = TVector3 {}; + + if constexpr (L1DetectorID::kTof == DetID) { + pPoint->Position(posIn); + pPoint->Position(posOut); + } + else { + pPoint->Position(posIn); + pPoint->PositionOut(posOut); + } + auto pos = 0.5 * (posIn + posOut); + double x = pos.X(); + double y = pos.Y(); + double z = pos.Z(); + double t = pPoint->GetTime() + fpMCEventList->GetEventTime(link); + + if (setting.fbSmear) { + auto dt = pHit->GetTimeError(); + auto du = pHit->GetDx(); + auto dv = pHit->GetDy(); + if constexpr (L1DetectorID::kSts == DetID) { + du = pHit->GetDu(); + dv = pHit->GetDv(); + } + auto [u, v] = fpDetInterface->ConvertXYtoUV(iSt, x, y); + u += fRandom.BoundedGaus<double>(0., du, kHitWSize); + v += fRandom.BoundedGaus<double>(0., dv, kHitWSize); + t += fRandom.BoundedGaus<double>(0., dt, kHitWSize); + std::tie(x, y) = fpDetInterface->ConvertUVtoXY(iSt, u, v); + } + pHit->SetX(x); + pHit->SetY(y); + pHit->SetZ(z); + pHit->SetTime(t); + } // IF setting.fMode == 1 + + // Check out the index of cluster + if constexpr (L1DetectorID::kSts == DetID) { + iCluster = std::max(pHit->GetFrontClusterId(), iCluster); + iCluster = std::max(pHit->GetBackClusterId(), iCluster); + } + + PushBackHit(pHit); + auto* pHitMatchNew = new ((*fpBrHitMatches)[fHitCounter]) CbmMatch(); + pHitMatchNew->AddLinks(*static_cast<CbmMatch*>(fpBrHitMatchesTmp->At(iH))); + ++fHitCounter; + } + + // ------ Create hits from points + int nEvents = fpMCEventList->GetNofEvents(); + for (int iE = 0; iE < nEvents; ++iE) { + int fileId = fpMCEventList->GetFileIdByIndex(iE); + int eventId = fpMCEventList->GetEventIdByIndex(iE); + int nPoints = fpBrPoints->Size(fileId, eventId); + + // TODO: Write point selector for TOF + // NOTE: does not implemented for MVD + for (int iP = 0; iP < nPoints; ++iP) { + auto* pPoint = static_cast<Point_t*>(fpBrPoints->Get(fileId, eventId, iP)); + int iSt = fpDetInterface->GetTrackingStationIndex(pPoint->GetDetectorID()); + const auto& setting = fmStationPars.at(iSt); + if (setting.fMode != 2) { continue; } + double du = setting.fDu; + double dv = setting.fDv; + double dt = setting.fDt; + double dz = 0.; + auto [dx2, dxy, dy2] = fpDetInterface->ConvertCovMatrixUVtoXY(iSt, du * du, dv * dv); + + // Get point position and time + auto posIn = TVector3 {}; + auto posOut = TVector3 {}; + + if constexpr (L1DetectorID::kTof == DetID) { + pPoint->Position(posIn); + pPoint->Position(posOut); + } + else { + pPoint->Position(posIn); + pPoint->PositionOut(posOut); + } + auto pos = 0.5 * (posIn + posOut); + double x = pos.X(); + double y = pos.Y(); + double z = pos.Z(); + double t = pPoint->GetTime() + fpMCEventList->GetEventTime(eventId, fileId); + + if (setting.fbSmear) { + auto [u, v] = fpDetInterface->ConvertXYtoUV(iSt, x, y); + // TODO: Provide more realistic profiles for TRD + u += fRandom.BoundedGaus<double>(0., du, kHitWSize); + v += fRandom.BoundedGaus<double>(0., dv, kHitWSize); + t += fRandom.BoundedGaus<double>(0., dt, kHitWSize); + std::tie(x, y) = fpDetInterface->ConvertUVtoXY(iSt, u, v); + } + pos.SetX(x); + pos.SetY(y); + pos.SetZ(z); + auto dpos = TVector3 {std::sqrt(dx2), std::sqrt(dy2), dz}; + + // Push back a new artificial hit + if constexpr (L1DetectorID::kMvd == DetID) { new ((*fpBrHits)[fHitCounter]) CbmMvdHit(); } + else { + auto address = pPoint->GetDetectorID(); + if constexpr (L1DetectorID::kSts == DetID) { + int32_t iClustF = iCluster; + int32_t iClustB = iCluster + 1; + new ((*fpBrHits)[fHitCounter]) CbmStsHit(address, pos, dpos, dxy, iClustF, iClustB, t, dt, du, dv); + iCluster += 2; + } + else if constexpr (L1DetectorID::kMuch == DetID) { + // TODO: What is planeID in MuCh? + int32_t refId = -1; + int32_t planeId = -1; + new ((*fpBrHits)[fHitCounter]) CbmMuchPixelHit(address, pos, dpos, dxy, refId, planeId, t, dt); + } + else if constexpr (L1DetectorID::kTrd == DetID) { + int32_t refId = -1; + auto eLoss = pPoint->GetEnergyLoss(); + new ((*fpBrHits)[fHitCounter]) CbmTrdHit(address, pos, dpos, dxy, refId, eLoss, t, dt); + } + else if constexpr (L1DetectorID::kTof == DetID) { + int32_t refId = 0; + int32_t flag = 0; + int32_t channel = 0; + new ((*fpBrHits)[fHitCounter]) CbmTofHit(address, pos, dpos, refId, t, dt, flag, channel); + } + } + + // Update hit match + auto* pHitMatchNew = new ((*fpBrHitMatches)[fHitCounter]) CbmMatch(); + pHitMatchNew->AddLink(1., iP, eventId, fileId); + ++fHitCounter; + } // iP + } + + // ------ Clear temporary hits array + fpBrHitsTmp->Delete(); + delete fpBrHitsTmp; + fpBrHitsTmp = nullptr; + + fpBrHitMatchesTmp->Delete(); + delete fpBrHitMatchesTmp; + fpBrHitMatchesTmp = nullptr; + } + +} // namespace cbm::ca + +#endif // CbmCaIdealHitProducerDet_h