[PWGEM] Please consider the following formatting changes to #8484

.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/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)

Common/DataModel/PIDResponseITS.h

@@ -14,6 +14,7 @@
 /// \since  2024-11-12
 /// \author Nicolò Jacazio [email protected]
 /// \author Francesco Mazzaschi [email protected]
+/// \author Giorgio Alberto Lucia [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
 ///
@@ -54,19 +55,33 @@ struct ITSResponse {
     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);
+    if (id == o2::track::PID::Helium3 || id == o2::track::PID::Alpha) {
+      return (mITSRespParamsZ2[0] / (std::pow(bg, mITSRespParamsZ2[1])) + mITSRespParamsZ2[2]);
+    }
+    return (mITSRespParams[0] / (std::pow(bg, mITSRespParams[1])) + mITSRespParams[2]);
+  }
+
+  template <o2::track::PID::ID id>
+  static float expResolution(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;
+    float relRes = mResolutionParams[0] * std::erf((bg - mResolutionParams[1]) / mResolutionParams[2]);
+    return relRes;
   }
 
   template <o2::track::PID::ID id>
-  static float nSigmaITS(uint32_t itsClusterSizes, float momentum)
+  static float nSigmaITS(uint32_t itsClusterSizes, float momentum, float eta)
   {
     const float exp = expSignal<id>(momentum);
     const float average = averageClusterSize(itsClusterSizes);
-    const float resolution = mResolution * exp;
-    return (average - exp) / resolution;
+    const float coslInv = 1. / std::cosh(eta);
+    const float resolution = expResolution<id>(momentum) * exp;
+    return (average * coslInv - exp) / resolution;
   };
 
-  static void setParameters(float p0, float p1, float p2, float chargeFactor, float resolution)
+  static void setParameters(float p0, float p1, float p2, float p0_Z2, float p1_Z2, float p2_Z2, float p0_res, float p1_res, float p2_res)
   {
     if (mIsInitialized) {
       LOG(fatal) << "ITSResponse parameters already initialized";
@@ -75,79 +90,84 @@ struct ITSResponse {
     mITSRespParams[0] = p0;
     mITSRespParams[1] = p1;
     mITSRespParams[2] = p2;
-    mChargeFactor = chargeFactor;
-    mResolution = resolution;
+    mITSRespParamsZ2[0] = p0_Z2;
+    mITSRespParamsZ2[1] = p1_Z2;
+    mITSRespParamsZ2[2] = p2_Z2;
+    mResolutionParams[0] = p0_res;
+    mResolutionParams[1] = p1_res;
+    mResolutionParams[2] = p2_res;
   }
 
  private:
   static std::array<float, 3> mITSRespParams;
-  static float mChargeFactor;
-  static float mResolution;
+  static std::array<float, 3> mITSRespParamsZ2;
+  static std::array<float, 3> mResolutionParams;
   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;
+std::array<float, 3> ITSResponse::mITSRespParams = {1.1576, 1.684, 1.9453};
+std::array<float, 3> ITSResponse::mITSRespParamsZ2 = {2.8752, 1.1246, 5.0259};
+// relative resolution is modelled with an erf function: [0]*TMath::Erf((x-[1])/[2])
+std::array<float, 3> ITSResponse::mResolutionParams = {0.2431, -0.3293, 1.533};
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Electron>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Muon>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Pion>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Kaon>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Proton>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Deuteron>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Triton>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Helium3>(itsClusterSizes, momentum, eta);
                            });
 
 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);
+                           [](uint32_t itsClusterSizes, float momentum, float eta) -> float {
+                             return ITSResponse::nSigmaITS<o2::track::PID::Alpha>(itsClusterSizes, momentum, eta);
                            });
 
 // 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>;
+using ITSNSigmaEl = ITSNSigmaElImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaMu = ITSNSigmaMuImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaPi = ITSNSigmaPiImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaKa = ITSNSigmaKaImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaPr = ITSNSigmaPrImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaDe = ITSNSigmaDeImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaTr = ITSNSigmaTrImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaHe = ITSNSigmaHeImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
+using ITSNSigmaAl = ITSNSigmaAlImp<o2::aod::track::ITSClusterSizes, o2::aod::track::P, o2::aod::track::Eta>;
 
 } // namespace pidits
 } // namespace o2::aod

Common/TableProducer/CMakeLists.txt

@@ -129,3 +129,8 @@ o2physics_add_dpl_workflow(ese-table-producer
                     SOURCES eseTableProducer.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(mftmch-matching-data
+                    SOURCES match-mft-mch-data.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCCDB O2Physics::PWGDQCore O2Physics::EventFilteringUtils
+                    COMPONENT_NAME Analysis)