[PWGCF] FemtoDream: added event sphericity calculation

PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h

@@ -1,4 +1,4 @@
 // Copyright 2019-2025 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.
 //
@@ -11,21 +11,21 @@
 
 /// \file FemtoDreamCollisionSelection.h
 /// \brief FemtoDreamCollisionSelection - event selection within the o2femtodream framework
 /// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de
 
 #ifndef PWGCF_FEMTODREAM_CORE_FEMTODREAMCOLLISIONSELECTION_H_
 #define PWGCF_FEMTODREAM_CORE_FEMTODREAMCOLLISIONSELECTION_H_
 
 #include <string>
 #include <iostream>
 #include "Common/CCDB/TriggerAliases.h"
 #include "Common/DataModel/EventSelection.h"
 #include "Framework/HistogramRegistry.h"
 #include "Framework/Logger.h"
 
 using namespace o2::framework;
 
 namespace o2::analysis::femtoDream
 {
 
 /// \class FemtoDreamCollisionSelection
@@ -41,7 +41,7 @@
   /// \param checkTrigger whether or not to check for the trigger alias
   /// \param trig Requested trigger alias
   /// \param checkOffline whether or not to check for offline selection criteria
-  void setCuts(float zvtxMax, bool checkTrigger, int trig, bool checkOffline, bool addCheckOffline, bool checkRun3)
+  void setCuts(float zvtxMax, bool checkTrigger, int trig, bool checkOffline, bool addCheckOffline, bool checkRun3, float minSphericity, float sphericityPtmin)
   {
     mCutsSet = true;
     mZvtxMax = zvtxMax;
@@ -50,6 +50,8 @@
     mCheckOffline = checkOffline;
     mAddCheckOffline = addCheckOffline;
     mCheckIsRun3 = checkRun3;
+    mMinSphericity = minSphericity;
+    mSphericityPtmin = sphericityPtmin;
   }
 
   /// Initializes histograms for the task
@@ -66,6 +68,7 @@
     mHistogramRegistry->add("Event/MultNTracksPV", "; MultNTracksPV; Entries", kTH1F, {{200, 0, 200}});
     mHistogramRegistry->add("Event/MultNTracklets", "; MultNTrackslets; Entries", kTH1F, {{300, 0, 300}});
     mHistogramRegistry->add("Event/MultTPC", "; MultTPC; Entries", kTH1F, {{600, 0, 600}});
+    mHistogramRegistry->add("Event/Sphericity", "; Sphericity; Entries", kTH1F, {{100, 0, 1}});
   }
 
   /// Print some debug information
@@ -76,6 +79,8 @@
     LOG(info) << "Check trigger: " << mCheckTrigger;
     LOG(info) << "Trigger: " << mTrigger;
     LOG(info) << " Check offline: " << mCheckOffline;
+    LOG(info) << "Min sphericity: " << mMinSphericity;
+    LOG(info) << "Min Pt (sphericity): " << mSphericityPtmin;
   }
 
   /// Check whether the collisions fulfills the specified selections
@@ -180,9 +185,65 @@
   /// \param tracks All tracks
   /// \return value of the sphericity of the event
   template <typename T1, typename T2>
-  float computeSphericity(T1 const& /*col*/, T2 const& /*tracks*/)
+  float computeSphericity(T1 const& col, T2 const& tracks)
   {
-    return 2.f;
+    double ptTot = 0.;
+    double s00 = 0.; // elements of the sphericity matrix taken form EPJC72:2124
+    double s01 = 0.;
+    // double s10 = 0.;
+    double s11 = 0.;
+
+    int numOfTracks = col.numContrib();
+    if (numOfTracks < 3)
+      return -9999.;
+
+    for (auto const& track : tracks) {
+      double pt = track.pt();
+      double eta = track.eta();
+      double px = track.px();
+      double py = track.py();
+      if (TMath::Abs(pt) < mSphericityPtmin || TMath::Abs(eta) > 0.8) {
+        continue;
+      }
+
+      ptTot += pt;
+
+      s00 += px * px / pt;
+      s01 += px * py / pt;
+      // s10 = s01;
+      s11 += py * py / pt;
+    }
+
+    // normalize to total Pt to obtain a linear form:
+    if (ptTot == 0.)
+      return -9999.;
+    s00 /= ptTot;
+    s11 /= ptTot;
+    s01 /= ptTot;
+
+    // Calculate the trace of the sphericity matrix:
+    double T = s00 + s11;
+    // Calculate the determinant of the sphericity matrix:
+    double D = s00 * s11 - s01 * s01; // S10 = S01
+
+    // Calculate the eigenvalues of the sphericity matrix:
+    double lambda1 = 0.5 * (T + std::sqrt(T * T - 4. * D));
+    double lambda2 = 0.5 * (T - std::sqrt(T * T - 4. * D));
+
+    if ((lambda1 + lambda2) == 0.)
+      return -9999.;
+
+    double spt = -1.;
+
+    if (lambda2 > lambda1) {
+      spt = 2. * lambda1 / (lambda1 + lambda2);
+    } else {
+      spt = 2. * lambda2 / (lambda1 + lambda2);
+    }
+
+    mHistogramRegistry->fill(HIST("Event/Sphericity"), spt);
+
+    return spt;
   }
 
  private:
@@ -194,6 +255,8 @@
   bool mCheckIsRun3 = false;                       ///< Check if running on Pilot Beam
   triggerAliases mTrigger = kINT7;                 ///< Trigger to check for
   float mZvtxMax = 999.f;                          ///< Maximal deviation from nominal z-vertex (cm)
+  float mMinSphericity = 0.f;
+  float mSphericityPtmin = 0.f;
 };
 } // namespace o2::analysis::femtoDream
 

PWGCF/FemtoDream/TableProducer/CMakeLists.txt

@@ -14,11 +14,6 @@ o2physics_add_dpl_workflow(femtodream-producer
           PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils
           COMPONENT_NAME Analysis)
 
-o2physics_add_dpl_workflow(femtodream-producer-withcascades
-          SOURCES femtoDreamProducerTaskWithCascades.cxx
-          PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils
-          COMPONENT_NAME Analysis)
-
 o2physics_add_dpl_workflow(femtodream-producer-reduced
           SOURCES femtoDreamProducerReducedTask.cxx
           PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore

PWGCF/FemtoDream/TableProducer/femtoDreamProducerReducedTask.cxx

@@ -79,6 +79,8 @@ struct femtoDreamProducerReducedTask {
   Configurable<int> ConfEvtTriggerSel{"ConfEvtTriggerSel", kINT7, "Evt sel: trigger"};
   Configurable<bool> ConfEvtOfflineCheck{"ConfEvtOfflineCheck", false, "Evt sel: check for offline selection"};
   Configurable<bool> ConfEvtAddOfflineCheck{"ConfEvtAddOfflineCheck", false, "Evt sel: additional checks for offline selection (not part of sel8 yet)"};
+  Configurable<float> ConfEvtMinSphericity{"ConfEvtMinSphericity", 0.0f, "Evt sel: Min. sphericity of event"};
+  Configurable<float> ConfEvtSphericityPtmin{"ConfEvtSphericityPtmin", 0.0f, "Evt sel: Min. Pt for sphericity calculation"};
 
   Configurable<bool> ConfTrkRejectNotPropagated{"ConfTrkRejectNotPropagated", false, "True: reject not propagated tracks"};
 
@@ -116,7 +118,7 @@ struct femtoDreamProducerReducedTask {
     int CutBits = 8 * sizeof(o2::aod::femtodreamparticle::cutContainerType);
     Registry.add("AnalysisQA/CutCounter", "; Bit; Counter", kTH1F, {{CutBits + 1, -0.5, CutBits + 0.5}});
 
-    colCuts.setCuts(ConfEvtZvtx.value, ConfEvtTriggerCheck.value, ConfEvtTriggerSel.value, ConfEvtOfflineCheck.value, ConfEvtAddOfflineCheck.value, ConfIsRun3.value);
+    colCuts.setCuts(ConfEvtZvtx.value, ConfEvtTriggerCheck.value, ConfEvtTriggerSel.value, ConfEvtOfflineCheck.value, ConfEvtAddOfflineCheck.value, ConfIsRun3.value, ConfEvtMinSphericity.value, ConfEvtSphericityPtmin.value);
     colCuts.init(&qaRegistry);
 
     trackCuts.setSelection(ConfTrkCharge, femtoDreamTrackSelection::kSign, femtoDreamSelection::kEqual);

PWGCF/FemtoDream/TableProducer/femtoDreamProducerTask.cxx

@@ -131,6 +131,8 @@ struct femtoDreamProducerTask {
   Configurable<bool> ConfEvtAddOfflineCheck{"ConfEvtAddOfflineCheck", false, "Evt sel: additional checks for offline selection (not part of sel8 yet)"};
   Configurable<bool> ConfIsActivateV0{"ConfIsActivateV0", true, "Activate filling of V0 into femtodream tables"};
   Configurable<bool> ConfIsActivateReso{"ConfIsActivateReso", false, "Activate filling of sl Resonances into femtodream tables"};
+  Configurable<float> ConfEvtMinSphericity{"ConfEvtMinSphericity", 0.0f, "Evt sel: Min. sphericity of event"};
+  Configurable<float> ConfEvtSphericityPtmin{"ConfEvtSphericityPtmin", 0.0f, "Evt sel: Min. Pt for sphericity calculation"};
 
   Configurable<bool> ConfTrkRejectNotPropagated{"ConfTrkRejectNotPropagated", false, "True: reject not propagated tracks"};
   // Configurable<bool> ConfRejectITSHitandTOFMissing{ "ConfRejectITSHitandTOFMissing", false, "True: reject if neither ITS hit nor TOF timing satisfied"};
@@ -270,7 +272,7 @@ struct femtoDreamProducerTask {
 
     rctChecker.init(rctCut.cfgEvtRCTFlagCheckerLabel, false, rctCut.cfgEvtRCTFlagCheckerLimitAcceptAsBad);
 
-    colCuts.setCuts(ConfEvtZvtx.value, ConfEvtTriggerCheck.value, ConfEvtTriggerSel.value, ConfEvtOfflineCheck.value, ConfEvtAddOfflineCheck.value, ConfIsRun3.value);
+    colCuts.setCuts(ConfEvtZvtx.value, ConfEvtTriggerCheck.value, ConfEvtTriggerSel.value, ConfEvtOfflineCheck.value, ConfEvtAddOfflineCheck.value, ConfIsRun3.value, ConfEvtMinSphericity.value, ConfEvtSphericityPtmin.value);
     colCuts.init(&qaRegistry);
 
     trackCuts.setSelection(ConfTrkCharge, femtoDreamTrackSelection::kSign, femtoDreamSelection::kEqual);