diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 571ed8bbb5c204d902b168c9e7f285444f761e6a..b08d0a4afe680ba8b8ef108727469c5c9871de09 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -48,7 +48,7 @@ set(SRCS L1Algo/L1Grid.cxx CbmL1Performance.cxx CbmL1ReadEvent.cxx - CbmCaPerformance.cxx + CbmCaMCModule.cxx L1Algo/L1Station.cxx L1Algo/L1Fit.cxx L1Algo/L1MCEvent.cxx @@ -69,8 +69,8 @@ set(SRCS L1Algo/utils/L1AlgoDraw.cxx L1Algo/utils/L1AlgoEfficiencyPerformance.cxx L1Algo/utils/L1AlgoPulls.cxx - catools/CaToolsMcData.cxx - catools/CaToolsMcDataManager.cxx + catools/CaToolsMCData.cxx + catools/CaToolsMCPoint.cxx catools/CaToolsPerformance.cxx catools/CaToolsDebugger.cxx catools/CaToolsWindowFinder.cxx @@ -96,6 +96,7 @@ set(NO_DICT_SRCS set(HEADERS CbmL1Constants.h CbmL1CATrdTrackFinderSA.h + CbmL1DetectorID.h CbmL1MCPoint.h CbmL1Hit.h CbmL1Track.h diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx new file mode 100644 index 0000000000000000000000000000000000000000..001fe6a53cbed5db6300fc8e7d0fe40d02483ee1 --- /dev/null +++ b/reco/L1/CbmCaMCModule.cxx @@ -0,0 +1,446 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmCaMCModule.cxx +/// \brief CA Tracking performance interface for CBM (implementation) +/// \since 23.09.2022 +/// \author S.Zharko <s.zharko@gsi.de> + +#include "CbmCaMCModule.h" + +#include "CbmEvent.h" +#include "CbmL1.h" // for L1DetectorID +#include "CbmMCDataManager.h" +#include "CbmMCDataObject.h" +#include "CbmMCEventList.h" +#include "CbmMCTrack.h" +#include "CbmTimeSlice.h" + +#include "FairEventHeader.h" +#include "FairMCEventHeader.h" +#include "FairRootManager.h" +#include "FairRunAna.h" +#include "Logger.h" + +#include "TDatabasePDG.h" +#include "TLorentzVector.h" +#include "TVector3.h" + +#include <algorithm> +#include <cassert> +#include <fstream> // TODO: SZh 07.12.2022: For debug, should be removed! +#include <stdexcept> // for std::logic_error + + +// ********************************* +// ** Action definition functions ** +// ********************************* + +// --------------------------------------------------------------------------------------------------------------------- +// +bool CbmCaMCModule::InitRun() +try { + LOG(info) << "CA MC Module: initializing CA tracking Monte-Carlo module... "; + + auto fairManager = FairRootManager::Instance(); + assert(fairManager); + + auto mcManager = dynamic_cast<CbmMCDataManager*>(fairManager->GetObject("MCDataManager")); + assert(mcManager); + + if (!fbLegacyEventMode) { + fpTimeSlice = dynamic_cast<CbmTimeSlice*>(fairManager->GetObject("TimeSlice.")); + fpMCEventList = dynamic_cast<CbmMCEventList*>(fairManager->GetObject("MCEventList.")); + } + + fpMCEventHeader = mcManager->GetObject("MCEventHeader."); + fpMCTracks = mcManager->InitBranch("MCTrack"); + + fpMvdPoints = nullptr; + fpStsPoints = nullptr; + fpMuchPoints = nullptr; + fpTrdPoints = nullptr; + fpTofPoints = nullptr; + + fpMvdHitMatches = nullptr; + fpStsHitMatches = nullptr; + fpMuchHitMatches = nullptr; + fpTrdHitMatches = nullptr; + fpTofHitMatches = nullptr; + + fpStsClusterMatches = nullptr; + + fFileEventIDs.clear(); + + if (fbUseMvd) { + LOG(info) << "CA MC Module: initializing branches for MVD"; + fpMvdPoints = mcManager->InitBranch("MvdPoint"); + fpMvdHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("MvdHitMatch")); + } + + if (fbUseSts) { + LOG(info) << "CA MC Module: initializing branches for STS"; + fpStsPoints = mcManager->InitBranch("StsPoint"); + fpStsHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsHitMatch")); + fpStsClusterMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsClusterMatch")); + } + + if (fbUseMuch) { + LOG(info) << "CA MC Module: initializing branches for MuCh"; + fpMuchPoints = mcManager->InitBranch("MuchPoint"); + fpMuchHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("MuchPixelHitMatch")); + } + + if (fbUseTrd) { + LOG(info) << "CA MC Module: initializing branches for TRD"; + fpTrdPoints = mcManager->InitBranch("TrdPoint"); + fpTrdHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("TrdHitMatch")); + } + + if (fbUseTof) { + LOG(info) << "CA MC Module: initializing branches for TOF"; + fpTofPoints = mcManager->InitBranch("TofPoint"); + fpTofHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("TofHitMatch")); + } + + // Check initialization + this->CheckInit(); + + LOG(info) << "CA MC Module: initializing CA tracking Monte-Carlo module... \033[1;32mDone!\033[0m"; + + LOG(info) << "Detector subsystems used:"; + LOG(info) << "\tMVD: " << (fbUseMvd ? "ON" : "OFF"); + LOG(info) << "\tSTS: " << (fbUseSts ? "ON" : "OFF"); + LOG(info) << "\tMuCh: " << (fbUseMuch ? "ON" : "OFF"); + LOG(info) << "\tTRD: " << (fbUseTrd ? "ON" : "OFF"); + LOG(info) << "\tTOF: " << (fbUseTof ? "ON" : "OFF"); + + return true; +} +catch (const std::logic_error& error) { + LOG(info) << "CA MC Module: initializing CA tracking Monte-Carlo module... \033[1;31mFailed\033[0m\nReason: " + << error.what(); + return false; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::InitEvent(CbmEvent* pEvent) +{ + std::cout << "\033[1;32mCALL CbmCAMCModule::InitEvent\033[0m\n"; + + // ------ Fill a set of file and event indexes + fFileEventIDs.clear(); + if (fbLegacyEventMode) { + int iFile = FairRunAna::Instance()->GetEventHeader()->GetInputFileId(); + int iEvent = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber(); + fFileEventIDs.insert({iFile, iEvent}); + } + else { + int nEvents = fpMCEventList->GetNofEvents(); + for (int iE = 0; iE < nEvents; ++iE) { + int iFile = fpMCEventList->GetFileIdByIndex(iE); + int iEvent = fpMCEventList->GetEventIdByIndex(iE); + fFileEventIDs.insert({iFile, iEvent}); + } + } + + // ------ Read data + fMCData.Clear(); + this->ReadMCTracks(); + this->ReadMCPoints(); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::Finish() { std::cout << "\033[1;32mFinishing performance\033[0m\n"; } + + +// *********************** +// ** Accessors ** +// *********************** + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::SetDetector(L1DetectorID detID, bool flag) +{ + switch (detID) { + case L1DetectorID::kMvd: fbUseMvd = flag; break; + case L1DetectorID::kSts: fbUseSts = flag; break; + case L1DetectorID::kMuch: fbUseMuch = flag; break; + case L1DetectorID::kTrd: fbUseTrd = flag; break; + case L1DetectorID::kTof: fbUseTof = flag; break; + } +} + + +// ******************************* +// ** Utility functions ** +// ******************************* + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::CheckInit() const +{ + // ----- Check parameters + if (!fpParameters) { throw std::logic_error("Tracking parameters object was not defined"); } + + // Check event list + if (!fbLegacyEventMode && !fpMCEventList) { throw std::logic_error("MC event list was not found"); } + + if (!fbLegacyEventMode && !fpTimeSlice) { throw std::logic_error("Time slice was not found"); } + + // Tracks branch + if (!fpMCTracks) { throw std::logic_error("MC tracks branch is unavailable"); } + + // Event header + if (!fpMCEventHeader) { throw std::logic_error("MC event header is unavailable"); } + + // Check detectors initialization + if (fbUseMvd) { + if (!fpMvdPoints) { throw std::logic_error("MC points unavailable for MVD"); } + if (!fpMvdHitMatches) { throw std::logic_error("Hit matches unavailable for MVD"); } + } + + if (fbUseSts) { + if (!fpStsPoints) { throw std::logic_error("MC points unavailable for STS"); } + if (!fpStsHitMatches) { throw std::logic_error("Hit matches unavailable for STS"); } + if (!fpStsClusterMatches) { throw std::logic_error("Cluster matches unavailable for STS"); } + } + + if (fbUseMuch) { + if (!fpMuchPoints) { throw std::logic_error("MC points unavailable for MuCh"); } + if (!fpMuchHitMatches) { throw std::logic_error("Hit matches unavailable for MuCh"); } + } + + if (fbUseTrd) { + if (!fpTrdPoints) { throw std::logic_error("MC points unavailable for TRD"); } + if (!fpTrdHitMatches) { throw std::logic_error("Hit matches unavailable for TRD"); } + } + + if (fbUseTof) { + if (!fpTofPoints) { throw std::logic_error("MC points unavailable for TOF"); } + if (!fpTofHitMatches) { throw std::logic_error("Hit matches unavailable for TOF"); } + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +template<> +void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* pPoints) +{ + fvNofPointsUsed[int(L1DetectorID::kTof)] = 0; + for (const auto& key : fFileEventIDs) { + int iFile = key.first; + int iEvent = key.second; + + // Prepare TODO + int nTofStationsGeo = fpParameters->GetNstationsGeometry(L1DetectorID::kTof); + std::vector<std::vector<int>> vTofPointToTrack(nTofStationsGeo); // TODO + std::vector<std::vector<float>> vTofPointToTrackDZ(nTofStationsGeo); // TODO + for (int iStLocGeo = 0; iStLocGeo < nTofStationsGeo; ++iStLocGeo) { + vTofPointToTrack[iStLocGeo].resize(fMCData.GetNofTracks(), -1); + vTofPointToTrackDZ[iStLocGeo].resize(fMCData.GetNofTracks(), 1.e+5f); + } + + // Array of flags, if a given TOF point is matched + std::vector<char> vbTofPointMatched(pPoints->Size(iFile, iEvent), 0); + + // Check whether a particular TOF point is matched to a hit + auto* pTofHits = dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("TofHit")); + for (int iHit = 0; iHit < pTofHits->GetEntriesFast(); ++iHit) { + auto* pHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fpTofHitMatches->At(iHit)); + LOG_IF(fatal, !pHitMatch) << "CA MC Module: hit match was not found for TOF hit " << iHit; + for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) { + const auto& link = pHitMatch->GetLink(iLink); + if (link.GetFile() != iFile) { continue; } + if (link.GetEntry() != iEvent) { continue; } + int iP = link.GetIndex(); + if (iP < 0) { continue; } + vbTofPointMatched[iP] = 1; + } + } + + int nPointsEvent = pPoints->Size(iFile, iEvent); + int counter = 0; + int counter1 = 0; + // loop over all TOF points in event and file + for (int iP = 0; iP < nPointsEvent; ++iP) { + if (!vbTofPointMatched[iP]) { + counter1++; + continue; + } + + auto* pExtPoint = L1_DYNAMIC_CAST<CbmTofPoint*>(pPoints->Get(iFile, iEvent, iP)); + LOG_IF(fatal, !pExtPoint) << "NO POINT: TOF: file, event, index = " << iFile << ' ' << iEvent << ' ' << iP + << '\n'; + if (!pExtPoint) { continue; } + + if (pExtPoint->GetTrackID() < 0) { continue; } + + // Cut on time-slice time + if (!fbLegacyEventMode) { + double pointT = pExtPoint->GetTime() + fpMCEventList->GetEventTime(iEvent, iFile); + double startT = fpTimeSlice->GetStartTime(); + double endT = fpTimeSlice->GetEndTime(); + if ((startT > 0. && pointT < startT) || (endT > 0. && pointT > endT)) { continue; } + } + + // Select station index for a point + double zPos = pExtPoint->GetZ(); + float bestDist = 1000.; // arbitrary large length [cm] + int iStSelected = -1; // local geometry index of TOF station + for (int iStLocGeo = 0; iStLocGeo < nTofStationsGeo; ++iStLocGeo) { + int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, L1DetectorID::kTof); + if (iStActive < 0) { continue; } + const L1Station& station = fpParameters->GetStation(iStActive); + if (fabs(zPos - station.fZ[0]) < bestDist) { + bestDist = fabs(zPos - station.fZ[0]); + iStSelected = iStLocGeo; + } + } + + int iTrack = fMCData.FindGlobalTrackIndex(pExtPoint->GetTrackID(), iEvent, iFile); + assert(iTrack > -1); + if (iStSelected != -1) { + int iStActive = fpParameters->GetStationIndexActive(iStSelected, L1DetectorID::kTof); + const L1Station& station = fpParameters->GetStation(iStActive); + if (fabs(station.fZ[0] - zPos) < vTofPointToTrackDZ[iStSelected][iTrack]) { + vTofPointToTrack[iStSelected][iTrack] = iP; + vTofPointToTrackDZ[iStSelected][iTrack] = fabs(station.fZ[0] - zPos); + } + } + } // loop over all TOF points in event and file: end + + for (int iStLocGeo = 0; iStLocGeo < nTofStationsGeo; ++iStLocGeo) { + int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, L1DetectorID::kTof); + if (iStActive < 0) { continue; } + for (int iTrk = 0; iTrk < (int) vTofPointToTrack[iStLocGeo].size(); ++iTrk) { + if (vTofPointToTrack[iStLocGeo][iTrk] < 0) { continue; } + ca::tools::MCPoint aPoint; + if (FillMCPoint<L1DetectorID::kTof>(vTofPointToTrack[iStLocGeo][iTrk], iEvent, iFile, aPoint)) { + aPoint.SetStationID(iStActive); + int iGlobP = this->CalcGlobMCPointIndex(vTofPointToTrack[iStLocGeo][iTrk], L1DetectorID::kTof); + fMCData.AddPoint(aPoint, iGlobP, iEvent, iFile); + ++fvNofPointsUsed[int(L1DetectorID::kTof)]; + } + } // loop over tracks: end + } // loop over geometry stations: end + } // loop over key: end +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::ReadMCPoints() +{ + int nPointsEstimated = 5 * fMCData.GetNofTracks() * fpParameters->GetNstationsActive(); + fMCData.ReserveNofPoints(nPointsEstimated); + + // ----- Initialise the number of points + std::fill(fvNofPointsOrig.begin(), fvNofPointsOrig.end(), 0); + + for (const auto& key : fFileEventIDs) { + int iFile = key.first; + int iEvent = key.second; + if (fbUseMvd) { fvNofPointsOrig[int(L1DetectorID::kMvd)] += fpMvdPoints->Size(iFile, iEvent); } + if (fbUseSts) { fvNofPointsOrig[int(L1DetectorID::kSts)] += fpStsPoints->Size(iFile, iEvent); } + if (fbUseMuch) { fvNofPointsOrig[int(L1DetectorID::kMuch)] += fpMuchPoints->Size(iFile, iEvent); } + if (fbUseTrd) { fvNofPointsOrig[int(L1DetectorID::kTrd)] += fpTrdPoints->Size(iFile, iEvent); } + if (fbUseTof) { fvNofPointsOrig[int(L1DetectorID::kTof)] += fpTofPoints->Size(iFile, iEvent); } + } + + // ----- Read MC points in MVD, STS, MuCh, TRD and TOF + if (fbUseMvd) { this->ReadMCPointsForDetector<L1DetectorID::kMvd>(fpMvdPoints); } + if (fbUseSts) { this->ReadMCPointsForDetector<L1DetectorID::kSts>(fpStsPoints); } + if (fbUseMuch) { this->ReadMCPointsForDetector<L1DetectorID::kMuch>(fpMuchPoints); } + if (fbUseTrd) { this->ReadMCPointsForDetector<L1DetectorID::kTrd>(fpTrdPoints); } + if (fbUseTof) { this->ReadMCPointsForDetector<L1DetectorID::kTof>(fpTofPoints); } + + + // ******* DEBUG: START + static int nCallsDBG = 0; + if (nCallsDBG < 1) { + auto openMode = nCallsDBG ? std::ios_base::app : std::ios_base::out; + std::ofstream out("mc_points_dbg_new.txt", openMode); + out << "N of points orig:\n"; + out << "\tMVD: " << fvNofPointsOrig[0] << '\n'; + out << "\tSTS: " << fvNofPointsOrig[1] << '\n'; + out << "\tMuCh: " << fvNofPointsOrig[2] << '\n'; + out << "\tTRD: " << fvNofPointsOrig[3] << '\n'; + out << "\tTOF: " << fvNofPointsOrig[4] << '\n'; + out << "N of points used:\n"; + out << "\tMVD: " << fvNofPointsUsed[0] << '\n'; + out << "\tSTS: " << fvNofPointsUsed[1] << '\n'; + out << "\tMuCh: " << fvNofPointsUsed[2] << '\n'; + out << "\tTRD: " << fvNofPointsUsed[3] << '\n'; + out << "\tTOF: " << fvNofPointsUsed[4] << '\n'; + out << "\ttotal: " << std::accumulate(fvNofPointsUsed.begin(), fvNofPointsUsed.end(), 0) << '\n'; + out << "POINTS:\n"; + + out << fMCData.GetPoint(0).ToString(2, true) << '\n'; + for (const auto& point : fMCData.GetPointContainer()) { + out << point.ToString(2) << '\n'; + } + out.close(); + ++nCallsDBG; + } + // ******* DEBUG: END +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::ReadMCTracks() +{ + // ----- Total number of tracks + int nTracksTot = 0; + for (const auto& key : fFileEventIDs) { + nTracksTot += fpMCTracks->Size(key.first, key.second); /// iFile, iEvent + } + fMCData.ReserveNofTracks(nTracksTot); + + // ----- Loop over MC events + for (const auto& key : fFileEventIDs) { + int iFile = key.first; + int iEvent = key.second; + auto pEvtHeader = dynamic_cast<FairMCEventHeader*>(fpMCEventHeader->Get(iFile, iEvent)); + if (!pEvtHeader) { + LOG(fatal) << "CbmCaMCModule: event header is not found for file " << iFile << " and event " << iEvent; + } + + // ----- Loop over MC tracks + int nTracks = fpMCTracks->Size(iFile, iEvent); + for (int iTrkExt = 0; iTrkExt < nTracks; ++iTrkExt) { + CbmMCTrack* pExtMCTrk = dynamic_cast<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTrkExt)); + if (!pExtMCTrk) { + LOG(warn) << "CbmCaMCModule: track with (iF, iE, iT) = " << iFile << ' ' << iEvent << ' ' << iTrkExt + << " not found"; + } + TVector3 vtx; + TLorentzVector mom4; + pExtMCTrk->GetStartVertex(vtx); + pExtMCTrk->Get4Momentum(mom4); + int pdg = pExtMCTrk->GetPdgCode(); + unsigned procID = pExtMCTrk->GetGeantProcessId(); + int motherID = pExtMCTrk->GetMotherId(); + double charge = 0.; // in [e] + double mass = 0.; // in [GeV/c2] + + // Initialize charge and mass from PDG + auto pPdgBase = TDatabasePDG::Instance()->GetParticle(pdg); + if (pPdgBase) { + charge = pPdgBase->Charge() / 3.; + mass = pPdgBase->Mass(); + } + + // Create a CbmL1MCTrack + int iTrkInt = fMCData.GetNofTracks(); // index of track in the MCData container + CbmL1MCTrack aTrk(mass, charge, vtx, mom4, iTrkInt, motherID, pdg, procID); + aTrk.time = pExtMCTrk->GetStartT(); + aTrk.iFile = iFile; + aTrk.iEvent = iEvent; + aTrk.isSignal = aTrk.IsPrimary() && (aTrk.z > pEvtHeader->GetZ() + 1.e-10); + fMCData.AddTrack(aTrk, iTrkExt, iEvent, iFile); + } // Loop over MC tracks: end + } // Loop over MC events: end +} diff --git a/reco/L1/CbmCaPerformance.cxx b/reco/L1/CbmCaPerformance.cxx deleted file mode 100644 index f9e4110cc5761d4245ed6e8c522d14be71234d71..0000000000000000000000000000000000000000 --- a/reco/L1/CbmCaPerformance.cxx +++ /dev/null @@ -1,143 +0,0 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/// \file CbmCaPerformance.cxx -/// \brief CA Tracking performance interface for CBM (implementation) -/// \since 23.09.2022 -/// \author S.Zharko <s.zharko@gsi.de> - -#include "CbmCaPerformance.h" - -#include "CbmEvent.h" -#include "CbmL1.h" // for L1DetectorID -#include "CbmMCDataManager.h" - -#include "FairRootManager.h" -#include "Logger.h" - -#include <cassert> -#include <stdexcept> // for std::logic_error - - -// ********************************* -// ** Action definition functions ** -// ********************************* - -// --------------------------------------------------------------------------------------------------------------------- -// -bool CbmCaPerformance::Init() -try { - LOG(info) << "Initializing CA tracking Monte-Carlo module... "; - - auto fairManager = FairRootManager::Instance(); - assert(fairManager); - - auto mcManager = dynamic_cast<CbmMCDataManager*>(fairManager->GetObject("MCDataManager")); - assert(mcManager); - - fpMvdPoints = nullptr; - fpStsPoints = nullptr; - fpMuchPoints = nullptr; - fpTrdPoints = nullptr; - fpTofPoints = nullptr; - - if (fbUseMvd) { - LOG(info) << "CA MC Module: initializing branches for MVD"; - fpMvdPoints = mcManager->InitBranch("MvdPoint"); - } - - if (fbUseSts) { - LOG(info) << "CA MC Module: initializing branches for STS"; - fpStsPoints = mcManager->InitBranch("StsPoint"); - } - - if (fbUseMuch) { - LOG(info) << "CA MC Module: initializing branches for MuCh"; - fpMuchPoints = mcManager->InitBranch("MuchPoint"); - } - - if (fbUseTrd) { - LOG(info) << "CA MC Module: initializing branches for TRD"; - fpTrdPoints = mcManager->InitBranch("TrdPoint"); - } - - if (fbUseTof) { - LOG(info) << "CA MC Module: initializing branches for TOF"; - fpTofPoints = mcManager->InitBranch("TofPoint"); - } - - // Check initialization - this->CheckInit(); - - LOG(info) << "Initializing CA tracking Monte-Carlo module... \033[1;32mDone!\033[0m"; - return true; -} -catch (const std::logic_error& error) { - LOG(info) << "Initializing CA tracking Monte-Carlo module... \033[1;31mFailed\033[0m\nReason: " << error.what(); - return false; -} - -// --------------------------------------------------------------------------------------------------------------------- -// -void CbmCaPerformance::ProcessEvent(CbmEvent* pEvent) -{ - assert(pEvent); - std::cout << "\033[1;32mProcessing performance event\033[0m\n"; -} - -// --------------------------------------------------------------------------------------------------------------------- -// -void CbmCaPerformance::Finish() { std::cout << "\033[1;32mFinishing performance\033[0m\n"; } - - -// *********************** -// ** Accessors ** -// *********************** - -// --------------------------------------------------------------------------------------------------------------------- -// -void CbmCaPerformance::SetDetector(L1DetectorID detID, bool flag) -{ - switch (detID) { - case L1DetectorID::kMvd: fbUseMvd = flag; break; - case L1DetectorID::kSts: fbUseSts = flag; break; - case L1DetectorID::kMuch: fbUseMuch = flag; break; - case L1DetectorID::kTrd: fbUseTrd = flag; break; - case L1DetectorID::kTof: fbUseTof = flag; break; - } -} - - -// ******************************* -// ** Utility functions ** -// ******************************* - -// --------------------------------------------------------------------------------------------------------------------- -// -void CbmCaPerformance::CheckInit() const -{ - // Check parameters - if (!fpParameters) { throw std::logic_error("Tracking parameters object was not defined"); } - - // Check detectors initialization - if (fbUseMvd) { - if (!fpMvdPoints) { throw std::logic_error("MC points unavailable for MVD"); } - } - - if (fbUseSts) { - if (!fpStsPoints) { throw std::logic_error("MC points unavailable for STS"); } - } - - if (fbUseMuch) { - if (!fpMuchPoints) { throw std::logic_error("MC points unavailable for MuCh"); } - } - - if (fbUseTrd) { - if (!fpTrdPoints) { throw std::logic_error("MC points unavailable for TRD"); } - } - - if (fbUseTof) { - if (!fpTofPoints) { throw std::logic_error("MC points unavailable for TOF"); } - } -} diff --git a/reco/L1/CbmCaPerformance.h b/reco/L1/CbmCaPerformance.h deleted file mode 100644 index 2f721ecf870dc706fbda19548d160c41071c9199..0000000000000000000000000000000000000000 --- a/reco/L1/CbmCaPerformance.h +++ /dev/null @@ -1,123 +0,0 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/// \file CbmCaPerformance.h -/// \brief CA Tracking performance interface for CBM (header) -/// \since 23.09.2022 -/// \author S.Zharko <s.zharko@gsi.de> - -#ifndef CbmCaPerformance_h -#define CbmCaPerformance_h 1 - -#include "CaToolsMcDataManager.h" -#include "CaToolsPerformance.h" - - -class CbmEvent; -class CbmMCDataObject; -class CbmMCDataArray; -class CbmMCEventList; -enum class L1DetectorID; - - -/// Class CbmCaPerformcance is an interface to communicate between -class CbmCaPerformance { -public: - // ***************************************** - // ** Constructors and destructor ** - // ***************************************** - - /// Default constructor - CbmCaPerformance() = default; - - /// Destructor - ~CbmCaPerformance() = default; - - /// Copy constructor - CbmCaPerformance(const CbmCaPerformance&) = delete; - - /// Move constructor - CbmCaPerformance(CbmCaPerformance&&) = delete; - - /// Copy assignment operator - CbmCaPerformance& operator=(const CbmCaPerformance&) = delete; - - /// Move assignment operator - CbmCaPerformance& operator=(CbmCaPerformance&&) = delete; - - - // ***************************************** - // ** Action definition functions ** - // ***************************************** - - /// Defines performance action in the beginning of the run - /// \return Success flag - bool Init(); - - /// Defines performance action in each event or timeslice - /// \param pEvent Pointer to a current CbmEvent - void ProcessEvent(CbmEvent* pEvent); - - /// Defines performance action in the end of the run - void Finish(); - - - // *********************** - // ** Accessors ** - // *********************** - - /// Registers pointer to the tracking parameters object - void SetParameters(const L1Parameters* pParameters) { fpParameters = pParameters; } - - /// Sets used detector subsystems - /// \param detID Id of detector - /// \param flag Flag: true - detector is used - void SetDetector(L1DetectorID detID, bool flag); - - -private: - // ******************************* - // ** Utility functions ** - // ******************************* - - /// Checks class initialization. Throws std::logic_error, if initialization is incomplete at initialization call - void CheckInit() const; - - - // ******************************* - // ** Utility variables ** - // ******************************* - - ca::tools::Performance fPerformance = {}; ///< Instance of the internal performance object - ca::tools::McDataManager fMcDataManager = {}; ///< Instance of the MC data manager - const L1Parameters* fpParameters = nullptr; ///< Pointer to tracking parameters object - - bool fbUseMvd = false; - bool fbUseSts = false; - bool fbUseMuch = false; - bool fbUseTrd = false; - bool fbUseTof = false; - - - // ********************************* - // ** Input data branches ** - // ********************************* - - // Mc-event - // const CbmMCEventList* fpEventList = nullptr; ///< MC event list - // const CbmMCDataObject* fpMcEventHeader = nullptr; ///< MC event header - // const CbmMCDataArray* fpMcTracks = nullptr; ///< MC tracks - - // Mc-points - const CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points container - const CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points container - const CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points container - const CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points container - const CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points container - - - // Matching information -}; - -#endif // CbmCaPerformance_h diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index 5ba72c87a67b1fb605924171cc610844a6b3d51e..66926522b6242d1bcc763c5b5b58c11b1cb092b8 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -1004,27 +1004,11 @@ void CbmL1::Reconstruct(CbmEvent* event) if ((fPerformance) && (fSTAPDataMode < 2)) { InputPerformance(); } - // TODO: Remove this code (S.Zharko) - // for (unsigned int iH = 0; iH < (*fpAlgo->vHits).size(); ++iH) { - //#ifdef USE_EVENT_NUMBER - // L1Hit& h = const_cast<L1Hit&>((*fpAlgo->vHits)[iH]); - // h.n = -1; - //#endif - // if (fvExternalHits[iH].mcPointIds.size() == 0) continue; - //#ifdef USE_EVENT_NUMBER - // const CbmL1MCPoint& mcp = fvMCPoints[fvExternalHits[iH].mcPointIds[0]]; - // h.n = mcp.event; - //#endif - // } // output performance if (fPerformance) { if (fVerbose > 1) { cout << "Performance..." << endl; } TrackMatch(); - } - - - if (fPerformance) { EfficienciesPerformance(); HistoPerformance(); TrackFitPerformance(); @@ -1032,20 +1016,8 @@ void CbmL1::Reconstruct(CbmEvent* event) // TimeHist(); /// WriteSIMDKFData(); } - if (fVerbose > 1) { cout << "End of L1" << endl; } - - static bool ask = 0; - char symbol; - if (ask) { - std::cout << std::endl << "L1 run (any/r/q) > "; - do { - std::cin.get(symbol); - if (symbol == 'r') ask = false; - if ((symbol == 'e') || (symbol == 'q')) exit(0); - } while (symbol != '\n'); - } - - // } + if (fVerbose > 1) { LOG(info) << "Tracking performance... done"; } + if (fVerbose > 1) { LOG(info) << "End of L1"; } } // ----- Finish CbmStsFitPerformanceTask task ----------------------------- diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h index 35c65a33b224b8f335411db907a6f93ae7512503..5cd6d31c20b35fa42c6b0bb4b557ea9efc42ef88 100644 --- a/reco/L1/CbmL1.h +++ b/reco/L1/CbmL1.h @@ -22,6 +22,8 @@ #define _CbmL1_h_ +#include "CbmCaMCModule.h" +#include "CbmL1DetectorID.h" #include "CbmL1Hit.h" #include "CbmL1MCPoint.h" #include "CbmL1MCTrack.h" @@ -66,7 +68,6 @@ class CbmL1Track; class CbmL1MCTrack; class KFTopoPerformance; class CbmMCDataObject; - class CbmEvent; class TProfile2D; class TNtuple; @@ -103,18 +104,6 @@ namespace std } // namespace std -/// Enumeration for the detector subsystems used in L1 tracking -/// It is important for the subsystems to be specified in the actual order. The order is used -/// for the L1Station array filling. -/// Note: L1DetectorID has a forward declaration in L1InitManager.h and L1BaseStationInfo.h -enum class L1DetectorID -{ - kMvd, - kSts, - kMuch, - kTrd, - kTof -}; // TODO: insert documentation! (S.Zh.) // @@ -129,7 +118,7 @@ public: // ** Types definition ** // ********************** - using DFSET = std::set<std::pair<int, int>>; // Why std::set<std::pair<>> instead of std::map<>? + using DFSET = std::set<std::pair<int, int>>; // ************************** // ** Friends classes list ** @@ -158,20 +147,20 @@ public: /// \param name Name of the task /// \param verbose Verbosity level /// \param performance Performance run flag: - /// 0: run without performance measurement - /// 1: standard efficiency definition - /// 2: QA efficiency definition + /// - #0 run without performance measurement + /// - #1 standard efficiency definition + /// - #2 QA efficiency definition /// \param dataMode Option to work with files for the standalone mode - /// 0: standalone mode is not used - /// 1: data for standalone mode is written to configuration file (currently does not work) - /// 2: tracking runs in standalone mode using configuration file (currently does not work) - /// 3: data is written and read (currently does not work) + /// - #0 standalone mode is not used + /// - #1 data for standalone mode is written to configuration file (currently does not work) + /// - #2 tracking runs in standalone mode using configuration file (currently does not work) + /// - #3 data is written and read (currently does not work) /// \param dataDir Name of directory for configuration file /// \param findParticleMode Find particle utility mode - /// 0: FindParticles is not used - /// 1: All MC particles are reconstructable - /// 2: MC particles are reconstructable if created from reconstructable tracks - /// 3: MC particles are reconstructable if created from reconstructed tracks + /// - #0 FindParticles is not used + /// - #1 All MC particles are reconstructable + /// - #2 MC particles are reconstructable if created from reconstructable tracks + /// - #3: MC particles are reconstructable if created from reconstructed tracks CbmL1(const char* name, Int_t verbose = 1, Int_t performance = 0, int dataMode = 0, const TString& dataDir = "./", int findParticleMode = 0); @@ -302,12 +291,10 @@ public: // void SetGhostSuppression( Bool_t b ){ fGhostSuppression= b; } // void SetDetectorEfficiency( Double_t eff ){ fDetectorEfficiency = eff; } - /// Reconstructs an event + /// Reconstructs tracks in a given event /// \param event Pointer to current CbmEvent object void Reconstruct(CbmEvent* event = nullptr); - // bool ReadMCDataFromFile(const char work_dir[100], const int maxNEvent, const int iVerbose); - // vector<CbmL1MCPoint> fvMCPoints; // static bool compareZ(const int &a, const int &b ); inline double Get_Z_vMCPoint(int a) const { return fvMCPoints[a].z; } @@ -369,7 +356,7 @@ private: /// Matches hit with MC point /// \tparam DetId Detector ID /// \param iHit External index of hit - /// \return tuple of MC point index in fvMCPoints array + /// \return MC-point index in fvMCPoints array template<L1DetectorID DetId> int MatchHitWithMc(int iHit) const; @@ -458,14 +445,18 @@ private: // ** Member variables list ** // *************************** + + L1InitManager fInitManager; ///< Tracking parameters data manager + L1IODataManager fIODataManager; ///< Input-output data manager + + std::unique_ptr<CbmCaMCModule> fpMCModule = nullptr; ///< MC-module for tracking + + public: // ** Basic data members ** L1Algo* fpAlgo = nullptr; ///< Pointer to the L1 track finder algorithm - L1InitManager fInitManager; ///< Tracking parameters data manager - L1IODataManager fIODataManager; ///< Input-output data manager - bool fMissingHits = false; ///< Turns on several ad-hock settings for "mcbm_beam_2021_07_surveyed.100ev" setup L1Algo::TrackingMode fTrackingMode = L1Algo::TrackingMode::kSts; ///< Tracking mode: kSts, kMcbm or kGlobal @@ -491,7 +482,6 @@ private: int fNpointsTofAll = 0; ///< Number of MC points for TOF L1Vector<CbmL1MCPoint> fvMCPoints = {"CbmL1::fvMCPoints"}; ///< Container of MC points - L1Vector<int> fvMCPointIndexesTs = {"CbmL1::fvMCPointIndexesTs"}; ///< Indexes of MC points in TS int fNStations = 0; ///< number of total active detector stations int fNMvdStations = 0; ///< number of active MVD stations @@ -603,7 +593,6 @@ private: TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches TClonesArray* fpMuchDigiMatches = nullptr; ///< Array of MuCh digi matches (NOTE: currently unsused) - vector<vector<int>> fTofPointToTrack; ///< L1Vector<CbmL1MCTrack> fvMCTracks = {"CbmL1::fvMCTracks"}; ///< Array of MC tracks L1Vector<int> fvHitPointIndexes = {"CbmL1::fvHitPointIndexes"}; ///< Indexes of MC points vs. hit index diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h new file mode 100644 index 0000000000000000000000000000000000000000..6bffea7026eb95166b225f8b28b44914a3df189f --- /dev/null +++ b/reco/L1/CbmL1DetectorID.h @@ -0,0 +1,26 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmL1DetectorID.h +/// \brief Implementation of L1DetectorID enum class for CBM +/// \author S.Zharko +/// \data 01.12.2022 + +#ifndef CbmL1DetectorID_h +#define CbmL1DetectorID_h 1 + +/// Enumeration for the detector subsystems used in L1 tracking +/// It is important for the subsystems to be specified in the actual order. The order is used +/// for the L1Station array filling. +/// Note: L1DetectorID has a forward declaration in L1InitManager.h and L1BaseStationInfo.h +enum class L1DetectorID +{ + kMvd, + kSts, + kMuch, + kTrd, + kTof +}; + +#endif // CbmL1DetectorID_h diff --git a/reco/L1/CbmL1MCPoint.h b/reco/L1/CbmL1MCPoint.h index 660d1635bbe02ab4f7926ed97967282127468a11..6cab8ceb58ea1c6448433fd29c8a32846f46c3ac 100644 --- a/reco/L1/CbmL1MCPoint.h +++ b/reco/L1/CbmL1MCPoint.h @@ -21,6 +21,10 @@ #ifndef CbmL1MCPoint_H #define CbmL1MCPoint_H +#include <iomanip> +#include <sstream> +#include <string> + #include "L1Vector.h" struct CbmL1MCPoint { @@ -66,8 +70,60 @@ struct CbmL1MCPoint { int pointId = -1; int file = -1; int event = -1; - int trackId = -1; L1Vector<int> hitIds {"CbmL1MCPoint::hitIds"}; // indices of CbmL1Hits in L1->vStsHits array + + /// Temporary log function for debugging + std::string ToString(int verbose, bool printHeader = false) const + { + if (verbose < 1) { return std::string(); } + + std::stringstream msg; + if (printHeader) { + if (verbose > 0) { + msg << std::setw(10) << std::setfill(' ') << "track ID" << ' '; + msg << std::setw(10) << std::setfill(' ') << "mother ID" << '|'; + msg << std::setw(10) << std::setfill(' ') << "station ID" << '|'; + msg << std::setw(10) << std::setfill(' ') << "PDG" << ' '; + msg << std::setw(10) << std::setfill(' ') << "m [GeV/c2]" << ' '; + msg << std::setw(10) << std::setfill(' ') << "q [e]" << '|'; + msg << std::setw(14) << std::setfill(' ') << "t [ns]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "x [cm]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "y [cm]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "z [cm]" << '|'; + if (verbose > 1) { + msg << std::setw(14) << std::setfill(' ') << "zIn [cm]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "zOut [cm]" << '|'; + msg << std::setw(14) << std::setfill(' ') << "p [GeV/c]" << '|'; + msg << std::setw(10) << std::setfill(' ') << "point ID" << ' '; + msg << std::setw(10) << std::setfill(' ') << "event ID" << ' '; + msg << std::setw(10) << std::setfill(' ') << "file ID" << ' '; + } + } + } + else { + if (verbose > 0) { + msg << std::setw(10) << std::setfill(' ') << ID << ' '; + msg << std::setw(10) << std::setfill(' ') << mother_ID << '|'; + msg << std::setw(10) << std::setfill(' ') << iStation << '|'; + msg << std::setw(10) << std::setfill(' ') << pdg << ' '; + msg << std::setw(10) << std::setfill(' ') << mass << ' '; + msg << std::setw(10) << std::setfill(' ') << q << '|'; + msg << std::setw(14) << std::setfill(' ') << time << ' '; + msg << std::setw(14) << std::setfill(' ') << x << ' '; + msg << std::setw(14) << std::setfill(' ') << y << ' '; + msg << std::setw(14) << std::setfill(' ') << z << '|'; + if (verbose > 1) { + msg << std::setw(14) << std::setfill(' ') << zIn << ' '; + msg << std::setw(14) << std::setfill(' ') << zOut << '|'; + msg << std::setw(14) << std::setfill(' ') << p << '|'; + msg << std::setw(10) << std::setfill(' ') << pointId << ' '; + msg << std::setw(10) << std::setfill(' ') << event << ' '; + msg << std::setw(10) << std::setfill(' ') << file << ' '; + } + } + } + return msg.str(); + } }; #endif diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx index 6e52c5c939b322100184e09ecd1dc6b28542a01b..e3330b8a100821c1e190b53f7cf970337780652e 100644 --- a/reco/L1/CbmL1MCTrack.cxx +++ b/reco/L1/CbmL1MCTrack.cxx @@ -79,20 +79,6 @@ void CbmL1MCTrack::Init() CalculateIsReconstructable(); } - -float CbmL1MCTrack::Fraction_MC() -{ - // if ( Points.size() == 0 ) return 0; - - CbmL1* L1 = CbmL1::Instance(); - int counter = 0; - for (unsigned int iP = 0; iP < Points.size(); iP++) { - if (L1->fvMCPointIndexesTs[Points[iP]] > 0) counter++; - }; - return (Points.size() > 0) ? float(counter) / float(Points.size()) : 0; -} - - void CbmL1MCTrack::CalculateMCCont() { CbmL1* L1 = CbmL1::Instance(); @@ -194,8 +180,6 @@ void CbmL1MCTrack::CalculateIsReconstructable() // f &= (maxNStaHits <= 4); f &= (maxNStaMC <= 4); // f &= (maxNSensorMC <= 1); - if (L1->fPerformance == 4) - isReconstructable = f & (nMCContStations >= CbmL1Constants::MinNStations) & (Fraction_MC() > 0.5); if (L1->fPerformance == 3) isReconstructable = f & (nMCContStations >= CbmL1Constants::MinNStations); // L1-MC if (L1->fPerformance == 2) isReconstructable = f & (nStations >= CbmL1Constants::MinNStations); // QA definition if (L1->fPerformance == 1) diff --git a/reco/L1/CbmL1MCTrack.h b/reco/L1/CbmL1MCTrack.h index 330a4a20a6a10baeb8f03afa3c780b10be133b8e..cd11625bd325f91a6920507e38f9f51d84c005ed 100644 --- a/reco/L1/CbmL1MCTrack.h +++ b/reco/L1/CbmL1MCTrack.h @@ -49,7 +49,6 @@ public: int NHitContStations() const { return nHitContStations; } int NMCStations() const { return nMCStations; } int NMCContStations() const { return nMCContStations; } - float Fraction_MC(); void Init(); diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index 04ed4ea15da821acc5b6607b69cbc9fa17f8649c..817f91c73aa6b7fb85d5542479baafa0e65c44c0 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -256,7 +256,6 @@ int CbmL1::MatchHitWithMc<L1DetectorID::kTrd>(int iHit) const auto itPoint = fmMCPointsLinksMap.find(CbmL1LinkKey(iIndex, iEvent, iFile)); assert(itPoint != fmMCPointsLinksMap.cend()); - if (itPoint == fmMCPointsLinksMap.cend()) return iPoint; iPoint = itPoint->second; } } @@ -367,8 +366,6 @@ void CbmL1::ReadEvent(CbmEvent* event) // Fill fvMCPoints and fmMCPointsLinksMap fvMCPoints.clear(); fvMCPoints.reserve(5 * fvMCTracks.size() * fNStations); - fvMCPointIndexesTs.clear(); - fvMCPointIndexesTs.reserve(fvMCPoints.capacity()); fNpointsMvdAll = 0; fNpointsStsAll = 0; @@ -420,7 +417,6 @@ void CbmL1::ReadEvent(CbmEvent* event) fmMCPointsLinksMap[CbmL1LinkKey(iMC, iEvent, iFile)] = fvMCPoints.size(); fvMCPoints.push_back(MC); - fvMCPointIndexesTs.push_back(0); fNpointsMvd++; } } @@ -459,7 +455,6 @@ void CbmL1::ReadEvent(CbmEvent* event) fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size()); fmMCPointsLinksMap[CbmL1LinkKey(iMC + fNpointsMvdAll, iEvent, iFile)] = fvMCPoints.size(); fvMCPoints.push_back(MC); - fvMCPointIndexesTs.push_back(0); fNpointsSts++; } } @@ -478,17 +473,15 @@ void CbmL1::ReadEvent(CbmEvent* event) const L1Station* sta = fpAlgo->GetParameters()->GetStations().begin(); for (Int_t iSt = 0; iSt < fNMuchStationsGeom; iSt++) { int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kMuch); - assert(iStActive != -1); if (MC.z > sta[iStActive].fZ[0] - 2.5) { MC.iStation = iStActive; } } - assert(MC.iStation >= 0); + if (MC.iStation < 0) { continue; } // Reject MC points from inactive stations; auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile)); assert(itTrack != fmMCTracksLinksMap.cend()); MC.ID = itTrack->second; fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size()); fmMCPointsLinksMap[CbmL1LinkKey(iMC + fNpointsMvdAll + +fNpointsStsAll, iEvent, iFile)] = fvMCPoints.size(); fvMCPoints.push_back(MC); - fvMCPointIndexesTs.push_back(0); fNpointsMuch++; } } @@ -504,11 +497,10 @@ void CbmL1::ReadEvent(CbmEvent* event) const L1Station* sta = fpAlgo->GetParameters()->GetStations().begin(); for (Int_t iSt = 0; iSt < fNTrdStationsGeom; iSt++) { int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kTrd); - - assert(iStActive != -1); + if (iStActive < 0) { continue; } if (MC.z > sta[iStActive].fZ[0] - 20.0) { MC.iStation = iStActive; } } - assert(MC.iStation >= 0); + if (MC.iStation < 0) { continue; } // Reject MC points from inactive stations auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile)); assert(itTrack != fmMCTracksLinksMap.cend()); MC.ID = itTrack->second; @@ -516,7 +508,6 @@ void CbmL1::ReadEvent(CbmEvent* event) fmMCPointsLinksMap[CbmL1LinkKey(iMC + fNpointsMvdAll + fNpointsStsAll + fNpointsMuchAll, iEvent, iFile)] = fvMCPoints.size(); fvMCPoints.push_back(MC); - fvMCPointIndexesTs.push_back(0); fNpointsTrd++; } } @@ -526,20 +517,20 @@ void CbmL1::ReadEvent(CbmEvent* event) firstTofPoint = fvMCPoints.size(); if (fUseTOF && fpTofPoints) { - vector<float> TofPointToTrackdZ[fNTofStations]; + vector<float> vTofPointToTrackdZ[fNTofStationsGeom]; // active stations! - fTofPointToTrack.resize(fNTofStations); + vector<vector<int>> vTofPointToTrack(fNTofStationsGeom); // active stations! - for (Int_t i = 0; i < fNTofStations; i++) { + for (Int_t i = 0; i < fNTofStationsGeom; i++) { - fTofPointToTrack[i].resize(fvMCTracks.size(), 0); - TofPointToTrackdZ[i].resize(fvMCTracks.size()); + vTofPointToTrack[i].resize(fvMCTracks.size(), 0); + vTofPointToTrackdZ[i].resize(fvMCTracks.size()); } - for (int iSt = 0; iSt < fNTofStations; iSt++) - for (unsigned int i = 0; i < TofPointToTrackdZ[iSt].size(); i++) { - TofPointToTrackdZ[iSt][i] = 100000; - fTofPointToTrack[iSt][i] = -1; + for (int iSt = 0; iSt < fNTofStationsGeom; iSt++) + for (unsigned int i = 0; i < vTofPointToTrackdZ[iSt].size(); i++) { + vTofPointToTrackdZ[iSt][i] = 100000; + vTofPointToTrack[iSt][i] = -1; } @@ -557,56 +548,60 @@ void CbmL1::ReadEvent(CbmEvent* event) } } + for (Int_t iMC = 0; iMC < fpTofPoints->Size(iFile, iEvent); iMC++) { - if (isTofPointMatched[iMC] == 0) continue; + if (isTofPointMatched[iMC] == 0) { continue; } CbmL1MCPoint MC; if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 4)) { - auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile)); - assert(itTrack != fmMCTracksLinksMap.cend()); - int iTrack = itTrack->second; - MC.iStation = -1; const L1Station* sta = fpAlgo->GetParameters()->GetStations().begin(); - float dist = 1000; - int iSta = -1; + float dist = 1000; + int iSta = -1; + int iStaGlob = -1; for (int iSt = 0; iSt < fNTofStationsGeom; iSt++) { int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kTof); if (iStActive == -1) { continue; } if (fabs(MC.z - sta[iStActive].fZ[0]) < dist) { - dist = fabs(MC.z - sta[iStActive].fZ[0]); - iSta = iSt; + dist = fabs(MC.z - sta[iStActive].fZ[0]); + iSta = iSt; + iStaGlob = iStActive; } } + MC.iStation = fpAlgo->GetParameters()->GetStationIndexActive(iSta, L1DetectorID::kTof); assert(MC.iStation >= 0); + + auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile)); + assert(itTrack != fmMCTracksLinksMap.cend()); + int iTrack = itTrack->second; + if (iSta >= 0) { - if (fabs(sta[iSta].fZ[0] - MC.z) < TofPointToTrackdZ[iSta][iTrack]) { - fTofPointToTrack[iSta][iTrack] = iMC; - TofPointToTrackdZ[iSta][iTrack] = fabs(sta[iSta].fZ[0] - MC.z); + if (fabs(sta[iStaGlob].fZ[0] - MC.z) < vTofPointToTrackdZ[iSta][iTrack]) { + vTofPointToTrack[iSta][iTrack] = iMC; + vTofPointToTrackdZ[iSta][iTrack] = fabs(sta[iStaGlob].fZ[0] - MC.z); } } } } - for (int iTofSta = 0; iTofSta < fNTofStations; iTofSta++) - for (unsigned int iTrack = 0; iTrack < fTofPointToTrack[iTofSta].size(); iTrack++) { + for (int iTofSta = 0; iTofSta < fNTofStationsGeom; iTofSta++) + for (unsigned int iTrack = 0; iTrack < vTofPointToTrack[iTofSta].size(); iTrack++) { - if (fTofPointToTrack[iTofSta][iTrack] < 0) continue; + if (vTofPointToTrack[iTofSta][iTrack] < 0) continue; CbmL1MCPoint MC; - if (!ReadMCPoint(&MC, fTofPointToTrack[iTofSta][iTrack], iFile, iEvent, 4)) { + if (!ReadMCPoint(&MC, vTofPointToTrack[iTofSta][iTrack], iFile, iEvent, 4)) { - MC.iStation = (fNMvdStations + fNStsStations + fNMuchStations + fNTrdStations + iTofSta); + MC.iStation = fpAlgo->GetParameters()->GetStationIndexActive(iTofSta, L1DetectorID::kTof); fvMCTracks[iTrack].Points.push_back_no_warning(fvMCPoints.size()); MC.ID = iTrack; int iMC = - fTofPointToTrack[iTofSta][iTrack] + fNpointsMvdAll + +fNpointsStsAll + fNpointsMuchAll + fNpointsTrdAll; + vTofPointToTrack[iTofSta][iTrack] + fNpointsMvdAll + +fNpointsStsAll + fNpointsMuchAll + fNpointsTrdAll; fmMCPointsLinksMap[CbmL1LinkKey(iMC, iEvent, iFile)] = fvMCPoints.size(); fvMCPoints.push_back(MC); - fvMCPointIndexesTs.push_back(0); fNpointsTof++; } } @@ -626,18 +621,17 @@ void CbmL1::ReadEvent(CbmEvent* event) } } //iTr if (fVerbose >= 10) cout << "Points in fvMCTracks are sorted." << endl; - } //fPerformance /* - * MC hits and tracks gathering: END + * MC points and tracks gathering: END */ /* * Measured hits gathering: START * * In this section the measured hits from different detector subsystems are reformatted according to TmpHit structure - * (NOTE: independent from particular detector design) and then pushed to the tmpHit vector. In the performance study mode + * (NOTE: independently from particular detector design) and then pushed to the tmpHit vector. In the performance study mode * matching with MC points is done */ diff --git a/reco/L1/L1Algo/L1Constants.h b/reco/L1/L1Algo/L1Constants.h index c11681ad6a5874a46f9cc4d00324a33b9a5d9b99..4917146e0b47ff7fcc5d11c618f8be9d26e29d19 100644 --- a/reco/L1/L1Algo/L1Constants.h +++ b/reco/L1/L1Algo/L1Constants.h @@ -17,6 +17,7 @@ /// Namespace contains compile-time constants definition for the L1 tracking algorithm namespace L1Constants { + using cbm::algo::ca::fvec; /// Array sizes namespace size @@ -88,16 +89,6 @@ namespace L1Constants constexpr int kAlignment = 16; ///< Default alignment of data (bytes) } // namespace misc - /// NoInit constants (aliases) - namespace noin - { - constexpr float kF = L1NaN::SetNaN<float>(); - constexpr double kD = L1NaN::SetNaN<double>(); - constexpr int k32I = L1NaN::SetNaN<int>(); - constexpr unsigned int k32U = L1NaN::SetNaN<unsigned int>(); - } // namespace noin - - // Colors of terminal log namespace clrs { diff --git a/reco/L1/L1Algo/L1Material.h b/reco/L1/L1Algo/L1Material.h index 319be73ebd7e2462e27a82afe939d9b565369f82..1727d3c363d3ef73de54541b53e7bc941501e40f 100644 --- a/reco/L1/L1Algo/L1Material.h +++ b/reco/L1/L1Algo/L1Material.h @@ -15,6 +15,7 @@ #include "L1Def.h" #include "L1NaN.h" #include "L1SimdSerializer.h" +#include "L1Undef.h" /// Class L1Material describes a map of station thickness in units of radiation length (X0) to the specific point in XY plane class L1Material { @@ -85,11 +86,10 @@ public: void Repare(); private: - int fNbins {L1NaN::SetNaN<decltype(fNbins)>()}; ///< Number of rows (columns) in the material budget table - float fRmax {L1NaN::SetNaN<decltype(fRmax)>()}; ///< Size of the station in x and y dimensions [cm] - float fFactor { - L1NaN::SetNaN<decltype(fFactor)>()}; ///< Factor used in the recalculation of point coordinates to row/column id - std::vector<float> fTable {}; ///< Material budget table + int fNbins = undef::kI32; ///< Number of rows (columns) in the material budget table + float fRmax = undef::kF; ///< Size of the station in x and y dimensions [cm] + float fFactor = undef::kF; ///< Factor used in the recalculation of point coordinates to row/column id + std::vector<float> fTable {}; ///< Material budget table /// Serialization function friend class boost::serialization::access; diff --git a/reco/L1/L1Algo/L1NaN.h b/reco/L1/L1Algo/L1NaN.h index 81b25073adb1fdbc74e0426f581b7d7c490f9637..129cd32ca05884708537cac7cb38eefab1ea9a77 100644 --- a/reco/L1/L1Algo/L1NaN.h +++ b/reco/L1/L1Algo/L1NaN.h @@ -19,7 +19,10 @@ #include "CaSimd.h" #include "L1Def.h" -using namespace cbm::algo::ca; // TODO: remove "using" from headers + +using cbm::algo::ca::fmask; +using cbm::algo::ca::fscal; // TODO: remove "using" from headers +using cbm::algo::ca::fvec; // TODO: remove "using" from headers /// Namespace L1NaN defines functions to set variables to NaN and check wether they are NaN or not /// diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h index baa4638624718b46d2bc39c52cce1f1302587999..0b86797812e65485c35e5981bab2afaff63d73dd 100644 --- a/reco/L1/L1Algo/L1Station.h +++ b/reco/L1/L1Algo/L1Station.h @@ -10,9 +10,9 @@ #include "L1Field.h" #include "L1Material.h" -#include "L1NaN.h" #include "L1SimdSerializer.h" #include "L1UMeasurementInfo.h" +#include "L1Undef.h" #include "L1Utils.h" #include "L1XYMeasurementInfo.h" @@ -21,14 +21,13 @@ /// class L1Station { public: - int type {L1NaN::SetNaN<decltype(type)>()}; - int timeInfo {L1NaN::SetNaN<decltype(timeInfo)>()}; ///< flag: if time information can be used - int fieldStatus {L1NaN::SetNaN<decltype( - fieldStatus)>()}; ///< flag: 1 - station is INSIDE the field, 0 - station is OUTSIDE the field (replace with enum) - fvec fZ {L1NaN::SetNaN<decltype(fZ)>()}; ///< z position of station [cm] - fvec fZthick {L1NaN::SetNaN<decltype(fZthick)>()}; ///< z thickness of the station [cm] - fvec Rmin {L1NaN::SetNaN<decltype(Rmin)>()}; ///< min radius of the station [cm] - fvec Rmax {L1NaN::SetNaN<decltype(Rmax)>()}; ///< max radius of the station [cm] + int type = undef::kI32; + int timeInfo = undef::kI32; ///< flag: if time information can be used + int fieldStatus = undef::kI32; ///< flag: 1 - station is INSIDE the field, 0 - station is OUTSIDE the field + fvec fZ = undef::kFvc; ///< z position of station [cm] + fvec fZthick = undef::kFvc; ///< z thickness of the station [cm] + fvec Rmin = undef::kFvc; ///< min radius of the station [cm] + fvec Rmax = undef::kFvc; ///< max radius of the station [cm] L1FieldSlice fieldSlice {}; L1UMeasurementInfo frontInfo {}; diff --git a/reco/L1/L1Algo/L1UMeasurementInfo.h b/reco/L1/L1Algo/L1UMeasurementInfo.h index 5405719c6036efe7ca197fc98a1a125d761658f4..51e05f7040ec1dbdef434e9321fcececa4c98091 100644 --- a/reco/L1/L1Algo/L1UMeasurementInfo.h +++ b/reco/L1/L1Algo/L1UMeasurementInfo.h @@ -10,13 +10,14 @@ #include "L1Def.h" #include "L1NaN.h" #include "L1SimdSerializer.h" +#include "L1Undef.h" #include "L1Utils.h" class L1UMeasurementInfo { public: - fvec cos_phi {L1NaN::SetNaN<decltype(cos_phi)>()}; - fvec sin_phi {L1NaN::SetNaN<decltype(sin_phi)>()}; + fvec cos_phi = undef::kFvc; + fvec sin_phi = undef::kFvc; /// String representation of class contents /// \param indentLevel number of indent characters in the output @@ -26,7 +27,7 @@ public: void CheckConsistency() const; /// Checks, if the fields are NaN - bool IsNaN() const { return L1NaN::IsNaN(cos_phi) || L1NaN::IsNaN(sin_phi); } + bool IsNaN() const { return isnan(cos_phi).isNotEmpty() || isnan(sin_phi).isNotEmpty(); } /// Serialization function friend class boost::serialization::access; diff --git a/reco/L1/L1Algo/L1XYMeasurementInfo.h b/reco/L1/L1Algo/L1XYMeasurementInfo.h index 35cfca20f85a9995c86b9eb428e1719b85908fa8..aa48af6501128dbfb15611d2920e571e5b26bfab 100644 --- a/reco/L1/L1Algo/L1XYMeasurementInfo.h +++ b/reco/L1/L1Algo/L1XYMeasurementInfo.h @@ -10,12 +10,13 @@ #include "L1Def.h" #include "L1NaN.h" #include "L1SimdSerializer.h" +#include "L1Undef.h" class L1XYMeasurementInfo { public: - fvec C00 {L1NaN::SetNaN<decltype(C00)>()}; - fvec C10 {L1NaN::SetNaN<decltype(C10)>()}; - fvec C11 {L1NaN::SetNaN<decltype(C11)>()}; + fvec C00 = undef::kFvc; + fvec C10 = undef::kFvc; + fvec C11 = undef::kFvc; /// Consistency checker void CheckConsistency() const; diff --git a/reco/L1/catools/CaToolsLinkKey.h b/reco/L1/catools/CaToolsLinkKey.h index aecc629ec4ed9dfed9a0a7ae438e9dec158c599b..78c91dea0531b9029ce6ff6c117d9f388b572901 100644 --- a/reco/L1/catools/CaToolsLinkKey.h +++ b/reco/L1/catools/CaToolsLinkKey.h @@ -12,6 +12,8 @@ #include <boost/functional/hash.hpp> +#include "L1Undef.h" + namespace ca::tools { struct LinkKey { @@ -27,9 +29,9 @@ namespace ca::tools return lhs.fFile == rhs.fFile && lhs.fEvent == rhs.fEvent && lhs.fIndex == rhs.fIndex; } - int fIndex = -1; ///< Index of MC point/track in external data structures - int fEvent = -1; ///< Index of MC event - int fFile = -1; ///< Index of MC file + int fIndex = undef::kI32; + int fEvent = undef::kI32; + int fFile = undef::kI32; }; } // namespace ca::tools diff --git a/reco/L1/catools/CaToolsMCData.cxx b/reco/L1/catools/CaToolsMCData.cxx new file mode 100644 index 0000000000000000000000000000000000000000..865d66103265b57ab917dc316ebfdfde8a18fbf8 --- /dev/null +++ b/reco/L1/catools/CaToolsMCData.cxx @@ -0,0 +1,115 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CaToolsMCData.cxx +/// \brief Data structure for internal tracking MC-information (implementation) +/// \since 23.09.2022 +/// \author S.Zharko <s.zharko@gsi.de> + +#include "CaToolsMCData.h" + +#include <iomanip> +#include <sstream> +#include <utility> // for std::move + +using namespace ca::tools; + +// --------------------------------------------------------------------------------------------------------------------- +// +MCData::MCData() {} + +// --------------------------------------------------------------------------------------------------------------------- +// +MCData::MCData(const MCData& other) + : fvPoints(other.fvPoints) + , fvTracks(other.fvTracks) + , fmPointsLinksMap(other.fmPointsLinksMap) + , fmTracksLinksMap(other.fmPointsLinksMap) +{ +} + +// --------------------------------------------------------------------------------------------------------------------- +// +MCData::MCData(MCData&& other) noexcept { this->Swap(other); } + +// --------------------------------------------------------------------------------------------------------------------- +// +MCData& MCData::operator=(const MCData& other) +{ + if (this != &other) { MCData(other).Swap(*this); } + return *this; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +MCData& MCData::operator=(MCData&& other) noexcept +{ + if (this != &other) { + MCData tmp(std::move(other)); + this->Swap(tmp); + } + return *this; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCData::Swap(MCData& other) noexcept +{ + std::swap(fvPoints, other.fvPoints); + std::swap(fvTracks, other.fvTracks); + std::swap(fmPointsLinksMap, other.fmPointsLinksMap); + std::swap(fmTracksLinksMap, other.fmTracksLinksMap); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCData::AddPoint(const MCPoint& point, int index, int event, int file) +{ + fmPointsLinksMap[LinkKey(index, event, file)] = static_cast<int>(fvPoints.size()); + fvPoints.push_back(point); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCData::AddTrack(const CbmL1MCTrack& track, int index, int event, int file) +{ + fmTracksLinksMap[LinkKey(index, event, file)] = static_cast<int>(fvTracks.size()); + fvTracks.push_back(track); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCData::Clear() +{ + fvPoints.clear(); + fvTracks.clear(); + fmPointsLinksMap.clear(); + fmTracksLinksMap.clear(); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +std::string MCData::ToString(int verbose) const +{ + if (verbose < 1) { return std::string(); } + std::stringstream msg; + msg << "MCData: " << fvTracks.size() << " tracks, " << fvPoints.size() << " points"; + if (verbose > 1) { + using std::setfill; + using std::setw; + const int kNofTracksToPrint = 5; + msg << "\n Track sample (first 5 tracks):\n"; + msg << setw(10) << setfill(' ') << "ID" << ' ' << setw(10) << setfill(' ') << "motherID" << ' ' << setw(10) + << setfill(' ') << "pdg" << ' ' << setw(10) << setfill(' ') << "zVTX [cm]" << ' ' << setw(10) << setfill(' ') + << "px [GeV/c]" << ' ' << setw(10) << setfill(' ') << "py [GeV/c]" << ' ' << setw(10) << setfill(' ') + << "pz [GeV/c]" << '\n'; + for (int i = 0; i < kNofTracksToPrint; ++i) { + msg << setw(10) << setfill(' ') << fvTracks[i].ID << ' ' << setw(10) << setfill(' ') << fvTracks[i].mother_ID + << ' ' << setw(10) << setfill(' ') << fvTracks[i].pdg << ' ' << setw(10) << setfill(' ') << fvTracks[i].z + << ' ' << setw(10) << setfill(' ') << fvTracks[i].px << ' ' << setw(10) << setfill(' ') << fvTracks[i].py + << ' ' << setw(10) << setfill(' ') << fvTracks[i].pz << '\n'; + } + } + return msg.str(); +} diff --git a/reco/L1/catools/CaToolsMCData.h b/reco/L1/catools/CaToolsMCData.h new file mode 100644 index 0000000000000000000000000000000000000000..fba046cbf60b35046b3ce04d0fa2a6c28a0ba5b3 --- /dev/null +++ b/reco/L1/catools/CaToolsMCData.h @@ -0,0 +1,148 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CaToolsMCData.h +/// \brief Data structure for internal tracking MC-information (header) +/// \since 23.09.2022 +/// \author S.Zharko <s.zharko@gsi.de> + +#ifndef CaToolsMCData_h +#define CaToolsMCData_h 1 + +#include "CbmL1MCTrack.h" + +#include <numeric> +#include <string> + +#include "CaToolsLinkKey.h" +#include "CaToolsMCPoint.h" +#include "L1Constants.h" +#include "L1Vector.h" + +enum class L1DetectorID; + +namespace ca::tools +{ + /// This class represents a package for tracking-related data + class MCData { + public: + // ********************************* + // ** Constructors and destructor ** + // ********************************* + + /// Default constructor + MCData(); + + /// Destructor + ~MCData() = default; + + /// Copy constructor + MCData(const MCData& other); + + /// Move constructor + MCData(MCData&& other) noexcept; + + /// Copy assignment operator + MCData& operator=(const MCData& other); + + /// Move assignment operator + MCData& operator=(MCData&& other) noexcept; + + /// Swap method + void Swap(MCData& other) noexcept; + + + /// Adds an MC point + /// \param point MC point object + /// \param index Index of point in a link + /// \param event Event index of a link + /// \param file File index of a link + void AddPoint(const MCPoint& point, int index, int event, int file); + + /// Adds an MC track + /// \param track MC track object + /// \param index Index of point in a link + /// \param event Event index of a link + /// \param file File index of a link + void AddTrack(const CbmL1MCTrack& track, int index, int event, int file); + + /// Clears contents + void Clear(); + + /// Reserves memory for tracks to avoid extra allocations + void ReserveNofTracks(int nTracks) { fvTracks.reserve(nTracks); } + + /// Reserves memory for points to avoid extra allocations + void ReserveNofPoints(int nPoints) { fvPoints.reserve(nPoints); } + + + // ********************** + // ** Getters ** + // ********************** + + /// Gets a reference to MC point by its index + const auto& GetPoint(int idx) const { return fvPoints[idx]; } + + /// Gets a reference to the vector of points + const auto& GetPointContainer() const { return fvPoints; } + + /// Gets a reference to MC track by its index + const auto& GetTrack(int idx) const { return fvTracks[idx]; } + + /// Gets a reference to the vector of tracks + const auto& GetTrackContainer() const { return fvTracks; } + + + /// Finds an MC point global index by given local index, event and file of a link + /// \param index Local index of MC point + /// \param event Event of link + /// \param file File of link + int FindGlobalPointIndex(int index, int event, int file) + { + auto it = fmPointsLinksMap.find(LinkKey(index, event, file)); + return (it != fmPointsLinksMap.cend()) ? it->second : -1; + } + + /// Finds an MC track global index by given local index, event and file of a link + /// \param index Local index of track + /// \param event Index of event + /// \param file Index of file + int FindGlobalTrackIndex(int index, int event, int file) + { + auto it = fmTracksLinksMap.find(LinkKey(index, event, file)); + return (it != fmTracksLinksMap.cend()) ? it->second : -1; + } + + /// Gets number of tracks + int GetNofTracks() const { return fvTracks.size(); } + + /// Gets number of points + int GetNofPoints() { return fvPoints.size(); } + + + // ********************* + // ** Setters ** + // ********************* + + /// Prints an example of tracks and points + /// \param verbose Verbose level: + /// - #0: Nothing is printed + /// - #1: Only numbers of tracks and points are printed + /// - #2: First five tracks and points are printed (partially) + std::string ToString(int verbose = 1) const; + + private: + // ********************** + // ** Member variables ** + // ********************** + + L1Vector<MCPoint> fvPoints = {"ca::tools::MCData::fvMCPoints"}; ///< Container of points + L1Vector<CbmL1MCTrack> fvTracks = {"ca::tools::MCData::fvMCTracks"}; ///< Container of tracks + + std::unordered_map<LinkKey, int> fmPointsLinksMap = {}; ///< MC point internal index vs. link + std::unordered_map<LinkKey, int> fmTracksLinksMap = {}; ///< MC track internal index vs. link + }; +} // namespace ca::tools + +#endif // CaToolsMCData_h diff --git a/reco/L1/catools/CaToolsMCPoint.cxx b/reco/L1/catools/CaToolsMCPoint.cxx new file mode 100644 index 0000000000000000000000000000000000000000..8c53161671075a21b8c65ad282154b7f92f1a1e0 --- /dev/null +++ b/reco/L1/catools/CaToolsMCPoint.cxx @@ -0,0 +1,80 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CaToolsMCPoint.cxx +/// \brief Internal class describing a MC point for CA tracking QA and performance (implementation) +/// \date 18.11.2022 +/// \author S.Zharko <s.zharko@gsi.de> + + +#include "CaToolsMCPoint.h" + +#include <iomanip> +#include <sstream> + +using namespace ca::tools; + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCPoint::SetLinkAddress(int pointID, int eventID, int fileID) +{ + // TODO: SZh 21.11.2022: Test uniqueness of the pointID among different subsystems in link definition + fPointID = pointID; + fEventID = eventID; + fFileID = fileID; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +std::string MCPoint::ToString(int verbose, bool printHeader) const +{ + if (verbose < 1) { return std::string(); } + + std::stringstream msg; + if (printHeader) { + if (verbose > 0) { + msg << std::setw(10) << std::setfill(' ') << "track ID" << ' '; + msg << std::setw(10) << std::setfill(' ') << "mother ID" << '|'; + msg << std::setw(10) << std::setfill(' ') << "station ID" << '|'; + msg << std::setw(10) << std::setfill(' ') << "PDG" << ' '; + msg << std::setw(10) << std::setfill(' ') << "m [GeV/c2]" << ' '; + msg << std::setw(10) << std::setfill(' ') << "q [e]" << '|'; + msg << std::setw(14) << std::setfill(' ') << "t [ns]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "x [cm]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "y [cm]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "z [cm]" << '|'; + if (verbose > 1) { + msg << std::setw(14) << std::setfill(' ') << "zIn [cm]" << ' '; + msg << std::setw(14) << std::setfill(' ') << "zOut [cm]" << '|'; + msg << std::setw(14) << std::setfill(' ') << "p [GeV/c]" << '|'; + msg << std::setw(10) << std::setfill(' ') << "point ID" << ' '; + msg << std::setw(10) << std::setfill(' ') << "event ID" << ' '; + msg << std::setw(10) << std::setfill(' ') << "file ID" << ' '; + } + } + } + else { + if (verbose > 0) { + msg << std::setw(10) << std::setfill(' ') << fTrackID << ' '; + msg << std::setw(10) << std::setfill(' ') << fMotherID << '|'; + msg << std::setw(10) << std::setfill(' ') << fStationID << '|'; + msg << std::setw(10) << std::setfill(' ') << fPdgCode << ' '; + msg << std::setw(10) << std::setfill(' ') << fMass << ' '; + msg << std::setw(10) << std::setfill(' ') << fCharge << '|'; + msg << std::setw(14) << std::setfill(' ') << fTime << ' '; + msg << std::setw(14) << std::setfill(' ') << fPos[0] << ' '; + msg << std::setw(14) << std::setfill(' ') << fPos[1] << ' '; + msg << std::setw(14) << std::setfill(' ') << fPos[2] << '|'; + if (verbose > 1) { + msg << std::setw(14) << std::setfill(' ') << fPosIn[2] << ' '; + msg << std::setw(14) << std::setfill(' ') << fPosOut[2] << '|'; + msg << std::setw(14) << std::setfill(' ') << this->GetP() << '|'; + msg << std::setw(10) << std::setfill(' ') << fPointID << ' '; + msg << std::setw(10) << std::setfill(' ') << fEventID << ' '; + msg << std::setw(10) << std::setfill(' ') << fFileID << ' '; + } + } + } + return msg.str(); +} diff --git a/reco/L1/catools/CaToolsMCPoint.h b/reco/L1/catools/CaToolsMCPoint.h new file mode 100644 index 0000000000000000000000000000000000000000..dedc3ddff25a9a09dc851e5dbfa32223029741a0 --- /dev/null +++ b/reco/L1/catools/CaToolsMCPoint.h @@ -0,0 +1,287 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CaToolsMCPoint.h +/// \brief Internal class describing a MC point for CA tracking QA and performance (header) +/// \date 17.11.2022 +/// \author S.Zharko <s.zharko@gsi.de> + +#ifndef CaToolsMCPoint_h +#define CaToolsMCPoint_h 1 + +#include <string> + +#include "L1Undef.h" +#include "L1Vector.h" + +enum class L1DetectorID; + +/// Class describes a Monte-Carlo point used in CA tracking QA analysis +namespace ca::tools +{ + class MCTrack; + + class MCPoint { + public: + /// Default constructor + MCPoint() = default; + + /// Destructor + ~MCPoint() = default; + + /// Copy constructor + MCPoint(const MCPoint&) = default; + + /// Move constructor + MCPoint(MCPoint&&) = default; + + /// Copy assignment operator + MCPoint& operator=(const MCPoint&) = default; + + /// Move assignment operator + MCPoint& operator=(MCPoint&&) = default; + + + /// Adds hit ID + /// \param hitID A hit index in the external hits array + void AddHitID(int hitID) { fHitIDs.push_back_no_warning(hitID); } + + + // ********************* + // ** Getters ** + // ********************* + + /// Gets charge of the particle [e] + double GetCharge() const { return fCharge; } + + /// Gets detector ID + L1DetectorID GetDetectorID() const { return fDetectorID; } + + /// Gets event ID + int GetEventID() const { return fEventID; } + + /// Gets file ID + int GetFileID() const { return fFileID; } + + /// Gets mass of the particle [GeV/c2] + double GetMass() const { return fMass; } + + /// Gets mother ID of the track + int GetMotherID() const { return fMotherID; } + + /// Gets track momentum absolute value at reference z of station [GeV/c] + double GetP() const { return std::sqrt(fMom[0] * fMom[0] + fMom[1] * fMom[1] + fMom[2] * fMom[2]); } + + /// Gets PDG code of the particle + int GetPdgCode() const { return fPdgCode; } + + /// Gets track momentum absolute value at entrance to station [GeV/c] + double GetPIn() const { return std::sqrt(fMomIn[0] * fMomIn[0] + fMomIn[1] * fMomIn[1] + fMomIn[2] * fMomIn[2]); } + + /// Gets point ID + int GetPointID() const { return fPointID; } + + /// Gets track momentum absolute value at exit of station [GeV/c] + double GetPOut() const + { + return std::sqrt(fMomOut[0] * fMomOut[0] + fMomOut[1] * fMomOut[1] + fMomOut[2] * fMomOut[2]); + } + + /// Gets track momentum x component at reference z of station [GeV/c] + double GetPx() const { return fMom[0]; } + + /// Gets track momentum x component at entrance to station [GeV/c] + double GetPxIn() const { return fMomIn[0]; } + + /// Gets track momentum x component at exit of station [GeV/c] + double GetPxOut() const { return fMomOut[0]; } + + /// Gets track momentum y component at reference z of station [GeV/c] + double GetPy() const { return fMom[1]; } + + /// Gets track momentum y component at entrance to station [GeV/c] + double GetPyIn() const { return fMomIn[1]; } + + /// Gets track momentum y component at exit of station [GeV/c] + double GetPyOut() const { return fMomOut[1]; } + + /// Gets track momentum z component at reference z of station [GeV/c] + double GetPz() const { return fMom[2]; } + + /// Gets track momentum z component at entrance to station [GeV/c] + double GetPzIn() const { return fMomIn[2]; } + + /// Gets track momentum z component at exit of station [GeV/c] + double GetPzOut() const { return fMomOut[2]; } + + /// Gets global ID of the active tracking station + int GetStationID() const { return fStationID; } + + /// Gets time [ns] + double GetTime() const { return fTime; } + + /// Gets ID of track + int GetTrackID() const { return fTrackID; } + + /// Gets x coordinate at reference z of station [cm] + double GetX() const { return fPos[0]; } + + /// Gets x coordinate at entrance to station [cm] + double GetXIn() const { return fPosIn[0]; } + + /// Gets x coordinate at exit of station [cm] + double GetXOut() const { return fPosOut[0]; } + + /// Gets y coordinate at reference z of station [cm] + double GetY() const { return fPos[1]; } + + /// Gets y coordinate at entrance to station [cm] + double GetYIn() const { return fPosIn[1]; } + + /// Gets y coordinate at exit of station [cm] + double GetYOut() const { return fPosOut[1]; } + + /// Gets z coordinate at reference z of station [cm] + double GetZ() const { return fPos[2]; } + + /// Gets z coordinate at entrance to station [cm] + double GetZIn() const { return fPosIn[2]; } + + /// Gets z coordinate at exit of station [cm] + double GetZOut() const { return fPosOut[2]; } + + + /// Gets Hit ID + int MapHitID(int iHitID) { return fHitIDs[iHitID]; } + + + // ********************* + // ** Setters ** + // ********************* + + /// Sets particle charge [e] + void SetCharge(double charge) { fCharge = charge; } + + /// Sets detector ID + void SetDetectorID(L1DetectorID detID) { fDetectorID = detID; } + + /// Sets particle mass [GeV/c2] + void SetMass(double mass) { fMass = mass; } + + /// Sets mother ID + void SetMotherID(int motherID) { fMotherID = motherID; } + + /// Sets link address + /// \param pointID Index of + /// \param eventID + /// \param fileID + void SetLinkAddress(int pointID, int eventID, int fileID); + + /// Sets PDG code + void SetPdgCode(int pdg) { fPdgCode = pdg; } + + /// Sets track momentum x component at reference z of station [GeV/c] + void SetPx(double px) { fMom[0] = px; } + + /// Sets track momentum x component at entrance to station [GeV/c] + void SetPxIn(double px) { fMomIn[0] = px; } + + /// Sets track momentum x component at exit of station [GeV/c] + void SetPxOut(double px) { fMomOut[0] = px; } + + /// Sets track momentum y component at reference z of station [GeV/c] + void SetPy(double py) { fMom[1] = py; } + + /// Sets track momentum y component at entrance to station [GeV/c] + void SetPyIn(double py) { fMomIn[1] = py; } + + /// Sets track momentum y component at exit of station [GeV/c] + void SetPyOut(double py) { fMomOut[1] = py; } + + /// Sets track momentum z component at reference z of station [GeV/c] + void SetPz(double pz) { fMom[2] = pz; } + + /// Sets track momentum z component at entrance to station [GeV/c] + void SetPzIn(double pz) { fMomIn[2] = pz; } + + /// Sets track momentum z component at exit of station [GeV/c] + void SetPzOut(double pz) { fMomOut[2] = pz; } + + /// Sets global index of active station + void SetStationID(int stationID) { fStationID = stationID; } + + /// Sets time [ns] + void SetTime(double time) { fTime = time; } + + /// Sets track ID + void SetTrackID(int trackID) { fTrackID = trackID; } + + /// Sets x coordinate at reference z of station [cm] + void SetX(double x) { fPos[0] = x; } + + /// Sets x coordinate at entrance to station [cm] + void SetXIn(double x) { fPosIn[0] = x; } + + /// Sets x coordinate at exit of station [cm] + void SetXOut(double x) { fPosOut[0] = x; } + + /// Sets y coordinate at reference z of station [cm] + void SetY(double y) { fPos[1] = y; } + + /// Sets y coordinate at entrance to station [cm] + void SetYIn(double y) { fPosIn[1] = y; } + + /// Sets x coordinate at exit of station [cm] + void SetYOut(double y) { fPosOut[1] = y; } + + /// Sets z coordinate at reference z of station [cm] + void SetZ(double z) { fPos[2] = z; } + + /// Sets z coordinate at entrance to station [cm] + void SetZIn(double z) { fPosIn[2] = z; } + + /// Sets z coordinate at exit of station [cm] + void SetZOut(double z) { fPosOut[2] = z; } + + /// Prints content for a given verbosity level + /// \param verbose Verbosity level: + /// -#0: Prints nothing + /// -#1: Prints track ID, station ID, time and position + /// -#2: Also prints zIn, zOut, absolute momentum, as well as point, event and file IDs + /// \param printHeader If true, parameter names will be printed instead of the parameters themselves + std::string ToString(int verbose, bool printHeader = false) const; + + // **************************** + // ** Data variables ** + // **************************** + + private: + std::array<double, 3> fPos = {undef::kD, undef::kD, undef::kD}; ///< Position at reference z of station [cm] + std::array<double, 3> fPosIn = {undef::kD, undef::kD, undef::kD}; ///< Position at entrance to station [cm] + std::array<double, 3> fPosOut = {undef::kD, undef::kD, undef::kD}; ///< Position at exit of station [cm] + std::array<double, 3> fMom = {undef::kD, undef::kD, undef::kD}; ///< Momentum at reference z of station [cm] + std::array<double, 3> fMomIn = {undef::kD, undef::kD, undef::kD}; ///< Momentum at entrance to station [cm] + std::array<double, 3> fMomOut = {undef::kD, undef::kD, undef::kD}; ///< Momentum at exit of station [cm] + + int fPdgCode = undef::kI32; ///< Particle PDG code + double fMass = undef::kD; ///< Particle mass [GeV/c2] + double fCharge = undef::kD; ///< Particle charge [e] + double fTime = undef::kD; ///< Point time [ns] + + int fTrackID = undef::kI32; ///< Track ID + int fMotherID = undef::kI32; ///< Mother particle ID (NOTE: probably is redundant) + int fStationID = undef::kI32; ///< Global index of active tracking station + + L1DetectorID fDetectorID; ///< Detector ID of MC point + + int fPointID = undef::kI32; ///< TODO + int fEventID = undef::kI32; + int fFileID = undef::kI32; + + L1Vector<int> fHitIDs {"ca::tools::MCPoint::fHitIDs"}; + }; +} // namespace ca::tools + + +#endif // CaToolsMCPoint_h diff --git a/reco/L1/catools/CaToolsMcData.cxx b/reco/L1/catools/CaToolsMcData.cxx deleted file mode 100644 index 9b534c2ba0a8ff129ae7102d6697df74d179f1a3..0000000000000000000000000000000000000000 --- a/reco/L1/catools/CaToolsMcData.cxx +++ /dev/null @@ -1,49 +0,0 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/// \file CaToolsMcData.cxx -/// \brief Data structure for internal tracking MC-information (implementation) -/// \since 23.09.2022 -/// \author S.Zharko <s.zharko@gsi.de> - -#include "CaToolsMcData.h" - -#include <utility> // for std::move - -using namespace ca::tools; - -// --------------------------------------------------------------------------------------------------------------------- -// -McData::McData() {} - -// --------------------------------------------------------------------------------------------------------------------- -// -McData::McData(const McData& /*other*/) {} - -// --------------------------------------------------------------------------------------------------------------------- -// -McData::McData(McData&& other) noexcept { this->Swap(other); } - -// --------------------------------------------------------------------------------------------------------------------- -// -McData& McData::operator=(const McData& other) -{ - if (this != &other) { McData(other).Swap(*this); } - return *this; -} - -// --------------------------------------------------------------------------------------------------------------------- -// -McData& McData::operator=(McData&& other) noexcept -{ - if (this != &other) { - McData tmp(std::move(other)); - this->Swap(tmp); - } - return *this; -} - -// --------------------------------------------------------------------------------------------------------------------- -// -void McData::Swap(McData& /*other*/) noexcept {} diff --git a/reco/L1/catools/CaToolsMcData.h b/reco/L1/catools/CaToolsMcData.h deleted file mode 100644 index c677e72a773d632b5c5301c86b75835d3ad57242..0000000000000000000000000000000000000000 --- a/reco/L1/catools/CaToolsMcData.h +++ /dev/null @@ -1,57 +0,0 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/// \file CaToolsMcData.h -/// \brief Data structure for internal tracking MC-information (header) -/// \since 23.09.2022 -/// \author S.Zharko <s.zharko@gsi.de> - -#ifndef CaToolsMcData_h -#define CaToolsMcData_h 1 - -namespace ca -{ - namespace tools - { - class McData { - public: - // ********************************* - // ** Constructors and destructor ** - // ********************************* - - /// Default constructor - McData(); - - /// Destructor - ~McData() = default; - - /// Copy constructor - McData(const McData& other); - - /// Move constructor - McData(McData&& other) noexcept; - - /// Copy assignment operator - McData& operator=(const McData& other); - - /// Move assignment operator - McData& operator=(McData&& other) noexcept; - - /// Swap method - void Swap(McData& other) noexcept; - - - // ********************** - // ** Access functions ** - // ********************** - - private: - // ********************** - // ** Member variables ** - // ********************** - }; - } // namespace tools -} // namespace ca - -#endif // CaToolsMcData_h diff --git a/reco/L1/catools/CaToolsMcDataManager.cxx b/reco/L1/catools/CaToolsMcDataManager.cxx deleted file mode 100644 index 3fae609bb4a6289077b67c7f735a9adce764f7af..0000000000000000000000000000000000000000 --- a/reco/L1/catools/CaToolsMcDataManager.cxx +++ /dev/null @@ -1,26 +0,0 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/// \file CaToolsMcDataManager.cxx -/// \brief Manager class for handling McData structure (implementation) -/// \since 23.09.2022 -/// \author S.Zharko <s.zharko@gsi.de> - -#include "CaToolsMcDataManager.h" - -#include <cassert> -#include <utility> - -#include "CaToolsPerformance.h" - -using namespace ca::tools; - -// --------------------------------------------------------------------------------------------------------------------- -// -bool McDataManager::SendMcData(Performance* pPerformance) -{ - assert(pPerformance); - pPerformance->ReceiveMcData(std::move(fMcData)); - return true; -} diff --git a/reco/L1/catools/CaToolsMcDataManager.h b/reco/L1/catools/CaToolsMcDataManager.h deleted file mode 100644 index 0e6cb9ebf604c9ee4272b65788f625c97e42b2f5..0000000000000000000000000000000000000000 --- a/reco/L1/catools/CaToolsMcDataManager.h +++ /dev/null @@ -1,58 +0,0 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/// \file CaToolsMcDataManager.h -/// \brief Manager class for handling McData structure (header) -/// \since 23.09.2022 -/// \author S.Zharko <s.zharko@gsi.de> - -#ifndef CaToolsMcDataManager_h -#define CaToolsMcDataManager_h 1 - -#include "CaToolsMcData.h" - -namespace ca -{ - namespace tools - { - class Performance; - class McDataManager { - public: - // ********************************* - // ** Constructors and destructor ** - // ********************************* - - /// Default constructor - McDataManager() {} - - /// Destructor - ~McDataManager() = default; - - /// Copy constructor - McDataManager(const McDataManager&) = delete; - - /// Move constructor - McDataManager(McDataManager&&) = delete; - - /// Copy assignment operator - McDataManager& operator=(const McDataManager&) = delete; - - /// Move assignment operator - McDataManager& operator=(McDataManager&&) = delete; - - /// Sends (moves) MC data object to the destination Performance instance - /// \param pPerformance Pointer to the destination Performance instance - /// \return Success flag - bool SendMcData(Performance* pPerformance); - - private: - // ********************** - // ** Member variables ** - // ********************** - McData fMcData {}; ///< Object with internal MC data - }; - } // namespace tools -} // namespace ca - -#endif // CaToolsMcDataManager_h diff --git a/reco/L1/catools/CaToolsPerformance.cxx b/reco/L1/catools/CaToolsPerformance.cxx index 37fc72ebee3144345880339736901868b8d3690a..f814df8d3ec5ee83be98842d235e711d677ee146 100644 --- a/reco/L1/catools/CaToolsPerformance.cxx +++ b/reco/L1/catools/CaToolsPerformance.cxx @@ -15,7 +15,7 @@ using namespace ca::tools; // --------------------------------------------------------------------------------------------------------------------- // -void Performance::ReceiveMcData(McData&& mcData) { fMcData = std::move(mcData); } +void Performance::ReceiveMCData(MCData&& mcData) { fMCData = std::move(mcData); } // --------------------------------------------------------------------------------------------------------------------- // diff --git a/reco/L1/catools/CaToolsPerformance.h b/reco/L1/catools/CaToolsPerformance.h index cb5c4730ca9ccb2384355430a6349edca7d4c518..4f2cb15858cfbd4034459d1ea92b29b1807c483f 100644 --- a/reco/L1/catools/CaToolsPerformance.h +++ b/reco/L1/catools/CaToolsPerformance.h @@ -10,7 +10,7 @@ #ifndef CaToolsPerformance_h #define CaToolsPerformance_h 1 -#include "CaToolsMcData.h" +#include "CaToolsMCData.h" class L1Parameters; @@ -50,13 +50,13 @@ namespace ca // ** Accessors ** // *********************** - /// Receives MC data object from McDataManager - void ReceiveMcData(McData&& mcData); + /// Receives MC data object from MCDataManager + void ReceiveMCData(MCData&& mcData); private: // const L1Parameters* fpParameters = nullptr; ///< Instance of the tracking core class parameters - McData fMcData = {}; ///< Object of MC data + MCData fMCData = {}; ///< Object of MC data }; } // namespace tools } // namespace ca