diff --git a/PWGCF/Femto/Core/baseSelection.h b/PWGCF/Femto/Core/baseSelection.h index 5e5b0deab97..3e58511914b 100644 --- a/PWGCF/Femto/Core/baseSelection.h +++ b/PWGCF/Femto/Core/baseSelection.h @@ -78,7 +78,7 @@ class BaseSelection { // check index if (static_cast(observableIndex) >= NumObservables) { - LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables; + LOG(fatal) << "Selection '" << selectionName << "': observable is not valid. Observable (index) has to be smaller than " << NumObservables; } // init selection container for selection at given index mSelectionContainers.at(observableIndex) = selectioncontainer::SelectionContainer(selectionName, @@ -112,7 +112,7 @@ class BaseSelection bool isOptionalCut) { if (static_cast(observableIndex) >= NumObservables) { - LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables; + LOG(fatal) << "Selection '" << selectionName << "': observable is not valid. Observable (index) has to be smaller than " << NumObservables; } mSelectionContainers.at(observableIndex) = selectioncontainer::SelectionContainer(selectionName, lowerLimit, @@ -146,7 +146,7 @@ class BaseSelection { // check index if (static_cast(observableIndex) >= NumObservables) { - LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables; + LOG(fatal) << "Selection '" << selectionName << "': observable is not valid. Observable (index) has to be smaller than " << NumObservables; } // init selection container for selection at given index mSelectionContainers.at(observableIndex) = selectioncontainer::SelectionContainer(selectionName, @@ -191,7 +191,7 @@ class BaseSelection mSelectionContainers.at(observableIndex) = selectioncontainer::SelectionContainer(selectionName, std::vector{1}, limits::LimitType::kEqual, false, false, false); break; default: - LOG(fatal) << "Invalid switch for boolean selection"; + LOG(fatal) << "Selection '" << selectionName << "': invalid switch for boolean selection"; } addSelection(observableIndex); } @@ -577,7 +577,7 @@ class BaseSelection } if (mNSelectionBits > sizeof(BitmaskType) * CHAR_BIT) { - LOG(fatal) << "Too many selections. At most " << sizeof(BitmaskType) * CHAR_BIT << " number of bits are supported"; + LOG(fatal) << "Selection '" << mSelectionContainers.at(observableIndex).getSelectionName() << "': too many selections. At most " << sizeof(BitmaskType) * CHAR_BIT << " number of bits are supported"; } } diff --git a/PWGCF/Femto/Core/selectionContainer.h b/PWGCF/Femto/Core/selectionContainer.h index 351ee34b674..0017215b718 100644 --- a/PWGCF/Femto/Core/selectionContainer.h +++ b/PWGCF/Femto/Core/selectionContainer.h @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -122,10 +123,10 @@ class SelectionContainer mIsOptionalCut(isOptionalCut) { if (mSelectionValues.size() > sizeof(BitmaskType) * CHAR_BIT) { - LOG(fatal) << "Too many selections for single a observable. Limit is " << sizeof(BitmaskType) * CHAR_BIT; + LOG(fatal) << "Selection '" << mSelectionName << "': 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"; + LOG(fatal) << "Selection '" << mSelectionName << "': a selection cannot be both minimal and optional"; } // values for selection are not necessarily ordered correctly sortSelections(); @@ -157,13 +158,13 @@ class SelectionContainer mIsOptionalCut(isOptionalCut) { if (mLimitType != limits::kInsideRange && mLimitType != limits::kOutsideRange) { - LOG(fatal) << "Range-string selections only support kInsideRange or kOutsideRange"; + LOG(fatal) << "Selection '" << mSelectionName << "': range-string selections only support kInsideRange or kOutsideRange"; } if (rangeStrings.size() > sizeof(BitmaskType) * CHAR_BIT) { - LOG(fatal) << "Too many selections for a single observable. Limit is " << sizeof(BitmaskType) * CHAR_BIT; + LOG(fatal) << "Selection '" << mSelectionName << "': 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"; + LOG(fatal) << "Selection '" << mSelectionName << "': a selection cannot be both minimal and optional"; } for (auto const& rangeStr : rangeStrings) { T lower, upper; @@ -196,7 +197,7 @@ class SelectionContainer : (mSelectionRanges[i].first <= mSelectionRanges[i + 1].first && mSelectionRanges[i].second >= mSelectionRanges[i + 1].second); if (!nested) { - LOG(fatal) << "Ranges for minimal cut " << selectionName + LOG(fatal) << "Selection '" << selectionName << "': ranges for minimal cut" << " are not nested. Range [" << mSelectionRanges[i].first << ";" << mSelectionRanges[i].second << "] is not " << (mLimitType == limits::kOutsideRange ? "contained in" : "containing") << " [" << mSelectionRanges[i + 1].first << ";" << mSelectionRanges[i + 1].second << "]"; @@ -229,10 +230,10 @@ class SelectionContainer mIsOptionalCut(isOptionalCut) { if (functions.size() > sizeof(BitmaskType) * CHAR_BIT) { - LOG(fatal) << "Too many selections for single a observable. Limit is " << sizeof(BitmaskType) * CHAR_BIT; + LOG(fatal) << "Selection '" << mSelectionName << "': 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"; + LOG(fatal) << "Selection '" << mSelectionName << "': 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); @@ -255,7 +256,7 @@ class SelectionContainer // 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"; + LOG(fatal) << "Selection '" << mSelectionName << "': number of comments and number of selections are inconsistent"; } mComments = comments; } @@ -382,7 +383,7 @@ class SelectionContainer void evaluate(std::vector const& values) { if (values.size() != mSelectionValues.size()) { - LOG(fatal) << "Wrong number of values have been passed"; + LOG(fatal) << "Selection '" << mSelectionName << "': wrong number of values have been passed"; } // reset the bitmask before evaluating mBitmask.reset(); @@ -462,7 +463,7 @@ class SelectionContainer T getLoosestSelection() const { if (mSelectionValues.empty()) { - LOG(fatal) << "No selections configured"; + LOG(fatal) << "Selection '" << mSelectionName << "': no selections configured"; } return mSelectionValues.at(0); } @@ -605,7 +606,7 @@ class SelectionContainer [[nodiscard]] int getBitPosition(int selectionIndex) const { if (selectionIndex == 0 && mSkipMostPermissiveBit) { - LOG(fatal) << "Trying to accessed the bit position of a skipped selection. Breaking..."; + LOG(fatal) << "Selection '" << mSelectionName << "': trying to access the bit position of a skipped selection. Breaking..."; return -1; } if (mSkipMostPermissiveBit) { @@ -706,7 +707,10 @@ class SelectionContainer { auto pos = rangeStr.find(';'); if (pos == std::string::npos) { - LOG(fatal) << "Range string '" << rangeStr << "' is missing ';' separator. Expected format: 'lower;upper'"; + LOG(fatal) << "Selection '" << mSelectionName << "': range string '" << rangeStr << "' is missing ';' separator. Expected format: 'lower;upper'"; + } + if (rangeStr.find(';', pos + 1) != std::string::npos) { + LOG(fatal) << "Selection '" << mSelectionName << "': range string '" << rangeStr << "' has more than one ';' separator. Expected format: 'lower;upper'"; } std::string lowerStr = rangeStr.substr(0, pos); std::string upperStr = rangeStr.substr(pos + 1); @@ -723,14 +727,36 @@ class SelectionContainer return; } - lower = lowerStr.empty() ? -std::numeric_limits::max() : static_cast(std::stod(lowerStr)); - upper = upperStr.empty() ? std::numeric_limits::max() : static_cast(std::stod(upperStr)); + lower = lowerStr.empty() ? -std::numeric_limits::max() : parseBound(rangeStr, lowerStr); + upper = upperStr.empty() ? std::numeric_limits::max() : parseBound(rangeStr, upperStr); if (lower >= upper) { - LOG(fatal) << "Range string '" << rangeStr << "' has lower bound >= upper bound"; + LOG(fatal) << "Selection '" << mSelectionName << "': range string '" << rangeStr << "' has lower bound >= upper bound"; } } + /// \brief Parse a single range bound and make sure the whole string was consumed, so that + /// malformed input such as "3;5" (extra ';') or "-3,5" (comma used as decimal + /// separator) is rejected instead of silently truncated by std::stod (which happily + /// parses just the numeric prefix "3" or "-3" and ignores the rest). + /// \param rangeStr Full "lower;upper" string, only used for the error message. + /// \param boundStr The individual bound substring (lower or upper) to parse. + T parseBound(std::string const& rangeStr, std::string const& boundStr) const + { + std::size_t idx = 0; + double value = 0.; + try { + value = std::stod(boundStr, &idx); + } catch (std::exception const& e) { + LOG(fatal) << "Selection '" << mSelectionName << "': range string '" << rangeStr << "' contains bound '" << boundStr << "' that is not a valid number (" << e.what() << ")"; + } + if (idx != boundStr.size()) { + LOG(fatal) << "Selection '" << mSelectionName << "': range string '" << rangeStr << "' contains bound '" << boundStr + << "' with trailing characters ('" << boundStr.substr(idx) << "') that could not be parsed"; + } + return static_cast(value); + } + 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) diff --git a/PWGCF/Femto/Core/twoTrackResonanceBuilder.h b/PWGCF/Femto/Core/twoTrackResonanceBuilder.h index fece3d0e06f..26228021eea 100644 --- a/PWGCF/Femto/Core/twoTrackResonanceBuilder.h +++ b/PWGCF/Femto/Core/twoTrackResonanceBuilder.h @@ -33,6 +33,7 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include #include namespace o2::analysis::femto::twotrackresonancebuilder @@ -341,5 +342,120 @@ class TwoTrackResonanceBuilder }; // namespace twotrackresonancebuilder +struct ConfTwoTrackResonanceTablesDerivedToDerived : o2::framework::ConfigurableGroup { + std::string prefix = std::string("TwoTrackResonanceTables"); + o2::framework::Configurable limitRho0{"limitRho0", 1, "Require at least this many Rho0s in the collision. Set to 0 to skip Rho0s entirely (not written to output)."}; + o2::framework::Configurable limitPhi{"limitPhi", 0, "Require at least this many Phis in the collision. Set to 0 to skip Phis entirely (not written to output)."}; +}; + +struct TwoTrackResonanceBuilderDerivedToDerivedProducts : o2::framework::ProducesGroup { + o2::framework::Produces producedRho0s; + o2::framework::Produces producedRho0Masks; + o2::framework::Produces producedPhis; + o2::framework::Produces producedPhiMasks; +}; + +class TwoTrackResonanceBuilderDerivedToDerived +{ + public: + template + void init(T& config) + { + mLimitRho0 = config.limitRho0.value; + mLimitPhi = config.limitPhi.value; + + if (mLimitRho0 < 0 || mLimitPhi < 0) { + LOG(fatal) << "Two-track resonance limits must be non-negative (got " << mLimitRho0 << " and " << mLimitPhi << "). Breaking..."; + } + if (mLimitRho0 == 0 && mLimitPhi == 0) { + LOG(fatal) << "Both Rho0 limit and Phi limit are 0. Breaking..."; + } + } + + template + bool collisionHasTooFewRho0s(T1 const& col, T2 const& /*rho0Table*/, T3& partitionRho0, T4& cache) + { + if (mLimitRho0 == 0) { // Rho0s disabled, cannot reject on them + return false; + } + auto rho0Slice = partitionRho0->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + return rho0Slice.size() < static_cast(mLimitRho0); + } + + template + bool collisionHasTooFewPhis(T1 const& col, T2 const& /*phiTable*/, T3& partitionPhi, T4& cache) + { + if (mLimitPhi == 0) { // Phis disabled, cannot reject on them + return false; + } + auto phiSlice = partitionPhi->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + return phiSlice.size() < static_cast(mLimitPhi); + } + + template + void processRho0s(T1 const& col, T2 const& /*rho0Table*/, T3 const& oldTrackTable, T4& partitionRho0, T5& trackBuilder, T6& cache, T7& newRho0Table, T8& newTrackTable, T9& newCollisionTable) + { + if (mLimitRho0 == 0) { // Rho0s disabled + return; + } + + auto rho0Slice = partitionRho0->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + + for (auto const& rho0 : rho0Slice) { + auto posDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(rho0.posDauId(), oldTrackTable)); + auto negDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(rho0.negDauId(), oldTrackTable)); + + int64_t posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable); + int64_t negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable); + + newRho0Table.producedRho0s(newCollisionTable.producedCollision.lastIndex(), + rho0.pt(), + rho0.eta(), + rho0.phi(), + rho0.mass(), + posDaughterIndex, + negDaughterIndex); + newRho0Table.producedRho0Masks(rho0.maskPosDau(), + rho0.posDauHasHighMomentum(), + rho0.maskNegDau(), + rho0.negDauHasHighMomentum()); + } + } + + template + void processPhis(T1 const& col, T2 const& /*phiTable*/, T3 const& oldTrackTable, T4& partitionPhi, T5& trackBuilder, T6& cache, T7& newPhiTable, T8& newTrackTable, T9& newCollisionTable) + { + if (mLimitPhi == 0) { // Phis disabled + return; + } + + auto phiSlice = partitionPhi->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + + for (auto const& phi : phiSlice) { + auto posDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(phi.posDauId(), oldTrackTable)); + auto negDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(phi.negDauId(), oldTrackTable)); + + int64_t posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable); + int64_t negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable); + + newPhiTable.producedPhis(newCollisionTable.producedCollision.lastIndex(), + phi.pt(), + phi.eta(), + phi.phi(), + phi.mass(), + posDaughterIndex, + negDaughterIndex); + newPhiTable.producedPhiMasks(phi.maskPosDau(), + phi.posDauHasHighMomentum(), + phi.maskNegDau(), + phi.negDauHasHighMomentum()); + } + } + + private: + int mLimitRho0 = 0; + int mLimitPhi = 0; +}; + } // namespace o2::analysis::femto::twotrackresonancebuilder #endif // PWGCF_FEMTO_CORE_TWOTRACKRESONANCEBUILDER_H_ diff --git a/PWGCF/Femto/DataModel/FemtoTables.h b/PWGCF/Femto/DataModel/FemtoTables.h index 029e5d3632e..60e4527d93a 100644 --- a/PWGCF/Femto/DataModel/FemtoTables.h +++ b/PWGCF/Femto/DataModel/FemtoTables.h @@ -552,12 +552,14 @@ DECLARE_SOA_TABLE_STAGED_VERSIONED(FPhis_001, "FPHI", 1, //! femto phis femtobase::dynamic::Pz, femtobase::dynamic::Theta); using FPhis = FPhis_001; +using StoredFPhis = StoredFPhis_001; DECLARE_SOA_TABLE_STAGED_VERSIONED(FPhiMasks_001, "FPHIMASK", 1, //! mask for phis femtotwotrackresonances::MaskPosDau, femtotwotrackresonances::PosDauHasHighMomentum, femtotwotrackresonances::MaskNegDau, femtotwotrackresonances::NegDauHasHighMomentum); using FPhiMasks = FPhiMasks_001; +using StoredFPhiMasks = StoredFPhiMasks_001; // table for kstars DECLARE_SOA_TABLE_STAGED_VERSIONED(FKstar0s_001, "FKSTAR0", 1, //! femto k0star @@ -599,12 +601,14 @@ DECLARE_SOA_TABLE_STAGED_VERSIONED(FRho0s_001, "FRHO0", 1, //! femto rho0s femtobase::dynamic::Pz, femtobase::dynamic::Theta); using FRho0s = FRho0s_001; +using StoredFRho0s = StoredFRho0s_001; DECLARE_SOA_TABLE_STAGED_VERSIONED(FRho0Masks_001, "FRHO0MASK", 1, //! rho0s masks femtotwotrackresonances::MaskPosDau, femtotwotrackresonances::PosDauHasHighMomentum, femtotwotrackresonances::MaskNegDau, femtotwotrackresonances::NegDauHasHighMomentum); using FRho0Masks = FRho0Masks_001; +using StoredFRho0Masks = StoredFRho0Masks_001; namespace femtov0s { diff --git a/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx b/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx index 3ca85915f55..3ce6b84aca3 100644 --- a/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx +++ b/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx @@ -19,6 +19,7 @@ #include "PWGCF/Femto/Core/mcBuilder.h" #include "PWGCF/Femto/Core/partitions.h" #include "PWGCF/Femto/Core/trackBuilder.h" +#include "PWGCF/Femto/Core/twoTrackResonanceBuilder.h" #include "PWGCF/Femto/Core/v0Builder.h" #include "PWGCF/Femto/DataModel/FemtoTables.h" @@ -47,6 +48,8 @@ struct FemtoProducerDerivedToDerived { using FemtoTracksWithMass = o2::soa::Join; using FemtoLambdas = o2::soa::Join; using FemtoK0shorts = o2::soa::Join; + using FemtoRho0s = o2::soa::Join; + using FemtoPhis = o2::soa::Join; using FemtoXis = o2::soa::Join; using FemtoOmegas = o2::soa::Join; using FemtoSigma = o2::soa::Join; @@ -93,6 +96,19 @@ struct FemtoProducerDerivedToDerived { o2::framework::Partition k0shortPartition = MAKE_K0SHORT_PARTITION(k0shortSelection1); o2::framework::Preslice perColK0shorts = o2::aod::femtobase::stored::fColId; + // two-track resonance builder + twotrackresonancebuilder::TwoTrackResonanceBuilderDerivedToDerived resonanceBuilder; + twotrackresonancebuilder::TwoTrackResonanceBuilderDerivedToDerivedProducts resonanceBuilderProducts; + twotrackresonancebuilder::ConfTwoTrackResonanceTablesDerivedToDerived confResonanceBuilder; + + twotrackresonancebuilder::ConfRho0Selection rho0Selection; + o2::framework::Partition rho0Partition = MAKE_RESONANCE_0_PARTITON(rho0Selection); + o2::framework::Preslice perColRho0s = o2::aod::femtobase::stored::fColId; + + twotrackresonancebuilder::ConfPhiSelection phiSelection; + o2::framework::Partition phiPartition = MAKE_RESONANCE_0_PARTITON(phiSelection); + o2::framework::Preslice perColPhis = o2::aod::femtobase::stored::fColId; + // kink builder kinkbuilder::KinkBuilderDerivedToDerived kinkBuilder; kinkbuilder::KinkBuilderDerivedToDerivedProducts kinkBuilderProducts; @@ -143,6 +159,7 @@ struct FemtoProducerDerivedToDerived { static_cast(doprocessTracksXis) + static_cast(doprocessTracksOmegas) + static_cast(doprocessTracksK0shorts) + static_cast(doprocessTracksSigma) + static_cast(doprocessTracksSigmaPlus) + + static_cast(doprocessTracksRho0) + static_cast(doprocessTracksPhi) + static_cast(doprocessTracksMc) + static_cast(doprocessTracksSigmaMc) + static_cast(doprocessTracksSigmaPlusMc); if (activeProcesses != 1) { @@ -153,6 +170,7 @@ struct FemtoProducerDerivedToDerived { v0Builder.init(confV0Builder); cascadeBuilder.init(confCascadeBuilder); kinkBuilder.init(confKinkBuilder); + resonanceBuilder.init(confResonanceBuilder); mcBuilder.init(confMcBuilder); } @@ -274,6 +292,30 @@ struct FemtoProducerDerivedToDerived { } PROCESS_SWITCH(FemtoProducerDerivedToDerived, processTracksSigmaPlus, "Process sigmaPlus and tracks", false); + void processTracksRho0(FilteredFemtoCollision const& col, FemtoTracks const& tracks, FemtoRho0s const& rho0s) + { + if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache) || resonanceBuilder.collisionHasTooFewRho0s(col, rho0s, rho0Partition, cache)) { + return; + } + trackBuilder.reset(tracks); + collisionBuilder.processCollision(col, collisionBuilderProducts); + trackBuilder.processTracks(col, tracks, trackPartition1, trackPartition2, cache, trackBuilderProducts, collisionBuilderProducts); + resonanceBuilder.processRho0s(col, rho0s, tracks, rho0Partition, trackBuilder, cache, resonanceBuilderProducts, trackBuilderProducts, collisionBuilderProducts); + } + PROCESS_SWITCH(FemtoProducerDerivedToDerived, processTracksRho0, "Process rho0 and tracks", false); + + void processTracksPhi(FilteredFemtoCollision const& col, FemtoTracks const& tracks, FemtoPhis const& phis) + { + if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache) || resonanceBuilder.collisionHasTooFewPhis(col, phis, phiPartition, cache)) { + return; + } + trackBuilder.reset(tracks); + collisionBuilder.processCollision(col, collisionBuilderProducts); + trackBuilder.processTracks(col, tracks, trackPartition1, trackPartition2, cache, trackBuilderProducts, collisionBuilderProducts); + resonanceBuilder.processPhis(col, phis, tracks, phiPartition, trackBuilder, cache, resonanceBuilderProducts, trackBuilderProducts, collisionBuilderProducts); + } + PROCESS_SWITCH(FemtoProducerDerivedToDerived, processTracksPhi, "Process phi and tracks", false); + // mc processing // // NOTE: unlike the other process functions, this one subscribes to the full // collision table and iterates manually. The MC index maps are keyed by source