[PWGLF] He3-L derived data analysis

PWGLF/DataModel/LFSlimHeLambda.h

@@ -0,0 +1,82 @@
+// Copyright 2019-2020 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 LFSlimNucleiTables.h
+/// \brief Slim nuclei tables
+///
+
+#ifndef PWGLF_DATAMODEL_LFSLIMNUCLEITABLES_H_
+#define PWGLF_DATAMODEL_LFSLIMNUCLEITABLES_H_
+
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+
+#include <Math/Vector4D.h>
+
+namespace o2::aod
+{
+namespace lfv0he3
+{
+DECLARE_SOA_COLUMN(Z, z, float);
+DECLARE_SOA_COLUMN(CentT0C, centT0C, float);
+} // namespace lfv0he3
+DECLARE_SOA_TABLE(LFEvents, "AOD", "LFEVENT", o2::soa::Index<>, lfv0he3::Z, lfv0he3::CentT0C);
+
+namespace lfv0he3
+{
+DECLARE_SOA_INDEX_COLUMN(LFEvent, lfEvent); // Collision ID for the event
+DECLARE_SOA_COLUMN(Pt, pt, float);
+DECLARE_SOA_COLUMN(Eta, eta, float);
+DECLARE_SOA_COLUMN(Phi, phi, float);
+DECLARE_SOA_COLUMN(Mass, mass, float);
+DECLARE_SOA_COLUMN(CosPA, cosPA, float);
+DECLARE_SOA_COLUMN(DCAxy, dcaXY, float);
+DECLARE_SOA_COLUMN(DCAz, dcaZ, float);
+DECLARE_SOA_COLUMN(TPCnCls, tpcNCls, int);
+DECLARE_SOA_COLUMN(ITSClusterSizes, itsClusterSizes, uint32_t);
+DECLARE_SOA_COLUMN(NsigmaTPC, nSigmaTPC, float);
+// DECLARE_SOA_COLUMN(NsigmaTPCproton, nSigmaTPCproton, float);
+DECLARE_SOA_COLUMN(DCAdaughters, dcaDaughters, float);
+DECLARE_SOA_COLUMN(DCAPVProton, dcaPVProton, float);
+DECLARE_SOA_COLUMN(DCAPVPion, dcaPVPion, float);
+DECLARE_SOA_COLUMN(V0Radius, v0Radius, float);
+DECLARE_SOA_COLUMN(Sign, sign, int8_t);
+} // namespace lfv0he3
+DECLARE_SOA_TABLE(LFHe3, "AOD", "LFHE3V0", lfv0he3::LFEventId, lfv0he3::Pt, lfv0he3::Eta, lfv0he3::Phi, lfv0he3::DCAxy, lfv0he3::DCAz, lfv0he3::TPCnCls, lfv0he3::ITSClusterSizes, lfv0he3::NsigmaTPC, lfv0he3::Sign);
+DECLARE_SOA_TABLE(LFLambda, "AOD", "LFLAMBDA", lfv0he3::LFEventId, lfv0he3::Pt, lfv0he3::Eta, lfv0he3::Phi, lfv0he3::Mass, lfv0he3::CosPA, lfv0he3::DCAdaughters, lfv0he3::DCAPVProton, lfv0he3::DCAPVPion, lfv0he3::V0Radius, lfv0he3::Sign);
+} // namespace o2::aod
+
+struct he3Candidate {
+  ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double>> momentum; // 4-momentum of the He3 candidate
+  float nSigmaTPC = -999.f;                                            // TPC nSigma for He3
+  float dcaXY = -999.f;
+  float dcaZ = -999.f;
+  int tpcNClsFound = 0;         // Number of TPC clusters found
+  int itsNCls = 0;              // Number of ITS clusters
+  uint32_t itsClusterSizes = 0; // ITS cluster sizes
+  int8_t sign = 0;              // Charge sign of the He3 candidate
+};
+
+struct lambdaCandidate {
+  ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double>> momentum;
+  float mass = -1.f;              // Lambda mass
+  float cosPA = -2.f;             // Cosine of pointing angle
+  float dcaV0Daughters = -999.f;  // DCA between V0 daughters
+  float dcaProtonToPV = -999.f;   // DCA of the proton to primary vertex
+  float dcaPionToPV = -999.f;     // DCA of the pion to primary vertex
+  float v0Radius = -1.f;          // V0 radius
+  float protonNSigmaTPC = -999.f; // Proton TPC nSigma
+  float pionNSigmaTPC = -999.f;
+  int8_t sign = 0; // Charge sign of the Lambda candidate
+};
+
+#endif // PWGLF_DATAMODEL_LFSLIMNUCLEITABLES_H_

PWGLF/TableProducer/Nuspex/he3LambdaAnalysis.cxx

@@ -9,6 +9,8 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
+#include "PWGLF/DataModel/LFSlimHeLambda.h"
+
 #include "Common/Core/trackUtilities.h"
 #include "Common/DataModel/Centrality.h"
 #include "Common/DataModel/EventSelection.h"
@@ -45,39 +47,6 @@ using namespace o2;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
 using namespace o2::constants::physics;
-
-namespace o2::aod
-{
-namespace lfv0he3
-{
-DECLARE_SOA_COLUMN(Z, z, float);
-DECLARE_SOA_COLUMN(CentT0C, centT0C, float);
-} // namespace lfv0he3
-DECLARE_SOA_TABLE(LFEvents, "AOD", "LFEVENT", o2::soa::Index<>, lfv0he3::Z, lfv0he3::CentT0C);
-
-namespace lfv0he3
-{
-DECLARE_SOA_INDEX_COLUMN(LFEvent, lfEvent); // Collision ID for the event
-DECLARE_SOA_COLUMN(Pt, pt, float);
-DECLARE_SOA_COLUMN(Eta, eta, float);
-DECLARE_SOA_COLUMN(Phi, phi, float);
-DECLARE_SOA_COLUMN(Mass, mass, float);
-DECLARE_SOA_COLUMN(CosPA, cosPA, float);
-DECLARE_SOA_COLUMN(DCAxy, dcaXY, float);
-DECLARE_SOA_COLUMN(DCAz, dcaZ, float);
-DECLARE_SOA_COLUMN(TPCnCls, tpcNCls, int);
-DECLARE_SOA_COLUMN(ITSClusterSizes, itsClusterSizes, uint32_t);
-DECLARE_SOA_COLUMN(NsigmaTPC, nSigmaTPC, float);
-DECLARE_SOA_COLUMN(DCAdaughters, dcaDaughters, float);
-DECLARE_SOA_COLUMN(DCAPVProton, dcaPVProton, float);
-DECLARE_SOA_COLUMN(DCAPVPion, dcaPVPion, float);
-DECLARE_SOA_COLUMN(V0Radius, v0Radius, float);
-DECLARE_SOA_COLUMN(Sign, sign, int8_t);
-} // namespace lfv0he3
-DECLARE_SOA_TABLE(LFHe3, "AOD", "LFHE3V0", lfv0he3::LFEventId, lfv0he3::Pt, lfv0he3::Eta, lfv0he3::Phi, lfv0he3::DCAxy, lfv0he3::DCAz, lfv0he3::TPCnCls, lfv0he3::ITSClusterSizes, lfv0he3::NsigmaTPC, lfv0he3::Sign);
-DECLARE_SOA_TABLE(LFLambda, "AOD", "LFLAMBDA", lfv0he3::LFEventId, lfv0he3::Pt, lfv0he3::Eta, lfv0he3::Phi, lfv0he3::Mass, lfv0he3::CosPA, lfv0he3::DCAdaughters, lfv0he3::DCAPVProton, lfv0he3::DCAPVPion, lfv0he3::V0Radius, lfv0he3::Sign);
-} // namespace o2::aod
-
 namespace
 {
 constexpr double betheBlochDefault[1][6]{{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}};
@@ -106,37 +75,14 @@ std::shared_ptr<TH2> hTPCnSigmaAll;
 std::shared_ptr<TH2> hTPCnSigmaHe3;
 std::shared_ptr<TH2> hArmenterosPodolanskiAll;
 std::shared_ptr<TH2> hArmenterosPodolanskiSelected;
-std::shared_ptr<TH2> hInvariantMass;
+std::shared_ptr<TH2> hInvariantMassUS;
+std::shared_ptr<TH2> hInvariantMassLS;
 
 }; // namespace
 
 using TracksFull = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::pidTPCFullPr, aod::pidTPCFullPi>;
 using CollisionsFull = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms>;
 
-struct he3Candidate {
-  ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double>> momentum; // 4-momentum of the He3 candidate
-  float nSigmaTPC;                                                     // TPC nSigma for He3
-  float dcaXY;
-  float dcaZ;
-  int tpcNClsFound;         // Number of TPC clusters found
-  int itsNCls;              // Number of ITS clusters
-  uint32_t itsClusterSizes; // ITS cluster sizes
-  int8_t sign;              // Charge sign of the He3 candidate
-};
-
-struct lambdaCandidate {
-  ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double>> momentum;
-  float mass;           // Lambda mass
-  float cosPA;          // Cosine of pointing angle
-  float dcaV0Daughters; // DCA between V0 daughters
-  float dcaProtonToPV;  // DCA of the proton to primary vertex
-  float dcaPionToPV;    // DCA of the pion to primary vertex
-  float v0Radius;
-  float protonNSigmaTPC; // Proton TPC nSigma
-  float pionNSigmaTPC;
-  int8_t sign; // Charge sign of the Lambda candidate
-};
-
 struct he3LambdaAnalysis {
 
   // Services
@@ -251,7 +197,8 @@ struct he3LambdaAnalysis {
     hArmenterosPodolanskiSelected = mRegistry.add<TH2>("hArmenterosPodolanskiSelected", "Armenteros-Podolanski Selected", {HistType::kTH2D, {{100, -1., 1.}, {100, 0., 0.5}}});
 
     constexpr double ConstituentsMass = o2::constants::physics::MassProton + o2::constants::physics::MassNeutron * 2 + o2::constants::physics::MassSigmaPlus;
-    hInvariantMass = mRegistry.add<TH2>("hInvariantMass", "Invariant Mass", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
+    hInvariantMassUS = mRegistry.add<TH2>("hInvariantMassUS", "Invariant Mass", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
+    hInvariantMassLS = mRegistry.add<TH2>("hInvariantMassLS", "Invariant Mass", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
 
     LOGF(info, "He3-Lambda analysis initialized");
   }
@@ -439,7 +386,7 @@ struct he3LambdaAnalysis {
     for (const auto& he3 : he3Candidates) {
       for (const auto& lambda : lambdaCandidates) {
         auto pairMomentum = lambda.momentum + he3.momentum; // Calculate invariant mass
-        hInvariantMass->Fill(pairMomentum.Pt(), pairMomentum.M());
+        (he3.sign * lambda.sign > 0 ? hInvariantMassLS : hInvariantMassUS)->Fill(pairMomentum.Pt(), pairMomentum.M());
       }
     }
   }

PWGLF/Tasks/Nuspex/CMakeLists.txt

@@ -145,4 +145,9 @@ o2physics_add_dpl_workflow(kink-pika
     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(he3-lambda-derived-analysis
+    SOURCES he3LambdaDerivedAnalysis.cxx
+    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+    COMPONENT_NAME Analysis)
+
 endif()

PWGLF/Tasks/Nuspex/he3LambdaDerivedAnalysis.cxx

@@ -0,0 +1,140 @@
+// Copyright 2019-2020 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.
+
+#include "PWGLF/DataModel/LFSlimHeLambda.h"
+
+#include <Framework/ASoAHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/runDataProcessing.h>
+
+#include <Math/Vector4D.h>
+
+#include <memory>
+
+namespace
+{
+std::shared_ptr<TH2> hInvariantMassUS[2];
+std::shared_ptr<TH2> hInvariantMassLS[2];
+std::shared_ptr<TH2> hRotationInvariantMassUS[2];
+std::shared_ptr<TH2> hRotationInvariantMassLS[2];
+std::shared_ptr<TH2> hRotationInvariantMassAntiLSeta[2];
+std::shared_ptr<TH2> hInvariantMassLambda[2];
+std::shared_ptr<TH2> hCosPALambda;
+std::shared_ptr<TH2> hNsigmaHe3;
+std::shared_ptr<TH2> hNsigmaProton;
+}; // namespace
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::constants::physics;
+
+struct he3LambdaDerivedAnalysis {
+  HistogramRegistry mRegistry{"He3LambdaDerivedAnalysis"};
+
+  Configurable<int> cfgNrotations{"cfgNrotations", 7, "Number of rotations for He3 candidates"};
+  Configurable<bool> cfgMirrorEta{"cfgMirrorEta", true, "Mirror eta for He3 candidates"};
+  Configurable<float> cfgMinCosPA{"cfgMinCosPA", 0.99, "Minimum cosPA for Lambda candidates"};
+  Configurable<float> cfgMaxNSigmaTPCHe3{"cfgMaxNSigmaTPCHe3", 2.0, "Maximum nSigmaTPC for He3 candidates"};
+  Configurable<float> cfgMinLambdaPt{"cfgMinLambdaPt", 0.4, "Minimum pT for Lambda candidates"};
+  Configurable<float> cfgMaxLambdaDeltaM{"cfgMaxLambdaDeltaM", 10.0e-3, "Maximum deltaM for Lambda candidates"};
+
+  void init(InitContext const&)
+  {
+    constexpr double ConstituentsMass = o2::constants::physics::MassProton + o2::constants::physics::MassNeutron * 2 + o2::constants::physics::MassSigmaPlus;
+    for (int i = 0; i < 2; ++i) {
+      hInvariantMassUS[i] = mRegistry.add<TH2>(Form("hInvariantMassUS%i", i), "Invariant Mass", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
+      hInvariantMassLS[i] = mRegistry.add<TH2>(Form("hInvariantMassLS%i", i), "Invariant Mass", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
+      hRotationInvariantMassUS[i] = mRegistry.add<TH2>(Form("hRotationInvariantMassUS%i", i), "Rotation Invariant Mass", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
+      hRotationInvariantMassLS[i] = mRegistry.add<TH2>(Form("hRotationInvariantMassLS%i", i), "Rotation Invariant Mass", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
+      hInvariantMassLambda[i] = mRegistry.add<TH2>(Form("hInvariantMassLambda%i", i), "Invariant Mass Lambda", {HistType::kTH2D, {{50, 0., 10.}, {30, o2::constants::physics::MassLambda0 - 0.015, o2::constants::physics::MassLambda0 + 0.015}}});
+      hRotationInvariantMassAntiLSeta[i] = mRegistry.add<TH2>(Form("hRotationInvariantMassAntiLSeta%i", i), "Rotation Invariant Mass Anti-Lambda", {HistType::kTH2D, {{45, 1., 10}, {100, ConstituentsMass - 0.05, ConstituentsMass + 0.05}}});
+    }
+    hCosPALambda = mRegistry.add<TH2>("hCosPALambda", "Cosine of Pointing Angle for Lambda", {HistType::kTH2D, {{50, 0., 10.}, {500, 0.9, 1.}}});
+    hNsigmaHe3 = mRegistry.add<TH2>("hNsigmaHe3", "nSigma TPC for He3", {HistType::kTH2D, {{100, -10., 10.}, {200, -5, 5.}}});
+    hNsigmaProton = mRegistry.add<TH2>("hNsigmaProton", "nSigma TPC for Proton", {HistType::kTH2D, {{100, -10., 10.}, {200, -5, 5.}}});
+  }
+
+  void processSameEvent(o2::aod::LFEvents::iterator const& collision, o2::aod::LFHe3 const& he3s, o2::aod::LFLambda const& lambdas)
+  {
+    std::vector<he3Candidate> he3Candidates;
+    he3Candidates.reserve(he3s.size());
+    std::vector<lambdaCandidate> lambdaCandidates;
+    lambdaCandidates.reserve(lambdas.size());
+    for (const auto& he3 : he3s) {
+      if (he3.lfEventId() != collision.globalIndex()) {
+        std::cout << "He3 candidate does not match event index, skipping." << std::endl;
+        return;
+      }
+      he3Candidate candidate;
+      candidate.momentum = ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double>>(he3.pt(), he3.eta(), he3.phi(), o2::constants::physics::MassHelium3);
+      candidate.nSigmaTPC = he3.nSigmaTPC();
+      candidate.dcaXY = he3.dcaXY();
+      candidate.dcaZ = he3.dcaZ();
+      candidate.tpcNClsFound = he3.tpcNCls();
+      candidate.itsNCls = he3.itsClusterSizes();
+      candidate.itsClusterSizes = he3.itsClusterSizes();
+      candidate.sign = he3.sign();
+      hNsigmaHe3->Fill(he3.pt() * he3.sign(), he3.nSigmaTPC());
+      if (std::abs(he3.nSigmaTPC()) > cfgMaxNSigmaTPCHe3) {
+        continue; // Skip candidates with nSigmaTPC outside range
+      }
+      he3Candidates.push_back(candidate);
+    }
+    for (const auto& lambda : lambdas) {
+      if (lambda.lfEventId() != collision.globalIndex()) {
+        std::cout << "Lambda candidate does not match event index, skipping." << std::endl;
+        return;
+      }
+      lambdaCandidate candidate;
+      candidate.momentum.SetCoordinates(lambda.pt(), lambda.eta(), lambda.phi(), o2::constants::physics::MassLambda0);
+      candidate.mass = lambda.mass();
+      candidate.cosPA = lambda.cosPA();
+      candidate.dcaV0Daughters = lambda.dcaDaughters();
+      hCosPALambda->Fill(lambda.pt(), candidate.cosPA);
+      // hNsigmaProton->Fill(lambda.pt() * lambda.sign(), lambda.protonNSigmaTPC());
+      hInvariantMassLambda[0]->Fill(lambda.pt(), lambda.mass());
+      if (candidate.cosPA < cfgMinCosPA || lambda.pt() < cfgMinLambdaPt ||
+          std::abs(lambda.mass() - o2::constants::physics::MassLambda0) > cfgMaxLambdaDeltaM) {
+        continue; // Skip candidates with low cosPA
+      }
+      hInvariantMassLambda[1]->Fill(lambda.pt(), lambda.mass());
+      lambdaCandidates.push_back(candidate);
+    }
+
+    for (const auto& he3 : he3Candidates) {
+      for (const auto& lambda : lambdaCandidates) {
+        auto pairMomentum = lambda.momentum + he3.momentum; // Calculate invariant mass
+        (he3.sign * lambda.sign > 0 ? hInvariantMassLS : hInvariantMassUS)[he3.sign > 0]->Fill(pairMomentum.Pt(), pairMomentum.M());
+      }
+      for (int iEta{0}; iEta <= cfgMirrorEta; ++iEta) {
+        for (int iR{0}; iR <= cfgNrotations; ++iR) {
+          auto he3Momentum = ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double>>(he3.momentum.Pt(), (1. - iEta * 2.) * he3.momentum.Eta(), he3.momentum.Phi() + TMath::Pi() * (0.75 + 0.5 * iR / cfgNrotations), he3.momentum.M());
+          for (const auto& lambda : lambdaCandidates) {
+            auto pairMomentum = lambda.momentum + he3Momentum; // Calculate invariant mass
+            (he3.sign * lambda.sign > 0 ? hRotationInvariantMassLS : hRotationInvariantMassUS)[he3.sign > 0]->Fill(pairMomentum.Pt(), pairMomentum.M());
+            if (he3.sign < 0 && lambda.sign < 0) {
+              hRotationInvariantMassAntiLSeta[iEta]->Fill(pairMomentum.Pt(), pairMomentum.M());
+            }
+          }
+        }
+      }
+    }
+  }
+  PROCESS_SWITCH(he3LambdaDerivedAnalysis, processSameEvent, "Process same event", true);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<he3LambdaDerivedAnalysis>(cfgc)};
+}