Skip to content

Commit

Permalink
feat: Add angles of incidence in RootMeasurementWriter (#2493)
Browse files Browse the repository at this point in the history 
For measurement calibrations, the angles of incidence on a surface are very useful. Previously, the MeasurementReader only saved the truth particle direction in the global frame, so this PR adds the angle of incidence.

To validate, here are the angles of incidence for 1-pixel clusters on the Generic detector barrel:

![f_phi](https://github.com/acts-project/acts/assets/2743680/c2e2b0b6-2d8a-4388-a63b-61aa2490f05e)

![f_theta](https://github.com/acts-project/acts/assets/2743680/bac7f163-8980-4523-a9cd-2fe018c6db85)

(presumably, the phi distribution is asymetric because of the Lorentz angle)

closes #2484
  • Loading branch information
@gagnonlg
gagnonlg authored Oct 2, 2023
1 parent 363d111 commit dd3348c
Show file tree
Hide file tree
Showing 5 changed files with 116 additions and 7 deletions.
16 changes: 16 additions & 0 deletions Core/include/Acts/Utilities/VectorHelpers.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,5 +216,21 @@ inline auto makeVector4(const Eigen::MatrixBase<vector3_t>& vec3,
return vec4;
}

/// Calculate the incident angles of a vector with in a given reference frame
/// @tparam Derived Eigen derived concrete type
/// @param direction The crossing direction in the global frame
/// @param globalToLocal Rotation from global to local frame
/// @return The angles of incidence in the two normal planes
inline std::pair<double, double> incidentAngles(
const Acts::Vector3& direction,
const Acts::RotationMatrix3& globalToLocal) {
Acts::Vector3 trfDir = globalToLocal * direction;
// The angles are defined with respect to the measurement axis
// i.e. "head-on" == pi/2, parallel = 0
double phi = std::atan2(trfDir[2], trfDir[0]);
double theta = std::atan2(trfDir[2], trfDir[1]);
return {phi, theta};
}

} // namespace VectorHelpers
} // namespace Acts
10 changes: 9 additions & 1 deletion Examples/Io/Root/include/ActsExamples/Io/Root/RootMeasurementWriter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,8 @@ class RootMeasurementWriter final : public WriterT<MeasurementContainer> {
float trueGx = 0.;
float trueGy = 0.;
float trueGz = 0.;
float incidentPhi = 0.;
float incidentTheta = 0.;

/// Reconstruction information
float recBound[Acts::eBoundSize] = {};
Expand Down Expand Up @@ -130,6 +132,8 @@ class RootMeasurementWriter final : public WriterT<MeasurementContainer> {
tree->Branch("true_x", &trueGx);
tree->Branch("true_y", &trueGy);
tree->Branch("true_z", &trueGz);
tree->Branch("true_incident_phi", &incidentPhi);
tree->Branch("true_incident_theta", &incidentTheta);
}

/// Constructor from GeometryIdentifier
Expand Down Expand Up @@ -198,7 +202,8 @@ class RootMeasurementWriter final : public WriterT<MeasurementContainer> {
/// @param xt The true 4D global position
/// @param dir The true particle direction
void fillTruthParameters(const Acts::Vector2& lp, const Acts::Vector4& xt,
const Acts::Vector3& dir) {
const Acts::Vector3& dir,
const std::pair<double, double> angles) {
trueBound[Acts::eBoundLoc0] = lp[Acts::eBoundLoc0];
trueBound[Acts::eBoundLoc1] = lp[Acts::eBoundLoc1];
trueBound[Acts::eBoundPhi] = Acts::VectorHelpers::phi(dir);
Expand All @@ -208,6 +213,9 @@ class RootMeasurementWriter final : public WriterT<MeasurementContainer> {
trueGx = xt[Acts::ePos0];
trueGy = xt[Acts::ePos1];
trueGz = xt[Acts::ePos2];

incidentPhi = angles.first;
incidentTheta = angles.second;
}

/// Convenience function to fill bound parameters
Expand Down
9 changes: 8 additions & 1 deletion Examples/Io/Root/src/RootMeasurementWriter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,14 @@ ActsExamples::ProcessCode ActsExamples::RootMeasurementWriter::writeT(
// Use average truth in the case of multiple contributing sim hits
auto [local, pos4, dir] = averageSimHits(ctx.geoContext, surface,
simHits, indices, logger());
dTree->fillTruthParameters(local, pos4, dir);
Acts::RotationMatrix3 rot =
surface
.referenceFrame(ctx.geoContext, pos4.segment<3>(Acts::ePos0),
dir)
.inverse();
std::pair<double, double> angles =
Acts::VectorHelpers::incidentAngles(dir, rot);
dTree->fillTruthParameters(local, pos4, dir, angles);
dTree->fillBoundMeasurement(m);
if (not clusters.empty()) {
const auto& c = clusters[hitIdx];
Expand Down
10 changes: 5 additions & 5 deletions Examples/Python/tests/root_file_hashes.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,12 +43,12 @@ test_material_mapping__material-map_tracks.root: 4e1c866038f0c06b099aa74fd01c3d8
test_material_mapping__propagation-material.root: 646b8e2bbacec40d0bc4132236f9ab3f03b088e656e6e9b80c47ae03eaf6eab5
test_volume_material_mapping__material-map-volume_tracks.root: b95561a6247df9e3599a997daa6c1d76461e58f83059b82f2ec27229c9b35e6c
test_volume_material_mapping__propagation-volume-material.root: b7597dada372d1b4aaec2c4fc3c0db830ce147ecf515c367ac6ba8ffc2708302
test_digitization_example[smeared]__measurements.root: 8dea40cccf3acaab2792f8994701e8baafb8b334277a7fdf2c7ff45fdbae776b
test_digitization_example[geometric]__measurements.root: 8b45ed4abc5d7d6dc87de1921df68981465d5917b1d5d458fb5d105700f13402
test_digitization_example[smeared]__measurements.root: 0888b7bb4da63d15be18c14c5a419c8813f9d76731da94751247dc4a6772c895
test_digitization_example[geometric]__measurements.root: bade4e094a950f4ecf523b192294e4336da14b4ba6dc5164160cc8e3308cdcd2
test_digitization_example_input[smeared]__particles.root: 8549ba6e20338004ab8ba299fc65e1ee5071985b46df8f77f887cb6fef56a8ec
test_digitization_example_input[smeared]__measurements.root: df0e31749f64f52b7917e44c47fa1eda9e72ef53c44c0e09f3c75b7143b99d24
test_digitization_example_input[smeared]__measurements.root: dd6d1af353a676d28ac398d015a20db736dccec0fb104ea909d81c9156072e18
test_digitization_example_input[geometric]__particles.root: 8549ba6e20338004ab8ba299fc65e1ee5071985b46df8f77f887cb6fef56a8ec
test_digitization_example_input[geometric]__measurements.root: b4dba5bb5e60d4540c1d07fb297b2add19d779dfcd29737c8d44f208ff0c0522
test_digitization_example_input[geometric]__measurements.root: 56a64b3df65fcd04d677fac6b6ee3f1978fe6d7777599bb918d40e514ef45a1f
test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: 1786b1f1e90a8fe2ad6f534072506bacf8dbcf1a7e25c283a5ab261a319479fb
test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: 359be063e44247439490d4ae99fa3c9c240e8c0777b27421b26b63f62c85e2f0
test_ckf_tracks_example[generic-full_seeding]__performance_seeding_trees.root: 0e0676ffafdb27112fbda50d1cf627859fa745760f98073261dcf6db3f2f991e
Expand Down Expand Up @@ -78,7 +78,7 @@ test_root_prop_step_writer[kwargsConstructor]__prop_steps.root: 53adaabc7856c888
test_root_particle_writer[configPosConstructor]__particles.root: 8ca3987523360e4cd52f9bc138e82230996baa105c77ddb2cde775733d8acc4c
test_root_particle_writer[configKwConstructor]__particles.root: 8ca3987523360e4cd52f9bc138e82230996baa105c77ddb2cde775733d8acc4c
test_root_particle_writer[kwargsConstructor]__particles.root: 8ca3987523360e4cd52f9bc138e82230996baa105c77ddb2cde775733d8acc4c
test_root_meas_writer__meas.root: 283ccfbbe10bf612f499a9d5c4ca7ac1e1dd4fad7f1629342813ecd36ec06acd
test_root_meas_writer__meas.root: 1e4138b38f7dfe612faef042aed0e9de872ffdaf89c1a386052d937276c35c30
test_root_simhits_writer[configPosConstructor]__meas.root: 26b51bbd97cd5f915628e150aa286a766a30cecdd511f8476e0c6b2ea62227fa
test_root_simhits_writer[configKwConstructor]__meas.root: 26b51bbd97cd5f915628e150aa286a766a30cecdd511f8476e0c6b2ea62227fa
test_root_simhits_writer[kwargsConstructor]__meas.root: 26b51bbd97cd5f915628e150aa286a766a30cecdd511f8476e0c6b2ea62227fa
Expand Down
78 changes: 78 additions & 0 deletions Tests/UnitTests/Core/Utilities/HelpersTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,84 @@ BOOST_AUTO_TEST_CASE(safeInverse) {
}
}

BOOST_AUTO_TEST_CASE(incidentAnglesTest) {
RotationMatrix3 ref = RotationMatrix3::Identity();

// Right angle in both planes
for (size_t i = 0; i < 3; i++) {
Vector3 dir = Vector3::Zero();
dir[i] = 1;
std::pair<double, double> angles = incidentAngles(dir, ref);
double expect = (i < 2) ? 0 : M_PI_2;
CHECK_CLOSE_ABS(angles.first, expect,
std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(angles.second, expect,
std::numeric_limits<ActsScalar>::epsilon());
}

// 45 degree on both axes
{
Vector3 dir = Vector3({1, 1, 1}).normalized();
auto [a0, a1] = incidentAngles(dir, ref);
CHECK_CLOSE_ABS(a0, M_PI_4, std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(a1, M_PI_4, std::numeric_limits<ActsScalar>::epsilon());
}

// 45 degree on first axis
{
Vector3 dir = Vector3({1, 0, 1}).normalized();
auto [a0, a1] = incidentAngles(dir, ref);
CHECK_CLOSE_ABS(a0, M_PI_4, std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(a1, M_PI_2, std::numeric_limits<ActsScalar>::epsilon());
}

// 45 degree on second axis
{
Vector3 dir = Vector3({0, 1, 1}).normalized();
auto [a0, a1] = incidentAngles(dir, ref);
CHECK_CLOSE_ABS(a0, M_PI_2, std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(a1, M_PI_4, std::numeric_limits<ActsScalar>::epsilon());
}

// Reverse crossing
{
Vector3 dir = {0, 0, -1};
auto [a0, a1] = incidentAngles(dir, ref);
CHECK_CLOSE_ABS(a0, -M_PI_2, std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(a1, -M_PI_2, std::numeric_limits<ActsScalar>::epsilon());
}

// 45 degree but different quadrant
{
Vector3 dir = {-1, -1, 1};
auto [a0, a1] = incidentAngles(dir, ref);
CHECK_CLOSE_ABS(a0, 3 * M_PI_4, std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(a1, 3 * M_PI_4, std::numeric_limits<ActsScalar>::epsilon());
}

// 45 degree but different quadrant & other side
{
Vector3 dir = {-1, -1, -1};
auto [a0, a1] = incidentAngles(dir, ref);
CHECK_CLOSE_ABS(a0, -3 * M_PI_4,
std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(a1, -3 * M_PI_4,
std::numeric_limits<ActsScalar>::epsilon());
}

// Rotate the reference frame instead
{
double s45 = std::sin(M_PI_4);
double c45 = std::cos(M_PI_4);
RotationMatrix3 ref45;
ref45 << c45, 0, s45, 0, 1, 0, -s45, 0, c45;
Vector3 dir = {0, 0, 1};
auto [a0, a1] = incidentAngles(dir, ref45);
CHECK_CLOSE_ABS(a0, M_PI_4, std::numeric_limits<ActsScalar>::epsilon());
CHECK_CLOSE_ABS(a1, M_PI_2, std::numeric_limits<ActsScalar>::epsilon());
}
}

BOOST_AUTO_TEST_SUITE_END()

} // namespace Test
Expand Down

0 comments on commit dd3348c

Please sign in to comment.