Skip to content

Commit

Permalink
Merge branch 'ACTS-385_Refine_BField_Handling' into 'master'
Browse files Browse the repository at this point in the history 
BField caching and handling Updates

Closes ACTS-385 and ACTS-379

See merge request !353
  • Loading branch information
@asalzburger
asalzburger committed Oct 23, 2017
2 parents d668c25 + 0f30822 commit 95dd052
Show file tree
Hide file tree
Showing 8 changed files with 321 additions and 126 deletions.
80 changes: 38 additions & 42 deletions Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -500,7 +500,8 @@ Acts::RungeKuttaEngine<MagneticField>::propagateWithJacobian(
SA[0] = SA[1] = SA[2] = 0.;
pCache.maxPathLimit = false;

if (pCache.mcondition && std::abs(pCache.pVector[6]) > 100.) return false;
// @todo the inverse momentum condition is hard-coded
if (pCache.mcondition && std::abs(pCache.pVector[6]) > 100000.) return false;

// Step estimation until surface
bool Q;
Expand Down Expand Up @@ -646,20 +647,18 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStep(int navigationStep,
double* R = &(pCache.pVector[0]); // Coordinates
double* A = &(pCache.pVector[3]); // Directions
double* sA = &(pCache.pVector[42]);

// @todo bring the numbers into initialize
double scaleM = 1e-3 / units::_T;
double Pi = 149.89626 * pCache.pVector[6] * units::_MeV; // Invert mometum/2.
// invert momentum in right units
double Pi = 0.5 / units::Nat2SI<units::MOMENTUM>(1. / pCache.pVector[6]);
double dltm = m_cfg.dlt * .03;

double f0[3], f[3];

// if new field is required get it
if (pCache.newfield) {
Vector3D bField = m_cfg.fieldService->getField(Vector3D(R[0], R[1], R[2]));
f0[0] = scaleM * bField.x();
f0[1] = scaleM * bField.y();
f0[2] = scaleM * bField.z();
f0[0] = bField.x();
f0[1] = bField.y();
f0[2] = bField.z();
} else {
f0[0] = pCache.field[0];
f0[1] = pCache.field[1];
Expand Down Expand Up @@ -690,9 +689,9 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStep(int navigationStep,
double gP[3] = {R[0] + A1 * S4, R[1] + B1 * S4, R[2] + C1 * S4};
Vector3D bField
= m_cfg.fieldService->getField(Vector3D(gP[0], gP[1], gP[2]));
f[0] = scaleM * bField.x();
f[1] = scaleM * bField.y();
f[2] = scaleM * bField.z();
f[0] = bField.x();
f[1] = bField.y();
f[2] = bField.z();
} else {
f[0] = f0[0];
f[1] = f0[1];
Expand All @@ -716,9 +715,9 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStep(int navigationStep,
double gP[3] = {R[0] + S * A4, R[1] + S * B4, R[2] + S * C4};
Vector3D bField
= m_cfg.fieldService->getField(Vector3D(gP[0], gP[1], gP[2]));
f[0] = scaleM * bField.x();
f[1] = scaleM * bField.y();
f[2] = scaleM * bField.z();
f[0] = bField.x();
f[1] = bField.y();
f[2] = bField.z();
} else {
f[0] = f0[0];
f[1] = f0[1];
Expand Down Expand Up @@ -890,10 +889,7 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStepWithGradient(
double* A = &(pCache.pVector[3]); // Directions
double* sA = &(pCache.pVector[42]);
// express the conversion factor with units
// @todo bring the conversion into initialize
double scaleM = 1e-3 / units::_T;
double Pi = 149.89626 * pCache.pVector[6] * units::_MeV; // Invert mometum/2.

double Pi = 0.5 / units::Nat2SI<units::MOMENTUM>(1. / pCache.pVector[6]);
double dltm = m_cfg.dlt * .03;

double f0[3], f1[3], f2[3], g0[9], g1[9], g2[9], H0[12], H1[12], H2[12];
Expand All @@ -902,18 +898,18 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStepWithGradient(
deriv0 << 0., 0., 0., 0., 0., 0., 0., 0., 0.;
Vector3D bField0 = m_cfg.fieldService->getFieldGradient(
Vector3D(R[0], R[1], R[2]), deriv0);
f0[0] = scaleM * bField0.x();
f0[1] = scaleM * bField0.y();
f0[2] = scaleM * bField0.z();
g0[0] = scaleM * deriv0(0, 0);
g0[1] = scaleM * deriv0(0, 1);
g0[2] = scaleM * deriv0(0, 2);
g0[3] = scaleM * deriv0(1, 0);
g0[4] = scaleM * deriv0(1, 1);
g0[5] = scaleM * deriv0(1, 2);
g0[6] = scaleM * deriv0(2, 0);
g0[7] = scaleM * deriv0(2, 1);
g0[8] = scaleM * deriv0(2, 2);
f0[0] = bField0.x();
f0[1] = bField0.y();
f0[2] = bField0.z();
g0[0] = deriv0(0, 0);
g0[1] = deriv0(0, 1);
g0[2] = deriv0(0, 2);
g0[3] = deriv0(1, 0);
g0[4] = deriv0(1, 1);
g0[5] = deriv0(1, 2);
g0[6] = deriv0(2, 0);
g0[7] = deriv0(2, 1);
g0[8] = deriv0(2, 2);

while (S != 0.) {
double S3 = C33 * S, S4 = .25 * S, PS2 = Pi * S;
Expand All @@ -940,18 +936,18 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStepWithGradient(
deriv1 << 0., 0., 0., 0., 0., 0., 0., 0., 0.;
Vector3D bField1 = m_cfg.fieldService->getFieldGradient(
Vector3D(gP1[0], gP1[1], gP1[2]), deriv1);
f1[0] = scaleM * bField1.x();
f1[1] = scaleM * bField1.y();
f1[2] = scaleM * bField1.z();
g1[0] = scaleM * deriv1(0, 0);
g1[1] = scaleM * deriv1(0, 1);
g1[2] = scaleM * deriv1(0, 2);
g1[3] = scaleM * deriv1(1, 0);
g1[4] = scaleM * deriv1(1, 1);
g1[5] = scaleM * deriv1(1, 2);
g1[6] = scaleM * deriv1(2, 0);
g1[7] = scaleM * deriv1(2, 1);
g1[8] = scaleM * deriv1(2, 2);
f1[0] = bField1.x();
f1[1] = bField1.y();
f1[2] = bField1.z();
g1[0] = deriv1(0, 0);
g1[1] = deriv1(0, 1);
g1[2] = deriv1(0, 2);
g1[3] = deriv1(1, 0);
g1[4] = deriv1(1, 1);
g1[5] = deriv1(1, 2);
g1[6] = deriv1(2, 0);
g1[7] = deriv1(2, 1);
g1[8] = deriv1(2, 2);
H1[0] = f1[0] * PS2;
H1[1] = f1[1] * PS2;
H1[2] = f1[2] * PS2;
Expand Down
76 changes: 76 additions & 0 deletions Core/include/ACTS/MagneticField/SharedBFieldMap.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
// This file is part of the ACTS project.
//
// Copyright (C) 2016 ACTS project team
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef ACTS_MAGNETICFIELD_SHAREDBFIELD_H
#define ACTS_MAGNETICFIELD_SHAREDBFIELD_H

#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp"
#include "ACTS/Utilities/Definitions.hpp"

namespace Acts {

/// @ingroup MagneticField
///
/// @brief allows to use a shared magnetic field
/// in several places and with multiple steppers
/// mainly targeted to save memory
template <typename BField>
class SharedBField
{
public:
/// @brief the constructur with a shared pointer
/// @note since it is a shared field, we enforce it to be const
/// @tparam bField is the shared BField to be stored
SharedBField(std::shared_ptr<const BField> bField) : m_bField(bField) {}

/// @brief retrieve magnetic field value
///
/// @param [in] position global 3D position
///
/// @return magnetic field vector at given position
Vector3D
getField(const Vector3D& position) const
{
return m_bField->getField(position);
}

/// @brief retrieve field cell for given position
///
/// @param [in] position global 3D position
/// @return field cell containing the given global position
///
/// @pre The given @c position must lie within the range of the underlying
/// magnetic field map.
concept::AnyFieldCell<>
getFieldCell(const Vector3D& position) const
{
return m_bField->getFieldCell(position);
}

/// @brief retrieve magnetic field value & its gradient
///
/// @param [in] position global 3D position
/// @param [out] derivative gradient of magnetic field vector as (3x3) matrix
/// @return magnetic field vector
///
/// @note currently the derivative is not calculated
/// @todo return derivative
Vector3D
getFieldGradient(const Vector3D& position,
ActsMatrixD<3, 3>& derivative) const
{
return m_bField->getField(position);
}

private:
std::shared_ptr<const BField> m_bField;
};

} // namespace Acts

#endif // ACTS_MAGNETICFIELD_SHAREDBFIELD_H
35 changes: 32 additions & 3 deletions Core/include/ACTS/Propagator/AtlasStepper.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@

#include <cmath>
#include "ACTS/EventData/TrackParameters.hpp"
#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp"
#include "ACTS/Surfaces/Surface.hpp"
#include "ACTS/Utilities/Units.hpp"

Expand All @@ -37,6 +38,12 @@ class AtlasStepper
const ActsSymMatrixD<NGlobalPars>* covariance;
double jacobian[NGlobalPars * NGlobalPars];

/// Lazily initialized cache
/// It caches the current magneticl field cell and stays interpolates within
/// as long as this is valid. See step() code for details.
bool field_cache_ready = false;
concept::AnyFieldCell<> field_cache;

Vector3D
position() const
{
Expand Down Expand Up @@ -287,6 +294,9 @@ class AtlasStepper
return CurvilinearParameters(std::move(cov), gp, mom, charge);
}

/// convert method into bound parameters
/// @param cache is the propagation cache to be converted
/// @param s the target surface
static BoundParameters
convert(Cache& cache, const Surface& s)
{
Expand Down Expand Up @@ -446,6 +456,25 @@ class AtlasStepper
return i.pathLength;
}

/// Get the field for the stepping
/// It checks first if the access is still within the Cell,
/// and updates the cell if necessary, then it takes the field
/// from the cell
/// @param [in,out] cache is the propagation cache associated with the track
/// the magnetic field cell is used (and potentially updated)
/// @param [in] pos is the field position
Vector3D
getField(Cache& cache, const Vector3D& pos) const
{
if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)) {
cache.field_cache_ready = true;
cache.field_cache = m_bField.getFieldCell(pos);
}
// get the field from the cell
cache.field = cache.field_cache.getField(pos);
return cache.field;
}

double
step(Cache& cache, double& h) const
{
Expand All @@ -463,7 +492,7 @@ class AtlasStepper
// if new field is required get it
if (cache.newfield) {
const Vector3D pos(R[0], R[1], R[2]);
f0 = m_bField.getField(pos);
f0 = getField(cache, pos);
} else {
f0 = cache.field;
}
Expand Down Expand Up @@ -491,7 +520,7 @@ class AtlasStepper
//
if (!Helix) {
const Vector3D pos(R[0] + A1 * S4, R[1] + B1 * S4, R[2] + C1 * S4);
f = m_bField.getField(pos);
f = getField(cache, pos);
} else {
f = f0;
}
Expand All @@ -511,7 +540,7 @@ class AtlasStepper
//
if (!Helix) {
const Vector3D pos(R[0] + h * A4, R[1] + h * B4, R[2] + h * C4);
f = m_bField.getField(pos);
f = getField(cache, pos);
} else {
f = f0;
}
Expand Down
Loading

0 comments on commit 95dd052

Please sign in to comment.