Feat: finish QA tasks and tools

Author: ariedel-cern

PWGCF/FemtoUnited/Core/BaseSelection.h

@@ -55,7 +55,29 @@ class BaseSelection
     mSelections.at(Observable) = SelectionContainer<T>(configSelections, limitType, SkipLastBit);
   }
 
-  void resetMinimalSelection() { mMinimalSelected = true; }
+  /// Pass the Configurable of selection values in the analysis task to the selection class
+  /// \param configSelection Vector from configurable containing the values employed for the selection
+  /// \param observableType Observable to be employed for the selection
+  /// \param limitType Type of the selection limit
+  void addSelection(std::string baseName, T lowerLimit, T upperLimit, std::vector<std::string>& configSelections, int Observable, limits::LimitType limitType, bool SkipLastBit)
+  {
+    if (static_cast<size_t>(Observable) >= NObservables) {
+      LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NObservables;
+    }
+    mNSelections += configSelections.size();
+    if (mNSelections >= 8 * sizeof(BitmaskType)) {
+      LOG(fatal) << "Too many selections. At most " << 8 * sizeof(BitmaskType) << " are supported";
+    }
+    mSelections.at(Observable) = SelectionContainer<T>(baseName, lowerLimit, upperLimit, configSelections, limitType, SkipLastBit);
+  }
+
+  void updateLimits(int observable, T value) { mSelections.at(observable).updateLimits(value); }
+
+  void reset()
+  {
+    mBitmask.reset();
+    mMinimalSelected = true;
+  }
 
   void setCheckMinimalSelection(bool checkMinimalSelection) { mCheckMinimalSelection = checkMinimalSelection; }
 
@@ -64,7 +86,11 @@ class BaseSelection
   /// \param value Value of the observable
   void setBitmaskForObservable(int observable, T value)
   {
-    // if any object did not pass minimal selections, there is no point in setting bitmask for other observables
+    // if there are no values set, bail out
+    if (mSelections.at(observable).empty()) {
+      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
     // can be deactivate by setting mCheckMinimalSelection to false
     if (mCheckMinimalSelection == true && mMinimalSelected == false) {
@@ -80,7 +106,11 @@ class BaseSelection
   }
 
   /// check if minimal Selections are passed
-  bool getMinimalSelection() { return mMinimalSelected; }
+  bool
+    getMinimalSelection()
+  {
+    return mMinimalSelected;
+  }
 
   /// check if any Selections are passed
   bool getAnySelection()
@@ -96,9 +126,11 @@ class BaseSelection
     // if minimal selections are not passed, just set bitmask to 0
     if (mCheckMinimalSelection == true && mMinimalSelected == false) {
       mBitmask.reset();
+      return;
     }
 
     // to assemble bitmask, convert all bitmask into integers
+    // shift the current one and add the new bits
     uint64_t result = 0u;
     int shift = 0;
     uint64_t value = 0u;

PWGCF/FemtoUnited/Core/DataTypes.h

@@ -23,13 +23,10 @@ namespace femtodatatypes
 // Bitmaks for tracks
 using TrackMaskType = uint32_t;
 using TrackPidMaskType = uint32_t;
-// using TrackTPCMaskType = uint16_t;
-// using TrackTOFMaskType = uint16_t;
-// using TrackTPCTOFMaskType = uint16_t;
 
 // Bitmaks for vzeros and daughters
 using VzeroMaskType = uint32_t;
-// using VzeroDauTrackMaskType = uint8_t;
+// using VzeroDaughterMaskType = uint8_t;
 // using VzeroDauTPCMaskType = uint8_t;
 
 // Bitmaks for cascades, vzero daughter and bachelor

PWGCF/FemtoUnited/Core/SelectionContainer.h

@@ -21,6 +21,8 @@
 #include <vector>
 #include <algorithm>
 
+#include <TF1.h>
+
 #include "fairlogger/Logger.h"
 #include "CommonConstants/MathConstants.h"
 
@@ -30,13 +32,15 @@ namespace o2::analysis::femtounited
 /// 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)
-                 kEqual ///< values need to be equal, e.g. sign = 1
+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
+                 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 upper limit of a function value, e.g. DCA_xy > f(pt)
+                 kAbsLowerFunctionLimit  ///< upper limit of an absolute value given by a function, e.g. |DCA_xy| > f(pt)
 };
 }
 // bitsets need number of bits at compile time. Set reasonable limit here
@@ -70,6 +74,26 @@ class SelectionContainer
     sortSelections();
   }
 
+  /// Constructor
+  /// \param values Values for the selection
+  /// \param limitType Type of limit of the selection
+  SelectionContainer(std::string baseName, T lowerLimit, T upperLimit, std::vector<std::string>& functions, limits::LimitType limitType, bool SkipLastBit)
+    : mLimitType(limitType), mSkipLastBit(SkipLastBit)
+  {
+    if (mValues.size() > BitmaskMaxSize) {
+      LOG(fatal) << "Too many selections for single a observable. Current limit is " << BitmaskMaxSize;
+    }
+    for (std::size_t i = 0; i < functions.size(); i++) {
+      mFunctions.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
+    // here we rely on the user that the functions can be ordered like this over the whole interval
+    sortFunctions((lowerLimit + upperLimit) / 2);
+    // initialize the values also to the midpoint
+    this->updateLimits((lowerLimit + upperLimit) / 2);
+  }
+
   /// Destructor
   virtual ~SelectionContainer() = default;
 
@@ -86,6 +110,34 @@ class SelectionContainer
       case (limits::kEqual):
         std::sort(mValues.begin(), mValues.end(), [](T a, T b) { return a <= b; });
         break;
+      default:
+        break;
+    }
+  }
+
+  // sort limit functions
+  void sortFunctions(T value)
+  {
+    switch (mLimitType) {
+      case (limits::kUpperFunctionLimit):
+      case (limits::kAbsUpperFunctionLimit):
+        std::sort(mFunctions.begin(), mFunctions.end(), [value](TF1 a, TF1 b) { return a.Eval(value) >= b.Eval(value); });
+        break;
+      case (limits::kLowerFunctionLimit):
+      case (limits::kAbsLowerFunctionLimit):
+        std::sort(mFunctions.begin(), mFunctions.end(), [value](TF1 a, TF1 b) { return a.Eval(value) <= b.Eval(value); });
+        break;
+      default:
+        break;
+    }
+  }
+
+  // update the selection limits depending on the passed function
+  void updateLimits(T value)
+  {
+    // functions are ordered so just add the values in the same order
+    for (std::size_t i = 0; i < mValues.size(); i++) {
+      mValues.at(i) = mFunctions.at(i).Eval(value);
     }
   }
 
@@ -95,7 +147,7 @@ class SelectionContainer
   {
     // better safe than sorry and reset the bitmask before you evaluate a new observable
     mBitmask.reset();
-    // the values are order, as soon as one comparison is not true, we can break out of the loop
+    // the values are ordered, for most loost 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
     for (size_t i = 0; i < mValues.size(); i++) {
@@ -177,6 +229,9 @@ class SelectionContainer
   /// Get shift for final bitmask
   int getShift()
   {
+    if (mValues.empty()) {
+      return 0;
+    }
     if (mSkipLastBit) {
       return static_cast<int>(mValues.size() - 1);
     } else {
@@ -186,6 +241,7 @@ class SelectionContainer
 
  private:
   std::vector<T> mValues{};             ///< Values used for the selection
+  std::vector<TF1> mFunctions{};        ///< Values used for the selection
   limits::LimitType mLimitType;         ///< Limit type of selection
   std::bitset<BitmaskMaxSize> mBitmask; ///< bitmask for a given observable
   bool mSkipLastBit = false;

PWGCF/FemtoUnited/Core/TrackHistManager.h

@@ -1,3 +1,4 @@
+//
 // Copyright 2019-2022 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
@@ -23,6 +24,7 @@
 #include "Framework/HistogramRegistry.h"
 #include "PWGCF/FemtoUnited/Core/HistManager.h"
 #include "PWGCF/FemtoUnited/Core/Modes.h"
+#include "PWGCF/FemtoUnited/Utils/FemtoUtils.h"
 
 using namespace ::o2::framework;
 
@@ -42,6 +44,9 @@ enum TrackHist {
   kItsClusterIb,
   kTpcCluster,
   kTpcClusterShared,
+  // kDcaxy,
+  // kDcaz,
+  // kDca,
   // 2d qa
   kPtVsEta,
   kPtVsPhi,
@@ -50,7 +55,11 @@ enum TrackHist {
   kPtVsTpcCluster,
   kPtVsTpcClusterShared,
   kTpcClusterVsTpcClusterShared,
+  kPtVsDcaxy,
+  kPtVsDcaz,
+  kPtVsDca,
   // its pid
+  kItsSignal,
   kItsElectron,
   kItsPion,
   kItsKaon,
@@ -109,6 +118,10 @@ constexpr std::array<Histmanager::HistInfo<TrackHist>, kTrackHistogramLast> Hist
    {kPtVsTpcCluster, kTH2F, "hPtVsTpcCluster", "p_{T} vs TPC cluster found; p_{T} (GeV/#it{c}) ; TPC cluster found"},
    {kPtVsTpcClusterShared, kTH2F, "hPtVsTpcClusterShared", "p_{T} vs TPC cluster shared; p_{T} (GeV/#it{c}) ; TPC cluster shared"},
    {kTpcClusterVsTpcClusterShared, kTH2F, "hTpcClusterVsTpcClusterShared", "TPC cluster found vs TPC cluster shared; TPC cluster found; TPC cluster shared"},
+   {kPtVsDcaxy, kTH2F, "hPtVsDcaxy", "p_{T} vs DCA_{XY}; p_{T} (GeV/#it{c}); DCA_{XY} (cm)"},
+   {kPtVsDcaz, kTH2F, "hPtVsDcaz", "p_{T} vs DCA_{Z}; p_{T} (GeV/#it{c}); DCA_{Z} (cm)"},
+   {kPtVsDca, kTH2F, "hPtVsDca", "p_{T} vs DCA; p_{T} (GeV/#it{c}); DCA (cm)"},
+   {kItsSignal, kTH2F, "hItsSignal", "ITS Signal; p (GeV/#it{c}) ; <ITS Cluster Size> x <cos #lambda>"},
    {kItsElectron, kTH2F, "hItsPidElectron", "TPC PID Electron; p (GeV/#it{c}) ; n#sigma_{TPC,el}"},
    {kItsPion, kTH2F, "hItsPidPion", "ITS PID Pion; p (GeV/#it{c}) ; n#sigma_{ITS,pi}"},
    {kItsKaon, kTH2F, "hItsPidKaon", "ITS PID Kaon; p (GeV/#it{c}) ; n#sigma_{ITS,ka}"},
@@ -182,8 +195,14 @@ class TrackHistManager
       mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsTpcClusterShared, HistTable), GetHistDesc(kPtVsTpcClusterShared, HistTable), GetHistType(kPtVsTpcClusterShared, HistTable), {Specs[kPtVsTpcClusterShared]});
       mHistogramRegistry->add(qaDir + GetHistNamev2(kTpcClusterVsTpcClusterShared, HistTable), GetHistDesc(kTpcClusterVsTpcClusterShared, HistTable), GetHistType(kTpcClusterVsTpcClusterShared, HistTable), {Specs[kTpcClusterVsTpcClusterShared]});
 
+      // dca
+      mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsDcaxy, HistTable), GetHistDesc(kPtVsDcaxy, HistTable), GetHistType(kPtVsDcaxy, HistTable), {Specs[kPtVsDcaxy]});
+      mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsDcaz, HistTable), GetHistDesc(kPtVsDcaz, HistTable), GetHistType(kPtVsDcaz, HistTable), {Specs[kPtVsDcaz]});
+      mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsDca, HistTable), GetHistDesc(kPtVsDca, HistTable), GetHistType(kPtVsDca, HistTable), {Specs[kPtVsDca]});
+
       std::string pidDir = std::string(PidDir);
 
+      mHistogramRegistry->add(pidDir + GetHistNamev2(kItsSignal, HistTable), GetHistDesc(kItsSignal, HistTable), GetHistType(kItsSignal, HistTable), {Specs[kItsSignal]});
       mHistogramRegistry->add(pidDir + GetHistNamev2(kItsElectron, HistTable), GetHistDesc(kItsElectron, HistTable), GetHistType(kItsElectron, HistTable), {Specs[kItsElectron]});
       mHistogramRegistry->add(pidDir + GetHistNamev2(kItsPion, HistTable), GetHistDesc(kItsPion, HistTable), GetHistType(kItsPion, HistTable), {Specs[kItsPion]});
       mHistogramRegistry->add(pidDir + GetHistNamev2(kItsKaon, HistTable), GetHistDesc(kItsKaon, HistTable), GetHistType(kItsKaon, HistTable), {Specs[kItsKaon]});
@@ -244,6 +263,11 @@ class TrackHistManager
       mHistogramRegistry->fill(HIST(QaDir) + HIST(GetHistName(kPtVsTpcClusterShared, HistTable)), track.pt(), static_cast<float>(track.tpcNClsShared()));
       mHistogramRegistry->fill(HIST(QaDir) + HIST(GetHistName(kTpcClusterVsTpcClusterShared, HistTable)), static_cast<float>(track.tpcNClsFound()), static_cast<float>(track.tpcNClsShared()));
 
+      mHistogramRegistry->fill(HIST(QaDir) + HIST(GetHistName(kPtVsDcaxy, HistTable)), track.pt(), track.dcaXY());
+      mHistogramRegistry->fill(HIST(QaDir) + HIST(GetHistName(kPtVsDcaz, HistTable)), track.pt(), track.dcaZ());
+      mHistogramRegistry->fill(HIST(QaDir) + HIST(GetHistName(kPtVsDca, HistTable)), track.pt(), utils::geometricMean(track.dcaXY(), track.dcaZ()));
+
+      mHistogramRegistry->fill(HIST(PidDir) + HIST(GetHistName(kItsSignal, HistTable)), track.p(), o2::analysis::femtounited::utils::itsSignal(track));
       mHistogramRegistry->fill(HIST(PidDir) + HIST(GetHistName(kItsElectron, HistTable)), track.p(), track.itsNSigmaEl());
       mHistogramRegistry->fill(HIST(PidDir) + HIST(GetHistName(kItsPion, HistTable)), track.p(), track.itsNSigmaPi());
       mHistogramRegistry->fill(HIST(PidDir) + HIST(GetHistName(kItsKaon, HistTable)), track.p(), track.itsNSigmaKa());

PWGCF/FemtoUnited/Core/TrackPidSelection.h

@@ -17,6 +17,7 @@
 #define PWGCF_FEMTOUNITED_CORE_TRACKPIDSELECTION_H_
 
 #include <cmath>
+
 #include "PWGCF/FemtoUnited/Core/BaseSelection.h"
 #include "PWGCF/FemtoUnited/Core/DataTypes.h"
 #include "PWGCF/FemtoUnited/Utils/FemtoUtils.h"
@@ -73,6 +74,8 @@ class TrackPidSelection : public BaseSelection<float, o2::aod::femtodatatypes::T
   template <class track>
   void ApplySelections(track const& Track)
   {
+    this->reset();
+
     // its pid
     this->setBitmaskForObservable(TrackPidSels::kItsElectron, Track.itsNSigmaEl());
     this->setBitmaskForObservable(TrackPidSels::kItsPion, Track.itsNSigmaPi());

PWGCF/FemtoUnited/Core/TrackSelection.h

@@ -33,8 +33,8 @@ enum TrackSels {
   kTPCsClsMax,   ///< Max. number of shared TPC clusters
   kITSnClsMin,   ///< Min. number of ITS clusters
   kITSnClsIbMin, ///< Min. number of ITS clusters in the inner barrel
-  kDCAxyMax,     ///< Max. DCA_xy (cm)
-  kDCAzMax,      ///< Max. DCA_z (cm)
+  kDCAxyMax,     ///< Max. DCA_xy (cm) as a function of pT
+  kDCAzMax,      ///< Max. DCA_z (cm) as a function of pT
   kTrackselMax
 };
 
@@ -48,13 +48,21 @@ class TrackSelection : public BaseSelection<float, o2::aod::femtodatatypes::Trac
   template <class track>
   void ApplySelections(track const& Track)
   {
-    this->resetMinimalSelection();
+    this->reset();
     this->setBitmaskForObservable(TrackSels::kSign, Track.sign());
     this->setBitmaskForObservable(TrackSels::kTPCnClsMin, Track.tpcNClsFound());
     this->setBitmaskForObservable(TrackSels::kTPCcRowsMin, Track.tpcNClsCrossedRows());
     this->setBitmaskForObservable(TrackSels::kTPCsClsMax, Track.tpcNClsShared());
     this->setBitmaskForObservable(TrackSels::kITSnClsMin, Track.itsNCls());
     this->setBitmaskForObservable(TrackSels::kITSnClsIbMin, Track.itsNClsInnerBarrel());
+
+    // evalue bitmask for pt dependent dca cuts
+    this->updateLimits(TrackSels::kDCAxyMax, Track.pt());
+    this->setBitmaskForObservable(TrackSels::kDCAxyMax, Track.dcaXY());
+
+    this->updateLimits(TrackSels::kDCAzMax, Track.pt());
+    this->setBitmaskForObservable(TrackSels::kDCAzMax, Track.dcaZ());
+
     this->assembleBismask();
   };
 }; // namespace femtoDream