diff --git a/DataFormats/Reconstruction/src/TrackLTIntegral.cxx b/DataFormats/Reconstruction/src/TrackLTIntegral.cxx index 2a680eea1293f..3f7b39624ab8b 100644 --- a/DataFormats/Reconstruction/src/TrackLTIntegral.cxx +++ b/DataFormats/Reconstruction/src/TrackLTIntegral.cxx @@ -10,12 +10,10 @@ #include "ReconstructionDataFormats/TrackLTIntegral.h" #include "CommonConstants/PhysicsConstants.h" - -#ifndef GPUCA_GPUCODE_DEVICE -#include -#endif +#include "GPUCommonMath.h" using namespace o2::track; +using namespace o2::gpu; //_____________________________________________________ GPUd() void TrackLTIntegral::print() const @@ -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; } } diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt index 5d0b031415a4c..dda6c34f843a1 100644 --- a/Detectors/Base/CMakeLists.txt +++ b/Detectors/Base/CMakeLists.txt @@ -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 @@ -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 diff --git a/Detectors/Base/include/DetectorsBase/Propagator.h b/Detectors/Base/include/DetectorsBase/Propagator.h index c882bbb127b9e..8a85455d251bb 100644 --- a/Detectors/Base/include/DetectorsBase/Propagator.h +++ b/Detectors/Base/include/DetectorsBase/Propagator.h @@ -44,6 +44,11 @@ namespace field class MagFieldFast; } +namespace gpu +{ +class GPUTPCGMPolynomialField; +} + namespace base { class Propagator @@ -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() @@ -123,11 +130,13 @@ class Propagator #endif GPUd() MatBudget getMatBudget(MatCorrType corrType, const o2::math_utils::Point3D& p0, const o2::math_utils::Point3D& p1) const; + GPUd() void getFiedXYZ(const math_utils::Point3D 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); }; diff --git a/Detectors/Base/src/Propagator.cxx b/Detectors/Base/src/Propagator.cxx index 2ad4ebb335d34..2e26615b1c9e2 100644 --- a/Detectors/Base/src/Propagator.cxx +++ b/Detectors/Base/src/Propagator.cxx @@ -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" @@ -118,19 +124,19 @@ GPUd() bool Propagator::PropagateToXBxByBz(o2::track::TrackParCov& track, float } gpu::gpustd::array 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) { @@ -177,19 +183,19 @@ GPUd() bool Propagator::PropagateToXBxByBz(o2::track::TrackPar& track, float xTo } gpu::gpustd::array 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) { @@ -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; } @@ -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) { @@ -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; } @@ -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) { @@ -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(); @@ -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(); @@ -427,27 +433,27 @@ GPUd() bool Propagator::propagateToDCA(const math_utils::Point3D& 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: "; @@ -471,27 +477,27 @@ GPUd() bool Propagator::propagateToDCABxByBz(const math_utils::Point3D& 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: "; @@ -509,10 +515,26 @@ GPUd() bool Propagator::propagateToDCABxByBz(const math_utils::Point3D& v //____________________________________________________________ GPUd() MatBudget Propagator::getMatBudget(Propagator::MatCorrType corrType, const math_utils::Point3D& p0, const math_utils::Point3D& 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 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 + } +} diff --git a/GPU/GPUTracking/Base/GPUReconstructionIncludesDevice.h b/GPU/GPUTracking/Base/GPUReconstructionIncludesDevice.h index ece919daa8a90..6ae8c9119c467 100644 --- a/GPU/GPUTracking/Base/GPUReconstructionIncludesDevice.h +++ b/GPU/GPUTracking/Base/GPUReconstructionIncludesDevice.h @@ -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" diff --git a/GPU/GPUTracking/Merger/GPUTPCGMPolynomialField.h b/GPU/GPUTracking/Merger/GPUTPCGMPolynomialField.h index 674bd8092a789..d334cc7ca2563 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMPolynomialField.h +++ b/GPU/GPUTracking/Merger/GPUTPCGMPolynomialField.h @@ -14,7 +14,7 @@ #ifndef GPUTPCGMPOLYNOMIALFIELD_H #define GPUTPCGMPOLYNOMIALFIELD_H -#include "GPUTPCDef.h" +#include "GPUCommonDef.h" namespace GPUCA_NAMESPACE { @@ -44,13 +44,13 @@ class GPUTPCGMPolynomialField GPUdi() float GetNominalBz() const { return mNominalBz; } - GPUd() void GetField(float x, float y, float z, float B[3]) const; + GPUd() void GetField(float x, float y, float z, float* B) const; GPUd() float GetFieldBz(float x, float y, float z) const; - GPUd() void GetFieldTrd(float x, float y, float z, float B[3]) const; + GPUd() void GetFieldTrd(float x, float y, float z, float* B) const; GPUd() float GetFieldTrdBz(float x, float y, float z) const; - GPUd() void GetFieldIts(float x, float y, float z, float B[3]) const; + GPUd() void GetFieldIts(float x, float y, float z, float* B) const; GPUd() float GetFieldItsBz(float x, float y, float z) const; void Print() const; @@ -164,7 +164,7 @@ GPUdi() void GPUTPCGMPolynomialField::GetPolynomsTpc(float x, float y, float z, f[9] = z * z; } -GPUdi() void GPUTPCGMPolynomialField::GetField(float x, float y, float z, float B[3]) const +GPUdi() void GPUTPCGMPolynomialField::GetField(float x, float y, float z, float* B) const { const float* fBxS = &mTpcBx[1]; const float* fByS = &mTpcBy[1]; @@ -220,7 +220,7 @@ GPUdi() void GPUTPCGMPolynomialField::GetPolynomsTrd(float x, float y, float z, f[19] = z * zz; } -GPUdi() void GPUTPCGMPolynomialField::GetFieldTrd(float x, float y, float z, float B[3]) const +GPUdi() void GPUTPCGMPolynomialField::GetFieldTrd(float x, float y, float z, float* B) const { float f[NTRDM]; GetPolynomsTrd(x, y, z, f); @@ -266,7 +266,7 @@ GPUdi() void GPUTPCGMPolynomialField::GetPolynomsIts(float x, float y, float z, */ } -GPUdi() void GPUTPCGMPolynomialField::GetFieldIts(float x, float y, float z, float B[3]) const +GPUdi() void GPUTPCGMPolynomialField::GetFieldIts(float x, float y, float z, float* B) const { const float* fBxS = &mItsBx[1]; const float* fByS = &mItsBy[1];