diff --git a/PWGCF/Femto/FemtoNuclei/TableProducer/PiNucleiFemto.cxx b/PWGCF/Femto/FemtoNuclei/TableProducer/PiNucleiFemto.cxx index c7272ad3e30..7cb1864593e 100644 --- a/PWGCF/Femto/FemtoNuclei/TableProducer/PiNucleiFemto.cxx +++ b/PWGCF/Femto/FemtoNuclei/TableProducer/PiNucleiFemto.cxx @@ -63,6 +63,7 @@ #include #include #include +#include #include // std::prev #include #include @@ -204,12 +205,13 @@ struct PiNucleiFemto { PresliceUnsorted hypPerCol = o2::aod::hyperrec::collisionId; // binning for EM background - ConfigurableAxis axisVertex{"axisVertex", {30, -10, 10}, "Binning for multiplicity"}; + ConfigurableAxis axisVertex{"axisVertex", {30, -10, 10}, "Binning for vtxz"}; ConfigurableAxis axisCentrality{"axisCentrality", {40, 0, 100}, "Binning for centrality"}; using BinningType = ColumnBinningPolicy; BinningType binningPolicy{{axisVertex, axisCentrality}, true}; SliceCache cache; SameKindPair mPair{binningPolicy, settingNoMixedEvents, -1, &cache}; + // Pair hyperPair{binningPolicy, settingNoMixedEvents, -1, &cache}; std::array mBBparamsDe; @@ -229,6 +231,7 @@ struct PiNucleiFemto { "QA", {{"hVtxZ", "Vertex distribution in Z;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}}, {"hNcontributor", "Number of primary vertex contributor", {HistType::kTH1F, {{2000, 0.0f, 2000.0f}}}}, + {"hCentrality", "Centrality", {HistType::kTH1F, {{200, 0.0f, 100.0f}}}}, {"hTrackSel", "Accepted tracks", {HistType::kTH1F, {{Selections::kAll, -0.5, static_cast(Selections::kAll) - 0.5}}}}, {"hEvents", "; Events;", {HistType::kTH1F, {{3, -0.5, 2.5}}}}, {"hEmptyPool", "svPoolCreator did not find track pairs false/true", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, @@ -240,6 +243,9 @@ struct PiNucleiFemto { {"hNuPitInvMass", "; M(Nu + p) (GeV/#it{c}^{2})", {HistType::kTH1F, {{300, 3.74f, 4.34f}}}}, {"hNuPt", "#it{p}_{T} distribution; #it{p}_{T} (GeV/#it{c})", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, {"hPiPt", "Pt distribution; #it{p}_{T} (GeV/#it{c})", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}}, + {"hHe3TPCnsigma", "NsigmaHe3 TPC distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TPC}(He3)", {HistType::kTH2F, {{100, -2.0f, 2.0f}, {200, -5.0f, 5.0f}}}}, + {"hHe3P", "Pin distribution; p (GeV/#it{c})", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}}, + {"hHe3P_preselected", "Pin distribution_preselected; p (GeV/#it{c})", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}}, {"hNuEta", "eta distribution; #eta(Nu)", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}}, {"hPiEta", "eta distribution; #eta(#pi)", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}}, {"hNuPhi", "phi distribution; phi(Nu)", {HistType::kTH1F, {{600, -4.0f, 4.0f}}}}, @@ -268,6 +274,52 @@ struct PiNucleiFemto { false, true}; + int numOfCentBins = 40; + int numOfVertexZBins = 30; + float Vz_low = -10.0f; + float Vz_high = 10.0f; + float Vz_step = (Vz_high - Vz_low) / numOfVertexZBins; + + struct EventRef { + int64_t collisionId; + }; + + struct PoolBin { + std::deque events; + }; + + std::vector All_Event_pool; + bool isInitialized = false; + + int nPoolBins() const { return numOfVertexZBins * numOfCentBins; } + + void initializePools() + { + All_Event_pool.clear(); + All_Event_pool.resize(nPoolBins()); + isInitialized = true; + } + + int where_pool(float vz, float v0Centr) const + { + float CentBinWidth = 100.0 / numOfCentBins; // = 2.5 + + int iy = static_cast(std::floor(v0Centr / CentBinWidth)); + if (iy < 0) + iy = 0; + if (iy >= numOfCentBins) + iy = numOfCentBins - 1; + + int ix = static_cast(std::floor((vz - Vz_low) / Vz_step)); + if (ix < 0) + ix = 0; + if (ix >= numOfVertexZBins) + ix = numOfVertexZBins - 1; + + int bin = ix + numOfVertexZBins * iy; + return bin; + } + void init(o2::framework::InitContext&) { mZorroSummary.setObject(mZorro.getZorroSummary()); @@ -488,6 +540,41 @@ struct PiNucleiFemto { return false; } + float averageClusterSizeCosl(uint32_t itsClusterSizes, float eta) + { + float average = 0; + int nclusters = 0; + const float cosl = 1. / std::cosh(eta); + const int nlayerITS = 7; + + for (int layer = 0; layer < nlayerITS; layer++) { + if ((itsClusterSizes >> (layer * 4)) & 0xf) { + nclusters++; + average += (itsClusterSizes >> (layer * 4)) & 0xf; + } + } + if (nclusters == 0) { + return 0; + } + return average * cosl / nclusters; + }; + + bool selectionPIDHyper(const aod::DataHypCandsWColl::iterator& V0Hyper) + { + mQaRegistry.fill(HIST("hHe3P_preselected"), V0Hyper.tpcMomHe()); + float averClusSizeHe = averageClusterSizeCosl(V0Hyper.itsClusterSizesHe(), V0Hyper.etaHe3()); + if (averClusSizeHe <= 4) { + return false; + } + if (V0Hyper.tpcChi2He() <= 0.5) { + return false; + } + mQaRegistry.fill(HIST("hHe3P"), V0Hyper.tpcMomHe()); + mQaRegistry.fill(HIST("hHe3TPCnsigma"), V0Hyper.nSigmaHe()); + + return true; + } + // ================================================================================================================== template @@ -720,21 +807,23 @@ struct PiNucleiFemto { template void pairTracksSameEventHyper(const Ttrack& piTracks, const Thypers& V0Hypers) { - for (const auto& piTrack : piTracks) { - - mQaRegistry.fill(HIST("hTrackSel"), Selections::kNoCuts); - - if (!selectTrack(piTrack)) { + for (const auto& V0Hyper : V0Hypers) { + if (!selectionPIDHyper(V0Hyper)) { continue; } - mQaRegistry.fill(HIST("hTrackSel"), Selections::kTrackCuts); + for (const auto& piTrack : piTracks) { - if (!selectionPIDPion(piTrack)) { - continue; - } - mQaRegistry.fill(HIST("hTrackSel"), Selections::kPID); + mQaRegistry.fill(HIST("hTrackSel"), Selections::kNoCuts); - for (const auto& V0Hyper : V0Hypers) { + if (!selectTrack(piTrack)) { + continue; + } + mQaRegistry.fill(HIST("hTrackSel"), Selections::kTrackCuts); + + if (!selectionPIDPion(piTrack)) { + continue; + } + mQaRegistry.fill(HIST("hTrackSel"), Selections::kPID); SVCand pair; pair.tr0Idx = piTrack.globalIndex(); @@ -770,6 +859,29 @@ struct PiNucleiFemto { } } + template + void pairHyperEventMixing(T1& pionCands, T2& hypCands) + { + for (const auto& hypCand : hypCands) { + if (!selectionPIDHyper(hypCand)) { + continue; + } + for (const auto& pionCand : pionCands) { + if (!selectTrack(pionCand) || !selectionPIDPion(pionCand)) { + continue; + } + + SVCand pair; + pair.tr0Idx = hypCand.globalIndex(); + pair.tr1Idx = pionCand.globalIndex(); + const int collIdx = hypCand.collisionId(); + CollBracket collBracket{collIdx, collIdx}; + pair.collBracket = collBracket; + mTrackHypPairs.push_back(pair); + } + } + } + template void fillTable(const PiNucandidate& piNucand, const Tcoll& collision) { @@ -977,6 +1089,15 @@ struct PiNucleiFemto { if (!fillCandidateInfoHyper(v0hyper, piTrack, piNucand, isMixedEvent)) { continue; } + mQaRegistry.fill(HIST("hNuPt"), piNucand.recoPtNu()); + mQaRegistry.fill(HIST("hPiPt"), piNucand.recoPtPi()); + mQaRegistry.fill(HIST("hNuEta"), piNucand.recoEtaNu()); + mQaRegistry.fill(HIST("hPiEta"), piNucand.recoEtaPi()); + mQaRegistry.fill(HIST("hNuPhi"), piNucand.recoPhiNu()); + mQaRegistry.fill(HIST("hPiPhi"), piNucand.recoPhiPi()); + mQaRegistry.fill(HIST("hNuPitInvMass"), piNucand.invMass); + mQaRegistry.fill(HIST("hNClsPiITS"), piNucand.nClsItsPi); + mQaRegistry.fill(HIST("hisBkgEM"), piNucand.isBkgEM); auto collision = collisions.rawIteratorAt(piNucand.collisionID); @@ -1021,6 +1142,7 @@ struct PiNucleiFemto { { mGoodCollisions.clear(); mGoodCollisions.resize(collisions.size(), false); + // LOG(info) << "Number of hyperCandidates read = " << V0Hypers.size(); for (const auto& collision : collisions) { @@ -1068,6 +1190,76 @@ struct PiNucleiFemto { fillPairs(collisions, tracks, /*isMixedEvent*/ true); } PROCESS_SWITCH(PiNucleiFemto, processMixedEvent, "Process Mixed event", false); + + /*void processMixedEventHyper(const CollisionsFull& collisions, o2::aod::DataHypCandsWColl const& V0Hypers, const TrackCandidates& pitracks) + { + LOG(debug) << "Processing mixed event for hypertriton"; + mTrackHypPairs.clear(); + + for (const auto& [c1, tracks1, c2, V0Hypers2] : hyperPair) { + if (!c1.sel8() || !c2.sel8()) { + continue; + } + + mQaRegistry.fill(HIST("hNcontributor"), c2.numContrib()); + //mQaRegistry.fill(HIST("hCentrality"), c2.centFT0C()); + mQaRegistry.fill(HIST("hVtxZ"), c2.posZ()); + + pairHyperEventMixing(tracks1, V0Hypers2); + } +} +PROCESS_SWITCH(PiNucleiFemto, processMixedEventHyper, "Process Mixed event", false);*/ + + void processMixedEventHyperPool(const CollisionsFull& collisions, o2::aod::DataHypCandsWColl const& V0Hypers, const TrackCandidates& pitracks) + { + LOG(info) << "Processing mixed event for hypertriton"; + + mTrackHypPairs.clear(); + if (!isInitialized) { + initializePools(); + LOG(info) << "Initialized event pool with size = " << All_Event_pool.size(); + } + for (auto const& collision : collisions) { + int poolIndexPi = where_pool(collision.posZ(), collision.centFT0C()); + auto& pool = All_Event_pool[poolIndexPi]; + + if (poolIndexPi < 0 || static_cast(poolIndexPi) >= All_Event_pool.size()) { + continue; + } + + for (auto const& storedEvent : pool.events) { + auto c1 = collisions.iteratorAt(storedEvent.collisionId); + const auto& c2 = collision; + if (!c1.sel8() || !c2.sel8()) + continue; + + std::vector tracks1; + for (auto const& t : pitracks) { + if (t.collisionId() == c1.globalIndex()) { + tracks1.push_back(t); + } + } + + std::vector hypers2; + for (auto const& h : V0Hypers) { + if (h.collisionId() != c2.globalIndex()) + continue; + int poolIndexHyp = where_pool(h.zPrimVtx(), h.centralityFT0C()); + if (poolIndexHyp != poolIndexPi) + continue; + hypers2.push_back(h); + } + pairHyperEventMixing(tracks1, hypers2); + } + fillPairsHyper(collisions, pitracks, V0Hypers, /*isMixedEvent*/ true); + + if (static_cast(pool.events.size()) >= settingNoMixedEvents) { + pool.events.pop_front(); + } + pool.events.push_back({collision.globalIndex()}); + } + } + PROCESS_SWITCH(PiNucleiFemto, processMixedEventHyperPool, "Process Mixed event", false); }; WorkflowSpec defineDataProcessing(const ConfigContext& cfgc)