[PWGJE,EMCAL-689] Add DeltaEta and DeltaPhi values to `EMCALMatchedTr…

PWGJE/Core/CMakeLists.txt

@@ -25,5 +25,6 @@ o2physics_target_root_dictionary(PWGJECore
                       JetBkgSubUtils.h
                       JetDerivedDataUtilities.h
                       emcalCrossTalkEmulation.h
+                      utilsTrackMatchingEMC.h
               LINKDEF PWGJECoreLinkDef.h)
 endif()

PWGJE/Core/JetUtilities.h

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 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.
 //
@@ -52,7 +52,7 @@
  * @returns (cluster to track index map, track to cluster index map)
  */
 template <typename T>
 std::tuple<std::vector<std::vector<int>>, std::vector<std::vector<int>>> MatchClustersAndTracks(
   std::vector<T>& clusterPhi,
   std::vector<T>& clusterEta,
   std::vector<T>& trackPhi,
@@ -76,11 +76,6 @@
     throw std::invalid_argument("track collection eta and phi sizes don't match. Check the inputs.");
   }
 
-  // for (std::size_t iTrack = 0; iTrack < nTracks; iTrack++) {
-  //   if (trackEta[iTrack] == 0)
-  //     LOG(warning) << "Track eta is 0!";
-  // }
-
   // Build the KD-trees using vectors
   // We build two trees:
   // treeBase, which contains the base collection.
@@ -144,7 +139,7 @@
 template <typename T, typename U>
 float deltaR(T const& A, U const& B)
 {
   float dPhi = RecoDecay::constrainAngle(A.phi() - B.phi(), -M_PI);
   float dEta = A.eta() - B.eta();
 
   return std::sqrt(dEta * dEta + dPhi * dPhi);
@@ -153,7 +148,7 @@
 template <typename T, typename U, typename V, typename W>
 float deltaR(T const& eta1, U const& phi1, V const& eta2, W const& phi2)
 {
   float dPhi = RecoDecay::constrainAngle(phi1 - phi2, -M_PI);
   float dEta = eta1 - eta2;
 
   return std::sqrt(dEta * dEta + dPhi * dPhi);

PWGJE/Core/utilsTrackMatchingEMC.h

@@ -0,0 +1,119 @@
+// Copyright 2019-2020 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.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file utilsTrackMatchingEMC.h
+/// \brief EMCal track matching related utils
+/// \author Marvin Hemmer <marvin.hemmer@cern.ch>
+
+#ifndef PWGJE_CORE_UTILSTRACKMATCHINGEMC_H_
+#define PWGJE_CORE_UTILSTRACKMATCHINGEMC_H_
+
+#include <TKDTree.h>
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <span>
+#include <stdexcept>
+#include <vector>
+
+namespace tmemcutilities
+{
+
+struct MatchResult {
+  std::vector<std::vector<int>> matchIndexTrack;
+  std::vector<std::vector<float>> matchDeltaPhi;
+  std::vector<std::vector<float>> matchDeltaEta;
+};
+
+/**
+ * Match clusters and tracks.
+ *
+ * Match cluster with tracks, where maxNumberMatches are considered in dR=maxMatchingDistance.
+ * If no unique match was found for a jet, an index of -1 is stored.
+ * The same map is created for clusters matched to tracks e.g. for electron analyses.
+ *
+ * @param clusterPhi cluster collection phi.
+ * @param clusterEta cluster collection eta.
+ * @param trackPhi track collection phi.
+ * @param trackEta track collection eta.
+ * @param maxMatchingDistance Maximum matching distance.
+ * @param maxNumberMatches Maximum number of matches (e.g. 5 closest).
+ *
+ * @returns (cluster to track index map, track to cluster index map)
+ */
+MatchResult matchTracksToCluster(
+  std::span<float> clusterPhi,
+  std::span<float> clusterEta,
+  std::span<float> trackPhi,
+  std::span<float> trackEta,
+  double maxMatchingDistance,
+  int maxNumberMatches)
+{
+  const std::size_t nClusters = clusterEta.size();
+  const std::size_t nTracks = trackEta.size();
+  MatchResult result;
+
+  result.matchIndexTrack.resize(nClusters);
+  result.matchDeltaPhi.resize(nClusters);
+  result.matchDeltaEta.resize(nClusters);
+
+  if (nClusters == 0 || nTracks == 0) {
+    // There are no jets, so nothing to be done.
+    return result;
+  }
+  // Input sizes must match
+  if (clusterPhi.size() != clusterEta.size()) {
+    throw std::invalid_argument("cluster collection eta and phi sizes don't match. Check the inputs.");
+  }
+  if (trackPhi.size() != trackEta.size()) {
+    throw std::invalid_argument("track collection eta and phi sizes don't match. Check the inputs.");
+  }
+
+  // Build the KD-trees using vectors
+  // We build two trees:
+  // treeBase, which contains the base collection.
+  // treeTag, which contains the tag collection.
+  // The trees are built to match in two dimensions (eta, phi)
+  TKDTree<int, float> treeTrack(trackEta.size(), 2, 1);
+  treeTrack.SetData(0, trackEta.data());
+  treeTrack.SetData(1, trackPhi.data());
+  treeTrack.Build();
+
+  // Find the track closest to each cluster.
+  for (std::size_t iCluster = 0; iCluster < nClusters; iCluster++) {
+    float point[2] = {clusterEta[iCluster], clusterPhi[iCluster]};
+    int index[50];      // size 50 for safety
+    float distance[50]; // size 50 for safery
+    std::fill_n(index, 50, -1);
+    std::fill_n(distance, 50, std::numeric_limits<float>::max());
+    treeTrack.FindNearestNeighbors(point, maxNumberMatches, index, distance);
+
+    // allocate enough memory
+    result.matchIndexTrack[iCluster].reserve(maxNumberMatches);
+    result.matchDeltaPhi[iCluster].reserve(maxNumberMatches);
+    result.matchDeltaEta[iCluster].reserve(maxNumberMatches);
+
+    // test whether indices are matching:
+    for (int m = 0; m < maxNumberMatches; m++) {
+      if (index[m] >= 0 && distance[m] < maxMatchingDistance) {
+        result.matchIndexTrack[iCluster].push_back(index[m]);
+        result.matchDeltaPhi[iCluster].push_back(trackPhi[index[m]] - clusterPhi[iCluster]);
+        result.matchDeltaEta[iCluster].push_back(trackEta[index[m]] - clusterEta[iCluster]);
+      }
+    }
+  }
+  return result;
+}
+}; // namespace tmemcutilities
+
+#endif // PWGJE_CORE_UTILSTRACKMATCHINGEMC_H_

PWGJE/DataModel/EMCALClusters.h

@@ -95,8 +95,8 @@
 };
 
 DECLARE_SOA_INDEX_COLUMN(Collision, collision);                        //! collisionID used as index for matched clusters
 DECLARE_SOA_INDEX_COLUMN(BC, bc);                                      //! bunch crossing ID used as index for ambiguous clusters
 DECLARE_SOA_COLUMN(ID, id, int);                                       //! cluster ID identifying cluster in event
 DECLARE_SOA_COLUMN(Energy, energy, float);                             //! cluster energy (GeV)
 DECLARE_SOA_COLUMN(CoreEnergy, coreEnergy, float);                     //! cluster core energy (GeV)
 DECLARE_SOA_COLUMN(RawEnergy, rawEnergy, float);                       //! raw cluster energy (GeV)
@@ -108,7 +108,7 @@
 DECLARE_SOA_COLUMN(Time, time, float);                                 //! cluster time (ns)
 DECLARE_SOA_COLUMN(IsExotic, isExotic, bool);                          //! flag to mark cluster as exotic
 DECLARE_SOA_COLUMN(DistanceToBadChannel, distanceToBadChannel, float); //! distance to bad channel
 DECLARE_SOA_COLUMN(NLM, nlm, int);                                     //! number of local maxima
 DECLARE_SOA_COLUMN(Definition, definition, int);                       //! cluster definition, see EMCALClusterDefinition.h
 
 } // namespace emcalcluster
@@ -148,11 +148,11 @@
 namespace emcalclustercell
 {
 // declare index column pointing to cluster table
 DECLARE_SOA_INDEX_COLUMN(EMCALCluster, emcalcluster); //! linked to EMCalClusters table
 DECLARE_SOA_INDEX_COLUMN(Calo, calo);                 //! linked to calo cells
 
 // declare index column pointing to ambiguous cluster table
 DECLARE_SOA_INDEX_COLUMN(EMCALAmbiguousCluster, emcalambiguouscluster); //! linked to EMCalAmbiguousClusters table
 } // namespace emcalclustercell
 DECLARE_SOA_TABLE(EMCALClusterCells, "AOD", "EMCCLUSCELLS",                                               //!
                   o2::soa::Index<>, emcalclustercell::EMCALClusterId, emcalclustercell::CaloId);          //!
@@ -162,10 +162,13 @@
 using EMCALAmbiguousClusterCell = EMCALAmbiguousClusterCells::iterator;
 namespace emcalmatchedtrack
 {
-DECLARE_SOA_INDEX_COLUMN(Track, track); //! linked to Track table only for tracks that were matched
+DECLARE_SOA_INDEX_COLUMN(Track, track);        //! linked to Track table only for tracks that were matched
+DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); //! difference between matched track and cluster azimuthal angle
+DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); //! difference between matched track and cluster pseudorapidity
 } // namespace emcalmatchedtrack
-DECLARE_SOA_TABLE(EMCALMatchedTracks, "AOD", "EMCMATCHTRACKS",                                     //!
-                  o2::soa::Index<>, emcalclustercell::EMCALClusterId, emcalmatchedtrack::TrackId); //!
+DECLARE_SOA_TABLE(EMCALMatchedTracks, "AOD", "EMCMATCHTRACKS", //!
+                  o2::soa::Index<>, emcalclustercell::EMCALClusterId, emcalmatchedtrack::TrackId,
+                  emcalmatchedtrack::DeltaPhi, emcalmatchedtrack::DeltaEta); //!
 using EMCALMatchedTrack = EMCALMatchedTracks::iterator;
 } // namespace o2::aod
 #endif // PWGJE_DATAMODEL_EMCALCLUSTERS_H_