Merge branch 'AliceO2Group:master' into master

Author: Shunsuke-Kurita

.github/labeler.yml

@@ -23,6 +23,10 @@ infrastructure:
           - 'packaging/**'
           - 'pyproject.toml'
 
+datamodel:
+  - changed-files:
+      - any-glob-to-any-file: ['DataModel/**', '*/DataModel/**']
+
 dpg:
   - changed-files:
       - any-glob-to-any-file: ['DPG/**']

CODEOWNERS

@@ -15,7 +15,7 @@
 /Common/CCDB @alibuild @jgrosseo @iarsene @ekryshen @ddobrigk
 /Common/Tools/Multiplicity @alibuild @ddobrigk @victor-gonzalez
 /ALICE3 @alibuild @njacazio @hscheid
-/DPG @alibuild @chiarazampolli @noferini
+/DPG @alibuild @chiarazampolli @alcaliva @catalinristea
 /DPG/Tasks/AOTEvent @alibuild @ekryshen @strogolo @altsybee
 /DPG/Tasks/AOTTrack @alibuild @mfaggin @iouribelikov @njacazio
 /DPG/Tasks/TOF @alibuild @noferini @njacazio

Common/Core/PID/PIDTOF.h

@@ -28,17 +28,19 @@
 #include "TMath.h"
 #include "TGraph.h"
 #include "TFile.h"
+#include "TF2.h"
 
 // O2 includes
 #include "DataFormatsTOF/ParameterContainers.h"
 #include "Framework/Logger.h"
 #include "ReconstructionDataFormats/PID.h"
 #include "Framework/DataTypes.h"
+#include "CommonConstants/PhysicsConstants.h"
 
 namespace o2::pid::tof
 {
 
-// Utility values
+// Utility values (to remove!)
 static constexpr float kCSPEED = TMath::C() * 1.0e2f * 1.0e-12f; /// Speed of light in TOF units (cm/ps)
 static constexpr float kCSPEDDInv = 1.f / kCSPEED;               /// Inverse of the Speed of light in TOF units (ps/cm)
 static constexpr float defaultReturnValue = -999.f;              /// Default return value in case TOF measurement is not available
@@ -55,7 +57,7 @@ class Beta
   /// \param length Length in cm of the track
   /// \param tofSignal TOF signal in ps for the track
   /// \param collisionTime collision time in ps for the event of the track
-  static float GetBeta(const float length, const float tofSignal, const float collisionTime) { return length / (tofSignal - collisionTime) * kCSPEDDInv; }
+  static float GetBeta(const float length, const float tofSignal, const float collisionTime) { return length / (tofSignal - collisionTime) * o2::constants::physics::invLightSpeedCm2PS; }
 
   /// Gets the beta for the track of interest
   /// \param track Track of interest
@@ -139,6 +141,11 @@ class TOFResoParamsV2 : public o2::tof::Parameters<13>
 
   ~TOFResoParamsV2() = default;
 
+  template <o2::track::PID::ID pid>
+  float getResolution(const float, const float) const
+  {
+    return -1.f;
+  }
   // Momentum shift for charge calibration
   void setMomentumChargeShiftParameters(std::unordered_map<std::string, float> const& pars)
   {
@@ -254,14 +261,12 @@ class TOFResoParamsV2 : public o2::tof::Parameters<13>
 class TOFResoParamsV3 : public o2::tof::Parameters<13>
 {
  public:
-  TOFResoParamsV3() : Parameters(std::array<std::string, 13>{"TrkRes.Pi.P0", "TrkRes.Pi.P1", "TrkRes.Pi.P2", "TrkRes.Pi.P3", "time_resolution",
-                                                             "TrkRes.Ka.P0", "TrkRes.Ka.P1", "TrkRes.Ka.P2", "TrkRes.Ka.P3",
-                                                             "TrkRes.Pr.P0", "TrkRes.Pr.P1", "TrkRes.Pr.P2", "TrkRes.Pr.P3"},
+  TOFResoParamsV3() : Parameters(std::array<std::string, 13>{"time_resolution", "time_resolution", "time_resolution", "time_resolution", "time_resolution",
+                                                             "time_resolution", "time_resolution", "time_resolution", "time_resolution",
+                                                             "time_resolution", "time_resolution", "time_resolution", "time_resolution"},
                                  "TOFResoParamsV3")
   {
-    setParameters(std::array<float, 13>{0.008, 0.008, 0.002, 40.0, 60.0,
-                                        0.008, 0.008, 0.002, 40.0,
-                                        0.008, 0.008, 0.002, 40.0});
+    setParameters(std::array<float, 13>{60.0});
   } // Default constructor with default parameters
 
   ~TOFResoParamsV3() = default;
@@ -313,7 +318,7 @@ class TOFResoParamsV3 : public o2::tof::Parameters<13>
     }
   }
 
-  // Time shift for post calibration
+  // Time shift for post calibration to realign as a function of eta
   void setTimeShiftParameters(std::unordered_map<std::string, float> const& pars, bool positive)
   {
     std::string baseOpt = positive ? "TimeShift.Pos." : "TimeShift.Neg.";
@@ -367,13 +372,55 @@ class TOFResoParamsV3 : public o2::tof::Parameters<13>
     return gNegEtaTimeCorr->Eval(eta);
   }
 
+  void setResolutionParametrization(std::unordered_map<std::string, float> const& pars)
+  {
+    static constexpr std::array<const char*, 9> particleNames = {"El", "Mu", "Pi", "Ka", "Pr", "De", "Tr", "He", "Al"};
+    for (int i = 0; i < 9; ++i) {
+      const std::string baseOpt = Form("tofResTrack.%s_", particleNames[i]);
+      // Check if a key begins with a string
+      for (const auto& [key, value] : pars) {
+        if (key.find(baseOpt) == 0) {
+          // Remove from the key the baseOpt
+          const std::string fun = key.substr(baseOpt.size());
+          mResolution[i] = new TF2(baseOpt.c_str(), fun.c_str(), 0., 20, -1, 1.);
+          LOG(info) << "Set the resolution function for " << particleNames[i] << " with formula " << mResolution[i]->GetFormula()->GetExpFormula();
+          break;
+        }
+      }
+    }
+    // Print a summary
+    for (int i = 0; i < 9; ++i) {
+      if (!mResolution[i]) {
+        LOG(info) << "Resolution function for " << particleNames[i] << " not provided, using default " << mDefaultResoParams[i];
+        mResolution[i] = new TF2(Form("tofResTrack.%s_Default", particleNames[i]), mDefaultResoParams[i], 0., 20, -1, 1.);
+      }
+      LOG(info) << "Resolution function for " << particleNames[i] << " is " << mResolution[i]->GetName() << " with formula " << mResolution[i]->GetFormula()->GetExpFormula();
+    }
+  }
+
+  template <o2::track::PID::ID pid>
+  float getResolution(const float p, const float eta) const
+  {
+    return mResolution[pid]->Eval(p, eta);
+  }
+
  private:
   // Charge calibration
   int mEtaN = 0; // Number of eta bins, 0 means no correction
   float mEtaStart = 0.f;
   float mEtaStop = 0.f;
   float mInvEtaWidth = 9999.f;
   std::vector<float> mContent;
+  std::array<TF2*, 9> mResolution{nullptr};
+  static constexpr std::array<const char*, 9> mDefaultResoParams{"14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)",
+                                                                 "14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)",
+                                                                 "14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)",
+                                                                 "42.66*TMath::Power((TMath::Max(x-0.417,0.1))*(1-0.4235*y*y),-0.7145)",
+                                                                 "99.46*TMath::Power((TMath::Max(x-0.447,0.1))*(1-0.4235*y*y),-0.8094)",
+                                                                 "216*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)",
+                                                                 "315*TMath::Power((TMath::Max(x-0.811,0.1))*(1-0.4235*y*y),-0.783)",
+                                                                 "157*TMath::Power((TMath::Max(x-0.556,0.1))*(1-0.4235*y*y),-0.783)",
+                                                                 "216*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)"};
 
   // Time shift for post calibration
   TGraph* gPosEtaTimeCorr = nullptr; /// Time shift correction for positive tracks
@@ -391,7 +438,7 @@ class ExpTimes
   static constexpr float mMassZSqared = mMassZ * mMassZ;                   /// (M/z)^2
 
   /// Computes the expected time of a track, given it TOF expected momentum
-  static float ComputeExpectedTime(const float tofExpMom, const float length) { return length * sqrt((mMassZSqared) + (tofExpMom * tofExpMom)) / (kCSPEED * tofExpMom); }
+  static float ComputeExpectedTime(const float tofExpMom, const float length) { return length * sqrt((mMassZSqared) + (tofExpMom * tofExpMom)) / (o2::constants::physics::LightSpeedCm2PS * tofExpMom); }
 
   /// Gets the expected signal of the track of interest under the PID assumption
   /// \param track Track of interest
@@ -401,7 +448,7 @@ class ExpTimes
       return defaultReturnValue;
     }
     if (track.trackType() == o2::aod::track::Run2Track) {
-      return ComputeExpectedTime(track.tofExpMom() * kCSPEDDInv, track.length());
+      return ComputeExpectedTime(track.tofExpMom() * o2::constants::physics::invLightSpeedCm2PS, track.length());
     }
     return ComputeExpectedTime(track.tofExpMom(), track.length());
   }
@@ -416,7 +463,7 @@ class ExpTimes
       return defaultReturnValue;
     }
     if (track.trackType() == o2::aod::track::Run2Track) {
-      return ComputeExpectedTime(track.tofExpMom() * kCSPEDDInv / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta())), track.length());
+      return ComputeExpectedTime(track.tofExpMom() * o2::constants::physics::invLightSpeedCm2PS / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta())), track.length());
     }
     LOG(debug) << "TOF exp. mom. " << track.tofExpMom() << " shifted = " << track.tofExpMom() / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta()));
     return ComputeExpectedTime(track.tofExpMom() / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta())), track.length()) + parameters.getTimeShift(track.eta(), track.sign());
@@ -432,9 +479,14 @@ class ExpTimes
   static float GetExpectedSigma(const ParamType& parameters, const TrackType& track, const float tofSignal, const float collisionTimeRes)
   {
     const float& mom = track.p();
+    const float& eta = track.eta();
     if (mom <= 0) {
       return -999.f;
     }
+    const float reso = parameters.template getResolution<id>(mom, eta);
+    if (reso > 0) {
+      return std::sqrt(reso * reso + parameters[4] * parameters[4] + collisionTimeRes * collisionTimeRes);
+    }
     if constexpr (id <= o2::track::PID::Pion) {
       LOG(debug) << "Using parameters for the pion hypothesis and ID " << id;
       const float dpp = parameters[0] + parameters[1] * mom + parameters[2] * mMassZ / mom; // mean relative pt resolution;

Common/Core/RecoDecay.h

@@ -1035,6 +1035,119 @@ struct RecoDecay {
     }
     return OriginType::None;
   }
+
+  /// based on getCharmHardronOrigin in order to extend general particle
+  /// Finding the origin (from charm hadronisation or beauty-hadron decay) of paritcle (b, c and others)
+  /// \param particlesMC  table with MC particles
+  /// \param particle  MC particle
+  /// \param searchUpToQuark if true tag origin based on charm/beauty quark otherwise on the presence of a b-hadron or c-hadron
+  /// \param idxBhadMothers optional vector of b-hadron indices (might be more than one in case of searchUpToQuark in case of beauty resonances)
+  /// \return an integer corresponding to the origin (0: none(others), 1: charm, 2: beauty) as in OriginType
+  template <typename T>
+  static int getParticleOrigin(const T& particlesMC,
+                               const typename T::iterator& particle,
+                               const bool searchUpToQuark = false,
+                               std::vector<int>* idxBhadMothers = nullptr)
+  {
+    int stage = 0; // mother tree level (just for debugging)
+
+    // vector of vectors with mother indices; each line corresponds to a "stage"
+    std::vector<std::vector<int64_t>> arrayIds{};
+    std::vector<int64_t> initVec{particle.globalIndex()};
+    arrayIds.push_back(initVec); // the first vector contains the index of the original particle
+    auto PDGParticle = std::abs(particle.pdgCode());
+    bool couldBeCharm = false;
+    if (PDGParticle / 100 == 4 || PDGParticle / 1000 == 4) {
+      couldBeCharm = true;
+    }
+    while (arrayIds[-stage].size() > 0) {
+      // vector of mother indices for the current stage
+      std::vector<int64_t> arrayIdsStage{};
+      for (auto& iPart : arrayIds[-stage]) { // check all the particles that were the mothers at the previous stage
+        auto particleMother = particlesMC.rawIteratorAt(iPart - particlesMC.offset());
+        if (particleMother.has_mothers()) {
+
+          // we break immediately if searchUpToQuark is false and the first mother is a parton (an hadron should never be the mother of a parton)
+          if (!searchUpToQuark) {
+            auto mother = particlesMC.rawIteratorAt(particleMother.mothersIds().front() - particlesMC.offset());
+            auto PDGParticleIMother = std::abs(mother.pdgCode()); // PDG code of the mother
+            if (PDGParticleIMother < 9 || (PDGParticleIMother > 20 && PDGParticleIMother < 38)) {
+              auto PDGPaticle = std::abs(particleMother.pdgCode());
+              if (
+                (PDGParticleIMother / 100 == 5 || // b mesons
+                 PDGParticleIMother / 1000 == 5)  // b baryons
+              ) {
+                return OriginType::NonPrompt; // beauty
+              }
+              if (
+                (PDGParticleIMother / 100 == 4 || // c mesons
+                 PDGParticleIMother / 1000 == 4)  // c baryons
+              ) {
+                return OriginType::Prompt; // charm
+              }
+              break;
+            }
+          }
+
+          for (auto iMother = particleMother.mothersIds().front(); iMother <= particleMother.mothersIds().back(); ++iMother) { // loop over the mother particles of the analysed particle
+            if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) {                       // if a mother is still present in the vector, do not check it again
+              continue;
+            }
+            auto mother = particlesMC.rawIteratorAt(iMother - particlesMC.offset());
+            // Check status code
+            auto motherStatusCode = std::abs(mother.getGenStatusCode());
+            auto PDGParticleIMother = std::abs(mother.pdgCode()); // PDG code of the mother
+            // Check mother's PDG code.
+            // printf("getMother: ");
+            // for (int i = stage; i < 0; i++) // Indent to make the tree look nice.
+            // printf(" ");
+            // printf("Stage %d: Mother PDG: %d, status: %d, Index: %d\n", stage, PDGParticleIMother, motherStatusCode, iMother);
+
+            if (searchUpToQuark) {
+              if (idxBhadMothers) {
+                if (PDGParticleIMother / 100 == 5 || // b mesons
+                    PDGParticleIMother / 1000 == 5)  // b baryons
+                {
+                  idxBhadMothers->push_back(iMother);
+                }
+              }
+              if (PDGParticleIMother == 5) {  // b quark
+                return OriginType::NonPrompt; // beauty
+              }
+              if (PDGParticleIMother == 4) { // c quark
+                return OriginType::Prompt;   // charm
+              }
+            } else {
+              if (
+                (PDGParticleIMother / 100 == 5 || // b mesons
+                 PDGParticleIMother / 1000 == 5)  // b baryons
+              ) {
+                if (idxBhadMothers) {
+                  idxBhadMothers->push_back(iMother);
+                }
+                return OriginType::NonPrompt; // beauty
+              }
+              if (
+                (PDGParticleIMother / 100 == 4 || // c mesons
+                 PDGParticleIMother / 1000 == 4)  // c baryons
+              ) {
+                couldBeCharm = true;
+              }
+            }
+            // add mother index in the vector for the current stage
+            arrayIdsStage.push_back(iMother);
+          }
+        }
+      }
+      // add vector of mother indices for the current stage
+      arrayIds.push_back(arrayIdsStage);
+      stage--;
+    }
+    if (couldBeCharm) {
+      return OriginType::Prompt; // charm
+    }
+    return OriginType::None;
+  }
 };
 
 /// Calculations using (pT, η, φ) coordinates, aka (transverse momentum, pseudorapidity, azimuth)