diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 5f69db7bb0d801d803294a6e4f4c7544a249caa7..14f11c0ee287ae3b9db5552d43c82b6aceecdeca 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -24,6 +24,8 @@ set(SRCS
#### Tracker interfaces (will be moved to core/detector/*) ################
CbmTrackingDetectorInterfaceInit.cxx
###########################################################################
+ CbmCaMCModule.cxx
+ CbmCaTimeSliceReader.cxx
CbmL1.cxx
OffLineInterface/CbmL1StsTrackFinder.cxx
@@ -32,9 +34,6 @@ set(SRCS
CbmL1Util.cxx
CbmL1Performance.cxx
- CbmL1ReadEvent.cxx
- CbmCaMCModule.cxx
- CbmCaTimeSliceReader.cxx
CbmL1MCTrack.cxx
CbmL1Track.cxx
@@ -84,6 +83,7 @@ set(HEADERS
CbmL1Hit.h
CbmL1Track.h
CbmL1Vtx.h
+ CbmCaTimeSliceReader.h
L1Algo/utils/CaUvConverter.h
catools/CaToolsWindowFinder.h
catools/CaToolsLinkKey.h
@@ -135,6 +135,11 @@ ENDIF(SSE_FOUND)
set(LIBRARY_NAME L1)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
set(PUBLIC_DEPENDENCIES
+ CbmMuchBase
+ CbmMvdBase
+ CbmStsBase
+ CbmTofBase
+ CbmTrdBase
CaCore
CbmBase
CbmData
@@ -149,13 +154,8 @@ set(PUBLIC_DEPENDENCIES
)
set(PRIVATE_DEPENDENCIES
- CbmMuchBase
- CbmMvdBase
CbmRecoSts
CbmSimSteer
- CbmStsBase
- CbmTofBase
- CbmTrdBase
CbmRecoBase
KFParticle
external::yaml-cpp
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index 6195dbc621a372810448e7c341a493fd4b271955..437f1beb1dffe661f941d42b8600d0ebe01b1185 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -252,20 +252,24 @@ void MCModule::Finish() { LOG(info) << '\n' << fMonitor.ToString(); }
// ** Reco and MC matching **
// **********************************
+
// ---------------------------------------------------------------------------------------------------------------------
//
-void MCModule::MatchRecoAndMC()
+void MCModule::MatchHits()
{
this->MatchPointsAndHits<ca::EDetectorID::kMvd>();
this->MatchPointsAndHits<ca::EDetectorID::kSts>();
this->MatchPointsAndHits<ca::EDetectorID::kMuch>();
this->MatchPointsAndHits<ca::EDetectorID::kTrd>();
this->MatchPointsAndHits<ca::EDetectorID::kTof>();
- LOG(info) << "DBG 1";
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void MCModule::MatchTracks()
+{
this->MatchRecoAndMCTracks();
- LOG(info) << "DBG 2";
this->InitTrackInfo();
- LOG(info) << "DBG 3";
for (const auto& trkMC : fpMCData->GetTrackContainer()) {
if (trkMC.IsReconstructable()) {
fMonitor.IncrementCounter(EMonitorKey::kMcTrackReconstructable);
@@ -273,6 +277,7 @@ void MCModule::MatchRecoAndMC()
}
}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void MCModule::MatchRecoAndMCTracks()
@@ -579,19 +584,30 @@ void MCModule::ReadMCTracks()
aTrk.SetCharge(pExtMCTrk->GetCharge());
// Set index of mother track. We assume, that mother track was recorded into the internal array before
- int extMotherId = pExtMCTrk->GetMotherId();
+ int extMotherId = pExtMCTrk->GetMotherId();
+ int extChainId = iTrkExt;
+ int extChainParent = extMotherId;
+ while (extChainParent >= 0) {
+ extChainId = extChainParent;
+ const auto* pExtParentTrk = dynamic_cast<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, extChainParent));
+ extChainParent = pExtParentTrk->GetMotherId();
+ }
+
if (extMotherId < 0) {
// This is a primary track, and it's mother ID equals -1 or -2. This value is taken also as an internal mother
- // ID.
+ // ID. The same rules here and below are applied to the chainId
aTrk.SetMotherId(extMotherId);
+ aTrk.SetChainId(extChainId);
}
else {
// This is a secondary track, mother ID should be recalculated for the internal track array.
int motherId = fpMCData->FindInternalTrackIndex(extMotherId, iEvent, iFile);
+ int chainId = fpMCData->FindInternalTrackIndex(extChainId, iEvent, iFile);
if (motherId == -1) {
motherId = -3;
} // Mother is neutral particle, which is rejected
aTrk.SetMotherId(motherId);
+ aTrk.SetChainId(chainId);
}
aTrk.SetProcessId(pExtMCTrk->GetGeantProcessId());
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
index 9d1216fadc7ffb94ada6a4cdc9209f26e5e261f6..cff89752fc81f97107413e7f5be7bec3f24378b4 100644
--- a/reco/L1/CbmCaMCModule.h
+++ b/reco/L1/CbmCaMCModule.h
@@ -107,11 +107,13 @@ namespace cbm::ca
template<ca::EDetectorID DetId>
std::tuple<int, std::vector<int>> MatchHitWithMc(int iHitExt);
+ /// @brief Match reconstructed hits and MC points
+ void MatchHits();
+
/// @brief Match reconstructed and MC data
///
- /// Runs matching of hits with points and reconstructed tracks with MC ones. Reconstructed tracks are updated with
- /// true information.
- void MatchRecoAndMC();
+ /// Runs matching of reconstructed tracks with MC ones. Reconstructed tracks are updated with true information.
+ void MatchTracks();
/// @brief Processes event
///
@@ -154,7 +156,6 @@ namespace cbm::ca
/// @brief Gets verbosity level
int GetVerbosity() const { return fVerbose; }
- private:
/// @brief Check class initialization
/// @note The function throws std::logic_error, if initialization is incomplete
void CheckInit() const;
@@ -183,6 +184,7 @@ namespace cbm::ca
/// will be added as an index of "touched" track.
void MatchRecoAndMCTracks();
+ private:
/// @brief Reads MC tracks from external trees and saves them to MCDataObject
void ReadMCTracks();
@@ -467,9 +469,6 @@ namespace cbm::ca
return;
}
- if ((*fpvFstHitId)[static_cast<int>(DetID)] == (*fpvFstHitId)[static_cast<int>(DetID) + 1]) {
- return;
- }
for (int iH = (*fpvFstHitId)[static_cast<int>(DetID)]; iH < (*fpvFstHitId)[static_cast<int>(DetID) + 1]; ++iH) {
auto& hit = (*fpvQaHits)[iH];
auto [iBestP, vAllP] = MatchHitWithMc<DetID>(hit.ExtIndex);
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index cc0c32551fe203ffa64d882465026adf2df6fb54..c0ba1b46bd8b03d1eb788a0221939701f3a5aef2 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -188,8 +188,7 @@ catch (const std::logic_error& error) {
//
void TimeSliceReader::ReadEvent(CbmEvent* pEvent)
{
- static int iEvt = 0;
- fpEvent = pEvent;
+ fpEvent = pEvent;
this->ReadHits();
this->ReadRecoTracks();
}
@@ -405,7 +404,7 @@ void TimeSliceReader::ReadHits()
}
// Sort debug hits
- if (fpvQaHits) {
+ if (fpvQaHits && fbSortQaHits) {
this->SortQaHits();
}
diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h
index 2b8df810e2bb7255bd5829f392018253bbc86028..e9b3b1c66227783f59568dab5ad6b4088c3dc817 100644
--- a/reco/L1/CbmCaTimeSliceReader.h
+++ b/reco/L1/CbmCaTimeSliceReader.h
@@ -10,11 +10,11 @@
#pragma once
#include "CaConstants.h"
+#include "CaDataManager.h"
#include "CaToolsHitRecord.h"
#include "CaVector.h"
#include "CbmCaMCModule.h"
#include "CbmEvent.h"
-#include "CbmL1.h" // TMP: for CbmL1HitStore
#include "CbmL1DetectorID.h"
#include "CbmL1Hit.h"
#include "CbmL1Track.h"
@@ -130,6 +130,10 @@ namespace cbm::ca
/// @param mode Tracking mode (from ECbmTrackingMode)
void SetTrackingMode(ECbmCaTrackingMode mode) { fTrackingMode = mode; }
+ /// @brief Sets flag to additionally sort QA hits by stations
+ /// @param doSortQaHits true: hits are sorted
+ void SetSortQaHits(bool doSortQaHits) { fbSortQaHits = doSortQaHits; }
+
private:
/// @brief Check class initialization
/// @note The function throws std::logic_error, if initialization is incomplete
@@ -191,6 +195,8 @@ namespace cbm::ca
ECbmCaTrackingMode fTrackingMode = ECbmCaTrackingMode::kSTS; ///< Tracking mode
+ bool fbSortQaHits = false; ///< Flag, if the QA hits must be sorted after reading
+
// Variables for storing cache
int fNofHits = 0; ///< Stored number of hits
int fNofHitKeys = 0; ///< Recorded number of hit keys
@@ -266,6 +272,11 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector()
hitRecord.fT = pHit->GetTime();
hitRecord.fDt2 = pHit->GetTimeError() * pHit->GetTimeError();
+ // Apply hit error smearing according to misalignment
+ hitRecord.fDx2 += fpParameters->GetMisalignmentXsq(DetID);
+ hitRecord.fDy2 += fpParameters->GetMisalignmentYsq(DetID);
+ hitRecord.fDt2 += fpParameters->GetMisalignmentTsq(DetID);
+
std::tie(hitRecord.fRangeX, hitRecord.fRangeY, hitRecord.fRangeT) = pDetInterface->GetHitRanges(*pHit);
hitRecord.fDet = static_cast<int>(DetID);
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 80564cb6fe3e5168427cb86d9080f97fe95aed86..b1b8ab00483a7e3828679ca44894560961a17a53 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -73,6 +73,8 @@ using std::endl;
using std::ios;
using CaTrack = cbm::algo::ca::Track;
using cbm::algo::ca::Parameters;
+using cbm::ca::MCModule;
+using cbm::ca::TimeSliceReader;
using cbm::ca::tools::MaterialHelper;
ClassImp(CbmL1);
@@ -259,8 +261,6 @@ InitStatus CbmL1::Init()
fpStsClusters = nullptr;
- fvSelectedMcEvents.clear();
-
fTimeSlice = (CbmTimeSlice*) fairManager->GetObject("TimeSlice.");
if (!fTimeSlice) {
LOG(fatal) << GetName() << ": No time slice branch in the tree!";
@@ -706,20 +706,57 @@ InitStatus CbmL1::Init()
//
// ** Send formed parameters object to ca::Framework instance **
//
- fpAlgo->ReceiveParameters(fInitManager.TakeParameters());
-
- /*** Get numbers of active stations ***/
+ {
+ auto parameters = fInitManager.TakeParameters();
+ fpParameters = std::make_shared<ca::Parameters>(parameters);
+ fpAlgo->ReceiveParameters(std::move(parameters));
+ }
- fNMvdStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kMvd);
- fNStsStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kSts);
- fNTrdStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kTrd);
- fNMuchStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kMuch);
- fNTofStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kTof);
- fNStations = fpAlgo->GetParameters().GetNstationsActive();
+ // Initialize time-slice reader
+ fpTSReader = std::make_unique<TimeSliceReader>();
+ fpTSReader->SetDetector(ca::EDetectorID::kMvd, fUseMVD);
+ fpTSReader->SetDetector(ca::EDetectorID::kSts, fUseSTS);
+ fpTSReader->SetDetector(ca::EDetectorID::kMuch, fUseMUCH);
+ fpTSReader->SetDetector(ca::EDetectorID::kTrd, fUseTRD);
+ fpTSReader->SetDetector(ca::EDetectorID::kTof, fUseTOF);
+
+ fpTSReader->RegisterParameters(fpParameters);
+ fpTSReader->RegisterHitIndexContainer(fvExternalHits);
+ fpTSReader->RegisterIODataManager(fpIODataManager);
+ fpTSReader->RegisterQaHitContainer(fvHitDebugInfo);
+ if (!fpTSReader->InitRun()) {
+ return kFATAL;
+ }
- LOG(info) << fpAlgo->GetParameters().ToString(1);
+ if (fPerformance) {
+ fpMCModule = std::make_unique<MCModule>(fVerbose, fPerformance);
+ fpMCModule->SetDetector(ca::EDetectorID::kMvd, fUseMVD);
+ fpMCModule->SetDetector(ca::EDetectorID::kSts, fUseSTS);
+ fpMCModule->SetDetector(ca::EDetectorID::kMuch, fUseMUCH);
+ fpMCModule->SetDetector(ca::EDetectorID::kTrd, fUseTRD);
+ fpMCModule->SetDetector(ca::EDetectorID::kTof, fUseTOF);
+
+ fpMCModule->RegisterMCData(fMCData);
+ fpMCModule->RegisterRecoTrackContainer(fvRecoTracks);
+ fpMCModule->RegisterHitIndexContainer(fvExternalHits);
+ fpMCModule->RegisterParameters(fpParameters);
+ fpMCModule->RegisterQaHitContainer(fvHitDebugInfo);
+ fpMCModule->RegisterFirstHitIndexes(fpTSReader->GetHitFirstIndexDet());
+ if (!fpMCModule->InitRun()) {
+ return kFATAL;
+ }
+ }
+ LOG(info) << fpParameters->ToString(1);
LOG(info) << "----- Numbers of stations active in tracking -----";
+ fNMvdStations = fpParameters->GetNstationsActive(ca::EDetectorID::kMvd);
+ fNStsStations = fpParameters->GetNstationsActive(ca::EDetectorID::kSts);
+ fNTrdStations = fpParameters->GetNstationsActive(ca::EDetectorID::kTrd);
+ fNMuchStations = fpParameters->GetNstationsActive(ca::EDetectorID::kMuch);
+ fNTofStations = fpParameters->GetNstationsActive(ca::EDetectorID::kTof);
+ fNStations = fpParameters->GetNstationsActive();
+
+
LOG(info) << " MVD: " << fNMvdStations;
LOG(info) << " STS: " << fNStsStations;
LOG(info) << " MuCh: " << fNMuchStations;
@@ -750,67 +787,22 @@ InitStatus CbmL1::Init()
void CbmL1::Reconstruct(CbmEvent* event)
{
- LOG(info) << "\n\n\n\n\n !!!!! EVENT = " << event << "\n\n\n\n\n";
+ fpTSReader->ReadEvent(event);
+ LOG(info) << "READ: hit ids = " << fvExternalHits.size() << ", hits = " << fpIODataManager->GetNofHits()
+ << ", dbg hits = " << fvHitDebugInfo.size();
ca::TrackingMonitorData monitorData;
- fvSelectedMcEvents.clear();
-
- int bestMcFile = -1;
- int bestMcEvent = -1;
-
- if (fPerformance) {
- if (event) {
- CbmMatch* match = event->GetMatch();
- if (!match) {
- LOG(error) << "CbmL1: match between reco and MC events missing!! Performance can not be evaluated!!";
- fPerformance = 0;
- }
- else {
- cout << "CbmL1: mc events all " << fpMcEventList->GetNofEvents() << " mc events linked " << match->GetNofLinks()
- << std::endl;
- for (int iLink = 0; iLink < match->GetNofLinks(); iLink++) {
- const CbmLink& link = match->GetLink(iLink);
- fvSelectedMcEvents.insert(DFSET::value_type(link.GetFile(), link.GetEntry()));
- cout << "CbmL1: linked mc event file " << link.GetFile() << " event " << link.GetEntry() << std::endl;
- }
- if (match->GetNofLinks() > 1) {
- const CbmLink& link = match->GetMatchedLink();
- bestMcFile = link.GetFile();
- bestMcEvent = link.GetEntry();
- }
- }
- }
- else {
- int nofEvents = fpMcEventList->GetNofEvents();
- for (int iE = 0; iE < nofEvents; iE++) {
- int fileId = fpMcEventList->GetFileIdByIndex(iE);
- int eventId = fpMcEventList->GetEventIdByIndex(iE);
- fvSelectedMcEvents.insert(DFSET::value_type(fileId, eventId));
- }
- }
- }
-
- if (fVerbose > 1) {
- if (event) {
- LOG(info) << "\nCbmL1::Exec event " << fEventNo << " ...\n";
- }
- else {
- LOG(info) << "\nCbmL1::Exec time slice " << fEventNo << " ...\n";
- }
- }
+ LOG_IF(info, fVerbose > 1) << "\nCbmL1::Exec " << (event ? "event " : "time slice ") << fEventNo << " ...\n";
// ----- Read data from branches and send data from IODataManager to ca::Framework ----------------------------------------
- ReadEvent(event);
- //if (fPerformance) {
- // LOG(info) << "**** DEBUG: BEGIN";
- // LOG(info) << fvMCPoints[0].ToString(10, true); // header
- // for (const auto& point: fvMCPoints) { LOG(info) << point.ToString(10); }
- // LOG(info) << "**** DEBUG: END";
- //}
+ if (fPerformance) {
+ fpMCModule->InitEvent(event); // reads points and MC tracks
+ fpMCModule->MatchHits(); // matches points with hits
+ }
// Material monitoring: mark active areas
{
@@ -819,37 +811,19 @@ void CbmL1::Reconstruct(CbmEvent* event)
fMaterialMonitor[h.Station()].MarkActiveBin(h.X(), h.Y());
}
}
- //
-
- if (fPerformance) {
- HitMatch();
- // calculate the max number of Hits\mcPoints on continuous(consecutive) stations
- for (auto it = fvMCTracks.begin(); it != fvMCTracks.end(); ++it) {
- it->Init();
- }
- if (bestMcFile >= 0) { // suppress mc tracks from complementary mc events
- for (auto it = fvMCTracks.begin(); it != fvMCTracks.end(); ++it) {
- if (it->iFile != bestMcFile || it->iEvent != bestMcEvent) {
- it->SetIsReconstructable(false);
- }
- }
- }
- }
+ LOG(info) << "CHECK: hit ids = " << fvExternalHits.size() << ", hits = " << fpIODataManager->GetNofHits()
+ << ", dbg hits = " << fvHitDebugInfo.size();
// FieldApproxCheck();
// FieldIntegralCheck();
-
+ fpAlgo->ReceiveInputData(fpIODataManager->TakeInputData());
fpAlgo->SetMonitorData(monitorData);
- if (fVerbose > 1) {
- LOG(info) << "L1 Track finder...";
- }
+ LOG_IF(info, fVerbose > 1) << "L1 Track finder...";
fpAlgo->fTrackFinder.FindTracks();
// IdealTrackFinder();
fTrackingTime = fpAlgo->fCaRecoTime; // TODO: remove (not used)
- if (fVerbose > 1) {
- LOG(info) << "L1 Track finder ok";
- }
+ LOG_IF(info, fVerbose > 1) << "L1 Track finder ok";
// Update monitor data after the actual tracking
monitorData = fpAlgo->GetMonitorData();
@@ -861,6 +835,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
int trackFirstHit = 0;
+ // FIXME: Rewrite
for (auto it = fpAlgo->fRecoTracks.begin(); it != fpAlgo->fRecoTracks.end(); trackFirstHit += it->fNofHits, it++) {
CbmL1Track t;
t.Set(it->fParFirst);
@@ -886,6 +861,8 @@ void CbmL1::Reconstruct(CbmEvent* event)
LOG(info) << "Performance...";
}
+ fpMCModule->MatchTracks(); // matches reco and MC tracks, fills the MC-truth fields of reco tracks
+
//
// tracker input performance is moved to external QA tasks.
// InputPerformance() method is not optimised to run with the event builder
@@ -893,54 +870,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
// InputPerformance();
//
- TrackMatch();
- // ****** DEBUG: BEGIN
- if constexpr (0) {
- using std::setw;
- {
- int id = 0;
- if (fvRecoTracks.size()) {
- LOG(info) << setw(4) << "No." << ' ' << fvRecoTracks[0].ToString(3, true);
- }
- for (const auto& trk : fvRecoTracks) {
- LOG(info) << setw(4) << (id++) << ' ' << trk.ToString(3);
- }
- }
- {
- LOG(info) << "---------- Reco hit sample";
- int id = 0;
- if (fvHitDebugInfo.size()) {
- LOG(info) << setw(4) << "No." << ' ' << fvHitDebugInfo[0].ToString(3, true);
- }
- for (const auto& h : fvHitDebugInfo) {
- LOG(info) << setw(4) << (id++) << ' ' << h.ToString(3);
- }
- }
- if (fPerformance) {
- LOG(info) << "---------- MC tracks sample";
- int id = 0;
- if (fvMCTracks.size()) {
- LOG(info) << setw(4) << "No." << ' ' << fvMCTracks[0].ToString(3, true);
- }
- for (const auto& trk : fvMCTracks) {
- if (trk.Hits.size() || trk.Points.size()) {
- LOG(info) << setw(4) << (id++) << ' ' << trk.ToString(3);
- }
- }
- }
- if (fPerformance) {
- LOG(info) << "---------- MC points sample";
- int id = 0;
- if (fvMCPoints.size()) {
- LOG(info) << setw(4) << "No." << ' ' << fvMCPoints[0].ToString(3, true);
- }
- for (const auto& p : fvMCPoints) {
- LOG(info) << setw(4) << (id++) << ' ' << p.ToString(3);
- }
- }
- }
- // ****** DEBUG: END
EfficienciesPerformance();
HistoPerformance();
TrackFitPerformance();
@@ -950,23 +880,8 @@ void CbmL1::Reconstruct(CbmEvent* event)
// TimeHist();
/// WriteSIMDKFData();
}
- if (fVerbose > 1) {
- LOG(info) << "Tracking performance... done";
- }
- if (fVerbose > 1) {
- LOG(info) << "End of CA";
- }
-
- static bool ask = 0;
- char symbol;
- if (ask) {
- LOG(info) << "\nCA run (any/r/q) > ";
- do {
- std::cin.get(symbol);
- if (symbol == 'r') ask = false;
- if ((symbol == 'e') || (symbol == 'q')) exit(0);
- } while (symbol != '\n');
- }
+ LOG_IF(info, fVerbose > 1) << "Tracking performance... done";
+ LOG_IF(info, fVerbose > 1) << "End of CA";
++fEventNo;
fMonitor.AddMonitorData(monitorData);
@@ -1112,22 +1027,21 @@ void CbmL1::IdealTrackFinder()
fpAlgo->fRecoTracks.clear();
fpAlgo->fRecoHits.clear();
- for (auto& MC : fvMCTracks) {
+ for (auto& MC : fMCData.GetTrackContainer()) {
if (!MC.IsReconstructable()) continue;
- if (!(MC.ID >= 0)) continue;
+ if (!(MC.GetId() >= 0)) continue;
- if (MC.Hits.size() < 4) continue;
+ if (MC.GetNofHits() < 4) continue;
CaTrack algoTr;
algoTr.fNofHits = 0;
ca::Vector<int> hitIndices("CbmL1::hitIndices", fpAlgo->GetParameters().GetNstationsActive(), -1);
- for (unsigned int iH = 0; iH < MC.Hits.size(); iH++) {
- const int hitI = MC.Hits[iH];
- const CbmL1HitDebugInfo& hit = fvHitDebugInfo[hitI];
+ for (unsigned int iH : MC.GetHitIndexes()) {
+ const CbmL1HitDebugInfo& hit = fvHitDebugInfo[iH];
const int iStation = hit.GetStationId();
- if (iStation >= 0) hitIndices[iStation] = hitI;
+ if (iStation >= 0) hitIndices[iStation] = iH;
}
@@ -1148,13 +1062,13 @@ void CbmL1::IdealTrackFinder()
fpAlgo->fRecoHits.push_back(hitI);
}
- algoTr.fParFirst.X() = MC.x;
- algoTr.fParFirst.Y() = MC.y;
- algoTr.fParFirst.Z() = MC.z;
- algoTr.fParFirst.Tx() = MC.px / MC.pz;
- algoTr.fParFirst.Ty() = MC.py / MC.pz;
- algoTr.fParFirst.Qp() = MC.q / MC.p;
- algoTr.fParFirst.Time() = MC.time;
+ algoTr.fParFirst.X() = MC.GetStartX();
+ algoTr.fParFirst.Y() = MC.GetStartY();
+ algoTr.fParFirst.Z() = MC.GetStartZ();
+ algoTr.fParFirst.Tx() = MC.GetTx();
+ algoTr.fParFirst.Ty() = MC.GetTy();
+ algoTr.fParFirst.Qp() = MC.GetCharge() / MC.GetP();
+ algoTr.fParFirst.Time() = MC.GetStartT();
fpAlgo->fRecoTracks.push_back(algoTr);
}
} // void CbmL1::IdealTrackFinder()
@@ -1286,9 +1200,8 @@ void CbmL1::WriteSIMDKFData()
for (auto RecTrack = fvRecoTracks.begin(); RecTrack != fvRecoTracks.end(); ++RecTrack) {
if (RecTrack->IsGhost()) continue;
-
- CbmL1MCTrack* MCTrack = RecTrack->GetMCTrack();
- if (!(MCTrack->IsPrimary())) continue;
+ const auto& mcTrk = fMCData.GetTrack(RecTrack->GetMatchedMCTrackIndex());
+ if (!(mcTrk.IsPrimary())) continue;
int NHits = (RecTrack->Hits).size();
float x[20], y[20], z[20];
@@ -1313,8 +1226,9 @@ void CbmL1::WriteSIMDKFData()
}
}
- Tracks << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px << " "
- << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NHits << endl;
+ Tracks << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
+ << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " " << NHits
+ << endl;
float AngleXAxis = 0, AngleYAxis = 0;
for (int i = 0; i < NHits; i++)
@@ -1322,23 +1236,30 @@ void CbmL1::WriteSIMDKFData()
<< endl;
Tracks << endl;
- int NMCPoints = (MCTrack->Points).size();
-
- McTracksIn << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px << " "
- << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl;
- McTracksOut << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px << " "
- << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl;
- McTracksCentr << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px
- << " " << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl;
-
- for (int iPoint = 0; iPoint < NMCPoints; iPoint++) {
- CbmL1MCPoint& MCPoint = fvMCPoints[MCTrack->Points[iPoint]];
- McTracksIn << " " << MCPoint.iStation << " " << MCPoint.xIn << " " << MCPoint.yIn << " " << MCPoint.zIn << " "
- << MCPoint.pxIn << " " << MCPoint.pyIn << " " << MCPoint.pzIn << endl;
- McTracksOut << " " << MCPoint.iStation << " " << MCPoint.xOut << " " << MCPoint.yOut << " " << MCPoint.zOut
- << " " << MCPoint.pxOut << " " << MCPoint.pyOut << " " << MCPoint.pzOut << endl;
- McTracksCentr << " " << MCPoint.iStation << " " << MCPoint.x << " " << MCPoint.y << " " << MCPoint.z << " "
- << MCPoint.px << " " << MCPoint.py << " " << MCPoint.pz << endl;
+ int NMCPoints = mcTrk.GetNofPoints();
+
+ McTracksIn << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
+ << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " "
+ << NMCPoints << endl;
+ McTracksOut << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
+ << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " "
+ << NMCPoints << endl;
+ McTracksCentr << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
+ << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " "
+ << NMCPoints << endl;
+
+ for (int iPoint : mcTrk.GetPointIndexes()) {
+ const auto& mcPoint = fMCData.GetPoint(iPoint);
+ McTracksIn << " " << mcPoint.GetStationId() << " " << mcPoint.GetXIn() << " " << mcPoint.GetYIn() << " "
+ << mcPoint.GetZIn() << " " << mcPoint.GetPxIn() << " " << mcPoint.GetPyIn() << " " << mcPoint.GetPzIn()
+ << endl;
+ McTracksOut << " " << mcPoint.GetStationId() << " " << mcPoint.GetXOut() << " " << mcPoint.GetYOut() << " "
+ << mcPoint.GetZOut() << " " << mcPoint.GetPxOut() << " " << mcPoint.GetPyOut() << " "
+ << mcPoint.GetPzOut() << endl;
+ McTracksCentr << " " << mcPoint.GetStationId() << " " << mcPoint.GetX() << " " << mcPoint.GetY() << " "
+ << mcPoint.GetZ() << " " << mcPoint.GetPx() << " " << mcPoint.GetPy() << " " << mcPoint.GetPz()
+ << endl;
+ ;
}
McTracksIn << endl;
McTracksOut << endl;
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 1ee02e02c5c1f6b9836f4fa19c4c86ffa1cc5cae..077a6c24ce4e1e2142e5ae65edbcbd81428dfd6e 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -21,7 +21,6 @@
#ifndef _CbmL1_h_
#define _CbmL1_h_
-
#include "AlgoFairloggerCompat.h"
#include "CaDataManager.h"
#include "CaFramework.h"
@@ -30,6 +29,8 @@
#include "CaMonitor.h"
#include "CaTrackParam.h"
#include "CaVector.h"
+#include "CbmCaMCModule.h"
+#include "CbmCaTimeSliceReader.h"
#include "CbmL1DetectorID.h"
#include "CbmL1Hit.h"
#include "CbmL1MCPoint.h"
@@ -65,6 +66,7 @@ class TProfile2D;
class TNtuple;
class TGeoNode;
+
namespace
{
namespace ca = cbm::algo::ca;
@@ -167,7 +169,6 @@ class CbmL1 : public FairTask {
/// Pointer to CbmL1 instance
static CbmL1* Instance() { return fpInstance; }
-
// ** Member functions override from FairTask **
/// Defines action in the beginning of the run (initialization)
@@ -179,6 +180,8 @@ class CbmL1 : public FairTask {
/// Defines action in the end of the run (saves results)
void Finish();
+ /// \brief Gets reference to MC data object
+ const cbm::ca::tools::MCData& GetMCData() const { return fMCData; }
// ** Specific member functions **
@@ -236,6 +239,10 @@ class CbmL1 : public FairTask {
/// Gets flag: to correct input hits on MC or not
bool GetCorrectHitsOnMC() const { return fIfCorrectHitsOnMC; }
+ /// Gets vector of the extended QA hits
+ /// TODO: It is a temporary function, do not rely on it
+ const auto& GetQaHits() const { return fvHitDebugInfo; }
+
/// Utility to map the L1DetectorID items into detector names
static constexpr const char* GetDetectorName(ca::EDetectorID detectorID)
{
@@ -275,17 +282,6 @@ class CbmL1 : public FairTask {
/// \param event Pointer to current CbmEvent object
void Reconstruct(CbmEvent* event = nullptr);
-
- // static bool compareZ(const int &a, const int &b );
- inline double Get_Z_vMCPoint(int a) const { return fvMCPoints[a].z; }
-
- const ca::Vector<CbmL1MCPoint>& GetMcPoints() const { return fvMCPoints; }
-
- const ca::Vector<CbmL1MCTrack>& GetMcTracks() const { return fvMCTracks; }
-
- const ca::Vector<std::vector<int>>& GetHitMcRefs() const { return fvHitPointIndices; }
- const ca::Vector<int>& GetHitBestMcRefs() const { return fvHitBestPointIndices; }
-
static double boundedGaus(double sigma);
/// Obsolete setters to be removed
@@ -338,20 +334,6 @@ class CbmL1 : public FairTask {
/// @brief Runs ideal track finder: copies all MC-tracks into reconstructed tracks
void IdealTrackFinder();
- /// @brief Reads input data in the event
- /// @param event Pointer to the current CbmEvent object
- ///
- /// Reads information about hits, mcPoints and mcTracks into L1 classes. Repacks data from the external data branches
- /// to the internal data arrays.
- void ReadEvent(CbmEvent* event = nullptr);
-
- /// Converts data from generic FairMCPoint based class to the CbmL1MCPoint
- /// \param iDet Index of the detector subsystem
- /// \param file Index of the input file
- /// \param event Index of the input event
- /// \param iPoint Index of the point into the input MC points CbmMCDataArray object for the particular detector
- void ReadMCPoint(ca::EDetectorID det, int iPoint, int file, int event);
-
// static bool compareZ(const int &a, const int &b );
// bool compareZ(const int &a, const int &b );
@@ -369,11 +351,6 @@ class CbmL1 : public FairTask {
template<ca::EDetectorID DetId>
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
- /// should be called after reading the event
- void HitMatch();
-
void FieldApproxCheck(); // Build histograms with difference between Field map and approximated field
void FieldIntegralCheck(); // Build 2D histogram: dependence of the field integral on phi and theta
void InputPerformance(); // Build histograms about input data, like hit pulls, etc.
@@ -387,10 +364,6 @@ class CbmL1 : public FairTask {
/// @param seed Random seed
void SetRandomSeed(unsigned int seed) { fInitManager.SetRandomSeed(seed); }
- /// Matches reconstructed and MC tracks
- /// \note Should be called before Performances
- void TrackMatch();
-
/// Calculates tracking efficiencies (counters)
void EfficienciesPerformance();
@@ -398,7 +371,7 @@ class CbmL1 : public FairTask {
/// \note Should be called only after CbmL1::Performance()
void TrackFitPerformance();
- void FillFitHistos(cbm::algo::ca::TrackParamV& tr, const CbmL1MCPoint& mc, bool isTimeFitted, TH1F* h[]);
+ void FillFitHistos(cbm::algo::ca::TrackParamV& tr, const cbm::ca::tools::MCPoint& mc, bool isTimeFitted, TH1F* h[]);
/// Fills performance histograms
void HistoPerformance();
@@ -461,13 +434,17 @@ class CbmL1 : public FairTask {
void DumpMaterialToFile(TString fileName);
private:
- std::string fInputDataFilename = ""; ///< File name to read/write input hits
- std::string fsInputSearchWindowsFilename = ""; ///< File name to read search windows
-
// ***************************
// ** Member variables list **
// ***************************
+ std::string fInputDataFilename = ""; ///< File name to read/write input hits
+ std::string fsInputSearchWindowsFilename = ""; ///< File name to read search windows
+
+ std::unique_ptr<cbm::ca::TimeSliceReader> fpTSReader = nullptr; ///< event/TS reader
+ std::unique_ptr<cbm::ca::MCModule> fpMCModule = nullptr; ///< MC module
+ std::shared_ptr<ca::Parameters> fpParameters = nullptr; ///< External parameters object
+ cbm::ca::tools::MCData fMCData; ///< MC Data object
ca::InitManager fInitManager; ///< Tracking parameters data manager
std::shared_ptr<ca::DataManager> fpIODataManager = nullptr; ///< Input-output data manager
@@ -481,11 +458,7 @@ class CbmL1 : public FairTask {
ca::Framework* fpAlgo = nullptr; ///< Pointer to the L1 track finder algorithm
- ca::Framework::TrackingMode fTrackingMode =
- ca::Framework::TrackingMode::kSts; ///< Tracking mode: kSts, kMcbm or kGlobal
-
- DFSET fvSelectedMcEvents{}; ///< Set of selected MC events with fileID and eventId
-
+ ca::Framework::TrackingMode fTrackingMode = ca::Framework::TrackingMode::kSts; ///< Tracking mode
ca::Vector<CbmL1Track> fvRecoTracks = {"CbmL1::fvRecoTracks"}; ///< Reconstructed tracks container
@@ -499,20 +472,6 @@ class CbmL1 : public FairTask {
double fTargetY{1.e10}; ///< target position Y
double fTargetZ{1.e10}; ///< target position Z
- int fNpointsMvd = 0; ///< Number of MC points for MVD
- int fNpointsSts = 0; ///< Number of MC points for STS
- int fNpointsMuch = 0; ///< Number of MC points for MuCh
- int fNpointsTrd = 0; ///< Number of MC points for TRD
- int fNpointsTof = 0; ///< Number of MC points for TOF
-
- int fNpointsMvdAll = 0; ///< Number of MC points for MVD
- int fNpointsStsAll = 0; ///< Number of MC points for STS
- int fNpointsMuchAll = 0; ///< Number of MC points for MuCh
- int fNpointsTrdAll = 0; ///< Number of MC points for TRD
- int fNpointsTofAll = 0; ///< Number of MC points for TOF
-
- ca::Vector<CbmL1MCPoint> fvMCPoints = {"CbmL1::fvMCPoints"}; ///< Container of MC points
-
int fNStations = 0; ///< number of total active detector stations
int fNMvdStations = 0; ///< number of active MVD stations
int fNStsStations = 0; ///< number of active STS stations
@@ -620,6 +579,8 @@ class CbmL1 : public FairTask {
TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches
TClonesArray* fpMuchDigiMatches = nullptr; ///< Array of MuCh digi matches (NOTE: currently unsused)
+ ca::Vector<CbmL1MCPoint> fvMCPoints = {"CbmL1::fvMCPoints"}; ///< Container of MC points
+
ca::Vector<CbmL1MCTrack> fvMCTracks = {"CbmL1::fvMCTracks"}; ///< Array of MC tracks
ca::Vector<std::vector<int>> fvHitPointIndices = {
@@ -627,6 +588,7 @@ class CbmL1 : public FairTask {
ca::Vector<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
@@ -645,6 +607,7 @@ class CbmL1 : public FairTask {
std::unordered_map<CbmL1LinkKey, int> fmMCPointsLinksMap = {}; /// Internal MC point index vs. link
std::unordered_map<CbmL1LinkKey, int> fmMCTracksLinksMap = {}; /// Internal MC track index vs. link
+
cbm::ca::DetIdArr_t<std::set<int>> fvmDisabledStationIDs; /// Array of local indices of disabled tracking stations
// *****************************
@@ -679,66 +642,7 @@ class CbmL1 : public FairTask {
ClassDef(CbmL1, 0);
};
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<ca::EDetectorID DetID>
-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
- if constexpr (ca::EDetectorID::kMvd == DetID) {
- pHitMatch = dynamic_cast<CbmMatch*>(fpMvdHitMatches->At(iHitExt));
- indexShift = 0;
- }
- else if constexpr (ca::EDetectorID::kSts == DetID) {
- pHitMatch = dynamic_cast<CbmMatch*>(fpStsHitMatches->At(iHitExt));
- indexShift = fNpointsMvdAll;
- }
- else if constexpr (ca::EDetectorID::kMuch == DetID) {
- pHitMatch = dynamic_cast<CbmMatch*>(fpMuchHitMatches->At(iHitExt));
- indexShift = fNpointsMvdAll + fNpointsStsAll;
- }
- else if constexpr (ca::EDetectorID::kTrd == DetID) {
- pHitMatch = dynamic_cast<CbmMatch*>(fpTrdHitMatches->At(iHitExt));
- indexShift = fNpointsMvdAll + fNpointsStsAll + fNpointsMuchAll;
- }
- else if constexpr (ca::EDetectorID::kTof == DetID) {
- pHitMatch = dynamic_cast<CbmMatch*>(fpTofHitMatches->At(iHitExt));
- indexShift = fNpointsMvdAll + fNpointsStsAll + fNpointsMuchAll + fNpointsTrdAll;
- }
-
- 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) {
- 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()) {
- iPointInt = itPoint->second;
- }
- }
- }
- return std::tuple(iPointInt, vPoints);
-}
#endif //_CbmL1_h_
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 51e05993e5367363e50abba3271f69db63f35b38..1b2d80943d2110d1125ae6b88555eb3b8a76a98a 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -67,96 +67,6 @@ using std::pair;
using std::setw;
using std::vector;
-void CbmL1::TrackMatch()
-{
- map<int, CbmL1MCTrack*> pMCTrackMap;
- pMCTrackMap.clear();
-
- // ***** DEBUG: Start
- for (int iTrk = 0; iTrk < int(fvMCTracks.size()); ++iTrk) {
- assert(iTrk == fvMCTracks[iTrk].ID);
- }
- // ***** DEBUG: End
-
- // fill pMCTrackMap
- for (vector<CbmL1MCTrack>::iterator i = fvMCTracks.begin(); i != fvMCTracks.end(); ++i) {
- CbmL1MCTrack& MC = *i;
- if (pMCTrackMap.find(MC.ID) == pMCTrackMap.end()) {
- pMCTrackMap.insert(pair<int, CbmL1MCTrack*>(MC.ID, &MC));
- }
- else {
- cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl;
- }
- }
-
- // -- prepare information about reconstructed tracks
- const int nRTracks = fvRecoTracks.size();
- for (int iR = 0; iR < nRTracks; iR++) {
- CbmL1Track* prtra = &(fvRecoTracks[iR]);
-
- // cout<<iR<<" iR"<<endl;
-
- int hitsum = prtra->Hits.size(); // number of hits in track
-
- // count how many hits from each mcTrack belong to current recoTrack
- 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) {
- 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
-
- // RTrack <-> MCTrack identification
- double max_percent = 0.0; // [%]. maximum persent of hits, which belong to one mcTrack.
- for (map<int, int>::iterator posIt = hitmap.begin(); posIt != hitmap.end();
- posIt++) { // loop over all touched MCTracks
-
- 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 (max_percent < purity) {
- max_percent = purity;
- }
-
- // set relation to the mcTrack
- assert(pMCTrackMap.find(iMcTrack) != pMCTrackMap.end());
-
- CbmL1MCTrack* pmtra = pMCTrackMap[iMcTrack];
-
- if (purity >= CbmL1Constants::MinPurity) { // found correspondent MCTrack
- pmtra->AddRecoTrack(prtra);
- pmtra->AddRecoTrackIndex(iR);
- prtra->AddMCTrack(pmtra);
- }
- else {
- pmtra->AddTouchTrack(prtra);
- pmtra->AddTouchTrackIndex(iR);
- }
- } // for hitmap
- prtra->SetMaxPurity(max_percent);
-
- } // for reco tracks
-} //
-
-
struct TL1PerfEfficiencies : public TL1Efficiencies {
TL1PerfEfficiencies()
: TL1Efficiencies()
@@ -335,9 +245,10 @@ void CbmL1::EfficienciesPerformance()
}
cout << endl;
for (map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++) {
- CbmL1MCTrack& t = fvMCTracks[posIt->first];
- cout << "mc " << posIt->first << " pdg " << t.pdg << " mother: " << t.mother_ID << " chain " << t.chainID;
- cout << " n mc stations: " << t.NMCStations() << endl;
+ const auto& mcTrk = fMCData.GetTrack(posIt->first);
+ cout << "mc " << posIt->first << " pdg " << mcTrk.GetPdgCode() << " mother: " << mcTrk.GetMotherId()
+ << " chain " << mcTrk.GetChainId();
+ cout << " n mc stations: " << mcTrk.GetNofConsStationsWithPoint() << endl;
}
for (unsigned int i = 0; i < rtraIt->Hits.size(); i++) {
const ca::Hit& h = fpAlgo->GetInputData().GetHit(rtraIt->Hits[i]);
@@ -359,36 +270,34 @@ void CbmL1::EfficienciesPerformance()
sta_nfakes[i] = 0;
}
- for (auto mtraIt = fvMCTracks.begin(); mtraIt != fvMCTracks.end(); mtraIt++) {
- CbmL1MCTrack& mtra = *(mtraIt);
-
- // if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
- if (!mtra.IsReconstructable() && !mtra.IsAdditional()) continue;
+ for (const auto& mcTrk : fMCData.GetTrackContainer()) {
+ // if( !( mcTrk.GetPdgCode() == -11 && mcTrk.GetMotherId() == -1 ) ) continue; // electrons only
+ if (!mcTrk.IsReconstructable() && !mcTrk.IsAdditional()) continue;
// -- find used constans --
// is track reconstructed
- const bool reco = (mtra.IsReconstructed());
+ const bool reco = (mcTrk.IsReconstructed());
// is track killed. At least one hit of it belong to some recoTrack
- const bool killed = !reco && mtra.IsDisturbed();
+ const bool killed = !reco && mcTrk.IsDisturbed();
// ration length for current mcTrack
- auto& rTracks = mtra.GetRecoTracks(); // for length calculations
double ratio_length = 0;
double ratio_fakes = 0;
- double mc_length_hits = mtra.NStations();
+ double mc_length_hits = mcTrk.GetTotNofStationsWithHit();
- int mc_length = mtra.NMCStations();
+ int mc_length = mcTrk.GetTotNofStationsWithPoint();
if (reco) {
- for (unsigned int irt = 0; irt < rTracks.size(); irt++) {
- ratio_length += static_cast<double>(rTracks[irt]->GetNofHits()) * rTracks[irt]->GetMaxPurity() / mc_length_hits;
- ratio_fakes += 1 - rTracks[irt]->GetMaxPurity();
+ for (unsigned int irt : mcTrk.GetRecoTrackIndexes()) {
+ const auto& rTrk = fvRecoTracks[irt];
+ ratio_length += static_cast<double>(rTrk.GetNofHits()) * rTrk.GetMaxPurity() / mc_length_hits;
+ ratio_fakes += 1 - rTrk.GetMaxPurity();
}
}
// number of clones
int nclones = 0;
- if (reco) nclones = mtra.GetNClones();
+ if (reco) nclones = mcTrk.GetNofClones();
// if (nclones){ // Debug. Look at clones
// for (int irt = 0; irt < rTracks.size(); irt++){
// const int ista = fvHitDebugInfo[rTracks[irt]->Hits[0]].iStation;
@@ -397,9 +306,9 @@ void CbmL1::EfficienciesPerformance()
// cout << mtra.NStations() << endl;
// }
- if (mtra.IsAdditional()) { // short
+ if (mcTrk.IsAdditional()) { // short
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "short");
- switch (mtra.pdg) {
+ switch (mcTrk.GetPdgCode()) {
case 11:
case -11:
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortE");
@@ -424,15 +333,15 @@ void CbmL1::EfficienciesPerformance()
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total");
- if (mtra.isSignal) { // D0
+ if (mcTrk.IsSignal()) { // D0
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "d0");
}
- if (mtra.p > CbmL1Constants::MinFastMom) { // fast tracks
+ if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { // fast tracks
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast");
- if (mtra.IsPrimary()) { // fast primary
- if (mtra.NStations() == fNStations) { // long fast primary
+ if (mcTrk.IsPrimary()) { // fast primary
+ if (mcTrk.GetTotNofStationsWithHit() == fNStations) { // long fast primary
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "long_fast_prim");
}
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_prim");
@@ -444,7 +353,7 @@ void CbmL1::EfficienciesPerformance()
else { // slow set of tracks
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow");
- if (mtra.IsPrimary()) { // slow primary
+ if (mcTrk.IsPrimary()) { // slow primary
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_prim");
}
else {
@@ -452,13 +361,13 @@ void CbmL1::EfficienciesPerformance()
}
} // if slow
- if (mtra.pdg == 11 || mtra.pdg == -11) {
+ if (mcTrk.GetPdgCode() == 11 || mcTrk.GetPdgCode() == -11) {
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total_e");
- if (mtra.p > CbmL1Constants::MinFastMom) { // fast tracks
+ if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { // fast tracks
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_e");
- if (mtra.IsPrimary()) { // fast primary
+ if (mcTrk.IsPrimary()) { // fast primary
}
else { // fast secondary
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec_e");
@@ -467,7 +376,7 @@ void CbmL1::EfficienciesPerformance()
else { // slow set of tracks
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_e");
- if (mtra.IsPrimary()) { // slow primary
+ if (mcTrk.IsPrimary()) { // slow primary
}
else {
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec_e");
@@ -866,8 +775,9 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
h_ghost_chi2->Fill(prtra->GetChiSq() / prtra->GetNdf());
h_ghost_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
}
- else if (prtra->GetNMCTracks() > 0) {
- if (prtra->GetMCTrack()[0].IsReconstructable()) {
+ else if (prtra->GetNofMCTracks() > 0) {
+ const auto& mcTrk = fMCData.GetTrack(prtra->GetMatchedMCTrackIndex());
+ if (mcTrk.IsReconstructable()) {
h_reco_chi2->Fill(prtra->GetChiSq() / prtra->GetNdf());
h_reco_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
}
@@ -912,29 +822,28 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
} // for reco tracks
int mc_total = 0; // total amount of reconstructable mcTracks
- for (vector<CbmL1MCTrack>::iterator mtraIt = fvMCTracks.begin(); mtraIt != fvMCTracks.end(); mtraIt++) {
- CbmL1MCTrack& mtra = *(mtraIt);
- // if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
+ for (const auto& mcTrk : fMCData.GetTrackContainer()) {
+ // if( !( mcTrk.GetPdgCode() == -11 && mcTrk.GetPdgCode() == -1 ) ) continue; // electrons only
// No Sts hits?
- int nmchits = mtra.Hits.size();
+ int nmchits = mcTrk.GetNofHits();
if (nmchits == 0) continue;
- double momentum = mtra.p;
- double pt = sqrt(mtra.px * mtra.px + mtra.py * mtra.py);
- double theta = acos(mtra.pz / mtra.p) * 180 / 3.1415;
- double phi = TMath::ATan2(-mtra.py, -mtra.px);
+ double momentum = mcTrk.GetP();
+ double pt = mcTrk.GetPt();
+ double theta = mcTrk.GetTheta() * 180 / 3.1415;
+ double phi = mcTrk.GetPhi();
h_mcp->Fill(momentum);
h_nmchits->Fill(nmchits);
- int nSta = mtra.NStations();
+ int nSta = mcTrk.GetTotNofStationsWithHit();
- CbmL1HitDebugInfo& fh = fvHitDebugInfo[*(mtra.Hits.begin())];
- CbmL1HitDebugInfo& lh = fvHitDebugInfo[*(mtra.Hits.rbegin())];
+ CbmL1HitDebugInfo& fh = fvHitDebugInfo[*(mcTrk.GetHitIndexes().begin())];
+ CbmL1HitDebugInfo& lh = fvHitDebugInfo[*(mcTrk.GetHitIndexes().rbegin())];
h_reg_MCmom->Fill(momentum);
- if (mtra.IsPrimary()) {
+ if (mcTrk.IsPrimary()) {
h_reg_mom_prim->Fill(momentum);
h_reg_phi_prim->Fill(phi);
h_reg_prim_MCmom->Fill(momentum);
@@ -955,15 +864,15 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
}
}
- if (mtra.IsAdditional()) {
+ if (mcTrk.IsAdditional()) {
h_acc_mom_short123s->Fill(momentum);
}
- if (!mtra.IsReconstructable()) continue;
+ if (!mcTrk.IsReconstructable()) continue;
mc_total++;
h_acc_MCmom->Fill(momentum);
- if (mtra.IsPrimary()) {
+ if (mcTrk.IsPrimary()) {
h_acc_mom_prim->Fill(momentum);
h_acc_phi_prim->Fill(phi);
h_acc_prim_MCmom->Fill(momentum);
@@ -986,37 +895,37 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
// vertex distribution of primary and secondary tracks
- h2_vertex->Fill(mtra.z, mtra.y);
- //h3_vertex->Fill(mtra.z, mtra.x, mtra.y);
- if (mtra.mother_ID < 0) { // primary
- h2_prim_vertex->Fill(mtra.z, mtra.y);
- //h3_prim_vertex->Fill(mtra.z, mtra.x, mtra.y);
+ h2_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
+ //h3_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
+ if (mcTrk.GetMotherId() < 0) { // primary
+ h2_prim_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
+ //h3_prim_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
}
else {
- h2_sec_vertex->Fill(mtra.z, mtra.y);
- //h3_sec_vertex->Fill(mtra.z, mtra.x, mtra.y);
+ h2_sec_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
+ //h3_sec_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
}
- int iph = mtra.Hits[0];
+ int iph = mcTrk.GetHitIndexes()[0];
CbmL1HitDebugInfo& ph = fvHitDebugInfo[iph];
h_sec_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
- if (fabs(mtra.pz) > 1.e-8) {
- h_tx->Fill(mtra.px / mtra.pz);
- h_ty->Fill(mtra.py / mtra.pz);
+ if (fabs(mcTrk.GetPz()) > 1.e-8) {
+ h_tx->Fill(mcTrk.GetTx());
+ h_ty->Fill(mcTrk.GetTy());
}
// reconstructed tracks
- bool reco = (mtra.IsReconstructed());
- int nMCHits = mtra.NStations();
+ bool reco = (mcTrk.IsReconstructed());
+ int nMCHits = mcTrk.GetTotNofStationsWithHit();
if (reco) {
p_eff_all_vs_mom->Fill(momentum, 100.0);
p_eff_all_vs_nhits->Fill(nMCHits, 100.0);
p_eff_all_vs_pt->Fill(pt, 100.0);
h_reco_MCmom->Fill(momentum);
- if (mtra.mother_ID < 0) { // primary
+ if (mcTrk.GetMotherId() < 0) { // primary
p_eff_prim_vs_mom->Fill(momentum, 100.0);
p_eff_prim_vs_nhits->Fill(nMCHits, 100.0);
p_eff_prim_vs_pt->Fill(pt, 100.0);
@@ -1026,7 +935,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
p_eff_sec_vs_mom->Fill(momentum, 100.0);
p_eff_sec_vs_nhits->Fill(nMCHits, 100.0);
}
- if (mtra.mother_ID < 0) { // primary
+ if (mcTrk.GetMotherId() < 0) { // primary
h_reco_mom_prim->Fill(momentum);
h_reco_phi_prim->Fill(phi);
h_reco_prim_MCmom->Fill(momentum);
@@ -1053,7 +962,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
p_eff_all_vs_mom->Fill(momentum, 0.0);
p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
p_eff_all_vs_pt->Fill(pt, 0.0);
- if (mtra.mother_ID < 0) { // primary
+ if (mcTrk.GetMotherId() < 0) { // primary
p_eff_prim_vs_mom->Fill(momentum, 0.0);
p_eff_prim_vs_nhits->Fill(nMCHits, 0.0);
p_eff_prim_vs_pt->Fill(pt, 0.0);
@@ -1063,7 +972,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
p_eff_sec_vs_mom->Fill(momentum, 0.0);
p_eff_sec_vs_nhits->Fill(nMCHits, 0.0);
}
- if (mtra.mother_ID < 0) { // primary
+ if (mcTrk.GetMotherId() < 0) { // primary
h_rest_mom_prim->Fill(momentum);
h_rest_phi_prim->Fill(phi);
h_rest_nhits_prim->Fill(nSta);
@@ -1090,10 +999,10 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
h_notfound_station->Fill(ph.iStation);
h_notfound_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
- if (mtra.Points.size() > 0) {
- CbmL1MCPoint& pMC = fvMCPoints[mtra.Points[0]];
- h_notfound_tx->Fill(pMC.px / pMC.pz);
- h_notfound_ty->Fill(pMC.py / pMC.pz);
+ if (mcTrk.GetNofPoints() > 0) {
+ const auto& pMC = fMCData.GetPoint(mcTrk.GetPointIndexes()[0]);
+ h_notfound_tx->Fill(pMC.GetTx());
+ h_notfound_ty->Fill(pMC.GetTy());
}
// CbmL1HitDebugInfo &ph21 = fvHitDebugInfo[mtra.Hits[0]];
@@ -1110,7 +1019,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
}
- if (mtra.isSignal) { // D0
+ if (mcTrk.IsSignal()) { // D0
h_reco_d0_mom->Fill(momentum);
if (reco)
p_eff_d0_vs_mom->Fill(momentum, 100.0);
@@ -1122,7 +1031,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 = fvHitBestPointIndices[ih]; // TODO2: adapt to linking
+ int iMC = fvHitDebugInfo[ih].GetBestMcPointId(); // TODO2: adapt to linking
if (iMC < 0) NFakes++;
}
@@ -1134,6 +1043,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
h_reco_fakeNhit->Fill(fpAlgo->GetInputData().GetNhits() - NFakes, NFakes);
+ // NOTE: Must be called once
h_reg_MCmom->Scale(1.f / NEvents);
h_acc_MCmom->Scale(1.f / NEvents);
h_reco_MCmom->Scale(1.f / NEvents);
@@ -1313,76 +1223,74 @@ void CbmL1::TrackFitPerformance()
do { // first hit
- if (it->GetNMCTracks() < 1) {
+ if (it->GetNofMCTracks() < 1) {
break;
}
-
- const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]);
-
- int imcPoint = fvHitBestPointIndices[it->Hits.front()];
+ const auto& mcTrk = fMCData.GetTrack(it->GetMCTrackIndexes()[0]);
+ int imcPoint = fvHitDebugInfo[it->Hits.front()].GetBestMcPointId();
// extrapolate to the first mc point of the mc track !
- imcPoint = mcTrack.Points[0];
+ imcPoint = mcTrk.GetPointIndexes()[0];
if (imcPoint < 0) break;
- const CbmL1MCPoint& mcP = fvMCPoints[imcPoint];
+ const auto& mcP = fMCData.GetPoint(imcPoint);
TrackParamV tr(*it);
FillFitHistos(tr, mcP, isTimeFitted, h_fit);
- double dx = tr.GetX()[0] - mcP.xIn;
- double dy = tr.GetY()[0] - mcP.yIn;
+ double dx = tr.GetX()[0] - mcP.GetXIn();
+ double dy = tr.GetY()[0] - mcP.GetYIn();
double dca = sqrt(dx * dx + dy * dy);
// make dca distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
- if (mcTrack.ID % 2) dca = -dca;
- double pt = sqrt(mcP.px * mcP.px + mcP.py * mcP.py);
- double dP = (fabs(1. / tr.GetQp()[0]) - mcP.p) / mcP.p;
+ if (mcTrk.GetId() % 2) dca = -dca;
+ double pt = mcP.GetPt();
+ double dP = (fabs(1. / tr.GetQp()[0]) - mcP.GetP()) / mcP.GetP();
- PRes2D->Fill(mcP.p, 100. * dP);
+ PRes2D->Fill(mcP.GetP(), 100. * dP);
- if (mcTrack.IsPrimary()) {
+ if (mcTrk.IsPrimary()) {
h_fit[4]->Fill(100. * dP);
- PRes2DPrim->Fill(mcP.p, 100. * dP);
+ PRes2DPrim->Fill(mcP.GetP(), 100. * dP);
- if (abs(mcTrack.pdg) == 211) {
- pion_res_p_fstDca->Fill(mcP.p, dca * 1.e4);
+ if (abs(mcTrk.GetPdgCode()) == 211) {
+ pion_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
pion_res_pt_fstDca->Fill(pt, dca * 1.e4);
}
- if (abs(mcTrack.pdg) == 321) {
- kaon_res_p_fstDca->Fill(mcP.p, dca * 1.e4);
+ if (abs(mcTrk.GetPdgCode()) == 321) {
+ kaon_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
kaon_res_pt_fstDca->Fill(pt, dca * 1.e4);
}
- if (abs(mcTrack.pdg) == 2212) {
- pton_res_p_fstDca->Fill(mcP.p, dca * 1.e4);
+ if (abs(mcTrk.GetPdgCode()) == 2212) {
+ pton_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
pton_res_pt_fstDca->Fill(pt, dca * 1.e4);
}
}
else {
- PRes2DSec->Fill(mcP.p, 100. * dP);
+ PRes2DSec->Fill(mcP.GetP(), 100. * dP);
}
} while (0);
{ // last hit
- int iMC = fvHitBestPointIndices[it->Hits.back()]; // TODO2: adapt to linking
+ int iMC = fvHitDebugInfo[it->Hits.back()].GetBestMcPointId(); // TODO2: adapt to linking
if (iMC < 0) continue;
- CbmL1MCPoint& mcP = fvMCPoints[iMC];
+ const auto& mcP = fMCData.GetPoint(iMC);
TrackParamV tr(it->TLast);
FillFitHistos(tr, mcP, isTimeFitted, h_fitL);
}
do { // vertex
- if (it->GetNMCTracks() < 1) {
+ if (it->GetNofMCTracks() < 1) {
break;
}
- CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
+ const auto& mc = fMCData.GetTrack(it->GetMatchedMCTrackIndex());
fit.SetTrack(*it);
const TrackParamV& tr = fit.Tr();
@@ -1393,14 +1301,15 @@ void CbmL1::TrackFitPerformance()
{ // extrapolate to vertex
ca::FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
- fit.Extrapolate(mc.z, fld);
+ fit.Extrapolate(mc.GetStartZ(), fld);
// add material
const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
- const int dir = int((mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
- / fabs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]));
+ const int dir = int((mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
+ / fabs(mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0]));
// if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]) > 10.) continue; // can't extrapolate on large distance
- for (int iSta = fSta /*+dir*/; (iSta >= 0) && (iSta < fNStations)
- && (dir * (mc.z - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
+ for (int iSta = fSta /*+dir*/;
+ (iSta >= 0) && (iSta < fNStations)
+ && (dir * (mc.GetStartZ() - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
// cout << iSta << " " << dir << endl;
auto radThick = fpAlgo->GetParameters().GetMaterialThickness(iSta, fit.Tr().GetX(), fit.Tr().GetY());
@@ -1408,7 +1317,7 @@ void CbmL1::TrackFitPerformance()
fit.EnergyLossCorrection(radThick, fvec::One());
}
}
- if (mc.z != tr.GetZ()[0]) continue;
+ if (mc.GetStartZ() != tr.GetZ()[0]) continue;
// static int good = 0;
// static int bad = 0;
@@ -1422,30 +1331,32 @@ void CbmL1::TrackFitPerformance()
// calculate pulls
//h_fitSV[0]->Fill( (mc.x-tr.GetX()[0]) *1.e4);
//h_fitSV[1]->Fill( (mc.y-tr.GetY()[0]) *1.e4);
- h_fitSV[0]->Fill((tr.GetX()[0] - mc.x));
- h_fitSV[1]->Fill((tr.GetY()[0] - mc.y));
- h_fitSV[2]->Fill((tr.GetTx()[0] - mc.px / mc.pz) * 1.e3);
- h_fitSV[3]->Fill((tr.GetTy()[0] - mc.py / mc.pz) * 1.e3);
- h_fitSV[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mc.p - 1.));
+ h_fitSV[0]->Fill((tr.GetX()[0] - mc.GetStartX()));
+ h_fitSV[1]->Fill((tr.GetY()[0] - mc.GetStartY()));
+ h_fitSV[2]->Fill((tr.GetTx()[0] - mc.GetTx()) * 1.e3);
+ h_fitSV[3]->Fill((tr.GetTy()[0] - mc.GetTy()) * 1.e3);
+ h_fitSV[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mc.GetP() - 1.));
if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0) {
- h_fitSV[5]->Fill((tr.GetX()[0] - mc.x) / sqrt(tr.C00()[0]));
+ h_fitSV[5]->Fill((tr.GetX()[0] - mc.GetStartX()) / sqrt(tr.C00()[0]));
}
if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
- h_fitSV[6]->Fill((tr.GetY()[0] - mc.y) / sqrt(tr.C11()[0]));
+ h_fitSV[6]->Fill((tr.GetY()[0] - mc.GetStartY()) / sqrt(tr.C11()[0]));
if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
- h_fitSV[7]->Fill((tr.GetTx()[0] - mc.px / mc.pz) / sqrt(tr.C22()[0]));
+ h_fitSV[7]->Fill((tr.GetTx()[0] - mc.GetTx()) / sqrt(tr.C22()[0]));
if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
- h_fitSV[8]->Fill((tr.GetTy()[0] - mc.py / mc.pz) / sqrt(tr.C33()[0]));
+ h_fitSV[8]->Fill((tr.GetTy()[0] - mc.GetTy()) / sqrt(tr.C33()[0]));
if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
- h_fitSV[9]->Fill((tr.GetQp()[0] - mc.q / mc.p) / sqrt(tr.C44()[0]));
+ h_fitSV[9]->Fill((tr.GetQp()[0] - mc.GetCharge() / mc.GetP()) / sqrt(tr.C44()[0]));
h_fitSV[10]->Fill(tr.GetQp()[0]);
- h_fitSV[11]->Fill(mc.q / mc.p);
+ h_fitSV[11]->Fill(mc.GetCharge() / mc.GetP());
if (isTimeFitted) {
- h_fitSV[12]->Fill(tr.GetTime()[0] - mc.time);
+ h_fitSV[12]->Fill(tr.GetTime()[0] - mc.GetStartT());
if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
- h_fitSV[13]->Fill((tr.GetTime()[0] - mc.time) / sqrt(tr.C55()[0]));
+ h_fitSV[13]->Fill((tr.GetTime()[0] - mc.GetStartT()) / sqrt(tr.C55()[0]));
}
- double dvi = tr.GetVi()[0] - sqrt(1. + mc.mass * mc.mass / mc.p / mc.p) * constants::phys::SpeedOfLightInvD;
+ double dvi =
+ tr.GetVi()[0]
+ - sqrt(1. + mc.GetMass() * mc.GetMass() / mc.GetP() / mc.GetP()) * constants::phys::SpeedOfLightInvD;
h_fitSV[14]->Fill(dvi * constants::phys::SpeedOfLightD);
if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
h_fitSV[15]->Fill(dvi / sqrt(tr.C66()[0]));
@@ -1462,12 +1373,12 @@ void CbmL1::TrackFitPerformance()
// add material
const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
- const int dir = (mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
- / abs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]);
+ const int dir = (mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
+ / abs(mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0]);
// if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta].fZ[0]) > 10.) continue; // can't extrapolate on large distance
for (int iSta = fSta + dir; (iSta >= 0) && (iSta < fNStations)
- && (dir * (mc.z - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
+ && (dir * (mc.GetStartZ() - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
fit.Extrapolate(fpAlgo->GetParameters().GetStation(iSta).fZ, fld);
@@ -1475,55 +1386,57 @@ void CbmL1::TrackFitPerformance()
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec::One());
}
- fit.Extrapolate(mc.z, fld);
+ fit.Extrapolate(mc.GetStartZ(), fld);
}
- if (mc.z != tr.GetZ()[0]) continue;
+ if (mc.GetStartZ() != tr.GetZ()[0]) continue;
- double dx = tr.GetX()[0] - mc.x;
- double dy = tr.GetY()[0] - mc.y;
+ double dx = tr.GetX()[0] - mc.GetStartX();
+ double dy = tr.GetY()[0] - mc.GetStartY();
double dt = sqrt(dx * dx + dy * dy);
// make dt distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
- if (mc.ID % 2) dt = -dt;
+ if (mc.GetId() % 2) dt = -dt;
- double pt = sqrt(mc.px * mc.px + mc.py * mc.py);
+ double pt = mc.GetPt();
- if (abs(mc.pdg) == 211) {
- pion_res_p_vtxDca->Fill(mc.p, dt * 1.e4);
+ if (abs(mc.GetPdgCode()) == 211) {
+ pion_res_p_vtxDca->Fill(mc.GetP(), dt * 1.e4);
pion_res_pt_vtxDca->Fill(pt, dt * 1.e4);
}
- if (abs(mc.pdg) == 321) {
- kaon_res_p_vtxDca->Fill(mc.p, dt * 1.e4);
- kaon_res_pt_vtxDca->Fill(pt, dt * 1.e4);
+ if (abs(mc.GetPdgCode()) == 321) {
+ kaon_res_p_vtxDca->Fill(mc.GetPdgCode(), dt * 1.e4);
+ kaon_res_pt_vtxDca->Fill(mc.GetPt(), dt * 1.e4);
}
- if (abs(mc.pdg) == 2212) {
- pton_res_p_vtxDca->Fill(mc.p, dt * 1.e4);
+ if (abs(mc.GetPdgCode()) == 2212) {
+ pton_res_p_vtxDca->Fill(mc.GetP(), dt * 1.e4);
pton_res_pt_vtxDca->Fill(pt, dt * 1.e4);
}
// calculate pulls
- h_fitPV[0]->Fill((mc.x - tr.GetX()[0]));
- h_fitPV[1]->Fill((mc.y - tr.GetY()[0]));
- h_fitPV[2]->Fill((mc.px / mc.pz - tr.GetTx()[0]) * 1.e3);
- h_fitPV[3]->Fill((mc.py / mc.pz - tr.GetTy()[0]) * 1.e3);
- h_fitPV[4]->Fill(100. * (fabs(1 / tr.GetQp()[0]) / mc.p - 1.));
+ h_fitPV[0]->Fill((mc.GetStartX() - tr.GetX()[0]));
+ h_fitPV[1]->Fill((mc.GetStartY() - tr.GetY()[0]));
+ h_fitPV[2]->Fill((mc.GetTx() - tr.GetTx()[0]) * 1.e3);
+ h_fitPV[3]->Fill((mc.GetTy() - tr.GetTy()[0]) * 1.e3);
+ h_fitPV[4]->Fill(100. * (fabs(1 / tr.GetQp()[0]) / mc.GetP() - 1.));
if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0)
- h_fitPV[5]->Fill((mc.x - tr.GetX()[0]) / sqrt(tr.C00()[0]));
+ h_fitPV[5]->Fill((mc.GetStartX() - tr.GetX()[0]) / sqrt(tr.C00()[0]));
if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
- h_fitPV[6]->Fill((mc.y - tr.GetY()[0]) / sqrt(tr.C11()[0]));
+ h_fitPV[6]->Fill((mc.GetStartY() - tr.GetY()[0]) / sqrt(tr.C11()[0]));
if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
- h_fitPV[7]->Fill((mc.px / mc.pz - tr.GetTx()[0]) / sqrt(tr.C22()[0]));
+ h_fitPV[7]->Fill((mc.GetTx() - tr.GetTx()[0]) / sqrt(tr.C22()[0]));
if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
- h_fitPV[8]->Fill((mc.py / mc.pz - tr.GetTy()[0]) / sqrt(tr.C33()[0]));
+ h_fitPV[8]->Fill((mc.GetTy() - tr.GetTy()[0]) / sqrt(tr.C33()[0]));
if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
- h_fitPV[9]->Fill((mc.q / mc.p - tr.GetQp()[0]) / sqrt(tr.C44()[0]));
+ h_fitPV[9]->Fill((mc.GetCharge() / mc.GetP() - tr.GetQp()[0]) / sqrt(tr.C44()[0]));
h_fitPV[10]->Fill(tr.GetQp()[0]);
- h_fitPV[11]->Fill(mc.q / mc.p);
+ h_fitPV[11]->Fill(mc.GetCharge() / mc.GetP());
if (isTimeFitted) {
- h_fitPV[12]->Fill(tr.GetTime()[0] - mc.time);
+ h_fitPV[12]->Fill(tr.GetTime()[0] - mc.GetStartT());
if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
- h_fitPV[13]->Fill((tr.GetTime()[0] - mc.time) / sqrt(tr.C55()[0]));
+ h_fitPV[13]->Fill((tr.GetTime()[0] - mc.GetStartT()) / sqrt(tr.C55()[0]));
}
- double dvi = tr.GetVi()[0] - sqrt(1. + mc.mass * mc.mass / mc.p / mc.p) * constants::phys::SpeedOfLightInvD;
+ double dvi =
+ tr.GetVi()[0]
+ - sqrt(1. + mc.GetMass() * mc.GetMass() / mc.GetP() / mc.GetP()) * constants::phys::SpeedOfLightInvD;
h_fitPV[14]->Fill(dvi * constants::phys::SpeedOfLightD);
if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
h_fitPV[15]->Fill(dvi / sqrt(tr.C66()[0]));
@@ -1595,7 +1508,7 @@ void CbmL1::TrackFitPerformance()
} // void CbmL1::TrackFitPerformance()
-void CbmL1::FillFitHistos(TrackParamV& track, const CbmL1MCPoint& mcP, bool isTimeFitted, TH1F* h[])
+void CbmL1::FillFitHistos(TrackParamV& track, const cbm::ca::tools::MCPoint& mcP, bool isTimeFitted, TH1F* h[])
{
ca::TrackFit fit;
fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
@@ -1605,33 +1518,36 @@ void CbmL1::FillFitHistos(TrackParamV& track, const CbmL1MCPoint& mcP, bool isTi
fit.SetTrack(track);
ca::FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
- fit.Extrapolate(mcP.zOut, fld);
+ fit.Extrapolate(mcP.GetZOut(), fld);
track = fit.Tr();
const TrackParamV& tr = track;
- h[0]->Fill((tr.GetX()[0] - mcP.xOut) * 1.e4);
- h[1]->Fill((tr.GetY()[0] - mcP.yOut) * 1.e4);
- h[2]->Fill((tr.GetTx()[0] - mcP.pxOut / mcP.pzOut) * 1.e3);
- h[3]->Fill((tr.GetTy()[0] - mcP.pyOut / mcP.pzOut) * 1.e3);
- h[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mcP.p - 1.));
+ h[0]->Fill((tr.GetX()[0] - mcP.GetXOut()) * 1.e4);
+ h[1]->Fill((tr.GetY()[0] - mcP.GetYOut()) * 1.e4);
+ h[2]->Fill((tr.GetTx()[0] - mcP.GetTxOut()) * 1.e3);
+ h[3]->Fill((tr.GetTy()[0] - mcP.GetTyOut()) * 1.e3);
+ h[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mcP.GetP() - 1.));
- if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0) h[5]->Fill((tr.GetX()[0] - mcP.xOut) / sqrt(tr.C00()[0]));
- if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0) h[6]->Fill((tr.GetY()[0] - mcP.yOut) / sqrt(tr.C11()[0]));
+ if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0)
+ h[5]->Fill((tr.GetX()[0] - mcP.GetXOut()) / sqrt(tr.C00()[0]));
+ if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
+ h[6]->Fill((tr.GetY()[0] - mcP.GetYOut()) / sqrt(tr.C11()[0]));
if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
- h[7]->Fill((tr.GetTx()[0] - mcP.pxOut / mcP.pzOut) / sqrt(tr.C22()[0]));
+ h[7]->Fill((tr.GetTx()[0] - mcP.GetTxOut()) / sqrt(tr.C22()[0]));
if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
- h[8]->Fill((tr.GetTy()[0] - mcP.pyOut / mcP.pzOut) / sqrt(tr.C33()[0]));
+ h[8]->Fill((tr.GetTy()[0] - mcP.GetTyOut()) / sqrt(tr.C33()[0]));
if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
- h[9]->Fill((tr.GetQp()[0] - mcP.q / mcP.p) / sqrt(tr.C44()[0]));
+ h[9]->Fill((tr.GetQp()[0] - mcP.GetQp()) / sqrt(tr.C44()[0]));
h[10]->Fill(tr.GetQp()[0]);
- h[11]->Fill(mcP.q / mcP.p);
+ h[11]->Fill(mcP.GetQp());
if (isTimeFitted) {
- h[12]->Fill(tr.GetTime()[0] - mcP.time);
+ h[12]->Fill(tr.GetTime()[0] - mcP.GetTime());
if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
- h[13]->Fill((tr.GetTime()[0] - mcP.time) / sqrt(tr.C55()[0]));
+ h[13]->Fill((tr.GetTime()[0] - mcP.GetTime()) / sqrt(tr.C55()[0]));
}
- double dvi = tr.GetVi()[0] - sqrt(1. + mcP.mass * mcP.mass / mcP.p / mcP.p) * constants::phys::SpeedOfLightInvD;
+
+ double dvi = tr.GetVi()[0] - mcP.GetInvSpeed();
h[14]->Fill(dvi * constants::phys::SpeedOfLightD);
if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
h[15]->Fill(dvi / sqrt(tr.C66()[0]));
@@ -2394,17 +2310,18 @@ void CbmL1::DumpMCTripletsToTree()
float z; ///< z-component of point position [cm]
};
- for (const auto& tr : fvMCTracks) {
+ for (const auto& tr : fMCData.GetTrackContainer()) {
// Use only reconstructable tracks
if (!tr.IsReconstructable()) {
continue;
}
std::vector<ReducedMcPoint> vPoints;
- vPoints.reserve(tr.Points.size());
- for (unsigned int iP = 0; iP < tr.Points.size(); ++iP) {
- const auto& point = fvMCPoints[tr.Points[iP]];
- vPoints.emplace_back(ReducedMcPoint{point.iStation, float(point.x), float(point.y), float(point.z)});
+ vPoints.reserve(tr.GetNofPoints());
+ for (unsigned int iP : tr.GetPointIndexes()) {
+ const auto& point = fMCData.GetPoint(iP);
+ vPoints.emplace_back(
+ ReducedMcPoint{point.GetStationId(), float(point.GetX()), float(point.GetY()), float(point.GetZ())});
}
std::sort(vPoints.begin(), vPoints.end(),
@@ -2415,12 +2332,12 @@ void CbmL1::DumpMCTripletsToTree()
// TODO: SZh 20.10.2022 Add cases for jump iterations
if (vPoints[i + 1].s == vPoints[i].s + 1 && vPoints[i + 2].s == vPoints[i].s + 2) {
// Fill MC-triplets tree
- brMotherId = tr.mother_ID;
- brPdg = tr.pdg;
- brProcId = tr.process_ID;
- brP = tr.p;
- brQ = tr.q;
- brVertexZ = tr.z;
+ brMotherId = tr.GetMotherId();
+ brPdg = tr.GetPdgCode();
+ brProcId = tr.GetProcessId();
+ brP = tr.GetP();
+ brQ = tr.GetCharge();
+ brVertexZ = tr.GetStartZ();
brStation = vPoints[i].s;
brX0 = vPoints[i].x;
brY0 = vPoints[i].y;
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
deleted file mode 100644
index 930ff3cb3ae4658129e5eb40a2e31be41c8aeb2b..0000000000000000000000000000000000000000
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ /dev/null
@@ -1,1130 +0,0 @@
-/* Copyright (C) 2006-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Ivan Kisel, Sergey Gorbunov [committer], Irina Rostovtseva, Valentina Akishina, Maksym Zyzak */
-
-/*
- *====================================================================
- *
- * CBM Level 1 Reconstruction
- *
- * Authors: I.Kisel, S.Gorbunov, I. Rostovtseva (2016)
- *
- *
- * e-mail : ikisel@kip.uni-heidelberg.de
- *
- *====================================================================
- *
- * Read Event information to L1
- *
- *====================================================================
- */
-
-#include "CaFramework.h"
-#include "CaToolsDebugger.h"
-#include "CbmEvent.h"
-#include "CbmL1.h"
-#include "CbmMCDataObject.h"
-#include "CbmMatch.h"
-#include "CbmMuchPixelHit.h"
-#include "CbmMuchPoint.h"
-#include "CbmMuchTrackingInterface.h"
-#include "CbmMvdTrackingInterface.h"
-#include "CbmStsAddress.h"
-#include "CbmStsCluster.h"
-#include "CbmStsDigi.h"
-#include "CbmStsHit.h"
-#include "CbmStsPoint.h"
-#include "CbmStsSetup.h"
-#include "CbmStsTrackingInterface.h"
-#include "CbmTofHit.h"
-#include "CbmTofPoint.h"
-#include "CbmTofTrackingInterface.h"
-#include "CbmTrdHit.h"
-#include "CbmTrdPoint.h"
-#include "CbmTrdTrackingInterface.h"
-#include "FairMCEventHeader.h"
-#include "TDatabasePDG.h"
-#include "TRandom.h"
-
-#include <iostream>
-
-using std::cout;
-using std::endl;
-using std::find;
-using namespace cbm::algo;
-using cbm::ca::tools::Debugger;
-
-static bool compareMcPointZ(const int& a, const int& b)
-{
- // return (CbmL1::Instance()->fvMCPoints[a].z < CbmL1::Instance()->fvMCPoints[b].z);
- const CbmL1* l1 = CbmL1::Instance();
- return l1->Get_Z_vMCPoint(a) < l1->Get_Z_vMCPoint(b);
-}
-
-/// Local structure for a hit, containing both measured and MC information.
-/// The structure is used to sort hits before writing them into L1 input arrays
-///
-struct TmpHit {
- int iStripF; ///< index of front strip
- int iStripB; ///< index of back strip
- int iStation; ///< index of station
- int64_t fDataStream; ///< global index of detector module (4 left bits -- station number, other - address)
- int ExtIndex; ///< index of hit in the external TClonesArray array (NOTE: negative for MVD)
- double x; ///< position of hit along x axis [cm]
- double y; ///< position of hit along y axis [cm]
- double z; ///< position of hit along z axis [cm]
- 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]
- 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;
-
- static constexpr bool kVerboseAccept = true;
- /// @brief Tests quality of the hit (eg., if one of its quantity is nan)
- /// @return true Hit is passed and send to the CaCore
- /// @return false Hit is discarded due to corrupted quantity values
- bool IfAccept() const
- {
- if constexpr (kVerboseAccept) {
- std::stringstream msg;
- if (std::isnan(x) || std::isinf(x)) {
- msg << ", x = " << x;
- }
- if (std::isnan(y) || std::isinf(y)) {
- msg << ", y = " << y;
- }
- if (std::isnan(z) || std::isinf(z)) {
- msg << ", z = " << z;
- }
- if (std::isnan(time) || std::isinf(time)) {
- msg << ", t = " << time;
- }
- if (std::isnan(dx2) || std::isinf(dx2)) {
- msg << ", dx2 = " << dx2;
- }
- if (std::isnan(dy2) || std::isinf(dy2)) {
- msg << ", dy2 = " << dy2;
- }
- if (std::isnan(dxy) || std::isinf(dxy)) {
- msg << ", dxy = " << dxy;
- }
- if (std::isnan(dt2) || std::isinf(dt2)) {
- msg << ", dt2 = " << dt2;
- }
- const auto& sMsg = msg.str();
- if (sMsg.size()) {
- LOG(warn) << "CbmL1: Discarding hit " << ExtIndex << ": det = " << Det
- << ", addr = " << static_cast<int32_t>(fDataStream) << ", " << sMsg;
- return false;
- }
- }
- else {
- bool res = true;
- res = res && !(std::isnan(x) || std::isinf(x));
- res = res && !(std::isnan(y) || std::isinf(y));
- res = res && !(std::isnan(z) || std::isinf(z));
- res = res && !(std::isnan(time) || std::isinf(time));
- res = res && !(std::isnan(dx2) || std::isinf(dx2));
- res = res && !(std::isnan(dy2) || std::isinf(dy2));
- res = res && !(std::isnan(dt2) || std::isinf(dt2));
- res = res && !(std::isnan(dxy) || std::isinf(dxy));
- return res;
- }
- return true;
- }
-};
-
-
-// **************************
-// ** Event reader **
-// **************************
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void CbmL1::ReadEvent(CbmEvent* event)
-{
- static int nCalls = 0;
-
- if (fVerbose >= 10) cout << "ReadEvent: start." << endl;
-
- // clear arrays for next event
- 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();
-
-
- if (fVerbose >= 10) cout << "ReadEvent: clear is done." << endl;
-
- ca::Vector<TmpHit> tmpHits("CbmL1ReadEvent::tmpHits");
-
- { // reserve enough space for hits
- int nHitsTotal = 0;
- if (event) { // be careful: GetNofData() can return -1
- if (fUseMVD && fpMvdHits) {
- nHitsTotal += event->GetNofData(ECbmDataType::kMvdHit);
- }
- if (fUseSTS && fpStsHits) {
- nHitsTotal += event->GetNofData(ECbmDataType::kStsHit);
- }
- if (fUseMUCH && fpMuchPixelHits) {
- nHitsTotal += event->GetNofData(ECbmDataType::kMuchPixelHit);
- }
- if (fUseTRD && fpTrdHits) {
- nHitsTotal += event->GetNofData(ECbmDataType::kTrdHit);
- }
- if (fUseTOF && fpTofHits) {
- nHitsTotal += event->GetNofData(ECbmDataType::kTofHit);
- }
- }
- else {
- if (fUseMVD && fpMvdHits) {
- nHitsTotal += fpMvdHits->GetEntriesFast();
- }
- if (fUseSTS && fpStsHits) {
- nHitsTotal += fpStsHits->GetEntriesFast();
- }
- if (fUseMUCH && fpMuchPixelHits) {
- nHitsTotal += fpMuchPixelHits->GetEntriesFast();
- }
- if (fUseTRD && fpTrdHits) {
- nHitsTotal += fpTrdHits->GetEntriesFast();
- }
- if (fUseTOF && fpTofHits) {
- nHitsTotal += fpTofHits->GetEntriesFast();
- }
- }
- tmpHits.reserve(nHitsTotal);
- }
-
- fNpointsMvd = 0;
- fNpointsSts = 0;
- fNpointsTrd = 0;
- fNpointsMuch = 0;
- fNpointsTof = 0;
-
- // get MVD hits
- int nMvdHits = 0;
- int nStsHits = 0;
- int nMuchHits = 0;
- int nTrdHits = 0;
- int nTofHits = 0;
-
- /*
- * MC hits and tracks gathering: START
- *
- * If performance is studied, i.e. fPerformanse > 0, MC-information is used for
- * efficiencies calculation. Otherwise this
- */
-
- if (fPerformance) {
- // Fill fvMCTracks and fmMCTracksLinksMap
- Fill_vMCTracks();
-
- // Fill fvMCPoints and fmMCPointsLinksMap
- fvMCPoints.clear();
- fvMCPoints.reserve(5 * fvMCTracks.size() * fNStations);
-
- fNpointsMvdAll = 0;
- fNpointsStsAll = 0;
- fNpointsMuchAll = 0;
- fNpointsTrdAll = 0;
- fNpointsTofAll = 0;
-
- // count total amount of mc points to set up the global indexing
- for (DFSET::iterator set_it = fvSelectedMcEvents.begin(); set_it != fvSelectedMcEvents.end(); ++set_it) {
- Int_t iFile = set_it->first;
- Int_t iEvent = set_it->second;
- if (fUseMVD) fNpointsMvdAll += fpMvdPoints->Size(iFile, iEvent);
- if (fUseSTS) fNpointsStsAll += fpStsPoints->Size(iFile, iEvent);
- if (fUseMUCH) fNpointsMuchAll += fpMuchPoints->Size(iFile, iEvent);
- if (fUseTRD) fNpointsTrdAll += fpTrdPoints->Size(iFile, iEvent);
- if (fUseTOF) fNpointsTofAll += fpTofPoints->Size(iFile, iEvent);
- }
-
- for (DFSET::iterator set_it = fvSelectedMcEvents.begin(); set_it != fvSelectedMcEvents.end(); ++set_it) {
- Int_t iFile = set_it->first;
- Int_t iEvent = set_it->second;
-
- if (fUseMVD && fpMvdPoints) {
- Int_t nMC = fpMvdPoints->Size(iFile, iEvent);
- for (Int_t iMC = 0; iMC < nMC; iMC++) {
- ReadMCPoint(ca::EDetectorID::kMvd, iFile, iEvent, iMC);
- }
- }
-
- if (fUseSTS && fpStsPoints) {
- Int_t nMC = fpStsPoints->Size(iFile, iEvent);
- for (Int_t iMC = 0; iMC < nMC; iMC++) {
- ReadMCPoint(ca::EDetectorID::kSts, iFile, iEvent, iMC);
- }
- }
-
- if (fUseMUCH && fpMuchPoints) {
- Int_t nMC = fpMuchPoints->Size(iFile, iEvent);
- for (Int_t iMC = 0; iMC < nMC; iMC++) {
- ReadMCPoint(ca::EDetectorID::kMuch, iFile, iEvent, iMC);
- }
- }
-
- if (fUseTRD && fpTrdPoints) {
- for (Int_t iMC = 0; iMC < fpTrdPoints->Size(iFile, iEvent); iMC++) {
- ReadMCPoint(ca::EDetectorID::kTrd, iFile, iEvent, iMC);
- }
- }
-
-
- if (fUseTOF && fpTofPoints) {
-
- // FIXME: Use mask from CbmTofAddress (but test before!!)
- constexpr int kNofBitsRpcAddress = 11;
-
- // Fill map of flags: the external index of the matched point vs. iTrExt and the RPC address
- int nPoints = fpTofPoints->Size(iFile, iEvent);
- std::map<std::pair<int, int>, int> mMatchedPointId; // map (iTr, addr) -> is matched
- for (int iH = 0; iH < fpTofHitMatches->GetEntriesFast(); ++iH) {
- auto* pHitMatch = dynamic_cast<CbmMatch*>(fpTofHitMatches->At(iH));
- LOG_IF(fatal, !pHitMatch) << "CA MC Module: hit match was not found for TOF hit " << iH;
- double bestWeight = 0;
- for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) {
- const auto& link = pHitMatch->GetLink(iLink);
- if (link.GetFile() != iFile || link.GetEntry() != iEvent) {
- continue;
- }
- int iPext = link.GetIndex();
- if (iPext < 0) {
- continue;
- }
- auto* pExtPoint = static_cast<CbmTofPoint*>(fpTofPoints->Get(link));
- int trkId = pExtPoint->GetTrackID();
- int rpcAddr = pExtPoint->GetDetectorID() << kNofBitsRpcAddress; // FIXME:
- auto key = std::make_pair(trkId, rpcAddr);
- auto prev = mMatchedPointId.find(key);
- if (prev == mMatchedPointId.end()) {
- bestWeight = link.GetWeight();
- mMatchedPointId[key] = iPext;
- }
- else { // If we find two links for the same interaction, we select the one with the largest weight
- if (bestWeight < link.GetWeight()) {
- mMatchedPointId[key] = iPext;
- }
- }
- }
- } // iH
-
- // Select one point for the given track ID and RPC address. If at least one of the points for (trackID, RPC address)
- // produced a hit, the point from the best link is selected. Otherwise the point closest to the RPC z-center is selected
- {
- int iPointSelected = -1;
- int iTmcCurr = -1;
- int rpcAddrCurr = -1;
- bool bTrkHasHits = false;
- double zDist = std::numeric_limits<double>::max();
- double zCell = std::numeric_limits<double>::signaling_NaN();
- for (int iP = 0; iP < nPoints; ++iP) {
- auto* pExtPoint = dynamic_cast<CbmTofPoint*>(fpTofPoints->Get(iFile, iEvent, iP));
- LOG_IF(fatal, !pExtPoint) << "NO POINT: TOF: file, event, index = " << iFile << ' ' << iEvent << ' ' << iP;
-
- int iTmc = pExtPoint->GetTrackID();
- int rpcAddr = pExtPoint->GetDetectorID() << kNofBitsRpcAddress;
- // New point
- if (rpcAddrCurr != rpcAddr || iTmcCurr != iTmc) {
- if (iPointSelected != -1) {
- ReadMCPoint(ca::EDetectorID::kTof, iFile, iEvent, iPointSelected);
- }
- iTmcCurr = iTmc;
- rpcAddrCurr = rpcAddr;
- auto key = std::make_pair(iTmc, rpcAddr);
- auto found = mMatchedPointId.find(key);
- bTrkHasHits = found != mMatchedPointId.end();
- if (bTrkHasHits) {
- iPointSelected = found->second;
- }
- else {
- // First iteration
- zCell = CbmTofTrackingInterface::Instance()->GetZrefModule(pExtPoint->GetDetectorID());
- zDist = std::fabs(zCell - pExtPoint->GetZ());
- iPointSelected = iP;
- }
- }
- else {
- if (!bTrkHasHits) {
- auto newDist = std::fabs(pExtPoint->GetZ() - zCell);
- if (newDist < zDist) {
- zDist = newDist;
- iPointSelected = iP;
- }
- }
- }
- }
- // Add the last point
- if (iPointSelected != -1) {
- ReadMCPoint(ca::EDetectorID::kTof, iFile, iEvent, iPointSelected);
- }
- } // iP
- }
- } //time_slice
-
- for (unsigned int iTr = 0; iTr < fvMCTracks.size(); iTr++) {
-
- sort(fvMCTracks[iTr].Points.begin(), fvMCTracks[iTr].Points.end(), compareMcPointZ);
-
- if (fvMCTracks[iTr].mother_ID >= 0) {
- auto iFile = fvMCTracks[iTr].iFile;
- auto iEvent = fvMCTracks[iTr].iEvent;
- auto iMother = fvMCTracks[iTr].mother_ID;
- auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(iMother, iEvent, iFile));
- fvMCTracks[iTr].mother_ID = ((itTrack != fmMCTracksLinksMap.cend()) ? itTrack->second : -2);
- }
- } //iTr
- if (fVerbose >= 10) cout << "Points in fvMCTracks are sorted." << endl;
- } //fPerformance
-
- /*
- * 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: independently from particular detector design) and then pushed to the tmpHit vector. In the performance study mode
- * matching with MC points is done
- */
-
- // TODO: Refactor this code: MC information and matches should be collected separately in a separate class (S.Zharko)
-
- int NStrips = 0;
-
- //
- // get MVD hits
- if (fUseMVD && fpMvdHits) {
-
- int firstDetStrip = NStrips;
-
- // MVD hits are not yet fully integrated to the event builder
- // catch the situation when MVD hits are present in data, but there is no corresp. branch in the event
- // read hits directly from the time slice
-
- Int_t nEntries = (event ? event->GetNofData(ECbmDataType::kMvdHit) : fpMvdHits->GetEntriesFast());
-
- for (int j = 0; j < nEntries; j++) {
-
- Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kMvdHit, j) : j);
-
- TmpHit th;
- CbmMvdHit* h = dynamic_cast<CbmMvdHit*>(fpMvdHits->At(hitIndex));
- {
- th.ExtIndex = hitIndex;
- th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
- CbmMvdTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kMvd);
-
- if (th.iStation < 0) continue;
-
- th.iStripF = firstDetStrip + j;
- th.iStripB = th.iStripF;
- if (NStrips <= th.iStripF) {
- NStrips = th.iStripF + 1;
- }
-
- TVector3 pos, err;
- h->Position(pos);
- h->PositionError(err);
-
- th.dx2 = h->GetDx() * h->GetDx() + fpAlgo->GetParameters().GetMisalignmentXsq(ca::EDetectorID::kMvd);
- th.dy2 = h->GetDy() * h->GetDy() + fpAlgo->GetParameters().GetMisalignmentYsq(ca::EDetectorID::kMvd);
- th.dxy = h->GetDxy();
-
- th.x = pos.X();
- th.y = pos.Y();
- th.z = pos.Z();
-
- // Get time
- th.time = h->GetTime(); // currently ignored by the tracking
- th.dt2 =
- h->GetTimeError() * h->GetTimeError()
- + fpAlgo->GetParameters().GetMisalignmentTsq(ca::EDetectorID::kMvd); // currently ignored by the tracking
-
- std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmMvdTrackingInterface::Instance()->GetHitRanges(*h);
- }
- th.Det = 0;
- if (fPerformance) {
- std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<ca::EDetectorID::kMvd>(hitIndex);
- }
- th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
- //if( th.iMC >=0 ) // DEBUG !!!!
- if (th.IfAccept()) {
- tmpHits.push_back(th);
- nMvdHits++;
- }
- else {
- //fMonitor.IncrementCounter(EMonitorKey::kDiscardedHitMvd);
- }
- } // for j
- } // if fpMvdHits
-
- //
- // Get STS hits
- //
- if (fUseSTS && fpStsHits) {
-
- Int_t nEntSts = (event ? event->GetNofData(ECbmDataType::kStsHit) : fpStsHits->GetEntriesFast());
-
- int firstDetStrip = NStrips;
-
- for (Int_t j = 0; j < nEntSts; j++) {
-
- Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kStsHit, j) : j);
-
- // ***********************************
- // ** Fill the temporary hit object **
- // ***********************************
-
- TmpHit th;
- CbmStsHit* h = dynamic_cast<CbmStsHit*>(fpStsHits->At(hitIndex));
-
- // Fill reconstructed information
- {
- th.ExtIndex = hitIndex;
- th.Det = 1;
-
- th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
- CbmStsTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kSts);
-
- if (th.iStation < 0) continue;
-
- th.iStripF = firstDetStrip + h->GetFrontClusterId();
- th.iStripB = firstDetStrip + h->GetBackClusterId();
-
- if (NStrips <= th.iStripF) {
- NStrips = th.iStripF + 1;
- }
- if (NStrips <= th.iStripB) {
- NStrips = th.iStripB + 1;
- }
-
- //Get time
- th.time = h->GetTime();
- th.dt2 =
- h->GetTimeError() * h->GetTimeError() + fpAlgo->GetParameters().GetMisalignmentTsq(ca::EDetectorID::kSts);
-
- TVector3 pos, err;
- h->Position(pos);
- h->PositionError(err);
-
- th.x = pos.X();
- th.y = pos.Y();
- th.z = pos.Z();
-
- th.dx2 = h->GetDx() * h->GetDx() + fpAlgo->GetParameters().GetMisalignmentXsq(ca::EDetectorID::kSts);
- th.dy2 = h->GetDy() * h->GetDy() + fpAlgo->GetParameters().GetMisalignmentYsq(ca::EDetectorID::kSts);
- th.dxy = h->GetDxy();
-
- std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmStsTrackingInterface::Instance()->GetHitRanges(*h);
- }
-
- if (fPerformance) {
- std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<ca::EDetectorID::kSts>(hitIndex);
- }
-
- th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
- if (th.IfAccept()) {
- tmpHits.push_back(th);
- nStsHits++;
- }
- else {
- //fMonitor.IncrementCounter(EMonitorKey::kDiscardedHitSts);
- }
- } // for j
- } // if fpStsHits
-
-
- //
- // Get MuCh hits
- //
-
- if (fUseMUCH && fpMuchPixelHits) {
-
- Int_t nEnt = (event ? event->GetNofData(ECbmDataType::kMuchPixelHit) : fpMuchPixelHits->GetEntriesFast());
-
- int firstDetStrip = NStrips;
-
- for (int j = 0; j < nEnt; j++) {
- TmpHit th;
- Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kMuchPixelHit, j) : j);
- CbmMuchPixelHit* h = static_cast<CbmMuchPixelHit*>(fpMuchPixelHits->At(hitIndex));
- {
- th.ExtIndex = hitIndex;
- th.Det = 2;
-
- th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
- CbmMuchTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kMuch);
-
- if (th.iStation < 0) continue;
-
- //Get time
- th.time = h->GetTime();
- th.dt2 =
- h->GetTimeError() * h->GetTimeError() + fpAlgo->GetParameters().GetMisalignmentTsq(ca::EDetectorID::kMuch);
-
- // th.iSector = 0;
- th.iStripF = firstDetStrip + hitIndex;
- th.iStripB = th.iStripF;
- if (NStrips <= th.iStripF) {
- NStrips = th.iStripF + 1;
- }
-
- th.x = h->GetX();
- th.y = h->GetY();
- th.z = h->GetZ();
-
- th.dx2 = h->GetDx() * h->GetDx() + fpAlgo->GetParameters().GetMisalignmentXsq(ca::EDetectorID::kMuch);
- th.dy2 = h->GetDy() * h->GetDy() + fpAlgo->GetParameters().GetMisalignmentYsq(ca::EDetectorID::kMuch);
- th.dxy = 0; /// FIXME: ???
-
- std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmMuchTrackingInterface::Instance()->GetHitRanges(*h);
- }
- if (fPerformance) {
- std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<ca::EDetectorID::kMuch>(hitIndex);
- }
- th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
- if (th.IfAccept()) {
- tmpHits.push_back(th);
- nMuchHits++;
- }
- else {
- //fMonitor.IncrementCounter(EMonitorKey::kDiscardedHitMuch);
- }
-
- } // for j
- } // if listMuchHits
-
-
- //
- // Get TRD hits
- //
-
- if (fUseTRD && fpTrdHits) {
-
- assert(fpTrdHits);
-
- std::vector<char> mcUsed(fNpointsTrd, 0);
-
- Int_t nEntTrd = (event ? event->GetNofData(ECbmDataType::kTrdHit) : fpTrdHits->GetEntriesFast());
-
- for (int iHit = 0; iHit < nEntTrd; iHit++) {
- TmpHit th;
-
- Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kTrdHit, iHit) : iHit);
- CbmTrdHit* h = dynamic_cast<CbmTrdHit*>(fpTrdHits->At(hitIndex));
-
- if ((ca::Framework::TrackingMode::kGlobal == fTrackingMode) && (int) h->GetClassType() != 1) {
- // SGtrd2d!! skip TRD-1D hit
- continue;
- }
-
- th.ExtIndex = hitIndex;
- th.Det = 3;
-
- th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
- CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kTrd);
-
- if (th.iStation < 0) continue;
-
- // if (h->GetPlaneId()==0) continue;
- // if (h->GetPlaneId()==1) continue;
- // if (h->GetPlaneId()==2) continue;
- // if (h->GetPlaneId()==3) continue;
-
- th.time = h->GetTime();
- th.dt2 =
- h->GetTimeError() * h->GetTimeError() + fpAlgo->GetParameters().GetMisalignmentTsq(ca::EDetectorID::kTrd);
-
- // th.iSector = 0;
- th.iStripF = NStrips;
- th.iStripB = th.iStripF; //TrdHitsOnStationIndex[num+1][k];
- NStrips++;
-
- TVector3 pos, err;
- h->Position(pos);
- h->PositionError(err);
-
- th.x = pos.X();
- th.y = pos.Y();
- th.z = pos.Z();
-
- th.dx2 = h->GetDx() * h->GetDx() + fpAlgo->GetParameters().GetMisalignmentXsq(ca::EDetectorID::kTrd);
- th.dy2 = h->GetDy() * h->GetDy() + fpAlgo->GetParameters().GetMisalignmentYsq(ca::EDetectorID::kTrd);
- th.dxy = 0;
-
- std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmTrdTrackingInterface::Instance()->GetHitRanges(*h);
-
- if (fPerformance) {
- std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<ca::EDetectorID::kTrd>(hitIndex);
- }
-
- th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
- if (th.IfAccept()) {
- tmpHits.push_back(th);
- nTrdHits++;
- }
- else {
- //fMonitor.IncrementCounter(EMonitorKey::kDiscardedHitTrd);
- }
- } // for fpTrdHits
- } // read TRD hits
-
- //
- // Get ToF hits
- //
-
- if (fUseTOF && fpTofHits) {
- int firstDetStrip = NStrips;
-
- Int_t nEntTof = (event ? event->GetNofData(ECbmDataType::kTofHit) : fpTofHits->GetEntriesFast());
-
- for (int j = 0; j < nEntTof; j++) {
-
- Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kTofHit, j) : j);
-
- TmpHit th;
-
- CbmTofHit* h = dynamic_cast<CbmTofHit*>(fpTofHits->At(hitIndex));
-
-
- th.ExtIndex = hitIndex;
- th.Det = 4;
-
- if (0x00202806 == h->GetAddress() || 0x00002806 == h->GetAddress()) continue; // TODO: Why? (S.Zharko)
-
- if (5 == CbmTofAddress::GetSmType(h->GetAddress())) continue; // Skip Bmon hits from TOF
-
- th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
- CbmTofTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kTof);
-
- if (th.iStation < 0) continue;
-
- th.time = h->GetTime();
- th.dt2 =
- h->GetTimeError() * h->GetTimeError() + fpAlgo->GetParameters().GetMisalignmentTsq(ca::EDetectorID::kTof);
-
- th.dx2 = h->GetDx() * h->GetDx() + fpAlgo->GetParameters().GetMisalignmentXsq(ca::EDetectorID::kTof);
- th.dy2 = h->GetDy() * h->GetDy() + fpAlgo->GetParameters().GetMisalignmentYsq(ca::EDetectorID::kTof);
- th.dxy = h->GetDxy();
-
- // th.iSector = 0;
- th.iStripF = firstDetStrip + j;
- th.iStripB = th.iStripF;
- if (NStrips <= th.iStripF) {
- NStrips = th.iStripF + 1;
- }
-
- TVector3 pos, err;
- h->Position(pos);
- h->PositionError(err);
-
- th.x = pos.X();
- th.y = pos.Y();
- th.z = pos.Z();
-
- std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmTofTrackingInterface::Instance()->GetHitRanges(*h);
-
- //TODO: is it still needed here? (S.Zharko)
- if (ca::Framework::TrackingMode::kMcbm == fTrackingMode && th.z > 400) continue;
-
- if (fPerformance) {
- std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<ca::EDetectorID::kTof>(hitIndex);
- }
-
- th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
- if (th.IfAccept()) {
- tmpHits.push_back(th);
- nTofHits++;
- }
- else {
- //fMonitor.IncrementCounter(EMonitorKey::kDiscardedHitTof);
- }
- } // for j
- } // if listTofHits
-
-
- if (fVerbose > 1) {
- LOG(info) << "L1 ReadEvent: nhits mvd " << nMvdHits << " sts " << nStsHits << " much " << nMuchHits << " trd "
- << nTrdHits << " tof " << nTofHits << endl;
- }
-
- // total number of hits
- int nHits = nMvdHits + nStsHits + nMuchHits + nTrdHits + nTofHits;
-
- /*
- * Measured hits gathering: END
- */
-
- if (fVerbose >= 10) {
- cout << "ReadEvent: strips are read." << endl;
- }
-
-
- // ----- Fill and save fvExternalHits, fvHitDebugInfo and fvHitPointIndices vectors as well as fpData->vHits
-
- fvExternalHits.reserve(nHits);
- fvHitDebugInfo.reserve(nHits);
- fvHitPointIndices.reserve(nHits);
- fvHitBestPointIndices.reserve(nHits);
- fpIODataManager->ResetInputData(nHits);
- fpIODataManager->SetNhitKeys(NStrips);
-
- if (0) {
- Debugger::Instance().AddNtuple("hits", "ev:mcpoint:mc:mcMother:nmc:sta:x:y:z");
- for (int iHit = 0; iHit < nHits; ++iHit) {
- TmpHit& h = tmpHits[iHit];
- int mc = -1;
- int mcMother = -1;
- if (h.iBestMc >= 0) {
- mc = fvMCPoints[h.iBestMc].ID;
- mcMother = fvMCPoints[h.iBestMc].mother_ID;
- }
- Debugger::Instance().FillNtuple("hits", nCalls, h.iBestMc, mc, mcMother, h.vMc.size(), h.iStation, h.x, h.y, h.z);
- }
- }
-
- // ----- Fill
- for (int iHit = 0; iHit < nHits; ++iHit) {
- TmpHit& th = tmpHits[iHit];
-
- CbmL1HitDebugInfo s;
- s.Det = th.Det;
- s.ExtIndex = th.ExtIndex;
- s.iStation = th.iStation;
- s.x = th.x;
- s.y = th.y;
- s.z = th.z;
- s.dx = sqrt(th.dx2);
- s.dy = sqrt(th.dy2);
- s.dxy = th.dxy;
- s.time = th.time;
-
- assert(th.iStripF >= 0 || th.iStripF < NStrips);
- assert(th.iStripB >= 0 || th.iStripB < NStrips);
-
- ca::Hit h;
-
- h.SetFrontKey(th.iStripF);
- h.SetBackKey(th.iStripB);
- h.SetX(th.x);
- h.SetY(th.y);
- h.SetZ(th.z);
- h.SetT(th.time);
- h.SetDx2(th.dx2);
- h.SetDy2(th.dy2);
- h.SetDxy(th.dxy);
- h.SetDt2(th.dt2);
-
- h.SetRangeX(th.rangeX);
- h.SetRangeY(th.rangeY);
- h.SetRangeT(th.rangeT);
-
- h.SetId(iHit);
- h.SetStation(th.iStation);
-
- // save hit
- fvExternalHits.push_back(CbmL1HitId(th.Det, th.ExtIndex));
-
- // TODO: Here one should fill in the fvExternalHits[iHit].mcPointIds
-
- fpIODataManager->PushBackHit(h, th.fDataStream);
-
- fvHitDebugInfo.push_back(s);
- fvHitPointIndices.push_back(th.vMc);
- fvHitBestPointIndices.push_back(th.iBestMc);
- }
- if (fPerformance) {
- HitMatch();
- } /// OLD
-
- if (fVerbose >= 2) cout << "ReadEvent: mvd and sts are saved." << endl;
-
- // ----- Send data from IODataManager to ca::Framework --------------------------------------------------------------------
- if (1 == fSTAPDataMode) {
- WriteSTAPAlgoInputData(nCalls);
- }
- if (2 == fSTAPDataMode) {
- ReadSTAPAlgoInputData(nCalls);
- }
- // TODO: SZh: If we read data from file, we don't need to collect them above. This should be addressed
- fpAlgo->ReceiveInputData(fpIODataManager->TakeInputData());
-
- if (fPerformance) {
- if (fVerbose >= 10) cout << "HitMatch is done." << endl;
- if (fVerbose >= 10) cout << "MCPoints and MCTracks are saved." << endl;
- }
- if (fVerbose >= 10) cout << "ReadEvent is done." << endl;
-
- ++nCalls;
-} // void CbmL1::ReadEvent()
-//
-//--------------------------------------------------------------------------------------------------
-//
-void CbmL1::Fill_vMCTracks()
-{
- fvMCTracks.clear();
-
- // Count the total number of tracks in this event and reserve memory
- {
- Int_t nMCTracks = 0;
- for (DFSET::iterator set_it = fvSelectedMcEvents.begin(); set_it != fvSelectedMcEvents.end(); ++set_it) {
- Int_t iFile = set_it->first;
- Int_t iEvent = set_it->second;
- nMCTracks += fpMcTracks->Size(iFile, iEvent);
- }
- fvMCTracks.reserve(nMCTracks);
- }
-
- /* Loop over MC event */
- for (DFSET::iterator set_it = fvSelectedMcEvents.begin(); set_it != fvSelectedMcEvents.end(); ++set_it) {
- Int_t iFile = set_it->first;
- Int_t iEvent = set_it->second;
-
- auto header = dynamic_cast<FairMCEventHeader*>(fpMcEventHeader->Get(iFile, iEvent));
- assert(header);
-
- double eventTime = fpMcEventList->GetEventTime(iEvent, iFile);
-
- Int_t nMCTrack = fpMcTracks->Size(iFile, iEvent);
- /* Loop over MC tracks */
- for (Int_t iMCTrack = 0; iMCTrack < nMCTrack; iMCTrack++) {
- CbmMCTrack* MCTrack = dynamic_cast<CbmMCTrack*>(fpMcTracks->Get(iFile, iEvent, iMCTrack));
- if (!MCTrack) {
- continue;
- }
-
- int mother_ID = MCTrack->GetMotherId();
-
- int chainID = iMCTrack;
- int chainParent = mother_ID;
- while (chainParent >= 0) {
- chainID = chainParent;
- CbmMCTrack* parent = dynamic_cast<CbmMCTrack*>(fpMcTracks->Get(iFile, iEvent, chainParent));
- chainParent = parent->GetMotherId();
- }
-
- if (mother_ID < 0 && mother_ID != -2) mother_ID = -iEvent - 1;
- TVector3 vr;
- TLorentzVector vp;
- MCTrack->GetStartVertex(vr);
- MCTrack->Get4Momentum(vp);
- Int_t pdg = MCTrack->GetPdgCode();
- unsigned int processID = MCTrack->GetGeantProcessId();
- Double_t q = 0, mass = 0.;
- //if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
- // q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0;
- // mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
- //}
- // TODO: Add light nuclei (d, t, He3, He4): they are common in tracking but not accounted in TDatabasePDG (S.Zharko)
-
- q = MCTrack->GetCharge();
- mass = MCTrack->GetMass();
-
- Int_t iTrack = fvMCTracks.size(); //or iMCTrack?
- CbmL1MCTrack T(mass, q, vr, vp, iTrack, mother_ID, pdg, processID);
- T.chainID = chainID;
-
- // CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg);
- T.time = eventTime + MCTrack->GetStartT();
- T.iFile = iFile;
- T.iEvent = iEvent;
- // signal: primary tracks, displaced from the primary vertex
- T.isSignal = T.IsPrimary() && (T.z > header->GetZ() + 1.e-10);
-
- fvMCTracks.push_back(T);
- fmMCTracksLinksMap[CbmL1LinkKey(iMCTrack, iEvent, iFile)] = fvMCTracks.size() - 1;
- } //iTracks
- } //Links
-} //Fill_vMCTracks
-
-//
-//--------------------------------------------------------------------------------------------------
-//
-// TODO: Probably, we should reduce code here, rewriting this funciton as a template from CbmMvdPoint (S.Zharko)
-void CbmL1::ReadMCPoint(ca::EDetectorID iDet, int file, int event, int iPoint)
-{
- FairMCPoint* fairPoint = nullptr;
- TVector3 xyzO, PO;
-
- int iStLoc = -1;
- if (ca::EDetectorID::kMvd == iDet) {
- CbmMvdPoint* pt = dynamic_cast<CbmMvdPoint*>(fpMvdPoints->Get(file, event, iPoint));
- assert(pt);
- fairPoint = pt;
- pt->PositionOut(xyzO);
- pt->MomentumOut(PO);
- iStLoc = CbmMvdTrackingInterface::Instance()->GetTrackingStationIndex(pt);
- }
-
- if (ca::EDetectorID::kSts == iDet) {
- CbmStsPoint* pt = dynamic_cast<CbmStsPoint*>(fpStsPoints->Get(file, event, iPoint));
- assert(pt);
- fairPoint = pt;
- pt->PositionOut(xyzO);
- pt->MomentumOut(PO);
- iStLoc = CbmStsTrackingInterface::Instance()->GetTrackingStationIndex(pt);
- }
-
-
- if (ca::EDetectorID::kMuch == iDet) {
- CbmMuchPoint* pt = dynamic_cast<CbmMuchPoint*>(fpMuchPoints->Get(file, event, iPoint));
- assert(pt);
- fairPoint = pt;
- pt->PositionOut(xyzO);
- pt->MomentumOut(PO);
- iStLoc = CbmMuchTrackingInterface::Instance()->GetTrackingStationIndex(pt);
- }
-
-
- if (ca::EDetectorID::kTrd == iDet) {
- CbmTrdPoint* pt = dynamic_cast<CbmTrdPoint*>(fpTrdPoints->Get(file, event, iPoint));
- assert(pt);
- fairPoint = pt;
- pt->PositionOut(xyzO);
- pt->MomentumOut(PO);
- iStLoc = CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(pt);
- }
-
- if (ca::EDetectorID::kTof == iDet) {
- CbmTofPoint* pt = dynamic_cast<CbmTofPoint*>(fpTofPoints->Get(file, event, iPoint));
- assert(pt);
- fairPoint = pt;
- pt->Position(xyzO);
- pt->Momentum(PO);
- iStLoc = CbmTofTrackingInterface::Instance()->GetTrackingStationIndex(pt);
- }
-
- CbmL1MCPoint MC;
-
- MC.iStation = fpAlgo->GetParameters().GetStationIndexActive(iStLoc, iDet);
- if (MC.iStation < 0) {
- return;
- }
-
- TVector3 xyzI, PI;
- fairPoint->Position(xyzI);
- fairPoint->Momentum(PI);
- TVector3 xyz = .5 * (xyzI + xyzO);
- TVector3 P = .5 * (PI + PO);
-
- MC.x = xyz.X();
- MC.y = xyz.Y();
- MC.z = xyz.Z();
- MC.px = P.X();
- MC.py = P.Y();
- MC.pz = P.Z();
- MC.xIn = xyzI.X();
- MC.yIn = xyzI.Y();
- MC.zIn = xyzI.Z();
- MC.pxIn = PI.X();
- MC.pyIn = PI.Y();
- MC.pzIn = PI.Z();
- MC.xOut = xyzO.X();
- MC.yOut = xyzO.Y();
- MC.zOut = xyzO.Z();
- MC.pxOut = PO.X();
- MC.pyOut = PO.Y();
- MC.pzOut = PO.Z();
- MC.p = sqrt(fabs(MC.px * MC.px + MC.py * MC.py + MC.pz * MC.pz));
-
- CbmMCTrack* MCTrack = dynamic_cast<CbmMCTrack*>(fpMcTracks->Get(file, event, fairPoint->GetTrackID()));
-
- auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(fairPoint->GetTrackID(), event, file));
- assert(itTrack != fmMCTracksLinksMap.cend());
- MC.ID = itTrack->second;
-
- MC.pointId = iPoint;
- MC.file = file;
- MC.event = event;
-
- MC.time = fpMcEventList->GetEventTime(event, file) + fairPoint->GetTime();
-
- assert(MC.ID >= 0);
- if (MC.ID < 0) return;
-
-
- if (!MCTrack) {
- return;
- }
- MC.pdg = MCTrack->GetPdgCode();
- MC.mother_ID = MCTrack->GetMotherId();
-
- TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(MC.pdg);
- if (particlePDG) {
- MC.q = particlePDG->Charge() / 3.0;
- MC.mass = particlePDG->Mass();
- }
-
- // END OF CUTS. Please, don't put any cuts below!
-
- int iPointInt = iPoint;
- switch (iDet) {
- case ca::EDetectorID::kMvd: {
- ++fNpointsMvd;
- break;
- }
- case ca::EDetectorID::kSts:
- iPointInt += fNpointsMvdAll;
- ++fNpointsSts;
- break;
- case ca::EDetectorID::kMuch:
- iPointInt += fNpointsMvdAll + fNpointsStsAll;
- ++fNpointsMuch;
- break;
- case ca::EDetectorID::kTrd:
- iPointInt += fNpointsMvdAll + fNpointsStsAll + fNpointsMuchAll;
- ++fNpointsTrd;
- break;
- case ca::EDetectorID::kTof:
- iPointInt += fNpointsMvdAll + fNpointsStsAll + fNpointsMuchAll + fNpointsTrdAll;
- ++fNpointsTof;
- break;
- case ca::EDetectorID::kEND: break;
- }
-
- fmMCPointsLinksMap[CbmL1LinkKey(iPointInt, event, file)] = fvMCPoints.size();
- fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size());
- fvMCPoints.push_back(MC);
-}
-
-
-//--------------------------------------------------------------------------------------------------
-
-// TODO: Duplicated code (from CbmL1::ReadEvent)
-void CbmL1::HitMatch()
-{
- // Clear contents
- 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];
- 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/CbmL1Track.cxx b/reco/L1/CbmL1Track.cxx
index d7aec0a2f65f2d179bec9cc0e5b03a85eeec13dc..8b92c6cdc4bd7d3a35cdd9fc4c4216611c8605c3 100644
--- a/reco/L1/CbmL1Track.cxx
+++ b/reco/L1/CbmL1Track.cxx
@@ -23,8 +23,7 @@ std::string CbmL1Track::ToString(int verbose, bool header) const
msg << setw(12) << setfill(' ') << "Tx" << ' ';
msg << setw(8) << setfill(' ') << "N hits" << ' ';
msg << setw(8) << setfill(' ') << "N MC tracks" << ' ';
- msg << setw(12) << setfill(' ') << "MC tr (old)" << ' ';
- msg << setw(12) << setfill(' ') << "MC tr (new)" << ' ';
+ msg << setw(12) << setfill(' ') << "MC tr" << ' ';
msg << setw(8) << setfill(' ') << "IsGhost" << ' ';
msg << setw(8) << setfill(' ') << "Max.pur." << ' ';
msg << '\n';
@@ -36,7 +35,6 @@ std::string CbmL1Track::ToString(int verbose, bool header) const
msg << setw(12) << setfill(' ') << GetTx() << ' ';
msg << setw(8) << setfill(' ') << GetNofHits() << ' ';
msg << setw(8) << setfill(' ') << GetNofMCTracks() << ' ';
- msg << setw(12) << setfill(' ') << (mcTracks.size() ? mcTracks[0]->ID : -1) << ' ';
msg << setw(12) << setfill(' ') << GetMatchedMCTrackIndex() << ' ';
msg << setw(8) << setfill(' ') << IsGhost() << ' ';
msg << setw(8) << setfill(' ') << GetMaxPurity() << ' ';
@@ -46,14 +44,10 @@ std::string CbmL1Track::ToString(int verbose, bool header) const
for (int iH : GetHitIndexes()) {
msg << iH << ' ';
}
- msg << "\n\tmc tracks: (NEW)";
+ msg << "\n\tmc tracks:";
for (int iT : GetMCTrackIndexes()) {
msg << iT << ' ';
}
- msg << "\n\tmc tracks: (Old)";
- for (auto* iT : mcTracks) {
- msg << iT->ID << ' ';
- }
}
}
return msg.str();
diff --git a/reco/L1/CbmL1Track.h b/reco/L1/CbmL1Track.h
index beeea5ebcc4d01b1b8709c523e29fe2adeba8332..a89093b129a2e0cac2adb3b48ed210ef035ad0c2 100644
--- a/reco/L1/CbmL1Track.h
+++ b/reco/L1/CbmL1Track.h
@@ -44,21 +44,15 @@ class CbmL1Track : public cbm::algo::ca::TrackParamD {
public:
CbmL1Track() = default;
- /// Adds pointer to MC track (TODO: remove this and related methods)
- void AddMCTrack(CbmL1MCTrack* mcTr) { mcTracks.push_back(mcTr); }
-
/// Adds an index of MC track index
void AddMCTrackIndex(int iMT) { fvMcTrackIndexes.push_back_no_warning(iMT); }
/// Clears the contents of matched MC track indexes (and pointers)
void ClearMatchedMCTracks()
{
- mcTracks.clear();
fvMcTrackIndexes.clear();
}
- int GetNMCTracks() { return mcTracks.size(); }
-
bool IsGhost() const { return (maxPurity < CbmL1Constants::MinPurity); }
/// @brief Gets first hit index
@@ -71,23 +65,12 @@ class CbmL1Track : public cbm::algo::ca::TrackParamD {
const auto& GetHitIndexes() const { return Hits; }
- /// Gets number of MC tracks
- // TODO: Remove this method
- int GetNMCTracks() const { return mcTracks.size(); }
-
/// Gets a reference to MC track indexes
const auto& GetMCTrackIndexes() const { return fvMcTrackIndexes; }
/// Gets max purity
double GetMaxPurity() const { return maxPurity; }
- /// Gets container of pointers to MC tracks
- // TODO: this function is to be replaced with GetMCTrackIndexes()
- std::vector<CbmL1MCTrack*>& GetMCTracks() { return mcTracks; }
-
- /// Gets pointer of the associated MC track
- CbmL1MCTrack* GetMCTrack() { return mcTracks[0]; }
-
/// @brief Gets index of matched MC track
/// @note If two tracks are matched (should not happen, if purity > 50%), the first
///
@@ -128,11 +111,7 @@ class CbmL1Track : public cbm::algo::ca::TrackParamD {
std::map<int, int> hitMap; // N hits (second) from each mcTrack (first is a MC track ID) belong to current recoTrack
// FIXME: SZh 14.12.2022: map => unordered_map
-
private:
- // next members filled and used in Performance
- std::vector<CbmL1MCTrack*> mcTracks; // array of assosiated recoTracks. Should be only one.
-
ca::Vector<int> fvMcTrackIndexes = {"CbmL1Track::fvMcTrackIndexes"}; // global indexes of MC tracks
// NOTE: mcTracks should be replaced with fvMcTrackIndexes
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index cc0b4f001ec5f9c72f9bdbcbab6dda51dfa33e0c..3321b5fea2abfaab488db4da7e7e3472ca2fb204 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -78,9 +78,13 @@ void L1AlgoDraw::InitL1Draw(ca::Framework* algo_)
algo = algo_;
vHits.clear();
+ vHitsQa.clear();
vHits.reserve(algo->fWindowHits.size());
+ vHitsQa.reserve(algo->fWindowHits.size());
for (unsigned int i = 0; i < algo->fWindowHits.size(); i++) {
vHits.push_back(algo->fWindowHits[i]);
+ int iQaHit = algo->fWindowHits[i].Id();
+ vHitsQa.push_back(CbmL1::Instance()->GetQaHits()[iQaHit]);
}
NStations = algo->GetParameters().GetNstationsActive();
for (int i = 0; i < NStations; i++) {
@@ -103,32 +107,32 @@ void L1AlgoDraw::DrawMCTracks()
cout << "Gray - secondary p < 0.5 - (third\\second iteration)" << endl;
};
- for (vector<CbmL1MCTrack>::iterator it = L1.fvMCTracks.begin(); it != L1.fvMCTracks.end(); ++it) {
- CbmL1MCTrack& T = *it;
+ const auto& mcData = L1.GetMCData();
+ for (const auto& mcTrk : mcData.GetTrackContainer()) {
//draw reconstructable tracks only
- if (!T.IsReconstructable()) continue;
- //if (( T.mother_ID< 0 )&&( T.nContStations<3 )) continue;
- //if (( T.mother_ID>=0 )&&( T.nContStations<4 )) continue;
- //if( T.p<.2 ) continue;
- if (T.p < 0.1) continue;
+ if (!mcTrk.IsReconstructable()) continue;
+ //if (( mcTrk.GetMotherId() < 0 )&&( mcTrk.GetNofConsStationsWithHit() < 3 )) continue;
+ //if (( mcTrk.GetMotherId()>=0 )&&( mcTrk.GetNofConsStationsWithHit() <4 )) continue;
+ //if( mcTrk.GetP()<.2 ) continue;
+ if (mcTrk.GetP() < 0.1) continue;
pline.SetLineColor(kRed);
- if (T.p < 0.5) pline.SetLineColor(kBlue);
- if (T.mother_ID != -1) pline.SetLineColor(8);
- if ((T.mother_ID != -1) && (T.p < 0.5)) pline.SetLineColor(12);
- if (fVerbose >= 1)
- cout << "MC Track: p = " << T.p << " mother_ID = " << T.mother_ID << " PDG = " << T.pdg << " x,y,z = (" << T.x
- << ", " << T.y << ", " << T.z << ")" << endl;
+ if (mcTrk.GetP() < 0.5) pline.SetLineColor(kBlue);
+ if (mcTrk.GetMotherId() != -1) pline.SetLineColor(8);
+ if ((mcTrk.GetMotherId() != -1) && (mcTrk.GetMotherId() < 0.5)) pline.SetLineColor(12);
+ LOG_IF(info, fVerbose > 0) << "MC Track: p = " << mcTrk.GetP() << " mother_ID = " << mcTrk.GetMotherId()
+ << " PDG = " << mcTrk.GetPdgCode() << " x,y,z = (" << mcTrk.GetStartX() << ", "
+ << mcTrk.GetStartY() << ", " << mcTrk.GetStartZ() << ")";
double par[6];
- par[0] = T.x;
- par[1] = T.y;
- //if( fabs(T.pz)<0.1 ) continue;
- par[2] = T.px / T.pz;
- par[3] = T.py / T.pz;
- par[4] = T.q / T.p;
- par[5] = T.z;
-
- int npoints = T.Points.size();
- if (fVerbose >= 10) cout << " NMCPoints = " << npoints << endl;
+ par[0] = mcTrk.GetStartX();
+ par[1] = mcTrk.GetStartY();
+ //if( fabs(mcTrk.GetPz())<0.1 ) continue;
+ par[2] = mcTrk.GetTx();
+ par[3] = mcTrk.GetTy();
+ par[4] = mcTrk.GetCharge() / mcTrk.GetP();
+ par[5] = mcTrk.GetStartZ();
+
+ int npoints = mcTrk.GetNofPoints();
+ LOG_IF(info, fVerbose >= 10) << " NMCPoints = " << npoints;
if (npoints < 1) continue;
vector<double> lx, ly, lz;
@@ -139,26 +143,26 @@ void L1AlgoDraw::DrawMCTracks()
bool ok = true;
if (fVerbose >= 4) {
- cout << "hits = ";
- for (unsigned int ih = 0; ih < T.Hits.size(); ih++)
- cout << T.Hits[ih] << " ";
- cout << endl;
+ std::stringstream msg;
+ msg << "hits = ";
+ for (unsigned int ih : mcTrk.GetHitIndexes())
+ msg << ih << " ";
+ LOG(info) << msg.str();
}
- for (int ip = 0; ip < npoints; ip++) {
- CbmL1MCPoint& p = L1.fvMCPoints[T.Points[ip]];
+ for (int ip : mcTrk.GetPointIndexes()) {
+ auto& point = mcData.GetPoint(ip);
double par1[6];
//if( fabs(p.pz)<0.05 ) continue;
- par1[0] = p.x;
- par1[1] = p.y;
- par1[2] = p.px / p.pz;
- par1[3] = p.py / p.pz;
- par1[4] = p.q / p.p;
- par1[5] = p.z;
+ par1[0] = point.GetX();
+ par1[1] = point.GetY();
+ par1[2] = point.GetTx();
+ par1[3] = point.GetTy();
+ par1[4] = point.GetQp();
+ par1[5] = point.GetZ();
if (fVerbose >= 5) {
static fscal pz = -1;
- if (fabs(pz - p.z) > 1.0) cout << "-- ";
- cout << "point.z = " << p.z << endl;
- pz = p.z;
+ LOG(info) << (fabs(pz - point.GetZ()) > 1.0 ? "-- " : "") << "point z = " << point.GetZ();
+ pz = point.GetZ();
}
double Zfrst = par[5];
@@ -191,12 +195,12 @@ void L1AlgoDraw::DrawMCTracks()
lz.push_back(zl);
}
- par[0] = p.x;
- par[1] = p.y;
- par[2] = p.px / p.pz;
- par[3] = p.py / p.pz;
- par[4] = p.q / p.p;
- par[5] = p.z;
+ par[0] = point.GetX();
+ par[1] = point.GetY();
+ par[2] = point.GetTx();
+ par[3] = point.GetTy();
+ par[4] = point.GetQp();
+ par[5] = point.GetZ();
lx.push_back(par[0]);
ly.push_back(par[1]);
lz.push_back(par[5]);
@@ -569,7 +573,8 @@ void L1AlgoDraw::DrawInputHits()
Int_t n_poly_fake = 0;
for (int ih = HitsStartIndex[ista]; ih < HitsStopIndex[ista]; ih++) {
ca::Hit& h = vHits[ih];
- int iMC = CbmL1::Instance()->GetHitBestMcRefs()[ih];
+ const auto& hQa = vHitsQa[ih];
+ int iMC = hQa.GetBestMcPointId();
//if( (vSFlag[h.f] | vSFlagB[h.b] )&0x02 ) continue; // if used
fscal x = h.X(), y = h.Y(), z = h.Z();
@@ -702,7 +707,8 @@ void L1AlgoDraw::DrawRestHits(ca::HitIndex_t* StsRestHitsStartIndex, ca::HitInde
for (ca::HitIndex_t iRestHit = StsRestHitsStartIndex[ista]; iRestHit < StsRestHitsStopIndex[ista]; iRestHit++) {
int ih = realIHit[iRestHit];
ca::Hit& h = vHits[ih];
- int iMC = CbmL1::Instance()->GetHitBestMcRefs()[ih];
+ const auto& qaHit = vHitsQa[ih];
+ int iMC = qaHit.GetBestMcPointId();
//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/L1AlgoDraw.h b/reco/L1/L1Algo/utils/L1AlgoDraw.h
index d753fcb8eb0efe45487f4d15cfea2144fadca31c..881c4f17189e9aed2913c6bcf67858732c4ac01a 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.h
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.h
@@ -8,6 +8,7 @@
#include "CaHit.h"
#include "CaStation.h"
#include "CaTriplet.h"
+#include "CbmL1Hit.h"
#include <TString.h>
@@ -69,6 +70,7 @@ class L1AlgoDraw {
ca::Framework* algo{nullptr};
+ std::vector<CbmL1HitDebugInfo> vHitsQa{};
std::vector<ca::Hit> vHits{};
int HitsStartIndex[20]{0};
int HitsStopIndex[20]{0};
diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
index b405f3d44d0e7293d206476827e40e8b4ad16749..b8388d28d2ecd519b5d534995288d589918b17e2 100644
--- a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
+++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
@@ -164,16 +164,15 @@ void L1AlgoEfficiencyPerformance<NHits>::Init()
template<int NHits>
void L1AlgoEfficiencyPerformance<NHits>::FillMC()
{
- for (vector<CbmL1MCTrack>::iterator mtraIt = fL1->fvMCTracks.begin(); mtraIt != fL1->fvMCTracks.end(); mtraIt++) {
- CbmL1MCTrack& mtra = *(mtraIt);
- // if( ! mtra.IsReconstructable() ) continue;
-
- const int NMCPoints = mtra.Points.size();
+ const auto& mcData = fL1->GetMCData();
+ for (const auto& mcTrk : mcData.GetTrackContainer()) {
+ const int NMCPoints = mcData.GetNofPoints();
// const int NIter = mtra.NStations()-NHits+1; // number of possible tracklets
int lastIterSta = -1;
for (int iterOffset = 0; iterOffset < NMCPoints; iterOffset++) { // first mcPoint on the station
// const int iterMcId = mtra.Points[iterOffset];
- int iterSta = fL1->fvMCPoints[mtra.Points[iterOffset]].iStation;
+ int iterSta = mcData.GetPoint(mcTrk.GetPointIndexes()[iterOffset]).GetStationId();
+
if (iterSta == lastIterSta) continue; // find offset for next station
lastIterSta = iterSta;
@@ -182,14 +181,14 @@ void L1AlgoEfficiencyPerformance<NHits>::FillMC()
// TODO: Should we create all possible tracklets?
L1MCTracklet trlet; // TODO: don't use hits in mcTracklet
for (int is = 0, offset = iterOffset; ((offset < NMCPoints) && (is < NHits)); offset++) {
- const int mcId = mtra.Points[offset];
- CbmL1MCPoint* mcPs = &(fL1->fvMCPoints[mcId]);
- is = mcPs->iStation - iterSta;
+ const int mcId = mcTrk.GetPointIndexes()[offset];
+ const auto& mcPoint = mcData.GetPoint(mcId);
+ is = mcPoint.GetStationId() - iterSta;
if (is < NHits) {
- const int NPointHits = mcPs->hitIds.size();
+ const int NPointHits = mcPoint.GetHitIndexes().size();
for (int ih = 0; ih < NPointHits; ih++) {
- trlet.hitIds[is].push_back(mcPs->hitIds[ih]);
+ trlet.hitIds[is].push_back(mcPoint.GetHitIndexes()[ih]);
}
}
} // for is
@@ -201,10 +200,10 @@ void L1AlgoEfficiencyPerformance<NHits>::FillMC()
if (!good) continue; // can't create tracklet
trlet.iStation = iterSta;
- trlet.mcTrackId = mtra.ID;
- trlet.mcMotherTrackId = mtra.mother_ID;
- trlet.p = mtra.p;
- if (mtra.p > 1. / fL1->fpAlgo->GetMaxInvMom()) trlet.SetAsReconstructable();
+ trlet.mcTrackId = mcTrk.GetId();
+ trlet.mcMotherTrackId = mcTrk.GetMotherId();
+ trlet.p = mcTrk.GetP();
+ if (mcTrk.GetP() > 1. / fL1->fpAlgo->GetMaxInvMom()) trlet.SetAsReconstructable();
mcTracklets.push_back(trlet);
} // for Iter = stations
@@ -215,6 +214,7 @@ void L1AlgoEfficiencyPerformance<NHits>::FillMC()
template<int NHits>
bool L1AlgoEfficiencyPerformance<NHits>::AddOne(ca::HitIndex_t* iHits)
{
+ const auto& mcData = fL1->GetMCData();
L1RecoTracklet trlet(iHits);
// --- check if all hits belong to one mcTrack ---
@@ -224,7 +224,7 @@ bool L1AlgoEfficiencyPerformance<NHits>::AddOne(ca::HitIndex_t* iHits)
for (int iih = 0; iih < NHits; iih++) {
const int iMCP = fL1->fvHitDebugInfo[iHits[iih]].GetBestMcPointId();
if (iMCP > -1) {
- int mcId = fL1->fvMCPoints[iMCP].ID;
+ int mcId = mcData.GetPoint(iMCP).GetTrackId();
mcIds[iih].push_back(mcId);
}
} // for hits
@@ -254,7 +254,7 @@ bool L1AlgoEfficiencyPerformance<NHits>::AddOne(ca::HitIndex_t* iHits)
trlet.mcTrackId = mcsN[i];
int iMCP = fL1->fvHitDebugInfo[iHits[0]].GetBestMcPointId();
assert(iMCP > -1);
- trlet.iStation = fL1->fvMCPoints[fL1->fvHitDebugInfo[iHits[0]].GetBestMcPointId()].iStation;
+ trlet.iStation = mcData.GetPoint(fL1->fvHitDebugInfo[iHits[0]].GetBestMcPointId()).GetStationId();
break;
}
}
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
index cf6133e9654b2a7d914a4530e7b428f5dd3dfa47..9c4064ce91f1207ad7a0871223cab0707aeaaa2a 100644
--- a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
@@ -130,7 +130,7 @@ void L1AlgoPulls::AddOne(TrackParamV& T_, int i, ca::HitIndex_t ih)
err.qp = tr.GetQpError()[i];
// mc data
- int iMCP = fL1->GetHitBestMcRefs()[ih];
+ int iMCP = fL1->GetQaHits()[ih].GetBestMcPointId();
if (iMCP < 0) return;
CbmL1MCPoint& mcP = fL1->fvMCPoints[iMCP];
TL1TrackParameters mc(mcP);
diff --git a/reco/L1/catools/CaToolsMCPoint.h b/reco/L1/catools/CaToolsMCPoint.h
index b757c75d45cb586cd2b3b9c28b76c9fcbeb65af6..202375ab182262085c6db171cd7ede9f7e6ecd57 100644
--- a/reco/L1/catools/CaToolsMCPoint.h
+++ b/reco/L1/catools/CaToolsMCPoint.h
@@ -131,6 +131,24 @@ namespace cbm::ca::tools
return std::sqrt(fMomOut[0] * fMomOut[0] + fMomOut[1] * fMomOut[1] + fMomOut[2] * fMomOut[2]);
}
+ /// @brief Gets track transverse momentum at reference z of station [GeV/c]
+ double GetPt() const { return sqrt(fMom[0] * fMom[0] + fMom[1] * fMom[1]); }
+
+ /// @brief Gets track transverse momentum at entrance to station [GeV/c]
+ double GetPtIn() const
+ {
+ return sqrt(fMomIn[0] * fMomIn[0] + fMomIn[1] * fMomIn[1]);
+ ;
+ }
+
+ /// @brief Gets track transverse momentum at exit of station [GeV/c]
+ double GetPtOut() const
+ {
+ return sqrt(fMomOut[0] * fMomOut[0] + fMomOut[1] * fMomOut[1]);
+ ;
+ }
+
+
/// @brief Gets track momentum x component at reference z of station [GeV/c]
double GetPx() const { return fMom[0]; }
diff --git a/reco/L1/catools/CaToolsMCTrack.h b/reco/L1/catools/CaToolsMCTrack.h
index 059d012bf16352c939bf0dc447a7532a8c86ce28..b545e8c5e079ddfc6ff6bf098b3f05167c1f937a 100644
--- a/reco/L1/catools/CaToolsMCTrack.h
+++ b/reco/L1/catools/CaToolsMCTrack.h
@@ -83,6 +83,9 @@ namespace cbm::ca::tools
// ** Getters **
// *********************
+ /// Gets index of the first particle in the decay chain
+ int GetChainId() const { return fChainId; }
+
/// Gets charge [e]
double GetCharge() const { return fCharge; }
@@ -261,6 +264,9 @@ namespace cbm::ca::tools
// ** Setters **
// *********************
+ /// Sets index of the first particle in the decay chain
+ void SetChainId(int chainId) { fChainId = chainId; }
+
/// Sets charge [e]
void SetCharge(double q) { fCharge = q; }
@@ -350,6 +356,7 @@ namespace cbm::ca::tools
// Track address
int fId = constants::Undef<int>; ///< Index of MC track in internal container for TS/event
int fMotherId = constants::Undef<int>; ///< Index of mother MC track in the external tracks container
+ int fChainId = constants::Undef<int>; ///< Index of the first particle in the decay chain
LinkKey fLinkKey = {constants::Undef<int>, constants::Undef<int>,
constants::Undef<int>}; ///< A link key of this track in the external data structures
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index 198937bdc459290d11a8e684e189c43a6379da49..0db50630c9a4aa93f4342300964ca15aadfe87b3 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -470,6 +470,7 @@ InitStatus OutputQa::InitDataBranches()
fpTSReader->RegisterTracksContainer(fvRecoTracks);
fpTSReader->RegisterQaHitContainer(fvHits);
fpTSReader->RegisterHitIndexContainer(fvHitIds);
+ fpTSReader->SetSortQaHits(true);
if (!fpTSReader->InitRun()) {
return kFATAL;
}
@@ -932,8 +933,6 @@ InitStatus OutputQa::InitTimeSlice()
int nHits = 0;
int nRecoTracks = 0;
- LOG(info) << "\n\n\n\n\n !!!!! EVENT = " << this->GetCurrentEvent() << "\n\n\n\n\n";
-
// Read reconstructed input
fpTSReader->ReadEvent(this->GetCurrentEvent());
nHits = fvHits.size();
@@ -950,7 +949,8 @@ InitStatus OutputQa::InitTimeSlice()
fpMCModule->InitEvent(this->GetCurrentEvent());
nMCPoints = fMCData.GetNofPoints();
nMCTracks = fMCData.GetNofTracks();
- fpMCModule->MatchRecoAndMC();
+ fpMCModule->MatchHits();
+ fpMCModule->MatchTracks();
fMonitor.IncrementCounter(EMonitorKey::kMcPoint, nMCPoints);
fMonitor.IncrementCounter(EMonitorKey::kMcTrack, nMCTracks);