Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add possibility to use the GPUPolynomialField in the O2 propagator and make the O2 propagator compile on GPU #5096

Merged
merged 3 commits into from
Dec 25, 2020
Merged

Add possibility to use the GPUPolynomialField in the O2 propagator and make the O2 propagator compile on GPU #5096

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9299733
GPU: Add option to o2 propagator to optionally use GPU PolynomialFiel…
davidrohr Dec 24, 2020
2a7ab61
GPU: Fix the compilation of the o2 propagator on GPU
davidrohr Dec 24, 2020
64fb6f5
GPU: Enable the o2 propagator when compiling for GPU
davidrohr Dec 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
8 changes: 3 additions & 5 deletions DataFormats/Reconstruction/src/TrackLTIntegral.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,10 @@

#include "ReconstructionDataFormats/TrackLTIntegral.h"
#include "CommonConstants/PhysicsConstants.h"

#ifndef GPUCA_GPUCODE_DEVICE
#include <cmath>
#endif
#include "GPUCommonMath.h"

using namespace o2::track;
using namespace o2::gpu;

//_____________________________________________________
GPUd() void TrackLTIntegral::print() const
Expand All @@ -38,7 +36,7 @@ GPUd() void TrackLTIntegral::addStep(float dL, const TrackPar& track)
float dTns = dL * 1000.f / o2::constants::physics::LightSpeedCm2NS; // time change in ps for beta = 1 particle
for (int id = 0; id < getNTOFs(); id++) {
float m2z = o2::track::PID::getMass2Z(id);
float betaInv = std::sqrt(1.f + m2z * m2z * p2);
float betaInv = CAMath::Sqrt(1.f + m2z * m2z * p2);
mT[id] += dTns * betaInv;
}
}
8 changes: 5 additions & 3 deletions Detectors/Base/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ o2_add_library(DetectorsBase
src/MatLayerCyl.cxx
src/MatLayerCylSet.cxx
src/Ray.cxx
src/DCAFitter.cxx
src/DCAFitter.cxx
src/BaseDPLDigitizer.cxx
src/CTFCoderBase.cxx
PUBLIC_LINK_LIBRARIES FairRoot::Base
Expand All @@ -29,8 +29,10 @@ o2_add_library(DetectorsBase
O2::Framework
FairMQ::FairMQ
O2::DataFormatsParameters
O2::SimConfig
ROOT::VMC)
O2::SimConfig
ROOT::VMC
PRIVATE_INCLUDE_DIRECTORIES ${CMAKE_SOURCE_DIR}/GPU/GPUTracking/Merger # Must not link to avoid cyclic dependency
)

o2_target_root_dictionary(DetectorsBase
HEADERS include/DetectorsBase/Detector.h
Expand Down
11 changes: 10 additions & 1 deletion Detectors/Base/include/DetectorsBase/Propagator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,11 @@ namespace field
class MagFieldFast;
}

namespace gpu
{
class GPUTPCGMPolynomialField;
}

namespace base
{
class Propagator
Expand Down Expand Up @@ -104,6 +109,8 @@ class Propagator

GPUd() void setMatLUT(const o2::base::MatLayerCylSet* lut) { mMatLUT = lut; }
GPUd() const o2::base::MatLayerCylSet* getMatLUT() const { return mMatLUT; }
GPUd() void setGPUField(const o2::gpu::GPUTPCGMPolynomialField* field) { mGPUField = field; }
GPUd() const o2::gpu::GPUTPCGMPolynomialField* getGPUField() const { return mGPUField; }

#ifndef GPUCA_GPUCODE
static Propagator* Instance()
Expand All @@ -123,11 +130,13 @@ class Propagator
#endif

GPUd() MatBudget getMatBudget(MatCorrType corrType, const o2::math_utils::Point3D<float>& p0, const o2::math_utils::Point3D<float>& p1) const;
GPUd() void getFiedXYZ(const math_utils::Point3D<float> xyz, float* bxyz) const;

const o2::field::MagFieldFast* mField = nullptr; ///< External fast field (barrel only for the moment)
float mBz = 0; // nominal field

const o2::base::MatLayerCylSet* mMatLUT = nullptr; // externally set LUT
const o2::base::MatLayerCylSet* mMatLUT = nullptr; // externally set LUT
const o2::gpu::GPUTPCGMPolynomialField* mGPUField = nullptr; // externally set GPU Field

ClassDef(Propagator, 0);
};
Expand Down
104 changes: 63 additions & 41 deletions Detectors/Base/src/Propagator.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,13 +10,19 @@

#include "DetectorsBase/Propagator.h"
#include "GPUCommonLogger.h"
#include "Field/MagFieldFast.h"
#include "GPUCommonMath.h"
#include "GPUTPCGMPolynomialField.h"
#include "MathUtils/Utils.h"
#include "ReconstructionDataFormats/Vertex.h"

using namespace o2::base;
using namespace o2::gpu;

#ifndef GPUCA_STANDALONE
#if !defined(GPUCA_GPUCODE)
#include "Field/MagFieldFast.h" // Don't use this on the GPU
#endif

#if !defined(GPUCA_STANDALONE) && !defined(GPUCA_GPUCODE)
#include "Field/MagneticField.h"
#include "DataFormatsParameters/GRPObject.h"
#include "DetectorsBase/GeometryManager.h"
Expand Down Expand Up @@ -118,19 +124,19 @@ GPUd() bool Propagator::PropagateToXBxByBz(o2::track::TrackParCov& track, float
}

gpu::gpustd::array<float, 3> b;
while (std::abs(dx) > Epsilon) {
auto step = std::min(std::abs(dx), maxStep);
while (CAMath::Abs(dx) > Epsilon) {
auto step = CAMath::Min(CAMath::Abs(dx), maxStep);
if (dir < 0) {
step = -step;
}
auto x = track.getX() + step;
auto xyz0 = track.getXYZGlo();
mField->Field(xyz0, &b[0]);
getFiedXYZ(xyz0, &b[0]);

if (!track.propagateTo(x, b)) {
return false;
}
if (maxSnp > 0 && std::abs(track.getSnp()) >= maxSnp) {
if (maxSnp > 0 && CAMath::Abs(track.getSnp()) >= maxSnp) {
return false;
}
if (matCorr != MatCorrType::USEMatCorrNONE) {
Expand Down Expand Up @@ -177,19 +183,19 @@ GPUd() bool Propagator::PropagateToXBxByBz(o2::track::TrackPar& track, float xTo
}

gpu::gpustd::array<float, 3> b;
while (std::abs(dx) > Epsilon) {
auto step = std::min(std::abs(dx), maxStep);
while (CAMath::Abs(dx) > Epsilon) {
auto step = CAMath::Min(CAMath::Abs(dx), maxStep);
if (dir < 0) {
step = -step;
}
auto x = track.getX() + step;
auto xyz0 = track.getXYZGlo();
mField->Field(xyz0, &b[0]);
getFiedXYZ(xyz0, &b[0]);

if (!track.propagateParamTo(x, b)) {
return false;
}
if (maxSnp > 0 && std::abs(track.getSnp()) >= maxSnp) {
if (maxSnp > 0 && CAMath::Abs(track.getSnp()) >= maxSnp) {
return false;
}
if (matCorr != MatCorrType::USEMatCorrNONE) {
Expand Down Expand Up @@ -234,8 +240,8 @@ GPUd() bool Propagator::propagateToX(o2::track::TrackParCov& track, float xToGo,
signCorr = -dir; // sign of eloss correction is not imposed
}

while (std::abs(dx) > Epsilon) {
auto step = std::min(std::abs(dx), maxStep);
while (CAMath::Abs(dx) > Epsilon) {
auto step = CAMath::Min(CAMath::Abs(dx), maxStep);
if (dir < 0) {
step = -step;
}
Expand All @@ -245,7 +251,7 @@ GPUd() bool Propagator::propagateToX(o2::track::TrackParCov& track, float xToGo,
if (!track.propagateTo(x, bZ)) {
return false;
}
if (maxSnp > 0 && std::abs(track.getSnp()) >= maxSnp) {
if (maxSnp > 0 && CAMath::Abs(track.getSnp()) >= maxSnp) {
return false;
}
if (matCorr != MatCorrType::USEMatCorrNONE) {
Expand Down Expand Up @@ -292,8 +298,8 @@ GPUd() bool Propagator::propagateToX(o2::track::TrackPar& track, float xToGo, fl
signCorr = -dir; // sign of eloss correction is not imposed
}

while (std::abs(dx) > Epsilon) {
auto step = std::min(std::abs(dx), maxStep);
while (CAMath::Abs(dx) > Epsilon) {
auto step = CAMath::Min(CAMath::Abs(dx), maxStep);
if (dir < 0) {
step = -step;
}
Expand All @@ -303,7 +309,7 @@ GPUd() bool Propagator::propagateToX(o2::track::TrackPar& track, float xToGo, fl
if (!track.propagateParamTo(x, bZ)) {
return false;
}
if (maxSnp > 0 && std::abs(track.getSnp()) >= maxSnp) {
if (maxSnp > 0 && CAMath::Abs(track.getSnp()) >= maxSnp) {
return false;
}
if (matCorr != MatCorrType::USEMatCorrNONE) {
Expand Down Expand Up @@ -337,27 +343,27 @@ GPUd() bool Propagator::propagateToDCA(const o2::dataformats::VertexBase& vtx, o
// propagate track to DCA to the vertex
float sn, cs, alp = track.getAlpha();
o2::math_utils::sincos(alp, sn, cs);
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = std::sqrt((1.f - snp) * (1.f + snp));
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = CAMath::Sqrt((1.f - snp) * (1.f + snp));
float xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
x -= xv;
y -= yv;
//Estimate the impact parameter neglecting the track curvature
float d = std::abs(x * snp - y * csp);
float d = CAMath::Abs(x * snp - y * csp);
if (d > maxD) {
return false;
}
float crv = track.getCurvature(bZ);
float tgfv = -(crv * x - snp) / (crv * y + csp);
sn = tgfv / std::sqrt(1.f + tgfv * tgfv);
cs = std::sqrt((1. - sn) * (1. + sn));
cs = (std::abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
sn = tgfv / CAMath::Sqrt(1.f + tgfv * tgfv);
cs = CAMath::Sqrt((1. - sn) * (1. + sn));
cs = (CAMath::Abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;

x = xv * cs + yv * sn;
yv = -xv * sn + yv * cs;
xv = x;

auto tmpT(track); // operate on the copy to recover after the failure
alp += std::asin(sn);
alp += CAMath::ASin(sn);
if (!tmpT.rotate(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
LOG(WARNING) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx << " | Track is: ";
tmpT.print();
Expand All @@ -382,27 +388,27 @@ GPUd() bool Propagator::propagateToDCABxByBz(const o2::dataformats::VertexBase&
// propagate track to DCA to the vertex
float sn, cs, alp = track.getAlpha();
o2::math_utils::sincos(alp, sn, cs);
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = std::sqrt((1.f - snp) * (1.f + snp));
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = CAMath::Sqrt((1.f - snp) * (1.f + snp));
float xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
x -= xv;
y -= yv;
//Estimate the impact parameter neglecting the track curvature
float d = std::abs(x * snp - y * csp);
float d = CAMath::Abs(x * snp - y * csp);
if (d > maxD) {
return false;
}
float crv = track.getCurvature(mBz);
float tgfv = -(crv * x - snp) / (crv * y + csp);
sn = tgfv / std::sqrt(1.f + tgfv * tgfv);
cs = std::sqrt((1. - sn) * (1. + sn));
cs = (std::abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
sn = tgfv / CAMath::Sqrt(1.f + tgfv * tgfv);
cs = CAMath::Sqrt((1. - sn) * (1. + sn));
cs = (CAMath::Abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;

x = xv * cs + yv * sn;
yv = -xv * sn + yv * cs;
xv = x;

auto tmpT(track); // operate on the copy to recover after the failure
alp += std::asin(sn);
alp += CAMath::ASin(sn);
if (!tmpT.rotate(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
LOG(WARNING) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx << " | Track is: ";
tmpT.print();
Expand All @@ -427,27 +433,27 @@ GPUd() bool Propagator::propagateToDCA(const math_utils::Point3D<float>& vtx, o2
// propagate track to DCA to the vertex
float sn, cs, alp = track.getAlpha();
o2::math_utils::sincos(alp, sn, cs);
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = std::sqrt((1.f - snp) * (1.f + snp));
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = CAMath::Sqrt((1.f - snp) * (1.f + snp));
float xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
x -= xv;
y -= yv;
//Estimate the impact parameter neglecting the track curvature
float d = std::abs(x * snp - y * csp);
float d = CAMath::Abs(x * snp - y * csp);
if (d > maxD) {
return false;
}
float crv = track.getCurvature(bZ);
float tgfv = -(crv * x - snp) / (crv * y + csp);
sn = tgfv / std::sqrt(1.f + tgfv * tgfv);
cs = std::sqrt((1. - sn) * (1. + sn));
cs = (std::abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
sn = tgfv / CAMath::Sqrt(1.f + tgfv * tgfv);
cs = CAMath::Sqrt((1. - sn) * (1. + sn));
cs = (CAMath::Abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;

x = xv * cs + yv * sn;
yv = -xv * sn + yv * cs;
xv = x;

auto tmpT(track); // operate on the copy to recover after the failure
alp += std::asin(sn);
alp += CAMath::ASin(sn);
if (!tmpT.rotateParam(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
LOG(WARNING) << "failed to propagate to alpha=" << alp << " X=" << xv << " for vertex "
<< vtx.X() << ' ' << vtx.Y() << ' ' << vtx.Z() << " | Track is: ";
Expand All @@ -471,27 +477,27 @@ GPUd() bool Propagator::propagateToDCABxByBz(const math_utils::Point3D<float>& v
// propagate track to DCA to the vertex
float sn, cs, alp = track.getAlpha();
o2::math_utils::sincos(alp, sn, cs);
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = std::sqrt((1.f - snp) * (1.f + snp));
float x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = CAMath::Sqrt((1.f - snp) * (1.f + snp));
float xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
x -= xv;
y -= yv;
//Estimate the impact parameter neglecting the track curvature
float d = std::abs(x * snp - y * csp);
float d = CAMath::Abs(x * snp - y * csp);
if (d > maxD) {
return false;
}
float crv = track.getCurvature(mBz);
float tgfv = -(crv * x - snp) / (crv * y + csp);
sn = tgfv / std::sqrt(1.f + tgfv * tgfv);
cs = std::sqrt((1. - sn) * (1. + sn));
cs = (std::abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
sn = tgfv / CAMath::Sqrt(1.f + tgfv * tgfv);
cs = CAMath::Sqrt((1. - sn) * (1. + sn));
cs = (CAMath::Abs(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;

x = xv * cs + yv * sn;
yv = -xv * sn + yv * cs;
xv = x;

auto tmpT(track); // operate on the copy to recover after the failure
alp += std::asin(sn);
alp += CAMath::ASin(sn);
if (!tmpT.rotateParam(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
LOG(WARNING) << "failed to propagate to alpha=" << alp << " X=" << xv << " for vertex "
<< vtx.X() << ' ' << vtx.Y() << ' ' << vtx.Z() << " | Track is: ";
Expand All @@ -509,10 +515,26 @@ GPUd() bool Propagator::propagateToDCABxByBz(const math_utils::Point3D<float>& v
//____________________________________________________________
GPUd() MatBudget Propagator::getMatBudget(Propagator::MatCorrType corrType, const math_utils::Point3D<float>& p0, const math_utils::Point3D<float>& p1) const
{
#ifndef GPUCA_STANDALONE
#if !defined(GPUCA_STANDALONE) && !defined(GPUCA_GPUCODE)
if (corrType == MatCorrType::USEMatCorrTGeo) {
return GeometryManager::meanMaterialBudget(p0, p1);
}
#endif
return mMatLUT->getMatBudget(p0.X(), p0.Y(), p0.Z(), p1.X(), p1.Y(), p1.Z());
}

GPUd() void Propagator::getFiedXYZ(const math_utils::Point3D<float> xyz, float* bxyz) const
{
if (mGPUField) {
#if defined(GPUCA_GPUCODE_DEVICE) && defined(GPUCA_HAS_GLOBAL_SYMBOL_CONSTANT_MEM)
const auto* f = &GPUCA_CONSMEM.param.polynomialField; // Access directly from constant memory on GPU (copied here to avoid complicated header dependencies)
#else
const auto* f = mGPUField;
#endif
f->GetField(xyz.X(), xyz.Y(), xyz.Z(), bxyz);
} else {
#ifndef GPUCA_GPUCODE
mField->Field(xyz, bxyz); // Must not call the host-only function in GPU compilation
#endif
}
}
2 changes: 2 additions & 0 deletions GPU/GPUTracking/Base/GPUReconstructionIncludesDevice.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,8 @@ using namespace GPUCA_NAMESPACE::gpu;
// O2 track model
#include "TrackParametrization.cxx"
#include "TrackParametrizationWithError.cxx"
#include "Propagator.cxx"
#include "TrackLTIntegral.cxx"

// Files for GPU dEdx
#include "GPUdEdx.cxx"
Expand Down