[PWGCF] Add Track-TwoPhotonResonance pair task to femto framework

PWGCF/Femto/Core/closePairRejection.h

@@ -95,6 +95,7 @@ constexpr const char PrefixCprTrackLcPion[] = "CprTrackLcPion";
 constexpr const char PrefixCprD0DaughterD0DaughterPos[] = "CprD0DaughterD0DaughterPos";
 constexpr const char PrefixCprD0DaughterD0DaughterNeg[] = "CprD0DaughterD0DaughterNeg";
 constexpr const char PrefixCprTrackResonanceDaughter[] = "CprTrackResonanceDaughter";
+constexpr const char PrefixCprTrackTwoPhotonResonanceDaughter[] = "CprTrackTwoPhotonResonanceDaughter";
 constexpr const char PrefixCprTrackKinkDaughter[] = "CprTrackKinkDaughter";
 constexpr const char PrefixCprV0DaughterV0DaughterPos[] = "CprV0DaughterV0DaughterPos";
 constexpr const char PrefixCprV0DaughterV0DaughterNeg[] = "CprV0DaughterV0DaughterNeg";
@@ -113,6 +114,8 @@ using ConfCprTrackLcPion = ConfCpr<PrefixCprTrackLcPion>;
 using ConfCprD0DaugherD0DaughterPos = ConfCpr<PrefixCprD0DaughterD0DaughterPos>;
 using ConfCprD0DaugherD0DaughterNeg = ConfCpr<PrefixCprD0DaughterD0DaughterNeg>;
 using ConfCprTrackResonanceDaughter = ConfCpr<PrefixCprTrackResonanceDaughter>;
+// no cut by default: CPR is a no-op for two-photon resonances (see ClosePairRejectionTrackTwoPhotonResonance)
+using ConfCprTrackTwoPhotonResonanceDaughter = ConfCpr<PrefixCprTrackTwoPhotonResonanceDaughter, false>;
 using ConfCprTrackKinkDaughter = ConfCpr<PrefixCprTrackKinkDaughter>;
 using ConfCprV0DaugherV0DaughterPos = ConfCpr<PrefixCprV0DaughterV0DaughterPos>;
 using ConfCprV0DaugherV0DaughterNeg = ConfCpr<PrefixCprV0DaughterV0DaughterNeg>;
@@ -152,6 +155,8 @@ constexpr char PrefixV0TwoTrackResonancePosMe[] = "CPR_V0Resonance_PosDau/ME/";
 constexpr char PrefixV0TwoTrackResonanceNegMe[] = "CPR_V0Resonance_NegDau/ME/";
 constexpr char PrefixTrackTwoTrackResonanceSe[] = "CPR_TrackResonanceDau/SE/";
 constexpr char PrefixTrackTwoTrackResonanceMe[] = "CPR_TrackResonanceDau/ME/";
+constexpr char PrefixTrackTwoPhotonResonanceSe[] = "CPR_TrackTwoPhotonResonanceDau/SE/";
+constexpr char PrefixTrackTwoPhotonResonanceMe[] = "CPR_TrackTwoPhotonResonanceDau/ME/";
 constexpr char PrefixTrackCascadeBachelorSe[] = "CPR_TrackCascadeBachelor/SE/";
 constexpr char PrefixTrackCascadeBachelorMe[] = "CPR_TrackCascadeBachelor/ME/";
 constexpr char PrefixTrackKinkSe[] = "CPR_TrackKink/SE/";
@@ -320,14 +325,22 @@ class CloseTrackRejection
       swapTracks = (mSwapDist(mRng) == 1);
     }
 
-    auto t1 = swapTracks ? track2 : track1;
-    auto t2 = swapTracks ? track1 : track2;
+    // get kinematics of the particles
+    // if configured, randomly swap particle 1 and particle 2
+    float eta1 = swapTracks ? track2.eta() : track1.eta();
+    float eta2 = swapTracks ? track1.eta() : track2.eta();
+    float phi1 = swapTracks ? track2.phi() : track1.phi();
+    float phi2 = swapTracks ? track1.phi() : track2.phi();
+    float signedPt1 = swapTracks ? track2.signedPt() : track1.signedPt();
+    float signedPt2 = swapTracks ? track1.signedPt() : track2.signedPt();
+    int chargeAbs1 = swapTracks ? mChargeAbsTrack2 : mChargeAbsTrack1;
+    int chargeAbs2 = swapTracks ? mChargeAbsTrack1 : mChargeAbsTrack2;
 
-    mDeta = t1.eta() - t2.eta();
+    mDeta = eta1 - eta2;
 
     for (size_t i = 0; i < TpcRadii.size(); i++) {
-      auto phistar1 = phistar(mMagField, TpcRadii.at(i), mChargeAbsTrack1 * t1.signedPt(), t1.phi());
-      auto phistar2 = phistar(mMagField, TpcRadii.at(i), mChargeAbsTrack2 * t2.signedPt(), t2.phi());
+      auto phistar1 = phistar(mMagField, TpcRadii.at(i), chargeAbs1 * signedPt1, phi1);
+      auto phistar2 = phistar(mMagField, TpcRadii.at(i), chargeAbs2 * signedPt2, phi2);
       if (phistar1 && phistar2) {
         mDphistar.at(i) = RecoDecay::constrainAngle(phistar1.value() - phistar2.value(), -o2::constants::math::PI); // constrain angular difference between -pi and pi
         mDphistarMask.at(i) = true;
@@ -338,14 +351,14 @@ class CloseTrackRejection
     if (count > 0) {
       mAverageDphistar = std::accumulate(mDphistar.begin(), mDphistar.end(), 0.f) / count; // only average values if phistar could be computed
     } else {
-      mAverageDphistar = 0.f; // if computation at all radii fail, set it 0
+      mAverageDphistar = 0.f; // if computation at all radii failed, set it 0
     }
 
     if (mPlotAngularCorrelation) {
-      mPhi1 = t1.phi();
-      mPhi2 = t2.phi();
-      mEta1 = t1.eta();
-      mEta2 = t2.eta();
+      mPhi1 = phi1;
+      mPhi2 = phi2;
+      mEta1 = eta1;
+      mEta2 = eta2;
     }
   }
 
@@ -731,6 +744,32 @@ class ClosePairRejectionTrackKink
   CloseTrackRejection<prefix> mCtr;
 };
 
+// PCM photon daughters are, for now, not indexed into the track table (see photonBuilder.h):
+// they are TPC-only tracks, so they cannot be confused for a global track the way two V0
+// daughters or a track and a V0 daughter can -- there is no detector-overlap effect to correct
+// for. This is a no-op placeholder so Track-TwoPhotonResonance pairs can use the same generic
+// pair-building code path.
+template <auto& prefix>
+class ClosePairRejectionTrackTwoPhotonResonance
+{
+ public:
+  template <typename T>
+  void init(o2::framework::HistogramRegistry* /*registry*/,
+            std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& /*specs*/,
+            T const& /*confCpr*/,
+            int /*absChargeTrack*/)
+  {
+  }
+
+  void setMagField(float /*magField*/) {}
+
+  template <typename T1, typename T2, typename T3, typename T4>
+  [[nodiscard]] bool isClosePair(T1 const& /*track*/, T2 const& /*resonance*/, T3 const& /*photonTable*/, T4 const& /*pairHistManager*/)
+  {
+    return false;
+  }
+};
+
 template <auto& prefix>
 class ClosePairRejectionMcParticleMcParticle
 {

PWGCF/Femto/Core/femtoUtils.h

@@ -141,6 +141,12 @@ inline double getPdgMass(int pdgCode)
     case kOmegaMinus:
       mass = o2::constants::physics::MassOmegaMinus;
       break;
+    case kPi0:
+      mass = o2::constants::physics::MassPi0;
+      break;
+    case o2::constants::physics::kEta:
+      mass = o2::constants::physics::MassEta;
+      break;
     default:
       LOG(warn) << "PDG code " << pdgCode << " is not suppored. Return 0... ";
   }

PWGCF/Femto/Core/modes.h

@@ -113,6 +113,7 @@ enum class Particle : o2::analysis::femto::datatypes::ParticleType {
   kKink = 4,
   kCascade = 5,
   kCharmHadron = 6,
+  kTwoPhotonResonance = 7,
 };
 
 enum class McOrigin : o2::analysis::femto::datatypes::McOriginType {

PWGCF/Femto/Core/pairBuilder.h

@@ -29,6 +29,7 @@
 #include "PWGCF/Femto/Core/pairProcessHelpers.h"
 #include "PWGCF/Femto/Core/particleCleaner.h"
 #include "PWGCF/Femto/Core/trackHistManager.h"
+#include "PWGCF/Femto/Core/twoPhotonResonanceHistManager.h"
 #include "PWGCF/Femto/Core/twoTrackResonanceHistManager.h"
 #include "PWGCF/Femto/Core/v0HistManager.h"
 #include "PWGCF/Femto/DataModel/FemtoTables.h"
@@ -1384,6 +1385,129 @@ class PairTrackV0Builder
   int mMixingDepth = 5;
 };
 
+// Track-TwoPhotonResonance pair builder.
+// Note: close pair rejection is a no-op (ClosePairRejectionTrackTwoPhotonResonance) and pair cleaning
+// reuses the plain TrackTrackPairCleaner (globalIndex-only), since photon daughters are, for
+// now, not indexed into FTracks: they are TPC-only tracks, so they cannot be confused for a
+// global track the way two V0 daughters or a track and a V0 daughter can.
+template <auto& prefixTrack,
+          auto& prefixResonance,
+          auto& prefixDau1,
+          auto& prefixDau2,
+          auto& prefixSe,
+          auto& prefixMe,
+          auto& prefixCprSe,
+          auto& prefixCprMe,
+          modes::TwoPhotonResonance resonanceType>
+class PairTrackTwoPhotonResonanceBuilder
+{
+ public:
+  PairTrackTwoPhotonResonanceBuilder() = default;
+  ~PairTrackTwoPhotonResonanceBuilder() = default;
+
+  template <modes::Mode modeSe,
+            modes::Mode modeMe,
+            typename T1,
+            typename T2,
+            typename T3,
+            typename T4,
+            typename T5,
+            typename T6,
+            typename T7,
+            typename T8,
+            typename T9,
+            typename T10,
+            typename T11,
+            typename T12,
+            typename T13,
+            typename T14,
+            typename T15>
+  void init(o2::framework::HistogramRegistry* registry,
+            T1 const& confCollisionBinning,
+            T2 const& confTrackSelection,
+            T3 const& confResonanceSelection,
+            T4 const& confCpr,
+            T5 const& confMixing,
+            T6 const& confPairBinning,
+            T7 const& confPairCuts,
+            T8 const& confDauQaBinning,
+            std::map<T9, std::vector<o2::framework::AxisSpec>> const& colHistSpec,
+            std::map<T10, std::vector<o2::framework::AxisSpec>> const& trackHistSpec,
+            std::map<T11, std::vector<o2::framework::AxisSpec>> const& resonanceHistSpec,
+            std::map<T12, std::vector<o2::framework::AxisSpec>> const& dauHistSpec,
+            std::map<T13, std::vector<o2::framework::AxisSpec>> const& pairHistSpec,
+            std::map<T14, std::vector<o2::framework::AxisSpec>> const& cprHistSpec,
+            std::map<T15, std::vector<o2::framework::AxisSpec>> const& pairCleanerHistSpec)
+  {
+    mColHistManager.template init<modeSe>(registry, colHistSpec, confCollisionBinning);
+
+    mTrackHistManager.template init<modeSe>(registry, trackHistSpec, confTrackSelection);
+    mResonanceHistManager.template init<modeSe>(registry, resonanceHistSpec, dauHistSpec, confDauQaBinning);
+
+    mPairHistManagerSe.template init<modeSe>(registry, pairHistSpec, confPairBinning, confPairCuts, confMixing);
+    mPairHistManagerSe.setMass(confTrackSelection.pdgCodeAbs.value, confResonanceSelection.pdgCodeAbs.value);
+    mPairHistManagerSe.setCharge(confTrackSelection.chargeAbs.value, 1);
+    mCprSe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value);
+
+    mPairHistManagerMe.template init<modeMe>(registry, pairHistSpec, confPairBinning, confPairCuts, confMixing);
+    mPairHistManagerMe.setMass(confTrackSelection.pdgCodeAbs.value, confResonanceSelection.pdgCodeAbs.value);
+    mPairHistManagerMe.setCharge(confTrackSelection.chargeAbs.value, 1);
+    mCprMe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value);
+
+    mPcSe.template init<modeSe>(registry, pairCleanerHistSpec, confPairCuts);
+    mPcMe.template init<modeMe>(registry, pairCleanerHistSpec, confPairCuts);
+
+    // setup mixing
+    mMixingPolicy = static_cast<pairhistmanager::MixingPolicy>(confMixing.policy.value);
+    mMixingDepth = confMixing.depth.value;
+  }
+
+  /// \return true if at least one pair passed all pair selections (usable as pair trigger)
+  template <modes::Mode mode, typename T1, typename T2, typename T3, typename T4, typename T5>
+  bool processSameEvent(T1 const& col, T2& trackPartition, T3& resonancePartition, T4& photonTable, T5& cache)
+  {
+    auto trackSlice = trackPartition->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
+    auto resonanceSlice = resonancePartition->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
+    if (trackSlice.size() < nLimitPartitionParticles || resonanceSlice.size() < nLimitPartitionParticles) {
+      return false;
+    }
+    mColHistManager.template fill<mode>(col);
+    mCprSe.setMagField(col.magField());
+    return pairprocesshelpers::processSameEvent<mode>(trackSlice, resonanceSlice, photonTable, col, mTrackHistManager, mResonanceHistManager, mPairHistManagerSe, mCprSe, mPcSe);
+  }
+
+  template <modes::Mode mode, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8>
+  void processMixedEvent(T1 const& cols, T2& trackPartition, T3& resonancePartition, T4& photonTable, T5& cache, T6& binsVtxMult, T7& binsVtxCent, T8& binsVtxMultCent)
+  {
+    switch (mMixingPolicy) {
+      case static_cast<int>(pairhistmanager::kVtxMult):
+        pairprocesshelpers::processMixedEvent<mode>(cols, trackPartition, resonancePartition, photonTable, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPcMe);
+        break;
+      case static_cast<int>(pairhistmanager::kVtxCent):
+        pairprocesshelpers::processMixedEvent<mode>(cols, trackPartition, resonancePartition, photonTable, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPcMe);
+        break;
+      case static_cast<int>(pairhistmanager::kVtxMultCent):
+        pairprocesshelpers::processMixedEvent<mode>(cols, trackPartition, resonancePartition, photonTable, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPcMe);
+        break;
+      default:
+        LOG(fatal) << "Invalid binning policiy specifed. Breaking...";
+    }
+  }
+
+ private:
+  colhistmanager::CollisionHistManager mColHistManager;
+  trackhistmanager::TrackHistManager<prefixTrack> mTrackHistManager;
+  twophotonresonancehistmanager::TwoPhotonResonanceHistManager<prefixResonance, prefixDau1, prefixDau2, resonanceType> mResonanceHistManager;
+  pairhistmanager::PairHistManager<prefixSe, modes::Particle::kTrack, modes::Particle::kTwoPhotonResonance> mPairHistManagerSe;
+  pairhistmanager::PairHistManager<prefixMe, modes::Particle::kTrack, modes::Particle::kTwoPhotonResonance> mPairHistManagerMe;
+  closepairrejection::ClosePairRejectionTrackTwoPhotonResonance<prefixCprSe> mCprSe;
+  closepairrejection::ClosePairRejectionTrackTwoPhotonResonance<prefixCprMe> mCprMe;
+  paircleaner::TrackTrackPairCleaner<paircleaner::PrefixPairCleanerTrackTwoPhotonResonanceSe> mPcSe;
+  paircleaner::TrackTrackPairCleaner<paircleaner::PrefixPairCleanerTrackTwoPhotonResonanceMe> mPcMe;
+  pairhistmanager::MixingPolicy mMixingPolicy = pairhistmanager::MixingPolicy::kVtxMult;
+  int mMixingDepth = 5;
+};
+
 template <auto& prefixTrack,
           auto& prefixResonance,
           auto& prefixPosDau,

PWGCF/Femto/Core/pairCleaner.h

@@ -84,6 +84,8 @@ constexpr char PrefixPairCleanerTrackLcSe[] = "TrackLcCleaner/SE/";
 constexpr char PrefixPairCleanerTrackLcMe[] = "TrackLcCleaner/ME/";
 constexpr char PrefixPairCleanerTrackResonanceSe[] = "TrackResonanceCleaner/SE/";
 constexpr char PrefixPairCleanerTrackResonanceMe[] = "TrackResonanceCleaner/ME/";
+constexpr char PrefixPairCleanerTrackTwoPhotonResonanceSe[] = "TrackTwoPhotonResonanceCleaner/SE/";
+constexpr char PrefixPairCleanerTrackTwoPhotonResonanceMe[] = "TrackTwoPhotonResonanceCleaner/ME/";
 constexpr char PrefixPairCleanerV0ResonanceSe[] = "V0ResonanceCleaner/SE/";
 constexpr char PrefixPairCleanerV0ResonanceMe[] = "V0ResonanceCleaner/ME/";
 constexpr char PrefixPairCleanerTrackKinkSe[] = "TrackKinkCleaner/SE/";

PWGCF/Femto/Core/pairHistManager.h

@@ -570,6 +570,11 @@ constexpr char PrefixV0V0Me[] = "V0V0/ME/";
 constexpr char PrefixTrackResonanceSe[] = "TrackResonance/SE/";
 constexpr char PrefixTrackResonanceMe[] = "TrackResonance/ME/";
 
+constexpr char PrefixTrackPi0Se[] = "TrackPi0/SE/";
+constexpr char PrefixTrackPi0Me[] = "TrackPi0/ME/";
+constexpr char PrefixTrackEtaSe[] = "TrackEta/SE/";
+constexpr char PrefixTrackEtaMe[] = "TrackEta/ME/";
+
 constexpr char PrefixV0ResonanceSe[] = "V0Resonance/SE/";
 constexpr char PrefixV0ResonanceMe[] = "V0Resonance/ME/";
 

PWGCF/Femto/Core/photonBuilder.h

@@ -304,9 +304,11 @@ struct ConfPhotonTables : o2::framework::ConfigurableGroup {
 
 /// \brief Builder for PCM photons (V0PhotonsKF + V0Legs, from PWGEM/PhotonMeson)
 ///
-/// Unlike the strangeness V0 builder, daughter tracks are not indexed into FTracks: the
-/// conversion legs (V0Legs) carry a reduced set of columns and are not meant to be paired
-/// individually in femtoscopy, so only photon-level kinematics plus daughter PID/quality
+/// Unlike the strangeness V0 builder, daughter (conversion electron) tracks are, for now, not
+/// indexed into FTracks: the legs (V0Legs) are TPC-only tracks, so they cannot be confused for
+/// (i.e. share hits with/be the same reconstructed object as) a global track the way two V0
+/// daughters or a track and a V0 daughter can -- there is no detector-overlap effect for close
+/// pair rejection to correct for. Only photon-level kinematics plus daughter PID/quality
 /// diagnostics (FPhotonExtras) are stored.
 class PhotonBuilder
 {