Phi expression column, use MathConstants.h (#7583)

jgrosseo · web-flow · commit 8dc418b7a6ae · 2021-11-12T08:24:43.000+01:00
diff --git a/Common/Constants/include/CommonConstants/MathConstants.h b/Common/Constants/include/CommonConstants/MathConstants.h
@@ -29,6 +29,8 @@ constexpr float VeryBig = 1.f / Almost0;
 constexpr float PI = 3.14159274101257324e+00f;
 constexpr float TwoPI = 2.f * PI;
 constexpr float PIHalf = 0.5f * PI;
+constexpr float PIThird = PI / 3.0f;
+constexpr float PIQuarter = 0.25f * PI;
 constexpr float Rad2Deg = 180.f / PI;
 constexpr float Deg2Rad = PI / 180.f;
 
diff --git a/Detectors/AOD/src/StandaloneAODProducer.cxx b/Detectors/AOD/src/StandaloneAODProducer.cxx
@@ -138,7 +138,7 @@ void fillCollisionAndTrackTable()
       auto collCursor = collBuilder.cursor<o2::aod::Collisions>();
 
       TableBuilder trackBuilder;
-      auto trackCursor = trackBuilder.cursor<o2::aod::Tracks>();
+      auto trackCursor = trackBuilder.cursor<o2::aod::StoredTracks>();
 
       int index = 0;
       for (auto& v : *vertices) {
@@ -190,7 +190,7 @@ void fillCollisionAndTrackTable()
           std::array<float, 3> pxpypz;
           track->getPxPyPzGlo(pxpypz);
           trackCursor(0, index, 0 /* CORRECT THIS */, track->getX(), track->getAlpha(), track->getY(), track->getZ(), track->getSnp(), track->getTgl(),
-                      track->getPt() /*CHECK!!*/, track->getPhi(), pxpypz[0], pxpypz[1], pxpypz[2]);
+                      track->getPt() /*CHECK!!*/);
         }
         index++;
       }
diff --git a/Framework/Core/include/Framework/AnalysisDataModel.h b/Framework/Core/include/Framework/AnalysisDataModel.h
@@ -15,8 +15,11 @@
 #include "MathUtils/Utils.h"
 #include <cmath>
 #include "Framework/DataTypes.h"
+#include "CommonConstants/MathConstants.h"
 #include "CommonConstants/PhysicsConstants.h"
 
+using namespace o2::constants::math;
+
 namespace o2
 {
 namespace aod
@@ -89,29 +92,21 @@ using Collision = Collisions::iterator;
 namespace track
 {
 // TRACKPAR TABLE definition
-DECLARE_SOA_INDEX_COLUMN(Collision, collision);      //! Collision to which this track belongs
-DECLARE_SOA_COLUMN(TrackType, trackType, uint8_t);   //! Type of track. See enum TrackTypeEnum
-DECLARE_SOA_COLUMN(X, x, float);                     //!
-DECLARE_SOA_COLUMN(Alpha, alpha, float);             //!
-DECLARE_SOA_COLUMN(Y, y, float);                     //!
-DECLARE_SOA_COLUMN(Z, z, float);                     //!
-DECLARE_SOA_COLUMN(Snp, snp, float);                 //!
-DECLARE_SOA_COLUMN(Tgl, tgl, float);                 //!
-DECLARE_SOA_COLUMN(Signed1Pt, signed1Pt, float);     //! (sign of charge)/Pt in c/GeV. Use pt() and sign() instead
-DECLARE_SOA_EXPRESSION_COLUMN(RawPhi, phiraw, float, //! Raw Phi (not folded onto [0, 2pi)). Use phi() instead
-                              nasin(aod::track::snp) + aod::track::alpha);
-// FIXME: make expression column when conditional nodes are supported in Gandiva
-DECLARE_SOA_DYNAMIC_COLUMN(NormalizedPhi, phi, //! Phi of the track, in radians within [0, 2pi)
-                           [](float phi) -> float {
-                             constexpr float twopi = 2.0f * static_cast<float>(M_PI);
-                             if (phi < 0)
-                               phi += twopi;
-                             if (phi > twopi)
-                               phi -= twopi;
-                             return phi;
-                           });
+DECLARE_SOA_INDEX_COLUMN(Collision, collision);    //! Collision to which this track belongs
+DECLARE_SOA_COLUMN(TrackType, trackType, uint8_t); //! Type of track. See enum TrackTypeEnum
+DECLARE_SOA_COLUMN(X, x, float);                   //!
+DECLARE_SOA_COLUMN(Alpha, alpha, float);           //!
+DECLARE_SOA_COLUMN(Y, y, float);                   //!
+DECLARE_SOA_COLUMN(Z, z, float);                   //!
+DECLARE_SOA_COLUMN(Snp, snp, float);               //!
+DECLARE_SOA_COLUMN(Tgl, tgl, float);               //!
+DECLARE_SOA_COLUMN(Signed1Pt, signed1Pt, float);   //! (sign of charge)/Pt in c/GeV. Use pt() and sign() instead
+DECLARE_SOA_EXPRESSION_COLUMN(Phi, phi, float,     //! Phi of the track, in radians within [0, 2pi)
+                              ifnode(nasin(aod::track::snp) + aod::track::alpha < 0.0f, nasin(aod::track::snp) + aod::track::alpha + TwoPI,
+                                     ifnode(nasin(aod::track::snp) + aod::track::alpha > TwoPI, nasin(aod::track::snp) + aod::track::alpha - TwoPI,
+                                            nasin(aod::track::snp) + aod::track::alpha)));
 DECLARE_SOA_EXPRESSION_COLUMN(Eta, eta, float, //! Pseudorapidity
-                              -1.f * nlog(ntan(0.25f * static_cast<float>(M_PI) - 0.5f * natan(aod::track::tgl))));
+                              -1.f * nlog(ntan(PIQuarter - 0.5f * natan(aod::track::tgl))));
 DECLARE_SOA_EXPRESSION_COLUMN(Pt, pt, float, //! Transverse momentum of the track in GeV/c
                               nabs(1.f / aod::track::signed1Pt));
 
@@ -140,7 +135,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, //! Momentum in z-direction in GeV/c
                            });
 
 DECLARE_SOA_EXPRESSION_COLUMN(P, p, float, //! Momentum in Gev/c
-                              0.5f * (ntan(0.25f * static_cast<float>(M_PI) - 0.5f * natan(aod::track::tgl)) + 1.f / ntan(0.25f * static_cast<float>(M_PI) - 0.5f * natan(aod::track::tgl))) / nabs(aod::track::signed1Pt));
+                              0.5f * (ntan(PIQuarter - 0.5f * natan(aod::track::tgl)) + 1.f / ntan(PIQuarter - 0.5f * natan(aod::track::tgl))) / nabs(aod::track::signed1Pt));
 
 // TRACKPARCOV TABLE definition
 DECLARE_SOA_COLUMN(SigmaY, sigmaY, float);        //! Covariance matrix
@@ -275,7 +270,6 @@ DECLARE_SOA_TABLE_FULL(StoredTracks, "Tracks", "AOD", "TRACK", //! On disk versi
                        track::X, track::Alpha,
                        track::Y, track::Z, track::Snp, track::Tgl,
                        track::Signed1Pt,
-                       track::NormalizedPhi<track::RawPhi>,
                        track::Px<track::Signed1Pt, track::Snp, track::Alpha>,
                        track::Py<track::Signed1Pt, track::Snp, track::Alpha>,
                        track::Pz<track::Signed1Pt, track::Tgl>,
@@ -285,7 +279,7 @@ DECLARE_SOA_EXTENDED_TABLE(Tracks, StoredTracks, "TRACK", //! Basic track proper
                            aod::track::Pt,
                            aod::track::P,
                            aod::track::Eta,
-                           aod::track::RawPhi);
+                           aod::track::Phi);
 
 DECLARE_SOA_TABLE_FULL(StoredTracksCov, "TracksCov", "AOD", "TRACKCOV", //! On disk version of the TracksCov table
                        track::SigmaY, track::SigmaZ, track::SigmaSnp, track::SigmaTgl, track::Sigma1Pt,
@@ -372,11 +366,11 @@ DECLARE_SOA_COLUMN(TrackTimeRes, trackTimeRes, float);
 DECLARE_SOA_DYNAMIC_COLUMN(Sign, sign,                                                       //! Sign of the track eletric charge
                            [](float signed1Pt) -> short { return (signed1Pt > 0) ? 1 : -1; });
 DECLARE_SOA_EXPRESSION_COLUMN(Eta, eta, float, //!
-                              -1.f * nlog(ntan(0.25f * static_cast<float>(M_PI) - 0.5f * natan(aod::fwdtrack::tgl))));
+                              -1.f * nlog(ntan(PIQuarter - 0.5f * natan(aod::fwdtrack::tgl))));
 DECLARE_SOA_EXPRESSION_COLUMN(Pt, pt, float, //!
                               nabs(1.f / aod::fwdtrack::signed1Pt));
 DECLARE_SOA_EXPRESSION_COLUMN(P, p, float, //!
-                              0.5f * (ntan(0.25f * static_cast<float>(M_PI) - 0.5f * natan(aod::fwdtrack::tgl)) + 1.f / ntan(0.25f * static_cast<float>(M_PI) - 0.5f * natan(aod::fwdtrack::tgl))) / nabs(aod::fwdtrack::signed1Pt));
+                              0.5f * (ntan(PIQuarter - 0.5f * natan(aod::fwdtrack::tgl)) + 1.f / ntan(PIQuarter - 0.5f * natan(aod::fwdtrack::tgl))) / nabs(aod::fwdtrack::signed1Pt));
 DECLARE_SOA_DYNAMIC_COLUMN(Px, px, //!
                            [](float pt, float phi) -> float {
                              return pt * std::cos(phi);
@@ -807,7 +801,7 @@ DECLARE_SOA_COLUMN(Vy, vy, float);
 DECLARE_SOA_COLUMN(Vz, vz, float);                                                      //! Z production vertex in cm
 DECLARE_SOA_COLUMN(Vt, vt, float);                                                      //! Production time
 DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi,                                                    //! Phi
-                           [](float px, float py) -> float { return static_cast<float>(M_PI) + std::atan2(-py, -px); });
+                           [](float px, float py) -> float { return PI + std::atan2(-py, -px); });
 DECLARE_SOA_DYNAMIC_COLUMN(ProducedByGenerator, producedByGenerator, //! True if particle produced by the generator (==TMCProcess::kPrimary); False if by the transport code
                            [](uint8_t flags) -> bool { return (flags & o2::aod::mcparticle::enums::ProducedByTransport) == 0x0; });
 DECLARE_SOA_DYNAMIC_COLUMN(FromBackgroundEvent, fromBackgroundEvent, //! Particle from background event
@@ -820,7 +814,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, //! True if par
                            [](uint8_t flags) -> bool { return (flags & o2::aod::mcparticle::enums::PhysicalPrimary) == o2::aod::mcparticle::enums::PhysicalPrimary; });
 
 // DECLARE_SOA_EXPRESSION_COLUMN(Phi, phi, float, //! Phi: NOTE this waits that the atan2 function is defined for expression columns
-//                               static_cast<float>(M_PI) + natan2(-aod::mcparticle::py, -aod::mcparticle::px));
+//                               PI + natan2(-aod::mcparticle::py, -aod::mcparticle::px));
 DECLARE_SOA_EXPRESSION_COLUMN(Eta, eta, float, //! Pseudorapidity, conditionally defined to avoid FPEs
                               ifnode((nsqrt(aod::mcparticle::px * aod::mcparticle::px +
                                             aod::mcparticle::py * aod::mcparticle::py +
diff --git a/Framework/Core/test/test_Expressions.cxx b/Framework/Core/test/test_Expressions.cxx
@@ -204,7 +204,7 @@ BOOST_AUTO_TEST_CASE(TestGandivaTreeCreation)
 BOOST_AUTO_TEST_CASE(TestConditionalExpressions)
 {
   // simple conditional
-  Filter cf = nabs(o2::aod::track::eta) < 1.0f && ifnode((o2::aod::track::pt < 1.0f), (o2::aod::track::phiraw > (float)(M_PI / 2.)), (o2::aod::track::phiraw < (float)(M_PI / 2.)));
+  Filter cf = nabs(o2::aod::track::eta) < 1.0f && ifnode((o2::aod::track::pt < 1.0f), (o2::aod::track::phi > (float)(M_PI / 2.)), (o2::aod::track::phi < (float)(M_PI / 2.)));
   auto cfspecs = createOperations(cf);
   BOOST_REQUIRE_EQUAL(cfspecs[0].left, (DatumSpec{1u, atype::BOOL}));
   BOOST_REQUIRE_EQUAL(cfspecs[0].right, (DatumSpec{2u, atype::BOOL}));
@@ -219,11 +219,11 @@ BOOST_AUTO_TEST_CASE(TestConditionalExpressions)
   BOOST_REQUIRE_EQUAL(cfspecs[2].right, (DatumSpec{LiteralNode::var_t{1.0f}, atype::FLOAT}));
   BOOST_REQUIRE_EQUAL(cfspecs[2].result, (DatumSpec{5u, atype::BOOL}));
 
-  BOOST_REQUIRE_EQUAL(cfspecs[3].left, (DatumSpec{std::string{"fRawPhi"}, typeid(o2::aod::track::RawPhi).hash_code(), atype::FLOAT}));
+  BOOST_REQUIRE_EQUAL(cfspecs[3].left, (DatumSpec{std::string{"fPhi"}, typeid(o2::aod::track::Phi).hash_code(), atype::FLOAT}));
   BOOST_REQUIRE_EQUAL(cfspecs[3].right, (DatumSpec{LiteralNode::var_t{(float)(M_PI / 2.)}, atype::FLOAT}));
   BOOST_REQUIRE_EQUAL(cfspecs[3].result, (DatumSpec{4u, atype::BOOL}));
 
-  BOOST_REQUIRE_EQUAL(cfspecs[4].left, (DatumSpec{std::string{"fRawPhi"}, typeid(o2::aod::track::RawPhi).hash_code(), atype::FLOAT}));
+  BOOST_REQUIRE_EQUAL(cfspecs[4].left, (DatumSpec{std::string{"fPhi"}, typeid(o2::aod::track::Phi).hash_code(), atype::FLOAT}));
   BOOST_REQUIRE_EQUAL(cfspecs[4].right, (DatumSpec{LiteralNode::var_t{(float)(M_PI / 2.)}, atype::FLOAT}));
   BOOST_REQUIRE_EQUAL(cfspecs[4].result, (DatumSpec{3u, atype::BOOL}));
 
@@ -237,13 +237,13 @@ BOOST_AUTO_TEST_CASE(TestConditionalExpressions)
 
   auto infield1 = o2::aod::track::Pt::asArrowField();
   auto infield2 = o2::aod::track::Eta::asArrowField();
-  auto infield3 = o2::aod::track::RawPhi::asArrowField();
+  auto infield3 = o2::aod::track::Phi::asArrowField();
   auto schema = std::make_shared<arrow::Schema>(std::vector{infield1, infield2, infield3});
   auto gandiva_tree = createExpressionTree(cfspecs, schema);
   auto gandiva_condition = makeCondition(gandiva_tree);
   auto gandiva_filter = createFilter(schema, gandiva_condition);
 
-  BOOST_CHECK_EQUAL(gandiva_tree->ToString(), "bool less_than(float absf((float) fEta), (const float) 1 raw(3f800000)) && if (bool less_than((float) fPt, (const float) 1 raw(3f800000))) { bool greater_than((float) fRawPhi, (const float) 1.5708 raw(3fc90fdb)) } else { bool less_than((float) fRawPhi, (const float) 1.5708 raw(3fc90fdb)) }");
+  BOOST_CHECK_EQUAL(gandiva_tree->ToString(), "bool less_than(float absf((float) fEta), (const float) 1 raw(3f800000)) && if (bool less_than((float) fPt, (const float) 1 raw(3f800000))) { bool greater_than((float) fPhi, (const float) 1.5708 raw(3fc90fdb)) } else { bool less_than((float) fPhi, (const float) 1.5708 raw(3fc90fdb)) }");
 
   // nested conditional
   Filter cfn = o2::aod::track::signed1Pt > 0.f && ifnode(std::move(*cf.node), nabs(o2::aod::track::x) > 1.0f, nabs(o2::aod::track::y) > 1.0f);
@@ -256,5 +256,5 @@ BOOST_AUTO_TEST_CASE(TestConditionalExpressions)
   auto gandiva_condition2 = makeCondition(gandiva_tree2);
   auto gandiva_filter2 = createFilter(schema2, gandiva_condition2);
   BOOST_REQUIRE_EQUAL(gandiva_tree2->ToString(),
-                      "bool greater_than((float) fSigned1Pt, (const float) 0 raw(0)) && if (bool less_than(float absf((float) fEta), (const float) 1 raw(3f800000)) && if (bool less_than((float) fPt, (const float) 1 raw(3f800000))) { bool greater_than((float) fRawPhi, (const float) 1.5708 raw(3fc90fdb)) } else { bool less_than((float) fRawPhi, (const float) 1.5708 raw(3fc90fdb)) }) { bool greater_than(float absf((float) fX), (const float) 1 raw(3f800000)) } else { bool greater_than(float absf((float) fY), (const float) 1 raw(3f800000)) }");
+                      "bool greater_than((float) fSigned1Pt, (const float) 0 raw(0)) && if (bool less_than(float absf((float) fEta), (const float) 1 raw(3f800000)) && if (bool less_than((float) fPt, (const float) 1 raw(3f800000))) { bool greater_than((float) fPhi, (const float) 1.5708 raw(3fc90fdb)) } else { bool less_than((float) fPhi, (const float) 1.5708 raw(3fc90fdb)) }) { bool greater_than(float absf((float) fX), (const float) 1 raw(3f800000)) } else { bool greater_than(float absf((float) fY), (const float) 1 raw(3f800000)) }");
 }
