Adding MFT cluster sizes and tracks flags to AO2D (#12358)

* Including MFT cluster size in AO2D

* Fix clang

* Minor fix in a comment

* Fix clang

* Fixing layer limit in setClusterSize

---------

Co-authored-by: Lucamicheletti93 <luca.mike93@gmail.com>

Author: lucamicheletti93

DataFormats/Detectors/ITSMFT/MFT/include/DataFormatsMFT/TrackMFT.h

@@ -70,6 +70,24 @@ class TrackMFT : public o2::track::TrackParCovFwd
   const o2::track::TrackParCovFwd& getOutParam() const { return mOutParameters; }       ///< Returns track parameters fitted outwards
   void setOutParam(const o2::track::TrackParCovFwd parcov) { mOutParameters = parcov; } ///< Set track out parameters
 
+  void setClusterSize(int l, int size)
+  {
+    if (l >= 10) {
+      return;
+    }
+    if (size > 63) {
+      size = 63;
+    }
+
+    mClusterSizes &= ~(0x3fULL << (l * 6));
+    mClusterSizes |= (static_cast<uint64_t>(size) << (l * 6));
+  }
+
+  uint64_t getClusterSizes() const
+  {
+    return mClusterSizes;
+  }
+
  private:
   Bool_t mIsCA = false; ///< Track finding method CA vs. LTF
 
@@ -79,8 +97,9 @@ class TrackMFT : public o2::track::TrackParCovFwd
 
   Double_t mSeedinvQPtFitChi2 = 0.; ///< Seed InvQPt Chi2 from FCF clusters X,Y positions
   Double_t mInvQPtSeed;             ///< Seed InvQPt from FCF clusters X,Y positions
+  uint64_t mClusterSizes = 0;       ///< MFT cluster sizes per track
 
-  ClassDefNV(TrackMFT, 2);
+  ClassDefNV(TrackMFT, 3);
 };
 
 class TrackMFTExt : public TrackMFT
@@ -94,14 +113,28 @@ class TrackMFTExt : public TrackMFT
   using TrackMFT::TrackMFT; // inherit base constructors
 
   int getExternalClusterIndex(int i) const { return mExtClsIndex[i]; }
+  int getExternalClusterSize(int i) const { return mExtClsSize[i]; }
+  int getExternalClusterLayer(int i) const { return mExtClsLayer[i]; }
 
   void setExternalClusterIndex(int np, int idx)
   {
     mExtClsIndex[np] = idx;
   }
 
+  void setExternalClusterSize(int np, int clsSize)
+  {
+    mExtClsSize[np] = clsSize;
+  }
+
+  void setExternalClusterLayer(int np, int clsLayer)
+  {
+    mExtClsLayer[np] = clsLayer;
+  }
+
  protected:
   std::array<int, MaxClusters> mExtClsIndex = {-1}; ///< External indices of associated clusters
+  std::array<int, MaxClusters> mExtClsSize = {-1};  ///< Cluster size
+  std::array<int, MaxClusters> mExtClsLayer = {-1}; ///< Cluster layer
 
   ClassDefNV(TrackMFTExt, 1);
 };

Detectors/AOD/src/AODProducerWorkflowSpec.cxx

@@ -361,14 +361,18 @@ void AODProducerWorkflowDPL::addToMFTTracksTable(mftTracksCursorType& mftTracksC
   }
   trackTime -= bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS;
 
+  // the Cellular Automaton track-finding algorithm flag is stored in first of the 4 bits not used for the cluster size
+  uint64_t mftClusterSizesAndTrackFlags = track.getClusterSizes();
+  mftClusterSizesAndTrackFlags |= (track.isCA()) ? (1ULL << (60)) : 0;
+
   mftTracksCursor(collisionID,
                   track.getX(),
                   track.getY(),
                   truncateFloatFraction(track.getZ(), mTrackX), // for the forward tracks Z has the same role as X in barrel
                   truncateFloatFraction(track.getPhi(), mTrackAlpha),
                   truncateFloatFraction(track.getTanl(), mTrackTgl),
                   truncateFloatFraction(track.getInvQPt(), mTrack1Pt),
-                  track.getNumberOfPoints(),
+                  mftClusterSizesAndTrackFlags,
                   truncateFloatFraction(track.getTrackChi2(), mTrackChi2),
                   truncateFloatFraction(trackTime, mTrackTime),
                   truncateFloatFraction(trackTimeRes, mTrackTimeError));

Detectors/ITSMFT/MFT/tracking/include/MFTTracking/ROframe.h

@@ -47,6 +47,7 @@ class ROframe
     while (layer--) {
       mClusters[layer].reserve(nClusters * fraction);
       mClusterExternalIndices[layer].reserve(nClusters * fraction);
+      mClusterSizes[layer].reserve(nClusters * fraction);
     }
     mTracks.reserve(nClusters * fraction);
   }
@@ -58,6 +59,8 @@ class ROframe
 
   const Int_t getClusterExternalIndex(Int_t layerId, const Int_t clusterId) const { return mClusterExternalIndices[layerId][clusterId]; }
 
+  const Int_t getClusterSize(Int_t layerId, const Int_t clusterId) const { return mClusterSizes[layerId][clusterId]; }
+
   std::vector<T>& getTracks() { return mTracks; }
   T& getCurrentTrack() { return mTracks.back(); }
 
@@ -70,6 +73,7 @@ class ROframe
   }
   void addClusterLabelToLayer(Int_t layer, const MCCompLabel label) { mClusterLabels[layer].emplace_back(label); }
   void addClusterExternalIndexToLayer(Int_t layer, const Int_t idx) { mClusterExternalIndices[layer].push_back(idx); }
+  void addClusterSizeToLayer(Int_t layer, const Int_t clusterSize) { mClusterSizes[layer].push_back(clusterSize); }
 
   void addTrack(bool isCA = false)
   {
@@ -87,6 +91,7 @@ class ROframe
       mClusters[iLayer].clear();
       mClusterLabels[iLayer].clear();
       mClusterExternalIndices[iLayer].clear();
+      mClusterSizes[iLayer].clear();
     }
     mTracks.clear();
     mRoads.clear();
@@ -98,6 +103,7 @@ class ROframe
   std::array<std::vector<Cluster>, constants::mft::LayersNumber> mClusters;
   std::array<std::vector<MCCompLabel>, constants::mft::LayersNumber> mClusterLabels;
   std::array<std::vector<Int_t>, constants::mft::LayersNumber> mClusterExternalIndices;
+  std::array<std::vector<Int_t>, constants::mft::LayersNumber> mClusterSizes;
   std::vector<T> mTracks;
   std::vector<Road> mRoads;
 };

Detectors/ITSMFT/MFT/tracking/include/MFTTracking/TrackCA.h

@@ -45,7 +45,7 @@ class TrackLTF : public TrackMFTExt
   const std::array<Int_t, constants::mft::LayersNumber>& getLayers() const { return mLayer; }
   const std::array<Int_t, constants::mft::LayersNumber>& getClustersId() const { return mClusterId; }
   const std::array<MCCompLabel, constants::mft::LayersNumber>& getMCCompLabels() const { return mMCCompLabels; }
-  void setPoint(const Cluster& cl, const Int_t layer, const Int_t clusterId, const MCCompLabel label, const Int_t extClsIndex);
+  void setPoint(const Cluster& cl, const Int_t layer, const Int_t clusterId, const MCCompLabel label, const Int_t extClsIndex, const Int_t clsSize);
 
   void sort();
 
@@ -84,7 +84,7 @@ class TrackLTFL : public TrackLTF // A track model for B=0
 };
 
 //_________________________________________________________________________________________________
-inline void TrackLTF::setPoint(const Cluster& cl, const Int_t layer, const Int_t clusterId, const MCCompLabel label, const Int_t extClsIndex)
+inline void TrackLTF::setPoint(const Cluster& cl, const Int_t layer, const Int_t clusterId, const MCCompLabel label, const Int_t extClsIndex, const Int_t clsSize)
 {
   auto nPoints = getNumberOfPoints();
   if (nPoints > 0) {
@@ -106,6 +106,8 @@ inline void TrackLTF::setPoint(const Cluster& cl, const Int_t layer, const Int_t
   mClusterId[nPoints] = clusterId;
   mMCCompLabels[nPoints] = label;
   setExternalClusterIndex(nPoints, extClsIndex);
+  setExternalClusterSize(nPoints, clsSize);
+  setExternalClusterLayer(nPoints, layer);
   setNumberOfPoints(nPoints + 1);
 }
 

Detectors/ITSMFT/MFT/tracking/src/IOUtils.cxx

@@ -47,6 +47,7 @@ int ioutils::loadROFrameData(const o2::itsmft::ROFRecord& rof, ROframe<T>& event
   auto first = rof.getFirstEntry();
   auto clusters_in_frame = rof.getROFData(clusters);
   auto nClusters = clusters_in_frame.size();
+  int clusterSize = -999;
 
   bool skip_ROF = true;
   auto& trackingParam = MFTTrackingParam::Instance();
@@ -69,13 +70,16 @@ int ioutils::loadROFrameData(const o2::itsmft::ROFRecord& rof, ROframe<T>& event
       sigmaY2 = dict->getErr2Z(pattID); // ALPIDE local Z coordinate => MFT global Y coordinate (ALPIDE columns)
       if (!dict->isGroup(pattID)) {
         locXYZ = dict->getClusterCoordinates(c);
+        clusterSize = dict->getNpixels(pattID);
       } else {
         o2::itsmft::ClusterPattern patt(pattIt);
         locXYZ = dict->getClusterCoordinates(c, patt);
+        clusterSize = patt.getNPixels();
       }
     } else {
       o2::itsmft::ClusterPattern patt(pattIt);
       locXYZ = dict->getClusterCoordinates(c, patt, false);
+      clusterSize = patt.getNPixels();
     }
     if (skip_ROF) { // Skip filtered-out ROFs after processing pattIt
       clusterId++;
@@ -96,6 +100,7 @@ int ioutils::loadROFrameData(const o2::itsmft::ROFRecord& rof, ROframe<T>& event
       event.addClusterLabelToLayer(layer, *(mcLabels->getLabels(first + clusterId).begin()));
     }
     event.addClusterExternalIndexToLayer(layer, first + clusterId);
+    event.addClusterSizeToLayer(layer, clusterSize);
     clusterId++;
   }
   return nClusters;

Detectors/ITSMFT/MFT/tracking/src/Tracker.cxx

@@ -268,6 +268,7 @@ void Tracker<T>::findTracksLTF(ROframe<T>& event)
   Int_t binR_proj, binPhi_proj, bin;
   Int_t binIndex, clsMinIndex, clsMaxIndex, clsMinIndexS, clsMaxIndexS;
   Int_t extClsIndex;
+  Int_t clsSize;
   Float_t dz = 0., dRCone = 1.;
   Float_t dR2, dR2min, dR2cut = mLTFclsR2Cut;
   Bool_t hasDisk[constants::mft::DisksNumber], newPoint;
@@ -395,7 +396,8 @@ void Tracker<T>::findTracksLTF(ROframe<T>& event)
             Cluster& cluster = event.getClustersInLayer(layer)[clsInLayer];
             mcCompLabel = mUseMC ? event.getClusterLabels(layer, cluster.clusterId) : MCCompLabel();
             extClsIndex = event.getClusterExternalIndex(layer, cluster.clusterId);
-            event.getCurrentTrack().setPoint(cluster, layer, clsInLayer, mcCompLabel, extClsIndex);
+            clsSize = event.getClusterSize(layer, cluster.clusterId);
+            event.getCurrentTrack().setPoint(cluster, layer, clsInLayer, mcCompLabel, extClsIndex, clsSize);
             // mark the used clusters
             cluster.setUsed(true);
           }
@@ -418,6 +420,7 @@ void Tracker<T>::findTracksLTFfcs(ROframe<T>& event)
   Int_t binR_proj, binPhi_proj, bin;
   Int_t binIndex, clsMinIndex, clsMaxIndex, clsMinIndexS, clsMaxIndexS;
   Int_t extClsIndex;
+  Int_t clsSize;
   Float_t dR2, dR2min, dR2cut = mLTFclsR2Cut;
   Bool_t hasDisk[constants::mft::DisksNumber], newPoint;
 
@@ -530,7 +533,8 @@ void Tracker<T>::findTracksLTFfcs(ROframe<T>& event)
           Cluster& cluster = event.getClustersInLayer(layer)[clsInLayer];
           mcCompLabel = mUseMC ? event.getClusterLabels(layer, cluster.clusterId) : MCCompLabel();
           extClsIndex = event.getClusterExternalIndex(layer, cluster.clusterId);
-          event.getCurrentTrack().setPoint(cluster, layer, clsInLayer, mcCompLabel, extClsIndex);
+          clsSize = event.getClusterSize(layer, cluster.clusterId);
+          event.getCurrentTrack().setPoint(cluster, layer, clsInLayer, mcCompLabel, extClsIndex, clsSize);
           // mark the used clusters
           cluster.setUsed(true);
         }
@@ -1069,14 +1073,17 @@ void Tracker<T>::addCellToCurrentTrackCA(const Int_t layer1, const Int_t cellId,
   MCCompLabel mcCompLabel2 = mUseMC ? event.getClusterLabels(layer2, cluster2.clusterId) : MCCompLabel();
 
   Int_t extClsIndex;
+  Int_t clsSize;
 
   if (trackCA.getNumberOfPoints() == 0) {
     extClsIndex = event.getClusterExternalIndex(layer2, cluster2.clusterId);
-    trackCA.setPoint(cluster2, layer2, clsInLayer2, mcCompLabel2, extClsIndex);
+    clsSize = event.getClusterSize(layer2, cluster2.clusterId);
+    trackCA.setPoint(cluster2, layer2, clsInLayer2, mcCompLabel2, extClsIndex, clsSize);
   }
 
   extClsIndex = event.getClusterExternalIndex(layer1, cluster1.clusterId);
-  trackCA.setPoint(cluster1, layer1, clsInLayer1, mcCompLabel1, extClsIndex);
+  clsSize = event.getClusterSize(layer2, cluster2.clusterId);
+  trackCA.setPoint(cluster1, layer1, clsInLayer1, mcCompLabel1, extClsIndex, clsSize);
 }
 
 //_________________________________________________________________________________________________

Detectors/ITSMFT/MFT/workflow/src/TrackerSpec.cxx

@@ -197,6 +197,9 @@ void TrackerDPL::run(ProcessingContext& pc)
       int ncl = trc.getNumberOfPoints();
       for (int ic = 0; ic < ncl; ic++) {
         auto externalClusterID = trc.getExternalClusterIndex(ic);
+        auto clusterSize = trc.getExternalClusterSize(ic);
+        auto clusterLayer = trc.getExternalClusterLayer(ic);
+        trc.setClusterSize(clusterLayer, clusterSize);
         allClusIdx.push_back(externalClusterID);
       }
       allTracks.emplace_back(trc);

Framework/Core/include/Framework/AnalysisDataModel.h

@@ -506,7 +506,7 @@ DECLARE_SOA_COLUMN(Phi, phi, float);
 DECLARE_SOA_COLUMN(Tgl, tgl, float);                                                         //! TrackParFwd parameter tan(\lamba); (\lambda = 90 - \theta_{polar})
 DECLARE_SOA_COLUMN(Signed1Pt, signed1Pt, float);                                             //! TrackParFwd parameter: charged inverse transverse momentum; (q/pt)
 DECLARE_SOA_COLUMN(NClusters, nClusters, int8_t);                                            //! Number of clusters
-DECLARE_SOA_COLUMN(MFTClusterSizesAndTracksFlags, mftClusterSizesAndTrackFlags, uint64_t);   //! Cluster sizes per track, stored per layer (each 6 bits). Remaining 4 bits for MFT flags
+DECLARE_SOA_COLUMN(MFTClusterSizesAndTrackFlags, mftClusterSizesAndTrackFlags, uint64_t);    //! Cluster sizes per track, stored per layer (each 6 bits). Remaining 4 bits for MFT flags
 DECLARE_SOA_COLUMN(Chi2, chi2, float);                                                       //! Track chi^2
 DECLARE_SOA_COLUMN(PDca, pDca, float);                                                       //! PDca for MUONStandalone
 DECLARE_SOA_COLUMN(RAtAbsorberEnd, rAtAbsorberEnd, float);                                   //! RAtAbsorberEnd for MUONStandalone tracks and GlobalMuonTrackstracks
@@ -563,7 +563,7 @@ namespace v001
 DECLARE_SOA_DYNAMIC_COLUMN(NClusters, nClusters, //! Number of MFT clusters
                            [](uint64_t mftClusterSizesAndTrackFlags) -> int8_t {
                              int8_t nClusters = 0;
-                             for (int layer = 0; layer < 11; layer++) {
+                             for (int layer = 0; layer < 10; layer++) {
                                if ((mftClusterSizesAndTrackFlags >> (layer * 6)) & 0x3F) {
                                  nClusters++;
                                }
@@ -635,7 +635,7 @@ DECLARE_SOA_TABLE_FULL(StoredMFTTracks_000, "MFTTracks", "AOD", "MFTTRACK", //!
 DECLARE_SOA_TABLE_FULL_VERSIONED(StoredMFTTracks_001, "MFTTracks", "AOD", "MFTTRACK", 1, //! On disk version of MFTTracks, version 1
                                  o2::soa::Index<>, fwdtrack::CollisionId,
                                  fwdtrack::X, fwdtrack::Y, fwdtrack::Z, fwdtrack::Phi, fwdtrack::Tgl,
-                                 fwdtrack::Signed1Pt, fwdtrack::v001::NClusters<fwdtrack::MFTClusterSizesAndTracksFlags>, fwdtrack::MFTClusterSizesAndTracksFlags, fwdtrack::IsCA<fwdtrack::MFTClusterSizesAndTracksFlags>,
+                                 fwdtrack::Signed1Pt, fwdtrack::v001::NClusters<fwdtrack::MFTClusterSizesAndTrackFlags>, fwdtrack::MFTClusterSizesAndTrackFlags, fwdtrack::IsCA<fwdtrack::MFTClusterSizesAndTrackFlags>,
                                  fwdtrack::Px<fwdtrack::Pt, fwdtrack::Phi>,
                                  fwdtrack::Py<fwdtrack::Pt, fwdtrack::Phi>,
                                  fwdtrack::Pz<fwdtrack::Pt, fwdtrack::Tgl>,
@@ -652,8 +652,8 @@ DECLARE_SOA_EXTENDED_TABLE(MFTTracks_001, StoredMFTTracks_001, "MFTTRACK", //! A
                            aod::fwdtrack::Eta,
                            aod::fwdtrack::P);
 
-using MFTTracks = MFTTracks_000;
-using StoredMFTTracks = StoredMFTTracks_000;
+using MFTTracks = MFTTracks_001;
+using StoredMFTTracks = StoredMFTTracks_001;
 
 using MFTTrack = MFTTracks::iterator;
 
@@ -1549,6 +1549,7 @@ DECLARE_EQUIVALENT_FOR_INDEX(aod::StoredTracksExtra_000, aod::StoredTracksExtra_
 DECLARE_EQUIVALENT_FOR_INDEX(aod::HMPID_000, aod::HMPID_001);
 DECLARE_EQUIVALENT_FOR_INDEX(aod::StoredMFTTracks, aod::StoredMFTTracks_000);
 DECLARE_EQUIVALENT_FOR_INDEX(aod::StoredMFTTracks, aod::StoredMFTTracks_001);
+DECLARE_EQUIVALENT_FOR_INDEX(aod::StoredMFTTracks_000, aod::StoredMFTTracks_001);
 } // namespace soa
 
 namespace aod