[PWGLF] Add debug table to kink builder + mc info for kink sigma

PWGLF/DataModel/LFKinkDecayTables.h

@@ -15,10 +15,11 @@
 /// \author Francesco Mazzaschi <francesco.mazzaschi@cern.ch>
 ///
 
-#include "Framework/AnalysisDataModel.h"
-#include "Framework/ASoAHelpers.h"
 #include "Common/Core/RecoDecay.h"
 
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+
 #ifndef PWGLF_DATAMODEL_LFKINKDECAYTABLES_H_
 #define PWGLF_DATAMODEL_LFKINKDECAYTABLES_H_
 
@@ -77,6 +78,13 @@ DECLARE_SOA_DYNAMIC_COLUMN(MSigmaPlus, mSigmaPlus, //! mass under sigma plus hyp
                             float pzneut = pzmoth - pzch;
                             return RecoDecay::m(std::array{std::array{pxch, pych, pzch}, std::array{pxneut, pyneut, pzneut}}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionNeutral}); });
 
+DECLARE_SOA_DYNAMIC_COLUMN(MXiMinus, mXiMinus, //! mass under Xi minus hypothesis
+                           [](float pxmoth, float pymoth, float pzmoth, float pxch, float pych, float pzch) -> float {
+                            float pxneut = pxmoth - pxch;
+                            float pyneut = pymoth - pych;
+                            float pzneut = pzmoth - pzch;
+                            return RecoDecay::m(std::array{std::array{pxch, pych, pzch}, std::array{pxneut, pyneut, pzneut}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); });
+
 } // namespace kinkcand
 
 DECLARE_SOA_TABLE(KinkCands, "AOD", "KINKCANDS",
@@ -93,7 +101,24 @@ DECLARE_SOA_TABLE(KinkCands, "AOD", "KINKCANDS",
                   kinkcand::PtMoth<kinkcand::PxMoth, kinkcand::PyMoth>,
                   kinkcand::PtDaug<kinkcand::PxDaug, kinkcand::PyDaug>,
                   kinkcand::MSigmaMinus<kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug>,
-                  kinkcand::MSigmaPlus<kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug>);
+                  kinkcand::MSigmaPlus<kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug>,
+                  kinkcand::MXiMinus<kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug>);
+
+DECLARE_SOA_TABLE(KinkCandsUnbound, "AOD", "UBKINKCANDS",
+                  o2::soa::Index<>, kinkcand::XDecVtx, kinkcand::YDecVtx, kinkcand::ZDecVtx,
+                  kinkcand::MothSign, kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth,
+                  kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug,
+                  kinkcand::DcaMothPv, kinkcand::DcaDaugPv, kinkcand::DcaKinkTopo,
+
+                  // dynamic columns
+                  kinkcand::PxDaugNeut<kinkcand::PxMoth, kinkcand::PxDaug>,
+                  kinkcand::PyDaugNeut<kinkcand::PyMoth, kinkcand::PyDaug>,
+                  kinkcand::PzDaugNeut<kinkcand::PzMoth, kinkcand::PzDaug>,
+                  kinkcand::PtMoth<kinkcand::PxMoth, kinkcand::PyMoth>,
+                  kinkcand::PtDaug<kinkcand::PxDaug, kinkcand::PyDaug>,
+                  kinkcand::MSigmaMinus<kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug>,
+                  kinkcand::MSigmaPlus<kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug>,
+                  kinkcand::MXiMinus<kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug>);
 
 } // namespace o2::aod
 

PWGLF/TableProducer/Common/kinkBuilder.cxx

@@ -13,30 +13,31 @@
 /// \brief Builder task for kink decay topologies using ITS standalone tracks for the mother
 /// \author Francesco Mazzaschi <francesco.mazzaschi@cern.ch>
 
-#include <memory>
-#include <string>
-#include <array>
-#include <vector>
-#include <algorithm>
+#include "PWGLF/DataModel/LFKinkDecayTables.h"
+#include "PWGLF/DataModel/LFParticleIdentification.h"
+#include "PWGLF/Utils/svPoolCreator.h"
 
-#include "Framework/runDataProcessing.h"
-#include "Framework/AnalysisTask.h"
-#include "Framework/AnalysisDataModel.h"
-#include "Framework/ASoAHelpers.h"
-#include "ReconstructionDataFormats/Track.h"
+#include "Common/DataModel/Centrality.h"
 #include "Common/DataModel/EventSelection.h"
 #include "Common/DataModel/Multiplicity.h"
-#include "Common/DataModel/Centrality.h"
-#include "DetectorsBase/Propagator.h"
-#include "DetectorsBase/GeometryManager.h"
-#include "DataFormatsParameters/GRPObject.h"
-#include "DataFormatsParameters/GRPMagField.h"
+
 #include "CCDB/BasicCCDBManager.h"
 #include "DCAFitter/DCAFitterN.h"
+#include "DataFormatsParameters/GRPMagField.h"
+#include "DataFormatsParameters/GRPObject.h"
+#include "DetectorsBase/GeometryManager.h"
+#include "DetectorsBase/Propagator.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+#include "ReconstructionDataFormats/Track.h"
 
-#include "PWGLF/DataModel/LFParticleIdentification.h"
-#include "PWGLF/Utils/svPoolCreator.h"
-#include "PWGLF/DataModel/LFKinkDecayTables.h"
+#include <algorithm>
+#include <array>
+#include <memory>
+#include <string>
+#include <vector>
 
 using namespace o2;
 using namespace o2::framework;
@@ -93,9 +94,12 @@
                   kHyperhelium4sigma };
 
   Produces<aod::KinkCands> outputDataTable;
+  Produces<aod::KinkCandsUnbound> outputDataTableUB;
+
   Service<o2::ccdb::BasicCCDBManager> ccdb;
 
   Configurable<int> hypoMoth{"hypoMoth", kSigmaMinus, "Mother particle hypothesis"};
+  Configurable<bool> fillDebugTable{"fillDebugTable", false, "If true, fill the debug table with all candidates unbound"};
   // Selection criteria
   Configurable<float> maxDCAMothToPV{"maxDCAMothToPV", 0.1, "Max DCA of the mother to the PV"};
   Configurable<float> minDCADaugToPV{"minDCADaugToPV", 0., "Min DCA of the daughter to the PV"};
@@ -204,7 +208,7 @@
     h2ClsMapPtMoth = qaRegistry.add<TH2>("h2ClsMapPtMoth", "; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis});
     h2ClsMapPtDaug = qaRegistry.add<TH2>("h2ClsMapPtDaug", "; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis});
 
     for (int i = 0; i < 5; i++) {
       mBBparamsDaug[i] = cfgBetheBlochParams->get("Daughter", Form("p%i", i));
     }
     mBBparamsDaug[5] = cfgBetheBlochParams->get("Daughter", "resolution");
@@ -213,8 +217,8 @@
   template <typename T>
   bool selectMothTrack(const T& candidate)
   {
     if (candidate.has_collision() && candidate.hasITS() && !candidate.hasTPC() && !candidate.hasTOF() && candidate.itsNCls() < 6 &&
         candidate.itsNClsInnerBarrel() == 3 && candidate.itsChi2NCl() < 36 && candidate.pt() > minPtMoth) {
       return true;
     }
     return false;
@@ -232,8 +236,8 @@
     }
 
     bool isGoodTPCCand = false;
     if (candidate.itsNClsInnerBarrel() == 0 && candidate.itsNCls() < 4 &&
         candidate.tpcNClsCrossedRows() > 0.8 * candidate.tpcNClsFindable() && candidate.tpcNClsFound() > nTPCClusMinDaug) {
       isGoodTPCCand = true;
     }
 
@@ -336,13 +340,13 @@
 
       // get last layer hitted by the mother and the first layer hitted by the daughter
       int lastLayerMoth = 0, firstLayerDaug = 0;
       for (int i = 0; i < 7; i++) {
         if (trackMoth.itsClusterMap() & (1 << i)) {
           lastLayerMoth = i;
         }
       }
 
       for (int i = 0; i < 7; i++) {
         if (trackDaug.itsClusterMap() & (1 << i)) {
           firstLayerDaug = i;
           break;
@@ -357,13 +361,13 @@
         continue;
       }
 
       for (int i = 0; i < 3; i++) {
         kinkCand.decVtx[i] -= kinkCand.primVtx[i];
       }
 
       propMothTrack.getPxPyPzGlo(kinkCand.momMoth);
       propDaugTrack.getPxPyPzGlo(kinkCand.momDaug);
       for (int i = 0; i < 3; i++) {
         kinkCand.momMoth[i] *= charge;
         kinkCand.momDaug[i] *= charge;
       }
@@ -373,7 +377,7 @@
       kinkCand.kinkAngle = std::acos(spKink / (pMoth * pDaug));
 
       std::array<float, 3> neutDauMom{0.f, 0.f, 0.f};
       for (int i = 0; i < 3; i++) {
         neutDauMom[i] = kinkCand.momMoth[i] - kinkCand.momDaug[i];
       }
 
@@ -430,11 +434,18 @@
     std::sort(kinkCandidates.begin(), kinkCandidates.end(), [](const kinkCandidate& a, const kinkCandidate& b) { return a.collisionID < b.collisionID; });
 
     for (const auto& kinkCand : kinkCandidates) {
-      outputDataTable(kinkCand.collisionID, kinkCand.mothTrackID, kinkCand.daugTrackID,
-                      kinkCand.decVtx[0], kinkCand.decVtx[1], kinkCand.decVtx[2],
-                      kinkCand.mothSign, kinkCand.momMoth[0], kinkCand.momMoth[1], kinkCand.momMoth[2],
-                      kinkCand.momDaug[0], kinkCand.momDaug[1], kinkCand.momDaug[2],
-                      kinkCand.dcaXYmoth, kinkCand.dcaXYdaug, kinkCand.dcaKinkTopo);
+      if (fillDebugTable) {
+        outputDataTableUB(kinkCand.decVtx[0], kinkCand.decVtx[1], kinkCand.decVtx[2],
+                          kinkCand.mothSign, kinkCand.momMoth[0], kinkCand.momMoth[1], kinkCand.momMoth[2],
+                          kinkCand.momDaug[0], kinkCand.momDaug[1], kinkCand.momDaug[2],
+                          kinkCand.dcaXYmoth, kinkCand.dcaXYdaug, kinkCand.dcaKinkTopo);
+      } else {
+        outputDataTable(kinkCand.collisionID, kinkCand.mothTrackID, kinkCand.daugTrackID,
+                        kinkCand.decVtx[0], kinkCand.decVtx[1], kinkCand.decVtx[2],
+                        kinkCand.mothSign, kinkCand.momMoth[0], kinkCand.momMoth[1], kinkCand.momMoth[2],
+                        kinkCand.momDaug[0], kinkCand.momDaug[1], kinkCand.momDaug[2],
+                        kinkCand.dcaXYmoth, kinkCand.dcaXYdaug, kinkCand.dcaKinkTopo);
+      }
     }
   }
 };

PWGLF/Tasks/Strangeness/sigmaminustask.cxx

@@ -13,11 +13,13 @@
 /// \brief Example of a simple task for the analysis of the Sigma-minus
 /// \author Francesco Mazzaschi <francesco.mazzaschi@cern.ch>
 
-#include "Framework/runDataProcessing.h"
-#include "Framework/AnalysisTask.h"
-#include "Common/DataModel/EventSelection.h"
 #include "PWGLF/DataModel/LFKinkDecayTables.h"
+
+#include "Common/DataModel/EventSelection.h"
 #include "Common/DataModel/PIDResponse.h"
+
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
 #include "ReconstructionDataFormats/PID.h"
 
 using namespace o2;
@@ -26,6 +28,7 @@ using namespace o2::framework::expressions;
 
 using TracksFull = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::pidTPCPi>;
 using CollisionsFull = soa::Join<aod::Collisions, aod::EvSel>;
+using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels>;
 
 struct sigmaminustask {
   // Histograms are defined with HistogramRegistry
@@ -42,18 +45,24 @@ struct sigmaminustask {
     const AxisSpec ptAxis{50, -10, 10, "#it{p}_{T} (GeV/#it{c})"};
     const AxisSpec nSigmaPiAxis{100, -5, 5, "n#sigma_{#pi}"};
     const AxisSpec sigmaMassAxis{100, 1.1, 1.4, "m (GeV/#it{c}^{2})"};
+    const AxisSpec xiMassAxis{100, 1.2, 1.6, "m_{#Xi} (GeV/#it{c}^{2})"};
+    const AxisSpec pdgAxis{10001, -5000, 5000, "PDG code"};
     const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"};
 
     // Event selection
     rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}});
-
     // Sigma-minus reconstruction
     rSigmaMinus.add("h2MassSigmaMinusPt", "h2MassSigmaMinusPt", {HistType::kTH2F, {ptAxis, sigmaMassAxis}});
+    rSigmaMinus.add("h2SigmaMassVsXiMass", "h2SigmaMassVsXiMass", {HistType::kTH2F, {xiMassAxis, sigmaMassAxis}});
     rSigmaMinus.add("h2NSigmaPiPt", "h2NSigmaPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}});
+
+    if (doprocessMC) {
+      // Add MC histograms if needed
+      rSigmaMinus.add("h2MassSigmaMinusPtMC", "h2MassSigmaMinusPtMC", {HistType::kTH2F, {ptAxis, sigmaMassAxis}});
+    }
   }
 
-  void process(soa::Join<aod::Collisions, aod::EvSels>::iterator const& collision,
-               aod::KinkCands const& KinkCands, TracksFull const&)
+  void processData(CollisionsFull::iterator const& collision, aod::KinkCands const& KinkCands, TracksFull const&)
   {
     if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) {
       return;
@@ -65,9 +74,54 @@ struct sigmaminustask {
         continue;
       }
       rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus());
+      rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus());
       rSigmaMinus.fill(HIST("h2NSigmaPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi());
     }
   }
+  PROCESS_SWITCH(sigmaminustask, processData, "Data processing", true);
+
+  void processMC(CollisionsFullMC::iterator const& collision, aod::KinkCands const& KinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const&, TracksFull const&)
+  {
+    if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) {
+      return;
+    }
+    rEventSelection.fill(HIST("hVertexZRec"), collision.posZ());
+    for (const auto& kinkCand : KinkCands) {
+      auto dauTrack = kinkCand.trackDaug_as<TracksFull>();
+      auto mothTrack = kinkCand.trackMoth_as<TracksFull>();
+      if (dauTrack.sign() != mothTrack.sign()) {
+        LOG(info) << "Skipping kink candidate with opposite sign daughter and mother: " << kinkCand.globalIndex();
+        continue; // Skip if the daughter has the opposite sign as the mother
+      }
+      if (abs(dauTrack.tpcNSigmaPi()) > cutNSigmaPi) {
+        continue;
+      }
+
+      rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus());
+      rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus());
+      rSigmaMinus.fill(HIST("h2NSigmaPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi());
+      // do MC association
+      auto mcLabSigma = trackLabelsMC.rawIteratorAt(mothTrack.globalIndex());
+      auto mcLabPiDau = trackLabelsMC.rawIteratorAt(dauTrack.globalIndex());
+      if (mcLabSigma.has_mcParticle() && mcLabPiDau.has_mcParticle()) {
+        auto mcTrackSigma = mcLabSigma.mcParticle_as<aod::McParticles>();
+        auto mcTrackPiDau = mcLabPiDau.mcParticle_as<aod::McParticles>();
+        if (!mcTrackPiDau.has_mothers()) {
+          continue;
+        }
+        for (auto& piMother : mcTrackPiDau.mothers_as<aod::McParticles>()) {
+          if (piMother.globalIndex() != mcTrackSigma.globalIndex()) {
+            continue;
+          }
+          if (std::abs(mcTrackSigma.pdgCode()) != 3112 || std::abs(mcTrackPiDau.pdgCode()) != 211) {
+            continue;
+          }
+          rSigmaMinus.fill(HIST("h2MassSigmaMinusPtMC"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus());
+        }
+      }
+    }
+  }
+  PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)