diff --git a/reco/qa/CbmRecoQaTask.cxx b/reco/qa/CbmRecoQaTask.cxx
index 57ddfe9a6b4afccf381f9a32d012186c5e1adfb7..bbc881af9db2537f22baca29e8003a6390e5b33c 100644
--- a/reco/qa/CbmRecoQaTask.cxx
+++ b/reco/qa/CbmRecoQaTask.cxx
@@ -1,8 +1,8 @@
/* Copyright (C) 2024 Hulubei National Institute of Physics and Nuclear Engineering - Horia Hulubei, Bucharest
SPDX-License-Identifier: GPL-3.0-only
- Authors: Alexandru Bercuci [committer] */
+ Authors: Alexandru Bercuci [committer], Omveer Singh */
-// CBM headers
+//CBM headers
#include "CbmRecoQaTask.h"
#include "../../core/detectors/rich/utils/CbmRichUtil.h"
@@ -10,18 +10,23 @@
#include "CbmEvent.h"
#include "CbmFsdHit.h"
#include "CbmGlobalTrack.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCDataManager.h"
#include "CbmMvdHit.h"
#include "CbmPsdHit.h"
#include "CbmRichHit.h"
#include "CbmSetup.h"
#include "CbmStsAddress.h"
#include "CbmStsHit.h"
+#include "CbmStsSetup.h"
#include "CbmTimeSlice.h"
#include "CbmTofAddress.h"
#include "CbmTofHit.h"
#include "CbmTrdAddress.h"
#include "CbmTrdHit.h"
// FAIR headers
+#include "FairMCPoint.h"
+
#include <FairRootManager.h>
#include <FairSink.h>
// ROOT headers
@@ -30,6 +35,9 @@
#include <TGeoNode.h>
#include <TH2D.h>
+#include <iterator>
+#include <regex>
+
using namespace std;
std::bitset<kRecoQaNConfigs> CbmRecoQaTask::fuRecoConfig = {};
@@ -39,8 +47,9 @@ CbmRecoQaTask::CbmRecoQaTask() : FairTask("RecoQA", 0) {}
//_____________________________________________________________________
CbmRecoQaTask::Detector* CbmRecoQaTask::AddDetector(ECbmModuleId id)
{
- /* Interface to the user. Check that the det defined in the outside world fulfills the quality criteria.
- */
+ /* Interface to the user. Check that the det defined in the outside world
+ * fulfills the quality criteria.
+ */
switch (id) {
case ECbmModuleId::kMvd:
case ECbmModuleId::kSts:
@@ -53,7 +62,8 @@ CbmRecoQaTask::Detector* CbmRecoQaTask::AddDetector(ECbmModuleId id)
case ECbmModuleId::kPsd:
if (GetDetector(id)) {
LOG(warn) << GetName() << "::AddDetector : det " << ToString(id)
- << " already registered. Using \"CbmRecoQaTask::GetDetector(ECbmModuleId).\"";
+ << " already registered. Using "
+ "\"CbmRecoQaTask::GetDetector(ECbmModuleId).\"";
return GetDetector(id);
}
LOG(debug) << GetName() << "::AddDetector(" << ToString(id) << ").";
@@ -97,7 +107,7 @@ InitStatus CbmRecoQaTask::Init()
fFitter.Init();
fFitter.SetSkipUnmeasuredCoordinates(true);
- //Get ROOT Manager
+ // Get ROOT Manager
FairRootManager* ioman = FairRootManager::Instance();
if (!ioman) {
@@ -133,19 +143,38 @@ InitStatus CbmRecoQaTask::Init()
LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2024\". ";
InitMcbm24();
break;
- case kDefault: LOG(info) << GetName() << "::Init: Setup config for \"DEFAULT\". Not implemented."; break;
+ case kDefault:
+ LOG(info) << GetName() << "::Init: Setup config for \"Default\". ";
+ InitDefault();
+ break;
}
fOutFolder.Clear();
for (auto& detp : fDetQa) { // SYSTEM LOOP
auto& det = detp.second;
+ cbm_mc_manager = dynamic_cast<CbmMCDataManager*>(ioman->GetObject("MCDataManager"));
+
+ if (cbm_mc_manager) {
+ fPoints[det.id] = cbm_mc_manager->InitBranch(det.point.name.data());
+ if (!fPoints[det.id])
+ LOG(warn) << GetName() << "::Init: MC Point array for " << ToString(det.id) << " not found!";
+ }
+ else {
+ LOG(warn) << GetName() << "::Init: MCDataManager not found in Tree.";
+ }
+
fHits[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.hit.name.data()));
if (!fHits[det.id]) {
LOG(warn) << GetName() << "::Init: Hit array for " << ToString(det.id) << " not found!";
continue;
}
+ fHitMatch[det.id] = static_cast<TClonesArray*>(ioman->GetObject((det.hit.name + "Match").c_str()));
+
+ if (!fHitMatch[det.id])
+ LOG(warn) << GetName() << "::Init: Hit Match array for " << ToString(det.id) << " not found!";
+
if (fGTracks && det.trk.id != ECbmDataType::kUnknown) {
fTracks[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.trk.name.data()));
if (!fTracks[det.id]) LOG(warn) << GetName() << "::Init: Track array for " << ToString(det.id) << " not found!";
@@ -153,8 +182,8 @@ InitStatus CbmRecoQaTask::Init()
if (!det.Init(&fOutFolder)) continue;
}
- // TDirectoryFile* trkDir = (TDirectoryFile*)fOutFolder.mkdir("track", "CA QA");
- // for (auto& trkly : fTrkView) {
+ // TDirectoryFile* trkDir = (TDirectoryFile*)fOutFolder.mkdir("track", "CA
+ // QA"); for (auto& trkly : fTrkView) {
// for (auto& view : trkly.second) view.Register(trkDir);
// }
@@ -167,7 +196,7 @@ InitStatus CbmRecoQaTask::Init()
/// Hit classification on system and view
template<class Hit>
-bool CbmRecoQaTask::Detector::View::HasAddress(const CbmHit* h, double&, double&) const
+bool CbmRecoQaTask::Detector::View::HasAddress(const CbmHit* h, double&, double&, double&, double&) const
{
LOG(error) << "Unprocessed hit in view " << name;
cout << h->ToString();
@@ -176,7 +205,8 @@ bool CbmRecoQaTask::Detector::View::HasAddress(const CbmHit* h, double&, double&
}
// STS specialization
template<>
-bool CbmRecoQaTask::Detector::View::HasAddress<CbmStsHit>(const CbmHit* h, double& x, double& y) const
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmStsHit>(const CbmHit* h, double& x, double& y, double& dx,
+ double& dy) const
{
int32_t a = h->GetAddress();
uint32_t sel[] = {CbmStsAddress::GetElementId(a, EStsElementLevel::kStsUnit),
@@ -187,8 +217,12 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmStsHit>(const CbmHit* h, doubl
if (uint(ii) != sel[idx]) break;
idx++;
}
- if (idx != fSelector.size())
+ if (fSetup == eSetup::kDefault && sel[0] != uint(fSelector[0])) {
return false;
+ }
+ else if (fSetup != eSetup::kDefault && idx != fSelector.size()) {
+ return false;
+ }
else
LOG(debug4) << "Accept Sts hit for " << sel[0] << " " << sel[1] << " " << sel[2];
@@ -197,51 +231,78 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmStsHit>(const CbmHit* h, doubl
LOG(error) << "Failed loading STS hit in view " << name;
return false;
}
- x = h0->GetX();
- y = h0->GetY();
+ x = h0->GetX();
+ dx = h0->GetDx();
+ y = h0->GetY();
+ dy = h0->GetDy();
return true;
}
// Rich specialization
template<>
-bool CbmRecoQaTask::Detector::View::HasAddress<CbmRichHit>(const CbmHit* h, double& x, double& y) const
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmRichHit>(const CbmHit* h, double& x, double& y, double& dx,
+ double& dy) const
{
int16_t uId = CbmRichUtil::GetDirichId(h->GetAddress());
// TODO implement at RichUtil level
int8_t modId = (uId >> 8) & 0xF;
- if (modId >= fSelector[0]) return false;
+
+ if (fSetup != eSetup::kDefault) {
+ if (modId >= fSelector[0]) return false;
+ }
const CbmRichHit* h0 = dynamic_cast<const CbmRichHit*>(h);
if (!h0) {
LOG(error) << "Failed loading RICH hit in view " << name;
return false;
}
- x = h0->GetX();
- y = h0->GetY();
+ x = h0->GetX();
+ dx = h0->GetDx();
+ y = h0->GetY();
+ dy = h0->GetDy();
return true;
}
// TRD specialization
template<>
-bool CbmRecoQaTask::Detector::View::HasAddress<CbmTrdHit>(const CbmHit* h, double& x, double& y) const
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmTrdHit>(const CbmHit* h, double& x, double& y, double& dx,
+ double& dy) const
{
- int16_t uId = CbmTrdAddress::GetModuleAddress(h->GetAddress());
- if (uId != fSelector[0])
+ std::vector<int> sel;
+
+ if (fSetup == eSetup::kDefault) {
+ sel = {static_cast<int>(CbmTrdAddress::GetLayerId(CbmTrdAddress::GetModuleAddress(h->GetAddress()))),
+ static_cast<int>(CbmTrdAddress::GetModuleId(h->GetAddress())), -1};
+ }
+ else {
+ sel = {static_cast<int>(CbmTrdAddress::GetModuleAddress(h->GetAddress())), -1, -1};
+ }
+
+ uint idx = 0;
+ for (auto ii : fSelector) {
+ if (ii != sel[idx]) break;
+ idx++;
+ }
+
+ if (idx != fSelector.size()) {
return false;
- else
- LOG(debug4) << "Accept Trd hit for " << uId;
+ }
+ LOG(debug4) << "Accept Sts hit for " << sel[0] << " " << sel[1] << " " << sel[2];
const CbmTrdHit* h0 = dynamic_cast<const CbmTrdHit*>(h);
if (!h0) {
LOG(error) << "Failed loading TRD hit in view " << name;
return false;
}
- x = h0->GetX();
- y = h0->GetY();
+ x = h0->GetX();
+ dx = h0->GetDx();
+ y = h0->GetY();
+ dy = h0->GetDy();
return true;
}
// ToF specialization
template<>
-bool CbmRecoQaTask::Detector::View::HasAddress<CbmTofHit>(const CbmHit* h, double& x, double& y) const
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmTofHit>(const CbmHit* h, double& x, double& y, double& dx,
+ double& dy) const
{
int32_t a = h->GetAddress();
int32_t sel[] = {CbmTofAddress::GetSmId(a), CbmTofAddress::GetSmType(a), CbmTofAddress::GetRpcId(a)};
@@ -250,8 +311,12 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmTofHit>(const CbmHit* h, doubl
if (ii != sel[idx]) break;
idx++;
}
- if (idx != fSelector.size())
+ if (fSetup == eSetup::kDefault && idx != 2) {
return false;
+ }
+ else if (fSetup != eSetup::kDefault && idx != fSelector.size()) {
+ return false;
+ }
else
LOG(debug4) << "Accept Tof hit for " << sel[0] << " " << sel[1] << " " << sel[2];
@@ -260,38 +325,55 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmTofHit>(const CbmHit* h, doubl
LOG(error) << "Failed loading ToF hit in view " << name;
return false;
}
- x = h0->GetX();
- y = h0->GetY();
+ x = h0->GetX();
+ dx = h0->GetDx();
+ y = h0->GetY();
+ dy = h0->GetDy();
return true;
}
template<class Hit>
-bool CbmRecoQaTask::Detector::View::Load(const CbmHit* h, const CbmEvent* ev)
+bool CbmRecoQaTask::Detector::View::Load(const CbmHit* h, const FairMCPoint* point, const CbmEvent* ev)
{
- double x(0), y(0);
- if (!HasAddress<Hit>(h, x, y)) return false;
+ double x(0), y(0), dx(0), dy(0);
+ if (!HasAddress<Hit>(h, x, y, dx, dy)) {
+ LOG(debug1) << "view " << name << " does not own hit " << h->ToString();
+ return false;
+ }
- //printf("View[%s] z(%d)[%f]\n", name.data(), (int)h->GetType(), h->GetZ());
+ // printf("View[%s] z(%d)[%f]\n", name.data(), (int)h->GetType(), h->GetZ());
int scale(0);
TH2* hh(nullptr);
- if (ev) {
- if (fProjection.find(eViewProjection::kChdT) != fProjection.end()) {
- scale = get<0>(fProjection[eViewProjection::kChdT]);
- hh = get<2>(fProjection[eViewProjection::kChdT]);
- hh->Fill(x, scale * (h->GetTime() - ev->GetStartTime()));
- }
- if (fProjection.find(eViewProjection::kXYh) != fProjection.end()) {
- hh = get<2>(fProjection[eViewProjection::kXYh]);
- hh->Fill(x, y);
+ auto fillView = [&](eViewProjection proj, double xx, double yy, bool scaleY = true) {
+ auto it = fProjection.find(proj);
+ if (it != fProjection.end()) {
+ scale = std::get<0>(it->second);
+ hh = std::get<2>(it->second);
+ hh->Fill(xx, scaleY ? yy * scale : yy);
}
+ };
+
+ if (ev) {
+ fillView(eViewProjection::kChdT, x, h->GetTime() - ev->GetStartTime());
+ fillView(eViewProjection::kXYh, x, y, false);
}
- if (fProjection.find(eViewProjection::kXYs) != fProjection.end()) {
- hh = get<2>(fProjection[eViewProjection::kXYs]);
- hh->Fill(x, y);
+
+ fillView(eViewProjection::kXYs, x, y, false);
+
+ double event_time = ev ? ev->GetStartTime() : 1000;
+
+ if (point) {
+ fillView(eViewProjection::kXYhMC, point->GetX(), point->GetY(), false);
+ fillView(eViewProjection::kResidualX, point->GetX(), x - point->GetX());
+ fillView(eViewProjection::kResidualY, point->GetY(), y - point->GetY());
+ fillView(eViewProjection::kResidualTX, point->GetX(), h->GetTime() - event_time - point->GetTime());
+ fillView(eViewProjection::kResidualTY, point->GetY(), h->GetTime() - event_time - point->GetTime());
+ fillView(eViewProjection::kPullX, point->GetX(), (x - point->GetX()) / dx);
+ fillView(eViewProjection::kPullY, point->GetY(), (y - point->GetY()) / dy);
}
return true;
}
-bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::FitNode* n)
+bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::FitNode* n, const FairMCPoint* point)
{
double dx = n->fMxy.X() - n->fTrack.X(), dy = n->fMxy.Y() - n->fTrack.Y(), dt = n->fMt.T() - n->fTrack.Time(),
pullx = dx / sqrt(n->fMxy.Dx2() + n->fTrack.GetCovariance(0, 0)),
@@ -313,11 +395,22 @@ bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::FitNode* n)
case eViewProjection::kYpY: hh->Fill(n->fMxy.Y(), pully); break;
default: break;
}
+
+ if (point) {
+ double dxMC = point->GetX() - n->fTrack.X(), dyMC = point->GetY() - n->fTrack.Y();
+
+ switch (projection.first) {
+ case eViewProjection::kXdXMC: hh->Fill(point->GetX(), scale * dxMC); break;
+ case eViewProjection::kYdYMC: hh->Fill(point->GetY(), scale * dyMC); break;
+ default: break;
+ }
+ }
}
return true;
}
//_____________________________________________________________________
+
void CbmRecoQaTask::Exec(Option_t*)
{
LOG(info) << GetName() << "::Exec : Evs[" << (fEvents ? fEvents->GetEntriesFast() : 0) << "] Trks["
@@ -326,228 +419,190 @@ void CbmRecoQaTask::Exec(Option_t*)
fFitter.SetDefaultMomentumForMs(0.5);
fFitter.FixMomentumForMs(true); // consider only mcbm case for a moment
- int iev(0), itrack(0), nnodes(0);
- if (fEvents) {
- for (auto evObj : *fEvents) { // EVENT LOOP
- auto ev = dynamic_cast<CbmEvent*>(evObj);
- // track QA filtering
- if (!FilterEvent(ev)) continue;
-
- for (auto detp : fDetQa) { // SYSTEM LOOP
- auto det = detp.second;
- // check SYSTEM @ hit level
- int nh = max(int(0), int(ev->GetNofData(det.hit.id)));
-
- if (!fHits[det.id]) {
- LOG(error) << GetName() << "::Exec() : Hits for " << ToString(det.id) << " not available. Skip.";
- continue;
- }
- for (int ih(0); ih < nh; ih++) { // HIT LOOP
- const int jh = ev->GetIndex(det.hit.id, ih);
- const CbmHit* hit = dynamic_cast<CbmHit*>(fHits[det.id]->At(jh));
- // printf("%s hit_%2d dt=%f time_event=%f time_hit=%f\n", ToString(det.id).data(), ih, ev->GetStartTime() - hit->GetTime(), ev->GetStartTime(), hit->GetTime());
-
- bool ret(false);
- for (auto view : det.fViews) { // VIEW LOOP
- switch (det.id) {
- case ECbmModuleId::kMvd: ret = view.Load<CbmMvdHit>(hit, ev); break;
- case ECbmModuleId::kSts: ret = view.Load<CbmStsHit>(hit, ev); break;
- case ECbmModuleId::kMuch: ret = view.Load<CbmMuchPixelHit>(hit, ev); break;
- case ECbmModuleId::kRich: ret = view.Load<CbmRichHit>(hit, ev); break;
- case ECbmModuleId::kTrd: ret = view.Load<CbmTrdHit>(hit, ev); break;
- case ECbmModuleId::kTrd2d: ret = view.Load<CbmTrdHit>(hit, ev); break;
- case ECbmModuleId::kTof: ret = view.Load<CbmTofHit>(hit, ev); break;
- case ECbmModuleId::kFsd: ret = view.Load<CbmFsdHit>(hit, ev); break;
- case ECbmModuleId::kPsd: ret = view.Load<CbmPsdHit>(hit, ev); break;
- default: LOG(fatal) << GetName() << "::Exec : unsupported det " << ToString(det.id); break;
- }
- if (ret) break;
- } // END VIEW LOOP
- } // END HIT LOOP
- } // END SYSTEM LOOP
- if (!fGTracks) continue;
+ int iev = 0, itrack = 0, nnodes = 0;
- int ntrk = max(int(0), int(ev->GetNofData(ECbmDataType::kGlobalTrack)));
- for (int itrk(0); itrk < ntrk; itrk++) { // TRACK LOOP
- auto trackObj = (CbmGlobalTrack*) (*fGTracks)[ev->GetIndex(ECbmDataType::kGlobalTrack, itrk)];
+ auto processHits = [&](CbmEvent* ev) {
+ for (auto& detp : fDetQa) {
+ auto& det = detp.second;
+ if (!fHits[det.id]) {
+ LOG(error) << GetName() << "::Exec() : Hits for " << ToString(det.id) << " not available. Skip.";
+ continue;
+ }
- // track QA filtering
- if (!FilterTrack(trackObj)) continue;
+ const int nh = (ev) ? max(int(0), int(ev->GetNofData(det.hit.id))) : fHits[det.id]->GetEntriesFast();
+ for (int ih = 0; ih < nh; ++ih) {
+ const int jh = (ev) ? ev->GetIndex(det.hit.id, ih) : ih;
- auto track = dynamic_cast<const CbmGlobalTrack*>(trackObj);
- CbmKfTrackFitter::Track trkKf;
- if (!fFitter.CreateGlobalTrack(trkKf, *track)) {
- LOG(fatal) << GetName() << "::Exec: can not create the track for the fit! ";
- break;
- }
+ const FairMCPoint* mcpoint = nullptr;
+ const CbmHit* hit = dynamic_cast<CbmHit*>(fHits[det.id]->At(jh));
+ if (!hit) continue;
- // minimalistic QA tool for tracks used for target projection
- int nhit[(size_t) ECbmModuleId::kLastModule] = {0};
- for (auto& n : trkKf.fNodes) {
- // check if hit is well defined
- if (n.fHitSystemId == ECbmModuleId::kNotExist) continue;
- // check if detector is registered
- if (fDetQa.find(n.fHitSystemId) == fDetQa.end()) continue;
-
- auto det = fDetQa[n.fHitSystemId];
- // det.Print();
- auto view = det.FindView(n.fMxy.X(), n.fMxy.Y(), n.fZ);
- if (!view) {
- LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
- continue;
+ if (fHitMatch[det.id] && fPoints[det.id]) {
+ if (det.id == ECbmModuleId::kRich) {
+ Int_t pointID = hit->GetRefId();
+ if (pointID >= 0) {
+ // mcpoint = dynamic_cast<FairMCPoint*>(fPoints[det.id]->At(pointID));
+ if (!mcpoint) continue;
+ }
+ }
+ else {
+ const CbmMatch* match = dynamic_cast<CbmMatch*>(fHitMatch[det.id]->At(jh));
+
+ if (match && match->GetNofLinks() > 0) {
+ const auto& link = match->GetMatchedLink();
+
+ if (ev) {
+ int file_id{0}, event_id{0};
+
+ if (ev && ev->GetMatch() && ev->GetMatch()->GetNofLinks() > 0) {
+ file_id = ev->GetMatch()->GetMatchedLink().GetFile();
+ event_id = ev->GetMatch()->GetMatchedLink().GetEntry();
+ }
+ else {
+ event_id = FairRootManager::Instance()->GetEntryNr();
+ }
+ if (link.GetFile() != file_id || link.GetEntry() != event_id) {
+ LOG(warn) << "match from different event";
+ }
+ }
+ mcpoint = dynamic_cast<FairMCPoint*>(fPoints[det.id]->Get(link));
+ }
}
-
- n.fIsTimeSet = false;
- n.fIsXySet = false;
- fFitter.FitTrack(trkKf);
- view->Load(&n);
- nhit[(int) n.fHitSystemId]++;
- n.fIsTimeSet = true;
- n.fIsXySet = true;
}
- // Fit track for intermediate projections
- int iprj(0);
- if (!nnodes) nnodes = (int) trkKf.fNodes.size();
- for (auto zprj : fProjs) {
- CbmKfTrackFitter::FitNode n = CbmKfTrackFitter::FitNode();
- n.fZ = zprj.first;
- n.fIsXySet = true;
- n.fReference1 = iprj++;
- trkKf.fNodes.push_back(n);
- }
- trkKf.OrderNodesInZ();
- fFitter.FitTrack(trkKf);
- for (auto& n : trkKf.fNodes) {
- if (n.fReference1 < 0) continue;
- if (fProjs.find(n.fReference1 + nnodes) == fProjs.end()) {
- LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fReference1 + nnodes << " not defined.";
- continue;
+ bool ret(false);
+ for (auto& view : det.fViews) {
+ switch (det.id) {
+ case ECbmModuleId::kMvd: ret = view.Load<CbmMvdHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kSts: ret = view.Load<CbmStsHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kMuch: ret = view.Load<CbmMuchPixelHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kRich: ret = view.Load<CbmRichHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kTrd: ret = view.Load<CbmTrdHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kTrd2d: ret = view.Load<CbmTrdHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kTof: ret = view.Load<CbmTofHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kFsd: ret = view.Load<CbmFsdHit>(hit, mcpoint, ev); break;
+ case ECbmModuleId::kPsd: ret = view.Load<CbmPsdHit>(hit, mcpoint, ev); break;
+ default: LOG(fatal) << GetName() << "::Exec : unsupported det " << ToString(det.id); break;
}
-
- // View* mod = &fProjs[n.fReference1 + nnodes][0];
- // mod->Load(&n);
- //
- // if (mod->hdXx.second)
- // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X(), mod->hdXx.first * n.fTrack.Y());
- // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X(), n.fTrack.Time());
+ if (ret) break;
}
- itrack++;
- } // END TRACK LOOP
- iev++;
- } // END EVENT LOOP
-
- LOG(info) << GetName() << "::Exec : Evs(%)[" << 1.e2 * iev / fEvents->GetEntriesFast() << "] Trks(%)["
- << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
- return; // END EbyE QA
- }
-
- // FULL TS LOOP
- for (auto detp : fDetQa) { // SYSTEM LOOP
- auto det = detp.second;
- // check SYSTEM @ hit level
- if (!fHits[det.id]) {
- LOG(error) << GetName() << "::Exec() : Hits for " << ToString(det.id) << " not available. Skip.";
- continue;
+ }
}
- int nh = fHits[det.id]->GetEntriesFast();
-
- for (int ih(0); ih < nh; ih++) { // HIT LOOP
- const CbmHit* hit = dynamic_cast<CbmHit*>(fHits[det.id]->At(ih));
- // printf("%s hit_%2d dt=%f time_event=%f time_hit=%f\n", ToString(det.id).data(), ih, ev->GetStartTime() - hit->GetTime(), ev->GetStartTime(), hit->GetTime());
-
- bool ret(false);
- for (auto view : det.fViews) { // VIEW LOOP
- switch (det.id) {
- case ECbmModuleId::kMvd: ret = view.Load<CbmMvdHit>(hit); break;
- case ECbmModuleId::kSts: ret = view.Load<CbmStsHit>(hit); break;
- case ECbmModuleId::kMuch: ret = view.Load<CbmMuchPixelHit>(hit); break;
- case ECbmModuleId::kRich: ret = view.Load<CbmRichHit>(hit); break;
- case ECbmModuleId::kTrd: ret = view.Load<CbmTrdHit>(hit); break;
- case ECbmModuleId::kTrd2d: ret = view.Load<CbmTrdHit>(hit); break;
- case ECbmModuleId::kTof: ret = view.Load<CbmTofHit>(hit); break;
- case ECbmModuleId::kFsd: ret = view.Load<CbmFsdHit>(hit); break;
- case ECbmModuleId::kPsd: ret = view.Load<CbmPsdHit>(hit); break;
- default: LOG(fatal) << GetName() << "::Exec : unsupported det " << ToString(det.id); break;
+ };
+
+ auto processTracks = [&](CbmEvent* ev) {
+ const int ntrk = (ev) ? ev->GetNofData(ECbmDataType::kGlobalTrack) : fGTracks->GetEntriesFast();
+ for (int itrk = 0; itrk < ntrk; ++itrk) {
+ auto trackObj = (ev) ? (CbmGlobalTrack*) (*fGTracks)[ev->GetIndex(ECbmDataType::kGlobalTrack, itrk)]
+ : static_cast<CbmGlobalTrack*>((*fGTracks)[itrk]);
+ if (!FilterTrack(trackObj)) continue;
+
+ auto track = dynamic_cast<const CbmGlobalTrack*>(trackObj);
+ CbmKfTrackFitter::Track trkKf;
+ if (!fFitter.CreateGlobalTrack(trkKf, *track)) {
+ LOG(fatal) << GetName() << "::Exec: can not create the track for the fit!";
+ break;
+ }
+
+ if (!nnodes) nnodes = trkKf.fNodes.size();
+
+ // minimalistic QA tool for tracks used for target projection
+ int nhit[(size_t) ECbmModuleId::kLastModule] = {0};
+ for (auto& n : trkKf.fNodes) {
+ // check if hit is well defined and detector is registered
+ if (n.fHitSystemId == ECbmModuleId::kNotExist || fDetQa.find(n.fHitSystemId) == fDetQa.end()) continue;
+ auto& det = fDetQa[n.fHitSystemId];
+ // det.Print();
+ auto view = det.FindView(n.fMxy.X(), n.fMxy.Y(), n.fZ);
+ if (!view) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
+ continue;
}
- if (ret) break;
- } // END VIEW LOOP
- } // END HIT LOOP
- } // END SYSTEM LOOP
- if (!fGTracks) {
- LOG(info) << GetName() << "::Exec : TS local reco only.";
- return; // END TS QA (no tracking)
- }
- for (auto trackObj : *fGTracks) { // TRACK LOOP
- auto track = dynamic_cast<const CbmGlobalTrack*>(trackObj);
+ n.fIsTimeSet = n.fIsXySet = false;
+ fFitter.FitTrack(trkKf);
- // track QA filtering
- if (!FilterTrack(track)) continue;
+ // match to MC
+ const CbmHit* hit = dynamic_cast<CbmHit*>(fHits[n.fHitSystemId]->At(n.fHitIndex));
+ if (!hit) continue;
+ const FairMCPoint* mcpoint = nullptr;
- CbmKfTrackFitter::Track trkKf;
- if (!fFitter.CreateGlobalTrack(trkKf, *track)) {
- LOG(fatal) << GetName() << "::Exec: can not create the track for the fit! ";
- break;
- }
+ if (fHitMatch[n.fHitSystemId] && fPoints[n.fHitSystemId]) {
+ if (n.fHitSystemId == ECbmModuleId::kRich) {
+ Int_t pointID = hit->GetRefId();
+ if (pointID >= 0) {
+ // mcpoint = dynamic_cast<FairMCPoint*>(fPoints[n.fHitSystemId]->At(pointID));
+ }
+ }
+ else {
+ const CbmMatch* match = dynamic_cast<CbmMatch*>(fHitMatch[n.fHitSystemId]->At(n.fHitIndex));
+ if (match && match->GetNofLinks() > 0) {
+ const auto& link = match->GetMatchedLink();
+ mcpoint = dynamic_cast<FairMCPoint*>(fPoints[n.fHitSystemId]->Get(link));
+ }
+ }
+ }
- if (!nnodes) nnodes = (int) trkKf.fNodes.size();
-
- // minimalistic QA tool for tracks used for target projection
- int nhit[(size_t) ECbmModuleId::kLastModule] = {0};
- for (auto& n : trkKf.fNodes) {
- // check if hit is well defined
- if (n.fHitSystemId == ECbmModuleId::kNotExist) continue;
- // check if detector is registered
- if (fDetQa.find(n.fHitSystemId) == fDetQa.end()) continue;
-
- auto det = fDetQa[n.fHitSystemId];
- // det.Print();
- auto view = det.FindView(n.fMxy.X(), n.fMxy.Y(), n.fZ);
- if (!view) {
- LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
- continue;
+ view->Load(&n, mcpoint);
+ nhit[(int) n.fHitSystemId]++;
+ n.fIsTimeSet = n.fIsXySet = true;
}
- n.fIsTimeSet = false;
- n.fIsXySet = false;
+ // Fit track with all hits ON for vx projection
+ int iprj(0);
+ for (auto& zprj : fProjs) {
+ CbmKfTrackFitter::FitNode n = CbmKfTrackFitter::FitNode();
+ n.fZ = zprj.first;
+ n.fIsXySet = true;
+ n.fReference1 = iprj++;
+ trkKf.fNodes.push_back(n);
+ }
+ trkKf.OrderNodesInZ();
fFitter.FitTrack(trkKf);
- view->Load(&n);
- nhit[(int) n.fHitSystemId]++;
- n.fIsTimeSet = true;
- n.fIsXySet = true;
+ for (auto& n : trkKf.fNodes) {
+ if (n.fReference1 < 0) continue;
+ if (fProjs.find(n.fReference1 + nnodes) == fProjs.end()) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fReference1 + nnodes << " not defined.";
+ continue;
+ }
+ // View* mod = &fProjs[n.fReference1 + nnodes][0];
+ // mod->Load(&n);
+ // if (mod->hdXx.second)
+ // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X()[0],
+ // mod->hdXx.first * n.fTrack.Y()[0]);
+ // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first *
+ // n.fTrack.X()[0], n.fTrack.Time()[0]);
+ }
+ itrack++;
}
+ };
- // Fit track with all hits ON for vx projection
- int iprj(0);
- for (auto zprj : fProjs) {
- CbmKfTrackFitter::FitNode n = CbmKfTrackFitter::FitNode();
- n.fZ = zprj.first;
- n.fIsXySet = true;
- n.fReference1 = iprj++;
- trkKf.fNodes.push_back(n);
+ if (fEvents) {
+ for (auto evObj : *fEvents) {
+ auto ev = dynamic_cast<CbmEvent*>(evObj);
+ if (!FilterEvent(ev)) continue;
+ processHits(ev);
+ if (fGTracks) processTracks(ev);
+ iev++;
}
- trkKf.OrderNodesInZ();
- fFitter.FitTrack(trkKf);
- for (auto& n : trkKf.fNodes) {
- if (n.fReference1 < 0) continue;
- if (fProjs.find(n.fReference1 + nnodes) == fProjs.end()) {
- LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fReference1 + nnodes << " not defined.";
- continue;
- }
+ LOG(info) << GetName() << "::Exec : Evs(%)[" << 1.e2 * iev / fEvents->GetEntriesFast() << "] Trks(%)["
+ << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
+ return;
+ }
+ else {
+
+ processHits(nullptr);
- // View* mod = &fProjs[n.fReference1 + nnodes][0];
- // mod->Load(&n);
- // if (mod->hdXx.second)
- // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X(), mod->hdXx.first * n.fTrack.Y());
- // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X(), n.fTrack.Time());
+ if (!fGTracks) {
+ LOG(info) << GetName() << "::Exec : TS local reco only.";
+ return; // END TS QA (no tracking)
}
- itrack++;
- } // END TRACK LOOP
- LOG(info) << GetName() << "::Exec : Trks(%)[" << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
- return; // END TS QA (tracking)
+ processTracks(nullptr);
+ }
+ LOG(info) << GetName() << "::Exec : Trks(%)[" << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
+ return;
}
//_____________________________________________________________________
@@ -555,7 +610,8 @@ void CbmRecoQaTask::SetProjections(std::vector<double> vzpj)
{
fProjs.clear();
for (auto zprj : vzpj) {
- // fProjView[zprj] = new TH2D(Form("hxx"), Form("; X (cm); Y (cm)"), 450, -25, 20, 450, -25, 20);
+ // fProjView[zprj] = new TH2D(Form("hxx"), Form("; X (cm); Y (cm)"), 450,
+ // -25, 20, 450, -25, 20);
fProjs.emplace(zprj, Detector::View(Form("z%.2f", zprj), "", {}));
}
}
@@ -585,6 +641,57 @@ bool CbmRecoQaTask::FilterTrack(const CbmGlobalTrack* ptr)
return true;
}
+//_____________________________________________________________________
+TString CbmRecoQaTask::GetGeoTagForDetector(const TString& detector)
+{
+ gGeoManager->CdTop();
+ TGeoNode* cave = gGeoManager->GetCurrentNode();
+ if (!cave) {
+ LOG(error) << "Error: Could not get the top node in the geometry manager." << std::endl;
+ return "";
+ }
+
+ for (Int_t iNode = 0; iNode < cave->GetNdaughters(); iNode++) {
+ TString name = cave->GetDaughter(iNode)->GetVolume()->GetName();
+ if (name.Contains(detector, TString::kIgnoreCase)) {
+ return name.Contains("mcbm", TString::kIgnoreCase) ? TString(name(5, name.Length()))
+ : TString(name(5, name.Length() - 5));
+ }
+ }
+ return "";
+}
+
+//_____________________________________________________________________
+std::vector<TString> CbmRecoQaTask::GetPath(TGeoNode* node, TString detector, TString activeNodeName, int depth,
+ const TString& path)
+{
+ std::vector<TString> nodePaths;
+
+ if (!node) {
+ return nodePaths;
+ }
+
+ TString nodePath = path;
+ if (!nodePath.IsNull()) {
+ nodePath += "/";
+ }
+ nodePath += node->GetName();
+
+ if (TString(node->GetName()).Contains(activeNodeName)) {
+ if (nodePath.Contains(detector)) nodePaths.push_back(nodePath);
+ }
+
+ // Recursively traverse the daughters
+ Int_t numDaughters = node->GetNdaughters();
+ for (Int_t i = 0; i < numDaughters; ++i) {
+ TGeoNode* daughterNode = node->GetDaughter(i);
+
+ std::vector<TString> result = GetPath(daughterNode, detector, activeNodeName, depth + 1, nodePath);
+ nodePaths.insert(nodePaths.end(), result.begin(), result.end());
+ }
+
+ return nodePaths;
+}
//_____________________________________________________________________
void CbmRecoQaTask::InitMcbm22()
{
@@ -594,130 +701,113 @@ void CbmRecoQaTask::InitMcbm22()
return;
}
- TString dtag;
- Detector::View* v(nullptr);
+ Detector::View* v = nullptr;
+
+ std::vector<TString> path;
+ TGeoNode* topNode = gGeoManager->GetTopNode();
+ if (!topNode) {
+ LOG(error) << "Error: Top node not found.";
+ return;
+ }
+
+ auto processDetector = [&](const std::string& detector, const std::string& component) {
+ TString geoTag = GetGeoTagForDetector(detector);
+
+ if (geoTag.Length() > 0) {
+ LOG(info) << detector << ": geometry tag is " << geoTag;
+ }
+ else {
+ LOG(warn) << "Warning: No geometry tag found for detector " << detector;
+ }
+
+ path = GetPath(topNode, detector, component, 0);
+ };
+
if (setup->IsActive(ECbmModuleId::kSts)) {
- dtag = "v22f_mcbm";
+
+ processDetector("sts", "Sensor");
+
+ std::regex pattern("/Station(\\d+)_(\\d+)/Ladder(\\d+)_(\\d+)/HalfLadder\\d+d_(\\d+)/"
+ "HalfLadder\\d+d_Module(\\d+)_(\\d+)/Sensor(\\d+)_(\\d+)");
+
Detector* sts = AddDetector(ECbmModuleId::kSts);
- // U0 L0 M0
- v = sts->AddView(
- "U0L0M0",
- Form("/cave_1/sts_%s_0/Station01_1/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
- {0, 0, 0});
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v->SetProjection(eViewProjection::kYdY, 4, "mm");
- v = sts->AddView(
- "U0L0M1",
- Form("/cave_1/sts_%s_0/Station01_1/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
- {0, 0, 1});
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v->SetProjection(eViewProjection::kYdY, 4, "mm");
- // U0 L1
- v = sts->AddView(
- "U0L1M0",
- Form("/cave_1/sts_%s_0/Station01_1/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
- {0, 1, 0});
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v->SetProjection(eViewProjection::kYdY, 4, "mm");
- v = sts->AddView(
- "U0L1M1",
- Form("/cave_1/sts_%s_0/Station01_1/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
- {0, 1, 1});
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v->SetProjection(eViewProjection::kYdY, 4, "mm");
- // U1 L0
- v = sts->AddView(
- "U1L0M0",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module03_1/Sensor03_1", dtag.Data()),
- {1, 0, 0});
- v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
- v->SetProjection(eViewProjection::kYdY, 3, "mm");
- v = sts->AddView(
- "U1L0M1",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module04_2/Sensor04_1", dtag.Data()),
- {1, 0, 1});
- v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
- v->SetProjection(eViewProjection::kYdY, 3, "mm");
- // U1 L1
- v = sts->AddView(
- "U1L1M0",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module03_1/Sensor03_1", dtag.Data()),
- {1, 1, 0});
- v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
- v->SetProjection(eViewProjection::kYdY, 3, "mm");
- v = sts->AddView(
- "U1L1M1",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module04_2/Sensor04_1", dtag.Data()),
- {1, 1, 1});
- v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
- v->SetProjection(eViewProjection::kYdY, 3, "mm");
- // U1 L2
- v = sts->AddView(
- "U1L2M0",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_1/Sensor03_1", dtag.Data()),
- {1, 2, 0});
- v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
- v->SetProjection(eViewProjection::kYdY, 3, "mm");
- v = sts->AddView(
- "U1L2M1",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1", dtag.Data()),
- {1, 2, 1});
- v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
- v->SetProjection(eViewProjection::kYdY, 3, "mm");
- v = sts->AddView(
- "U1L2M2",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1", dtag.Data()),
- {1, 2, 2});
- v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
- v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ for (const auto& str : path) {
+ std::string Str(str.Data());
+ std::smatch match;
+ if (std::regex_search(Str, match, pattern)) {
+ int station = std::stoi(match[2]);
+ int ladder = std::stoi(match[4]);
+ int module = std::stoi(match[7]);
+
+ const char* fType = Form("U%dL%dM%d", station - 1, ladder - 1, module - 1);
+
+ v = sts->AddView(fType, str.Data(), {station - 1, ladder - 1, module - 1});
+ v->SetProjection(eViewProjection::kXdX, (station == 1) ? 1 : 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, (station == 1) ? 4 : 3, "mm");
+ }
+ else {
+ std::cout << "No match found in string: " << str << std::endl;
+ }
+ }
}
+
if (setup->IsActive(ECbmModuleId::kTrd)) {
- dtag = "v22h_mcbm";
+ processDetector("trd", "module");
+
+ // Trd2d
Detector* trd = AddDetector(ECbmModuleId::kTrd2d);
- // Trd2D
- v = trd->AddView("2D", Form("/cave_1/trd_%s_0/layer01_20101/module9_101001001/gas_0", dtag.Data()), {5});
- v->SetProjection(eViewProjection::kChdT, 400, "ns");
- v->SetProjection(eViewProjection::kXdX, 10, "mm");
- v->SetProjection(eViewProjection::kYdY, 20, "mm");
- // Trd1Dx
+
+ for (const auto& str : path) {
+ if (!str.Contains("module9")) continue;
+ v = trd->AddView("2D", str.Data(), {5});
+ v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ v->SetProjection(eViewProjection::kXdX, 10, "mm");
+ v->SetProjection(eViewProjection::kYdY, 20, "mm");
+ }
+ // Trd1DxDy
trd = AddDetector(ECbmModuleId::kTrd);
- v = trd->AddView("1Dx", Form("/cave_1/trd_%s_0/layer02_10202/module8_101002001", dtag.Data()), {21});
- v->SetProjection(eViewProjection::kChdT, 400, "ns");
- v->SetProjection(eViewProjection::kXdX, 1.5, "cm");
- // Trd1Dy
- v = trd->AddView("1Dy", Form("/cave_1/trd_%s_0/layer03_11303/module8_101303001", dtag.Data()), {37});
- v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ for (const auto& str : path) {
+ LOG(info) << str;
+ if (str.Contains("layer02")) {
+ v = trd->AddView("1Dx", str.Data(), {21});
+ v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ v->SetProjection(eViewProjection::kXdX, 1.5, "cm");
+ }
+ if (str.Contains("layer03")) {
+ v = trd->AddView("1Dy", str.Data(), {37});
+ v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ }
+ }
}
+
if (setup->IsActive(ECbmModuleId::kTof)) {
- dtag = "v21k_mcbm";
+ processDetector("tof", "counter");
+
Detector* tof = AddDetector(ECbmModuleId::kTof);
- tof->hit.name = "TofCalHit";
- vector<int> tofSelect(3);
- // add type 0
- tofSelect[1] = 0; // ToF type
- for (int ism(0); ism < 4; ism++) {
- tofSelect[0] = ism;
- for (int irpc(0); irpc < 5; irpc++) {
- tofSelect[2] = irpc;
- v = tof->AddView(Form("Sm%dRpc%d", tofSelect[0], tofSelect[2]),
- Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(),
- dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
- tofSelect);
- v->SetProjection(eViewProjection::kChdT, 15, "ns");
+ tof->hit.name = "TofHit";
+ std::regex pattern("module_(\\d+)_(\\d+)/gas_box_(\\d+)/counter_(\\d+)");
+
+ for (const auto& str : path) {
+ std::string Str(str.Data());
+ std::smatch match;
+ if (std::regex_search(Str, match, pattern)) {
+ int type = std::stoi(match[1]);
+ int smid = std::stoi(match[2]);
+ int rpc = std::stoi(match[4]);
+ const char* name = Form("Sm%d_%dRpc%d", type, smid, rpc);
+ if (type == 0) {
+ v = tof->AddView(name, str.Data(), {smid, type, rpc});
+ v->SetProjection(eViewProjection::kChdT, 15, "ns");
+ }
+ else {
+ v = tof->AddView(name, str.Data(), {smid, type, rpc});
+ }
+ }
+ else {
+ std::cout << "No match found in string: " << str << std::endl;
}
- }
- // add type 2
- tofSelect[1] = 2; // ToF type
- tofSelect[0] = 0; // Tof SM id
- for (int irpc(0); irpc < 5; irpc++) {
- tofSelect[2] = irpc;
- tof->AddView(Form("Sm_2%dRpc%d", tofSelect[0], tofSelect[2]),
- Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
- tofSelect[1], tofSelect[0], tofSelect[2]),
- tofSelect);
}
}
-
// ===============================================================================
// TRG - upstream projections
// int iprjLy = 0;
@@ -726,11 +816,13 @@ void CbmRecoQaTask::InitMcbm22()
// prj.second.name = Form("Z=%+4.1f", prj.first);
// prj.second.hdXx.first = 1;
// prj.second.hdXx.second =
- // new TH2D(Form("Hxy%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk} (cm); Y_{trk} (cm)", prj.first),
+ // new TH2D(Form("Hxy%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk}
+ // (cm); Y_{trk} (cm)", prj.first),
// 700, -25, 10, 100, -5, 5);
// prj.second.hdT.first = 1;
// prj.second.hdT.second =
- // new TH2D(Form("Hxt%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk} (cm); T_{trk} - T_{ev} (ns)", prj.first),
+ // new TH2D(Form("Hxt%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk}
+ // (cm); T_{trk} - T_{ev} (ns)", prj.first),
// 700, -25, 10, 350, -15.5, 5.5);
// iprjLy++;
// }
@@ -752,79 +844,91 @@ void CbmRecoQaTask::InitMcbm24()
dtag = "v24c_mcbm";
Detector* sts = AddDetector(ECbmModuleId::kSts);
// U0 L0 M0
- v = sts->AddView(
- "U0L0M0",
- Form("/cave_1/sts_%s_0/Station01_1/Ladder13_1/HalfLadder13u_1/HalfLadder13u_Module03_1/Sensor03_1", dtag.Data()),
- {0, 0, 0});
+ v = sts->AddView("U0L0M0",
+ Form("/cave_1/sts_%s_0/Station01_1/Ladder13_1/"
+ "HalfLadder13u_1/HalfLadder13u_Module03_1/Sensor03_1",
+ dtag.Data()),
+ {0, 0, 0});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
// U1 L0 M0
- v = sts->AddView(
- "U1L0M0",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
- {1, 0, 0});
+ v = sts->AddView("U1L0M0",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_1/"
+ "HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1",
+ dtag.Data()),
+ {1, 0, 0});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v = sts->AddView(
- "U1L0M1",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
- {1, 0, 1});
+ v = sts->AddView("U1L0M1",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_1/"
+ "HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1",
+ dtag.Data()),
+ {1, 0, 1});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
// U1 L1
- v = sts->AddView(
- "U1L1M0",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
- {1, 1, 0});
+ v = sts->AddView("U1L1M0",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_2/"
+ "HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1",
+ dtag.Data()),
+ {1, 1, 0});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v = sts->AddView(
- "U1L1M1",
- Form("/cave_1/sts_%s_0/Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
- {1, 1, 1});
+ v = sts->AddView("U1L1M1",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_2/"
+ "HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1",
+ dtag.Data()),
+ {1, 1, 1});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
// U2 L0
- v = sts->AddView(
- "U2L0M0",
- Form("/cave_1/sts_%s_0/Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module03_1/Sensor03_1", dtag.Data()),
- {2, 0, 0});
+ v = sts->AddView("U2L0M0",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder10_1/"
+ "HalfLadder10d_2/HalfLadder10d_Module03_1/Sensor03_1",
+ dtag.Data()),
+ {2, 0, 0});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v = sts->AddView(
- "U2L0M1",
- Form("/cave_1/sts_%s_0/Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module04_2/Sensor04_1", dtag.Data()),
- {2, 0, 1});
+ v = sts->AddView("U2L0M1",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder10_1/"
+ "HalfLadder10d_2/HalfLadder10d_Module04_2/Sensor04_1",
+ dtag.Data()),
+ {2, 0, 1});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
// U2 L1
- v = sts->AddView(
- "U2L1M0",
- Form("/cave_1/sts_%s_0/Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module03_1/Sensor03_1", dtag.Data()),
- {2, 1, 0});
+ v = sts->AddView("U2L1M0",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder12_2/"
+ "HalfLadder12d_2/HalfLadder12d_Module03_1/Sensor03_1",
+ dtag.Data()),
+ {2, 1, 0});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v = sts->AddView(
- "U2L1M1",
- Form("/cave_1/sts_%s_0/Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module04_2/Sensor04_1", dtag.Data()),
- {2, 1, 1});
+ v = sts->AddView("U2L1M1",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder12_2/"
+ "HalfLadder12d_2/HalfLadder12d_Module04_2/Sensor04_1",
+ dtag.Data()),
+ {2, 1, 1});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
// U2 L2
- v = sts->AddView(
- "U2L2M0",
- Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_1/Sensor03_1", dtag.Data()),
- {2, 2, 0});
+ v = sts->AddView("U2L2M0",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/"
+ "HalfLadder11d_2/HalfLadder11d_Module03_1/Sensor03_1",
+ dtag.Data()),
+ {2, 2, 0});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
- v = sts->AddView(
- "U2L2M1",
- Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1", dtag.Data()),
- {2, 2, 1});
- v = sts->AddView(
- "U2L2M2",
- Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1", dtag.Data()),
- {2, 2, 2});
+ v = sts->AddView("U2L2M1",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/"
+ "HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1",
+ dtag.Data()),
+ {2, 2, 1});
+ v = sts->AddView("U2L2M2",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/"
+ "HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1",
+ dtag.Data()),
+ {2, 2, 2});
v->SetProjection(eViewProjection::kYdY, 5, "mm");
v->SetProjection(eViewProjection::kXdX, 1, "mm");
}
@@ -858,8 +962,9 @@ void CbmRecoQaTask::InitMcbm24()
for (int irpc(0); irpc < 5; irpc++) {
tofSelect[2] = irpc;
v = tof->AddView(Form("Sm%dRpc%d", ism, irpc),
- Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(),
- dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/"
+ "gas_box_0/counter_%d",
+ dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
tofSelect);
v->SetProjection(eViewProjection::kXdX, 7, "cm");
v->SetProjection(eViewProjection::kYdY, 7, "cm");
@@ -871,8 +976,9 @@ void CbmRecoQaTask::InitMcbm24()
for (int irpc(0); irpc < 2; irpc++) {
tofSelect[2] = irpc;
tof->AddView(Form("BuchRpc%d", irpc),
- Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
- tofSelect[1], tofSelect[0], tofSelect[2]),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/"
+ "counter_%d",
+ dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
tofSelect);
}
// add type 9
@@ -882,8 +988,9 @@ void CbmRecoQaTask::InitMcbm24()
for (int irpc(0); irpc < 2; irpc++) {
tofSelect[2] = irpc;
tof->AddView(Form("Test%dRpc%d", ism, irpc),
- Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
- tofSelect[1], tofSelect[0], tofSelect[2]),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/"
+ "gas_box_0/counter_%d",
+ dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
tofSelect);
}
}
@@ -894,8 +1001,9 @@ void CbmRecoQaTask::InitMcbm24()
for (int irpc(0); irpc < 5; irpc++) {
tofSelect[2] = irpc;
tof->AddView(Form("Sm2%dRpc%d", ism, irpc),
- Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
- tofSelect[1], tofSelect[0], tofSelect[2]),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/"
+ "gas_box_0/counter_%d",
+ dtag.Data(), dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
tofSelect);
}
}
@@ -909,30 +1017,165 @@ void CbmRecoQaTask::InitMcbm24()
rich->AddView("Aerogel", Form("/cave_1/rich_%s_0/box_1/Gas_1", dtag.Data()), {4});
}
}
+//_____________________________________________________________________
+void CbmRecoQaTask::InitDefault()
+{
+ LOG(info) << "Init Default ....";
+ CbmSetup* setup = CbmSetup::Instance();
+ if (!setup) {
+ LOG(fatal) << GetName() << "::InitDefault() : Missing setup definition.";
+ return;
+ }
+
+ Detector::View* v = nullptr;
+ std::vector<TString> path;
+ TGeoNode* topNode = gGeoManager->GetTopNode();
+ if (!topNode) {
+ LOG(error) << "Error: Top node not found.";
+ return;
+ }
+
+ TString geoTag;
+ auto processDetector = [&](const std::string& detector, const std::string& component) {
+ geoTag = GetGeoTagForDetector(detector);
+
+ if (geoTag.Length() > 0) {
+ LOG(info) << detector << ": geometry tag is " << geoTag;
+ }
+ else {
+ LOG(warn) << "Warning: No geometry tag found for detector " << detector;
+ }
+ path = GetPath(topNode, detector, component, 0);
+ };
+
+ if (setup->IsActive(ECbmModuleId::kSts)) {
+ processDetector("sts", "Unit");
+
+ // std::regex pattern("Unit\\d{2}[LR]_\\d+");
+ std::regex pattern(R"(Unit(\d{2}[LR])_(\d+))");
+
+ Detector* sts = AddDetector(ECbmModuleId::kSts);
+
+ for (const auto& str : path) {
+ std::string Str(str.Data());
+ std::smatch match;
+
+ if (std::regex_search(Str, match, pattern)) {
+ int unitid = std::stoi(match[2].str());
+ std::string unitname = (Str.find("Unit") != std::string::npos) ? Str.substr(Str.find("Unit")) : "";
+
+ v = sts->AddView(unitname.c_str(), str.Data(), {unitid - 1, -1, -1});
+ v->SetMode(eSetup::kDefault);
+ }
+ else {
+ std::cout << "No match found in string: " << str << std::endl;
+ }
+ }
+ }
+
+ if (setup->IsActive(ECbmModuleId::kTrd)) {
+ processDetector("trd", "module");
+ std::regex pattern("layer(\\d+)_(\\d+)/module(\\d+)_(\\d+)");
+ Detector* trd = AddDetector(ECbmModuleId::kTrd);
+ char name[256];
+ for (const auto& str : path) {
+ std::string Str(str.Data());
+ std::smatch match;
+ if (std::regex_search(Str, match, pattern)) {
+ int layer = std::stoi(match[1]);
+ int layercopyNr = std::stoi(match[2]);
+
+ int module = std::stoi(match[3]);
+ int modulecopyNr = std::stoi(match[4]);
+
+ int fLayer = ((layercopyNr / 100) % 10);
+ int fModuleCopy = (modulecopyNr % 100);
+
+ sprintf(name, "layer%d_%d_module_%d_%d", layer, layercopyNr, module, modulecopyNr);
+
+ v = trd->AddView(name, str.Data(), {fLayer - 1, fModuleCopy - 1, -1});
+ v->SetMode(eSetup::kDefault);
+ }
+ else {
+ std::cout << "No match found in string: " << str << std::endl;
+ }
+ }
+ }
+ if (setup->IsActive(ECbmModuleId::kRich)) {
+ processDetector("rich", "sens");
+ Detector* rich = AddDetector(ECbmModuleId::kRich);
+ for (const auto& str : path) {
+ v = rich->AddView("rich", str.Data(), {-1, -1, -1});
+ v->SetMode(eSetup::kDefault);
+ }
+ }
+
+ if (setup->IsActive(ECbmModuleId::kTof)) {
+ processDetector("tof", "module");
+
+ Detector* tof = AddDetector(ECbmModuleId::kTof);
+ tof->hit.name = "TofHit";
+ std::regex pattern("module_(\\d+)_(\\d+)");
+ for (const auto& str : path) {
+ std::string Str(str.Data());
+ std::smatch match;
+ if (std::regex_search(Str, match, pattern)) {
+ int type = std::stoi(match[1]);
+ int smid = std::stoi(match[2]);
+
+ std::string modulename = (Str.find("module") != std::string::npos) ? Str.substr(Str.find("module")) : "";
+
+ v = tof->AddView(modulename.c_str(), str.Data(), {smid, type, -1});
+ v->SetMode(eSetup::kDefault);
+ }
+ else {
+ std::cout << "No match found in string: " << str << std::endl;
+ }
+ }
+ }
+}
//========= DETECTOR ======================
CbmRecoQaTask::Detector::Detector(ECbmModuleId did)
{
id = did;
switch (id) {
- case ECbmModuleId::kMvd: new (&hit) Data(ECbmDataType::kMvdHit, "MvdHit"); break;
+ case ECbmModuleId::kMvd:
+ new (&hit) Data(ECbmDataType::kMvdHit, "MvdHit");
+ new (&point) Data(ECbmDataType::kMvdPoint, "MvdPoint");
+ break;
case ECbmModuleId::kSts:
new (&hit) Data(ECbmDataType::kStsHit, "StsHit");
+ new (&point) Data(ECbmDataType::kStsPoint, "StsPoint");
new (&trk) Data(ECbmDataType::kStsTrack, "StsTrack");
break;
- case ECbmModuleId::kMuch: new (&hit) Data(ECbmDataType::kMuchPixelHit, "MuchHit"); break;
- case ECbmModuleId::kRich: new (&hit) Data(ECbmDataType::kRichHit, "RichHit"); break;
+ case ECbmModuleId::kMuch:
+ new (&hit) Data(ECbmDataType::kMuchPixelHit, "MuchHit");
+ new (&point) Data(ECbmDataType::kMuchPoint, "MuchPoint");
+ break;
+ case ECbmModuleId::kRich:
+ new (&hit) Data(ECbmDataType::kRichHit, "RichHit");
+ new (&point) Data(ECbmDataType::kRichPoint, "RichPoint");
+ break;
case ECbmModuleId::kTrd:
case ECbmModuleId::kTrd2d:
new (&hit) Data(ECbmDataType::kTrdHit, "TrdHit");
+ new (&point) Data(ECbmDataType::kTrdPoint, "TrdPoint");
new (&trk) Data(ECbmDataType::kTrdTrack, "TrdTrack");
break;
case ECbmModuleId::kTof:
new (&hit) Data(ECbmDataType::kTofHit, "TofHit");
+ new (&point) Data(ECbmDataType::kTofPoint, "TofPoint");
new (&trk) Data(ECbmDataType::kTofTrack, "TofTrack");
break;
- case ECbmModuleId::kFsd: new (&hit) Data(ECbmDataType::kFsdHit, "FsdHit"); break;
- case ECbmModuleId::kPsd: new (&hit) Data(ECbmDataType::kPsdHit, "PsdHit"); break;
+ case ECbmModuleId::kFsd:
+ new (&hit) Data(ECbmDataType::kFsdHit, "FsdHit");
+ new (&point) Data(ECbmDataType::kFsdHit, "FsdPoint");
+ break;
+ case ECbmModuleId::kPsd:
+ new (&hit) Data(ECbmDataType::kPsdHit, "PsdHit");
+ new (&point) Data(ECbmDataType::kPsdPoint, "PsdPoint");
+ break;
default:
LOG(warn) << "QA unsupported for Detector=" << ToString(did);
id = ECbmModuleId::kNotExist;
@@ -951,11 +1194,20 @@ CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::AddView(const char* n, c
v->AddProjection(CbmRecoQaTask::eViewProjection::kChdT);
v->AddProjection(CbmRecoQaTask::eViewProjection::kXYh);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXYhMC);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kPullX);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kPullY);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kResidualX);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kResidualY);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kResidualTX);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kResidualTY);
v->AddProjection(CbmRecoQaTask::eViewProjection::kXYa);
v->AddProjection(CbmRecoQaTask::eViewProjection::kXYp);
v->AddProjection(CbmRecoQaTask::eViewProjection::kXdX);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXdXMC);
v->AddProjection(CbmRecoQaTask::eViewProjection::kYdY);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kYdYMC);
v->AddProjection(CbmRecoQaTask::eViewProjection::kWdT);
v->AddProjection(CbmRecoQaTask::eViewProjection::kXpX);
v->AddProjection(CbmRecoQaTask::eViewProjection::kYpY);
@@ -974,7 +1226,9 @@ CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::GetView(const char* n)
CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::FindView(double x, double y, double z)
{
for (auto& v : fViews) {
- // printf("Looking for p[%f %f %f] in V[%s-%s] @ %f %f %f (%f %f %f)\n", x, y, z, ToString(id).data(), v.name.data(), v.pos[0], v.pos[1], v.pos[2], v.size[0], v.size[1], v.size[2]);
+ // printf("Looking for p[%f %f %f] in V[%s-%s] @ %f %f %f (%f %f %f)\n", x,
+ // y, z, ToString(id).data(), v.name.data(), v.pos[0], v.pos[1], v.pos[2],
+ // v.size[0], v.size[1], v.size[2]);
if (abs(v.pos[0] - x) > 0.5 * v.size[0]) continue;
if (abs(v.pos[1] - y) > 0.5 * v.size[1]) continue;
if (abs(v.pos[2] - z) > 0.5 * v.size[2]) continue;
@@ -987,7 +1241,7 @@ CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::FindView(double x, doubl
bool CbmRecoQaTask::Detector::Init(TDirectoryFile* fOut)
{
/* Query the geometry and trigger detector views init
- */
+ */
if (!gGeoManager) {
LOG(fatal) << "CbmRecoQaTask::Detector::Init() " << ToString(id) << " missing geometry.";
return false;
@@ -1050,7 +1304,8 @@ bool CbmRecoQaTask::Detector::View::SetProjection(eViewProjection prj, float ran
{
if (fProjection.find(prj) == fProjection.end()) {
LOG(warn) << "Projection " << ToString(prj)
- << " not initialized. Calling \"CbmRecoQaTask::Detector::View::AddProjection()\"";
+ << " not initialized. Calling "
+ "\"CbmRecoQaTask::Detector::View::AddProjection()\"";
return AddProjection(prj, range, unit);
}
int scale(1);
@@ -1107,26 +1362,36 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
double xlo = pos[0] - 0.5 * size[0], xhi = pos[0] + 0.5 * size[0], ylo = pos[1] - 0.5 * size[1],
yhi = pos[1] + 0.5 * size[1];
switch (projection.first) {
+
+ case eViewProjection::kXdXMC: xy_id = 'M';
+ // fall through
case eViewProjection::kXdX:
- nbinsx = 10 * ceil(size[0]); // mm binning
- unit = makeYrange(scale, yrange);
- nbinsy = 200; //* ceil(yrange);
- get<2>(projection.second) = new TH2D(Form("hxx_%s%s", dname, name.data()),
- Form("X resolution %s [%s]; X_{%s-%s} (cm); #Delta X (%s)", name.data(),
- dname, dname, name.data(), unit.data()),
- nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
+ if (!xy_id) xy_id = 't';
+ nbinsx = 10 * ceil(size[0]); // mm binning
+ unit = makeYrange(scale, yrange);
+ nbinsy = 200; //* ceil(yrange);
+ get<2>(projection.second) =
+ new TH2D(Form("hxx%s_%s%s", (xy_id == 'M' ? "MC" : ""), dname, name.data()),
+ Form("X resolution %s %s [%s]; X_{%s-%s} (cm); #Delta X (%s)", (xy_id == 'M' ? "MC" : ""),
+ name.data(), dname, dname, name.data(), unit.data()),
+ nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
break;
+ case eViewProjection::kYdYMC: xy_id = 'M';
+ // fall through
case eViewProjection::kYdY:
- nbinsx = 10 * ceil(size[1]); // mm binning
- unit = makeYrange(scale, yrange);
- nbinsy = 200; //* ceil(yrange);
- get<2>(projection.second) = new TH2D(Form("hyy_%s%s", dname, name.data()),
- Form("Y resolution %s [%s]; Y_{%s-%s} (cm); #Delta Y (%s)", name.data(),
- dname, dname, name.data(), unit.data()),
- nbinsx, ylo, yhi, nbinsy, -yrange, yrange);
+ if (!xy_id) xy_id = 't';
+ nbinsx = 10 * ceil(size[1]); // mm binning
+ unit = makeYrange(scale, yrange);
+ nbinsy = 200; //* ceil(yrange);
+ get<2>(projection.second) =
+ new TH2D(Form("hyy%s_%s%s", (xy_id == 'M' ? "MC" : ""), dname, name.data()),
+ Form("Y resolution %s %s [%s]; Y_{%s-%s} (cm); #Delta Y (%s)", (xy_id == 'M' ? "MC" : ""),
+ name.data(), dname, dname, name.data(), unit.data()),
+ nbinsx, ylo, yhi, nbinsy, -yrange, yrange);
break;
+
case eViewProjection::kXpX:
nbinsx = 10 * ceil(size[0]); // mm binning
if (yrange < 0) {
@@ -1189,6 +1454,68 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
name.data(), dname, dname, name.data(), dname, name.data()),
nbinsx, xlo, xhi, nbinsy, ylo, yhi);
break;
+
+ case eViewProjection::kXYhMC:
+ if (!xy_id) xy_id = 'h';
+ nbinsx = dscale * ceil(size[0]); // mm binning
+ nbinsy = dscale * ceil(size[1]);
+ get<2>(projection.second) =
+ new TH2D(Form("hxymc%c_%s%s", xy_id, dname, name.data()),
+ Form("Hit_{%s} %s [%s]; X_{%s-%s} (cm); Y_{%s-%s} (cm)", (xy_id == 'h' ? "mc" : "attach"),
+ name.data(), dname, dname, name.data(), dname, name.data()),
+ nbinsx, xlo, xhi, nbinsy, ylo, yhi);
+ break;
+
+ case eViewProjection::kResidualX:
+ // fall through
+ case eViewProjection::kResidualY: {
+ bool isResidualX = (projection.first == eViewProjection::kResidualX);
+
+ nbinsx = dscale * ceil(size[0]); // mm binning
+ nbinsy = dscale * ceil(size[1]);
+ unit = makeYrange(scale, yrange);
+
+ get<2>(projection.second) = new TH2D(
+ Form("h%sR_%s%s", (isResidualX ? "xx" : "yy"), dname, name.data()),
+ Form("%s resolution %s [%s]; %s_{%s-%s} (cm); #Delta %s (%s)", (isResidualX ? "X" : "Y"), name.data(), dname,
+ (isResidualX ? "X" : "Y"), dname, name.data(), (isResidualX ? "X" : "Y"), unit.data()),
+ nbinsx, (isResidualX ? xlo : ylo), (isResidualX ? xhi : yhi), nbinsy, -yrange, yrange);
+ break;
+ }
+
+ case eViewProjection::kResidualTX:
+ // fall through
+ case eViewProjection::kResidualTY: {
+ bool isResidualTX = (projection.first == eViewProjection::kResidualTX);
+
+ nbinsx = dscale * ceil(size[0]); // mm binning
+ nbinsy = dscale * ceil(size[1]);
+
+ get<2>(projection.second) =
+ new TH2D(Form("h%sR_%s%s", (isResidualTX ? "tx" : "ty"), dname, name.data()),
+ Form("%s resolution %s [%s]; %s_{%s-%s} (cm); #Delta %s (ns)", (isResidualTX ? "T" : "T"),
+ name.data(), dname, (isResidualTX ? "X" : "Y"), dname, name.data(), (isResidualTX ? "T" : "T")),
+ nbinsx, (isResidualTX ? xlo : ylo), (isResidualTX ? xhi : yhi), nbinsy, -yrange, yrange);
+ break;
+ }
+
+ case eViewProjection::kPullX:
+ // fall through
+ case eViewProjection::kPullY: {
+ bool isPullX = (projection.first == eViewProjection::kPullX);
+
+ nbinsx = dscale * ceil(size[0]); // mm binning
+ nbinsy = dscale * ceil(size[1]);
+ unit = makeYrange(scale, yrange);
+
+ get<2>(projection.second) =
+ new TH2D(Form("h%sP_%s%s", (isPullX ? "xx" : "yy"), dname, name.data()),
+ Form("%s pull %s [%s]; %s_{%s-%s} (cm); Pull %s (%s)", (isPullX ? "X" : "Y"), name.data(), dname,
+ (isPullX ? "X" : "Y"), dname, name.data(), (isPullX ? "X" : "Y"), unit.data()),
+ nbinsx, (isPullX ? xlo : ylo), (isPullX ? xhi : yhi), nbinsy, -yrange, yrange);
+ break;
+ }
+
case eViewProjection::kXYp:
nbinsx = dscale * ceil(size[0] * 2.); // mm binning
nbinsy = dscale * ceil(size[1] * 2.);
@@ -1239,7 +1566,9 @@ uint CbmRecoQaTask::Detector::View::Register(TDirectoryFile* fOut)
case eViewProjection::kXYa:
case eViewProjection::kXYp:
case eViewProjection::kXdX:
+ case eViewProjection::kXdXMC:
case eViewProjection::kYdY:
+ case eViewProjection::kYdYMC:
case eViewProjection::kWdT:
case eViewProjection::kXpX:
case eViewProjection::kYpY:
@@ -1247,7 +1576,15 @@ uint CbmRecoQaTask::Detector::View::Register(TDirectoryFile* fOut)
break;
case eViewProjection::kChdT:
case eViewProjection::kXYh:
+ case eViewProjection::kXYhMC:
+ case eViewProjection::kResidualX:
+ case eViewProjection::kResidualY:
+ case eViewProjection::kResidualTX:
+ case eViewProjection::kResidualTY:
+ case eViewProjection::kPullX:
+ case eViewProjection::kPullY:
if (eDir) eDir->Add(hh);
+ if (sDir) sDir->Add(hh);
break;
case eViewProjection::kXYs:
if (sDir) sDir->Add(hh);
@@ -1267,12 +1604,22 @@ string CbmRecoQaTask::Detector::View::ToString(eViewProjection prj)
case eViewProjection::kXYa: s = "x-y (attach)"; break;
case eViewProjection::kXYp: s = "x-y (track)"; break;
case eViewProjection::kXdX: s = "dx-x"; break;
+ case eViewProjection::kXdXMC: s = "dx-x (MC)"; break;
+ case eViewProjection::kResidualX: s = "residual(x)-x:MC"; break;
+ case eViewProjection::kResidualY: s = "residual(y)-y:MC"; break;
+ case eViewProjection::kResidualTX: s = "residual(t)-x:MC"; break;
+ case eViewProjection::kResidualTY: s = "residual(t)-y:MC"; break;
+ case eViewProjection::kPullX: s = "pull(x)-x:MC"; break;
+ case eViewProjection::kPullY: s = "pull(y)-y:MC"; break;
case eViewProjection::kYdY: s = "dy-y"; break;
+ case eViewProjection::kYdYMC: s = "dy-y (MC)"; break;
case eViewProjection::kXpX: s = "pull(x)-x"; break;
case eViewProjection::kYpY: s = "pull(y)-y"; break;
case eViewProjection::kWdT: s = "dt(trk)-w"; break;
case eViewProjection::kChdT: s = "dt(ev)-w"; break;
case eViewProjection::kXYh: s = "x-y (hit)"; break;
+ case eViewProjection::kXYhMC: s = "x-y (point)"; break;
+
case eViewProjection::kXYs: s = "x-y (TS)"; break;
default: LOG(error) << "Detector::View::ToString() : Unknown projection " << int(prj); break;
}
@@ -1428,12 +1775,14 @@ void CbmRecoQaTask::EventFilter::HelpMess() const
switch (fType) {
case eEventCut::kMultTrk:
LOG(info) << ToString();
- LOG(info) << "\tDepends one two values : min and max of the no of global tracks / event.";
+ LOG(info) << "\tDepends one two values : min and max of the no of global "
+ "tracks / event.";
LOG(info) << "\tValue should follow the relation max >= min.";
break;
case eEventCut::kMultHit:
LOG(info) << ToString();
- LOG(info) << "\tDepends one three values : total hits in the following systems (in this order) : STS TRD TOF.";
+ LOG(info) << "\tDepends one three values : total hits in the following "
+ "systems (in this order) : STS TRD TOF.";
LOG(info) << "\tNegative values are interpreted as no-cut.";
break;
default:
@@ -1525,7 +1874,7 @@ bool CbmRecoQaTask::TrackFilter::SetFilter(std::vector<float> cuts)
case eTrackCut::kSts:
if (cuts.size() < 1) {
LOG(warning) << "Improper definition for track filter :\n\t" << ToString() << endl;
- //HelpMess();
+ // HelpMess();
return false;
}
fNSts = int(cuts[0]);
@@ -1533,7 +1882,7 @@ bool CbmRecoQaTask::TrackFilter::SetFilter(std::vector<float> cuts)
case eTrackCut::kTrd:
if (cuts.size() < 1) {
LOG(warning) << "Improper definition for track filter :\n\t" << ToString() << endl;
- //HelpMess();
+ // HelpMess();
return false;
}
fNTrd = int(cuts[0]);
@@ -1541,7 +1890,7 @@ bool CbmRecoQaTask::TrackFilter::SetFilter(std::vector<float> cuts)
case eTrackCut::kTof:
if (cuts.size() < 1) {
LOG(warning) << "Improper definition for track filter :\n\t" << ToString() << endl;
- //HelpMess();
+ // HelpMess();
return false;
}
fNTof = int(cuts[0]);
diff --git a/reco/qa/CbmRecoQaTask.h b/reco/qa/CbmRecoQaTask.h
index 2f219191cea9137862edb379716d0a8924bdfdf6..197264f9d4b87cac65928d5df4ab15aefcab621f 100644
--- a/reco/qa/CbmRecoQaTask.h
+++ b/reco/qa/CbmRecoQaTask.h
@@ -1,6 +1,6 @@
/* Copyright (C) 2024 Hulubei National Institute of Physics and Nuclear Engineering - Horia Hulubei, Bucharest
SPDX-License-Identifier: GPL-3.0-only
- Authors: Alexandru Bercuci [committer] */
+ Authors: Alexandru Bercuci [committer], Omveer Singh */
#ifndef CBMRECOQATASK_H
#define CBMRECOQATASK_H 1
@@ -11,15 +11,19 @@
#include <FairTask.h>
#include <TDirectoryFile.h>
+#include <TGeoNode.h>
#include <bitset>
#include <map>
#include <vector>
#define kRecoQaNConfigs 8
+class FairMCPoint;
class CbmEvent;
class CbmHit;
+class CbmMCDataManager;
+class CbmMCDataArray;
class CbmTimeSlice;
class TClonesArray;
class TH2;
@@ -44,16 +48,26 @@ class CbmRecoQaTask : public FairTask {
};
enum class eViewProjection : int
{
- kXYa = 0, /// X-Y hit coorelation in track filtered data
- kXYp, /// X-Y track projections
- kXdX, /// X to TRK residuals as function of local X in view
- kYdY, /// Y to TRK residuals as function of local Y in view
- kWdT, /// Time to TRK residuals as function of local high resolution coordinate in view
- kXpX, /// X to TRK pulls as function of local X in view
- kYpY, /// Y to TRK pulls as function of local Y in view
- kChdT, /// Time to EV residuals as function of coordinate in view
- kXYh, /// X-Y hit coorelation in local view
- kXYs /// X-Y hit coorelation in time slice
+ kXYa = 0, /// X-Y hit coorelation in track filtered data
+ kXYp, /// X-Y track projections
+ kXdX, /// X to TRK residuals as function of local X in view
+ kXdXMC, /// X to TRK residuals w.r.t MC points
+ kYdY, /// Y to TRK residuals as function of local Y in view
+ kYdYMC, /// Y to TRK residuals w.r.t MC points
+ kWdT, /// Time to TRK residuals as function of local high resolution
+ /// coordinate in view
+ kXpX, /// X to TRK pulls as function of local X in view
+ kYpY, /// Y to TRK pulls as function of local Y in view
+ kChdT, /// Time to EV residuals as function of coordinate in view
+ kXYh, /// X-Y hit coorelation in local view
+ kXYhMC, /// X-Y MC point coorelation in local view (using HitMatch)
+ kPullX, /// Pull distribution X: (RC - MC) / dx_RC
+ kPullY, /// Pull distribution Y:
+ kResidualX, /// Residual distribution X: (x_RC - x_MC) in cm
+ kResidualY, /// Residual distribution Y:
+ kResidualTX, /// Residual distribution T:
+ kResidualTY, /// Residual distribution T:
+ kXYs /// X-Y hit coorelation in time slice
};
struct Detector {
struct Data {
@@ -71,10 +85,10 @@ class CbmRecoQaTask : public FairTask {
virtual ~View() = default;
bool SetProjection(eViewProjection prj, float range, const char* unit);
template<class Hit>
- bool HasAddress(const CbmHit* h, double& x, double& y) const;
+ bool HasAddress(const CbmHit* h, double& x, double& y, double& dx, double& dy) const;
template<class Hit>
- bool Load(const CbmHit* h, const CbmEvent* ev = nullptr);
- bool Load(const CbmKfTrackFitter::FitNode* n);
+ bool Load(const CbmHit* h, const FairMCPoint* point, const CbmEvent* ev);
+ bool Load(const CbmKfTrackFitter::FitNode* n, const FairMCPoint* point);
std::string ToString() const;
static std::string ToString(eViewProjection prj);
@@ -84,14 +98,20 @@ class CbmRecoQaTask : public FairTask {
double pos[3] = {0.}; /// detection element center x0 y0 z0
std::vector<int> fSelector = {}; /// defining subset of the address set for this view
std::map<eViewProjection, std::tuple<int, float, TH2*>> fProjection =
- {}; /// map of projections indexed by their type. Each projection contains also the relative scale [int] wrt to default units (ns, cm, keV) and the range [float].
+ {}; /// map of projections indexed by their type. Each projection
+ /// contains also the relative scale [int] wrt to default units
+ /// (ns, cm, keV) and the range [float].
+ void SetMode(CbmRecoQaTask::eSetup setup) { fSetup = setup; }
+ mutable eSetup fSetup = eSetup::kMcbm22;
+
protected:
bool AddProjection(eViewProjection prj, float range = -1, const char* unit = "cm");
- /** \brief convert geo address into pointer to geometry.
- * By the time of the call the geometry has to be available in memory.
+ /** \brief convert geo address into pointer to geometry.
+ * By the time of the call the geometry has to be available in memory.
* Failing to identify the named physical node will rezult in error */
bool Init(const char* dname);
- /** \brief build directory structure for all projections of current view.*/
+ /** \brief build directory structure for all projections of current
+ * view.*/
uint Register(TDirectoryFile* f);
/** helper functions to estimate the representation (y) axis
* \param[in] scale read-only unit defining parameter
@@ -101,25 +121,28 @@ class CbmRecoQaTask : public FairTask {
std::string makeTrange(const int scale, float& range);
std::string makeYrange(const int scale, float& range);
- ClassDef(CbmRecoQaTask::Detector::View, 1); // A QA significant correlation
+ ClassDef(CbmRecoQaTask::Detector::View,
+ 1); // A QA significant correlation
};
Detector(ECbmModuleId did = ECbmModuleId::kNotExist);
virtual ~Detector() = default;
View* AddView(const char* n, const char* p, std::vector<int> set);
View* GetView(const char* n);
View* FindView(double x, double y, double z);
- /** \brief Check geometry and trigger Init() for all registered views. Build main directory outut for the current detector.
- * Failing to identify the geometry will rezult in fatal error.
- * \return true If ALL the subsequent calls to Init result in a true */
+ /** \brief Check geometry and trigger Init() for all registered views. Build
+ * main directory outut for the current detector. Failing to identify the
+ * geometry will rezult in fatal error. \return true If ALL the subsequent
+ * calls to Init result in a true */
bool Init(TDirectoryFile* f);
void Print() const;
ECbmModuleId id;
- Data hit;
+ Data hit, point;
Data trk;
std::vector<View> fViews = {};
- ClassDef(CbmRecoQaTask::Detector, 1); // QA representation of a detector unit
+ ClassDef(CbmRecoQaTask::Detector,
+ 1); // QA representation of a detector unit
};
struct TrackFilter {
@@ -188,7 +211,8 @@ class CbmRecoQaTask : public FairTask {
public:
CbmRecoQaTask();
virtual ~CbmRecoQaTask() { ; }
- /** Copy the qa test defined for detector det from the steering macro to the current class */
+ /** Copy the qa test defined for detector det from the steering macro to the
+ * current class */
virtual Detector* AddDetector(ECbmModuleId did);
virtual EventFilter* AddEventFilter(EventFilter::eEventCut cut);
virtual TrackFilter* AddTrackFilter(TrackFilter::eTrackCut cut);
@@ -206,9 +230,13 @@ class CbmRecoQaTask : public FairTask {
virtual void Finish();
- /** \brief Define the set of extra z positions where the track should be projected in the x-y plane */
+ /** \brief Define the set of extra z positions where the track should be
+ * projected in the x-y plane */
void SetProjections(std::vector<double> vzpj);
void SetSetupClass(CbmRecoQaTask::eSetup setup) { fSetupClass = setup; }
+ TString GetGeoTagForDetector(const TString& detector);
+ std::vector<TString> GetPath(TGeoNode* node, TString, TString activeNodeName, int depth = 0,
+ const TString& path = "");
protected:
static std::bitset<kRecoQaNConfigs> fuRecoConfig;
@@ -217,29 +245,33 @@ class CbmRecoQaTask : public FairTask {
CbmRecoQaTask(const CbmRecoQaTask&);
CbmRecoQaTask& operator=(const CbmRecoQaTask&);
- /** \brief Filter events for QA use (e.g. event multiplicity)
+ /** \brief Filter events for QA use (e.g. event multiplicity)
* \param[in] ptr cbm event
* \return true if event accepted
*/
virtual bool FilterEvent(const CbmEvent* ptr);
- /** \brief Filter tracks for further use (e.g. track projections)
+ /** \brief Filter tracks for further use (e.g. track projections)
* \param[in] ptr global track
* \return true if track accepted
*/
virtual bool FilterTrack(const CbmGlobalTrack* ptr);
void InitMcbm22();
void InitMcbm24();
+ void InitDefault();
CbmKfTrackFitter fFitter;
- TClonesArray* fGTracks = nullptr; //! reconstructed global tracks / event
- std::map<ECbmModuleId, TClonesArray*> fTracks = {}; //! reconstructed global tracks / event
- TClonesArray* fTrackMatches = nullptr; //! MC info for the global tracks
- TClonesArray* fEvents = nullptr; //! reconstructed events
- CbmTimeSlice* fTimeSlice = nullptr; //! Time slice info
- std::map<ECbmModuleId, TClonesArray*> fHits = {}; //! reconstructed hits
- std::map<ECbmModuleId, Detector> fDetQa = {};
- std::vector<EventFilter> fFilterEv = {};
- std::vector<TrackFilter> fFilterTrk = {};
+ TClonesArray* fGTracks = nullptr; //! reconstructed global tracks / event
+ std::map<ECbmModuleId, TClonesArray*> fTracks = {}; //! reconstructed global tracks / event
+ TClonesArray* fTrackMatches = nullptr; //! MC info for the global tracks
+ TClonesArray* fEvents = nullptr; //! reconstructed events
+ CbmTimeSlice* fTimeSlice = nullptr; //! Time slice info
+ std::map<ECbmModuleId, TClonesArray*> fHits = {}; //! reconstructed hits
+ std::map<ECbmModuleId, CbmMCDataArray*> fPoints = {}; //! mc points
+ std::map<ECbmModuleId, TClonesArray*> fHitMatch = {}; //! reconstructed hits
+ std::map<ECbmModuleId, Detector> fDetQa = {};
+ std::vector<EventFilter> fFilterEv = {};
+ std::vector<TrackFilter> fFilterTrk = {};
+ CbmMCDataManager* cbm_mc_manager{nullptr};
TDirectoryFile fOutFolder = {"RecoQA", "CA track driven reco QA"}; //!
eSetup fSetupClass = eSetup::kMcbm24;