diff --git a/reco/qa/CbmRecoQaTask.cxx b/reco/qa/CbmRecoQaTask.cxx
index cbca9816885c980447e8e80557bb27bca04909fc..259210dc0bc8d04b1caff856e418143974fb52bb 100644
--- a/reco/qa/CbmRecoQaTask.cxx
+++ b/reco/qa/CbmRecoQaTask.cxx
@@ -53,17 +53,15 @@ InitStatus CbmRecoQaTask::Init()
}
LOG(info) << GetName() << "::Init: Registered projections " << fProjView.size();
- for (auto zprj : fProjView)
- LOG(debug) << GetName() << "::Init: Project trks @ " << zprj.first;
- int itrkLy(0);
+ int itrkLy(0), iprjLy(0);
TFolder *modDir(nullptr), *smDir(nullptr);
fOutFolder.Clear();
switch (fSetupClass) {
case kMcbm22:
LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2022\".";
- modDir = fOutFolder.AddFolder("STS", "STS reziduals to CA");
// STS U0 - station 0
+ modDir = fOutFolder.AddFolder("STS", "STS reziduals to CA");
fModView[itrkLy].push_back(ModuleView());
fModView[itrkLy][0].hdXx.first = 10;
fModView[itrkLy][0].hdXx.second =
@@ -105,6 +103,7 @@ InitStatus CbmRecoQaTask::Init()
new TH2D(Form("hpully%d", itrkLy), Form("STS U1 pully-tht; dy/dz_{trk}; pull Y"), 160, -0.8, 0.8, 100, -5, 5);
modDir->Add(fModView[itrkLy][0].hpullYtht);
itrkLy++;
+
// ===============================================================================
modDir = fOutFolder.AddFolder("TRD", "TRD reziduals to CA");
// TRD 2D - station 2
@@ -170,10 +169,11 @@ InitStatus CbmRecoQaTask::Init()
new TH2D(Form("hpully%d", itrkLy), Form("TRD1Dy pully-tht; dy/dz_{trk}; pull Y"), 100, -0.5, 0.5, 100, -5, 5);
modDir->Add(fModView[itrkLy][0].hpullYtht);
itrkLy++;
+
// ===============================================================================
- modDir = fOutFolder.AddFolder("ToF", "ToF reziduals to CA");
// ToF Sm0 & Sm3 (Typ0) - station 5
- smDir = modDir->AddFolder("SM0_0", "ToF reziduals to CA for Sm0 and Sm3 type 0");
+ modDir = fOutFolder.AddFolder("ToF", "ToF reziduals to CA");
+ smDir = modDir->AddFolder("SM0_0", "ToF reziduals to CA for Sm0 and Sm3 type 0");
for (int irpc(0); irpc < 5; irpc++) {
fModView[itrkLy].push_back(ModuleView());
ModuleView* cdet = &fModView[itrkLy][irpc];
@@ -202,7 +202,6 @@ InitStatus CbmRecoQaTask::Init()
smDir->Add(cdet->hpullYtht);
}
itrkLy++;
-
// ToF Sm1 (Typ0) - station 6
smDir = modDir->AddFolder("SM1_0", "ToF reziduals to CA for Sm1 type 0");
for (int irpc(0); irpc < 5; irpc++) {
@@ -233,7 +232,6 @@ InitStatus CbmRecoQaTask::Init()
smDir->Add(cdet->hpullYtht);
}
itrkLy++;
-
// ToF Sm0 (Typ2) - station 7
smDir = modDir->AddFolder("SM0_2", "ToF reziduals to CA for Sm0 type 2");
for (int irpc(0); irpc < 5; irpc++) {
@@ -264,7 +262,6 @@ InitStatus CbmRecoQaTask::Init()
smDir->Add(cdet->hpullYtht);
}
itrkLy++;
-
// ToF Sm2 (Typ0) - station 8
smDir = modDir->AddFolder("SM2_0", "ToF reziduals to CA for Sm2 type 0");
for (int irpc(0); irpc < 5; irpc++) {
@@ -294,6 +291,27 @@ InitStatus CbmRecoQaTask::Init()
-0.45, 0.45, 100, -5, 5);
smDir->Add(cdet->hpullYtht);
}
+ itrkLy++;
+
+ // ===============================================================================
+ // TRG - upstream projections
+ modDir = fOutFolder.AddFolder("Target", "Target tomography with CA");
+ iprjLy = itrkLy;
+ for (auto zprj : fProjView) {
+ LOG(debug) << GetName() << "::Init: Project trks @ " << zprj.first;
+ fModView[iprjLy].push_back(ModuleView());
+ fModView[iprjLy][0].hdXx.first = 1;
+ fModView[iprjLy][0].hdXx.second =
+ new TH2D(Form("hxy%d", iprjLy - itrkLy), Form("Project @ z = %+5.2f; X_{trk} (cm); Y_{trk} (cm)", zprj.first),
+ 700, -25, 10, 100, -5, 5);
+ modDir->Add(fModView[iprjLy][0].hdXx.second);
+ fModView[iprjLy][0].hdYy.first = 1;
+ fModView[iprjLy][0].hdYy.second =
+ new TH2D(Form("hxt%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk} (cm); T_{trk} - T_{ev} (ns)", zprj.first),
+ 700, -25, 10, 350, -15.5, 5.5);
+ modDir->Add(fModView[iprjLy][0].hdYy.second);
+ iprjLy++;
+ }
break;
case kMcbm24: LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2024\". Not implemented."; break;
case kDefault: LOG(info) << GetName() << "::Init: Setup config for \"DEFAULT\". Not implemented."; break;
@@ -307,7 +325,7 @@ void CbmRecoQaTask::Exec(Option_t*)
{
LOG(debug) << GetName() << "::Exec : Evs[" << (fEvents ? fEvents->GetEntriesFast() : 0) << "] Trks["
<< fTracks->GetEntriesFast() << "]";
- int itrack(0);
+ int itrack(0), nnodes(0);
for (auto trackObj : *fTracks) {
//printf("GlobalTrack %d\n", itrack);
auto track = dynamic_cast<const CbmGlobalTrack*>(trackObj);
@@ -318,15 +336,12 @@ void CbmRecoQaTask::Exec(Option_t*)
}
trkKf.fMsQp0 = 1. / 0.5;
trkKf.fIsMsQp0Set = true;
+ 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) {
if (n.fHitSystemId == ECbmModuleId::kNotExist) continue;
- //printf(" %d z[%6.2f]\n", n.fMaterialLayer, n.fZ);
- //printf(" D[%d] hit[%d %d] Flags[%c %c %c %c]\n",
- // (int)n.fHitSystemId, n.fHitAddress, n.fHitIndex,
- // (n.fIsXySet ? 'y' : 'n'), (n.fIsTimeSet ? 'y' : 'n'), (n.fIsFitted ? 'y' : 'n'), (n.fIsRadThickSet ? 'y' : 'n'));
- //printf(" hit x[%+6.2f] y[%+6.2f] t[%f]\n", (double)n.fMxy.X()[0], (double)n.fMxy.Y()[0], (double)n.fMt.T()[0]);
-
n.fIsTimeSet = false;
n.fIsXySet = false;
trkKf.MakeConsistent();
@@ -334,20 +349,33 @@ void CbmRecoQaTask::Exec(Option_t*)
double dx = n.fMxy.X()[0] - n.fTrack.X()[0], dy = n.fMxy.Y()[0] - n.fTrack.Y()[0],
dt = n.fMt.T()[0] - n.fTrack.Time()[0],
pullx = dx / sqrt(n.fMxy.Dx2()[0] + n.fTrack.GetCovariance(0, 0)[0]),
- pully = dy / sqrt(n.fMxy.Dy2()[0] + n.fTrack.GetCovariance(1, 1)[0]),
- pullt = dt / sqrt(n.fMt.Dt2()[0] + n.fTrack.GetCovariance(5, 5)[0]);
- //printf(" res dx[%+6.4f] dy[%+6.4f] dt[%+6.4f]\n pull x[%+6.4f] y[%+6.4f] t[%+6.4f]\n", dx, dy, dt, pullx, pully, pullt);
+ pully = dy / sqrt(n.fMxy.Dy2()[0] + n.fTrack.GetCovariance(1, 1)[0]);
+ // pullt = dt / sqrt(n.fMt.Dt2()[0] + n.fTrack.GetCovariance(5, 5)[0]);
// add granularity for specific systems
int modId = 0;
- if (n.fHitSystemId == ECbmModuleId::kTof) modId = CbmTofAddress::GetRpcId(n.fHitAddress);
+ switch (n.fHitSystemId) {
+ case ECbmModuleId::kTof:
+ modId = CbmTofAddress::GetRpcId(n.fHitAddress);
+ nhit[(int) n.fHitSystemId]++;
+ break;
+ default: nhit[(int) n.fHitSystemId]++; break;
+ }
// check that we have correctly created all views
- if (fModView.find(n.fMaterialLayer) == fModView.end()) continue;
- if (fModView[n.fMaterialLayer].size() <= (size_t) modId) continue;
+ if (fModView.find(n.fMaterialLayer) == fModView.end()) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fMaterialLayer << " not defined.";
+ continue;
+ }
+ if (fModView[n.fMaterialLayer].size() <= (size_t) modId) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fMaterialLayer << " @ " << modId
+ << " not defined.";
+ continue;
+ }
ModuleView* mod = &fModView[n.fMaterialLayer][modId];
if (!mod) {
- LOG(warn) << GetName() << "::Exec: can not access mod view for " << modId;
+ LOG(warn) << GetName() << "::Exec: can not access view for tracking layer " << n.fMaterialLayer << " @ "
+ << modId;
continue;
}
if (mod->hdXx.second) mod->hdXx.second->Fill(n.fMxy.X()[0], mod->hdXx.first * dx);
@@ -367,6 +395,9 @@ void CbmRecoQaTask::Exec(Option_t*)
n.fIsXySet = true;
}
+ // track QA filtering
+ if (!FilterTrack(nhit)) continue;
+
// Fit track with all hits ON for vx projection
int iprj(0);
for (auto zprj : fProjView) {
@@ -380,7 +411,15 @@ void CbmRecoQaTask::Exec(Option_t*)
fFitter.FitTrack(trkKf);
for (auto& n : trkKf.fNodes) {
if (n.fReference1 < 0) continue;
- //printf(" prj[%d] z[%+6.4f] x[%+6.4f] y[%+6.4f] t[%+6.4f]\n", n.fReference1, n.fZ, (double)n.fTrack.X()[0], (double)n.fTrack.Y()[0], (double)n.fTrack.Time()[0]);
+ if (fModView.find(n.fReference1 + nnodes) == fModView.end()) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fReference1 + nnodes << " not defined.";
+ continue;
+ }
+ ModuleView* mod = &fModView[n.fReference1 + nnodes][0];
+
+ 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++;
}
@@ -395,6 +434,23 @@ void CbmRecoQaTask::SetProjections(std::vector<double> vzpj)
}
}
+//_____________________________________________________________________
+bool CbmRecoQaTask::FilterTrack(int* ptr)
+{
+ bool ret = false;
+ switch (fSetupClass) {
+ case kMcbm22: {
+ int* nhit = ptr;
+ if (nhit[int(ECbmModuleId::kSts)] == 2 && nhit[int(ECbmModuleId::kTrd2d)] && nhit[int(ECbmModuleId::kTof)])
+ ret = true;
+ break;
+ }
+ case kMcbm24: break;
+ case kDefault: break;
+ }
+ return ret;
+}
+
//_____________________________________________________________________
void CbmRecoQaTask::Finish()
{
diff --git a/reco/qa/CbmRecoQaTask.h b/reco/qa/CbmRecoQaTask.h
index 156088158fc38da5371fe2c4017b688caf6cf121..44dca4e5d08b3c6959cb9e3de45e1d81aa2747cc 100644
--- a/reco/qa/CbmRecoQaTask.h
+++ b/reco/qa/CbmRecoQaTask.h
@@ -34,6 +34,11 @@ class CbmRecoQaTask : public FairTask {
/** \brief Executed task **/
virtual void Exec(Option_t* option);
+ /** \brief Filter tracks for further use (e.g. track projections)
+ * \param[in] ptr information on which to filter (e.g. nhits from detectors)
+ * \return true if track accepted
+ */
+ virtual bool FilterTrack(int* ptr);
virtual void Finish();
/** \brief Define the set of extra z positions where the track should be projected in the x-y plane */