From 374b0d0edeaa0443bd12dc73a314207841ff054e Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 8 Oct 2026 10:19:38 +0200 Subject: [PATCH 1/2] Stream the MC kinematics in the ITS and MFT MC track QC tasks This reduces the memory of ITSTrackSimTask and QcMFTTrackMCTask from the size of all kinematics of the TF to that of the largest event. - The tasks cached all MC tracks of all events and kept a per-track record for each of them. - Events are now read and released one by one, and records are kept only for the tracks that are used. - The cluster-label container of ITSTrackSimTask is no longer leaked. - The histograms are unchanged (45 ITS and 9 MFT histograms compared bin by bin, same entries). Peak memory (sum of PSS) on one PbPb timeframe of 8 orbits (28 signal events, 75.6M MC tracks): | | before | after | |---|---:|---:| | ITS task total | 8.44 GB | 2.40 GB | | MFT task total | 7.80 GB | 2.16 GB | | ITS `o2-qc` device | 7.66 GB | 1.57 GB | | MFT `o2-qc` device | 7.04 GB | 1.40 GB | | wall time ITS / MFT | 33 s / 40 s | 36 s / 40 s | Co-Authored-By: Claude Sonnet 5.5 --- Modules/ITS/src/ITSTrackSimTask.cxx | 287 ++++++++++++++------------- Modules/MFT/src/QcMFTTrackMCTask.cxx | 106 ++++++---- 2 files changed, 214 insertions(+), 179 deletions(-) diff --git a/Modules/ITS/src/ITSTrackSimTask.cxx b/Modules/ITS/src/ITSTrackSimTask.cxx index 8d55234ddd..e2ef2d8f37 100644 --- a/Modules/ITS/src/ITSTrackSimTask.cxx +++ b/Modules/ITS/src/ITSTrackSimTask.cxx @@ -42,6 +42,8 @@ #include "TFile.h" #include "TTree.h" #include +#include +#include using namespace o2::constants::math; using namespace o2::itsmft; using namespace o2::its; @@ -142,83 +144,84 @@ void ITSTrackSimTask::monitorData(o2::framework::ProcessingContext& ctx) { ILOG(Debug, Devel) << "START DOING QC General" << ENDM; o2::steer::MCKinematicsReader reader(mCollisionsContextPath.c_str()); - info.resize(reader.getNEvents(0)); - for (int i = 0; i < reader.getNEvents(0); ++i) { - std::vector const& mcArr = reader.getTracks(i); - info[i].resize(mcArr.size()); - } + const int nEvents = reader.getNEvents(0); auto clusArr = ctx.inputs().get>("compclus"); // used to get hit information - auto clusLabArr = ctx.inputs().get*>("mcclustruth").release(); - - for (int iCluster = 0; iCluster < clusArr.size(); iCluster++) { - + auto clusLabArr = ctx.inputs().get*>("mcclustruth"); + + // layers with a cluster of an MC track, as a bit mask; indexed [event][track], sized from the labels only + std::vector> layerMask(nEvents); + { + std::vector nTracks(nEvents, 0); + for (int iCluster = 0; iCluster < (int)clusArr.size(); iCluster++) { + auto lab = (clusLabArr->getLabels(iCluster))[0]; + if (!lab.isValid() || lab.getSourceID() != 0 || lab.getTrackID() < 0 || !lab.isCorrect()) { + continue; + } + nTracks[lab.getEventID()] = std::max(nTracks[lab.getEventID()], lab.getTrackID() + 1); + } + for (int i = 0; i < nEvents; ++i) { + layerMask[i].assign(nTracks[i], 0); + } + } + for (int iCluster = 0; iCluster < (int)clusArr.size(); iCluster++) { auto lab = (clusLabArr->getLabels(iCluster))[0]; - if (!lab.isValid() || lab.getSourceID() != 0) - continue; - - int TrackID = lab.getTrackID(); - - if (TrackID < 0) { + if (!lab.isValid() || lab.getSourceID() != 0 || lab.getTrackID() < 0 || !lab.isCorrect()) { continue; } - - if (!lab.isCorrect()) - continue; - const auto& Cluster = (clusArr)[iCluster]; - unsigned short& ok = info[lab.getEventID()][lab.getTrackID()].clusters; // bitmask with track hits at each layer - auto layer = mGeom->getLayer(Cluster.getSensorID()); - float r = 0.f; - if (layer == 0) - ok |= 0b1; - if (layer == 1) - ok |= 0b10; - if (layer == 2) - ok |= 0b100; - if (layer == 3) - ok |= 0b1000; - if (layer == 4) - ok |= 0b10000; - if (layer == 5) - ok |= 0b100000; - if (layer == 6) - ok |= 0b1000000; + const auto layer = mGeom->getLayer(clusArr[iCluster].getSensorID()); + if (layer >= 0 && layer <= 6) { + layerMask[lab.getEventID()][lab.getTrackID()] |= (1 << layer); // bitmask with track hits at each layer + } } - for (int i = 0; i < reader.getNEvents(0); ++i) { - std::vector const& mcArr = reader.getTracks(i); - auto mcHeader = reader.getMCEventHeader(0, i); // SourceID=0 for ITS - - for (int mc = 0; mc < mcArr.size(); mc++) { - const auto& mcTrack = (mcArr)[mc]; - - info[i][mc].isFilled = false; - if (mcTrack.Vx() * mcTrack.Vx() + mcTrack.Vy() * mcTrack.Vy() > 1) - continue; - if (TMath::Abs(mcTrack.GetPdgCode()) != 211) - continue; // Select pions - if (TMath::Abs(mcTrack.GetEta()) > 1.2) - continue; - if (info[i][mc].clusters != 0b1111111) - continue; - Double_t distance = sqrt(pow(mcHeader.GetX() - mcTrack.Vx(), 2) + pow(mcHeader.GetY() - mcTrack.Vy(), 2) + pow(mcHeader.GetZ() - mcTrack.Vz(), 2)); - info[i][mc].isFilled = true; - info[i][mc].r = distance; - info[i][mc].pt = mcTrack.GetPt(); - info[i][mc].eta = mcTrack.GetEta(); - info[i][mc].phi = mcTrack.GetPhi(); - info[i][mc].z = mcTrack.Vz(); - info[i][mc].isPrimary = mcTrack.isPrimary(); - if (mcTrack.isPrimary()) { - hPrimaryGen_pt->Fill(mcTrack.GetPt()); - // True Generated primaries: denominator of the efficiency plots - hDenTrue_r[4]->Fill(distance); - hDenTrue_pt[4]->Fill(mcTrack.GetPt()); - hDenTrue_eta[4]->Fill(mcTrack.GetEta()); - hDenTrue_phi[4]->Fill(mcTrack.GetPhi()); - hDenTrue_z[4]->Fill(mcTrack.Vz()); + // MC truth of the tracks passing the selection, keyed by (event, track); the kinematics are read one event at a time + struct SelectedTrack { + float r, pt, eta, phi, z; + bool isPrimary; + int isReco = 0; + }; + std::unordered_map selected; + auto keyOf = [](int event, int track) { return (uint64_t(uint32_t(event)) << 32) | uint32_t(track); }; + + for (int i = 0; i < nEvents; ++i) { + { + std::vector const& mcArr = reader.getTracks(i); + auto mcHeader = reader.getMCEventHeader(0, i); // SourceID=0 for ITS + const auto& mask = layerMask[i]; + + for (int mc = 0; mc < (int)mcArr.size(); mc++) { + const auto& mcTrack = (mcArr)[mc]; + + if (mc >= (int)mask.size() || mask[mc] != 0b1111111) + continue; + if (mcTrack.Vx() * mcTrack.Vx() + mcTrack.Vy() * mcTrack.Vy() > 1) + continue; + if (TMath::Abs(mcTrack.GetPdgCode()) != 211) + continue; // Select pions + if (TMath::Abs(mcTrack.GetEta()) > 1.2) + continue; + Double_t distance = sqrt(pow(mcHeader.GetX() - mcTrack.Vx(), 2) + pow(mcHeader.GetY() - mcTrack.Vy(), 2) + pow(mcHeader.GetZ() - mcTrack.Vz(), 2)); + SelectedTrack& sel = selected[keyOf(i, mc)]; + sel.r = distance; + sel.pt = mcTrack.GetPt(); + sel.eta = mcTrack.GetEta(); + sel.phi = mcTrack.GetPhi(); + sel.z = mcTrack.Vz(); + sel.isPrimary = mcTrack.isPrimary(); + if (mcTrack.isPrimary()) { + hPrimaryGen_pt->Fill(mcTrack.GetPt()); + // True Generated primaries: denominator of the efficiency plots + hDenTrue_r[4]->Fill(distance); + hDenTrue_pt[4]->Fill(mcTrack.GetPt()); + hDenTrue_eta[4]->Fill(mcTrack.GetEta()); + hDenTrue_phi[4]->Fill(mcTrack.GetPhi()); + hDenTrue_z[4]->Fill(mcTrack.Vz()); + } } } + reader.releaseTracksForSourceAndEvent(0, i); // the tracks of this event are not needed any more + std::vector().swap(layerMask[i]); } auto trackArr = ctx.inputs().get>("tracks"); // MC Tracks @@ -236,104 +239,106 @@ void ITSTrackSimTask::monitorData(o2::framework::ProcessingContext& ctx) hAngularDistribution->Fill(track.getEta(), track.getPhi()); - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isFilled) { + auto selIt = selected.find(keyOf(MCinfo.getEventID(), MCinfo.getTrackID())); + if (selIt != selected.end()) { + SelectedTrack& rec = selIt->second; - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isPrimary) { + if (rec.isPrimary) { // True Generated primaries for QoverPt plot, because MCTrack does not have charge function - hDenTrue_QoverPt[4]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hDenTrue_QoverPt[4]->Fill(track.getSign() / rec.pt); if (iNClusters == 4) { - hDenTrue_QoverPt[0]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_r[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hDenTrue_pt[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_eta[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hDenTrue_phi[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hDenTrue_z[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); + hDenTrue_QoverPt[0]->Fill(track.getSign() / rec.pt); + hDenTrue_r[0]->Fill(rec.r); + hDenTrue_pt[0]->Fill(rec.pt); + hDenTrue_eta[0]->Fill(rec.eta); + hDenTrue_phi[0]->Fill(rec.phi); + hDenTrue_z[0]->Fill(rec.z); } else if (iNClusters == 5) { - hDenTrue_QoverPt[1]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_r[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hDenTrue_pt[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_eta[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hDenTrue_phi[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hDenTrue_z[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); + hDenTrue_QoverPt[1]->Fill(track.getSign() / rec.pt); + hDenTrue_r[1]->Fill(rec.r); + hDenTrue_pt[1]->Fill(rec.pt); + hDenTrue_eta[1]->Fill(rec.eta); + hDenTrue_phi[1]->Fill(rec.phi); + hDenTrue_z[1]->Fill(rec.z); } else if (iNClusters == 6) { - hDenTrue_QoverPt[2]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_r[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hDenTrue_pt[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_eta[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hDenTrue_phi[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hDenTrue_z[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); + hDenTrue_QoverPt[2]->Fill(track.getSign() / rec.pt); + hDenTrue_r[2]->Fill(rec.r); + hDenTrue_pt[2]->Fill(rec.pt); + hDenTrue_eta[2]->Fill(rec.eta); + hDenTrue_phi[2]->Fill(rec.phi); + hDenTrue_z[2]->Fill(rec.z); } else if (iNClusters == 7) { - hDenTrue_QoverPt[3]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_r[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hDenTrue_pt[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hDenTrue_eta[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hDenTrue_phi[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hDenTrue_z[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); + hDenTrue_QoverPt[3]->Fill(track.getSign() / rec.pt); + hDenTrue_r[3]->Fill(rec.r); + hDenTrue_pt[3]->Fill(rec.pt); + hDenTrue_eta[3]->Fill(rec.eta); + hDenTrue_phi[3]->Fill(rec.phi); + hDenTrue_z[3]->Fill(rec.z); } } - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isReco != 0) { - hNumDuplicate_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hNumDuplicate_phi->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hNumDuplicate_eta->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hNumDuplicate_z->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); - hNumDuplicate_r->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); + if (rec.isReco != 0) { + hNumDuplicate_pt->Fill(rec.pt); + hNumDuplicate_phi->Fill(rec.phi); + hNumDuplicate_eta->Fill(rec.eta); + hNumDuplicate_z->Fill(rec.z); + hNumDuplicate_r->Fill(rec.r); continue; } - info[MCinfo.getEventID()][MCinfo.getTrackID()].isReco++; + rec.isReco++; if (MCinfo.isFake()) { - hNumRecoFake_pt[4]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hNumRecoFake_phi[4]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hNumRecoFake_eta[4]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hNumRecoFake_z[4]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); - hNumRecoFake_r[4]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hNumRecoFake_QoverPt[4]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hNumRecoFake_pt[4]->Fill(rec.pt); + hNumRecoFake_phi[4]->Fill(rec.phi); + hNumRecoFake_eta[4]->Fill(rec.eta); + hNumRecoFake_z[4]->Fill(rec.z); + hNumRecoFake_r[4]->Fill(rec.r); + hNumRecoFake_QoverPt[4]->Fill(track.getSign() / rec.pt); if (iNClusters == 4) { - hNumRecoFake_pt[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hNumRecoFake_phi[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hNumRecoFake_eta[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hNumRecoFake_z[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); - hNumRecoFake_r[0]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hNumRecoFake_QoverPt[0]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hNumRecoFake_pt[0]->Fill(rec.pt); + hNumRecoFake_phi[0]->Fill(rec.phi); + hNumRecoFake_eta[0]->Fill(rec.eta); + hNumRecoFake_z[0]->Fill(rec.z); + hNumRecoFake_r[0]->Fill(rec.r); + hNumRecoFake_QoverPt[0]->Fill(track.getSign() / rec.pt); } else if (iNClusters == 5) { - hNumRecoFake_pt[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hNumRecoFake_phi[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hNumRecoFake_eta[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hNumRecoFake_z[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); - hNumRecoFake_r[1]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hNumRecoFake_QoverPt[1]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hNumRecoFake_pt[1]->Fill(rec.pt); + hNumRecoFake_phi[1]->Fill(rec.phi); + hNumRecoFake_eta[1]->Fill(rec.eta); + hNumRecoFake_z[1]->Fill(rec.z); + hNumRecoFake_r[1]->Fill(rec.r); + hNumRecoFake_QoverPt[1]->Fill(track.getSign() / rec.pt); } else if (iNClusters == 6) { - hNumRecoFake_pt[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hNumRecoFake_phi[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hNumRecoFake_eta[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hNumRecoFake_z[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); - hNumRecoFake_r[2]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hNumRecoFake_QoverPt[2]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hNumRecoFake_pt[2]->Fill(rec.pt); + hNumRecoFake_phi[2]->Fill(rec.phi); + hNumRecoFake_eta[2]->Fill(rec.eta); + hNumRecoFake_z[2]->Fill(rec.z); + hNumRecoFake_r[2]->Fill(rec.r); + hNumRecoFake_QoverPt[2]->Fill(track.getSign() / rec.pt); } else if (iNClusters == 7) { - hNumRecoFake_pt[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hNumRecoFake_phi[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hNumRecoFake_eta[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hNumRecoFake_z[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); - hNumRecoFake_r[3]->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); - hNumRecoFake_QoverPt[3]->Fill(track.getSign() / info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hNumRecoFake_pt[3]->Fill(rec.pt); + hNumRecoFake_phi[3]->Fill(rec.phi); + hNumRecoFake_eta[3]->Fill(rec.eta); + hNumRecoFake_z[3]->Fill(rec.z); + hNumRecoFake_r[3]->Fill(rec.r); + hNumRecoFake_QoverPt[3]->Fill(track.getSign() / rec.pt); } - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isPrimary) { + if (rec.isPrimary) { hTrackImpactTransvFake->Fill(ip[0]); - hPrimaryReco_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hPrimaryReco_pt->Fill(rec.pt); } } else { - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isPrimary) { + if (rec.isPrimary) { // True primaries reconstructed: numerator of efficiency plots - hNumRecoValid_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hNumRecoValid_phi->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hNumRecoValid_eta->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - hNumRecoValid_z->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].z); - hNumRecoValid_r->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].r); + hNumRecoValid_pt->Fill(rec.pt); + hNumRecoValid_phi->Fill(rec.phi); + hNumRecoValid_eta->Fill(rec.eta); + hNumRecoValid_z->Fill(rec.z); + hNumRecoValid_r->Fill(rec.r); hTrackImpactTransvValid->Fill(ip[0]); - hPrimaryReco_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hPrimaryReco_pt->Fill(rec.pt); } } } diff --git a/Modules/MFT/src/QcMFTTrackMCTask.cxx b/Modules/MFT/src/QcMFTTrackMCTask.cxx index 7f062ce81b..3c48891cab 100644 --- a/Modules/MFT/src/QcMFTTrackMCTask.cxx +++ b/Modules/MFT/src/QcMFTTrackMCTask.cxx @@ -26,6 +26,7 @@ #include #include #include +#include #include //ADDED FOR MC FILE READING #include #include @@ -122,57 +123,86 @@ void QcMFTTrackMCTask::monitorData(o2::framework::ProcessingContext& ctx) { ILOG(Debug, Devel) << "START DOING QC General" << ENDM; o2::steer::MCKinematicsReader reader(mCollisionsContextPath.c_str()); - info.resize(reader.getNEvents(0)); - for (int i = 0; i < reader.getNEvents(0); ++i) { - std::vector const& mcArr = reader.getTracks(i); - info[i].resize(mcArr.size()); + const int nEvents = reader.getNEvents(0); + + auto trackArr = ctx.inputs().get>("tracks"); // MFT Tracks + auto MCTruth = ctx.inputs().get>("mctruth"); // MC track label, contains info about EventID, TrackID, SourceID etc + + // the MC tracks referenced by the reconstructed tracks, per event and sorted by track ID + std::vector> wanted(nEvents); + for (int itrack = 0; itrack < (int)trackArr.size(); itrack++) { + const auto& MCinfo = MCTruth[itrack]; + if (MCinfo.isNoise() || !MCinfo.isValid()) { + continue; + } + wanted[MCinfo.getEventID()].push_back(MCinfo.getTrackID()); + } + for (auto& w : wanted) { + std::sort(w.begin(), w.end()); + w.erase(std::unique(w.begin(), w.end()), w.end()); } - for (int i = 0; i < reader.getNEvents(0); ++i) { - std::vector const& mcArr = reader.getTracks(i); - auto mcHeader = reader.getMCEventHeader(0, i); - for (int mc = 0; mc < mcArr.size(); mc++) { - const auto& mcTrack = (mcArr)[mc]; - info[i][mc].isFilled = false; - info[i][mc].isFilled = true; - info[i][mc].pt = mcTrack.GetPt(); - info[i][mc].eta = mcTrack.GetEta(); - info[i][mc].phi = TMath::ATan2(mcTrack.Py(), mcTrack.Px()); - info[i][mc].isPrimary = mcTrack.isPrimary(); - if (mcTrack.isPrimary()) { - hPrimaryGen_pt->Fill(mcTrack.GetPt()); + // MC truth of the referenced tracks, parallel to 'wanted'; the kinematics are read one event at a time + struct MCTruthInfo { + float pt = 0, eta = 0, phi = 0; + bool isPrimary = false; + bool isFilled = false; + int isReco = 0; + }; + std::vector> truth(nEvents); + + for (int i = 0; i < nEvents; ++i) { + { + std::vector const& mcArr = reader.getTracks(i); + truth[i].resize(wanted[i].size()); + size_t iw = 0; + for (int mc = 0; mc < (int)mcArr.size(); mc++) { + const auto& mcTrack = (mcArr)[mc]; + const double phi = TMath::ATan2(mcTrack.Py(), mcTrack.Px()); + if (iw < wanted[i].size() && wanted[i][iw] == mc) { + auto& t = truth[i][iw++]; + t.isFilled = true; + t.pt = mcTrack.GetPt(); + t.eta = mcTrack.GetEta(); + t.phi = phi; + t.isPrimary = mcTrack.isPrimary(); + } + if (mcTrack.isPrimary()) { + hPrimaryGen_pt->Fill(mcTrack.GetPt()); + } + hTrue_pt->Fill(mcTrack.GetPt()); + hTrue_eta->Fill(mcTrack.GetEta()); + hTrue_phi->Fill(phi); } - hTrue_pt->Fill(mcTrack.GetPt()); - hTrue_eta->Fill(mcTrack.GetEta()); - hTrue_phi->Fill(TMath::ATan2(mcTrack.Py(), mcTrack.Px())); } + reader.releaseTracksForSourceAndEvent(0, i); // the tracks of this event are not needed any more } - auto trackArr = ctx.inputs().get>("tracks"); // MFT Tracks - auto MCTruth = ctx.inputs().get>("mctruth"); // MC track label, contains info about EventID, TrackID, SourceID etc - - for (int itrack = 0; itrack < trackArr.size(); itrack++) { + for (int itrack = 0; itrack < (int)trackArr.size(); itrack++) { const auto& track = trackArr[itrack]; const auto& MCinfo = MCTruth[itrack]; - if (MCinfo.isNoise()) + if (MCinfo.isNoise() || !MCinfo.isValid()) continue; - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isFilled) { - info[MCinfo.getEventID()][MCinfo.getTrackID()].isReco++; + const auto& w = wanted[MCinfo.getEventID()]; + const auto iw = std::lower_bound(w.begin(), w.end(), MCinfo.getTrackID()) - w.begin(); + auto& rec = truth[MCinfo.getEventID()][iw]; + if (rec.isFilled) { + rec.isReco++; if (MCinfo.isFake()) { - hRecoFake_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hRecoFake_phi->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hRecoFake_eta->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isPrimary) { - hPrimaryReco_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); + hRecoFake_pt->Fill(rec.pt); + hRecoFake_phi->Fill(rec.phi); + hRecoFake_eta->Fill(rec.eta); + if (rec.isPrimary) { + hPrimaryReco_pt->Fill(rec.pt); } } else { - hRecoValid_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hRecoValid_phi->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].phi); - hRecoValid_eta->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].eta); - if (info[MCinfo.getEventID()][MCinfo.getTrackID()].isPrimary) { - hPrimaryReco_pt->Fill(info[MCinfo.getEventID()][MCinfo.getTrackID()].pt); - hResolution_pt->Fill(1 / (track.getPt()) - 1 / (info[MCinfo.getEventID()][MCinfo.getTrackID()].pt)); + hRecoValid_pt->Fill(rec.pt); + hRecoValid_phi->Fill(rec.phi); + hRecoValid_eta->Fill(rec.eta); + if (rec.isPrimary) { + hPrimaryReco_pt->Fill(rec.pt); + hResolution_pt->Fill(1 / (track.getPt()) - 1 / (rec.pt)); } } } From f51d2515ffb8a0abfd23c47b9d0bd662b7c2ef84 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 8 Oct 2026 12:38:40 +0200 Subject: [PATCH 2/2] Bound the event ID in the MC track QC tasks This skips MC labels whose event ID is outside the event range of the kinematics in ITSTrackSimTask and QcMFTTrackMCTask. - Such labels would index the per-event arrays out of range. --- Modules/ITS/src/ITSTrackSimTask.cxx | 4 ++-- Modules/MFT/src/QcMFTTrackMCTask.cxx | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Modules/ITS/src/ITSTrackSimTask.cxx b/Modules/ITS/src/ITSTrackSimTask.cxx index e2ef2d8f37..7feafd39f1 100644 --- a/Modules/ITS/src/ITSTrackSimTask.cxx +++ b/Modules/ITS/src/ITSTrackSimTask.cxx @@ -155,7 +155,7 @@ void ITSTrackSimTask::monitorData(o2::framework::ProcessingContext& ctx) std::vector nTracks(nEvents, 0); for (int iCluster = 0; iCluster < (int)clusArr.size(); iCluster++) { auto lab = (clusLabArr->getLabels(iCluster))[0]; - if (!lab.isValid() || lab.getSourceID() != 0 || lab.getTrackID() < 0 || !lab.isCorrect()) { + if (!lab.isValid() || lab.getSourceID() != 0 || lab.getTrackID() < 0 || lab.getEventID() >= nEvents || !lab.isCorrect()) { continue; } nTracks[lab.getEventID()] = std::max(nTracks[lab.getEventID()], lab.getTrackID() + 1); @@ -166,7 +166,7 @@ void ITSTrackSimTask::monitorData(o2::framework::ProcessingContext& ctx) } for (int iCluster = 0; iCluster < (int)clusArr.size(); iCluster++) { auto lab = (clusLabArr->getLabels(iCluster))[0]; - if (!lab.isValid() || lab.getSourceID() != 0 || lab.getTrackID() < 0 || !lab.isCorrect()) { + if (!lab.isValid() || lab.getSourceID() != 0 || lab.getTrackID() < 0 || lab.getEventID() >= nEvents || !lab.isCorrect()) { continue; } const auto layer = mGeom->getLayer(clusArr[iCluster].getSensorID()); diff --git a/Modules/MFT/src/QcMFTTrackMCTask.cxx b/Modules/MFT/src/QcMFTTrackMCTask.cxx index 3c48891cab..6a131224cb 100644 --- a/Modules/MFT/src/QcMFTTrackMCTask.cxx +++ b/Modules/MFT/src/QcMFTTrackMCTask.cxx @@ -132,7 +132,7 @@ void QcMFTTrackMCTask::monitorData(o2::framework::ProcessingContext& ctx) std::vector> wanted(nEvents); for (int itrack = 0; itrack < (int)trackArr.size(); itrack++) { const auto& MCinfo = MCTruth[itrack]; - if (MCinfo.isNoise() || !MCinfo.isValid()) { + if (MCinfo.isNoise() || !MCinfo.isValid() || MCinfo.getEventID() >= nEvents) { continue; } wanted[MCinfo.getEventID()].push_back(MCinfo.getTrackID()); @@ -181,7 +181,7 @@ void QcMFTTrackMCTask::monitorData(o2::framework::ProcessingContext& ctx) for (int itrack = 0; itrack < (int)trackArr.size(); itrack++) { const auto& track = trackArr[itrack]; const auto& MCinfo = MCTruth[itrack]; - if (MCinfo.isNoise() || !MCinfo.isValid()) + if (MCinfo.isNoise() || !MCinfo.isValid() || MCinfo.getEventID() >= nEvents) continue; const auto& w = wanted[MCinfo.getEventID()];