Corrected table filling for MC reconstructed particles, added mc decay radius column, added few resolution histograms

Author: morgmatt

PWGLF/DataModel/LFKinkDecayTables.h

@@ -59,7 +59,7 @@ DECLARE_SOA_COLUMN(MothPdgCode, mothPdgCode, int);   //! PDG code of the Sigma d
 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
 
@@ -152,7 +152,7 @@ DECLARE_SOA_TABLE(SlimKinkCandsMC, "AOD", "SLIMKINKCANDSMC",
                   kinkcand::NSigmaTPCPi, kinkcand::NSigmaTPCPr, kinkcand::NSigmaTPCKa,
                   kinkcand::NSigmaTOFPi, kinkcand::NSigmaTOFPr, kinkcand::NSigmaTOFKa,
                   kinkcand::MothPdgCode, kinkcand::DaugPdgCode,
-                  kinkcand::PtMC, kinkcand::MassMC);
+                  kinkcand::PtMC, kinkcand::MassMC, kinkcand::DecayRadiusMC);
 
 
 

PWGLF/Tasks/Strangeness/sigmaminustask.cxx

@@ -53,13 +53,16 @@ struct sigmaminustask {
   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
@@ -71,6 +74,11 @@ struct sigmaminustask {
       // 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}});
     }
   }
 
@@ -151,7 +159,19 @@ struct sigmaminustask {
             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) {
@@ -163,7 +183,7 @@ struct sigmaminustask {
                               dauTrack.tpcNSigmaPi(), dauTrack.tpcNSigmaPr(), dauTrack.tpcNSigmaKa(),
                               dauTrack.tofNSigmaPi(), dauTrack.tofNSigmaPr(), dauTrack.tofNSigmaKa(),
                               mcTrackSigma.pdgCode(), mcTrackPiDau.pdgCode(),
-                              kinkCand.ptMoth(), kinkCand.mSigmaMinus());
+                              MotherpTMC, MotherMassMC, decayRadiusMC);
             }
           }
         } // MC association and selection
@@ -195,6 +215,7 @@ struct sigmaminustask {
         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);
 
@@ -208,7 +229,7 @@ struct sigmaminustask {
                           -999, -999, -999,
                           -999, -999, -999,
                           mcPart.pdgCode(), daug_pdg,
-                          mcPart.pt(), mcMass);
+                          mcPart.pt(), mcMass, mcDecayRadius);
       }
     }
   }