diff --git a/PWGDQ/Core/MixingHandler.cxx b/PWGDQ/Core/MixingHandler.cxx index 7d3a8e49f62..7bdaa62f072 100644 --- a/PWGDQ/Core/MixingHandler.cxx +++ b/PWGDQ/Core/MixingHandler.cxx @@ -62,6 +62,8 @@ void MixingHandler::AddMixingVariable(int var, const std::vector& binLims { fVariables[var] = fVariableLimits.size(); fVariableLimits.push_back(binLims); + // FillEvent() only fills variables marked as used + VarManager::SetUseVariable(var); } /* @@ -125,6 +127,8 @@ int MixingHandler::FindEventCategory(float* values) // loop over the variables and find out in which bin the value of the variable for the event is located std::vector bin; + // number of bins per variable in the iteration order of fVariables (fVariableLimits is in insertion order) + std::vector nBins; for (auto [var, pos] : fVariables) { // check that the value is within limits, if not return -1 to exclude the event from mixing size_t binValue = std::distance(fVariableLimits[pos].begin(), std::upper_bound(fVariableLimits[pos].begin(), fVariableLimits[pos].end(), values[var])); @@ -132,6 +136,7 @@ int MixingHandler::FindEventCategory(float* values) return -1; // all variables must be inside limits } bin.push_back(binValue - 1); + nBins.push_back(fVariableLimits[pos].size() - 1); } // Hash the bin values to define a unique category @@ -149,7 +154,7 @@ int MixingHandler::FindEventCategory(float* values) if (iv2 == iv1) { tempCategory *= bin[iv2]; } else { - tempCategory *= (fVariableLimits[iv2].size() - 1); + tempCategory *= nBins[iv2]; } } category += tempCategory; @@ -167,15 +172,25 @@ int MixingHandler::GetBinFromCategory(VarManager::Variables var, int category) c return -1; } - // Search for the position of the variable "var" in the internal variable list of the handler - int ivar = fVariables.at(var); + // number of bins and position of var in the iteration order of fVariables, as used by FindEventCategory() + std::vector nBins; + int ivar = -1; + for (auto const& [v, pos] : fVariables) { + if (v == var) { + ivar = static_cast(nBins.size()); + } + nBins.push_back(fVariableLimits[pos].size() - 1); + } + if (ivar < 0) { + return -1; + } // extract the bin position in variable "var" from the category int norm = 1; - for (int i = fVariables.size() - 1; i > ivar; --i) { - norm *= (fVariableLimits[i].size() - 1); + for (size_t i = nBins.size() - 1; i > static_cast(ivar); --i) { + norm *= nBins[i]; } int truncatedCategory = category - (category % norm); truncatedCategory /= norm; - return truncatedCategory % (fVariableLimits[ivar].size() - 1); + return truncatedCategory % nBins[ivar]; } diff --git a/PWGDQ/Core/MixingHandler.h b/PWGDQ/Core/MixingHandler.h index 931db50f8a2..839666f36ec 100644 --- a/PWGDQ/Core/MixingHandler.h +++ b/PWGDQ/Core/MixingHandler.h @@ -33,17 +33,31 @@ class MixingHandler : public TNamed { public: + // number of track cuts which fit in the 32-bit filtering masks + static constexpr int NMaxCuts = 32; + // Struct to define track properties relevant for mixing and few utility functions struct MixingTrack { float pt; float eta; float phi; uint32_t filteringFlags; + // globalIndex is unique only within a dataframe, so the dataframe sequence is part of the track identity + uint64_t dataFrameSequence = 0; + uint64_t trackGlobalIndex = 0; + // electric charge of the track; 0 means "not set" and disables the charge dependent pair variables + int8_t sign = 0; + bool IsSamePhysicalTrack(const MixingTrack& other) const + { + return dataFrameSequence == other.dataFrameSequence && trackGlobalIndex == other.trackGlobalIndex; + } // Clear a bit once the track was used in mixing for that bit for the required pool depth. void ClearBit(uint32_t mask) { filteringFlags &= ~mask; } void Print() const { - std::cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", filteringFlags: " << filteringFlags << std::endl; + std::cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", sign: " << static_cast(sign) + << ", filteringFlags: " << filteringFlags + << ", dataframe: " << dataFrameSequence << ", track: " << trackGlobalIndex << std::endl; } }; @@ -157,6 +171,17 @@ class MixingHandler : public TNamed CleanPool(); events.push_back(event); } + // Same, but the stored events are aged only for the cuts in agingMask. Passing the filtering mask of the + // incoming event ages an event only for the cuts for which a mixed pair was actually produced, so that the + // pool depth is a number of mixed partners and not a number of arrivals. + void UpdatePool(const MixingEvent& event, int16_t poolDepth, uint32_t agingMask) + { + for (auto& poolEvent : events) { + poolEvent.IncrementCounters(agingMask, poolDepth); + } + CleanPool(); + events.push_back(event); + } // getter for the events in the pool const std::vector& GetEvents() const { return events; } @@ -176,6 +201,8 @@ class MixingHandler : public TNamed // setters void AddMixingVariable(int var, const std::vector& binLims); void SetPoolDepth(int16_t depth) { fPoolDepth = depth; } + // remove all pools (e.g. at a run change) + void ClearPools() { fPools.clear(); } // getters // int GetNMixingVariables() const { return fVariables.size(); } diff --git a/PWGDQ/Tasks/tableReader.cxx b/PWGDQ/Tasks/tableReader.cxx index 8501e0792e5..ef3e6713b2f 100644 --- a/PWGDQ/Tasks/tableReader.cxx +++ b/PWGDQ/Tasks/tableReader.cxx @@ -273,6 +273,10 @@ struct AnalysisEventSelection { if (fMixHandler != nullptr) { int hh = fMixHandler->FindEventCategory(VarManager::fgValues); + // events outside the mixing limits (-1) get a distinct negative hash so that they are not mixed with each other + if (hh < 0) { + hh = -1 - static_cast(event.globalIndex()); + } hash(hh); } } diff --git a/PWGDQ/Tasks/tableReader_withAssoc.h b/PWGDQ/Tasks/tableReader_withAssoc.h index 541a8f039cb..b8028f2b773 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.h +++ b/PWGDQ/Tasks/tableReader_withAssoc.h @@ -494,6 +494,10 @@ struct AnalysisEventSelection { // create the mixing hash and publish it into the hash table if (fMixHandler != nullptr) { int hh = fMixHandler->FindEventCategory(dqtablereader_helpers::varValues()); + // events outside the mixing limits (-1) get a distinct negative hash so that they are not mixed with each other + if (hh < 0) { + hh = -1 - static_cast(event.globalIndex()); + } hash(hh); } } @@ -1427,6 +1431,8 @@ struct AnalysisSameEventPairing { HistogramManager* fHistMan = nullptr; MixingHandler fMixingHandler; + // dataframe counter, part of the track identity in the mixing pools + uint64_t fMixingDataFrameSequence = 0; o2::analysis::DQMlResponse fDQMlResponse; std::vector fOutputMlPsi2ee; // TODO: check this is needed or not @@ -1714,6 +1720,9 @@ struct AnalysisSameEventPairing { } if (fConfigRunMixingAcrossTFs) { + if (fNCutsBarrel > MixingHandler::NMaxCuts) { + LOGF(fatal, "Across-TF mixing supports at most %d barrel track-cut bits, got %d", MixingHandler::NMaxCuts, fNCutsBarrel); + } TString mixVarsString = fConfigMixingVariables.value; TString mixVarsJsonString = fConfigMixingVariablesJson.value; std::unique_ptr objArray(mixVarsString.Tokenize(",")); @@ -1899,6 +1908,10 @@ struct AnalysisSameEventPairing { { if (events.size() > 0) { // Additional protection to avoid crashing of events.begin().runNumber() if (fCurrentRun != events.begin().runNumber()) { + if (fConfigRunMixingAcrossTFs) { + // do not mix events from different runs + fMixingHandler.ClearPools(); + } initParamsFromCCDB(events.begin().timestamp(), events.begin().runNumber(), TTwoProngFitter); fCurrentRun = events.begin().runNumber(); } @@ -1981,6 +1994,10 @@ struct AnalysisSameEventPairing { // constexpr bool fillFlowReso = eventHasQvector || eventHasQvectorCentr; bool isSelectedBDT = false; fNPairPerEvent = 0; + uint64_t currentMixingDataFrameSequence = 0; + if (fConfigRunMixingAcrossTFs) { + currentMixingDataFrameSequence = ++fMixingDataFrameSequence; + } for (auto const& event : events) { if (!event.isEventSelected_bit(0)) { @@ -2006,6 +2023,10 @@ struct AnalysisSameEventPairing { } VarManager::FillEventFlowResoFactor(ResoFlowSP, ResoFlowEP); } + int mixingCategory = -1; + if (fConfigRunMixingAcrossTFs) { + mixingCategory = fMixingHandler.FindEventCategory(dqtablereader_helpers::varValues()); + } bool isFirst = true; for (auto const& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { @@ -2478,6 +2499,9 @@ struct AnalysisSameEventPairing { if (fConfigRunMixingAcrossTFs) { // run event mixing across TFs + if (mixingCategory < 0) { + continue; + } // 1) create a MixingEvent and fill it with the relevant tracks MixingHandler::MixingEvent mixingEvent; uint32_t trackFilterForMixing = 0; @@ -2488,7 +2512,7 @@ struct AnalysisSameEventPairing { continue; } auto t1 = assoc.template reducedtrack_as(); - MixingHandler::MixingTrack mixingTrack(t1.pt(), t1.eta(), t1.phi(), trackFilterForMixing); + MixingHandler::MixingTrack mixingTrack(t1.pt(), t1.eta(), t1.phi(), trackFilterForMixing, currentMixingDataFrameSequence, static_cast(assoc.reducedtrackId()), static_cast(t1.sign())); if (t1.sign() > 0) { mixingEvent.AddTrack1(mixingTrack); } else { @@ -2496,8 +2520,11 @@ struct AnalysisSameEventPairing { } } } + if (mixingEvent.tracks1.empty() && mixingEvent.tracks2.empty()) { + continue; + } // 2) run the mixing with the events in the pool corresponding to this event - auto& pool = fMixingHandler.GetPool(fMixingHandler.FindEventCategory(dqtablereader_helpers::varValues())); + auto& pool = fMixingHandler.GetPool(mixingCategory); for (auto const& poolEvent : pool.GetEvents()) { for (auto const& t1 : mixingEvent.tracks1) { // run +- pairing @@ -2516,9 +2543,9 @@ struct AnalysisSameEventPairing { } // run ++ pairing for (auto const& t2 : poolEvent.tracks1) { - // check the two-track filter for the mixed pair + // check the two-track filter for the mixed pair and skip the same track associated to both collisions uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags; - if (!mixedTwoTrackFilter) { + if (!mixedTwoTrackFilter || t1.IsSamePhysicalTrack(t2)) { continue; } VarManager::FillPairMEAcrossTFs(t1, t2); @@ -2546,9 +2573,9 @@ struct AnalysisSameEventPairing { } // run -- pairing for (auto const& t2 : poolEvent.tracks2) { - // check the two-track filter for the mixed pair + // check the two-track filter for the mixed pair and skip the same track associated to both collisions uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags; - if (!mixedTwoTrackFilter) { + if (!mixedTwoTrackFilter || t1.IsSamePhysicalTrack(t2)) { continue; } VarManager::FillPairMEAcrossTFs(t1, t2); @@ -2561,7 +2588,7 @@ struct AnalysisSameEventPairing { } } // 3) add the current event to the pool - pool.UpdatePool(mixingEvent, fMixingHandler.GetPoolDepth()); + pool.UpdatePool(mixingEvent, fMixingHandler.GetPoolDepth(), mixingEvent.filteringMask); // pool.Print(); } } // end loop over events diff --git a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx index 781930985c7..95fceb1e889 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx @@ -494,6 +494,16 @@ struct AnalysisEventSelection { VarManager::FillBC(bc); VarManager::FillEvent(event); + // the hash table is joined to the events by row order, so publish one row per event before any event selection + if (fMixHandler != nullptr) { + int hh = fMixHandler->FindEventCategory(VarManager::fgValues); + // events outside the mixing limits (-1) get a distinct negative hash so that they are not mixed with each other + if (hh < 0) { + hh = -1 - static_cast(event.globalIndex()); + } + hash(hh); + } + bool decision = false; if (fConfigQA) { fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); @@ -533,10 +543,6 @@ struct AnalysisEventSelection { auto& evIndices = fBCCollMap[bc.globalBC()]; evIndices.push_back(event.globalIndex()); } - if (fMixHandler != nullptr) { - int hh = fMixHandler->FindEventCategory(VarManager::fgValues); - hash(hh); - } } }