Merge branch 'AliceO2Group:master' into master

CODEOWNERS

@@ -55,7 +55,7 @@
 /PWGMM/UE   @alibuild @aalkin @aortizve
 
 /PWGUD @alibuild @pbuehler @abylinkin @rolavick
-/PWGJE @alibuild @lhavener @maoyx @nzardosh @fjonasALICE @mfasDa
+/PWGJE @alibuild @lhavener @maoyx @nzardosh @fjonasALICE @mfasDa @mhemmer-cern
 /Tools/PIDML @alibuild @saganatt
 /Tools/ML @alibuild @fcatalan92 @fmazzasc
 /Tutorials/PWGCF @alibuild @jgrosseo @saganatt @victor-gonzalez @zchochul

Common/DataModel/PIDResponseITS.h

@@ -0,0 +1,155 @@
+// 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   PIDResponseITS.h
+/// \since  2024-11-12
+/// \author Nicolò Jacazio [email protected]
+/// \author Francesco Mazzaschi [email protected]
+/// \brief  Set of tables, tasks and utilities to provide the interface between
+///         the analysis data model and the PID response of the ITS
+///
+
+#ifndef COMMON_DATAMODEL_PIDRESPONSEITS_H_
+#define COMMON_DATAMODEL_PIDRESPONSEITS_H_
+
+// O2 includes
+#include "Framework/ASoA.h"
+#include "Framework/AnalysisDataModel.h"
+#include "ReconstructionDataFormats/PID.h"
+#include "Framework/Logger.h"
+
+namespace o2::aod
+{
+
+struct ITSResponse {
+  static float averageClusterSize(uint32_t itsClusterSizes)
+  {
+    float average = 0;
+    int nclusters = 0;
+
+    for (int layer = 0; layer < 7; layer++) {
+      if ((itsClusterSizes >> (layer * 4)) & 0xf) {
+        nclusters++;
+        average += (itsClusterSizes >> (layer * 4)) & 0xf;
+      }
+    }
+    if (nclusters == 0) {
+      return 0;
+    }
+    return average / nclusters;
+  };
+
+  template <o2::track::PID::ID id>
+  static float expSignal(const float momentum)
+  {
+    static constexpr float inverseMass = 1. / o2::track::pid_constants::sMasses[id];
+    static constexpr float charge = static_cast<float>(o2::track::pid_constants::sCharges[id]);
+    const float bg = momentum * inverseMass;
+    return (mITSRespParams[0] / (std::pow(bg, mITSRespParams[1])) + mITSRespParams[2]) * std::pow(charge, mChargeFactor);
+  }
+
+  template <o2::track::PID::ID id>
+  static float nSigmaITS(uint32_t itsClusterSizes, float momentum)
+  {
+    const float exp = expSignal<id>(momentum);
+    const float average = averageClusterSize(itsClusterSizes);
+    const float resolution = mResolution * exp;
+    return (average - exp) / resolution;
+  };
+
+  static void setParameters(float p0, float p1, float p2, float chargeFactor, float resolution)
+  {
+    if (mIsInitialized) {
+      LOG(fatal) << "ITSResponse parameters already initialized";
+    }
+    mIsInitialized = true;
+    mITSRespParams[0] = p0;
+    mITSRespParams[1] = p1;
+    mITSRespParams[2] = p2;
+    mChargeFactor = chargeFactor;
+    mResolution = resolution;
+  }
+
+ private:
+  static std::array<float, 3> mITSRespParams;
+  static float mChargeFactor;
+  static float mResolution;
+  static bool mIsInitialized;
+};
+
+std::array<float, 3> ITSResponse::mITSRespParams = {0.903, 2.014, 2.440};
+float ITSResponse::mChargeFactor = 2.299999952316284f;
+float ITSResponse::mResolution = 0.15f;
+bool ITSResponse::mIsInitialized = false;
+
+namespace pidits
+{
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaElImp, itsNSigmaEl, //! Nsigma separation with the ITS detector for electrons
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Electron>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaMuImp, itsNSigmaMu, //! Nsigma separation with the ITS detector for muons
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Muon>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaPiImp, itsNSigmaPi, //! Nsigma separation with the ITS detector for pions
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Pion>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaKaImp, itsNSigmaKa, //! Nsigma separation with the ITS detector for kaons
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Kaon>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaPrImp, itsNSigmaPr, //! Nsigma separation with the ITS detector for protons
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Proton>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaDeImp, itsNSigmaDe, //! Nsigma separation with the ITS detector for deuterons
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Deuteron>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaTrImp, itsNSigmaTr, //! Nsigma separation with the ITS detector for tritons
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Triton>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaHeImp, itsNSigmaHe, //! Nsigma separation with the ITS detector for helium3
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Helium3>(itsClusterSizes, momentum);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(ITSNSigmaAlImp, itsNSigmaAl, //! Nsigma separation with the ITS detector for alphas
+                           [](uint32_t itsClusterSizes, float momentum) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Alpha>(itsClusterSizes, momentum);
+                           });
+
+// Define user friendly names for the columns to join with the tracks
+using ITSNSigmaEl = ITSNSigmaElImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaMu = ITSNSigmaMuImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaPi = ITSNSigmaPiImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaKa = ITSNSigmaKaImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaPr = ITSNSigmaPrImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaDe = ITSNSigmaDeImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaTr = ITSNSigmaTrImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaHe = ITSNSigmaHeImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+using ITSNSigmaAl = ITSNSigmaAlImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P>;
+
+} // namespace pidits
+} // namespace o2::aod
+
+#endif // COMMON_DATAMODEL_PIDRESPONSEITS_H_

Common/TableProducer/PID/CMakeLists.txt

@@ -9,8 +9,13 @@
 # granted to it by virtue of its status as an Intergovernmental Organization
 # or submit itself to any jurisdiction.
 
-# TOF
+# ITS
+o2physics_add_dpl_workflow(pid-its
+    SOURCES pidITS.cxx
+    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+    COMPONENT_NAME Analysis)
 
+# TOF
 o2physics_add_dpl_workflow(pid-tof-base
                     SOURCES pidTOFBase.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::TOFBase

Common/TableProducer/PID/pidITS.cxx

@@ -0,0 +1,110 @@
+// 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   pidITS.cxx
+/// \since  2024-11-12
+/// \author Nicolò Jacazio [email protected]
+/// \author Francesco Mazzaschi [email protected]
+/// \brief  Task to produce PID tables for ITS split for each particle.
+///         Only the tables for the mass hypotheses requested are filled, the others are sent empty.
+///
+
+#include <utility>
+#include <vector>
+#include <string>
+
+// O2 includes
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "ReconstructionDataFormats/Track.h"
+#include "CCDB/BasicCCDBManager.h"
+#include "TOFBase/EventTimeMaker.h"
+
+// O2Physics includes
+#include "Common/DataModel/PIDResponseITS.h"
+#include "MetadataHelper.h"
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::track;
+
+MetadataHelper metadataInfo;
+
+static constexpr int nCases = 2;
+static constexpr int nParameters = 5;
+static const std::vector<std::string> casesNames{"Data", "MC"};
+static const std::vector<std::string> parameterNames{"bb1", "bb2", "bb3", "Charge exponent", "Resolution"};
+static constexpr float defaultParameters[nCases][nParameters]{{0.903, 2.014, 2.440, 2.299999952316284f, 0.15f},
+                                                              {0.903, 2.014, 2.440, 2.299999952316284f, 0.15f}};
+
+/// Task to produce the ITS PID information for each particle species
+/// The parametrization is: [p0/(bg)**p1 + p2] * pow(q, p3), being bg = p/m and q the charge
+struct itsPid {
+
+  Configurable<LabeledArray<float>> itsParams{"itsParams",
+                                              {defaultParameters[0], nCases, nParameters, casesNames, parameterNames},
+                                              "Response parameters"};
+  Configurable<bool> getFromCCDB{"getFromCCDB", false, "Get the parameters from CCDB"};
+
+  Service<o2::ccdb::BasicCCDBManager> ccdb;
+  Configurable<std::string> paramfile{"param-file", "", "Path to the parametrization object, if empty the parametrization is not taken from file"};
+  Configurable<std::string> url{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
+  Configurable<std::string> ccdbPath{"ccdbPath", "Analysis/PID/TPC/Response", "Path of the TPC parametrization on the CCDB"};
+  Configurable<std::string> recoPass{"recoPass", "", "Reconstruction pass name for CCDB query (automatically takes latest object for timestamp if blank)"};
+  Configurable<int64_t> ccdbTimestamp{"ccdb-timestamp", 0, "timestamp of the object used to query in CCDB the detector response. Exceptions: -1 gets the latest object, 0 gets the run dependent timestamp"};
+
+  void init(o2::framework::InitContext&)
+  {
+    if (getFromCCDB) {
+      ccdb->setURL(url.value);
+      ccdb->setCaching(true);
+      ccdb->setLocalObjectValidityChecking();
+      ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
+      LOG(fatal) << "Not implemented yet";
+    } else {
+      const char* key = metadataInfo.isMC() ? "MC" : "Data";
+      o2::aod::ITSResponse::setParameters(itsParams->get(key, "bb1"),
+                                          itsParams->get(key, "bb2"),
+                                          itsParams->get(key, "bb3"),
+                                          itsParams->get(key, "Charge exponent"),
+                                          itsParams->get(key, "Resolution"));
+    }
+  }
+
+  /// Dummy process function for BCs, needed in case both Run2 and Run3 process functions are disabled
+  void process(aod::Timestamps const&) {}
+
+  void processTest(o2::soa::Join<aod::TracksIU, aod::TracksExtra> const& tracks)
+  {
+    auto tracksWithPid = soa::Attach<o2::soa::Join<aod::TracksIU, aod::TracksExtra>,
+                                     aod::pidits::ITSNSigmaEl, aod::pidits::ITSNSigmaMu, aod::pidits::ITSNSigmaPi,
+                                     aod::pidits::ITSNSigmaKa, aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe,
+                                     aod::pidits::ITSNSigmaTr, aod::pidits::ITSNSigmaHe, aod::pidits::ITSNSigmaAl>(tracks);
+
+    for (const auto& track : tracksWithPid) {
+      LOG(info) << track.itsNSigmaEl();
+      LOG(info) << track.itsNSigmaPi();
+      LOG(info) << track.itsNSigmaPr();
+    }
+  }
+  PROCESS_SWITCH(itsPid, processTest, "Produce a test", false);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  // Parse the metadata
+  metadataInfo.initMetadata(cfgc);
+  auto workflow = WorkflowSpec{adaptAnalysisTask<itsPid>(cfgc)};
+  return workflow;
+}

EventFiltering/PWGHF/HFFilterHelpers.h

@@ -410,6 +410,8 @@ class HfFilterHelper
   int8_t isBDTSelected(const T& scores, const U& thresholdBDTScores);
   template <bool isKaonTrack, typename T>
   bool isSelectedKaonFromXicResoToSigmaC(const T& track);
+  template <typename T1>
+  inline bool isCharmHadronMassInSbRegions(T1 const& massHypo1, T1 const& massHypo2, const float& lowLimitSB, const float& upLimitSB);
 
   // helpers
   template <typename T>
@@ -1615,6 +1617,23 @@ inline bool HfFilterHelper::isSelectedKaon4Charm3Prong(const T& track)
   return true;
 }
 
+/// Method to check if charm candidates has mass between sideband limits
+/// \param massHypo1 is the array for the candidate D daughter momentum after reconstruction of secondary vertex
+/// \param massHypo2 is the array for the candidate bachelor pion momentum after reconstruction of secondary vertex
+/// \param lowLimitSB  is the dca of the D daughter track
+/// \param upLimitSB  is the dca of the pion daughter track
+/// \return true if the candidate passes the mass selection.
+template <typename T1>
+inline bool isCharmHadronMassInSbRegions(T1 const& massHypo1, T1 const& massHypo2, const float& lowLimitSB, const float& upLimitSB)
+{
+
+  if ((massHypo1 < lowLimitSB || massHypo1 > upLimitSB) && (massHypo2 < lowLimitSB || massHypo2 > upLimitSB)) {
+    return false;
+  }
+
+  return true;
+}
+
 /// Update the TPC PID baesd on the spline of particles
 /// \param track is a track parameter
 /// \param pidSpecies is the particle species to be considered