Acts/Utilities/BinAdjustment.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Surfaces/CylinderBounds.hpp"
0013 #include "Acts/Surfaces/RadialBounds.hpp"
0014 #include "Acts/Surfaces/RectangleBounds.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0017 #include "Acts/Utilities/AxisDefinitions.hpp"
0018 #include "Acts/Utilities/BinUtility.hpp"
0019
0020 #include <stdexcept>
0021
0022 namespace Acts {
0023
0024 /// @brief adjust the BinUtility bu to the dimensions of radial bounds
0025 ///
0026 /// @param bu BinUtility at source
0027 /// @param rBounds the Radial bounds to adjust to
0028 /// @param transform Transform for the adjusted @c BinUtility
0029 ///
0030 /// @return new updated BinUtiltiy
0031 static inline BinUtility adjustBinUtility(const BinUtility& bu,
0032                                           const RadialBounds& rBounds,
0033                                           const Transform3& transform) {
0034   // Default constructor
0035   BinUtility uBinUtil(transform);
0036
0037   // The parameters from the cylinder bounds
0038   double minR = rBounds.get(RadialBounds::eMinR);
0039   double maxR = rBounds.get(RadialBounds::eMaxR);
0040   double minPhi = rBounds.get(RadialBounds::eAveragePhi) -
0041                   rBounds.get(RadialBounds::eHalfPhiSector);
0042   double maxPhi = rBounds.get(RadialBounds::eAveragePhi) +
0043                   rBounds.get(RadialBounds::eHalfPhiSector);
0044   // Retrieve the binning data
0045   const std::vector<BinningData>& bData = bu.binningData();
0046   // Loop over the binning data and adjust the dimensions
0047   for (auto& bd : bData) {
0048     // The binning value
0049     AxisDirection bval = bd.binvalue;
0050     // Throw exceptions is stuff doesn't make sense:
0051     // - not the right binning value
0052     // - not equidistant
0053     if (bd.type == arbitrary) {
0054       throw std::invalid_argument("Arbitrary binning can not be adjusted.");
0055     } else if (bval != AxisDirection::AxisR && bval != AxisDirection::AxisPhi) {
0056       throw std::invalid_argument("Disc binning must be: phi, r");
0057     }
0058     float min = 0., max = 0.;
0059     // Perform the value adjustment
0060     if (bval == AxisDirection::AxisPhi) {
0061       min = minPhi;
0062       max = maxPhi;
0063     } else {
0064       min = minR;
0065       max = maxR;
0066     }
0067     // Create the updated BinningData
0068     BinningData uBinData(bd.option, bval, bd.bins(), min, max);
0069     uBinUtil += BinUtility(uBinData);
0070   }
0071   return uBinUtil;
0072 }
0073
0074 /// @brief adjust the BinUtility bu to the dimensions of cylinder bounds
0075 ///
0076 /// @param bu BinUtility at source
0077 /// @param cBounds the Cylinder bounds to adjust to
0078 /// @param transform Transform for the adjusted @c BinUtility
0079 ///
0080 /// @return new updated BinUtiltiy
0081 static inline BinUtility adjustBinUtility(const BinUtility& bu,
0082                                           const CylinderBounds& cBounds,
0083                                           const Transform3& transform) {
0084   // Default constructor
0085   BinUtility uBinUtil(transform);
0086
0087   // The parameters from the cylinder bounds
0088   double cR = cBounds.get(CylinderBounds::eR);
0089   double cHz = cBounds.get(CylinderBounds::eHalfLengthZ);
0090   double avgPhi = cBounds.get(CylinderBounds::eAveragePhi);
0091   double halfPhi = cBounds.get(CylinderBounds::eHalfPhiSector);
0092   double minPhi = avgPhi - halfPhi;
0093   double maxPhi = avgPhi + halfPhi;
0094
0095   // Retrieve the binning data
0096   const std::vector<BinningData>& bData = bu.binningData();
0097   // Loop over the binning data and adjust the dimensions
0098   for (auto& bd : bData) {
0099     // The binning value
0100     AxisDirection bval = bd.binvalue;
0101     // Throw exceptions if stuff doesn't make sense:
0102     // - not the right binning value
0103     // - not equidistant
0104     if (bd.type == arbitrary) {
0105       throw std::invalid_argument("Arbitrary binning can not be adjusted.");
0106     } else if (bval != AxisDirection::AxisRPhi &&
0107                bval != AxisDirection::AxisPhi && bval != AxisDirection::AxisZ) {
0108       throw std::invalid_argument("Cylinder binning must be: rphi, phi, z");
0109     }
0110     float min = 0., max = 0.;
0111     // Perform the value adjustment
0112     if (bval == AxisDirection::AxisPhi) {
0113       min = minPhi;
0114       max = maxPhi;
0115     } else if (bval == AxisDirection::AxisRPhi) {
0116       min = cR * minPhi;
0117       max = cR * maxPhi;
0118     } else {
0119       min = -cHz;
0120       max = cHz;
0121     }
0122     // Create the updated BinningData
0123     BinningData uBinData(bd.option, bval, bd.bins(), min, max);
0124     uBinUtil += BinUtility(uBinData);
0125   }
0126   return uBinUtil;
0127 }
0128
0129 /// @brief adjust the BinUtility bu to the dimensions of plane bounds
0130 ///
0131 /// @param bu BinUtility at source
0132 /// @param pBounds the Rectangle bounds to adjust to
0133 /// @param transform Transform for the adjusted @c BinUtility
0134 ///
0135 /// @return new updated BinUtiltiy
0136 static inline BinUtility adjustBinUtility(const BinUtility& bu,
0137                                           const RectangleBounds& pBounds,
0138                                           const Transform3& transform) {
0139   // Default constructor
0140   BinUtility uBinUtil(transform);
0141
0142   // The parameters from the cylinder bounds
0143   double minX = pBounds.get(RectangleBounds::eMinX);
0144   double minY = pBounds.get(RectangleBounds::eMinY);
0145   double maxX = pBounds.get(RectangleBounds::eMaxX);
0146   double maxY = pBounds.get(RectangleBounds::eMaxY);
0147
0148   // Retrieve the binning data
0149   const std::vector<BinningData>& bData = bu.binningData();
0150   // Loop over the binning data and adjust the dimensions
0151   for (auto& bd : bData) {
0152     // The binning value
0153     AxisDirection bval = bd.binvalue;
0154     // Throw exceptions if stuff doesn't make sense:
0155     // - not the right binning value
0156     // - not equidistant
0157     if (bd.type == arbitrary) {
0158       throw std::invalid_argument("Arbitrary binning can not be adjusted.");
0159     } else if (bval != AxisDirection::AxisX && bval != AxisDirection::AxisY) {
0160       throw std::invalid_argument("Rectangle binning must be: x, y. ");
0161     }
0162     float min = 0., max = 0.;
0163     // Perform the value adjustment
0164     if (bval == AxisDirection::AxisX) {
0165       min = minX;
0166       max = maxX;
0167     } else {
0168       min = minY;
0169       max = maxY;
0170     }
0171     // Create the updated BinningData
0172     BinningData uBinData(bd.option, bval, bd.bins(), min, max);
0173     uBinUtil += BinUtility(uBinData);
0174   }
0175
0176   return uBinUtil;
0177 }
0178
0179 /// @brief adjust the BinUtility bu to the dimensions of plane bounds
0180 ///
0181 /// @param bu BinUtility at source
0182 /// @param pBounds the Trapezoid bounds to adjust to
0183 /// @param transform Transform for the adjusted @c BinUtility
0184 ///
0185 /// @return new updated BinUtiltiy
0186 static inline BinUtility adjustBinUtility(const BinUtility& bu,
0187                                           const TrapezoidBounds& pBounds,
0188                                           const Transform3& transform) {
0189   // Default constructor
0190   BinUtility uBinUtil(transform);
0191
0192   // The parameters from the cylinder bounds
0193
0194   double halfX = std::max(pBounds.get(Acts::TrapezoidBounds::eHalfLengthXnegY),
0195                           pBounds.get(Acts::TrapezoidBounds::eHalfLengthXposY));
0196   double halfY = pBounds.get(Acts::TrapezoidBounds::eHalfLengthY);
0197
0198   // Retrieve the binning data
0199   const std::vector<BinningData>& bData = bu.binningData();
0200   // Loop over the binning data and adjust the dimensions
0201   for (auto& bd : bData) {
0202     // The binning value
0203     AxisDirection bval = bd.binvalue;
0204     // Throw exceptions if stuff doesn't make sense:
0205     // - not the right binning value
0206     // - not equidistant
0207     if (bd.type == arbitrary) {
0208       throw std::invalid_argument("Arbitrary binning can not be adjusted.");
0209     } else if (bval != AxisDirection::AxisX && bval != AxisDirection::AxisY) {
0210       throw std::invalid_argument("Rectangle binning must be: x, y. ");
0211     }
0212     float min = 0., max = 0.;
0213     // Perform the value adjustment
0214     if (bval == AxisDirection::AxisX) {
0215       min = -1 * halfX;
0216       max = halfX;
0217     } else {
0218       min = -1 * halfY;
0219       max = halfY;
0220     }
0221     // Create the updated BinningData
0222     BinningData uBinData(bd.option, bval, bd.bins(), min, max);
0223     uBinUtil += BinUtility(uBinData);
0224   }
0225
0226   return uBinUtil;
0227 }
0228
0229 /// @brief adjust the BinUtility bu to a surface
0230 ///
0231 /// @param bu BinUtility at source
0232 /// @param surface Surface to which the adjustment is being done
0233 /// @param gctx Geometry context to get the surfaces transform
0234 ///
0235 /// @return new updated BinUtiltiy
0236 static inline BinUtility adjustBinUtility(const BinUtility& bu,
0237                                           const Surface& surface,
0238                                           const GeometryContext& gctx) {
0239   if (auto b = dynamic_cast<const CylinderBounds*>(&(surface.bounds()));
0240       b != nullptr) {
0241     return adjustBinUtility(bu, *b, surface.transform(gctx));
0242   }
0243   if (auto b = dynamic_cast<const RadialBounds*>(&(surface.bounds()));
0244       b != nullptr) {
0245     return adjustBinUtility(bu, *b, surface.transform(gctx));
0246   }
0247   if (surface.type() == Surface::Plane) {
0248     if (auto b = dynamic_cast<const RectangleBounds*>(&(surface.bounds()));
0249         b != nullptr) {
0250       return adjustBinUtility(bu, *b, surface.transform(gctx));
0251     }
0252     if (auto b = dynamic_cast<const TrapezoidBounds*>(&(surface.bounds()));
0253         b != nullptr) {
0254       return adjustBinUtility(bu, *b, surface.transform(gctx));
0255     }
0256   }
0257
0258   std::stringstream ss;
0259   ss << surface.toStream({});
0260   throw std::invalid_argument(
0261       "Bin adjustment not implemented for this surface yet:\n" + ss.str());
0262 }
0263
0264 }  // namespace Acts