TPC Splines: init inverse from the inverse voxel map; rebase

Author: sgorbuno

Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h

@@ -86,15 +86,14 @@ class TPCFastSpaceChargeCorrectionHelper
 
   /// Create SpaceCharge correction out of the voxel tree
   std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> createFromTrackResiduals(
-    const o2::tpc::TrackResiduals& trackResiduals, TTree* voxResTree, bool useSmoothed = false, bool invertSigns = false);
+    const o2::tpc::TrackResiduals& trackResiduals, TTree* voxResTree, TTree* voxResTreeInverse, bool useSmoothed, bool invertSigns);
+
   /// _______________  Utilities   ________________________
 
   const TPCFastTransformGeo& getGeometry() { return mGeo; }
 
   TPCFastSpaceChargeCorrectionMap& getCorrectionMap() { return mCorrectionMap; }
 
-  void fillSpaceChargeCorrectionFromMap(TPCFastSpaceChargeCorrection& correction);
-
   void testGeometry(const TPCFastTransformGeo& geo) const;
 
   /// initialise inverse transformation
@@ -103,15 +102,13 @@ class TPCFastSpaceChargeCorrectionHelper
   /// initialise inverse transformation from linear combination of several input corrections
   void initInverse(std::vector<o2::gpu::TPCFastSpaceChargeCorrection*>& corrections, const std::vector<float>& scaling, bool prn);
 
+  void MergeCorrections(std::vector<o2::gpu::TPCFastSpaceChargeCorrection*>& corrections, const std::vector<float>& scaling, bool prn);
+
  private:
   /// geometry initialization
   void initGeometry();
 
-  /// get space charge correction in internal TPCFastTransform coordinates u,v->dx,du,dv
-  void getSpaceChargeCorrection(const TPCFastSpaceChargeCorrection& correction, int slice, int row, o2::gpu::TPCFastSpaceChargeCorrectionMap::CorrectionPoint p, double& su, double& sv, double& dx, double& du, double& dv);
-
-  /// initialise max drift length
-  void initMaxDriftLength(o2::gpu::TPCFastSpaceChargeCorrection& correction, bool prn);
+  void fillSpaceChargeCorrectionFromMap(TPCFastSpaceChargeCorrection& correction, bool processingInverseCorrection);
 
   static TPCFastSpaceChargeCorrectionHelper* sInstance; ///< singleton instance
   bool mIsInitialized = 0;                              ///< initialization flag

Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx

@@ -177,8 +177,8 @@ void TPCFastTransformHelperO2::testGeometry(const TPCFastTransformGeo& geo) cons
 {
   const Mapper& mapper = Mapper::instance();
 
-  if (geo.getNumberOfRocs() != Sector::MAXSECTOR) {
-    LOG(fatal) << "Wrong number of sectors :" << geo.getNumberOfRocs() << " instead of " << Sector::MAXSECTOR << std::endl;
+  if (geo.getNumberOfSectors() != Sector::MAXSECTOR) {
+    LOG(fatal) << "Wrong number of sectors :" << geo.getNumberOfSectors() << " instead of " << Sector::MAXSECTOR << std::endl;
   }
 
   if (geo.getNumberOfRows() != mapper.getNumberOfRows()) {

Detectors/TPC/reconstruction/src/TPCFastTransformHelperO2.cxx

@@ -53,7 +53,7 @@ BOOST_AUTO_TEST_CASE(FastTransform_test1)
 
   BOOST_CHECK_EQUAL(geo.test(), 0);
 
-  BOOST_CHECK_EQUAL(geo.getNumberOfRocs(), Sector::MAXSECTOR);
+  BOOST_CHECK_EQUAL(geo.getNumberOfSectors(), Sector::MAXSECTOR);
   BOOST_CHECK_EQUAL(geo.getNumberOfRows(), mapper.getNumberOfRows());
 
   double maxDx = 0, maxDy = 0;
@@ -71,15 +71,16 @@ BOOST_AUTO_TEST_CASE(FastTransform_test1)
     for (int pad = 0; pad < nPads; pad++) {
       const GlobalPadNumber p = mapper.globalPadNumber(PadPos(row, pad));
       const PadCentre& c = mapper.padCentre(p);
-      float u = 0, v = 0;
-      fastTransform.convPadTimeToUV(row, pad, 0, u, v, 0.);
-
+      float y = 0, z = 0;
+      int sector = 0;
+      float time = 0.;
+      fastTransform.convPadTimeToLocal(sector, row, pad, time, y, z, 0.);
       double dx = x - c.X();
-      double dy = u - (-c.Y()); // diferent sign convention for Y coordinate in the map
+      double dy = y - (-c.Y()); // diferent sign convention for Y coordinate in the map
       BOOST_CHECK(fabs(dx) < 1.e-6);
       BOOST_CHECK(fabs(dy) < 1.e-5);
       if (fabs(dy) >= 1.e-5) {
-        std::cout << "row " << row << " pad " << pad << " y calc " << u << " y in map " << -c.Y() << " dy " << dy << std::endl;
+        std::cout << "row " << row << " pad " << pad << " y calc " << y << " y in map " << -c.Y() << " dy " << dy << std::endl;
       }
       if (fabs(maxDx) < fabs(dx)) {
         maxDx = dx;
@@ -104,46 +105,46 @@ BOOST_AUTO_TEST_CASE(FastTransform_test_setSpaceChargeCorrection)
   std::unique_ptr<TPCFastTransform> fastTransform0(TPCFastTransformHelperO2::instance()->create(0));
   const TPCFastTransformGeo& geo = fastTransform0->getGeometry();
 
-  auto correctionUV = [&](int roc, int /*row*/, const double u, const double v, double& dX, double& dU, double& dV) {
+  auto correctionUV = [&](int sector, int /*row*/, const double u, const double v, double& dX, double& dU, double& dV) {
     // float lx = geo.getRowInfo(row).x;
     dX = 1. + 1 * u + 0.1 * u * u;
     dU = 2. + 0.2 * u + 0.002 * u * u; // + 0.001 * u * u * u;
     dV = 3. + 0.1 * v + 0.01 * v * v;  //+ 0.0001 * v * v * v;
   };
 
-  auto correctionLocal = [&](int roc, int row, double ly, double lz,
+  auto correctionLocal = [&](int sector, int row, double ly, double lz,
                              double& dx, double& dly, double& dlz) {
     float u, v;
-    geo.convLocalToUV(roc, ly, lz, u, v);
+    geo.convLocalToUV(sector, ly, lz, u, v);
     double du, dv;
-    correctionUV(roc, row, u, v, dx, du, dv);
+    correctionUV(sector, row, u, v, dx, du, dv);
     float ly1, lz1;
-    geo.convUVtoLocal(roc, u + du, v + dv, ly1, lz1);
+    geo.convUVtoLocal(sector, u + du, v + dv, ly1, lz1);
     dly = ly1 - ly;
     dlz = lz1 - lz;
   };
 
-  int nRocs = geo.getNumberOfRocs();
+  int nSectors = geo.getNumberOfSectors();
   int nRows = geo.getNumberOfRows();
   TPCFastSpaceChargeCorrectionMap& scData = TPCFastTransformHelperO2::instance()->getCorrectionMap();
-  scData.init(nRocs, nRows);
+  scData.init(nSectors, nRows);
 
-  for (int iRoc = 0; iRoc < nRocs; iRoc++) {
+  for (int iSector = 0; iSector < nSectors; iSector++) {
     for (int iRow = 0; iRow < nRows; iRow++) {
       double dsu = 1. / (3 * 8 - 3);
       double dsv = 1. / (3 * 20 - 3);
       for (double su = 0.f; su < 1.f + .5 * dsu; su += dsv) {
         for (double sv = 0.f; sv < 1.f + .5 * dsv; sv += dsv) {
           float ly = 0.f, lz = 0.f;
-          geo.convScaledUVtoLocal(iRoc, iRow, su, sv, ly, lz);
+          geo.convScaledUVtoLocal(iSector, iRow, su, sv, ly, lz);
           double dx, dy, dz;
-          correctionLocal(iRoc, iRow, ly, lz, dx, dy, dz);
-          scData.addCorrectionPoint(iRoc, iRow,
+          correctionLocal(iSector, iRow, ly, lz, dx, dy, dz);
+          scData.addCorrectionPoint(iSector, iRow,
                                     ly, lz, dx, dy, dz);
         }
       }
     } // row
-  } // roc
+  } // sector
 
   std::unique_ptr<TPCFastTransform> fastTransform(TPCFastTransformHelperO2::instance()->create(0));
 
@@ -158,12 +159,12 @@ BOOST_AUTO_TEST_CASE(FastTransform_test_setSpaceChargeCorrection)
   double statDiff = 0., statN = 0.;
   double statDiffFile = 0., statNFile = 0.;
 
-  for (int roc = 0; roc < geo.getNumberOfRocs(); roc += 1) {
-    // std::cout << "roc " << roc << " ... " << std::endl;
+  for (int sector = 0; sector < geo.getNumberOfSectors(); sector += 1) {
+    // std::cout << "sector " << sector << " ... " << std::endl;
 
-    const TPCFastTransformGeo::RocInfo& rocInfo = geo.getRocInfo(roc);
+    const TPCFastTransformGeo::SectorInfo& sectorInfo = geo.getSectorInfo(sector);
 
-    float lastTimeBin = fastTransform->getMaxDriftTime(roc, 0.f);
+    float lastTimeBin = fastTransform->getMaxDriftTime(sector, 0.f);
 
     for (int row = 0; row < geo.getNumberOfRows(); row++) {
 
@@ -172,31 +173,31 @@ BOOST_AUTO_TEST_CASE(FastTransform_test_setSpaceChargeCorrection)
       for (int pad = 0; pad < nPads; pad += 10) {
 
         for (float time = 0; time < lastTimeBin; time += 30) {
-          // std::cout<<"roc "<<roc<<" row "<<row<<" pad "<<pad<<" time "<<time<<std::endl;
+          // std::cout<<"sector "<<sector<<" row "<<row<<" pad "<<pad<<" time "<<time<<std::endl;
 
           fastTransform->setApplyCorrectionOff();
           float x0, y0, z0;
-          fastTransform->Transform(roc, row, pad, time, x0, y0, z0);
+          fastTransform->Transform(sector, row, pad, time, x0, y0, z0);
 
-          BOOST_CHECK_EQUAL(geo.test(roc, row, y0, z0), 0);
+          BOOST_CHECK_EQUAL(geo.test(sector, row, y0, z0), 0);
 
           fastTransform->setApplyCorrectionOn();
           float x1, y1, z1;
-          fastTransform->Transform(roc, row, pad, time, x1, y1, z1);
+          fastTransform->Transform(sector, row, pad, time, x1, y1, z1);
 
           // local to UV
           float u0, v0, u1, v1;
-          geo.convLocalToUV(roc, y0, z0, u0, v0);
-          geo.convLocalToUV(roc, y1, z1, u1, v1);
+          geo.convLocalToUV(sector, y0, z0, u0, v0);
+          geo.convLocalToUV(sector, y1, z1, u1, v1);
           double dx, du, dv;
-          correctionUV(roc, row, u0, v0, dx, du, dv);
+          correctionUV(sector, row, u0, v0, dx, du, dv);
           statDiff += fabs((x1 - x0) - dx) + fabs((u1 - u0) - du) + fabs((v1 - v0) - dv);
           statN += 3;
           // std::cout << (x1 - x0) - dx << " " << (u1 - u0) - du << " " << (v1 - v0) - dv << std::endl; //": v0 " << v0 <<" z0 "<<z0<<" v1 "<< v1<<" z1 "<<z1 << std::endl;
           // BOOST_CHECK_MESSAGE(0, "SG");
 
           float x1f, y1f, z1f;
-          fromFile->Transform(roc, row, pad, time, x1f, y1f, z1f);
+          fromFile->Transform(sector, row, pad, time, x1f, y1f, z1f);
           statDiffFile += fabs(x1f - x1) + fabs(y1f - y1) + fabs(z1f - z1);
           statNFile += 3;
         }

Detectors/TPC/reconstruction/test/testTPCFastTransform.cxx

@@ -564,7 +564,7 @@ GPUd() int32_t GPUTPCGMMerger::RefitSectorTrack(GPUTPCGMSectorTrack& sectorTrack
   trk.SinPhi() = inTrack->Param().GetSinPhi();
   trk.DzDs() = inTrack->Param().GetDzDs();
   trk.QPt() = inTrack->Param().GetQPt();
-  trk.TOffset() = Param().par.continuousTracking ? GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convZOffsetToVertexTime(sector, inTrack->Param().GetZOffset(), Param().continuousMaxTimeBin) : 0;
+  trk.TOffset() = Param().par.continuousTracking ? GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convZOffsetToVertexTime(inTrack->Param().GetZOffset(), Param().continuousMaxTimeBin) : 0;
   const auto tmp = sectorTrack.ClusterTN() > sectorTrack.ClusterT0() ? std::array<float, 2>{sectorTrack.ClusterTN(), sectorTrack.ClusterT0()} : std::array<float, 2>{sectorTrack.ClusterT0(), sectorTrack.ClusterTN()};
   trk.ShiftZ(this, sector, tmp[0], tmp[1], inTrack->Param().GetX()); // We do not store the inner / outer cluster X, so we just use the track X instead
   sectorTrack.SetX2(0.f);
@@ -1951,7 +1951,7 @@ GPUd() void GPUTPCGMMerger::MergeLoopersInit(int32_t nBlocks, int32_t nThreads,
     const float qptabs = CAMath::Abs(p.GetQPt());
     if (trk.OK() && trk.NClusters() && trk.Leg() == 0 && qptabs * Param().qptB5Scaler > 5.f && qptabs * Param().qptB5Scaler <= lowPtThresh) {
       const int32_t sector = mClusters[trk.FirstClusterRef() + trk.NClusters() - 1].sector;
-      const float refz = p.GetZ() + (Param().par.continuousTracking ? GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, p.GetTOffset(), Param().continuousMaxTimeBin) : 0) + (trk.CSide() ? -100 : 100);
+      const float refz = p.GetZ() + (Param().par.continuousTracking ? GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(p.GetTOffset(), Param().continuousMaxTimeBin) : 0) + (trk.CSide() ? -100 : 100);
       float sinA, cosA;
       CAMath::SinCos(trk.GetAlpha(), sinA, cosA);
       float gx = cosA * p.GetX() - sinA * p.GetY();

GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx

@@ -502,7 +502,7 @@ GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestric
     return -1e6f;
   }
 
-  const float zOffset = param.par.continuousTracking ? Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, mTOffset, param.continuousMaxTimeBin) : 0; // TODO: do some validatiomns for the transform conv functions...
+  const float zOffset = param.par.continuousTracking ? Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(mTOffset, param.continuousMaxTimeBin) : 0; // TODO: do some validatiomns for the transform conv functions...
   const float y0 = row.Grid().YMin();
   const float stepY = row.HstepY();
   const float z0 = row.Grid().ZMin() - zOffset; // We can use our own ZOffset, since this is only used temporarily anyway

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx

@@ -555,7 +555,7 @@ void GPUDisplay::DrawFinal(int32_t iSector, int32_t /*iCol*/, const GPUTPCGMProp
           if constexpr (std::is_same_v<T, GPUTPCGMMergedTrack>) {
             auto cl = mIOPtrs->mergedTrackHits[track->FirstClusterRef() + lastCluster];
             const auto& cln = mIOPtrs->clustersNative->clustersLinear[cl.num];
-            GPUTPCConvertImpl::convert(*mCalib->fastTransform, *mParam, cl.slice, cl.row, cln.getPad(), cln.getTime(), x, y, z);
+            GPUTPCConvertImpl::convert(*mCalib->fastTransform, *mParam, cl.sector, cl.row, cln.getPad(), cln.getTime(), x, y, z);
             ZOffset = mCalib->fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(track->GetParam().GetTOffset(), mParam->continuousMaxTimeBin);
           } else {
             uint8_t sector, row;

GPU/GPUTracking/display/render/GPUDisplayDraw.cxx

@@ -369,7 +369,41 @@ class Spline1DSpec<DataT, YdimT, 0> : public Spline1DContainer<DataT>
     dDr = v * a;
     // F(u) = dSl * Sl + dSr * Sr + dDl * Dl + dDr * Dr;
   }
-
+  /*
+    template <typename T>
+    GPUd() void getUsecondDerivatives(const Knot& knotL, DataT u,
+                                      T& dSl, T& dDl, T& dSr, T& dDr,
+                                      T& dSl2, T& dDl2, T& dSr2, T& dDr2) const
+    {
+      /// Get derivatives of the interpolated value {S(u): 1D -> nYdim} at the segment [knotL, next knotR]
+      /// over the spline values Sl, Sr and the slopes Dl, Dr
+
+      if (u < (DataT)0) {
+        u = (DataT)0;
+      }
+      if (u > (DataT)TBase::getUmax()) {
+        u = (DataT)TBase::getUmax();
+      }
+
+      u = u - knotL.u;
+      T v = u * T(knotL.Li); // scaled u
+      T vm1 = v - 1.;
+      T a = u * vm1;
+      T v2 = v * v;
+      dSr = v2 * (3. - 2 * v);
+      dSl = 1. - dSr;
+      dDl = vm1 * a;
+      dDr = v * a;
+      T dv = T(knotL.Li);
+      dSr2 = 6. * v * (1. - v) * dv;
+      dSl2 = -dSr2;
+      dDl2 = (v - 1) * (3 * v - 1);
+      dDr = u * (v * v - v);
+      dDr2 = 3.f * v * v - 2.f * v;
+      // F(u) = dSl * Sl + dSr * Sr + dDl * Dl + dDr * Dr;
+      // dF(u)/du = dSl2 * Sl + dSr2 * Sr + dDl2 * Dl + dDr2 * Dr;
+    }
+  */
   using TBase::convXtoU;
   using TBase::getKnot;
   using TBase::getKnots;