From 9db9c6aae811f470830042305cd7d8edea7f6e3b Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Wed, 10 May 2023 23:10:05 +0200
Subject: [PATCH] CA: QA updates
- Kaniadakis gaussian distribution is introduced to fit residuals of tracking parameters
- Reorganisation of MC module
---
core/base/CbmTrackingDetectorInterfaceBase.h | 3 +-
core/qa/CMakeLists.txt | 11 +-
core/qa/CbmQaCmpDrawer.cxx | 2 +
core/qa/CbmQaCmpDrawer.h | 15 +
core/qa/CbmQaTask.h | 2 +-
reco/L1/CMakeLists.txt | 2 +
reco/L1/CbmCaMCModule.cxx | 101 ++--
reco/L1/CbmCaMCModule.h | 458 ++++++++++---------
reco/L1/CbmCaTimeSliceReader.cxx | 142 +++---
reco/L1/CbmCaTimeSliceReader.h | 18 +-
reco/L1/CbmL1.h | 3 +-
reco/L1/CbmL1DetectorID.h | 19 +-
reco/L1/L1Algo/L1EArray.h | 43 ++
reco/L1/L1Algo/L1TrackPar.h | 8 +-
reco/L1/catools/CaToolsMCData.cxx | 19 +-
reco/L1/catools/CaToolsMCData.h | 95 +++-
reco/L1/qa/CbmCaInputQaSts.cxx | 44 +-
reco/L1/qa/CbmCaOutputQa.cxx | 66 ++-
reco/L1/qa/CbmCaOutputQa.h | 4 +-
reco/L1/qa/CbmCaTrackFitQa.cxx | 248 +++++++---
reco/L1/qa/CbmCaTrackFitQa.h | 150 ++++--
reco/L1/qa/CbmCaTrackTypeQa.cxx | 45 +-
reco/L1/qa/CbmCaTrackTypeQa.h | 26 +-
23 files changed, 959 insertions(+), 565 deletions(-)
create mode 100644 reco/L1/L1Algo/L1EArray.h
diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h
index d14c936430..e6cb0c9102 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.h
+++ b/core/base/CbmTrackingDetectorInterfaceBase.h
@@ -95,8 +95,7 @@ public:
/// Technical flag: true - all downcasts are done with dynamic_cast followed by checks for nullptr (increases
// computing time, better for debug), false - all downcasts are done with static_cast without sanity checks
- // (decreases computting time, but can cause bugs)
-
+ // (decreases computing time, but can cause bugs)
static constexpr bool kUseDynamicCast {true};
diff --git a/core/qa/CMakeLists.txt b/core/qa/CMakeLists.txt
index 0974dd5734..5c4eeb7907 100644
--- a/core/qa/CMakeLists.txt
+++ b/core/qa/CMakeLists.txt
@@ -45,7 +45,16 @@ set(INTERFACE_DEPENDENCIES
generate_cbm_library()
-Install(FILES CbmQaTask.h CbmQaCanvas.h CbmQaTable.h CbmQaHist.h CbmQaEff.h CbmQaPie.h CbmQaConstants.h CbmQaCmpDrawer.h checker/CbmQaCheckerTypedefs.h
+Install(FILES
+ CbmQaTask.h
+ CbmQaCanvas.h
+ CbmQaTable.h
+ CbmQaHist.h
+ CbmQaEff.h
+ CbmQaPie.h
+ CbmQaConstants.h
+ CbmQaCmpDrawer.h
+ checker/CbmQaCheckerTypedefs.h
DESTINATION include
)
diff --git a/core/qa/CbmQaCmpDrawer.cxx b/core/qa/CbmQaCmpDrawer.cxx
index 1e117bf50e..ed2b2305db 100644
--- a/core/qa/CbmQaCmpDrawer.cxx
+++ b/core/qa/CbmQaCmpDrawer.cxx
@@ -16,3 +16,5 @@ templateClassImp(CbmQaCmpDrawer);
template class CbmQaCmpDrawer<TH1F>;
template class CbmQaCmpDrawer<TProfile>;
+template class CbmQaCmpDrawer<TH2F>;
+template class CbmQaCmpDrawer<TProfile2D>;
diff --git a/core/qa/CbmQaCmpDrawer.h b/core/qa/CbmQaCmpDrawer.h
index d25fce7180..14832b880a 100644
--- a/core/qa/CbmQaCmpDrawer.h
+++ b/core/qa/CbmQaCmpDrawer.h
@@ -160,6 +160,21 @@ void CbmQaCmpDrawer<Obj>::Draw(Option_t* /*opt*/) const
(*it)->Draw(it == fvpObjects.begin() ? "" : "same");
}
}
+ else if constexpr (std::is_base_of_v<TH2, Obj> || std::is_base_of_v<TH3, Obj>) {
+ int nObjects = fvpObjects.size();
+ int nRows = static_cast<int>(std::sqrt(nObjects));
+ int nCols = nObjects / nRows;
+ if (nCols * nRows < nObjects) { ++nCols; }
+ gPad->Divide(nCols, nRows);
+
+ for (auto it = fvpObjects.begin(); it != fvpObjects.end(); ++it) {
+ if constexpr (std::is_base_of_v<TH2, Obj>) { (*it)->Draw("colz"); }
+ else {
+ (*it)->Draw();
+ }
+ }
+ return; // do not draw legend
+ }
// Draw legend
double legX0 = kLegRightBound - kLegEntryWidth;
diff --git a/core/qa/CbmQaTask.h b/core/qa/CbmQaTask.h
index 7f9e39e71b..da033d11dd 100644
--- a/core/qa/CbmQaTask.h
+++ b/core/qa/CbmQaTask.h
@@ -241,7 +241,7 @@ T* CbmQaTask::MakeHistoFromConfig(const char* nameBase, int id0, int id1, int id
double minZ = node["minz"].as<double>();
double maxZ = node["maxz"].as<double>();
std::string opt = node["opt"].as<std::string>("");
- pHist = new T(name.data(), title.data(), nBinsX, minX, maxX, nBinsY, minY, maxY, minZ, maxZ, opt.data());
+ pHist = new T(name.data(), title.data(), nBinsX, minX, maxX, nBinsY, minY, maxY, nBinsZ, minZ, maxZ, opt.data());
}
// 2D-histograms
else if constexpr (std::is_base_of_v<TH2, T>) {
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 36e1375d6d..53086f8c77 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -111,6 +111,7 @@ set(HEADERS
CbmL1Vtx.h
L1Algo/L1Def.h
L1Algo/L1Vector.h
+ L1Algo/L1EArray.h
L1Algo/L1Undef.h
catools/CaToolsWindowFinder.h
catools/CaToolsLinkKey.h
@@ -224,5 +225,6 @@ install(FILES L1Algo/L1Algo.h
L1Algo/L1Field.h
L1Algo/L1Hit.h
L1Algo/L1Vector.h
+ L1Algo/L1EArray.h
DESTINATION include/L1Algo
)
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index cfaea58318..9d55256837 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -41,9 +41,11 @@
// ** Action definition functions **
// *********************************
+using cbm::ca::MCModule;
+
// ---------------------------------------------------------------------------------------------------------------------
//
-bool CbmCaMCModule::InitRun()
+bool MCModule::InitRun()
try {
LOG(info) << "CA MC Module: initializing CA tracking Monte-Carlo module... ";
@@ -125,7 +127,7 @@ catch (const std::logic_error& error) {
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::InitEvent(CbmEvent* /*pEvent*/)
+void MCModule::InitEvent(CbmEvent* /*pEvent*/)
{
// Fill a set of file and event indexes
fFileEventIDs.clear();
@@ -153,18 +155,15 @@ void CbmCaMCModule::InitEvent(CbmEvent* /*pEvent*/)
aTrk.SortPointIndexes(
[&](const int& iPl, const int& iPr) { return fpMCData->GetPoint(iPl).GetZ() < fpMCData->GetPoint(iPr).GetZ(); });
}
-
- // ------ Match reconstructed and MC data
- this->MatchRecoAndMC();
}
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::ProcessEvent(CbmEvent*) {}
+void MCModule::ProcessEvent(CbmEvent*) {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::InitTrackInfo()
+void MCModule::InitTrackInfo()
{
// ----- Initialize stations arrangement and hit indexes
fpMCData->InitTrackInfo(*fpvQaHits);
@@ -207,7 +206,7 @@ void CbmCaMCModule::InitTrackInfo()
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::Finish() {}
+void MCModule::Finish() {}
// **********************************
@@ -216,7 +215,7 @@ void CbmCaMCModule::Finish() {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::MatchRecoAndMC()
+void MCModule::MatchRecoAndMC()
{
this->MatchPointsAndHits();
this->MatchRecoAndMCTracks();
@@ -227,13 +226,14 @@ void CbmCaMCModule::MatchRecoAndMC()
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMvd>(int iHitExt) const
+int MCModule::MatchHitWithMc<L1DetectorID::kMvd>(int iHitExt) const
{
int iPoint = -1;
if (fpMvdHitMatches) {
const auto* hitMatch = dynamic_cast<CbmMatch*>(fpMvdHitMatches->At(iHitExt));
assert(hitMatch);
- if (hitMatch->GetNofLinks() > 0 && hitMatch->GetLink(0).GetIndex() < fvNofPointsUsed[int(L1DetectorID::kMvd)]) {
+ if (hitMatch->GetNofLinks() > 0
+ && hitMatch->GetLink(0).GetIndex() < fpMCData->GetNofPointsUsed(L1DetectorID::kMvd)) {
iPoint = hitMatch->GetLink(0).GetIndex();
}
}
@@ -243,7 +243,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMvd>(int iHitExt) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHitExt) const
+int MCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHitExt) const
{
int iPoint = -1;
const auto* sh = dynamic_cast<CbmStsHit*>(fpStsHits->At(iHitExt));
@@ -288,7 +288,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHitExt) const
if (link.GetWeight() > bestWeight) {
bestWeight = link.GetWeight();
- iPoint = fpMCData->FindInternalPointIndex(CalcGlobMCPointIndex(iIndex, L1DetectorID::kSts), iEvent, iFile);
+ iPoint = fpMCData->FindInternalPointIndex(L1DetectorID::kSts, iIndex, iEvent, iFile);
assert(iPoint != -1);
}
}
@@ -299,13 +299,13 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHitExt) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHitExt) const
+int MCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHitExt) const
{
int iPoint = -1;
const auto* hitMatchMuch = dynamic_cast<CbmMatch*>(fpMuchHitMatches->At(iHitExt));
if (hitMatchMuch) {
for (int iLink = 0; iLink < hitMatchMuch->GetNofLinks(); ++iLink) {
- if (hitMatchMuch->GetLink(iLink).GetIndex() < fvNofPointsUsed[int(L1DetectorID::kMuch)]) {
+ if (hitMatchMuch->GetLink(iLink).GetIndex() < fpMCData->GetNofPointsUsed(L1DetectorID::kMuch)) {
int iMc = hitMatchMuch->GetLink(iLink).GetIndex();
if (iMc < 0) {
iPoint = -1;
@@ -313,7 +313,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHitExt) const
}
int iFile = hitMatchMuch->GetLink(iLink).GetFile();
int iEvent = hitMatchMuch->GetLink(iLink).GetEntry();
- iPoint = fpMCData->FindInternalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kMuch), iEvent, iFile);
+ iPoint = fpMCData->FindInternalPointIndex(L1DetectorID::kMuch, iMc, iEvent, iFile);
assert(iPoint != -1);
}
}
@@ -324,7 +324,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHitExt) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHitExt) const
+int MCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHitExt) const
{
int iPoint = -1;
const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTrdHitMatches->At(iHitExt));
@@ -332,14 +332,14 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHitExt) const
int iMC = -1;
if (hitMatch->GetNofLinks() > 0) {
iMC = hitMatch->GetLink(0).GetIndex();
- assert(iMC >= 0 && iMC < fvNofPointsUsed[int(L1DetectorID::kTrd)]);
+ assert(iMC >= 0 && iMC < fpMCData->GetNofPointsUsed(L1DetectorID::kTrd));
int iMc = hitMatch->GetLink(0).GetIndex();
if (iMc < 0) return iPoint;
int iFile = hitMatch->GetLink(0).GetFile();
int iEvent = hitMatch->GetLink(0).GetEntry();
- iPoint = fpMCData->FindInternalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kTrd), iEvent, iFile);
+ iPoint = fpMCData->FindInternalPointIndex(L1DetectorID::kTrd, iMc, iEvent, iFile);
assert(iPoint != -1);
}
}
@@ -349,7 +349,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHitExt) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
+int MCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
{
int iPoint = -1;
const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTofHitMatches->At(iHitExt));
@@ -361,7 +361,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
if (iMc < 0) return iPoint;
assert(iMc >= 0 && iMc < fpTofPoints->Size(iFile, iEvent));
- int iPointPrelim = fpMCData->FindInternalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kTof), iEvent, iFile);
+ int iPointPrelim = fpMCData->FindInternalPointIndex(L1DetectorID::kTof, iMc, iEvent, iFile);
if (iPointPrelim == -1) { continue; }
iPoint = iPointPrelim;
}
@@ -371,7 +371,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::MatchPointsAndHits()
+void MCModule::MatchPointsAndHits()
{
// Clear point indexes of each hit
// NOTE: at the moment only the best point is stored to the hit, but we keep the possibility for keeping many
@@ -420,7 +420,7 @@ void CbmCaMCModule::MatchPointsAndHits()
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::MatchRecoAndMCTracks()
+void MCModule::MatchRecoAndMCTracks()
{
for (int iTre = 0; iTre < int(fpvRecoTracks->size()); ++iTre) {
auto& trkRe = (*fpvRecoTracks)[iTre];
@@ -454,7 +454,7 @@ void CbmCaMCModule::MatchRecoAndMCTracks()
if (double(nHitsTrkMc) > double(nHitsTrkRe) * maxPurity) { maxPurity = double(nHitsTrkMc) / double(nHitsTrkRe); }
- ca::tools::MCTrack& trkMc = fpMCData->GetTrack(iTmc);
+ ::ca::tools::MCTrack& trkMc = fpMCData->GetTrack(iTmc);
// Match reconstructed and MC tracks, if purity with this MC track passes the threshold
if (double(nHitsTrkMc) >= CbmL1Constants::MinPurity * double(nHitsTrkRe)) {
@@ -480,7 +480,7 @@ void CbmCaMCModule::MatchRecoAndMCTracks()
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::SetDetector(L1DetectorID detID, bool flag)
+void MCModule::SetDetector(L1DetectorID detID, bool flag)
{
switch (detID) {
case L1DetectorID::kMvd: fbUseMvd = flag; break;
@@ -488,6 +488,7 @@ void CbmCaMCModule::SetDetector(L1DetectorID detID, bool flag)
case L1DetectorID::kMuch: fbUseMuch = flag; break;
case L1DetectorID::kTrd: fbUseTrd = flag; break;
case L1DetectorID::kTof: fbUseTof = flag; break;
+ case L1DetectorID::kEND: break;
}
}
@@ -497,7 +498,7 @@ void CbmCaMCModule::SetDetector(L1DetectorID detID, bool flag)
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::CheckInit() const
+void MCModule::CheckInit() const
{
// Check parameters
if (!fpParameters.get()) { throw std::logic_error("Tracking parameters object was not defined"); }
@@ -551,9 +552,8 @@ void CbmCaMCModule::CheckInit() const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* pPoints)
+void MCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* pPoints)
{
- fvNofPointsUsed[int(L1DetectorID::kTof)] = 0;
for (const auto& key : fFileEventIDs) {
int iFile = key.first;
int iEvent = key.second;
@@ -586,7 +586,7 @@ void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray*
}
int nPointsEvent = pPoints->Size(iFile, iEvent);
- // loop over all TOF points in event and ffvNofPointsUsed[int(L1DetectorID::kMvd)]ile
+ // loop over all TOF points in event and file
for (int iP = 0; iP < nPointsEvent; ++iP) {
if (!vbTofPointMatched[iP]) { continue; } // iP was not matched
@@ -627,14 +627,14 @@ void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray*
if (iStActive < 0) { continue; }
for (int iTrk = 0; iTrk < (int) vTofPointToTrack[iStLocGeo].size(); ++iTrk) {
if (vTofPointToTrack[iStLocGeo][iTrk] < 0) { continue; }
- ca::tools::MCPoint aPoint;
+ ::ca::tools::MCPoint aPoint;
if (FillMCPoint<L1DetectorID::kTof>(vTofPointToTrack[iStLocGeo][iTrk], iEvent, iFile, aPoint)) {
aPoint.SetStationId(iStActive);
- aPoint.SetExternalId(CalcGlobMCPointIndex(vTofPointToTrack[iStLocGeo][iTrk], L1DetectorID::kTof));
+ int iExtGlob = fpMCData->GetPointGlobExtIndex(L1DetectorID::kTof, vTofPointToTrack[iStLocGeo][iTrk]);
+ aPoint.SetExternalId(iExtGlob);
int iPInt = fpMCData->GetNofPoints(); // assign an index of point in internal container
if (aPoint.GetTrackId() > -1) { fpMCData->GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); }
fpMCData->AddPoint(aPoint);
- ++fvNofPointsUsed[int(L1DetectorID::kTof)];
}
} // loop over tracks: end
} // loop over geometry stations: end
@@ -643,24 +643,33 @@ void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray*
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::ReadMCPoints()
+void MCModule::ReadMCPoints()
{
int nPointsEstimated = 5 * fpMCData->GetNofTracks() * fpParameters->GetNstationsActive();
fpMCData->ReserveNofPoints(nPointsEstimated);
- // ----- Initialise the number of points
- std::fill(fvNofPointsOrig.begin(), fvNofPointsOrig.end(), 0);
-
+ int nPointsMvd = 0;
+ int nPointsSts = 0;
+ int nPointsMuch = 0;
+ int nPointsTrd = 0;
+ int nPointsTof = 0;
for (const auto& key : fFileEventIDs) {
int iFile = key.first;
int iEvent = key.second;
- if (fbUseMvd) { fvNofPointsOrig[int(L1DetectorID::kMvd)] += fpMvdPoints->Size(iFile, iEvent); }
- if (fbUseSts) { fvNofPointsOrig[int(L1DetectorID::kSts)] += fpStsPoints->Size(iFile, iEvent); }
- if (fbUseMuch) { fvNofPointsOrig[int(L1DetectorID::kMuch)] += fpMuchPoints->Size(iFile, iEvent); }
- if (fbUseTrd) { fvNofPointsOrig[int(L1DetectorID::kTrd)] += fpTrdPoints->Size(iFile, iEvent); }
- if (fbUseTof) { fvNofPointsOrig[int(L1DetectorID::kTof)] += fpTofPoints->Size(iFile, iEvent); }
+
+ if (fbUseMvd) { nPointsMvd += fpMvdPoints->Size(iFile, iEvent); }
+ if (fbUseSts) { nPointsSts += fpStsPoints->Size(iFile, iEvent); }
+ if (fbUseMuch) { nPointsMuch += fpMuchPoints->Size(iFile, iEvent); }
+ if (fbUseTrd) { nPointsTrd += fpTrdPoints->Size(iFile, iEvent); }
+ if (fbUseTof) { nPointsTof += fpTofPoints->Size(iFile, iEvent); }
}
+ fpMCData->SetNofPointsOrig(L1DetectorID::kMvd, nPointsMvd);
+ fpMCData->SetNofPointsOrig(L1DetectorID::kSts, nPointsSts);
+ fpMCData->SetNofPointsOrig(L1DetectorID::kMuch, nPointsMuch);
+ fpMCData->SetNofPointsOrig(L1DetectorID::kTrd, nPointsTrd);
+ fpMCData->SetNofPointsOrig(L1DetectorID::kTof, nPointsTof);
+
// ----- Read MC points in MVD, STS, MuCh, TRD and TOF
if (fbUseMvd) { this->ReadMCPointsForDetector<L1DetectorID::kMvd>(fpMvdPoints); }
if (fbUseSts) { this->ReadMCPointsForDetector<L1DetectorID::kSts>(fpStsPoints); }
@@ -671,7 +680,7 @@ void CbmCaMCModule::ReadMCPoints()
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::ReadMCTracks()
+void MCModule::ReadMCTracks()
{
// ----- Total number of tracks
int nTracksTot = 0;
@@ -686,7 +695,7 @@ void CbmCaMCModule::ReadMCTracks()
int iEvent = key.second;
auto pEvtHeader = dynamic_cast<FairMCEventHeader*>(fpMCEventHeader->Get(iFile, iEvent));
if (!pEvtHeader) {
- LOG(fatal) << "CbmCaMCModule: event header is not found for file " << iFile << " and event " << iEvent;
+ LOG(fatal) << "cbm::ca::MCModule: event header is not found for file " << iFile << " and event " << iEvent;
}
double eventTime = fpMCEventList->GetEventTime(iEvent, iFile);
@@ -695,11 +704,11 @@ void CbmCaMCModule::ReadMCTracks()
for (int iTrkExt = 0; iTrkExt < nTracks; ++iTrkExt) {
CbmMCTrack* pExtMCTrk = dynamic_cast<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTrkExt));
if (!pExtMCTrk) {
- LOG(warn) << "CbmCaMCModule: track with (iF, iE, iT) = " << iFile << ' ' << iEvent << ' ' << iTrkExt
+ LOG(warn) << "cbm::ca::MCModule: track with (iF, iE, iT) = " << iFile << ' ' << iEvent << ' ' << iTrkExt
<< " not found";
}
// Create a CbmL1MCTrack
- ca::tools::MCTrack aTrk;
+ ::ca::tools::MCTrack aTrk;
aTrk.SetId(fpMCData->GetNofTracks()); // assign current number of tracks read so far as an ID of a new track
aTrk.SetExternalId(iTrkExt); // external index of track is its index from CbmMCTrack objects container
@@ -722,7 +731,7 @@ void CbmCaMCModule::ReadMCTracks()
aTrk.SetMass(pExtMCTrk->GetMass());
// The charge in CbmMCTrack is similarly to mass defined from ROOT PDG data base. The value of charge there is
- // given in the units of 1/3e (absolute value of d-quark charge). In ca::tools::MCTrack we recalculate this
+ // given in the units of 1/3e (absolute value of d-quark charge). In ::ca::tools::MCTrack we recalculate this
// value to the units of e.
aTrk.SetCharge(pExtMCTrk->GetCharge() / 3.);
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
index c5ceef1cab..95e5ecf167 100644
--- a/reco/L1/CbmCaMCModule.h
+++ b/reco/L1/CbmCaMCModule.h
@@ -47,229 +47,237 @@ class CbmL1Track;
enum class L1DetectorID;
-/// Class CbmCaPerformance is an interface to communicate between
-///
-class CbmCaMCModule {
-public:
- /// @brief Constructor
- /// @param verbosity Verbosity level
- /// @param perfMode Performance mode (defines cut on number of consecutive stations with hit or point)
- CbmCaMCModule(int verb = 0, int perfMode = 3) : fVerbose(verb), fPerformanceMode(perfMode) {}
-
- /// @brief Destructor
- ~CbmCaMCModule() = default;
-
- /// @brief Copy constructor
- CbmCaMCModule(const CbmCaMCModule&) = delete;
-
- /// @brief Move constructor
- CbmCaMCModule(CbmCaMCModule&&) = delete;
-
- /// @brief Copy assignment operator
- CbmCaMCModule& operator=(const CbmCaMCModule&) = delete;
-
- /// @brief Move assignment operator
- CbmCaMCModule& operator=(CbmCaMCModule&&) = delete;
-
- /// @brief Defines performance action in the end of the run
- void Finish();
-
- /// @brief Gets the first point index for a given detector subsystem
- int GetFirstPointIndex(L1DetectorID detID) const
- {
- return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + int(detID), 0);
- }
-
- /// @brief Gets the first point index for a given detector subsystem
- int GetLastPointIndex(L1DetectorID detID) const
- {
- return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + int(detID) + 1, 0) - 1;
- }
-
- /// @brief Defines performance action in the beginning of each event or time slice
- /// @note This function should be called before hits initialization
- /// @param pEvent Pointer to a current CbmEvent
- void InitEvent(CbmEvent* pEvent);
-
- /// @brief Defines action on the module in the beginning of the run
- /// @return Success flag
- bool InitRun();
-
- /// @brief Initializes MC track
- ///
- /// Initializes information about arrangement of points and hits of MC tracks within stations and the status
- /// of track ability to be reconstructed, calculates max number of points and hits on a station, number of
- /// consecutive stations containing a hit or point and number of stations and points with hits.
- void InitTrackInfo();
-
- /// @brief Match reconstructed and MC data
- ///
- /// Runs matching of hits with points and reconstructed tracks with
- void MatchRecoAndMC();
-
- /// @brief Processes event
- ///
- /// Fills histograms and tables, should be called after the tracking done
- void ProcessEvent(CbmEvent* pEvent);
-
- /// @brief Sets first hit indexes container in different detectors
- /// @param source Array of indexes
- void RegisterFirstHitIndexes(const std::array<int, L1Constants::size::kMaxNdetectors + 1>& source)
- {
- fpvFstHitId = &source;
- }
-
- /// @brief Sets used detector subsystems
- /// @param detID Id of detector
- /// @param flag Flag: true - detector is used
- /// @note Should be called before this->Init()
- void SetDetector(L1DetectorID detID, bool flag);
-
- /// @brief Sets legacy event mode:
- /// @param flag Flag:
- /// - true: runs on events base,
- /// - false: runs on time-slice base
- void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; }
-
- /// @brief Registers MC data object
- /// @param mcData Instance of MC data
- void RegisterMCData(ca::tools::MCData& mcData) { fpMCData = &mcData; }
-
- /// @brief Registers reconstructed track container
- /// @param vRecoTrack Reference to reconstructed track container
- void RegisterRecoTrackContainer(L1Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; }
-
- /// @brief Registers hit index container
- /// @param vHitIds Reference to hit index container
- void RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; }
-
- /// @brief Registers CA parameters object
- /// @param pParameters A shared pointer to the parameters object
- void RegisterParameters(std::shared_ptr<L1Parameters>& pParameters) { fpParameters = pParameters; }
-
- /// @brief Registers debug hit container
- /// @param vQaHits Reference to debug hit container
- void RegisterQaHitContainer(L1Vector<CbmL1HitDebugInfo>& vQaHits) { fpvQaHits = &vQaHits; }
-
-private:
- /// @brief Calculates global index of MC point
- /// @param iPointLocal Local index of MC point
- /// @param detID Detector ID
- ///
- /// The function calculates global index of MC point as a sum of a given local index and total provided number of
- /// points in previous detector subsystem.
- int CalcGlobMCPointIndex(int iPointLocal, L1DetectorID detID) const
- {
- return iPointLocal + std::accumulate(fvNofPointsOrig.cbegin(), fvNofPointsOrig.cbegin() + int(detID), 0);
- }
-
- /// @brief Check class initialization
- /// @note The function throws std::logic_error, if initialization is incomplete
- void CheckInit() const;
-
- /// @brief Matches hit with MC point
- /// @tparam DetId Detector ID
- /// @param iHitExt External index of hit
- /// @return MC-point index in fvMCPoints array
+namespace cbm::ca
+{
+ /// @brief Class CbmCaPerformance is an interface to communicate between
///
- /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
- /// link in the match and returns the corresponding global ID of the MC points.
- template<L1DetectorID DetId>
- int MatchHitWithMc(int iHitExt) const;
+ class MCModule {
+ public:
+ /// @brief Constructor
+ /// @param verbosity Verbosity level
+ /// @param perfMode Performance mode (defines cut on number of consecutive stations with hit or point)
+ MCModule(int verb = 0, int perfMode = 3) : fVerbose(verb), fPerformanceMode(perfMode) {}
+
+ /// @brief Destructor
+ ~MCModule() = default;
+
+ /// @brief Copy constructor
+ MCModule(const MCModule&) = delete;
+
+ /// @brief Move constructor
+ MCModule(MCModule&&) = delete;
+
+ /// @brief Copy assignment operator
+ MCModule& operator=(const MCModule&) = delete;
+
+ /// @brief Move assignment operator
+ MCModule& operator=(MCModule&&) = delete;
+
+ /// @brief Creates hits from MC points
+ ///
+ /// This function creates hits from MC points in selected tracking stations and re-fills hit arrays.
+ void CreateHitsFromPoints() {}
+
+ /// @brief Sets hit parameters from the matched MC point
+ void AdjustHitsToPoints();
+
+ /// @brief Defines performance action in the end of the run
+ void Finish();
+
+ /// @brief Defines performance action in the beginning of each event or time slice
+ /// @note This function should be called before hits initialization
+ /// @param pEvent Pointer to a current CbmEvent
+ void InitEvent(CbmEvent* pEvent);
+
+ /// @brief Defines action on the module in the beginning of the run
+ /// @return Success flag
+ bool InitRun();
+
+ /// @brief Match reconstructed and MC data
+ ///
+ /// Runs matching of hits with points and reconstructed tracks with MC ones. Reconstructed tracks are updated with
+ /// true information.
+ void MatchRecoAndMC();
+
+ /// @brief Processes event
+ ///
+ /// Fills histograms and tables, should be called after the tracking done
+ void ProcessEvent(CbmEvent* pEvent);
+
+ /// @brief Sets first hit indexes container in different detectors
+ /// @param source Array of indexes
+ void RegisterFirstHitIndexes(const std::array<int, L1Constants::size::kMaxNdetectors + 1>& source)
+ {
+ fpvFstHitId = &source;
+ }
- /// @brief Match MC points and reconstructed hits
- ///
- /// Writes indexes of MC points to each hit and indexes of hits to each MC point.
- void MatchPointsAndHits();
+ /// @brief Sets used detector subsystems
+ /// @param detID Id of detector
+ /// @param flag Flag: true - detector is used
+ /// @note Should be called before this->Init()
+ void SetDetector(L1DetectorID detID, bool flag);
+
+ /// @brief Sets legacy event mode:
+ /// @param flag Flag:
+ /// - true: runs on events base,
+ /// - false: runs on time-slice base
+ void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; }
+
+ /// @brief Creates hits from MC points for a particular detector and station
+ /// @param detID DetectorID
+ /// @param iStLoc Local index of station, provided by geometry
+ /// @param du Error u-coordinate [cm]
+ /// @param dv Error v-coordinate [cm]
+ /// @param dt Error of time [ns]
+ /// @param ifSmear Flag:
+ /// - true: MC point quantities are smeared along the error by gaussian,
+ /// - false: Precise point quantities are used
+ ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, double du, double dv, double dt, bool ifSmear) {}
+
+ /// @brief Adjusts existing reconstructed hits to MC points for a particular detector and station
+ /// @param detID DetectorID
+ /// @param iStLoc Local index of station, provided by geometry
+ /// @param ifSmear Flag:
+ /// - true: MC point quantities are smeared along the error by gaussian,
+ /// - false: Precise point quantities are used
+ ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, bool ifSmear) {}
+
+
+ /// @brief Registers MC data object
+ /// @param mcData Instance of MC data
+ void RegisterMCData(::ca::tools::MCData& mcData) { fpMCData = &mcData; }
+
+ /// @brief Registers reconstructed track container
+ /// @param vRecoTrack Reference to reconstructed track container
+ void RegisterRecoTrackContainer(L1Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; }
+
+ /// @brief Registers hit index container
+ /// @param vHitIds Reference to hit index container
+ void RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; }
+
+ /// @brief Registers CA parameters object
+ /// @param pParameters A shared pointer to the parameters object
+ void RegisterParameters(std::shared_ptr<L1Parameters>& pParameters) { fpParameters = pParameters; }
+
+ /// @brief Registers debug hit container
+ /// @param vQaHits Reference to debug hit container
+ void RegisterQaHitContainer(L1Vector<CbmL1HitDebugInfo>& vQaHits) { fpvQaHits = &vQaHits; }
+
+ private:
+ /// @brief Check class initialization
+ /// @note The function throws std::logic_error, if initialization is incomplete
+ void CheckInit() const;
+
+ /// @brief Initializes MC track
+ ///
+ /// Initializes information about arrangement of points and hits of MC tracks within stations and the status
+ /// of track ability to be reconstructed, calculates max number of points and hits on a station, number of
+ /// consecutive stations containing a hit or point and number of stations and points with hits.
+ void InitTrackInfo();
+
+
+ /// @brief Matches hit with MC point
+ /// @tparam DetId Detector ID
+ /// @param iHitExt External index of hit
+ /// @return MC-point index in fvMCPoints array
+ ///
+ /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
+ /// link in the match and returns the corresponding global ID of the MC points.
+ template<L1DetectorID DetId>
+ int MatchHitWithMc(int iHitExt) const;
+
+ /// @brief Match MC points and reconstructed hits
+ ///
+ /// Writes indexes of MC points to each hit and indexes of hits to each MC point.
+ void MatchPointsAndHits();
+
+ /// @brief Matches reconstructed tracks with MC tracks
+ ///
+ /// In the procedure, the maps of associated MC track indexes to corresponding number of hits are filled out and the
+ /// max purity is calculated for each reconstructed track in the TS/event. If the value of purity for a given MC track
+ /// is larger then a threshold, the corresponding track index is saved to the reconstructed track object, in the same
+ /// time the index of the reconstructed track object is saved to this MC track object. If purity for the MC track does
+ /// not pass the threshold, its index will not be accounted as a matched track, and the index of reconstructed track
+ /// will be added as an index of "touched" track.
+ void MatchRecoAndMCTracks();
+
+ /// @brief Reads MC tracks from external trees and saves them to MCDataObject
+ void ReadMCTracks();
+
+ /// @brief Reads MC points from external trees and saves them to MCDataObject
+ void ReadMCPoints();
+
+ /// @brief Reads MC points in particular detector
+ template<L1DetectorID DetID>
+ void ReadMCPointsForDetector(CbmMCDataArray* pPoints);
+
+ /// @brief Fills a single detector-specific MC point
+ /// @tparam DetID Detector subsystem ID
+ /// @param[in] iExtId Index of point in the external points container
+ /// @param[in] iEvent EventID of point in the external points container
+ /// @param[in] iFile FileID of point int the external points container
+ /// @param[out] intMCPoint Reference to the internal tracking MC point object (ca::tools::MCData)
+ /// @return true Point is correct and is to be saved to the MCData object
+ /// @return false Point is incorrect and will be ignored
+ template<L1DetectorID DetID>
+ bool FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point);
+
+ std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
+
+ // ------ Flags
+ bool fbUseMvd = false; ///< is MVD used
+ bool fbUseSts = false; ///< is STS used
+ bool fbUseMuch = false; ///< is MuCh used
+ bool fbUseTrd = false; ///< is TRD used
+ bool fbUseTof = false; ///< is TOF used
+
+ bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
+ int fVerbose = 0; ///< Verbosity level
+ int fPerformanceMode = -1; ///< Mode of performance
+
+ // ------ Input data branches
+ const CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice
+
+ // Mc-event
+ CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
+ CbmMCDataObject* fpMCEventHeader = nullptr; ///< MC event header
+ CbmMCDataArray* fpMCTracks = nullptr; ///< MC tracks input
+
+ // Mc-points
+ CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points input container
+ CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points input container
+ CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points input container
+ CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container
+ CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container
+
+ // Matches
+ TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
+ TClonesArray* fpStsHitMatches = nullptr; ///< Array of STS hit matches
+ TClonesArray* fpMuchHitMatches = nullptr; ///< Array of MuCh hit matches
+ TClonesArray* fpTrdHitMatches = nullptr; ///< Array of TRD hit matches
+ TClonesArray* fpTofHitMatches = nullptr; ///< Array of TOF hit matches
+
+ TClonesArray* fpStsHits = nullptr; ///< Array of STS hits (currently needed for matching)
+ TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches
+
+ // Matching information
+ std::set<std::pair<int, int>> fFileEventIDs; ///< Set of file and event indexes: first - iFile, second - iEvent
+
+ // MC data
+ ::ca::tools::MCData* fpMCData = nullptr; ///< MC information (hits and tracks) instance
+
+ // Reconstructed data
+ L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
+ L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to hit index container
+ L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to QA hit container
+
+ /// @brief Pointer to array of first hit indexes in the detector subsystem
+ ///
+ /// This array must be initialized in the run initialization function.
+ const std::array<int, L1Constants::size::kMaxNdetectors + 1>* fpvFstHitId = nullptr;
+ };
+} // namespace cbm::ca
- /// @brief Matches reconstructed tracks with MC tracks
- ///
- /// In the procedure, the maps of associated MC track indexes to corresponding number of hits are filled out and the
- /// max purity is calculated for each reconstructed track in the TS/event. If the value of purity for a given MC track
- /// is larger then a threshold, the corresponding track index is saved to the reconstructed track object, in the same
- /// time the index of the reconstructed track object is saved to this MC track object. If purity for the MC track does
- /// not pass the threshold, its index will not be accounted as a matched track, and the index of reconstructed track
- /// will be added as an index of "touched" track.
- void MatchRecoAndMCTracks();
-
- /// @brief Reads MC tracks from external trees and saves them to MCDataObject
- void ReadMCTracks();
-
- /// @brief Reads MC points from external trees and saves them to MCDataObject
- void ReadMCPoints();
-
- /// @brief Reads MC points in particular detector
- template<L1DetectorID DetID>
- void ReadMCPointsForDetector(CbmMCDataArray* pPoints);
-
- /// @brief Fills a single detector-specific MC point
- /// @tparam DetID Detector subsystem ID
- /// @param[in] iExtId Index of point in the external points container
- /// @param[in] iEvent EventID of point in the external points container
- /// @param[in] iFile FileID of point int the external points container
- /// @param[out] intMCPoint Reference to the internal tracking MC point object (ca::tools::MCData)
- /// @return true Point is correct and is to be saved to the MCData object
- /// @return false Point is incorrect and will be ignored
- template<L1DetectorID DetID>
- bool FillMCPoint(int iExtId, int iEvent, int iFile, ca::tools::MCPoint& point);
-
- std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
-
- // ------ Flags
- bool fbUseMvd = false; ///< is MVD used
- bool fbUseSts = false; ///< is STS used
- bool fbUseMuch = false; ///< is MuCh used
- bool fbUseTrd = false; ///< is TRD used
- bool fbUseTof = false; ///< is TOF used
-
- bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
- int fVerbose = 0; ///< Verbosity level
- int fPerformanceMode = -1; ///< Mode of performance
-
- // ------ Input data branches
- const CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice
-
- // Mc-event
- CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
- CbmMCDataObject* fpMCEventHeader = nullptr; ///< MC event header
- CbmMCDataArray* fpMCTracks = nullptr; ///< MC tracks input
-
- // Mc-points
- CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points input container
- CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points input container
- CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points input container
- CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container
- CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container
-
- std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsOrig = {0}; ///< Number of points by detector provided
- std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsUsed = {0}; ///< Number of points used in performance
-
- // Matches
- TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
- TClonesArray* fpStsHitMatches = nullptr; ///< Array of STS hit matches
- TClonesArray* fpMuchHitMatches = nullptr; ///< Array of MuCh hit matches
- TClonesArray* fpTrdHitMatches = nullptr; ///< Array of TRD hit matches
- TClonesArray* fpTofHitMatches = nullptr; ///< Array of TOF hit matches
-
- TClonesArray* fpStsHits = nullptr; ///< Array of STS hits (currently needed for matching)
- TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches
-
- // Matching information
- std::set<std::pair<int, int>> fFileEventIDs; ///< Set of file and event indexes: first - iFile, second - iEvent
-
- // MC data
- ca::tools::MCData* fpMCData = nullptr; ///< MC information (hits and tracks) instance
-
- // Reconstructed data
- L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
- L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to hit index container
- L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to QA hit container
-
- /// @brief Pointer to array of first hit indexes in the detector subsystem
- ///
- /// This array must be initialized in the run initialization function.
- const std::array<int, L1Constants::size::kMaxNdetectors + 1>* fpvFstHitId = nullptr;
-};
// **********************************************
// ** Template function implementation **
@@ -278,7 +286,7 @@ private:
// ---------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
-bool CbmCaMCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ca::tools::MCPoint& point)
+bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point)
{
// ----- Get positions, momenta, time and track ID
TVector3 posIn; // Position at entrance to station [cm]
@@ -467,22 +475,20 @@ bool CbmCaMCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ca::tools::MC
// ---------------------------------------------------------------------------------------------------------------------
// NOTE: template is used, because another template function FillMCPoint is used inside
template<L1DetectorID DetID>
-void CbmCaMCModule::ReadMCPointsForDetector(CbmMCDataArray* pPoints)
+void cbm::ca::MCModule::ReadMCPointsForDetector(CbmMCDataArray* pPoints)
{
if constexpr (L1DetectorID::kTof != DetID) {
- fvNofPointsUsed[int(DetID)] = 0;
for (const auto& key : fFileEventIDs) {
int iFile = key.first;
int iEvent = key.second;
int nPointsEvent = pPoints->Size(iFile, iEvent);
for (int iP = 0; iP < nPointsEvent; ++iP) {
- ca::tools::MCPoint aPoint;
+ ::ca::tools::MCPoint aPoint;
if (FillMCPoint<DetID>(iP, iEvent, iFile, aPoint)) {
- aPoint.SetExternalId(CalcGlobMCPointIndex(iP, DetID));
+ aPoint.SetExternalId(fpMCData->GetPointGlobExtIndex(DetID, iP));
int iPInt = fpMCData->GetNofPoints(); // assign an index of point in internal container
if (aPoint.GetTrackId() > -1) { fpMCData->GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); }
fpMCData->AddPoint(aPoint);
- ++fvNofPointsUsed[int(DetID)];
}
} // iP: end
} // key: end
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index 64dc3dd247..85e0ec20d3 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -39,11 +39,7 @@ TimeSliceReader::TimeSliceReader(ECbmTrackingMode mode) : fTrackingMode(mode) {}
void TimeSliceReader::Clear()
{
// Detector used
- fbUseMvd = false;
- fbUseSts = false;
- fbUseMuch = false;
- fbUseTrd = false;
- fbUseTof = false;
+ fvbUseDet.fill(false);
// Input branches
fpBrTimeSlice = nullptr;
@@ -84,11 +80,11 @@ void TimeSliceReader::CheckInit() const
if (!fpBrTimeSlice) { throw std::logic_error("Time slice branch is unavailable"); }
- if (fbUseMvd && !fpBrMvdHits) { throw std::logic_error("MVD hits branch is unavailable"); }
- if (fbUseSts && !fpBrStsHits) { throw std::logic_error("STS hits branch is unavailable"); }
- if (fbUseMuch && !fpBrMuchHits) { throw std::logic_error("MuCh hits branch is unavailable"); }
- if (fbUseTrd && !fpBrTrdHits) { throw std::logic_error("TRD hits branch is unavailable"); }
- if (fbUseTof && !fpBrTofHits) { throw std::logic_error("TOF hits branch is unavailable"); }
+ if (fvbUseDet[L1DetectorID::kMvd] && !fpBrMvdHits) { throw std::logic_error("MVD hits branch is unavailable"); }
+ if (fvbUseDet[L1DetectorID::kSts] && !fpBrStsHits) { throw std::logic_error("STS hits branch is unavailable"); }
+ if (fvbUseDet[L1DetectorID::kMuch] && !fpBrMuchHits) { throw std::logic_error("MuCh hits branch is unavailable"); }
+ if (fvbUseDet[L1DetectorID::kTrd] && !fpBrTrdHits) { throw std::logic_error("TRD hits branch is unavailable"); }
+ if (fvbUseDet[L1DetectorID::kTof] && !fpBrTofHits) { throw std::logic_error("TOF hits branch is unavailable"); }
if (fpvTracks) {
if (ECbmTrackingMode::kSTS == fTrackingMode) {
@@ -96,10 +92,12 @@ void TimeSliceReader::CheckInit() const
}
else if (ECbmTrackingMode::kMCBM == fTrackingMode) {
if (!fpBrRecoTracks) { throw std::logic_error("GlobalTrack branch is unavailable"); }
- if (fbUseSts && !fpBrStsTracks) { throw std::logic_error("StsTrack branch is unavailable"); }
- if (fbUseMuch && !fpBrMuchTracks) { throw std::logic_error("MuchTrack branch is unavailable"); }
- if (fbUseTrd && !fpBrRecoTracks) { throw std::logic_error("TrdTrack branch is unavailable"); }
- if (fbUseTof && !fpBrRecoTracks) { throw std::logic_error("TofTrack branch is unavailable"); }
+ if (fvbUseDet[L1DetectorID::kSts] && !fpBrStsTracks) { throw std::logic_error("StsTrack branch is not found"); }
+ if (fvbUseDet[L1DetectorID::kMuch] && !fpBrMuchTracks) {
+ throw std::logic_error("MuchTrack branch is not found");
+ }
+ if (fvbUseDet[L1DetectorID::kTrd] && !fpBrRecoTracks) { throw std::logic_error("TrdTrack branch is not found"); }
+ if (fvbUseDet[L1DetectorID::kTof] && !fpBrRecoTracks) { throw std::logic_error("TofTrack branch is not found"); }
}
}
}
@@ -125,11 +123,13 @@ try {
fpBrTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairManager->GetObject("TimeSlice."));
// Init hit branches
- if (fbUseMvd) { fpBrMvdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MvdHit")); }
- if (fbUseSts) { fpBrStsHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsHit")); }
- if (fbUseMuch) { fpBrMuchHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchPixelHit")); }
- if (fbUseTrd) { fpBrTrdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdHit")); }
- if (fbUseTof) { fpBrTofHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofHit")); }
+ if (fvbUseDet[L1DetectorID::kMvd]) { fpBrMvdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MvdHit")); }
+ if (fvbUseDet[L1DetectorID::kSts]) { fpBrStsHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsHit")); }
+ if (fvbUseDet[L1DetectorID::kMuch]) {
+ fpBrMuchHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchPixelHit"));
+ }
+ if (fvbUseDet[L1DetectorID::kTrd]) { fpBrTrdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdHit")); }
+ if (fvbUseDet[L1DetectorID::kTof]) { fpBrTofHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofHit")); }
// Init track branches
if (fpvTracks) {
@@ -139,10 +139,18 @@ try {
break;
case ECbmTrackingMode::kMCBM:
fpBrRecoTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("GlobalTrack"));
- if (fbUseSts) { fpBrStsTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsTrack")); }
- if (fbUseMuch) { fpBrMuchTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchTrack")); }
- if (fbUseTrd) { fpBrTrdTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdTrack")); }
- if (fbUseTof) { fpBrTofTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofTrack")); }
+ if (fvbUseDet[L1DetectorID::kSts]) {
+ fpBrStsTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsTrack"));
+ }
+ if (fvbUseDet[L1DetectorID::kMuch]) {
+ fpBrMuchTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchTrack"));
+ }
+ if (fvbUseDet[L1DetectorID::kTrd]) {
+ fpBrTrdTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdTrack"));
+ }
+ if (fvbUseDet[L1DetectorID::kTof]) {
+ fpBrTofTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofTrack"));
+ }
break;
}
}
@@ -223,11 +231,11 @@ void TimeSliceReader::ReadRecoTracks()
// ** Fill information from local tracks **
// STS tracks (+ MVD)
- if (fbUseSts) {
+ if (fvbUseDet[L1DetectorID::kSts]) {
int iStsTrkId = pInputTrack->GetStsTrackIndex();
if (iStsTrkId > -1) {
auto* pStsTrack = static_cast<CbmStsTrack*>(fpBrStsTracks->At(iStsTrkId));
- if (fbUseMvd) {
+ if (fvbUseDet[L1DetectorID::kMvd]) {
for (int iH = 0; iH < pStsTrack->GetNofMvdHits(); ++iH) {
int iHext = pStsTrack->GetMvdHitIndex(iH);
int iHint = vHitMvdIds[iHext];
@@ -243,7 +251,7 @@ void TimeSliceReader::ReadRecoTracks()
}
// MUCH tracks
- if (fbUseMuch) {
+ if (fvbUseDet[L1DetectorID::kMuch]) {
int iMuchTrkId = pInputTrack->GetMuchTrackIndex();
if (iMuchTrkId > -1) {
auto* pMuchTrack = static_cast<CbmMuchTrack*>(fpBrMuchTracks->At(iMuchTrkId));
@@ -256,7 +264,7 @@ void TimeSliceReader::ReadRecoTracks()
}
// TRD tracks
- if (fbUseTrd) {
+ if (fvbUseDet[L1DetectorID::kTrd]) {
int iTrdTrkId = pInputTrack->GetTrdTrackIndex();
if (iTrdTrkId > -1) {
const auto* pTrdTrack = static_cast<const CbmTrdTrack*>(fpBrTrdTracks->At(iTrdTrkId));
@@ -269,7 +277,7 @@ void TimeSliceReader::ReadRecoTracks()
}
// TOF tracks
- if (fbUseTof) {
+ if (fvbUseDet[L1DetectorID::kTof]) {
int iTofTrkId = pInputTrack->GetTofTrackIndex();
if (iTofTrkId > -1) {
const auto* pTofTrack = static_cast<const CbmTofTrack*>(fpBrTofTracks->At(iTofTrkId));
@@ -279,7 +287,7 @@ void TimeSliceReader::ReadRecoTracks()
track.Hits.push_back(iHint);
} // iH
} // if iTofTrkId > -1
- } // if fbUseTof
+ } // if fvbUseDet[L1DetectorID::kTof]
} // iT
break;
}
@@ -293,19 +301,6 @@ void TimeSliceReader::RegisterIODataManager(std::shared_ptr<L1IODataManager>& pI
fpIODataManager = pIODataManager;
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void TimeSliceReader::SetDetector(L1DetectorID detID, bool flag)
-{
- switch (detID) {
- case L1DetectorID::kMvd: fbUseMvd = flag; break;
- case L1DetectorID::kSts: fbUseSts = flag; break;
- case L1DetectorID::kMuch: fbUseMuch = flag; break;
- case L1DetectorID::kTrd: fbUseTrd = flag; break;
- case L1DetectorID::kTof: fbUseTof = flag; break;
- }
-}
-
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
@@ -352,15 +347,13 @@ void TimeSliceReader::ReadHits()
fNofHitKeys = 0;
fFirstHitKey = 0;
- // Get total number of hits
- std::array<int, L1Constants::size::kMaxNdetectors> nHitsDet = {0};
- if (fbUseMvd) { nHitsDet[int(L1DetectorID::kMvd)] = fpBrMvdHits->GetEntriesFast(); }
- if (fbUseSts) { nHitsDet[int(L1DetectorID::kSts)] = fpBrStsHits->GetEntriesFast(); }
- if (fbUseMuch) { nHitsDet[int(L1DetectorID::kMuch)] = fpBrMuchHits->GetEntriesFast(); }
- if (fbUseTrd) { nHitsDet[int(L1DetectorID::kTrd)] = fpBrTrdHits->GetEntriesFast(); }
- if (fbUseTof) { nHitsDet[int(L1DetectorID::kTof)] = fpBrTofHits->GetEntriesFast(); }
+ if (fvbUseDet[L1DetectorID::kMvd]) { fvNofHitsTotal[L1DetectorID::kMvd] = fpBrMvdHits->GetEntriesFast(); }
+ if (fvbUseDet[L1DetectorID::kSts]) { fvNofHitsTotal[L1DetectorID::kSts] = fpBrStsHits->GetEntriesFast(); }
+ if (fvbUseDet[L1DetectorID::kMuch]) { fvNofHitsTotal[L1DetectorID::kMuch] = fpBrMuchHits->GetEntriesFast(); }
+ if (fvbUseDet[L1DetectorID::kTrd]) { fvNofHitsTotal[L1DetectorID::kTrd] = fpBrTrdHits->GetEntriesFast(); }
+ if (fvbUseDet[L1DetectorID::kTof]) { fvNofHitsTotal[L1DetectorID::kTof] = fpBrTofHits->GetEntriesFast(); }
- int nHitsTot = std::accumulate(nHitsDet.begin(), nHitsDet.end(), 0);
+ int nHitsTot = std::accumulate(fvNofHitsTotal.begin(), fvNofHitsTotal.end(), 0);
// Resize the containers
if (fpvHitIds) {
@@ -378,22 +371,30 @@ void TimeSliceReader::ReadHits()
std::fill(fvHitFirstIndexDet.begin(), fvHitFirstIndexDet.end(), 0);
- // Save number of hits stored
- std::array<int, L1Constants::size::kMaxNdetectors> vNofHits = {0};
-
// Read hits for different detectors
- if (fbUseMvd) { vNofHits[int(L1DetectorID::kMvd)] += ReadHitsForDetector<L1DetectorID::kMvd>(fpBrMvdHits); }
- if (fbUseSts) { vNofHits[int(L1DetectorID::kSts)] += ReadHitsForDetector<L1DetectorID::kSts>(fpBrStsHits); }
- if (fbUseMuch) { vNofHits[int(L1DetectorID::kMuch)] += ReadHitsForDetector<L1DetectorID::kMuch>(fpBrMuchHits); }
- if (fbUseTrd) { vNofHits[int(L1DetectorID::kTrd)] += ReadHitsForDetector<L1DetectorID::kTrd>(fpBrTrdHits); }
- if (fbUseTof) { vNofHits[int(L1DetectorID::kTof)] += ReadHitsForDetector<L1DetectorID::kTof>(fpBrTofHits); }
+ if (fvbUseDet[L1DetectorID::kMvd]) {
+ fvNofHitsUsed[L1DetectorID::kMvd] = ReadHitsForDetector<L1DetectorID::kMvd>(fpBrMvdHits);
+ }
+ if (fvbUseDet[L1DetectorID::kSts]) {
+ fvNofHitsUsed[L1DetectorID::kSts] = ReadHitsForDetector<L1DetectorID::kSts>(fpBrStsHits);
+ }
+ if (fvbUseDet[L1DetectorID::kMuch]) {
+ fvNofHitsUsed[L1DetectorID::kMuch] = ReadHitsForDetector<L1DetectorID::kMuch>(fpBrMuchHits);
+ }
+ if (fvbUseDet[L1DetectorID::kTrd]) {
+ fvNofHitsUsed[L1DetectorID::kTrd] = ReadHitsForDetector<L1DetectorID::kTrd>(fpBrTrdHits);
+ }
+ if (fvbUseDet[L1DetectorID::kTof]) {
+ fvNofHitsUsed[L1DetectorID::kTof] = ReadHitsForDetector<L1DetectorID::kTof>(fpBrTofHits);
+ }
// Save first hit index for different detector subsystems
- for (uint32_t iDet = 0; iDet < vNofHits.size(); ++iDet) {
- fvHitFirstIndexDet[iDet + 1] = fvHitFirstIndexDet[iDet] + vNofHits[iDet];
- fNofHits += vNofHits[iDet];
+ for (uint32_t iDet = 0; iDet < fvNofHitsUsed.size(); ++iDet) {
+ fvHitFirstIndexDet[iDet + 1] = fvHitFirstIndexDet[iDet] + fvNofHitsUsed[iDet];
}
+ fNofHits = std::accumulate(fvNofHitsUsed.cbegin(), fvNofHitsUsed.cend(), 0);
+
// Update number of hit keys in input data object
if (fpIODataManager) { fpIODataManager->SetNhitKeys(fNofHitKeys); }
@@ -405,16 +406,19 @@ void TimeSliceReader::ReadHits()
if (fpvHitIds) {
auto UpdateHitIndexMap = [&](L1Vector<int> & vIds, L1DetectorID detID) constexpr
{
- vIds.reset(nHitsDet[int(detID)]);
- for (int iH = fvHitFirstIndexDet[int(detID)]; iH < fvHitFirstIndexDet[int(detID) + 1]; ++iH) {
- vIds[(*fpvQaHits)[iH].ExtIndex] = iH;
+ if (fvbUseDet[detID]) {
+ vIds.reset(fvNofHitsTotal[detID]);
+ for (int iH = fvHitFirstIndexDet[int(detID)]; iH < fvHitFirstIndexDet[int(detID) + 1]; ++iH) {
+ vIds[(*fpvQaHits)[iH].ExtIndex] = iH;
+ }
}
};
- if (fbUseMvd) { UpdateHitIndexMap(vHitMvdIds, L1DetectorID::kMvd); }
- if (fbUseSts) { UpdateHitIndexMap(vHitStsIds, L1DetectorID::kSts); }
- if (fbUseMuch) { UpdateHitIndexMap(vHitMuchIds, L1DetectorID::kMuch); }
- if (fbUseTrd) { UpdateHitIndexMap(vHitTrdIds, L1DetectorID::kTrd); }
- if (fbUseTof) { UpdateHitIndexMap(vHitTofIds, L1DetectorID::kTof); }
+
+ UpdateHitIndexMap(vHitMvdIds, L1DetectorID::kMvd);
+ UpdateHitIndexMap(vHitStsIds, L1DetectorID::kSts);
+ UpdateHitIndexMap(vHitMuchIds, L1DetectorID::kMuch);
+ UpdateHitIndexMap(vHitTrdIds, L1DetectorID::kTrd);
+ UpdateHitIndexMap(vHitTofIds, L1DetectorID::kTof);
}
}
diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h
index 32668cab1f..1b9e60dd05 100644
--- a/reco/L1/CbmCaTimeSliceReader.h
+++ b/reco/L1/CbmCaTimeSliceReader.h
@@ -151,7 +151,7 @@ namespace cbm::ca
/// @param detID Id of detector
/// @param flag Flag: true - detector is used
/// @note Should be called before this->Init()
- void SetDetector(L1DetectorID detID, bool flag = true);
+ void SetDetector(L1DetectorID detID, bool flag = true) { fvbUseDet[detID] = flag; }
private:
/// @brief Reads hits for a given detector subsystem
@@ -173,14 +173,9 @@ namespace cbm::ca
void StoreHitRecord(const HitRecord& hitRecord);
- // Flags for detector subsystems being used
- bool fbUseMvd = false; ///< is MVD used
- bool fbUseSts = false; ///< is STS used
- bool fbUseMuch = false; ///< is MuCh used
- bool fbUseTrd = false; ///< is TRD used
- bool fbUseTof = false; ///< is TOF used
bool fbReadTracks = true; ///< flag to read reconstructed tracks from reco.root
+
// Pointers to the tracking detector interfaces
CbmTrackingDetectorInterfaceBase* fpMvdInterface = nullptr;
CbmTrackingDetectorInterfaceBase* fpStsInterface = nullptr;
@@ -197,6 +192,8 @@ namespace cbm::ca
TClonesArray* fpBrTrdHits = nullptr; ///< Input branch for TRD hits ("TrdHit")
TClonesArray* fpBrTofHits = nullptr; ///< Input branch for TOF hits ("TofHit")
+ CbmCaDetIdArr_t<TClonesArray*> fvpBrHits = {nullptr}; ///< Input branch for hits
+
// Branches for reconstructed tracks. The input at the moment (as for 27.02.2023) depends on the selected
// tracking mode. For simulations in CBM, the CA tracking is used only in STS + MVD detectors. In this case
// the reconstructed tracks are saved to the "StsTrack" branch as CbmStsTrack objects. For mCBM, the tracks from
@@ -207,7 +204,6 @@ namespace cbm::ca
TClonesArray* fpBrTrdTracks = nullptr; ///< Input branch for reconstructed TRD tracks ("TrdTrack")
TClonesArray* fpBrTofTracks = nullptr; ///< Input branch for reconstructed TOF tracks ("TofTrack")
-
// Pointers to output data containers
L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to array of hit index objects
L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to array of debug hits
@@ -223,6 +219,9 @@ namespace cbm::ca
L1Vector<int> vHitTrdIds {"CbmCaTimeSliceReader::vHitTrdIds"};
L1Vector<int> vHitTofIds {"CbmCaTimeSliceReader::vHitTofIds"};
+ CbmCaDetIdArr_t<int> fvNofHitsTotal = {0}; ///< Total hit number in detector
+ CbmCaDetIdArr_t<int> fvNofHitsUsed = {0}; ///< Number of used hits in detector
+ CbmCaDetIdArr_t<bool> fvbUseDet = {false}; ///< Flag: is detector subsystem used
// Additional
ECbmTrackingMode fTrackingMode; ///< Tracking mode
@@ -232,8 +231,7 @@ namespace cbm::ca
int fNofHitKeys = 0; ///< Recorded number of hit keys
int fFirstHitKey = 0; ///< First index of hit key for the detector subsystem
- std::array<int, L1Constants::size::kMaxNdetectors + 1> fvHitFirstIndexDet = {
- 0}; ///< First index of hit in detector
+ std::array<int, L1Constants::size::kMaxNdetectors + 1> fvHitFirstIndexDet = {0}; ///< First hit index in detector
};
} // namespace cbm::ca
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 9c6633e1f6..c7d8f1de20 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -38,6 +38,7 @@
#include "FairDetector.h"
#include "FairRootManager.h"
#include "FairTask.h"
+#include "Logger.h"
#include "TClonesArray.h"
#include "TH1.h"
@@ -269,8 +270,8 @@ public:
case L1DetectorID::kMuch: return "MuCh";
case L1DetectorID::kTrd: return "TRD";
case L1DetectorID::kTof: return "TOF";
+ case L1DetectorID::kEND: break;
}
- // TODO: Probably, we should provide default with throwing exception here... (S.Zharko)
return "";
}
diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h
index c63df00909..1ec1ae7eb2 100644
--- a/reco/L1/CbmL1DetectorID.h
+++ b/reco/L1/CbmL1DetectorID.h
@@ -2,14 +2,16 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergei Zharko [committer] */
-/// \file CbmL1DetectorID.h
-/// \brief Implementation of L1DetectorID enum class for CBM
-/// \author S.Zharko
-/// \data 01.12.2022
+/// @file CbmL1DetectorID.h
+/// @brief Implementation of L1DetectorID enum class for CBM
+/// @author S.Zharko
+/// @since 01.12.2022
#ifndef CbmL1DetectorID_h
#define CbmL1DetectorID_h 1
+#include "L1EArray.h"
+
/// Enumeration for the detector subsystems used in L1 tracking
/// It is important for the subsystems to be specified in the actual order. The order is used
/// for the L1Station array filling.
@@ -20,9 +22,16 @@ enum class L1DetectorID
kSts,
kMuch,
kTrd,
- kTof
+ kTof,
+ kEND ///< End of enumeration
};
+/// @brief Alias to array, indexed by L1DetectorID enum
+/// @note To be used only in CBM-specific code
+template<typename T>
+using CbmCaDetIdArr_t = L1EArray<L1DetectorID, T>;
+
+
/// @brief Enumeration for different tracking running modes
enum class ECbmTrackingMode
{
diff --git a/reco/L1/L1Algo/L1EArray.h b/reco/L1/L1Algo/L1EArray.h
new file mode 100644
index 0000000000..272238f731
--- /dev/null
+++ b/reco/L1/L1Algo/L1EArray.h
@@ -0,0 +1,43 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file L1EArray.h
+/// @brief Implementation of L1EArray class
+/// @since 02.05.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#ifndef L1EArray_h
+#define L1EArray_h 1
+
+#include <array>
+
+/// @brief Class of arrays, which can be accessed by an enum class entry as an index
+/// @note The enum-array must contain an entry kEND, which represents the number of the enumeration entries and
+/// is used, as an array size
+/// @tparam E The enum class type
+/// @tparam T Type of data in the underlying array
+template<class E, class T>
+class L1EArray : public std::array<T, static_cast<std::size_t>(E::kEND)> {
+ using U = typename std::underlying_type<E>::type; ///< Underlying type of enumeration
+public:
+ /// @brief Mutable access operator, indexed by enum entry
+ T& operator[](const E& entry)
+ {
+ return std::array<T, static_cast<std::size_t>(E::kEND)>::operator[](static_cast<U>(entry));
+ }
+
+ /// @brief Mutable access operator, indexed by underlying index type
+ T& operator[](U index) { return std::array<T, static_cast<std::size_t>(E::kEND)>::operator[](index); }
+
+ /// @brief Constant access operator, indexed by enum entry
+ const T& operator[](const E& entry) const
+ {
+ return std::array<T, static_cast<std::size_t>(E::kEND)>::operator[](static_cast<U>(entry));
+ }
+
+ /// @brief Constant access operator, indexed by underlying index type
+ const T& operator[](U index) const { return std::array<T, static_cast<std::size_t>(E::kEND)>::operator[](index); }
+};
+
+#endif // L1EArray_h
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index af9eca6e0e..46e2ec7e73 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -173,10 +173,10 @@ inline fscal L1TrackPar::GetPhiErr() const
fscal phiDFactor = 1. / (GetTx() * GetTx() + GetTy() * GetTy());
fscal phiDTx = -phiDFactor * GetTy(); // partial derivative of phi over Tx
fscal phiDTy = +phiDFactor * GetTx(); // partial derivative of phi over Ty
-
+
fscal varTx = (std::isfinite(C22[0]) && C22[0] > 0) ? C22[0] : undef::kF; // variance of Tx
fscal varTy = (std::isfinite(C33[0]) && C33[0] > 0) ? C33[0] : undef::kF; // variance of Ty
- fscal covTxTy = std::isfinite(C32[0]) ? C32[0] : undef::kF; // covariance of Tx and Ty
+ fscal covTxTy = std::isfinite(C32[0]) ? C32[0] : undef::kF; // covariance of Tx and Ty
fscal varPhi = phiDTx * phiDTx * varTx + phiDTy * phiDTy * varTy + 2 * phiDTx * phiDTy * covTxTy;
return std::sqrt(varPhi);
@@ -186,14 +186,14 @@ inline fscal L1TrackPar::GetPhiErr() const
//
inline fscal L1TrackPar::GetThetaErr() const
{
- fscal sumSqSlopes = GetTx() * GetTx() + GetTy() * GetTy();
+ fscal sumSqSlopes = GetTx() * GetTx() + GetTy() * GetTy();
fscal thetaDFactor = 1. / ((sumSqSlopes + 1) * sqrt(sumSqSlopes));
fscal thetaDTx = thetaDFactor * GetTx();
fscal thetaDTy = thetaDFactor * GetTy();
fscal varTx = (std::isfinite(C22[0]) && C22[0] > 0) ? C22[0] : undef::kF; // variance of Tx
fscal varTy = (std::isfinite(C33[0]) && C33[0] > 0) ? C33[0] : undef::kF; // variance of Ty
- fscal covTxTy = std::isfinite(C32[0]) ? C32[0] : undef::kF; // covariance of Tx and Ty
+ fscal covTxTy = std::isfinite(C32[0]) ? C32[0] : undef::kF; // covariance of Tx and Ty
fscal varTheta = thetaDTx * thetaDTx * varTx + thetaDTy * thetaDTy * varTy + 2 * thetaDTx * thetaDTy * covTxTy;
return std::sqrt(varTheta);
diff --git a/reco/L1/catools/CaToolsMCData.cxx b/reco/L1/catools/CaToolsMCData.cxx
index 33213add26..70d56bf99c 100644
--- a/reco/L1/catools/CaToolsMCData.cxx
+++ b/reco/L1/catools/CaToolsMCData.cxx
@@ -26,6 +26,8 @@ MCData::MCData() {}
MCData::MCData(const MCData& other)
: fvPoints(other.fvPoints)
, fvTracks(other.fvTracks)
+ , fvNofPointsOrig(other.fvNofPointsOrig)
+ , fvNofPointsUsed(other.fvNofPointsUsed)
, fmPointLinkMap(other.fmPointLinkMap)
, fmTrackLinkMap(other.fmTrackLinkMap)
{
@@ -60,25 +62,12 @@ void MCData::Swap(MCData& other) noexcept
{
std::swap(fvPoints, other.fvPoints);
std::swap(fvTracks, other.fvTracks);
+ std::swap(fvNofPointsOrig, other.fvNofPointsOrig);
+ std::swap(fvNofPointsUsed, other.fvNofPointsUsed);
std::swap(fmPointLinkMap, other.fmPointLinkMap);
std::swap(fmTrackLinkMap, other.fmTrackLinkMap);
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void MCData::AddPoint(const MCPoint& point)
-{
- fmPointLinkMap[point.GetLinkKey()] = point.GetId();
- fvPoints.push_back(point);
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void MCData::AddTrack(const MCTrack& track)
-{
- fmTrackLinkMap[track.GetLinkKey()] = track.GetId();
- fvTracks.push_back(track);
-}
// ---------------------------------------------------------------------------------------------------------------------
//
diff --git a/reco/L1/catools/CaToolsMCData.h b/reco/L1/catools/CaToolsMCData.h
index d0f6ec431d..72c3d95e07 100644
--- a/reco/L1/catools/CaToolsMCData.h
+++ b/reco/L1/catools/CaToolsMCData.h
@@ -17,6 +17,7 @@
#include "CaToolsMCPoint.h"
#include "CaToolsMCTrack.h"
#include "L1Constants.h"
+#include "L1EArray.h"
#include "L1Vector.h"
enum class L1DetectorID;
@@ -51,10 +52,6 @@ namespace ca::tools
/// Swap method
void Swap(MCData& other) noexcept;
- /// Adds hit index to MC point
- /// \param iPoint Index of MC point
- /// \param
-
/// Adds an MC point to points container and a corresponding link key to the point index map
/// \param point MC point object
void AddPoint(const MCPoint& point);
@@ -66,23 +63,25 @@ namespace ca::tools
/// Clears contents
void Clear();
- /// Finds an index of MC point in internal point container
- /// \param index Index of MC point in external point container
- /// \param event Index of MC event
- /// \param file Index of MC file
- /// \return Index of MC point in internal point container within event/TS
+ /// @brief Finds an index of MC point in internal point container
+ /// @param detID Detector ID
+ /// @param index Index of MC point in external point container
+ /// @param event Index of MC event
+ /// @param file Index of MC file
+ /// @return Index of MC point in internal point container within event/TS
/// If the point is not found the function returns -1
- int FindInternalPointIndex(int index, int event, int file) const
+ int FindInternalPointIndex(L1DetectorID detID, int index, int event, int file) const
{
- auto it = fmPointLinkMap.find(LinkKey(index, event, file));
+ int indexGlob = GetPointGlobExtIndex(detID, index);
+ auto it = fmPointLinkMap.find(LinkKey(indexGlob, event, file));
return (it != fmPointLinkMap.cend()) ? it->second : -1;
}
- /// Finds an index of MC track in internal track container
- /// \param index Index of MC track in external track container
- /// \param event Index of MC event
- /// \param file Index of MC file
- /// \return Index of MC track in internal track container within event/TS
+ /// @brief Finds an index of MC track in internal track container
+ /// @param index Index of MC track in external track container
+ /// @param event Index of MC event
+ /// @param file Index of MC file
+ /// @return Index of MC track in internal track container within event/TS
/// If the track is not found, the function returns -1
int FindInternalTrackIndex(int index, int event, int file) const
{
@@ -90,12 +89,46 @@ namespace ca::tools
return (it != fmTrackLinkMap.cend()) ? it->second : -1;
}
+ /// @brief Gets the first point index for a given detector subsystem
+ /// @param detID Detector ID
+ int GetFirstPointIndex(L1DetectorID detID) const
+ {
+ return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + static_cast<int>(detID), 0);
+ }
+
+ /// @brief Gets the first point index for a given detector subsystem
+ /// @param detID Detector ID
+ int GetLastPointIndex(L1DetectorID detID) const
+ {
+ return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + static_cast<int>(detID) + 1, 0) - 1;
+ }
+
+ /// @brief Calculates global index of MC point
+ /// @param iPointLocal Local index of MC point
+ /// @param detID Detector ID
+ ///
+ /// The function calculates global external index of MC point as a sum of a given local index and total provided
+ /// number of points in previous detector subsystem.
+ int GetPointGlobExtIndex(L1DetectorID detID, int iPointLocal) const
+ {
+ return iPointLocal + std::accumulate(fvNofPointsOrig.cbegin(), fvNofPointsOrig.cbegin() + int(detID), 0);
+ }
+
+
/// Gets number of tracks in this event/TS
int GetNofTracks() const { return fvTracks.size(); }
/// Gets number of points in this event/TS
int GetNofPoints() const { return fvPoints.size(); }
+ /// @brief Gets original number of MC points in different detectors
+ /// @param detID Detector ID
+ int GetNofPointsOrig(L1DetectorID detID) const { return fvNofPointsOrig[static_cast<int>(detID)]; }
+
+ /// @brief Gets used number of MC points in different detectors
+ /// @param detID Detector ID
+ int GetNofPointsUsed(L1DetectorID detID) const { return fvNofPointsUsed[static_cast<int>(detID)]; }
+
/// Gets a reference to MC point by its index
const auto& GetPoint(int idx) const { return fvPoints[idx]; }
@@ -132,6 +165,11 @@ namespace ca::tools
/// Reserves memory for points to avoid extra allocations
void ReserveNofPoints(int nPoints) { fvPoints.reserve(nPoints); }
+ /// @brief Sets original number of MC points in different detectors
+ /// @param detID Detector ID
+ /// @param nPoints Number of points
+ void SetNofPointsOrig(L1DetectorID detID, int nPoints) { fvNofPointsOrig[static_cast<int>(detID)] = nPoints; }
+
/// Prints an example of tracks and points
/// \param verbose Verbose level:
/// - #0: Nothing is printed
@@ -147,9 +185,34 @@ namespace ca::tools
L1Vector<MCPoint> fvPoints = {"ca::tools::MCData::fvPoints"}; ///< Container of points
L1Vector<MCTrack> fvTracks = {"ca::tools::MCData::fvTracks"}; ///< Container of tracks
+ std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsOrig = {0}; ///< Total number of points by detector
+ std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsUsed = {0}; ///< Number of points used vs. detector
+
std::unordered_map<LinkKey, int> fmPointLinkMap = {}; ///< MC point internal index vs. link
std::unordered_map<LinkKey, int> fmTrackLinkMap = {}; ///< MC track internal index vs. link
};
} // namespace ca::tools
+
+// *********************************************************
+// ** Template and inline function implementation **
+// *********************************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+inline void ca::tools::MCData::AddPoint(const MCPoint& point)
+{
+ fmPointLinkMap[point.GetLinkKey()] = point.GetId();
+ fvPoints.push_back(point);
+ ++fvNofPointsUsed[static_cast<int>(point.GetDetectorId())];
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+inline void ca::tools::MCData::AddTrack(const MCTrack& track)
+{
+ fmTrackLinkMap[track.GetLinkKey()] = track.GetId();
+ fvTracks.push_back(track);
+}
+
#endif // CaToolsMCData_h
diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx
index 181950dd31..40ad8b5e7c 100644
--- a/reco/L1/qa/CbmCaInputQaSts.cxx
+++ b/reco/L1/qa/CbmCaInputQaSts.cxx
@@ -145,18 +145,24 @@ bool CbmCaInputQaSts::Check()
// Function to fit a residual distribution, returns a structure
auto FitResidualsAndGetMean = [&](TH1* pHist) {
- auto fit = TF1("fitres", "gausn");
- double statMean = pHist->GetMean();
- double statSigm = pHist->GetStdDev();
- fit.SetParameters(100., statMean, statSigm);
- pHist->Fit("fitres", "MQ");
- pHist->Fit("fitres", "MQ");
- pHist->Fit("fitres", "MQ");
+ // Fit function - Kaniadakis Gaussian Distribution:
+ // [0] - Scale (correlated with [3])
+ // [1] - Mean
+ // [2] - Sigma (biased from standard deviation)
+ // [3] - Kaniadakis parameter (0 - 2)
+ auto fit =
+ TF1("fitres", "[0]*TMath::Exp(TMath::ASinH(-0.5*[3]*((x-[1])/[2])**2)/[3])/([2]*TMath::Sqrt(2*TMath::Pi()))",
+ -10., 10.);
+ fit.SetParameters(100, 0., 1., .3);
+ fit.SetParLimits(3, 0, 2);
+ pHist->Fit("fitres", "Q");
// NOTE: Several fit procedures are made to avoid empty fit results
std::array<double, 3> result;
result[0] = fit.GetParameter(1);
- result[1] = -fit.GetParameter(2) * fResMeanThrsh;
- result[2] = +fit.GetParameter(2) * fResMeanThrsh;
+ double fitK = fit.GetParameter(3);
+ double fitHWHM = fit.GetParameter(2) * std::sqrt((1 - std::pow(2, -2 * fitK)) / fitK);
+ result[1] = -fitHWHM * fResMeanThrsh;
+ result[2] = +fitHWHM * fResMeanThrsh;
return result;
};
@@ -190,17 +196,19 @@ bool CbmCaInputQaSts::Check()
// Fit pull distributions
- // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF
- //if (!gROOT->FindObject("Expk")) {
- // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])");
- //}
- //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
-
auto FitPullsAndGetSigma = [&](TH1* pHist) {
- auto fit = TF1("fitpull", "gausn(0)");
- fit.SetParameters(100, 0.001, 1.000);
+ // Fit function - Kaniadakis Gaussian Distribution:
+ // [0] - Scale (correlated with [3])
+ // [1] - Mean
+ // [2] - Sigma (biased from standard deviation)
+ // [3] - Kaniadakis parameter (0 - 2)
+ auto fit = TF1("fitpull", "[0] * TMath::Exp(TMath::ASinH(-0.5*[3]*((x-[1])/[2])**2)/[3])", -10., 10.);
+ fit.SetParameters(100, 0., 1., .3);
+ fit.SetParLimits(3, 0, 2);
pHist->Fit("fitpull", "Q");
- return fit.GetParameter(2);
+ double fitK = fit.GetParameter(3);
+ double fitHWHM = fit.GetParameter(2) * std::sqrt((1 - std::pow(2, -2 * fitK)) / fitK);
+ return fitHWHM;
};
auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index 220f69a167..69f1a33db5 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -9,6 +9,7 @@
#include "CbmCaOutputQa.h"
+#include "CbmCaMCModule.h"
#include "CbmQaCanvas.h"
#include "FairRootManager.h"
@@ -21,6 +22,7 @@
using ca::tools::Debugger;
using ::ca::tools::MCTrack;
+using cbm::ca::MCModule;
using cbm::ca::OutputQa;
// ---------------------------------------------------------------------------------------------------------------------
@@ -377,6 +379,10 @@ InitStatus OutputQa::InitCanvases()
MakeCanvas<CbmQaCanvas>("mc_yMC", "MC reconstructable track MC rapidity", kCXSIZEPX * 3, kCYSIZEPX * 2);
DrawTrackDistributions(pc_mc_yMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_mc_yMC; });
+ // MC rapidity vs. MC momentum
+ auto* pc_mc_pMC_yMC =
+ MakeCanvas<CbmQaCanvas>("mc_pMC_yMC", "MC track MC mom. vs. rapidity ", kCXSIZEPX * 3, kCYSIZEPX * 2);
+ DrawSetOf<TH2F>(vCmpTypesGeneral, [&](ETrackType t) -> TH2F* { return fvpTrackHistograms[t]->fph_reco_pMC_yMC; });
// ** Efficiencies **
@@ -386,6 +392,31 @@ InitStatus OutputQa::InitCanvases()
auto* pc_eff_yMC = MakeCanvas<CbmQaCanvas>("eff_yMC", "Tracking Eff. vs. MC rapidity", kCXSIZEPX * 3, kCYSIZEPX);
DrawTrackEfficiens(pc_eff_yMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_yMC; });
+
+ auto* pc_eff_thetaMC =
+ MakeCanvas<CbmQaCanvas>("eff_thetaMC", "Tracking Eff. vs. MC polar angle", kCXSIZEPX * 3, kCYSIZEPX);
+ DrawTrackEfficiens(pc_eff_thetaMC,
+ [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_thetaMC; });
+
+ auto* pc_eff_phiMC =
+ MakeCanvas<CbmQaCanvas>("eff_phiMC", "Tracking Eff. vs. MC azimuthal angle", kCXSIZEPX * 3, kCYSIZEPX);
+ DrawTrackEfficiens(pc_eff_phiMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_phiMC; });
+
+ auto* pc_eff_etaMC =
+ MakeCanvas<CbmQaCanvas>("eff_etaMC", "Tracking Eff. vs. MC pseudorapidity", kCXSIZEPX * 3, kCYSIZEPX);
+ DrawTrackEfficiens(pc_eff_etaMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_etaMC; });
+
+
+ // ** Pulls and residuals **
+ // NOTE: stored in a subdirectory for a given track type and point type
+ for (int iType = 0; iType < ETrackType::kEND; ++iType) {
+ if (fvbTrackTypeOn[iType] && fvpTrackHistograms[iType]->IsMCUsed()) {
+ fvpTrackHistograms[iType]->fpFitQaFirstHit->CreateResidualPlot();
+ fvpTrackHistograms[iType]->fpFitQaFirstHit->CreatePullPlot();
+ fvpTrackHistograms[iType]->fpFitQaLastHit->CreateResidualPlot();
+ fvpTrackHistograms[iType]->fpFitQaLastHit->CreatePullPlot();
+ }
+ }
}
return kSUCCESS;
@@ -436,7 +467,7 @@ InitStatus OutputQa::InitDataBranches()
// Initialize MC module
if (IsMCUsed()) {
- if (!fpMCModule.get()) { fpMCModule = std::make_unique<CbmCaMCModule>(fVerbose, fPerformanceMode); }
+ if (!fpMCModule.get()) { fpMCModule = std::make_unique<MCModule>(fVerbose, fPerformanceMode); }
fpMCModule->SetDetector(L1DetectorID::kMvd, fbUseMvd);
fpMCModule->SetDetector(L1DetectorID::kSts, fbUseSts);
fpMCModule->SetDetector(L1DetectorID::kMuch, fbUseMuch);
@@ -699,18 +730,39 @@ InitStatus OutputQa::InitTimeSlice()
int nHits = 0;
int nRecoTracks = 0;
+ // Read MC tracks and points
+ if (IsMCUsed()) {
+ fpMCModule->InitEvent(nullptr);
+ nMCPoints = fMCData.GetNofPoints();
+ nMCTracks = fMCData.GetNofTracks();
+ }
+
// Read reconstructed input
fpTSReader->InitTimeSlice();
nHits = fvHits.size();
nRecoTracks = fvRecoTracks.size();
- if (IsMCUsed()) {
- // Read MC information
- fpMCModule->InitEvent(nullptr);
- nMCPoints = fMCData.GetNofPoints();
- nMCTracks = fMCData.GetNofTracks();
+ static bool bDo = true;
+ if (bDo) {
+ for (int iP = 0; iP < nMCPoints; ++iP) {
+ const auto& point = fMCData.GetPoint(iP);
+ LOG(info) << iP << ' ' << (int) point.GetDetectorId() << ' ' << point.GetStationId();
+ }
+
+ for (int iH = 0; iH < nHits; ++iH) {
+ const auto& hit = fvHits[iH];
+ LOG(info) << iH << ' ' << (int) hit.GetDetectorType() << ' ' << hit.GetStationId();
+ }
+ bDo = false;
}
+ for (const auto& hit : fvHits) {
+ assert(hit.GetMCPointIndexes().size() < 2);
+ }
+
+ // Match tracks and points
+ if (IsMCUsed()) { fpMCModule->MatchRecoAndMC(); }
+
LOG_IF(info, fVerbose > 1) << fName << ": Data sample consists of " << nHits << " hits, " << nRecoTracks
<< " reco tracks, " << nMCTracks << " MC tracks, " << nMCPoints << " MC points";
@@ -728,7 +780,7 @@ void OutputQa::ReadParameters(const char* filename)
manager.ReadParametersObject(filename);
manager.SendParameters(*fpParameters);
- LOG(info) << fpParameters->ToString(0);
+ LOG(info) << fpParameters->ToString(5);
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/qa/CbmCaOutputQa.h b/reco/L1/qa/CbmCaOutputQa.h
index 95e8466986..1cc6216224 100644
--- a/reco/L1/qa/CbmCaOutputQa.h
+++ b/reco/L1/qa/CbmCaOutputQa.h
@@ -26,7 +26,7 @@
namespace cbm::ca
{
- /// @brief Enumeration of track category
+ /// @brief Enumeration fors track category
enum ETrackType
{
kAll, ///< all tracks
@@ -235,7 +235,7 @@ namespace cbm::ca
ECbmTrackingMode fTrackingMode = ECbmTrackingMode::kSTS; ///< Tracking mode
std::unique_ptr<TimeSliceReader> fpTSReader = nullptr; ///< Reader of the time slice
- std::unique_ptr<CbmCaMCModule> fpMCModule = nullptr; ///< MC module
+ std::unique_ptr<cbm::ca::MCModule> fpMCModule = nullptr; ///< MC module
std::shared_ptr<L1IODataManager> fpDataManager = nullptr; ///< Data manager
std::shared_ptr<::ca::tools::Debugger> fpDebugger = nullptr; ///< Debugger
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Tracking parameters object
diff --git a/reco/L1/qa/CbmCaTrackFitQa.cxx b/reco/L1/qa/CbmCaTrackFitQa.cxx
index 21b23c0c4c..75bf0ff43c 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.cxx
+++ b/reco/L1/qa/CbmCaTrackFitQa.cxx
@@ -10,12 +10,25 @@
#include "CbmCaTrackFitQa.h"
#include "CbmL1Track.h"
+#include "CbmQaCanvas.h"
+
+#include "TF1.h"
+#include "TFormula.h"
+#include "TH1.h"
+#include "TProfile.h"
+#include "TString.h"
+
+#include <algorithm>
#include "CaToolsMCData.h"
#include "L1Field.h"
#include "L1Fit.h"
-#include "TProfile.h"
-#include "TH1.h"
+#include "L1Utils.h"
+
+
+// *******************************************************
+// ** Implementation of cbm::ca::TrackFitQa functions **
+// *******************************************************
using cbm::ca::TrackFitQa;
@@ -27,52 +40,90 @@ TrackFitQa::TrackFitQa(const char* pointTag, const char* prefixName, TFolder* pP
fStoringMode = EStoringMode::kSAMEDIR;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaCanvas* TrackFitQa::CreateResidualPlot()
+{
+ auto* pCanv = MakeCanvas<CbmQaCanvas>("canv_residuals", " residuals", kCXSIZEPX * 4, kCYSIZEPX * 2);
+ pCanv->Divide2D(7);
+
+
+ for (int iType = static_cast<int>(ETrackParType::kBEGIN); iType < static_cast<int>(ETrackParType::kEND); ++iType) {
+ ETrackParType type = static_cast<ETrackParType>(iType);
+ if (fvbParIgnored[type]) { continue; }
+ pCanv->cd(iType + 1);
+ fvphResiduals[type]->Draw();
+ }
+
+ return pCanv;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaCanvas* TrackFitQa::CreatePullPlot()
+{
+ auto* pCanv = MakeCanvas<CbmQaCanvas>(fsPrefix + "_canv_pulls", " pulls", kCXSIZEPX * 4, kCYSIZEPX * 2);
+ pCanv->Divide2D(7);
+
+ for (int iType = static_cast<int>(ETrackParType::kBEGIN); iType < static_cast<int>(ETrackParType::kEND); ++iType) {
+ ETrackParType type = static_cast<ETrackParType>(iType);
+ if (fvbParIgnored[type]) { continue; }
+ auto fit = TF1("fitpull", "[0] * TMath::Exp(TMath::ASinH(-0.5*[3]*((x-[1])/[2])**2)/[3])", -10., 10.);
+ fit.SetParameters(100, 0., 1., .3);
+ fit.SetParLimits(3, 0., 2.);
+ pCanv->cd(iType + 1);
+ fvphPulls[type]->Draw();
+ fvphPulls[type]->Fit("fitpull", "Q");
+ }
+
+ return pCanv;
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void TrackFitQa::Init()
{
- fph_res_x = MakeHisto<TH1F>("res_x", "", fBinsRESX, fLoRESX, fUpRESX);
- fph_res_y = MakeHisto<TH1F>("res_y", "", fBinsRESY, fLoRESY, fUpRESY);
- fph_res_t = MakeHisto<TH1F>("res_t", "", fBinsREST, fLoREST, fUpREST);
- fph_res_tx = MakeHisto<TH1F>("res_tx", "", fBinsRESTX, fLoRESTX, fUpRESTX);
- fph_res_ty = MakeHisto<TH1F>("res_ty", "", fBinsRESTY, fLoRESTY, fUpRESTY);
- fph_res_qp = MakeHisto<TH1F>("res_qp", "", fBinsRESQP, fLoRESQP, fUpRESQP);
- fph_res_vi = MakeHisto<TH1F>("res_vi", "", fBinsRESVI, fLoRESVI, fUpRESVI);
+ // Init default distribution properties
+ SetDefaultProperties();
+
+ auto CreateResidualHisto = [&](ETrackParType t, const char* name, const char* title) {
+ if (fvbParIgnored[t]) { return; }
+ TString sPrefix = (fsTitle.Length() > 0) ? fsTitle + " point residual for " : "residual for ";
+ fvphResiduals[t] = MakeHisto<TH1F>(name, sPrefix + title, fvRBins[t], fvRLo[t], fvRUp[t]);
+ };
+
+ auto CreatePullHisto = [&](ETrackParType t, const char* name, const char* title) {
+ if (fvbParIgnored[t]) { return; }
+ TString sPrefix = (fsTitle.Length() > 0) ? fsTitle + " point pull for " : "pull for ";
+ fvphPulls[t] = MakeHisto<TH1F>(name, sPrefix + title, fvPBins[t], fvPLo[t], fvPUp[t]);
+ };
+
+ CreateResidualHisto(ETrackParType::kX, "res_x", "x-coordinate;x^{reco} - x^{MC} [cm]");
+ CreateResidualHisto(ETrackParType::kY, "res_y", "y-coordinate;y^{reco} - y^{MC} [cm]");
+ CreateResidualHisto(ETrackParType::kTX, "res_tx", "slope along x-axis;T_{x}^{reco} - T_{x}^{MC}");
+ CreateResidualHisto(ETrackParType::kTY, "res_ty", "slope along y-axis;T_{y}^{reco} - T_{y}^{MC}");
+ CreateResidualHisto(ETrackParType::kQP, "res_qp", "charge over mom.;(q/p)^{reco} - (q/p)^{MC} [ec/GeV]");
+ CreateResidualHisto(ETrackParType::kTIME, "res_t", "time;t^{reco} - t^{MC} [ns]");
+ CreateResidualHisto(ETrackParType::kVI, "res_vi", "inverse speed;v^{-1}_{reco} - v^{-1}_{MC} [c^{-1}]");
+
+ CreatePullHisto(ETrackParType::kX, "pull_x", "x-coordinate;(x^{reco} - x^{MC})/#sigma_{x}");
+ CreatePullHisto(ETrackParType::kY, "pull_y", "y-coordinate;(y^{reco} - y^{MC})/#sigma_{y}");
+ CreatePullHisto(ETrackParType::kTX, "pull_tx", "slope along x-axis;(T_{x}^{reco} - T_{x}^{MC})/#sigma_{T_{x}}");
+ CreatePullHisto(ETrackParType::kTY, "pull_ty", "slope along y-axis;(T_{y}^{reco} - T_{y}^{MC})/#sigma_{T_{y}}");
+ CreatePullHisto(ETrackParType::kQP, "pull_qp", "charge over mom.;((q/p)^{reco} - (q/p)^{MC})/#sigma_{q/p}");
+ CreatePullHisto(ETrackParType::kTIME, "pull_t", "time;(t^{reco} - t^{MC})/#sigma_{t}");
+ CreatePullHisto(ETrackParType::kVI, "pull_vi", "inverse speed;(v^{-1}_{reco} - v^{-1}_{MC})/#sigma_{v^{-1}}");
// FIXME: Replace hardcoded parameters with variables
- fph_res_p_vs_pMC = MakeHisto<TProfile>("res_p_vs_pMC", "", 100, 0.0, 10.0, -2., 2.);
- fph_res_phi_vs_phiMC = MakeHisto<TProfile>("res_phi_vs_phiMC", "", 100, -3.2, 3.2, -2., 2.);
- fph_res_theta_vs_thetaMC = MakeHisto<TProfile>("res_theta_vs_phiMC", "", 100, 0., 3.2, -2., 2.);
-
- fph_pull_x = MakeHisto<TH1F>("pull_x", "", fBinsPULLX, fLoPULLX, fUpPULLX);
- fph_pull_y = MakeHisto<TH1F>("pull_y", "", fBinsPULLY, fLoPULLY, fUpPULLY);
- fph_pull_t = MakeHisto<TH1F>("pull_t", "", fBinsPULLT, fLoPULLT, fUpPULLT);
- fph_pull_tx = MakeHisto<TH1F>("pull_tx", "", fBinsPULLTX, fLoPULLTX, fUpPULLTX);
- fph_pull_ty = MakeHisto<TH1F>("pull_ty", "", fBinsPULLTY, fLoPULLTY, fUpPULLTY);
- fph_pull_qp = MakeHisto<TH1F>("pull_qp", "", fBinsPULLQP, fLoPULLQP, fUpPULLQP);
- fph_pull_vi = MakeHisto<TH1F>("pull_vi", "", fBinsPULLVI, fLoPULLVI, fUpPULLVI);
+ fph_res_p_pMC = MakeHisto<TProfile>("res_p_vs_pMC", "", 100, 0.0, 10.0, -2., 2.);
+ fph_res_phi_phiMC = MakeHisto<TProfile>("res_phi_vs_phiMC", "", 100, -3.2, 3.2, -2., 2.);
+ fph_res_theta_thetaMC = MakeHisto<TProfile>("res_theta_vs_phiMC", "", 100, 0., 3.2, -2., 2.);
// Set histogram titles
TString sPrefix = (fsTitle.Length() > 0) ? fsTitle + " point " : "";
- fph_res_x->SetTitle(sPrefix + " residual for x-coordinate;x_{reco} - x_{MC} [cm]");
- fph_res_y->SetTitle(sPrefix + " residual for y-coordinate;y_{reco} - y_{MC} [cm]");
- fph_res_t->SetTitle(sPrefix + " residual for time;t_{reco} - t_{MC} [ns]");
- fph_res_tx->SetTitle(sPrefix + " residual for slope along x-axis;T_{x}^{reco} - T_{x}^{MC}");
- fph_res_ty->SetTitle(sPrefix + " residual for slope along y-axis;T_{y}^{reco} - T_{y}^{MC}");
- fph_res_qp->SetTitle(sPrefix + " residual for q/p;(q/p)_{reco} - (q/p)_{MC} [ec/GeV]");
- fph_res_vi->SetTitle(sPrefix + " residual for inverse speed;v^{-1}_{reco} - v^{-1}_{MC} [c^{-1}]");
-
- fph_res_p_vs_pMC->SetTitle(sPrefix + " resolution of momentum;p^{MC} [GeV/c];#delta p [GeV/c]");
- fph_res_phi_vs_phiMC->SetTitle(sPrefix + " resolution of polar angle;#phi^{MC} [rad];#delta#phi [rad]");
- fph_res_theta_vs_thetaMC->SetTitle(sPrefix + " resolution of polar angle;#theta^{MC} [rad];#delta#theta [rad]");
-
- fph_pull_x->SetTitle(sPrefix + " pull for x-coordinate;(x_{reco} - x_{MC}) / #sigma_{x}");
- fph_pull_y->SetTitle(sPrefix + " pull for y-coordinate;(y_{reco} - y_{MC}) / #sigma_{y}");
- fph_pull_t->SetTitle(sPrefix + " pull for time;(t_{reco} - t_{MC}) / #sigma_{t}");
- fph_pull_tx->SetTitle(sPrefix + " pull for slope along x-axis;(T_{x}^{reco} - T_{x}^{MC}) / #sigma_{T_{x}}");
- fph_pull_ty->SetTitle(sPrefix + " pull for slope along y-axis;(T_{y}^{reco} - T_{y}^{MC}) / #sigma_{T_{y}}");
- fph_pull_qp->SetTitle(sPrefix + " pull for q/p;((q/p)_{reco} - (q/p)_{MC}) / #sigma_{q/p}");
- fph_pull_vi->SetTitle(sPrefix + " pull for inverse speed;(v^{-1}_{reco} - v^{-1}_{MC}) / #sigma_{v^{-1}}");
+ fph_res_p_pMC->SetTitle(sPrefix + " resolution of momentum;p^{MC} [GeV/c];#delta p [GeV/c]");
+ fph_res_phi_phiMC->SetTitle(sPrefix + " resolution of polar angle;#phi^{MC} [rad];#delta#phi [rad]");
+ fph_res_theta_thetaMC->SetTitle(sPrefix + " resolution of polar angle;#theta^{MC} [rad];#delta#theta [rad]");
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -92,52 +143,97 @@ void TrackFitQa::Fill(const L1TrackPar& trPar, const ::ca::tools::MCPoint& mcPoi
const L1TrackPar& trParExtr = fitter.Tr(); // Track parameters extrapolated to given MC point
// ** Time-independent measurements **
+ FillResAndPull(ETrackParType::kX, trParExtr.GetX(), trParExtr.GetXErr(), mcPoint.GetXOut(), weight);
+ FillResAndPull(ETrackParType::kY, trParExtr.GetY(), trParExtr.GetYErr(), mcPoint.GetYOut(), weight);
+ FillResAndPull(ETrackParType::kTX, trParExtr.GetTx(), trParExtr.GetTxErr(), mcPoint.GetTxOut(), weight);
+ FillResAndPull(ETrackParType::kTY, trParExtr.GetTy(), trParExtr.GetTyErr(), mcPoint.GetTyOut(), weight);
+ FillResAndPull(ETrackParType::kQP, trParExtr.GetQp(), trParExtr.GetQpErr(), mcPoint.GetQpOut(), weight);
- double resX = trParExtr.GetX() - mcPoint.GetXOut(); // residual of x-position [cm]
- double resY = trParExtr.GetY() - mcPoint.GetYOut(); // residual of y-position [cm]
- double resTx = trParExtr.GetTx() - mcPoint.GetTxOut(); // residual of slope along x-axis [1]
- double resTy = trParExtr.GetTy() - mcPoint.GetTyOut(); // residual of slope along y-axis [1]
- double resQp = trParExtr.GetQp() - mcPoint.GetQpOut(); // residual of q/p [ec/GeV]
-
- // TODO: in which point do we calculate MC parameters (center, in, out)??
+ // TODO: in which point do we calculate MC parameters (central, entrance, exit)??
double recoP = std::fabs(mcPoint.GetCharge() / trParExtr.GetQp()); // reco mom. (with MC charge)
double resP = recoP - mcPoint.GetPOut(); // residual of total momentum
double resPhi = trParExtr.GetPhi() - mcPoint.GetPhiOut(); // residual of azimuthal angle
double resTheta = trParExtr.GetTheta() - mcPoint.GetThetaOut(); // residual of polar angle
- double pullX = resX / trParExtr.GetXErr(); // pull of x-position
- double pullY = resY / trParExtr.GetYErr(); // pull of y-position
- double pullTx = resTx / trParExtr.GetTxErr(); // pull of slope along x-axis
- double pullTy = resTy / trParExtr.GetTyErr(); // pull of slope along y-axis
- double pullQp = resQp / trParExtr.GetQpErr(); // pull of q/p
-
- fph_res_x->Fill(resX, weight);
- fph_res_y->Fill(resY, weight);
- fph_res_tx->Fill(resTx, weight);
- fph_res_ty->Fill(resTy, weight);
- fph_res_qp->Fill(resQp, weight);
-
- fph_res_p_vs_pMC->Fill(mcPoint.GetPOut(), resP);
- fph_res_phi_vs_phiMC->Fill(mcPoint.GetPhiOut(), resPhi);
- fph_res_theta_vs_thetaMC->Fill(mcPoint.GetThetaOut(), resTheta);
-
- fph_pull_x->Fill(pullX, weight);
- fph_pull_y->Fill(pullY, weight);
- fph_pull_tx->Fill(pullTx, weight);
- fph_pull_ty->Fill(pullTy, weight);
- fph_pull_qp->Fill(pullQp, weight);
+ resPhi = std::atan2(std::sin(resPhi), std::cos(resPhi)); // overflow over (-pi, pi] protection
+ fph_res_p_pMC->Fill(mcPoint.GetPOut(), resP);
+ fph_res_phi_phiMC->Fill(mcPoint.GetPhiOut(), resPhi);
+ fph_res_theta_thetaMC->Fill(mcPoint.GetThetaOut(), resTheta);
// ** Time-dependent measurements **
if (bTimeMeasured) {
- double resT = trParExtr.GetTime() - mcPoint.GetTime(); // residual of time [ns]
- double resVi = (trParExtr.GetInvSpeed() - mcPoint.GetInvSpeedOut()) * L1Constants::phys::kSpeedOfLightInv;
- double pullT = resT / trParExtr.GetTimeErr(); // pull of time
- double pullVi = resVi / trParExtr.GetInvSpeedErr(); // pull of inverse speed
-
- fph_res_t->Fill(resT, weight);
- fph_res_vi->Fill(resVi, weight);
- fph_pull_t->Fill(pullT, weight);
- fph_pull_vi->Fill(pullVi, weight);
+ FillResAndPull(ETrackParType::kTIME, trParExtr.GetTime(), trParExtr.GetTimeErr(), mcPoint.GetTime(), weight);
+ FillResAndPull(ETrackParType::kVI, trParExtr.GetInvSpeed(), trParExtr.GetInvSpeedErr(), mcPoint.GetInvSpeedOut(),
+ weight);
}
}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFitQa::SetDefaultProperties()
+{
+ // ** Residual distribution parameters **
+ fvRBins[ETrackParType::kX] = 200; ///< Number of bins, residual of x
+ fvRLo[ETrackParType::kX] = -4.e-3; ///< Lower boundary, residual of x [cm]
+ fvRUp[ETrackParType::kX] = +4.e-3; ///< Upper boundary, residual of x [cm]
+ fvRBins[ETrackParType::kY] = 200; ///< Number of bins, residual of y
+ fvRLo[ETrackParType::kY] = -4.e-2; ///< Lower boundary, residual of y [cm]
+ fvRUp[ETrackParType::kY] = +4.e-2; ///< Upper boundary, residual of y [cm]
+ fvRBins[ETrackParType::kTX] = 200; ///< Number of bins, residual of slope along x-axis
+ fvRLo[ETrackParType::kTX] = -4.e-3; ///< Lower boundary, residual of slope along x-axis
+ fvRUp[ETrackParType::kTX] = +4.e-3; ///< Upper boundary, residual of slope along x-axis
+ fvRBins[ETrackParType::kTY] = 200; ///< Number of bins, residual of slope along y-axis
+ fvRLo[ETrackParType::kTY] = -4.e-3; ///< Lower boundary, residual of slope along y-axis
+ fvRUp[ETrackParType::kTY] = +4.e-3; ///< Upper boundary, residual of slope along y-axis
+ fvRBins[ETrackParType::kQP] = 200; ///< Number of bins, residual of q/p
+ fvRLo[ETrackParType::kQP] = -1.; ///< Lower boundary, residual of q/p [ec/GeV]
+ fvRUp[ETrackParType::kQP] = +1.; ///< Upper boundary, residual of q/p [ec/GeV]
+ fvRBins[ETrackParType::kTIME] = 200; ///< Number of bins, residual of time
+ fvRLo[ETrackParType::kTIME] = -1.; ///< Lower boundary, residual of time [ns]
+ fvRUp[ETrackParType::kTIME] = +1.; ///< Upper boundary, residual of time [ns]
+ fvRBins[ETrackParType::kVI] = 200; ///< Number of bins, residual of inverse speed
+ fvRLo[ETrackParType::kVI] = -2.; ///< Lower boundary, residual of inverse speed [1/c]
+ fvRUp[ETrackParType::kVI] = +2.; ///< Upper boundary, residual of inverse speed [1/c]
+
+ // ** Pulls distribution parameters **
+ fvPBins[ETrackParType::kX] = 200; ///< Number of bins, pull of x
+ fvPLo[ETrackParType::kX] = -4.; ///< Lower boundary, pull of x [cm]
+ fvPUp[ETrackParType::kX] = +4.; ///< Upper boundary, pull of x [cm]
+ fvPBins[ETrackParType::kY] = 200; ///< Number of bins, pull of y
+ fvPLo[ETrackParType::kY] = -4.; ///< Lower boundary, pull of y [cm]
+ fvPUp[ETrackParType::kY] = +4.; ///< Upper boundary, pull of y [cm]
+ fvPBins[ETrackParType::kTX] = 200; ///< Number of bins, pull of slope along x-axis
+ fvPLo[ETrackParType::kTX] = -4.; ///< Lower boundary, pull of slope along x-axis
+ fvPUp[ETrackParType::kTX] = +4.; ///< Upper boundary, pull of slope along x-axis
+ fvPBins[ETrackParType::kTY] = 200; ///< Number of bins, pull of slope along y-axis
+ fvPLo[ETrackParType::kTY] = -4.; ///< Lower boundary, pull of slope along y-axis
+ fvPUp[ETrackParType::kTY] = +4.; ///< Upper boundary, pull of slope along y-axis
+ fvPBins[ETrackParType::kQP] = 200; ///< Number of bins, pull of q/p
+ fvPLo[ETrackParType::kQP] = -10.; ///< Lower boundary, pull of q/p [ec/GeV]
+ fvPUp[ETrackParType::kQP] = +10.; ///< Upper boundary, pull of q/p [ec/GeV]
+ fvPBins[ETrackParType::kTIME] = 200; ///< Number of bins, pull of time
+ fvPLo[ETrackParType::kTIME] = -10.; ///< Lower boundary, pull of time [ns]
+ fvPUp[ETrackParType::kTIME] = +10.; ///< Upper boundary, pull of time [ns]
+ fvPBins[ETrackParType::kVI] = 200; ///< Number of bins, pull of inverse speed
+ fvPLo[ETrackParType::kVI] = -2.; ///< Lower boundary, pull of inverse speed [1/c]
+ fvPUp[ETrackParType::kVI] = +2.; ///< Upper boundary, pull of inverse speed [1/c]
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFitQa::SetResidualHistoProperties(ETrackParType type, int nBins, double lo, double up)
+{
+ fvRBins[type] = nBins;
+ fvRLo[type] = lo;
+ fvRUp[type] = up;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFitQa::SetPullHistoProperties(ETrackParType type, int nBins, double lo, double up)
+{
+ fvPBins[type] = nBins;
+ fvPLo[type] = lo;
+ fvPUp[type] = up;
+}
diff --git a/reco/L1/qa/CbmCaTrackFitQa.h b/reco/L1/qa/CbmCaTrackFitQa.h
index 837c7e3eee..91686b14ea 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.h
+++ b/reco/L1/qa/CbmCaTrackFitQa.h
@@ -12,7 +12,10 @@
#include "CbmQaIO.h"
+#include <array>
+
#include "L1Constants.h"
+#include "L1EArray.h"
#include "L1Field.h"
#include "L1Fit.h"
#include "L1Vector.h"
@@ -29,9 +32,36 @@ class CbmL1HitDebugInfo;
class TFolder;
class TH1F;
class TProfile;
+class CbmQaCanvas;
namespace cbm::ca
{
+ /// TODO: SZh 02.05.2023: Move ETrackParType to another class and apply it everywhere to track parameters
+ /// @brief Enumeration for track parameter type
+ enum class ETrackParType
+ {
+ kX, ///< x-position
+ kY, ///< y-position
+ kTX, ///< slope along x-axis
+ kTY, ///< slope along y-axis
+ kQP, ///< charge over total momentum
+ kTIME, ///< time
+ kVI, ///< inverse speed
+ kEND, ///< end of enumeration
+ kBEGIN = 0 ///< begin of enumeration
+ };
+
+ /// @brief Prefix increment operator for ETrackParType
+ ETrackParType operator++(ETrackParType& type);
+
+ /// @brief Postfix increment operator for ETrackParType
+ ETrackParType operator++(ETrackParType& type, int);
+
+
+ /// @brief Alias to array, indexed with ETrackParType enumeration
+ template<typename T>
+ using TrackParArray_t = L1EArray<ETrackParType, T>;
+
/// @brief Set of histograms to monitor track parameters
///
/// Class provides histograms to monitor track fit parameters at a selected z-coordinate.
@@ -79,11 +109,37 @@ namespace cbm::ca
/// @param title Title of fit distributions
void SetTitle(const char* title) { fsTitle = title; }
+ /// @brief Fit histograms
+ void FitHistograms() {}
+
+ /// @brief Creates residuals plot
+ CbmQaCanvas* CreateResidualPlot();
+
+ /// @brief Creates pulls plot
+ CbmQaCanvas* CreatePullPlot();
+
+ /// @brief Creates resolutionis plot
+ CbmQaCanvas* CreateResolutionPlot() { return nullptr; }
+
+ /// @brief Sets properties for a residual histogram
+ /// @param type Type of track parameter
+ /// @param nBins Number of bins
+ /// @param lo Lower boundary
+ /// @param up Upper boundary
+ void SetResidualHistoProperties(ETrackParType type, int nBins, double lo, double up);
+
+ /// @brief Sets properties for a pull histogram
+ /// @param type Type of track parameter
+ /// @param nBins Number of bins
+ /// @param lo Lower boundary
+ /// @param up Upper boundary
+ void SetPullHistoProperties(ETrackParType type, int nBins, double lo, double up);
+
// ************************
// ** List of histograms **
// ************************
- // ** Residuals **
+
TH1F* fph_res_x = nullptr; ///< Residual of x-coordinate [cm]
TH1F* fph_res_y = nullptr; ///< Residual of y-coordinate [cm]
TH1F* fph_res_tx = nullptr; ///< Residual of slope along x-axis
@@ -102,66 +158,64 @@ namespace cbm::ca
TH1F* fph_pull_vi = nullptr; ///< Pull of inverse speed
// ** Resolution profiles **
- TProfile* fph_res_p_vs_pMC = nullptr; ///< Resolution of momentum [GeV/c]
- TProfile* fph_res_theta_vs_thetaMC = nullptr; ///< Resolution of polar angle [rad]
- TProfile* fph_res_phi_vs_phiMC = nullptr; ///< Resolution of azimuthal angle [rad]
+ TProfile* fph_res_p_pMC = nullptr; ///< Resolution of momentum [GeV/c]
+ TProfile* fph_res_phi_phiMC = nullptr; ///< Resolution of azimuthal angle [rad]
+ TProfile* fph_res_theta_thetaMC = nullptr; ///< Resolution of polar angle [rad]
// **************************
// ** Histogram properties **
// **************************
// ** Binning **
- int fBinsRESX = 200; ///< Number of bins, residual of x
- double fLoRESX = -4.e-3; ///< Lower boundary, residual of x [cm]
- double fUpRESX = +4.e-3; ///< Upper boundary, residual of x [cm]
- int fBinsRESY = 200; ///< Number of bins, residual of y
- double fLoRESY = -4.e-2; ///< Lower boundary, residual of y [cm]
- double fUpRESY = +4.e-2; ///< Upper boundary, residual of y [cm]
- int fBinsRESTX = 200; ///< Number of bins, residual of slope along x-axis
- double fLoRESTX = -4.e-3; ///< Lower boundary, residual of slope along x-axis
- double fUpRESTX = +4.e-3; ///< Upper boundary, residual of slope along x-axis
- int fBinsRESTY = 200; ///< Number of bins, residual of slope along y-axis
- double fLoRESTY = -4.e-3; ///< Lower boundary, residual of slope along y-axis
- double fUpRESTY = +4.e-3; ///< Upper boundary, residual of slope along y-axis
- int fBinsRESQP = 200; ///< Number of bins, residual of q/p
- double fLoRESQP = -10.; ///< Lower boundary, residual of q/p [ec/GeV]
- double fUpRESQP = +10.; ///< Upper boundary, residual of q/p [ec/GeV]
- int fBinsREST = 200; ///< Number of bins, residual of time
- double fLoREST = -10.; ///< Lower boundary, residual of time [ns]
- double fUpREST = +10.; ///< Upper boundary, residual of time [ns]
- int fBinsRESVI = 200; ///< Number of bins, residual of inverse speed
- double fLoRESVI = -2.; ///< Lower boundary, residual of inverse speed [1/c]
- double fUpRESVI = +2.; ///< Upper boundary, residual of inverse speed [1/c]
-
- int fBinsPULLX = 200; ///< Number of bins, pull of x
- double fLoPULLX = -4.; ///< Lower boundary, pull of x [cm]
- double fUpPULLX = +4.; ///< Upper boundary, pull of x [cm]
- int fBinsPULLY = 200; ///< Number of bins, pull of y
- double fLoPULLY = -4.; ///< Lower boundary, pull of y [cm]
- double fUpPULLY = +4.; ///< Upper boundary, pull of y [cm]
- int fBinsPULLTX = 200; ///< Number of bins, pull of slope along x-axis
- double fLoPULLTX = -4.; ///< Lower boundary, pull of slope along x-axis
- double fUpPULLTX = +4.; ///< Upper boundary, pull of slope along x-axis
- int fBinsPULLTY = 200; ///< Number of bins, pull of slope along y-axis
- double fLoPULLTY = -4.; ///< Lower boundary, pull of slope along y-axis
- double fUpPULLTY = +4.; ///< Upper boundary, pull of slope along y-axis
- int fBinsPULLQP = 200; ///< Number of bins, pull of q/p
- double fLoPULLQP = -10.; ///< Lower boundary, pull of q/p [ec/GeV]
- double fUpPULLQP = +10.; ///< Upper boundary, pull of q/p [ec/GeV]
- int fBinsPULLT = 200; ///< Number of bins, pull of time
- double fLoPULLT = -10.; ///< Lower boundary, pull of time [ns]
- double fUpPULLT = +10.; ///< Upper boundary, pull of time [ns]
- int fBinsPULLVI = 200; ///< Number of bins, pull of inverse speed
- double fLoPULLVI = -2.; ///< Lower boundary, pull of inverse speed [1/c]
- double fUpPULLVI = +2.; ///< Upper boundary, pull of inverse speed [1/c]
+ static constexpr int kCXSIZEPX = 600; ///< Canvas size along x-axis [px]
+ static constexpr int kCYSIZEPX = 600; ///< Canvas size along y-axis [px]
private:
+ /// @brief Sets default histogram and track parameter properties
+ void SetDefaultProperties();
+
+ using FnVal_t = std::function<double()>;
+ /// @brief Fills residual and pull for a given track parameter
+ /// @param type Type of the track parameter
+ /// @param recoVal Reconstructed error of quantity
+ /// @param recoErr Error of quantity
+ /// @param trueVal True value of quantity
+ /// @param weight Weight
+ void FillResAndPull(ETrackParType type, double recoVal, double recoErr, double trueVal, double w);
+
+ TrackParArray_t<TH1F*> fvphResiduals = {0}; ///< Residuals for different track parameters
+ TrackParArray_t<TH1F*> fvphPulls = {0}; ///< Pulls for different track parameters
+
+ TrackParArray_t<bool> fvbParIgnored = {0}; ///< Flag: true - parameter is ignored
+
+ // ** Distribution properties **
+ TrackParArray_t<int> fvRBins = {0}; ///< Number of bins, residuals
+ TrackParArray_t<double> fvRLo = {0}; ///< Lower boundary, residuals
+ TrackParArray_t<double> fvRUp = {0}; ///< Upper boundary, residuals
+
+ TrackParArray_t<int> fvPBins = {0}; ///< Number of bins, pulls
+ TrackParArray_t<double> fvPLo = {0}; ///< Lower boundary, pulls
+ TrackParArray_t<double> fvPUp = {0}; ///< Upper boundary, pulls
+
TString fsTitle = ""; ///< Title of the point
double fMass = L1Constants::phys::kMuonMass; /// Mass of particle
};
-
+ // *****************************
+ // ** Inline implementation **
+ // *****************************
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ // TODO: Test this function for performance penalties
+ inline void TrackFitQa::FillResAndPull(ETrackParType type, double recoVal, double recoErr, double trueVal, double w)
+ {
+ if (fvbParIgnored[type]) { return; }
+ double res = recoVal - trueVal;
+ double pull = res / recoErr;
+ fvphResiduals[type]->Fill(res, w);
+ fvphPulls[type]->Fill(pull, w);
+ }
} // namespace cbm::ca
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.cxx b/reco/L1/qa/CbmCaTrackTypeQa.cxx
index 8e2514d3dc..232e5a2d6b 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.cxx
+++ b/reco/L1/qa/CbmCaTrackTypeQa.cxx
@@ -45,8 +45,10 @@ void TrackTypeQa::Init()
fph_reco_eta = MakeHisto<TH1F>("reco_eta", "", kBinsETA, kLoETA, kUpETA);
fph_reco_phi = MakeHisto<TH1F>("reco_phi", "", kBinsPHI, kLoPHI, kUpPHI);
fph_reco_theta = MakeHisto<TH1F>("reco_theta", "", kBinsTHETA, kLoTHETA, kUpTHETA);
+ fph_reco_theta_phi = MakeHisto<TH2F>("reco_theta_phi", "", kBinsPHI, kLoPHI, kUpPHI, kBinsTHETA, kLoTHETA, kUpTHETA);
fph_reco_tx = MakeHisto<TH1F>("reco_tx", "", kBinsTX, kLoTX, kUpTX);
fph_reco_ty = MakeHisto<TH1F>("reco_ty", "", kBinsTY, kLoTY, kUpTY);
+ fph_reco_ty_tx = MakeHisto<TH2F>("reco_ty_tx", "", kBinsTX, kLoTX, kUpTX, kBinsTY, kLoTY, kUpTY);
fph_reco_fhitR = MakeHisto<TH1F>("reco_fhitR", "", kBinsFHITR, kLoFHITR, kUpFHITR);
// TODO: ...
@@ -55,8 +57,11 @@ void TrackTypeQa::Init()
fph_reco_phi->SetTitle("Azimuthal angle of reconstructed track;#phi^{reco} [rad];Counts");
fph_reco_eta->SetTitle("Pseudorapidity of reconstructed track;#eta^{reco};Counts");
fph_reco_theta->SetTitle("Polar angle of reconstructed track;#theta^{reco} [rad];Counts");
+ fph_reco_theta_phi->SetTitle(
+ "Polar angle vs. azimuthal angle of reconstructed track;#phi^{reco} [rad];#theta^{reco} [rad];Counts");
fph_reco_tx->SetTitle("Slope along x-axis of reconstructed tracks;t_{x}^{reco};Counts");
fph_reco_ty->SetTitle("Slope along y-axis of reconstructed tracks;t_{y}^{reco};Counts");
+ fph_reco_ty_tx->SetTitle("Slope along y-axis vs. x-axis of reconstructed tracks;t_{x}^{reco};t_{y}^{reco};");
fph_reco_fhitR->SetTitle("Distance of the first hit from z-axis for reconstructed tracks;R^{reco} [cm];Counts");
if (fbUseMC) {
@@ -80,15 +85,24 @@ void TrackTypeQa::Init()
fph_reco_nhits->SetTitle("Hit number of reconstructed tracks;N^{MC}_{hits};Counts");
- fph_mc_pMC = MakeHisto<TH1F>("mc_pMC", "", kBinsP, kLoP, kUpP);
- fph_mc_etaMC = MakeHisto<TH1F>("mc_etaMC", "", kBinsETA, kLoETA, kUpETA);
- fph_mc_yMC = MakeHisto<TH1F>("mc_yMC", "", kBinsY, kLoY, kUpY);
- fph_mc_pMC_yMC = MakeHisto<TH2F>("mc_pMC_yMC", "", kBinsY, kLoY, kUpY, kBinsP, kLoP, kUpP);
+ fph_mc_pMC = MakeHisto<TH1F>("mc_pMC", "", kBinsP, kLoP, kUpP);
+ fph_mc_etaMC = MakeHisto<TH1F>("mc_etaMC", "", kBinsETA, kLoETA, kUpETA);
+ fph_mc_yMC = MakeHisto<TH1F>("mc_yMC", "", kBinsY, kLoY, kUpY);
+ fph_mc_pMC_yMC = MakeHisto<TH2F>("mc_pMC_yMC", "", kBinsY, kLoY, kUpY, kBinsP, kLoP, kUpP);
+ fph_mc_txMC = MakeHisto<TH1F>("mc_txMC", "", kBinsTX, kLoTX, kUpTX);
+ fph_mc_tyMC = MakeHisto<TH1F>("mc_tyMC", "", kBinsTY, kLoTY, kUpTY);
+ fph_mc_tyMC_txMC = MakeHisto<TH2F>("mc_tyMC_txMC", "", kBinsTX, kLoTX, kUpTX, kBinsTY, kLoTY, kUpTY);
+ fph_mc_thetaMC_phiMC =
+ MakeHisto<TH2F>("mc_thetaMC_phiMC", "", kBinsPHI, kLoPHI, kUpPHI, kBinsTHETA, kLoTHETA, kUpTHETA);
fph_mc_pMC->SetTitle("MC total momentum of MC tracks;p^{MC} [GeV/c];Counts");
fph_mc_yMC->SetTitle("MC rapidity of MC tracks;y^{MC};Counts");
fph_mc_etaMC->SetTitle("MC pseudorapidity of MC track;#eta^{MC};Counts");
fph_mc_pMC_yMC->SetTitle("MC total momentum vs. MC rapidity of MC tracks;y^{MC};p^{MC} [GeV/c];Counts");
+ fph_mc_txMC->SetTitle("Slope along x-axis of MC tracks;t_{x}^{MC};Counts");
+ fph_mc_tyMC->SetTitle("Slope along y-axis of MC tracks;t_{y}^{MC};Counts");
+ fph_mc_tyMC_txMC->SetTitle("Slope along y-axis vs. x-axis of MC tracks;t_{x}^{MC};t_{y}^{MC};");
+ fph_mc_thetaMC_phiMC->SetTitle("Polar angle vs. azimuthal angle of MC track;#phi^{MC} [rad];#theta^{MC} [rad]");
//
// ** Efficiencies **
@@ -98,16 +112,22 @@ void TrackTypeQa::Init()
fph_eff_yMC = MakeHisto<TProfile>("eff_yMC", "", kBinsY, kLoY, kUpY, 0., 1.);
fph_eff_ptMC = MakeHisto<TProfile>("eff_ptMC", "", kBinsPT, kLoPT, kUpPT, 0., 1.);
fph_eff_thetaMC = MakeHisto<TProfile>("eff_thetaMC", "", kBinsTHETA, kLoTHETA, kUpTHETA, 0., 1.);
+ fph_eff_etaMC = MakeHisto<TProfile>("eff_etaMC", "", kBinsTHETA, kLoTHETA, kUpTHETA, 0., 1.);
fph_eff_phiMC = MakeHisto<TProfile>("eff_phiMC", "", kBinsPHI, kLoPHI, kUpPHI, 0., 1.);
fph_eff_nhitsMC = MakeHisto<TProfile>("eff_nhitsMC", "", kBinsNHITS, kLoNHITS, kUpNHITS, 0., 1.);
+ fph_eff_txMC = MakeHisto<TProfile>("eff_txMC", "", kBinsTX, kLoTX, kUpTX, 0., 1.);
+ fph_eff_tyMC = MakeHisto<TProfile>("eff_tyMC", "", kBinsTY, kLoTY, kUpTX, 0., 1.);
fph_eff_int->SetTitle("Integrated efficiency;;#epsilon_{CA}");
fph_eff_pMC->SetTitle("Efficiency vs. MC total momentum;p_{MC} [GeV/c];#epsilon_{CA}");
fph_eff_yMC->SetTitle("Efficiency vs. MC rapidity;y_{MC};#epsilon");
fph_eff_ptMC->SetTitle("Efficiency vs. MC total momentum;p_{T}^{MC} [GeV/c];#epsilon_{CA}");
fph_eff_thetaMC->SetTitle("Efficiency vs. MC polar angle;#theta^{MC};#epsilon_{CA}");
+ fph_eff_etaMC->SetTitle("Efficiency vs. MC pseudorapidity;#eta^{MC};#epsilon_{CA}");
fph_eff_phiMC->SetTitle("Efficiency vs. MC azimuthal angle;#phi^{MC};#epsilon_{CA}");
fph_eff_nhitsMC->SetTitle("Efficiency vs. MC number of stations with hits;N_{hit}^{MC};#epsilon_{CA}");
+ fph_eff_txMC->SetTitle("Efficiency vs. MC slope along x-axis;t_{x}^{MC};#epsilon_{CA}");
+ fph_eff_tyMC->SetTitle("Efficiency vs. MC slope along y-axis;t_{y}^{MC};#epsilon_{CA}");
fph_eff_pMC_yMC = MakeHisto<TProfile2D>("eff_pMC_yMC", "", kBinsY, kLoY, kUpY, kBinsP, kLoP, kUpP, 0., 1.);
fph_eff_thetaMC_phiMC =
@@ -125,10 +145,8 @@ void TrackTypeQa::Init()
fpFitQaLastHit = std::make_unique<TrackFitQa>("lst_hit", fsPrefix, fpFolderRoot);
fpFitQaLastHit->SetTitle("Last hit");
- fpFitQaLastHit->fLoRESX = -0.4;
- fpFitQaLastHit->fUpRESX = +0.4;
- fpFitQaLastHit->fLoRESY = -0.8;
- fpFitQaLastHit->fUpRESY = +0.8;
+ fpFitQaLastHit->SetResidualHistoProperties(ETrackParType::kX, 200, -0.4, +0.4);
+ fpFitQaLastHit->SetResidualHistoProperties(ETrackParType::kY, 200, -0.8, +0.8);
fpFitQaLastHit->Init();
fpFitQaVertex = std::make_unique<TrackFitQa>("vertex", fsPrefix, fpFolderRoot);
@@ -149,6 +167,8 @@ void TrackTypeQa::FillRecoTrack(int iTrkReco, double weight)
fph_reco_theta->Fill(recoTrack.GetTheta(), weight);
fph_reco_tx->Fill(recoTrack.GetTx(), weight);
fph_reco_ty->Fill(recoTrack.GetTy(), weight);
+ fph_reco_ty_tx->Fill(recoTrack.GetTx(), recoTrack.GetTy(), weight);
+ fph_reco_theta_phi->Fill(recoTrack.GetPhi(), recoTrack.GetTheta(), weight);
if (fbUseMC) {
int iTrkMC = recoTrack.GetMatchedMCTrackIndex();
if (iTrkMC > -1) {
@@ -234,6 +254,10 @@ void TrackTypeQa::FillMCTrack(int iTrkMC, double weight)
fph_mc_etaMC->Fill(mcTrack.GetEta(), weight);
fph_mc_yMC->Fill(mcTrack.GetRapidity(), weight);
fph_mc_pMC_yMC->Fill(mcTrack.GetRapidity(), mcTrack.GetP(), weight);
+ fph_mc_txMC->Fill(mcTrack.GetTx(), weight);
+ fph_mc_tyMC->Fill(mcTrack.GetTy(), weight);
+ fph_mc_tyMC_txMC->Fill(mcTrack.GetTx(), mcTrack.GetTy(), weight);
+ fph_mc_thetaMC_phiMC->Fill(mcTrack.GetPhi(), mcTrack.GetTheta(), weight);
// ** Efficiencies **
bool bReco = mcTrack.IsReconstructed();
@@ -241,14 +265,17 @@ void TrackTypeQa::FillMCTrack(int iTrkMC, double weight)
// NOTE: Weight is ignored in efficiencies
fph_eff_int->Fill(0., bReco);
fph_eff_pMC->Fill(mcTrack.GetP(), bReco);
+ fph_eff_etaMC->Fill(mcTrack.GetEta(), bReco);
fph_eff_yMC->Fill(mcTrack.GetRapidity(), bReco);
fph_eff_ptMC->Fill(mcTrack.GetPt(), bReco);
fph_eff_thetaMC->Fill(mcTrack.GetTheta(), bReco);
fph_eff_phiMC->Fill(mcTrack.GetPhi(), bReco);
fph_eff_nhitsMC->Fill(mcTrack.GetTotNofStationsWithHit(), bReco);
+ fph_eff_txMC->Fill(mcTrack.GetTx(), bReco);
+ fph_eff_tyMC->Fill(mcTrack.GetTy(), bReco);
fph_eff_pMC_yMC->Fill(mcTrack.GetRapidity(), mcTrack.GetP(), bReco);
- fph_eff_thetaMC_phiMC->Fill(mcTrack.GetPhi(), mcTrack.GetPhi(), bReco);
+ fph_eff_thetaMC_phiMC->Fill(mcTrack.GetPhi(), mcTrack.GetTheta(), bReco);
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.h b/reco/L1/qa/CbmCaTrackTypeQa.h
index a933122d45..c2ddfa617a 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.h
+++ b/reco/L1/qa/CbmCaTrackTypeQa.h
@@ -149,8 +149,10 @@ namespace cbm::ca
TH1F* fph_reco_pt = nullptr; ///< Transverse momentum over charge of reconstructed tracks
TH1F* fph_reco_phi = nullptr; ///< Azimuthal angle of reconstructed tracks
TH1F* fph_reco_theta = nullptr; ///< Polar angle of reconstructed tracks
+ TH2F* fph_reco_theta_phi = nullptr; ///< Polar angle vs. azimuthal angle of reconstructed tracks
TH1F* fph_reco_tx = nullptr; ///< Slope along x-axis of reconstructed tracks
TH1F* fph_reco_ty = nullptr; ///< Slope along y-axis of reconstructed tracks
+ TH2F* fph_reco_ty_tx = nullptr; ///< Slope along x-axis vs y-axis of reconstructed tracks
TH1F* fph_reco_eta = nullptr; ///< Pseudo-rapidity of reconstructed tracks
TH1F* fph_reco_fhitR = nullptr; ///< Distance of the first hit from z-axis for reconstructed tracks
TH1F* fph_reco_nhits = nullptr; ///< Hit number of reconstructed tracks
@@ -170,15 +172,17 @@ namespace cbm::ca
TH1F* fph_reco_tyMC = nullptr; ///< Slope along y-axis of reconstructed tracks
// ** MC track distributions **
- TH1F* fph_mc_pMC = nullptr; ///< MC total momentum over charge of MC tracks
- TH1F* fph_mc_yMC = nullptr; ///< MC rapidity of MC tracks
- TH1F* fph_mc_etaMC = nullptr; ///< MC pseudo-rapidity of MC tracks
- TH2F* fph_mc_pMC_yMC = nullptr; ///< MC total momentum vs. MC rapidity of MC tracks
- TH1F* fph_mc_ptMC = nullptr; ///< Transverse momentum over charge of reconstructed tracks
- TH1F* fph_mc_phiMC = nullptr; ///< Azimuthal angle of reconstructed tracks
- TH1F* fph_mc_thetaMC = nullptr; ///< Polar angle of reconstructed tracks
- TH1F* fph_mc_txMC = nullptr; ///< Slope along x-axis of reconstructed tracks
- TH1F* fph_mc_tyMC = nullptr; ///< Slope along y-axis of reconstructed tracks
+ TH1F* fph_mc_pMC = nullptr; ///< MC total momentum over charge of MC tracks
+ TH1F* fph_mc_yMC = nullptr; ///< MC rapidity of MC tracks
+ TH1F* fph_mc_etaMC = nullptr; ///< MC pseudo-rapidity of MC tracks
+ TH2F* fph_mc_pMC_yMC = nullptr; ///< MC total momentum vs. MC rapidity of MC tracks
+ TH1F* fph_mc_ptMC = nullptr; ///< Transverse momentum over charge of MC tracks
+ TH1F* fph_mc_phiMC = nullptr; ///< Azimuthal angle of MC tracks
+ TH1F* fph_mc_thetaMC = nullptr; ///< Polar angle of MC tracks
+ TH2F* fph_mc_thetaMC_phiMC = nullptr; ///< Polar angle vs. azimuthal angle of MC tracks
+ TH1F* fph_mc_txMC = nullptr; ///< Slope along x-axis of MC tracks
+ TH1F* fph_mc_tyMC = nullptr; ///< Slope along y-axis of MC tracks
+ TH2F* fph_mc_tyMC_txMC = nullptr; ///< Slope along x-axis vs y-axis of MC tracks
// ** Efficiencies **
TProfile* fph_eff_int = nullptr; ///< Integrated efficiency
@@ -186,8 +190,12 @@ namespace cbm::ca
TProfile* fph_eff_yMC = nullptr; ///< Efficiency vs. MC rapidity
TProfile* fph_eff_ptMC = nullptr; ///< Efficiency vs. MC transverse momentum
TProfile* fph_eff_thetaMC = nullptr; ///< Efficiency vs. MC polar angle
+ TProfile* fph_eff_etaMC = nullptr; ///< Efficiency vs. MC pseudorapidity
TProfile* fph_eff_phiMC = nullptr; ///< Efficiency vs. MC azimuthal angle
TProfile* fph_eff_nhitsMC = nullptr; ///< Efficiency vs. MC number of hits
+ TProfile* fph_eff_txMC = nullptr; ///< Efficiency vs. MC slope along x-axis
+ TProfile* fph_eff_tyMC = nullptr; ///< Efficiency vs. MC slope along y-axis
+
TProfile2D* fph_eff_thetaMC_phiMC = nullptr; ///< Efficiency vs. MC theta and MC phi
TProfile2D* fph_eff_pMC_yMC = nullptr; ///< Efficiency vs. MC momentum and MC rapidity
--
GitLab