Add helicity angle analysis

Author: a-m-andrushko

PWGCF/FemtoUniverse/Core/FemtoUniverseMath.h

@@ -64,13 +64,14 @@ class FemtoUniverseMath
     return 0.5 * trackRelK.P();
   }
 
+  /// Boost particles from LAB Frame to Pair Rest Frame (for lambda daughters)
   /// \tparam T type of tracks
   /// \param part1 Particle 1
   /// \param mass1 Mass of particle 1
   /// \param part2 Particle 2
   /// \param mass2 Mass of particle 2
   template <typename T>
-  static float getthetastar(const T& part1, const float mass1, const T& part2, const float mass2)
+  static ROOT::Math::PxPyPzMVector boostPRF(const T& part1, const float mass1, const T& part2, const float mass2)
   {
     const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1);
     const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2);
@@ -88,10 +89,7 @@ class FemtoUniverseMath
     partOneCMS = boostPRF(partOneCMS);
     partTwoCMS = boostPRF(partTwoCMS);
 
-    const ROOT::Math::PtEtaPhiMVector partOneCMSGeo(partOneCMS);
-    const ROOT::Math::PtEtaPhiMVector partTwoCMSGeo(partTwoCMS);
-
-    return (partOneCMSGeo.Theta() - partTwoCMSGeo.Theta());
+    return partOneCMS;
   }
 
   /// Compute the qij of a pair of particles

PWGCF/FemtoUniverse/Tasks/femtoUniverseDebugV0.cxx

@@ -24,12 +24,14 @@
 #include "Framework/ASoAHelpers.h"
 #include "Framework/RunningWorkflowInfo.h"
 #include "Framework/StepTHn.h"
+#include "Framework/O2DatabasePDGPlugin.h"
 #include "DataFormatsParameters/GRPObject.h"
 
 #include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"
 #include "PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h"
 #include "PWGCF/FemtoUniverse/Core/FemtoUniverseEventHisto.h"
 #include "PWGCF/FemtoUniverse/Core/femtoUtils.h"
+#include "PWGCF/FemtoUniverse/Core/FemtoUniverseMath.h"
 
 using namespace o2;
 using namespace o2::analysis::femto_universe;
@@ -38,6 +40,9 @@ using namespace o2::framework::expressions;
 using namespace o2::soa;
 
 struct femtoUniverseDebugV0 {
+
+  Service<o2::framework::O2DatabasePDG> pdg;
+
   SliceCache cache;
 
   Configurable<int> ConfPDGCodeV0{"ConfPDGCodePartOne", 3122, "V0 - PDG code"};
@@ -71,21 +76,24 @@ struct femtoUniverseDebugV0 {
   /// Histogram output
   HistogramRegistry EventRegistry{"Event", {}, OutputObjHandlingPolicy::AnalysisObject};
   HistogramRegistry V0Registry{"FullV0QA", {}, OutputObjHandlingPolicy::AnalysisObject};
+  HistogramRegistry thetaRegistry{"ThetaQA", {}, OutputObjHandlingPolicy::AnalysisObject};
 
   void init(InitContext&)
   {
     eventHisto.init(&EventRegistry);
     posChildHistos.init(&V0Registry, ConfChildTempFitVarpTBins, ConfChildTempFitVarBins, false, ConfPDGCodeChildPos.value, true);
     negChildHistos.init(&V0Registry, ConfChildTempFitVarpTBins, ConfChildTempFitVarBins, false, ConfPDGCodeChildNeg, true);
     V0Histos.init(&V0Registry, ConfV0TempFitVarpTBins, ConfV0TempFitVarBins, false, ConfPDGCodeV0.value, true);
+
+    thetaRegistry.add("Theta/hTheta", " ; p (GeV/#it{c}); cos(#theta)", kTH2F, {{100, 0, 10}, {50, -5, 5}});
   }
 
-  /// Porduce QA plots for V0 selection in FemtoUniverse framework
+  /// Produce QA plots for V0 selection in FemtoUniverse framework
   void process(o2::aod::FdCollision const& col, FemtoFullParticles const& parts)
   {
     auto groupPartsOne = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
     eventHisto.fillQA(col);
-    for (auto& part : groupPartsOne) {
+    for (const auto& part : groupPartsOne) {
       if (!part.has_children()) {
         continue;
       }
@@ -95,14 +103,21 @@ struct femtoUniverseDebugV0 {
         LOG(warn) << "Indices of V0 children do not match";
         continue;
       }
-      // check cuts on V0 children
-      if ((posChild.partType() == uint8_t(aod::femtouniverseparticle::ParticleType::kV0Child) && (posChild.cut() & ConfCutChildPos) == ConfCutChildPos) &&
-          (negChild.partType() == uint8_t(aod::femtouniverseparticle::ParticleType::kV0Child) && (negChild.cut() & ConfCutChildNeg) == ConfCutChildNeg) &&
+
+      // Check cuts on V0 children
+      if (posChild.partType() == uint8_t(aod::femtouniverseparticle::ParticleType::kV0Child) &&
+          negChild.partType() == uint8_t(aod::femtouniverseparticle::ParticleType::kV0Child) &&
           isFullPIDSelected(posChild.pidCut(), posChild.p(), 999.f, ConfChildPosIndex.value, ConfChildnSpecies.value, ConfChildPIDnSigmaMax.value, ConfChildPosPidnSigmaMax.value, 1.f) &&
           isFullPIDSelected(negChild.pidCut(), negChild.p(), 999.f, ConfChildNegIndex.value, ConfChildnSpecies.value, ConfChildPIDnSigmaMax.value, ConfChildNegPidnSigmaMax.value, 1.f)) {
+        auto protonMass = pdg->Mass(2212);
+        auto pionMass = pdg->Mass(211);
+        auto protonBoosted = FemtoUniverseMath::boostPRF<decltype(posChild)>(posChild, protonMass, negChild, pionMass);
+        auto cosineTheta = (protonBoosted.Px() * part.px() + protonBoosted.Py() * part.py() + protonBoosted.Pz() * part.pz()) / (protonBoosted.P() * part.p());
+
         V0Histos.fillQA<false, true>(part);
         posChildHistos.fillQA<false, true>(posChild);
         negChildHistos.fillQA<false, true>(negChild);
+        thetaRegistry.fill(HIST("Theta/hTheta"), part.p(), cosineTheta);
       }
     }
   }