GPU: Process dEdx with full qTot range

Author: davidrohr

DataFormats/Detectors/TPC/include/DataFormatsTPC/ClusterNative.h

@@ -73,28 +73,21 @@ struct ClusterNative {
   uint8_t sigmaTimePacked;  //< Sigma of the time in packed format
   uint8_t sigmaPadPacked;   //< Sigma of the pad in packed format
   uint16_t qMax;            //< QMax of the cluster
-  uint16_t qTot;            //< Total charge of the cluster
+  uint16_t qTotPacked;      //< Total charge of the cluster
 
   GPUd() static uint16_t packPad(float pad) { return (uint16_t)(pad * scalePadPacked + 0.5); }
   GPUd() static uint32_t packTime(float time) { return (uint32_t)(time * scaleTimePacked + 0.5); }
   GPUd() static float unpackPad(uint16_t pad) { return float(pad) * (1.f / scalePadPacked); }
   GPUd() static float unpackTime(uint32_t time) { return float(time) * (1.f / scaleTimePacked); }
 
   GPUdDefault() ClusterNative() = default;
-  GPUd() ClusterNative(uint32_t time, uint8_t flags, uint16_t pad, uint8_t sigmaTime, uint8_t sigmaPad, uint16_t qmax, uint16_t qtot) : padPacked(pad), sigmaTimePacked(sigmaTime), sigmaPadPacked(sigmaPad), qMax(qmax), qTot(qtot)
+  GPUd() ClusterNative(uint32_t time, uint8_t flags, uint16_t pad, uint8_t sigmaTime, uint8_t sigmaPad, uint16_t qmax, uint16_t qtotPacked) : padPacked(pad), sigmaTimePacked(sigmaTime), sigmaPadPacked(sigmaPad), qMax(qmax), qTotPacked(qtotPacked)
   {
     setTimePackedFlags(time, flags);
   }
 
   GPUd() uint16_t getQmax() const { return qMax; }
-  GPUd() uint16_t getQtot() const
-  {
-    if (isSaturated()) [[unlikely]] {
-      auto sQtot = getSaturatedQtot();
-      return sQtot < USHRT_MAX ? sQtot : USHRT_MAX;
-    }
-    return qTot;
-  }
+  GPUd() uint32_t getQtot() const { return isSaturated() ? getSaturatedQtot() : (uint32_t)qTotPacked; }
   GPUd() uint8_t getFlags() const { return timeFlagsPacked >> 24; }
   GPUd() uint32_t getTimePacked() const { return timeFlagsPacked & 0xFFFFFF; }
   GPUd() void setTimePackedFlags(uint32_t timePacked, uint8_t flags)
@@ -155,19 +148,19 @@ struct ClusterNative {
     sigmaPadPacked = tmp;
   }
 
-  GPUd() bool isSaturated() const { return qTot > maxRegularQtot; }
+  GPUd() bool isSaturated() const { return qTotPacked > maxRegularQtot; }
 
   GPUd() void setSaturatedQtot(uint32_t qtot)
   {
-    this->qTot = USHRT_MAX;
+    this->qTotPacked = USHRT_MAX;
     if (qtot < maxSaturatedQtot) {
-      this->qTot = ((qtot + scaleSaturatedQtot / 2) / scaleSaturatedQtot) + maxRegularQtot;
+      this->qTotPacked = ((qtot + scaleSaturatedQtot / 2) / scaleSaturatedQtot) + maxRegularQtot;
     }
   }
 
   GPUd() uint32_t getSaturatedQtot() const
   {
-    return uint32_t(qTot - maxRegularQtot) * scaleSaturatedQtot;
+    return uint32_t(qTotPacked - maxRegularQtot) * scaleSaturatedQtot;
   }
 
   GPUd() void setSaturatedTailLength(uint32_t tail)
@@ -192,8 +185,8 @@ struct ClusterNative {
       return (this->sigmaPadPacked < rhs.sigmaPadPacked);
     } else if (this->qMax != rhs.qMax) {
       return (this->qMax < rhs.qMax);
-    } else if (this->qTot != rhs.qTot) {
-      return (this->qTot < rhs.qTot);
+    } else if (this->qTotPacked != rhs.qTotPacked) {
+      return (this->getQtot() < rhs.getQtot());
     } else {
       return (this->getFlags() < rhs.getFlags());
     }
@@ -206,7 +199,7 @@ struct ClusterNative {
            this->sigmaTimePacked == rhs.sigmaTimePacked &&
            this->sigmaPadPacked == rhs.sigmaPadPacked &&
            this->qMax == rhs.qMax &&
-           this->qTot == rhs.qTot &&
+           this->qTotPacked == rhs.qTotPacked &&
            this->getFlags() == rhs.getFlags();
   }
 

DataFormats/Detectors/TPC/include/DataFormatsTPC/ClusterNativeHelper.h

@@ -312,7 +312,7 @@ class ClusterNativeHelper
         sigmaTime = rhs.getSigmaTime();
         sigmaPad = rhs.getSigmaPad();
         qMax = rhs.qMax;
-        qTot = rhs.qTot;
+        qTot = rhs.qTotPacked;
         flags = rhs.getFlags();
         return *this;
       }

Detectors/Align/src/AlignableDetectorTPC.cxx

@@ -166,18 +166,18 @@ int AlignableDetectorTPC::processPoints(GIndex gid, int npntCut, bool inv)
           // mController->getTPCCorrMaps()->Transform(sector, row, cl->getPad(), cl->getTime(), x, y, z, tOffset);
           currentRow = row;
           currentSector = sector;
-          charge = cl->qTot;
+          charge = cl->getQtot();
           clusterState = nextState;
           combRow = row;
           LOGP(debug, "starting a supercluster at row {} of sector {} -> {},{},{}", currentRow, currentSector, x, y, z);
         } else {
           // float xx, yy, zz;
           // mController->getTPCCorrMaps()->Transform(sector, row, cl->getPad(), cl->getTime(), xx, yy, zz, tOffset);
-          x += xTmp * cl->qTot;
-          y += yTmp * cl->qTot;
-          z += zTmp * cl->qTot;
-          combRow += row * cl->qTot;
-          charge += cl->qTot;
+          x += xTmp * cl->getQtot();
+          y += yTmp * cl->getQtot();
+          z += zTmp * cl->getQtot();
+          combRow += row * cl->getQtot();
+          charge += cl->getQtot();
           clusterState |= nextState;
           npntCut--;
           LOGP(debug, "merging cluster #{} at row {} to a supercluster starting at row {} ", clusters + 1, row, currentRow);

Detectors/TPC/calibration/src/CalculatedEdx.cxx

@@ -245,7 +245,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
     }
 
     // get charge values
-    float chargeTot = cl.qTot;
+    float chargeTot = cl.getQtot();
     float chargeMax = cl.qMax;
 
     // get threshold

Detectors/TPC/calibration/src/CalibPadGainTracks.cxx

@@ -100,7 +100,7 @@ void CalibPadGainTracks::processTrack(o2::tpc::TrackTPC track, o2::gpu::GPUO2Int
     }
 
     const int region = Mapper::REGION[rowIndex];
-    const float charge = (mChargeType == ChargeType::Max) ? cl.qMax : cl.qTot;
+    const float charge = (mChargeType == ChargeType::Max) ? cl.qMax : cl.getQtot();
     const float effectiveLength = mCalibTrackTopologyPol ? getTrackTopologyCorrectionPol(track, cl, region, charge) : getTrackTopologyCorrection(track, region);
 
     const unsigned char pad = std::clamp(static_cast<unsigned int>(cl.getPad() + 0.5f), static_cast<unsigned int>(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); // the left side of the pad is defined at e.g. 3.5 and the right side at 4.5

Detectors/TPC/calibration/src/TrackDump.cxx

@@ -73,7 +73,7 @@ void TrackDump::filter(const gsl::span<const TrackTPC> tracks, ClusterNativeAcce
       excludes[sector][padrow].emplace_back(clusterIndexInRow);
 
       if (clustersGlobal) {
-        auto& clGlobal = clustersGlobal->emplace_back(ClusterGlobal{clInfo.gx(), clInfo.gy(), cl.qMax, cl.qTot, sector, padrow});
+        auto& clGlobal = clustersGlobal->emplace_back(ClusterGlobal{clInfo.gx(), clInfo.gy(), cl.qMax, cl.getQtot(), sector, padrow});
       }
     }
   }

Detectors/TPC/reconstruction/src/HardwareClusterDecoder.cxx

@@ -87,7 +87,7 @@ int HardwareClusterDecoder::decodeClusters(std::vector<std::pair<const ClusterHa
             cOut.setSigmaPad(std::sqrt(cIn.getSigmaPad2()));
             cOut.setSigmaTime(std::sqrt(cIn.getSigmaTime2()));
             cOut.qMax = cIn.getQMax();
-            cOut.qTot = cIn.getQTot();
+            cOut.qTotPacked = cIn.getQTot();
             mIntegrator->integrateCluster(sector, padRowGlobal, pad, cIn.getQTot());
             if (outMCLabels) {
               auto& mcOut = outMCLabelContainers[containerRowCluster[sector][padRowGlobal]];

Detectors/TPC/reconstruction/test/testGPUCATracking.cxx

@@ -95,7 +95,7 @@ BOOST_AUTO_TEST_CASE(CATracking_test1)
     cont[i].clusters[0].setSigmaTime(1);
     cont[i].clusters[0].setSigmaPad(1);
     cont[i].clusters[0].qMax = 10;
-    cont[i].clusters[0].qTot = 50;
+    cont[i].clusters[0].qTotPacked = 50;
   }
   std::unique_ptr<ClusterNative[]> clusterBuffer;
   std::unique_ptr<ClusterNativeAccess> clusters = ClusterNativeHelper::createClusterNativeIndex(clusterBuffer, cont, nullptr, nullptr);

GPU/GPUTracking/Base/GPUReconstructionConvert.cxx

@@ -48,43 +48,6 @@ using namespace o2::tpc;
 using namespace o2::tpc::constants;
 using namespace std::string_literals;
 
-void GPUReconstructionConvert::ConvertNativeToClusterData(o2::tpc::ClusterNativeAccess* native, std::unique_ptr<GPUTPCClusterData[]>* clusters, uint32_t* nClusters, const TPCFastTransformPOD* transform, int32_t continuousMaxTimeBin)
-{
-  memset(nClusters, 0, NSECTORS * sizeof(nClusters[0]));
-  uint32_t offset = 0;
-  for (uint32_t i = 0; i < NSECTORS; i++) {
-    uint32_t nClSector = 0;
-    for (uint32_t j = 0; j < GPUTPCGeometry::NROWS; j++) {
-      nClSector += native->nClusters[i][j];
-    }
-    nClusters[i] = nClSector;
-    clusters[i].reset(new GPUTPCClusterData[nClSector]);
-    nClSector = 0;
-    for (uint32_t j = 0; j < GPUTPCGeometry::NROWS; j++) {
-      for (uint32_t k = 0; k < native->nClusters[i][j]; k++) {
-        const auto& clin = native->clusters[i][j][k];
-        float x = 0, y = 0, z = 0;
-        if (continuousMaxTimeBin == 0) {
-          transform->Transform(i, j, clin.getPad(), clin.getTime(), x, y, z);
-        } else {
-          transform->TransformInTimeFrame(i, j, clin.getPad(), clin.getTime(), x, y, z, continuousMaxTimeBin);
-        }
-        auto& clout = clusters[i].get()[nClSector];
-        clout.x = x;
-        clout.y = y;
-        clout.z = z;
-        clout.row = j;
-        clout.amp = clin.qTot;
-        clout.flags = clin.getFlags();
-        clout.id = offset + k;
-        nClSector++;
-      }
-      native->clusterOffset[i][j] = offset;
-      offset += native->nClusters[i][j];
-    }
-  }
-}
-
 void GPUReconstructionConvert::ConvertRun2RawToNative(o2::tpc::ClusterNativeAccess& native, std::unique_ptr<ClusterNative[]>& nativeBuffer, const AliHLTTPCRawCluster** rawClusters, uint32_t* nRawClusters)
 {
   memset((void*)&native, 0, sizeof(native));
@@ -110,7 +73,7 @@ void GPUReconstructionConvert::ConvertRun2RawToNative(o2::tpc::ClusterNativeAcce
       c.setSigmaTime(CAMath::Sqrt(org.GetSigmaTime2()));
       c.setSigmaPad(CAMath::Sqrt(org.GetSigmaPad2()));
       c.qMax = org.GetQMax();
-      c.qTot = org.GetCharge();
+      c.qTotPacked = org.GetCharge();
     }
   }
 }

GPU/GPUTracking/Base/GPUReconstructionConvert.h

@@ -50,7 +50,6 @@ class GPUReconstructionConvert
 {
  public:
   constexpr static uint32_t NSECTORS = o2::tpc::constants::MAXSECTOR;
-  static void ConvertNativeToClusterData(o2::tpc::ClusterNativeAccess* native, std::unique_ptr<GPUTPCClusterData[]>* clusters, uint32_t* nClusters, const TPCFastTransformPOD* transform, int32_t continuousMaxTimeBin = 0);
   static void ConvertRun2RawToNative(o2::tpc::ClusterNativeAccess& native, std::unique_ptr<o2::tpc::ClusterNative[]>& nativeBuffer, const AliHLTTPCRawCluster** rawClusters, uint32_t* nRawClusters);
   template <class S>
   static void RunZSEncoder(const S& in, std::unique_ptr<uint64_t[]>* outBuffer, uint32_t* outSizes, o2::raw::RawFileWriter* raw, const o2::InteractionRecord* ir, const GPUParam& param, int32_t version, bool verify, float threshold = 0.f, bool padding = false, std::function<void(std::vector<o2::tpc::Digit>&)> digitsFilter = nullptr);