From dd43ca5fc2bc46ac05c31a302fb135f91769702d Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Wed, 20 Sep 2023 16:48:30 +0000
Subject: [PATCH] CA: qa: account for several MC points per hit
---
reco/L1/CbmCaMCModule.cxx | 18 +++--
reco/L1/CbmCaMCModule.h | 42 ++++++++----
reco/L1/CbmCaTimeSliceReader.cxx | 1 +
reco/L1/CbmL1.cxx | 5 +-
reco/L1/CbmL1.h | 30 ++++++---
reco/L1/CbmL1Hit.h | 38 ++++++-----
reco/L1/CbmL1Performance.cxx | 66 ++++++++-----------
reco/L1/CbmL1ReadEvent.cxx | 54 ++++++++-------
reco/L1/L1Algo/L1Algo.cxx | 4 +-
reco/L1/L1Algo/utils/L1AlgoDraw.cxx | 4 +-
.../utils/L1AlgoEfficiencyPerformance.h | 8 +--
reco/L1/L1Algo/utils/L1AlgoPulls.cxx | 2 +-
reco/L1/qa/CbmCaTrackTypeQa.cxx | 4 +-
13 files changed, 151 insertions(+), 125 deletions(-)
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index 44c3642310..b296ec7821 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -19,6 +19,7 @@
#include "CbmMCTrack.h"
#include "CbmStsHit.h"
#include "CbmTimeSlice.h"
+
#include "FairEventHeader.h"
#include "FairMCEventHeader.h"
#include "FairRootManager.h"
@@ -240,17 +241,15 @@ void MCModule::MatchRecoAndMCTracks()
// ----- Count number of hits from each MC track belonging to this reconstructed track
auto& mNofHitsVsMCTrkID = trkRe.hitMap;
- //mNofHitsVsMCTrkID.clear();
+ mNofHitsVsMCTrkID.clear();
for (int iH : trkRe.Hits) {
- int iP = (*fpvQaHits)[iH].GetMCPointIndex();
- if (iP > -1) {
+ auto& vP = (*fpvQaHits)[iH].GetMcPointIds();
+ for (int iP : vP) {
+ if (iP < 0) { continue; }
int iTmc = fpMCData->GetPoint(iP).GetTrackId();
- if (mNofHitsVsMCTrkID.find(iTmc) == mNofHitsVsMCTrkID.end()) { mNofHitsVsMCTrkID[iTmc] = 1; }
- else {
- mNofHitsVsMCTrkID[iTmc] += 1;
- }
+ mNofHitsVsMCTrkID[iTmc]++;
}
- } // loop over hit ids stored for a reconstructed track: end
+ } // loop over hit ids stored for a reconstructed track: end
// ----- Calculate track max purity
@@ -282,8 +281,7 @@ void MCModule::MatchRecoAndMCTracks()
// Update max purity of the reconstructed track
trkRe.SetMaxPurity(maxPurity);
} // loop over hit map: end
- assert(trkRe.GetNofMCTracks() < 2);
- } // loop over reconstructed tracks: end
+ } // loop over reconstructed tracks: end
}
// *******************************
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
index c83315b159..d817e97e15 100644
--- a/reco/L1/CbmCaMCModule.h
+++ b/reco/L1/CbmCaMCModule.h
@@ -109,7 +109,7 @@ namespace cbm::ca
/// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
/// link in the match and returns the corresponding global ID of the MC points.
template<L1DetectorID DetId>
- int MatchHitWithMc(int iHitExt);
+ std::tuple<int, vector<int>> MatchHitWithMc(int iHitExt);
/// @brief Match reconstructed and MC data
///
@@ -229,10 +229,10 @@ namespace cbm::ca
::ca::Monitor<EMonitorKey> fMonitor {"CA MC Module"}; ///< Monitor
// ------ Flags
- DetIdArr_t<bool> fvbUseDet = {{false}}; ///< Flag: is detector subsystem used
- bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
- int fVerbose = 0; ///< Verbosity level
- int fPerformanceMode = -1; ///< Mode of performance
+ DetIdArr_t<bool> fvbUseDet = {{false}}; ///< Flag: is detector subsystem used
+ bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
+ int fVerbose = 0; ///< Verbosity level
+ int fPerformanceMode = -1; ///< Mode of performance
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
// ------ Input data branches
@@ -273,9 +273,10 @@ namespace cbm::ca
// ---------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
-int cbm::ca::MCModule::MatchHitWithMc(int iHitExt)
+std::tuple<int, vector<int>> cbm::ca::MCModule::MatchHitWithMc(int iHitExt)
{
- int iPoint = -1;
+ int iPoint = -1;
+ std::vector<int> vPoints;
const auto* pHitMatch = dynamic_cast<CbmMatch*>(fvpBrHitMatches[DetID]->At(iHitExt));
if (!pHitMatch) {
LOG(warn) << "Hit match with index " << iHitExt << " is missing for " << kDetName[DetID];
@@ -292,7 +293,17 @@ int cbm::ca::MCModule::MatchHitWithMc(int iHitExt)
else if constexpr (L1DetectorID::kTof == DetID) {
fMonitor.Increment(EMonitorKey::kMissedMatchesTof);
}
- return iPoint;
+ return std::tuple(iPoint, vPoints);
+ }
+
+ for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) {
+ const auto& link = pHitMatch->GetLink(iLink);
+ int iPointExt = link.GetIndex();
+ int iEvent = link.GetEntry();
+ int iFile = link.GetFile();
+ if (iPointExt < 0) continue;
+ int id = fpMCData->FindInternalPointIndex(DetID, iPointExt, iEvent, iFile);
+ vPoints.push_back(id);
}
if (pHitMatch->GetNofLinks() > 0) {
@@ -305,7 +316,7 @@ int cbm::ca::MCModule::MatchHitWithMc(int iHitExt)
}
}
- return iPoint;
+ return std::tuple(iPoint, vPoints);
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -456,11 +467,14 @@ void cbm::ca::MCModule::MatchPointsAndHits()
if (!fvbUseDet[DetID]) { return; }
for (int iH = (*fpvFstHitId)[static_cast<int>(DetID)]; iH < (*fpvFstHitId)[static_cast<int>(DetID) + 1]; ++iH) {
- auto& hit = (*fpvQaHits)[iH];
- int iP = MatchHitWithMc<DetID>(hit.ExtIndex);
- if (iP > -1) {
- hit.SetMCPointIndex(iP);
- fpMCData->GetPoint(iP).AddHitID(iH);
+ auto& hit = (*fpvQaHits)[iH];
+ auto [iBestP, vAllP] = MatchHitWithMc<DetID>(hit.ExtIndex);
+ if (iBestP >= 0) { hit.SetBestMcPointId(iBestP); }
+ for (auto iP : vAllP) {
+ if (iP >= 0) {
+ hit.AddMcPointId(iP);
+ fpMCData->GetPoint(iP).AddHitID(iH);
+ }
}
}
}
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index 94d7ff9af8..ffa58bc595 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -427,6 +427,7 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
aHitQa.dy = sqrt(hitRecord.fDy2);
aHitQa.dxy = hitRecord.fDxy;
aHitQa.time = hitRecord.fT;
+ aHitQa.dt = sqrt(hitRecord.fDt2);
fpvQaHits->push_back(aHitQa);
}
}
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 3b57facba8..61a2f307a9 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -1027,10 +1027,7 @@ void CbmL1::IdealTrackFinder()
for (unsigned int iH = 0; iH < MC.Hits.size(); iH++) {
const int hitI = MC.Hits[iH];
const CbmL1HitDebugInfo& hit = fvHitDebugInfo[hitI];
- const int iP = hit.GetMCPointIndex();
- assert(iP >= 0);
- const int iStation = fvMCPoints[iP].iStation;
-
+ const int iStation = hit.GetStationId();
if (iStation >= 0) hitIndices[iStation] = hitI;
}
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index f2283fb1b3..8ed28ff156 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -292,7 +292,8 @@ public:
const L1Vector<CbmL1MCTrack>& GetMcTracks() const { return fvMCTracks; }
- const L1Vector<int>& GetHitMCRefs() const { return fvHitPointIndexes; }
+ const L1Vector<vector<int>>& GetHitMcRefs() const { return fvHitPointIndices; }
+ const L1Vector<int>& GetHitBestMcRefs() const { return fvHitBestPointIndices; }
static double boundedGaus(double sigma);
@@ -378,7 +379,7 @@ private:
/// \param iHit External index of hit
/// \return MC-point index in fvMCPoints array
template<L1DetectorID DetId>
- int MatchHitWithMc(int iHit) const;
+ std::tuple<int, std::vector<int>> MatchHitWithMc(int iHit) const;
/// Procedure for match hits and MCPoints.
/// Reads information about correspondence between hits and MC-points and fill CbmL1MCPoint::hitIds and CbmL1Hit::mcPointIds arrays
@@ -629,8 +630,9 @@ private:
TClonesArray* fpMuchDigiMatches = nullptr; ///< Array of MuCh digi matches (NOTE: currently unsused)
- L1Vector<CbmL1MCTrack> fvMCTracks = {"CbmL1::fvMCTracks"}; ///< Array of MC tracks
- L1Vector<int> fvHitPointIndexes = {"CbmL1::fvHitPointIndexes"}; ///< Indexes of MC points vs. hit index
+ L1Vector<CbmL1MCTrack> fvMCTracks = {"CbmL1::fvMCTracks"}; ///< Array of MC tracks
+ L1Vector<vector<int>> fvHitPointIndices = {"CbmL1::fvHitPointIndices"}; ///< Indices of MC points vs. hit index
+ L1Vector<int> fvHitBestPointIndices = {"CbmL1::fvHitBestPointIndices"}; ///< Indices of best MC points vs. hit index
int fEventNo = 0; ///< Current number of event/TS
int fNofRecoTracks = 0; ///< Total number of reconstructed tracks
@@ -690,7 +692,7 @@ private:
// ---------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
-int CbmL1::MatchHitWithMc(int iHitExt) const
+std::tuple<int, std::vector<int>> CbmL1::MatchHitWithMc(int iHitExt) const
{
const CbmMatch* pHitMatch = nullptr;
int indexShift = 0; // Number of points in detectors, standing in front
@@ -718,6 +720,18 @@ int CbmL1::MatchHitWithMc(int iHitExt) const
assert(pHitMatch);
int iPointInt = -1;
+ std::vector<int> vPoints;
+ for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) {
+ const auto& link = pHitMatch->GetLink(iLink);
+ if (link.GetIndex() > -1) {
+ int iPointExt = link.GetIndex();
+ int iEvent = link.GetEntry();
+ int iFile = link.GetFile();
+ auto itPoint = fmMCPointsLinksMap.find(CbmL1LinkKey(iPointExt + indexShift, iEvent, iFile));
+ if (itPoint == fmMCPointsLinksMap.cend()) { continue; }
+ vPoints.push_back(itPoint->second);
+ }
+ }
if (pHitMatch->GetNofLinks() > 0) {
const auto& link = pHitMatch->GetMatchedLink();
if (link.GetIndex() > -1) {
@@ -725,11 +739,11 @@ int CbmL1::MatchHitWithMc(int iHitExt) const
int iEvent = link.GetEntry();
int iFile = link.GetFile();
auto itPoint = fmMCPointsLinksMap.find(CbmL1LinkKey(iPointExt + indexShift, iEvent, iFile));
- if (itPoint == fmMCPointsLinksMap.cend()) { return -1; }
- iPointInt = itPoint->second;
+ if (itPoint != fmMCPointsLinksMap.cend()) { iPointInt = itPoint->second; }
}
}
- return iPointInt;
+
+ return std::tuple(iPointInt, vPoints);
}
diff --git a/reco/L1/CbmL1Hit.h b/reco/L1/CbmL1Hit.h
index d6090ad65b..8441879aa2 100644
--- a/reco/L1/CbmL1Hit.h
+++ b/reco/L1/CbmL1Hit.h
@@ -65,7 +65,9 @@ public:
int GetExternalHitId() const { return ExtIndex; }
/// @brief Gets index of matched MC point
- int GetMCPointIndex() const { return fPointID; }
+ int GetBestMcPointId() const { return fBestMcPointId; }
+
+ const std::vector<int>& GetMcPointIds() const { return fMcPointIds; }
/// @brief Gets global index of active tracking station
int GetStationId() const { return iStation; }
@@ -84,8 +86,11 @@ public:
/// @brief Sets index of matched MC point
/// @param pointID Internal index of MC point
- void SetMCPointIndex(int pointID) { fPointID = pointID; }
+ void SetBestMcPointId(int pointID) { fBestMcPointId = pointID; }
+ /// @brief Sets index of matched MC point
+ /// @param pointID Internal index of MC point
+ void AddMcPointId(int pointID) { fMcPointIds.push_back(pointID); }
/// @brief String representation of the object
/// @param verbose Verbosity level
@@ -118,7 +123,7 @@ public:
msg << setw(8) << setfill(' ') << IntIndex << ' ';
msg << setw(8) << setfill(' ') << iStation << ' ';
msg << setw(8) << setfill(' ') << Det << ' ';
- msg << setw(8) << setfill(' ') << GetMCPointIndex() << ' ';
+ msg << setw(8) << setfill(' ') << GetBestMcPointId() << ' ';
msg << setw(14) << setfill(' ') << x << ' ';
msg << setw(14) << setfill(' ') << y << ' ';
msg << setw(14) << setfill(' ') << z << ' ';
@@ -134,19 +139,20 @@ public:
}
// TODO: SZh 2.03.2023: make the variables private
- int ExtIndex; ///< index of hit in the external branch
- int IntIndex; ///< index of hit in the internal array
- int iStation; ///< index of station in active stations array
- int Det; ///< detector subsystem ID
- double x; ///< x coordinate of position [cm]
- double y; ///< y coordinate of position [cm]
- double z; ///< z coordinate of position [cm]
- double time; ///< hit time [ns]
- double dx; ///< x coordinate error [cm]
- double dy; ///< y coordinate error [cm]
- double dt; ///< time error [ns]
- double dxy; ///< covariance between x and y [cm2]
- int fPointID = -1; ///< index of matched MC point
+ int ExtIndex; ///< index of hit in the external branch
+ int IntIndex; ///< index of hit in the internal array
+ int iStation; ///< index of station in active stations array
+ int Det; ///< detector subsystem ID
+ double x; ///< x coordinate of position [cm]
+ double y; ///< y coordinate of position [cm]
+ double z; ///< z coordinate of position [cm]
+ double time; ///< hit time [ns]
+ double dx; ///< x coordinate error [cm]
+ double dy; ///< y coordinate error [cm]
+ double dt; ///< time error [ns]
+ double dxy; ///< covariance between x and y [cm2]
+ int fBestMcPointId {-1}; ///< index of best matched MC point
+ std::vector<int> fMcPointIds {}; ///< indices of all matched MC points
};
#endif
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index e490df8e6e..0e73bf2bdc 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -99,22 +99,16 @@ void CbmL1::TrackMatch()
map<int, int>& hitmap = prtra->hitMap; // how many hits from each mcTrack belong to current recoTrack
map<int, int> hitmapChain;
for (vector<int>::iterator ih = (prtra->Hits).begin(); ih != (prtra->Hits).end(); ++ih) {
- int iMC = fvHitDebugInfo[*ih].GetMCPointIndex();
-
- // cout<<iMC<<" iMC"<<endl;
- int ID = -1;
- int chainID = -1;
- if (iMC >= 0) {
- ID = fvMCPoints[iMC].ID;
- chainID = fvMCTracks[ID].chainID;
- }
- if (hitmap.find(ID) == hitmap.end()) hitmap[ID] = 1;
- else {
- hitmap[ID] += 1;
- }
- if (hitmapChain.find(chainID) == hitmapChain.end()) hitmapChain[chainID] = 1;
- else {
- hitmapChain[chainID] += 1;
+ auto& vP = fvHitDebugInfo[*ih].GetMcPointIds();
+ for (int iP : vP) {
+ int ID = -1;
+ int chainID = -1;
+ if (iP >= 0) {
+ ID = fvMCPoints[iP].ID;
+ chainID = fvMCTracks[ID].chainID;
+ }
+ hitmap[ID]++;
+ hitmapChain[chainID]++;
}
} // for iHit
@@ -123,16 +117,26 @@ void CbmL1::TrackMatch()
for (map<int, int>::iterator posIt = hitmap.begin(); posIt != hitmap.end();
posIt++) { // loop over all touched MCTracks
- if (posIt->first < 0) continue; // not a MC track - based on fake hits
+ int iMcTrack = posIt->first;
+ int nMcTrackHits = posIt->second;
+
+ if (iMcTrack < 0) continue; // not a MC track - based on fake hits
+
+ if (0) { // debug: consider tracks from decay chains as one MC track
+ nMcTrackHits = hitmapChain[fvMCTracks[iMcTrack].chainID];
+ }
+
+ double purity = double(nMcTrackHits) / double(hitsum);
// count max-purity
- if (double(posIt->second) > max_percent * double(hitsum)) max_percent = double(posIt->second) / double(hitsum);
+ if (max_percent < purity) { max_percent = purity; }
// set relation to the mcTrack
- if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) { continue; }
- CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
+ assert(pMCTrackMap.find(iMcTrack) != pMCTrackMap.end());
- if (double(posIt->second) >= CbmL1Constants::MinPurity * double(hitsum)) { // found correspondent MCTrack
+ CbmL1MCTrack* pmtra = pMCTrackMap[iMcTrack];
+
+ if (purity >= CbmL1Constants::MinPurity) { // found correspondent MCTrack
pmtra->AddRecoTrack(prtra);
pmtra->AddRecoTrackIndex(iR);
prtra->AddMCTrack(pmtra);
@@ -144,20 +148,6 @@ void CbmL1::TrackMatch()
} // for hitmap
prtra->SetMaxPurity(max_percent);
- if (0) { // debug
-
- double max_percentChain = 0.0; // [%]. maximum persent of hits, which belong to one mcTrack.
- for (map<int, int>::iterator posIt = hitmapChain.begin(); posIt != hitmapChain.end();
- posIt++) { // loop over all touched MCTracks
- if (posIt->first < 0) continue; // not a MC track - based on fake hits
- // count max-purity
- if (double(posIt->second) > max_percentChain * double(hitsum)) {
- max_percentChain = double(posIt->second) / double(hitsum);
- }
- } // for hitmap
- prtra->SetMaxPurity(max_percentChain);
- }
-
} // for reco tracks
} //
@@ -1120,7 +1110,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
int NFakes = 0;
for (unsigned int ih = 0; ih < fpAlgo->GetInputData().GetNhits(); ih++) {
- int iMC = fvHitPointIndexes[ih]; // TODO2: adapt to linking
+ int iMC = fvHitBestPointIndices[ih]; // TODO2: adapt to linking
if (iMC < 0) NFakes++;
}
@@ -1334,7 +1324,7 @@ void CbmL1::TrackFitPerformance()
const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]);
- int imcPoint = fvHitPointIndexes[it->Hits.front()];
+ int imcPoint = fvHitBestPointIndices[it->Hits.front()];
// extrapolate to the first mc point of the mc track !
@@ -1385,7 +1375,7 @@ void CbmL1::TrackFitPerformance()
{ // last hit
- int iMC = fvHitPointIndexes[it->Hits.back()]; // TODO2: adapt to linking
+ int iMC = fvHitBestPointIndices[it->Hits.back()]; // TODO2: adapt to linking
if (iMC < 0) continue;
CbmL1MCPoint& mcP = fvMCPoints[iMC];
L1TrackPar tr;
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index f92fefcdc5..36f6276a45 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -77,12 +77,13 @@ struct TmpHit {
double dx2; ///< hit position uncertainty along x axis [cm]
double dy2; ///< hit position uncertainty along y axis [cm]
double dxy; ///< hit position covariance along x and y axes [cm2]
- int iMC; ///< index of MCPoint in the fvMCPoints array
- double time = 0.; ///< time of the hit [ns]
- double dt2 = 1.e8; ///< time error of the hit [ns]
- double rangeX; ///< hit range in X [cm]
- double rangeY; ///< hit range in Y [cm]
- double rangeT; ///< hit range in T [ns]
+ std::vector<int> vMc {};
+ int iBestMc {-1}; ///< index of MCPoint in the fvMCPoints array
+ double time = 0.; ///< time of the hit [ns]
+ double dt2 = 1.e8; ///< time error of the hit [ns]
+ double rangeX; ///< hit range in X [cm]
+ double rangeY; ///< hit range in Y [cm]
+ double rangeT; ///< hit range in T [ns]
int Det;
@@ -139,11 +140,12 @@ void CbmL1::ReadEvent(CbmEvent* event)
if (fVerbose >= 10) cout << "ReadEvent: start." << endl;
// clear arrays for next event
- fvMCPoints.clear(); /* <CbmL1MCPoint> */
- fvMCTracks.clear(); /* <CbmL1MCTrack> */
- fvExternalHits.clear(); /* <CbmL1Hit> */
- fvHitPointIndexes.clear(); /* <int>: indexes of MC-points in fvMCPoints (by index of fpAlgo->vHits) */
- fvHitDebugInfo.clear(); /* <CbmL1HitDebugInfo> */
+ fvMCPoints.clear(); /* <CbmL1MCPoint> */
+ fvMCTracks.clear(); /* <CbmL1MCTrack> */
+ fvExternalHits.clear(); /* <CbmL1Hit> */
+ fvHitPointIndices.clear(); /* <int>: indices of MC-points in fvMCPoints (by index of fpAlgo->vHits) */
+ fvHitBestPointIndices.clear(); /* <int>: indices of MC-points in fvMCPoints (by index of fpAlgo->vHits) */
+ fvHitDebugInfo.clear(); /* <CbmL1HitDebugInfo> */
fmMCPointsLinksMap.clear();
fmMCTracksLinksMap.clear();
@@ -473,8 +475,8 @@ void CbmL1::ReadEvent(CbmEvent* event)
std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmMvdTrackingInterface::Instance()->GetHitRanges(*h);
}
- th.Det = 0;
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMvd>(hitIndex) : -1;
+ th.Det = 0;
+ if (fPerformance) { std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<L1DetectorID::kMvd>(hitIndex); }
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
//if( th.iMC >=0 ) // DEBUG !!!!
if (th.IfAccept()) {
@@ -542,7 +544,8 @@ void CbmL1::ReadEvent(CbmEvent* event)
std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmStsTrackingInterface::Instance()->GetHitRanges(*h);
}
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kSts>(hitIndex) : -1;
+ if (fPerformance) { std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<L1DetectorID::kSts>(hitIndex); }
+
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
if (th.IfAccept()) {
tmpHits.push_back(th);
@@ -597,7 +600,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmMuchTrackingInterface::Instance()->GetHitRanges(*h);
}
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMuch>(th.ExtIndex) : -1;
+ if (fPerformance) { std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<L1DetectorID::kMuch>(hitIndex); }
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
if (th.IfAccept()) {
tmpHits.push_back(th);
@@ -669,7 +672,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmTrdTrackingInterface::Instance()->GetHitRanges(*h);
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTrd>(hitIndex) : -1;
+ if (fPerformance) { std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<L1DetectorID::kTrd>(hitIndex); }
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
if (th.IfAccept()) {
@@ -737,7 +740,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
//TODO: is it still needed here? (S.Zharko)
if (L1Algo::TrackingMode::kMcbm == fTrackingMode && th.z > 400) continue;
- th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTof>(hitIndex) : -1;
+ if (fPerformance) { std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<L1DetectorID::kTof>(hitIndex); }
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
if (th.IfAccept()) {
@@ -766,11 +769,12 @@ void CbmL1::ReadEvent(CbmEvent* event)
if (fVerbose >= 10) { cout << "ReadEvent: strips are read." << endl; }
- // ----- Fill and save fvExternalHits, fvHitDebugInfo and fvHitPointIndexes vectors as well as fpData->vHits
+ // ----- Fill and save fvExternalHits, fvHitDebugInfo and fvHitPointIndices vectors as well as fpData->vHits
fvExternalHits.reserve(nHits);
fvHitDebugInfo.reserve(nHits);
- fvHitPointIndexes.reserve(nHits);
+ fvHitPointIndices.reserve(nHits);
+ fvHitBestPointIndices.reserve(nHits);
fpIODataManager->ResetInputData(nHits);
fpIODataManager->SetNhitKeys(NStrips);
@@ -821,7 +825,8 @@ void CbmL1::ReadEvent(CbmEvent* event)
fpIODataManager->PushBackHit(h, th.fDataStream);
fvHitDebugInfo.push_back(s);
- fvHitPointIndexes.push_back(th.iMC);
+ fvHitPointIndices.push_back(th.vMc);
+ fvHitBestPointIndices.push_back(th.iBestMc);
}
if (fPerformance) { HitMatch(); } /// OLD
@@ -1092,16 +1097,17 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i
void CbmL1::HitMatch()
{
// Clear contents
- std::for_each(fvHitDebugInfo.begin(), fvHitDebugInfo.end(), [&](auto& hit) { hit.SetMCPointIndex(-1); });
std::for_each(fvMCPoints.begin(), fvMCPoints.end(), [&](auto& point) { point.hitIds.clear(); });
// Fill new contents
const int NHits = fvHitDebugInfo.size();
for (int iH = 0; iH < NHits; iH++) {
CbmL1HitDebugInfo& hit = fvHitDebugInfo[iH];
- int iP = fvHitPointIndexes[iH];
- if (iP >= 0) {
- hit.SetMCPointIndex(iP);
+ hit.SetBestMcPointId(fvHitBestPointIndices[iH]);
+ hit.fMcPointIds.clear();
+ for (int iP : fvHitPointIndices[iH]) {
+ if (iP < 0) { continue; }
+ hit.AddMcPointId(iP);
fvMCPoints[iP].hitIds.push_back_no_warning(iH);
}
}
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 7e4e734b46..c1beb70966 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -103,7 +103,7 @@ std::pair<fscal, fscal> L1Algo::GetHitCoorOnGrid(const L1Hit& h)
int L1Algo::GetMcTrackIdForCaHit(int iHit) const
{
int hitId = iHit;
- int iMcPoint = CbmL1::Instance()->GetHitMCRefs()[hitId];
+ int iMcPoint = CbmL1::Instance()->GetHitBestMcRefs()[hitId];
if (iMcPoint < 0) return -1;
return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
}
@@ -111,7 +111,7 @@ int L1Algo::GetMcTrackIdForCaHit(int iHit) const
int L1Algo::GetMcTrackIdForGridHit(int iHit) const
{
int hitId = fGridHitIds[iHit];
- int iMcPoint = CbmL1::Instance()->GetHitMCRefs()[hitId];
+ int iMcPoint = CbmL1::Instance()->GetHitBestMcRefs()[hitId];
if (iMcPoint < 0) return -1;
return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
}
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index 400df85bac..6b7b4134d1 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -558,7 +558,7 @@ void L1AlgoDraw::DrawInputHits()
Int_t n_poly_fake = 0;
for (int ih = HitsStartIndex[ista]; ih < HitsStopIndex[ista]; ih++) {
L1Hit& h = vHits[ih];
- int iMC = CbmL1::Instance()->fvHitPointIndexes[ih];
+ int iMC = CbmL1::Instance()->GetHitBestMcRefs()[ih];
//if( (vSFlag[h.f] | vSFlagB[h.b] )&0x02 ) continue; // if used
fscal x = h.x, y = h.y, z = h.z;
@@ -691,7 +691,7 @@ void L1AlgoDraw::DrawRestHits(L1HitIndex_t* StsRestHitsStartIndex, L1HitIndex_t*
for (L1HitIndex_t iRestHit = StsRestHitsStartIndex[ista]; iRestHit < StsRestHitsStopIndex[ista]; iRestHit++) {
int ih = realIHit[iRestHit];
L1Hit& h = vHits[ih];
- int iMC = CbmL1::Instance()->fvHitPointIndexes[ih];
+ int iMC = CbmL1::Instance()->GetHitBestMcRefs()[ih];
//if( (vSFlag[h.f] | vSFlagB[h.b] )&0x02 ) continue; // if used
fscal x = h.x, y = h.y, z = h.z;
diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
index 6d767fe104..5dc83b177a 100644
--- a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
+++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
@@ -219,12 +219,12 @@ bool L1AlgoEfficiencyPerformance<NHits>::AddOne(L1HitIndex_t* iHits)
{
L1RecoTracklet trlet(iHits);
- // --- check is all hits belong to one mcTrack ---
+ // --- check if all hits belong to one mcTrack ---
vector<int> mcIds[NHits];
// obtain mc data
for (int iih = 0; iih < NHits; iih++) {
- const int iMCP = fL1->fvHitDebugInfo[iHits[iih]].GetMCPointIndex();
+ const int iMCP = fL1->fvHitDebugInfo[iHits[iih]].GetBestMcPointId();
if (iMCP > -1) {
int mcId = fL1->fvMCPoints[iMCP].ID;
mcIds[iih].push_back(mcId);
@@ -254,9 +254,9 @@ bool L1AlgoEfficiencyPerformance<NHits>::AddOne(L1HitIndex_t* iHits)
for (unsigned int i = 0; i < mcsN.size(); i++) {
if (mcsN[i] >= 0) {
trlet.mcTrackId = mcsN[i];
- int iMCP = fL1->fvHitDebugInfo[iHits[0]].GetMCPointIndex();
+ int iMCP = fL1->fvHitDebugInfo[iHits[0]].GetBestMcPointId();
assert(iMCP > -1);
- trlet.iStation = fL1->fvMCPoints[fL1->fvHitDebugInfo[iHits[0]].GetMCPointIndex()].iStation;
+ trlet.iStation = fL1->fvMCPoints[fL1->fvHitDebugInfo[iHits[0]].GetBestMcPointId()].iStation;
break;
}
}
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
index 7a84778507..c5ac52ef40 100644
--- a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
@@ -127,7 +127,7 @@ void L1AlgoPulls::AddOne(L1TrackPar& T_, int i, L1HitIndex_t ih)
err.qp = sqrt(tr.C44[i]);
// mc data
- int iMCP = fL1->fvHitPointIndexes[ih];
+ int iMCP = fL1->GetHitBestMcRefs()[ih];
if (iMCP < 0) return;
CbmL1MCPoint& mcP = fL1->fvMCPoints[iMCP];
TL1TrackParameters mc(mcP);
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.cxx b/reco/L1/qa/CbmCaTrackTypeQa.cxx
index 2aad2a6c14..8241263cab 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.cxx
+++ b/reco/L1/qa/CbmCaTrackTypeQa.cxx
@@ -212,7 +212,7 @@ void TrackTypeQa::FillRecoTrack(int iTrkReco, double weight)
// ** First hit **
int iHfst = recoTrack.GetFirstHitIndex();
- int iPfst = (*fpvHits)[iHfst].GetMCPointIndex();
+ int iPfst = (*fpvHits)[iHfst].GetBestMcPointId();
if (iPfst > -1) {
const auto& mcPoint = fpMCData->GetPoint(iPfst);
L1TrackPar trPar(recoTrack.T, recoTrack.C);
@@ -221,7 +221,7 @@ void TrackTypeQa::FillRecoTrack(int iTrkReco, double weight)
// ** Last hit **
int iHlst = recoTrack.GetLastHitIndex();
- int iPlst = (*fpvHits)[iHlst].GetMCPointIndex();
+ int iPlst = (*fpvHits)[iHlst].GetBestMcPointId();
if (iPlst > -1) {
const auto& mcPoint = fpMCData->GetPoint(iPlst);
L1TrackPar trPar(recoTrack.TLast, recoTrack.CLast);
--
GitLab