diff --git a/PWGCF/Femto/Core/baseSelection.h b/PWGCF/Femto/Core/baseSelection.h index c2280dc21eb..d85829beb7a 100644 --- a/PWGCF/Femto/Core/baseSelection.h +++ b/PWGCF/Femto/Core/baseSelection.h @@ -170,15 +170,46 @@ class BaseSelection // 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 (mSelectionContainers.at(observableIndex).passesAsMinimalCut() == false) { mPassesMinimalSelections = false; } // check if any optional selection holds + // if one optional selection is fullfilled, the condition succeeds if (mSelectionContainers.at(observableIndex).passesAsOptionalCut() == true) { mPassesOptionalSelections = true; } } + /// \brief Evaluate a single observable against its configured selections. + /// \param observableIndex Index 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 (mPassesMinimalSelections == false) { + return; + } + // set bitmask for given observable + mSelectionContainers.at(observableIndex).evaluate(values); + // check if minimal selction for this observable holds + if (mSelectionContainers.at(observableIndex).passesAsMinimalCut() == false) { + mPassesMinimalSelections = false; + } + // check if any optional selection holds + if (mSelectionContainers.at(observableIndex).passesAsOptionalCut() == true) { + mPassesOptionalSelections = true; + } + } + + /// \brief Add comments to specific observabel + 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. bool passesAllRequiredSelections() const @@ -223,20 +254,38 @@ class BaseSelection /// \return The combined selection bitmask. BitmaskType getBitmask(int observableIndex) const { return static_cast(mSelectionContainers.at(observableIndex).getBitmask().to_ullong()); } - /// \brief Retrieve the assembled bitmask as an integer value. + /// \brief Set the assembled bitmask for on observable /// \return The combined selection bitmask. template void setBitmask(int observableIndex, R bitmask) { + // 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 (mPassesMinimalSelections == false) { + return; + } + // set bitmask for given observable mSelectionContainers.at(observableIndex).setBitmask(bitmask); + // check if minimal selction for this observable holds + if (mSelectionContainers.at(observableIndex).passesAsMinimalCut() == false) { + mPassesMinimalSelections = false; + } + // check if any optional selection holds + if (mSelectionContainers.at(observableIndex).passesAsOptionalCut() == true) { + mPassesOptionalSelections = true; + } } /// \brief Print detailed information about all configured selections. /// \tparam MapType Type used in the observable name map (usually an enum or int). /// \param objectName Name of the current object (e.g. particle species). /// \param observableNames Map from observable index to human-readable names. - template - void printSelections(const std::string& objectName, const std::unordered_map& observableNames) const + template + void printSelections(const std::string& objectName, const std::unordered_map& observableNames) const { LOG(info) << "Printing Configuration of " << objectName; @@ -248,7 +297,7 @@ class BaseSelection continue; } - const MapType key = static_cast(idx); + R key = static_cast(idx); const std::string& name = observableNames.count(key) ? observableNames.at(key) : "[Unknown]"; LOG(info) << "Observable: " << name << " (index " << idx << ")"; @@ -260,9 +309,11 @@ class BaseSelection const auto& values = container.getSelectionValues(); const auto& functions = container.getSelectionFunction(); - const bool useFunctions = !functions.empty(); - const size_t numSelections = useFunctions ? functions.size() : values.size(); - const bool skipMostPermissive = container.skipMostPermissiveBit(); + bool useFunctions = !functions.empty(); + size_t numSelections = useFunctions ? functions.size() : values.size(); + bool skipMostPermissive = container.skipMostPermissiveBit(); + const auto& comments = container.getComments(); + bool hasComments = !comments.empty(); int valWidth = 20; int bitWidth = 30; @@ -273,7 +324,9 @@ class BaseSelection // Selection string (either value or function) const std::string& sel = useFunctions ? std::string(functions[j].GetFormula()->GetExpFormula().Data()) : std::to_string(values[j]); line << " " << std::left << std::setw(valWidth) << sel; - + if (hasComments) { + line << "(" << comments.at(j) << ")"; + } // Bitmask if (skipMostPermissive && j == 0) { line << std::setw(bitWidth) << "-> loosest minimal selection, no bit saved"; diff --git a/PWGCF/Femto/Core/collisionBuilder.h b/PWGCF/Femto/Core/collisionBuilder.h index adcfd40b5ee..71f9bc48028 100644 --- a/PWGCF/Femto/Core/collisionBuilder.h +++ b/PWGCF/Femto/Core/collisionBuilder.h @@ -25,14 +25,17 @@ #include "Common/CCDB/EventSelectionParams.h" #include "EventFiltering/Zorro.h" +#include "DataFormatsParameters/GRPMagField.h" #include "Framework/AnalysisHelpers.h" #include "Framework/Configurable.h" #include "fairlogger/Logger.h" +#include #include #include #include +#include namespace o2::analysis::femto { @@ -73,13 +76,14 @@ struct ConfCollisionBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> occupancyMax{"occupancyMax", {}, "Maximum occpancy"}; o2::framework::Configurable> sphericityMin{"sphericityMin", {}, "Minimum sphericity"}; o2::framework::Configurable> sphericityMax{"sphericityMax", {}, "Maximum sphericity"}; + o2::framework::Configurable> triggers{"triggers", {}, "List of all triggers to be used"}; }; -struct ConfCollisionTriggers : o2::framework::ConfigurableGroup { - std::string prefix = std::string("CollisionTriggers"); - o2::framework::Configurable useTrigger{"useTrigger", false, "Set to true to only selected triggered collisions"}; - o2::framework::Configurable ccdbPath{"ccdbPath", std::string("EventFiltering/Zorro/"), "CCDB path for trigger information"}; - o2::framework::Configurable triggers{"triggers", std::string("fPPP,fPPL"), "Comma seperated list of all triggers to be used"}; +struct ConfCcdb : o2::framework::ConfigurableGroup { + std::string prefix = std::string("ConfCcdb"); + o2::framework::Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "URL to ccdb"}; + o2::framework::Configurable grpPath{"grpPath", "GLO/Config/GRPMagField", "Path to GRP object (Run3 -> GLO/Config/GRPMagField/Run2 -> GLO/GRP/GRP"}; + o2::framework::Configurable triggerPath{"triggerPath", "EventFiltering/Zorro/", "CCDB path for trigger information"}; }; struct ConfCollisionRctFlags : o2::framework::ConfigurableGroup { @@ -125,6 +129,8 @@ enum CollisionSels { kSphericityMin, ///< Min. sphericity kSphericityMax, ///< Max. sphericity + kTriggers, + kCollisionSelsMax }; @@ -146,14 +152,16 @@ const std::unordered_map colSelsToString = { {kOccupancyMin, "Minimum Occupancy"}, {kOccupancyMax, "Maximum Occupancy"}, {kSphericityMin, "Minimum Sphericity"}, - {kSphericityMax, "Maximum Sphericity"} + {kSphericityMax, "Maximum Sphericity"}, + + {kTriggers, "Triggers"} }; class CollisionSelection : public BaseSelection { public: - CollisionSelection() {} + CollisionSelection() = default; virtual ~CollisionSelection() = default; template @@ -188,6 +196,10 @@ class CollisionSelection : public BaseSelectionaddSelection(config.occupancyMax.value, kOccupancyMax, limits::kUpperLimit, true, true); this->addSelection(config.sphericityMin.value, kSphericityMin, limits::kLowerLimit, true, true); this->addSelection(config.sphericityMax.value, kSphericityMax, limits::kUpperLimit, true, true); + + std::vector triggerValues(config.triggers.value.size(), 1.f); + this->addSelection(triggerValues, kTriggers, limits::kEqualArray, false, true); + this->addComments(kTriggers, config.triggers.value); }; void setMagneticField(int MagField) @@ -255,13 +267,12 @@ class CollisionSelection : public BaseSelection - void applySelections(T const& col) + void applySelections(T const& col, std::vector const& triggerDecisions) { this->reset(); - // casting bool to float gurantees - // false -> 0 - // true -> 1 + // casting bool to float gurantees false -> 0 and true -> 1 + // and we check for equality to 1, so evaluation succeeds if the selection bit is true this->evaluateObservable(kSel8, static_cast(col.sel8())); this->evaluateObservable(kNoSameBunchPileUp, static_cast(col.selection_bit(o2::aod::evsel::kNoSameBunchPileup))); this->evaluateObservable(kIsVertexItsTpc, static_cast(col.selection_bit(o2::aod::evsel::kIsVertexITSTPC))); @@ -280,6 +291,13 @@ class CollisionSelection : public BaseSelectionevaluateObservable(kSphericityMin, mSphericity); this->evaluateObservable(kSphericityMax, mSphericity); + // for the trigger we need to pass an vector of 0 (false) and 1 (true) for all configured trigger selections + if (!triggerDecisions.empty()) { + std::vector trigger(triggerDecisions.size()); + std::transform(triggerDecisions.begin(), triggerDecisions.end(), trigger.begin(), [](bool b) { return b ? 1.0f : 0.0f; }); + this->evaluateObservable(kTriggers, trigger); + } + this->assembleBitmask(); }; @@ -330,18 +348,24 @@ class CollisionBuilder virtual ~CollisionBuilder() = default; template - void init(T1& confFilter, T2& confBits, T3& confRct, T4& confTrigger, T5& confTable, T6& initContext) + void init(T1& confFilter, T2& confBits, T3& confRct, T4& confCcdb, T5& confTable, T6& initContext) { mCollisionSelection.configure(confFilter, confBits); - if (confTrigger.useTrigger.value) { + if (!confBits.triggers.value.empty()) { mUseTrigger = true; - mTriggerNames = confTrigger.triggers.value; - mZorro.setBaseCCDBPath(confTrigger.ccdbPath.value); + for (size_t i = 0; i < confBits.triggers.value.size(); ++i) { + mTriggerNames += confBits.triggers.value[i]; + if (i != confBits.triggers.value.size() - 1) { + mTriggerNames += ","; + } + } + mZorro.setBaseCCDBPath(confCcdb.triggerPath.value); } if (confRct.useRctFlags.value) { mUseRctFlags = true; mRctFlagsChecker.init(confRct.label.value, confRct.useZdc.value, confRct.treatLimitedAcceptanceAsBad.value); } + mGrpPath = confCcdb.grpPath.value; LOG(info) << "Initialize femto collision builder..."; mProducedCollisions = utils::enableTable("FCols_001", confTable.produceCollisions.value, initContext); @@ -360,31 +384,45 @@ class CollisionBuilder LOG(info) << "Initialization done..."; } - template - void buildCollision(T1& bc, T2& col, T3& tracks, T4& ccdb, int magField) + template + void initCollision(T1& bc, T2& col, T3& tracks, T4& ccdb, T5& histRegistry) { - if (mUseTrigger) { - mZorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), mTriggerNames); + if (mRunNumber != bc.runNumber()) { + mRunNumber = bc.runNumber(); + static o2::parameters::GRPMagField* grpo = nullptr; + grpo = ccdb->template getForRun(mGrpPath, mRunNumber); + if (grpo == nullptr) { + LOG(fatal) << "GRP object not found for Run " << mRunNumber; + } + mMagField = static_cast(grpo->getNominalL3Field()); // get magnetic field in kG + + if (mUseTrigger) { + mZorro.initCCDB(ccdb.service, mRunNumber, bc.timestamp(), mTriggerNames); + mZorro.populateHistRegistry(histRegistry, mRunNumber); + } } - mCollisionSelection.setMagneticField(magField); + + mCollisionSelection.setMagneticField(mMagField); mCollisionSelection.setSphericity(tracks); mCollisionSelection.setMultiplicity(col); mCollisionSelection.setCentrality(col); - mCollisionSelection.applySelections(col); + + std::vector triggerDecisions = {}; + if (mUseTrigger) { + triggerDecisions = mZorro.getTriggerOfInterestResults(bc.globalBC()); + } + + mCollisionSelection.applySelections(col, triggerDecisions); } - template - bool checkCollision(T1 const& bc, T2 const& col) + template + bool checkCollision(T1 const& col) { - // First: if triggers are enabled, the object must be selected - if (mUseTrigger && !mZorro.isSelected(bc.globalBC())) { - return false; - } - // Then: if RCT flags are enabled, check them + // check RCT flags first if (mUseRctFlags && !mRctFlagsChecker(col)) { return false; } - // Finally: do the expensive checks + // make other checks return mCollisionSelection.checkFilters(col) && mCollisionSelection.passesAllRequiredSelections(); } @@ -438,6 +476,9 @@ class CollisionBuilder CollisionSelection mCollisionSelection; Zorro mZorro; bool mUseTrigger = false; + int mRunNumber = -1; + std::string mGrpPath = std::string(""); + int mMagField = 0; aod::rctsel::RCTFlagsChecker mRctFlagsChecker; bool mUseRctFlags = false; std::string mTriggerNames = std::string(""); diff --git a/PWGCF/Femto/Core/selectionContainer.h b/PWGCF/Femto/Core/selectionContainer.h index 420af6554d6..0ebeb803b0e 100644 --- a/PWGCF/Femto/Core/selectionContainer.h +++ b/PWGCF/Femto/Core/selectionContainer.h @@ -41,11 +41,12 @@ enum LimitType { kUpperLimit, ///< simple upper limit for the value, 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 - kEqual, ///< values need to be equal, e.g. sign = 1 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) + 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 }; std::unordered_map limitTypeAsStrings = { @@ -53,11 +54,14 @@ std::unordered_map limitTypeAsStrings = { {kAbsUpperLimit, "Absolute Upper Limit"}, {kLowerLimit, "Lower Limit"}, {kAbsLowerLimit, "Absolute Lower Limit"}, - {kEqual, "Equal"}, {kUpperFunctionLimit, "Upper Function Limit"}, {kAbsUpperFunctionLimit, "Absolute Upper Function Limit"}, {kLowerFunctionLimit, "Lower Function Limit"}, - {kAbsLowerFunctionLimit, "Absolute Lower Function Limit"}}; + {kAbsLowerFunctionLimit, "Absolute Lower Function Limit"}, + {kEqual, "Equal"}, + {kEqualArray, "EqualArray"}, + +}; }; // namespace limits @@ -70,7 +74,8 @@ class SelectionContainer { public: /// Default constructor - SelectionContainer() {} + SelectionContainer() = default; + ~SelectionContainer() = default; /// \brief Constructor for static value-based selection. /// \param SelectionValues Vector of values for the selection. @@ -113,9 +118,7 @@ class SelectionContainer LOG(fatal) << "Too many selections for single a observable. Limit is " << sizeof(BitmaskType) * CHAR_BIT; } for (std::size_t i = 0; i < functions.size(); i++) { - const std::string& func = functions.at(i); - const std::string& safeFunc = func.empty() ? "0.1" : func; // in case string is empty, set to constant value of 0.1 - mSelectionFunctions.emplace_back((baseName + std::to_string(i)).c_str(), safeFunc.c_str(), lowerLimit, upperLimit); + mSelectionFunctions.emplace_back((baseName + 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 @@ -128,8 +131,6 @@ class SelectionContainer } } - virtual ~SelectionContainer() = default; - /// \brief Sort static selection values based on the limit type. void sortSelections() { @@ -140,7 +141,6 @@ class SelectionContainer break; case (limits::kLowerLimit): case (limits::kAbsLowerLimit): - case (limits::kEqual): std::sort(mSelectionValues.begin(), mSelectionValues.end(), [](T a, T b) { return a < b; }); break; default: @@ -166,6 +166,17 @@ class SelectionContainer } } + /// \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 getComments() const { return mComments; } + /// \brief Update selection limits using internal functions evaluated at a given value. /// \param value Input value to evaluate functions at. void updateLimits(T value) @@ -237,6 +248,26 @@ class SelectionContainer } } + /// \brief Evaluate which selection criteria are fulfilled for a given value. + /// \param values Values of the observable to evaluate + void evaluate(std::vector& 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): + if (std::fabs(values.at(i) - mSelectionValues.at(i)) < constants::math::Epsilon) { + mBitmask.set(i); + } + break; + default: + continue; + } + } + } + /// \brief Retrieve the bitmask indicating which selections were passed. /// \return Bitset representing passed selections. std::bitset getBitmask() const @@ -260,8 +291,9 @@ class SelectionContainer bool passesAsMinimalCut() const { if (mIsMinimalCut) { - // check if loosest bit is set - return mBitmask.test(0); + // 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(); } else { // if selection is not marked as a minimal cut, we return true by default return true; @@ -325,6 +357,7 @@ class SelectionContainer private: std::vector mSelectionValues = {}; ///< Values used for the selection + std::vector mComments = {}; ///< Comments for the values std::vector mSelectionFunctions = {}; ///< Function used for the selection limits::LimitType mLimitType; ///< Limit type of selection std::bitset mBitmask = {}; ///< bitmask for the observable diff --git a/PWGCF/Femto/TableProducer/femtoProducer.cxx b/PWGCF/Femto/TableProducer/femtoProducer.cxx index 207309a5781..c21281d0616 100644 --- a/PWGCF/Femto/TableProducer/femtoProducer.cxx +++ b/PWGCF/Femto/TableProducer/femtoProducer.cxx @@ -32,7 +32,6 @@ #include "Common/DataModel/TrackSelectionTables.h" #include "CCDB/BasicCCDBManager.h" -#include "DataFormatsParameters/GRPMagField.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisHelpers.h" #include "Framework/AnalysisTask.h" @@ -59,8 +58,6 @@ namespace o2::analysis::femto namespace consumeddata { using Run3PpCollisions = soa::Join; -// FIXME: sometimes people want to run analyis even though centrality calibration is not available yet -// using Run3PpWithoutCentCollisions = soa::Join; using Run3FullPidTracks = soa::Join ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "URL to ccdb"}; - Configurable grpPath{"grpPath", "GLO/Config/GRPMagField", "Path to GRP object (Run3 -> GLO/Config/GRPMagField/Run2 -> GLO/GRP/GRP"}; - } ConfOptions; + // ccdb + collisionbuilder::ConfCcdb confCcdb; // collision builder collisionbuilder::CollisionBuilderProducts collisionBuilderProducts; @@ -92,7 +85,6 @@ struct FemtoProducer { collisionbuilder::ConfCollisionFilters confCollisionFilters; collisionbuilder::ConfCollisionBits confCollisionBits; collisionbuilder::ConfCollisionRctFlags confCollisionRctFlags; - collisionbuilder::ConfCollisionTriggers confCollisionTriggers; collisionbuilder::CollisionBuilder collisionBuilder; // track builder @@ -160,40 +152,23 @@ struct FemtoProducer { // histogramming // add histograms in next iteration - // HistogramRegistry hRegistry{"FemtoProducer", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry hRegistry{"FemtoProducer", {}, OutputObjHandlingPolicy::AnalysisObject}; // data members - int runNumber = -1; - int magField = 0.f; Service ccdb; /// Accessing the CCDB std::unordered_map indexMapTracks; // for mapping tracks to lambdas, cascades and resonances - void initFromCcdb(o2::aod::BCsWithTimestamps::iterator const& bc) - { - if (runNumber == bc.runNumber()) - return; - auto timestamp = bc.timestamp(); - static o2::parameters::GRPMagField* grpo = nullptr; - grpo = ccdb->getForTimeStamp(ConfOptions.grpPath.value, timestamp); - if (grpo == nullptr) { - LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); - return; - } - magField = static_cast(grpo->getNominalL3Field()); // get magnetic field in kG - runNumber = bc.runNumber(); - }; - void init(InitContext& context) { // init ccdb - ccdb->setURL(ConfOptions.ccdbUrl.value); + ccdb->setURL(confCcdb.ccdbUrl.value); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); int64_t now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); ccdb->setCreatedNotAfter(now); // collision selection - collisionBuilder.init(confCollisionFilters, confCollisionBits, confCollisionRctFlags, confCollisionTriggers, confCollisionTables, context); + collisionBuilder.init(confCollisionFilters, confCollisionBits, confCollisionRctFlags, confCcdb, confCollisionTables, context); // configure track builder trackBuilder.init(confTrackBits, confTrackFilters, confTrackTables, context); @@ -232,9 +207,8 @@ struct FemtoProducer { void processTracks(T1 const& col, T2 const& /* bcs*/, T3 const& tracks, T4 const& tracksWithItsPid) { auto bc = col.template bc_as(); - initFromCcdb(bc); - collisionBuilder.buildCollision(bc, col, tracks, ccdb, magField); - if (!collisionBuilder.checkCollision(bc, col)) { + collisionBuilder.initCollision(bc, col, tracks, ccdb, hRegistry); + if (!collisionBuilder.checkCollision(col)) { return; } collisionBuilder.fillCollision(collisionBuilderProducts, col);