diff --git a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx index 47417fa070f..ec6aaffd8cf 100644 --- a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx +++ b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx @@ -37,7 +37,9 @@ #include "Math/Vector4D.h" #include "TDatabasePDG.h" +#include "THnSparse.h" #include "TParticlePDG.h" +#include "TTree.h" #include #include @@ -185,6 +187,7 @@ struct NonPromptCascadeTask { using CollisionCandidatesRun3 = soa::Join; using CollisionCandidatesRun3MC = soa::Join; using CollisionsWithLabel = soa::Join; + using TracksWithLabel = soa::Join; Preslice perCollision = aod::track::collisionId; Preslice perCollisionMC = aod::track::collisionId; @@ -212,8 +215,7 @@ struct NonPromptCascadeTask { Configurable cfgMaxMult{"cfgMaxMult", 8000.f, "Upper range of multiplicty histo"}; Configurable cfgMaxMultFV0{"cfgMaxMultFV0", 10000.f, "Upper range of multiplicty FV0 histo"}; - Configurable cfgMinMult{"cfgMinMult", 3000.f, "Lower range of FT0M histo in zoomed histo"}; - Configurable cfgMaxCent{"cfgMaxCent", 8.0025f, "Upper range of FT0M histo"}; + Configurable cfgPtEdgesdNdeta{"ptEdges", "0,0.2,0.4,0.6,0.8,1,1.2,1.6,2.0,2.4,2.8,3.2,3.6,4,4.5,5,5.5,6,7,8,10", "Pt bin edges (comma-separated)"}; Zorro mZorro; OutputObj mZorroSummary{"ZorroSummary"}; @@ -225,23 +227,11 @@ struct NonPromptCascadeTask { o2::vertexing::DCAFitterN<2> mDCAFitter; std::array mProcessCounter = {0, 0}; // {Tracked, All} std::map mToiMap; - std::unordered_map> mHistsPerRunMultVsCent; - std::unordered_map> mHistsPerRunMultVsCentZoom; - std::unordered_map> mHistsPerRunNtracktVsCent; - std::unordered_map> mHistsPerRunNtracktVsCentZoom; - - int nBinsMult = cfgMaxMult; - int nBinsMultFV0 = cfgMaxMultFV0; - int nBinsMultZoom = cfgMaxMult - cfgMinMult; - int nBinsCentZoom = (cfgMaxCent + 0.0025) / 0.005; - - AxisSpec multAxis = {nBinsMult, 0, cfgMaxMult, "Multiplicity FT0M"}; - AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FT0M"}; - AxisSpec centAxis = {101, -0.025, 101.025, "Centrality"}; - AxisSpec centAxisZoom = {nBinsCentZoom, -0.0025, cfgMaxCent, "Centrality"}; - AxisSpec multAxisZoom = {nBinsMultZoom, cfgMinMult, cfgMaxMult, "Multiplicity FT0M"}; - AxisSpec nTracksAxis = {100, 0., 100., "NTracksGlobal"}; - AxisSpec nTracksAxisMC = {100, 0., 100., "NTracksMC"}; + // + HistogramRegistry mRegistryMults{"Multhistos"}; + HistogramRegistry mRegistrydNdeta{"dNdetahistos"}; + + // void initCCDB(aod::BCsWithTimestamps::iterator const& bc) { @@ -262,7 +252,7 @@ struct NonPromptCascadeTask { } } - void init(InitContext const&) + void init(InitContext& context) { mZorroSummary.setObject(mZorro.getZorroSummary()); mCCDB->setURL(ccdbUrl); @@ -283,18 +273,77 @@ struct NonPromptCascadeTask { std::array cutsNames{"# candidates", "hasTOF", "nClusTPC", "nSigmaTPCbach", "nSigmaTPCprotontrack", "nSigmaTPCpiontrack", "cosPA"}; auto cutsOmega{std::get>(mRegistry.add("h_PIDcutsOmega", ";;Invariant mass (GeV/#it{c}^{2})", HistType::kTH2D, {{cutsNames.size(), -0.5, -0.5 + cutsNames.size()}, {125, 1.650, 1.700}}))}; auto cutsXi{std::get>(mRegistry.add("h_PIDcutsXi", ";;Invariant mass (GeV/#it{c}^{2})", HistType::kTH2D, {{6, -0.5, 5.5}, {125, 1.296, 1.346}}))}; - mRegistry.add("hMultVsCent", "hMultVsCent", HistType::kTH2F, {centAxis, multAxis}); - mRegistry.add("hMultVsCentZoom", "hMultVsCentZoom", HistType::kTH2F, {centAxisZoom, multAxisZoom}); - mRegistry.add("hNTracksVsCent", "hNTracksVsCent", HistType::kTH2F, {centAxis, nTracksAxis}); - mRegistry.add("hNTracksVsCentZoom", "hNTracksVsCentZoom", HistType::kTH2F, {centAxisZoom, nTracksAxis}); - mRegistry.add("hMultFV0VshNTracks", "hMultFV0VshNTracks", HistType::kTH2F, {nTracksAxis, multAxisFV0}); - mRegistry.add("hNTracksVsCentFV0A", "hNTracksVsCentFV0A", HistType::kTH2F, {nTracksAxis, centAxis}); - mRegistry.add("hNTracksMCVsTracksReco", "hNTracksMCVsTracksReco", HistType::kTH2F, {nTracksAxisMC, nTracksAxis}); - mRegistry.add("hNTracksMCNotInReco", "hNTracksMCNotInReco", HistType::kTH1F, {nTracksAxisMC}); for (size_t iBin{0}; iBin < cutsNames.size(); ++iBin) { cutsOmega->GetYaxis()->SetBinLabel(iBin + 1, cutsNames[iBin].c_str()); cutsXi->GetYaxis()->SetBinLabel(iBin + 1, cutsNames[iBin].c_str()); } + // + int nBinsMult = cfgMaxMult; + int nBinsMultFV0 = cfgMaxMultFV0; + + AxisSpec multAxis = {nBinsMult, 0, cfgMaxMult, "Multiplicity FT0M"}; + AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FT0M"}; + AxisSpec nTracksAxis = {100, 0., 100., "NTracksGlobal"}; + AxisSpec nTracksAxisMC = {100, 0., 100., "NTracksMC"}; + std::vector centBinning; + std::vector multBinning; + std::vector trackBinning; + std::vector runsBinning; + double cent = -0.0; + while (cent < 100) { + if (cent < 0.1) { + centBinning.push_back(cent); + cent += 0.01; + } else if (cent < 1) { + centBinning.push_back(cent); + cent += 0.1; + } else { + centBinning.push_back(cent); + cent += 1; + } + } + double ntracks = 0; + while (ntracks < 100) { + trackBinning.push_back(ntracks); + ntracks++; + } + double run = 550367 - 0.5; + for (int i = 0; i < 9000; i++) { + runsBinning.push_back(run); + run++; + } + AxisSpec centAxis{centBinning, "Centrality (%)"}; + // AxisSpec multAxis{multBinning, "Multiplicity"}; + AxisSpec trackAxisMC{trackBinning, "NTracks MC"}; + AxisSpec trackAxis{trackBinning, "NTracks Global Reco"}; + AxisSpec runsAxis{runsBinning, "Run Number"}; + + mRegistryMults.add("hCentMultsRuns", "hCentMultsRuns", HistType::kTHnSparseF, {centAxis, multAxis, centAxis, multAxisFV0, nTracksAxis, runsAxis}); + // + // dN/deta + // + bool runMCdNdeta = context.options().get("processdNdetaMC"); + // std::cout << "runMCdNdeta: " << runMCdNdeta << std::endl; + if (runMCdNdeta) { + std::vector ptBins; + std::vector tokens = o2::utils::Str::tokenize(cfgPtEdgesdNdeta, ','); + for (auto const& pts : tokens) { + double pt = 0; + try { + pt = std::stof(pts); + } catch (...) { + LOG(error) << "Wrong cfgPtEdgesdNdeta string:" << cfgPtEdgesdNdeta << std::endl; + } + ptBins.push_back(pt); + } + AxisSpec ptAxisMC{ptBins, "pT MC"}; + AxisSpec ptAxisReco{ptBins, "pT Reco"}; + + // multMeasured, multMC, ptMeasured, ptMC + mRegistrydNdeta.add("hdNdetaRM/hdNdetaRM", "hdNdetaRM", HistType::kTHnSparseF, {nTracksAxisMC, nTracksAxis, ptAxisMC, ptAxisReco}); + mRegistrydNdeta.add("hdNdetaRM/hdNdetaRMNotInRecoCol", "hdNdetaRMNotInRecoCol", HistType::kTHnSparseF, {nTracksAxisMC, ptAxisMC}); + mRegistrydNdeta.add("hdNdetaRM/hdNdetaRMNotInRecoTrk", "hdNdetaRMNotInRecoTrk", HistType::kTHnSparseF, {nTracksAxisMC, ptAxisMC}); + } } template @@ -361,27 +410,13 @@ struct NonPromptCascadeTask { { // std::cout << "Filling mult histos" << std::endl; for (const auto& coll : collisions) { - std::string histNameMvC = "mult/hMultVsCent_run" + std::to_string(mRunNumber); - std::string histNameMvCZ = "mult/hMultVsCentZoom_run" + std::to_string(mRunNumber); - std::string histNameTvC = "mult/hNTracksVsCent_run" + std::to_string(mRunNumber); - std::string histNameTvCZ = "mult/hNTracksVsCentZoom_run" + std::to_string(mRunNumber); - if (!mHistsPerRunMultVsCent.contains(histNameMvC)) { - mHistsPerRunMultVsCent[histNameMvC] = std::get>(mRegistry.add(histNameMvC.c_str(), histNameMvC.c_str(), HistType::kTH2F, {centAxis, multAxis})); - mHistsPerRunMultVsCentZoom[histNameMvCZ] = std::get>(mRegistry.add(histNameMvCZ.c_str(), histNameMvCZ.c_str(), HistType::kTH2F, {centAxisZoom, multAxisZoom})); - mHistsPerRunNtracktVsCent[histNameTvC] = std::get>(mRegistry.add(histNameTvC.c_str(), histNameTvC.c_str(), HistType::kTH2F, {centAxis, nTracksAxis})); - mHistsPerRunNtracktVsCentZoom[histNameTvCZ] = std::get>(mRegistry.add(histNameTvCZ.c_str(), histNameTvCZ.c_str(), HistType::kTH2F, {centAxisZoom, nTracksAxis})); - } - mHistsPerRunMultVsCent[histNameMvC]->Fill(coll.centFT0M(), coll.multFT0M()); - mHistsPerRunMultVsCentZoom[histNameMvCZ]->Fill(coll.centFT0M(), coll.multFT0M()); - mHistsPerRunNtracktVsCent[histNameTvC]->Fill(coll.centFT0M(), coll.multNTracksGlobal()); - mHistsPerRunNtracktVsCentZoom[histNameTvCZ]->Fill(coll.centFT0M(), coll.multNTracksGlobal()); - // run integrated histos - mRegistry.fill(HIST("hMultVsCent"), coll.centFT0M(), coll.multFT0M()); - mRegistry.fill(HIST("hMultVsCentZoom"), coll.centFT0M(), coll.multFT0M()); - mRegistry.fill(HIST("hNTracksVsCent"), coll.centFT0M(), (float)coll.multNTracksGlobal()); - mRegistry.fill(HIST("hNTracksVsCentZoom"), coll.centFT0M(), coll.multNTracksGlobal()); - mRegistry.fill(HIST("hMultFV0VshNTracks"), coll.multNTracksGlobal(), coll.multFV0A()); - mRegistry.fill(HIST("hNTracksVsCentFV0A"), coll.multNTracksGlobal(), coll.centFV0A()); + float centFT0M = coll.centFT0M(); + float multFT0M = coll.multFT0M(); + float centFV0A = coll.centFV0A(); + float multFV0A = coll.multFV0A(); + float multNTracks = coll.multNTracksGlobal(); + float run = mRunNumber; + mRegistryMults.fill(HIST("hCentMultsRuns"), centFT0M, multFT0M, centFV0A, multFV0A, multNTracks, run); } }; @@ -758,7 +793,7 @@ struct NonPromptCascadeTask { } PROCESS_SWITCH(NonPromptCascadeTask, processCascadesData, "process cascades: Data analysis", false); - int getMCMult(aod::McParticles const& mcParticles, int mcCollId) + int getMCMult(aod::McParticles const& mcParticles, int mcCollId, std::vector& ptMCvec) { int mult = 0; for (auto const& mcp : mcParticles) { @@ -776,35 +811,97 @@ struct NonPromptCascadeTask { accept = accept && (q != 0); if (accept) { ++mult; + ptMCvec.push_back(mcp.pt()); } } } return mult; } - void processNegMC(CollisionsWithLabel const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles) + void processdNdetaMC(CollisionsWithLabel const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TracksWithLabel const& tracks) { // std::cout << "ProcNegMC" << std::endl; + // Map: collision index -> reco multiplicity + std::vector recoMult(colls.size(), 0); + for (auto const& trk : tracks) { + int collId = trk.collisionId(); + // Here you can impose same track cuts as for MC (eta, primaries, etc.) + if (std::abs(trk.eta()) > 0.5f) { + continue; + } + ++recoMult[collId]; + } + std::vector isReco(mcParticles.size(), 0); + std::vector isRecoMult(mcParticles.size(), 0); std::vector mcReconstructed(mcCollisions.size(), 0); + // std::cout << " reco cols with mc:" << colls.size() << " tracks:" << tracks.size() << " mccols:" << mcCollisions.size() << "mcParticles:" << mcParticles.size() << std::endl; for (auto const& col : colls) { int mcCollId = col.mcCollisionId(); // col.template mcCollision_as(); + int collId = col.globalIndex(); // auto mc = col.mcCollision(); // int mcId = mc.globalIndex(); // std::cout << "globalIndex:" << mcId << " colID:" << mcCollId << std::endl; - int mult = getMCMult(mcParticles, mcCollId); + std::vector mcptvec; + float mult = getMCMult(mcParticles, mcCollId, mcptvec); mcReconstructed[mcCollId] = 1; - mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal()); + for (auto const& trk : tracks) { + int mcPid = trk.mcParticleId(); // depends on your label table + if (mcPid >= 0 && mcPid < (int)mcParticles.size()) { + isReco[mcPid] = 1; + isRecoMult[mcPid] = mult; + } else { + continue; + } + if (trk.collisionId() != collId) { + continue; // different event + } + auto mcPar = mcParticles.rawIteratorAt(mcPid); + // Apply same acceptance as in MC multiplicity + if (!mcPar.isPhysicalPrimary()) { + continue; + } + if (std::abs(mcPar.eta()) > 0.5f) { + continue; + } + int q = 0; + auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode()); + if (pdgEntry) { + q = int(std::round(pdgEntry->Charge() / 3.0)); + } + if (q == 0) { + continue; + } + // float multReco = recoMult[collId]; + float multReco = col.multNTracksGlobal(); + float ptReco = trk.pt(); + float ptMC = mcPar.pt(); + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco); + } + + // mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal()); } + // count mc particles with no reco tracks + for (auto const& mcp : mcParticles) { + int mcPidG = mcp.globalIndex(); + // std::cout << "mcPidG:" << mcPidG << std::endl; + if (!isReco[mcPidG]) { + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[mcPidG], mcp.pt()); + } + } + // count unreconstructed mc collisions for (auto const& mc : mcCollisions) { int gindex = mc.globalIndex(); // std::cout << "mc globalIndex:" << gindex << std::endl; if (!mcReconstructed[gindex]) { - int mult = getMCMult(mcParticles, gindex); + std::vector mcptvec; + int mult = getMCMult(mcParticles, gindex, mcptvec); // std::cout << "===> unreconstructed:" << mult << std::endl; - mRegistry.fill(HIST("hNTracksMCNotInReco"), mult); + for (auto const& pt : mcptvec) { + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt); + } } } } - PROCESS_SWITCH(NonPromptCascadeTask, processNegMC, "process mc", false); + PROCESS_SWITCH(NonPromptCascadeTask, processdNdetaMC, "process mc dN/deta", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)