Skip to content

Commit 5faf8b0

Browse files
lietavalietava
committed
dev: fixing mem leaks and optimising speed
1 parent 4eb919f commit 5faf8b0

File tree

1 file changed

+170
-104
lines changed

1 file changed

+170
-104
lines changed

‎PWGLF/Tasks/Strangeness/nonPromptCascade.cxx‎

Lines changed: 170 additions & 104 deletions
Original file line numberDiff line numberDiff line change
@@ -281,7 +281,7 @@ struct NonPromptCascadeTask {
281281
int nBinsMultFV0 = cfgMaxMultFV0;
282282

283283
AxisSpec multAxis = {nBinsMult, 0, cfgMaxMult, "Multiplicity FT0M"};
284-
AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FT0M"};
284+
AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FV0"};
285285
AxisSpec nTracksAxis = {100, 0., 100., "NTracksGlobal"};
286286
AxisSpec nTracksAxisMC = {100, 0., 100., "NTracksMC"};
287287
std::vector<double> centBinning;
@@ -311,13 +311,13 @@ struct NonPromptCascadeTask {
311311
runsBinning.push_back(run);
312312
run++;
313313
}
314-
AxisSpec centAxis{centBinning, "Centrality (%)"};
315-
// AxisSpec multAxis{multBinning, "Multiplicity"};
314+
AxisSpec centAxisFT0M{centBinning, "Centrality FT0M (%)"};
315+
AxisSpec centAxisFV0{centBinning, "Centrality FV0 (%)"};
316316
AxisSpec trackAxisMC{trackBinning, "NTracks MC"};
317317
AxisSpec trackAxis{trackBinning, "NTracks Global Reco"};
318318
AxisSpec runsAxis{runsBinning, "Run Number"};
319319

320-
mRegistryMults.add("hCentMultsRuns", "hCentMultsRuns", HistType::kTHnSparseF, {centAxis, multAxis, centAxis, multAxisFV0, nTracksAxis, runsAxis});
320+
mRegistryMults.add("hCentMultsRuns", "hCentMultsRuns", HistType::kTHnSparseF, {centAxisFT0M, multAxis, centAxisFV0, multAxisFV0, nTracksAxis, runsAxis});
321321
//
322322
// dN/deta
323323
//
@@ -792,114 +792,180 @@ struct NonPromptCascadeTask {
792792
}
793793
PROCESS_SWITCH(NonPromptCascadeTask, processCascadesData, "process cascades: Data analysis", false);
794794

795-
int getMCMult(aod::McParticles const& mcParticles, int mcCollId, std::vector<float>& ptMCvec)
796-
{
797-
int mult = 0;
798-
for (auto const& mcp : mcParticles) {
799-
if (mcp.mcCollisionId() == mcCollId) {
800-
// multiplicity definition:
801-
bool accept = mcp.isPhysicalPrimary();
802-
accept = accept && (mcp.eta() < 0.5) && (mcp.eta() > -0.5);
803-
int q = 0;
804-
auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcp.pdgCode());
805-
if (pdgEntry) {
806-
q = int(std::round(pdgEntry->Charge() / 3.0));
807-
} else {
808-
// LOG(warn) << "No pdg assuming neutral";
809-
}
810-
accept = accept && (q != 0);
811-
if (accept) {
812-
++mult;
813-
ptMCvec.push_back(mcp.pt());
814-
}
815-
}
795+
// colls : Join<aod::Collisions, ...>
796+
// tracks: Join<aod::Tracks, aod::McTrackLabels>
797+
// mcCollisions: aod::McCollisions
798+
// mcParticles : aod::McParticles
799+
800+
void processdNdetaMC(CollisionCandidatesRun3MC const& colls,
801+
aod::McCollisions const& mcCollisions,
802+
aod::McParticles const& mcParticles,
803+
TracksWithLabel const& tracks)
804+
{
805+
//-------------------------------------------------------------
806+
// MC mult for all MC coll
807+
//--------------------------------------------------------------
808+
std::vector<int> mcMult(mcCollisions.size(), 0);
809+
for (auto const& mcp : mcParticles) {
810+
int mcid = mcp.mcCollisionId();
811+
if (mcid < 0 || mcid >= (int)mcMult.size()) continue;
812+
813+
// apply your primary/eta/charge definition here
814+
if (!mcp.isPhysicalPrimary()) continue;
815+
if (std::abs(mcp.eta()) > 0.5f) continue;
816+
int q = 0;
817+
if (auto pdg = TDatabasePDG::Instance()->GetParticle(mcp.pdgCode())) {
818+
q = int(std::round(pdg->Charge() / 3.0));
816819
}
817-
return mult;
820+
if (q == 0) continue;
821+
822+
++mcMult[mcid];
818823
}
819-
void processdNdetaMC(CollisionCandidatesRun3MC const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TracksWithLabel const& tracks)
820-
{
821-
// std::cout << "ProcNegMC" << std::endl;
822-
// Map: collision index -> reco multiplicity
823-
std::vector<int> recoMult(colls.size(), 0);
824-
for (auto const& trk : tracks) {
825-
int collId = trk.collisionId();
826-
// Here you can impose same track cuts as for MC (eta, primaries, etc.)
827-
if (std::abs(trk.eta()) > 0.5f) {
828-
continue;
829-
}
830-
++recoMult[collId];
831-
}
832-
std::vector<int> isReco(mcParticles.size(), 0);
833-
std::vector<int> isRecoMult(mcParticles.size(), 0);
834-
std::vector<int> mcReconstructed(mcCollisions.size(), 0);
835-
// std::cout << " reco cols with mc:" << colls.size() << " tracks:" << tracks.size() << " mccols:" << mcCollisions.size() << "mcParticles:" << mcParticles.size() << std::endl;
836-
for (auto const& col : colls) {
837-
int mcCollId = col.mcCollisionId(); // col.template mcCollision_as<aod::McCollisions>();
838-
int collId = col.globalIndex();
839-
// auto mc = col.mcCollision();
840-
// int mcId = mc.globalIndex();
841-
// std::cout << "globalIndex:" << mcId << " colID:" << mcCollId << std::endl;
842-
std::vector<float> mcptvec;
843-
float mult = getMCMult(mcParticles, mcCollId, mcptvec);
844-
mcReconstructed[mcCollId] = 1;
845-
for (auto const& trk : tracks) {
846-
int mcPid = trk.mcParticleId(); // depends on your label table
847-
if (mcPid >= 0 && mcPid < (int)mcParticles.size()) {
848-
isReco[mcPid] = 1;
849-
isRecoMult[mcPid] = mult;
850-
} else {
851-
continue;
852-
}
853-
if (trk.collisionId() != collId) {
854-
continue; // different event
855-
}
856-
auto mcPar = mcParticles.rawIteratorAt(mcPid);
857-
// Apply same acceptance as in MC multiplicity
858-
if (!mcPar.isPhysicalPrimary()) {
859-
continue;
860-
}
861-
if (std::abs(mcPar.eta()) > 0.5f) {
862-
continue;
863-
}
864-
int q = 0;
865-
auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode());
866-
if (pdgEntry) {
867-
q = int(std::round(pdgEntry->Charge() / 3.0));
868-
}
869-
if (q == 0) {
870-
continue;
871-
}
872-
// float multReco = recoMult[collId];
873-
float multReco = col.multNTracksGlobal();
874-
float ptReco = trk.pt();
875-
float ptMC = mcPar.pt();
876-
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco);
877-
}
878824

879-
// mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal());
825+
// ------------------------------------------------------------
826+
// Build mapping: (aod::Collisions row id used by tracks.collisionId())
827+
// -> dense index in 'colls' (0..colls.size()-1)
828+
// We assume col.globalIndex() refers to the original aod::Collisions row.
829+
// ------------------------------------------------------------
830+
int maxCollRowId = -1;
831+
for (auto const& trk : tracks) {
832+
maxCollRowId = std::max(maxCollRowId, (int)trk.collisionId());
833+
}
834+
std::vector<int> collRowIdToDense(maxCollRowId + 1, -1);
835+
836+
int dense = 0;
837+
for (auto const& col : colls) {
838+
const int collRowId = col.globalIndex(); // row id in aod::Collisions
839+
if (collRowId >= 0 && collRowId < (int)collRowIdToDense.size()) {
840+
collRowIdToDense[collRowId] = dense;
880841
}
881-
// count mc particles with no reco tracks
882-
for (auto const& mcp : mcParticles) {
883-
int mcPidG = mcp.globalIndex();
884-
// std::cout << "mcPidG:" << mcPidG << std::endl;
885-
if (!isReco[mcPidG]) {
886-
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[mcPidG], mcp.pt());
887-
}
842+
++dense;
843+
}
844+
845+
// ------------------------------------------------------------
846+
// Reco multiplicity per *dense collision index in colls*
847+
// ------------------------------------------------------------
848+
std::vector<int> recoMultDense(colls.size(), 0);
849+
for (auto const& trk : tracks) {
850+
if (std::abs(trk.eta()) > 0.5f) {
851+
continue;
888852
}
889-
// count unreconstructed mc collisions
890-
for (auto const& mc : mcCollisions) {
891-
int gindex = mc.globalIndex();
892-
// std::cout << "mc globalIndex:" << gindex << std::endl;
893-
if (!mcReconstructed[gindex]) {
894-
std::vector<float> mcptvec;
895-
int mult = getMCMult(mcParticles, gindex, mcptvec);
896-
// std::cout << "===> unreconstructed:" << mult << std::endl;
897-
for (auto const& pt : mcptvec) {
898-
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt);
899-
}
853+
const int collRowId = trk.collisionId();
854+
if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) {
855+
continue;
856+
}
857+
const int dIdx = collRowIdToDense[collRowId];
858+
if (dIdx >= 0) {
859+
++recoMultDense[dIdx];
860+
}
861+
}
862+
863+
// ------------------------------------------------------------
864+
// MC bookkeeping: index by ROW INDEX (0..size-1), not globalIndex()
865+
// ------------------------------------------------------------
866+
std::vector<char> isReco(mcParticles.size(), 0);
867+
std::vector<float> isRecoMult(mcParticles.size(), 0.f);
868+
std::vector<char> mcReconstructed(mcCollisions.size(), 0);
869+
870+
// Optional cache of MC multiplicity per MC collision
871+
std::vector<float> mcMultCache(mcCollisions.size(), -1.f);
872+
873+
// ------------------------------------------------------------
874+
// Single pass over tracks: fill RM for tracks whose collision is in colls
875+
// ------------------------------------------------------------
876+
for (auto const& trk : tracks) {
877+
// Accept reco track
878+
if (std::abs(trk.eta()) > 0.5f) {
879+
continue;
880+
}
881+
882+
// Map track's collision row id -> dense colls index
883+
const int collRowId = trk.collisionId();
884+
if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) {
885+
continue;
886+
}
887+
const int dIdx = collRowIdToDense[collRowId];
888+
if (dIdx < 0) {
889+
continue; // this track's collision is not in our 'colls' view
890+
}
891+
892+
// Get the collision row (dense index in colls view)
893+
auto col = colls.rawIteratorAt(dIdx);
894+
895+
// MC collision id (row index in aod::McCollisions)
896+
const int mcCollId = col.mcCollisionId();
897+
if (mcCollId < 0 || mcCollId >= (int)mcCollisions.size()) {
898+
continue;
899+
}
900+
mcReconstructed[mcCollId] = 1;
901+
902+
// MC particle id (row index in aod::McParticles)
903+
const int mcPid = trk.mcParticleId();
904+
if (mcPid < 0 || mcPid >= (int)mcParticles.size()) {
905+
continue;
906+
}
907+
908+
// MC multiplicity for that MC collision (cache)
909+
float mult = mcMultCache[mcCollId];
910+
if (mult < 0.f) {
911+
std::vector<float> tmp;
912+
mult = mcMult[mcCollId];
913+
mcMultCache[mcCollId] = mult;
914+
}
915+
916+
auto mcPar = mcParticles.rawIteratorAt(mcPid);
917+
918+
// Apply the same acceptance as in MC multiplicity definition
919+
if (!mcPar.isPhysicalPrimary()) {
920+
continue;
921+
}
922+
if (std::abs(mcPar.eta()) > 0.5f) {
923+
continue;
924+
}
925+
926+
int q = 0;
927+
if (auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode())) {
928+
q = int(std::round(pdgEntry->Charge() / 3.0));
929+
}
930+
if (q == 0) {
931+
continue;
932+
}
933+
934+
// Mark reconstructed MC particle (now that it truly passed & matched)
935+
isReco[mcPid] = 1;
936+
isRecoMult[mcPid] = mult;
937+
938+
const float multReco = col.multNTracksGlobal(); // or recoMultDense[dIdx]
939+
const float ptReco = trk.pt();
940+
const float ptMC = mcPar.pt();
941+
942+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco);
943+
}
944+
945+
// ------------------------------------------------------------
946+
// MC particles with no reco track (iterate by row index)
947+
// ------------------------------------------------------------
948+
for (int pid = 0; pid < (int)mcParticles.size(); ++pid) {
949+
if (!isReco[pid]) {
950+
auto mcp = mcParticles.rawIteratorAt(pid);
951+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[pid], mcp.pt());
952+
}
953+
}
954+
955+
// ------------------------------------------------------------
956+
// Unreconstructed MC collisions (iterate by row index)
957+
// ------------------------------------------------------------
958+
for (int mcid = 0; mcid < (int)mcCollisions.size(); ++mcid) {
959+
if (!mcReconstructed[mcid]) {
960+
std::vector<float> mcptvec;
961+
const int mult = mcMult[mcid];
962+
for (auto const& pt : mcptvec) {
963+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt);
900964
}
901965
}
902966
}
967+
}
968+
903969
PROCESS_SWITCH(NonPromptCascadeTask, processdNdetaMC, "process mc dN/deta", false);
904970
};
905971

0 commit comments

Comments
 (0)