diff --git a/PWGCF/Femto/Core/collisionBuilder.h b/PWGCF/Femto/Core/collisionBuilder.h index 6b2660ae66b..9ffa9848882 100644 --- a/PWGCF/Femto/Core/collisionBuilder.h +++ b/PWGCF/Femto/Core/collisionBuilder.h @@ -81,7 +81,8 @@ struct ConfCollisionBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> sphericityMin{"sphericityMin", {}, "Minimum sphericity"}; o2::framework::Configurable> sphericityMax{"sphericityMax", {}, "Maximum sphericity"}; o2::framework::Configurable> triggers{"triggers", {}, "List of all triggers to be used"}; - o2::framework::Configurable qvecDetector{"qvecDetector", 0, "Detector used to estimate the Q-vector: 0 -> FT0C, 1 -> FT0A"}; + o2::framework::Configurable eventPlaneAngleDetector{"eventPlaneAngleDetector", 0, "Detector used to estimate the event plane angle: 0 -> FT0C, 1 -> FT0A"}; + o2::framework::Configurable qvecDetector{"qvecDetector", 0, "Detector used to estimate the Q-vector: 0 -> FT0C, 1 -> FT0A"}; o2::framework::Configurable qvecHarmonic{"qvecHarmonic", 2, "Harmonic n of the Q-vector and event plane angle Psi_n: 2 -> elliptic, 3 -> triangular"}; }; @@ -204,7 +205,7 @@ class CollisionSelection : public baseselection::BaseSelection @@ -244,13 +245,19 @@ class CollisionSelection : public baseselection::BaseSelection(config.qvecDetector.value); - if (mQvecDetector >= modes::QvecDetector::kQvecDetectorLast) { - LOG(fatal) << "Qvector Detector is not supported"; + mEventPlaneAngleDetector = static_cast(config.eventPlaneAngleDetector.value); + if (mEventPlaneAngleDetector >= modes::EventShapeDetector::kEventShapeDetectorLast) { + LOG(fatal) << "Detector " << static_cast(mEventPlaneAngleDetector) << " for event plane angle is not supported"; } + + mQvecDetector = static_cast(config.qvecDetector.value); + if (mQvecDetector >= modes::EventShapeDetector::kEventShapeDetectorLast) { + LOG(fatal) << "Detector " << static_cast(mQvecDetector) << " for q-vector is not supported"; + } + mQvecHarmonic = static_cast(config.qvecHarmonic.value); if (mQvecHarmonic < modes::QvecHarmonic::kN2 || mQvecHarmonic >= modes::QvecHarmonic::kQvecHarmonicLast) { - LOG(fatal) << "Qvector Harmonic is not supported"; + LOG(fatal) << "Harmonic " << static_cast(mQvecHarmonic) << " is not supported"; } this->addSelection(kSel8, collisionSelectionNames.at(kSel8), config.sel8.value); @@ -358,43 +365,40 @@ class CollisionSelection : public baseselection::BaseSelection void setQvector(T const& col) { + const auto index = static_cast(static_cast(mQvecHarmonic) - 2); // n=2 -> 0, n=3 -> 1 switch (mQvecDetector) { - case modes::QvecDetector::kFT0C: - mQvec = std::hypot(col.qvecFT0CReVec()[0], col.qvecFT0CImVec()[0]) * std::sqrt(col.sumAmplFT0C()); - break; - case modes::QvecDetector::kFT0A: - mQvec = std::hypot(col.qvecFT0AReVec()[0], col.qvecFT0AImVec()[0]) * std::sqrt(col.sumAmplFT0A()); + case modes::EventShapeDetector::kFT0C: + mQvec = computeReducedQvec(col.qvecFT0CReVec(), col.qvecFT0CImVec(), col.sumAmplFT0C(), index); break; - case modes::QvecDetector::kQvecDetectorLast: - LOG(fatal) << "Invalid Q-vector detector"; + case modes::EventShapeDetector::kFT0A: + mQvec = computeReducedQvec(col.qvecFT0AReVec(), col.qvecFT0AImVec(), col.sumAmplFT0A(), index); break; default: - LOG(fatal) << "Invalid Q-vector detector"; + LOG(fatal) << "Invalid detector for q-vector"; break; } } [[nodiscard]] float getQvector() const { return mQvec; } + /// \brief Event plane angle Psi_n in [0, 2pi/n) for the configured detector and harmonic template void setEventPlane(T const& col) { - auto harmonic = static_cast(mQvecHarmonic); - int index = static_cast(mQvecHarmonic) - 2; // get index in the qvector vector - switch (mQvecDetector) { - case modes::QvecDetector::kFT0C: - mEventPlane = RecoDecay::constrainAngle((std::atan2(col.qvecFT0CImVec()[index], col.qvecFT0CReVec()[index])) / harmonic, 0, harmonic); // constrain between 0 and 2pi/harmonic - break; - case modes::QvecDetector::kFT0A: - mEventPlane = RecoDecay::constrainAngle((std::atan2(col.qvecFT0AImVec()[index], col.qvecFT0AReVec()[index])) / harmonic, 0, harmonic); // constrain between 0 and 2pi/harmonic + const int harmonic = static_cast(mQvecHarmonic); + const auto index = static_cast(harmonic - 2); // n=2 -> 0, n=3 -> 1 + switch (mEventPlaneAngleDetector) { + case modes::EventShapeDetector::kFT0C: + mEventPlane = computeEventPlane(col.qvecFT0CReVec(), col.qvecFT0CImVec(), index, harmonic); break; - case modes::QvecDetector::kQvecDetectorLast: - LOG(fatal) << "Invalid Q-vector detector"; + case modes::EventShapeDetector::kFT0A: + mEventPlane = computeEventPlane(col.qvecFT0AReVec(), col.qvecFT0AImVec(), index, harmonic); break; default: - LOG(fatal) << "Invalid Q-vector detector"; + LOG(fatal) << "Invalid detector for event plane angle"; break; } } @@ -527,6 +531,30 @@ class CollisionSelection : public baseselection::BaseSelection(2. * lambda2 / (lambda1 + lambda2)); } + /// \brief |q_n| = |Q_n| * sqrt(M); returns 0 for non-positive amplitude sum to avoid NaN + template + static float computeReducedQvec(TRe const& re, TIm const& im, float sumAmpl, std::size_t index) + { + if (index >= re.size() || index >= im.size()) { + LOG(fatal) << "Requested Q-vector harmonic index " << index << " but only " << re.size() << " harmonics are stored"; + } + if (!(sumAmpl > 0.f)) { + return 0.f; + } + return static_cast(std::hypot(re[index], im[index]) * std::sqrt(sumAmpl)); + } + + /// \brief Psi_n = atan2(Im Q_n, Re Q_n) / n, constrained to [0, 2pi/n) + template + static float computeEventPlane(TRe const& re, TIm const& im, std::size_t index, int harmonic) + { + if (index >= re.size() || index >= im.size()) { + LOG(fatal) << "Requested Q-vector harmonic index " << index << " but only " << re.size() << " harmonics are stored"; + } + const double psi = std::atan2(static_cast(im[index]), static_cast(re[index])) / static_cast(harmonic); + return static_cast(RecoDecay::constrainAngle(psi, 0., static_cast(harmonic))); + } + // filter cuts float mVtxZMin = -12.f; float mVtxZMax = 12.f; @@ -546,7 +574,9 @@ class CollisionSelection : public baseselection::BaseSelection, kColHistLast> HistTable = { {kCentVsSphericity, o2::framework::HistType::kTH2F, "hCentVsSphericity", "Centrality vs Sphericity; Centrality (%); Sphericity"}, {kFT0AvsFT0C, o2::framework::HistType::kTH2F, "hFT0AvsFT0C", "FT0A centrality vs FT0C centrality; Centrality_{FT0A} (%); Centrality_{FT0C}"}, // event shape - {kQvector, o2::framework::HistType::kTH1F, "hQvector", "Q-vector; Q-vector; Entries"}, - {kEventPlaneAngle, o2::framework::HistType::kTH1F, "hEventPlaneAngle", "Event Plane angle; #Psi_{n}; Entries"}, + {kQvector, o2::framework::HistType::kTH1F, "hQvector", "q-vector; q-vector; Entries"}, + {kEventPlaneAngle, o2::framework::HistType::kTH1F, "hEventPlaneAngle", "Event Plane angle; #Psi_{EP}; Entries"}, // mc {kTruePosZ, o2::framework::HistType::kTH1F, "hTruePosZ", "True vertex Z (mc-truth collision); V_{Z,True} (cm); Entries"}, {kTrueCent, o2::framework::HistType::kTH1F, "hTrueCent", "True centrality (mc-truth collision); Centrality_{True} (%); Entries"}, diff --git a/PWGCF/Femto/Core/dataTypes.h b/PWGCF/Femto/Core/dataTypes.h index 0ed2fb12ec7..25c5d84147c 100644 --- a/PWGCF/Femto/Core/dataTypes.h +++ b/PWGCF/Femto/Core/dataTypes.h @@ -64,7 +64,7 @@ using CharmHadronMaskType = uint32_t; using CharmHadronType = uint16_t; // datatypes for event shape enums -using QvecDetectorType = uint8_t; +using EventShapeDetectorType = uint8_t; using QvecHarmonicType = uint8_t; // datatype for kinematic variable diff --git a/PWGCF/Femto/Core/modes.h b/PWGCF/Femto/Core/modes.h index bb665fe1057..f55830b958a 100644 --- a/PWGCF/Femto/Core/modes.h +++ b/PWGCF/Femto/Core/modes.h @@ -16,6 +16,8 @@ #ifndef PWGCF_FEMTO_CORE_MODES_H_ #define PWGCF_FEMTO_CORE_MODES_H_ +#include "dataTypes.h" + #include "PWGCF/Femto/Core/dataTypes.h" #include @@ -178,10 +180,10 @@ enum class CharmHadron : o2::analysis::femto::datatypes::CharmHadronType { kLcBar }; -enum class QvecDetector : o2::analysis::femto::datatypes::QvecDetectorType { +enum class EventShapeDetector : o2::analysis::femto::datatypes::EventShapeDetectorType { kFT0C = 0, kFT0A = 1, - kQvecDetectorLast = 2 + kEventShapeDetectorLast = 2 }; enum class QvecHarmonic : o2::analysis::femto::datatypes::QvecHarmonicType { diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index 08a40744013..5f2a285f3ec 100644 --- a/PWGCF/Femto/Core/pairHistManager.h +++ b/PWGCF/Femto/Core/pairHistManager.h @@ -363,7 +363,7 @@ constexpr std::array, kPairHistogramLast> {kTrueQoutVsQout, o2::framework::HistType::kTH2F, "hTrueQoutVsQout", "q_{out,True} vs q_{out}; q_{out,True} (GeV/#it{c}); q_{out} (GeV/#it{c})"}, {kTrueQsideVsQside, o2::framework::HistType::kTH2F, "hTrueQsideVsQside", "q_{side,True} vs q_{side}; q_{side,True} (GeV/#it{c}); q_{side} (GeV/#it{c})"}, {kTrueQlongVsQlong, o2::framework::HistType::kTH2F, "hTrueQlongVsQlong", "q_{long,True} vs q_{long}; q_{long,True} (GeV/#it{c}); q_{long} (GeV/#it{c})"}, - {kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, o2::framework::HistType::kTHnSparseF, "hQoutQsideQlongEventPlaneAngleQvector", "Event shape enginering; q_{out} (GeV/#it{c}); q_{side} (GeV/#it{c}); q_{long} (GeV/#it{c}); #varphi_{EP}; q-vector;"}, + {kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, o2::framework::HistType::kTHnSparseF, "hQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector", "Event shape enginering; q_{out} (GeV/#it{c}); q_{side} (GeV/#it{c}); q_{long} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); centrality (%); #varphi_{pair} - #Psi_{EP}; q-vector;"}, }}; // NOLINTNEXTLINE(cppcoreguidelines-macro-usage)