Skip to content

Commit fe37353

Browse files
lietavalietava
authored
[PWGLF] Omegahm (#14279)
1 parent 286e851 commit fe37353

File tree

1 file changed

+160
-90
lines changed

1 file changed

+160
-90
lines changed

‎PWGLF/Tasks/Strangeness/nonPromptCascade.cxx‎

Lines changed: 160 additions & 90 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,184 @@ 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)
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)
796804
{
797-
int mult = 0;
805+
//-------------------------------------------------------------
806+
// MC mult for all MC coll
807+
//--------------------------------------------------------------
808+
std::vector<int> mcMult(mcCollisions.size(), 0);
798809
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-
}
810+
int mcid = mcp.mcCollisionId();
811+
if (mcid < 0 || mcid >= (int)mcMult.size())
812+
continue;
813+
814+
// apply your primary/eta/charge definition here
815+
if (!mcp.isPhysicalPrimary())
816+
continue;
817+
if (std::abs(mcp.eta()) > 0.5f)
818+
continue;
819+
int q = 0;
820+
if (auto pdg = TDatabasePDG::Instance()->GetParticle(mcp.pdgCode())) {
821+
q = int(std::round(pdg->Charge() / 3.0));
815822
}
823+
if (q == 0)
824+
continue;
825+
826+
++mcMult[mcid];
816827
}
817-
return mult;
818-
}
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);
828+
829+
// ------------------------------------------------------------
830+
// Build mapping: (aod::Collisions row id used by tracks.collisionId())
831+
// -> dense index in 'colls' (0..colls.size()-1)
832+
// We assume col.globalIndex() refers to the original aod::Collisions row.
833+
// ------------------------------------------------------------
834+
int maxCollRowId = -1;
835+
for (auto const& trk : tracks) {
836+
maxCollRowId = std::max(maxCollRowId, (int)trk.collisionId());
837+
}
838+
std::vector<int> collRowIdToDense(maxCollRowId + 1, -1);
839+
840+
int dense = 0;
841+
for (auto const& col : colls) {
842+
const int collRowId = col.globalIndex(); // row id in aod::Collisions
843+
if (collRowId >= 0 && collRowId < (int)collRowIdToDense.size()) {
844+
collRowIdToDense[collRowId] = dense;
845+
}
846+
++dense;
847+
}
848+
849+
// ------------------------------------------------------------
850+
// Reco multiplicity per *dense collision index in colls*
851+
// ------------------------------------------------------------
852+
std::vector<int> recoMultDense(colls.size(), 0);
824853
for (auto const& trk : tracks) {
825-
int collId = trk.collisionId();
826-
// Here you can impose same track cuts as for MC (eta, primaries, etc.)
827854
if (std::abs(trk.eta()) > 0.5f) {
828855
continue;
829856
}
830-
++recoMult[collId];
857+
const int collRowId = trk.collisionId();
858+
if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) {
859+
continue;
860+
}
861+
const int dIdx = collRowIdToDense[collRowId];
862+
if (dIdx >= 0) {
863+
++recoMultDense[dIdx];
864+
}
831865
}
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);
866+
867+
// ------------------------------------------------------------
868+
// MC bookkeeping: index by ROW INDEX (0..size-1), not globalIndex()
869+
// ------------------------------------------------------------
870+
std::vector<char> isReco(mcParticles.size(), 0);
871+
std::vector<float> isRecoMult(mcParticles.size(), 0.f);
872+
std::vector<char> mcReconstructed(mcCollisions.size(), 0);
873+
874+
// Optional cache of MC multiplicity per MC collision
875+
std::vector<float> mcMultCache(mcCollisions.size(), -1.f);
876+
877+
// ------------------------------------------------------------
878+
// Single pass over tracks: fill RM for tracks whose collision is in colls
879+
// ------------------------------------------------------------
880+
for (auto const& trk : tracks) {
881+
// Accept reco track
882+
if (std::abs(trk.eta()) > 0.5f) {
883+
continue;
884+
}
885+
886+
// Map track's collision row id -> dense colls index
887+
const int collRowId = trk.collisionId();
888+
if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) {
889+
continue;
890+
}
891+
const int dIdx = collRowIdToDense[collRowId];
892+
if (dIdx < 0) {
893+
continue; // this track's collision is not in our 'colls' view
894+
}
895+
896+
// Get the collision row (dense index in colls view)
897+
auto col = colls.rawIteratorAt(dIdx);
898+
899+
// MC collision id (row index in aod::McCollisions)
900+
const int mcCollId = col.mcCollisionId();
901+
if (mcCollId < 0 || mcCollId >= (int)mcCollisions.size()) {
902+
continue;
903+
}
844904
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);
905+
906+
// MC particle id (row index in aod::McParticles)
907+
const int mcPid = trk.mcParticleId();
908+
if (mcPid < 0 || mcPid >= (int)mcParticles.size()) {
909+
continue;
877910
}
878911

879-
// mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal());
912+
// MC multiplicity for that MC collision (cache)
913+
float mult = mcMultCache[mcCollId];
914+
if (mult < 0.f) {
915+
std::vector<float> tmp;
916+
mult = mcMult[mcCollId];
917+
mcMultCache[mcCollId] = mult;
918+
}
919+
920+
auto mcPar = mcParticles.rawIteratorAt(mcPid);
921+
922+
// Apply the same acceptance as in MC multiplicity definition
923+
if (!mcPar.isPhysicalPrimary()) {
924+
continue;
925+
}
926+
if (std::abs(mcPar.eta()) > 0.5f) {
927+
continue;
928+
}
929+
930+
int q = 0;
931+
if (auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode())) {
932+
q = int(std::round(pdgEntry->Charge() / 3.0));
933+
}
934+
if (q == 0) {
935+
continue;
936+
}
937+
938+
// Mark reconstructed MC particle (now that it truly passed & matched)
939+
isReco[mcPid] = 1;
940+
isRecoMult[mcPid] = mult;
941+
942+
const float multReco = col.multNTracksGlobal(); // or recoMultDense[dIdx]
943+
const float ptReco = trk.pt();
944+
const float ptMC = mcPar.pt();
945+
946+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco);
880947
}
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());
948+
949+
// ------------------------------------------------------------
950+
// MC particles with no reco track (iterate by row index)
951+
// ------------------------------------------------------------
952+
for (int pid = 0; pid < (int)mcParticles.size(); ++pid) {
953+
if (!isReco[pid]) {
954+
auto mcp = mcParticles.rawIteratorAt(pid);
955+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[pid], mcp.pt());
887956
}
888957
}
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]) {
958+
959+
// ------------------------------------------------------------
960+
// Unreconstructed MC collisions (iterate by row index)
961+
// ------------------------------------------------------------
962+
for (int mcid = 0; mcid < (int)mcCollisions.size(); ++mcid) {
963+
if (!mcReconstructed[mcid]) {
894964
std::vector<float> mcptvec;
895-
int mult = getMCMult(mcParticles, gindex, mcptvec);
896-
// std::cout << "===> unreconstructed:" << mult << std::endl;
965+
const int mult = mcMult[mcid];
897966
for (auto const& pt : mcptvec) {
898967
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt);
899968
}
900969
}
901970
}
902971
}
972+
903973
PROCESS_SWITCH(NonPromptCascadeTask, processdNdetaMC, "process mc dN/deta", false);
904974
};
905975

0 commit comments

Comments
 (0)