ClusterLines.h

Go to the documentation of this file.

// 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.
 
#ifndef O2_ITS_CLUSTERLINES_H
#define O2_ITS_CLUSTERLINES_H
 
#include <gsl/span>
#include <vector>
#include <array>
#ifndef GPUCA_GPUCODE
#include <cmath>
#else
#include "GPUCommonAlgorithm.h" // warp reductions of WarpFit
#endif
#include "ITStracking/Cluster.h"
#include "ITStracking/Cell.h"
#include "ITSMFTTracking/Constants.h"
#include "ITStracking/Tracklet.h"
#include "GPUCommonDef.h"
#include "GPUCommonMath.h"
 
namespace o2::its
{
 
struct Line final {
  GPUhdDefault() Line() = default;
  Line(const Tracklet&, const Cluster*, const Cluster*);
  GPUhdi() Line(const float origin[3], const float direction[3], const TimeEstBC& time) : mTime(time)
  {
    float norm2 = 0.f;
    for (int i = 0; i < 3; ++i) {
      norm2 += direction[i] * direction[i];
    }
    const float norm = o2::gpu::GPUCommonMath::Sqrt(norm2);
    for (int i = 0; i < 3; ++i) {
      originPoint[i] = origin[i];
      cosinesDirector[i] = norm2 > 0.f ? direction[i] / norm : 0.f;
    }
  }
  Line(const std::array<float, 3>& origin, const std::array<float, 3>& direction, const TimeEstBC& time)
    : Line(origin.data(), direction.data(), time) {}
  bool operator==(const Line&) const = default;
 
  GPUhdi() static float getDistance2FromPoint(const Line& line, const float point[3])
  {
    float delta[3], proj = 0.f;
    for (int i = 0; i < 3; ++i) {
      delta[i] = point[i] - line.originPoint[i];
      proj += delta[i] * line.cosinesDirector[i];
    }
    float d2 = 0.f;
    for (int i = 0; i < 3; ++i) {
      const float residual = delta[i] - proj * line.cosinesDirector[i];
      d2 += residual * residual;
    }
    return d2;
  }
 
  GPUhdi() static void getDCAComponents(const Line& line, const float point[3], float out[6])
  {
    float delta[3], proj = 0.f;
    for (int i = 0; i < 3; ++i) {
      delta[i] = line.originPoint[i] - point[i];
      proj += delta[i] * line.cosinesDirector[i];
    }
    float r[3];
    for (int i = 0; i < 3; ++i) {
      r[i] = delta[i] - proj * line.cosinesDirector[i];
    }
    out[0] = r[0];
    out[2] = r[1];
    out[5] = r[2];
#ifndef GPUCA_GPUCODE
    out[1] = std::hypot(r[0], r[1]);
    out[3] = std::hypot(r[0], r[2]);
    out[4] = std::hypot(r[1], r[2]);
#else
    out[1] = o2::gpu::GPUCommonMath::Hypot(r[0], r[1]);
    out[3] = o2::gpu::GPUCommonMath::Hypot(r[0], r[2]);
    out[4] = o2::gpu::GPUCommonMath::Hypot(r[1], r[2]);
#endif
  }
 
  GPUhdi() static bool getZAtBeam(const Line& line, const float beamX, const float beamY, float& z)
  {
    const float ux = line.cosinesDirector[0];
    const float uy = line.cosinesDirector[1];
    const float den = ux * ux + uy * uy;
    if (den <= constants::Tolerance) {
      return false;
    }
    const float dx = line.originPoint[0] - beamX;
    const float dy = line.originPoint[1] - beamY;
    const float s0 = -(dx * ux + dy * uy) / den;
    z = line.originPoint[2] + s0 * line.cosinesDirector[2];
    return true;
  }
 
  static float getDistance2FromPoint(const Line& line, const std::array<float, 3>& point)
  {
    return getDistance2FromPoint(line, point.data());
  }
  static float getDistanceFromPoint(const Line& line, const std::array<float, 3>& point);
  static float getDCA2(const Line&, const Line&, const float precision = constants::Tolerance);
  static float getDCA(const Line&, const Line&, const float precision = constants::Tolerance);
  GPUhdi() bool isEmpty() const noexcept
  {
    return originPoint[0] == 0.f && originPoint[1] == 0.f && originPoint[2] == 0.f &&
           cosinesDirector[0] == 0.f && cosinesDirector[1] == 0.f && cosinesDirector[2] == 0.f;
  }
  void print() const;
 
  float originPoint[3] = {0.f, 0.f, 0.f};
  float cosinesDirector[3] = {0.f, 0.f, 0.f};
  TimeEstBC mTime;
};
 
GPUdi() bool getCellLineAtBeam(const CellSeed& cell, const float beamX, const float beamY, const float bz, const float curvatureScale, std::array<float, 3>& origin, std::array<float, 3>& direction)
{
  cell.getXYZGlo(origin);
  if (!cell.getPxPyPzGlo(direction)) {
    return false;
  }
  o2::track::TrackParametrization<float> par{cell};
  par.setQ2Pt(par.getQ2Pt() * curvatureScale);
  std::array<float, 3> dcaOrigin, dcaDirection;
  if (par.propagateParamToDCA({beamX, beamY, 0.f}, bz)) {
    par.getXYZGlo(dcaOrigin);
    if (par.getPxPyPzGlo(dcaDirection)) {
      origin = dcaOrigin;
      direction = dcaDirection;
    }
  }
  return true;
}
 
class LineFit
{
 public:
  GPUhdi() bool add(const Line& line)
  {
    const double d0 = line.cosinesDirector[0], d1 = line.cosinesDirector[1], d2 = line.cosinesDirector[2];
    const double o0 = line.originPoint[0], o1 = line.originPoint[1], o2 = line.originPoint[2];
    const double det = d0 * d0 + d1 * d1 + d2 * d2; // == 1 for a normalised director
    if (det <= 0.) {
      return false;
    }
    mA[0] += (det - d0 * d0) / det;
    mA[1] += (-d0 * d1) / det;
    mA[2] += (-d0 * d2) / det;
    mA[3] += (det - d1 * d1) / det;
    mA[4] += (-d1 * d2) / det;
    mA[5] += (det - d2 * d2) / det;
    const double dDotO = d0 * o0 + d1 * o1 + d2 * o2;
    mB[0] += (d0 * dDotO - det * o0) / det;
    mB[1] += (d1 * dDotO - det * o1) / det;
    mB[2] += (d2 * dDotO - det * o2) / det;
    return true;
  }
 
  GPUhdi() bool solve(float vertex[3]) const
  {
    const double c00 = mA[3] * mA[5] - mA[4] * mA[4];
    const double c01 = mA[2] * mA[4] - mA[1] * mA[5];
    const double c02 = mA[1] * mA[4] - mA[2] * mA[3];
    const double c11 = mA[0] * mA[5] - mA[2] * mA[2];
    const double c12 = mA[1] * mA[2] - mA[0] * mA[4];
    const double c22 = mA[0] * mA[3] - mA[1] * mA[1];
    const double det = mA[0] * c00 + mA[1] * c01 + mA[2] * c02;
    if (o2::gpu::GPUCommonMath::Abs(det) < 1.e-12) {
      return false;
    }
    const double invDet = 1. / det;
    const double x0 = (c00 * mB[0] + c01 * mB[1] + c02 * mB[2]) * invDet;
    const double x1 = (c01 * mB[0] + c11 * mB[1] + c12 * mB[2]) * invDet;
    const double x2 = (c02 * mB[0] + c12 * mB[1] + c22 * mB[2]) * invDet;
    vertex[0] = static_cast<float>(-x0);
    vertex[1] = static_cast<float>(-x1);
    vertex[2] = static_cast<float>(-x2);
    return true;
  }
 
  GPUhdi() void computeClusterCentroid() { mIsValid = solve(mVertex); }
 
  GPUhdi() bool isValid() const noexcept { return mIsValid; }
  GPUhdi() const float* getVertex() const noexcept { return mVertex; }
  GPUhdi() int getSize() const noexcept { return mNContributors; }
  GPUhdi() int getNContributors() const noexcept { return mNContributors; }
  GPUhdi() const TimeEstBC& getTimeStamp() const noexcept { return mTime; }
  GPUhdi() float getR2() const noexcept { return (mVertex[0] * mVertex[0]) + (mVertex[1] * mVertex[1]); }
  GPUhdi() float getR() const noexcept { return o2::gpu::GPUCommonMath::Sqrt(getR2()); }
  GPUhdi() const float* getRMS2() const noexcept { return mRMS2; } // {XX, XY, YY, XZ, YZ, ZZ}
  GPUhdi() float getAvgDistance2() const noexcept { return mAvgDistance2; }
 
 protected:
  double mA[6] = {0., 0., 0., 0., 0., 0.};         // AX=B, packed symmetric normal matrix
  double mB[3] = {0., 0., 0.};                     // AX=B, right-hand side
  float mVertex[3] = {0.f, 0.f, 0.f};              // cluster centroid position
  float mRMS2[6] = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f}; // mean DCA components about the vertex
  float mAvgDistance2 = 0.f;                       // mean squared distance about the vertex
  int mNContributors = 0;                          // lines accumulated into the normal equations
  bool mIsValid = false;                           // true if the linear system was solved successfully
  TimeEstBC mTime;                                 // intersection of the contributors' time intervals
};
 
class ClusterLines final : public LineFit
{
 public:
  GPUhdDefault() ClusterLines() = default;
  ClusterLines(gsl::span<const int> lineIndices, gsl::span<const Line> lines);
 
  GPUhdi() void add(const Line& line)
  {
    if (!LineFit::add(line)) {
      return;
    }
    if (mNContributors <= 0) {
      mTime = line.mTime;
    } else {
      mTime += line.mTime;
    }
    ++mNContributors;
  }
 
  GPUhdi() void addResidual(const Line& line) { addResidual(line, mVertex); }
 
  GPUhdi() void addResidual(const Line& line, const float vertex[3])
  {
    float dca[6];
    Line::getDCAComponents(line, vertex, dca);
    const float d2 = Line::getDistance2FromPoint(line, vertex);
    ++mResidualCount;
    const float inv = 1.f / static_cast<float>(mResidualCount);
    for (int i = 0; i < 6; ++i) {
      mRMS2[i] += (dca[i] - mRMS2[i]) * inv;
    }
    mAvgDistance2 += (d2 - mAvgDistance2) / static_cast<float>(mResidualCount);
  }
 
  bool operator==(const ClusterLines& rhs) const noexcept;
 
 private:
  int mResidualCount = 0; // lines folded into the running means
};
 
#ifdef GPUCA_GPUCODE
struct WarpFit : public LineFit {
  uint32_t tLo = 0, tHi = 0xffffffffu;
 
  GPUdi() void add(const Line& l)
  {
    if (LineFit::add(l)) {
      tLo = o2::gpu::GPUCommonMath::Max(tLo, l.mTime.lower());
      tHi = o2::gpu::GPUCommonMath::Min(tHi, l.mTime.upper());
      ++mNContributors;
    }
  }
 
  template <int Width>
  GPUdi() void reduce()
  {
    for (int i = 0; i < 6; ++i) {
      mA[i] = warp_reduce_add<Width>(mA[i]);
    }
    for (int i = 0; i < 3; ++i) {
      mB[i] = warp_reduce_add<Width>(mB[i]);
    }
    tLo = warp_reduce_max<Width>(tLo);
    tHi = warp_reduce_min<Width>(tHi);
    mNContributors = warp_reduce_add<Width>(mNContributors);
    mTime = TimeEstBC{tLo, tHi - tLo};
  }
 
  GPUdi() void addResidual(const Line& l, const float vertex[3])
  {
    float dca[6];
    Line::getDCAComponents(l, vertex, dca);
    for (int i = 0; i < 6; ++i) {
      rms2Sum[i] += dca[i];
    }
    dist2Sum += Line::getDistance2FromPoint(l, vertex);
  }
 
  template <int Width>
  GPUdi() void reduceResiduals(const int n)
  {
    for (int i = 0; i < 6; ++i) {
      mRMS2[i] = warp_reduce_add<Width>(rms2Sum[i]) / n;
    }
    mAvgDistance2 = warp_reduce_add<Width>(dist2Sum) / n;
  }
 
  double rms2Sum[6] = {0., 0., 0., 0., 0., 0.};
  double dist2Sum = 0.;
};
#endif
 
} // namespace o2::its
#endif /* O2_ITS_CLUSTERLINES_H */