[PWGCF] Add new process function in femto producer (#17388)

Author: ariedel-cern

PWGCF/Femto/Core/collisionBuilder.h

@@ -34,10 +34,6 @@
 #include <Framework/HistogramRegistry.h>
 #include <Framework/Logger.h>
 
-#include <sys/stat.h>
-
-#include <Rtypes.h>
-
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
@@ -86,7 +82,7 @@ struct ConfCollisionBits : o2::framework::ConfigurableGroup {
   o2::framework::Configurable<std::vector<float>> sphericityMax{"sphericityMax", {}, "Maximum sphericity"};
   o2::framework::Configurable<std::vector<std::string>> triggers{"triggers", {}, "List of all triggers to be used"};
   o2::framework::Configurable<datatypes::QvecDetectorType> qvecDetector{"qvecDetector", 0, "Detector used to estimate the Q-vector: 0 -> FT0C, 1 -> FT0A"};
-  o2::framework::Configurable<datatypes::QvecHarmonicType> qvecHarmonic{"qvecHarmonic", 2, "Harmonic n of the Q-vector and event plane angle Psi_n: 1 -> direct, 2 -> elliptic, 3 -> triangular"};
+  o2::framework::Configurable<datatypes::QvecHarmonicType> qvecHarmonic{"qvecHarmonic", 2, "Harmonic n of the Q-vector and event plane angle Psi_n: 2 -> elliptic, 3 -> triangular"};
 };
 
 struct ConfCcdb : o2::framework::ConfigurableGroup {
@@ -249,7 +245,13 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
 
     // event shape
     mQvecDetector = static_cast<modes::QvecDetector>(config.qvecDetector.value);
+    if (mQvecDetector >= modes::QvecDetector::kQvecDetectorLast) {
+      LOG(fatal) << "Qvector Detector is not supported";
+    }
     mQvecHarmonic = static_cast<modes::QvecHarmonic>(config.qvecHarmonic.value);
+    if (mQvecHarmonic < modes::QvecHarmonic::kN2 || mQvecHarmonic >= modes::QvecHarmonic::kQvecHarmonicLast) {
+      LOG(fatal) << "Qvector Harmonic is not supported";
+    }
 
     this->addSelection(kSel8, collisionSelectionNames.at(kSel8), config.sel8.value);
     this->addSelection(kNoSameBunchPileUp, collisionSelectionNames.at(kNoSameBunchPileUp), config.noSameBunchPileup.value);
@@ -366,20 +368,33 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
       case modes::QvecDetector::kFT0A:
         mQvec = std::hypot(col.qvecFT0AReVec()[0], col.qvecFT0AImVec()[0]) * std::sqrt(col.sumAmplFT0A());
         break;
+      case modes::QvecDetector::kQvecDetectorLast:
+        LOG(fatal) << "Invalid Q-vector detector";
+        break;
+      default:
+        LOG(fatal) << "Invalid Q-vector detector";
+        break;
     }
   }
   [[nodiscard]] float getQvector() const { return mQvec; }
 
   template <modes::System system, typename T>
   void setEventPlane(T const& col)
   {
-    float harmonic = static_cast<float>(mQvecHarmonic);
+    auto harmonic = static_cast<float>(mQvecHarmonic);
+    int index = static_cast<int>(mQvecHarmonic) - 2; // get index in the qvector vector
     switch (mQvecDetector) {
       case modes::QvecDetector::kFT0C:
-        mEventPlane = RecoDecay::constrainAngle((std::atan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])) / harmonic, 0, harmonic); // constrain between 0 and 2pi/harmonic
+        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()[0], col.qvecFT0AReVec()[0])) / harmonic, 0, harmonic); // constrain between 0 and 2pi/harmonic
+        mEventPlane = RecoDecay::constrainAngle((std::atan2(col.qvecFT0AImVec()[index], col.qvecFT0AReVec()[index])) / harmonic, 0, harmonic); // constrain between 0 and 2pi/harmonic
+        break;
+      case modes::QvecDetector::kQvecDetectorLast:
+        LOG(fatal) << "Invalid Q-vector detector";
+        break;
+      default:
+        LOG(fatal) << "Invalid Q-vector detector";
         break;
     }
   }

PWGCF/Femto/Core/modes.h

@@ -180,13 +180,14 @@ enum class CharmHadron : o2::analysis::femto::datatypes::CharmHadronType {
 
 enum class QvecDetector : o2::analysis::femto::datatypes::QvecDetectorType {
   kFT0C = 0,
-  kFT0A = 1
+  kFT0A = 1,
+  kQvecDetectorLast = 2
 };
 
 enum class QvecHarmonic : o2::analysis::femto::datatypes::QvecHarmonicType {
-  kN1 = 1,
   kN2 = 2,
-  kN3 = 3
+  kN3 = 3,
+  kQvecHarmonicLast = 4
 };
 
 }; // namespace o2::analysis::femto::modes

PWGCF/Femto/Core/v0HistManager.h

@@ -165,7 +165,7 @@ constexpr std::array<histmanager::HistInfo<V0Hist>, kV0HistLast> HistTable = {
    {kCosPa, o2::framework::HistType::kTH1F, "hCosPa", "Cosine of pointing angle; cos(#alpha); Entries"},
    {kDecayDauDca, o2::framework::HistType::kTH1F, "hDauDca", "Daughter DCA at decay vertex ; DCA_{Decay vertex} (cm); Entries"},
    {kStrangeTofPosDau, o2::framework::HistType::kTH1F, "hStrangeTofPosDau", "Strange TOF of positive Daughter ; n#sigma_{TOF, strange}; Entries"},
-   {kStrangeTofNegDau, o2::framework::HistType::kTH1F, "hStrangeTofNegDau", "Strange TOF of negative Daughter ; n#sigma+{TOF, strange}; Entries"},
+   {kStrangeTofNegDau, o2::framework::HistType::kTH1F, "hStrangeTofNegDau", "Strange TOF of negative Daughter ; n#sigma_{TOF, strange}; Entries"},
    {kDecayVtxX, o2::framework::HistType::kTH1F, "hDecayVtxX", "X coordinate of decay vertex ; DV_{X} (cm); Entries"},
    {kDecayVtxY, o2::framework::HistType::kTH1F, "hDecayVtxY", "Y coordinate of decay vertex ; DV_{Y} (cm); Entries"},
    {kDecayVtxZ, o2::framework::HistType::kTH1F, "hDecayVtxZ", "Z coordinate of decay vertex ; DV_{Z} (cm); Entries"},

PWGCF/Femto/TableProducer/femtoProducer.cxx

@@ -34,6 +34,7 @@
 #include "Common/DataModel/PIDResponseITS.h"
 #include "Common/DataModel/PIDResponseTOF.h"
 #include "Common/DataModel/PIDResponseTPC.h"
+#include "Common/DataModel/Qvectors.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 
 #include <CCDB/BasicCCDBManager.h>
@@ -77,7 +78,7 @@ using Run3FullPidTracks =
             o2::aod::pidTOFbeta, o2::aod::pidTOFmass>;
 using Run3McRecoTracks = soa::Join<Run3FullPidTracks, o2::aod::McTrackLabels>;
 
-using Run3Vzeros = o2::soa::Join<o2::aod::V0Datas, o2::aod::V0TOFPIDs, o2::aod::V0TOFNSigmas>;
+using Run3Vzeros = o2::soa::Join<o2::aod::V0Datas, o2::aod::V0TOFNSigmas>;
 using Run3RecoVzeros = o2::soa::Join<Run3Vzeros, o2::aod::McV0Labels>;
 
 using Run3D0s = soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfMlD0>;
@@ -204,13 +205,15 @@ struct FemtoProducer {
     // ---- guard: enabled tables vs. enabled process function ------------------
     if ((xiBuilder.fillAnyTable() || omegaBuilder.fillAnyTable()) &&
         (!doprocessTracksV0sCascadesRun3pp && !doprocessTracksV0sCascadesRun3PbPb &&
+         !doprocessTracksV0sCascadesRun3PbPbWithEventShape &&
          !doprocessTracksV0sCascadesKinksRun3pp && !doprocessTracksV0sCascadesRun3ppMc &&
          !doprocessTracksV0sCascadesRun3PbPbMc)) {
       LOG(fatal) << "At least one cascade table is enabled, but wrong process function is enabled. Breaking...";
     }
     if ((lambdaBuilder.fillAnyTable() || antilambdaBuilder.fillAnyTable() || k0shortBuilder.fillAnyTable()) &&
         (!doprocessTracksV0sCascadesRun3pp && !doprocessTracksV0sCascadesRun3PbPb &&
          !doprocessTracksV0sRun3pp && !doprocessTracksV0sCascadesKinksRun3pp &&
+         !doprocessTracksV0sCascadesRun3PbPbWithEventShape &&
          !doprocessTracksV0sRun3ppMc && !doprocessTracksV0sRun3PbPb && !doprocessTracksV0sRun3PbPbMc &&
          !doprocessTracksV0sCascadesRun3ppMc && !doprocessTracksV0sCascadesRun3PbPbMc &&
          !doprocessTracksV0sKinksRun3ppMc)) {
@@ -285,6 +288,7 @@ struct FemtoProducer {
       static_cast<int>(doprocessTracksV0sRun3PbPb) +
       static_cast<int>(doprocessTracksV0sCascadesRun3pp) +
       static_cast<int>(doprocessTracksV0sCascadesRun3PbPb) +
+      static_cast<int>(doprocessTracksV0sCascadesRun3PbPbWithEventShape) +
       static_cast<int>(doprocessTracksKinksRun3pp) +
       static_cast<int>(doprocessTracksV0sCascadesKinksRun3pp) +
       static_cast<int>(doprocessTracksD0sRun3pp) +
@@ -531,6 +535,23 @@ struct FemtoProducer {
   }
   PROCESS_SWITCH(FemtoProducer, processTracksV0sCascadesRun3PbPb, "Provide tracks, v0s and cascades in PbPb collisions", false);
 
+  void processTracksV0sCascadesRun3PbPbWithEventShape(rawinputs::Run3PbPbCollisionsWithEventShape::iterator const& col,
+                                                      o2::aod::BCsWithTimestamps const& bcs,
+                                                      rawinputs::Run3FullPidTracks const& tracks,
+                                                      rawinputs::Run3Vzeros const& v0s,
+                                                      rawinputs::Run3Cascades const& cascades)
+  {
+    if (!processCollisions<modes::System::kPbPb_Run3>(col, bcs, tracks)) {
+      return;
+    }
+    auto tracksWithItsPid = o2::soa::Attach<rawinputs::Run3FullPidTracks, o2::aod::pidits::ITSNSigmaEl, o2::aod::pidits::ITSNSigmaPi, o2::aod::pidits::ITSNSigmaKa,
+                                            o2::aod::pidits::ITSNSigmaPr, o2::aod::pidits::ITSNSigmaDe, o2::aod::pidits::ITSNSigmaTr, o2::aod::pidits::ITSNSigmaHe>(tracks);
+    processTracks<modes::System::kPbPb_Run3>(col, tracksWithItsPid);
+    processV0s<modes::System::kPbPb_Run3>(col, tracks, v0s);
+    processCascades<modes::System::kPbPb_Run3>(col, tracks, cascades);
+  }
+  PROCESS_SWITCH(FemtoProducer, processTracksV0sCascadesRun3PbPbWithEventShape, "Provide tracks, v0s and cascades in PbPb collisions with event shape information", false);
+
   void processTracksKinksRun3pp(rawinputs::Run3PpCollisions::iterator const& col,
                                 o2::aod::BCsWithTimestamps const& bcs,
                                 rawinputs::Run3FullPidTracks const& tracks,