Committing older uncom changes

Author: Aswathy Menon

PWGDQ/Core/VarManager.cxx

@@ -816,6 +816,10 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kMultDimuonsME] = "";
   fgVariableNames[kCentFT0C] = "Centrality FT0C";
   fgVariableUnits[kCentFT0C] = "%";
+  fgVariableNames[kSpherocity] = "Spherocity"; // transverse spherocity
+  fgVariableUnits[kSpherocity] = "";
+  fgVariableNames[kSpherocityPtWeighted] = "SpherocityPtWeighted"; // pT-weighted transverse spherocity
+  fgVariableUnits[kSpherocityPtWeighted] = "";
   fgVariableNames[kMCEventGeneratorId] = "MC Generator ID";
   fgVariableNames[kMCEventSubGeneratorId] = "MC SubGenerator ID";
   fgVariableNames[kMCVtxX] = "MC Vtx X";
@@ -2096,6 +2100,8 @@ void VarManager::SetDefaultVarNames()
   fgVarNamesMap["kVtxChi2"] = kVtxChi2;
   fgVarNamesMap["kCentVZERO"] = kCentVZERO;
   fgVarNamesMap["kCentFT0C"] = kCentFT0C;
+  fgVarNamesMap["kSpherocity"] = kSpherocity; // transverse spherocity
+  fgVarNamesMap["kSpherocityPtWeighted"] = kSpherocityPtWeighted; // pT-weighted transverse spherocity
   fgVarNamesMap["kCentFT0A"] = kCentFT0A;
   fgVarNamesMap["kCentFT0M"] = kCentFT0M;
   fgVarNamesMap["kMultTPC"] = kMultTPC;

PWGDQ/Core/VarManager.h

@@ -118,6 +118,7 @@ class VarManager : public TObject
     RapidityGapFilter = BIT(21),
     Fit = BIT(22),
     ReducedFit = BIT(23),
+    ReducedEventSpherocity = BIT(24),  // Achu-Changes
     Track = BIT(0),
     TrackCov = BIT(1),
     TrackExtra = BIT(2),
@@ -283,6 +284,8 @@ class VarManager : public TObject
     kMultNTracksITSTPC,
     kMultNTracksPVeta1,
     kMultNTracksPVetaHalf,
+    kSpherocity,         // transverse spherocity
+    kSpherocityPtWeighted, // pT-weighted transverse spherocity
     kTrackOccupancyInTimeRange,
     kFT0COccupancyInTimeRange,
     kNoCollInTimeRangeStandard,
@@ -1438,6 +1441,12 @@ class VarManager : public TObject
   {
     return RecoDecay::constrainAngle(psi1 - psi2, -o2::constants::math::PI / harmonic, harmonic);
   }
+
+  // Spherocity function declaration--Start
+  template <typename T>
+  static float CalculateSpherocity(T const & tracks, float ptMin = 0.15, float ptMax = 10.0, float etaMin = -0.8, float etaMax = 0.8, int minMult = 10, bool usePtWeight = false);
+  //------------------------------------End
+
   template <typename T, typename T1>
   static o2::dataformats::VertexBase RecalculatePrimaryVertex(T const& track0, T const& track1, const T1& collision);
 
@@ -2316,6 +2325,11 @@ void VarManager::FillEvent(T const& event, float* values)
     }
   }
 
+  if constexpr ((fillMap & ObjTypes::ReducedEventSpherocity) > 0) {
+    values[kSpherocity] = event.spherocity();
+    values[kSpherocityPtWeighted] = event.spherocityPtWeighted();
+  }
+
   if constexpr ((fillMap & ReducedEventVtxCov) > 0) {
     values[kVtxCovXX] = event.covXX();
     values[kVtxCovXY] = event.covXY();
@@ -7996,5 +8010,83 @@ void VarManager::FillTripletVertexingALICE3(C const& collision, T const& t1, T c
     values[kVertexingTauxyzProjected] = values[kVertexingLxyzProjected] * v123.M() / (v123.P());
   }
 }
+template <typename T>
+float VarManager::CalculateSpherocity(T const& tracks, float ptMin, float ptMax, float etaMin, float etaMax, int minMult, bool usePtWeight) 
+{
+  LOGF(info, "    >>> CalculateSpherocity CALLED: ptMin=%.2f, ptMax=%.2f, etaMin=%.2f, etaMax=%.2f, minMult=%d, usePtWeight=%s",
+    ptMin, ptMax, etaMin, etaMax, minMult, usePtWeight ? "true" : "false");
+  
+  std::vector<float> pxNorm;
+  std::vector<float> pyNorm;
+  std::vector<float> ptValues;
+  float sumPt = 0.0;
+  int nTracks = 0;
+
+  // Collect valid tracks
+  for (const auto& track : tracks) {
+    float pt = track.pt();
+    float eta = track.eta();
+    
+    if (pt < ptMin || pt > ptMax) continue;
+    if (eta < etaMin || eta > etaMax) continue;
+    
+    float px = track.px();
+    float py = track.py();
+
+    if (pt <= 0.f) {
+      continue;
+    }
+    pxNorm.push_back(px / pt);
+    pyNorm.push_back(py / pt);
+    ptValues.push_back(pt);
+    sumPt += pt;
+    nTracks++;
+  }
+
+  // Check minimum multiplicity
+  if (nTracks < minMult) {
+    return -1.0f;
+  }
+
+  // Minimize over azimuthal angle
+  float spherocity = 5.0f; //float prevSpherocity = 5.0f; const float convergenceThreshold = 1e-4f;
+  const int nSteps = 360;
+
+  for (int iStep = 0; iStep < nSteps; iStep++) {
+    float phi = 2.0f * M_PI * static_cast<float>(iStep) / static_cast<float>(nSteps);
+    float nx = std::cos(phi);
+    float ny = std::sin(phi);
+
+    float numerator = 0.0f;
+    float denominator = 0.0f;
+
+    for (int iTrk = 0; iTrk < nTracks; iTrk++) {
+      // Cross product magnitude: |p_T x n| = |px*ny - py*nx|
+      float crossProduct = std::abs(pxNorm[iTrk] * ny - pyNorm[iTrk] * nx);
+      
+      if (usePtWeight) {
+        numerator += ptValues[iTrk] * crossProduct;
+        denominator += ptValues[iTrk];
+      } else {
+        numerator += crossProduct;
+        denominator += 1.0f;
+      }
+    }
+
+    float sph = std::pow(numerator / denominator, 2);
+    if (sph < spherocity) {
+      spherocity = sph;
+    }
+    // // Early exit if converged
+    // if (iStep > 10 && std::abs(spherocity - prevSpherocity) < convergenceThreshold) {
+    //   break;
+    // }
+    // prevSpherocity = spherocity;
+  }
+
+  spherocity *= (M_PI * M_PI / 4.0f);
+  return spherocity;
+}
+
 
 #endif // PWGDQ_CORE_VARMANAGER_H_

PWGDQ/DataModel/ReducedInfoTables.h

@@ -69,6 +69,8 @@ DECLARE_SOA_COLUMN(NTPCoccupMeanTimeShortA, nTPCoccupMeanTimeShortA, float);
 DECLARE_SOA_COLUMN(NTPCoccupMeanTimeShortC, nTPCoccupMeanTimeShortC, float);                                 //!  TPC pileup mean time on C side (short time range)
 DECLARE_SOA_COLUMN(NTPCoccupMedianTimeShortA, nTPCoccupMedianTimeShortA, float);                             //!  TPC pileup median time on A side (short time range)
 DECLARE_SOA_COLUMN(NTPCoccupMedianTimeShortC, nTPCoccupMedianTimeShortC, float);                             //!  TPC pileup median time on C side (short time range)
+DECLARE_SOA_COLUMN(Spherocity, spherocity, float);                                                           //! transverse spherocity
+DECLARE_SOA_COLUMN(SpherocityPtWeighted, spherocityPtWeighted, float);                                       //! pT-weighted transverse spherocity// namespace reducedevent
 DECLARE_SOA_COLUMN(DCAzBimodalityCoefficient, dcazBimodalityCoefficient, float);                             //!  Bimodality coefficient of the DCAz distribution of the tracks in the event
 DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinned, dcazBimodalityCoefficientBinned, float);                 //!  Bimodality coefficient of the DCAz distribution of the tracks in the event, binned
 DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinnedTrimmed1, dcazBimodalityCoefficientBinnedTrimmed1, float); //!  Bimodality coefficient of the DCAz distribution of the tracks in the event, binned and trimmed 1
@@ -177,7 +179,23 @@ DECLARE_SOA_TABLE_VERSIONED(ReducedEventsExtended_001, "AOD", "REEXTENDED", 1, /
                             mult::MultFDDA, mult::MultFDDC, mult::MultZNA, mult::MultZNC, mult::MultTracklets, mult::MultNTracksPV,
                             cent::CentFT0C, cent::CentFT0A, cent::CentFT0M);
 
+DECLARE_SOA_TABLE_VERSIONED(ReducedEventsExtended_002, "AOD", "REEXTENDED", 2, //!  Extended event information
+                            bc::GlobalBC, evsel::Alias, evsel::Selection, timestamp::Timestamp, cent::CentRun2V0M,
+                            mult::MultTPC, mult::MultFV0A, mult::MultFV0C, mult::MultFT0A, mult::MultFT0C,
+                            mult::MultFDDA, mult::MultFDDC, mult::MultZNA, mult::MultZNC, mult::MultTracklets, mult::MultNTracksPV,
+                            cent::CentFT0C, cent::CentFT0A, cent::CentFT0M,
+                            reducedevent::Spherocity,
+                            reducedevent::SpherocityPtWeighted);
+
 using ReducedEventsExtended = ReducedEventsExtended_001;
+// using ReducedEventsExtended = ReducedEventsExtended_002;
+
+// // In ReducedInfoTables.h, add this after the ReducedEventsExtended declaration
+// DECLARE_SOA_TABLE(ReducedEventsSpherocity, "AOD", "RESPHEROCITY",
+//                   reducedevent::Spherocity,
+//                   reducedevent::SpherocityPtWeighted);
+
+// using ReducedEventSphero = ReducedEventsSpherocity::iterator;
 
 DECLARE_SOA_TABLE(ReducedEventsMultPV_000, "AOD", "REMULTPV", //!  Multiplicity information for primary vertex
                   mult::MultNTracksHasITS, mult::MultNTracksHasTPC, mult::MultNTracksHasTOF, mult::MultNTracksHasTRD,

PWGDQ/TableProducer/tableMaker_withAssoc.cxx

@@ -222,6 +222,7 @@ struct TableMaker {
     // Q-vector related tables, to be filled only if the user selects the corresponding option
     Produces<ReducedEventsQvectorCentr> eventQvectorCentr;
     Produces<ReducedEventsQvectorCentrExtra> eventQvectorCentrExtra;
+    Produces<ReducedEventsSpherocity> eventSpherocity; // Achu-Changes
   } outTables;
 
   OutputObj<THashList> fOutputList{"output"}; //! the histogram manager output list
@@ -333,6 +334,16 @@ struct TableMaker {
     Configurable<bool> fExcludeShort{"cfgTPCExcludeShort", true, "Exclude short term from long term occupancy (micro-seconds)"};
   } fConfigVariousOptions;
 
+  // Spherocity configurables
+  struct : ConfigurableGroup {
+    Configurable<float> fConfigSpherocityPtMin{"cfgSpherocityPtMin", 0.15f, "Min pT for spherocity calculation"};
+    Configurable<float> fConfigSpherocityPtMax{"cfgSpherocityPtMax", 10.0f, "Max pT for spherocity calculation"};
+    Configurable<float> fConfigSpherocityEtaMin{"cfgSpherocityEtaMin", -0.8f, "Min eta for spherocity calculation"};
+    Configurable<float> fConfigSpherocityEtaMax{"cfgSpherocityEtaMax", 0.8f, "Max eta for spherocity calculation"};
+    Configurable<int> fConfigSpherocityMinMult{"cfgSpherocityMinMult", 5, "Min multiplicity for spherocity calculation"};
+    Configurable<bool> fConfigSpherocityUsePtWeight{"cfgSpherocityUsePtWeight", false, "Use pT weighting for spherocity"};
+  } fConfigSpherocity;
+
   Service<o2::ccdb::BasicCCDBManager> fCCDB{};
   o2::ccdb::CcdbApi fCCDBApi;
 
@@ -1033,15 +1044,22 @@ struct TableMaker {
     }
   }
 
-  template <uint32_t TEventFillMap, typename TEvents, typename TBCs,
-            typename TZdcs, typename TFt0s, typename TFv0as, typename TFdds>
+  template <uint32_t TEventFillMap, uint32_t TTrackFillMap, typename TEvents, typename TBCs,
+            typename TZdcs, typename TFt0s, typename TFv0as, typename TFdds, typename TTrackAssoc, typename TTracks>
   void skimCollisions(TEvents const& collisions, TBCs const& bcs, TZdcs const& /*zdcs*/,
-                      TFt0s const& ft0s, TFv0as const& fv0as, TFdds const& fdds)
+                      TFt0s const& ft0s, TFv0as const& fv0as, TFdds const& fdds,
+                      TTrackAssoc const& trackAssocs, TTracks const& tracksBarrel)
   {
     // Skim collisions
     // NOTE: So far, collisions are filtered based on the user specified analysis cuts AND the filterPP or Zorro event filter.
     //      The collision-track associations which point to an event that is not selected for writing are discarded!
 
+    LOGF(info, "=== skimCollisions CALLED ===");
+    LOGF(info, "TEventFillMap = %u", TEventFillMap);
+    LOGF(info, "TTrackFillMap = %u", TTrackFillMap);
+    LOGF(info, "TTrackFillMap is zero: %s", TTrackFillMap == 0 ? "TRUE" : "FALSE");
+    LOGF(info, "Number of collisions = %d", static_cast<int>(collisions.size()));
+
     VarManager::FillTimeFrame(collisions);
     fCollIndexMap.clear();
     int multTPC = -1.0;
@@ -1222,6 +1240,78 @@ struct TableMaker {
       }
       (dynamic_cast<TH2D*>(fStatsList->At(kStatsEvent)))->Fill(3.0, static_cast<float>(o2::aod::evsel::kNsel));
 
+      // Calculate spherocity for this collision
+      float spherocityValue = -1.0f;
+      float spherocityPtWeighted = -1.0f;
+
+      if constexpr (static_cast<bool>(TTrackFillMap)) {
+        // DIAGNOSTIC: Check table sizes
+        LOGF(info, "=== skimCollisions DIAGNOSTIC ===");
+        LOGF(info, "Total collisions: %d", static_cast<int>(collisions.size()));
+        LOGF(info, "Total tracksBarrel: %d", static_cast<int>(tracksBarrel.size()));
+        LOGF(info, "Total trackAssocs: %d", static_cast<int>(trackAssocs.size()));
+        LOGF(info, "================================");
+
+        LOGF(info, ">>> Collision %d: TTrackFillMap check PASSED", collision.globalIndex());
+
+        auto groupedTrackIndices = trackAssocs.sliceBy(trackIndicesPerCollision, collision.globalIndex());
+        LOGF(info, ">>> Collision %d: groupedTrackIndices.size() = %d",
+             collision.globalIndex(), static_cast<int>(groupedTrackIndices.size()));
+
+        if (static_cast<int>(groupedTrackIndices.size()) >= fConfigSpherocity.fConfigSpherocityMinMult) {
+          LOGF(info, ">>> Collision %d: PASSED multiplicity check, collecting tracks...",
+               collision.globalIndex());
+
+          std::vector<typename std::decay_t<decltype(tracksBarrel.begin())>> trackVec;
+          trackVec.reserve(groupedTrackIndices.size());
+
+          int validTracks = 0;
+          for (auto const& assoc : groupedTrackIndices) {
+            if (assoc.trackId() < tracksBarrel.size()) {
+              auto track = tracksBarrel.rawIteratorAt(assoc.trackId());
+              trackVec.push_back(track);
+              validTracks++;
+            }
+          }
+
+          LOGF(info, ">>> Collision %d: Collected %d valid tracks",
+               collision.globalIndex(), validTracks);
+
+          if (validTracks >= fConfigSpherocity.fConfigSpherocityMinMult) {
+            // Calculate unweighted spherocity
+            spherocityValue = VarManager::CalculateSpherocity(
+              trackVec,
+              fConfigSpherocity.fConfigSpherocityPtMin,
+              fConfigSpherocity.fConfigSpherocityPtMax,
+              fConfigSpherocity.fConfigSpherocityEtaMin,
+              fConfigSpherocity.fConfigSpherocityEtaMax,
+              fConfigSpherocity.fConfigSpherocityMinMult,
+              false);
+
+            // Calculate pT-weighted spherocity
+            spherocityPtWeighted = VarManager::CalculateSpherocity(
+              trackVec,
+              fConfigSpherocity.fConfigSpherocityPtMin,
+              fConfigSpherocity.fConfigSpherocityPtMax,
+              fConfigSpherocity.fConfigSpherocityEtaMin,
+              fConfigSpherocity.fConfigSpherocityEtaMax,
+              fConfigSpherocity.fConfigSpherocityMinMult,
+              true);
+
+            LOGF(info, "Event %d: Spherocity = %.3f, pT-weighted = %.3f, nTracks = %d",
+                 collision.globalIndex(), spherocityValue, spherocityPtWeighted, validTracks);
+          }
+        }
+      }
+
+      // CRITICAL: Fill VarManager::fgValues so histograms can be filled!
+      VarManager::fgValues[VarManager::kSpherocity] = spherocityValue;
+      VarManager::fgValues[VarManager::kSpherocityPtWeighted] = spherocityPtWeighted;
+      LOGF(info, "Filled VarManager: Spherocity = %.3f, pT-weighted = %.3f",
+           VarManager::fgValues[VarManager::kSpherocity],
+           VarManager::fgValues[VarManager::kSpherocityPtWeighted]);
+
+      // fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues);
       fHistMan->FillHistClass("Event_AfterCuts", dqtablemaker_helpers::varValues());
 
       // create the event tables
@@ -1245,9 +1335,10 @@ struct TableMaker {
         centFT0M = collision.centFT0M();
       }
       outTables.eventExtended(bc.globalBC(), collision.alias_raw(), collision.selection_raw(), bc.timestamp(), VarManager::fgValues[VarManager::kCentVZERO],
-                              multTPC, multFV0A, multFV0C, multFT0A, multFT0C, multFDDA, multFDDC, multZNA, multZNC, multTracklets, multTracksPV, centFT0C, centFT0A, centFT0M);
+                              multTPC, multFV0A, multFV0C, multFT0A, multFT0C, multFDDA, multFDDC, multZNA, multZNC, multTracklets, multTracksPV, centFT0C, centFT0A, centFT0M, spherocityValue, spherocityPtWeighted);
       outTables.eventVtxCov(collision.covXX(), collision.covXY(), collision.covXZ(), collision.covYY(), collision.covYZ(), collision.covZZ(), collision.chi2());
       outTables.eventInfo(collision.globalIndex());
+      outTables.eventSpherocity(spherocityValue, spherocityPtWeighted);
 
       if constexpr ((TEventFillMap & VarManager::ObjTypes::CollisionQvectCentr) > 0) {
         if (fConfigQvector.fConfigQvectCalibAvailable) {
@@ -1914,8 +2005,9 @@ struct TableMaker {
     outTables.event.reserve(collisions.size());
     outTables.eventExtended.reserve(collisions.size());
     outTables.eventVtxCov.reserve(collisions.size());
+    outTables.eventSpherocity.reserve(collisions.size()); // ADD THIS LINE//Achu-Changes
 
-    skimCollisions<TEventFillMap>(collisions, bcs, zdcs, ft0s, fv0as, fdds);
+    skimCollisions<TEventFillMap, TTrackFillMap>(collisions, bcs, zdcs, ft0s, fv0as, fdds, trackAssocs, tracksBarrel);
     if (fCollIndexMap.empty()) {
       return;
     }