[PWGCF] Update in D0 analysis (MC,ML, SB)

PWGCF/FemtoUniverse/Core/FemtoUniverseSoftPionRemoval.h

@@ -0,0 +1,120 @@
+// 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.
+//
+// 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 FemtoUniverseSoftPionRemoval.h
+/// \brief FemtoUniverseSoftPionRemoval - Checking the soft pions from D* decay and removing them
+/// \author Katarzyna Gwiździel, WUT, katarzyna.gwizdziel@cern.ch
+
+#ifndef PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSESOFTPIONREMOVAL_H_
+#define PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSESOFTPIONREMOVAL_H_
+
+#include <vector>
+
+#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"
+#include "Framework/HistogramRegistry.h"
+
+namespace o2::analysis::femto_universe
+{
+
+/// \class FemtoUniverseSoftPionRemoval
+/// \brief Class taking care of removing soft pions from D* decays
+/// \tparam partOne Type of particle 1 (Track/D0/...)
+/// \tparam partTwo Type of particle 2 (Track/D0/...)
+template <o2::aod::femtouniverseparticle::ParticleType partOne, o2::aod::femtouniverseparticle::ParticleType partTwo>
+class FemtoUniverseSoftPionRemoval
+{
+ public:
+  /// Destructor
+  virtual ~FemtoUniverseSoftPionRemoval() = default;
+
+  /// Initialization of the QA histograms
+  /// \param registry HistogramRegistry
+  void init(HistogramRegistry* registry)
+  {
+    if (registry) {
+      mHistogramRegistry = registry;
+      mHistogramRegistry->add("SoftPion/softPionMassVsPt", "; M(K#pi#pi-K#pi); p_{T}", kTH2F, {{200, 0.0, 0.2}, {36, 0., 36.}});
+    }
+  }
+
+  /// Check whether a track is a soft pion from D* decay
+  /// \tparam Part Data type of the particle
+  /// \tparam Parts Data type of the collection of all particles
+  /// \param part1 Particle 1
+  /// \param part2 Particle 2
+  /// \param particles Collection of all particles passed to the task
+  /// \return Whether the track is a soft pion
+  template <typename Part, typename Parts>
+  bool isSoftPion(Part const& part1, Part const& part2, Parts const& particles, bool isD0Cand, bool isD0barCand, double sigma)
+  {
+    if constexpr (kPartOneType == o2::aod::femtouniverseparticle::ParticleType::kTrack && kPartTwoType == o2::aod::femtouniverseparticle::ParticleType::kD0) {
+      /// Track-D0 combination part1 is hadron and part2 is D0
+      if (part2.partType() != o2::aod::femtouniverseparticle::ParticleType::kD0) {
+        LOG(fatal) << "FemtoUniverseSoftPionRemoval: passed arguments don't agree with FemtoUniverseSoftPionRemoval instantiation! Please provide second argument kD0 candidate.";
+        return false;
+      }
+      // Getting D0 (part2) children
+      const auto& posChild = particles.iteratorAt(part2.index() - 2);
+      const auto& negChild = particles.iteratorAt(part2.index() - 1);
+      // Pion and kaon mass
+      double massPion = o2::constants::physics::MassPiPlus;
+      double massKaon = o2::constants::physics::MassKPlus;
+      // D* reconstruction
+      double pSum2 = std::pow(posChild.px() + negChild.px() + part1.px(), 2.0) + std::pow(posChild.py() + negChild.py() + part1.py(), 2.0) + std::pow(posChild.pz() + negChild.pz() + part1.pz(), 2.0);
+      // Energies of the daughters -> D0->K-pi+
+      double e1Pi = std::sqrt(std::pow(massPion, 2.0) + std::pow(posChild.px(), 2.0) + std::pow(posChild.py(), 2.0) + std::pow(posChild.pz(), 2.0));
+      double e1K = std::sqrt(std::pow(massKaon, 2.0) + std::pow(negChild.px(), 2.0) + std::pow(negChild.py(), 2.0) + std::pow(negChild.pz(), 2.0));
+      // Energies of the daughters -> D0bar->K+pi-
+      double e2Pi = std::sqrt(std::pow(massPion, 2.0) + std::pow(negChild.px(), 2.0) + std::pow(negChild.py(), 2.0) + std::pow(negChild.pz(), 2.0));
+      double e2K = std::sqrt(std::pow(massKaon, 2.0) + std::pow(posChild.px(), 2.0) + std::pow(posChild.py(), 2.0) + std::pow(posChild.pz(), 2.0));
+      // Soft pion energy
+      auto ePion = RecoDecay::e(massPion, part1.p());
+      // D* masses
+      double mDstar1 = std::sqrt(std::pow(e1Pi + e1K + ePion, 2.0) - pSum2);
+      double mDstar2 = std::sqrt(std::pow(e2Pi + e2K + ePion, 2.0) - pSum2);
+
+      bool isSoftPion = false;
+      double softPiMass = 0.14542; // pion mass in D*->D0pi decay
+      double lowMassLimitSoftPion = softPiMass - 3.0 * sigma;
+      double highMassLimitSoftPion = softPiMass + 3.0 * sigma;
+
+      if (isD0Cand) {
+        if (mDstar1 - part2.mLambda() > 0.) {
+          mHistogramRegistry->fill(HIST("SoftPion/softPionMassVsPt"), mDstar1 - part2.mLambda(), part2.pt());
+        }
+        if ((std::abs(mDstar1 - part2.mLambda()) > lowMassLimitSoftPion) && (std::abs(mDstar1 - part2.mLambda()) < highMassLimitSoftPion)) {
+          isSoftPion = true;
+        }
+      }
+
+      if (isD0barCand) {
+        if (mDstar2 - part2.mAntiLambda() > 0.) {
+          mHistogramRegistry->fill(HIST("SoftPion/softPionMassVsPt"), mDstar2 - part2.mAntiLambda(), part2.pt());
+        }
+        if ((std::abs(mDstar2 - part2.mAntiLambda()) > lowMassLimitSoftPion) && (std::abs(mDstar2 - part2.mAntiLambda()) < highMassLimitSoftPion)) {
+          isSoftPion = true;
+        }
+      }
+      return isSoftPion;
+    } else {
+      LOG(fatal) << "FemtoUniverseSoftPionRemoval: Combination of objects not defined - quitting!";
+      return false;
+    }
+  }
+
+ private:
+  HistogramRegistry* mHistogramRegistry;                                                ///< For QA output
+  static constexpr o2::aod::femtouniverseparticle::ParticleType kPartOneType = partOne; ///< Type of particle 1
+  static constexpr o2::aod::femtouniverseparticle::ParticleType kPartTwoType = partTwo; ///< Type of particle 2
+};
+} // namespace o2::analysis::femto_universe
+
+#endif // PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSESOFTPIONREMOVAL_H_

PWGCF/FemtoUniverse/DataModel/FemtoDerived.h

@@ -105,11 +105,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(Theta, theta, //! Compute the theta of the track
                            });
 DECLARE_SOA_DYNAMIC_COLUMN(Px, px, //! Compute the momentum in x in GeV/c
                            [](float pt, float phi) -> float {
-                             return pt * std::sin(phi);
+                             return pt * std::cos(phi);
                            });
 DECLARE_SOA_DYNAMIC_COLUMN(Py, py, //! Compute the momentum in y in GeV/c
                            [](float pt, float phi) -> float {
-                             return pt * std::cos(phi);
+                             return pt * std::sin(phi);
                            });
 DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, //! Compute the momentum in z in GeV/c
                            [](float pt, float eta) -> float {
@@ -123,7 +123,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! Compute the overall momentum in GeV/c
 DECLARE_SOA_COLUMN(Sign, sign, int8_t);                                                  //! Sign of the track charge
 DECLARE_SOA_COLUMN(TpcNClsFound, tpcNClsFound, uint8_t);                                 //! Number of TPC clusters
 DECLARE_SOA_COLUMN(TpcNClsCrossedRows, tpcNClsCrossedRows, uint8_t);                     //! Number of TPC crossed rows
-DECLARE_SOA_COLUMN(TPCFractionSharedCls, tpcFractionSharedCls, float);                   //! Number of TPC crossed rows
+DECLARE_SOA_COLUMN(TpcFractionSharedCls, tpcFractionSharedCls, float);                   //! Number of TPC crossed rows
 DECLARE_SOA_COLUMN(ItsNCls, itsNCls, uint8_t);                                           //! Number of ITS clusters
 DECLARE_SOA_COLUMN(ItsNClsInnerBarrel, itsNClsInnerBarrel, uint8_t);                     //! Number of ITS clusters in the inner barrel                             //! TPC signal
 DECLARE_SOA_DYNAMIC_COLUMN(TpcCrossedRowsOverFindableCls, tpcCrossedRowsOverFindableCls, //! Compute the number of crossed rows over findable TPC clusters
@@ -182,7 +182,7 @@ DECLARE_SOA_TABLE(FDExtParticles, "AOD", "FDEXTPARTICLE",
                   track::TPCNClsFindable,
                   femtouniverseparticle::TpcNClsCrossedRows,
                   track::TPCNClsShared,
-                  femtouniverseparticle::TPCFractionSharedCls,
+                  femtouniverseparticle::TpcFractionSharedCls,
                   track::TPCInnerParam,
                   femtouniverseparticle::ItsNCls,
                   femtouniverseparticle::ItsNClsInnerBarrel,
@@ -260,6 +260,8 @@ enum ParticleOriginMCTruth {
   kFake,              //! particle, that has NOT the PDG code of the current analysed particle
   kDaughterLambda,    //! Daughter from a Lambda decay
   kDaughterSigmaplus, //! Daughter from a Sigma^plus decay
+  kPrompt,            //! Orgin for D0/D0bar mesons
+  kNonPrompt,         //! Orgin for D0/D0bar mesons
   kNOriginMCTruthTypes
 };
 
@@ -271,7 +273,9 @@ static constexpr std::string_view ParticleOriginMCTruthName[kNOriginMCTruthTypes
   "_NotPrimary",
   "_Fake",
   "_DaughterLambda",
-  "DaughterSigmaPlus"};
+  "DaughterSigmaPlus",
+  "_Prompt",
+  "_NonPrompt"};
 
 /// Distinguished between reconstructed and truth
 enum MCType {