[PWGLF] Improvements in data model for strangeness, strangeness tof pid cleanup

PWGLF/DataModel/LFStrangenessPIDTables.h

@@ -1,4 +1,4 @@
 // 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.
 //
@@ -8,43 +8,163 @@
 // 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.
 
 // Defines TOF PID tables for strangeness.
 // Entries calculated per candidate, tables are joinable with v0/cascdata tables.
 
 #ifndef PWGLF_DATAMODEL_LFSTRANGENESSPIDTABLES_H_
 #define PWGLF_DATAMODEL_LFSTRANGENESSPIDTABLES_H_
 
-#include <cmath>
-#include "Framework/AnalysisDataModel.h"
+#include "PWGLF/DataModel/LFStrangenessTables.h"
+
 #include "Common/Core/RecoDecay.h"
+
 #include "CommonConstants/PhysicsConstants.h"
+#include "Framework/AnalysisDataModel.h"
+
+#include <cmath>
 
 namespace o2::aod
 {
 namespace dautrack
 {
-// ==== TPC INFORMATION ===
-DECLARE_SOA_COLUMN(TPCSignal, tpcSignal, float);     //! track TPC signal
+// ==== define packing helpers ===
+namespace packing
+{
+// define variables for packing
+static constexpr int nbins = (1 << 8 * sizeof(int8_t)) - 2;
+static constexpr int8_t overflowBin = nbins >> 1;
+static constexpr int8_t underflowBin = -(nbins >> 1);
+static constexpr float binned_max = 6.35;
+static constexpr float binned_min = -6.35;
+static constexpr float bin_width = (binned_max - binned_min) / nbins;
+static constexpr float underflow_return = -100.0f;
+static constexpr float overflow_return = +100.0f;
+
+// define helper function to do packing
+int8_t packInInt8(float nSigma)
+{
+  // calculate
+  if (nSigma <= binned_min)
+    return underflowBin;
+  if (nSigma >= binned_max)
+    return overflowBin;
+  if (nSigma >= 0) {
+    return static_cast<int8_t>((nSigma / bin_width) + 0.5f);
+  }
+  // automatic: this is the case in which nSigma < 0
+  return static_cast<int8_t>((nSigma / bin_width) - 0.5f);
+}
+
+// define helper function to do unpacking
+float unpackInt8(int8_t nSigma)
+{
+  if (nSigma == underflowBin) {
+    return underflow_return;
+  }
+  if (nSigma == overflowBin) {
+    return overflow_return;
+  }
+  return bin_width * nSigma;
+}
+
+} // namespace packing
+} // namespace dautrack
+
+namespace dautrack_legacy
+{
+// ==== LEGACY TPC INFORMATION (full size tables) ===
 DECLARE_SOA_COLUMN(TPCNSigmaEl, tpcNSigmaEl, float); //! Nsigma proton
 DECLARE_SOA_COLUMN(TPCNSigmaPi, tpcNSigmaPi, float); //! Nsigma proton
 DECLARE_SOA_COLUMN(TPCNSigmaKa, tpcNSigmaKa, float); //! Nsigma proton
 DECLARE_SOA_COLUMN(TPCNSigmaPr, tpcNSigmaPr, float); //! Nsigma proton
 DECLARE_SOA_COLUMN(TPCNSigmaHe, tpcNSigmaHe, float); //! Nsigma proton
+} // namespace dautrack_legacy
+
+namespace dautrack
+{
+// ==== COMPACT TPC INFORMATION (full size tables) ===
+DECLARE_SOA_COLUMN(TPCSignal, tpcSignal, float);                  //! track TPC signal
+DECLARE_SOA_COLUMN(PackedTPCNSigmaEl, packedTpcNSigmaEl, int8_t); //! Nsigma proton
+DECLARE_SOA_COLUMN(PackedTPCNSigmaPi, packedTpcNSigmaPi, int8_t); //! Nsigma proton
+DECLARE_SOA_COLUMN(PackedTPCNSigmaKa, packedTpcNSigmaKa, int8_t); //! Nsigma proton
+DECLARE_SOA_COLUMN(PackedTPCNSigmaPr, packedTpcNSigmaPr, int8_t); //! Nsigma proton
+
+DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaEl, tpcNSigmaEl, //! unpacked TPC nsigma
+                           [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); });
+DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaPi, tpcNSigmaPi, //! unpacked TPC nsigma
+                           [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); });
+DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaKa, tpcNSigmaKa, //! unpacked TPC nsigma
+                           [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); });
+DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaPr, tpcNSigmaPr, //! unpacked TPC nsigma
+                           [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); });
 
 // ==== TOF INFORMATION ===
+DECLARE_SOA_INDEX_COLUMN(DauTrackExtra, dauTrackExtra); //! point to daughter this TOF info belongs to
+DECLARE_SOA_INDEX_COLUMN(StraCollision, straCollision); //! point to collision associated with this track (not the V0/Casc)
 DECLARE_SOA_COLUMN(TOFSignal, tofSignal, float); //! track TOF signal
-DECLARE_SOA_COLUMN(TOFEvTime, tofEvTime, float); //! track TOF signal
-DECLARE_SOA_COLUMN(Length, length, float);       //! track TOF signal
+DECLARE_SOA_COLUMN(TOFEvTime, tofEvTime, float); //! event time
+DECLARE_SOA_COLUMN(Length, length, float);       //! track length (to assigned PV)
+DECLARE_SOA_COLUMN(TOFExpMom, tofExpMom, float); //! tof Exp Mom (to assigned PV)
+
+// dynamics with expected times
+DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimeEl, tofExpTimeEl, //! Expected time for the track to reach the TOF under the electron hypothesis
+                           [](float length, float tofExpMom) -> float {
+                             constexpr float massSquared = o2::constants::physics::MassElectron * o2::constants::physics::MassElectron;
+                             return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimePi, tofExpTimePi, //! Expected time for the track to reach the TOF under the pion hypothesis
+                           [](float length, float tofExpMom) -> float {
+                             constexpr float massSquared = o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged;
+                             return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimeKa, tofExpTimeKa, //! Expected time for the track to reach the TOF under the kaon hypothesis
+                           [](float length, float tofExpMom) -> float {
+                             constexpr float massSquared = o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged;
+                             return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared);
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimePr, tofExpTimePr, //! Expected time for the track to reach the TOF under the proton hypothesis
+                           [](float length, float tofExpMom) -> float {
+                             constexpr float massSquared = o2::constants::physics::MassProton * o2::constants::physics::MassProton;
+                             return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared);
+                           });
+
 } // namespace dautrack
 
-DECLARE_SOA_TABLE(DauTrackTPCPIDs, "AOD", "DAUTRACKTPCPID", // nsigma table (for analysis)
-                  dautrack::TPCSignal, dautrack::TPCNSigmaEl,
-                  dautrack::TPCNSigmaPi, dautrack::TPCNSigmaKa,
-                  dautrack::TPCNSigmaPr, dautrack::TPCNSigmaHe);
-DECLARE_SOA_TABLE(DauTrackTOFPIDs, "AOD", "DAUTRACKTOFPID", // raw table (for posterior TOF calculation)
+DECLARE_SOA_TABLE(DauTrackTPCPIDs_000, "AOD", "DAUTRACKTPCPID", // nsigma table (for analysis)
+                  dautrack::TPCSignal, dautrack_legacy::TPCNSigmaEl,
+                  dautrack_legacy::TPCNSigmaPi, dautrack_legacy::TPCNSigmaKa,
+                  dautrack_legacy::TPCNSigmaPr, dautrack_legacy::TPCNSigmaHe);
+
+DECLARE_SOA_TABLE_VERSIONED(DauTrackTPCPIDs_001, "AOD", "DAUTRACKTPCPID", 1, // nsigma table (for analysis)
+                            dautrack::TPCSignal,
+                            dautrack::PackedTPCNSigmaEl, dautrack::PackedTPCNSigmaPi,
+                            dautrack::PackedTPCNSigmaKa, dautrack::PackedTPCNSigmaPr,
+                            dautrack::TPCNSigmaEl<dautrack::PackedTPCNSigmaEl>,
+                            dautrack::TPCNSigmaPi<dautrack::PackedTPCNSigmaPi>,
+                            dautrack::TPCNSigmaKa<dautrack::PackedTPCNSigmaKa>,
+                            dautrack::TPCNSigmaPr<dautrack::PackedTPCNSigmaPr>);
+
+using DauTrackTPCPIDs = DauTrackTPCPIDs_001; // second gen: packed Nsigma, no He
+
+DECLARE_SOA_TABLE(DauTrackTOFPIDs_000, "AOD", "DAUTRACKTOFPID", // raw table (for posterior TOF calculation)
                   dautrack::TOFSignal, dautrack::TOFEvTime, dautrack::Length);
 
+DECLARE_SOA_TABLE_VERSIONED(DauTrackTOFPIDs_001, "AOD", "DAUTRACKTOFPID", 1, // raw table (for posterior TOF calculation)
+                            o2::soa::Index<>,
+                            dautrack::StraCollisionId, dautrack::DauTrackExtraId,
+                            dautrack::TOFSignal, dautrack::TOFEvTime,
+                            dautrack::Length, dautrack::TOFExpMom,
+                            dautrack::TOFExpTimeEl<dautrack::Length, dautrack::TOFExpMom>,
+                            dautrack::TOFExpTimePi<dautrack::Length, dautrack::TOFExpMom>,
+                            dautrack::TOFExpTimeKa<dautrack::Length, dautrack::TOFExpMom>,
+                            dautrack::TOFExpTimePr<dautrack::Length, dautrack::TOFExpMom>);
+
+using DauTrackTOFPIDs = DauTrackTOFPIDs_001; // second gen: with collision Id, with TOFExpMom
+
 namespace v0data
 {
 // define constants for NSigma operation

PWGLF/DataModel/LFStrangenessTables.h

@@ -63,6 +63,9 @@ DECLARE_SOA_DYNAMIC_COLUMN(EnergyCommonZNA, energyCommonZNA, //! get the total s
                            [](float value) -> float { return value; });
 DECLARE_SOA_DYNAMIC_COLUMN(EnergyCommonZNC, energyCommonZNC, //! get the total sum of the ZN A amplitudes
                            [](float value) -> float { return value; });
+
+// event time
+DECLARE_SOA_COLUMN(EventTime, eventTime, float); //! event time (FT0, TOF) for TOF PID (stored once per event)
 } // namespace stracollision
 
 //______________________________________________________
@@ -308,6 +311,8 @@ DECLARE_SOA_TABLE(StraStamps_000, "AOD", "STRASTAMPS", //! information for ID-in
                   bc::RunNumber, timestamp::Timestamp);
 DECLARE_SOA_TABLE_VERSIONED(StraStamps_001, "AOD", "STRASTAMPS", 1, //! information for ID-ing mag field if needed
                             bc::RunNumber, timestamp::Timestamp, bc::GlobalBC);
+DECLARE_SOA_TABLE(StraEvTimes, "AOD", "STRAEVTIMES", //! event time (FT0, TOF)
+                  stracollision::EventTime);
 
 using StraRawCents = StraRawCents_004;
 using StraCents = StraCents_001;

PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt

@@ -9,11 +9,21 @@
 # granted to it by virtue of its status as an Intergovernmental Organization
 # or submit itself to any jurisdiction.
 
+o2physics_add_dpl_workflow(stradautrackstpcpidconverter
+                    SOURCES stradautrackstpcpidconverter.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(stradautrackstofpidconverter
                     SOURCES stradautrackstofpidconverter.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(stradautrackstofpidconverter2
+                    SOURCES stradautrackstofpidconverter2.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(stradautracksextraconverter2
                     SOURCES stradautracksextraconverter2.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore

PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter.cxx

@@ -53,7 +53,7 @@ struct stradautrackstofpidconverter {
       lTOFEvTimes[casc.bachTrackExtraId()] = casc.bachTOFEventTime();
     }
     for (int ii = 0; ii < dauTracks.size(); ii++) {
-      dautracktofpids(lTOFSignals[ii], lTOFEvTimes[ii], lLengths[ii]);
+      dautracktofpids(-1, -1, lTOFSignals[ii], lTOFEvTimes[ii], lLengths[ii], 0.0f);
     }
   }
 };

PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter2.cxx

@@ -0,0 +1,66 @@
+// 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/LFStrangenessPIDTables.h"
+#include "PWGLF/DataModel/LFStrangenessTables.h"
+
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+
+using namespace o2;
+using namespace o2::framework;
+
+// converts DauTrackTOFPIDs_000 to _001
+struct stradautrackstofpidconverter2 {
+  Produces<aod::DauTrackTOFPIDs_001> dautracktofpids;
+  Produces<aod::StraEvTimes> straEvTimes;
+
+  void process(aod::StraCollisions const& collisions, soa::Join<aod::DauTrackExtras, aod::DauTrackTOFPIDs_000> const& dauTracks, soa::Join<aod::V0CollRefs, aod::V0Cores, aod::V0Extras> const& v0s)
+  {
+    // create new TOFPIDs
+    for (int ii = 0; ii < dauTracks.size(); ii++) {
+      auto dauTrack = dauTracks.rawIteratorAt(ii);
+      dautracktofpids(-1, -1, dauTrack.tofSignal(), dauTrack.tofEvTime(), dauTrack.length(), 0.0f);
+    }
+
+    // fill EvTimes (created simultaneously, so done in the same converter)
+    // recover as much as possible from the previous format
+    std::vector<double> collisionEventTime(collisions.size(), 0.0);
+    std::vector<int> collisionNtracks(collisions.size(), 0);
+
+    for (const auto& v0 : v0s) {
+      auto posTrackTOF = dauTracks.rawIteratorAt(v0.posTrackExtraId());
+      auto negTrackTOF = dauTracks.rawIteratorAt(v0.negTrackExtraId());
+      if (posTrackTOF.hasTOF()) {
+        collisionEventTime[v0.straCollisionId()] += posTrackTOF.tofEvTime();
+        collisionNtracks[v0.straCollisionId()]++;
+      }
+      if (negTrackTOF.hasTOF()) {
+        collisionEventTime[v0.straCollisionId()] += negTrackTOF.tofEvTime();
+        collisionNtracks[v0.straCollisionId()]++;
+      }
+    }
+    for (const auto& collision : collisions) {
+      if (collisionNtracks[collision.globalIndex()] > 0) {
+        collisionEventTime[collision.globalIndex()] /= static_cast<double>(collisionNtracks[collision.globalIndex()]);
+      } else {
+        collisionEventTime[collision.globalIndex()] = -1e+6; // undefined
+      }
+      straEvTimes(collisionEventTime[collision.globalIndex()]);
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<stradautrackstofpidconverter2>(cfgc)};
+}

PWGLF/TableProducer/Strangeness/Converters/stradautrackstpcpidconverter.cxx

@@ -0,0 +1,42 @@
+// 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/LFStrangenessPIDTables.h"
+#include "PWGLF/DataModel/LFStrangenessTables.h"
+
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+
+using namespace o2;
+using namespace o2::framework;
+
+// converts DauTrackTOFPIDs_000 to _001
+struct stradautrackstpcpidconverter {
+  Produces<aod::DauTrackTPCPIDs_001> dautrackpcpids;
+
+  void process(aod::DauTrackTPCPIDs_000 const& v000s)
+  {
+    for (int ii = 0; ii < v000s.size(); ii++) {
+      auto dauTrackTPCPID = v000s.rawIteratorAt(ii);
+      dautrackpcpids(dauTrackTPCPID.tpcSignal(),
+                     aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaEl()),
+                     aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaPi()),
+                     aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaKa()),
+                     aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaPr()));
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<stradautrackstpcpidconverter>(cfgc)};
+}