Add Bs workflow

Author: fchinu

PWGHF/Core/HfHelper.h

@@ -663,6 +663,12 @@ class HfHelper
     return candidate.m(std::array{o2::constants::physics::MassDSBar, o2::constants::physics::MassPiPlus});
   }
 
+  template <typename T>
+  auto invMassBsToJPsiPhi(const T& candidate)
+  {
+    return candidate.m(std::array{o2::constants::physics::MassJPsi, o2::constants::physics::MassKPlus});
+  }
+
   template <typename T>
   auto cosThetaStarBs(const T& candidate)
   {
@@ -879,7 +885,7 @@ class HfHelper
 
   // Apply topological cuts as defined in SelectorCuts.h
   /// \param candBp B+ candidate
-  /// \param cuts B+ candidate selection per pT bin"
+  /// \param cuts B+ candidate selection per pT bin
   /// \param binsPt pT bin limits
   /// \return true if candidate passes all selections
   template <typename T1, typename T2, typename T3>
@@ -1064,6 +1070,114 @@ class HfHelper
     return true;
   }
 
+  // Apply topological cuts as defined in SelectorCuts.h
+  /// \param candBs Bs candidate
+  /// \param cuts Bs candidate selection per pT bin
+  /// \param binsPt pT bin limits
+  /// \return true if candidate passes all selections
+  template <typename T1, typename T2, typename T3, typename T4, typename T5>
+  bool selectionBsToJPsiPhiTopol(const T1& candBs, const T2& candKa0, const T3& candKa1, const T4& cuts, const T5& binsPt)
+  {
+    auto ptcandBs = candBs.pt();
+    auto mcandBs = invMassBsToJPsiPhi(candBs);
+    std::array<float, 3> pVecKa0 = {candKa0.px(), candKa0.py(), candKa0.pz()};
+    std::array<float, 3> pVecKa1 = {candKa1.px(), candKa1.py(), candKa1.pz()};
+    auto mcandPhi = RecoDecay::m(std::array{pVecKa0, pVecKa1}, std::array{o2::constants::physics::MassKPlus, o2::constants::physics::MassKPlus});
+    auto ptJPsi = RecoDecay::pt(candBs.pxProng0(), candBs.pyProng0());
+    auto candJPsi = candBs.jPsi();
+    float pseudoPropDecLen = candBs.decayLengthXY() * mcandBs / ptcandBs;
+
+    int pTBin = o2::analysis::findBin(binsPt, ptcandBs);
+    if (pTBin == -1) {
+      return false;
+    }
+
+    // Bs mass cut
+    if (std::abs(mcandBs - o2::constants::physics::MassBPlus) > cuts->get(pTBin, "m")) {
+      return false;
+    }
+
+    // kaon pt
+    if (candKa0.pt() < cuts->get(pTBin, "pT K") &&
+        candKa1.pt() < cuts->get(pTBin, "pT K")) {
+      return false;
+    }
+
+    // J/Psi pt
+    if (ptJPsi < cuts->get(pTBin, "pT J/Psi")) {
+      return false;
+    }
+
+    // phi mass
+    if (std::abs(mcandPhi - o2::constants::physics::MassPhi) < cuts->get(pTBin, "DeltaM phi")) {
+      return false;
+    }
+
+    // J/Psi mass
+    if (std::abs(candJPsi.m() - o2::constants::physics::MassJPsi) < cuts->get(pTBin, "DeltaM J/Psi")) {
+      return false;
+    }
+
+    // d0(J/Psi)xd0(phi)
+    if (candBs.impactParameterProduct() > cuts->get(pTBin, "B Imp. Par. Product")) {
+      return false;
+    }
+
+    // Bs Decay length
+    if (candBs.decayLength() < cuts->get(pTBin, "B decLen")) {
+      return false;
+    }
+
+    // Bs Decay length XY
+    if (candBs.decayLengthXY() < cuts->get(pTBin, "B decLenXY")) {
+      return false;
+    }
+
+    // Bs CPA cut
+    if (candBs.cpa() < cuts->get(pTBin, "CPA")) {
+      return false;
+    }
+
+    // Bs CPAXY cut
+    if (candBs.cpaXY() < cuts->get(pTBin, "CPAXY")) {
+      return false;
+    }
+
+    // d0 of phi
+    if (std::abs(candBs.impactParameter1()) < cuts->get(pTBin, "d0 phi")) {
+      return false;
+    }
+
+    // d0 of J/Psi
+    if (std::abs(candBs.impactParameter0()) < cuts->get(pTBin, "d0 J/Psi")) {
+      return false;
+    }
+
+    // B pseudoproper decay length
+    if (pseudoPropDecLen < cuts->get(pTBin, "B pseudoprop. decLen")) {
+      return false;
+    }
+
+    return true;
+  }
+
+  /// Apply PID selection
+  /// \param pidTrackKa PID status of trackKa (prong1 of B+ candidate)
+  /// \param acceptPIDNotApplicable switch to accept Status::NotApplicable
+  /// \return true if prong1 of B+ candidate passes all selections
+  template <typename T1 = int, typename T2 = bool>
+  bool selectionBsToJPsiPhiPid(const T1& pidTrackKa, const T2& acceptPIDNotApplicable)
+  {
+    if (!acceptPIDNotApplicable && pidTrackKa != TrackSelectorPID::Accepted) {
+      return false;
+    }
+    if (acceptPIDNotApplicable && pidTrackKa == TrackSelectorPID::Rejected) {
+      return false;
+    }
+
+    return true;
+  }
+
   /// Apply topological cuts as defined in SelectorCuts.h
   /// \param candLb Lb candidate
   /// \param cuts Lb candidate selection per pT bin"

PWGHF/Core/HfMlResponseBsToJPsiPhiReduced.h

@@ -0,0 +1,180 @@
+// Copyright 2019-2020 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.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file HfMlResponseBsToJPsiPhiReduced.h
+/// \brief Class to compute the ML response for Bs0 → J/Psi phi analysis selections in the reduced format
+/// \author Fabrizio Chinu <fabrizio.chinu@cern.ch>, Università degli Studi and INFN Torino
+
+#ifndef PWGHF_CORE_HFMLRESPONSEBSTOJPSIPHIREDUCED_H_
+#define PWGHF_CORE_HFMLRESPONSEBSTOJPSIPHIREDUCED_H_
+
+#include <map>
+#include <string>
+#include <vector>
+
+#include "PWGHF/Core/HfMlResponse.h"
+#include "PWGHF/D2H/Utils/utilsRedDataFormat.h"
+
+// Fill the map of available input features
+// the key is the feature's name (std::string)
+// the value is the corresponding value in EnumInputFeatures
+#define FILL_MAP_BS(FEATURE)                                                 \
+  {                                                                          \
+    #FEATURE, static_cast<uint8_t>(InputFeaturesBsToJPsiPhiReduced::FEATURE) \
+  }
+
+// Check if the index of mCachedIndices (index associated to a FEATURE)
+// matches the entry in EnumInputFeatures associated to this FEATURE
+// if so, the inputFeatures vector is filled with the FEATURE's value
+// by calling the corresponding GETTER from OBJECT
+#define CHECK_AND_FILL_VEC_BS_FULL(OBJECT, FEATURE, GETTER)              \
+  case static_cast<uint8_t>(InputFeaturesBsToJPsiPhiReduced::FEATURE): { \
+    inputFeatures.emplace_back(OBJECT.GETTER());                         \
+    break;                                                               \
+  }
+
+// Check if the index of mCachedIndices (index associated to a FEATURE)
+// matches the entry in EnumInputFeatures associated to this FEATURE
+// if so, the inputFeatures vector is filled with the FEATURE's value
+// by calling the GETTER function taking OBJECT in argument
+#define CHECK_AND_FILL_VEC_BS_FUNC(OBJECT, FEATURE, GETTER)              \
+  case static_cast<uint8_t>(InputFeaturesBsToJPsiPhiReduced::FEATURE): { \
+    inputFeatures.emplace_back(GETTER(OBJECT));                          \
+    break;                                                               \
+  }
+
+// Specific case of CHECK_AND_FILL_VEC_BS_FULL(OBJECT, FEATURE, GETTER)
+// where OBJECT is named candidate and FEATURE = GETTER
+#define CHECK_AND_FILL_VEC_BS(GETTER)                                   \
+  case static_cast<uint8_t>(InputFeaturesBsToJPsiPhiReduced::GETTER): { \
+    inputFeatures.emplace_back(candidate.GETTER());                     \
+    break;                                                              \
+  }
+
+namespace o2::analysis
+{
+
+enum class InputFeaturesBsToJPsiPhiReduced : uint8_t {
+  ptProng0 = 0,
+  ptProng1,
+  impactParameter0,
+  impactParameter1,
+  impactParameterProduct,
+  chi2PCA,
+  decayLength,
+  decayLengthXY,
+  decayLengthNormalised,
+  decayLengthXYNormalised,
+  cpa,
+  cpaXY,
+  maxNormalisedDeltaIP,
+  tpcNSigmaKa0,
+  tofNSigmaKa0,
+  tpcTofNSigmaKa0,
+  tpcNSigmaKa1,
+  tofNSigmaKa1,
+  tpcTofNSigmaKa1
+};
+
+template <typename TypeOutputScore = float>
+class HfMlResponseBsToJPsiPhiReduced : public HfMlResponse<TypeOutputScore>
+{
+ public:
+  /// Default constructor
+  HfMlResponseBsToJPsiPhiReduced() = default;
+  /// Default destructor
+  virtual ~HfMlResponseBsToJPsiPhiReduced() = default;
+
+  /// Method to get the input features vector needed for ML inference
+  /// \param candidate is the Bs candidate
+  /// \param prong1 is the candidate's prong1
+  /// \return inputFeatures vector
+  template <typename T1, typename T2, typename T3>
+  std::vector<float> getInputFeatures(T1 const& candidate,
+                                      T2 const& prong1,
+                                      T3 const& prong2)
+  {
+    std::vector<float> inputFeatures;
+
+    for (const auto& idx : MlResponse<TypeOutputScore>::mCachedIndices) {
+      switch (idx) {
+        CHECK_AND_FILL_VEC_BS(ptProng0);
+        CHECK_AND_FILL_VEC_BS(ptProng1);
+        CHECK_AND_FILL_VEC_BS(impactParameter0);
+        CHECK_AND_FILL_VEC_BS(impactParameter1);
+        CHECK_AND_FILL_VEC_BS(impactParameterProduct);
+        CHECK_AND_FILL_VEC_BS(chi2PCA);
+        CHECK_AND_FILL_VEC_BS(decayLength);
+        CHECK_AND_FILL_VEC_BS(decayLengthXY);
+        CHECK_AND_FILL_VEC_BS(decayLengthNormalised);
+        CHECK_AND_FILL_VEC_BS(decayLengthXYNormalised);
+        CHECK_AND_FILL_VEC_BS(cpa);
+        CHECK_AND_FILL_VEC_BS(cpaXY);
+        CHECK_AND_FILL_VEC_BS(maxNormalisedDeltaIP);
+        // TPC PID variable
+        CHECK_AND_FILL_VEC_BS_FULL(prong1, tpcNSigmaKa0, tpcNSigmaKa);
+        // TOF PID variable
+        CHECK_AND_FILL_VEC_BS_FULL(prong1, tofNSigmaKa0, tofNSigmaKa);
+        // Combined PID variables
+        CHECK_AND_FILL_VEC_BS_FUNC(prong1, tpcTofNSigmaKa0, o2::pid_tpc_tof_utils::getTpcTofNSigmaKa1);
+        // TPC PID variable
+        CHECK_AND_FILL_VEC_BS_FULL(prong2, tpcNSigmaKa1, tpcNSigmaKa);
+        // TOF PID variable
+        CHECK_AND_FILL_VEC_BS_FULL(prong2, tofNSigmaKa1, tofNSigmaKa);
+        // Combined PID variables
+        CHECK_AND_FILL_VEC_BS_FUNC(prong2, tpcTofNSigmaKa1, o2::pid_tpc_tof_utils::getTpcTofNSigmaKa1);
+      }
+    }
+
+    return inputFeatures;
+  }
+
+ protected:
+  /// Method to fill the map of available input features
+  void setAvailableInputFeatures()
+  {
+    MlResponse<TypeOutputScore>::mAvailableInputFeatures = {
+      FILL_MAP_BS(ptProng0),
+      FILL_MAP_BS(ptProng1),
+      FILL_MAP_BS(impactParameter0),
+      FILL_MAP_BS(impactParameter1),
+      FILL_MAP_BS(impactParameterProduct),
+      FILL_MAP_BS(chi2PCA),
+      FILL_MAP_BS(decayLength),
+      FILL_MAP_BS(decayLengthXY),
+      FILL_MAP_BS(decayLengthNormalised),
+      FILL_MAP_BS(decayLengthXYNormalised),
+      FILL_MAP_BS(cpa),
+      FILL_MAP_BS(cpaXY),
+      FILL_MAP_BS(maxNormalisedDeltaIP),
+      // TPC PID variable
+      FILL_MAP_BS(tpcNSigmaKa0),
+      // TOF PID variable
+      FILL_MAP_BS(tofNSigmaKa0),
+      // Combined PID variable
+      FILL_MAP_BS(tpcTofNSigmaKa0),
+      // TPC PID variable
+      FILL_MAP_BS(tpcNSigmaKa1),
+      // TOF PID variable
+      FILL_MAP_BS(tofNSigmaKa1),
+      // Combined PID variable
+      FILL_MAP_BS(tpcTofNSigmaKa1)};
+  }
+};
+
+} // namespace o2::analysis
+
+#undef FILL_MAP_BS
+#undef CHECK_AND_FILL_VEC_BS_FULL
+#undef CHECK_AND_FILL_VEC_BS_FUNC
+#undef CHECK_AND_FILL_VEC_BS
+
+#endif // PWGHF_CORE_HFMLRESPONSEBSTOJPSIPHIREDUCED_H_

PWGHF/Core/SelectorCuts.h

@@ -1050,6 +1050,62 @@ static const std::vector<std::string> labelsPt = {
 static const std::vector<std::string> labelsCutVar = {"m", "CPA", "Chi2PCA", "d0 Ds", "d0 Pi", "pT Ds", "pT Pi", "Bs decLen", "Bs decLenXY", "Imp. Par. Product"};
 } // namespace hf_cuts_bs_to_ds_pi
 
+namespace hf_cuts_bs_to_jpsi_phi
+{
+static constexpr int NBinsPt = 12;
+static constexpr int NCutVars = 13;
+// default values for the pT bin edges (can be used to configure histogram axis)
+// offset by 1 from the bin numbers in cuts array
+constexpr double BinsPt[NBinsPt + 1] = {
+  0,
+  0.5,
+  1.0,
+  2.0,
+  3.0,
+  4.0,
+  5.0,
+  7.0,
+  10.0,
+  13.0,
+  16.0,
+  20.0,
+  24.0};
+
+auto vecBinsPt = std::vector<double>{BinsPt, BinsPt + NBinsPt + 1};
+
+// default values for the cuts
+// DeltaM CPA d0JPsi d0K pTJPsi pTK BDecayLength BDecayLengthXY BIPProd DeltaMJPsi JPsiIPProd
+constexpr double Cuts[NBinsPt][NCutVars] = {{1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 0 < pt < 0.5 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 0.5 < pt < 1 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 1 < pt < 2 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 2 < pt < 3 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 3 < pt < 4 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 4 < pt < 5 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 5 < pt < 7 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 7 < pt < 10 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 10 < pt < 13 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 13 < pt < 16 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.},  /* 16 < pt < 20 */
+                                            {1., 0.8, 0.8, 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.02, 0.}}; /* 20 < pt < 24 */
+// row labels
+static const std::vector<std::string> labelsPt = {
+  "pT bin 0",
+  "pT bin 1",
+  "pT bin 2",
+  "pT bin 3",
+  "pT bin 4",
+  "pT bin 5",
+  "pT bin 6",
+  "pT bin 7",
+  "pT bin 8",
+  "pT bin 9",
+  "pT bin 10",
+  "pT bin 11"};
+
+// column labels
+static const std::vector<std::string> labelsCutVar = {"m", "CPA", "CPAXY", "d0 J/Psi", "d0 phi", "pT J/Psi", "pT K", "B decLen", "B decLenXY", "B Imp. Par. Product", "DeltaM J/Psi", "DeltaM phi", "B pseudoprop. decLen"};
+} // namespace hf_cuts_bs_to_jpsi_phi
+
 namespace hf_cuts_bplus_to_d0_pi
 {
 static constexpr int NBinsPt = 12;

PWGHF/D2H/DataModel/ReducedDataModel.h

@@ -722,6 +722,7 @@ DECLARE_SOA_TABLE(HfRedBsDsMls, "AOD", "HFREDBSDSML", //! Table with ML scores f
                   o2::soa::Marker<1>);
 
 using HfRedCandBs = soa::Join<HfCandBsExt, HfRedBsProngs>;
+using HfRedCandBsToJPsiPhi = soa::Join<HfCandBsExt, HfRedBs2JPsiDaus>;
 
 namespace hf_cand_lb_reduced
 {