[PWGLF] Added output table to sigmaminustask

PWGLF/DataModel/LFKinkDecayTables.h

@@ -47,6 +47,20 @@ DECLARE_SOA_COLUMN(DcaMothPv, dcaMothPv, float);     //! DCA of the mother to th
 DECLARE_SOA_COLUMN(DcaDaugPv, dcaDaugPv, float);     //! DCA of the daughter kink to the primary vertex
 DECLARE_SOA_COLUMN(DcaKinkTopo, dcaKinkTopo, float); //! DCA of the kink topology
 
+DECLARE_SOA_COLUMN(NSigmaTPCPi, nSigmaTPCPi, float); //! Number of sigmas for the pion candidate from Sigma kink in TPC
+DECLARE_SOA_COLUMN(NSigmaTPCPr, nSigmaTPCPr, float); //! Number of sigmas for the proton candidate from Sigma kink in TPC
+DECLARE_SOA_COLUMN(NSigmaTPCKa, nSigmaTPCKa, float); //! Number of sigmas for the kaon candidate from Sigma kink in TPC
+DECLARE_SOA_COLUMN(NSigmaTOFPi, nSigmaTOFPi, float); //! Number of sigmas for the pion candidate from Sigma kink in TOF
+DECLARE_SOA_COLUMN(NSigmaTOFPr, nSigmaTOFPr, float); //! Number of sigmas for the proton candidate from Sigma kink in TOF
+DECLARE_SOA_COLUMN(NSigmaTOFKa, nSigmaTOFKa, float); //! Number of sigmas for the kaon candidate from Sigma kink in TOF
+
+// MC Columns
+DECLARE_SOA_COLUMN(MothPdgCode, mothPdgCode, int);       //! PDG code of the Sigma daughter
+DECLARE_SOA_COLUMN(DaugPdgCode, daugPdgCode, int);       //! PDG code of the kink daughter
+DECLARE_SOA_COLUMN(PtMC, ptMC, float);                   //! pT of the candidate in MC
+DECLARE_SOA_COLUMN(MassMC, massMC, float);               //! Invariant mass of the candidate in MC
+DECLARE_SOA_COLUMN(DecayRadiusMC, decayRadiusMC, float); //! Decay radius of the candidate in MC
+
 // DYNAMIC COLUMNS
 
 DECLARE_SOA_DYNAMIC_COLUMN(PxDaugNeut, pxDaugNeut, //! Px of the daughter neutral particle
@@ -120,6 +134,26 @@ DECLARE_SOA_TABLE(KinkCandsUnbound, "AOD", "UBKINKCANDS",
                   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(SlimKinkCands, "AOD", "SLIMKINKCANDS",
+                  kinkcand::XDecVtx, kinkcand::YDecVtx, kinkcand::ZDecVtx,
+                  kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth,
+                  kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug,
+                  kinkcand::DcaMothPv, kinkcand::DcaDaugPv, kinkcand::DcaKinkTopo,
+                  kinkcand::MothSign,
+                  kinkcand::NSigmaTPCPi, kinkcand::NSigmaTPCPr, kinkcand::NSigmaTPCKa,
+                  kinkcand::NSigmaTOFPi, kinkcand::NSigmaTOFPr, kinkcand::NSigmaTOFKa);
+
+DECLARE_SOA_TABLE(SlimKinkCandsMC, "AOD", "SLIMKINKCANDSMC",
+                  kinkcand::XDecVtx, kinkcand::YDecVtx, kinkcand::ZDecVtx,
+                  kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth,
+                  kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug,
+                  kinkcand::DcaMothPv, kinkcand::DcaDaugPv, kinkcand::DcaKinkTopo,
+                  kinkcand::MothSign,
+                  kinkcand::NSigmaTPCPi, kinkcand::NSigmaTPCPr, kinkcand::NSigmaTPCKa,
+                  kinkcand::NSigmaTOFPi, kinkcand::NSigmaTOFPr, kinkcand::NSigmaTOFKa,
+                  kinkcand::MothPdgCode, kinkcand::DaugPdgCode,
+                  kinkcand::PtMC, kinkcand::MassMC, kinkcand::DecayRadiusMC);
+
 } // namespace o2::aod
 
 #endif // PWGLF_DATAMODEL_LFKINKDECAYTABLES_H_

PWGLF/Tasks/Strangeness/sigmaminustask.cxx

@@ -26,11 +26,18 @@
 using namespace o2::framework;
 using namespace o2::framework::expressions;
 
-using TracksFull = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::pidTPCPi>;
+using TracksFull = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU,
+                             aod::pidTPCFullPi, aod::pidTPCFullPr, aod::pidTPCFullKa,
+                             aod::pidTOFFullPi, aod::pidTOFFullPr, aod::pidTOFFullKa>;
 using CollisionsFull = soa::Join<aod::Collisions, aod::EvSel>;
 using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSel>;
 
 struct sigmaminustask {
+
+  // Output Tables
+  Produces<aod::SlimKinkCands> outputDataTable;
+  Produces<aod::SlimKinkCandsMC> outputDataTableMC;
+
   // Histograms are defined with HistogramRegistry
   HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
   HistogramRegistry rSigmaMinus{"sigmaminus", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
@@ -39,18 +46,23 @@
   Configurable<float> cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"};
   Configurable<float> cutNSigmaPi{"cutNSigmaPi", 4, "NSigmaTPCPion"};
 
+  Configurable<bool> fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Kink candidates"};
+
   Preslice<aod::KinkCands> mPerCol = aod::track::collisionId;
 
   void init(InitContext const&)
   {
     // Axes
-    const AxisSpec ptAxis{50, -10, 10, "#it{p}_{T} (GeV/#it{c})"};
+    const AxisSpec ptAxis{100, -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]"};
 
+    const AxisSpec ptResolutionAxis{100, -0.5, 0.5, "#it{p}_{T}^{rec} - #it{p}_{T}^{gen} (GeV/#it{c})"};
+    const AxisSpec massResolutionAxis{100, -0.1, 0.1, "m_{rec} - m_{gen} (GeV/#it{c}^{2})"};
+
     // Event selection
     rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}});
     // Sigma-minus reconstruction
@@ -62,6 +74,11 @@
       // Add MC histograms if needed
       rSigmaMinus.add("h2MassPtMCRec", "h2MassPtMCRec", {HistType::kTH2F, {ptAxis, sigmaMassAxis}});
       rSigmaMinus.add("h2MassPtMCGen", "h2MassPtMCGen", {HistType::kTH2F, {ptAxis, sigmaMassAxis}});
+
+      rSigmaMinus.add("h2MassResolution_minus", "h2MassResolution_minus", {HistType::kTH2F, {sigmaMassAxis, massResolutionAxis}});
+      rSigmaMinus.add("h2PtResolution_minus", "h2PtResolution_minus", {HistType::kTH2F, {ptAxis, ptResolutionAxis}});
+      rSigmaMinus.add("h2MassResolution_plus", "h2MassResolution_plus", {HistType::kTH2F, {sigmaMassAxis, massResolutionAxis}});
+      rSigmaMinus.add("h2PtResolution_plus", "h2PtResolution_plus", {HistType::kTH2F, {ptAxis, ptResolutionAxis}});
     }
   }
 
@@ -71,14 +88,27 @@
       return;
     }
     rEventSelection.fill(HIST("hVertexZRec"), collision.posZ());
+
     for (const auto& kinkCand : KinkCands) {
       auto dauTrack = kinkCand.trackDaug_as<TracksFull>();
+
       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());
+
+      if (fillOutputTree) {
+        outputDataTable(kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx(),
+                        kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth(),
+                        kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug(),
+                        kinkCand.dcaMothPv(), kinkCand.dcaDaugPv(), kinkCand.dcaKinkTopo(),
+                        kinkCand.mothSign(),
+                        dauTrack.tpcNSigmaPi(), dauTrack.tpcNSigmaPr(), dauTrack.tpcNSigmaKa(),
+                        dauTrack.tofNSigmaPi(), dauTrack.tofNSigmaPr(), dauTrack.tofNSigmaKa());
+      }
     }
   }
   PROCESS_SWITCH(sigmaminustask, processData, "Data processing", true);
@@ -92,20 +122,23 @@
 
       rEventSelection.fill(HIST("hVertexZRec"), collision.posZ());
       auto kinkCandPerColl = KinkCands.sliceBy(mPerCol, collision.globalIndex());
+
       for (const auto& kinkCand : kinkCandPerColl) {
+
         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());
@@ -115,40 +148,91 @@
           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) {
+            if (std::abs(mcTrackSigma.pdgCode()) != 3112) {
               continue;
             }
+            if (std::abs(mcTrackPiDau.pdgCode()) != 211 && std::abs(mcTrackPiDau.pdgCode()) != 2212) {
+              continue;
+            }
+
+            float MotherMassMC = std::sqrt(piMother.e() * piMother.e() - piMother.p() * piMother.p());
+            float MotherpTMC = piMother.pt();
+            float decayRadiusMC = std::sqrt(mcTrackPiDau.vx() * mcTrackPiDau.vx() + mcTrackPiDau.vy() * mcTrackPiDau.vy());
+
             rSigmaMinus.fill(HIST("h2MassPtMCRec"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus());
+            if (mcTrackSigma.pdgCode() > 0) {
+              rSigmaMinus.fill(HIST("h2MassResolution_plus"), kinkCand.mSigmaMinus(), kinkCand.mSigmaMinus() - MotherMassMC);
+              rSigmaMinus.fill(HIST("h2PtResolution_plus"), kinkCand.ptMoth(), kinkCand.ptMoth() - MotherpTMC);
+            } else {
+              rSigmaMinus.fill(HIST("h2MassResolution_minus"), kinkCand.mSigmaMinus(), kinkCand.mSigmaMinus() - MotherMassMC);
+              rSigmaMinus.fill(HIST("h2PtResolution_minus"), kinkCand.ptMoth(), kinkCand.ptMoth() - MotherpTMC);
+            }
+
+            // fill the output table with Mc information
+            if (fillOutputTree) {
+              outputDataTableMC(kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx(),
+                                kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth(),
+                                kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug(),
+                                kinkCand.dcaMothPv(), kinkCand.dcaDaugPv(), kinkCand.dcaKinkTopo(),
+                                kinkCand.mothSign(),
+                                dauTrack.tpcNSigmaPi(), dauTrack.tpcNSigmaPr(), dauTrack.tpcNSigmaKa(),
+                                dauTrack.tofNSigmaPi(), dauTrack.tofNSigmaPr(), dauTrack.tofNSigmaKa(),
+                                mcTrackSigma.pdgCode(), mcTrackPiDau.pdgCode(),
+                                MotherpTMC, MotherMassMC, decayRadiusMC);
+            }
           }
-        }
-      }
-    }
+        } // MC association and selection
+      } // kink cand loop
+    } // collision loop
+
+    // Loop over all generated particles to fill MC histograms
     for (const auto& mcPart : particlesMC) {
-      if (std::abs(mcPart.pdgCode()) != 3112 || std::abs(mcPart.y()) > 0.5) {
+      if (std::abs(mcPart.pdgCode()) != 3112 || std::abs(mcPart.y()) > 1.0) { // only sigma mothers and rapidity cut
         continue;
       }
       if (!mcPart.has_daughters()) {
         continue; // Skip if no daughters
       }
       bool hasSigmaDaughter = false;
+      int daug_pdg = 0;
+      std::array<float, 3> secVtx;
+      std::array<float, 3> momDaug;
       for (const auto& daughter : mcPart.daughters_as<aod::McParticles>()) {
-        if (std::abs(daughter.pdgCode()) == 211) { // Sigma PDG code
+        if (std::abs(daughter.pdgCode()) == 211 || std::abs(daughter.pdgCode()) == 2212) { // Pi or proton daughter
           hasSigmaDaughter = true;
-          break; // Found a pi daughter, exit loop
+          secVtx = {daughter.vx(), daughter.vy(), daughter.vz()};
+          momDaug = {daughter.px(), daughter.py(), daughter.pz()};
+          daug_pdg = daughter.pdgCode();
+          break; // Found a daughter, exit loop
         }
       }
       if (!hasSigmaDaughter) {
-        continue; // Skip if no pi daughter found
+        continue; // Skip if no pi/proton daughter found
       }
       float mcMass = std::sqrt(mcPart.e() * mcPart.e() - mcPart.p() * mcPart.p());
+      float mcDecayRadius = std::sqrt(secVtx[0] * secVtx[0] + secVtx[1] * secVtx[1]);
       int sigmaSign = mcPart.pdgCode() > 0 ? 1 : -1; // Determine the sign of the Sigma
       rSigmaMinus.fill(HIST("h2MassPtMCGen"), sigmaSign * mcPart.pt(), mcMass);
+
+      // Fill output table with non reconstructed MC candidates
+      if (fillOutputTree) {
+        outputDataTableMC(-999, -999, -999,
+                          -999, -999, -999,
+                          -999, -999, -999,
+                          -999, -999, -999,
+                          sigmaSign,
+                          -999, -999, -999,
+                          -999, -999, -999,
+                          mcPart.pdgCode(), daug_pdg,
+                          mcPart.pt(), mcMass, mcDecayRadius);
+      }
     }
   }
+
   PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false);
 };