diff --git a/PWGCF/Femto/Core/baseSelection.h b/PWGCF/Femto/Core/baseSelection.h index ceb2184a96a..fcce0c11c13 100644 --- a/PWGCF/Femto/Core/baseSelection.h +++ b/PWGCF/Femto/Core/baseSelection.h @@ -40,9 +40,9 @@ namespace o2::analysis::femto /// It evaluates which selections are fulfilled, assembles a final bitmask, and tracks required vs. optional cuts. /// /// \tparam T Type of observable values (mostly floats). -/// \tparam BitmaskType Type used for internal bitmask operations (e.g., uint32_t, uint64_t). +/// \tparam BitmaskType Integer type used for bitmask operations (e.g., uint32_t, uint64_t). /// \tparam NumObservables Total number of observables handled. -template +template class BaseSelection { public: @@ -53,57 +53,41 @@ class BaseSelection virtual ~BaseSelection() = default; /// \brief Add a static-value based selection for a specific observable. - /// \param selectionValues Vector of threshold values. /// \param observableIndex Index of the observable. + /// \param selectionName Name of the selection. + /// \param selectionValues Vector of threshold values. /// \param limitType Type of limit (from limits::LimitType). - /// \param skipMostPermissiveBit Whether to skip the loosest threshold in the bitmask. - /// \param isMinimalCut Whether this cut is mandatory or optional. + /// \param skipMostPermissiveBit Whether to skip the loosest threshold when assembling the bitmask. + /// \param isMinimalCut Whether this cut is mandatory (must be passed for the candidate to be accepted). + /// \param isOptionalCut Whether this cut is optional (candidate is accepted if any optional cut passes). void addSelection(int observableIndex, std::string const& selectionName, std::vector const& selectionValues, limits::LimitType limitType, bool skipMostPermissiveBit, bool isMinimalCut, - bool isOptionCut) + bool isOptionalCut) { // check index - if (static_cast(observableIndex) >= NumObservables) { + if (static_cast(observableIndex) >= NumObservables) { LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables; } // init selection container for selection at given index - mSelectionContainers.at(observableIndex) = SelectionContainer(selectionName, selectionValues, limitType, skipMostPermissiveBit, isMinimalCut, isOptionCut); - - // check if any selections are configured - if (mSelectionContainers.at(observableIndex).isEmpty()) { - return; - } - - // keep track of selections and bits - mNSelectionBits += mSelectionContainers.at(observableIndex).getShift(); - mNSelection += mSelectionContainers.at(observableIndex).getNSelections(); + mSelectionContainers.at(observableIndex) = SelectionContainer(selectionName, selectionValues, limitType, skipMostPermissiveBit, isMinimalCut, isOptionalCut); - if (mNSelectionBits > sizeof(BitmaskType) * CHAR_BIT) { - LOG(fatal) << "Too many selections. At most " << sizeof(BitmaskType) * CHAR_BIT << " number of bits are supported"; - } - // check if any selection is minimal - if (mSelectionContainers.at(observableIndex).isMinimalCut()) { - mHasMinimalSelection = true; - } - // check if selection is optional - if (mSelectionContainers.at(observableIndex).isOptionalCut()) { - mHasOptionalSelection = true; - } + init(observableIndex); } /// \brief Add a function-based selection for a specific observable. - /// \param baseName Base name for TF1 functions. - /// \param lowerLimit Lower bound for the TF1 domain. - /// \param upperLimit Upper bound for the TF1 domain. - /// \param selectionValues Function definitions as strings. /// \param observableIndex Index of the observable. - /// \param limitType Type of limit. - /// \param skipMostPermissiveBit Whether to skip the loosest threshold in the bitmask. - /// \param isMinimalCut Whether this cut is mandatory or optional. + /// \param selectionName Name of the selection. + /// \param lowerLimit Lower bound of the TF1 domain. + /// \param upperLimit Upper bound of the TF1 domain. + /// \param functions Selection threshold functions as strings (parsed as TF1 expressions). + /// \param limitType Type of limit (from limits::LimitType). + /// \param skipMostPermissiveBit Whether to skip the loosest threshold when assembling the bitmask. + /// \param isMinimalCut Whether this cut is mandatory (must be passed for the candidate to be accepted). + /// \param isOptionalCut Whether this cut is optional (candidate is accepted if any optional cut passes). void addSelection(int observableIndex, std::string const& selectionName, T lowerLimit, @@ -114,73 +98,75 @@ class BaseSelection bool isMinimalCut, bool isOptionalCut) { - if (static_cast(observableIndex) >= NumObservables) { + if (static_cast(observableIndex) >= NumObservables) { LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables; } mSelectionContainers.at(observableIndex) = SelectionContainer(selectionName, lowerLimit, upperLimit, functions, limitType, skipMostPermissiveBit, isMinimalCut, isOptionalCut); - // check if any selections are configured - if (mSelectionContainers.at(observableIndex).isEmpty()) { - return; - } - - // advance mNSelections so we can use it as offset for next selection - mNSelectionBits += mSelectionContainers.at(observableIndex).getShift(); - mNSelection += mSelectionContainers.at(observableIndex).getNSelections(); + init(observableIndex); + } - if (mNSelectionBits > sizeof(BitmaskType) * CHAR_BIT) { - LOG(fatal) << "Too many selections. At most " << sizeof(BitmaskType) * CHAR_BIT << " are supported"; - } - // keep track of selection selections - // check if any cut is minimal - if (mSelectionContainers.at(observableIndex).isMinimalCut()) { - mHasMinimalSelection = true; - } - // check if any selection is optional - if (mSelectionContainers.at(observableIndex).isOptionalCut()) { - mHasOptionalSelection = true; + /// \brief Add a static-value based selection for a specific observable. + /// \param observableIndex Index of the observable. + /// \param selectionName Name of the selection. + /// \param selectionValues Vector of threshold values. + /// \param limitType Type of limit (from limits::LimitType). + /// \param skipMostPermissiveBit Whether to skip the loosest threshold when assembling the bitmask. + /// \param isMinimalCut Whether this cut is mandatory (must be passed for the candidate to be accepted). + /// \param isOptionalCut Whether this cut is optional (candidate is accepted if any optional cut passes). + void addSelection(int observableIndex, + std::string const& selectionName, + std::vector const& selectionRanges, + bool skipMostPermissiveBit, + bool isMinimalCut, + bool isOptionalCut) + { + // check index + if (static_cast(observableIndex) >= NumObservables) { + LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables; } + // init selection container for selection at given index + mSelectionContainers.at(observableIndex) = SelectionContainer(selectionName, selectionRanges, skipMostPermissiveBit, isMinimalCut, isOptionalCut); + + init(observableIndex); } - /// \brief Add a boolean based selection for a specific observable. - /// \param mode Whether the selection is not applied, minimal or optional cut + /// \brief Add a boolean-based selection for a specific observable. /// \param observableIndex Index of the observable. + /// \param selectionName Name of the selection. + /// \param mode Controls how the cut is applied: + /// -1 = optional cut, bit is stored in bitmask; + /// 0 = cut is disabled, no bit stored; + /// 1 = minimal (mandatory) cut, no extra bit stored since only one threshold exists. void addSelection(int observableIndex, std::string const& selectionName, int mode) { switch (mode) { - case -1: // cut is optional and we store bit for the cut + case -1: // cut is optional and we store a bit for it mSelectionContainers.at(observableIndex) = SelectionContainer(selectionName, std::vector{1}, limits::LimitType::kEqual, false, false, true); - mHasOptionalSelection = true; - mNSelectionBits += 1; - mNSelection += 1; break; - case 0: // cut is not applied, initalize with empty vector, so we bail out later + case 0: // cut is disabled; initialize with empty vector so evaluation bails out early mSelectionContainers.at(observableIndex) = SelectionContainer(selectionName, std::vector{}, limits::LimitType::kEqual, false, false, false); break; - case 1: // cut is added as mininal selection (since it is only one value, no extra bit is stored) + case 1: // mandatory cut; only one threshold so the most permissive bit is skipped and no extra bit is stored mSelectionContainers.at(observableIndex) = SelectionContainer(selectionName, std::vector{1}, limits::LimitType::kEqual, true, true, false); - mHasMinimalSelection = true; - mNSelection += 1; break; default: LOG(fatal) << "Invalid switch for boolean selection"; } - if (mNSelectionBits > sizeof(BitmaskType) * CHAR_BIT) { - LOG(fatal) << "Too many selections. At most " << sizeof(BitmaskType) * CHAR_BIT << " are supported"; - } + init(observableIndex); } /// \brief Update the limits of a function-based selection for a specific observable. /// \param observable Index of the observable. - /// \param value Value at which to evaluate the selection functions. + /// \param value Value at which to re-evaluate the selection functions. void updateLimits(int observable, T value) { mSelectionContainers.at(observable).updateLimits(value); } - /// \brief Reset the internal bitmask and evaluation flags before evaluating a new event. + /// \brief Reset the internal bitmask and evaluation flags before processing a new candidate. void reset() { mFinalBitmask.reset(); @@ -195,6 +181,8 @@ class BaseSelection } } + /// \brief Reset the selection container for a single observable. + /// \param observableIndex Index of the observable to reset. void reset(int observableIndex) { mSelectionContainers.at(observableIndex).reset(); } /// \brief Evaluate a single observable against its configured selections. @@ -206,22 +194,20 @@ class BaseSelection if (mSelectionContainers.at(observableIndex).isEmpty()) { return; } - // if any previous observable did not pass minimal selections, there is no point in setting bitmask for other observables - // minimal selection for each observable is computed after adding it + // if a previous observable already failed a minimal selection, + // there is no point in evaluating further observables if (!mPassesMinimalSelections) { return; } // set bitmask for given observable mSelectionContainers.at(observableIndex).evaluate(value); - // check if minimal selction for this observable holds - // if one minimal selection is not fullfilled, the condition failes + // if any minimal selection is not fulfilled, the candidate is rejected if (mHasMinimalSelection) { if (!mSelectionContainers.at(observableIndex).passesAsMinimalCut()) { mPassesMinimalSelections = false; } } - // check if any optional selection holds - // if one optional selection is fullfilled, the condition succeeds + // if any optional selection is fulfilled, the candidate is accepted if (mHasOptionalSelection) { if (mSelectionContainers.at(observableIndex).passesAsOptionalCut()) { mPassesOptionalSelections = true; @@ -231,39 +217,41 @@ class BaseSelection /// \brief Evaluate a single observable against its configured selections. /// \param observableIndex Index of the observable. - /// \param values vector of values of the observable. + /// \param values Vector of values of the observable. void evaluateObservable(int observableIndex, std::vector values) { // if there are no selections configured, bail out if (mSelectionContainers.at(observableIndex).isEmpty()) { return; } - // if any previous observable did not pass minimal selections, there is no point in setting bitmask for other observables - // minimal selection for each observable is computed after adding it + // if a previous observable already failed a minimal selection, + // there is no point in evaluating further observables if (!mPassesMinimalSelections) { return; } // set bitmask for given observable mSelectionContainers.at(observableIndex).evaluate(values); - // check if minimal selction for this observable holds + // if any minimal selection is not fulfilled, the candidate is rejected if (mHasMinimalSelection) { - if (mSelectionContainers.at(observableIndex).passesAsMinimalCut() == false) { + if (!mSelectionContainers.at(observableIndex).passesAsMinimalCut()) { mPassesMinimalSelections = false; } } - // check if any optional selection holds + // if any optional selection is fulfilled, the candidate is accepted if (mHasOptionalSelection) { - if (mSelectionContainers.at(observableIndex).passesAsOptionalCut() == true) { + if (mSelectionContainers.at(observableIndex).passesAsOptionalCut()) { mPassesOptionalSelections = true; } } } - /// \brief Add comments to specific observabel + /// \brief Add comments to the selections of a specific observable. + /// \param observableIndex Index of the observable. + /// \param comments Vector of comment strings, one per selection threshold. void addComments(int observableIndex, std::vector const& comments) { mSelectionContainers.at(observableIndex).addComments(comments); } - /// \brief Check if all required (minimal) and optional cuts are passed. - /// \return True if all required and at least one optional cut (if present) is passed. + /// \brief Check whether all required and optional cuts are passed. + /// \return True if all minimal cuts pass and, if optional cuts are present, at least one of them passes. bool passesAllRequiredSelections() const { if (mHasMinimalSelection && !mHasOptionalSelection) { @@ -275,23 +263,25 @@ class BaseSelection if (mHasMinimalSelection && mHasOptionalSelection) { return mPassesMinimalSelections && mPassesOptionalSelections; } + // if there are no minimal or optional selections, we pass let it pass with true return true; } - /// \brief Check if the optional selection for a specific observable is passed. + /// \brief Check whether the optional selection for a specific observable is passed. /// \param observableIndex Index of the observable. - /// \return True if at least one optional selection is fulfilled. + /// \return True if at least one optional selection for this observable is fulfilled. bool passesOptionalSelection(int observableIndex) const { return mSelectionContainers.at(observableIndex).passesAsOptionalCut(); } - /// \brief Assemble the global selection bitmask from individual observable selections. + /// \brief Assemble the global selection bitmask from all individual observable selections. + /// \tparam HistName Name of the histogram used to track selection statistics. template void assembleBitmask() { mHistRegistry->fill(HIST(HistName), mNSelection); - // if the required selections are not passed, we can break early + // if the required selections are not passed, clear the bitmask and return early if (!this->passesAllRequiredSelections()) { mFinalBitmask.reset(); return; @@ -299,17 +289,16 @@ class BaseSelection mHistRegistry->fill(HIST(HistName), mNSelection + 1); int binCenter = 0; - // to assemble bitmask, convert all bitmask into integers - // shift the current one and add the new bits + // assemble the final bitmask by shifting each container's bitmask to its offset and OR-ing it in for (auto const& selectionContainer : mSelectionContainers) { - // if there are no selections for a certain observable, skip + // skip observables with no configured selections if (selectionContainer.isEmpty()) { continue; } - // Shift the result to its offset and add the new values + // shift the container's bitmask to its offset and merge mFinalBitmask |= (selectionContainer.getBitmask() << selectionContainer.getOffset()); - for (int j = 0; j < selectionContainer.getNSelections(); ++j) { + for (std::size_t j = 0; j < selectionContainer.getNSelections(); ++j) { if (j == 0 && selectionContainer.isMinimalCut()) { // minimal cuts are always filled mHistRegistry->fill(HIST(HistName), binCenter); @@ -325,28 +314,36 @@ class BaseSelection } /// \brief Retrieve the assembled bitmask as an integer value. - /// \return The combined selection bitmask. + /// \return The combined selection bitmask for all observables. BitmaskType getBitmask() const { return static_cast(mFinalBitmask.to_ullong()); } - /// \brief Retrieve the assembled bitmask as an integer value. - /// \return The combined selection bitmask. + /// \brief Retrieve the bitmask for a single observable as an integer value. + /// \param observableIndex Index of the observable. + /// \return The selection bitmask for the specified observable. BitmaskType getBitmask(int observableIndex) const { return static_cast(mSelectionContainers.at(observableIndex).getBitmask().to_ullong()); } - /// \brief Set the assembled bitmask for on observable - /// \return The combined selection bitmask. + /// \brief Manually set the bitmask for a specific observable. + /// \tparam R Integer type of the bitmask value. + /// \param observableIndex Index of the observable. + /// \param bitmask Bitmask value to set. template void setBitmask(int observableIndex, R bitmask) { mSelectionContainers.at(observableIndex).setBitmask(bitmask); } + /// \brief Retrieve the loosest (most permissive) selection threshold for a specific observable. + /// \param observableIndex Index of the observable. + /// \return The loosest threshold value configured for this observable. T getLoosestSelection(int observableIndex) const { return mSelectionContainers.at(observableIndex).getLoosestSelection(); } + /// \brief Print the full configuration of all selections to the log. + /// \param objectName Name of the object owning this selection (used as label in the log output). void printSelections(const std::string& objectName) const { LOG(info) << "Printing Configuration of " << objectName; for (size_t idx = 0; idx < mSelectionContainers.size(); ++idx) { - const auto& container = mSelectionContainers[idx]; + const auto& container = mSelectionContainers.at(idx); if (container.isEmpty()) { continue; } @@ -360,15 +357,11 @@ class BaseSelection LOG(info) << " Bitmask shift : " << container.getShift(); LOG(info) << " Selections:"; - const bool useFunctions = container.isUsingFunctions(); - const auto& values = container.getSelectionValues(); - const auto& functions = container.getSelectionFunction(); - const auto& comments = container.getComments(); - - for (int j = 0; j < container.getNSelections(); ++j) { + for (std::size_t j = 0; j < container.getNSelections(); ++j) { std::stringstream line; - std::string sel = useFunctions ? std::string(functions[j].GetExpFormula().Data()) : std::to_string(values[j]); + std::string sel = container.getValueAsString(j); + std::string comment = container.getComment(j); line << " " << std::left << std::setw(25) << sel; @@ -379,8 +372,8 @@ class BaseSelection line << "-> Bit: 0x" << std::hex << std::uppercase << (1ULL << bit) << std::dec; } - if (!comments.empty()) { - line << " (" << comments.at(j) << ")"; + if (!comment.empty()) { + line << " (" << comment << ")"; } LOG(info) << line.str(); } @@ -390,24 +383,27 @@ class BaseSelection LOG(info) << "Printing done"; } + /// \brief Initialize histograms and set bitmask offsets for all configured observables. + /// \tparam HistName Name of the histogram to create in the registry. + /// \param registry Pointer to the histogram registry. template void setupContainers(o2::framework::HistogramRegistry* registry) { mHistRegistry = registry; - // Create histogram with correct number of bins + // create histogram with one bin per selection, plus two summary bins (all analyzed, all passed) int nBins = mNSelection + 2; mHistRegistry->add(HistName, "; Selection Bits; Entries", o2::framework::kTH1F, {{nBins, -0.5, nBins - 0.5}}); - size_t binIndex = 0; + int binIndex = 0; int offset = 0; - for (size_t idx = 0; idx < mSelectionContainers.size(); ++idx) { + for (std::size_t idx = 0; idx < mSelectionContainers.size(); ++idx) { auto& container = mSelectionContainers[idx]; if (container.isEmpty()) { continue; } container.setOffset(offset); offset += container.getShift(); - for (int j = 0; j < container.getNSelections(); j++) { + for (std::size_t j = 0; j < container.getNSelections(); j++) { std::string label = container.getBinLabel(j); mHistRegistry->get(HIST(HistName))->GetXaxis()->SetBinLabel(binIndex + 1, label.c_str()); binIndex++; @@ -418,15 +414,40 @@ class BaseSelection } protected: + void init(int observableIndex) + { + // check if any selections are configured + if (mSelectionContainers.at(observableIndex).isEmpty()) { + return; + } + + // track the number of occupied bits and total selections + mNSelectionBits += mSelectionContainers.at(observableIndex).getShift(); + mNSelection += mSelectionContainers.at(observableIndex).getNSelections(); + + // check if any selection is minimal + if (mSelectionContainers.at(observableIndex).isMinimalCut()) { + mHasMinimalSelection = true; + } + // check if selection is optional + if (mSelectionContainers.at(observableIndex).isOptionalCut()) { + mHasOptionalSelection = true; + } + + if (mNSelectionBits > sizeof(BitmaskType) * CHAR_BIT) { + LOG(fatal) << "Too many selections. At most " << sizeof(BitmaskType) * CHAR_BIT << " number of bits are supported"; + } + } + o2::framework::HistogramRegistry* mHistRegistry = nullptr; - std::array, NumObservables> mSelectionContainers = {}; ///< Array containing all selections - std::bitset mFinalBitmask = {}; ///< final bitmaks - std::size_t mNSelectionBits = 0; ///< Number of selections (all - minimal selections) - int mNSelection = 0; ///< Number of selections all selections - bool mHasMinimalSelection = false; ///< Set to true if all minimal (mandatory) selections are passed - bool mPassesMinimalSelections = true; ///< Set to true if all minimal (mandatory) selections are passed - bool mHasOptionalSelection = false; ///< Set to true if at least one selections is optional - bool mPassesOptionalSelections = false; ///< Set to true if at least one optional (non-mandatory) selections is passed + std::array, NumObservables> mSelectionContainers = {}; ///< Array of selection containers, one per observable + std::bitset mFinalBitmask = {}; ///< Assembled bitmask combining all observable selections + std::size_t mNSelectionBits = 0; ///< Number of bits occupied in the bitmask (excludes skipped most-permissive bits) + std::size_t mNSelection = 0; ///< Total number of configured selection thresholds across all observables + bool mHasMinimalSelection = false; ///< True if at least one observable has a mandatory (minimal) cut configured + bool mPassesMinimalSelections = true; ///< True if all mandatory (minimal) cuts have been passed so far + bool mHasOptionalSelection = false; ///< True if at least one observable has an optional cut configured + bool mPassesOptionalSelections = false; ///< True if at least one optional cut has been passed }; } // namespace o2::analysis::femto diff --git a/PWGCF/Femto/Core/partitions.h b/PWGCF/Femto/Core/partitions.h index d14bbf76196..c4b039e8c05 100644 --- a/PWGCF/Femto/Core/partitions.h +++ b/PWGCF/Femto/Core/partitions.h @@ -24,18 +24,24 @@ (o2::aod::femtocollisions::magField >= static_cast(selection.magFieldMin) && o2::aod::femtocollisions::magField <= static_cast(selection.magFieldMax)) && \ ncheckbit(o2::aod::femtocollisions::mask, selection.collisionMask) +// macro for track momentum, i.e. ||q|*pT/q| * cosh(eta) +// there is no ncosh function, so we have to make our own, i.e. cosh(x) = (exp(x)+exp(-x))/2 +#define TRACK_MOMENTUM(chargeAbs, signedPt, eta) nabs((chargeAbs) * (signedPt)) * (nexp(eta) + nexp(-1.f * (eta))) / 2.f + // standard track partition -#define MAKE_TRACK_PARTITION(selection) \ - ifnode(selection.chargeSign.node() != 0, ifnode(selection.chargeSign.node() > 0, o2::aod::femtobase::stored::signedPt > 0.f, o2::aod::femtobase::stored::signedPt < 0.f), true) && \ - (nabs(selection.chargeAbs.node() * o2::aod::femtobase::stored::signedPt) > selection.ptMin) && \ - (nabs(selection.chargeAbs.node() * o2::aod::femtobase::stored::signedPt) < selection.ptMax) && \ - (o2::aod::femtobase::stored::eta > selection.etaMin) && \ - (o2::aod::femtobase::stored::eta < selection.etaMax) && \ - (o2::aod::femtobase::stored::phi > selection.phiMin) && \ - (o2::aod::femtobase::stored::phi < selection.phiMax) && \ - ifnode(nabs(selection.chargeAbs.node() * o2::aod::femtobase::stored::signedPt) * (nexp(o2::aod::femtobase::stored::eta) + nexp(-1.f * o2::aod::femtobase::stored::eta)) / (2.f) <= selection.pidThres, \ - ncheckbit(o2::aod::femtotracks::mask, selection.maskLowMomentum), \ - ncheckbit(o2::aod::femtotracks::mask, selection.maskHighMomentum)) +#define MAKE_TRACK_PARTITION(selection) \ + ifnode(selection.chargeSign.node() != 0, ifnode(selection.chargeSign.node() > 0, o2::aod::femtobase::stored::signedPt > 0.f, o2::aod::femtobase::stored::signedPt < 0.f), true) && \ + (nabs(selection.chargeAbs * o2::aod::femtobase::stored::signedPt) > selection.ptMin) && \ + (nabs(selection.chargeAbs * o2::aod::femtobase::stored::signedPt) < selection.ptMax) && \ + (o2::aod::femtobase::stored::eta > selection.etaMin) && \ + (o2::aod::femtobase::stored::eta < selection.etaMax) && \ + (o2::aod::femtobase::stored::phi > selection.phiMin) && \ + (o2::aod::femtobase::stored::phi < selection.phiMax) && \ + ifnode(TRACK_MOMENTUM(selection.chargeAbs, o2::aod::femtobase::stored::signedPt, o2::aod::femtobase::stored::eta) <= selection.pidThres, \ + ncheckbit(o2::aod::femtotracks::mask, selection.maskLowMomentum) && \ + (o2::aod::femtotracks::mask & selection.rejectionMaskLowMomentum) == static_cast(0), \ + ncheckbit(o2::aod::femtotracks::mask, selection.maskHighMomentum) && \ + (o2::aod::femtotracks::mask & selection.rejectionMaskHighMomentum) == static_cast(0)) // partition for phis and rhos, i.e. resonance that are their own antiparticle #define MAKE_RESONANCE_0_PARTITON(selection) \ diff --git a/PWGCF/Femto/Core/selectionContainer.h b/PWGCF/Femto/Core/selectionContainer.h index e9d67e6354a..d67819b6f8d 100644 --- a/PWGCF/Femto/Core/selectionContainer.h +++ b/PWGCF/Femto/Core/selectionContainer.h @@ -27,8 +27,10 @@ #include #include #include +#include #include #include +#include #include namespace o2::analysis::femto @@ -37,20 +39,21 @@ namespace o2::analysis::femto /// Limit type for selections namespace limits { -enum LimitType { kUpperLimit, ///< simple upper limit for the value, e.g. p_T < 1 GeV/c - kAbsUpperLimit, ///< upper limit of the absolute value, e.g. |eta| < 0.8 - kLowerLimit, ///< simple lower limit for the value, e.g. p_T > 0.2 GeV/c - kAbsLowerLimit, ///< lower limit of the absolute value, e.g. |DCA_xyz| > 0.05 cm - kUpperFunctionLimit, ///< simple upper limit of a function value, e.g. DCA_xy > f(pt) - kAbsUpperFunctionLimit, ///< upper limit of an absolute value given by a function, e.g. |DCA_xy| > f(pt) - kLowerFunctionLimit, ///< simple lower limit of a function value, e.g. DCA_xy < f(pt) - kAbsLowerFunctionLimit, ///< lower limit of an absolute value given by a function, e.g. |DCA_xy| < f(pt) - kEqual, ///< values need to be equal, e.g. sign = 1 - kEqualArray, ///< values inside an array need to be equal +enum LimitType { kUpperLimit, ///< simple upper limit, e.g. p_T < 1 GeV/c + kAbsUpperLimit, ///< upper limit on the absolute value, e.g. |eta| < 0.8 + kLowerLimit, ///< simple lower limit, e.g. p_T > 0.2 GeV/c + kAbsLowerLimit, ///< lower limit on the absolute value, e.g. |DCA_xy| > 0.05 cm + kUpperFunctionLimit, ///< upper limit given by a function, e.g. DCA_xy < f(pt) + kAbsUpperFunctionLimit, ///< upper limit on the absolute value given by a function, e.g. |DCA_xy| < f(pt) + kLowerFunctionLimit, ///< lower limit given by a function, e.g. DCA_xy > f(pt) + kAbsLowerFunctionLimit, ///< lower limit on the absolute value given by a function, e.g. |DCA_xy| > f(pt) + kEqual, ///< value must equal a fixed threshold, e.g. sign == 1 + kEqualArray, ///< each element of a value array must equal the corresponding threshold + kRange, ///< value must fall within [lower, upper]; either bound can be omitted (open interval) kLimitTypeLast }; -std::unordered_map limitTypeAsStrings = { +inline const std::unordered_map limitTypeAsStrings = { {kUpperLimit, "Upper Limit"}, {kAbsUpperLimit, "Absolute Upper Limit"}, {kLowerLimit, "Lower Limit"}, @@ -61,283 +64,358 @@ std::unordered_map limitTypeAsStrings = { {kAbsLowerFunctionLimit, "Absolute Lower Function Limit"}, {kEqual, "Equal"}, {kEqualArray, "EqualArray"}, - + {kRange, "Range"}, + {kLimitTypeLast, "Last Limit Type"}, }; }; // namespace limits /// \class SelectionContainer -/// \brief Class for storing and evaluating multiple selection thresholds for a single observable. -/// \tparam T Data type for selection values (mostly floats) -/// \tparam BitmaskType Type used for bitmask storage (e.g., uint8_t, uint32_t). +/// \brief Stores and evaluates multiple selection thresholds for a single observable. +/// +/// Selections can be based on static threshold values or dynamically evaluated TF1 functions. +/// Thresholds are sorted from most permissive to most restrictive so that evaluation can bail +/// out early once a threshold fails. +/// +/// \tparam T Data type for selection values (mostly floats). +/// \tparam BitmaskType Integer type used for bitmask storage (e.g., uint8_t, uint32_t). template class SelectionContainer { public: - /// Default constructor + /// \brief Default constructor. SelectionContainer() = default; + + /// \brief Destructor. ~SelectionContainer() = default; - /// \brief Constructor for static value-based selection. - /// \param SelectionValues Vector of values for the selection. + /// \brief Constructor for static value-based selections. + /// \param selectionName Name of the observable this container manages. + /// \param selectionValues Vector of threshold values. /// \param limitType Type of limit (from limits::LimitType). - /// \param SkipMostPermissiveBit Whether to skip the most permissive bit in the bitmask. - /// \param IsMinimalCut Whether this selection should be treated as a minimal required cut. - SelectionContainer(std::string const& SelectionName, - std::vector const& SelectionValues, + /// \param skipMostPermissiveBit Whether to skip the most permissive threshold when assembling the bitmask. + /// \param isMinimalCut Whether this selection is mandatory (candidate is rejected if it fails). + /// \param isOptionalCut Whether this selection is optional (candidate is accepted if any optional cut passes). + SelectionContainer(std::string const& selectionName, + std::vector const& selectionValues, limits::LimitType limitType, - bool SkipMostPermissiveBit, - bool IsMinimalCut, + bool skipMostPermissiveBit, + bool isMinimalCut, bool isOptionalCut) - : mSelectionName(SelectionName), - mSelectionValues(SelectionValues), + : mSelectionName(selectionName), + mSelectionValues(selectionValues), mLimitType(limitType), - mSkipMostPermissiveBit(SkipMostPermissiveBit), - mIsMinimalCut(IsMinimalCut), + mSkipMostPermissiveBit(skipMostPermissiveBit), + mIsMinimalCut(isMinimalCut), mIsOptionalCut(isOptionalCut) { if (mSelectionValues.size() > sizeof(BitmaskType) * CHAR_BIT) { LOG(fatal) << "Too many selections for single a observable. Limit is " << sizeof(BitmaskType) * CHAR_BIT; } + if (isMinimalCut && isOptionalCut) { + LOG(fatal) << "A selection cannot be both minimal and optional"; + } // values for selection are not necessarily ordered correctly sortSelections(); } - /// \brief Constructor for function-based dynamic selection. - /// \param baseName Base name for TF1 functions. - /// \param lowerLimit Lower bound for TF1 domain. - /// \param upperLimit Upper bound for TF1 domain. - /// \param functions Vector of strings defining TF1 functions. - /// \param limitType Type of limit. - /// \param skipMostPermissiveBit Whether to skip the most permissive bit in the bitmask. - /// \param IsMinimalCut Whether this selection should be treated as a minimal required cut. - SelectionContainer(std::string const& SelectionName, + /// \brief Constructor for range-based selections defined as "lower;upper" strings. + /// Either bound may be omitted to represent an open interval, e.g. ";1.0" means value < 1.0. + /// \param selectionName Name of the observable this container manages. + /// \param rangeStrings Vector of range strings, each of the form "lower;upper". + /// \param skipMostPermissiveBit Whether to skip the most permissive threshold when assembling the bitmask. + /// \param isMinimalCut Whether this selection is mandatory (candidate is rejected if it fails). + /// \param isOptionalCut Whether this selection is optional (candidate is accepted if any optional cut passes). + SelectionContainer(std::string const& selectionName, + std::vector const& rangeStrings, + bool skipMostPermissiveBit, + bool isMinimalCut, + bool isOptionalCut) + : mSelectionName(selectionName), + mLimitType(limits::kRange), + mSkipMostPermissiveBit(skipMostPermissiveBit), + mIsMinimalCut(isMinimalCut), + mIsOptionalCut(isOptionalCut) + { + if (rangeStrings.size() > sizeof(BitmaskType) * CHAR_BIT) { + LOG(fatal) << "Too many selections for a single observable. Limit is " << sizeof(BitmaskType) * CHAR_BIT; + } + if (isMinimalCut && isOptionalCut) { + LOG(fatal) << "A selection cannot be both minimal and optional"; + } + for (auto const& rangeStr : rangeStrings) { + T lower, upper; + parseRangeString(rangeStr, lower, upper); + mSelectionRanges.emplace_back(lower, upper); + } + // ranges are sorted by their lenghts, i.e. from widest range to tightest range + // in principle the ranges do not have to include each other, exepct this is configured as minimal cut, check is added at the end + // to cover both cases, this also means we always check all ranges, when assembling the bit mask + sortSelections(); + // init mSelectionValues to be the widths of the intervals + for (std::size_t i = 0; i < mSelectionRanges.size(); i++) { + mSelectionValues.push_back(mSelectionRanges[i].second - mSelectionRanges[i].first); + } + + // for minimal range cut, ranges must be strictly nested (each range contains all narrower ones) + // this is required for the early-exit logic in passesAsMinimalCut() to be correct + if (isMinimalCut) { + for (std::size_t i = 0; i + 1 < mSelectionRanges.size(); i++) { + if (mSelectionRanges[i].first > mSelectionRanges[i + 1].first || + mSelectionRanges[i].second < mSelectionRanges[i + 1].second) { + LOG(fatal) << "Ranges for minimal cut " << selectionName + << " are not nested. Range [" << mSelectionRanges[i].first << ";" << mSelectionRanges[i].second + << "] does not contain [" << mSelectionRanges[i + 1].first << ";" << mSelectionRanges[i + 1].second << "]"; + } + } + } + } + + /// \brief Constructor for function-based dynamic selections. + /// \param selectionName Name of the observable this container manages. + /// \param lowerLimit Lower bound of the TF1 domain. + /// \param upperLimit Upper bound of the TF1 domain. + /// \param functions Vector of strings defining TF1 threshold functions. + /// \param limitType Type of limit (from limits::LimitType). + /// \param skipMostPermissiveBit Whether to skip the most permissive threshold when assembling the bitmask. + /// \param isMinimalCut Whether this selection is mandatory (candidate is rejected if it fails). + /// \param isOptionalCut Whether this selection is optional (candidate is accepted if any optional cut passes). + SelectionContainer(std::string const& selectionName, T lowerLimit, T upperLimit, std::vector const& functions, limits::LimitType limitType, bool skipMostPermissiveBit, - bool IsMinimalCut, + bool isMinimalCut, bool isOptionalCut) - : mSelectionName(SelectionName), + : mSelectionName(selectionName), mLimitType(limitType), mSkipMostPermissiveBit(skipMostPermissiveBit), - mIsMinimalCut(IsMinimalCut), + mIsMinimalCut(isMinimalCut), mIsOptionalCut(isOptionalCut) { if (functions.size() > sizeof(BitmaskType) * CHAR_BIT) { LOG(fatal) << "Too many selections for single a observable. Limit is " << sizeof(BitmaskType) * CHAR_BIT; } + if (isMinimalCut && isOptionalCut) { + LOG(fatal) << "A selection cannot be both minimal and optional"; + } for (std::size_t i = 0; i < functions.size(); i++) { mSelectionFunctions.emplace_back((mSelectionName + std::to_string(i)).c_str(), functions.at(i).c_str(), lowerLimit, upperLimit); } - // functions for selection are not necessarily ordered correctly - // use value at midpoint to order them - // here we rely on the user that the functions can be ordered like this over the whole interval + // functions are not necessarily ordered correctly; + // use the midpoint of the domain to establish their order. + // this relies on the user ensuring the ordering is consistent across the whole interval. T midPoint = (lowerLimit + upperLimit) / 2.; sortFunctions(midPoint); - // initialize the values also to the midpoint + // initialize threshold values to the functions evaluated at the midpoint for (std::size_t i = 0; i < functions.size(); i++) { mSelectionValues.push_back(mSelectionFunctions.at(i).Eval(midPoint)); } } - /// \brief Sort static selection values based on the limit type. - void sortSelections() + /// \brief Attach comments to the selection thresholds, one per threshold in the same order as the input values. + /// \param comments Vector of comment strings. + void addComments(std::vector const& comments) { - switch (mLimitType) { - case (limits::kUpperLimit): - case (limits::kAbsUpperLimit): - std::sort(mSelectionValues.begin(), mSelectionValues.end(), [](T a, T b) { return a > b; }); - break; - case (limits::kLowerLimit): - case (limits::kAbsLowerLimit): - std::sort(mSelectionValues.begin(), mSelectionValues.end(), [](T a, T b) { return a < b; }); - break; - default: - break; + // note: threshold values may be reordered by sortSelections() or sortFunctions(), + // so comments must be provided in the already-sorted order + if (comments.size() != getNSelections()) { + LOG(fatal) << "Number of comments and number of selections are inconsistent"; } + mComments = comments; } - /// \brief Sort selection functions based on evaluation at a given point. - /// \param value Point at which to evaluate the functions for ordering. - void sortFunctions(T value) + /// \brief Get comments attached to the selection thresholds. + /// \return Vector of comment strings. + std::string getComment(int selectionIndex) const { - switch (mLimitType) { - case (limits::kUpperFunctionLimit): - case (limits::kAbsUpperFunctionLimit): - std::sort(mSelectionFunctions.begin(), mSelectionFunctions.end(), [value](TF1 const& a, TF1 const& b) { return a.Eval(value) > b.Eval(value); }); - break; - case (limits::kLowerFunctionLimit): - case (limits::kAbsLowerFunctionLimit): - std::sort(mSelectionFunctions.begin(), mSelectionFunctions.end(), [value](TF1 const& a, TF1 const& b) { return a.Eval(value) < b.Eval(value); }); - break; - default: - break; + if (mComments.empty()) { + return std::string(""); } + return mComments.at(selectionIndex); } - /// \brief Add comments to the selection values - /// \param comments Vector of comments - void addComments(std::vector const& comments) - { - // make sure that the comments are in correct order - // the values passed to the selection container can be reordered based on the limit type - mComments = comments; - } - - std::vector const& getComments() const { return mComments; } + /// \brief Get the name of this selection. + /// \return Selection name string. std::string const& getSelectionName() const { return mSelectionName; } - /// \brief Update selection limits using internal functions evaluated at a given value. - /// \param value Input value to evaluate functions at. + /// \brief Update threshold values by re-evaluating the internal TF1 functions at a given point. + /// \param value Input value at which to evaluate the functions. void updateLimits(T value) { - // functions are ordered so just add the values in the same order + // functions are already sorted, so evaluate in the same order as mSelectionValues for (std::size_t i = 0; i < mSelectionValues.size(); i++) { mSelectionValues.at(i) = mSelectionFunctions.at(i).Eval(value); } } - /// \brief Evaluate which selection criteria are fulfilled for a given value. + /// \brief Evaluate which selection thresholds are passed for a given observable value. /// \param value Value of the observable to evaluate. void evaluate(T value) { - // better safe than sorry and reset the bitmask before you evaluate and set minimal selection to true + // reset the bitmask before evaluating mBitmask.reset(); - // the values are ordered, from most loose to most tight, as soon as one comparison is not true, we can break out of the loop - bool breakLoop = false; - // iterate over all limits and set the corresponding bit if we pass the selection, otherwise break out as soon as we can - // only break if the observable is used for the minimal selection - for (size_t i = 0; i < mSelectionValues.size(); i++) { - switch (mLimitType) { - case (limits::kUpperLimit): - case (limits::kUpperFunctionLimit): + + // switch on limit type once outside the loop; + // thresholds are sorted from most permissive to most restrictive, + // so we can break early as soon as one comparison fails + switch (mLimitType) { + case (limits::kUpperLimit): + case (limits::kUpperFunctionLimit): + for (std::size_t i = 0; i < mSelectionValues.size(); i++) { if (value <= mSelectionValues.at(i)) { mBitmask.set(i); } else { - breakLoop = true; + break; } - break; - case (limits::kAbsUpperLimit): - case (limits::kAbsUpperFunctionLimit): + } + break; + case (limits::kAbsUpperLimit): + case (limits::kAbsUpperFunctionLimit): + for (std::size_t i = 0; i < mSelectionValues.size(); i++) { if (std::abs(value) <= mSelectionValues.at(i)) { mBitmask.set(i); } else { - breakLoop = true; + break; } - break; - case (limits::kLowerLimit): - case (limits::kLowerFunctionLimit): + } + break; + case (limits::kLowerLimit): + case (limits::kLowerFunctionLimit): + for (std::size_t i = 0; i < mSelectionValues.size(); i++) { if (value >= mSelectionValues.at(i)) { mBitmask.set(i); } else { - breakLoop = true; + break; } - break; - case (limits::kAbsLowerLimit): - case (limits::kAbsLowerFunctionLimit): + } + break; + case (limits::kAbsLowerLimit): + case (limits::kAbsLowerFunctionLimit): + for (std::size_t i = 0; i < mSelectionValues.size(); i++) { if (std::abs(value) >= mSelectionValues.at(i)) { mBitmask.set(i); } else { - breakLoop = true; + break; } - break; - case (limits::kEqual): - // special case for kEqual since here we cannot really establish an order so we need to check all cases explicitly and we cannot bail early + } + break; + case (limits::kRange): + // ranges are sorted widest-first but a narrower range passing does not imply a wider one fails (no check on boundries), + // so all ranges must be checked explicitly — no early exit + for (std::size_t i = 0; i < mSelectionRanges.size(); i++) { + if (value >= mSelectionRanges.at(i).first && value <= mSelectionRanges.at(i).second) { + mBitmask.set(i); + } + } + break; + case (limits::kEqual): + // kEqual has no natural ordering, so all thresholds must be checked and we cannot bail early + for (std::size_t i = 0; i < mSelectionValues.size(); i++) { if (std::fabs(value - mSelectionValues.at(i)) < constants::math::Epsilon) { mBitmask.set(i); } - break; - default: - breakLoop = true; - } - // bail early if a comparison fails - // the values are ordered, so all following we also fail, there there is no point in contiuing - if (breakLoop) { + } break; - } + default: + LOG(warn) << "Limit type not known, no selection is applied"; } } - /// \brief Evaluate which selection criteria are fulfilled for a given value. - /// \param values Values of the observable to evaluate + /// \brief Evaluate which selection thresholds are passed for a vector of observable values. + /// Only kEqualArray is supported for now; each element is compared against the corresponding threshold. + /// \param values Values of the observable to evaluate. void evaluate(std::vector const& values) { if (values.size() != mSelectionValues.size()) { LOG(fatal) << "Wrong number of values have been passed"; } - for (size_t i = 0; i < mSelectionValues.size(); i++) { - switch (mLimitType) { - case (limits::kEqualArray): + // reset the bitmask before evaluating + mBitmask.reset(); + + switch (mLimitType) { + case (limits::kEqualArray): + for (std::size_t i = 0; i < mSelectionValues.size(); i++) { if (std::fabs(values.at(i) - mSelectionValues.at(i)) < constants::math::Epsilon) { mBitmask.set(i); } - break; - default: - continue; - } + } + break; + default: + LOG(warn) << "Limit type not known, no selection is applied"; } } - /// \brief Retrieve the bitmask indicating which selections were passed. + /// \brief Retrieve the bitmask indicating which thresholds were passed. + /// If mSkipMostPermissiveBit is set, the bit for the loosest threshold is removed by shifting. /// \return Bitset representing passed selections. std::bitset getBitmask() const { - // if we do not skip the last bit, return full bitmask - if (mSkipMostPermissiveBit == false) { + if (!mSkipMostPermissiveBit) { return mBitmask; } else { - // for the other selections we can remove the first bit since it is the minimal selection and therefore always true + // remove the first (most permissive) bit since it corresponds to the minimal selection and is always true return mBitmask >> 1; } } + + /// \brief Manually set the internal bitmask. + /// \tparam R Integer type of the bitmask value. + /// \param bitmask Bitmask value to set. template void setBitmask(R bitmask) { mBitmask = std::bitset(bitmask); } - /// \brief Check whether the minimal cut condition is fulfilled. - /// \return True if minimal selection is fulfilled, false otherwise. + /// \brief Reset the internal bitmask to zero. + void reset() { mBitmask.reset(); } + + /// \brief Check whether the mandatory (minimal) cut condition is fulfilled. + /// \return True if the minimal selection passes or if this container is not marked as a minimal cut. bool passesAsMinimalCut() const { if (mIsMinimalCut) { - // check if any bit is set - // in case were bits are evaluted in order, if the loosests fails, all fail, so testing any is safe here - return mBitmask.any(); + // if any bit is set the loosest threshold passed; since thresholds are ordered, + // the loosest passing implies all looser ones also pass + return mBitmask.test(0); } - // if a selection is not marked as minimal cut we return true by default + // not a minimal cut — return true by default so it does not block the candidate return true; } - /// \brief Check whether any optional cuts are fulfilled. - /// \return True if at least one optional cut is passed. + /// \brief Check whether any optional cut is fulfilled. + /// \return True if at least one optional threshold is passed, false if this container is not marked as optional. bool passesAsOptionalCut() const { if (mIsOptionalCut) { - // check if any bit is set - // in case were bits are evaluted in order, if the loosests fails, all fail, so testing any is safe here + // if any bit is set the loosest threshold passed return mBitmask.any(); } - // if a selection is not marked as optional cut we return false by default + // not an optional cut — return false by default so it does not accidentally accept the candidate return false; } - /// \brief Get the loosest (most permissive) selection value. - /// \return First (loosest) selection value. - T getLoosestSelection() const { return mSelectionValues.at(0); } + /// \brief Get the loosest (most permissive) selection threshold value. + /// \return The first element of the sorted threshold vector. + T getLoosestSelection() const + { + if (mSelectionValues.empty()) { + LOG(fatal) << "No selections configured"; + } + return mSelectionValues.at(0); + } - /// \brief Check if there are any selection values configured. We also init values in case of function so this is safe - /// \return True if no selections are configured. + /// \brief Check whether any selection thresholds are configured. + /// For function-based selections, mSelectionValues is always populated (initialised at the midpoint), + /// so this check is safe for both static and function-based containers. + /// \return True if no thresholds are configured. bool isEmpty() const { return mSelectionValues.empty(); } - /// \brief Check if there are any selection values configured. - /// \return True if no selections are configured. - bool isUsingFunctions() const { return !mSelectionFunctions.empty(); } - - /// \brief Get the number of bits to shift for the final bitmask. - /// \return Number of bits to shift. + /// \brief Get the number of bits this container contributes to the global bitmask. + /// If the most permissive bit is skipped, the contribution is reduced by one. + /// \return Number of bits to add to the global bitmask offset. int getShift() const { if (mSelectionValues.empty()) { @@ -345,22 +423,47 @@ class SelectionContainer } if (mSkipMostPermissiveBit) { return static_cast(mSelectionValues.size() - 1); - } else { - return static_cast(mSelectionValues.size()); } + return static_cast(mSelectionValues.size()); } + /// \brief Set the bit offset of this container within the global bitmask. + /// \param offset Bit position at which this container's bits start. void setOffset(int offset) { mOffset = offset; } + + /// \brief Get the bit offset of this container within the global bitmask. + /// \return Bit offset. int getOffset() const { return mOffset; } - int getNSelections() const { return mSelectionValues.size(); } + /// \brief Get the total number of configured selection thresholds. + /// \return Number of thresholds. + std::size_t getNSelections() const { return mSelectionValues.size(); } + /// \brief Build a histogram bin label string encoding the full configuration of a single threshold. + /// \param selectionIndex Index of the threshold within this container. + /// \return Encoded label string. std::string getBinLabel(int selectionIndex) const { std::ostringstream oss; std::string sectionDelimiter = ":::"; std::string valueDelimiter = "___"; std::string noValue = "X"; + + // Determine value string + std::string valueStr; + + if (mLimitType == limits::kRange) { + // Print actual lower;upper interval + const auto& range = mSelectionRanges.at(selectionIndex); + std::ostringstream rangeStream; + rangeStream << range.first << ";" << range.second; + valueStr = rangeStream.str(); + } else if (mSelectionFunctions.empty()) { + valueStr = std::to_string(mSelectionValues.at(selectionIndex)); + } else { + valueStr = mSelectionFunctions.at(selectionIndex).GetExpFormula().Data(); + } + oss << "SelectionName" << valueDelimiter << mSelectionName << sectionDelimiter << "LimitType" << valueDelimiter << getLimitTypeAsString() << sectionDelimiter << "MinimalCut" << valueDelimiter << (mIsMinimalCut ? "1" : "0") << sectionDelimiter @@ -368,12 +471,32 @@ class SelectionContainer << "OptionalCut" << valueDelimiter << (mIsOptionalCut ? "1" : "0") << sectionDelimiter << "Shift" << valueDelimiter << getShift() << sectionDelimiter << "Offset" << valueDelimiter << mOffset << sectionDelimiter - << "Value" << valueDelimiter << (mSelectionFunctions.empty() ? std::to_string(mSelectionValues.at(selectionIndex)) : mSelectionFunctions.at(selectionIndex).GetExpFormula().Data()) << sectionDelimiter + << "Value" << valueDelimiter << valueStr << sectionDelimiter << "BitPosition" << valueDelimiter << (mSkipMostPermissiveBit ? (selectionIndex == 0 ? noValue : std::to_string(mOffset + selectionIndex - 1)) : std::to_string(mOffset + selectionIndex)) << sectionDelimiter << "Comment" << valueDelimiter << (mComments.empty() ? noValue : mComments.at(selectionIndex)); return oss.str(); } + std::string getValueAsString(int selectionIndex) const + { + if (this->isEmpty()) { + return std::string("No value configured"); + } + if (!mSelectionFunctions.empty()) { + return std::string(mSelectionFunctions.at(selectionIndex).GetExpFormula().Data()); + } + if (!mSelectionRanges.empty()) { + std::ostringstream oss; + oss << mSelectionRanges.at(selectionIndex).first << ";" << mSelectionRanges.at(selectionIndex).second; + return oss.str(); + } + return std::to_string(mSelectionValues.at(selectionIndex)); + } + + /// \brief Get the global bit position of a threshold within the final bitmask. + /// Calling this for the skipped most-permissive threshold is a fatal error. + /// \param selectionIndex Index of the threshold within this container. + /// \return Global bit position. int getBitPosition(int selectionIndex) const { if (selectionIndex == 0 && mSkipMostPermissiveBit) { @@ -387,43 +510,107 @@ class SelectionContainer } } - /// \brief Get string representation of the limit type. - /// \return String name of the limit type. - std::string getLimitTypeAsString() const { return limits::limitTypeAsStrings[mLimitType]; } + /// \brief Get the string representation of the configured limit type. + /// \return Human-readable limit type name. + std::string getLimitTypeAsString() const { return limits::limitTypeAsStrings.at(mLimitType); } - /// \brief Get a copy of all selection values. - /// \return Vector of selection values. + /// \brief Get the configured static threshold values. + /// \return Const reference to the vector of threshold values. std::vector const& getSelectionValues() const { return mSelectionValues; } - /// \brief Get a copy of all selection values. - /// \return Vector of selection values. + /// \brief Get the configured TF1 threshold functions. + /// \return Const reference to the vector of TF1 functions. std::vector const& getSelectionFunction() const { return mSelectionFunctions; } - /// \brief Check if this container is marked as minimal cut - /// \return True if minimal cut, false otherwise. + /// \brief Check whether this container is marked as a mandatory (minimal) cut. + /// \return True if this is a minimal cut. bool isMinimalCut() const { return mIsMinimalCut; } - /// \brief Check if this container is marked as optional cut - /// \return True if minimal cut, false otherwise. + /// \brief Check whether this container is marked as an optional cut. + /// \return True if this is an optional cut. bool isOptionalCut() const { return mIsOptionalCut; } - /// \brief Check whether the most permissive bit is skipped. - /// \return True if skipped, false otherwise. + /// \brief Check whether the most permissive threshold bit is skipped when assembling the bitmask. + /// \return True if the most permissive bit is skipped. bool skipMostPermissiveBit() const { return mSkipMostPermissiveBit; } - void reset() { mBitmask.reset(); } - private: - std::string mSelectionName = std::string(""); - std::vector mSelectionValues = {}; ///< Values used for the selection - std::vector mSelectionFunctions = {}; ///< Function used for the selection - limits::LimitType mLimitType = limits::kLimitTypeLast; ///< Limit type of selection - bool mSkipMostPermissiveBit = false; ///< whether to skip the last bit or not - bool mIsMinimalCut = false; ///< whether to use this observable for minimal selection or not - bool mIsOptionalCut = false; ///< whether to use this observable for minimal selection or not - std::vector mComments = {}; ///< Comments for the values - std::bitset mBitmask = {}; ///< bitmask for the observable - int mOffset = 0; + /// \brief Sort static threshold values from most permissive to most restrictive based on the limit type. + void sortSelections() + { + switch (mLimitType) { + case (limits::kUpperLimit): + case (limits::kAbsUpperLimit): + std::sort(mSelectionValues.begin(), mSelectionValues.end(), [](T a, T b) { return a > b; }); + break; + case (limits::kLowerLimit): + case (limits::kAbsLowerLimit): + std::sort(mSelectionValues.begin(), mSelectionValues.end(), [](T a, T b) { return a < b; }); + break; + case (limits::kRange): + // sort by range width descending so the most permissive (widest) range comes first + std::sort(mSelectionRanges.begin(), mSelectionRanges.end(), + [](std::pair const& a, std::pair const& b) { + return (a.second - a.first) > (b.second - b.first); + }); + break; + default: + break; + } + } + + /// \brief Sort selection functions from most permissive to most restrictive, evaluated at a given point. + /// \param value Point at which to evaluate the functions for ordering. + void sortFunctions(T value) + { + switch (mLimitType) { + case (limits::kUpperFunctionLimit): + case (limits::kAbsUpperFunctionLimit): + std::sort(mSelectionFunctions.begin(), mSelectionFunctions.end(), [value](TF1 const& a, TF1 const& b) { return a.Eval(value) > b.Eval(value); }); + break; + case (limits::kLowerFunctionLimit): + case (limits::kAbsLowerFunctionLimit): + std::sort(mSelectionFunctions.begin(), mSelectionFunctions.end(), [value](TF1 const& a, TF1 const& b) { return a.Eval(value) < b.Eval(value); }); + break; + default: + break; + } + } + + /// \brief Parse a range string of the form "lower;upper" into a lower and upper bound. + /// Either bound may be omitted (e.g. ";1.0" or "0.5;"), in which case + /// -/+ numeric_limits::max() is used respectively. + /// \param rangeStr Input string to parse. + /// \param lower Output lower bound. + /// \param upper Output upper bound. + void parseRangeString(std::string const& rangeStr, T& lower, T& upper) const + { + auto pos = rangeStr.find(';'); + if (pos == std::string::npos) { + LOG(fatal) << "Range string '" << rangeStr << "' is missing ';' separator. Expected format: 'lower;upper'"; + } + std::string lowerStr = rangeStr.substr(0, pos); + std::string upperStr = rangeStr.substr(pos + 1); + + lower = lowerStr.empty() ? -std::numeric_limits::max() : static_cast(std::stod(lowerStr)); + upper = upperStr.empty() ? std::numeric_limits::max() : static_cast(std::stod(upperStr)); + + if (lower >= upper) { + LOG(fatal) << "Range string '" << rangeStr << "' has lower bound >= upper bound"; + } + } + + std::string mSelectionName = ""; + std::vector mSelectionValues = {}; ///< Threshold values, sorted from most permissive to most restrictive + std::vector mSelectionFunctions = {}; ///< TF1 threshold functions (empty for static selections) + std::vector> mSelectionRanges = {}; ///< Lower and upper bounds for kRange selections, one pair per threshold + limits::LimitType mLimitType = limits::kLimitTypeLast; ///< Comparison type applied during evaluation + bool mSkipMostPermissiveBit = false; ///< If true, the most permissive threshold does not occupy a bit in the global bitmask + bool mIsMinimalCut = false; ///< If true, this selection is mandatory; failing it rejects the candidate + bool mIsOptionalCut = false; ///< If true, this selection is optional; passing it accepts the candidate + std::vector mComments = {}; ///< Optional comments per threshold, in the same order as mSelectionValues + std::bitset mBitmask = {}; ///< Bitmask indicating which thresholds were passed during the last evaluation + int mOffset = 0; ///< Bit offset of this container within the global bitmask }; } // namespace o2::analysis::femto diff --git a/PWGCF/Femto/Core/trackBuilder.h b/PWGCF/Femto/Core/trackBuilder.h index c33b5850468..ffb5df80c37 100644 --- a/PWGCF/Femto/Core/trackBuilder.h +++ b/PWGCF/Femto/Core/trackBuilder.h @@ -67,65 +67,65 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { // Electron PID cuts o2::framework::Configurable requirePidElectron{"requirePidElectron", false, "Make election PID optional"}; o2::framework::Configurable minMomTofElectron{"minMomTofElectron", 0.3, "Minimum momentum to required TOF PID for Electron"}; - o2::framework::Configurable> itsElectron{"itsElectron", {}, "Maximum |nsigma| for Electron PID"}; - o2::framework::Configurable> tpcElectron{"tpcElectron", {}, "Maximum |nsigma| for Electron PID"}; - o2::framework::Configurable> tofElectron{"tofElectron", {}, "Maximum |nsigma| for Electron PID"}; - o2::framework::Configurable> tpcitsElectron{"tpcitsElectron", {}, "Maximum |nsigma| for Electron PID"}; - o2::framework::Configurable> tpctofElectron{"tpctofElectron", {}, "Maximum |nsigma| for Electron PID"}; + o2::framework::Configurable> itsElectron{"itsElectron", {}, "Ranges LowerLimit;UpperLimit for nsigma_ITS for Electron PID"}; + o2::framework::Configurable> tpcElectron{"tpcElectron", {}, "Ranges LowerLimit;UpperLimit for nsigma_TPC for Electron PID"}; + o2::framework::Configurable> tofElectron{"tofElectron", {}, "Ranges LowerLimit;UpperLimit for nsigma_TOF for Electron PID"}; + o2::framework::Configurable> tpcitsElectron{"tpcitsElectron", {}, "Maximum nsigma_TPCITS for Electron PID"}; + o2::framework::Configurable> tpctofElectron{"tpctofElectron", {}, "Maximum nsigma_TPCTOF for Electron PID"}; // Pion PID cuts o2::framework::Configurable requirePidPion{"requirePidPion", false, "Make election PID optional"}; o2::framework::Configurable minMomTofPion{"minMomTofPion", 0.5, "Minimum momentum to required TOF PID for Pion"}; - o2::framework::Configurable> itsPion{"itsPion", {}, "Maximum |nsigma| for Pion PID"}; - o2::framework::Configurable> tpcPion{"tpcPion", {}, "Maximum |nsigma| for Pion PID"}; - o2::framework::Configurable> tofPion{"tofPion", {}, "Maximum |nsigma| for Pion PID"}; - o2::framework::Configurable> tpcitsPion{"tpcitsPion", {}, "Maximum |nsigma| for Pion PID"}; - o2::framework::Configurable> tpctofPion{"tpctofPion", {}, "Maximum |nsigma| for Pion PID"}; + o2::framework::Configurable> itsPion{"itsPion", {}, "Ranges LowerLimit;UpperLimit for nsigma_ITS for Pion PID"}; + o2::framework::Configurable> tpcPion{"tpcPion", {}, "Ranges LowerLimit;UpperLimit for nsigma_TPC for Pion PID"}; + o2::framework::Configurable> tofPion{"tofPion", {}, "Ranges LowerLimit;UpperLimit for nsigma_TOF for Pion PID"}; + o2::framework::Configurable> tpcitsPion{"tpcitsPion", {}, "Maximum nsigma_TPCITS for Pion PID"}; + o2::framework::Configurable> tpctofPion{"tpctofPion", {}, "Maximum nsigma_TPCTOF for Pion PID"}; // Kaon PID cuts o2::framework::Configurable requirePidKaon{"requirePidKaon", false, "Make election PID optional"}; o2::framework::Configurable minMomTofKaon{"minMomTofKaon", 0.4, "Minimum momentum to required TOF PID for Kaon"}; - o2::framework::Configurable> itsKaon{"itsKaon", {}, "Maximum |nsigma| for Kaon PID"}; - o2::framework::Configurable> tpcKaon{"tpcKaon", {}, "Maximum |nsigma| for Kaon PID"}; - o2::framework::Configurable> tofKaon{"tofKaon", {}, "Maximum |nsigma| for Kaon PID"}; - o2::framework::Configurable> tpcitsKaon{"tpcitsKaon", {}, "Maximum |nsigma| for Kaon PID"}; - o2::framework::Configurable> tpctofKaon{"tpctofKaon", {}, "Maximum |nsigma| for Kaon PID"}; + o2::framework::Configurable> itsKaon{"itsKaon", {}, "Ranges LowerLimit;UpperLimit for nsigma_ITS for Kaon PID"}; + o2::framework::Configurable> tpcKaon{"tpcKaon", {}, "Ranges LowerLimit;UpperLimit for nsigma_TPC for Kaon PID"}; + o2::framework::Configurable> tofKaon{"tofKaon", {}, "Ranges LowerLimit;UpperLimit for nsigma_TOF for Kaon PID"}; + o2::framework::Configurable> tpcitsKaon{"tpcitsKaon", {}, "Maximum nsigma_TPCITS for Kaon PID"}; + o2::framework::Configurable> tpctofKaon{"tpctofKaon", {}, "Maximum nsigma_TPCTOF for Kaon PID"}; // Proton PID cuts o2::framework::Configurable requirePidProton{"requirePidProton", true, "Make election PID optional"}; o2::framework::Configurable minMomTofProton{"minMomTofProton", 0.75, "Minimum momentum to required TOF PID for Proton"}; - o2::framework::Configurable> itsProton{"itsProton", {}, "Maximum |nsigma| for Proton PID"}; - o2::framework::Configurable> tpcProton{"tpcProton", {}, "Maximum |nsigma| for Proton PID"}; - o2::framework::Configurable> tofProton{"tofProton", {}, "Maximum |nsigma| for Proton PID"}; - o2::framework::Configurable> tpcitsProton{"tpcitsProton", {}, "Maximum |nsigma| for Proton PID"}; - o2::framework::Configurable> tpctofProton{"tpctofProton", {}, "Maximum |nsigma| for Proton PID"}; + o2::framework::Configurable> itsProton{"itsProton", {}, "Ranges LowerLimit;UpperLimit for nsigma_ITS for Proton PID"}; + o2::framework::Configurable> tpcProton{"tpcProton", {}, "Ranges LowerLimit;UpperLimit for nsigma_TPC for Proton PID"}; + o2::framework::Configurable> tofProton{"tofProton", {}, "Ranges LowerLimit;UpperLimit for nsigma_TOF for Proton PID"}; + o2::framework::Configurable> tpcitsProton{"tpcitsProton", {}, "Maximum nsigma_TPCITS for Proton PID"}; + o2::framework::Configurable> tpctofProton{"tpctofProton", {}, "Maximum nsigma_TPCTOF for Proton PID"}; // Deuteron PID cuts o2::framework::Configurable requirePidDeuteron{"requirePidDeuteron", false, "Make election PID optional"}; o2::framework::Configurable minMomTofDeuteron{"minMomTofDeuteron", 1.2, "Minimum momentum to required TOF PID for Deuteron"}; - o2::framework::Configurable> itsDeuteron{"itsDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; - o2::framework::Configurable> tpcDeuteron{"tpcDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; - o2::framework::Configurable> tofDeuteron{"tofDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; - o2::framework::Configurable> tpcitsDeuteron{"tpcitsDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; - o2::framework::Configurable> tpctofDeuteron{"tpctofDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; + o2::framework::Configurable> itsDeuteron{"itsDeuteron", {}, "Ranges LowerLimit;UpperLimit for nsigma_ITS for Deuteron PID"}; + o2::framework::Configurable> tpcDeuteron{"tpcDeuteron", {}, "Ranges LowerLimit;UpperLimit for nsigma_TPC for Deuteron PID"}; + o2::framework::Configurable> tofDeuteron{"tofDeuteron", {}, "Ranges LowerLimit;UpperLimit for nsigma_TOF for Deuteron PID"}; + o2::framework::Configurable> tpcitsDeuteron{"tpcitsDeuteron", {}, "Maximum nsigma_TPCITS for Deuteron PID"}; + o2::framework::Configurable> tpctofDeuteron{"tpctofDeuteron", {}, "Maximum nsigma_TPCTOF for Deuteron PID"}; // Triton PID cuts o2::framework::Configurable requirePidTriton{"requirePidTriton", false, "Make election PID optional"}; o2::framework::Configurable minMomTofTriton{"minMomTofTriton", 1.4, "Minimum momentum to required TOF PID for Triton"}; - o2::framework::Configurable> itsTriton{"itsTriton", {}, "Maximum |nsigma| for Triton PID"}; - o2::framework::Configurable> tpcTriton{"tpcTriton", {}, "Maximum |nsigma| for Triton PID"}; - o2::framework::Configurable> tofTriton{"tofTriton", {}, "Maximum |nsigma| for Triton PID"}; - o2::framework::Configurable> tpcitsTriton{"tpcitsTriton", {}, "Maximum |nsigma| for Triton PID"}; - o2::framework::Configurable> tpctofTriton{"tpctofTriton", {}, "Maximum |nsigma| for Triton PID"}; + o2::framework::Configurable> itsTriton{"itsTriton", {}, "Ranges LowerLimit;UpperLimit for nsigma_ITS for Triton PID"}; + o2::framework::Configurable> tpcTriton{"tpcTriton", {}, "Ranges LowerLimit;UpperLimit for nsigma_TPC for Triton PID"}; + o2::framework::Configurable> tofTriton{"tofTriton", {}, "Ranges LowerLimit;UpperLimit for nsigma_TOF for Triton PID"}; + o2::framework::Configurable> tpcitsTriton{"tpcitsTriton", {}, "Maximum nsigma_TPCITS for Triton PID"}; + o2::framework::Configurable> tpctofTriton{"tpctofTriton", {}, "Maximum nsigma_TPCTOF for Triton PID"}; // Helium PID cuts o2::framework::Configurable requirePidHelium{"requirePidHelium", false, "Make election PID optional"}; o2::framework::Configurable minMomTofHelium{"minMomTofHelium", 1.6, "Minimum momentum to required TOF PID for Helium"}; - o2::framework::Configurable> itsHelium{"itsHelium", {}, "Maximum |nsigma| for Helium PID"}; - o2::framework::Configurable> tpcHelium{"tpcHelium", {}, "Maximum |nsigma| for Helium PID"}; - o2::framework::Configurable> tofHelium{"tofHelium", {}, "Maximum |nsigma| for Helium PID"}; - o2::framework::Configurable> tpcitsHelium{"tpcitsHelium", {}, "Maximum |nsigma| for Helium PID"}; - o2::framework::Configurable> tpctofHelium{"tpctofHelium", {}, "Maximum |nsigma| for Helium PID"}; + o2::framework::Configurable> itsHelium{"itsHelium", {}, "Ranges LowerLimit;UpperLimit for nsigma_ITS for Helium PID"}; + o2::framework::Configurable> tpcHelium{"tpcHelium", {}, "Ranges LowerLimit;UpperLimit for nsigma_TPC for Helium PID"}; + o2::framework::Configurable> tofHelium{"tofHelium", {}, "Ranges LowerLimit;UpperLimit for nsigma_TOF for Helium PID"}; + o2::framework::Configurable> tpcitsHelium{"tpcitsHelium", {}, "Maximum nsigma_TPCITS for Helium PID"}; + o2::framework::Configurable> tpctofHelium{"tpctofHelium", {}, "Maximum nsigma_TPCTOF for Helium PID"}; }; // define the template structure for TrackSelection @@ -144,8 +144,11 @@ struct ConfTrackSelection : public o2::framework::ConfigurableGroup { o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum phi"}; o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; // track selection masks - o2::framework::Configurable maskLowMomentum{"maskLowMomentum", 0x2u, "Bitmask for selections below momentum threshold"}; - o2::framework::Configurable maskHighMomentum{"maskHighMomentum", 0x1u, "Bitmask for selections above momentum threshold"}; + o2::framework::Configurable maskLowMomentum{"maskLowMomentum", 1ul, "Bitmask for selections below momentum threshold"}; + o2::framework::Configurable maskHighMomentum{"maskHighMomentum", 2ul, "Bitmask for selections above momentum threshold"}; + // track rejection masks + o2::framework::Configurable rejectionMaskLowMomentum{"rejectionMaskLowMomentum", 0ul, "Bitmask for rejections below momentum threshold"}; + o2::framework::Configurable rejectionMaskHighMomentum{"rejectionMaskHighMomentum", 0ul, "Bitmask for rejections above momentum threshold"}; // momentum threshold for PID usage o2::framework::Configurable pidThres{"pidThres", 1.2f, "Momentum threshold for using TPCTOF/TOF pid for tracks with large momentum (GeV/c)"}; }; @@ -283,12 +286,12 @@ class TrackSelection : public BaseSelection void configure(o2::framework::HistogramRegistry* registry, T1& config, T2& filter) { - mPtMin = filter.ptMin; - mPtMax = filter.ptMax; - mEtaMin = filter.etaMin; - mEtaMax = filter.etaMax; - mPhiMin = filter.phiMin; - mPhiMax = filter.phiMax; + mPtMin = filter.ptMin.value; + mPtMax = filter.ptMax.value; + mEtaMin = filter.etaMin.value; + mEtaMax = filter.etaMax.value; + mPhiMin = filter.phiMin.value; + mPhiMax = filter.phiMax.value; // add selections for track quality this->addSelection(kTPCnClsMin, trackSelectionNames.at(kTPCnClsMin), config.tpcClustersMin.value, limits::kLowerLimit, true, true, false); @@ -298,61 +301,61 @@ class TrackSelection : public BaseSelectionaddSelection(kTPCsClsFracMax, trackSelectionNames.at(kTPCsClsFracMax), config.tpcSharedClusterFractionMax.value, limits::kUpperLimit, true, true, false); this->addSelection(kITSnClsMin, trackSelectionNames.at(kITSnClsMin), config.itsClustersMin.value, limits::kLowerLimit, true, true, false); this->addSelection(kITSnClsIbMin, trackSelectionNames.at(kITSnClsIbMin), config.itsIbClustersMin.value, limits::kLowerLimit, true, true, false); - this->addSelection(kDCAxyMax, trackSelectionNames.at(kDCAxyMax), filter.ptMin, filter.ptMax.value, config.dcaxyMax.value, limits::kAbsUpperFunctionLimit, true, true, false); + this->addSelection(kDCAxyMax, trackSelectionNames.at(kDCAxyMax), filter.ptMin.value, filter.ptMax.value, config.dcaxyMax.value, limits::kAbsUpperFunctionLimit, true, true, false); this->addSelection(kDCAzMax, trackSelectionNames.at(kDCAzMax), filter.ptMin.value, filter.ptMax.value, config.dcazMax.value, limits::kAbsUpperFunctionLimit, true, true, false); // add selections for Electron pid - this->addSelection(kItsElectron, trackSelectionNames.at(kItsElectron), config.itsElectron.value, limits::kAbsUpperLimit, false, false, config.requirePidElectron); - this->addSelection(kTpcElectron, trackSelectionNames.at(kTpcElectron), config.tpcElectron.value, limits::kAbsUpperLimit, false, false, config.requirePidElectron); - this->addSelection(kTofElectron, trackSelectionNames.at(kTofElectron), config.tofElectron.value, limits::kAbsUpperLimit, false, false, config.requirePidElectron); + this->addSelection(kItsElectron, trackSelectionNames.at(kItsElectron), config.itsElectron.value, false, false, config.requirePidElectron); + this->addSelection(kTpcElectron, trackSelectionNames.at(kTpcElectron), config.tpcElectron.value, false, false, config.requirePidElectron); + this->addSelection(kTofElectron, trackSelectionNames.at(kTofElectron), config.tofElectron.value, false, false, config.requirePidElectron); this->addSelection(kTpcitsElectron, trackSelectionNames.at(kTpcitsElectron), config.tpcitsElectron.value, limits::kUpperLimit, false, false, config.requirePidElectron); this->addSelection(kTpctofElectron, trackSelectionNames.at(kTpctofElectron), config.tpctofElectron.value, limits::kUpperLimit, false, false, config.requirePidElectron); mElectronTofThres = config.minMomTofElectron.value; // add selections for Pion pid - this->addSelection(kItsPion, trackSelectionNames.at(kItsPion), config.itsPion.value, limits::kAbsUpperLimit, false, false, config.requirePidPion); - this->addSelection(kTpcPion, trackSelectionNames.at(kTpcPion), config.tpcPion.value, limits::kAbsUpperLimit, false, false, config.requirePidPion); - this->addSelection(kTofPion, trackSelectionNames.at(kTofPion), config.tofPion.value, limits::kAbsUpperLimit, false, false, config.requirePidPion); + this->addSelection(kItsPion, trackSelectionNames.at(kItsPion), config.itsPion.value, false, false, config.requirePidPion); + this->addSelection(kTpcPion, trackSelectionNames.at(kTpcPion), config.tpcPion.value, false, false, config.requirePidPion); + this->addSelection(kTofPion, trackSelectionNames.at(kTofPion), config.tofPion.value, false, false, config.requirePidPion); this->addSelection(kTpcitsPion, trackSelectionNames.at(kTpcitsPion), config.tpcitsPion.value, limits::kUpperLimit, false, false, config.requirePidPion); this->addSelection(kTpctofPion, trackSelectionNames.at(kTpctofPion), config.tpctofPion.value, limits::kUpperLimit, false, false, config.requirePidPion); mPionTofThres = config.minMomTofPion.value; // add selections for Kaon pid - this->addSelection(kItsKaon, trackSelectionNames.at(kItsKaon), config.itsKaon.value, limits::kAbsUpperLimit, false, false, config.requirePidKaon); - this->addSelection(kTpcKaon, trackSelectionNames.at(kTpcKaon), config.tpcKaon.value, limits::kAbsUpperLimit, false, false, config.requirePidKaon); - this->addSelection(kTofKaon, trackSelectionNames.at(kTofKaon), config.tofKaon.value, limits::kAbsUpperLimit, false, false, config.requirePidKaon); + this->addSelection(kItsKaon, trackSelectionNames.at(kItsKaon), config.itsKaon.value, false, false, config.requirePidKaon); + this->addSelection(kTpcKaon, trackSelectionNames.at(kTpcKaon), config.tpcKaon.value, false, false, config.requirePidKaon); + this->addSelection(kTofKaon, trackSelectionNames.at(kTofKaon), config.tofKaon.value, false, false, config.requirePidKaon); this->addSelection(kTpcitsKaon, trackSelectionNames.at(kTpcitsKaon), config.tpcitsKaon.value, limits::kUpperLimit, false, false, config.requirePidKaon); this->addSelection(kTpctofKaon, trackSelectionNames.at(kTpctofKaon), config.tpctofKaon.value, limits::kUpperLimit, false, false, config.requirePidKaon); mKaonTofThres = config.minMomTofKaon.value; // add selections for Proton pid - this->addSelection(kItsProton, trackSelectionNames.at(kItsProton), config.itsProton.value, limits::kAbsUpperLimit, false, false, config.requirePidProton); - this->addSelection(kTpcProton, trackSelectionNames.at(kTpcProton), config.tpcProton.value, limits::kAbsUpperLimit, false, false, config.requirePidProton); - this->addSelection(kTofProton, trackSelectionNames.at(kTofProton), config.tofProton.value, limits::kAbsUpperLimit, false, false, config.requirePidProton); + this->addSelection(kItsProton, trackSelectionNames.at(kItsProton), config.itsProton.value, false, false, config.requirePidProton); + this->addSelection(kTpcProton, trackSelectionNames.at(kTpcProton), config.tpcProton.value, false, false, config.requirePidProton); + this->addSelection(kTofProton, trackSelectionNames.at(kTofProton), config.tofProton.value, false, false, config.requirePidProton); this->addSelection(kTpcitsProton, trackSelectionNames.at(kTpcitsProton), config.tpcitsProton.value, limits::kUpperLimit, false, false, config.requirePidProton); this->addSelection(kTpctofProton, trackSelectionNames.at(kTpctofProton), config.tpctofProton.value, limits::kUpperLimit, false, false, config.requirePidProton); mProtonTofThres = config.minMomTofProton.value; // add selections for Deuteron pid - this->addSelection(kItsDeuteron, trackSelectionNames.at(kItsDeuteron), config.itsDeuteron.value, limits::kAbsUpperLimit, false, false, config.requirePidDeuteron); - this->addSelection(kTpcDeuteron, trackSelectionNames.at(kTpcDeuteron), config.tpcDeuteron.value, limits::kAbsUpperLimit, false, false, config.requirePidDeuteron); - this->addSelection(kTofDeuteron, trackSelectionNames.at(kTofDeuteron), config.tofDeuteron.value, limits::kAbsUpperLimit, false, false, config.requirePidDeuteron); + this->addSelection(kItsDeuteron, trackSelectionNames.at(kItsDeuteron), config.itsDeuteron.value, false, false, config.requirePidDeuteron); + this->addSelection(kTpcDeuteron, trackSelectionNames.at(kTpcDeuteron), config.tpcDeuteron.value, false, false, config.requirePidDeuteron); + this->addSelection(kTofDeuteron, trackSelectionNames.at(kTofDeuteron), config.tofDeuteron.value, false, false, config.requirePidDeuteron); this->addSelection(kTpcitsDeuteron, trackSelectionNames.at(kTpcitsDeuteron), config.tpcitsDeuteron.value, limits::kUpperLimit, false, false, config.requirePidDeuteron); this->addSelection(kTpctofDeuteron, trackSelectionNames.at(kTpctofDeuteron), config.tpctofDeuteron.value, limits::kUpperLimit, false, false, config.requirePidDeuteron); mDeuteronTofThres = config.minMomTofDeuteron.value; // add selections for Triton pid - this->addSelection(kItsTriton, trackSelectionNames.at(kItsTriton), config.itsTriton.value, limits::kAbsUpperLimit, false, false, config.requirePidTriton); - this->addSelection(kTpcTriton, trackSelectionNames.at(kTpcTriton), config.tpcTriton.value, limits::kAbsUpperLimit, false, false, config.requirePidTriton); - this->addSelection(kTofTriton, trackSelectionNames.at(kTofTriton), config.tofTriton.value, limits::kAbsUpperLimit, false, false, config.requirePidTriton); + this->addSelection(kItsTriton, trackSelectionNames.at(kItsTriton), config.itsTriton.value, false, false, config.requirePidTriton); + this->addSelection(kTpcTriton, trackSelectionNames.at(kTpcTriton), config.tpcTriton.value, false, false, config.requirePidTriton); + this->addSelection(kTofTriton, trackSelectionNames.at(kTofTriton), config.tofTriton.value, false, false, config.requirePidTriton); this->addSelection(kTpcitsTriton, trackSelectionNames.at(kTpcitsTriton), config.tpcitsTriton.value, limits::kUpperLimit, false, false, config.requirePidTriton); this->addSelection(kTpctofTriton, trackSelectionNames.at(kTpctofTriton), config.tpctofTriton.value, limits::kUpperLimit, false, false, config.requirePidTriton); mTritonTofThres = config.minMomTofTriton.value; // add selections for Helium pid - this->addSelection(kItsHelium, trackSelectionNames.at(kItsHelium), config.itsHelium.value, limits::kAbsUpperLimit, false, false, config.requirePidHelium); - this->addSelection(kTpcHelium, trackSelectionNames.at(kTpcHelium), config.tpcHelium.value, limits::kAbsUpperLimit, false, false, config.requirePidHelium); - this->addSelection(kTofHelium, trackSelectionNames.at(kTofHelium), config.tofHelium.value, limits::kAbsUpperLimit, false, false, config.requirePidHelium); + this->addSelection(kItsHelium, trackSelectionNames.at(kItsHelium), config.itsHelium.value, false, false, config.requirePidHelium); + this->addSelection(kTpcHelium, trackSelectionNames.at(kTpcHelium), config.tpcHelium.value, false, false, config.requirePidHelium); + this->addSelection(kTofHelium, trackSelectionNames.at(kTofHelium), config.tofHelium.value, false, false, config.requirePidHelium); this->addSelection(kTpcitsHelium, trackSelectionNames.at(kTpcitsHelium), config.tpcitsHelium.value, limits::kUpperLimit, false, false, config.requirePidHelium); this->addSelection(kTpctofHelium, trackSelectionNames.at(kTpctofHelium), config.tpctofHelium.value, limits::kUpperLimit, false, false, config.requirePidHelium); mHeliumTofThres = config.minMomTofHelium.value; @@ -378,30 +381,20 @@ class TrackSelection : public BaseSelectionevaluateObservable(its, nsigmaIts); - this->evaluateObservable(tpc, nsigmaTpc); - this->evaluateObservable(tpcits, std::hypot(nsigmaTpc, nsigmaIts)); - this->evaluateObservable(tof, nsigmaTof); - this->evaluateObservable(tpctof, std::hypot(nsigmaTpc, nsigmaTof)); + // above threshold without TOF: skip entirely unless forced + // forced evaluation is used in the second pass to populate rejection bits + if (!ignoreThreshold && Track.p() >= tofThreshold && !Track.hasTOF()) { return; } - // if track is above threshold, check if TOF PID is available - // if not, we dont check any selection and they stay at reseted values, i.e. the cut fails + this->evaluateObservable(its, nsigmaIts); + this->evaluateObservable(tpc, nsigmaTpc); + this->evaluateObservable(tpcits, std::hypot(nsigmaTpc, nsigmaIts)); if (Track.hasTOF()) { - // if tof inforamtion is available, check them first this->evaluateObservable(tof, nsigmaTof); this->evaluateObservable(tpctof, std::hypot(nsigmaTpc, nsigmaTof)); - // if both failed, the bitmask will be 0 and there is no need to check tpc and its information since we do not want to have this track - // so if we just bail out here, the PID for this particle type will failed for its, tpc and tof - if (this->passesOptionalSelection(tof) || this->passesOptionalSelection(tpctof)) { - this->evaluateObservable(its, nsigmaIts); - this->evaluateObservable(tpc, nsigmaTpc); - this->evaluateObservable(tpcits, std::hypot(nsigmaTpc, nsigmaIts)); - } } } @@ -416,90 +409,34 @@ class TrackSelection : public BaseSelectionevaluateObservable(kTPCsClsFracMax, static_cast(Track.tpcNClsShared()) / static_cast(Track.tpcNClsFound())); this->evaluateObservable(kITSnClsMin, Track.itsNCls()); this->evaluateObservable(kITSnClsIbMin, Track.itsNClsInnerBarrel()); - - // evalue bitmask for pt dependent dca cuts this->updateLimits(kDCAxyMax, Track.pt()); this->evaluateObservable(kDCAxyMax, Track.dcaXY()); - this->updateLimits(kDCAzMax, Track.pt()); this->evaluateObservable(kDCAzMax, Track.dcaZ()); - this->evaluatePid(Track, - mElectronTofThres, - Track.itsNSigmaEl(), - Track.tpcNSigmaEl(), - Track.tofNSigmaEl(), - kItsElectron, - kTpcElectron, - kTofElectron, - kTpcitsElectron, - kTpctofElectron); - - this->evaluatePid(Track, - mPionTofThres, - Track.itsNSigmaPi(), - Track.tpcNSigmaPi(), - Track.tofNSigmaPi(), - kItsPion, - kTpcPion, - kTofPion, - kTpcitsPion, - kTpctofPion); - - this->evaluatePid(Track, - mKaonTofThres, - Track.itsNSigmaKa(), - Track.tpcNSigmaKa(), - Track.tofNSigmaKa(), - kItsKaon, - kTpcKaon, - kTofKaon, - kTpcitsKaon, - kTpctofKaon); - - this->evaluatePid(Track, - mProtonTofThres, - Track.itsNSigmaPr(), - Track.tpcNSigmaPr(), - Track.tofNSigmaPr(), - kItsProton, - kTpcProton, - kTofProton, - kTpcitsProton, - kTpctofProton); - - this->evaluatePid(Track, - mDeuteronTofThres, - Track.itsNSigmaDe(), - Track.tpcNSigmaDe(), - Track.tofNSigmaDe(), - kItsDeuteron, - kTpcDeuteron, - kTofDeuteron, - kTpcitsDeuteron, - kTpctofDeuteron); - - this->evaluatePid(Track, - mTritonTofThres, - Track.itsNSigmaTr(), - Track.tpcNSigmaTr(), - Track.tofNSigmaTr(), - kItsTriton, - kTpcTriton, - kTofTriton, - kTpcitsTriton, - kTpctofTriton); - - this->evaluatePid(Track, - mHeliumTofThres, - Track.itsNSigmaHe(), - Track.tpcNSigmaHe(), - Track.tofNSigmaHe(), - kItsHelium, - kTpcHelium, - kTofHelium, - kTpcitsHelium, - kTpctofHelium); + // first pass: threshold-aware PID evaluation + // determines if the track passes any optional selection and should be stored + this->evaluatePid(Track, mElectronTofThres, Track.itsNSigmaEl(), Track.tpcNSigmaEl(), Track.tofNSigmaEl(), kItsElectron, kTpcElectron, kTofElectron, kTpcitsElectron, kTpctofElectron); + this->evaluatePid(Track, mPionTofThres, Track.itsNSigmaPi(), Track.tpcNSigmaPi(), Track.tofNSigmaPi(), kItsPion, kTpcPion, kTofPion, kTpcitsPion, kTpctofPion); + this->evaluatePid(Track, mKaonTofThres, Track.itsNSigmaKa(), Track.tpcNSigmaKa(), Track.tofNSigmaKa(), kItsKaon, kTpcKaon, kTofKaon, kTpcitsKaon, kTpctofKaon); + this->evaluatePid(Track, mProtonTofThres, Track.itsNSigmaPr(), Track.tpcNSigmaPr(), Track.tofNSigmaPr(), kItsProton, kTpcProton, kTofProton, kTpcitsProton, kTpctofProton); + this->evaluatePid(Track, mDeuteronTofThres, Track.itsNSigmaDe(), Track.tpcNSigmaDe(), Track.tofNSigmaDe(), kItsDeuteron, kTpcDeuteron, kTofDeuteron, kTpcitsDeuteron, kTpctofDeuteron); + this->evaluatePid(Track, mTritonTofThres, Track.itsNSigmaTr(), Track.tpcNSigmaTr(), Track.tofNSigmaTr(), kItsTriton, kTpcTriton, kTofTriton, kTpcitsTriton, kTpctofTriton); + this->evaluatePid(Track, mHeliumTofThres, Track.itsNSigmaHe(), Track.tpcNSigmaHe(), Track.tofNSigmaHe(), kItsHelium, kTpcHelium, kTofHelium, kTpcitsHelium, kTpctofHelium); + + // second pass: if the track passes minimal and any optional selection, + // re-evaluate all species ignoring thresholds so rejection bits are fully + // populated for all competing hypotheses. evaluate() resets each container + // before writing, so no explicit reset is needed before this pass. + if (this->passesAllRequiredSelections()) { + this->evaluatePid(Track, mElectronTofThres, Track.itsNSigmaEl(), Track.tpcNSigmaEl(), Track.tofNSigmaEl(), kItsElectron, kTpcElectron, kTofElectron, kTpcitsElectron, kTpctofElectron, true); + this->evaluatePid(Track, mPionTofThres, Track.itsNSigmaPi(), Track.tpcNSigmaPi(), Track.tofNSigmaPi(), kItsPion, kTpcPion, kTofPion, kTpcitsPion, kTpctofPion, true); + this->evaluatePid(Track, mKaonTofThres, Track.itsNSigmaKa(), Track.tpcNSigmaKa(), Track.tofNSigmaKa(), kItsKaon, kTpcKaon, kTofKaon, kTpcitsKaon, kTpctofKaon, true); + this->evaluatePid(Track, mProtonTofThres, Track.itsNSigmaPr(), Track.tpcNSigmaPr(), Track.tofNSigmaPr(), kItsProton, kTpcProton, kTofProton, kTpcitsProton, kTpctofProton, true); + this->evaluatePid(Track, mDeuteronTofThres, Track.itsNSigmaDe(), Track.tpcNSigmaDe(), Track.tofNSigmaDe(), kItsDeuteron, kTpcDeuteron, kTofDeuteron, kTpcitsDeuteron, kTpctofDeuteron, true); + this->evaluatePid(Track, mTritonTofThres, Track.itsNSigmaTr(), Track.tpcNSigmaTr(), Track.tofNSigmaTr(), kItsTriton, kTpcTriton, kTofTriton, kTpcitsTriton, kTpctofTriton, true); + this->evaluatePid(Track, mHeliumTofThres, Track.itsNSigmaHe(), Track.tpcNSigmaHe(), Track.tofNSigmaHe(), kItsHelium, kTpcHelium, kTofHelium, kTpcitsHelium, kTpctofHelium, true); + } this->assembleBitmask(); }