[PWGLF] addition of derived data path for lambda polarization (#9100)

ddobrigk · alibuild · web-flow · commit 71dcb8e5caec · 2024-12-22T20:38:55.000+01:00
Co-authored-by: ALICE Builder &lt;alibuild@users.noreply.github.com&gt;
diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h
@@ -229,8 +229,9 @@ DECLARE_SOA_TABLE(StraTPCQVs, "AOD", "STRATPCQVS", //! tpc Qvec
                   qvec::QvecBPosRe, qvec::QvecBPosIm, epcalibrationtable::QTPCR);
 DECLARE_SOA_TABLE(StraFT0CQVsEv, "AOD", "STRAFT0CQVSEv", //! events used to compute t0c Qvec
                   epcalibrationtable::TriggerEventEP);
-DECLARE_SOA_TABLE(StraZDCSP, "AOD", "STRAZDCSP", //! events used to compute the ZDC spectator plane
-                  spcalibrationtable::TriggerEventSP, spcalibrationtable::PsiZDCA, spcalibrationtable::PsiZDCC);
+DECLARE_SOA_TABLE(StraZDCSP, "AOD", "STRAZDCSP", //! ZDC SP information
+                  spcalibrationtable::TriggerEventSP,
+                  spcalibrationtable::PsiZDCA, spcalibrationtable::PsiZDCC);
 DECLARE_SOA_TABLE(StraStamps_000, "AOD", "STRASTAMPS", //! information for ID-ing mag field if needed
                   bc::RunNumber, timestamp::Timestamp);
 DECLARE_SOA_TABLE_VERSIONED(StraStamps_001, "AOD", "STRASTAMPS", 1, //! information for ID-ing mag field if needed
@@ -262,6 +263,8 @@ namespace dautrack
 {
 //______________________________________________________
 // Daughter track declarations for derived data analysis
+// These definitions are for the first version of the table
+// The latest version will inherit most properties from TracksExtra
 DECLARE_SOA_COLUMN(ITSChi2PerNcl, itsChi2PerNcl, float);        //! ITS chi2 per N cluster
 DECLARE_SOA_COLUMN(DetectorMap, detectorMap, uint8_t);          //! detector map for reference (see DetectorMapEnum)
 DECLARE_SOA_COLUMN(ITSClusterSizes, itsClusterSizes, uint32_t); //! ITS cluster sizes per layer
@@ -306,6 +309,18 @@ DECLARE_SOA_DYNAMIC_COLUMN(HasITSTracker, hasITSTracker, //! Flag to check if tr
                            [](uint8_t detectorMap, float itsChi2PerNcl) -> bool { return (detectorMap & o2::aod::track::ITS) ? (itsChi2PerNcl > -1e-3f) : false; });
 DECLARE_SOA_DYNAMIC_COLUMN(HasITSAfterburner, hasITSAfterburner, //! Flag to check if track is from ITS AB
                            [](uint8_t detectorMap, float itsChi2PerNcl) -> bool { return (detectorMap & o2::aod::track::ITS) ? (itsChi2PerNcl < -1e-3f) : false; });
+
+// sub-namespace for compatibility purposes
+namespace compatibility
+{                                                    // adds dynamics that ensure full backwards compatibility with previous getters
+DECLARE_SOA_DYNAMIC_COLUMN(TPCClusters, tpcClusters, //! number of TPC clusters
+                           [](uint8_t tpcNClsFindable, int8_t tpcNClsFindableMinusFound) -> int16_t { return (int16_t)tpcNClsFindable - tpcNClsFindableMinusFound; });
+DECLARE_SOA_DYNAMIC_COLUMN(TPCCrossedRows, tpcCrossedRows, //! Number of crossed TPC Rows
+                           [](uint8_t tpcNClsFindable, int8_t TPCNClsFindableMinusCrossedRows) -> int16_t { return (int16_t)tpcNClsFindable - TPCNClsFindableMinusCrossedRows; });
+DECLARE_SOA_DYNAMIC_COLUMN(ITSChi2PerNcl, itsChi2PerNcl, //! simple equivalent return
+                           [](float itsChi2NCl) -> float { return (float)itsChi2NCl; });
+} // namespace compatibility
+
 } // namespace dautrack
 
 DECLARE_SOA_TABLE(DauTrackExtras_000, "AOD", "DAUTRACKEXTRA", //! detector properties of decay daughters
@@ -337,10 +352,42 @@ DECLARE_SOA_TABLE_VERSIONED(DauTrackExtras_001, "AOD", "DAUTRACKEXTRA", 1, //! d
                             dautrack::HasITSTracker<dautrack::DetectorMap, dautrack::ITSChi2PerNcl>,
                             dautrack::HasITSAfterburner<dautrack::DetectorMap, dautrack::ITSChi2PerNcl>);
 
+DECLARE_SOA_TABLE_VERSIONED(DauTrackExtras_002, "AOD", "DAUTRACKEXTRA", 2, //! detector properties of decay daughters
+                            track::ITSChi2NCl,
+                            dautrack::DetectorMap, // here we don´t save everything so we simplify this
+                            track::ITSClusterSizes,
+                            track::TPCNClsFindable,
+                            track::TPCNClsFindableMinusFound,
+                            track::TPCNClsFindableMinusCrossedRows,
+
+                            // Dynamics for ITS matching TracksExtra
+                            track::v001::ITSNClsInnerBarrel<track::ITSClusterSizes>,
+                            track::v001::ITSClsSizeInLayer<track::ITSClusterSizes>,
+                            track::v001::ITSClusterMap<track::ITSClusterSizes>,
+                            track::v001::ITSNCls<track::ITSClusterSizes>,
+                            track::v001::IsITSAfterburner<track::v001::DetectorMap, track::ITSChi2NCl>,
+                            /*compatibility*/ dautrack::HasITSTracker<dautrack::DetectorMap, track::ITSChi2NCl>,
+                            /*compatibility*/ dautrack::HasITSAfterburner<dautrack::DetectorMap, track::ITSChi2NCl>,
+
+                            // dynamics for TPC tracking properties matching main data model
+                            track::TPCCrossedRowsOverFindableCls<track::TPCNClsFindable, track::TPCNClsFindableMinusCrossedRows>,
+                            track::TPCFoundOverFindableCls<track::TPCNClsFindable, track::TPCNClsFindableMinusFound>,
+                            track::TPCNClsFound<track::TPCNClsFindable, track::TPCNClsFindableMinusFound>,
+                            track::TPCNClsCrossedRows<track::TPCNClsFindable, track::TPCNClsFindableMinusCrossedRows>,
+                            /*compatibility*/ dautrack::compatibility::TPCClusters<track::TPCNClsFindable, track::TPCNClsFindableMinusFound>,
+                            /*compatibility*/ dautrack::compatibility::TPCCrossedRows<track::TPCNClsFindable, track::TPCNClsFindableMinusCrossedRows>,
+                            /*compatibility*/ dautrack::compatibility::ITSChi2PerNcl<track::ITSChi2NCl>,
+
+                            // dynamics to identify detectors
+                            dautrack::HasITS<dautrack::DetectorMap>,
+                            dautrack::HasTPC<dautrack::DetectorMap>,
+                            dautrack::HasTRD<dautrack::DetectorMap>,
+                            dautrack::HasTOF<dautrack::DetectorMap>);
+
 DECLARE_SOA_TABLE(DauTrackMCIds, "AOD", "DAUTRACKMCID", // index table when using AO2Ds
                   dautrack::ParticleMCId);
 
-using DauTrackExtras = DauTrackExtras_001;
+using DauTrackExtras = DauTrackExtras_002;
 using DauTrackExtra = DauTrackExtras::iterator;
 
 namespace motherParticle
diff --git a/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt
@@ -14,6 +14,11 @@ o2physics_add_dpl_workflow(stradautrackstofpidconverter
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(stradautracksextraconverter2
+                    SOURCES stradautracksextraconverter2.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(strarawcentsconverter
                     SOURCES strarawcentsconverter.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
diff --git a/PWGLF/TableProducer/Strangeness/Converters/stradautracksextraconverter2.cxx b/PWGLF/TableProducer/Strangeness/Converters/stradautracksextraconverter2.cxx
@@ -0,0 +1,41 @@
+// 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 "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/AnalysisDataModel.h"
+#include "PWGLF/DataModel/LFStrangenessTables.h"
+#include "PWGLF/DataModel/LFStrangenessPIDTables.h"
+
+using namespace o2;
+using namespace o2::framework;
+
+// Converts daughter TracksExtra from 1 to 2
+struct stradautracksextraconverter2 {
+  Produces<aod::DauTrackExtras_002> dauTrackExtras_002;
+
+  void process(aod::DauTrackExtras_001 const& dauTrackExtras_001)
+  {
+    for (auto& values : dauTrackExtras_001) {
+      dauTrackExtras_002(values.itsChi2PerNcl(),
+                         values.detectorMap(),
+                         values.itsClusterSizes(),
+                         static_cast<uint8_t>(0),   // findable (unknown in old format)
+                         -values.tpcClusters(),     // findable minus found: we know found
+                         -values.tpcCrossedRows()); // findable minus crossed rows: we know crossed rows
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<stradautracksextraconverter2>(cfgc)};
+}
diff --git a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx
@@ -483,8 +483,11 @@ struct strangederivedbuilder {
     for (auto const& tr : tracksExtra) {
       if (trackMap[tr.globalIndex()] >= 0) {
         dauTrackExtras(tr.itsChi2NCl(),
-                       tr.detectorMap(), tr.itsClusterSizes(),
-                       tr.tpcNClsFound(), tr.tpcNClsCrossedRows());
+                       tr.detectorMap(),
+                       tr.itsClusterSizes(),
+                       tr.tpcNClsFindable(),
+                       tr.tpcNClsFindableMinusFound(),
+                       tr.tpcNClsFindableMinusCrossedRows());
       }
     }
     // done!
@@ -574,8 +577,11 @@ struct strangederivedbuilder {
     for (auto const& tr : tracksExtra) {
       if (trackMap[tr.globalIndex()] >= 0) {
         dauTrackExtras(tr.itsChi2NCl(),
-                       tr.detectorMap(), tr.itsClusterSizes(),
-                       tr.tpcNClsFound(), tr.tpcNClsCrossedRows());
+                       tr.detectorMap(),
+                       tr.itsClusterSizes(),
+                       tr.tpcNClsFindable(),
+                       tr.tpcNClsFindableMinusFound(),
+                       tr.tpcNClsFindableMinusCrossedRows());
 
         // if the table has MC info
         if constexpr (requires { tr.mcParticle(); }) {
@@ -807,7 +813,8 @@ struct strangederivedbuilder {
   }
   void processZDCSP(soa::Join<aod::Collisions, aod::SPCalibrationTables>::iterator const& collision)
   {
-    StraZDCSP(collision.triggereventsp(), collision.psiZDCA(), collision.psiZDCC());
+    StraZDCSP(collision.triggereventsp(),
+              collision.psiZDCA(), collision.psiZDCC());
   }
   void processFT0MQVectors(soa::Join<aod::Collisions, aod::QvectorFT0Ms>::iterator const& collision)
   {
diff --git a/PWGLF/Tasks/Strangeness/lambdapolsp.cxx b/PWGLF/Tasks/Strangeness/lambdapolsp.cxx
@@ -54,13 +54,17 @@
 #include "DataFormatsParameters/GRPMagField.h"
 #include "CCDB/BasicCCDBManager.h"
 #include "PWGLF/DataModel/LFStrangenessTables.h"
+#include "PWGLF/DataModel/LFStrangenessPIDTables.h"
 #include "Common/DataModel/FT0Corrected.h"
 
 using namespace o2;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
 using std::array;
 
+using dauTracks = soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>;
+using v0Candidates = soa::Join<aod::V0CollRefs, aod::V0Cores, aod::V0Extras>;
+
 struct lambdapolsp {
 
   int mRunNumber;
@@ -427,6 +431,49 @@ struct lambdapolsp {
     return true;
   }
 
+  template <typename TV0>
+  bool isCompatible(TV0 const& v0, int pid /*0: lambda, 1: antilambda*/)
+  {
+    // checks if this V0 is compatible with the requested hypothesis
+
+    // de-ref track extras
+    auto posTrackExtra = v0.template posTrackExtra_as<dauTracks>();
+    auto negTrackExtra = v0.template negTrackExtra_as<dauTracks>();
+
+    // check for desired kinematics
+    if (pid == 0 && (v0.positivept() < cfgDaughPrPt || v0.negativept() < cfgDaughPiPt)) {
+      return false; // doesn´t pass lambda pT sels
+    }
+    if (pid == 1 && (v0.positivept() < cfgDaughPiPt || v0.negativept() < cfgDaughPrPt)) {
+      return false; // doesn´t pass antilambda pT sels
+    }
+    if (std::abs(v0.positiveeta()) > ConfDaughEta || std::abs(v0.negativeeta()) > ConfDaughEta) {
+      return false;
+    }
+
+    // check TPC tracking properties
+    if (posTrackExtra.tpcNClsCrossedRows() < 70 || negTrackExtra.tpcNClsCrossedRows() < 70) {
+      return false;
+    }
+    if (posTrackExtra.tpcNClsFound() < ConfDaughTPCnclsMin || negTrackExtra.tpcNClsFound() < ConfDaughTPCnclsMin) {
+      return false;
+    }
+    if (posTrackExtra.tpcCrossedRowsOverFindableCls() < 0.8 || negTrackExtra.tpcCrossedRowsOverFindableCls() < 0.8) {
+      return false;
+    }
+
+    // check TPC PID
+    if (pid == 0 && ((std::abs(posTrackExtra.tpcNSigmaPr()) > ConfDaughPIDCuts) || (std::abs(negTrackExtra.tpcNSigmaPi()) > ConfDaughPIDCuts))) {
+      return false;
+    }
+    if (pid == 1 && ((std::abs(posTrackExtra.tpcNSigmaPi()) > ConfDaughPIDCuts) || (std::abs(negTrackExtra.tpcNSigmaPr()) > ConfDaughPIDCuts))) {
+      return false;
+    }
+
+    // if we made it this far, it's good
+    return true;
+  }
+
   double GetPhiInRange(double phi)
   {
     double result = phi;
@@ -743,6 +790,136 @@ struct lambdapolsp {
     }
   }
   PROCESS_SWITCH(lambdapolsp, processData, "Process data", true);
+
+  // process function for derived data - mimics the functionality of the original data
+  void processDerivedData(soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraStamps, aod::StraZDCSP>::iterator const& collision, v0Candidates const& V0s, dauTracks const&)
+  {
+    //___________________________________________________________________________________________________
+    // event selection
+    if (!collision.sel8()) {
+      return;
+    }
+    auto centrality = collision.centFT0C();
+    if (!collision.triggereventsp()) { // provided by StraZDCSP
+      return;
+    }
+
+    if (additionalEvSel && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) {
+      return;
+    }
+    // histos.fill(HIST("hCentrality2"), centrality);
+    //  if (additionalEvSel2 && (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard))) {
+    if (additionalEvSel2 && (collision.trackOccupancyInTimeRange() > cfgMaxOccupancy || collision.trackOccupancyInTimeRange() < cfgMinOccupancy)) {
+      return;
+    }
+    // histos.fill(HIST("hCentrality3"), centrality);
+    if (additionalEvSel3 && (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(aod::evsel::kNoITSROFrameBorder))) {
+      return;
+    }
+
+    //___________________________________________________________________________________________________
+    // retrieve further info provided by StraZDCSP
+    auto psiZDCC = collision.psiZDCC();
+    auto psiZDCA = collision.psiZDCA();
+
+    // fill histograms
+    histos.fill(HIST("hCentrality"), centrality);
+    if (!checkwithpub) {
+      // histos.fill(HIST("hVtxZ"), collision.posZ());
+      histos.fill(HIST("hpRes"), centrality, (TMath::Cos(GetPhiInRange(psiZDCA - psiZDCC))));
+      if (QA) {
+        histos.fill(HIST("hpResSin"), centrality, (TMath::Sin(GetPhiInRange(psiZDCA - psiZDCC))));
+        histos.fill(HIST("hpCosPsiA"), centrality, (TMath::Cos(GetPhiInRange(psiZDCA))));
+        histos.fill(HIST("hpCosPsiC"), centrality, (TMath::Cos(GetPhiInRange(psiZDCC))));
+        histos.fill(HIST("hpSinPsiA"), centrality, (TMath::Sin(GetPhiInRange(psiZDCA))));
+        histos.fill(HIST("hpSinPsiC"), centrality, (TMath::Sin(GetPhiInRange(psiZDCC))));
+      }
+    }
+
+    //___________________________________________________________________________________________________
+    // loop over V0s as necessary
+    for (auto v0 : V0s) {
+      bool LambdaTag = isCompatible(v0, 0);
+      bool aLambdaTag = isCompatible(v0, 1);
+
+      if (LambdaTag == aLambdaTag)
+        continue;
+
+      if (!SelectionV0(collision, v0)) {
+        continue;
+      }
+
+      if (LambdaTag) {
+        Proton = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPr);
+        Pion = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPi);
+      }
+      if (aLambdaTag) {
+        Proton = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPr);
+        Pion = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPi);
+      }
+      Lambda = Proton + Pion;
+
+      ROOT::Math::Boost boost{Lambda.BoostToCM()};
+      fourVecDauCM = boost(Proton);
+      threeVecDauCM = fourVecDauCM.Vect();
+      // beamvector = ROOT::Math::XYZVector(0, 0, 1);
+      // eventplaneVec = ROOT::Math::XYZVector(collision.qFT0C(), collision.qFT0A(), 0); //this needs to be changed
+      // eventplaneVecNorm = eventplaneVec.Cross(beamvector); //z'
+      phiangle = TMath::ATan2(fourVecDauCM.Py(), fourVecDauCM.Px());
+      // double phiangledir = fourVecDauCM.Phi();
+
+      auto phiminuspsiC = GetPhiInRange(phiangle - psiZDCC);
+      auto phiminuspsiA = GetPhiInRange(phiangle - psiZDCA);
+      // histos.fill(HIST("hpsiApsiC"), psiZDCA, psiZDCC);
+      // histos.fill(HIST("hpsiApsiC"), GetPhiInRange(GetPhiInRange(phiangle) - GetPhiInRange(psiZDCA)), phiminuspsiA);
+      // histos.fill(HIST("hphiminuspsiA"), (phiminuspsiA));
+      // histos.fill(HIST("hphiminuspsiC"), (phiminuspsiC));
+      // auto cosThetaStar = eventplaneVecNorm.Dot(threeVecDauCM) / std::sqrt(threeVecDauCM.Mag2()) / std::sqrt(eventplaneVecNorm.Mag2());
+      auto cosThetaStar = fourVecDauCM.Pz() / fourVecDauCM.P(); // A0 correction
+      auto sinThetaStar = TMath::Sqrt(1 - (cosThetaStar * cosThetaStar));
+      auto PolC = TMath::Sin(phiminuspsiC);
+      auto PolA = TMath::Sin(phiminuspsiA);
+
+      // needed for corrections
+      auto sinPhiStar = TMath::Sin(GetPhiInRange(phiangle));
+      auto cosPhiStar = TMath::Cos(GetPhiInRange(phiangle));
+      auto sinThetaStarcosphiphiStar = sinThetaStar * TMath::Cos(2 * GetPhiInRange(Lambda.Phi() - phiangle)); // A2 correction
+      auto phiphiStar = GetPhiInRange(Lambda.Phi() - phiangle);
+
+      auto candmass = 0.0;
+      auto candpt = 0.0;
+      auto candeta = 0.0;
+
+      if (LambdaTag) {
+        candmass = v0.mLambda();
+        candpt = v0.pt();
+        candeta = v0.eta();
+
+        histos.fill(HIST("hSparseLambdaPolA"), candmass, candpt, candeta, PolA, centrality);
+        histos.fill(HIST("hSparseLambdaPolC"), candmass, candpt, candeta, PolC, centrality);
+        histos.fill(HIST("hSparseLambda_corr1a"), candmass, candpt, candeta, sinPhiStar, centrality);
+        histos.fill(HIST("hSparseLambda_corr1b"), candmass, candpt, candeta, cosPhiStar, centrality);
+        histos.fill(HIST("hSparseLambda_corr1c"), candmass, candpt, candeta, phiphiStar, centrality);
+        histos.fill(HIST("hSparseLambda_corr2a"), candmass, candpt, candeta, sinThetaStar, centrality);
+        histos.fill(HIST("hSparseLambda_corr2b"), candmass, candpt, candeta, sinThetaStarcosphiphiStar, centrality);
+      }
+
+      if (aLambdaTag) {
+        candmass = v0.mAntiLambda();
+        candpt = v0.pt();
+        candeta = v0.eta();
+
+        histos.fill(HIST("hSparseAntiLambdaPolA"), candmass, candpt, candeta, PolA, centrality);
+        histos.fill(HIST("hSparseAntiLambdaPolC"), candmass, candpt, candeta, PolC, centrality);
+        histos.fill(HIST("hSparseAntiLambda_corr1a"), candmass, candpt, candeta, sinPhiStar, centrality);
+        histos.fill(HIST("hSparseAntiLambda_corr1b"), candmass, candpt, candeta, cosPhiStar, centrality);
+        histos.fill(HIST("hSparseAntiLambda_corr1c"), candmass, candpt, candeta, phiphiStar, centrality);
+        histos.fill(HIST("hSparseAntiLambda_corr2a"), candmass, candpt, candeta, sinThetaStar, centrality);
+        histos.fill(HIST("hSparseAntiLambda_corr2b"), candmass, candpt, candeta, sinThetaStarcosphiphiStar, centrality);
+      }
+    } // end loop over V0s
+  }
+  PROCESS_SWITCH(lambdapolsp, processDerivedData, "Process derived data", false);
 };
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
