LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PhotonVoxels.cxx
Go to the documentation of this file.
1 
7 // library header
10 
11 // C++ standard libraries
12 #include <algorithm> // std::min(), std::max()
13 #include <cmath> // std::abs(), std::floor()
14 #include <stdexcept> // std::runtime_error
15 #include <string>
16 
17 namespace sim {
18 
19  //----------------------------------------------------------------------------
20  // PhotonVoxelDef class
21  //----------------------------------------------------------------------------
23  double xMax,
24  int xN,
25  double yMin,
26  double yMax,
27  int yN,
28  double zMin,
29  double zMax,
30  int zN)
31  : fLowerCorner(xMin, yMin, zMin)
32  , fUpperCorner(xMax, yMax, zMax)
33  , fxSteps(xN)
34  , fySteps(yN)
35  , fzSteps(zN)
36  {}
37 
38  //----------------------------------------------------------------------------
39  std::array<unsigned int, 3U> PhotonVoxelDef::GetSteps() const
40  {
41  // BUG the double brace syntax is required to work around clang bug 21629
42  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
43  return {{fxSteps, fySteps, fzSteps}};
44  }
45 
46  //----------------------------------------------------------------------------
48  {
49  return ((GetRegionUpperCorner() == right.GetRegionUpperCorner()) &&
51  (GetSteps() == right.GetSteps()));
52  }
53 
54  //----------------------------------------------------------------------------
55  unsigned int PhotonVoxelDef::GetNVoxels() const
56  {
57  return fxSteps * fySteps * fzSteps;
58  }
59 
60  //----------------------------------------------------------------------------
61  int PhotonVoxelDef::GetVoxelID(double const* Position) const
62  {
63  return GetVoxelIDImpl(geo::vect::makeFromCoords<geo::Point_t>(Position));
64  }
65 
66  //----------------------------------------------------------------------------
67  std::optional<std::array<sim::PhotonVoxelDef::NeiInfo, 8U>>
69  {
70  if (!isInside(v)) return {};
71 
72  std::array<sim::PhotonVoxelDef::NeiInfo, 8U> ret;
73 
74  // Position in voxel coordinates including floating point part
75  auto const rStepD = GetVoxelStepCoordsUnchecked(v);
76 
77  // The neighbours are the 8 corners of a cube around this point
78  std::size_t iNeigh = 0U;
79  for (int dx : {0, 1}) {
80  for (int dy : {0, 1}) {
81  for (int dz : {0, 1}) {
82  // The full 3D step
83  const int dr[3] = {dx, dy, dz};
84 
85  // The integer-only position of the current corner
86  int rStepI[3];
87  for (int d = 0; d < 3; ++d) {
88  // Round down to get the "lower left" corner
89  rStepI[d] = int(rStepD[d]);
90  // Ensure we'll stay in-bounds
91  rStepI[d] = std::max(0, rStepI[d]);
92  rStepI[d] = std::min(rStepI[d], int(GetSteps()[d]) - 2);
93  // Adjust to the corner we're actually considering
94  rStepI[d] += dr[d];
95  }
96 
97  double w = 1;
98  for (int d = 0; d < 3; ++d) {
99  // These expressions will interpolate when between the 8 corners,
100  // and extrapolate in the half-voxel space around the edges.
101  if (dr[d] == 0)
102  w *= 1 + rStepI[d] - rStepD[d];
103  else
104  w *= 1 - rStepI[d] + rStepD[d];
105  }
106 
107  const int id = (rStepI[0] + rStepI[1] * (fxSteps) + rStepI[2] * (fxSteps * fySteps));
108 
109  ret[iNeigh++] = {id, w};
110  }
111  }
112  }
113 
114  // Sanity check the weights sum to 1
115  double wSum = 0;
116  for (const NeiInfo& n : ret)
117  wSum += n.weight;
118  if (std::abs(wSum - 1) > 1e-3) {
119  std::string msg = "PhotonVoxelDef::GetNeighboringVoxelIDs():"
120  " Weights sum to " +
121  std::to_string(wSum) +
122  " (should be 1)."
123  " Weights are:";
124  for (const NeiInfo& n : ret) {
125  msg += ' ';
126  msg += std::to_string(n.weight);
127  }
128  throw std::runtime_error(msg);
129  }
130  return {ret};
131  }
132 
133  //----------------------------------------------------------------------------
135  {
136  // float TempID = (float) ID;
137 
138  // Decompose ID into steps in each direction
139  int xStep = ID % fxSteps;
140  int yStep = ((ID - xStep) / fxSteps) % fySteps;
141  int zStep = ((ID - xStep - (yStep * fxSteps)) / (fySteps * fxSteps)) % fzSteps;
142 
143  auto const VoxelSize = GetVoxelSize<geo::Vector_t>();
144 
145  double const xMin = VoxelSize.X() * (xStep) + fLowerCorner.X();
146  double const xMax = VoxelSize.X() * (xStep + 1) + fLowerCorner.X();
147  double const yMin = VoxelSize.Y() * (yStep) + fLowerCorner.Y();
148  double const yMax = VoxelSize.Y() * (yStep + 1) + fLowerCorner.Y();
149  double const zMin = VoxelSize.Z() * (zStep) + fLowerCorner.Z();
150  double const zMax = VoxelSize.Z() * (zStep + 1) + fLowerCorner.Z();
151 
152  return PhotonVoxel(xMin, xMax, yMin, yMax, zMin, zMax);
153  }
154 
155  //----------------------------------------------------------------------------
157  {
158  return ((ID >= 0) && (static_cast<unsigned int>(ID) < GetNVoxels()));
159  }
160 
161  std::array<int, 3U> PhotonVoxelDef::GetVoxelCoords(int ID) const
162  {
163  std::array<int, 3U> ReturnVector;
164  ReturnVector[0] = ID % fxSteps;
165  ReturnVector[1] = ((ID - ReturnVector[0]) / fxSteps) % fySteps;
166  ReturnVector[2] =
167  ((ID - ReturnVector[0] - (ReturnVector[1] * fxSteps)) / (fySteps * fxSteps)) % fzSteps;
168  return ReturnVector;
169  }
170 
171  //----------------------------------------------------------------------------
172  std::array<double, 3U> PhotonVoxelDef::GetVoxelStepCoordsUnchecked(geo::Point_t const& p) const
173  {
174 
175  auto const span = fUpperCorner - fLowerCorner;
176  auto const relPos = p - fLowerCorner;
177 
178  // BUG the double brace syntax is required to work around clang bug 21629
179  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
180  return {{(relPos.X() / span.X()) * fxSteps,
181  (relPos.Y() / span.Y()) * fySteps,
182  (relPos.Z() / span.Z()) * fzSteps}};
183  } // PhotonVoxelDef::GetVoxelStepCoordsUnchecked()
184 
185  //----------------------------------------------------------------------------
187  {
188  if (!isInside(p)) return -1;
189 
190  auto const stepCoords = GetVoxelStepCoordsUnchecked(p);
191 
192  // figure out how many steps this point is in the x,y,z directions;
193  // `p` is guaranteed to be in the mapped volume by the previous check
194  int xStep = static_cast<int>(stepCoords[0]);
195  int yStep = static_cast<int>(stepCoords[1]);
196  int zStep = static_cast<int>(stepCoords[2]);
197 
198  // if within bounds, generate the voxel ID
199  return (xStep + yStep * (fxSteps) + zStep * (fxSteps * fySteps));
200  }
201 
202  //----------------------------------------------------------------------------
204  geo::Point_t const& lower,
205  geo::Point_t const& upper)
206  {
207  return isInsideRange(point.X(), lower.X(), upper.X()) &&
208  isInsideRange(point.Y(), lower.Y(), upper.Y()) &&
209  isInsideRange(point.Z(), lower.Z(), upper.Z());
210  }
211 
212  bool PhotonVoxelDef::isInsideRange(double value, double lower, double upper)
213  {
214 
215  return (value >= lower) && (value < upper);
216 
217  } // PhotonVoxelDef::isInsideRange()
218 
219  //----------------------------------------------------------------------------
220 
221  std::ostream& operator<<(std::ostream& out, sim::PhotonVoxelDef const& voxelDef)
222  {
223 
224  auto const& lower = voxelDef.GetRegionLowerCorner();
225  auto const& upper = voxelDef.GetRegionUpperCorner();
226  auto const& steps = voxelDef.GetSteps();
227  auto const& stepSize = voxelDef.GetVoxelSize();
228 
229  out << "Volume " << voxelDef.GetVolumeSize() << " cm^3 split in " << voxelDef.GetNVoxels()
230  << " voxels:"
231  << "\n - x axis: [ " << lower.X() << " ; " << upper.X() << " ] split in " << steps[0]
232  << "x " << stepSize.X() << " cm steps"
233  << "\n - y axis: [ " << lower.Y() << " ; " << upper.Y() << " ] split in " << steps[1]
234  << "x " << stepSize.Y() << " cm steps"
235  << "\n - z axis: [ " << lower.Z() << " ; " << upper.Z() << " ] split in " << steps[2]
236  << "x " << stepSize.Z() << " cm steps"
237  << "\n";
238 
239  return out;
240  } // operator<< (sim::PhotonVoxelDef)
241 
242  //----------------------------------------------------------------------------
243 
244 } // namespace sim
Definitions of voxel data structures.
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
unsigned int fySteps
Definition: PhotonVoxels.h:65
PhotonVoxelDef()=default
Vector GetVoxelSize() const
Returns a vector describing the span of a single voxel in x, y an z [cm].
Definition: PhotonVoxels.h:208
int GetVoxelIDImpl(geo::Point_t const &p) const
bool operator==(const PhotonVoxelDef &rhs) const
Representation of a region of space diced into voxels.
Definition: PhotonVoxels.h:58
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::array< double, 3U > GetVoxelStepCoordsUnchecked(geo::Point_t const &p) const
Returns the coordinates of the cvoxel containing p in step units.
int GetVoxelID(Point const &p) const
Returns the ID of the voxel containing p, or -1 if none.
Definition: PhotonVoxels.h:217
std::array< int, 3U > GetVoxelCoords(int ID) const
bool isInside(geo::Point_t const &p) const
Returns whether point p is inside the region (upper border excluded).
Definition: PhotonVoxels.h:137
static bool isInsideRange(double value, double lower, double upper)
unsigned int fzSteps
Definition: PhotonVoxels.h:66
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
decltype(auto) GetRegionUpperCorner() const
Returns the volume vertex (type Point) with the highest coordinates.
Utilities to extend the interface of geometry vectors.
Float_t d
Definition: plot.C:235
Monte Carlo Simulation.
bool IsLegalVoxelID(int) const
double value
Definition: spectrum.C:18
unsigned int fxSteps
Definition: PhotonVoxels.h:64
geo::Point_t fLowerCorner
Definition: PhotonVoxels.h:62
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
Representation of a single small volume (voxel).
Definition: PhotonVoxels.h:20
geo::Point_t fUpperCorner
Definition: PhotonVoxels.h:63
span(IterB &&b, IterE &&e, Adaptor &&adaptor) -> span< decltype(adaptor(std::forward< IterB >(b))), decltype(adaptor(std::forward< IterE >(e))) >
unsigned int GetNVoxels() const
Returns the total number of voxels in the volume.
std::array< unsigned int, 3U > GetSteps() const
Returns the number of voxels along each of the three dimensions.
Char_t n[5]
Vector GetVolumeSize() const
Returns a vector describing the full span in x, y an z [cm].
Definition: PhotonVoxels.h:97
Float_t e
Definition: plot.C:35
std::ostream & operator<<(std::ostream &output, const LArVoxelData &data)
static bool isInsideVolume(geo::Point_t const &point, geo::Point_t const &lower, geo::Point_t const &upper)
Float_t w
Definition: plot.C:20
std::optional< std::array< NeiInfo, 8U > > GetNeighboringVoxelIDsImpl(geo::Point_t const &v) const
decltype(auto) GetRegionLowerCorner() const
Returns the volume vertex (type Point) with the lowest coordinates.
PhotonVoxel GetPhotonVoxel(int ID) const