diff --git a/reco/qa/CbmRecoQaTask.cxx b/reco/qa/CbmRecoQaTask.cxx
index e5617fe623cca27f4f483a913d0ff366b28474f5..3dab7e450ea7fc160c5e54883d98c706517784d4 100644
--- a/reco/qa/CbmRecoQaTask.cxx
+++ b/reco/qa/CbmRecoQaTask.cxx
@@ -41,7 +41,6 @@
using namespace std;
using namespace cbm::algo;
-std::bitset<kRecoQaNConfigs> CbmRecoQaTask::fuRecoConfig = {};
//_____________________________________________________________________
CbmRecoQaTask::CbmRecoQaTask() : FairTask("RecoQA", 0) {}
@@ -52,6 +51,12 @@ CbmRecoQaTask::Detector* CbmRecoQaTask::AddDetector(ECbmModuleId id)
/* Interface to the user. Check that the det defined in the outside world
* fulfills the quality criteria.
*/
+ if (GetDetector(id)) {
+ LOG(warn) << GetName() << "::AddDetector(" << ToString(id) << ")."
+ << " already registered. Using "
+ << "\"CbmRecoQaTask::GetDetector(ECbmModuleId).\"";
+ return GetDetector(id);
+ }
switch (id) {
case ECbmModuleId::kMvd:
case ECbmModuleId::kSts:
@@ -62,12 +67,6 @@ CbmRecoQaTask::Detector* CbmRecoQaTask::AddDetector(ECbmModuleId id)
case ECbmModuleId::kTof:
case ECbmModuleId::kFsd:
case ECbmModuleId::kPsd:
- if (GetDetector(id)) {
- LOG(warn) << GetName() << "::AddDetector : det " << ToString(id)
- << " already registered. Using "
- "\"CbmRecoQaTask::GetDetector(ECbmModuleId).\"";
- return GetDetector(id);
- }
LOG(debug) << GetName() << "::AddDetector(" << ToString(id) << ").";
fDetQa.emplace(id, id);
break;
@@ -129,13 +128,6 @@ InitStatus CbmRecoQaTask::Init()
else
fuRecoConfig.set(kRecoEvents);
- fTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrackMatch"));
- if (!fTrackMatches) {
- LOG(warn) << GetName() << "::Init: MC info for Global track not available !";
- }
-
- LOG(info) << GetName() << "::Init: Registered projections " << fProjs.size();
-
switch (fSetupClass) {
case kMcbm22:
LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2022\".";
@@ -151,19 +143,32 @@ InitStatus CbmRecoQaTask::Init()
break;
}
+ if (fuRecoConfig[kUseMC]) {
+ cbm_mc_manager = dynamic_cast<CbmMCDataManager*>(ioman->GetObject("MCDataManager"));
+ if (!cbm_mc_manager) {
+ LOG(warn) << GetName() << "::Init: MC data manager not available even though asked by user !";
+ fuRecoConfig.reset(kUseMC);
+ }
+ else {
+ fTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrackMatch"));
+ if (!fTrackMatches) {
+ LOG(warn) << GetName() << "::Init: MC info for Global track not available !";
+ }
+ }
+ }
+
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.";
+ 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!";
}
fHits[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.hit.name.data()));
@@ -172,33 +177,54 @@ InitStatus CbmRecoQaTask::Init()
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!";
}
- if (!det.Init(&fOutFolder)) continue;
- }
-
- // TDirectoryFile* trkDir = (TDirectoryFile*)fOutFolder.mkdir("track", "CA
- // QA"); for (auto& trkly : fTrkView) {
- // for (auto& view : trkly.second) view.Register(trkDir);
- // }
+ if (!det.Init(&fOutFolder, (bool) fuRecoConfig[kUseMC])) continue;
+ }
+
+ TDirectoryFile* trgDir = (fuRecoConfig[kRecoEvents] || GetNviews(eViewType::kTrkProj))
+ ? (TDirectoryFile*) fOutFolder.mkdir("TRG", "Target tomography with CA")
+ : nullptr;
+
+ // register track projection on z planes
+ if (fPrjPlanes.size() && trgDir) {
+ fViews.emplace("prj", View("TrkProj", "", {}));
+ fViews["prj"].fType = eViewType::kTrkProj;
+ int i(0);
+ for (auto plane : fPrjPlanes) { // add projection planes according to user request
+ if (i >= kNtrkProjections) {
+ LOG(warn) << GetName() << "::Init: Only " << kNtrkProjections << " are supported. Skipping the rest.";
+ fPrjPlanes.erase(fPrjPlanes.begin() + int(kNtrkProjections), fPrjPlanes.end());
+ break;
+ }
- TDirectoryFile* trgDir = (TDirectoryFile*) fOutFolder.mkdir("target", "Target tomography with CA");
- for (auto& prj : fProjs)
- prj.second.Register(trgDir);
+ fViews["prj"].AddProjection(CbmRecoQaTask::eProjectionType(int(CbmRecoQaTask::eProjectionType::kXYt0) + i),
+ plane.Z());
+ i++;
+ }
+ fViews["prj"].Init("");
+ fViews["prj"].Register(trgDir);
+ }
+ // register primary vertex QA
+ if (fuRecoConfig[kRecoEvents] && trgDir) {
+ fViews.emplace("vx", View("Vx", "", {}));
+ fViews["vx"].fType = eViewType::kPV;
+ fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVxy);
+ fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVxz);
+ fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVyz);
+ fViews["vx"].AddProjection(CbmRecoQaTask::eProjectionType::kPVmult);
+ fViews["vx"].Init("Prim");
+ fViews["vx"].Register(trgDir);
+ }
return kSUCCESS;
}
/// Hit classification on system and view
template<class Hit>
-bool CbmRecoQaTask::Detector::View::HasAddress(const CbmHit* h, double&, double&, double&, double&) const
+bool CbmRecoQaTask::View::HasAddress(const CbmHit* h, double&, double&, double&, double&) const
{
LOG(error) << "Unprocessed hit in view " << name;
cout << h->ToString();
@@ -207,8 +233,7 @@ 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, double& dx,
- double& dy) const
+bool CbmRecoQaTask::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),
@@ -241,8 +266,7 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmStsHit>(const CbmHit* h, doubl
}
// Rich specialization
template<>
-bool CbmRecoQaTask::Detector::View::HasAddress<CbmRichHit>(const CbmHit* h, double& x, double& y, double& dx,
- double& dy) const
+bool CbmRecoQaTask::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
@@ -265,8 +289,7 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmRichHit>(const CbmHit* h, doub
}
// TRD specialization
template<>
-bool CbmRecoQaTask::Detector::View::HasAddress<CbmTrdHit>(const CbmHit* h, double& x, double& y, double& dx,
- double& dy) const
+bool CbmRecoQaTask::View::HasAddress<CbmTrdHit>(const CbmHit* h, double& x, double& y, double& dx, double& dy) const
{
std::vector<int> sel;
@@ -303,8 +326,7 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmTrdHit>(const CbmHit* h, doubl
}
// ToF specialization
template<>
-bool CbmRecoQaTask::Detector::View::HasAddress<CbmTofHit>(const CbmHit* h, double& x, double& y, double& dx,
- double& dy) const
+bool CbmRecoQaTask::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)};
@@ -334,7 +356,7 @@ bool CbmRecoQaTask::Detector::View::HasAddress<CbmTofHit>(const CbmHit* h, doubl
return true;
}
template<class Hit>
-bool CbmRecoQaTask::Detector::View::Load(const CbmHit* h, const FairMCPoint* point, const CbmEvent* ev)
+bool CbmRecoQaTask::View::Load(const CbmHit* h, const FairMCPoint* point, const CbmEvent* ev)
{
double x(0), y(0), dx(0), dy(0);
if (!HasAddress<Hit>(h, x, y, dx, dy)) {
@@ -345,33 +367,34 @@ bool CbmRecoQaTask::Detector::View::Load(const CbmHit* h, const FairMCPoint* poi
// printf("View[%s] z(%d)[%f]\n", name.data(), (int)h->GetType(), h->GetZ());
int scale(0);
TH2* hh(nullptr);
- auto fillView = [&](eViewProjection proj, double xx, double yy, bool scaleY = true) {
+ auto fillView = [&](eProjectionType 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 (hh) hh->Fill(xx, scaleY ? yy * scale : yy);
}
};
if (ev) {
- fillView(eViewProjection::kChdT, x, h->GetTime() - ev->GetStartTime());
- fillView(eViewProjection::kXYh, x, y, false);
+ fillView(eProjectionType::kChdT, x, h->GetTime() - ev->GetStartTime());
+ //fillView(eProjectionType::kXYh, x, y, false);
}
- fillView(eViewProjection::kXYs, x, y, false);
+ fillView(eProjectionType::kXYh, 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);
- }
+ fillView(eProjectionType::kXYhMC, point->GetX(), point->GetY(), false);
+ fillView(eProjectionType::kResidualX, point->GetX(), x - point->GetX());
+ fillView(eProjectionType::kResidualY, point->GetY(), y - point->GetY());
+ fillView(eProjectionType::kResidualTX, point->GetX(), h->GetTime() - event_time - point->GetTime());
+ fillView(eProjectionType::kResidualTY, point->GetY(), h->GetTime() - event_time - point->GetTime());
+ fillView(eProjectionType::kPullX, point->GetX(), (x - point->GetX()) / dx);
+ fillView(eProjectionType::kPullY, point->GetY(), (y - point->GetY()) / dy);
+ }
+ fMult++; // register hit/point multiplicity
return true;
}
@@ -389,13 +412,14 @@ bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::TrajectoryNode*
if (!hh) continue;
switch (projection.first) {
- case eViewProjection::kXYa: hh->Fill(n->fMxy.X(), n->fMxy.Y()); break;
+ case eProjectionType::kXYa: hh->Fill(n->fMxy.X(), n->fMxy.Y()); break;
case eViewProjection::kXYp: hh->Fill(t.X(), t.Y()); break;
- case eViewProjection::kXdX: hh->Fill(n->fMxy.X(), scale * dx); break;
- case eViewProjection::kYdY: hh->Fill(n->fMxy.Y(), scale * dy); break;
- case eViewProjection::kWdT: hh->Fill(n->fMxy.X(), scale * dt); break;
- case eViewProjection::kXpX: hh->Fill(n->fMxy.X(), pullx); break;
- case eViewProjection::kYpY: hh->Fill(n->fMxy.Y(), pully); break;
+ case eProjectionType::kXdX: hh->Fill(n->fMxy.X(), scale * dx); break;
+ case eProjectionType::kYdY: hh->Fill(n->fMxy.Y(), scale * dy); break;
+ case eProjectionType::kWdT: hh->Fill(n->fMxy.X(), scale * dt); break;
+ case eProjectionType::kXpX: hh->Fill(n->fMxy.X(), pullx); break;
+ case eProjectionType::kYpY: hh->Fill(n->fMxy.Y(), pully); break;
+
default: break;
}
@@ -403,8 +427,8 @@ bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::TrajectoryNode*
double dxMC = point->GetX() - t.X(), dyMC = point->GetY() - t.Y();
switch (projection.first) {
- case eViewProjection::kXdXMC: hh->Fill(point->GetX(), scale * dxMC); break;
- case eViewProjection::kYdYMC: hh->Fill(point->GetY(), scale * dyMC); break;
+ case eProjectionType::kXdXMC: hh->Fill(point->GetX(), scale * dxMC); break;
+ case eProjectionType::kYdYMC: hh->Fill(point->GetY(), scale * dyMC); break;
default: break;
}
}
@@ -412,19 +436,66 @@ bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::TrajectoryNode*
return true;
}
+bool CbmRecoQaTask::View::Load(TVector3* p)
+{
+ for (auto& projection : fProjection) {
+ int scale = get<0>(projection.second);
+ TH2* hh = get<2>(projection.second);
+ if (!hh) continue;
+ switch (projection.first) {
+ case eProjectionType::kDmult:
+ if (int(p->Z()) == -124) hh->Fill(p->X(), p->Y());
+ break;
+ case eProjectionType::kXYt0:
+ if (int(p->Z()) == 0) hh->Fill(scale * p->X(), scale * p->Y());
+ break;
+ // fall through
+ case eProjectionType::kXYt1:
+ if (int(p->Z()) == 1) hh->Fill(scale * p->X(), scale * p->Y());
+ break;
+ // fall through
+ case eProjectionType::kXYt2:
+ if (int(p->Z()) == 2) hh->Fill(scale * p->X(), scale * p->Y());
+ break;
+ // fall through
+ case eProjectionType::kXYt3:
+ if (int(p->Z()) == 3) hh->Fill(scale * p->X(), scale * p->Y());
+ break;
+ // fall through
+ case eProjectionType::kXYt4:
+ if (int(p->Z()) == 4) hh->Fill(scale * p->X(), scale * p->Y());
+ break;
+ // fall through
+ case eProjectionType::kXYt5:
+ if (int(p->Z()) == 5) hh->Fill(scale * p->X(), scale * p->Y());
+ break;
+ // fall through
+ case eProjectionType::kPVxy: hh->Fill(scale * p->X(), scale * p->Y()); break;
+ case eProjectionType::kPVxz: hh->Fill(p->Z(), scale * p->X()); break;
+ case eProjectionType::kPVyz: hh->Fill(p->Z(), scale * p->Y()); break;
+ case eProjectionType::kPVmult:
+ if (int(p->Z()) == -123) hh->Fill(scale * p->X(), scale * p->Y());
+ break;
+ default: break;
+ }
+ }
+ return true;
+}
+
//_____________________________________________________________________
void CbmRecoQaTask::Exec(Option_t*)
{
LOG(info) << GetName() << "::Exec : Evs[" << (fEvents ? fEvents->GetEntriesFast() : 0) << "] Trks["
- << (fGTracks ? fGTracks->GetEntriesFast() : 0) << "]";
+ << (fGTracks ? fGTracks->GetEntriesFast() : 0) << "] Hits["
+ << (fHits[ECbmModuleId::kSts] ? fHits[ECbmModuleId::kSts]->GetEntriesFast() : 0) << " "
+ << (fHits[ECbmModuleId::kTrd] ? fHits[ECbmModuleId::kTrd]->GetEntriesFast() : 0) << " "
+ << (fHits[ECbmModuleId::kTof] ? fHits[ECbmModuleId::kTof]->GetEntriesFast() : 0) << "]";
+ // fixed momentum no magentic field for mCBM
fFitter.SetDefaultMomentumForMs(0.5);
- fFitter.FixMomentumForMs(true); // consider only mcbm case for a moment
-
-
+ fFitter.FixMomentumForMs(true);
int iev = 0, itrack = 0, nnodes = 0;
-
auto processHits = [&](CbmEvent* ev) {
for (auto& detp : fDetQa) {
auto& det = detp.second;
@@ -439,7 +510,11 @@ void CbmRecoQaTask::Exec(Option_t*)
const FairMCPoint* mcpoint = nullptr;
const CbmHit* hit = dynamic_cast<CbmHit*>(fHits[det.id]->At(jh));
- if (!hit) continue;
+ if (!hit) {
+ LOG(warning) << GetName() << "::Exec() : Hit " << jh << " for " << ToString(det.id)
+ << " not available. Skip.";
+ continue;
+ }
if (fHitMatch[det.id] && fPoints[det.id]) {
if (det.id == ECbmModuleId::kRich) {
@@ -451,13 +526,10 @@ void CbmRecoQaTask::Exec(Option_t*)
}
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();
@@ -491,25 +563,58 @@ void CbmRecoQaTask::Exec(Option_t*)
if (ret) break;
}
}
+ TVector3 mult;
+ for (auto& v : det.fViews) {
+ mult.SetXYZ(nh, double(v.fMult) / nh, -124);
+ v.Load(&mult);
+ v.fMult = 0;
+ }
}
};
auto processTracks = [&](CbmEvent* ev) {
const int ntrk = (ev) ? ev->GetNofData(ECbmDataType::kGlobalTrack) : fGTracks->GetEntriesFast();
+ // read in the vertex and the list of tracks used for its definition
+ TVector3 pvx, evx;
+ int ntrkDet[3] = {0};
+
+ CbmVertex* vx(nullptr);
+ int nTrkVx(0);
+ if (ev && fViews.find("vx") != fViews.end()) {
+ vx = ev->GetVertex();
+ nTrkVx = vx->GetNTracks();
+ auto& v = fViews["vx"];
+ if (nTrkVx >= 2) {
+ vx->Position(pvx);
+ for (int i(0); i < 3; i++) {
+ evx[i] = vx->GetCovariance(i, i);
+ if (evx[i] > 0.) evx[i] = sqrt(evx[i]);
+ }
+ v.Load(&pvx);
+ }
+ pvx.SetXYZ(ntrk, nTrkVx, -123);
+ v.Load(&pvx);
+ }
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;
+ int trkIdx = ev ? ev->GetIndex(ECbmDataType::kGlobalTrack, itrk) : itrk;
+ auto track = dynamic_cast<const CbmGlobalTrack*>((*fGTracks)[trkIdx]);
+ // track QA filtering
+ if (!FilterTrack(track)) continue;
+
+ if (nTrkVx >= 2) {
+ if (vx->FindTrackByIndex(trkIdx)) {
+ if (track->GetStsTrackIndex() >= 0) ntrkDet[0]++;
+ if (track->GetTrdTrackIndex() >= 0) ntrkDet[1]++;
+ if (track->GetTofTrackIndex() >= 0) ntrkDet[2]++;
+ }
+ }
- auto track = dynamic_cast<const CbmGlobalTrack*>(trackObj);
CbmKfTrackFitter::Trajectory 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) {
@@ -522,7 +627,6 @@ void CbmRecoQaTask::Exec(Option_t*)
LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
continue;
}
-
n.fIsTimeSet = n.fIsXySet = false;
fFitter.FitTrajectory(trkKf);
@@ -552,28 +656,31 @@ void CbmRecoQaTask::Exec(Option_t*)
n.fIsTimeSet = n.fIsXySet = true;
}
- // Fit track with all hits ON for vx projection
- int iprj(0);
- for (auto& zprj : fProjs) {
+ // Fit track with all hits ON for track projections
+ if (!nnodes) nnodes = (int) trkKf.fNodes.size();
+ if (fViews.find("prj") != fViews.end()) {
+ auto v = fViews["prj"];
+ if (v.fType != eViewType::kTrkProj)
+ LOG(error) << GetName() << "::Exec: view for track projection with wrong type. Skipping.";
+ else {
+ int iprj(0);
+ for (auto plane : fPrjPlanes) {
CbmKfTrackFitter::TrajectoryNode n = CbmKfTrackFitter::TrajectoryNode();
- n.fZ = zprj.first;
- n.fIsXySet = true;
- n.fReference1 = iprj++;
- trkKf.fNodes.push_back(n);
- }
- trkKf.OrderNodesInZ();
+ n.fZ = plane.Z();
+ n.fIsXySet = true;
+ n.fReference1 = iprj++;
+ trkKf.fNodes.push_back(n);
+ }
+ trkKf.OrderNodesInZ();
fFitter.FitTrajectory(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;
+
+ TVector3 xyz;
+ for (auto& n : trkKf.fNodes) {
+ if (n.fReference1 < 0) continue;
+ xyz.SetXYZ(n.fTrack.X(), n.fTrack.Y(), n.fReference1); // to communicate to the Load function
+ v.Load(&xyz);
+ }
}
- // 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]);
}
@@ -589,9 +696,10 @@ void CbmRecoQaTask::Exec(Option_t*)
if (fGTracks) processTracks(ev);
iev++;
}
- LOG(info) << GetName() << "::Exec : Evs(%)[" << 1.e2 * iev / fEvents->GetEntriesFast() << "] Trks(%)["
- << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
- return;
+ LOG(info) << GetName() << "::Exec : Evs(%)["
+ << (fEvents->GetEntriesFast() ? 1.e2 * iev / fEvents->GetEntriesFast() : 0) << "] Trks(%)["
+ << (fGTracks && fGTracks->GetEntriesFast() ? 1.e2 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
+ return; // END EbyE QA
}
else {
@@ -604,21 +712,23 @@ void CbmRecoQaTask::Exec(Option_t*)
processTracks(nullptr);
}
- LOG(info) << GetName() << "::Exec : Trks(%)[" << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
+ LOG(info) << GetName() << "::Exec : Trks(%)["
+ << (fGTracks && fGTracks->GetEntriesFast() ? 1.e2 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
return;
}
//_____________________________________________________________________
-void CbmRecoQaTask::SetProjections(std::vector<double> vzpj)
+int CbmRecoQaTask::GetNviews(eViewType type) const
{
- fProjs.clear();
- for (auto zprj : vzpj) {
- // 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), "", {}));
+ int n(0);
+ for (auto v : fViews) {
+ if (v.second.fType != type) continue;
+ n++;
}
+ return n;
}
+
//_____________________________________________________________________
void CbmRecoQaTask::Finish()
{
@@ -629,6 +739,9 @@ void CbmRecoQaTask::Finish()
//_____________________________________________________________________
bool CbmRecoQaTask::FilterEvent(const CbmEvent* ptr)
{
+ // sanity checks
+ if (ptr->GetStartTime() < 0) return false; // this is found in MC
+
for (auto evCut : fFilterEv) {
if (!evCut.Accept(ptr, this)) return false;
}
@@ -703,8 +816,11 @@ void CbmRecoQaTask::InitMcbm22()
LOG(fatal) << GetName() << "::InitMcbm22() : Missing setup definition.";
return;
}
+ // fixed momentum no magentic field for mCBM
+ fFitter.SetDefaultMomentumForMs(0.5);
+ fFitter.FixMomentumForMs(true);
- Detector::View* v = nullptr;
+ View* v = nullptr;
std::vector<TString> path;
TGeoNode* topNode = gGeoManager->GetTopNode();
@@ -745,8 +861,14 @@ void CbmRecoQaTask::InitMcbm22()
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");
+ v->SetProjection(eProjectionType::kXdX, (station == 1) ? 1 : 0.75, "mm");
+ v->SetProjection(eProjectionType::kYdY, (station == 1) ? 4 : 3, "mm");
+ if (fuRecoConfig[kUseMC]) {
+ v->SetProjection(eProjectionType::kXdXMC, 1, "mm");
+ v->SetProjection(eProjectionType::kYdYMC, 3, "mm");
+ v->SetProjection(eProjectionType::kResidualX, 200, "um");
+ v->SetProjection(eProjectionType::kResidualY, 500, "um");
+ }
}
else {
std::cout << "No match found in string: " << str << std::endl;
@@ -763,9 +885,15 @@ void CbmRecoQaTask::InitMcbm22()
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");
+ v->SetProjection(eProjectionType::kChdT, 400, "ns");
+ v->SetProjection(eProjectionType::kXdX, 10, "mm");
+ v->SetProjection(eProjectionType::kYdY, 20, "mm");
+ if (fuRecoConfig[kUseMC]) {
+ v->SetProjection(eProjectionType::kXdXMC, 10, "mm");
+ v->SetProjection(eProjectionType::kYdYMC, 20, "mm");
+ v->SetProjection(eProjectionType::kResidualX, 1.5, "mm");
+ v->SetProjection(eProjectionType::kResidualY, 5.0, "mm");
+ }
}
// Trd1DxDy
trd = AddDetector(ECbmModuleId::kTrd);
@@ -773,16 +901,26 @@ void CbmRecoQaTask::InitMcbm22()
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");
+ v->SetProjection(eProjectionType::kChdT, 400, "ns");
+ v->SetProjection(eProjectionType::kXdX, 1.5, "cm");
+ if (fuRecoConfig[kUseMC]) {
+ v->SetProjection(eProjectionType::kXdXMC, 1.5, "cm");
+ v->SetProjection(eProjectionType::kResidualX, 1.5, "mm");
+ }
}
if (str.Contains("layer03")) {
v = trd->AddView("1Dy", str.Data(), {37});
- v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ v->SetProjection(eProjectionType::kChdT, 400, "ns");
+ v->SetProjection(eProjectionType::kYdY, 1.5, "cm");
+ if (fuRecoConfig[kUseMC]) {
+ v->SetProjection(eProjectionType::kYdYMC, 1.5, "cm");
+ v->SetProjection(eProjectionType::kResidualY, 1.5, "mm");
+ }
}
}
}
+
if (setup->IsActive(ECbmModuleId::kTof)) {
processDetector("tof", "counter");
@@ -800,7 +938,7 @@ void CbmRecoQaTask::InitMcbm22()
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");
+ v->SetProjection(eProjectionType::kChdT, 15, "ns");
}
else {
v = tof->AddView(name, str.Data(), {smid, type, rpc});
@@ -813,22 +951,13 @@ void CbmRecoQaTask::InitMcbm22()
}
// ===============================================================================
// TRG - upstream projections
- // int iprjLy = 0;
- // for (auto prj : fProjs) {
- // LOG(debug) << GetName() << "::Init: Project trks @ " << prj.first;
- // 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),
- // 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}
+ float angle = 25., L[] = {14.3, 0, -20, -38 /*-1*/, -50.5 /*-4.1*/};
+ //R[]= {2.5, 3, 0.5, 0.6, 0.6}, dx[5], dz[5];
+ for (int i(0); i < 5; i++) {
+ fPrjPlanes.emplace_back(L[i] * TMath::Sin(angle * TMath::DegToRad()), 0.,
+ L[i] * TMath::Cos(angle * TMath::DegToRad()));
+ }
// (cm); T_{trk} - T_{ev} (ns)", prj.first),
- // 700, -25, 10, 350, -15.5, 5.5);
- // iprjLy++;
- // }
}
//_____________________________________________________________________
@@ -839,8 +968,12 @@ void CbmRecoQaTask::InitMcbm24()
LOG(fatal) << GetName() << "::InitMcbm24() : Missing setup definition.";
return;
}
+ // fixed momentum no magentic field for mCBM
+ fFitter.SetDefaultMomentumForMs(0.5);
+ fFitter.FixMomentumForMs(true);
+
TString dtag;
- Detector::View* v(nullptr);
+ View* v(nullptr);
if (setup->IsActive(ECbmModuleId::kSts)) {
// setup->GetGeoTag(ECbmModuleId::kSts, dtag);
// LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
@@ -852,76 +985,76 @@ void CbmRecoQaTask::InitMcbm24()
"HalfLadder13u_1/HalfLadder13u_Module03_1/Sensor03_1",
dtag.Data()),
{0, 0, 0});
- v->SetProjection(eViewProjection::kYdY, 5, "mm");
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
// 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->SetProjection(eViewProjection::kYdY, 5, "mm");
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
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");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
// 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->SetProjection(eViewProjection::kYdY, 5, "mm");
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
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");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
// 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->SetProjection(eViewProjection::kYdY, 5, "mm");
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
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");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
// 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->SetProjection(eViewProjection::kYdY, 5, "mm");
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
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");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
// 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->SetProjection(eViewProjection::kYdY, 5, "mm");
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
v = sts->AddView("U2L2M1",
Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/"
"HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1",
@@ -932,8 +1065,8 @@ void CbmRecoQaTask::InitMcbm24()
"HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1",
dtag.Data()),
{2, 2, 2});
- v->SetProjection(eViewProjection::kYdY, 5, "mm");
- v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eProjectionType::kYdY, 2, "mm");
+ v->SetProjection(eProjectionType::kXdX, 500, "um");
}
if (setup->IsActive(ECbmModuleId::kTrd)) {
// setup->GetGeoTag(ECbmModuleId::kTrd, dtag);
@@ -942,16 +1075,18 @@ void CbmRecoQaTask::InitMcbm24()
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, 300, "ns");
+ v->SetProjection(eProjectionType::kChdT, 300, "ns");
// Trd1Dx
trd = AddDetector(ECbmModuleId::kTrd);
v = trd->AddView("1Dx", Form("/cave_1/trd_%s_0/layer02_10202/module5_101002001", dtag.Data()), {21});
- v->SetProjection(eViewProjection::kChdT, 300, "ns");
- v->SetProjection(eViewProjection::kXdX, 7, "cm");
+ v->SetProjection(eProjectionType::kChdT, 300, "ns");
+ v->SetProjection(eProjectionType::kXdX, 1.5, "cm");
+ v->SetProjection(eProjectionType::kYdY, 3.0, "cm");
// Trd1Dy
v = trd->AddView("1Dy", Form("/cave_1/trd_%s_0/layer03_11303/module5_101103001", dtag.Data()), {37});
- v->SetProjection(eViewProjection::kChdT, 300, "ns");
- v->SetProjection(eViewProjection::kXdX, 7, "cm");
+ v->SetProjection(eProjectionType::kChdT, 300, "ns");
+ v->SetProjection(eProjectionType::kXdX, 1.5, "cm");
+ v->SetProjection(eProjectionType::kYdY, 3.0, "cm");
}
if (setup->IsActive(ECbmModuleId::kTof)) {
// setup->GetGeoTag(ECbmModuleId::kTof, dtag);
@@ -969,8 +1104,8 @@ void CbmRecoQaTask::InitMcbm24()
"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");
+ v->SetProjection(eProjectionType::kXdX, 3, "cm");
+ v->SetProjection(eProjectionType::kYdY, 1, "cm");
}
}
// add type 6 (Buch)
@@ -1019,6 +1154,14 @@ void CbmRecoQaTask::InitMcbm24()
// Aerogel 1
rich->AddView("Aerogel", Form("/cave_1/rich_%s_0/box_1/Gas_1", dtag.Data()), {4});
}
+ // ===============================================================================
+ // TRG - upstream projections
+ float angle = 25., L[] = {14.3, 0, -20, -38, -50.5};
+ //R[]= {2.5, 3, 0.5, 0.6, 0.6}, dx[5], dz[5];
+ for (int i(0); i < 5; i++) {
+ fPrjPlanes.emplace_back(L[i] * TMath::Sin(angle * TMath::DegToRad()), 0.,
+ L[i] * TMath::Cos(angle * TMath::DegToRad()));
+ }
}
//_____________________________________________________________________
void CbmRecoQaTask::InitDefault()
@@ -1030,7 +1173,7 @@ void CbmRecoQaTask::InitDefault()
return;
}
- Detector::View* v = nullptr;
+ View* v = nullptr;
std::vector<TString> path;
TGeoNode* topNode = gGeoManager->GetTopNode();
if (!topNode) {
@@ -1068,7 +1211,7 @@ void CbmRecoQaTask::InitDefault()
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);
+ v->SetSetup(eSetup::kDefault);
}
else {
std::cout << "No match found in string: " << str << std::endl;
@@ -1097,7 +1240,7 @@ void CbmRecoQaTask::InitDefault()
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);
+ v->SetSetup(eSetup::kDefault);
}
else {
std::cout << "No match found in string: " << str << std::endl;
@@ -1110,7 +1253,7 @@ void CbmRecoQaTask::InitDefault()
Detector* rich = AddDetector(ECbmModuleId::kRich);
for (const auto& str : path) {
v = rich->AddView("rich", str.Data(), {-1, -1, -1});
- v->SetMode(eSetup::kDefault);
+ v->SetSetup(eSetup::kDefault);
}
}
@@ -1130,7 +1273,7 @@ void CbmRecoQaTask::InitDefault()
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);
+ v->SetSetup(eSetup::kDefault);
}
else {
std::cout << "No match found in string: " << str << std::endl;
@@ -1138,6 +1281,7 @@ void CbmRecoQaTask::InitDefault()
}
}
}
+
//========= DETECTOR ======================
CbmRecoQaTask::Detector::Detector(ECbmModuleId did)
{
@@ -1188,37 +1332,37 @@ CbmRecoQaTask::Detector::Detector(ECbmModuleId did)
//+++++++++++++++++++++ Detector ++++++++++++++++++++++++++++++++
//_________________________________________________________________
-CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::AddView(const char* n, const char* p, std::vector<int> set)
+CbmRecoQaTask::View* CbmRecoQaTask::Detector::AddView(const char* n, const char* p, std::vector<int> set)
{
View* v(nullptr);
fViews.emplace_back(n, p, set);
v = &fViews.back();
- v->AddProjection(CbmRecoQaTask::eViewProjection::kXYs);
-
- 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);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kChdT);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kXYh);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kXYhMC);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kDmult);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kDmultMC);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kPullX);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kPullY);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualX);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualY);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualTX);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kResidualTY);
+
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kXYa);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kXYp);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kXdX);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kXdXMC);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kYdY);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kYdYMC);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kWdT);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kXpX);
+ v->AddProjection(CbmRecoQaTask::eProjectionType::kYpY);
return v;
}
-CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::GetView(const char* n)
+CbmRecoQaTask::View* CbmRecoQaTask::Detector::GetView(const char* n)
{
for (auto& view : fViews)
if (view.name.compare(n) == 0) return &view;
@@ -1226,7 +1370,7 @@ CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::GetView(const char* n)
return nullptr;
}
-CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::FindView(double x, double y, double z)
+CbmRecoQaTask::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,
@@ -1241,7 +1385,7 @@ CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::FindView(double x, doubl
return nullptr;
}
-bool CbmRecoQaTask::Detector::Init(TDirectoryFile* fOut)
+bool CbmRecoQaTask::Detector::Init(TDirectoryFile* fOut, bool mc)
{
/* Query the geometry and trigger detector views init
*/
@@ -1249,13 +1393,14 @@ bool CbmRecoQaTask::Detector::Init(TDirectoryFile* fOut)
LOG(fatal) << "CbmRecoQaTask::Detector::Init() " << ToString(id) << " missing geometry.";
return false;
}
- LOG(debug) << "CbmRecoQaTask::Detector::Init() : " << ToString(id);
+ //LOG(debug) << "CbmRecoQaTask::Detector::Init() : " << ToString(id);
+ LOG(info) << "CbmRecoQaTask::Detector::Init() : " << ToString(id) << " MC = " << mc;
TDirectoryFile* modDir =
(TDirectoryFile*) fOut->mkdir(ToString(id).data(), Form("Reco QA for %s", ToString(id).data()));
bool ret(true);
for (auto& view : fViews) {
- bool vret = view.Init(ToString(id).data());
+ bool vret = view.Init(ToString(id).data(), mc);
view.Register(modDir);
ret = ret && vret;
}
@@ -1271,9 +1416,9 @@ void CbmRecoQaTask::Detector::Print() const
cout << s.str();
}
-//+++++++++++++++++++++ Detector::View ++++++++++++++++++++++++++++++++
+//+++++++++++++++++++++ View ++++++++++++++++++++++++++++++++
//_______________________________________________________________________
-bool CbmRecoQaTask::Detector::View::AddProjection(eViewProjection prj, float range, const char* unit)
+bool CbmRecoQaTask::View::AddProjection(eProjectionType prj, float range, const char* unit)
{
if (fProjection.find(prj) != fProjection.end()) {
LOG(warn) << "Projection " << ToString(prj) << " already registered";
@@ -1303,12 +1448,12 @@ bool CbmRecoQaTask::Detector::View::AddProjection(eViewProjection prj, float ran
return true;
}
-bool CbmRecoQaTask::Detector::View::SetProjection(eViewProjection prj, float range, const char* unit)
+bool CbmRecoQaTask::View::SetProjection(eProjectionType prj, float range, const char* unit)
{
if (fProjection.find(prj) == fProjection.end()) {
LOG(warn) << "Projection " << ToString(prj)
<< " not initialized. Calling "
- "\"CbmRecoQaTask::Detector::View::AddProjection()\"";
+ "\"CbmRecoQaTask::View::AddProjection()\"";
return AddProjection(prj, range, unit);
}
int scale(1);
@@ -1336,21 +1481,30 @@ bool CbmRecoQaTask::Detector::View::SetProjection(eViewProjection prj, float ran
return true;
}
-bool CbmRecoQaTask::Detector::View::Init(const char* dname)
+bool CbmRecoQaTask::View::Init(const char* dname, bool mc)
{
- if (!gGeoManager->cd(path.data())) return false;
- TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
- const double* tr(m->GetTranslation());
- TGeoVolume* v = gGeoManager->GetCurrentVolume();
- TGeoShape* bb = v->GetShape();
- double w_lo, w_hi;
- for (int i(0); i < 3; i++) { // loop over the spatial dimensions
- size[i] = bb->GetAxisRange(i + 1, w_lo, w_hi);
- pos[i] = 0.5 * (w_lo + w_hi) + tr[i];
+ /** All projections type according to enum eProjectionType are defined here.
+ * Additionally the linking with the geometry, in the case of detector view type is also done to check the alignment.
+ * If the mc flag is set, histograms relevant for MC truth are also build along side the reconstruction ones.
+ */
+ if (path.size() && fType == eViewType::kDetUnit) { // applies only for detection units
+ if (!gGeoManager->cd(path.data())) return false;
+ TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
+ const double* tr(m->GetTranslation());
+ TGeoVolume* v = gGeoManager->GetCurrentVolume();
+ TGeoShape* bb = v->GetShape();
+ double w_lo, w_hi;
+ for (int i(0); i < 3; i++) { // loop over the spatial dimensions
+ size[i] = bb->GetAxisRange(i + 1, w_lo, w_hi);
+ pos[i] = 0.5 * (w_lo + w_hi) + tr[i];
+ }
+
+ LOG(info) << "CbmRecoQaTask::Detector(" << dname << ")::View(" << name << ")::Init() : size [" << size[0] << " x "
+ << size[1] << " x " << size[2] << "]. mc[" << mc << "].";
}
+ else
+ LOG(info) << "CbmRecoQaTask::Detector(" << dname << ")::View(" << name << ")::Init() : mc[" << mc << "].";
- LOG(debug) << "CbmRecoQaTask::Detector(" << dname << ")::View(" << name << ")::Init() : size [" << size[0] << " x "
- << size[1] << " x " << size[2] << "].";
// TODO short-hand for XY display resolution
int dscale = 10;
@@ -1362,40 +1516,47 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
int scale = get<0>(projection.second);
float yrange = get<1>(projection.second);
char xy_id = 0;
+ char mc_id[2] = {' ', ' '};
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';
+ case eProjectionType::kXdXMC:
+ if (!mc) break;
+ xy_id = 'M';
+ mc_id[0] = 'M';
+ mc_id[1] = 'C';
// fall through
- case eViewProjection::kXdX:
+ case eProjectionType::kXdX:
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()),
+ Form("X resolution %s %s [%s]; X^{%s}_{%s-%s} (cm); #Delta X (%s)", (xy_id == 'M' ? "MC" : ""),
+ name.data(), dname, mc_id, dname, name.data(), unit.data()),
nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
break;
- case eViewProjection::kYdYMC: xy_id = 'M';
+ case eProjectionType::kYdYMC:
+ if (!mc) break;
+ xy_id = 'M';
+ mc_id[0] = 'M';
+ mc_id[1] = 'C';
// fall through
- case eViewProjection::kYdY:
+ case eProjectionType::kYdY:
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()),
+ Form("Y resolution %s %s [%s]; Y^{%s}_{%s-%s} (cm); #Delta Y (%s)", (xy_id == 'M' ? "MC" : ""),
+ name.data(), dname, mc_id, dname, name.data(), unit.data()),
nbinsx, ylo, yhi, nbinsy, -yrange, yrange);
break;
-
- case eViewProjection::kXpX:
+ case eProjectionType::kXpX:
nbinsx = 10 * ceil(size[0]); // mm binning
if (yrange < 0) {
LOG(debug) << "ProjectionP[" << name << "] using default range.";
@@ -1404,11 +1565,11 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
nbinsy = 200; //* ceil(yrange);
get<2>(projection.second) =
new TH2D(Form("hpx_%s%s", dname, name.data()),
- Form("X pulls %s [%s]; X_{%s-%s} (cm); pull(X)", name.data(), dname, dname, name.data()), nbinsx,
- xlo, xhi, nbinsy, -yrange, yrange);
+ Form("X pulls %s [%s]; X^{%s}_{%s-%s} (cm); pull(X)", name.data(), dname, mc_id, dname, name.data()),
+ nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
break;
- case eViewProjection::kYpY:
+ case eProjectionType::kYpY:
nbinsx = 10 * ceil(size[1]); // mm binning
if (yrange < 0) {
LOG(debug) << "ProjectionP[" << name << "] using default range.";
@@ -1417,11 +1578,11 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
nbinsy = 200; //* ceil(yrange);
get<2>(projection.second) =
new TH2D(Form("hpy_%s%s", dname, name.data()),
- Form("Y pulls %s [%s]; Y_{%s-%s} (cm); pull(Y)", name.data(), dname, dname, name.data()), nbinsx,
- ylo, yhi, nbinsy, -yrange, yrange);
+ Form("Y pulls %s [%s]; Y^{%s}_{%s-%s} (cm); pull(Y)", name.data(), dname, mc_id, dname, name.data()),
+ nbinsx, ylo, yhi, nbinsy, -yrange, yrange);
break;
- case eViewProjection::kWdT:
+ case eProjectionType::kWdT:
nbinsx = 10 * ceil(size[0]); // mm binning
unit = makeTrange(scale, yrange);
nbinsy = 2 * ceil(yrange);
@@ -1431,23 +1592,34 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
break;
- case eViewProjection::kChdT:
+ case eProjectionType::kChdT:
nbinsx = dscale * ceil(size[0]); // mm binning
unit = makeTrange(scale, yrange);
nbinsy = 10 * ceil(yrange);
get<2>(projection.second) = new TH2D(Form("hct_%s%s", dname, name.data()),
- Form("Hit_Event %s [%s]; X_{%s-%s} (cm); #Delta time_{EV} (%s)",
+ Form("Hit Event %s [%s]; X_{%s-%s} (cm); #Delta time_{EV} (%s)",
name.data(), dname, dname, name.data(), unit.data()),
nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
break;
+ case eProjectionType::kDmultMC:
+ if (!mc) break;
+ xy_id = 'm';
+ // fall through
+ case eProjectionType::kDmult:
+ nbinsx = 100;
+ xlo = -0.5;
+ xhi = xlo + nbinsx;
+ get<2>(projection.second) = new TH2D(
+ Form("hdm_%s%s", dname, name.data()),
+ Form("%s multiplicity [EbyE] %s [%s]; N_{%s}^{%s}; N_{%s-%s}^{%s} (%%)", (xy_id ? "point" : "hit"),
+ name.data(), dname, dname, (xy_id ? "point" : "hit"), dname, name.data(), (xy_id ? "point" : "hit")),
+ nbinsx, xlo, xhi, 200, 0, 1);
+ break;
- case eViewProjection::kXYa:
+ case eProjectionType::kXYa:
xy_id = 'a';
// fall through
- case eViewProjection::kXYs:
- if (!xy_id) xy_id = 's';
- // fall through
- case eViewProjection::kXYh:
+ case eProjectionType::kXYh:
if (!xy_id) xy_id = 'h';
nbinsx = dscale * ceil(size[0]); // mm binning
nbinsy = dscale * ceil(size[1]);
@@ -1458,8 +1630,11 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
nbinsx, xlo, xhi, nbinsy, ylo, yhi);
break;
- case eViewProjection::kXYhMC:
+ case eProjectionType::kXYhMC:
+ if (!mc) break;
if (!xy_id) xy_id = 'h';
+ mc_id[0] = 'M';
+ mc_id[1] = 'C';
nbinsx = dscale * ceil(size[0]); // mm binning
nbinsy = dscale * ceil(size[1]);
get<2>(projection.second) =
@@ -1469,74 +1644,166 @@ bool CbmRecoQaTask::Detector::View::Init(const char* dname)
nbinsx, xlo, xhi, nbinsy, ylo, yhi);
break;
- case eViewProjection::kResidualX:
+ case eProjectionType::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);
+ case eProjectionType::kResidualY:
+ if (!mc) break;
+ mc_id[0] = 'M';
+ mc_id[1] = 'C';
+ {
+ bool isResidualX = (projection.first == eProjectionType::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;
+ }
- 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 eProjectionType::kResidualTX:
+ // fall through
+ case eProjectionType::kResidualTY:
+ if (!mc) break;
+ {
+ bool isResidualTX = (projection.first == eProjectionType::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::kResidualTX:
+ case eProjectionType::kPullX:
// fall through
- case eViewProjection::kResidualTY: {
- bool isResidualTX = (projection.first == eViewProjection::kResidualTX);
+ case eProjectionType::kPullY:
+ if (!mc) break;
+ {
+ bool isPullX = (projection.first == eProjectionType::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)", (isPullX ? "X" : "Y"), name.data(), dname,
+ (isPullX ? "X" : "Y"), dname, name.data(), (isPullX ? "X" : "Y")),
+ nbinsx, (isPullX ? xlo : ylo), (isPullX ? xhi : yhi), nbinsy, -yrange, yrange);
+ break;
+ }
- nbinsx = dscale * ceil(size[0]); // mm binning
- nbinsy = dscale * ceil(size[1]);
+ case eProjectionType::kXYp:
+ nbinsx = dscale * ceil(size[0] * 2.); // mm binning
+ nbinsy = dscale * ceil(size[1] * 2.);
+ xlo = pos[0] - size[0];
+ xhi = pos[0] + size[0];
+ ylo = pos[1] - size[1];
+ yhi = pos[1] + 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);
+ get<2>(projection.second) = new TH2D(Form("hxyp_%s%s", dname, name.data()),
+ Form("Trk_{proj} %s [%s]; X_{%s-%s} (cm); Y_{%s-%s} (cm)", name.data(),
+ dname, dname, name.data(), dname, name.data()),
+ nbinsx, xlo, xhi, nbinsy, ylo, yhi);
break;
- }
-
- case eViewProjection::kPullX:
+ case eProjectionType::kXYt0: xy_id = '0';
+ // fall through
+ case eProjectionType::kXYt1:
+ if (!xy_id) xy_id = '1';
+ // fall through
+ case eProjectionType::kXYt2:
+ if (!xy_id) xy_id = '2';
+ // fall through
+ case eProjectionType::kXYt3:
+ if (!xy_id) xy_id = '3';
+ // fall through
+ case eProjectionType::kXYt4:
+ if (!xy_id) xy_id = '4';
// fall through
- case eViewProjection::kPullY: {
- bool isPullX = (projection.first == eViewProjection::kPullX);
+ case eProjectionType::kXYt5:
+ if (!xy_id) xy_id = '5';
+ nbinsx = 4500; // mm binning
+ nbinsy = 1000;
+ xlo = -25;
+ xhi = +20;
+ ylo = -5;
+ yhi = +5;
+
+ get<2>(projection.second) = new TH2D(Form("hxyt%c_%s%s", xy_id, dname, name.data()),
+ Form("Trk_{proj} z = %+.2f [%s]; X_{%s-%s} (cm); Y_{%s-%s} (cm)", yrange,
+ dname, dname, name.data(), dname, name.data()),
+ nbinsx, xlo, xhi, nbinsy, ylo, yhi);
+ break;
- nbinsx = dscale * ceil(size[0]); // mm binning
- nbinsy = dscale * ceil(size[1]);
- unit = makeYrange(scale, yrange);
+ // Primary vertex plots
+ case eProjectionType::kPVmult: // define primary vertex multiplicity
+ nbinsx = 50;
+ nbinsy = 20;
+ xlo = -0.5;
+ xhi = xlo + nbinsx;
+ ylo = -0.5;
+ yhi = ylo + nbinsy;
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);
+ new TH2D(Form("hMult_%s%s", dname, name.data()), "Vertex multiplicity; N_{trk}^{event}; N_{trk}^{PV}",
+ nbinsx, xlo, xhi, nbinsy, ylo, yhi);
break;
- }
-
- case eViewProjection::kXYp:
- nbinsx = dscale * ceil(size[0] * 2.); // mm binning
- nbinsy = dscale * ceil(size[1] * 2.);
- xlo = pos[0] - size[0];
- xhi = pos[0] + size[0];
- ylo = pos[1] - size[1];
- yhi = pos[1] + size[1];
- get<2>(projection.second) = new TH2D(Form("hxyp_%s%s", dname, name.data()),
- Form("Trk_{proj.} %s [%s]; X_{%s-%s} (cm); Y_{%s-%s} (cm)", name.data(),
- dname, dname, name.data(), dname, name.data()),
- nbinsx, xlo, xhi, nbinsy, ylo, yhi);
+ case eProjectionType::kPVxz: // define X-Z vertex projection
+ xy_id = 'y';
+ mc_id[0] = 'z'; // use to store x-axis title
+ mc_id[1] = 'x'; // use to store y-axis title
+ nbinsx = 100;
+ nbinsy = 200;
+ xlo = -10;
+ xhi = +10;
+ ylo = -2;
+ yhi = +2;
+ // fall through
+ case eProjectionType::kPVyz: // define Y-Z vertex projection
+ if (!xy_id) {
+ xy_id = 'x';
+ mc_id[0] = 'z'; // use to store x-axis title
+ mc_id[1] = 'y'; // use to store y-axis title
+ nbinsx = 100;
+ nbinsy = 200;
+ xlo = -10;
+ xhi = +10;
+ ylo = -2;
+ yhi = +2;
+ }
+ // fall through
+ case eProjectionType::kPVxy: // define X-Y vertex projection
+ if (!xy_id) {
+ xy_id = 'z';
+ mc_id[0] = 'x'; // use to store x-axis title
+ mc_id[1] = 'y'; // use to store y-axis title
+ nbinsx = 200;
+ nbinsy = 200;
+ xlo = -2;
+ xhi = +2;
+ ylo = -2;
+ yhi = +2;
+ }
+ get<2>(projection.second) =
+ new TH2D(Form("h%c_%s%s", xy_id, dname, name.data()),
+ Form("Vertex_{%c%c}; %c (cm); %c (cm)", mc_id[0], mc_id[1], mc_id[0], mc_id[1]), nbinsx, xlo, xhi,
+ nbinsy, ylo, yhi);
break;
}
}
return true;
}
-string CbmRecoQaTask::Detector::View::ToString() const
+string CbmRecoQaTask::View::ToString() const
{
stringstream s;
s << "V[" << name << "] path[" << path << "] @ [" << pos[0] << ", " << pos[1] << ", " << pos[2] << "] size["
@@ -1547,89 +1814,101 @@ string CbmRecoQaTask::Detector::View::ToString() const
return s.str();
}
-uint CbmRecoQaTask::Detector::View::Register(TDirectoryFile* fOut)
+uint CbmRecoQaTask::View::Register(TDirectoryFile* fOut)
{
uint n(0);
TDirectoryFile* lDir = (TDirectoryFile*) fOut->mkdir(name.data(), Form("QA for vol[%s]", path.data()));
- TDirectoryFile* sDir =
- (!CbmRecoQaTask::fuRecoConfig[kRecoEvents] ? (TDirectoryFile*) lDir->mkdir("TS", "TimeSlice QA") : nullptr);
- TDirectoryFile* eDir =
- (CbmRecoQaTask::fuRecoConfig[kRecoEvents] ? (TDirectoryFile*) lDir->mkdir("event", "Local Reco QA") : nullptr);
- TDirectoryFile* tDir =
- (CbmRecoQaTask::fuRecoConfig[kRecoTracks] ? (TDirectoryFile*) lDir->mkdir("track", "CA QA") : nullptr);
- LOG(debug) << "CbmRecoQaTask::Detector::View[" << name << "]::Register() : "
- << "\n\t level \"TS \" [" << (sDir ? "ON" : "OFF") << "]"
- << "\n\t level \"EVT\" [" << (eDir ? "ON" : "OFF") << "]"
- << "\n\t level \"TRK\" [" << (tDir ? "ON" : "OFF") << "]";
+ TDirectoryFile *sDir(nullptr), *tDir(nullptr);
+ // run extra directory structure for the case of detector view
+ if (fType == eViewType::kDetUnit) {
+ if (CbmRecoQaTask::fuRecoConfig[kRecoEvents])
+ sDir = (TDirectoryFile*) lDir->mkdir("event", "Local Reco QA");
+ else
+ sDir = (TDirectoryFile*) lDir->mkdir("TS", "TimeSlice QA");
+ tDir = (CbmRecoQaTask::fuRecoConfig[kRecoTracks] ? (TDirectoryFile*) lDir->mkdir("track", "CA QA") : nullptr);
+ }
for (auto& projection : fProjection) {
TH2* hh = get<2>(projection.second);
if (!hh) continue;
switch (projection.first) {
- 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:
+ case eProjectionType::kXYa:
+ case eProjectionType::kXYp:
+ case eProjectionType::kXdX:
+ case eProjectionType::kXdXMC:
+ case eProjectionType::kYdY:
+ case eProjectionType::kYdYMC:
+ case eProjectionType::kWdT:
+ case eProjectionType::kXpX:
+ case eProjectionType::kYpY:
if (tDir) tDir->Add(hh);
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);
+ case eProjectionType::kChdT:
+ case eProjectionType::kXYh:
+ case eProjectionType::kXYhMC:
+ case eProjectionType::kDmult:
+ case eProjectionType::kDmultMC:
+ case eProjectionType::kResidualX:
+ case eProjectionType::kResidualY:
+ case eProjectionType::kResidualTX:
+ case eProjectionType::kResidualTY:
+ case eProjectionType::kPullX:
+ case eProjectionType::kPullY:
if (sDir) sDir->Add(hh);
break;
- case eViewProjection::kXYs:
- if (sDir) sDir->Add(hh);
+ case eProjectionType::kXYt0:
+ case eProjectionType::kXYt1:
+ case eProjectionType::kXYt2:
+ case eProjectionType::kXYt3:
+ case eProjectionType::kXYt4:
+ case eProjectionType::kXYt5:
+ case eProjectionType::kPVxz:
+ case eProjectionType::kPVyz:
+ case eProjectionType::kPVxy:
+ case eProjectionType::kPVmult:
+ if (lDir) lDir->Add(hh);
break;
}
n++;
}
- LOG(debug) << "CbmRecoQaTask::Detector::View[" << name << "]::Register() : " << n << " projections for "
- << fOut->GetName();
+ LOG(debug) << "CbmRecoQaTask::View[" << name << "]::Register() : " << n << " projections for " << fOut->GetName();
return n;
}
-string CbmRecoQaTask::Detector::View::ToString(eViewProjection prj)
+string CbmRecoQaTask::View::ToString(eProjectionType prj)
{
string s;
switch (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;
+ case eProjectionType::kXYa: s = "x-y (attach)"; break;
+ case eProjectionType::kXYp:
+ case eProjectionType::kXYt0:
+ case eProjectionType::kXYt1:
+ case eProjectionType::kXYt2:
+ case eProjectionType::kXYt3:
+ case eProjectionType::kXYt4:
+ case eProjectionType::kXYt5: s = "x-y (track)"; break;
+ case eProjectionType::kXdX: s = "dx-x"; break;
+ case eProjectionType::kXdXMC: s = "dx-x (MC)"; break;
+ case eProjectionType::kResidualX: s = "residual(x)-x:MC"; break;
+ case eProjectionType::kResidualY: s = "residual(y)-y:MC"; break;
+ case eProjectionType::kResidualTX: s = "residual(t)-x:MC"; break;
+ case eProjectionType::kResidualTY: s = "residual(t)-y:MC"; break;
+ case eProjectionType::kPullX: s = "pull(x)-x:MC"; break;
+ case eProjectionType::kPullY: s = "pull(y)-y:MC"; break;
+ case eProjectionType::kYdY: s = "dy-y"; break;
+ case eProjectionType::kYdYMC: s = "dy-y (MC)"; break;
+ case eProjectionType::kXpX: s = "pull(x)-x"; break;
+ case eProjectionType::kYpY: s = "pull(y)-y"; break;
+ case eProjectionType::kWdT: s = "dt(trk)-w"; break;
+ case eProjectionType::kChdT: s = "dt(ev)-w"; break;
+ case eProjectionType::kXYh: s = "x-y (hit)"; break;
+ case eProjectionType::kXYhMC: s = "x-y (point)"; break;
+ default: LOG(error) << "View::ToString() : Unknown projection " << int(prj); break;
}
return s;
}
-string CbmRecoQaTask::Detector::View::makeYrange(const int scale, float& yrange)
+string CbmRecoQaTask::View::makeYrange(const int scale, float& yrange)
{
bool kDefaultRange = false;
string unit = "cm";
@@ -1649,7 +1928,7 @@ string CbmRecoQaTask::Detector::View::makeYrange(const int scale, float& yrange)
return unit;
}
-string CbmRecoQaTask::Detector::View::makeTrange(const int scale, float& yrange)
+string CbmRecoQaTask::View::makeTrange(const int scale, float& yrange)
{
bool kDefaultRange = false;
string unit = "ns";
@@ -1918,7 +2197,7 @@ std::string CbmRecoQaTask::TrackFilter::ToString() const
/* clang-format off */
ClassImp(CbmRecoQaTask)
ClassImp(CbmRecoQaTask::Detector)
-ClassImp(CbmRecoQaTask::Detector::View)
+ClassImp(CbmRecoQaTask::View)
ClassImp(CbmRecoQaTask::TrackFilter)
ClassImp(CbmRecoQaTask::EventFilter)
/* clang-format on */
diff --git a/reco/qa/CbmRecoQaTask.h b/reco/qa/CbmRecoQaTask.h
index f36c458da682e6125e4269081831c6ae655c1929..54495b14fd63b3a1a9d90bf286c4a9203c07990a 100644
--- a/reco/qa/CbmRecoQaTask.h
+++ b/reco/qa/CbmRecoQaTask.h
@@ -17,9 +17,7 @@
#include <map>
#include <vector>
-#define kRecoQaNConfigs 8
class FairMCPoint;
-
class CbmEvent;
class CbmHit;
class CbmMCDataManager;
@@ -28,6 +26,7 @@ class CbmTimeSlice;
class TClonesArray;
class TH2;
class TH2D;
+class TVector3;
class CbmRecoQaTask : public FairTask {
public:
enum eRecoConfig
@@ -39,6 +38,8 @@ class CbmRecoQaTask : public FairTask {
kTofHits, /// has ToF hits (TofHit branch)
kRichHits, /// has Rich hits (RichHit branch)
kMuchHits, /// has Much hits (MuchHit branch)
+ kUseMC, /// use MC even if available
+ kRecoQaNConfigs /// no of configuration flags
};
enum eSetup
{
@@ -46,10 +47,18 @@ class CbmRecoQaTask : public FairTask {
kMcbm24,
kDefault
};
- enum class eViewProjection : int
+ enum class eViewType : int
+ {
+ kDetUnit = 0, /// detector view
+ kTrkProj, /// set of track projection views
+ kPV /// primary vertex view
+ };
+
+#define kNtrkProjections 6
+ enum class eProjectionType : int
{
kXYa = 0, /// X-Y hit coorelation in track filtered data
- kXYp, /// X-Y track projections
+ kXYp, /// X-Y track projections on detection unit
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
@@ -60,6 +69,8 @@ class CbmRecoQaTask : public FairTask {
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
+ kDmult, /// local view hit multiplicity
+ kDmultMC, /// local view MC point multiplicity
kXYhMC, /// X-Y MC point coorelation in local view (using HitMatch)
kPullX, /// Pull distribution X: (RC - MC) / dx_RC
kPullY, /// Pull distribution Y:
@@ -67,63 +78,82 @@ class CbmRecoQaTask : public FairTask {
kResidualY, /// Residual distribution Y:
kResidualTX, /// Residual distribution T:
kResidualTY, /// Residual distribution T:
- kXYs /// X-Y hit coorelation in time slice
+ kPVxy, /// x-y projection of the primary vertex:
+ kPVxz, /// x-z projection of the primary vertex:
+ kPVyz, /// y-z projection of the primary vertex:
+ kPVmult, /// y-z projection of the primary vertex:
+ kXYt0, /// X-Y track projections on a random plane (value 0)
+ kXYt1, /// X-Y track projections on a random plane (value 1)
+ kXYt2, /// X-Y track projections on a random plane (value 2)
+ kXYt3, /// X-Y track projections on a random plane (value 3)
+ kXYt4, /// X-Y track projections on a random plane (value 4)
+ kXYt5 /// X-Y track projections on a random plane (value 5)
};
+
+ // Generic view definition
+ struct Detector;
+ struct View {
+ friend struct Detector;
+ friend class CbmRecoQaTask;
+
+ View() = default;
+ View(const char* n, const char* p, std::vector<int> set) : name(n), path(p), fSelector(set) { ; }
+ virtual ~View() = default;
+ bool SetProjection(eProjectionType prj, float range, const char* unit);
+ template<class Hit>
+ bool HasAddress(const CbmHit* h, double& x, double& y, double& dx, double& dy) const;
+ template<class Hit>
+ bool Load(const CbmHit* h, const FairMCPoint* point, const CbmEvent* ev);
+ bool Load(const CbmKfTrackFitter::TrajectoryNode* n, const FairMCPoint* point);
+ bool Load(TVector3* p);
+ std::string ToString() const;
+ static std::string ToString(eProjectionType prj);
+ void SetSetup(CbmRecoQaTask::eSetup setup) { fSetup = setup; }
+ void SetType(CbmRecoQaTask::eViewType type) { fType = type; }
+
+ std::string name = ""; /// name describing the module
+ std::string path = ""; /// path to the geo volume describing the module
+ double size[3] = {0.}; /// detection element geometrical dx dy dz dimmensions
+ 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<eProjectionType, 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].
+ int fMult = 0; /// multiplicity between 2 reset signals
+ mutable eSetup fSetup = eSetup::kMcbm22;
+ mutable eViewType fType = eViewType::kDetUnit;
+
+ protected:
+ bool AddProjection(eProjectionType prj, float range = -1, const char* unit = "cm");
+ /** \brief Define all type of QA histo known to the class.
+ * In case of detector view type, 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, bool mc = false);
+ /** \brief build directory structure for all projections of current
+ * view.*/
+ uint Register(TDirectoryFile* f);
+ /** \brief helper functions to estimate the representation (y) axis
+ * \param[in] scale read-only unit defining parameter
+ * \param[in] range read-write full range on the ordinate
+ * \return units for the ordinate and the range value
+ */
+ std::string makeTrange(const int scale, float& range);
+ std::string makeYrange(const int scale, float& range);
+ ClassDef(CbmRecoQaTask::View,
+ 1); // Stand-alone detection set to which QA is applied
+ }; // QA View definition
+
+ // Detector unit definition
struct Detector {
struct Data {
Data() = default;
Data(ECbmDataType i, const char* n) : id(i), name(n) { ; }
ECbmDataType id = ECbmDataType::kUnknown;
std::string name = "nn";
- };
- struct View {
- friend struct Detector;
- friend class CbmRecoQaTask;
-
- View() = default;
- View(const char* n, const char* p, std::vector<int> set) : name(n), path(p), fSelector(set) { ; }
- virtual ~View() = default;
- bool SetProjection(eViewProjection prj, float range, const char* unit);
- template<class Hit>
- bool HasAddress(const CbmHit* h, double& x, double& y, double& dx, double& dy) const;
- template<class Hit>
- bool Load(const CbmHit* h, const FairMCPoint* point, const CbmEvent* ev);
- bool Load(const CbmKfTrackFitter::TrajectoryNode* n, const FairMCPoint* point);
- std::string ToString() const;
- static std::string ToString(eViewProjection prj);
-
- std::string name = ""; /// name describing the module
- std::string path = ""; /// path to the geo volume describing the module
- double size[3] = {0.}; /// detection element geometrical dx dy dz dimmensions
- 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].
- 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.
- * 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.*/
- uint Register(TDirectoryFile* f);
- /** helper functions to estimate the representation (y) axis
- * \param[in] scale read-only unit defining parameter
- * \param[in] range read-write full range on the ordinate
- * \return units for the ordinate and the range value
- */
- 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
- };
+ }; // Identifier of data to be monitored
+
Detector(ECbmModuleId did = ECbmModuleId::kNotExist);
virtual ~Detector() = default;
View* AddView(const char* n, const char* p, std::vector<int> set);
@@ -133,7 +163,7 @@ class CbmRecoQaTask : public FairTask {
* 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);
+ bool Init(TDirectoryFile* f, bool mc = false);
void Print() const;
ECbmModuleId id;
@@ -143,7 +173,7 @@ class CbmRecoQaTask : public FairTask {
ClassDef(CbmRecoQaTask::Detector,
1); // QA representation of a detector unit
- };
+ }; // Detection system agregate
struct TrackFilter {
enum class eTrackCut : int
@@ -172,7 +202,7 @@ class CbmRecoQaTask : public FairTask {
// ToF definition of track
int fNTof = -1;
ClassDef(CbmRecoQaTask::TrackFilter, 1); // Track cut implementation
- };
+ }; // TrackFilter definition
struct EventFilter {
enum class eEventCut : int
@@ -206,7 +236,7 @@ class CbmRecoQaTask : public FairTask {
void HelpMess() const;
ClassDef(CbmRecoQaTask::EventFilter, 1); // Event cut implementation
- };
+ }; // EventFilter definition
public:
CbmRecoQaTask();
@@ -230,13 +260,17 @@ 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 */
- void SetProjections(std::vector<double> vzpj);
+ // /** \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 = "");
+ void UseMC(bool set = true)
+ {
+ if (set) fuRecoConfig.set(kUseMC);
+ }
protected:
static std::bitset<kRecoQaNConfigs> fuRecoConfig;
@@ -255,6 +289,9 @@ class CbmRecoQaTask : public FairTask {
* \return true if track accepted
*/
virtual bool FilterTrack(const CbmGlobalTrack* ptr);
+ /** \brief count views types registered with the task */
+ int GetNviews(eViewType type) const;
+ /** \brief build QA plots for particular setups */
void InitMcbm22();
void InitMcbm24();
void InitDefault();
@@ -268,14 +305,14 @@ class CbmRecoQaTask : public FairTask {
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};
-
+ CbmMCDataManager* cbm_mc_manager = nullptr;
TDirectoryFile fOutFolder = {"RecoQA", "CA track driven reco QA"}; //!
eSetup fSetupClass = eSetup::kMcbm24;
- std::map<double, Detector::View> fProjs = {}; //!
+ std::map<ECbmModuleId, Detector> fDetQa = {}; //! list of detector QA
+ std::map<const char*, View> fViews = {}; //! list of QA views
+ std::vector<TVector3> fPrjPlanes = {}; //! local storage for the z positions of track projection planes
ClassDef(CbmRecoQaTask, 1); // Reconstruction QA analyzed from CA perspective
};