Merge 0ea9033 into sapling-pr-archive-ktf

Author: ktf

Common/ML/include/ML/3rdparty/GPUORTFloat16.h

@@ -568,9 +568,11 @@ GPUdi() uint16_t BFloat16Impl<Derived>::ToUint16Impl(float v) noexcept
 template <class Derived>
 GPUdi() float BFloat16Impl<Derived>::ToFloatImpl() const noexcept
 {
+#ifndef __FAST_MATH__
   if (IsNaN()) {
     return o2::gpu::CAMath::QuietNaN();
   }
+#endif
   float result;
   char* const first = reinterpret_cast<char*>(&result);
   char* const second = first + sizeof(uint16_t);

DataFormats/Detectors/EMCAL/include/DataFormatsEMCAL/CellLabel.h

@@ -15,6 +15,7 @@
 #include <cstddef>
 #include <cstdint>
 #include <gsl/span>
+#include <vector>
 
 namespace o2
 {
@@ -34,10 +35,15 @@ class CellLabel
  public:
   // CellLabel() = default;
 
-  /// \brief Constructor
+  /// \brief Constructor using std::vector by moving NOT copying
   /// \param labels list of mc labels
   /// \param amplitudeFractions list of amplitude fractions
-  CellLabel(const gsl::span<const int> labels, const gsl::span<const float> amplitudeFractions);
+  CellLabel(std::vector<int> labels, std::vector<float> amplitudeFractions);
+
+  /// \brief Constructor using gsl::span
+  /// \param labels list of mc labels
+  /// \param amplitudeFractions list of amplitude fractions
+  CellLabel(gsl::span<const int> labels, gsl::span<const float> amplitudeFractions);
 
   // ~CellLabel() = default;
   // CellLabel(const CellLabel& clus) = default;
@@ -52,21 +58,21 @@ class CellLabel
   int32_t GetLabel(size_t index) const { return mLabels[index]; }
 
   /// \brief Getter for labels
-  gsl::span<const int32_t> GetLabels() const { return mLabels; }
+  std::vector<int32_t> GetLabels() const { return mLabels; }
 
   /// \brief Getter for amplitude fraction
   /// \param index index which amplitude fraction to get
   float GetAmplitudeFraction(size_t index) const { return mAmplitudeFraction[index]; }
 
   /// \brief Getter for amplitude fractions
-  gsl::span<const float> GetAmplitudeFractions() const { return mAmplitudeFraction; }
+  std::vector<float> GetAmplitudeFractions() const { return mAmplitudeFraction; }
 
   /// \brief Getter for label with leading amplitude fraction
   int32_t GetLeadingMCLabel() const;
 
  protected:
-  gsl::span<const int32_t> mLabels;          ///< List of MC particles that generated the cluster, ordered in deposited energy.
-  gsl::span<const float> mAmplitudeFraction; ///< List of the fraction of the cell energy coming from a MC particle. Index aligns with mLabels!
+  std::vector<int32_t> mLabels;          ///< List of MC particles that generated the cluster, ordered in deposited energy.
+  std::vector<float> mAmplitudeFraction; ///< List of the fraction of the cell energy coming from a MC particle. Index aligns with mLabels!
 };
 
 } // namespace emcal

DataFormats/Detectors/EMCAL/src/CellLabel.cxx

@@ -16,10 +16,19 @@
 #include <cstddef>
 #include <cstdint>
 #include <gsl/span>
+#include <vector>
+#include <utility>
 
 using namespace o2::emcal;
 
-CellLabel::CellLabel(const gsl::span<const int> labels, const gsl::span<const float> amplitudeFractions) : mLabels(labels), mAmplitudeFraction(amplitudeFractions)
+CellLabel::CellLabel(std::vector<int> labels, std::vector<float> amplitudeFractions) : mLabels(std::move(labels)), mAmplitudeFraction(std::move(amplitudeFractions))
+{
+  if (labels.size() != amplitudeFractions.size()) {
+    LOG(error) << "Size of labels " << labels.size() << " does not match size of amplitude fraction " << amplitudeFractions.size() << " !";
+  }
+}
+
+CellLabel::CellLabel(gsl::span<const int> labels, gsl::span<const float> amplitudeFractions) : mLabels(labels.begin(), labels.end()), mAmplitudeFraction(amplitudeFractions.begin(), amplitudeFractions.end())
 {
   if (labels.size() != amplitudeFractions.size()) {
     LOG(error) << "Size of labels " << labels.size() << " does not match size of amplitude fraction " << amplitudeFractions.size() << " !";

Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyConfig.h

@@ -28,7 +28,9 @@ struct TrackMCStudyConfig : o2::conf::ConfigurableParamHelper<TrackMCStudyConfig
   bool requireITSorTPCTrackRefs = true;
   bool requireTopBottomRefs = false;
   int minTPCRefsToExtractClRes = 2;
-  float rejectClustersResStat = 0.;
+  int nOccBinsDrift = 10; // number of bins for TPC max drift time, where we integrate the occupancies
+  int nTBPerOccBin = 48;  // number of TB per occ bin
+  float rejectClustersResStat = 0.1;
   float maxTPCRefExtrap = 2;                   // max dX to extrapolate the track ref when extrapolating track true posions
   int decayPDG[5] = {310, 3122, 411, 421, -1}; // decays to study, must end by -1
   O2ParamDef(TrackMCStudyConfig, "trmcconf");

Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h

@@ -32,7 +32,7 @@ struct MCTrackInfo {
   int getNITSClusForAB() const;
   int getLowestITSLayer() const;
   int getHighestITSLayer() const;
-
+  std::vector<float> occTPCV{};
   o2::track::TrackPar track{};
   o2::MCCompLabel label{};
   float occTPC = -1.f;
@@ -52,7 +52,24 @@ struct MCTrackInfo {
   uint8_t maxTPCRowSect = -1;
   int8_t nITSCl = 0;
   int8_t pattITSCl = 0;
-  ClassDefNV(MCTrackInfo, 4);
+  uint8_t flags = 0;
+
+  enum Flags : uint32_t { Primary = 0,
+                          AddedAtRecStage = 2,
+                          BitMask = 0xff };
+
+  bool isPrimary() const { return isBitSet(Primary); }
+  bool isAddedAtRecStage() const { return isBitSet(AddedAtRecStage); }
+  void setPrimary() { setBit(Primary); }
+  void setAddedAtRecStage() { setBit(AddedAtRecStage); }
+
+  uint8_t getBits() const { return flags; }
+  bool isBitSet(int bit) const { return flags & (0xff & (0x1 << bit)); }
+  void setBits(std::uint8_t b) { flags = b; }
+  void setBit(int bit) { flags |= BitMask & (0x1 << bit); }
+  void resetBit(int bit) { flags &= ~(BitMask & (0x1 << bit)); }
+
+  ClassDefNV(MCTrackInfo, 7);
 };
 
 struct RecTrack {
@@ -75,6 +92,7 @@ struct RecTrack {
   uint8_t nClTPC = 0;
   uint8_t pattITS = 0;
   int8_t lowestPadRow = -1;
+  int8_t padFromEdge = -1;
 
   bool isFakeGLO() const { return flags & FakeGLO; }
   bool isFakeITS() const { return flags & FakeITS; }
@@ -83,7 +101,7 @@ struct RecTrack {
   bool isFakeTOF() const { return flags & FakeTOF; }
   bool isFakeITSTPC() const { return flags & FakeITSTPC; }
 
-  ClassDefNV(RecTrack, 1);
+  ClassDefNV(RecTrack, 2);
 };
 
 struct TrackPairInfo {
@@ -270,7 +288,8 @@ struct MCVertex {
   int nTrackSel = 0; // number of selected MC charged tracks
   int ID = -1;
   std::vector<RecPV> recVtx{};
-  ClassDefNV(MCVertex, 1);
+  std::vector<float> occTPCV{};
+  ClassDefNV(MCVertex, 2);
 };
 
 } // namespace o2::trackstudy

Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx

@@ -109,7 +109,7 @@ class TrackMCStudy : public Task
   void updateTimeDependentParams(ProcessingContext& pc);
   float getDCAYCut(float pt) const;
 
-  gsl::span<const MCTrack> mCurrMCTracks;
+  const std::vector<o2::MCTrack>* mCurrMCTracks = nullptr;
   TVector3 mCurrMCVertex;
   o2::tpc::VDriftHelper mTPCVDriftHelper{};
   o2::tpc::CorrectionMapsLoader mTPCCorrMapsLoader{};
@@ -122,8 +122,9 @@ class TrackMCStudy : public Task
   std::vector<float> mTPCOcc;    ///< TPC occupancy for this interaction time
   std::vector<int> mITSOcc;      //< N ITS clusters in the ROF containing collision
   bool mCheckSV = false;         //< check SV binding (apart from prongs availability)
+  bool mRecProcStage = false;    //< flag that the MC particle was added only at the stage of reco tracks processing
   int mNTPCOccBinLength = 0;     ///< TPC occ. histo bin length in TBs
-  float mNTPCOccBinLengthInv;
+  float mNTPCOccBinLengthInv = -1.f;
   int mVerbose = 0;
   float mITSTimeBiasMUS = 0.f;
   float mITSROFrameLengthMUS = 0.f; ///< ITS RO frame in mus
@@ -181,10 +182,11 @@ void TrackMCStudy::run(ProcessingContext& pc)
   }
   mDecProdLblPool.clear();
   mMCVtVec.clear();
-  mCurrMCTracks = {};
+  mCurrMCTracks = nullptr;
 
   recoData.collectData(pc, *mDataRequest.get()); // select tracks of needed type, with minimal cuts, the real selected will be done in the vertexer
   updateTimeDependentParams(pc);                 // Make sure this is called after recoData.collectData, which may load some conditions
+  mRecProcStage = false;
   process(recoData);
 }
 
@@ -278,15 +280,21 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
     return patt;
   };
 
-  auto getLowestPadrow = [&recoData](const o2::tpc::TrackTPC& trc) {
+  auto getLowestPadrow = [&recoData](const o2::tpc::TrackTPC& trc, RecTrack& tref) {
     if (recoData.inputsTPCclusters) {
       uint8_t clSect = 0, clRow = 0;
       uint32_t clIdx = 0;
       const auto clRefs = recoData.getTPCTracksClusterRefs();
+      const auto tpcClusAcc = recoData.getTPCClusters();
       trc.getClusterReference(clRefs, trc.getNClusterReferences() - 1, clSect, clRow, clIdx);
-      return int(clRow);
+      const auto& clus = tpcClusAcc.clusters[clSect][clRow][clIdx];
+      int padFromEdge = int(clus.getPad()), npads = o2::gpu::GPUTPCGeometry::NPads(clRow);
+      if (padFromEdge > npads / 2) {
+        padFromEdge = npads - 1 - padFromEdge;
+      }
+      tref.padFromEdge = uint8_t(padFromEdge);
+      tref.lowestPadRow = clRow;
     }
-    return -1;
   };
 
   auto flagTPCClusters = [&recoData](const o2::tpc::TrackTPC& trc, o2::MCCompLabel lbTrc) {
@@ -338,6 +346,21 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
         }
         break;
       }
+      if (mNTPCOccBinLengthInv > 0.f) {
+        mcVtx.occTPCV.resize(params.nOccBinsDrift);
+        int grp = TMath::Max(1, TMath::Nint(params.nTBPerOccBin * mNTPCOccBinLengthInv));
+        for (int ib = 0; ib < params.nOccBinsDrift; ib++) {
+          float smb = 0;
+          int tbs = occBin + TMath::Nint(ib * params.nTBPerOccBin * mNTPCOccBinLengthInv);
+          for (int ig = 0; ig < grp; ig++) {
+            if (tbs >= 0 && tbs < int(mTBinClOccHist.size())) {
+              smb += mTBinClOccHist[tbs];
+            }
+            tbs++;
+          }
+          mcVtx.occTPCV[ib] = smb;
+        }
+      }
       if (rofCount >= ITSClusROFRec.size()) {
         mITSOcc.push_back(0); // IR after the last ROF
       }
@@ -352,15 +375,17 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
     int nev = mcReader.getNEvents(curSrcMC);
     bool okAccVtx = true;
     if (nev != (int)mMCVtVec.size()) {
-      LOGP(error, "source {} has {} events while {} MC vertices were booked", curSrcMC, nev, mMCVtVec.size());
+      LOGP(debug, "source {} has {} events while {} MC vertices were booked", curSrcMC, nev, mMCVtVec.size());
       okAccVtx = false;
+      if (nev > (int)mMCVtVec.size()) { // QED
+        continue;
+      }
     }
     for (curEvMC = 0; curEvMC < nev; curEvMC++) {
       if (mVerbose > 1) {
         LOGP(info, "Event {}", curEvMC);
       }
-      const auto& mt = mcReader.getTracks(curSrcMC, curEvMC);
-      mCurrMCTracks = gsl::span<const MCTrack>(mt.data(), mt.size());
+      mCurrMCTracks = &mcReader.getTracks(curSrcMC, curEvMC);
       const_cast<o2::dataformats::MCEventHeader&>(mcReader.getMCEventHeader(curSrcMC, curEvMC)).GetVertex(mCurrMCVertex);
       if (okAccVtx) {
         auto& pos = mMCVtVec[curEvMC].pos;
@@ -370,7 +395,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
           pos[2] = mCurrMCVertex.Z();
         }
       }
-      for (int itr = 0; itr < mCurrMCTracks.size(); itr++) {
+      for (int itr = 0; itr < mCurrMCTracks->size(); itr++) {
         processMCParticle(curSrcMC, curEvMC, itr);
       }
     }
@@ -382,6 +407,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
   }
 
   // add reconstruction info to MC particles. If MC particle was not selected before but was reconstrected, account MC info
+  mRecProcStage = true; // MC particles accepted only at this stage will be flagged
   for (int iv = 0; iv < nv; iv++) {
     if (mVerbose > 1) {
       LOGP(info, "processing PV {} of {}", iv, nv);
@@ -416,11 +442,10 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
             if (lbl.getSourceID() != curSrcMC || lbl.getEventID() != curEvMC) {
               curSrcMC = lbl.getSourceID();
               curEvMC = lbl.getEventID();
-              const auto& mt = mcReader.getTracks(curSrcMC, curEvMC);
-              mCurrMCTracks = gsl::span<const MCTrack>(mt.data(), mt.size());
+              mCurrMCTracks = &mcReader.getTracks(curSrcMC, curEvMC);
               const_cast<o2::dataformats::MCEventHeader&>(mcReader.getMCEventHeader(curSrcMC, curEvMC)).GetVertex(mCurrMCVertex);
             }
-            if (!acceptMCCharged(mCurrMCTracks[lbl.getTrackID()], lbl)) {
+            if (!acceptMCCharged((*mCurrMCTracks)[lbl.getTrackID()], lbl)) {
               continue;
             }
             entry = mSelMCTracks.find(lbl);
@@ -532,7 +557,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
         if (msk[DetID::TPC]) {
           const auto& trtpc = recoData.getTPCTrack(gidSet[GTrackID::TPC]);
           tref.nClTPC = trtpc.getNClusters();
-          tref.lowestPadRow = getLowestPadrow(trtpc);
+          getLowestPadrow(trtpc, tref);
           flagTPCClusters(trtpc, entry.first);
           if (trackFam.entTPC < 0) {
             trackFam.entTPC = tcnt;
@@ -968,7 +993,7 @@ float TrackMCStudy::getDCAYCut(float pt) const
 
 bool TrackMCStudy::processMCParticle(int src, int ev, int trid)
 {
-  const auto& mcPart = mCurrMCTracks[trid];
+  const auto& mcPart = (*mCurrMCTracks)[trid];
   int pdg = mcPart.GetPdgCode();
   bool res = false;
   while (true) {
@@ -990,7 +1015,7 @@ bool TrackMCStudy::processMCParticle(int src, int ev, int trid)
           break;
         }
         for (int idd = idd0; idd <= idd1; idd++) {
-          const auto& product = mCurrMCTracks[idd];
+          const auto& product = (*mCurrMCTracks)[idd];
           auto lbld = o2::MCCompLabel(idd, ev, src);
           if (!acceptMCCharged(product, lbld, decay)) {
             decay = -1; // discard decay
@@ -1088,10 +1113,17 @@ bool TrackMCStudy::addMCParticle(const MCTrack& mcPart, const o2::MCCompLabel& l
   mcEntry.mcTrackInfo.bcInTF = mIntBC[lb.getEventID()];
   mcEntry.mcTrackInfo.occTPC = mTPCOcc[lb.getEventID()];
   mcEntry.mcTrackInfo.occITS = mITSOcc[lb.getEventID()];
+  mcEntry.mcTrackInfo.occTPCV = mMCVtVec[lb.getEventID()].occTPCV;
+  if (mRecProcStage) {
+    mcEntry.mcTrackInfo.setAddedAtRecStage();
+  }
+  if (o2::mcutils::MCTrackNavigator::isPhysicalPrimary(mcPart, *mCurrMCTracks)) {
+    mcEntry.mcTrackInfo.setPrimary();
+  }
   int moth = -1;
   o2::MCCompLabel mclbPar;
   if ((moth = mcPart.getMotherTrackId()) >= 0) {
-    const auto& mcPartPar = mCurrMCTracks[moth];
+    const auto& mcPartPar = (*mCurrMCTracks)[moth];
     mcEntry.mcTrackInfo.pdgParent = mcPartPar.GetPdgCode();
   }
   if (mcPart.isPrimary() && mcReader.getNEvents(lb.getSourceID()) == mMCVtVec.size()) {

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h

@@ -91,6 +91,7 @@ struct TimeFrame {
                       gsl::span<const unsigned char>::iterator& pattIt,
                       const itsmft::TopologyDictionary* dict,
                       const dataformats::MCTruthContainer<MCCompLabel>* mcLabels = nullptr);
+  void resetROFrameData();
 
   int getTotalClusters() const;
   auto& getTotVertIteration() { return mTotVertPerIteration; }
@@ -138,7 +139,7 @@ struct TimeFrame {
   gsl::span<const MCCompLabel> getClusterLabels(int layerId, const int clId) const { return mClusterLabels->getLabels(mClusterExternalIndices[layerId][clId]); }
   int getClusterExternalIndex(int layerId, const int clId) const { return mClusterExternalIndices[layerId][clId]; }
   int getClusterSize(int clusterId) const { return mClusterSize[clusterId]; }
-  void setClusterSize(const bounded_vector<uint8_t>& v) { mClusterSize = v; }
+  void setClusterSize(bounded_vector<uint8_t>& v) { mClusterSize = std::move(v); }
 
   auto& getTrackletsLabel(int layer) { return mTrackletLabels[layer]; }
   auto& getCellsLabel(int layer) { return mCellLabels[layer]; }

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h

@@ -27,6 +27,8 @@
 #include <utility>
 #include <sstream>
 
+#include <oneapi/tbb/task_arena.h>
+
 #include "ITStracking/Configuration.h"
 #include "CommonConstants/MathConstants.h"
 #include "ITStracking/Definitions.h"
@@ -73,8 +75,7 @@ class Tracker
   void setBz(float bz) { mTraits->setBz(bz); }
   void setCorrType(const o2::base::PropagatorImpl<float>::MatCorrType type) { mTraits->setCorrType(type); }
   bool isMatLUT() const { return mTraits->isMatLUT(); }
-  void setNThreads(int n) { mTraits->setNThreads(n); }
-  int getNThreads() const { return mTraits->getNThreads(); }
+  void setNThreads(int n, std::shared_ptr<tbb::task_arena>& arena) { mTraits->setNThreads(n, arena); }
   void printSummary() const;
 
  private:

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h

@@ -80,8 +80,8 @@ class TrackerTraits
   void SetRecoChain(o2::gpu::GPUChainITS* chain) { mChain = chain; }
   void setSmoothing(bool v) { mApplySmoothing = v; }
   bool getSmoothing() const { return mApplySmoothing; }
-  void setNThreads(int n);
-  int getNThreads() const { return mNThreads; }
+  void setNThreads(int n, std::shared_ptr<tbb::task_arena>& arena);
+  int getNThreads() { return mTaskArena->max_concurrency(); }
 
   o2::gpu::GPUChainITS* getChain() const { return mChain; }
 
@@ -94,10 +94,9 @@ class TrackerTraits
   track::TrackParCov buildTrackSeed(const Cluster& cluster1, const Cluster& cluster2, const TrackingFrameInfo& tf3);
   bool fitTrack(TrackITSExt& track, int start, int end, int step, float chi2clcut = o2::constants::math::VeryBig, float chi2ndfcut = o2::constants::math::VeryBig, float maxQoverPt = o2::constants::math::VeryBig, int nCl = 0);
 
-  int mNThreads = 1;
   bool mApplySmoothing = false;
   std::shared_ptr<BoundedMemoryResource> mMemoryPool;
-  tbb::task_arena mTaskArena;
+  std::shared_ptr<tbb::task_arena> mTaskArena;
 
  protected:
   o2::base::PropagatorImpl<float>::MatCorrType mCorrType = o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrNONE;